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

#define NMAX 10004


/**************************************************************************
 **************************************************************************
 * This program calculates the prob(mu|x) for an array of mu's and an
 * array of sampled values of x.
 * The input information is contained in a file. The file contains the
 * following information:
 * 
 * 	1. The data file name.
 * 	2. The lower limit of possible data 
 * 	3. The upper limit of possible data 
 * 	4  The value of quantile (not used in this program)
 * 	5  The lower limit of prior (must be inside the data limits) 
 * 	6  The upper limit of prior (must be inside the data limits)
 * 	7  The number of mu values in the output (max value 10000)
 * 
 * The probability of mu at the upper and lower possible data limits is 
 * infinite. This is logical, since if mu is at one of these boundaries, 
 * The probability of mu is an impulse.  Thus, for numerical reasons, the 
 * lower and upper limits of the prior must be inside the lower and upper 
 * limits of possible data.  The amount these bounds must be inside is 
 * problem dependent, but a starting point is 0.1%.
 *
 * The output (probability of mu) is written to a file called "output.dat".
 * Also, the values of lambda, as a function of mu, are written to a file
 * called "mu-lambda.arr". The third column of "mu-lambda.arr" also contains
 * the values of f_0 for the root; i.e., how close is f_0 to 0.0,
 * 
 * 
 *************************************************************************
 **************************************************************************
 *
 *  Harold J. Price, March 27, 2001
 *
 ***************************************************************************
 **************************************************************************/

 /*
 * function f_0
 *
 * This function is used to find lambda, given a value for mu; i.e.,
 * lambda is a zero of this function.
 * 
 * */
double f_0(double a, double b, double lam, double mu)
{
    double temp1;
    double temp2;

    if(lam*(a-b) > 700.0)
        {
            return -b-1.0/lam-mu ;
        }
    else if(lam*(b-a) > 700.0)
        {
            return -a-1.0/lam-mu;
        }
    else if(fabs(lam*(b-a)) < 1.0e-3) 
        {
            return -mu-(b+a)/2.0+(b-a)*((b-a)/12.0)*lam;
        }
    else
        {
            temp1 = b*exp(-b*lam)-a*exp(-a*lam);
            temp2 = exp(-a*lam)-exp(-b*lam);
            return temp1/temp2 - 1.0/lam - mu;
        }

} /* end function f_0 */

/*
 * this function calculates the probability of mu, given the upper and lower
 * data limits (a and b) lambda (lam) and a data point (x).
 *
 */
double p_of_mu(double a,double b,double x,double lam)
{
/*
 * This correction is necessary, because of an inadvertant sign omission in my 
 * analytical work; i.e., lam need to be of the opposite sign
 */
    lam = -lam;

    if(x == b) return 0.0;
    if(x == a) return 0.0;
    if(fabs(lam*(b-a)) < 1.0e-4)
    {
        return 1/(b-a)+(a-2.0*a*x+b)*lam/(2.0*(b-a));
    }
    else
    {
        return lam*exp(-lam*(x-a))/(1.0-exp(-lam*(b-a)));
    }
} /* end function p_of_mu */

/*
 * main program
 */

int
main(void)
{
    double a,b,mu,musq,std_dev,exp_value;
    double lam_min,lam_max,low_lam,hi_lam,lam;
    double lambda[NMAX];
    double dummy,testarr[NMAX];
    double mu_array[NMAX],lam_array[NMAX];
    double low_mu,hi_mu,delmu;
    double mu_low,mu_hi,tr_mu_array,tr_delmu;
    double pofmu[NMAX];
    double xdata,tr_xdata;
    double integrate,average,variance;

    float xintermediary;

    int n,nmax,i;
    int test;
    int ndata;

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

    FILE *fpout,*fpout2,*fpin,*fpdata;

/*
 * read the input information from the information file.
 */
 
    printf("Please enter the information file name: ");
    scanf("%s",inputfile);
/*
 * open and read infofile
 */

    if((fpin = fopen(inputfile,"r")) == NULL)
    {
        printf("Can't open %s\n",inputfile);
        exit( 1);
    }
    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 = 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: %s\na: %e\tb: %e\nlow_mu: %e\thi_mu: %e\nOutput resolution: %d\n\n",
                                            datafile,a,b,low_mu,hi_mu,nmax);

/*
 * begin calculation of mu-lambda table
 */
    
    lam_min = -100.0;
    lam_max = 90.0;
   
    delmu = (hi_mu-low_mu)/(nmax-1);
    mu_low = -1.0 + 2.0*(low_mu-a)/(b-a);
    mu_hi = -1.0 + 2.0*(hi_mu-a)/(b-a);
    tr_delmu = (mu_hi-mu_low)/(nmax-1);
   
    for(n=1; n<=nmax; n++)
    {
        mu_array[n] = mu_low + (n-1)*tr_delmu;
/*
        printf("n,mu_array %d %e\n",n,mu_array[n]);
*/

/*
 * Here, we find the value of lambda (root) by a halfing procedure.
 *
 */
    }
    for(n=1; n<=nmax; n++)
    {
        test=0;
        low_lam = lam_min;
        hi_lam = lam_max;
        i=0;
        while(i<=64)
        {
            i += 1;
            lam = (low_lam + hi_lam)/2.0;
/*
            dummy = f_0(a,b,lam,mu_array[n]);
*/
            dummy = f_0(-1.0,1.0,lam,mu_array[n]);
/*
printf("i,lam,dummy: %d %e %e\n",i,lam,dummy);
printf("low_lam,hi_lam: %e %e\n",low_lam,hi_lam);
*/
            if(dummy < 0.0)
                {
                    low_lam = lam;
                }
            else
                {
                    hi_lam = lam;
                }
        } 
        lam_array[n] = (low_lam + hi_lam)/2.0;
        testarr[n] = f_0(-1.0,1.0,lam,mu_array[n]);

    }

/*
 * write mu-lambda data to file
 */
 
    if((fpout = fopen("mu_lambda.arr","w")) == NULL)
        {
            printf("Can't open %s: ","mu-lambda.arr");
            exit(1);
        }
    for(n=1 ;n<=nmax; n++)
    {
        fprintf(fpout,"%e %e %e\n",mu_array[n],lam_array[n],testarr[n]);
    }


    fclose(fpout);

/*
 * calculate probability of x, given mu vs mu.
 */

/*
 * open data file
 */

    if((fpdata=fopen(datafile,"r")) == NULL)
    {
        printf("Can't open %s\n",datafile);
        exit(1);
    }
/*
 * initialize pofmu array
 */
    for(n=1; n<=nmax; n++)
    {
        pofmu[n] = 1.0;
    }
/*
 * read in data and calulculate pofmu array;
 */
    i = 0;
    while((fscanf(fpdata,"%d %e",&ndata,&xintermediary)) != EOF)
    {
        xdata = xintermediary;
        tr_xdata = -1.0 + 2.0*(xdata-a)/(b-a);
/*
        printf("ndata,xdata %d %e\n",ndata,xdata);
*/
        i += 1;
        
        for(n=1; n<=nmax; n++)
        {
            pofmu[n] *= p_of_mu(-1.0,1.0,tr_xdata,lam_array[n]);
        }
/* 
 * normalize pofmu 
 */
        integrate = 0.0; 
        for(n=1; n<=nmax-1; n++)
        {
            integrate += ((pofmu[n] + pofmu[n+1]) * tr_delmu) / 2.0;
        }
/*
        printf("integrate %e\n",integrate);
*/
        for(n=1; n<=nmax; n++)
        {
            pofmu[n] /= integrate;
        }
/*
        printf("integrate 2 %e\n",integrate);
*/
    }

    fclose(fpdata);

/*
 * write pofmu array to file
 */
    if((fpout2=fopen("average.dat","w")) == NULL)
    {
            printf("Can't open %s\n","average.dat");
    }
/*
    fprintf(fpout2,"%e %e\n",a,pofmu[1]);
 */
    for(n=1; n<=nmax; n++)
    {
        tr_mu_array = a + ((b-a)/2.0)*(mu_array[n]+1.0);

        fprintf(fpout2,"%e %e\n",tr_mu_array,pofmu[n]*(2.0/(b-a)));

    }
    for(n=1; n<=nmax; n++)
    {
        pofmu[n] *= 2.0/(b-a);
    }

/*
 * check integral of pofmu and calculate average and variance (about average)
 */
    integrate = 0.0;
    for(n=1; n<=nmax-1; n++)

    {
        integrate += ((pofmu[n]+pofmu[n+1])/2.0)*(delmu);
    }

    printf("The integral of mu is: %e\n",integrate);
/*
 * calculate the average
 */
    printf("delmu %e\n",delmu);
    integrate = 0.0;
    mu = low_mu - delmu/2.0;
    for(n=1; n<=nmax-1; n++)
    {
        mu += delmu;
/*
        printf("n,mu %d %e %e \n",n,mu,mu_array[n]);
*/
        integrate += mu*((pofmu[n]+pofmu[n+1])/2.0)*(delmu);
    }
    exp_value = integrate;
    printf("The expected value of mu is: %e\n",exp_value);
/*
 * calculate the standard deviation
 */
    integrate = 0.0;
    mu = (low_mu) - delmu/2.0;

    for(n=1; n<=nmax-1; n++)
    {
        mu += delmu;
        musq = mu * mu;
/*
        printf("n,mu,musq,integrate %d %e %e %e\n",n,mu,musq,integrate);
*/
        integrate += musq*((pofmu[n]+pofmu[n+1])/2.0)*(delmu);
    }
    std_dev = sqrt(integrate - exp_value*exp_value);
    printf("The standard deviation of mu is: %e\n",std_dev);
    printf("The variance about the origin of mu is: %e\n",integrate);

    fclose(fpout);
        

} /* end main */