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

void CParam::init_logUnif_y_tilde(CData &Data,CFeasibilityMap &FM, Uniform &randUnif, int n_simul) {
  logUnif_y_tilde = ColumnVector(Data.n_faulty);
  for (int i = 1; i <= Data.n_faulty; i++) {
    ColumnVector s_i = S_Mat.column(i);
    int i_original = Data.Faulty2Original[i-1];
    // double log_f_y_tilde_q = 0;  Commented by HANG on 5/16/2015
    if (FM.useMap && Data.is_case(i_original,2)) {
      logUnif_y_tilde(i) = FM.Simulate_logUnif_case2(CFeasibilityMap::s_to_tau_fn(s_i),i_original,n_simul,randUnif);
    } else {
      logUnif_y_tilde(i) = Data.get_logUnif_y_tilde(s_i, CFeasibilityMap::s_to_tau_fn(s_i),i_original); 
    } 
    //logUnif_y_tilde(i) = Data.get_logUnif_y_tilde(s_i, CFeasibilityMap::s_to_tau_fn(s_i),i_original);
  }
}

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 CFeasibilityMap::initilize_D_and_S(CData &Data) {
	
  ColumnVector order_to_test = get_order_to_test(Data.n_tau, Data.n_var);
  ColumnVector list_feasible_type2 = get_feasible_tau(Data);
	
	for (int i_faulty = 1; i_faulty <= Data.n_faulty; i_faulty++) {  
		
    bool is_pass = false;
    int i_original = Data.Faulty2Original[i_faulty-1];
    ColumnVector x_tilde_i = (Data.D_Observed.row(i_original)).t(); 
    for (int i_order = 1; i_order <= order_to_test.nrows() && !is_pass; i_order++) {
      int i_tau = order_to_test(i_order);
      ColumnVector s_i = tau_to_s_fn(i_tau,Data.n_var);
    
      bool skip_for_type2 = false; 
      if (Data.is_case(i_original,2) && list_feasible_type2(i_tau) == 0) { skip_for_type2 = true;}
      if (!skip_for_type2){
        ColumnVector x_mean;
        int is_feasible = feasible_test_fn(Data, x_tilde_i, s_i, i_original, true, Data.epsilon, x_mean);
        
        if (is_feasible > 0) {
          //copy solution to a temp vector
          ColumnVector temp = x_tilde_i;
          for (int index = 1, count =0; index <= Data.n_var; index++){
            if (s_i(index) == 1) {
              count++;
              temp(index) = x_mean(count);
            }
          }
          Data.Debug = Debug;
          if (Data.PassEdits(temp)) { //then check if it satifies edits          
            is_pass = true;
            Data.initial_S_Mat.row(i_faulty) = s_i.t();   
            Data.D_initial.row(i_original) = temp.t();
          }
        }//is_feasible
        Debug = false;
      }
    } 
  
	}  
	
}

void CParam::init_sigmamu(CData &Data) {
  mu_bar = ColumnVector(n_var_independent);
  ColumnVector log_obs_mean(n_var), log_obs_sd(n_var) ;
  Matrix logD_editpassing(n_sample-Data.n_faulty,n_var) ;
  for (int i_sample=1, count = 1; i_sample<=n_sample; i_sample++){
    if (Data.is_case(i_sample,0)) {
      logD_editpassing.row(count++) = Data.log_D_Observed.row(i_sample);
    }
  }

  for (int i_var=1; i_var<=n_var; i_var++){
    ColumnVector temp_col = logD_editpassing.column(i_var) ;
    log_obs_mean(i_var) = 1.0/(n_sample-Data.n_faulty)*(temp_col.sum()) ;
    ColumnVector temp_mean(n_sample-Data.n_faulty) ; temp_mean = log_obs_mean(i_var) ;
    Matrix SumSq = (temp_col-temp_mean).t()*(temp_col-temp_mean) ;
    log_obs_sd(i_var) = sqrt( 1.0/(n_sample-Data.n_faulty-1)*SumSq(1,1) ) ;
  }

  Data.UpdateCompactVector(mu_bar,log_obs_mean);
  DiagonalMatrix LSigma_temp(n_var_independent);
  for (int i=1; i<= n_var_independent; i++){
    LSigma_temp(i) = log_obs_sd(i);
  }
  DiagonalMatrix LSigma_temp_i = LSigma_temp.i();

  LSIGMA = vector<LowerTriangularMatrix>(K);
  LSIGMA_i = vector<LowerTriangularMatrix>(K);
  SIGMA = vector<SymmetricMatrix>(K);
  for (int k=0 ; k<K; k++) {
    LSIGMA[k] = LSigma_temp;
    LSIGMA_i[k] = LSigma_temp_i;
    SIGMA[k] << LSigma_temp * LSigma_temp; //lazy
  }
  logdet_and_more = ColumnVector(K); logdet_and_more = -0.5*n_var* LOG_2_PI + logdet(LSIGMA_i[0]);

  Mu = Matrix((n_var_independent),K);      // Note that mu_k is Mu.column(k)
  for (int k=1; k<=K; k++){
    Mu.column(k) = mu_bar ;  // starting value
  }
}

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;
}