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

int CFeasibilityMap::count_x_out_fn(CData &Data,int i_tau, int i_original,int n_simul, Uniform &randUnif) {
	
  // double case2_count_out = 0;
	int case2_count_out = 0; // Changed by Hang on 5/16/2015
	
  ColumnVector s_i = tau_to_s_fn( i_tau, Data.n_var );   
  ColumnVector item_by_joint = Data.copy_non_balance_edit(s_i);
  ColumnVector tilde_y_i = Data.log_D_Observed.row(i_original).t();
  
	for (int i_simul=1; i_simul<=n_simul; i_simul++){
			//Generate from uniform distribution
			ColumnVector y_q = tilde_y_i;
			for ( int temp_j=1; temp_j<=Data.n_var; temp_j++ ){
				if ( item_by_joint(temp_j)==1 ){
					y_q(temp_j) = Data.logB_L(temp_j)+Data.logB_U_L(temp_j)*randUnif.Next(); 
				} 
			} 
	
			ColumnVector x_q = exp_ColumnVector(y_q) ;
      Data.update_full_x_for_balance_edit(x_q);
			// if (!Data.PassEdits(x_q)) { case2_count_out += 1.0;}
      if (!Data.PassEdits(x_q)) { case2_count_out += 1;}  // Changed by Hang on 5/16/2015
	} 
  if (case2_count_out ==0) {
    case2_count_out = 1;
  }
	return case2_count_out; // ADDED by Hang on 5/16/2015
}

void CFeasibilityMap::Simulate_logUnif_case2(int n_simul, Uniform &randUnif, CData &Data) {
  if (useMap) {return;}
  if (feasibleMap.nrows() == 0 || feasibleMap.maximum()==0) {
	  Rprintf( "Feasibility Map need to be set or computed first\n");

    return;
  }
  Data.logUnif_case2 =Matrix(Data.n_faulty,Data.n_tau); Data.logUnif_case2 = 0;
  for (int i_original=1; i_original<=Data.n_sample; i_original++){
    if (Data.is_case(i_original,2)) {
      int i_faulty = Data.Original2Faulty[i_original-1];
      ColumnVector tilde_y_i = Data.log_D_Observed.row(i_original).t();
      for (int i_tau=1; i_tau<=Data.n_tau; i_tau++){
        if ( feasibleMap(i_tau,i_faulty)==1 ){
          double case2_count_out = 0;
          ColumnVector s_i = tau_to_s_fn( i_tau, Data.n_var );
          ColumnVector item_by_joint = Data.copy_non_balance_edit(s_i);
          
    			for (int i_simul=1; i_simul<=n_simul; i_simul++){
    					//Generate from uniform distribution
    					ColumnVector y_q = tilde_y_i;
    					for ( int temp_j=1; temp_j<=Data.n_var; temp_j++ ){
    						if ( item_by_joint(temp_j)==1 ){
    							y_q(temp_j) = Data.logB_L(temp_j)+Data.logB_U_L(temp_j)*randUnif.Next(); 
    						} 
    					} 
    			
    					ColumnVector x_q = exp_ColumnVector(y_q) ;
              Data.update_full_x_for_balance_edit(x_q);
              if (!Data.PassEdits(x_q)) { case2_count_out += 1.0;} 
    			} 
    		
    			double Area = 1.0; 
    			for ( int temp_j=1; temp_j<=Data.n_var; temp_j++ ){
    				if ( item_by_joint(temp_j)==1 ){
    					Area = Area * Data.logB_U_L(temp_j) ; 
    				} 
    			} 
          Area = Area * case2_count_out / n_simul ;
          Data.set_logUnif_case2(i_original, i_tau, -log(Area));
     		} 
    	} 
    	
      if ( ((1.0*i_original/100)==(floor(1.0*i_original/100))) ){ 
		   Rprintf( "logUnif_y_tilde for i_sample= %d\n",i_original);
    	}
    }
  }
}

void CParam::S1(int iter, Uniform &randUnif, CData &Data, CFeasibilityMap &FM, int n_simul) {
  is_accept(1) = 0;
  for (int i_faulty=1; i_faulty<=Data.n_faulty; i_faulty++) {
    if (FM.useMap) {FM.pmm->iter = iter;}
    int i_original = Data.Faulty2Original[i_faulty-1];
    double g_option_q, g_mode_q;
    int what_type_move;
    ColumnVector s_i = S_Mat.column(i_faulty);
    int tau_q = FM.EvaluateMove(i_original, Data, s_i, iter, what_type_move, g_option_q, g_mode_q, true);
    // calculate g_mode_i & g_option_i
    ColumnVector s_q = CFeasibilityMap::tau_to_s_fn(tau_q, n_var);
    double g_option_i, g_mode_i;
    FM.EvaluateMove(i_original, Data, s_q, iter, what_type_move, g_option_i, g_mode_i, false);
    // calculate g_s_q & g_s_i
    double g_s_q = g_mode_q * g_option_q;
    double g_s_i = g_mode_i * g_option_i;

    // whether drawn by item_by_rnorm and balance edits
    ColumnVector item_by_bal;
    ColumnVector item_by_rnorm;
    if (Data.is_case(i_original,1)) {
      item_by_rnorm = Data.get_item_by_norm_indicator(s_q,item_by_bal);
    } else {
      item_by_rnorm = Data.copy_non_balance_edit(s_q);
    }
    
    //Generate from normal distribution
    ColumnVector mu_z_i = Mu.column(z_in(i_original)) ;
    ColumnVector tilde_y_i = Data.log_D_Observed.row(i_original).t();
    ColumnVector y_q(n_var);
    LowerTriangularMatrix LSigma_i;
    double log_cond_norm_q = calculate_log_cond_norm(Data, i_original, item_by_rnorm, tilde_y_i, y_q, true, LSigma_i, s_q); // MODIFIED 2015/02/16
    
    // Put values from balance edits
    ColumnVector x_q = exp_ColumnVector(y_q);
    if (Data.is_case(i_original,1)) {
      Data.set_balance_edit_values_for_x_q(s_q, x_q, item_by_bal); // CHANGED by Hang, 2014/12/29
    } else {
      Data.update_full_x_for_balance_edit(x_q);
    }

    // Acceptance/Rejection
    if (Data.PassEdits(x_q)) {  // Check constraints
      ColumnVector item_i_by_rnorm;
      if (Data.is_case(i_original,1)) {
        item_i_by_rnorm = Data.get_item_by_norm_indicator(s_i,item_by_bal);
      } else {
        item_i_by_rnorm = Data.copy_non_balance_edit(s_i);
      }
    
    double log_cond_norm_i = calculate_log_cond_norm(Data, i_original, item_i_by_rnorm, tilde_y_i, y_q, false, LSigma_i, s_i); // MODIFIED 2015/02/16

      y_q = log_ColumnVector(x_q) ;
      ColumnVector y_i = (Y_in.row(i_original)).t();
  
      // Calculate acceptance ratio
      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)));
      
    
      double log_f_y_tilde_q = 0;
      if (FM.useMap && Data.is_case(i_original,2)) {
        log_f_y_tilde_q = FM.Simulate_logUnif_case2(CFeasibilityMap::s_to_tau_fn(s_q),i_original,n_simul,randUnif);
      } else {
        log_f_y_tilde_q = Data.get_logUnif_y_tilde(s_q, CFeasibilityMap::s_to_tau_fn(s_q),i_original); 
      }
      double logNum = log_full_norm_q + log_f_y_tilde_q - log_cond_norm_q - log(g_s_q);
      double logDen = log_full_norm_i + logUnif_y_tilde(i_faulty) - log_cond_norm_i - log(g_s_i);
    
      accept_rate(1) = exp( logNum - logDen ) ;
      if ( randUnif.Next() < accept_rate(1) ){
        Y_in.row(i_original) = y_q.t() ;
        S_Mat.column(i_faulty) = s_q;
        logUnif_y_tilde(i_faulty) = log_f_y_tilde_q;
        is_accept(1)++;
      }
    }
    
  }
  is_accept(1) = is_accept(1) / Data.n_faulty;
}