//Programm zur Berechnung von Trend und Periodizitaet einer unregelmaessig beprobten Zeitreihe 
//Die Periodizität wird über die Korrelation mit einer (Co)Sinusschwingung ermittelt, indem die 
//"normierte" Zeitreihe (Mittelwert abgezogen, durch Standardabweichung dividiert) mit berechneten 
//Schwingungen verglichen wird. Die Korrelationskoeffizienten werden nach DAVIS ermittelt.
//Die Summe der Abweichungsquadrate (SAQ) ist besser geeignet! Die Signifikanzgrenzen wurden nach
//HARTUNG (1982) ermittelt. Es ergeben sich leichte Abweichungen zu den Tabellenwerken, da 
//aufgrund der numerischen Probleme auf eine Korrektur von der Normalverteilung auf die 
//STUDENT-Verteilung verzichtet wurde. Der Fehler wirkt sich bei den Korrelationskoeffizienten
//bei n>100 an der 2. Nachkommastelle aus.

//Program for the calculation of trend and periodicity of an unevenly sampled time series
//The periodicity is calculated by a correlation of the sampled time series with a sine-oscillation
//The sample time series is normed (average subtracted, divided by standard deviation). The 
//korrelation coefficient is calculated according to DAVIS. The sum of squared deviation may be more 
//efficient, therefore it is calculated too. In calculation od significance there are small 
//deviations to tables because of a missing correction to the STUDENT-distribution. This mistake
//is significant only to the 2nd position after the delimiter at n>100.


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

void welcome()
{
cout<<"Welcome to a special kind of time series analysis!\n";
cout<<"(C) Wolfgang Gossel, Muehlweg 49, D-06114 Halle\n";
cout<<"-------------------------------------------------------------\n";
cout<<"This time series analysis calculates\n";
cout<<"- statistical measures (average, variance, standarddeviation)\n";
cout<<"- a trend\n";
cout<<"- a periodogram\n";
cout<<"- the residuals of measurements, average, trend and periodicity analysis\n";
cout<<"of unevenly sampled datasets.\n";
cout<<"-------------------------------------------------------------\n";
}

int datenein1(char datei_ein[])
{
int a=0;
double b=0.0;

cout<<"Please write the name of the file\n";
cout<<"with the dataset to be analyzed!\n";
cout<<"The file must be a textfile with\n";
cout<<"- time (integer), ordered, and\n";
cout<<"- value (float)\n";
cout<<"and a blank between both columns.\n";
cout<<"File name:";
cin>>datei_ein;

int i=0;

ifstream fin(datei_ein);

while (!fin.eof())
  {
  fin>>a>>b;
  i++;
  }
fin.close();
i=i-1;
return i; 
}

void null_setzen(int wert[], int i)
{
int a=0;
for (a=0;a<=i-1;a++)
  {
  wert[a]=0;
  }
}

void null_setzen(double wert[], int i)
{
int a=0;
for (a=0;a<=i-1;a++)
  {
  wert[a]=0.0;
  }
}


void datenein2(char datei_ein[], int i, int twert[], double messwert[], double koeff[])
{
int a=0;
int *atwert = new int[i];
null_setzen(atwert,i);

ifstream fin(datei_ein);

while (!fin.eof())
  {
  fin>>atwert[a]>>messwert[a];
  a++;
  }

koeff[1]=atwert[0];
for (a=0;a<=i-1;a++)
  {
  twert[a]=atwert[a]-atwert[0]+1;
  }
fin.close();    
}


double mittelwert(double wert[], int i)
{
int n=0;
double sum=0.0,mittel=0.0;

for (n=0;n<=i-1;n++)
 {
 sum=sum+wert[n];
 }
return mittel=sum/(i+1);
}


double varianz(double wert[], double koeff[], int i)
{
int n;
double saq=0.0,var=0.0;

for (n=0;n<=i-1;n++)
  {
  saq=saq+((wert[n]-koeff[2])*(wert[n]-koeff[2]));
  }
return var=saq/i;
}

void aus_stat(double koeff[], int i)
{
cout<<"Average y:			"<<koeff[2]<<"\n";
cout<<"Variance y:			"<<koeff[3]<<"\n";
cout<<"Standard deviation y:		"<<koeff[4]<<"\n";
cout<<"Sum of squared deviations:	"<<koeff[3]*i<<"\n";
}

void test(double koeff[], int i)
{
int j=i-2;
double t9=0.0, t5=0.0, t0=0.0, 
       u9=0.0, u5=0.0, u0=0.0; 

 t9=sqrt(-2*log(1-0.99));
 t5=sqrt(-2*log(1-0.95));
 t0=sqrt(-2*log(1-0.90));

 u9=t9-((2.515517+0.802853*t9+0.010328*t9*t9)/(1+1.432788*t9+0.189269*t9*t9+0.001308*t9*t9*t9));
 u5=t5-((2.515517+0.802853*t5+0.010328*t5*t5)/(1+1.432788*t5+0.189269*t5*t5+0.001308*t5*t5*t5));
 u0=t0-((2.515517+0.802853*t0+0.010328*t0*t0)/(1+1.432788*t0+0.189269*t0*t0+0.001308*t0*t0*t0));

koeff[9]=u9/(sqrt((u9*u9)+j));
koeff[10]=u5/(sqrt((u5*u5)+j));
koeff[11]=u0/(sqrt((u0*u0)+j));
}

double kkk_ber(double messwert[], double y[], int i)
{
int a=0;
double sumx=0.0, sumy=0.0, sumx2=0.0, sumy2=0.0, sumxy=0.0, qxy=0.0, qx=0.0, qy=0.0, kkk_schw=0.0;

for (a=0;a<=i-1;a++)
{
sumx=sumx+messwert[a];
sumy=sumy+y[a];
sumx2=sumx2+(messwert[a]*messwert[a]);
sumy2=sumy2+(y[a]*y[a]);
sumxy=sumxy+(messwert[a]*y[a]);
}
qxy=sumxy-((sumx*sumy)/i);
qx=sumx2-((sumx*sumx)/i);
qy=sumy2-((sumy*sumy)/i);
return kkk_schw=qxy/(sqrt(qx*qy));
}


double saq_ber(double messwert[], double y[], int i)
{
int a=0;
double b=0.0;

for (a=0;a<=i-1;a++)
  {
  b=b+((messwert[a]-y[a])*(messwert[a]-y[a]));
  }
return b;  
}

//Allgemeine Form der Ausgleichsgeraden, keine regelmäßige Zeitreihe nötig
void trend(int twert[],double messwert[],int i,double koeff[])
{
int a=0;
double sxy=0.0, sx=0.0,sy=0.0,s2x=0.0,s2y=0.0, qxy=0.0, qx=0.0, qy=0.0;
double *y = new double[i];

for (a=0;a<=i-1;a++)
 {
 sxy=sxy+(twert[a]*messwert[a]);
 sx=sx+twert[a];
 sy=sy+messwert[a];
 s2x=s2x+(twert[a]*twert[a]);
 s2y=s2y+(messwert[a]*messwert[a]);
 }
koeff[6]=(sxy-(sx*sy/i))/(s2x-((sx*sx)/i));      //Steigung
koeff[5]=(sy/i)-(koeff[6]*(sx/i));               //Ordinate
qxy=sxy-(sx*sy/i);
qx=s2x-(sx*sx/i);
qy=s2y-(sy*sy/i);
koeff[7]=fabs(qxy/(sqrt(qx*qy)));

for (a=0;a<=i;a++)
  {
  y[a]=koeff[5]+koeff[6]*twert[a];
  }
koeff[15]=kkk_ber(messwert,y,i);  
koeff[8]=saq_ber(messwert,y,i);
}

void aus_trend(double koeff[])
{
cout<<"The term for the trend is given by:\n";
cout<<"y = "<<koeff[6]<<" * x + "<<koeff[5]<<"\n";
cout<<"The sum of squared deviations is "<<koeff[8]<<"\n";
cout<<"The correlation coefficient is "<<koeff[7]<<" or "<<koeff[15]<<"\n";

if (koeff[7]<=koeff[11])
{
cout<<"The correlation coefficient is lower than the 90%-level for significance of about "<<koeff[11]<<".\n";
}

else if (koeff[7]<=koeff[10])
{
cout<<"The correlation coefficient is between the 90%-level(r="<<koeff[11]<<") and 95%-level(r="<<koeff[10]<<") of significance.\n";
}

else if (koeff[7]<=koeff[9])
{
cout<<"The correlation coefficient is between the 95%-level(r="<<koeff[10]<<") and 99%-level(r="<<koeff[9]<<") of significance.\n";
}

else if (koeff[7]>koeff[9])
{
cout<<"The correlation coefficient is higher than the 99%-level of significance of about "<<koeff[9]<<".\n";
}
}


void dx(int twert[],int i,double koeff[])
{
int a=0, diffx=0;
double sdiffx=0.0;
koeff[16]=1000000; koeff[17]=0; koeff[18]=0;
for (a=0;a<=i-2;a++)
{
diffx=twert[a+1]-twert[a];
if (koeff[16]>diffx) koeff[16]=diffx; //dxmin wird Minimum
if (koeff[17]<diffx) koeff[17]=diffx; //dxmax wird Maximum
sdiffx=sdiffx+diffx;
}
koeff[18]=sdiffx/(i-1);
}


void dx_fest(int twert[], int i, double koeff[])
{

dx(twert,i,koeff);
cout<<"For the calculation of periodicity \n- the start and  \n- the increment\n";
cout<<"of the periodogram have to be set.\nFrom the time series the following values can be derived:\n";
cout<<"The minimum of the time differences is: "<<koeff[16]<<"\n";
cout<<"The maximum of the time differences is: "<<koeff[17]<<"\n";
cout<<"The average time difference between measurements is: "<<koeff[18]<<"\n";
cout<<"The proposed values are:\n- for the start "<<2*koeff[18]<<" and\n- for the increment "<<koeff[18]<<"\n";
cout<<"Please set the start value and increment:\n\n- Start:  ";
cin>>koeff[13];
cout<<"\n- Increment:  ";
cin>>koeff[12];
float p_max=twert[i-1]/2;
koeff[14]=(p_max-koeff[13])/koeff[12];
}



void per_korr(int twert[], double messwert[], int i, double koeff[],double phi[],double ampl[], double kkk_per[], double saq_per[])
{
int a=0,p=0,o=0,c=0,d=0,e=100;
int min_phi;
double omega=0.0, t=0.0, phi2=0.0, saq2=0.0, min_saq=0.0, min_ampl=0.0, 
       k_ampl=0.0,s_ampl=0.0,k_help=0.0,k_phi=0.0,s_phi=0.0,saq3=0.0;

double *y = new double[i];
double *saq = new double[e];
null_setzen(y,i);
null_setzen(saq,e);



for (a=0;a<=koeff[14]-1;a++)
{
  cout<<a<<" from "<<koeff[14]<<"\n";
  null_setzen(y,i);
  d=koeff[13]+(koeff[12]*a);
  omega=6.2831583/d;
  k_help=0.0;
  for (p=0;p<=e;p++)
  { 
  phi2=(6.2831583/e)*p;
  
    for (c=0;c<=i-1;c++)
    {
      t=koeff[5]+koeff[6]*twert[c];
      y[c]=t+koeff[4]*cos(omega*twert[c]+phi2);
      min_saq=min_saq+sqrt((y[c]-messwert[c])*(y[c]-messwert[c]));
    }
  k_phi=kkk_ber(y,messwert,i);
  s_phi=min_saq;

    if (k_phi>k_help)
    {
      phi[a]=phi2;
      k_help=k_phi;
    }

  }
  k_help=0.0;
  null_setzen(y,i);
  for (o=0;o<=koeff[3]+1;o++)
  {
    for (d=0;d<=i-1;d++)
    {  
      t=koeff[5]+koeff[6]*twert[d];
      y[d]=t+o*cos(omega*twert[d]+phi[a]);
      saq2=saq2+sqrt((y[d]-messwert[d])*(y[d]-messwert[d]));
    }
  k_ampl=kkk_ber(y,messwert,i);
  s_ampl=saq2;  

    if (k_ampl>k_help)
    {
      ampl[a]=o;
      k_help=k_ampl;
    }
  null_setzen(y,i);
  k_ampl=0.0;
  }
  for (d=0;d<=i-1;d++)
  {
  t=koeff[5]+koeff[6]*twert[d];
  y[d]=t+ampl[a]*cos(omega*twert[d]+phi[a]);
  saq3=saq3+((y[d]-messwert[d])*(y[d]-messwert[d]));
  }
kkk_per[a]=kkk_ber(y,messwert,i);
saq_per[a]=saq3;
saq3=0.0;
saq2=0.0;
phi2=0.0;
min_saq=0.0;  
}
} 


void aus_per(int twert[], double messwert[], int i, double koeff[], double phi[], double ampl[],double kkk_per[], double saq_per[])
{
int a=0, n=0, wahl=0;
double p=0.0,kkk_e_t=0.0,kkk_e_p=0.0;
double *t = new double[i];
double *fp = new double[i];

cout<<"In the following periodogram the intensities are given by correlation coefficients (cc) and\n";
cout<<"sums of quadratic deviations (qd).\n";
cout<<"The phase shift (Phi) and amplitudes (A) are calculated by optimization of these values.\n";

cout<<"Periode       Phi          A           cc          qd\n";

for (a=0;a<=koeff[14]-1;a++)
{
  cout.width(5);cout<<koeff[13]+koeff[12]*a;
  cout.width(15);cout<<phi[a];
  cout.width(15);cout<<ampl[a];
  cout.width(15);cout<<kkk_per[a];
  cout.width(15);cout<<saq_per[a]<<"\n";
}
cout<<"\nPlease choose one period for the calculation of function values!\n";
cout<<"Period:   ";
cin>>n;
wahl=(n-koeff[13])/koeff[12];

cout<<"time    measurement     FunctionValue Trend      Residuals Trend      FunctionValue period       Residuals period";

for (a=0;a<=i-1;a++)
{
  p=koeff[13]+koeff[12]*wahl;
  t[a]=koeff[5]+koeff[6]*twert[a];
  fp[a]=t[a]+ampl[wahl]*cos((6.2831583/p)*twert[a]+phi[wahl]);
  cout.width(5);cout<<twert[a];
  cout.width(15);cout<<messwert[a];
  cout.width(15);cout<<t[a];
  cout.width(15);cout<<messwert[a]-t[a];
  cout.width(15);cout<<fp[a];
  cout.width(15);cout<<messwert[a]-fp[a]<<"\n";
}
kkk_e_t=kkk_ber(messwert,t,i);
kkk_e_p=kkk_ber(messwert,fp,i);

cout<<"The correlation coefficient of the trend is:     "<<kkk_e_t<<"\n";
cout<<"The correlation coefficient of the oscillation is: "<<kkk_e_p<<"\n";

cout<<"The level of significance for 99% confidence interval is: "<<koeff[9]<<"\n";
cout<<"The level of significance for 95% confidence interval is: "<<koeff[10]<<"\n";
cout<<"The level of significance for 90% confidence interval is: "<<koeff[11]<<"\n";


char datei_aus[50];
cout<<"Please set the file name for output:";
cin>>datei_aus;

ofstream fout(datei_aus);

fout<<"Time    Measurement     FunctionValue Trend      Residuals Trend      FunctionValue Periods       Residuals Period";

for (a=0;a<=i-1;a++)
{
  p=koeff[13]+(koeff[12]*wahl);
  t[a]=koeff[5]+koeff[6]*twert[a];
  fp[a]=t[a]+koeff[4]*cos((6.2831583/p)*twert[a]);
  fout.width(5);fout<<twert[a];
  fout.width(15);fout<<messwert[a];
  fout.width(15);fout<<t[a];
  fout.width(15);fout<<messwert[a]-t[a];
  fout.width(15);fout<<fp[a];
  fout.width(15);fout<<messwert[a]-fp[a]<<"\n";
}
fout<<"\nThe correlation coefficient of the trend is:     "<<kkk_e_t<<"\n";
fout<<"The correlation coefficient of the oscillation is: "<<kkk_e_p<<"\n\n";
fout<<"The signifikance level for the 99% confidence interval is: "<<koeff[9]<<"\n";
fout<<"The signifikance level for the 95% confidence interval is: "<<koeff[10]<<"\n";
fout<<"The signifikance level for the 90% confidence interval is: "<<koeff[11]<<"\n";

fout<<"In the following periodogram the intensities are given by correlation coefficients (cc) and\n";
fout<<"sums of quadratic deviations (qd).\n";
fout<<"The phase shift (Phi) and amplitudes (A) are calculated by optimization of these values.\n\n";

fout<<"Period       Phi          A           cc          qd\n";

for (a=0;a<=koeff[14]-1;a++)
{
  fout.width(5);fout<<koeff[13]+koeff[12]*a;
  fout.width(15);fout<<phi[a];
  fout.width(15);fout<<ampl[a];
  fout.width(15);fout<<kkk_per[a];
  fout.width(15);fout<<saq_per[a]<<"\n";
}
fout.close();
}


int main(void)
{
  char datei_ein[50];

  welcome();
  
  int i=0, koeff_zahl=20;

  i=datenein1(datei_ein);

  int *twert = new int[i];
  double *messwert = new double[i];
  double *koeff = new double[koeff_zahl];

  null_setzen(twert,i);
  null_setzen(messwert,i);
  null_setzen(koeff,koeff_zahl);

  datenein2(datei_ein,i,twert,messwert,koeff);//koeff[1]=Anfangswert T
  
  koeff[2]=mittelwert(messwert,i);       //koeff[2]=Mittelwert
  koeff[3]=varianz(messwert,koeff,i);    //koeff[3]=Varianz
  koeff[4]=sqrt(koeff[3]);                   //koeff[4]=Standardabweichung

  aus_stat(koeff,i);
  
  trend(twert,messwert,i,koeff);         //koeff[5]=Ordinatenabschnitt
                                         //koeff[6]=Steigung
                                         //koeff[7]=KKK des Trends
                                         //koeff[8]=SAQ des Trends
                                          
  test(koeff,i);                         //koeff[9]=99%-Quantile
                                         //koeff[10]=95%-Quantile
                                         //koeff[11]=90%-Quantile
                                         
  aus_trend(koeff);

  dx_fest(twert,i,koeff);                //koeff[12]=Inkrement Periode
                                         //koeff[13]=Anfangswert Periode
                                         //koeff[14]=Anzahl Perioden

  int anz_p=koeff[14];
  
  double *phi = new double[anz_p];
  double *ampl = new double[anz_p];
  double *kkk_per = new double[anz_p];
  double *saq_per = new double[anz_p];

  null_setzen(kkk_per,anz_p);
  null_setzen(phi,anz_p);
  null_setzen(ampl,anz_p);
  null_setzen(saq_per,anz_p);
                                                                                
  per_korr(twert,messwert,i,koeff,phi,ampl,kkk_per,saq_per);

  aus_per(twert, messwert,i,koeff,phi,ampl,kkk_per,saq_per);

  
  return 0;
}