bindata2ipca.cc

  1 const char *help = "\
  2 progname: bindata2ipca.cc\n\
  3 code2html: This program a bindata file and a PCA model, projects patterns into it and reconstructs the patterns in the input space.\n\
  4 version: Torch3 vision2.0, 2003-2005\n\
  5 (c) Sebastien Marcel (marcel@idiap.ch)\n";
  6 
  7 #include "DiskXFile.h"
  8 #include "ImageGray.h"
  9 #include "PCAMachine.h"
 10 #include "matrix.h"
 11 #include "MyMeanVarNorm.h"
 12 #include "CmdLine.h"
 13 
 14 using namespace Torch;
 15 
 16 int main(int argc, char *argv[])
 17 {
 18 	char *bindata_filename;
 19 	char *model_filename;
 20 	int n_input;
 21 	char *filename_out;
 22 	real variance;
 23 	bool verbose;
 24 	int verbose_level;
 25 	bool normalize;
 26 	CmdLine cmd;
 27 	cmd.setBOption("write log", false);
 28 	
 29 	cmd.info(help);
 30 	cmd.addText("\nArguments:");
 31 	cmd.addSCmdArg("bindata filename", &bindata_filename, "bindata filename");
 32 	cmd.addSCmdArg("model filename", &model_filename, "PCA model filename");
 33 	cmd.addICmdArg("n_input", &n_input, "number of inputs");
 34 	cmd.addText("\nOptions:");
 35 	cmd.addSCmdOption("-o", &filename_out, "ipca.bindata", "bindata output file");
 36 	cmd.addRCmdOption("-variance", &variance, 0.95, "variance (-1 => 100\%)");
 37 	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
 38 	cmd.addICmdOption("-verbose_level", &verbose_level, 1, "level of verbose");
 39 	cmd.addBCmdOption("-normalize", &normalize, false, "normalize");
 40 	cmd.read(argc, argv);
 41 
 42 	if(verbose == false) verbose_level = 0;
 43 
 44 
 45 	//
 46 	PCAMachine *pca_machine = NULL;
 47 	pca_machine = new PCAMachine(n_input);
 48 
 49 	//
 50 	MyMeanVarNorm *mv_norm = NULL;
 51 
 52 	if(normalize)
 53 	{
 54 	   	if(verbose_level >= 1)
 55 			print("Loading PCA model and its normalisation: %s ...\n", model_filename);
 56 
 57 		mv_norm = new MyMeanVarNorm(n_input, 1);
 58 
 59 	   	pca_machine->load(model_filename, mv_norm);
 60 	}
 61 	else
 62 	{
 63 	   	if(verbose_level >= 1)
 64 			print("Loading PCA model: %s ...\n", model_filename);
 65 
 66 		DiskXFile *file = NULL;
 67 		file = new DiskXFile(model_filename, "r");
 68 		pca_machine->loadXFile(file);
 69 		delete file;
 70 	}
 71 	
 72 	//
 73 	pca_machine->setIOption("verbose_level", 1);
 74 	pca_machine->setROption("variance", variance);
 75 	pca_machine->init();
 76 
 77 	int n_output = pca_machine->n_outputs;
 78 
 79 	//
 80 	Mat *E_inv = NULL;
 81 	E_inv = new Mat(n_input, n_input);
 82 
 83 	print("Copying eigenvectors ...\n");
 84 	   
 85 	Mat *eigenvectors = new Mat(n_input, n_input);
 86 
 87 	for(int i = 0 ; i < n_input ; i++)
 88        		for(int j = 0 ; j < n_input ; j++)
 89 			eigenvectors->ptr[j][i] = pca_machine->eigenvectors[j][i];
 90 		
 91 	print("Transpose the eigenvectors matrix ...\n");
 92 	mxTrMat(eigenvectors, E_inv);
 93 	delete eigenvectors;
 94 
 95 	//
 96 	int dimIn;
 97 	int n_patterns;
 98 		
 99 	DiskXFile *pf = NULL;
100 	pf = new DiskXFile(bindata_filename, "r");
101 
102 	pf->read(&n_patterns, sizeof(int), 1);
103 	pf->read(&dimIn, sizeof(int), 1);
104 
105 	if(verbose_level >= 1)
106 	{
107 		print("n_inputs : %d\n", dimIn);
108 		print("n_patterns : %d\n", n_patterns);  
109 	}
110 
111 	if(n_input > dimIn)
112 	{
113 		error("Number of inputs specified (%d) bigger than into the file (%d)\n", n_input, dimIn);
114 	   
115 		delete pf;
116 
117 		return 0;
118 	}
119 
120 
121 	//
122 	float *realinput = NULL;
123         Sequence *seq;
124 	float *ipca = NULL;
125 
126         realinput = new float [n_input];
127         seq = new Sequence(&realinput, 1, n_input);	
128         ipca = new float [n_input];
129 
130 	//
131 	DiskXFile *pfOutput = NULL;
132 	pfOutput = new DiskXFile(filename_out, "w");
133 	
134 	//
135 	if(verbose_level >= 1)
136 		print("Projection bindata file into PCA space (%d -> %d) ...\n", n_input, n_output);
137 
138 	int P = n_patterns;
139 
140 	pfOutput->write(&P, sizeof(int), 1);
141 	pfOutput->write(&dimIn, sizeof(int), 1);
142 	
143 	real MSE = 0.0;
144 
145 	for(int p = 0 ; p < P ; p++)
146 	{
147 	   	//
148 	  	for(int i = 0 ; i < n_input ; i++)
149 			pf->read(&realinput[i], sizeof(float), 1);   
150 
151 	   	if(verbose_level >= 2)
152 			print(" Seq =     [%2.3f %2.3f %2.3f ...]\n", realinput[0], realinput[1], realinput[2]);
153 
154 		if(normalize)
155 			mv_norm->preProcessInputs(seq);
156 
157 		//
158 		pca_machine->forward(seq);
159 			
160 	   	if(verbose_level >= 2)
161 			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]);
162 
163 		float *machine_output = pca_machine->outputs->frames[0];
164 		
165 		for(int i = 0 ; i < n_input ; i++)
166 		{
167 			ipca[i] = 0.0;
168 
169 			for(int j = 0 ; j < n_output ; j++)
170 				ipca[i] += E_inv->ptr[j][i] * machine_output[j];
171 		}
172 		
173 		for(int i = 0 ; i < n_input ; i++) 
174 		{
175 		   	ipca[i] += pca_machine->Xm[i];
176 
177 			if(ipca[i] < 0.0) ipca[i] = 0.0;
178 			if(ipca[i] > 1.0) ipca[i] = 1.0;
179 		}
180 
181 		//
182 		float data;
183 		float mse = 0;
184 		for(int i = 0 ; i < n_input ; i++)
185 		{
186 		   	data = ipca[i];
187 
188 			real z = realinput[i] - data;
189 			mse += z*z;
190 			
191 			pfOutput->write(&data, sizeof(float), 1);   
192 		}
193 
194 		mse /= (float) n_input;
195 
196 		print(" mse = %g\n", mse);
197 		
198 		MSE += mse;
199 	}
200 
201 	MSE /= (real) P;
202 
203 	print("Total MSE = %g\n", MSE);
204 
205 	//
206 	delete pf;
207 	delete pfOutput;
208 	
209 	//
210 	delete seq;
211 	delete [] realinput;
212 	delete [] ipca;
213 	if(normalize) delete mv_norm;
214 	delete pca_machine;
215 	delete E_inv;
216 	
217   	return 0;
218 }
219