main_lasso.cpp

统计学习软件包

	  min_mean = meanE;
	  best_mu = mu_vec[ii];
	  best_std = spreadE;
	  best_rv = rv_vec[ii];
	  best_cv = ii;
      }
      cout <<"std. dev.: "<<sqrt(spreadE) << endl;
      cout <<"#RVs: " << rv_vec[ii] << endl;
      cout << "==========================" << endl;
  }

  cout << endl;
  cout << endl;
  spreadE = best_std;
  meanE = min_mean;
  mu = best_mu;

  cout << "=============================================" << endl;
  cout << "Optimal model: "<< endl;
  cout <<  "kappa: " << best_mu <<endl;
  cout <<"mean error: "<<meanE <<" (" << sqrt(spreadE)<<")"<< endl;
  cout << "#RVs: " << best_rv << endl;
  cout << "=============================================" << endl;


  if(retrain){
      cout << "*******************************************" << endl;
      cout <<"Re-training with optimal model:" <<endl;
      cout << "*******************************************" << endl;

       for(i=0; i<iter; i++){    
	   cout <<"Iteration: "<<i << endl;
	   cross_rand(fold_n);

	   for(j=0; j<XG->Ncols(); j++){
	       index[j] =0;
	       beta[j] = 0.0;
	   }

	   ri =  rand()%(XG->Ncols()-1);
	   index[ri] = 1;
	   beta[ri] = 1e-7;
	   
	   LM = new LogregLASSO(A, Y, const_val, mu, index, beta);
	   
	   LM->regress();
	   
	   erg->el(i,min_numb) = Predict(min_numb);    
	   rv_vec[min_numb] += rv_numb;
	   
       }

       meanE = 0.;
       spreadE = 0.;
       for(i=0; i<iter; i++)
	  meanE += erg->el(i,min_numb);
       meanE /= iter;
       for(i=0; i<iter; i++)
	   spreadE += pow(erg->el(i,min_numb)-meanE,2);
       spreadE /= iter;

       cout << endl << endl;
       cout <<"================================================== " << endl;
       cout <<"Error of optimal model: "<<meanE << " (" << sqrt(spreadE) <<")"<< endl;
       cout <<"================================================== " << endl;

      
      char line[10000];
      ii =0;
      

    
      
      ofstream ofile("relevant_genes.dat", ios::out);
      
      ofile << "# genes: "<< A->Ncols()-1<< endl;
      ofile << "kappa = " << mu << endl;
      ofile <<"mean error: "<<meanE << endl;
      ofile <<"std. dev: "<<sqrt(spreadE) << endl;
      ofile <<"#RVs: " << best_rv << endl;
      
      ifstream inf2;
      int thresh = int(iter/20.0);
      int* rel_rv_hist = new int[1000];
      char genes[1000][1000];

      ii =0;
      for(j=0; j< A->Ncols()-1; j++){
	 
	  if( rv_hist_mat[j][min_numb] > thresh){
	      ostrstream ossr(genes[ii],  1000);
	      ossr << rv_hist_mat[j][min_numb]/ (double)iter <<'\t'<< j 
		   <<'\t'<<  gene_names[j] <<'\0';
	      ii++;
	  }
      }
      inf2.close();
      int s_size =  sizeof genes[0];
      qsort(genes, ii, s_size, compare);
      int* ref_number = new int[ii];
      double score;
      for(i=0; i<ii; i++){
	  ofile << genes[i] <<endl;
	  istrstream issr(genes[i],  1000);
	  issr >> score >> ref_number[i] >> line;
      }
      ofile.close();
      
      char** gene_names_reduced = new char*[ii];
      for(i=0; i< ii; i++)
	  gene_names_reduced[i] = new char[1000];
      Relg->ReDimension(ii,ntotal);
      for(i=0; i<ii; i++){
	  for(k=0; k<ntotal; k++){
	      Relg->el(i,k) = XG->el(k,ref_number[i]);	    
	  }
	  strcpy(gene_names_reduced[i], gene_names[ref_number[i]]); 
	  rel_rv_hist[i] = rv_hist_mat[ref_number[i]][min_numb];
      }
      
      for(i=0; i<ii; i++){
	  strcpy(gene_names[i], gene_names_reduced[i]);
      }

      for(i=0; i< ii; i++)
	  delete []  gene_names_reduced[i];
      delete [] gene_names_reduced;

      *Relg_Buff = *Relg;
      
      int rel_num = ii;
      
      double max_r = -1e20;
      double min_r = 1e20;
      double mean_r =0.;
         
      for(i=0; i< rel_num; i++){
	  for(j=0; j<ntotal; j++){
	      if( Relg->el(i,j) > max_r)
		  max_r =  Relg->el(i,j);
	      if( Relg->el(i,j) < min_r)
		  min_r =  Relg->el(i,j);
	      
	  }
      }
      
      for(i=0; i< rel_num; i++)
	  for(j=0; j<ntotal; j++){
	      Relg->el(i,j) =  ((Relg->el(i,j) -min_r)/(max_r-min_r)) *255.0*1.0;
	  }
      
      mean_r =0.;
      for(i=0; i< rel_num; i++){
	  for(j=0; j<ntotal; j++){
	      mean_r +=  Relg->el(i,j);
	  }
      }
      mean_r /= (double)(rel_num*ntotal);
      
      for(i=0; i< rel_num; i++)
	  for(j=0; j<ntotal; j++){
	      Relg->el(i,j) =  (Relg->el(i,j)/ mean_r) *128.0;
	      if( Relg->el(i,j) <0)
		  Relg->el(i,j) =0;
	      if( Relg->el(i,j) >255)
		  Relg->el(i,j) =255;
	  }
      
      int* perm_local = new int[ntotal];
      
      int jj =0,a,b; 
      for(j=0; j<ntotal ;j++)
	  if(true_labs[j] < 0.5){
	      perm_local[jj] = j;
	      jj++;
	  }
      for(j=0; j<ntotal ;j++)
	  if(true_labs[j] > 0.5){
	      perm_local[jj] = j;
	      jj++;
	  }
      
      double sum_1 =0, sum_2 =0;
      int* rows_class = new int[rel_num];
      
      int anz_1 =0, anz_2 =0;
      for(i=0; i< rel_num; i++){
	  sum_1 =0;
	  sum_2 =0;
	  for(j=0; j<ntotal ;j++){
	      if(true_labs[j] < 0.5){ 
		  sum_1+=  Relg->el(i,j);
		  anz_1++;
	      }
	      else{
		  sum_2+=  Relg->el(i,j);
		  anz_2++;
	      }
	  }
	  if(sum_1/(double)anz_1 < sum_2/(double)anz_2 )
	      rows_class[i] = -1;
	  else
	      rows_class[i] = 1;
      }
      int* hist_norm = new int[rel_num];
      double max_h =0.;
      int scale_factor = 100;
      for(i=0; i< rel_num; i++)
	  if( rel_rv_hist[i] > max_h )
	      max_h =  rel_rv_hist[i];
      for(i=0; i< rel_num; i++){
	  hist_norm[i] = int(((double)rel_rv_hist[i]/max_h) * scale_factor);
      }
      ofstream of5("rg.pgm", ios::out);
      of5 <<  "P2" << endl;
      of5 <<"# rg.pgm " << endl;
      of5 << 5*ntotal + scale_factor+5<<' '<< 5*rel_num+ 10<< endl;
      of5 << 255 << endl; 
      // label indicator bar 
      for(b=0; b<5; b++){
	  for(j=0; j<ntotal ;j++)
	      for(a=0;a<5;a++)
		  if(true_labs[perm_local[j]] < 0.5)
		      of5 << 200 <<' ';
		  else
		      of5 << 0 <<' '; 
	  for(a=0;a<5;a++)
	      of5 << 255 <<' ';
	  for(j=0; j<scale_factor ;j++)
	      of5 << 255 <<' ';
	  
	  of5 << endl;
      }
      for(b=0; b<5; b++){
	  for(j=0; j<ntotal ;j++)
	      for(a=0;a<5;a++)
		  of5 << 255 <<' ';
	  for(a=0;a<5;a++)
	      of5 << 255 <<' ';
	  for(j=0; j<scale_factor ;j++)
	      of5 << 255 <<' ';
	  of5 << endl;
      }
      // main image
      for(i=0; i< rel_num; i++){
	  if( rows_class[i] == -1){
	      for(b=0; b<5; b++){
		  for(j=0; j<ntotal ;j++)
		      for(a=0;a<5;a++)
			  of5 << int(Relg->el(i,perm_local[j])) << ' ';
		  // historam
		  for(a=0;a<5;a++)
		      of5 << 255 <<' ';
		  for(j=0; j<scale_factor ;j++){
		      if(j+1 <  hist_norm[i])
			  of5 << 0<<' ';
		      else
			  of5 << 255 <<' ';
		  }
		  of5 << endl;
	      }
	  }
      }
      for(i=0; i< rel_num; i++){
	  if( rows_class[i] == 1){
	      for(b=0; b<5; b++){
		  for(j=0; j<ntotal ;j++)
		      for(a=0;a<5;a++)
			  of5 << int(Relg->el(i,perm_local[j])) << ' ';
		  for(a=0;a<5;a++)
		      of5 << 255 <<' ';
		  for(j=0; j<scale_factor ;j++){
		      if(j+1 <  hist_norm[i])
			  of5 << 0<<' ';
		      else
			  of5 << 255 <<' ';
		  }
		  of5 << endl;
	      }
	  }
      }
      of5.close();
      double sum_w =0.0, sum_c =0.0;
      for(i=8; i<20; i++){
	  sum_w += hist_wrong_mat->el(i,min_numb);
	  sum_c += hist_correct_mat->el(i,min_numb); 
      }
  
      
      
      double bd[20];
      double cd[20];
      double f[20];
  
      for(i=0; i< 20; i++){
	  f[i] = (i+1)/20.0-0.025;
	  bd[i] = (double)hist_wrong_mat->el(i,min_numb)/sum_w;
      }
      for(i=0; i< 20; i++){
	  cd[i] = (double)hist_correct_mat->el(i,min_numb)/sum_c;
      }
      double err[11];
      double unc[11];
      
      ofstream o2("doubt.dat" , ios::out);
      for(i=9; i< 20; i++){
	  err[i-9] =0.;
	  for(j=9; j<=i; j++){
	      err[i-9]+= bd[j];
	  }
	  err[i-9] = meanE*(1.0-err[i-9]);
	  
      }
      for(i=9; i< 20; i++){
	  unc[i-9] =0.;
	  for(j=9; j<=i; j++){
	      unc[i-9]+= cd[j];
	  }
	  unc[i-9] = (1.0-meanE)*(unc[i-9]);
	  o2  <<  f[i]+0.025 <<' '<<err[i-9] 
	      <<' '<<unc[i-9]+ meanE*(1.0-err[i-9]/meanE) << endl;
      }

      delete [] rel_rv_hist;
      delete [] ref_number;
      delete [] hist_norm;
      delete [] rows_class;
      delete [] perm_local;

  }


    delete erg;
 
    delete [] index;
    delete [] beta;
    delete [] mu_vec;
    delete [] rv_vec;

  
}
main(int argc, char** argv){

    srand(time(NULL));
    
    if(argc !=2)
	cout << "Usage: GenLASSO file_name" << endl;
    else{
	ifstream inf1(argv[1], ios::in);
	double cross_start, cross_add;
	int  ntotal, dim ,   cross_number, iter, reduce_flag;
	char data_file[500];
	char annotation_file[500];

	inf1 >> ntotal >> dim;
	inf1 >> cross_start >>  cross_add >> cross_number;
	inf1 >> iter;
	inf1 >> data_file;
	inf1 >> annotation_file;
	inf1 >> reduce_flag;
	inf1.close();
	
	main_LASSO lasso(data_file, annotation_file, ntotal, dim, cross_start, cross_add, cross_number, iter);
	lasso.ReadData();
	lasso.Classify();
	if(reduce_flag)
	    lasso.Classify_Reduced();
    }
}