//==============================================================================
//
// CorrMatr.c
//
// developer: Henry Guennadi Levkin
//
//==============================================================================

#include <stdio.h>

#include <math.h>

//==============================================================================
// print  vector
//==============================================================================
void VectorPrint(double *pVect, int nDim)
{
  int i;

  printf("--------------------------------------------------------------------\n");
  for(i=0; i<nDim; i++)
  {
    printf("%0.3lf ",pVect[i]);
  }
  printf("\n");  
}

//==============================================================================
//  Function for calculation of correlation matrix:
//
//  input:
//  nPoints    - number of data records for processing
//  nDim       - dimention of one data record
//  pfDataMatr - data matrix
//
//  output:
//  pfCorrMatr - matrix with coefficients of correlation
//
//  developer: Henry Guennadi Levkin
//
//==============================================================================
int CorrelationMatrix(int nPoints, int nDim, double* pfDataMatr, double* pfCorrMatr)
{
  int i, j, k;

  double fCorr;
  double fDenom;
  double fPoints = (double) nPoints;

  double* pfVectX;
  double* pfVectY;

  double fXSum  = 0., fYSum  = 0.;
  double fX2Sum = 0., fY2Sum = 0.;
  double fXYSum = 0.;

  if(nPoints < 2) return 1; // not enough data!
  if(nDim    < 2) return 2; // not enough parameters!
  
  pfVectX =(double*) malloc(nPoints * sizeof(double));
  pfVectY =(double*) malloc(nPoints * sizeof(double));

  for(i=0; i<nDim-1; i++) // first parameter
  {
    // read first vector
    for(k=0; k<nPoints; k++)
    {
      pfVectX[k] = pfDataMatr[k*nDim + i];
    }

    for(j=(i+1); j<nDim; j++) // second parameter
    {
      //read second vector
      for(k=0; k<nPoints; k++)
      {
        pfVectY[k] = pfDataMatr[k*nDim + j];
      }

      fXSum  = 0.0;   fYSum  = 0.0;
      fX2Sum = 0.0;   fY2Sum = 0.0;
      fXYSum = 0.0;

      // calculate correlation between parameters i and j
      for(k=0; k<nPoints; k++)
      {
        fXSum += pfVectX[k];
        fYSum += pfVectY[k];

        fX2Sum += pfVectX[k] * pfVectX[k];
        fY2Sum += pfVectY[k] * pfVectY[k];
    
        fXYSum += pfVectX[k] * pfVectY[k];
      }

      fDenom = (fPoints*fX2Sum - fXSum*fXSum)*(fPoints*fY2Sum - fYSum*fYSum);
      if(fDenom < 1.E-100) return 3;
  
      fCorr  = fPoints*fXYSum - fXSum*fYSum;
      fCorr  = fCorr / sqrt(fDenom);

      pfCorrMatr[j+i*nDim] = fCorr;
    }
  }
  
  for(i=0; i<nDim; i++)
  {
    pfCorrMatr[i+i*nDim] = 1.0;
  }

  for(i=0; i<nDim; i++)
  {
    for(j=(i+1); j<nDim; j++)
    {
      pfCorrMatr[i+j*nDim] = pfCorrMatr[j+i*nDim];
    }
  }

  return 0;
}

//==============================================================================
// read data matrix from text file
int MatrixReadFromFile(char *pFilename, int* pnRows, int* pnColumns, double** ppfMatr)
{
  FILE* pFile;
  int i;
  int nSize;
  
  pFile = fopen(pFilename, "r");
  if(pFile == NULL) 
  {
    return 1;
  }

  fscanf(pFile,"%d",pnRows);
  fscanf(pFile,"%d",pnColumns);

  nSize = (*pnRows) * (*pnColumns);
  *ppfMatr = (double*) malloc(nSize * sizeof(double));

  if(*ppfMatr == NULL)
  {
    return 2;
  }

  for(i=0; i<nSize; i++)
  {
    fscanf(pFile, "%lf", &((*ppfMatr)[i]));
  }

  return 0;
}

//==============================================================================
// print  matrix
void MatrixPrint(double *pMatr, int nRows, int nColumns)
{
  int i, j;

  printf("--------------------------------------------------------------------\n");
  for(i=0; i<nRows; i++)
  {
    for(j=0; j<nColumns; j++)
    {
      printf("%.4lf\t",pMatr[j+i*nColumns]);
    }
    printf("\n");
  }
}


//==============================================================================
//  TEST for function CorrelationMatrix :
//==============================================================================
int main(int nArgs, char** ppArgs)
{
  int nObjects;
  int nDim;
  double *pArr;
  double *pCorMatr;

  double fCorr;
  
  if(MatrixReadFromFile("datamatr.txt", &nObjects, &nDim, &pArr))
  {
    printf("File not found!\n");
    return 1;
  }
  printf("Matrix with data:\n");
  MatrixPrint(pArr, nObjects, nDim);

  pCorMatr = (double*) malloc(nDim * nDim * sizeof(double));

  printf("Correlation Matrix:\n");
  CorrelationMatrix(nObjects, nDim, pArr, pCorMatr);
  
  MatrixPrint(pCorMatr, nDim, nDim);
  
  free(pCorMatr);
  free(pArr);

  return 0;
}