/*
   run-analysis.c

   Driver program for the computation of piece-wise analytic fits to
   radial 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"

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

double rnd( void *s )
{
  double d;
  drand48_r(s,&d);
  return d;
}

void analyzerdf( char* filename, int nj, double qm, int km, int mm, int sm )
{
  int           j, m, n, nrj;
  double        y, y1, y2, z, q, rdf, rdfj, drdf, nrm, *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 pafit analysis\n");

  #pragma omp parallel for
  for (j = 0; j < nj; ++j) {
    struct drand48_data s;
    double zj, qj, *rj;
    srand48_r(101+j,&s);
    rj = (double*)malloc(nrj*sizeof(double));
    zj = pajbias(j,nj,n,r,bias,&s,rj);
    pasort(nrj,rj);
    fit[j] = paalloc(km, mm);
    qj = pafit(nrj, rj, qm, km, mm, sm, fit[j]);
    if (j==0) z = zj;
    if (j==0) q = qj;
    free(rj);
  }
  nrm = 2*z*y2*y2*y2/M_PI;

  printf("#Writing %d points of the distribution\n", m);
  for (y = y1; y <= y2; y += (y2-y1)/m) {
    rdf = drdf = 0;
    for (j = 0;j < nj; j++) {
      rdfj =  nrm*papd(y, fit[j]);
      rdf  += rdfj;
      drdf += rdfj*rdfj;
    }
    rdf /= nj;
    drdf = sqrt((nj-1)/(double)nj*(drdf/nj-rdf*rdf));
    printf("%f %f %f\n", y, rdf, 1.96*drdf);
  }

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

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

/* driver routine */
int main( int argc, char **argv )
{
  int     nj = 8;      
  double  qm = 0.6; 
  int     km = 21;     
  int     mm = 14;     

  if (argc <= 1) {
    fprintf(stderr,"%s\n%s\n%s%s\n",
            "run-analysis - not enough parameters!",
            "usage:",
            "run-analysis <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 );
    analyzerdf(argv[1], nj, qm, km, mm, argc<=6);
  }
}