#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <stdarg.h>

#define NMAX_1 101
#define NMAX_2 101

/*****************************************************************************
 *
 * Calculates p(mu,sigma) array, based on a data file.
 *
 * The program asks for the name of as input-information file containing
 * the data file name and the data bounds.
 *
 *
 *****************************************************************************
 *****************************************************************************
 *
 * Harold J. Price 
 * October 10, 2001
 *
 *****************************************************************************
 *****************************************************************************/

int
main(void)
{

    double a,b,low_mu,hi_mu;
    double xdata[NMAX_1],norm_xdata[NMAX_1];
    double sumx,sumxsq;
    double lam1[NMAX_2][NMAX_2],lam2[NMAX_2][NMAX_2];
    double f0_n[NMAX_2][NMAX_2];
    double mu,sigma,musq,sigsq,nusq,mu_peak,nusq_peak,sigma_peak;
    double pofmusig[NMAX_2][NMAX_2];
    double del_mu,del_sig,sum,sum_check;
    double marg_mu[NMAX_2],marg_sig[NMAX_2],sum_mu,sum_sig;

    float dummy_1,dummy_2,dummy_3;

    int nmax;
    int ndata,nx;
    int i,j,num_lam,num_mu,num_sig;
    int n_mu,n_sig;

    char inputfile[256],datafile[256];
    char *p,*delims={" /"},buffer[256];

    FILE *fpin,*fpdata,*fplambda_1,*fplambda_2,*fpf0_n,*fpout;
    FILE *fp_mu,*fp_sig;

/* Start calculation of p(mu,sigma) */

/* Read the input information from the information file */

    printf("please enter the input information file name: ");
    scanf("%s",inputfile);

    printf("%s\n",inputfile);


    if((fpin = fopen(inputfile,"r")) == NULL)
    {
        printf("Can't open %s\n",inputfile);
        exit(1);
    }

/* open and read infofile */

    fgets(datafile,255,fpin);
    strtok(datafile,delims);

    fgets(buffer,255,fpin);
    p = strtok(buffer,delims);
    a = atof(p);

    fgets(buffer,255,fpin);
    p = strtok(buffer,delims);
    b = atof(p);

    fgets(buffer,255,fpin);
    p = strtok(buffer,delims);
    dummy_1 = atof(p);
   
    fgets(buffer,255,fpin);
    p = strtok(buffer,delims);
    low_mu = atof(p);

    fgets(buffer,255,fpin);
    p = strtok(buffer,delims);
    hi_mu = atof(p);

    fgets(buffer,255,fpin);
    p = strtok(buffer,delims);
    nmax = atoi(p);

    fclose(fpin);
    printf("datafile,a,b,low_mu,hi_mu,nmax) %s %e %e %e %e %d\n",
                                     datafile,a,b,low_mu,hi_mu,nmax);

/* open data file */

    if((fpdata = fopen(datafile,"r")) == NULL)
    {
        printf("Can't open %s\n",datafile);
        exit(1);
    }

/* read in data */

    i = 1;
    while((fscanf(fpdata,"%d %e",&ndata,&dummy_1)) != EOF)
    {
        xdata[i] = dummy_1;
        norm_xdata[i] = -1 + 2.0 * (xdata[i] - a)/(b - a);
        
        printf("i,ndata,xdata,norm_xdata %3d %d %e %e\n",i,ndata,xdata[i],
                                                                     norm_xdata[i]);

        i += 1;

    }
    ndata = i - 1;
    fclose(fpdata);

/* calculate the peaks (mu_peak and nusq_peak) of the distribution */

    num_mu = NMAX_2 -1;
    num_sig = NMAX_2 - 1;

    del_mu = 2.0/(num_mu - 1);
    del_sig = 1.0/(num_sig - 1);

    nx = 3;
    sumx = 0.0;
    sumxsq = 0.0;
    for(i = 1;i <= nx;i++)
    { 
        sumx += norm_xdata[i];
        sumxsq += norm_xdata[i] * norm_xdata[i];
    }

    printf("sumx,sumxsq %e %e\n",sumx,sumxsq);

    mu_peak = sumx/nx;
    nusq_peak = sumxsq/nx;
    sigma_peak = sqrt(nusq_peak - (mu_peak * mu_peak));

    printf("mu_peak,sigma_peak : %e %e\n",mu_peak,sigma_peak);

    mu_peak = ((b+a) + mu_peak * (b - a))/2.0;
    sigma_peak = sigma_peak * (b-a)/2.0;

    printf("mu_peak,sigma_peak: %e %e\n",mu_peak,sigma_peak); 

/* read in the lambda table */
    
    if((fplambda_1 = fopen("lambda_1.fil","r")) == NULL)
    {
        printf("Can't open %s\n","lambda_1.fil");
        exit (1);
    }

    if((fplambda_2 = fopen("lambda_2.fil","r")) == NULL)
    {
        printf("Can't open %s\n","lambda_2.fil");
        exit (1);
    }

    if((fpf0_n = fopen("f0_n.fil","r")) == NULL)
    {
        printf("Can't open %s\n","f0_n.fil");
        exit (1);
    }

    for(n_mu = 1;n_mu <= num_mu;n_mu++)
    {
        for(n_sig = 1;n_sig <= num_sig;n_sig++)
        {
            fscanf(fplambda_1,"%e %e %e",&dummy_1,&dummy_2,&dummy_3);
            lam1[n_mu][n_sig] = dummy_3;
        }
    }

    printf(" %d %d\n",n_mu - 1,n_sig - 1);

    for(n_mu = 1;n_mu <= num_mu;n_mu++)
    {
        for(n_sig = 1;n_sig <= num_sig;n_sig++)
        {
            fscanf(fplambda_2,"%e %e %e",&dummy_1,&dummy_2,&dummy_3);
            lam2[n_mu][n_sig] = dummy_3;
        }
    }

    printf(" %d %d\n",n_mu - 1,n_sig - 1);

    for(n_mu = 1;n_mu <= num_mu;n_mu++)
    {
        for(n_sig = 1;n_sig <= num_sig;n_sig++)
        {
            fscanf(fpf0_n,"%e %e %e",&dummy_1,&dummy_2,&dummy_3);
            f0_n[n_mu][n_sig] = dummy_3;
        }
    }

    printf(" %d %d\n",n_mu - 1,n_sig - 1);

    fclose(fplambda_1);
    fclose(fplambda_2);
    fclose(fpf0_n);

/* calculate p(mu,sigma|data) */

    printf("del_mu,del_sig %e %e\n",del_mu,del_sig);

    sum = 0.0;
    for(n_mu = 1;n_mu <= num_mu;n_mu++)
    {
        mu = -1.0 + (n_mu-1) * del_mu;
        for(n_sig = 1;n_sig <= num_sig;n_sig++)
        {
            sigma = (n_sig - 1) * del_sig;
            nusq = mu * mu + sigma * sigma;
            if(f0_n[n_mu][n_sig] == 0)
            {
                pofmusig[n_mu][n_sig] = 0.0;
            }
            else
            {
                pofmusig[n_mu][n_sig] = exp(-(f0_n[n_mu][n_sig] * nx) - 
                                   lam1[n_mu][n_sig] * (sumx - nx * mu) -
                                   lam2[n_mu][n_sig] * (sumxsq - nx * nusq));
            }
            sum += pofmusig[n_mu][n_sig];
        }

    }

/* write the probability distribution to file */

    if((fpout = fopen("pofmusigma.fil","w")) == NULL)
    {
        printf("Can't open %s\n","pofmusig.fil");
    }

/* rescale variables */

    del_mu = (b - a)/(num_mu - 1);
    del_sig = 0.5 *((b - a)/(num_sig - 1));
    sum_check = 0.0;
    for(n_mu = 1;n_mu <= num_mu;n_mu++)
    {
        mu = a + (n_mu - 1) * del_mu;
        for(n_sig = 1;n_sig <= num_sig;n_sig++)
        {
            sigma = (n_sig - 1) * del_sig;
            pofmusig[n_mu][n_sig] /= (sum * del_mu * del_sig);   
            fprintf(fpout,"%e %e %e\n",mu,sigma,pofmusig[n_mu][n_sig]);
            sum_check += pofmusig[n_mu][n_sig];
        }
    }
    sum_check *= (del_mu * del_sig);

    printf("sum_check : %e\n",sum_check);

/* calculate the marginals and write them to files */

    if((fp_mu = fopen("marg_mu.fil","w")) == NULL)
    {
        printf("Can't open %s\n","marg_mu.fil");
    }
    
    if((fp_sig = fopen("marg_sig.fil","w")) == NULL)
    {
        printf("Can't open %s\n","marg_sig.fil");
    }

    sum_check = 0.0;
    for(n_mu = 1;n_mu <= num_mu;n_mu++)
    {
        mu = a + (n_mu - 1) * del_mu;
        sum_mu = 0.0;
        for(n_sig = 1;n_sig <= num_sig;n_sig++)
        {
            sum_mu += pofmusig[n_mu][n_sig];
        }
        sum_check += sum_mu;
        marg_mu[n_mu] = sum_mu * del_sig;
        fprintf(fp_mu,"%e %e\n",mu,marg_mu[n_mu]);
    }
    sum_check *= del_mu * del_sig;
    printf("mu sum_check : %e\n",sum_check);

    sum_check = 0.0;
    for(n_sig = 1;n_sig <= num_sig;n_sig++)
    {
        sigma = (n_sig - 1) * del_sig;
        sum_sig = 0.0;
        for(n_mu = 1;n_mu <= num_mu;n_mu++)
        {
            sum_sig += pofmusig[n_mu][n_sig];
        }
        sum_check += sum_sig;
        marg_sig[n_sig] = sum_sig * del_mu;
        fprintf(fp_sig,"%e %e\n",sigma,marg_sig[n_sig]);
    }
    sum_check *= del_sig * del_mu;
    printf("sigma sum_check : %e\n",sum_check);
    
    fclose(fp_mu);
    fclose(fp_sig);

    printf("end of program\n");


} /* end of main */