/*
 Copyright (C) 2000-2004

 Code contributed by Greg Collecutt, Joseph Hope and Paul Cochrane

 This file is part of xmds.
 
 This program is free software; you can redistribute it and/or
 modify it under the terms of the GNU General Public License
 as published by the Free Software Foundation; either version 2
 of the License, or (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/

/*
  $Id: xmdssimulation.cc,v 1.52 2005/10/27 00:58:39 joehope Exp $
*/

/*! @file xmdssimulation.cc
  @brief Simulation element parsing classes and methods

  More detailed explanation...
*/

#include<xmlbasics.h>
#include<dom3.h>
#include<xmdsutils.h>
#include<xmdsclasses.h>
#include<config.h>

// ******************************************************************************
// ******************************************************************************
//                              xmdsSimulation public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsSimulations=0;   //!< The number of xmds simulation objects

// ******************************************************************************
xmdsSimulation::xmdsSimulation(
                               const char *const yourRawFileName,
                               const bool& yourVerboseMode,
                               const bool& mpiAvailable) :
  xmdsElement(this,yourVerboseMode){
  if(debugFlag) {
    nxmdsSimulations++;
    printf("xmdsSimulation::xmdsSimulation\n");
    printf("nxmdsSimulations=%li\n",nxmdsSimulations);
  }

  myParameters.rawFileName = yourRawFileName;
  myParameters.stochastic = 0;
  myParameters.nThreads = 1;
  myParameters.nPaths = 1;
  myParameters.mpiAvailable = mpiAvailable;
  myParameters.usempi = 0;
  myParameters.bing = 0;
  myParameters.seed[0] = 1;
  myParameters.seed[1] = 2;
  myParameters.nNoises = 0;
  myParameters.noiseKind = "gaussian";
  myParameters.noiseMean = 5.0;
  myParameters.errorCheck = 1;
  myParameters.useWisdom = 0;
  myParameters.binaryOutput = 0;
  myParameters.useDouble = 1;
  myParameters.benchmark = 0;
  myParameters.version = VERSION;
  myParameters.release = RELEASE;
  myCurrentSegmentNumber = 0;
};

// ******************************************************************************
xmdsSimulation::~xmdsSimulation() {
  if(debugFlag) {
    nxmdsSimulations--;
    printf("xmdsSimulation::~xmdsSimulation\n");
    printf("nxmdsSimulations=%li\n",nxmdsSimulations);
  }
};

// ******************************************************************************
void xmdsSimulation::processElement(
                                    const Element *const yourElement) {
  if(debugFlag) {
    printf("xmdsSimulation::processElement\n");
  }

  // ************************************
  // process parameters
  // ************************************

  list<XMLString> myXMLStringList;
  list<bool> myBoolList;
  list<unsigned long> myULongList;

  // ************************************
  // find 'name'

  getAssignmentStrings(yourElement,"name",0,1,myXMLStringList);

  if(myXMLStringList.size()==1) {
    myParameters.simulationName=*myXMLStringList.begin();
    if(verbose()) {
      printf("simulation name = '%s'\n",myParameters.simulationName.c_str());
    }
  }
  else {
    // Create default simulation file name
    unsigned long lastdot=0;
    const XMLString tempName = myParameters.rawFileName;
    unsigned long i=0;
    while(i<tempName.length()) {
      if(tempName.data(i)=='.') {
        lastdot=i;
      }
      i++;
    }
    if(lastdot>0) {
      tempName.subString(myParameters.simulationName,0,lastdot);
    }
    else {
      myParameters.simulationName=tempName;
    }
    if(verbose()) {
      printf("simulation name defaulting to '%s'\n",myParameters.simulationName.c_str());
    }
  }

  // ************************************
  // find 'prop_dim'

  getAssignmentStrings(yourElement,"prop_dim",1,1,myXMLStringList);

  myParameters.propDimName=*myXMLStringList.begin();
  if(verbose()) {
    printf("simulation prop_dim = '%s'\n",myParameters.propDimName.c_str());
  }


  // ************************************
  // find 'author'

  getAssignmentStrings(yourElement,"author",NOT_REQD,0,myXMLStringList);

  if (myXMLStringList.size() > 0) {

    myParameters.authorName = *myXMLStringList.begin();
    myXMLStringList.pop_front();
    for(list<XMLString>::const_iterator pXMLString = myXMLStringList.begin(); pXMLString != myXMLStringList.end(); pXMLString++) {
      myParameters.authorName += " ";
      myParameters.authorName += *pXMLString;
    }
    if(verbose()) {
      printf("simulation author = '%s'\n",myParameters.authorName.c_str());
    }
  }
  else {
    // this warning may need to be taken out somehow, but I sort of want people
    // to be nice little coders and document their code nicely, and this is
    // one way to do it...  (PTC)
    printf("No <author> tag found.  It's not required, but it's a Good Idea.\n");
  }

  // ************************************
  // find 'description'

  getAssignmentStrings(yourElement,"description",NOT_REQD,0,myXMLStringList);

  // Storing the description can cause overflow errors if it is too long, so we'll
  // comment out the actual loading of the description, and put into the description
  // variable a note that it actually exists.
  if (myXMLStringList.size() > 0) {
    
    myParameters.description += "Description found.  See xmds file for the rest of it.";
    /*
    myParameters.description = *myXMLStringList.begin();
    myXMLStringList.pop_front();
    for(list<XMLString>::const_iterator pXMLString = myXMLStringList.begin(); pXMLString != myXMLStringList.end(); pXMLString++) {
      myParameters.description += " ";
      myParameters.description += *pXMLString;
    }
    */
    if(verbose()) {
      printf("simulation description = '%s'\n",myParameters.description.c_str());
    }
  }
  else {
    // this warning may need to be taken out somehow, but I sort of want people
    // to be nice little coders and document their code nicely, and this is
    // one way to do it...  (PTC)
    printf("No <description> tag found.  It's not required, but it's a Good Idea.\n");
  }

  // ************************************
  // find 'error_check'

  getAssignmentBools(yourElement,"error_check",0,1,myBoolList);

  if(myBoolList.size()==1) {

    myParameters.errorCheck=*myBoolList.begin();

    if(verbose()) {
      if(myParameters.errorCheck) {
        printf("simulation error_check is 'yes'\n");
      }
      else {
        printf("simulation error_check is 'no'\n");
      }
    }
  }
  else {

    if(verbose()) {
      printf("simulation error_check defaulting to 'yes'\n");
    }

    myParameters.errorCheck=1;  
  }

  // ************************************
  // find 'use_wisdom'

  getAssignmentBools(yourElement,"use_wisdom",0,1,myBoolList);

  if(myBoolList.size()==1) {

    myParameters.useWisdom=*myBoolList.begin();

    if(verbose()) {
      if(myParameters.useWisdom) {
        printf("simulation use_wisdom is 'yes'\n");
      }
      else {
        printf("simulation use_wisdom is 'no'\n");
      }
    }
  }
  else {

    if(verbose()) {
      printf("simulation use_wisdom defaulting to 'no'\n");
    }

    myParameters.useWisdom=0;
  }

  // ************************************
  // find 'use_prefs'

  getAssignmentBools(yourElement,"use_prefs",0,1,myBoolList);

  if(myBoolList.size()==1) {

    myParameters.usePrefs=*myBoolList.begin();

    if(verbose()) {
      if(myParameters.usePrefs) {
        printf("simulation use_prefs is 'yes'\n");
      }
      else {
        printf("simulation use_prefs is 'no'\n");
      }
    }
  }
  else {

    if(verbose()) {
      printf("simulation use_prefs defaulting to 'yes'\n");
    }

    myParameters.usePrefs=1;
  }

  // ***********************************
  // find binary output flag

  getAssignmentBools(yourElement,"binary_output",0,1,myBoolList);

  if(myBoolList.size()==1) {
    // now we need to print a warning and tell the user what they can do to fix the problem
    sprintf(errorMessage(),"Using the <binary_output> tag is obsolete.\n"
	    "  Please replace this tag and use the format=\"ascii|binary\"\n"
	    "  attribute of the <output> tag instead.");
    throw xmdsException(errorMessage());  

    // this is the old code
    /*
    myParameters.binaryOutput=*myBoolList.begin();
    
    if(verbose()) {
      if(myParameters.binaryOutput) {
        printf("simulation binary_output is 'yes'\n");
      }
      else
        printf("simulation binary_output is 'no'\n");
    }
    */
  }
  else {

    if(verbose()) {
      printf("simulation binary_output defaulting to 'no'\n");
    }
        
    myParameters.binaryOutput=0;
  }

  // ***********************************
  // find use_double element

  getAssignmentBools(yourElement,"use_double",0,1,myBoolList);

  if(myBoolList.size()==1) {
    // now we need to print a warning and tell the user what they can do to fix the problem
    sprintf(errorMessage(),"Using the <use_double> tag is obsolete.\n"
	    "  Please replace this tag and use the precision=\"double|single\"\n"
	    "  attribute of the <output> tag instead.");
    throw xmdsException(errorMessage());  

    // this is the old code
    /*
    myParameters.useDouble=*myBoolList.begin();

    if(verbose()) {
      if(myParameters.useDouble) {
        printf("simulation use_double is 'yes'\n");
	}
      else
        printf("simulation use_double is 'no'\n");
    }
    */
  }
  else {

    if(verbose()) {
      printf("simulation use_double defaulting to 'yes'\n");
    }

    myParameters.useDouble=1;
  }

  // ***********************************
  // find benchmark element
        
  getAssignmentBools(yourElement,"benchmark",0,1,myBoolList);
        
  if(myBoolList.size()==1) {
    myParameters.benchmark=*myBoolList.begin();

    if(verbose()) {
      if(myParameters.benchmark) {
        printf("simulation benchmark is 'yes'\n");
      }
      else {
        printf("simulation benchmark is 'no'\n");
      }
    }
  }
  else {

    if(verbose()) {
      printf("simulation benchmark defaulting to 'no'\n");
    }

    myParameters.benchmark=0;
  }


  // ************************************
  // find 'threads'

  getAssignmentULongs(yourElement,"threads",0,1,myULongList);

  if(myULongList.size()==1) {

    myParameters.nThreads = *myULongList.begin();

    if(myParameters.nThreads < 1) {
      throw xmdsException(yourElement,"number of threads must be > 0");
    }

    if(verbose()) {
      printf("number of threads = %li\n",myParameters.nThreads);
    }
  }
  else {

    if(verbose()) {
      printf("number of threads defaulting to 1\n");
    }

    myParameters.nThreads=1;  
  }

  // ************************************
  // find 'stochastic'

  getAssignmentBools(yourElement,"stochastic",0,1,myBoolList);

  if(myBoolList.size()==1) {

    myParameters.stochastic=*myBoolList.begin();

    if(verbose()) {
      if(myParameters.stochastic) {
        printf("simulation stochastic is 'yes'\n");
      }
      else {
        printf("simulation stochastic is 'no'\n");
      }
    }
  }
  else {

    if(verbose()) {
      printf("simulation stochastic defaulting to 'no'\n");
    }

    myParameters.stochastic=0;  
  }

  // ************************************
  // find 'use_mpi' assignment

  getAssignmentBools(yourElement,"use_mpi",0,1,myBoolList);

  if(myBoolList.size()==1) {

    myParameters.usempi=*myBoolList.begin();

    if(myParameters.usempi&!myParameters.mpiAvailable) {
      throw xmdsException(yourElement,"MPI was not available when xmds was installed");
    }

    if(verbose()) {
      if(myParameters.usempi) {
        printf("use_mpi is 'yes'\n");
      }
      else {
        printf("use_mpi is 'no'\n");
      }
    }
  }
  else {
    if(verbose()) {
      printf("use_mpi defaulting to 'no'\n");
    }

    myParameters.usempi=0;  
  }
        
  if((myParameters.usempi)&(myParameters.nThreads>1)&(!myParameters.stochastic)) {
    throw xmdsException(yourElement,"Non-determinstic MPI and threaded FFTs are not compatible");
  }
        
        
  // ************************************
  // find 'bing' assignment

  getAssignmentBools(yourElement,"bing",0,1,myBoolList);

  if(myBoolList.size()==1) {

    myParameters.bing=*myBoolList.begin();

    if(verbose()) {
      if(myParameters.bing) {
        printf("bing is 'yes'\n");
      }
      else {
        printf("bing is 'no'\n");
      }
    }
  }
  else {
    if(verbose()) {
      printf("bing defaulting to 'no'\n");
    }

    myParameters.bing=0;  
  }

  // ************************************
  // if stochastic get paths, and seed assignments
  if(myParameters.stochastic) {

    // ************************************
    // get paths

    getAssignmentULongs(yourElement,"paths",1,1,myULongList);
    myParameters.nPaths = *myULongList.begin();

    if(myParameters.nPaths < 1) {
      throw xmdsException(yourElement,"number of paths must be > 0");
    }

    if(verbose()) {
      printf("number of paths = %li\n",myParameters.nPaths);
    }


    // ************************************
    // get seeds

    getAssignmentULongs(yourElement,"seed",1,2,myULongList);

    list<unsigned long>::const_iterator pULong = myULongList.begin();
    myParameters.seed[0] = *pULong;
    pULong++;
    myParameters.seed[1] = *pULong;

    if((myParameters.seed[0] < 0)|(myParameters.seed[1] < 0)) {
      throw xmdsException(yourElement,"seeds must be positive");
    }

    if(myParameters.seed[0]==myParameters.seed[1]) {
      throw xmdsException(yourElement,"seeds are indentical");
    }

    if(verbose()) {
      printf("seeds are %li and %li\n",myParameters.seed[0],myParameters.seed[1]);
    }

    // ************************************
    // get noises

    getAssignmentULongs(yourElement,"noises",1,1,myULongList);

    myParameters.nNoises = 2*((*myULongList.begin()+1)/2);

    if(myParameters.nNoises < 2) {
      throw xmdsException(yourElement,"number of noises must be > 0");
    }

    if(verbose()) {
      printf("number of noises = %li\n",myParameters.nNoises);
    }

    // find out what kind of noise is specified, if any
    getAttributeStrings(yourElement, "noises", "kind", NOT_REQD, myParameters.noiseKind);
    // the kind has several values, current possible values are:
    // (this list may be augmented in time)
    // "gaussian"
    // "gaussFast"
    // "poissonian"
    // "uniform"
    if (debugFlag) {
      printf("After getAttributeStrings(), noiseKind = %s\n",myParameters.noiseKind.c_str()); 
    }
     
    if (myParameters.noiseKind == "") {
      if (verbose()) {
        printf("Setting the default noise kind (i.e. \"gaussian\")\n");
      }
      myParameters.noiseKind = "gaussian";  // this is the default option
    }
    else if (myParameters.noiseKind == "gaussian") {
      if (verbose()) {
        printf("Noise kind set to \"gaussian\"\n");
      }
    }
    else if (myParameters.noiseKind == "gaussFast") {
      if (verbose()) {
        printf("Noise kind set to \"gaussFast\"\n");
      }
    }
    else if (myParameters.noiseKind == "poissonian") {
      if (verbose()) {
        printf("Noise kind set to \"poissonian\"\n");
      }
      // if we have poissonian specified, we must have a (nonzero) mean for the
      // distribution specified as well
      XMLString noiseMeanString;
      getAttributeStrings(yourElement, "noises", "mean", REQD, noiseMeanString);
      if (debugFlag) {
	printf("After getAttributeStrings(), noiseMean = %s\n",noiseMeanString.c_str()); 
      }      
      // now convert the string to its relevant floating point value
      myParameters.noiseMean = atof(noiseMeanString.c_str());
      if (verbose()) {
	printf("Mean for Poissonian noise set to %g\n",myParameters.noiseMean);
      }
    }
    else if (myParameters.noiseKind == "uniform") {
      if (verbose()) {
        printf("Noise kind set to \"uniform\"\n");
      }
    }
    else {
      throw xmdsException(yourElement, 
                          "You must specify a correct noise kind.  Current choices are:\n gaussian\n gaussFast\n poissonian\n uniform");
    }
	
	// find out which scheduling method to use
	if(myParameters.usempi) {
		
		getAssignmentStrings(yourElement,"MPI_Method",0,1,myXMLStringList);
		
		if(myXMLStringList.size()==0) {
			// Default mpi method is "scheduling"
			myXMLStringList.push_back("Scheduling");
			myParameters.mpiMethod=*myXMLStringList.begin();
			if(verbose()) {
				printf("MPI method defaulting to '%s'\n",myParameters.mpiMethod.c_str());
			}
		}
		else {
			if(!((*myXMLStringList.begin()=="Scheduling")||(*myXMLStringList.begin()=="Uniform"))) {
				sprintf(errorMessage(),"MPI method '%s' unknown.\nCurrent methods are 'Scheduling' and 'Uniform'.",myXMLStringList.begin()->c_str());
				throw xmdsException(yourElement,errorMessage());
			}
			myParameters.mpiMethod=*myXMLStringList.begin();
			if(verbose()) {
				printf("MPI method = '%s'\n",myParameters.mpiMethod.c_str());
			}
		}
	}


  }

  // ************************************
  // find and process children
  // ************************************

  const NodeList* candidateElements;

  // ************************************
  // find argv element and process

  candidateElements = yourElement->getElementsByTagName("argv",0);

  if(candidateElements->length()>1) {
    throw xmdsException(yourElement,"Multiple <argv> elements!");
  }
  else if(candidateElements->length()==1) {
    myArgv = createxmdsArgv();
    const Element* yourElement = dynamic_cast<const Element*> (candidateElements->item(0));
    myArgv->processElement(yourElement);
  }

  // ************************************
  // find globals element and process if present

  candidateElements = yourElement->getElementsByTagName("globals",0);

  if(candidateElements->length()>1) {
    throw xmdsException(yourElement,"Multiple <globals> elements!");
  }
  else if(candidateElements->length()==1) {
    xmdsElement* newGlobals = createxmdsGlobals();
    const Element* yourElement = dynamic_cast<const Element*> (candidateElements->item(0));
    newGlobals->processElement(yourElement);
  }

  // ************************************
  // find field element and process

  candidateElements = yourElement->getElementsByTagName("field",0);

  if(candidateElements->length()>1) {
    throw xmdsException(yourElement,"Multiple <field> elements!");
  }
  else if(candidateElements->length()==1) {
    myField = createxmdsField();
    const Element* yourElement = dynamic_cast<const Element*> (candidateElements->item(0));
    myField->processElement(yourElement);
  }
  else
    throw xmdsException(yourElement,
                        "Cannot find <field> element");
                                
  // ************************************
  // find output element and process

  candidateElements = yourElement->getElementsByTagName("output",0);

  if(candidateElements->length()>1) {
    throw xmdsException(yourElement,"Multiple <output> elements!");
  }
  else if(candidateElements->length()==1) {
    // now grab the format attribute if available
    XMLString outputFormat = "ascii"; // set it to the default while declaring it
    const Node *thisElement = candidateElements->item(0);
    getAttributeStrings(yourElement, "output", "format", NOT_REQD, outputFormat);
    if (outputFormat == "") {
      outputFormat = "ascii";  // the default option
    }
    else {
      if (outputFormat == "ascii") {
        myParameters.binaryOutput = false;
        if (verbose()) {
          printf("Output format set to 'ascii' in output element\n");
        }
      }
      else if (outputFormat == "binary") {
        myParameters.binaryOutput = true;
        if (verbose()) {
          printf("Output format set to 'binary' in output element\n");
        }
      }
      else {
        throw xmdsException(thisElement, "Output format should be either 'ascii' or 'binary'");
      }
    }

    // now grab the precision attribute if available
    XMLString outputPrecision = "double";  // set it to the default while declaring it
    getAttributeStrings(yourElement, "output", "precision", NOT_REQD, outputPrecision);
    if (outputPrecision == "") {
      outputPrecision = "double";  // the default option
    }
    else {
      if (outputPrecision == "double") {
        myParameters.useDouble = true;
        if (verbose()) {
          printf("Output precision set to 'double' in output element\n");
        }
      }
      else if (outputPrecision == "single") {
        myParameters.useDouble = false;
        if (verbose()) {
          printf("Output precision set to 'single' in output element\n");
        }
      }
      else {
        throw xmdsException(thisElement, "Output precision should be either 'double' or 'single'");
      }
    }

    // now process the output element
    myOutput = createxmdsOutput();
    const Element* yourElement = dynamic_cast<const Element*> (candidateElements->item(0));
    myOutput->processElement(yourElement);
  }
  else
    throw xmdsException(yourElement,"Cannot find <output> element");

  // ************************************
  // find sequence element and process

  candidateElements = yourElement->getElementsByTagName("sequence",0);

  if(candidateElements->length()>1) {
    throw xmdsException(yourElement,"Multiple <sequence> elements at top level");
  }
  else if(candidateElements->length()==1) {
    mySequence = createxmdsSequence();
    const Element* yourElement = dynamic_cast<const Element*> (candidateElements->item(0));
    mySequence->processElement(yourElement);
  }
  else
    throw xmdsException(yourElement,"Cannot find <sequence> element");

  myField->outputSampleCount();
  
  mySequence->outputSampleCount();

  myOutput->finaliseGeometry();

  // ************************************
  // Allow non-stochastic mpi?

  if(myParameters.usempi) {
      if((myField->geometry()->nDims()<2)&(!myParameters.stochastic)) {
        myParameters.usempi=0;
        printf("MPI disabled.  MPI cannot be used for deterministic problems with transverse dimensions < 2\n");
      }
      if((myField->geometry()->nDims()>1)&(!myParameters.stochastic)){
        printf("MPI used for a deterministic problem. Ain't life great? \n");
      }
    }
  
};

// ******************************************************************************
const xmdsSimulation::simulationParametersStruct*
xmdsSimulation::parameters() const {
  return &myParameters;
}

// ******************************************************************************
xmdsField* xmdsSimulation::field() const {
  return myField;
};

// ******************************************************************************
xmdsSimulation::argvStruct* 
xmdsSimulation::argStruct() const {
  return &myArgStruct;
};

// ******************************************************************************
const xmdsOutput* xmdsSimulation::output() const {
  return myOutput;
};

// ******************************************************************************
const xmdsSequence* xmdsSimulation::sequence() const {
  return mySequence;
};

// ******************************************************************************
unsigned long xmdsSimulation::nextSegmentNumber() const {
  if(debugFlag) {
    printf("xmdsSimulation::nextSegmentNumber\n");
  }

  myCurrentSegmentNumber++;

  return myCurrentSegmentNumber-1;
};

// ******************************************************************************
void xmdsSimulation::makeCode(
                              const unsigned long& inFileSplitPoint) const {
  if(debugFlag) {
    printf("xmdsSimulation::makeCode\n");
  }

  if(myOutput != 0) {
    myOutput->setInFileSplitPoint(inFileSplitPoint);
  }

  XMLString fileName = myParameters.simulationName;
  fileName += ".cc";

  FILE* outfile = fopen(fileName.c_str(),"w");

  if(outfile==0) {
    sprintf(errorMessage(),"Unable to open file '%s' for write access",fileName.c_str());
    throw xmdsException(errorMessage());  
  }
  if(verbose()) {
    printf("Beginning to write code ...\n");
  }

  try {
    writeIncludes(outfile);
    writeDefines(outfile);
    writeGlobals(outfile);
    writePrototypes(outfile);
    writeRoutines(outfile);
  }
  catch (xmdsException xmdsExceptionErr) {
    printf("Error: simulation failed to write output code\n");
    printf("due to the following xmdsException:\n");
    printf("%s",xmdsExceptionErr.getError());
    fclose(outfile);
    throw xmdsException("Internal Error"); 
  }
  fclose(outfile);
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsSimulation private
// ******************************************************************************
// ******************************************************************************

// ******************************************************************************
void xmdsSimulation::writeIncludes(
                                   FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsSimulation::writeIncludes\n");
  }

  if(verbose()) {
    printf("Writing simulation includes ...\n");
  }

  fprintf(outfile,
          "// ********************************************************\n"
  "// simulation includes\n"
  "\n"
  "#include <stdio.h>\n"
  "#include <math.h>\n"
  "#include <string>\n"
  "#include <fstream>\n"
  "#include <iostream>\n"
  "#include <sstream>\n"
  "#include <cstdlib>\n"
  "#include <unistd.h>\n"
  "#include <sys/time.h>\n"
  "#include <time.h>\n");
  if(simulation()->argStruct()->nameList.size() != 0) {
    fprintf(outfile,"#include <getopt_xmds.h>\n");
  }
  fprintf(outfile,"#include <vector>\n");
  if(myParameters.usempi) {
    fprintf(outfile,"#include <mpi.h>\n");
  }
  if(myParameters.nThreads>1) {
    fprintf(outfile,"#include <fftw_threads.h>\n");
  }
  else if((!myParameters.stochastic)&(myParameters.usempi)) {
      fprintf(outfile,
              "#include <fftw_mpi.h>\n"
              "#include <fftw.h>\n");
    }
  else
    fprintf(outfile,"#include <fftw.h>\n");
  fprintf(outfile,"#include <xmdscomplex.h>\n");
  if (myParameters.binaryOutput) {
    fprintf(outfile,"#include <xmdsconfig.h>\n");
  }
  fprintf(outfile,
          "using namespace std;\n"
          "\n");
};

// ******************************************************************************
void xmdsSimulation::writeDefines(
                                  FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsSimulation::writeDefines\n");
  }

  if(verbose()) {
    printf("Writing simulation defines ...\n");
  }

  if(myParameters.mpiMethod=="Scheduling") {
  fprintf(outfile,
          "// ********************************************************\n"
		  "//Adaptive Scheduler defines\n"
		  "\n"
		  "#define MIN(a,b) (((a)<(b))?(a):(b))\n"
		  "#define MAX(a,b) (((a)<(b))?(b):(a))\n"
		  "#define DEBUG_MPI_Scheduling 1\n\n");
  }

  fprintf(outfile,
          "// ********************************************************\n"
          "// simulation defines\n"
          "\n");
  if(myParameters.nThreads>1) {
    fprintf(outfile,"#define _num_threads %li\n\n",myParameters.nThreads);
  }

  if(myParameters.stochastic) {
    fprintf(outfile,"#define _n_noises %li\n",myParameters.nNoises);
    for(long i=0;i<myParameters.nNoises;i++) {
      fprintf(outfile,"#define n_%li _noises[%li]\n",i+1,i);
    }
    if(myParameters.nPaths>1) {
      fprintf(outfile,"#define _n_paths %li\n",myParameters.nPaths);
    }
  }
  fprintf(outfile,"#define d%s _step\n",myParameters.propDimName.c_str());
  fprintf(outfile,"\n");

  xmdsElement::writeDefines(outfile);
  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsSimulation::writeGlobals(
                                  FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsSimulation::writeGlobals\n");
  }

  if(verbose()) {
    printf("Writing simulation globals ...\n");
  }

	if(myParameters.mpiMethod=="Scheduling") {
		fprintf(outfile,
				"// ********************************************************\n"
				"//Adaptive Scheduler globals\n"
				"\n"
				"pthread_mutex_t tasklock;    /*Ensures mutual exclusion when assigning tasks*/\n"
				"pthread_mutex_t finlock;     /*Lock to synchronize completion of thread and master*/\n"
				"\n"
				"int total_paths;             /*total number of paths of each type to be completed*/\n"
				"int *slave_stat;             /*Array of slaves indicating their present status IDLE(0) or BUSY(1)*/\n"
				"\n"
				"int *fschedcount;            /*Array of num full paths scheduled for each slave*/\n"
				"double *fcommave;            /*Array of full path communication latency for each slave*/\n"
				"double *fthroughput;         /*Array of full path throughput (task/sec) for each slave*/\n"
				"int fullpaths_assigned=0;    /*number of full paths assigned or completed*/\n"
				"static int fpathnum = 0;     /*for debugging only*/\n"
				"\n");
		if(myParameters.errorCheck) {
			fprintf(outfile,
					"int *hschedcount;            /*Array of num half paths scheduled for each slave*/\n"
					"double *hcommave;            /*Array of half path communication latency for each slave*/\n"
					"double *hthroughput;         /*Array of half path throughput (task/sec) for each slave*/\n"
					"int halfpaths_assigned=0;    /*number of half paths assigned or completed*/\n"
					"static int hpathnum = 0;     /*for debugging only*/\n\n");
		}
	}

  fprintf(outfile,
          "// ********************************************************\n"
          "// simulation globals\n"
          "\n");
  fprintf(outfile,"double %s;\n",myParameters.propDimName.c_str());
  fprintf(outfile,"\n");
  if(myParameters.errorCheck) {
    fprintf(outfile,"bool _half_step;\n\n");
  }
  if(myParameters.usempi){
    fprintf(outfile,"int rank, size;\n\n");
  }

  if(myParameters.stochastic) {
    if (myParameters.noiseKind == "poissonian") {
      fprintf(outfile, 
	      "double _noiseMean = %g;\n",myParameters.noiseMean);
    }
    if(myParameters.errorCheck) {
      fprintf(outfile,
              "unsigned short _gen1[3];\n"
              "unsigned short _gen2[3];\n");
    }
    else
      fprintf(outfile,"unsigned short _gen[3];\n");
    fprintf(outfile,"\n");
  }
  else if(myParameters.usempi){
    fprintf(outfile,
            "int local_nx, local_x_start, local_ny_after_transpose, local_y_start_after_transpose, total_local_size;\n\n");
  }

  xmdsElement::writeGlobals(outfile);
  fprintf(outfile,"\n");
};

// ******************************************************************************
void xmdsSimulation::writePrototypes(
                                     FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsSimulation::writePrototypes\n");
  }

  if(verbose()) {
    printf("Writing simulation prototypes ...\n");
  }

  if(myParameters.stochastic) {
  
	  if(myParameters.mpiMethod=="Scheduling") {
		  fprintf(outfile,
				  "// ********************************************************\n"
				  "//Adaptive Scheduler prototypes\n"
				  "\n"
				  "int print_array2(int *array, int size);           /*prints int arrays*/\n"
				  "int print_array(double *array, int size);         /*prints double arrays*/\n"
				  "void *slave(void *ptr);                           /*slave routine*/\n"
				  "int master();									  /*master routine*/\n"
				  "void *mslave(void *ptr);                          /*thread slave*/\n\n"
				  );
	  }

    if (debugFlag) {
	printf("xmdsSimulation::writePrototypes noiseKind = %s\n",myParameters.noiseKind.c_str()); 
    }
    if (myParameters.noiseKind == "gaussian") {
      fprintf(outfile,
	      "// ********************************************************\n"
	      "// simulation prototypes\n"
	      "\n"
	      "void _make_noises(\n"
	      "   unsigned short *const _generator,\n"
	      "   const double& _var,\n"
	      "   double *const _noise_vector,\n"
	      "   const unsigned long& _n);\n"
	      "\n");
    }
    else if (myParameters.noiseKind == "gaussFast") {
      fprintf(outfile,
	      "// ********************************************************\n"
	      "// simulation prototypes\n"
	      "\n"
	      "void _make_noises(\n"
	      "   unsigned short *const _generator,\n"
	      "   const double& _var,\n"
	      "   double *const _noise_vector,\n"
	      "   const unsigned long& _n);\n"
	      "\n");
    }
    else if (myParameters.noiseKind == "uniform") {
      fprintf(outfile,
	      "// ********************************************************\n"
	      "// simulation prototypes\n"
	      "\n"
	      "void _make_noises(\n"
	      "   unsigned short *const _generator,\n"
	      "   const double& _var,\n"
	      "   double *const _noise_vector,\n"
	      "   const unsigned long& _n);\n"
	      "\n");
    }
    else if (myParameters.noiseKind == "poissonian") {
      fprintf(outfile,
	      "// ********************************************************\n"
	      "// simulation prototypes\n"
	      "\n"
	      "void _make_noises(\n"
	      "   unsigned short *const _generator,\n"
	      "   const double& _var,\n"
	      "   double *const _noise_vector,\n"
	      "   const unsigned long& _n);\n"
	      "\n"
	      "double gammln(double _xx);\n"
	      "\n"
	      "double poidev(double xm, unsigned short *const _generator);\n"
	      "\n");
    }
    else {
      sprintf(errorMessage(),"Incorrectly specified noise kind (%s) in xmdsSimulation::writePrototypes",myParameters.noiseKind.c_str());
      throw xmdsException(errorMessage());
    }
  }

  xmdsElement::writePrototypes(outfile);
};

// ******************************************************************************
void xmdsSimulation::writeRoutines(
                                   FILE *const outfile) const {
  if(debugFlag) {
    printf("xmdsSimulation::writeRoutines\n");
  }

  if(verbose()) {
    printf("Writing simulation routines ...\n");
  }

  // NOTE:  it turns out the erand48() is more likely to be thread safe, and generate
  // decent random numbers.  rand() isn't necessarily thread safe, and the thread safe
  // version (rand_r()) isn't as good a random number generator.
  if(myParameters.stochastic) {
    if (debugFlag) {
	printf("xmdsSimulation::writeRoutines noiseKind = %s\n",myParameters.noiseKind.c_str());
    }
    if (myParameters.noiseKind == "gaussian") {
      fprintf(outfile,
              "// ********************************************************\n"
              "//               simulation routines\n"
              "// ********************************************************\n"
              "\n"
              "// *************************\n"
              "void _make_noises(\n"
              " unsigned short *const _generator,\n"
              " const double& _var,\n"
              " double *const _noise_vector,\n"
              " const unsigned long& _n) {\n"
              "\n"
              "for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
              " const double _mag = sqrt(-2*log(erand48(_generator))*_var);\n"
              " const double _arg = 2*M_PI*erand48(_generator);\n"
              " _noise_vector[_i0 + 0] = _mag*cos(_arg);\n"
              " _noise_vector[_i0 + 1] = _mag*sin(_arg);\n"
              " }\n"
              "}\n"
              "\n");
    }
    else if (myParameters.noiseKind == "gaussFast") {
      fprintf(outfile,
              "// ********************************************************\n"
              "//               simulation routines\n"
              "// ********************************************************\n"
              "\n"
              "// *************************\n"
              "void _make_noises(\n"
              " unsigned short *const _generator,\n"
              " const double& _var,\n"
              " double *const _noise_vector,\n"
              " const unsigned long& _n) {\n"
              "\n"
              "for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
              "  double _v1, _v2, _rsq;\n"
              "  do {\n"
              "    _v1 = 2.0*erand48(_generator)-1.0;\n"
              "    _v2 = 2.0*erand48(_generator)-1.0;\n"
              "    _rsq = _v1*_v1 + _v2*_v2;\n"
              "  } while(_rsq >= 1.0 || _rsq == 0.0);\n"
              "  const double _fac = sqrt(-2.0*_var*log(_rsq)/_rsq);\n"
              "  _noise_vector[_i0 + 0] = _v1*_fac;\n"
              "  _noise_vector[_i0 + 1] = _v2*_fac;\n"
              " }\n"
              "}\n"
              "\n");
    }
    else if (myParameters.noiseKind == "uniform") {
      fprintf(outfile,
              "// ********************************************************\n"
              "//               simulation routines\n"
              "// ********************************************************\n"
              "\n"
              "// *************************\n"
              "void _make_noises(\n"
              " unsigned short *const _generator,\n"
              " const double& _var,\n"
              " double *const _noise_vector,\n"
              " const unsigned long& _n) {\n"
              "\n"
              "for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
              "  const double _fac = sqrt(_var);\n"
              "  _noise_vector[_i0 + 0] = _fac*erand48(_generator);\n"
              "  _noise_vector[_i0 + 1] = _fac*erand48(_generator);\n"
              " }\n"
              "}\n"
              "\n");
    }
    else if (myParameters.noiseKind == "poissonian") {
      fprintf(outfile,
              "// ********************************************************\n"
              "//               simulation routines\n"
              "// ********************************************************\n"
              "\n"
              "// *************************\n"
              "double gammln(double _xx) {\n"
              "  double _x, _y, _tmp, _ser;\n"
              "  const double _cof[6] = {76.18009172947146,-86.50532032941677,\n"
              "                          24.01409824083091,-1.231739572450155,\n"
              "                          0.1208650973866179e-2,-0.5395239384953e-5};\n"
              "  int _j;\n"
              "  _y = _x = _xx;\n"
              "  _tmp = _x + 5.5;\n"
              "  _tmp -= (_x+0.5)*log(_tmp);\n"
              "  _ser = 1.000000000190015;\n"
              "  for (_j=0; _j<=5; _j++) _ser += _cof[_j]/++_y;\n"
              "  return -_tmp+log(2.5066282746310005*_ser/_x);\n"
              "}\n"
              "\n"
              "double poidev(double xm, unsigned short *const _generator) {\n"
              "  static double sq, alxm, g;\n"
              "  double em,t,y;\n"
              "  if (xm < 12.0) {        // Use direct method \n"
              "    g=exp(-xm);\n"
              "    em = -1;\n"
              "    t=1.0;\n"
              "    // Instead of adding exponential deviates it is equivalent\n"
              "    // to multiply uniform deviates.  We never actually have to\n"
              "    // take the log, merely compare to the pre-computed exponential\n"
              "    do {\n"
              "      ++em;\n"
              "      t *= erand48(_generator);\n"
              "    } while (t > g);\n"
              "  } else {                  // Use rejection method\n"
              "    sq=sqrt(2.0*xm);\n"
              "    alxm=log(xm);\n"
              "    g=xm*alxm-gammln(xm+1.0);\n"
              "    do {\n"
              "      do { // y is a deviate from a Lorenzian comparison function\n"
              "        y=tan(M_PI*erand48(_generator));\n"
              "        em=sq*y+xm;        // em is y, shifted and scaled\n"
              "      } while (em < 0.0);  // Reject if in regime of zero probability\n"
              "      em=floor(em);        // The trick for integer-valued distributions\n"
              "      t=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);\n"
              "      // The ratio of the desired distribution to the comparison\n"
              "      // function; we reject by comparing it to another uniform\n"
              "      // deviate. The factor 0.9 so that t never exceeds 1.\n"
              "    } while (erand48(_generator) > t);\n"
              "  }\n"
              "  return em;\n"
              "  }\n"
              "\n"
              "void _make_noises(\n"
              " unsigned short *const _generator,\n"
              " const double& _var,\n"
              " double *const _noise_vector,\n"
              " const unsigned long& _n) {\n"
              "\n"
              "for(unsigned long _i0=0; _i0<_n; _i0+=2) {\n"
              "  const double _fac = sqrt(_var);\n"
              "  _noise_vector[_i0 + 0] = _fac*poidev(_noiseMean, _generator);\n"
              "  _noise_vector[_i0 + 1] = _fac*poidev(_noiseMean, _generator);\n"
              " }\n"
              "}\n"
              "\n");
    }
    else {
      sprintf(errorMessage(),"Incorrectly specified noise kind (%s) in xmdsSimulation::writeRoutines",myParameters.noiseKind.c_str());
      throw xmdsException(errorMessage());
    }
  }

  xmdsElement::writeRoutines(outfile);
};

// ******************************************************************************
xmdsGlobals* xmdsSimulation::createxmdsGlobals() {
  if(debugFlag) {
    printf("xmdsSimulation::createxmdsGlobals\n");
  }

  xmdsGlobals* newxmdsGlobals =
    new xmdsGlobals(this,verbose());
  addChild((xmdsElement*) newxmdsGlobals);
  return newxmdsGlobals;
};

// ******************************************************************************
xmdsField* xmdsSimulation::createxmdsField() {
  if(debugFlag) {
    printf("xmdsSimulation::createxmdsField\n");
  }

  xmdsField* newxmdsField =
    new xmdsField(this,verbose());
  addChild((xmdsElement*) newxmdsField);
  return newxmdsField;
};

// ******************************************************************************
xmdsArgv* xmdsSimulation::createxmdsArgv() {
  if(debugFlag) {
    printf("xmdsSimulation::createxmdsArgv\n");
  }

  xmdsArgv* newxmdsArgv = new xmdsArgv(this,verbose());
  addChild((xmdsElement*) newxmdsArgv);
  return newxmdsArgv;
};

// ******************************************************************************
xmdsOutput* xmdsSimulation::createxmdsOutput() {
  if(debugFlag) {
    printf("xmdsSimulation::createxmdsOutput\n");
  }

  xmdsOutput* newxmdsOutput =
    new xmdsOutput(this,verbose());
  addChild((xmdsElement*) newxmdsOutput);
  return newxmdsOutput;
};

// ******************************************************************************
xmdsSequence* xmdsSimulation::createxmdsSequence() {
  if(debugFlag) {
    printf("xmdsSimulation::createxmdsSequence\n");
  }

  xmdsSequence* newxmdsSequence =
    new xmdsSequence(this,verbose());
  addChild((xmdsElement*) newxmdsSequence);
  return newxmdsSequence;
};