/*
   run-analysis.cc

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

   Ramses van Zon, December 14, 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 <iostream>
#include <cstdlib>
#include <cmath>
#include <algorithm>
#include "paxx.h"

using namespace std;
using namespace pa;

// biasing function 
double bias(double x)
{ 
  return 1.0/(x*x); 
}

// random number generator
double rnd(void *s)
{
  double d;
  drand48_r((struct drand48_data*)s,&d);
  return d;
}

// analysis routine
void analyzerdf(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) {
    struct drand48_data s;
    double zj, qj, *rj;
    srand48_r(101+j,&s);
    rj = new double[nrj];
    zj = jbias(j,nj,n,r,bias,rnd,&s,rj);
    sort(rj,rj+nrj);
    qj = fit[j].fit(nrj, rj, qm, km, mm, sm);
    if (j==0) z = zj;
    if (j==0) q = qj;
    delete[] rj;
  }
  delete[] r;
  double nrm = 2*z*y2*y2*y2/M_PI;
  cout << "#Writing " << m << " points of the distribution" << endl;
  double dy=(y2-y1)/m;
  for (double y = y1; y <= y2; y += dy) {
    double rdf = 0;
    double drdf = 0;
    for (int j = 0; j < nj; j++) {
      double rdfj = nrm*fit[j].pd(y);
      rdf += rdfj;
      drdf += rdfj*rdfj;
    }
    rdf /= nj;
    drdf = sqrt((nj-1)/(double)nj*(drdf/nj-rdf*rdf));
    cout << y << " " << rdf << " " << 1.96*drdf << endl;
  }
  cout << "#pafit parameters (q = " << q << ")" << endl;
  fit[0].write(stdout);
}

// driver routine 
int main(int argc, char **argv)
{
  if (argc <= 1) {
    cerr << "run-analysis - not enough parameters!\n"
         << "usage:\n"
         << "run-analysis <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;
    analyzerdf(argv[1], nj, qm, km, mm, argc<=6);
  }
}