/* FITTING ROUTINE -- WITH MODIFIED DATA TARGETS AND 4 INPUT PARAMETERS

This is a shell program which calls 'amebsa' from Numerical Recipes, a downhill
simplex method with simulated annealing. This shell program is tailored to the
TB individual-based model, where 'funx' and 'gof' were written by Adrienne Keen
and somewhat specific to the TB model. The rest of the code was adapted from
code for a testing program for amebsa fitting, 'atest' written by Clarence
Lehman in 2004, which in turn was based on the Numerical Recipes code for
amebsa.

Note: ptex is a batch file to convert these into printed page), calls tex and
then calls things to make dvi out of it.

Update: code modified to include parallel processing of separate calls to TB
program. This is done by integrating the MPI modules in mpi.c, written by
Clarence, based on his ideas and commands given to us by those at MSI.

Update2: This version reads in modified data targets and accepts only four
variable parameters, three disease risks and 'df'.
Example new program call:
     1    2           3             4
tb32i df=3 d1uk20=0.15 d2uk20=0.0004 d3uk20=0.09
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
#include "nrutil.h"
#include "mpi.h"
#include "fileio.h"

#define X              //Temporary marker on lines recently added or changed.
#define ITER   25      //Number of iterations at a given temperature.
#define W      .8      //Displacement in the starting simplex.
#define NPROC 500      //Maximum number of processors.
#define RED    .7      //Temperature reduction factor.
#define THRESH .01     //Convergence threshold for parameter convergence.

//Define statements specific to TB model
#define AC  4                   //Number of age classes for notification rates.
#define SEX 2                   //Number of sexes.
#define ROB 3                   //Number of regions of birth in the model.
#define YR  11                  //Number of years of model output for fitting
#define IND (AC*SEX*ROB*YR)     //to notification data (1999-2009).
#define NOTIF 0                 //Target data version, 0=rates, 1=numbers.
#define AR  (r==2? 1: a)        //Index forcing age class 1 for SSA.          X

//Function prototypes
double converge(int,int,double[]);     //Test for convergence of parameters.
float funx(float *);                   //Call TB model and evaluate fit of
                                       //output to calibration targets.
float gof(float[ROB][YR][SEX][AC]);    //Evaluate discrepancy between model
                                       //output and data (targets).
int mpOpen();                          //Open parallel processing.
int mpClose();                         //Close parallel processing.
int mpCommon(double[][],double[],int); //Gather output from multiple processors.
void simo(float[ROB][YR][SEX][AC],double[NPROC][ROB][YR][SEX][AC],double[ROB][YR][SEX][AC],double[ROB][YR][SEX][AC]);//Write simulation output to tex file.
int Error(double);
int Error1(double,char*,double);
int Error2(double,char*,double,char*,double);
int Error3(double,char*,double,char*,double,char*,double);

#define forROB for(r=0; r<ROB;   r++)  //Standard loops.                     X
#define forSEX for(s=0; s<SEX;   s++)                                        X
#define forAC  for(a=0; a<AC;    a++)                                        X
#define forYR  for(y=0; y<YR;    y++)                                        X
#define forN   for(n=0; n<nproc; n++)                                        X
#define forQ   for(q=0; q<QN;    q++)                                        X
                                                                             X
#define qS0  0    //Number of observations (power 0, becoming n).            X
#define qS1  1    //Sum of the observations (power 1, becoming the mu).      X
#define qS2  2    //Sum of squared observations (power 2, becoming sigma).   X
#define QN   3    //Size of statistics vector.                               X

//Global variables
double tar[ROB][YR][SEX][AC];   //Array which holds all calibration targets.
float tar2[ROB][YR][SEX][AC];   //Float array which holds all targets.
float sim4[ROB][SEX][AC][QN];   //Statistics for fitting, n, mu, sigma.      X


long idum;          //Random number working seed.

int ndim     = 4;    //Number of parameters to vary.
float temptr = 2.0;  //Current temperature.
int iter     = ITER; //Number of function evaluations at current temperature.

float **p;           //Vertices of the search simplex, a pointer to a list of
                     //pointers or array of pointers.
float  *y;           //Values of objective function at each vertex.

float ftol = .001;  //Fractional convergence tolerance for early return.
float *pb, yb;      //Best point and function value (returned by 'amebsa').
FILE  *pt;          //Pointer to tex file for printing progress of fitting.
extern int cproc;   //Current processor number.
extern int nproc;   //Total number of processors.
FILE *fout;         //File for checking simulated annealing progress.
FILE *plot;         //Simulation output in tex.

struct IO fmt2[] =  //Format statements for input/output.
{ /*00*/ { (double*)tar,
{-'r',ROB,-'y',YR,-'s',SEX,-'a',AC},{-'r',-'y',-'A',-'s'} } };

static char *macro[] = { "pictex", "graphs", 0 };

/*------------------------------------------------------------------------------
MAIN PROGRAM
*/

main(int argc, char *argv[])
{
  int i;
  ErrorInit();                            //Trap system failures.
  mpOpen();                               //Open parallel processing tools.
  if(nproc>NPROC) exit(3);                //Check for appropriate number
                                          //of processors.
  char fitout[100];                       //Create string for output file and
  sprintf(fitout,"fitR56-%02d.txt",cproc); //open file.
  fout=fopen(fitout,"w");
  if(fout==NULL)
  { printf("Error opening output file!"); exit(3); }
  if(cproc==0)
  { plot = fopen("plotR56.tex","w"); //Open file for simulation output
    if(plot==NULL)                        //and print headers.
    { printf("Error opening output file!"); exit(3); }
      for(i=0; macro[i]; i++)           //Opening all graphs.
        fprintf(plot, "\\input %s.tex\n", macro[i]);
      fprintf(plot, "\n"); }

fprintf(fout,"Starting fit5, cproc=%d and nproc=%d\n",cproc,nproc); fflush(fout);
fprintf(fout,"temptr = %f iter = %d and temptr red factor = %f\n",temptr,iter,RED); fflush(fout);
  if(NOTIF)                                  //Read in data targets, where
  { if(ROB==3)                               //appropriate file depends on
      FileIO("targc1n.txt",fmt2[0],"r|");    //model version (SSA/nonSSA) and
    else                                     //whether fitting to case numbers
      FileIO("targc0n.txt",fmt2[0],"r|"); }  //or notification rates (NOTIF=1
  else                                       //or NOTIF=0).
  { if(ROB==3)
      FileIO("targc1.txt",fmt2[0],"r|");
    else
      FileIO("targc0.txt",fmt2[0],"r|"); }

  char tout[100];                            //For testing: Write out targets.
  sprintf(tout,"tout%02d.txt",cproc);
  FileIO(tout,fmt2[0],"w|");

  if(argc>1) temptr = atof(argv[1]);         //Retrieve starting temperature.

  p = matrix((long)1, (long)ndim+1,          //Allocate and initialize
             (long)1, (long)ndim);           //the starting simplex.

p[1][1]=1.841155;  p[1][2]=0.163031;  p[1][3]=0.001048;  p[1][4]=0.0755;
p[2][1]=0.560994;  p[2][2]=0.106223;  p[2][3]=0.00038;  p[2][4]=0.169378;
p[3][1]=1.297929;  p[3][2]=0.07736;  p[3][3]=0.000933;  p[3][4]=0.089043;
p[4][1]=2.534156;  p[4][2]=0.099439;  p[4][3]=0.001034;  p[4][4]=0.027495;
p[5][1]=1.822337;  p[5][2]=0.094479;  p[5][3]=0.000948;  p[5][4]=0.086628;

  pb = vector((long)1, (long)ndim+1);

//- fprintf(fout,"About to get funx value at initial simplex points\n"); fflush(fout);

  y = vector((long)1, (long)ndim+1);         //Allocate and initialize the
  for(i=1; i<=ndim+1; i++)                   //objective function values at
    y[i] = funx(p[i]);                       //each vertex of the initial simplex.

//- fprintf(fout,"About to start sim annealing loop \n"); fflush(fout);

  for(yb=10000000; 1; temptr*=RED)           //Run through the annealing.
  { iter = ITER;
    amebsa(p, y, ndim, pb, &yb, ftol, funx, &iter, temptr);
    fprintf(fout,"! %% temptr=%f idum=%ld iter=%d yb=%f pb=%f,%f,%f,%f\n\n",
      temptr, idum, iter, yb, pb[1],pb[2],pb[3],pb[4]); fflush(fout);
    if(iter>0) break; }

  fclose(fout);                              //Close annealing output file.
  if(cproc==0)                               //If this is first processor,
    fclose(plot);                            //close simulation output file.
  mpClose();                                 //Close parallel processing.
  exit(0);                                   //Return to operating system.
}



/*-----------------------------------------------------------------------------
GOODNESS OF FIT

This function evaluates and returns the goodness-of-fit (GOF), using the sum of
least squares, or other GOF statistic, of two arrays. Originally, the sum of
least squares statistic was used for simplicity, though this is currently not
implemented. Instead, a related method is used, but for any given age/sex/rob
category, the target value (number of cases or notification rate) and simulation
output value are both divided by the target value. This is done because the
notification rates (and case numbers) are essentially on different scales for
different categories. For example, SSA notification rates are around 250 per
100,000 each year, for the age category used for fitting. In UK-born, for the
same age category, notification rates would be on the order of 6 per 100,000
each year. In addition, an R-squared option can be chosen or a log-likelihood
deviance.

ENTRY: 'si' contains values from simulation output.
       'tar2' contains the calibration targets or observed data.
       'STAT' contains the type of GOF statistic to be used.
         0 Sum of squared deviations, relative to the observed value.
         1 Log likelihood deviance function.
         2 R-squared, fraction of variance explained among categories.
X        3 Normalized squared deviation from the mean.
X        4 Normalized absolute deviation from the median (not implemented).

EXIT:  'gof' returns the goodness-of-fit of 'si' to 'tar2'.

*/

#define STAT 2
#define INF  1E100

float gof(float si[ROB][YR][SEX][AC])
{ float x,x1,x2,rs,rn; int r,y,s,a;       X

// --- SQUARED DEVIATION FITTING

//- fprintf(fout,"About to calculate GOF\n"); fflush(fout);
  if(STAT==0)                       //Get discrepancy between model and data
  { x=0;                            //using sum of least squares method.**
    for(r=0; r<2;     r++)          //First add discrepancy for for non-UK and
    for(y=0; y<YR;    y++)          //UK-born, with SSA below due to age classes.
    for(s=0; s<2;     s++)          //**NOTE: Not true LS method, see note
    for(a=1; a<AC;    a++)          //in text at beginning of function.
      x+=pow(1-si[r][y][s][a]/tar2[r][y][s][a],2);
    if(ROB==3)                     //Add Sub-Saharan African (one age class).
    { r=2; a=1;
      for(y=0; y<YR;  y++)
      for(s=0; s<2; s++)
        x+=pow(1-si[r][y][s][a]/tar2[r][y][s][a],2); }
    return x; }

// --- POISSON DEVIANCE FITTING

  else if(STAT==1)                 //Get discrepancy between model and data using
  { x=0;                           //loglikelihood 'deviance' method. First add
    for(r=0; r<2;     r++)         //for non-UK and UK-born.
    for(y=0; y<YR;    y++)
    for(s=0; s<2;     s++)
    for(a=1; a<AC;    a++)
      x+=tar2[r][y][s][a]-si[r][y][s][a]+tar2[r][y][s][a]*log(si[r][y][s][a])-tar2[r][y][s][a]*log(tar2[r][y][s][a]);
    if(ROB==3)                     //Add Sub-Saharan African (one age class).
    { r=2; a=1;
      for(y=0; y<YR;  y++)
      for(s=0; s<2;   s++)
        x+=tar2[r][y][s][a]-si[r][y][s][a]+tar2[r][y][s][a]*log(si[r][y][s][a])-tar2[r][y][s][a]*log(tar2[r][y][s][a]);
    }
    x=-2*x; return x; }             //Finish 'deviance' formula.

// --- R-SQUARED FITTING

X else if(STAT==2)
X { rs=rn=0; forROB forSEX forAC             //Loop over region, sex, and age.
X   {
X     x=0; forYR                             //Sum the squared deviations from
X       x += pow(si[r][y][s][AR]             //the observed values across
X            - tar2[r][y][s][AR],2);         //years.
X
X     rn += 1; x1 = sim4[r][s][AR][qS2];     //Accumulate the R^2 values,
X     rs += x==0? 1:                         //taking care to avoid dividing
X          x1==0? -INF:                      //when the overall variance is
X                 1-x/(x1*x1); }             //zero.
X
X   return -rs/rn; }                         //Return the average -R^2.

// --- NORMALIZED SQUARED DEVIATION FITTING
//
//     (This first rescales by subtracting mu and dividing by sigma.)

X else if(STAT==3)
X { x=0; forROB forSEX forAC forYR
X   { x1 = (  si[r][y][s][AR]-sim4[r][s][a][qS1]) / sim4[r][s][a][qS2];
X     x2 = (tar2[r][y][s][AR]-sim4[r][s][a][qS1]) / sim4[r][s][a][qS2];
X     x += pow(x1-x2, 2); }
X   return x; }

// --- OTHER FITTING, NOT SPECIFIED

X else return 0;                   //(Should never get here)
}



/*-----------------------------------------------------------------------------
OBTAIN DEVIANCE FROM DATA

This function will be minimized by amebsa routine. The function to be minimized
is some measure of the goodness-of-fit of TB simulation model output to the
observed data, or calibration targets.

ENTRY: 'tb32' is the TB IBM.
       'gof' calculates goodness-of-fit between model output and
         calibration targets (actually opposite, where the smaller value the
         better the fit to data).
       'tar' contains the calibration targets or observed data.
       'in' contains the input parameter values for this evaluation.

EXIT:  'funx' contains the goodness-of-fit between targets (data) and model
         output.
*/

char arg0[] = "                                            ";
char arg1[] = "                                            ";
char arg2[] = "                                            ";
char arg3[] = "                                            ";
char arg4[] = "                                            ";
char arg5[] = "                                            ";
char arg6[] = "                                            ";
char *argci[] = { arg0, arg1, arg2, arg3, arg4, arg5, arg6 };

static int step = 0;

float funx(float *in)
{
//- fprintf(fout,"Entering funx\n"); fflush(fout);
  int i,j,r,a,s,y,n,q; float x; double v,w;
  float sim[IND];                         //Holds average simulation output.
  float sim2[ROB][YR][SEX][AC];           //Average output in 4-d array.
  float sim3[ROB][YR][SEX][AC];           //Averages 'glob2' for testing.
  double global[NPROC][IND];              //Holds output from all runs of tb29.
  double glob2[NPROC][ROB][YR][SEX][AC];  //Holds global simulation output in
                                          //5-d array (for indexing).
  double local[IND];                      //Contains results from current
                                          //processor.
  double global2[NPROC][2];               //Holds GOF value from all processors.
  double local2[2];                       //Holds GOF value for one processor.
  double param[ndim+1];                   //Holds current parameters and GOF
                                          //for convergence test.
  double mini[ROB][YR][SEX][AC];          //Array of minimum values of model
                                          //output for one iteration.
  double maxi[ROB][YR][SEX][AC];          //Array of maximum values of model
                                          //output for one iteration.

  FILE *fp;                               //Create pointer to file pipe for
                                          //use with popen.

#define MAINIAC
#ifndef MAINIAC
  char call[1000];                         //String used for popen call.
  sprintf(call,                            //Call simulation program.
    "tb32 randseq=%d df=%.3f d1uk20=%.10f d2uk20=%.10f d3uk20=%.10f "
    "|grep ^! |sed -e 's/!/ /g'",
    -(cproc+1),in[1],in[2],in[3],in[4]);

//- fprintf(fout,"About to run popen call: %s\n",call); fflush(fout);

  fp = popen(call,"r");                    //Run simulation and save
  for(i=0; i<IND; i++)                     //output from the current
  { if(fscanf(fp,"%lf",&v)<1)              //processor to 'local'.
    { fprintf(fout, "*** Input failure on 'popen'\n"); exit(3); }
    local[i]=v; }
  pclose(fp);

#else /* MAINIAC */
  double *wl, *mainiac(int, char**);
  sprintf(argci[0], "fit5DummyProgramName");  //Create array of arguments to
  sprintf(argci[1], "randseq=%d", -(cproc+1));//pass to TB program.
  sprintf(argci[2], "df=%.3f",      in[1]);
  sprintf(argci[3], "d1uk20=%.10f", in[2]);
  sprintf(argci[4], "d2uk20=%.10f", in[3]);
  sprintf(argci[5], "d3uk20=%.10f", in[4]);
//- fprintf(fout,"About to make mainiac call\n"); fflush(fout);
  wl = mainiac(6, argci);                  //Save simulation results into array.
  for(i=0; i<IND; i++) local[i] = wl[i];
#endif /* MAINIAC */
//- fprintf(fout,"About to call mpCommon (1)\n"); fflush(fout);
  mpCommon(global,local,IND);              //Gather results from all processors.
//- fprintf(fout,"Just finished mpCommon call\n"); fflush(fout);

  for(i=0,n=0; n<nproc; n++)               //Put global results into properly
  { for(r=0; r<ROB; r++)                   //indexed array.
    for(y=0; y<YR;  y++)
    for(s=0; s<SEX; s++)
    for(a=0; a<AC;  a++)
      glob2[n][r][y][s][a]=global[n][i++];
    i=0; }

//- This should be the same as the above but simpler:
//-  for(i=0,n=0; i=0,n<nproc; n++)           //Put global results into properly
//-  for(r=0; r<ROB; r++)                     //indexed array.
//-  for(y=0; y<YR;  y++)
//-  for(s=0; s<SEX; s++)
//-  for(a=0; a<AC;  a++)
//-    glob2[n][r][y][s][a]=global[n][i++];

  for(r=0; r<ROB; r++)                     //Clear 'sim3' for averaging 'glob2'.
  for(y=0; y<YR;  y++)
  for(s=0; s<SEX; s++)
  for(a=0; a<AC;  a++)
    sim3[r][y][s][a]=0;

  for(n=0; n<nproc; n++)                   //Average 'glob2' for testing and put
  for(r=0; r<ROB; r++)                     //into 'sim3' for comparison with
  for(y=0; y<YR;  y++)                     //'sim2'.
  for(s=0; s<SEX; s++)
  for(a=0; a<AC;  a++)
    sim3[r][y][s][a]+=glob2[n][r][y][s][a]/nproc;

X #ifdef ACROSSYEARS
X forROB forSEX forAC forQ                        //Clear 'sim4', for mean and
X   sim4[r][s][a][q]=0;                           //variance across years.
X
X forROB forSEX forAC forYR forN                  //Sum powers 0, 1, and 2 of
X { sim4[r][s][a][qS0] += 1;                      //the data in each category,
X   sim4[r][s][a][qS1] += glob2[n][r][y][s][a];   //across years and processors.
X   sim4[r][s][a][qS2] += glob2[n][r][y][s][a]    //(SHOULD THIS BE THE VARIANCE
X                       * glob2[n][r][y][s][a]; } // IN THE ANNUAL AVERAGES?)
X
X forROB forSEX forAC                             //Convert the summed powers to
X { sim4[r][s][a][qS1] /= sim4[r][s][a][qS0];     //mean and standard deviation
X   sim4[r][s][a][qS2]  = sim4[r][s][a][qS2]      //(mu and sigma).
X                       / sim4[r][s][a][qS0]
X                       - sim4[r][s][a][qS1]
X                       * sim4[r][s][a][qS1];
X   sim4[r][s][a][qS2]  = sqrt(sim4[r][s][a][qS2]); }

X #else
X forROB forSEX forAC forQ                        //Clear 'sim4', for mean and
X   sim4[r][s][a][q]=0;                           //variance across years.
X
X forROB forSEX forAC                             //Step through all categories.
X {
X   vS0=vS1=vS2=0;
X   for(y=0; y<YR; y++)
X   {
X     wS0=wS1=wS2=0;                              //Sum powers 0, 1, and 2 of
X     for(n=0; n<nproc; n++)                      //the data in each category,
X     { wS0 += 1;                                 //for each year across
X       wS1 += glob2[n][r][y][s][a];              //processors.
X       wS2 += glob2[n][r][y][s][a]
X            * glob2[n][r][y][s][a]; }
X
X       wS1 /= wS0;                               //Convert the summed powers to
X       wS2  = sqrt(wS2/wS0 - wS1*wS1);           //mean and standard deviation,
X       vS0 += 1; vS1 += wS1; vS2 += wS2; }       //and accumulate them.
X
X   sim4[r][s][a][qS0]  = vS0;                    //Store away the overall mean
X   sim4[r][s][a][qS1]  = vS1/vS0;                //and the average standard
X   sim4[r][s][a][qS2]  = vS2/vS0; }              //deviation.
X #endif

  for(i=0; i<IND; i++)                     //Clear 'sim' before using it to
    sim[i]=0;                              //average simulation results.

  for(i=0; i<nproc; i++)                   //Average output and place inside
  for(j=0; j<IND;   j++)                   //'sim'.
    sim[j]+=global[i][j]/nproc;

  for(i=0,r=0; r<ROB; r++)                 //Transfer simulation results to
  for(y=0; y<YR;  y++)                     //four-dimensional array and transfer
  for(s=0; s<SEX; s++)                     //targets to float array.
  for(a=0; a<AC;  a++)
  { sim2[r][y][s][a]=sim[i++];
    tar2[r][y][s][a]=tar[r][y][s][a]; }

  for(r=0; r<ROB; r++)                     //Check that 'sim2' and 'sim3' are
  for(y=0; y<YR;  y++)                     //exactly the same.
  for(s=0; s<SEX; s++)
  for(a=0; a<AC;  a++)
    if(sim2[r][y][s][a]!=sim3[r][y][s][a])
    { fprintf(fout,"Error -- different averaged simulation results for 'sim2' & 'sim3'\n");
      exit(3); }

  for(r=0; r<ROB; r++)                     //Set 'mini' and 'maxi' to values of
  for(y=0; y<YR;  y++)                     //first processor, to begin with.
  for(s=0; s<SEX; s++)
  for(a=0; a<AC;  a++)
  { mini[r][y][s][a] = glob2[0][r][y][s][a];
    maxi[r][y][s][a] = glob2[0][r][y][s][a]; }

  for(n=1; n<nproc; n++)                   //Find minimum and maximum for each
  for(r=0; r<ROB;   r++)                   //value returned by the different
  for(y=0; y<YR;    y++)                   //processors.
  for(s=0; s<SEX;   s++)
  for(a=0; a<AC;    a++)
  { if(glob2[n][r][y][s][a]< mini[r][y][s][a])
      mini[r][y][s][a] = glob2[n][r][y][s][a];
    if(glob2[n][r][y][s][a]> maxi[r][y][s][a])
      maxi[r][y][s][a] = glob2[n][r][y][s][a]; }

  x=gof(sim2);                             //Obtain GOF and check GOF is the
  local2[0]=x; local2[1]=0;                //same for all processors.
  mpCommon(global2,local2,2);
  for(i=1; i<nproc; i++)
    if(global2[i-1][0]!=global2[i][0])
    { fprintf(fout,"Error -- different GOF results\n"); exit(3); }

  for(i=0; i<ndim; i++)                    //Copy input parameters and GOF
    param[i] = in[i+1];                    //together for a convergence test,
  param[ndim] = x;                         //and run the test.
  w = converge(step++, ndim+1, param);

  if(step-1==0)                            //Display Centinel headings.
    fprintf(fout, "|n |df |d1uk20 |d2uk20 |d3uk20 |gof |w\n");

  fprintf(fout,                            //Display the results.
         "|%3d |%.5f |%.14f |%.14f |%.14f |%.1f |%.4f\n",
        step-1, in[1],in[2], in[3], in[4], x,    w); fflush(fout);

//- fprintf(fout,"About to print to 'plot'\n"); fflush(fout);
  if(cproc==0)                             //If this is the first processor,
  { fprintf(plot,                          //write conditions for this step
          "step=%d df=%.3f d1uk10=%.12f "  //in tex file and write simulation
          "d2uk10=%.12f   d3uk10=%.12f GOF=%.3f\n",     //output using 'simo'.
             step-1,in[1],in[2],in[3],in[4],x); fflush(plot);
    simo(sim2,glob2,mini,maxi); }

  if(w>=0 && w<THRESH)                     //If convergence has been achieved,
  { fprintf(fout,                          //terminate processing.
           "Convergence achieved\n"
           "***BEST*** yb=%.3f pb=%.3f,%.14f,%.14f,%.14f\n",
           yb, pb[1],pb[2],pb[3],pb[4]);
    fflush(fout);
    fclose(fout);                          //(Close annealing output file.)
    if(cproc==0)                           //(Close simulation output file.)
    { fprintf(plot,"\\end \n");
      fprintf(plot,"***BEST*** yb=%f pb=%f,%f,%f,%f\n",
      yb, pb[1],pb[2],pb[3],pb[4]); fflush(plot);
      fclose(plot); }
    mpClose();                             //(Close parallel processing.)
    exit(0); }                             //(Return to operating system.)
//- fprintf(fout,"About to return GOF\n"); fflush(fout);
  return x;                                //Otherwise return goodness of fit.
}

#ifdef CLTESTS
/*-----------------------------------------------------------------------------
NEW MATERIAL FOR THE ABOVE ROUTINE

This section is what I am doing to see if I can piggyback some of my own work on
the large-scale runs, comparing a few with a Bayesian step at the end. It picks
up when the data are about to be collected from the global array and places in a
properly indexed array. It uses a new array 'glob3', ordered so that its
subarrays can be passed to routines like 'sort' and 'RandF'.  ---2011/09/24, CLL
*/

#define qN      0    //Index for the number of observations.
#define qS1     1    //Index for the sum of the observations.
#define qS2     2    //Index for the sum of the squares of the observations.
#define DSTAT 1000   //Number of bootstrapped fitness samples.

#define forROB  for(r=0; r<ROB;   r++)
#define forSEX  for(s=0; s<SEX;   s++)
#define forAC   for(a=0; a<AC;    a++)
#define forYR   for(y=1; y<YR+1;  y++)
#define forN    for(n=0; n<nproc; n++)
#define forQ    for(q=0; q<qS2+1; q++)

dec qstat[ROB][SEX][AC][qS2+1];      //Normalization statistics.
dec glob3[ROB][SEX][AC][YR][NPROC];  //Same as 'glob2' but ordered for calling.

  for(i=0,n=0; i=0,n<nproc; n++)                  //Put global results into
  forROB forYR forSEX forAC                       //a properly indexed array.
    glob3[r][s][a][y][n+1] = global[n][i++];

  forROB forSEX forAC forYR                       //Sort replicates across
    qmed[r][s][a][y] =                            //processors to ascending
      psort(glob3[r][s][a][y], nproc+2);          //order and record the median.

  forROB forSEX forAC forQ                        //Clear the array of
    qstat[r][s][a][q] = 0;                        //statistics.

  forROB forSEX forAC forYR forN                  //Sum the zeroth, first, and
    qstat[r][s][a][qN]  += 1,                     //second powers of the data
    qstat[r][s][a][qS1] += glob3[r][s][a][y][n],  //in each category, across
    qstat[r][s][a][qS2] += glob3[r][s][a][y][n]   //years and processors.
                         * glob3[r][s][a][y][n];

  forROB forSEX forAC                             //Convert the summed powers to
    qstat[r][s][a][qS1] /= qstat[r][s][a][qN],    //mean and standard deviation.
    qstat[r][s][a][qS2]  = sqrt(
                           qstat[r][s][a][qS2]/qstat[r][s][a][qN]
                         - qstat[r][s][a][qS1]*qstat[r][s][a][qS1]);

  forROB forSEX forAC forYR forN                  //Scale all the data by
    glob3[r][s][a][y][n] -= qstat[r][s][a][qS1],  //subtracting the mean and
    glob3[r][s][a][y][n] /= qstat[r][s][a][qS2];  //dividing by the dispersion.

  g0 = 0; forROB forSEX forAC                     //Compute the primary goodness
    g0 += gofc(glob3[r][s][a], tar1[r][s][a]);    //of fit for each category and
  g0 /= ROB*SEX*AC;                               //across all categories.
  if(DIRECT) return g0;

  forROB forSEX forAC forYR forN                  //Convert the replicates to
    glob3[r][s][a][y][n]+=glob3[r][s][a][y][n-1]; //cumulative distributions and
  forROB forSEX forAC forYR                       //balance the distributions.
    balance(glob3[r][s][a][y], nproc+2);

  if(glob3y[1]==0)                                //If this is the first time
  { for(i=1; i<nproc+1; i++)                      //through, create y-values for
      glob3y[i] = (dec)i/(nproc+2);               //all the distributions.
    glob3y[nproc+1] = 1; }

  for(i=0; i<DSTAT; i++)                          //Generate numerous random
  { gn = 0; forROB forSEX forAC forYR             //variants on this set of runs
      glob3b[r][s][a][y] = RandF(...),            //and calculate the fitness
      gn += gofc(glob3b[r][s][a], tar1[r][s][a]); //statistic for each.
    gv[i] = g0; gi[i+1] = i; }

  i = sort(gi, 1, DSTAT-1, order);                //Locate the prime statistic
  for(n=0; gi[i]; n+=1,i=gi[i])                   //within the distribution of
    if(g0<=gv[i]) break;                          //statistics.

  ....                                            //Convert its position to a
                                                  //probability.
#endif /* CLTESTS */

/*----------------------------------------------------------------------------
PLOT SIMULATION RESULTS

This routine is called from 'funx' and writes average simulation results for
each evaluation step to a Tex-formatted graphics file which can be typeset into
a PDF file.

ENTRY: 'step' contains the step number, starting with 0.
       'tar2' contains the calibration targets, from observed data.
       'sim' contains the simulation output data.
       'glob2' contains a five-dimensional array of all simulation
         results, indexed on 'rob', year, age class, and sex.
       'mini' contains the array of minimum values for each data point of
         simulation output.
       'maxi' contains the array of maximum values for each data point of
         simulation output.

EXIT:  'plot' contains graphics output in Pictex format.
*/
#define RANGE 0                         //0=No ranges printed, 1=ranges printed.
#define PLOTALL 1                       //0=Do not print all results, 1=do.
#define VR  (max>125?50: max>50?25: 10) //Vertical resolution on graph size.
#define HS   4.0                        //Horizontal graph size, inches.
#define HU  10.0                        //Horizontal graph units.
#define Phi  1.618                      //Horizontal to vertical graph ratio.
#define CBL  "{"                        //Left curly brace, for balancing.
#define CBR  "}"                        //Right curly brace, for balancing.

static int figure = 0;                  //Running figure number.

static char *wsex[] = { "Male", "Female" };
static char *wage[] = { "black ", "blue ", "green", "red " };
static char *wrob[] = { "Non-UK born", "UK-born",  "SSA-born" };

void simo(float sim[ROB][YR][SEX][AC], double glob2[NPROC][ROB][YR][SEX][AC],double mini[ROB][YR][SEX][AC],double maxi[ROB][YR][SEX][AC])
{
  int i,n,k,r,y,s,a; double max;

//- fprintf(fout,"Entering 'simo'\n"); fflush(fout);

  fprintf(plot,"\\overfullrule=0pt\n");

  for(r=ROB-1; r>=0; r--)               //Determine the vertical scale for the
  for(s=0;     s<=1; s++)               //current graph (auto-scaling).
  { max = -1E9;
    for(a=0; a<AC; a++)
    for(y=0; y<YR; y++)
    { if(max<sim[r][y][s][a]) max = sim[r][y][s][a];
      if(max<tar[r][y][s][a]) max = tar[r][y][s][a]; }
    k = (max+1.16*VR); k = k/VR*VR;

    fprintf(plot,"%s\\graph\n", CBL);   //Open the current graph.
    fprintf(plot,"$$ \\beginpicture\n");
    fprintf(plot,"\\setcoordinatesystem units < %g in, %g in>\n", HS/HU, (HS/Phi)/k);
    fprintf(plot,"\\setplotarea x from 0 to 10, y from 0 to %d \n", k);
    fprintf(plot,"\\plotheading {%s %s Notification Rates}\n", wrob[r], wsex[s]);
    fprintf(plot,"\\axis bottom %%label {Year}\n");
    fprintf(plot,"  ticks withvalues {1999} {2001} {2003} {2005} {2007} {2009} /\n");
    fprintf(plot,"  at 0 2 4 6 8 10 / /\n");
    fprintf(plot,"\\axis left label {\\vertical{\\lines{Cases per 100,000\\cr\\cr}}}\n");
    fprintf(plot,"  ticks numbered from 0 to %d by %d /\n", k, k<=10?2: k<=30?5: k<=150?25: 50);
    fprintf(plot,"\\axis right / \\axis top /\n");

    for(a=0; a<AC; a++)                 //Display observed values.
    { fprintf(plot,"\\%s%s\\plot ", wage[a], CBL);
      for(y=0; y<YR; y++) fprintf(plot,"%2d %6.2f\t ", y,tar[r][y][s][a]);
      fprintf(plot," / %s%%\n", CBR); }

    fprintf(plot, "\\setdashes\n");     //Display simulated values.
    for(a=0; a<AC; a++)
    { fprintf(plot,"\\%s%s\\plot ", wage[a], CBL);
      for(y=0; y<YR; y++) fprintf(plot,"%2d %6.2f\t ", y,sim[r][y][s][a]);
      fprintf(plot," / %s%%\n", CBR); }

    if(RANGE)
    { fprintf(plot,"\\setsolid\n");     //Display for minimum and maximum
      for(a=0; a<AC; a++)               //for multiple simulations.
      for(y=0; y<YR; y++)
      { fprintf(plot,"\\%s%s\\plot ", wage[a], CBL);
        fprintf(plot,"%2d %.2lf\t ", y,mini[r][y][s][a]);
        fprintf(plot,"%2d %.2lf\t ", y,maxi[r][y][s][a]);
        fprintf(plot," / %s%%\n", CBR); }
    }

    if(PLOTALL)
    { for(a=0; a<AC; a++)               //Display an 'x' for every individual
      for(y=0; y<YR; y++)               //simulation replicate.
      { fprintf(plot,"\\%s%s\\multiput {x} at ", wage[a], CBL);
        for(n=0; n<nproc; n++)
          fprintf(plot,"%2d %.2lf\t ", y,glob2[n][r][y][s][a]);
        fprintf(plot," / %s%%\n", CBR); }
    }

    fprintf(plot,"\\setsolid\n"         //Close the current graph.
                 "\\endpicture $$%s\n", CBR);

    fprintf(plot,"\\vskip-0pt \\noindent{\\bf Figure %d.} "
                 "Correspondence between simulated values (dashed) "
                 "versus observed data (solid) for age classes"
                 " 0--14 (black), "
                 "15--44 (blue), "
                 "45--64 (green), "
                 "and 65+ (red).\n"
                 "\\vskip 24pt\n\n", ++figure);

    if(figure%2==0)                     //Force a new page every other graph.
      fprintf(plot,"\\vfill\\eject\n"); }
}



/*------------------------------------------------------------------------------
DISPLAY INTERMEDIATE RESULTS

This routine is called from 'amebsa' to display progress. Results are put in a
form that gets typeset into an animated PDF file by TEX.
*/

simplex(int ndim, float temptr, float **p, float *y)
{
  int i;
 /*
  printf("\\vfill\\eject\n");
  printf("$$ \\beginpicture\n");
  printf("\\setcoordinatesystem units <.5in, .5in>\n");
  printf("\\setplotarea x from -6 to 6, y from -6 to 6\n");
  printf("\\plotheading {$T=%.4f$}\n", temptr);
  printf("\\putrule from -6 0 to 6 0\n");
  printf("\\putrule from 0 -6 to 0 6\n");
  printf("\\put {$\\bullet$} at 0 0\n");
  printf("\\put {$\\bullet$} at 2 2\n");

  printf("\\plot");                          //Note: plotting should be fixed
  for(i=1; i<=ndim+1; i++)                   //for the 4 parameters!!!
    printf(" %f %f ", p[i][1], p[i][2]);
  printf(" %f %f /\n", p[1][1], p[1][2]);

  for(i=1; i<=ndim+1; i++)
    printf("\\put {%.3g} at %f %f\n", y[i], p[i][1], p[i][2]);

  printf("\\endpicture $$\n\n");
*/
}

// ADRIENNE KEEN, APRIL 2011.