bindata2ipca.cc

 const char *help = "\
 progname: bindata2ipca.cc\n\
 code2html: This program a bindata file and a PCA model, projects patterns into it and reconstructs the patterns in the input space.\n\
 version: Torch3 vision2.0, 2003-2005\n\
 (c) Sebastien Marcel (marcel@idiap.ch)\n";
 
 #include "DiskXFile.h"
 #include "ImageGray.h"
 #include "PCAMachine.h"
 #include "matrix.h"
 #include "MyMeanVarNorm.h"
 #include "CmdLine.h"
 
 using namespace Torch;
 
 int main(int argc, char *argv[])
 {
 	char *bindata_filename;
 	char *model_filename;
 	int n_input;
 	char *filename_out;
 	real variance;
 	bool verbose;
 	int verbose_level;
 	bool normalize;
 	CmdLine cmd;
 	cmd.setBOption("write log", false);
 	
 	cmd.info(help);
 	cmd.addText("\nArguments:");
 	cmd.addSCmdArg("bindata filename", &bindata_filename, "bindata filename");
 	cmd.addSCmdArg("model filename", &model_filename, "PCA model filename");
 	cmd.addICmdArg("n_input", &n_input, "number of inputs");
 	cmd.addText("\nOptions:");
 	cmd.addSCmdOption("-o", &filename_out, "ipca.bindata", "bindata output file");
 	cmd.addRCmdOption("-variance", &variance, 0.95, "variance (-1 => 100\%)");
 	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
 	cmd.addICmdOption("-verbose_level", &verbose_level, 1, "level of verbose");
 	cmd.addBCmdOption("-normalize", &normalize, false, "normalize");
 	cmd.read(argc, argv);
 
 	if(verbose == false) verbose_level = 0;
 
 
 	//
 	PCAMachine *pca_machine = NULL;
 	pca_machine = new PCAMachine(n_input);
 
 	//
 	MyMeanVarNorm *mv_norm = NULL;
 
 	if(normalize)
 	{
 	   	if(verbose_level >= 1)
 			print("Loading PCA model and its normalisation: %s ...\n", model_filename);
 
 		mv_norm = new MyMeanVarNorm(n_input, 1);
 
 	   	pca_machine->load(model_filename, mv_norm);
 	}
 	else
 	{
 	   	if(verbose_level >= 1)
 			print("Loading PCA model: %s ...\n", model_filename);
 
 		DiskXFile *file = NULL;
 		file = new DiskXFile(model_filename, "r");
 		pca_machine->loadXFile(file);
 		delete file;
 	}
 	
 	//
 	pca_machine->setIOption("verbose_level", 1);
 	pca_machine->setROption("variance", variance);
 	pca_machine->init();
 
 	int n_output = pca_machine->n_outputs;
 
 	//
 	Mat *E_inv = NULL;
 	E_inv = new Mat(n_input, n_input);
 
 	print("Copying eigenvectors ...\n");
 	   
 	Mat *eigenvectors = new Mat(n_input, n_input);
 
 	for(int i = 0 ; i < n_input ; i++)
        		for(int j = 0 ; j < n_input ; j++)
 			eigenvectors->ptr[j][i] = pca_machine->eigenvectors[j][i];
 		
 	print("Transpose the eigenvectors matrix ...\n");
 	mxTrMat(eigenvectors, E_inv);
 	delete eigenvectors;
 
 	//
 	int dimIn;
 	int n_patterns;
 		
 	DiskXFile *pf = NULL;
 	pf = new DiskXFile(bindata_filename, "r");
 
 	pf->read(&n_patterns, sizeof(int), 1);
 	pf->read(&dimIn, sizeof(int), 1);
 
 	if(verbose_level >= 1)
 	{
 		print("n_inputs : %d\n", dimIn);
 		print("n_patterns : %d\n", n_patterns);  
 	}
 
 	if(n_input > dimIn)
 	{
 		error("Number of inputs specified (%d) bigger than into the file (%d)\n", n_input, dimIn);
 	   
 		delete pf;
 
 		return 0;
 	}
 
 
 	//
 	float *realinput = NULL;
         Sequence *seq;
 	float *ipca = NULL;
 
         realinput = new float [n_input];
         seq = new Sequence(&realinput, 1, n_input);	
         ipca = new float [n_input];
 
 	//
 	DiskXFile *pfOutput = NULL;
 	pfOutput = new DiskXFile(filename_out, "w");
 	
 	//
 	if(verbose_level >= 1)
 		print("Projection bindata file into PCA space (%d -> %d) ...\n", n_input, n_output);
 
 	int P = n_patterns;
 
 	pfOutput->write(&P, sizeof(int), 1);
 	pfOutput->write(&dimIn, sizeof(int), 1);
 	
 	real MSE = 0.0;
 
 	for(int p = 0 ; p < P ; p++)
 	{
 	   	//
 	  	for(int i = 0 ; i < n_input ; i++)
 			pf->read(&realinput[i], sizeof(float), 1);   
 
 	   	if(verbose_level >= 2)
 			print(" Seq =     [%2.3f %2.3f %2.3f ...]\n", realinput[0], realinput[1], realinput[2]);
 
 		if(normalize)
 			mv_norm->preProcessInputs(seq);
 
 		//
 		pca_machine->forward(seq);
 			
 	   	if(verbose_level >= 2)
 			print(" Output =   [%2.3f %2.3f %2.3f ...]\n", pca_machine->outputs->frames[0][0], pca_machine->outputs->frames[0][1], pca_machine->outputs->frames[0][2]);
 
 		float *machine_output = pca_machine->outputs->frames[0];
 		
 		for(int i = 0 ; i < n_input ; i++)
 		{
 			ipca[i] = 0.0;
 
 			for(int j = 0 ; j < n_output ; j++)
 				ipca[i] += E_inv->ptr[j][i] * machine_output[j];
 		}
 		
 		for(int i = 0 ; i < n_input ; i++) 
 		{
 		   	ipca[i] += pca_machine->Xm[i];
 
 			if(ipca[i] < 0.0) ipca[i] = 0.0;
 			if(ipca[i] > 1.0) ipca[i] = 1.0;
 		}
 
 		//
 		float data;
 		float mse = 0;
 		for(int i = 0 ; i < n_input ; i++)
 		{
 		   	data = ipca[i];
 
 			real z = realinput[i] - data;
 			mse += z*z;
 			
 			pfOutput->write(&data, sizeof(float), 1);   
 		}
 
 		mse /= (float) n_input;
 
 		print(" mse = %g\n", mse);
 		
 		MSE += mse;
 	}
 
 	MSE /= (real) P;
 
 	print("Total MSE = %g\n", MSE);
 
 	//
 	delete pf;
 	delete pfOutput;
 	
 	//
 	delete seq;
 	delete [] realinput;
 	delete [] ipca;
 	if(normalize) delete mv_norm;
 	delete pca_machine;
 	delete E_inv;
 	
   	return 0;
 }