/*
   analyze-pdf.c

   Driver program for the computation of piece-wise analytic fits to
   probability distribution functions, based on sample data.

   Ramses van Zon, December 5, 2009

   Reference:
   [] Ramses van Zon and Jeremy Schofield, "Constructing smooth
      potentials of mean force, radial distribution functions and
      probability densities from sampled data", 
      arxiv.0912.0465 [cond-mat.stat-mech].
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "pa.h"

int analyzepdf( char*filename, int nj, double qm, int km, int mm, int sm )
{
  int            i, j, m, n, nrj;;
  double         p, pj, dp, y, dy, y1, y2, q, *r;
  struct Pafit  *fit[nj];

  r   = readbinfile(filename,&y1,&y2,&m,&n);
  nrj = n-n/nj;

  printf("#Read %d numbers from file %s\n", n, filename);
  printf("#y1 = %6.3f y2 = %6.3f m = %d\n", y1, y2, m);
  printf("#Performing pa analysis\n");
  printf("# NJ = %d QCUT = %f KMAX = %d MPA = %d nosmooth = 0\n", nj, qm, km, mm);

  #pragma omp parallel for 
  for (j = 0; j < nj; ++j) {
    double *rj, qj;
    rj     = (double*)malloc(n*sizeof(double));
    fit[j] = paalloc(km, mm);
    pajack(j, nj, n, r, rj);
    pasort(nrj, rj);
    qj = pafit(nrj, rj, qm, km, mm, sm, fit[j]);
    if (j==0) q = qj;
    free(rj);
  }

  printf("#Writing distribution at %d points with jackknife errors\n", m);
  dy = (y2-y1)/m;
  for (y = y1; y < y2; y += dy) {
    p  = 0;
    dp = 0;
    for (j = 0; j < nj; ++j) {
      pj =  papd(y, fit[j]);
      p  += pj;
      dp += pj*pj;
    }
    p /= nj;
    dp = sqrt((nj-1)/(double)nj*(dp/nj-p*p));
    printf("%f %f %f\n", y, p, dp*1.96);
  }

  printf("#Parameters (q = %f):\n", q);
  pawrite(stdout, fit[0]);

  for (j = 0; j < nj; ++j)
    pafree(fit[j]);
}

int main( int argc, char **argv )
{
  int     nj = 8;      /* number of jackknifes to use */
  double  qm = 0.6;    /* satisfactory value for Q */
  int     km = 21;     /* maximum number of intervals */
  int     mm = 14;     /* maximum number of Fourier modes per interval */

  if (argc <= 1) {
    fprintf(stderr,"%s\n%s\n%s%s\n",
            "analyze-pdf - not enough parameters!",
            "usage:",
            "analyze-pdf <file> ",
            "[<#jacknifes> [<qcut> [<#fourier> [<#intervals> [-smooth]]]]]");
    exit(1);
  } else {
    if (argc > 2) nj = atoi(argv[2]);
    if (argc > 3) qm = atof(argv[3]);
    if (argc > 4) km = atoi(argv[4]);
    if (argc > 5) mm = atoi(argv[5]);      
    printf("# nj=%d qcut=%f km=%d mm=%d nosmooth=%d\n",
           nj, qm, km, mm, argc>6 );
    analyzepdf(argv[1], nj, qm, km, mm, argc<=6);
  }
}