bindata2ipcamse.cc

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