#include <stdio.h>
#include <math.h>
#include <time.h>

#define   EXCLUSION   15 /* data length not to be considered as neighbors */

main(argc,argv)
  int   argc;
  char  *argv[];
{
  int    i;
  int    step;
  int    length;    /* length of the data */ 
  int    dim;       /* embedding dimension */
  int    nnn;       /* number of the nearest neighbors */
  double p_nnn;     /* percentage of the number of nearest neighbors */
  double *xx;       /* variable for the time series data */ 
  double snr,SNR(); /* Signal-to-Noise ratio */ 
  double max_snr;   /* Maximam value of the SNR */ 
  FILE  *fp;

  /* reading parameters */ 
  fp=fopen("parameters.txt","r");
  if(fp==NULL){
    fprintf(stderr,"parameter file not found.\n"); 
    exit(0);
  }
  fscanf(fp,"%d %d",&length,&dim);
  fclose(fp);

  /* reading data from input file */ 
  xx = (double*)calloc(length,sizeof(double));
  fp=fopen("data.txt","r");
  if(fp==NULL){
    fprintf(stderr,"data file not found.\n"); 
    exit(0);
  }
  for(i=0;i<length;i++) fscanf(fp,"%lf",(xx+i));
  fclose(fp);

  step = (int)(0.025*((double)length));
  /* writing SNR values to output file */ 
  fp=fopen("dvs_plot.txt","w");
  for(nnn=dim+1;nnn<(length-dim-2*EXCLUSION);nnn+=step){
    /* computing the Signal-to-Noise ratio with a
       change in the number of nearest neighbors, nnn */ 
    snr = SNR(length,dim,nnn,xx);
    p_nnn = 100.0*((double)nnn)/((double)length);

    /* finding the maximum SNR */ 
    if(nnn==step || max_snr < snr) max_snr=snr; 

    /* printing SNR with percentage of the number of
       nearest neighbors */ 
    fprintf(fp,"%lf %lf\n",p_nnn,snr);
    fflush(fp);
  }
  fclose(fp);

  /* computing the Signal-to-Noise ratio by linear 
     prediction using all data points */ 
  nnn = length-dim-2*EXCLUSION-1;
  snr = SNR(length,dim,nnn,xx);

  /* nonlinear measure is finally computed by the difference 
     between maximal SNR and SNR using all data points */
  printf("Nonlinear Measure = %lf [dB]\n",max_snr-snr);

  fp=fopen("nlm_value.txt","w");
  fprintf(fp,"Nonlinear Measure = %lf [dB]\n",max_snr-snr);
  fclose(fp);
}    

double SNR(length,dim,nnn,xx)
  int    length,dim,nnn;
  double *xx; 
{
  int    i,j; 
  int    *nb;
  double *wd;
  double *rr,*ss;
  double Variance();
  double vs,vr;
  void   NearestNeighbors();
  void   LPC_Estimate();

  /* variable for linear predictive coefficients */
  wd = (double*)calloc(dim,sizeof(double)); 

  /* variable to locate original time series data */
  ss = (double*)calloc(length,sizeof(double));

  /* variable for residual data by local linear prediction */
  rr = (double*)calloc(length,sizeof(double));

  /* index for nearest neighbors */
  nb = (int*)calloc(length,sizeof(int));

  for(j=0;j<length;j++) *(ss+j)=*(xx+j);
  for(i=dim-1;(i+1)<length;i++){
    *(rr+i+1)=*(ss+i+1);

    /* finding the neighbors nearest to #i data */
    NearestNeighbors(length,dim,ss,i,nb);

    /* computing the linear coefficients by using 
       the nearest neighbors */
    LPC_Estimate(dim,nnn,ss,nb,wd);

    /* computing the residual signals by the local
       linear prediction */
    for(j=0;j<dim;j++) rr[i+1]-= (*(wd+j))*(*(ss+i-j));
  }

  /* computing variance of the original signal */
  vs=Variance(length,dim,ss);

  /* computing variance of the residual signal */
  vr=Variance(length,dim,rr);

  free(wd);
  free(rr);
  free(ss);
  free(nb);

  /* returning the Signal-to-Noise ratio */ 
  return(10.0*log10((vs/vr)));
}

double Variance(length,dim,x)
  int    length,dim;
  double *x;
{
  int    i;
  double norm,mean,variance;
  norm=mean=variance=0;
  for(i=dim;i<length;i++){ mean+=*(x+i); norm+=1.0; }
  mean=mean/norm;
  for(i=dim;i<length;i++) variance+=pow(*(x+i)-mean,2.0);
  variance=sqrt(variance/norm);
  return(variance);
}

void NearestNeighbors(length,dim,ss,c,nb)
  int    dim;
  double *ss;
  int    c;
  int    *nb;
{
  int    l,n,j; 
  double *distance;
  void Quick_Sort();
  distance = (double*)calloc(length,sizeof(double));
  for(n=0,l=dim-1;(l+1)<length;l++){
    if(EXCLUSION<abs(l-c)){
      for(*(distance+n)=0,j=0;j<dim;j++)
        *(distance+n)+=pow(*(ss+c-j)-*(ss+l-j),2.0);
      *(distance+n) = sqrt(*(distance+n));
      *(nb+n)=l; n++;
    }
  }
  Quick_Sort(nb,distance,0,n-1);
  free(distance);
}  

void Quick_Sort(order,x,left,right)
  int    *order; 
  double *x;
  int    left,right;
{
  int    i,j,k;
  int    v;
  double y,pivot;
  void Quick_Sort();
  if(left>=right) return;
  for(i=left-1,j=right,pivot=*(x+right);;){
    while(*(x+(++i)) < pivot);
    while(i < --j && pivot < *(x+j));
    if(i >= j) break;
    y=*(x+i); *(x+i)=*(x+j); *(x+j)=y;    
    k=*(order+i); *(order+i)=*(order+j); 
    *(order+j)=k;
  }
  y=*(x+i); *(x+i)=*(x+right); *(x+right)=y;
  k=*(order+i); *(order+i)=*(order+right); 
  *(order+right)=k; v=i;
  Quick_Sort(order,x,left,v-1);
  Quick_Sort(order,x,v+1,right);
}

void LPC_Estimate(dim,nnn,ss,nb,wd)
  int    dim,nnn;
  double *ss;
  int    *nb;
  double *wd;
{
  int    i,j,k,p;
  double Const;
  double Max;
  double *mx,*vc,*vt;

  vt=(double*)calloc(dim,sizeof(double));
  vc=(double*)calloc(2*dim,sizeof(double));
  mx=(double*)calloc(dim*2*dim,sizeof(double));

  for(i=0;i<dim;i++){ 
    *(vt+i)=0.0;
    for(k=0;k<nnn;k++) *(vt+i)+=(*(ss+*(nb+k)-i))*(*(ss+*(nb+k)+1));
  } 

  for(i=0;i<dim;i++)
    for(j=0;j<dim;j++){ 
      *(mx+2*dim*i+j) = 0.0;
      for(k=0;k<nnn;k++) 
        *(mx+2*dim*i+j) += (*(ss+*(nb+k)-i))*(*(ss+*(nb+k)-j));
    }

  for(i=0;i<dim;i++)
    for(j=0;j<dim;j++){
      if(i==j) *(mx+2*dim*i+j+dim)=1.0;
      else     *(mx+2*dim*i+j+dim)=0.0;
    }

  for(i=0;i<dim;i++){
    Max=fabs(*(mx+2*dim*i+i)); 
    p=i;
    for(j=i+1;j<dim;j++)
      if(Max<fabs(*(mx+2*dim*j+i))){ 
        p=j; Max=fabs(*(mx+2*dim*j+i));
      }
    if(i!=p){
      for(k=0;k<2*dim;k++) *(vc+k)         =*(mx+2*dim*i+k); 
      for(k=0;k<2*dim;k++) *(mx+2*dim*i+k) =*(mx+2*dim*p+k);
      for(k=0;k<2*dim;k++) *(mx+2*dim*p+k) =*(vc+k);
    }
    Const = *(mx+2*dim*i+i);
    for(j=0;j<2*dim;j++) *(mx+2*dim*i+j)/=Const;
    for(k=0;k<dim;k++){
      Const = *(mx+2*dim*k+i);
      for(j=0;j<2*dim;j++){
        if(k==i) break;
        *(mx+2*dim*k+j) -= Const**(mx+2*dim*i+j);
      }
    }
  }
  for(i=0;i<dim;i++)
    for(*(wd+i)=0,j=0;j<dim;j++)
      *(wd+i) += *(mx+2*dim*i+j+dim)*(*(vt+j)); 

  free(vt);
  free(vc);
  free(mx);
}