C++ CData::get_compact_vector方法代码示例

void CParam::S2_add(Uniform &randUnif,CData &Data) {
  int n_needtoupdate = 0;
  for (int i_faulty=1; i_faulty<=Data.n_faulty; i_faulty++){
    int i_original = Data.Faulty2Original[i_faulty-1];
    ColumnVector item_by_bal;
    ColumnVector s_i = S_Mat.column(i_faulty);
    ColumnVector item_by_rnorm = Data.get_item_by_norm_indicator(s_i,item_by_bal);

    //Generate from normal distribution
    if ( item_by_rnorm.sum() >= 1 ) { // if no random number, other values by balanc edits remain same
    n_needtoupdate++;
    ColumnVector mu_z_i = Mu.column(z_in(i_original)) ;
    ColumnVector tilde_y_i = Data.log_D_Observed.row(i_original).t();
    ColumnVector s_1_compact = Data.get_compact_vector(item_by_rnorm);
    ColumnVector Mu_1i = subvector(mu_z_i,s_1_compact);

    LowerTriangularMatrix LSigma_1i_i;
    ColumnVector y_q(n_var);
    double log_cond_norm_q = calculate_log_cond_norm(Data, i_original, item_by_rnorm, tilde_y_i, y_q, true, LSigma_1i_i, s_i); // MODIFIED 2015/02/16
    ColumnVector y_i = (Y_in.row(i_original)).t() ;
    // ColumnVector y_part_i = subvector(y_i,item_by_rnorm);

    // Put values from balance edits
    ColumnVector x_q = exp_ColumnVector(y_q) ;
    Data.set_balance_edit_values_for_x_q(s_i, x_q, item_by_bal); // CHANGED by Hang, 2014/12/29

    // double log_cond_norm_i = log_MVN_fn(y_part_i,Mu_1i,LSigma_1i_i);
    double log_cond_norm_i = calculate_log_cond_norm(Data, i_original, item_by_rnorm, tilde_y_i, y_q, false, LSigma_1i_i, s_i); // CHANGED 2015/01/27 , // MODIFIED 2015/02/16

    // Acceptance/Rejection
    if (Data.PassEdits(x_q)) {  // Check constraints
    y_q = log_ColumnVector(x_q) ;
    ColumnVector y_compact_q = Data.get_compact_vector(y_q);
    ColumnVector y_compact_i = Data.get_compact_vector(y_i);
    double log_full_norm_q = log_MVN_fn(y_compact_q,mu_z_i,LSIGMA_i[z_in(i_original)-1],logdet_and_more(z_in(i_original)));
    double log_full_norm_i = log_MVN_fn(y_compact_i,mu_z_i,LSIGMA_i[z_in(i_original)-1],logdet_and_more(z_in(i_original)));

    // Calculate acceptance ratio
    double logNum = log_full_norm_q - log_cond_norm_q;
    double logDen = log_full_norm_i - log_cond_norm_i;
    accept_rate(2) = exp( logNum - logDen );

    if (randUnif.Next() < accept_rate(2)){
      Y_in.row(i_original) = y_q.t();
      is_accept(2)++;
    }
    }
    }
  }
  is_accept(2) = is_accept(2) / n_needtoupdate;
}

double CParam::calculate_log_cond_norm(CData &Data, int i_original, ColumnVector &item_by_rnorm, ColumnVector &tilde_y_i, ColumnVector &y_q, bool is_q, LowerTriangularMatrix &LSigma_1_i, ColumnVector &s_q) { // MODIFIED 2015/02/16

  double log_cond_norm;

  if ( item_by_rnorm.sum() >= 1 ) {

    ColumnVector mu_z_i = Mu.column(z_in(i_original));
    ColumnVector s_1_compact = Data.get_compact_vector(item_by_rnorm); 
    ColumnVector Mu_1 = subvector(mu_z_i,s_1_compact);
    Matrix Sigma_1 = Submatrix_elem_2(SIGMA[z_in(i_original)-1],s_1_compact,s_1_compact);

    // ADDED 2015/01/27

    ColumnVector s_q_compact = Data.get_compact_vector(s_q) ; // MODIFIED 2015/02/16 
    ColumnVector VectorOne = s_q_compact ; VectorOne = 1 ; // MODIFIED 2015/02/16
    ColumnVector s_0_compact = VectorOne - s_q_compact ; // MODIFIED 2015/02/16
    int sum_s_0_comp = s_0_compact.sum() ;
    LowerTriangularMatrix LSigma_cond ;
    ColumnVector Mu_cond ;

    if ( sum_s_0_comp>0 ){

      ColumnVector Mu_0 = subvector(mu_z_i,s_0_compact); // (s_1_compact.sum()) vector
      Matrix Sigma_0 = Submatrix_elem_2(SIGMA[z_in(i_original)-1],s_0_compact,s_0_compact);
      Matrix Sigma_10 = Submatrix_elem_2(SIGMA[z_in(i_original)-1],s_1_compact,s_0_compact);
      ColumnVector y_tilde_compact = Data.get_compact_vector(tilde_y_i) ;
      ColumnVector y_tilde_0 = subvector(y_tilde_compact,s_0_compact) ;
      SymmetricMatrix Sigma_0_symm ; Sigma_0_symm << Sigma_0 ;

      LowerTriangularMatrix LSigma_0 = Cholesky(Sigma_0_symm) ;
      Mu_cond = Mu_1 + Sigma_10 * (LSigma_0.i()).t()*LSigma_0.i() * ( y_tilde_0-Mu_0 ) ;
      Matrix Sigma_cond = Sigma_1 - Sigma_10 * (LSigma_0.i()).t()*LSigma_0.i() * Sigma_10.t() ;
      SymmetricMatrix Sigma_cond_symm ; Sigma_cond_symm << Sigma_cond ;
      int sum_s_1_comp = s_1_compact.sum() ;
      DiagonalMatrix D(sum_s_1_comp) ; Matrix V(sum_s_1_comp,sum_s_1_comp) ;
      Jacobi(Sigma_cond_symm,D,V) ;
      int is_zero_exist = 0 ;
      for (int i_var=1; i_var<=sum_s_1_comp; i_var++){
    	if ( D(i_var) < 1e-9 ){
    	  D(i_var) = 1e-9 ;
    	  is_zero_exist = 1 ;
    	}
      } // for (int i_var=1; i_var<=sum_s_1_comp; i_var++)
      if ( is_zero_exist == 1 ){
    	Sigma_cond_symm << V * D * V.t() ;
    	if ( msg_level >= 1 ) {
		  Rprintf( "   Warning: When generating y_j from conditional normal(Mu_-j,Sigma_-j), Sigma_-j is non-positive definite because of computation precision. The eigenvalues D(j,j) smaller than 1e-9 is replaced with 1e-9, and let Sigma_-j = V D V.t().\n");
    	}
      } //
      LSigma_cond = Cholesky(Sigma_cond_symm);
      // y_part = rMVN_fn(Mu_cond,LSigma_cond);
      // log_cond_norm = log_MVN_fn(y_part,Mu_cond,LSigma_cond) ;

    } else {

      Mu_cond = Mu_1 ;
      SymmetricMatrix Sigma_1_symm = Submatrix_elem(SIGMA[z_in(i_original)-1],s_1_compact);
      LSigma_cond = Cholesky(Sigma_1_symm) ;
      // SymmetricMatrix Sigma_1_symm ; Sigma_1_symm << Sigma_1 ;
      // LowerTriangularMatrix LSigma_1 = Cholesky_Sigma_star_symm(Sigma_1_symm);
      // y_part = rMVN_fn(Mu_1,LSigma_1);
      // log_cond_norm = log_MVN_fn(y_part,Mu_1,LSigma_1) ;

    } // if ( sum_s_0_comp>0 ) else ...

    // ADDED 2015/01/26

    LowerTriangularMatrix LSigma_cond_i = LSigma_cond.i() ;
    // LowerTriangularMatrix LSigma_1 = Cholesky(Sigma_1);
    // LSigma_1_i = LSigma_1.i();

    ColumnVector y_part;
    if (is_q) {
      y_part = rMVN_fn(Mu_cond,LSigma_cond);
    } else {
      ColumnVector y_i = (Y_in.row(i_original)).t();
      y_part = subvector(y_i,item_by_rnorm);
    }
    
    log_cond_norm = log_MVN_fn(y_part,Mu_cond,LSigma_cond_i);

    if (is_q) {
      y_q = tilde_y_i;
      for ( int temp_j = 1,temp_count1 = 0; temp_j<=n_var; temp_j++ ){
        if ( item_by_rnorm(temp_j)==1 ){ 
          y_q(temp_j) = y_part(++temp_count1);
	}
      }
    } // if (is_q)

  } else {

    log_cond_norm = 0;
    if (is_q) { y_q = tilde_y_i;}

  } // if ( item_by_rnorm.sum() > = 1 ) else ..

  return log_cond_norm;
}