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

/**********************************************************
***********************************************************
This program calcuates a nonparametric distribution of the
quantile based on the data presented.

Input information is given in an info file in the
following format:

        1. Data file name.
        2. Lower limit on data distribution.
        3. Upper limit on data distribution.
        4. Quantile, expressed as a decimal fraction.
        5. Lower limit on constant prior.
        6. Upper limit on constant prior.
        7. Number of data points in the output array.

Harold J. Price, February 24, 2001.

***********************************************************
**********************************************************/

/**********************************************************
    prob of x1 given nu 
**********************************************************/
double prx1gv(int i,double a,double b,double data[],double v,
    double quantile){
    if(v < a || v > b) return 0.0;
    if(data[i] >= v){
        return ((1.-quantile))/(b-v);
    }
    else{
        return (quantile/(v-a));
    }
} /* end of prx1gv */
    
    
        
int
main(void)
{
    int i,j,ndata,nout;
    
    double a,b,quantile,aa,bb,delv,v,temp;
    double data[1000],dist[10000];
    
    char buffer[256],infofile[256],datafile[256];
    char *delims={" /"};
    char *p;
    
    FILE *fp_in,*fp_out,*fp_data;

    /**********************/
    /* get the input data */
    /**********************/

    printf("Please enter the information file name: ");
    scanf("%s",infofile);

    if((fp_in = fopen(infofile,"r")) == NULL){
        printf("Can't open %s\n",infofile);
        exit(1);
    }
    
    /* Open and read infofile */
    fgets(datafile,255,fp_in);
    strtok(datafile,delims);
    if((fp_data = fopen(datafile,"r")) == NULL){
        printf("Can't open %s\n",datafile);
        exit(1);
    }
    fgets(buffer,255,fp_in);
    p = strtok(buffer,delims);
    a = atof(p);
    fgets(buffer,255,fp_in);
    p=strtok(buffer,delims);
    b = atof(p);
    fgets(buffer,255,fp_in);
    p = strtok(buffer,delims);
    quantile = atof(p);
    fgets(buffer,255,fp_in);
    p = strtok(buffer,delims);
    aa = atof(p);
    fgets(buffer,255,fp_in);
    p = strtok(buffer,delims);
    bb = atof(p);
    fgets(buffer,255,fp_in);
    p = strtok(buffer,delims);
    nout = atoi(p);
    fclose(fp_in);

    printf("\n\n\ndata: %s\na: %e\tb: %e\nquantile: %e\noutput resolution: %d\n",datafile,a,b,quantile,nout);


    /*****************************/
    /* Calcuate the distribution */
    /*****************************/
     
    /* first read data input file */
    ndata = 1;
    while(fgets(buffer,255,fp_data) != NULL){
        p=strtok(buffer,delims); /* data point is 2nd */
        p=strtok(NULL,delims);
        data[ndata] = atof(p);
        ndata += 1;        
    }
    ndata -= 1;

    printf("\n\nThe Number of data points is:%d\n",ndata);

    /* Set the output array to the prior for v; i.e., */
    /* unity for aa >= v <= bb; 0.0, otherwise         */

    /* calculate delv */
    delv = (b-a)/(nout- 1);

    for(i = 1;i <= nout;i++){
        v = a + (i-1)*delv;
        if (v < aa) dist[i] = 0.0;
        if (v >= aa && v <= bb) dist[i] = 1./(bb-aa);
        if (v > bb) dist[i] = 0.0;
    }

    /* compute the distribution */

    
    for(i = 1;i <= ndata;i++){
        for(j = 1;j <= nout;j++){
            v = a + (j-1)*delv;
            dist[j] *= prx1gv(i,a,b,data,v,quantile);
        }
        
        /* Normalize dist */
        temp = 0.0;
        for(j = 1;j <= (nout-1);j++){
            temp += 0.5*(dist[j]+dist[j+1]);
        }
        for(j = 1;j <= nout;j++){
            dist[j] /= (temp*delv);
        }
    }
    fclose(fp_data);
    
    /*********************************/
    /* Write dist to the output file */
    /*********************************/   
            
    if((fp_out = fopen("output.dat","w")) == NULL){
        printf("Can't open %s\n","output.dat");
    }
    for(j = 1;j <= nout;j++){
        v = a + (j-1)*delv;
        fprintf(fp_out,"%e  %e\n",v,dist[j]);
    }
    fclose(fp_out);

}/* End of main */