//
// Purpose:
//    This library implements Yao's algorithm
//
// Authors:
//    Frederico R. B. Cruz
//    Rosangela H. Loschi
//    Departamento de Estatistica
//    Universidade Federal de Minas Gerais
//    E-mail: {fcruz,loschi}@est.ufmg.br
//
// Version:
//    1.00
//
// Date:
//    May/2003
//

#ifndef CHPOINT_H
#define CHPOINT_H
//
// inclusion of standard libraries
//
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
//
#ifndef PI
#define PI 3.14159265358979323846
#endif
//
// random number generators
//
#include "randx.c"
//
// seeds for the random number generators
//
unsigned long runifSeed1=362436069, runifSeed2=521288629;
int rnormSeed=13579;
int rbetaSeed=35791;
//
// gamma function
//
#include "gamma.cpp"
//
// time measurement tool
//
#include "clock.c"

#define float double

float NormEuclid(float *vec1, float *vec2, int size) {
  static int i;
  static float norm;
  norm=0.0;
  for (i=1;i<=size;i++) {
      norm+=pow((vec1[i]-vec2[i]),2);
//      fprintf(stdout,"NormEuclid: norm=%f\n",norm);
  }
  return sqrt(norm);
}

// class definition
class ChangePoints {
  FILE *fileData;
  FILE *fileBurnIn;
  FILE *fileMuSig;
  FILE *fileNumBlo;
  FILE *filePartic;
  // maximum number of partition to be printed
  int numPartic;
  // time series and its size
  float *X;
  int tamX;
  // probability p
  float pProb;
  // posterior parameters of mu
  float **mPost, **vPost;
  // posterior parameters of sigma^2
  float **dPost, **aPost;
  // predictive distribution
  float **fXij;
  // posterior cohesions
  float **cPost;
  // vector of u samples
  int  **vecU;
  // size of each sample u
  int tamU;
  // number of samples u
  int numAmost;
  // burn-in period
  int burnIn;
  // lag to be taken into account
  int lag;
  // "net" number of samples u which is int((numAmost-burnIn)/lag)
  int numAmostU;
  // number of occurences of the i-th partition
  int *partic;
  // number of blocks of the i-th sample
  int *numBlocks;
  // number of occurences of this block
  int *blocks;
  // vector lambda[0j]
  float *lambda0;
  // vector lambda[iN]
  float *lambdaN;
  // posterior relevance of the block [ij]
  float **relevPost;
  // posterior estimates of mu and sigma^2
  float *muPost, *sigma2Post;
  // private methods
  int ComputeParamPost(void);
  int ComputeCohesPost(void);
  int ComputeVecU(void);
  int ComputeNuBloc(void);
  int PrintNumBloc(void);
  int ComputeRelevYaoPost(void);
  int ComputeMuSigYaoPost(void);
  int PrintMuSigPost(void);
  int ComputeRhoPost(void);
  int ComputeRhoYaoPost(void);
  int PrintRhoPost(void);
  int ComputePbPost(void);
  int ComputePbYaoPost(void);
  int PrintPbPost(void);
public:
  ChangePoints(void);    // default constructor
  ~ChangePoints(void);   // destructor
  int ReadInputScript(char *scriptFileName);
  int Yao(void);
  int PrintResults();
};

// class implementation

ChangePoints::ChangePoints(void) {
    fprintf(stdout, "Object ChangePoints created\n");
    numPartic = 10;
}

ChangePoints::~ChangePoints(void) {
    fprintf(stdout, "Object ChangePoints destroied\n");
}

int ChangePoints::ReadInputScript(char *scriptFileName) {
    const int LLENGTH = 256;
    char Line[LLENGTH];
    int i;
    FILE *scriptFile;
    char fileName[LLENGTH];
    float dummy;

    fprintf(stdout, "ChangePoints::ReadData\n");
//    open script file
    fprintf(stdout, "Script file name\n\t%s\n", scriptFileName);
    scriptFile = fopen(scriptFileName, "r");
//    read script
    fgets(Line, LLENGTH, scriptFile);
    fgets(Line, LLENGTH, scriptFile);
    fscanf(scriptFile, "%d\n", &numAmost);
    fprintf(stdout, "numAmost\n\t%d\n", numAmost);
    fgets(Line, LLENGTH, scriptFile);
    fscanf(scriptFile, "%d\n", &burnIn);
    fprintf(stdout, "burnIn\n\t%d\n", burnIn);
    fgets(Line, LLENGTH, scriptFile);
    fscanf(scriptFile, "%d\n", &lag);
    fprintf(stdout, "Lag\n\t%d\n", lag);
    fgets(Line, LLENGTH, scriptFile);
    fscanf(scriptFile, "%lf\n", &pProb);
    fprintf(stdout, "probability p\n\t%f\n", pProb);
//
    fgets(Line, LLENGTH, scriptFile);
    fgets(Line, LLENGTH, scriptFile);
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ((fileName[i]=Line[i]) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for data\n\t%s\n",fileName);
    fileData = fopen(fileName, "r");
//
    fgets(Line, LLENGTH, scriptFile);
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ((fileName[i]=Line[i]) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for number of blocks\n\t%s\n",fileName);
    fileBurnIn = fopen(fileName, "w");
//
    fgets(Line, LLENGTH, scriptFile);
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ((fileName[i]=Line[i]) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for posterior estimates of mu & sigma^2\n\t%s\n",fileName);
    fileMuSig = fopen(fileName, "w");
//
    fgets(Line, LLENGTH, scriptFile);
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ((fileName[i]=Line[i]) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for posterior distribution of number of blocks\n\t%s\n",fileName);
    fileNumBlo = fopen(fileName, "w");
//
    fgets(Line, LLENGTH, scriptFile);
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ((fileName[i]=Line[i]) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for the most probable partitions\n\t%s\n",fileName);
    filePartic = fopen(fileName, "w");
    fclose(scriptFile);
//  count time series size
    tamX = 0;
    while (fscanf(fileData,"%lf\n", &dummy) == 1) {
          tamX++;
//          fprintf(stdout, "dummy = %lf\t tamX = %d\n", dummy, tamX);
    }
    // read time series
    rewind(fileData);
    fprintf(stdout, "Time series size\n\t%d\n", tamX);
    X = new float[tamX+1];
    tamU = tamX-1;
    for (i=1; i <= tamX; i++) {
        fscanf(fileData, "%lf\n", &X[i]);
//        fprintf(stdout, "X%d = %lf\n", i, X[i]);
    }
    return 0;
}

int ChangePoints::ComputeParamPost(void) {
    fprintf(stdout,"ChangePoints::ComputeParamPost\n");
    int i, j, tamBloc, k;
    float denom, XijBarr, Qij, Sum1, Sum2, Mij;
    mPost = new float*[tamX];
    vPost = new float*[tamX];
    aPost = new float*[tamX];
    dPost = new float*[tamX];
    fXij  = new float*[tamX];
    for (i=0; i<tamX; i++) {
      // create matrixes
      mPost[i] = new float[tamX+1];
      vPost[i] = new float[tamX+1];
      aPost[i] = new float[tamX+1];
      dPost[i] = new float[tamX+1];
      fXij[i] = new float[tamX+1];
      for (j=i+1; j<=tamX; j++) {
          // initialize constants
          tamBloc = j-i;
          denom = tamBloc*vPrior(i,j)+1;
          // compute average within block [ij]
          XijBarr = 0.0;
          for (k=i+1; k<=j; k++) {
              XijBarr += X[k];
          }
          XijBarr /= tamBloc;
          // compute Q for block [ij]
          Sum1 = 0.0;
          for (k=i+1; k<=j; k++) {
              Sum1 += pow( (X[k]-XijBarr),2 );
          }
          Qij = Sum1 + tamBloc * pow( (XijBarr-mPrior(i,j)), 2)/denom;
          // compute posterior parameters
          mPost[i][j] = (vPrior(i,j)*tamBloc*XijBarr+mPrior(i,j))/denom;
          vPost[i][j] =  vPrior(i,j) / denom;
          dPost[i][j] =  dPrior(i,j) + tamBloc;
          aPost[i][j] =  aPrior(i,j) + Qij;
          // compute predictive distribution
          Sum1 = 0.0;
          for (k=i+1; k<=j; k++) {
              Sum1 += pow( (X[k]-mPrior(i,j)),2 );
          }
          Sum2 = 0.0;
          for (k=i+1; k<=j; k++) {
              Sum2 += X[k]-mPrior(i,j);
          }
          Mij = Sum1 - pow(Sum2,2)*vPrior(i,j)/denom;
          fXij[i][j] = log(Gamma((dPrior(i,j)+tamBloc)/2)) +
                       dPrior(i,j)/2 * log(aPrior(i,j)) +
                       (-0.5) * log(1+tamBloc*vPrior(i,j)) +
                       (-(dPrior(i,j)+tamBloc)/2) * log(aPrior(i,j)+Mij) -
                       (tamBloc/2)*log(PI) - log(Gamma(dPrior(i,j)/2));
          fXij[i][j] = exp(fXij[i][j]);
//          fXij[i][j] = Gamma((dPrior(i,j)+tamBloc)/2) *
//                       pow(aPrior(i,j),dPrior(i,j)/2) *
//                       pow((1+tamBloc*vPrior(i,j)),-0.5) *
//                       pow((aPrior(i,j)+Mij),-(dPrior(i,j)+tamBloc)/2) /
//                       (pow(PI,tamBloc/2) * Gamma(dPrior(i,j)/2));
      }
    }
    // print results
//    fprintf(stdout, "mPost:\n");
//    for (i=0; i<tamX; i++) {
//        for (j=1; j<=i; j++)
//            fprintf(stdout, "         ");
//        for (j=i+1; j<=tamX; j++)
//            fprintf(stdout, "%f ", mPost[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "vPost:\n");
//    for (i=0; i<tamX; i++) {
//        for (j=1; j<=i; j++)
//            fprintf(stdout, "         ");
//        for (j=i+1; j<=tamX; j++)
//            fprintf(stdout, "%f ", vPost[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "aPost:\n");
//    for (i=0; i<tamX; i++) {
//        for (j=1; j<=i; j++)
//            fprintf(stdout, "         ");
//        for (j=i+1; j<=tamX; j++)
//            fprintf(stdout, "%f ", aPost[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "dPost:\n");
//    for (i=0; i<tamX; i++) {
//        for (j=1; j<=i; j++)
//            fprintf(stdout, "         ");
//        for (j=i+1; j<=tamX; j++)
//            fprintf(stdout, "%f ", dPost[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "fXij:\n");
//    for (i=0; i<tamX; i++) {
//        for (j=1; j<=i; j++)
//            fprintf(stdout, "         ");
//        for (j=i+1; j<=tamX; j++)
//            fprintf(stdout, "%f ", fXij[i][j]);
//        fprintf(stdout, "\n");
//    }
    return 0;
};

int ChangePoints::ComputeCohesPost(void) {
    fprintf(stdout,"ChangePoints::ComputeCohesPost\n");
    int i, j;
    cPost = new float*[tamX];
    for (i=0; i<tamX; i++) {
      // create matrixes
      cPost[i] = new float[tamX+1];
      for (j=i+1; j<tamX; j++) {
          // compute cohesion of block [ij]
          cPost[i][j] = log(pProb) + (j-i-1)*log(1-pProb) + log(fXij[i][j]);
          cPost[i][j] = exp(cPost[i][j]);
      }
      cPost[i][tamX] = (tamX-i-1)*log(1-pProb) + log(fXij[i][tamX]);
      cPost[i][tamX] = exp(cPost[i][tamX]);
    }
    // print results
//    fprintf(stdout, "cPost:\n");
//    for (i=0; i<tamX; i++) {
//        for (j=1; j<=i; j++)
//            fprintf(stdout, "         ");
//        for (j=i+1; j<=tamX; j++)
//            fprintf(stdout, "%f ", cPost[i][j]);
//        fprintf(stdout, "\n");
//    }
    return 0;
};

int ChangePoints::ComputeVecU(void) {
    fprintf(stdout, "ChangePoints::ComputeVecU\n");
    int k, i, l, x, y;
    float Ri, u;

    vecU = new int*[numAmost+1];
    for (k=0; k<=numAmost; k++) {
        vecU[k] = new int[tamU+1];
    }
    // create seed
    for (i=1; i<=tamU; i++) {
        vecU[0][i] = 0;
    }
//    fprintf(stdout, "vecU[%d]:\n", 0);
//    for (i=1; i<=tamU; i++) {
//        fprintf(stdout, "%d ", vecU[0][i]);
//    }
//    fprintf(stdout, "\n");
    // create kth sample
    for (k=1; k<=numAmost; k++) {
        for (i=1; i<=tamU; i++) {
            vecU[k][i] = vecU[k-1][i];
        }
        for (i=1; i<=tamU; i++) {
            // find x
            x = 0;
            for (l=1; l<i; l++) {
                if ( vecU[k][l] == 0 ) x = l;
            }
            // find y
            y = tamX;
            for (l=tamU; l>i; l--) {
                if ( vecU[k-1][l] == 0 ) y = l;
            }
            // compute Ri
            Ri = log(cPost[x][y]) - log(cPost[x][i]) - log(cPost[i][y]);
            Ri = exp(Ri);
            // generate Uki
//            u = UNI;
            u = rUnif2(&runifSeed1,&runifSeed2);
            if ( Ri >= (1-u)/u ) {
               vecU[k][i] = 1;
            } else {
               vecU[k][i] = 0;
            }
        }
//        fprintf(stdout, "vecU[%d]:\n", k);
//        for (i=1; i<=tamU; i++) {
//            fprintf(stdout, "%d ", vecU[k][i]);
//        }
//        fprintf(stdout, "\n");
    }
    return 0;
}

int ChangePoints::ComputeNuBloc(void) {
    fprintf(stdout, "ChangePoints::ComputeNuBloc\n");
    int k,i;
    numBlocks = new int[numAmost+1];
    for (k=1; k<=numAmost; k++) {
        numBlocks[k] = 1;
        for (i=1; i<=tamU; i++) {
            if (vecU[k][i] == 0) {
               numBlocks[k]++;
            }
        }
    }
   return 0;
}

int ChangePoints::PrintNumBloc() {
    fprintf(stdout, "ChangePoints::PrintNumBloc\n");
    for (int k=1; k<=numAmost; k++) {
        fprintf(fileBurnIn,"%d\n", numBlocks[k]);
    }
    return 0;
}

int ChangePoints::ComputeRelevYaoPost(void) {
    fprintf(stdout, "ChangePoints::ComputeRelevYaoPost\n");
    int i, j, k;
    lambda0 = new float[tamX+1];
    lambdaN = new float[tamX+1];
    relevPost = new float*[tamX];

    // compute lambda[0j]
    lambda0[0] = 1.0;
    lambda0[1] = cPost[0][1];
//    fprintf(stdout, "lambda0[0] is %f\n", lambda0[0]);
//    fprintf(stdout, "lambda0[1] is %f\n", lambda0[1]);
    for (j=2; j<=tamX; j++) {
        lambda0[j] = cPost[0][j];
//        fprintf(stdout, "lambda0[%d] is ", j);
//        fprintf(stdout, "%f ", cPost[0][j]);
        for (k=1; k<j; k++) {
            lambda0[j] += lambda0[k] * cPost[k][j];
//            fprintf(stdout, "+ (%f*%f) ", lambda0[k], cPost[k][j]);
        }
//        fprintf(stdout, " = %f\n", lambda0[j]);
    }
    // compute lambda[iN]
    lambdaN[tamX] = 1.0;
    lambdaN[tamX-1] = cPost[tamX-1][tamX];
//    fprintf(stdout, "lambdaN[%d] is %f\n", tamX, lambdaN[tamX]);
//    fprintf(stdout, "lambdaN[%d] is %f\n", tamX-1, lambdaN[tamX-1]);
    for (i=tamX-2; i>0; i--) {
        lambdaN[i] = cPost[i][tamX];
//        fprintf(stdout, "lambdaN[%d] is ", i);
//        fprintf(stdout, "%f ", cPost[i][tamX]);
        for (k=i+1; k<tamX; k++) {
            lambdaN[i] += cPost[i][k] * lambdaN[k];
//            fprintf(stdout, "+ (%f*%f) ", cPost[i][k], lambdaN[k]);
        }
//        fprintf(stdout, " = %f\n", lambdaN[i]);
    }
    // compute posterior relevance
    for (i=0; i<tamX; i++) {
        relevPost[i] = new float[tamX+1];
        for (j=i+1; j<=tamX; j++) {
//            relevPost[i][j] = log(lambda0[i]) + log(cPost[i][j]) +
//                              log(lambdaN[j]) - log(lambda0[tamX]);
//            relevPost[i][j] = exp(relevPost[i][j]);
            relevPost[i][j] = lambda0[i]*cPost[i][j]*lambdaN[j]/lambda0[tamX];
        }
    }
    // print results
//    for (i=0; i<tamX; i++) {
//        for (j=1; j<=i; j++)
//            fprintf(stdout, "  ");
//        for (j=i+1; j<=tamX; j++)
//            fprintf(stdout, "%f ", relevPost[i][j]);
//        fprintf(stdout, "\n");
//    }
    return 0;
}

int ChangePoints::ComputeMuSigYaoPost(void) {
    fprintf(stdout, "ChangePoints::ComputeMuSigYaoPost\n");
    int k, i, j;
//    float aux;
    muPost = new float[tamX+1];
    sigma2Post = new float[tamX+1];
    for (k=1; k<=tamX; k++) {
        muPost[k] = 0.0;
        sigma2Post[k] = 0.0;
        for (i=0; i<=k-1; i++) {
            for (j=k; j<=tamX; j++) {
//                aux = log(relevPost[i][j]) + log(mPost[i][j]);
                muPost[k] += relevPost[i][j] * mPost[i][j];
//                aux = log(relevPost[i][j]) + log(aPost[i][j]) -
//                      log(dPost[i][j] - 2);
                sigma2Post[k] += relevPost[i][j]*aPost[i][j]/
                                 (dPost[i][j] - 2);

            }
        }
    }
    return 0;
}

int ChangePoints::PrintMuSigPost(void) {
    fprintf(stdout, "ChangePoints::PrintMuSigPost\n");
    fprintf(fileMuSig, "muPost\t sigma2Post\n");
    for (int i=1; i<=tamX; i++) {
        fprintf(fileMuSig, "%f\t %f\n", muPost[i], sigma2Post[i]);
    }
    return 0;
}

int ChangePoints::ComputeRhoPost(void) {
    fprintf(stdout, "ChangePoints::ComputeRhoPost\n");
    int i, j, k;
    partic = new int[numAmost+1];
    for (i=0; i<=numAmost; i++) {
        partic[i] = 1;
    }
    for (i=burnIn+1; i<=numAmost; i+=lag) {
        if ( partic[i] != -1 ) {
           for (j=i+lag; j<=numAmost; j+=lag) {
               // comparar i-esima amostra com j-esima amostra
               k = 1;
               while ( (vecU[i][k]==vecU[j][k])&&(k<tamU)) {
                    k++;
               }
               if ( vecU[i][k]==vecU[j][k] ) {
                  partic[i]++;
                  partic[j] = -1;
               }
           }
        }
    }
    return 0;
}

int ChangePoints::ComputeRhoYaoPost(void) {
   fprintf(stdout, "ChangePoints::ComputeRhoYaoPost\n");
   // generate samples u
   ComputeVecU();
   // compute number of blocks for each sample
   ComputeNuBloc();
   // compute probability of each partition
   ComputeRhoPost();
   return 0;
}

int ChangePoints::PrintRhoPost(void) {
    fprintf(stdout, "ChangePoints::PrintRhoPost\n");
    int i, j, maxJ, max;
    fprintf(filePartic, "partition \t probability\t #\n");
    for (i=1; i<=numPartic; i++) {
        max = -1;
        for (j=burnIn+1; j<=numAmost; j+=lag) {
            if (partic[j]>max) {
               max = partic[j];
               maxJ = j;
            }
        }
//      imprima particoes mais provaveis (ignorando as de ocorrencia)
        if ( partic[maxJ] > 0 ) {
           for (j=1; j<=tamU; j++) {
               fprintf(filePartic, "%d ", vecU[maxJ][j]);
           }
           fprintf(filePartic, "\t %f\t %d\n",
                   float(partic[maxJ])/numAmostU, partic[maxJ] );
           partic[maxJ] = -1;
        }
    }
    return 0;
};

int ChangePoints::ComputePbPost(void) {
    fprintf(stdout, "ChangePoints::ComputePbPost\n");
    int i, j;
    blocks = new int[numAmost+1];
    for (i=1; i<=numAmost; i++) {
        blocks[i] = 1;
    }
    for (i=burnIn+1; i<=numAmost; i+=lag) {
        if ( blocks[i] != -1 ) {
           for (j=i+lag; j<=numAmost; j+=lag) {
               // comparar i-esima amostra com j-esima amostra
               if ( numBlocks[i]==numBlocks[j] ) {
                  blocks[i]++;
                  blocks[j] = -1;
               }
           }
        }
    }
//    for (i=burnIn+1; i<=numAmost; i+=lag) {
//        fprintf(arqNumBlo, "%d\n", blocks[i]);
//    }
    return 0;
};

int ChangePoints::ComputePbYaoPost(void) {
   fprintf(stdout, "ChangePoints::ComputePbYaoPost\n");
   // compute the posterior distribution of number of blocks
   ComputePbPost();
   return 0;
}

int ChangePoints::PrintPbPost(void) {
    fprintf(stdout, "ChangePoints::PrintPbPost\n");
    int i, j, minJ, min;
    fprintf(fileNumBlo, "#blocks\t #occur\t prob\n" );
    for (i=burnIn+1; i<=numAmost; i+=lag) {
        // find minimum
        min = tamX+1;
        for (j=burnIn+1; j<=numAmost; j+=lag) {
            if ( (numBlocks[j]<min)&&(blocks[j]!=-1)) {
               min = numBlocks[j];
               minJ = j;
            }
        }
        // print it
        if ((blocks[minJ]!=-1)&&(numBlocks[minJ]!=(tamX+1))) {
           fprintf(fileNumBlo, "%d\t %d\t %f\n",
                  numBlocks[minJ], blocks[minJ], float(blocks[minJ])/numAmostU );
        }
        numBlocks[minJ] = tamX+1;
    }
    return 0;
}

int ChangePoints::Yao(void) {
   fprintf(stdout, "ChangePoints::Yao\n");
   StartCronograph();
   // compute mij*, aij*, dij*, vij*, and fXij
   ComputeParamPost();
   // compute posterior cohesions
   ComputeCohesPost();
   // compute posterior relevance for each block [ij]
   ComputeRelevYaoPost();
   // compute muk and sigmak^2
   ComputeMuSigYaoPost();
   ElapsedTime();
   // compute probability of each partition
   ComputeRhoYaoPost();
   // compute the posterior distribution of number of blocks
   ComputePbYaoPost();
   // print results
   PrintResults();
   return 0;
}

int ChangePoints::PrintResults(void) {
   fprintf(stdout, "ChangePoints::PrintResults\n");
   // print results
   PrintMuSigPost();
   PrintNumBloc();
   PrintRhoPost();
   PrintPbPost();
   return 0;
}

#endif