/*===========================================================================
Name: conf_intervals_betafunction.c
===========================================================================*/

/*###########################################################################
This code calculates the 1-2-3 sigma (Gaussian) confidence intervals for
the posterior distribution of the IDC parameter following the prescriptions
of the reference paper reported below.
The same posterior distribution is calculated for a uniform sampling of the
[0,1] interval with 1e-3 resolution and printed to output.

Requires GSL libraries, see http://www.gnu.org/software/gsl/ 

Reference paper:  Romano et al. 2014, A&A
Title: "Constraining duty cycles through a Bayesian technique"

Date   : Oct 2014
Author : Cristiano Guidorzi (UNIFE)
###########################################################################*/
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_cdf.h>
#define  PSI  1001 // output posterior size
#define  BUF  200

main(int argc, char *argv[])
{
  FILE     *out;
  char     buf[BUF];
  double   IDCe, m, cl[3],di,p1[3],p2[3], b1[3], b2[3];
  double   x[PSI], y[PSI], a, b;
  int      N, i;

  if(argc!=4)
    {
      puts("-------------------------------------------------------------------");
      puts("This code calculates the 1-2-3 sigma (Gaussian) confidence intervals");
      puts("for IDC from its posterior distribution following the prescriptions");
      puts("of the reference paper reported below.");
      puts("The same posterior distribution is calculated for a uniform sampling of the");
      puts("[0,1] interval with 1e-3 resolution and printed to output.\n");
      puts("Reference paper:  Romano et al. 2014, A&A");
      puts("Title: \"Constraining duty cycles through a Bayesian technique\"\n");
      puts("Date   : Oct 2014");
      puts("Author : Cristiano Guidorzi (UNIFE)\n");

      puts("\nUsage> 3 mandatory args");
      puts("\t1. N");
      puts("\t2. IDCe");
      puts("\t3. output file\n");
      fprintf(stdout,"Example:\n%s 15 0.11 posterior.out\n\n",argv[0]);
      puts("-------------------------------------------------------------------\n");
      exit(1);
    }

  sscanf(argv[1], "%d", &N);
  sscanf(argv[2], "%lf", &IDCe);

  if((out=fopen(argv[3],"w"))==NULL)
    {
      printf("Cannot open file %s for writing\n",argv[3]);
      exit(1);
    }
  
  m = N*IDCe;

  fprintf(stdout, "# N    : %d\n",N);
  fprintf(stdout, "# IDCe : %lf\n",IDCe);
  fprintf(stdout, "# m    : %lf\n",m);
  fprintf(out, "# N    : %d\n",N);
  fprintf(out, "# IDCe : %lf\n",IDCe);
  fprintf(out, "# m    : %lf\n",m);

  // -----Uniform prior is default. Ask user -----
  fprintf(stdout, "Prior parameter a (default=1 for uniform): \n");
  fgets(buf,BUF,stdin);
  if(buf[0]=='\n') a = 1;
  else sscanf(buf, "%lg", &a);
  fprintf(stdout, "Prior parameter b (default=1 for uniform): \n");
  fgets(buf,BUF,stdin);
  if(buf[0]=='\n') b = 1;
  else sscanf(buf, "%lg", &b);

  fprintf(stdout, "# Prior parameter a  (default=1 for uniform prior) : %lg\n",a);
  fprintf(stdout, "# Prior parameter b  (default=1 for uniform prior) : %lg\n",b);
  fprintf(out, "# Prior parameter a  (default=1 for uniform prior) : %lg\n",a);
  fprintf(out, "# Prior parameter b  (default=1 for uniform prior) : %lg\n",b);
  //----------------------------------------------


  for(i=0; i<3; ++i) {
    di = (double) i+1;
    cl[i] = gsl_cdf_ugaussian_P(di) - gsl_cdf_ugaussian_P(-di);
    p1[i] = (1.-cl[i])/2.0;
    p2[i] = (1.+cl[i])/2.0;
    b1[i] = gsl_cdf_beta_Pinv(p1[i], m+a, N-m+b);
    b2[i] = gsl_cdf_beta_Pinv(p2[i], m+a, N-m+b);

    fprintf(stdout, "# Bounds at %d sig(Gauss): %lg\t%lg\n",i+1,b1[i],b2[i]);
    fprintf(out, "# Bounds at %d sig(Gauss): %lg\t%lg\n",i+1,b1[i],b2[i]);
  }

  for(i=0; i<PSI; ++i) {
    x[i] = ((double) i)/PSI;
    y[i] = gsl_ran_beta_pdf(x[i], m+a, N-m+b);

    //fprintf(stdout, "%5.3f\t%lg\n",x[i],y[i]);
    fprintf(out, "%5.3f\t%lg\n",x[i],y[i]);
  }

  fclose(out);

}