//
// Purpose:
//    This library implements a procedure
//    to obtain product estimates to a regression
//
// 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:
//    Jan/2003
//

#ifndef CHREG_H
#define CHREG_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"

#define float double

// prior parameters of alpha
#define aPrior(i,j) (0.0)
#define tau02 (1.0)

// prior parameters of psi
#define zPrior(i,j) (0.0)
#define gamma02 (1.0)

// prior parameters of sigma^2
#define nPrior(i,j) (1.00)
#define dPrior(i,j) (4.00)

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 ChangeReg {
  FILE *fileDataX;
  FILE *fileDataY;
  FILE *fileBurnIn;
  FILE *fileAlpPsiSig;
  FILE *fileNumBlo;
  FILE *filePartic;
  // maximum number of partition to be printed
  int numPartic;
  // time series and their sizes
  float *X;
  int tamX;
  float *Y;
  int tamY;
  // probability p
  float pProb;
  // posterior parameters of alpha
  float **aPost, **dPost, **vPostAlp, **nPostAlp;
  // posterior parameters of psi
  float **zPost, **vPostPsi, **nPostPsi;
  // posterior parameters of sigma^2
  float **nPostSig;
  // predictive distribution
  float **fYijXij, **vPostY, **nPostY;
  // 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]
  float **relevPost;
  // posterior estimates of alpha, psi, and sigma^2
  float *alpPost, *psiPost, *sigma2Post;
  // private methods
  int ComputeParamPost(void);
  int ComputeCohesPost(void);
  int ComputeVecU(void);
  int ComputeNuBloc(void);
  int PrintNumBloc(void);
  int ComputeRelevPost(void);
  int ComputeAlpPsiSigPost(void);
  int PrintAlpPsiSigPost(void);
  int ComputeRhoPost(void);
  int PrintRhoPost(void);
  int ComputePbPost(void);
  int PrintPbPost(void);
public:
  ChangeReg(void);    // default constructor
  ~ChangeReg(void);   // destructor
  int ReadInputScript(char *scriptFileName);
  int LoGibbs(void);
  int PrintResults();
};

// class implementation

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

ChangeReg::~ChangeReg(void) {
    fprintf(stdout, "Object ChangeReg destroyed()\n");
}

int ChangeReg::ReadInputScript(char *scriptFileName) {
    fprintf(stdout,"ChangeReg::ReadInputScript(char)\n");

    const int LLENGTH = 256;
    char Line[LLENGTH];
    int i;
    FILE *scriptFile;
    char fileName[LLENGTH];
    float dummy;

//    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 X\n\t%s\n",fileName);
    fileDataX = fopen(fileName, "r");
    fgets(Line, LLENGTH, scriptFile);
    i = 0;
    while ((fileName[i]=Line[i]) != '\n' ) i++;
    fileName[i] = '\0';
    fprintf(stdout, "File name for data Y\n\t%s\n",fileName);
    fileDataY = 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\n\t%s\n",fileName);
    fileAlpPsiSig = 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 X size
    tamX = 0;
    while (fscanf(fileDataX,"%lf\n", &dummy) == 1) {
          tamX++;
//          fprintf(stdout, "dummy = %lf\t tamX = %d\n", dummy, tamX);
    }
    // read time series
    rewind(fileDataX);
    fprintf(stdout, "Time series X size\n\t%d\n", tamX);
    X = new float[tamX+1];
    tamU = tamX-1;
    for (i=1; i <= tamX; i++) {
        fscanf(fileDataX, "%lf\n", &X[i]);
//        fprintf(stdout, "X%d = %lf\n", i, X[i]);
    }
//  count time series Y size
    tamY = 0;
    while (fscanf(fileDataY,"%lf\n", &dummy) == 1) {
          tamY++;
//          fprintf(stdout, "dummy = %lf\t tamY = %d\n", dummy, tamY);
    }
    // read time series
    rewind(fileDataY);
    fprintf(stdout, "Time series Y size\n\t%d\n", tamY);
    Y = new float[tamY+1];
    tamU = tamY-1;
    for (i=1; i <= tamY; i++) {
        fscanf(fileDataY, "%lf\n", &Y[i]);
//        fprintf(stdout, "Y%d = %lf\n", i, Y[i]);
    }
    return 0;
}

int ChangeReg::ComputeParamPost(void) {
    fprintf(stdout,"ChangeReg::ComputeParamPost()\n");

    int i, j, k;
    float sumX, sumY, sumXY, sumXX, sumYY, W;
    aPost = new float*[tamX];
    dPost = new float*[tamX];
    vPostAlp = new float*[tamX];
    nPostAlp = new float*[tamX];
    zPost = new float*[tamX];
    vPostPsi = new float*[tamX];
    nPostPsi = new float*[tamX];
    nPostSig = new float*[tamX];
    fYijXij = new float*[tamX];
    vPostY = new float*[tamX];
    nPostY = new float*[tamX];

    for (i=0; i<tamX; i++) {
      // create matrixes
      aPost[i] = new float[tamX+1];
      dPost[i] = new float[tamX+1];
      vPostAlp[i] = new float[tamX+1];
      nPostAlp[i] = new float[tamX+1];
      zPost[i] = new float[tamX+1];
      vPostPsi[i] = new float[tamX+1];
      nPostPsi[i] = new float[tamX+1];
      nPostSig[i] = new float[tamX+1];
      fYijXij[i] = new float[tamX+1];
      vPostY[i] = new float[tamX+1];
      nPostY[i] = new float[tamX+1];
      for (j=i+1; j<=tamX; j++) {
          // compute average within block [ij]
          sumX = 0.0; sumY = 0.0; sumXY = 0.0; sumXX = 0.0; sumYY = 0.0;
          for (k=i+1; k<=j; k++) {
              sumX += X[k];
              sumY += Y[k];
              sumXY += X[k]*Y[k];
              sumXX += X[k]*X[k];
              sumYY += Y[k]*Y[k];
          }
          // compute posterior parameters
          aPost[i][j] = (
                         (gamma02*sumXX+1)*(aPrior(i,j)+tau02*sumY)-
                         tau02*sumX*(zPrior(i,j)+gamma02*sumXY)
                        )/(
                          gamma02*sumXX+1+(j-i)*tau02*(gamma02*sumXX+1)-
                          tau02*gamma02*sumX*sumX
                        );
          dPost[i][j] = dPrior(i,j) + (j-i);
          vPostAlp[i][j] = (
                             tau02*(1+gamma02*sumXX)
                           )/(
                             (1+tau02*(j-i))*(gamma02*sumXX+1)-
                             gamma02*tau02*sumX*sumX
                           );
          nPostAlp[i][j] = nPrior(i,j) + (
                            (gamma02*sumXX+1)*(
                              aPrior(i,j)*aPrior(i,j)*gamma02 +
                              tau02*gamma02*sumYY +
                              zPrior(i,j)*zPrior(i,j)*tau02
                            )
                           )/(
                             tau02*gamma02*(gamma02*sumXX+1)) - (
                                gamma02*sumXY*sumXY +
                                2*zPrior(i,j)*sumXY +
                                zPrior(i,j)*zPrior(i,j)
                             )/(gamma02*sumXX+1) -
                             pow((
                               (gamma02*sumXX+1)*(aPrior(i,j)+tau02*sumY) -
                               tau02*sumX*(gamma02*sumXY+zPrior(i,j)
                             )),2)/(
                               ((gamma02*sumXX+1)*(1+tau02*(j-i)) -
                                 tau02*gamma02*sumX*sumX
                               )*
                               (tau02*(gamma02*sumXX+1))
                           );
          zPost[i][j] = ((
                          1+tau02*(j-i))*(gamma02*sumXY+zPrior(i,j)) -
                          gamma02*sumX*(tau02*sumY+aPrior(i,j))
                        )/(
                          (1+tau02*(j-i))*(gamma02*sumXX+1) -
                          gamma02*tau02*sumX*sumX
                        );
          vPostPsi[i][j] = (
                             gamma02*(1+tau02*(j-i))
                           )/(
                             (1+tau02*(j-i))*(gamma02*sumXX+1) -
                             gamma02*tau02*sumX*sumX
                           );
          W = (
               (1+tau02*(j-i))*(gamma02*sumXX+1)-gamma02*tau02*sumX*sumX
              )/(
               gamma02*(1+tau02*(j-i))
              );
          nPostPsi[i][j] = nPrior(i,j) +
                           W*(
                             (1+tau02*(j-i))*(tau02*gamma02*sumYY +
                              aPrior(i,j)*aPrior(i,j)*gamma02 +
                              zPrior(i,j)*zPrior(i,j)*tau02
                             ) - (
                              gamma02*pow((tau02*sumY+aPrior(i,j)),2)
                             )
                           )/(
                             tau02*(1+tau02*(j-i))*(gamma02*sumXX+1)-
                             gamma02*tau02*sumX*sumX
                           ) -
                           W*pow((
                              (
                                (1+tau02*(j-i))*(gamma02*sumXY+zPrior(i,j)) -
                                gamma02*sumX*(tau02*sumY+aPrior(i,j))
                              )/(
                                (1+tau02*(j-i))*(gamma02*sumXX+1)-
                                gamma02*tau02*sumX*sumX
                              )
                           ),2);
          nPostSig[i][j] = nPrior(i,j) + (
                            (1+tau02*(j-i))*(tau02*zPrior(i,j)*zPrior(i,j)+
                             tau02*gamma02*sumYY+
                             gamma02*aPrior(i,j)*aPrior(i,j)
                            )/(
                             tau02*gamma02*(1+tau02*(j-i))
                            )
                           ) - (
                             tau02*tau02*gamma02*sumY*sumY+
                             gamma02*aPrior(i,j)*aPrior(i,j)+
                             2*aPrior(i,j)*tau02*gamma02*sumY
                           )/(
                             tau02*gamma02*(1+tau02*(j-i))
                           ) - (
                             pow((
                                (1+tau02*(j-i))*(zPrior(i,j)+gamma02*sumXY)-
                                gamma02*tau02*sumY*sumX-
                                gamma02*aPrior(i,j)*sumX
                             ),2)
                           )/(
                             gamma02*(1+tau02*(j-i))*((1+tau02*(j-i))*
                               (gamma02*sumXX+1)-
                               gamma02*tau02*sumX*sumX
                             )
                           );
          // compute predictive distribution
          vPostY[i][j] = (
                           tau02*(1+gamma02*sumXX)*gamma02*
                           (1+(j-i)*tau02)-
                           tau02*tau02*gamma02*gamma02*sumX*sumX
                         )/(
                           pow((
                             (1+(j-i)*tau02)*(1+gamma02*sumXX)-
                             tau02*gamma02*sumX*sumX
                           ),2)
                         );
          nPostY[i][j] = nPrior(i,j) +
                         aPrior(i,j)*aPrior(i,j)/tau02 +
                         zPrior(i,j)*zPrior(i,j)/gamma02 +
                         sumYY -
                         aPost[i][j]*aPost[i][j]*(1+(j-i)*tau02)/tau02 -
                         2*aPost[i][j]*zPost[i][j]*sumX -
                         zPost[i][j]*zPost[i][j]*(1+gamma02*sumXX)/gamma02;
          fYijXij[i][j] = (
                           sqrt(vPostY[i][j])*pow(nPrior(i,j),(dPrior(i,j)/2))*
                           Gamma(dPost[i][j]/2)*
                           pow(nPostY[i][j],(-dPost[i][j]/2))
                          )/(
                            sqrt(gamma02*tau02)*pow(PI,((j-i)/2))*
                            Gamma(dPrior(i,j)/2)
                          );
      }
    }
    // print results
//    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, "vPostAlp:\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 ", vPostAlp[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "nPostAlp:\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 ", nPostAlp[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "zPost:\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 ", zPost[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "vPostPsi:\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 ", vPostPsi[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "nPostPsi:\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 ", nPostPsi[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "nPostSig:\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 ", nPostSig[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "vPosY:\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 ", vPostY[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "nPostY:\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 ", nPostY[i][j]);
//        fprintf(stdout, "\n");
//    }
//    fprintf(stdout, "fYijXij:\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 ", fYijXij[i][j]);
//        fprintf(stdout, "\n");
//    }
    return 0;
};

int ChangeReg::ComputeCohesPost(void) {
    fprintf(stdout,"ChangeReg::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(fYijXij[i][j]);
          cPost[i][j] = exp(cPost[i][j]);
      }
      cPost[i][tamX] = (tamX-i-1)*log(1-pProb) + log(fYijXij[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 ChangeReg::ComputeVecU(void) {
    fprintf(stdout, "ChangeReg::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 ChangeReg::ComputeNuBloc(void) {
   fprintf(stdout, "ChangeReg::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 ChangeReg::PrintNumBloc() {
    fprintf(stdout, "ChangeReg::PrintNumBloc()\n");
    for (int k=1; k<=numAmost; k++) {
        fprintf(fileBurnIn,"%d\n", numBlocks[k]);
    }
    return 0;
}

int ChangeReg::ComputeRelevPost(void) {
    fprintf(stdout, "ChangeReg::ComputeRelevPost()\n");
    int i, j, inicBloc, finalBloc;
    relevPost = new float*[tamX];
    // creat matrix
    for (i=0; i<tamX; i++) {
        relevPost[i] = new float[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);
    // 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 ChangeReg::ComputeAlpPsiSigPost(void) {
    fprintf(stdout, "ChangeReg::ComputeAlpPsiSigPost()\n");
    int k, i, j;
    alpPost = new float[tamX+1];
    psiPost = new float[tamX+1];
    sigma2Post = new float[tamX+1];
    for (k=1; k<=tamX; k++) {
        alpPost[k] = 0.0;
        psiPost[k] = 0.0;
        sigma2Post[k] = 0.0;
        for (i=0; i<=k-1; i++) {
            for (j=k; j<=tamX; j++) {
                alpPost[k] += relevPost[i][j] * aPost[i][j] / numAmostU;
                psiPost[k] += relevPost[i][j] * zPost[i][j] / numAmostU;
                sigma2Post[k] += relevPost[i][j]*nPostSig[i][j]/
                                 ((dPost[i][j] - 2)*numAmostU);
            }
        }
    }
    return 0;
}

int ChangeReg::PrintAlpPsiSigPost(void) {
    fprintf(stdout, "ChangeReg::PrintAlpPsiSigPost()\n");
    fprintf(fileAlpPsiSig, "alpPost\t psiPost\t sigma2Post\n");
    for (int i=1; i<=tamX; i++) {
        fprintf(fileAlpPsiSig, "%f\t %f\t %f\n", alpPost[i], psiPost[i], sigma2Post[i]);
    }
    return 0;
}

int ChangeReg::ComputeRhoPost(void) {
    fprintf(stdout, "ChangeReg::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 ChangeReg::PrintRhoPost(void) {
    fprintf(stdout, "ChangeReg::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 ChangeReg::ComputePbPost(void) {
    fprintf(stdout, "ChangeReg::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 ChangeReg::PrintPbPost(void) {
    fprintf(stdout, "ChangeReg::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 ChangeReg::LoGibbs(void) {
   fprintf(stdout, "ChangeReg::LoGibbs()\n");
   // compute aij*, dij*, VijAlp*, nijAlp*,
   // zij*, VijPsi*, nijPsi*, nijSig*, VijY*, nijY*, fYijXij
   ComputeParamPost();
   // compute posterior cohesions cij*
   ComputeCohesPost();
   // generate samples u
   ComputeVecU();
   // compute number of blocks for each sample
   ComputeNuBloc();
   // compute posterior relevance for each block [ij]
   ComputeRelevPost();
   // compute the product estimates of alpha, psi, sigma^2
   ComputeAlpPsiSigPost();
   // compute probability of each partition
   ComputeRhoPost();
   // compute the posterior distribution of number of blocks
   ComputePbPost();
   // print results
   PrintResults();
   return 0;
}

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

#endif