//
// This library implements several procedures
// to obtain product estimates
//
// Authors: Rosangela H. Loschi
//          Frederico R. B. Cruz
//          Departamento de Estatistica
//          Universidade Federal de Minas Gerais
//          E-mail: {loschi,fcruz}@est.ufmg.br
//
// Version: 1.00
//
// Date: Feb/2001
//

#ifndef CHPOINT_H
#define CHPOINT_H

// inclusion of standard libraries
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

// uniform generator
#include "uni.c"

// gamma function
#include "gamma.cpp"

#define PI (3.141592654)

// prior parameters of mu
#define mPrior(i,j) (0.0)
#define vPrior(i,j) (1.0)

// prior parameters of sqr(sigma)
#define dPrior(i,j) (4.0)
#define aPrior(i,j) (0.01)


// 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 probP;
  // prior parameters of probability p
  float alpha;
  float beta;
  // posterior parameters of mu
  float **mPost;
  float **vPost;
  // posterior parameters of sqr(sigma)
  float **dPost;
  float **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;
  // posterior relevance of the block [ij]
  int **relevPost;
  // posterior estimates of mu and sqr(sigma)
  float *muPost;
  float *sigma2Post;
  // private methods
  int ComputePostParam(void);
  int PrintPostParam(void);
  int ComputeVecUbeta(void);
  int ComputeNuBloc(void);
  int PrintNumBloc(void);
  int ComputePostRelev(void);
  int PrintPostRelev(void);
  int ComputePostMuSig(void);
  int PrintPostMuSig(void);
  int ComputePostRho(void);
  int PrintPostRho(void);
  int ComputePostPb(void);
  int PrintPostPb(void);
public:
  ChangePoints(void);    // default constructor
  ~ChangePoints(void);   // destructor
  int ReadInputScript(char *scriptFileName);
  int LoGibbsBeta(void);
};

// class implementation

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

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
    i = 0;
    while ( (scriptFileName[i] != '\n')&&(i<LLENGTH) ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "Script file name: %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: \t%d\n", numAmost);
    fgets(Line, LLENGTH, scriptFile);
    fscanf(scriptFile, "%d\n", &burnIn);
    fprintf(stdout, "burnIn: \t%d\n", burnIn);
    fgets(Line, LLENGTH, scriptFile);
    fscanf(scriptFile, "%d\n", &lag);
    fprintf(stdout, "Lag: \t%d\n", lag);
//
    fgets(Line, LLENGTH, scriptFile);
    fgets(Line, LLENGTH, scriptFile);
    fscanf(scriptFile, "%f\n", &alpha);
    fprintf(stdout, "Alpha: \t%f\n", alpha);
    fgets(Line, LLENGTH, scriptFile);
    fscanf(scriptFile, "%f\n", &beta);
    fprintf(stdout, "Beta: \t%f\n", beta);
//
    fgets(Line, LLENGTH, scriptFile);
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ( (fileName[i] = fgetc(scriptFile)) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for data: %s\n",fileName);
    fileData = fopen(fileName, "r");
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ( (fileName[i] = fgetc(scriptFile)) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for number of blocks: %s\n",fileName);
    fileBurnIn = fopen(fileName, "w");
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ( (fileName[i] = fgetc(scriptFile)) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for posterior estimates of mu & sqr(sigma): %s\n",fileName);
    fileMuSig = fopen(fileName, "w");
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ( (fileName[i] = fgetc(scriptFile)) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for posterior distribution of number of blocks: %s\n",fileName);
    fileNumBlo = fopen(fileName, "w");
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ( (fileName[i] = fgetc(scriptFile)) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for the most probable partitions: %s\n",fileName);
    filePartic = fopen(fileName, "w");
    fgets(Line, LLENGTH, scriptFile);
    fscanf(scriptFile, "%d/n", &numPartic);
    fclose(scriptFile);
    // count time series size
    tamX = 0;
    while (fscanf(fileData,"%f\n", &dummy) == 1) {
          tamX++;
//          fprintf(stdout, "dummy = %f\t tamX = %d\n", dummy, tamX);
    }
    // read time series
    rewind(fileData);
    fprintf(stdout, "Time series size: %d\n", tamX);
    X = new float[tamX+1];
    tamU = tamX-1;
    for (i=1; i <= tamX; i++) {
        fscanf(fileData, "%f\n", &X[i]);
//        fprintf(stdout, "X%d = %f\n", i, X[i]);
    }
    return 0;
}

int ChangePoints::ComputePostParam(void) {
    fprintf(stdout,"ChangePoints::ComputePostParam\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]);
      }
    }
    return 0;
};

int ChangePoints::PrintPostParam(void) {
    int i, j;
    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::ComputeVecUbeta(void) {
    fprintf(stdout, "ChangePoints::ComputeVecUbeta\n");
    int k, i, numBlocks, 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++) {
            // compute number of blocks
            numBlocks = 1;
            for (l=1; l<=tamU; l++) {
                if (vecU[k][l] == 0) numBlocks++;
            }
            // 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(fXij[x][y]) - log(fXij[x][i]) - log(fXij[i][y]) +
                 log(Gamma(numBlocks+alpha-2)) +
                 log(Gamma(tamX+beta-numBlocks+1)) -
                 log(Gamma(numBlocks+alpha-1)) -
                 log(Gamma(tamX+beta-numBlocks));
            Ri = exp(Ri);
//            fprintf(stdout, "fXij[%d][%d] = %f\n", x, y, fXij[x][y]);
//            fprintf(stdout, "fXij[%d][%d] = %f\n", x, i, fXij[x][i]);
//            fprintf(stdout, "fXij[%d][%d] = %f\n", i, y, fXij[i][y]);
//            fprintf(stdout, "Ri= (%f/(%f*%f))*(%f*%f)/(%f*%f)=%f\n",
//                    fXij[x][y], fXij[x][i], fXij[i][y],
//                    Gamma(numBlocks+alpha-2), Gamma(tamX+beta-numBlocks+1),
//                    Gamma(numBlocks+alpha-1), Gamma(tamX+beta-numBlocks), Ri);                    tamX, beta, numBlocks, fXij[x][i],
            // generate Uki
            u = UNI;
            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::ComputePostRelev(void) {
    fprintf(stdout, "ChangePoints::ComputePostRelev\n");
    int i, j, inicBloc, finalBloc;
    relevPost = new int*[tamX];
    // creat matrix
    for (i=0; i<tamX; i++) {
        relevPost[i] = new int[tamX+1];
        for (j=i+1; j<=tamX; j++)
            relevPost[i][j] = 0;
    }
    // make counts
    numAmostU = 0;
    for (i=burnIn+1; i<=numAmost; i+=lag) {
        numAmostU++;
        inicBloc = 0;
        for (j=1; j<tamX; j++) {
            if (vecU[i][j] == 0) {
               finalBloc = j;
               relevPost[inicBloc][finalBloc]++;
               inicBloc = j;
            }
        }
        relevPost[inicBloc][tamX]++;
    }
    fprintf(stdout, "Number of net samples: %d\n",numAmostU);
    return 0;
}

int ChangePoints::PrintPostRelev(void) {
    fprintf(stdout, "ChangePoints::PrintPostRelev\n");
    int i,j;
    for (i=0; i<tamX; i++) {
        for (j=1; j<=i; j++)
            fprintf(stdout, "  ");
        for (j=i+1; j<=tamX; j++)
            fprintf(stdout, "%d ", relevPost[i][j]);
        fprintf(stdout, "\n");
    }
    return 0;
};


int ChangePoints::ComputePostMuSig(void) {
    fprintf(stdout, "ChangePoints::ComputePostMuSig\n");
    int k, i, j;
    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++) {
                muPost[k] += relevPost[i][j] * mPost[i][j] / numAmostU;
                sigma2Post[k] += relevPost[i][j] * aPost[i][j] /
                                 ((dPost[i][j] - 2)*numAmostU);
            }
        }
    }
    return 0;
}

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

int ChangePoints::ComputePostRho(void) {
    fprintf(stdout, "ChangePoints::ComputePostRho\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::PrintPostRho(void) {
    fprintf(stdout, "ChangePoints::PrintPostRho\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::ComputePostPb(void) {
    fprintf(stdout, "ChangePoints::ComputePostPb\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::PrintPostPb(void) {
    fprintf(stdout, "ChangePoints::PrintPostPb\n");
    int i, j, minJ, min;
    fprintf(fileNumBlo, "#blocks\t #occur\t \tprob\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::LoGibbsBeta(void) {
    fprintf(stdout, "ChangePoints::LoGibbsBeta\n");
   // compute mij*, aij*, dij*, vij*, and fXij
   ComputePostParam();
//  PrintPostParam();
   // generate samples u
   ComputeVecUbeta();
   // compute number of blocks for each sample
   ComputeNuBloc();
   PrintNumBloc();
   // compute posterior relevance for each block [ij]
   ComputePostRelev();
//  PrintPostRelev();
   // compute muk and sqr(sigmak)
   ComputePostMuSig();
   PrintPostMuSig();
   // compute probability of each partition
   ComputePostRho();
   PrintPostRho();
   // compute the posterior distribution of number of blocks
   ComputePostPb();
   PrintPostPb();
   return 0;
}

#endif