// TO compile use: 
// g++ matrixalgebra.cpp -lm -o matal -Wno-deprecated
// To run use: :   ./matalg 10 15 f.txt m.txt out.xls 
// where 10 is the number of variables and 15 is the number of cases
// f.txt is the fingerprint files and m.txt is the output we are trying to
// determine the coefficients for.

// Builted by Peter A Noble at Alabama State University on November 2, 2013

#include <fstream>
#include <string>
#include <iostream>
#include <math.h>
#include "matrix_h"

using namespace std;

#define columns 20
#define rows 100000

double A[columns][rows];
double nA[columns][rows];
double sA[columns][rows];
double B[rows];
double AT[columns][rows];
double IAT[columns][rows];
double AB[rows];
double C[rows];

/// ******************************************************************///
/// ******************************************************************///

double max(int column,int x);						// DETERMINES THE MAX OF ARRAY
double max2(int column,int x);						// DETERMINES THE MAX OF ARRAY
double min(int column,int x);					    // DETERMINES THE MIN OF ARRAY
double sum(int column,int x);					    // DETERMINES THE SUM OF ARRAY
double sum2(int x);					    			// DETERMINES THE SUM OF ARRAY
double correl(int r, int c);						// DETERMINES THE CORREL OF ARRAY


/// ******************************************************************///
/// ******************************************************************///

double correl(int r, int c)
{
//
// This subroutine calculates the Pearson correlation between two
// strings of numbers

double y1=0.0,y2=0.0,sumsqy1=0.0,sumsqy2=0.0,Sumxy=0.0,SSY1=0.0,
	SSY2=0.0,SP=0.0;int y=0,i=r;

for (y=0;y<(r);y++) 	
	{
	y1=B[y]+y1;
	sumsqy1=(B[y]*B[y])+sumsqy1;
	y2=A[c][y]+y2;
	sumsqy2=(A[c][y]*A[c][y])+sumsqy2;
	Sumxy=Sumxy+(B[y]*A[c][y]);
	}
SSY1=sumsqy1-((y1*y1)/double(i));
SSY2=sumsqy2-((y2*y2)/double(i));
SP=Sumxy-((y1*y2)/double(i));
//cout << " i= " << i << " corr=  " << SP/sqroot(SSY1*SSY2) <<"\n" << flush; 

return SP/sqrt(SSY1*SSY2);

}

/// ******************************************************************///
/// ******************************************************************///

double sum2(int x)
{
int j=0; double temp=0;
	for (j=0;j<x;j++)
	{
	//cout << j << "\t" << C[j] << "\n" << flush;
	temp=C[j]+temp;
	}
return temp;
}

/// ******************************************************************///
/// ******************************************************************///

double sum(int column,int x)
{
int j=0; double temp=0;
	for (j=0;j<column;j++)
	{
	temp=A[j][x]+temp;
	}
return temp;
}

/// ******************************************************************///
/// ******************************************************************///

double max2(int column,int x)
{
int j=0; double maxnumber=sA[0][x];
	for (j=0;j<column;j++)
	{
	if (sA[j][x]>maxnumber) { maxnumber=sA[j][x];}
	}
return maxnumber;
}

/// ******************************************************************///
/// ******************************************************************///

double max(int column,int x)
{
int j=0; double maxnumber=A[0][x];
	for (j=0;j<column;j++)
	{
	if (A[j][x]>maxnumber) { maxnumber=A[j][x];}
	}
return maxnumber;
}

/// ******************************************************************///
/// ******************************************************************///

double min(int column,int x)
{
int j=0; double minnumber=A[0][x];
	for (j=0;j<column;j++)
	{
	if (A[j][x]<minnumber) { minnumber=A[j][x];}
	}
return minnumber;
}

/// ******************************************************************///
/// ******************************************************************///
/*
void dumpMatrixValues(matrix <double> M)  {
  bool xyz;
  double rv;
  for (int i=0; i < M.getactualsize(); i++)
    {
    cout << "i=" << i << ": ";
    for (int j=0; j<M.getactualsize(); j++)
      {
        M.getvalue(i,j,rv,xyz);
        cout << rv << " ";
      }
    cout << endl;
    }
};
*/
void put_backMatrixValues(matrix <double> M)  {
  bool xyz;
  double rv;
  for (int i=0; i < M.getactualsize(); i++)
    {
    for (int j=0; j<M.getactualsize(); j++)
      {
        M.getvalue(i,j,rv,xyz);
		IAT[i][j]=rv; 
	//	cout << "put_back_into_matrix_complete\n";
      }
    }
};
/// ******************************************************************///
/// ******************************************************************///

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

int c=0,r=0;

c=atoi(argv[1]);
r=atoi(argv[2]);

ifstream inA(argv[3]);
ifstream inB(argv[4]);
ofstream out(argv[5]); //proportions

double cor=0.0;
double r2=0.0;
//int c=3,r=86652;   /// <==== need to set columns and rows here.  86652

for (int y=0;y<r;y++) 	
    {
	for (int x1=0;x1<c;x1++) 	
		{
		inA >> A[x1][y]; // Fingerprints
		}
	}
	
for (int x1=0;x1<r;x1++) 	
    {
    inB >> B[x1]; /// Mixtures
	}

double temp=0;

for (int h=0;h<c;h++) 	
	{
	for (int b=0;b<c;b++) 	
		{
		temp=0;
		for (int y=0;y<r;y++) 	
			{	
			temp=(A[h][y]*A[b][y])+temp; 
			}
		AT[b][h]=temp;	
		}
	}

for (int x=0;x<c;x++) 	
	{	
	for (int y=0;y<r;y++) 	
		{	
		AB[x]=(A[x][y]*B[y])+AB[x]; 
		}
	}

matrix <double> M1(c,c); 
int k = 0;

  for (int i=0; i < M1.getactualsize(); i++)  
    for (int j=0; j<M1.getactualsize(); j++)
      {
        M1.setvalue(i,j,AT[i][j]);
        k++;
      }
 
M1.invert();  // invert the matrix

for (int j=0; j < c; j++)  // define random values for initial matrix
	{
	out << "Coef_" << j+1 << "\t" ;
	}
	out << "\tR2\n" ;
	

put_backMatrixValues(M1);
 
for (int j=0; j < c; j++)  // define random values for initial matrix
	{
	 temp=0; 
	 
	for (int i=0; i < c; i++)  // define random values for initial matrix
		{
		temp=IAT[j][i]*AB[i]+temp;
		}
		C[j]=temp;
		out << temp << "\t";
	}	
//	out << "\t";

double temp2=0;
//	out <<  r << "\n";

for (int j=0; j < c; j++)  // define random values for initial matrix
	{
	temp2=((double)C[j]/sum2(c));
	//out <<  temp2 << "\t";
	}
 
// Now go back and see how well it works!
	for (int y=0;y<r;y++) 	
		{
		temp=0;
		for (int x=0;x<c;x++) 	
			{
			temp=A[x][y]*C[x]+temp; // Fingerprints
			}
			A[c][y]=temp;
		}

		cor=correl(r,c);
		r2=cor*cor;
		out <<"\t" << r2 <<"\n" << flush;
		cout << "\nR2= " << r2 <<"\n\n" << flush;

exit(1);
}