bindata2ipcamse.cc

torch tracking code, it is a good code

#include "DiskXFile.h"
#include "ImageGray.h"

#include "PCAMachine.h"
#include "matrix.h"
#include "MyMeanVarNorm.h"

#include "CmdLine.h"

using namespace Torch;

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\
(c) Sebastien Marcel 2003-2004\n";

int main(int argc, char *argv[])
{
	char *bindata_filename;
	char *model_filename;
	char *ipca_filename;
	int n_input;
	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.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.addSCmdOption("-ipca", &ipca_filename, "", "iPCA model filename");
	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);
	DiskXFile *pfMxInverse = NULL;

	if(strcmp(ipca_filename, "") == 0)
	{
		print("Copying eigenvectors ...\n");
	   
		Mat *eigenvectors = NULL;
        	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("Inversion of eigenvectors matrix ...\n");
	
		mxInverse(eigenvectors, E_inv);

		print("Saving matrix inverse ...\n");
		pfMxInverse = new DiskXFile("ipca.matrix", "w");
		pfMxInverse->write(&n_input, sizeof(int), 1);
       		for(int j = 0 ; j < n_input ; j++)
			pfMxInverse->write(E_inv->ptr[j], sizeof(double), n_input);
	
		delete eigenvectors;
	}
	else
	{
		print("Loading matrix inverse %s ...\n", ipca_filename);
		pfMxInverse = new DiskXFile(ipca_filename, "r");
		int n_input_;
		pfMxInverse->read(&n_input_, sizeof(int), 1);
		if(n_input_ != n_input)
			error("Number of inputs %d != %d incorrect\n", n_input_, n_input);
       		for(int j = 0 ; j < n_input ; j++)
			pfMxInverse->read(E_inv->ptr[j], sizeof(double), n_input);
	}

	delete pfMxInverse;

	//
	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];

	//
	if(verbose_level >= 1)
		print("Projection bindata file into PCA space (%d -> %d) ...\n", n_input, n_output);

	int P = n_patterns;

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

		mse /= (float) n_input;

		//print(" mse = %g\n", mse);
		
		MSE += mse;
	}

	MSE /= (real) P;

	print("MSE = %g\n", MSE);

	//
	delete pf;
	
	//
	delete seq;
	delete [] realinput;
	delete [] ipca;
	if(normalize) delete mv_norm;
	delete pca_machine;
	delete E_inv;
	
  	return 0;
}