/*
   analyze-pdf.cc

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

   Ramses van Zon, December 14, 2009

   See: 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 <iostream>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include "paxx.h"

using namespace std;
using namespace pa;

/* analysis routine */
void analyzepdf(char*filename, int nj, double qm, int km, int mm, bool sm)
{
  int m, n;
  double y1, y2;
  double *r = newreadbinfile(filename, y1, y2, m, n);
  cout << "#Read " << n << " numbers from file " << filename << "\n"
       << "#y1 = " << y1 << " y2 = " << y2 << " m = " << m << "\n"
       << "#Performing pafit analysis" << endl;
  double z, q;
  int nrj = n-n/nj;
  Fit fit[nj];
  #pragma omp parallel for
  for (int j = 0; j < nj; ++j) {
    double *rj = new double[nrj];
    jack(j,nj,n,r,rj);  
    sort(rj,rj+nrj);
    double qj = fit[j].fit(nrj, rj, qm, km, mm, sm);
    if (j==0) q = qj;
    delete[] rj;
  }
  delete[] r;
  cout << "#Writing " << m
       << " points of the distribution with jackknife errors" << endl;
  double dy = (y2-y1)/m;
  for (double y = y1; y <= y2; y += dy) {
    double p = 0;
    double dp = 0;
    for (int j = 0; j < nj; j++) {
      double pj = fit[j].pd(y);
      p += pj;
      dp += pj*pj;
    }
    p /= nj;
    dp = sqrt((nj-1)/(double)nj*(dp/nj-p*p));
    cout << y << " " << p << " " << 1.96*dp << endl;
  }
  cout << "#pafit parameters (q = " << q << ")" << endl;
  fit[0].write(stdout);
}

/* driver routine */
int main(int argc, char **argv)
{
  if (argc <= 1) {
    cerr << "analyze-pdf - not enough parameters!\n"
         << "usage:\n"
         << "analyze-pdf <file> "
         << "[<#jacknifes> [<qcut> [<#fourier> [<#intervals> [-smooth]]]]]";
    exit(1);
  } else {
    int     nj = 8;      
    double  qm = 0.6; 
    int     km = 21;     
    int     mm = 14;     
    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]);
    cout << "# nj=" << nj << " qcut=" << qm << " km=" << km 
         << " mm=" << mm << " nosmooth=" << (int(argc>6)) << endl;
    analyzepdf(argv[1], nj, qm, km, mm, argc<=6);
  }
}