/* this is stand alone version of a program to do Bayesian
inference with the metropolis algorithm */

#include <math.h>
#include <iostream>
#include <fstream>
#include <string.h>
#include <stdlib.h>
using namespace std;


/* standard libraries to include */
//#include <math.h>
//#include <iostream.h>
//#include <fstream.h>
//#include <stdlib.h>
#pragma hdrstop

#define FREE_ARG char*

/* macros */
static double dsqrarg;
#define DSQR(a) ((dsqrarg=(a)) == 0.0 ? 0.0 : dsqrarg*dsqrarg)

static float sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)

#define SIGN(a,b) ((b)>= 0.0 ? fabs(a) : -fabs(a))

static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
   (dmaxarg1) : (dmaxarg2))

static double dminarg1,dminarg2;
#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
   (dminarg1) : (dminarg2))

static float maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
   (maxarg1) : (maxarg2))

static int iminarg1,iminarg2;
#define IMIN(a,b) (iminarg1=(a),iminarg2=(b),(iminarg1) < (iminarg2) ?\
   (iminarg1) : (iminarg2))

/* function declarations */
double minFcn(double *x);
void getNegHes(double *initValu,double **obsInfo,double *err,double h,
   int numpar);
void metrop(long setSeed,int lookInt,int noChn,int numIter,int numburn,
   int thinfac,int numpar,double **tryV,double **initValu,double ***simres);
void initlz(double **initValu,int noChn);
void sqrtR(double *stat,int noChn,int numiter,int numpar,double ***tens);
void dfunc(double *p,double *df);
void trnsfrm(double *theta,int n);
void untrnsfrm(double *theta,int n);
void powell(double *p,double **xi,int n,double ftol,int *iter,
   double *fret, double (*func)(double *));
void linmin(double *p,double *xi,int n,double *fret,
   double (*func)(double []));
void mnbrak(double *ax,double *bx,double *cx,double *fa,double *fb,
   double *fc, double (*func)(double));
double brent(double ax,double bx,double cx,double (*f)(double),
   double tol,double *xmin);
double f1dim(double x);
void frprmn(double *p,int n,double ftol,int *iter,double *fret,
   double (*func)(double *),void (*dfunc)(double *,double *));
double dbrent(double ax,double bx,double cx,double (*f)(double),
   double (*df)(double),double tol,double *xmin);
void dlinmin(double *p,double *xi,int n,double *fret,double (*func)(double *),
   void (*dfunc)(double *,double *));
double df1dim(double x);
void amebsa(double **p,double *y,int ndim,double *pb,double *yb,
   double ftol,double (*funk)(double *),int *iter,double temptr);
double amotsa(double **p,double *y,double *psum,int ndim,double *pb,
   double *yb,double (*funk)(double *),int ihi,double *yhi,double fac);
double gasdev(long *idum);
double ran1(long *idum);
double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl, long ndh);
void free_d3tensor(double ***t,long nrl,long nrh,long ncl,long nch,
   long ndl,long ndh);
double **dmatrix(long nrl,long nrh,long ncl,long nch);
void free_dmatrix(double **m, long nrl,long nrh,long ncl,long nch);
void free_matrix(float **m, long nrl,long nrh,long ncl,long nch);
double *dvector(long nl,long nh);
int *ivector(long nl,long nh);
void free_dvector(double *v,long nl,long nh);
void free_ivector(int *v,long nl, long nh);
void choldc(double **a,int n,double p[]);
void invert(double **a,double **res,int n);
void matprod(double **a,int anr,int anc,double **b,int bnr,int bnc,
   double **res);
double ldet(double **a,int n);
void nrerror(char error_text[]);

int nGlob=40,noPar=2;

double *yGlob;

int main()
{
/* this file holds data */
   ifstream filer1("ydata.dat");
   if(!filer1){
      cerr << "Failure to open filer1\n";
      exit(EXIT_FAILURE);
   }
   /* initial values for the maximization */
   ifstream filer2("initVal.dat");
   if(!filer2){
      cerr << "Failure to open filer2\n";
      exit(EXIT_FAILURE);
   }
   /* holds estimated asymptotic cov matrix */
   ifstream filer3("cov.dat");
   if(!filer3){
      cerr << "Failure to open filer3\n";
      exit(EXIT_FAILURE);
   }
   /* parameter value which maximizes posterior */
   ofstream filew1("postMode.dat");
   if(!filew1){
      cerr << "Failure to open filew1\n";
      exit(EXIT_FAILURE);
   }
   /* asymptotic cov matrix at post mode */
   ofstream filew2("asymCov.dat");
   if(!filew2){
      cerr << "Failure to open filew2\n";
      exit(EXIT_FAILURE);
   }
   /* samples from the posterior */
   ofstream filew3("postSmpl.dat");
   if(!filew3){
      cerr << "Failure to open filew3\n";
      exit(EXIT_FAILURE);
   }
   /* numIter is numiter per chain in metrop, numburn is number of burn in iterations,
   noChn is number of chains to use in MCMC, thinfac inicates how much thinning
   of chain output-use 1 to have none, and lookInt indicates how frequently
   convergence diagnostic statistics should be computed */
   int MAXX=1,SIMANL=0,GTVR=0,SMPL=0,numIter=200,numburn=100,noChn=3,lookInt=100,
      thinfac=1,ii=1,i,j,k,l,iter;
   long setSeed=-8623;
   double *initVal,**tryV,***simres;
   yGlob=dvector(1,nGlob);
   initVal=dvector(1,noPar);
   tryV=dmatrix(1,noPar,1,noPar);
   simres=d3tensor(1,noChn,1,(numIter-numburn)/thinfac,1,noPar);
   /* put in limits for censoring */
   /* read in sampled locations */
   filer1.seekg(0);
   for(i=1;i<=nGlob;i++){
      filer1 >> yGlob[i];
   }
   /* posterior mode so far-or initial values as means, cors and sds-read in and transform */
   filer2.seekg(0);
   for(i=1;i<=noPar;i++){
      filer2 >> initVal[i];
   }
   trnsfrm(initVal,noPar);
   cout << "result at initial value=" << minFcn(initVal) << endl;
   if(MAXX==1){
      int maxOpt;
      double retVal,**drctns;
      drctns=dmatrix(1,noPar,1,noPar);
      for(i=1;i<=noPar;i++){
         for(j=1;j<=noPar;j++){
            drctns[i][j]=0.0;
            if(i==j) drctns[i][j]=1.0;
         }
      }
      /* for direction set-no derivative */
      powell(initVal,drctns,noPar,1.0e-7,&iter,&retVal,minFcn);
      /* for conjugate gradient */
//      frprmn(initVal,noPar,1.0e-7,&iter,&retVal,minFcn,dfunc);
      cout << "Number of iterations in minimization=" << iter << endl;
      cout << "minimum value=" << retVal << endl;
      cout << "optimal value" << endl;
      ii=1;
      untrnsfrm(initVal,noPar);
      for(i=1;i<=noPar;i++){
         cout << initVal[i] << " ";
         filew1 << initVal[i] << " ";
      }
      free_dmatrix(drctns,1,noPar,1,noPar);
   }
   if(SIMANL==1){
      double temprtr,yb,*dev,*y,**p;
      y=dvector(1,noPar+1);
      dev=dvector(1,noPar);
      p=dmatrix(1,noPar+1,1,noPar);
      for(i=1;i<=noPar;i++) p[1][i]=initVal[i];
      y[1]=minFcn(initVal);
      cout << y[1] << " ";
      for(i=2;i<=noPar+1;i++){
         for(j=1;j<=noPar;j++){
            p[i][j]=p[1][j];
            if(i==(j+1)) p[i][j]+=1.0;
            dev[j]=p[i][j];
         }
         y[i]=minFcn(dev);
         cout << y[i] << " ";
      }
      yb=1.0e30;
      temprtr=100.0;
      for(i=1;i<=1000;i++){
//         cout << endl << "starting iteration " << i << " ";
         iter=500;
         temprtr*=0.9;
         amebsa(p,y,noPar,initVal,&yb,1.0e-7,minFcn,&iter,temprtr);
         if(iter>0) break;
//         cout << "minVal=" << yb << endl;
//         untrnsfrm(initVal,noPar);
//         for(j=1;j<=noPar;j++) cout << initVal[j] << " ";
//         trnsfrm(initVal,noPar);
      }
      cout << endl << "iter=" << iter << " min value=" << yb << endl;
      untrnsfrm(initVal,noPar);
      for(j=1;j<=noPar;j++){
         cout << initVal[j] << " ";
         filew1 << initVal[j] << " ";
      }
      free_dvector(y,1,noPar+1);
      free_dvector(dev,1,noPar);
      free_dmatrix(p,1,noPar+1,1,noPar);
   }
   if(GTVR==1){
      double hhLoc,maxerr=0.0,*drvErr,**nobsInfo;
      drvErr=dvector(1,noPar*noPar);
      nobsInfo=dmatrix(1,noPar,1,noPar);
      /* here you specify an initial step size */
      hhLoc=1.0e-4;
      getNegHes(initVal,nobsInfo,drvErr,hhLoc,noPar);
      cout << endl << "got nobsInfo-the errors in the hessian estimates are" << endl;
      for(i=1;i<=noPar*(noPar+1)/2;i++){
         cout << drvErr[i] << " ";
         if(drvErr[i]>maxerr) maxerr=drvErr[i];
      }
      cout << endl;
      cout << "largest error=" << maxerr << endl;
      /* nobsInfo is neg of obs Info-so now invert this-note that
      this does an extran chol decomp so a little redundant */
      invert(nobsInfo,tryV,noPar);
//      cout << "scaled variance matrix of MLE's" << endl;
      for(i=1;i<=noPar;i++){
         for(j=1;j<=noPar;j++){
//            cout << tryV[i][j] << " ";
            filew2 << tryV[i][j] << " ";
         }
//         cout << endl;
      }
      free_dvector(drvErr,1,noPar*noPar);
      free_dmatrix(nobsInfo,1,noPar,1,noPar);
   }
   if(SMPL==1){
      double *vec,**initValu;
      initValu=dmatrix(1,noPar,1,noChn);
      vec=dvector(1,noPar);
      /* here is cov matrix based on curvature at mode  */
      filer3.seekg(0);
      for(i=1;i<=noPar;i++){
         for(j=1;j<=noPar;j++){
            filer3 >> tryV[i][j];
//            if(i==j) tryV[i][j]=3.0e-3;
//           else tryV[i][j]=0.0;
         }
      }
      /* initial values moved to the end of the program-read in from a file */
      initlz(initValu,noChn);
      /* for scaling jumping cov matrix */
      for(i=1;i<=noPar;i++){
         for(j=1;j<=noPar;j++){
            tryV[i][j]*=(2.4/sqrt(noPar));
         }
      }
      /* do the metropolis step */
      metrop(setSeed,lookInt,noChn,numIter,numburn,thinfac,noPar,tryV,initValu,
         simres);
      /* write to file */
      for(i=1;i<=noChn;i++){
         for(j=1;j<=(numIter-numburn)/thinfac;j++){
//            if(j%thinfac==0){
               for(k=1;k<=noPar;k++) vec[k]=simres[i][j][k];
               untrnsfrm(vec,noPar);
               for(k=1;k<=noPar;k++){
                  filew3 << vec[k] << " ";
               }
//            }
         }
      }
      free_dvector(vec,1,noPar);
      free_dmatrix(initValu,1,noPar,1,noChn);
   }
   free_dvector(yGlob,1,nGlob);
   free_dvector(initVal,1,noPar);
   free_dmatrix(tryV,1,noPar,1,noPar);
   free_d3tensor(simres,1,noChn,1,(numIter-numburn)/thinfac,1,noPar);
}

double minFcn(double *x)
{
   /* must alter this function when you change the loglik */
   int i;
   double res=0.0;
   for(i=1;i<=nGlob;i++) res+=DSQR(yGlob[i]-x[1]);
   res/=exp(x[2]);
   res+=nGlob*x[2];
   return res;
}

void trnsfrm(double *theta,int n)
{
/* put paramter vector on unbounded scale */
   theta[2]=log(theta[2]);
}

void untrnsfrm(double *theta,int n)
{
/* put paramter vector on interpretable scale */
   theta[2]=exp(theta[2]);
}

void initlz(double **initValu,int noChn)
{
/* starting values for multiple chains-these initial values can be on
   the transformed scale (then set trans=1), e.g. log scale if they are var
   parameters */
   int i,j,k,trans=0;
   double *vec;
   vec=dvector(1,noPar);
   ifstream filer("initValu.dat");
   if(!filer){
      cerr << "Failure to open filer\n";
      exit(EXIT_FAILURE);
   }
   for(i=1;i<=noChn;i++){
      for(j=1;j<=noPar;j++){
         filer >> initValu[j][i];
      }
      if(trans==0){
         for(k=1;k<=noPar;k++) vec[k]=initValu[k][i];
         trnsfrm(vec,noPar);
         for(k=1;k<=noPar;k++) initValu[k][i]=vec[k];
      }
   }
   free_dvector(vec,1,noPar);
}


/* you probably don't want to alter anything below this line */

void metrop(long setSeed,int lookInt,int noChn,int numIter,int numburn,
   int thinfac,int numpar,double **tryV,double **initValu,double ***simres)
{
/* does the metropolis algorithm-lots of parameters on log scale-simres
   holds sim. results in form chain by iteration by parameter, tryV
   is the varaince matrix for the normal jumping kernel, and initValu
   is matrix of initial values for the chain this finds convergence
   diagnostic stat and reports it-stops if converrgence taken place */
   int ii,jj,i,j,k;
   double RSTOP=1.01,sum,*vec,*minVal,*accept,*temp1,*temp2,*Rstat,
      *pp,*mm,*canDev,**dev,**optDev;
   minVal=dvector(1,noChn);
   vec=dvector(1,numpar);
   accept=dvector(1,noChn);
   temp1=dvector(1,noChn);
   temp2=dvector(1,noChn);
   Rstat=dvector(1,numpar);
   pp=dvector(1,numpar);
   mm=dvector(1,numpar);
   canDev=dvector(1,numpar);
   dev=dmatrix(1,numpar,1,noChn);
   optDev=dmatrix(1,numpar,1,noChn);
   choldc(tryV,numpar,pp);
   for(ii=1;ii<=noChn;ii++){
      accept[ii]=0.0;
      for(i=1;i<=numpar;i++){
         dev[i][ii]=initValu[i][ii];
         optDev[i][ii]=initValu[i][ii];
         canDev[i]=initValu[i][ii];
      }
      temp1[ii]=minFcn(canDev);
      minVal[ii]=temp1[ii];
   }
   for(i=1;i<=numIter;i++){
//      cout << "starting iteration number " << i << endl;
      for(ii=1;ii<=noChn;ii++){
//         cout << "starting chain number " << ii << endl;
         for(j=1;j<=numpar;j++) mm[j]=gasdev(&setSeed);
         for(j=1;j<=numpar;j++){
            sum=0.0;
            for(k=1;k<j;k++) sum += tryV[j][k]*mm[k];
            canDev[j] = dev[j][ii]+sum+pp[j]*mm[j];
         }
         temp2[ii]=minFcn(canDev);
         if(temp2[ii] < minVal[ii]){
            for(j=1;j<=numpar;j++) optDev[j][ii]=canDev[j];
            minVal[ii]=temp2[ii];
         }
         if(log(ran1(&setSeed)) < temp1[ii]-temp2[ii]){
            accept[ii]+=1.0;
            temp1[ii]=temp2[ii];
            for(j=1;j<=numpar;j++) dev[j][ii]=canDev[j];
         }
         if(i > numburn){
            if((i-numburn)%thinfac==0){
               for(j=1;j<=numpar;j++){
                  simres[ii][(i-numburn)/thinfac][j]=dev[j][ii];
               }
            }
         }
//         cout << "the minimum neg loglik for chain " << ii << "="
//            << minVal[ii] << " acceptance rate=" << accept[ii]/i << endl;
      }
      if(i%lookInt==0 && i > numburn && noChn > 1){
         sqrtR(Rstat,noChn,(i-numburn)/thinfac,numpar,simres);
         cout << "convergence diagnostic statistics after " << i << " iterations" << endl;
         jj=0;
/* able or disable early return by commenting in/out */
         for(j=1;j<=numpar;j++){
//            if(Rstat[j] < RSTOP) jj+=1;
            cout << Rstat[j] << " ";
         }
         cout << endl;
      }
//      if(jj==numpar) break;
   }
   for(ii=1;ii<=noChn;ii++){
      cout << "for chain " << ii;
      cout << " the acceptance rate=" << accept[ii]/numIter << " ";
      cout << " the minimum neg loglik=" << minVal[ii] << endl;
      cout << "the parameter value for this is" << endl;
      for(j=1;j<=numpar;j++) vec[j]=optDev[j][ii];
      untrnsfrm(vec,numpar);
      for(j=1;j<=numpar;j++) optDev[j][ii]=vec[j];
      for(j=1;j<=numpar;j++) cout << optDev[j][ii] << " ";
      cout << endl;
   }
   free_dvector(minVal,1,noChn);
   free_dvector(vec,1,numpar);
   free_dvector(accept,1,noChn);
   free_dvector(temp1,1,noChn);
   free_dvector(temp2,1,noChn);
   free_dvector(Rstat,1,numpar);
   free_dvector(mm,1,numpar);
   free_dvector(canDev,1,numpar);
   free_dvector(pp,1,numpar);
   free_dmatrix(dev,1,numpar,1,noChn);
   free_dmatrix(optDev,1,numpar,1,noChn);
}

void getNegHes(double *initValu,double **obsInfo,double *err,double h,
   int numpar)
{
/* here we find the negative Hessian of the log lik */
   int NTABd=10,i,j,k,l,m;
   double ERRMAX=1.0e20,CON=1.4,CON2=1.4*1.4,BIG=1.0e30,SAFE=2.0,errt,
      fac,hh,ans,temp1,temp2,temp3,temp4,**a;
   a=dmatrix(1,NTABd,1,NTABd);
   m=1;
   for(i=1;i<=numpar;i++){
      for(j=i;j<=numpar;j++){
        hh=h;
        initValu[i] -= hh;
        initValu[j] -= hh;
        temp4= -minFcn(initValu);
        initValu[i] += 2.0*hh;
        temp2= -minFcn(initValu);
        initValu[j] += 2.0*hh;
        temp1= -minFcn(initValu);
        initValu[i] -= 2.0*hh;
        temp3= -minFcn(initValu);
        initValu[i] += hh;
        initValu[j] -= hh;
        a[1][1]=(temp1-temp2-temp3+temp4)/(4.0*hh*hh);
        err[m]=BIG;
        for(k=2;k<=NTABd;k++){
           hh /= CON;
           initValu[i] -= hh;
           initValu[j] -= hh;
           temp4= -minFcn(initValu);
           initValu[i] += 2.0*hh;
           temp2= -minFcn(initValu);
           initValu[j] += 2.0*hh;
           temp1= -minFcn(initValu);
           initValu[i] -= 2.0*hh;
           temp3= -minFcn(initValu);
           initValu[i] += hh;
           initValu[j] -= hh;
           a[1][k]=(temp1-temp2-temp3+temp4)/(4.0*hh*hh);
           fac=CON2;
           for(l=2;l<=k;l++){
              a[l][k]=(a[l-1][k]*fac-a[l-1][k-1])/(fac-1.0);
              fac=CON2*fac;
              errt=DMAX(fabs(a[l][k]-a[l-1][k]),
                 fabs(a[l][k]-a[l-1][k-1]));
              if(errt <= err[m]){
                 err[m]=errt;
                 ans=a[l][k];
              }
           }
           if(fabs(a[k][k]-a[k-1][k-1]) >= SAFE*(err[m])) break;
        }
      if(err[m] > ERRMAX){
         cout << "error on element " << m << " and error is " <<
            err[m] << endl;
         nrerror("error too big in getHes");
      }
      obsInfo[i][j]= -ans;
      obsInfo[j][i]= obsInfo[i][j];
      m+=1;
      }
   }
   free_dmatrix(a,1,NTABd,1,NTABd);
}

void sqrtR(double *stat,int noChn,int numiter,int numpar,double ***tens)
{
/* this calculates Gelman-Rubin sqrtR for assessing convergence
   for MCMC using overdispersed starting points-set trans=0 to
   calculate on original scale otherwise it calculates it on unbounded scale -
   if you want this yhou must modify below where the comments are */
   double mean(double *,int);
   double var(double *,int);
   int i,j,k,k1,trans=1;
   double *s,*m,*temp,W,inter,BonN;
   temp=dvector(1,numiter);
   s=dvector(1,noChn);
   m=dvector(1,noChn);
   for(i=1;i<=numpar;i++){
      for(j=1;j<=noChn;j++){
         for(k=1;k<=numiter;k++){
            /* can calculate on the original rather than the unbounded scale
            by changing the next few lines */
            if(trans==0 && i>noPar-1) temp[k]=exp(tens[j][k][i]);
            else temp[k]=tens[j][k][i];
         }
         m[j]=mean(temp,numiter);
         s[j]=var(temp,numiter);
      }
      W=mean(s,noChn);
      BonN=var(m,noChn);
      inter=W*(numiter-1.0)/numiter+BonN;
      if(inter<0.0 || W<=0.0) stat[i]=100.0;
      else stat[i]=sqrt((inter)/W);
   }
   free_dvector(temp,1,numiter);
   free_dvector(s,1,noChn);
   free_dvector(m,1,noChn);
}

double mean(double *data, int n)
{
   int j;
   double ave=0.0;
   for(j=1;j<=n;j++) ave += data[j];
   ave /= n;
   return ave;
}

double var(double *data, int n)
{
   double ave,dev,variance;
   double s=0.0;
   int j;
   ave=mean(data,n);
   for(j=1;j<=n;j++){
      dev = data[j]-ave;
      s += dev*dev;
   }
   variance=s/(n-1.0);
   return variance;
}

void dfunc(double *p,double *df)
{
   int i;
   double h=1.0e-4,x;
   x=minFcn(p);
   for(i=1;i<=noPar;i++){
      p[i]+=h;
      df[i]=(minFcn(p)-x)/h;
      p[i]-=h;
   }
}

#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

/* next several routines do direction set and conjugate gradient minimization */

void powell(double *p,double **xi,int n,double ftol,int *iter,
   double *fret, double (*func)(double *))
{
   int maxiter=200,i,ibig,j;
   double del,fp,fptt,t,*pt,*ptt,*xit;
   pt=dvector(1,n);
   ptt=dvector(1,n);
   xit=dvector(1,n);
   *fret=(*func)(p);
   for(j=1;j<=n;j++) pt[j]=p[j];
   for(*iter=1;;++(*iter)){
      fp=(*fret);
      ibig=0;
      del=0.0;
      for(i=1;i<=n;i++){
         for(j=1;j<=n;j++) xit[j]=xi[j][i];
         fptt=(*fret);
         linmin(p,xit,n,fret,func);
         if(fabs(fptt-(*fret)) > del){
            del=fabs(fptt-(*fret));
            ibig=i;
         }
      }
      if(2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))){
         free_dvector(xit,1,n);
         free_dvector(ptt,1,n);
         free_dvector(pt,1,n);
         return;
      }
      if(*iter==maxiter)
         nrerror("powell exceeding max number of iterations");
      for(j=1;j<=n;j++){
         ptt[j]=2.0*p[j]-pt[j];
         xit[j]=p[j]-pt[j];
         pt[j]=p[j];
      }
      fptt=(*func)(ptt);
      if(fptt < fp){
         t=2.0*(fp-2.0*(*fret)+fptt)*DSQR(fp-(*fret)-del)-
            del*DSQR(fp-fptt);
         if(t < 0.0){
            linmin(p,xit,n,fret,func);
            for(j=1;j<=n;j++){
               xi[j][ibig]=xi[j][n];
               xi[j][n]=xit[j];
            }
         }
      }
   }
}

int ncom;
double *pcom,*xicom,(*nrfuncmin)(double []);

void linmin(double *p,double *xi,int n,double *fret,
   double (*func)(double []))
{
   int j;
   double TOL=2.0e-4,xx,xmin,fx,fb,fa,bx,ax;
   ncom=n;
   pcom=dvector(1,n);
   xicom=dvector(1,n);
   nrfuncmin=func;
   for(j=1;j<=n;j++){
      pcom[j]=p[j];
      xicom[j]=xi[j];
   }
   ax=0.0;
   xx=1.0;
   mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
   *fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
   for(j=1;j<=n;j++){
      xi[j] *= xmin;
      p[j] += xi[j];
   }
   free_dvector(xicom,1,n);
   free_dvector(pcom,1,n);
}

extern int ncom;
extern double *pcom,*xicom,(*nrfuncmin)(double []);

double f1dim(double x)
{
   int j;
   double f,*xt;
   xt=dvector(1,ncom);
   for(j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];
   f=(*nrfuncmin)(xt);
   free_dvector(xt,1,ncom);
   return f;
}

void mnbrak(double *ax,double *bx,double *cx,double *fa,double *fb,
   double *fc, double (*func)(double))
{
   double GOLD=1.618034,GLIMIT=100.0,TINY=1.0e-20,ulim,u,r,q,fu,dum;
   *fa=(*func)(*ax);
   *fb=(*func)(*bx);
   if(*fb > *fa){
      SHFT(dum,*ax,*bx,dum)
      SHFT(dum,*fb,*fa,dum)
   }
   *cx=(*bx)+GOLD*(*bx-*ax);
   *fc=(*func)(*cx);
   while(*fb > *fc){
      r=(*bx-*ax)*(*fb-*fc);
      q=(*bx-*cx)*(*fb-*fa);
      u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
         (2.0*SIGN(DMAX(fabs(q-r),TINY),q-r));
      ulim=(*bx)+GLIMIT*(*cx-*bx);
      if((*bx-u)*(u-*cx) > 0.0){
         fu=(*func)(u);
         if(fu < *fc){
            *ax=(*bx);
            *bx=u;
            *fa=(*fb);
            *fb=fu;
            return;
         } else if(fu > *fb){
            *cx=u;
            *fc=fu;
            return;
         }
         u=(*cx)+GOLD*(*cx-*bx);
         fu=(*func)(u);
      } else if((*cx-u)*(u-ulim) > 0.0){
         fu=(*func)(u);
         if(fu < *fc){
            SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
            SHFT(*fb,*fc,fu,(*func)(u))
         }
      } else if((u-ulim)*(ulim-*cx) >= 0.0){
         u=ulim;
         fu=(*func)(u);
      } else {
         u=(*cx)+GOLD*(*cx-*bx);
         fu=(*func)(u);
      }
      SHFT(*ax,*bx,*cx,u)
      SHFT(*fa,*fb,*fc,fu)
   }
}

double brent(double ax,double bx,double cx,double (*f)(double),
   double tol,double *xmin)
{
   int maxiter=100,iter;
   double CGOLD=0.3819660,ZEPS=1.0e-10,a,b,d,etemp,fu,fv,fw,fx,p,q,r,
      tol1,tol2,u,v,w,x,xm;
   double e=0.0;
   a=(ax < cx ? ax : cx);
   b=(ax > cx ? ax : cx);
   x=w=v=bx;
   fw=fv=fx=(*f)(x);
   for(iter=1;iter<=maxiter;iter++){
      xm=0.5*(a+b);
      tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
      if(fabs(x-xm) <= (tol2=0.5*(b-a))){
         *xmin=x;
         return fx;
      }
      if(fabs(e) > tol1){
         r=(x-w)*(fx-fv);
         q=(x-v)*(fx-fw);
         p=(x-v)*q-(x-w)*r;
         q=2.0*(q-r);
         if(q > 0.0) p = -p;
         q=fabs(q);
         etemp=e;
         e=d;
         if(fabs(p)>=fabs(0.5*q*etemp) || p<=q*(a-x) || p>=q*(b-x))
            d=CGOLD*(e=(x >= xm ? a-x : b-x));
         else{
            d=p/q;
            u=x+d;
            if(u-a < tol2 || b-u < tol2)
               d=SIGN(tol1,xm-x);
         }
      } else {
         d=CGOLD*(e=(x >= xm ? a-x : b-x));
      }
      u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
      fu=(*f)(u);
      if(fu <= fx){
         if(u >= x) a=x; else b=x;
         SHFT(v,w,x,u)
         SHFT(fv,fw,fx,fu)
      } else {
         if(u < x) a=u; else b=u;
         if(fu <= fw || w == x){
            v=w;
            w=u;
            fv=fw;
            fw=fu;
         } else if(fu <= fv || v == x || v == w){
            v=u;
            fv=fu;
         }
      }
   }
   nrerror("Too many iterations in brent");
}

#define FREEALL free_dvector(xi,1,n);free_dvector(g,1,n); \
   free_dvector(h,1,n);

void frprmn(double *p,int n,double ftol,int *iter,double *fret,
   double (*func)(double *),void (*dfunc)(double *,double *))
{
   int j,its,itmax=200;
   double gg,gam,fp,dgg,EPS=1.0e-10;
   double *g,*h,*xi;
   g=dvector(1,n);
   h=dvector(1,n);
   xi=dvector(1,n);
   fp=(*func)(p);
   (*dfunc)(p,xi);
   for(j=1;j<=n;j++){
      g[j]= -xi[j];
      xi[j]=h[j]=g[j];
   }
   for(its=1;its<=itmax;its++){
      *iter=its;
      dlinmin(p,xi,n,fret,func,dfunc);
      if(2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)){
         FREEALL
         return;
      }
      fp=(*func)(p);
      (*dfunc)(p,xi);
      dgg=gg=0.0;
      for(j=1;j<=n;j++){
         gg+=g[j]*g[j];
         dgg+=(xi[j]+g[j])*xi[j];
      }
      if(gg==0.0){
         FREEALL
         return;
      }
      gam=dgg/gg;
      for(j=1;j<=n;j++){
         g[j]= -xi[j];
         xi[j]=h[j]=g[j]+gam*h[j];
      }
   }
   nrerror("too many iterations in frprmn");
}

#define MOV3(a,b,c,d,e,f) (a)=(d);(b)=(e);(c)=(f);

double dbrent(double ax,double bx,double cx,double (*f)(double),
   double (*df)(double),double tol,double *xmin)
{
   int iter,ok1,ok2,maxiter=100;
   double a,b,d,d1,d2,du,dv,dw,dx,e=0.0,zeps=1.0e-10;
   double fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,v,w,x,xm;
   a=(ax < cx ? ax : cx);
   b=(ax > cx ? ax : cx);
   x=w=v=bx;
   fw=fv=fx=(*f)(x);
   dw=dv=dx=(*df)(x);
   for(iter=1;iter<=maxiter;iter++){
      xm=0.5*(a+b);
      tol1=tol*fabs(x)+zeps;
      tol2=2.0*tol1;
      if(fabs(x-xm) <= (tol2-0.5*(b-a))){
         *xmin=x;
         return fx;
      }
      if(fabs(e) > tol1){
         d1=2.0*(b-a);
         d2=d1;
         if(dw != dx) d1=(w-x)*dx/(dx-dw);
         if(dv != dx) d2=(v-x)*dx/(dx-dv);
         u1=x+d1;
         u2=x+d2;
         ok1=(a-u1)*(u1-b) > 0.0 && dx*d1 <= 0.0;
         ok2=(a-u2)*(u2-b) > 0.0 && dx*d2 <= 0.0;
         olde=e;
         e=d;
         if(ok1 || ok2){
            if(ok1 && ok2)
               d=(fabs(d1) < fabs(d2) ? d1 : d2);
            else if(ok1)
               d=d1;
            else
               d=d2;
            if(fabs(d) <= fabs(0.5*olde)){
               u=x+d;
               if(u-a < tol2 || b-u < tol2)
                  d=SIGN(tol1,xm-x);
            } else{
               d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
            }
         } else{
            d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
         }
      } else {
         d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
      }
      if(fabs(d) >= tol1){
         u=x+d;
         fu=(*f)(u);
      } else {
         u=x+SIGN(tol1,d);
         fu=(*f)(u);
         if(fu > fx){
            *xmin=x;
            return fx;
         }
      }
      du=(*df)(u);
      if(fu <= fx){
         if(u >= x) a=x; else b=x;
         MOV3(v,fv,dv,w,fw,dw)
         MOV3(w,fw,dw,x,fx,dx)
         MOV3(x,fx,dx,u,fu,du)
      } else {
         if(u<x) a=u; else b=u;
         if(fu <= fw || w == x){
            MOV3(v,fv,dv,w,fw,dw)
            MOV3(w,fw,dw,u,fu,du)
         } else if(fu < fv || v==x || v==w){
            MOV3(v,fv,dv,u,fu,du)
         }
      }
   }
   nrerror("Too many iterations in routine dbrent");
   return 0.0;
}

void (*nrdfun)(double *,double *);

void dlinmin(double *p,double *xi,int n,double *fret,double (*func)(double *),
   void (*dfunc)(double *,double *))
{
   int j;
   double xx,xmin,fx,fb,fa,bx,ax,tol=2.0e-4;
   ncom=n;
   pcom=dvector(1,n);
   xicom=dvector(1,n);
   nrfuncmin=func;
   nrdfun=dfunc;
   for(j=1;j<=n;j++){
      pcom[j]=p[j];
      xicom[j]=xi[j];
   }
   ax=0.0;
   xx=1.0;
   mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
   *fret=dbrent(ax,xx,bx,f1dim,df1dim,tol,&xmin);
   for(j=1;j<=n;j++){
      xi[j] *= xmin;
      p[j] += xi[j];
   }
   free_dvector(xicom,1,n);
   free_dvector(pcom,1,n);
}

extern void (*nrdfun)(double *,double *);

double df1dim(double x)
{
   int j;
   double df1=0.0;
   double *xt,*df;
   xt=dvector(1,ncom);
   df=dvector(1,ncom);
   for(j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];
   (*nrdfun)(xt,df);
   for(j=1;j<=ncom;j++) df1 += df[j]*xicom[j];
   free_dvector(df,1,ncom);
   free_dvector(xt,1,ncom);
   return df1;
}

/* next two functions minimize via simulated annealing */

#define GET_PSUM \
                  for(n=1;n<=ndim;n++) {\
                  for(sum=0.0,m=1;m<=mpts;m++) sum += p[m][n];\
                  psum[n]=sum;}
long idum=-3856;
double tt;

void amebsa(double **p,double *y,int ndim,double *pb,double *yb,
   double ftol,double (*funk)(double *),int *iter,double temptr)
{
   int i,ihi,ilo,j,m,n,mpts=ndim+1;
   double rtol,sum,swap,yhi,ylo,ynhi,ysave,yt,ytry,*psum;
   psum=dvector(1,ndim);
   tt = -temptr;
   GET_PSUM
   for(;;){
      ilo=1;
      ihi=2;
      ynhi=ylo=y[1]+tt*log(ran1(&idum));
      yhi=y[2]+tt*log(ran1(&idum));
      if(ylo > yhi){
         ihi=1;
         ilo=2;
         ynhi=yhi;
         yhi=ylo;
         ylo=ynhi;
      }
      for(i=3;i<=mpts;i++){
         yt=y[i]+tt*log(ran1(&idum));
         if(yt <= ylo){
            ilo=i;
            ylo=yt;
         }
         if(yt > yhi){
            ynhi=yhi;
            ihi=i;
            yhi=yt;
         } else if(yt > ynhi){
            ynhi=yt;
         }
      }
      rtol=2.0*fabs(yhi-ylo)/(fabs(yhi)+fabs(ylo));
      if(rtol < ftol || *iter < 0){
         swap=y[1];
         y[1]=y[ilo];
         y[ilo]=swap;
         for(n=1;n<=ndim;n++){
            swap=p[1][n];
            p[1][n]=p[ilo][n];
            p[ilo][n]=swap;
         }
         break;
      }
      *iter -= 2;
      ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi,-1.0);
      if(ytry <= ylo){
         ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi,2.0);
      } else if(ytry >= ynhi){
         ysave=yhi;
         ytry=amotsa(p,y,psum,ndim,pb,yb,funk,ihi,&yhi,0.5);
         if(ytry >= ysave){
            for(i=1;i<=mpts;i++){
               if(i != ilo){
                  for(j=1;j<=ndim;j++){
                     psum[j]=0.5*(p[i][j]+p[ilo][j]);
                     p[i][j]=psum[j];
                  }
                  y[i]=(*funk)(psum);
               }
            }
            *iter -= ndim;
            GET_PSUM
         }
      } else ++(*iter);
   }
   free_dvector(psum,1,ndim);
}

extern long idum;
extern double tt;

double amotsa(double **p,double *y,double *psum,int ndim,double *pb,
   double *yb,double (*funk)(double *),int ihi,double *yhi,double fac)
{
   int j;
   double fac1,fac2,yflu,ytry,*ptry;
   ptry=dvector(1,ndim);
   fac1=(1.0-fac)/ndim;
   fac2=fac1-fac;
   for(j=1;j<=ndim;j++)
      ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
   ytry=(*funk)(ptry);
   if(ytry <= *yb){
      for(j=1;j<=ndim;j++) pb[j]=ptry[j];
      *yb=ytry;
   }
   yflu=ytry-tt*log(ran1(&idum));
   if(yflu < *yhi){
      y[ihi]=ytry;
      *yhi=yflu;
      for(j=1;j<=ndim;j++){
         psum[j] += ptry[j]-p[ihi][j];
         p[ihi][j]=ptry[j];
      }
   }
   free_dvector(ptry,1,ndim);
   return yflu;
}

/* this function draws a uniform [0,1] deviate */
#define IAran1 16807
#define IMran1 2147483647
#define AMran1 (1.0/IMran1)
#define IQran1 127773
#define IRran1 2836
#define NDIVran1 (1+(IMran1-1)/32)
#define EPSran1 1.2e-7
#define RNMXran1 (1.0-EPSran1)
double ran1(long *idum)
{
   const int NTABran1=32;
   int j;
   long k;
   static long iy=0;
   static long iv[NTABran1];
   double temp;
   if (*idum <= 0 || !iy) {
      if (-(*idum)<1) *idum=1;
      else *idum = -(*idum);
      for(j=NTABran1+7;j>=0;j--) {
         k=(*idum)/IQran1;
         *idum=IAran1*(*idum-k*IQran1)-IRran1*k;
         if (*idum<0) *idum += IMran1;
         if (j<NTABran1) iv[j] = *idum;
      }
      iy=iv[0];
   }
   k=(*idum)/IQran1;
   *idum=IAran1*(*idum-k*IQran1)-IRran1*k;
   if (*idum<0) *idum += IMran1;
   j=iy/NDIVran1;
   iy=iv[j];
   iv[j] = *idum;
   if ((temp=AMran1*iy)>RNMXran1) return RNMXran1;
   else return temp;
}

double gasdev(long *idum)
{
   /* draws a standard normal random variable */
   double ran1(long *idum);
   static int iset=0;
   static double gset;
   double fac,rsq,v1,v2;
   if(iset==0){
      do{
         v1=2.0*ran1(idum)-1.0;
         v2=2.0*ran1(idum)-1.0;
         rsq=v1*v1+v2*v2;
      } while (rsq>=1.0 || rsq==0.0);
      fac=sqrt(-2.0*log(rsq)/rsq);
      gset=v1*fac;
      iset=1;
      return v2*fac;
   } else{
      iset=0;
      return gset;
      }
}

/* routines for linear algebra */
void choldc(double **a, int n, double p[])
{
/* this finds the cholesky decomp. returns result in lower half
   of a and the diagonal in p */
   void nrerror(char error_text[]);
   int i,j,k;
   double sum;
   for(i=1;i<=n;i++){
      for(j=i;j<=n;j++){
         for(sum=a[i][j],k=i-1;k>=1;k--) sum -= a[i][k]*a[j][k];
         if (i==j){
             if (sum<=0.0){
                 nrerror("choldc failed");
              }
             p[i]=sqrt(sum);
         } else a[j][i]=sum/p[i];
      }
   }
}

void invert(double **a,double **res, int n)
/* this function inverts a covariance (sym pd) matrix using choldc
   be careful! a gets overwritten so don't try to use it again */
{
   double *p;
   p=dvector(1,n);
   for(int i=1;i<=n;i++){
      p[i]=1.0;
   }
   choldc(a,n,p);
   /* this loop calculates L^-1 */
   for(int i=1;i<=n;i++){
      a[i][i]=1.0/p[i];
      for(int j=i+1;j<=n;j++){
         double sum1=0.0;
         for(int k=i;k<j;k++) sum1 -= a[j][k]*a[k][i];
         a[j][i]=sum1/p[j];
      }
   }
   /* now find t(L^-1)%*%L^-1 */
   for(int i=1;i<=n;i++){
      for(int j=i;j<=n;j++){
         double sum2=0.0;
         for(int k=j;k<=n;k++) sum2 += a[k][i]*a[k][j];
         res[i][j]=sum2;
      }
   }
   for(int i=1;i<=n;i++){
      for(int j=1;j<=i;j++){
         res[i][j]=res[j][i];
      }
   }
   free_dvector(p,1,n);
}

void matprod(double **a,int anr,int anc,double **b,int bnr,int bnc,
   double **res)
{
/* multiplies two matrices and returns result in res */
   if(anc != bnr) nrerror("A and B can't be multiplied!");
   for(int i=1;i<=anr;i++){
      for(int j=1;j<=bnc;j++){
         double sum=0.0;
         for(int k=1;k<=anc;k++) sum += a[i][k]*b[k][j];
         res[i][j]=sum;
      }
   }
}

double ldet(double **a,int n)
{
   /* returns the log of the determinant of a pd matrix*/
   int i,j;
   double d=0.0,*p;
   double **b;
   b=dmatrix(1,n,1,n);
   p=dvector(1,n);
   for(i=1;i<=n;i++){
      for(j=1;j<=n;j++){
         b[i][j]=a[i][j];
      }
   }
   choldc(b,n,p);
   for(j=1;j<=n;j++) d += log(p[j]);
   free_dvector(p,1,n);
   free_dmatrix(b,1,n,1,n);
   return 2.0*d;
}

/* routines for dealing with memory and creating vectors and matrices */

double *dvector(long nl,long nh)
/* allocate a vector with index range v[nl...nh] */
{
   void nrerror(char error_text[]);
   double *v;
   v=(double *)malloc((size_t) ((nh-nl+2)*sizeof(double)));
   if (!v) nrerror("allocation failure in dvector()");
   return v-nl+1;
}

void free_dvector(double *v, long nl, long nh)
{
   free((FREE_ARG) (v+nl-1));
}

int *ivector(long nl,long nh)
/* allocate a vector with index range v[nl...nh] */
{
   void nrerror(char error_text[]);
   int *v;
   v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
   if (!v) nrerror("allocation failure in dvector()");
   return v-nl+1;
}

void free_ivector(int *v, long nl, long nh)
{
   free((FREE_ARG) (v+nl-1));
}

double **dmatrix(long nrl,long nrh,long ncl,long nch)
/* allocate a double matrix with subscript range m[nrl...nrh][ncl...nch] */
{
   long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
   double **m;
   void nrerror(char error_text[]);
   /* allocate pointers to rows */
   m=(double **) malloc((size_t)((nrow+1)*sizeof(double*)));
   if (!m) nrerror("allocation failure 1 in dmatrix()");
   m += 1;
   m -= nrl;
   /* allocate rows and set pointers to them */
   m[nrl]=(double *) malloc((size_t)((nrow*ncol+1)*sizeof(double)));
   if (!m[nrl]) nrerror("allocation failure 2 in dmatrix()");
   m[nrl] += 1;
   m[nrl] -= ncl;
   for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
   /* return pointer to array of pointers to rows */
   return m;
}

void free_dmatrix(double **m, long nrl,long nrh,long ncl,long nch)
/* free a double matrix allocated by dmatrix() */
{
   free((FREE_ARG) (m[nrl]+ncl-1));
   free((FREE_ARG) (m+nrl-1));
}

double ***d3tensor(long nrl,long nrh,long ncl,long nch,long ndl, long ndh)
/* allocate a double 3-tensor */
{
   long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1,ndep=ndh-ndl+1;
   double ***t;
   /* allocate pointers ot rows */
   t=(double ***) malloc((size_t)((nrow+1)*sizeof(double**)));
   if(!t) nrerror("allocation failure 1 in d3tensor()");
   t += 1;
   t -= nrl;
   /* allocate pointers to rows and set pointers to them */
   t[nrl]=(double **) malloc((size_t)((nrow*ncol+1)*sizeof(double*)));
   if(!t[nrl]) nrerror("allocation failure 2 in d3tensor()");
   t[nrl] += 1;
   t[nrl] -= ncl;
   /* allocate rows and set pointers to them */
   t[nrl][ncl]=(double *) malloc((size_t)((nrow*ncol*ndep+1)*sizeof(double)));
   if(!t[nrl][ncl]) nrerror("allocation failure 3 in d3tensor()");
   t[nrl][ncl] += 1;
   t[nrl][ncl] -= ndl;
   for(j=ncl+1;j<=nch;j++) t[nrl][j]=t[nrl][j-1]+ndep;
   for(i=nrl+1;i<=nrh;i++){
      t[i]=t[i-1]+ncol;
      t[i][ncl]=t[i-1][ncl]+ncol*ndep;
      for(j=ncl+1;j<=nch;j++) t[i][j]=t[i][j-1]+ndep;
   }
   /* return pointer to array of pointers to rows */
   return t;
}

void free_d3tensor(double ***t,long nrl,long nrh,long ncl,long nch,
   long ndl,long ndh)
/* free a double tensor declared by d3tensor */
{
   free((FREE_ARG) (t[nrl][ncl]+ndl-1));
   free((FREE_ARG) (t[nrl]+ncl-1));
   free((FREE_ARG) (t+nrl-1));
}

void nrerror(char error_text[])
/* Numerical Recipes standard error handler */
{
   cerr << "Run-time error...\n";
   cerr << error_text;
   cerr << "...now exiting to system...\n";
   exit(EXIT_FAILURE);
//   fprintf(stderr,"Run-time error...\n");
//   fprintf(stderr,"%s\n",error_text);
//   fprintf(stderr,"...now exiting to system...\n");
//   exit(1);
}