#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "nrutil.h"


void PCA(float **star,int Nstar,int Mp,double **basef,double **coff);

void psfPCA(float *tstar, int tNstar, int tMp, double *tbasef, double *tcoff)
{
  int i, j, k;
  float **star;
  double **basef;
  double **coff;

  fprintf(stdout, "libPCA.so: load libPCA for psf-Matrix-[%d, %d]\n", tNstar, tMp);
  fflush(stdout);

  star = matrix(0, tNstar-1, 0, tMp-1);
  for(i=0; i<tNstar; i++)
    for(j=0; j<tMp; j++)
    {
      k = j + i*tMp;
      star[i][j] = tstar[k];
    }

  int Nstar, Mp;
  Nstar = tNstar;
  Mp = tMp;

  basef= dmatrix(0, tNstar-1, 0, tMp-1);
  coff = dmatrix(0, tNstar-1, 0, tNstar-1);

  PCA(star, Nstar, Mp, basef, coff);

  for(i=0; i<tNstar; i++)
    for(j=0; j<tMp; j++)
    {
      k = j + i*tMp;
      tbasef[k] = basef[i][j];
    }
  for(i=0; i<tNstar; i++)
    for(j=0; j<tNstar; j++)
    {
      k = j + i*tNstar;
      tcoff[k] = coff[i][j];
    }
    fprintf(stdout, "libPCA.so: psf-PCAs have been found!\n");
    fflush(stdout);
}


void PCA(float **star,int Nstar,int Mp,double **basef,double **coff)
{
  
  void svdcmp(double **a, int m, int n, double *w, double **v);
  void indexx(int n, double arr[], int *indx);
  double **a,**v,*w,*check;
  int i,j,k,*indx;
  if(Mp>=Nstar){
    a=dmatrix(1,Mp,1,Nstar);
    v=dmatrix(1,Nstar,1,Nstar);
    w=dvector(1,Nstar);
    indx=ivector(1,Nstar);
    for(i=1;i<=Nstar;i++){for(j=1;j<=Mp;j++)a[j][i]=star[i-1][j-1];}
    svdcmp(a,Mp,Nstar,w,v);   
    indexx(Nstar,w,indx);    
    for(j=1;j<=Nstar;j++){
      k=indx[Nstar-j+1];
      for(i=1;i<=Mp;i++)basef[j-1][i-1]=a[i][k];
    }
    for(i=1;i<=Nstar;i++){
      for(j=1;j<=Nstar;j++){
	k=indx[Nstar-j+1];
	coff[i-1][j-1]= w[k]*v[i][k];
      }
    }
    //for(j=1;j<=Nstar;j++){for(i=1;i<=Mp;i++)basef[j-1][i-1]=a[i][j];}
    //for(i=1;i<=Nstar;i++){for(k=1;k<=Nstar;k++)coff[i-1][k-1]= w[k]*v[i][k];}
    
    free_dmatrix(a,1,Mp,1,Nstar);
    free_dmatrix(v,1,Nstar,1,Nstar);
    free_dvector(w,1,Nstar);
    free_ivector(indx,1,Nstar);
  }

  else{
    a=dmatrix(1,Nstar,1,Mp);
    v=dmatrix(1,Mp,1,Mp);
    w=dvector(1,Mp);
    indx=ivector(1,Mp);
    for(i=1;i<=Nstar;i++){for(j=1;j<=Mp;j++)a[i][j]=star[i-1][j-1];}
    svdcmp(a,Nstar,Mp,w,v);   
    indexx(Mp,w,indx);  
    for(j=1;j<=Mp;j++){    
      k=indx[Mp-j+1];
      for(i=1;i<=Mp;i++)basef[j-1][i-1]=v[i][k];
    }
    for(i=1;i<=Nstar;i++){
      for(j=1;j<=Mp;j++){
	k=indx[Mp-j+1];
	coff[i-1][j-1]= w[k]*a[i][k];
      }
    }  
    free_dmatrix(a,1,Nstar,1,Mp);
    free_dmatrix(v,1,Mp,1,Mp);
    free_dvector(w,1,Mp);
    free_ivector(indx,1,Mp);
  }
}
void PCA_err(float **err,int Nstar,int Mp,int nbf,double **basef,double **cofferr)
{
  float gasdev(long *iseed);
  long iseed=-1;
  int nc=100,i,j,k,ns;
  double sum,*mean,**coeffs,*start,*errt,det;
  FILE *fp;
  start=dvector(0,Mp);
  mean=dvector(0,nbf);
  coeffs=dmatrix(0,nc,0,nbf);
  for(ns=0;ns<Nstar;ns++){
    for(k=0;k<nbf;k++)mean[k]=0;
    for(k=0;k<nc;k++){
      for(j=0;j<Mp;j++)start[j]=gasdev(&iseed)*err[ns][j];
      for(i=0;i<nbf;i++){
	coeffs[k][i]=0;
	for(j=0;j<Mp;j++)coeffs[k][i]+=start[j]*basef[i][j];
	mean[i]+=coeffs[k][i];
      }
    }   
    for(k=0;k<nbf;k++)mean[k]/=nc;
    for(i=0;i<nbf;i++){
      sum=0;
      for(k=0;k<nc;k++){
	det=coeffs[k][i]-mean[i];
	sum+=det*det;
      }
      cofferr[ns][i]=sqrt(sum/(nc-1.));
    }
  }
  
  free_dvector(start,0,Mp);
  free_dvector(mean,0,nbf);
  free_dmatrix(coeffs,0,nc,0,nbf);
}