/* ---------------------------------------------------------------------------
MrAICM Utility.
written by C R Young
last updated 9/24/07

    Copyright (C) 2007  C R Young

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.


compile as:
gcc AICM.c -o AICM -lm    
 ---------------------------------------------------------------------------*/


#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#define MAXSTR 80
#define MAXLINE 2000
#define MAXMODELS 10
#define MAXFIELD 100
#define MAXFILES 20

struct Model_quantities {
    char name[MAXSTR];                   // name of the model (e.g., "unconstrained")
    char filelist[MAXFILES][MAXSTR];     // list of filenames for this model (e.g., "sym_16S_run1.p")
    int Numfiles;                        // number of files for this model
    int MCMCdraws[MAXFILES];             // number of MCMC draws in each file (calculated on file read)
    int burnin;							// burnin specified for this model
    int thin;							// thinning interval for this model
    double data[MAXFILES][MAXLINE];      // actual loglikelihood scores
    int totaldraws;                      // sum of the MCMC draws across all files minus the burnin
    int Alldraws;						// sum of the MCMC draws across all files
    double dhat;                         // effefctive number of model parameters
    double mean_loglik;                  // mean of the loglikelihood scores
    double var_loglik;                   // variance of the loglikelihood scores
    double skew_loglik;					// skewedness of the loglikelihood scores
    double kurt_loglik;					// kurtosis of the loglikelihood scores
    double AICM;                         // MCMC estimate of AIC
    double SE_AICM;                      // Standard error of AICM (assuming approximately independent MCMC draws!!!!)
    double harmonicmean;                 // Harmonic mean of the raw likelihood scores
    double AkaikeWeight;
    double MaxLogL;
    };

/* My error message */
void MyError(error_text)
char  error_text[];
{
	fprintf(stderr,"\nRY's run-time error...\n");
	fprintf(stderr,"\n\t%s\n",error_text);
	fprintf(stderr,"\n...now exiting to system...\n\n");
	exit(1);
}



void Moment(int model, struct Model_quantities ModelData[])
{
    //void MyError(char *error_text);
    int j, n, file, line;
    double ep=0.0, s, p, ave, var, adev, sdev, skew, kurt;
    
    n = ModelData[model].totaldraws;
    
    
    if(n <= 1) MyError("Error in Moment(): n can't be less than 2 to compute moments.\n");
    
    // Find the mean
    s = 0.0;
    for(file=0; file<ModelData[model].Numfiles; file++){
        for(line=0; line<ModelData[model].MCMCdraws[file]; line++){
            s += ModelData[model].data[file][line];  
        }
    }
    ave = s / n;
    ModelData[model].mean_loglik = ave;
    
    
    // Find the higher moments
    for(file=0; file<ModelData[model].Numfiles; file++){
        for(line=0; line<ModelData[model].MCMCdraws[file]; line++){
            adev += fabs(s=ModelData[model].data[file][line]-ave);
            ep += s;  
            var += (p=s*s);
            skew += (p *= s);
            kurt += (p *= s);
        }
    }
    adev /= n;
    var = (var - ep*ep/n)/(n-1);   // corrected two pass formula
    sdev = sqrt(var);
    if (var) {
        skew /= (n*var*sdev);
        kurt = kurt/(n*var*var) - 3.0;
    } else {
         MyError("Error in Moment(): variance = 0, can't compute higher moments.\n");
    }
    ModelData[model].var_loglik = var;
    ModelData[model].skew_loglik = skew;
    ModelData[model].kurt_loglik = kurt;

}




main (int argc, char *argv[])
{


        FILE *Files, *currentfile;
        char fname[MAXSTR];
        struct Model_quantities ModelData[MAXMODELS];
        int i, j, model, file, Nummodels, Numfiles[MAXFILES], BestModel, thin, Lastdraw, countsample;
        int end, countline, field, line, datasize[MAXMODELS][MAXFILES], linecount[MAXMODELS];
        double data[MAXMODELS][10000][2], dummynumber, SumX, Sum_Xsq, SumX_sq, firstterm, secondterm, thirdterm, shiftedgam; 
        char filelist[MAXMODELS][MAXFILES][MAXSTR];
        char dummystring[MAXLINE], check[2], strtest[MAXSTR];
        
        
/*********************/
/* Begin read files  */
/*********************/

        if( argc > 2 || argc < 2) {
            printf("usage: MrAICM [infile]\n");
            exit(1);
        }

        
        /* Read file list */
        
        Files = fopen(argv[1], "r");
        if (!Files)
        {
            printf ( "\nCould not find files.txt...\n\n");
            MyError("No files specified.\n");
        } else {
            printf("\nReading file list...\n");
            fscanf(Files,"%d", &Nummodels);  // Read the number of hypotheses to be compared
            printf("%d models to compare.\n", Nummodels);
            if(Nummodels > MAXMODELS) MyError("Too many models to compare. Increase MAXMODELS.\n");
            
            for(model=0; model<Nummodels; model++){
                printf("----------------------------------\n");
                printf("Informaiton for model %d:\n", model+1);
                printf("----------------------------------\n");
                fscanf(Files,"%d", &Numfiles[model]);  // Read a list of the numbers of MCMC files to open and read for each hypothesis
                ModelData[model].Numfiles = Numfiles[model]; 
                printf("%d files for model %d.\n", Numfiles[model], model+1);
                if(Numfiles[model] > MAXFILES) {
                    printf ( "\nToo many files for model %d\n\n", i);
                    MyError("Too many files. Increase MAXFILES.\n");
                }
                fscanf(Files,"%d %d", &ModelData[model].burnin, &ModelData[model].thin);  // Read burnin and thinning interval for each model
                printf("Burnin (applied to each file) is %d for model %d.\n", ModelData[model].burnin, model+1);
                printf("Thinning interval (applied to each file) is %d for model %d.\n\n", ModelData[model].thin, model+1);
            }
            printf ("\n");
            
            for(model=0; model<Nummodels; model++){
                printf ("Files for model %d:\n", model+1);
                for(file=0; file<Numfiles[model]; file++){
                    fscanf(Files,"%s", &ModelData[model].filelist[file]);
                    printf ("%s\n", ModelData[model].filelist[file]);
                }
                printf ("\n");
            }


        if(fclose(Files) != 0) MyError("Could not close file: files.txt.\n");
        
        } // END else condition

/*********************************************/        
/* Main loop reading through MCMC data       */
/*********************************************/ 
        
    for(model=0; model<Nummodels; model++){
                printf("------------------\n");
                printf("Scanning model %d:\n", model+1);
                printf("------------------\n");
       for(file=0; file<Numfiles[model]; file++){
            /* Open file for reading */
                currentfile =  fopen (ModelData[model].filelist[file], "r");
        
            /* Error trappiong for file opening */
            if (!currentfile)
            {
                printf ( "Error opening file for model %d file number %d!\n", model, file);
                MyError("Could not read MCMC file\n");
            }
        

            fgets(dummystring, MAXLINE, currentfile); // remove id header
            fgets(dummystring, MAXLINE, currentfile); // remove quantity header
            
            end = 0;
            countline = 1;
            countsample = 1;

            /* until you reach the end of this file... */
            while( end == 0){ 	
                /* Read line */
                if(end == 0) {
                
                        if(countline < ModelData[model].burnin + 1){
                            fgets(dummystring, MAXLINE, currentfile); // remove line without storing
                            countline += 1; // count this line that is not stored
                        } else {                                       
                                        
                            fscanf(currentfile, "%lf", &dummynumber); // remove generation
                            fscanf(currentfile, "%lf", &ModelData[model].data[file][countsample-1]); // store loglikelihood
                            fgets(dummystring, MAXLINE, currentfile); // remove rest of line
                            countline += 1; // count this line that is stored
                            countsample += 1; // count this sample
                            
                            // remove thinned samples
                            for(thin=1; thin<ModelData[model].thin; thin++){
                                fgets(dummystring, MAXLINE, currentfile); // remove lines for thinning
                                /* Force end of file read, then push chars back into stack */
                                for(i=0; i<2; i++){
                                    check[i] = getc(currentfile);
                                    end = feof(currentfile);
                                }
                                for(i=1; i>-1; i--){
                                    ungetc(check[i],currentfile); 
                                }
                                if(end == 0) {
                                    countline += 1;  // count these lines that are not stored
                                }
                            }

                        }
                    
                    
                    
                } /* END IF ENDFILE */
            } /* END WHILE ENDFILE */
            
            ModelData[model].MCMCdraws[file] = countsample - 1;
            printf("Number of draws in MCMC file: %d\n", countline - 1);
            printf("Number of draws after thinning and burnin: %d\n", ModelData[model].MCMCdraws[file]);

            /* Close all files */
            if(fclose(currentfile) != 0) MyError("Could not close file: INFILE.\n");

       } // END file loop
            
            ModelData[model].totaldraws = 0;
            for(file=0; file<ModelData[model].Numfiles; file++){
                ModelData[model].totaldraws += ModelData[model].MCMCdraws[file];
            }
            
            printf("Total MCMC draws in model %d (after burnin and thinning) = %d\n\n", model+1, ModelData[model].totaldraws);
            
    } // END model loop


/*********************************************/        
/* Compute AICM from the MCMC data           */
/*********************************************/ 

    // calculate moments of loglikelihood distribution
    for(model=0; model<Nummodels; model++){
        Moment(model, ModelData);
    }


    for(model=0; model<Nummodels; model++){
                        
         // calculate AICM
        ModelData[model].AICM = -2.0 * (ModelData[model].mean_loglik - ModelData[model].var_loglik);
        
        // calculate dhat (effective number of parameters)
        ModelData[model].dhat = 2.0 * ModelData[model].var_loglik;  
        
        // calculate SE of AICM : page 17 of Raferty et al 2007 Bayesian statistics 8, pp. 1-45
        firstterm = 4.0 * ModelData[model].dhat / (2.0 * ModelData[model].totaldraws);
        secondterm = 4.0 * ModelData[model].dhat * (11.0*ModelData[model].dhat/4.0 + 12.0) / ModelData[model].totaldraws;
        ModelData[model].SE_AICM = sqrt( firstterm + secondterm );  
        
        
    } // END model loop
    
    // rank the models
    BestModel = 0;
    for(model=1; model<Nummodels; model++){
        if(ModelData[model].AICM < ModelData[model-1].AICM)
        BestModel = model;
    }
    
    thirdterm = 0.0;
    for(model=0; model<Nummodels; model++){
        firstterm = ModelData[model].AICM - ModelData[BestModel].AICM;
        secondterm = -0.5 * firstterm;
        thirdterm += exp(secondterm);
    }
    
    for(model=0; model<Nummodels; model++){
        if(model != BestModel){
            ModelData[model].AkaikeWeight = exp(-0.5 * (ModelData[model].AICM - ModelData[BestModel].AICM)) / thirdterm; 
        } else {
            ModelData[model].AkaikeWeight = exp(-0.5 * (ModelData[model].AICM - ModelData[BestModel].AICM)) / thirdterm;
        }
    }
    


/*********************************************/        
/* Print results                             */
/*********************************************/ 
    
    printf("Model\tAkaike Weight\tAICM\t\tSE_AICM\t\td_hat\t\tE(logL)\t\tVar(logL)\n");
    printf("%d\t%.5f\t\t%.1f\t\t%.1f\t\t%.1f\t\t%.1f\t\t%.1f\n", BestModel+1, ModelData[BestModel].AkaikeWeight, ModelData[BestModel].AICM, ModelData[BestModel].SE_AICM, ModelData[BestModel].dhat, ModelData[BestModel].mean_loglik, ModelData[BestModel].var_loglik);

    
    for(model=0; model<Nummodels; model++){
        if(model != BestModel){
            printf("%d\t%.5f\t\t%.1f\t\t%.1f\t\t%.1f\t\t%.1f\t\t%.1f\n", model+1, ModelData[model].AkaikeWeight, ModelData[model].AICM, ModelData[model].SE_AICM, ModelData[model].dhat, ModelData[model].mean_loglik, ModelData[model].var_loglik);
        }
    }


}