/*
 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: xmdssequence.cc,v 1.43 2005/10/27 00:58:39 joehope Exp $
*/

/*! @file xmdssequence.cc
  @brief Sequence element parsing classes and methods

  More detailed explanation...
*/

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

// ******************************************************************************
// *****************************************************************************
//                             xmdsSequence public
// *****************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsSequences=0;  //!< Number of xmds sequence objects

// ******************************************************************************

xmdsSequence::xmdsSequence(
			   const xmdsSimulation *const yourSimulation,
			   const bool& yourVerboseMode) :
  xmdsSegment(yourSimulation,yourVerboseMode) {
  if(debugFlag) {
    nxmdsSequences++;
    printf("xmdsSequence::xmdsSequence\n");
    printf("nxmdsSequences=%li\n",nxmdsSequences);
  }
  nCycles=1;
};

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

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

  unsigned long i;
  const NodeList* myElements;
  const Element* nextElement;
  xmdsSegment* newSegment;

  if(verbose()) {
    printf("Processing sequence element ...\n");
  }

  // ************************************
  // process elements in order

  myElements = yourElement->getElementsByTagName("*",0);
  i=0;
  nextElement = dynamic_cast<const Element*> (myElements->item(0));
  if(*nextElement->nodeName() == "cycles") {
    if(segmentNumber==0) {
      throw xmdsException(nextElement,"Cannot have cycles in top level sequence");
    }
    if(!nextElement->textContent(0)->asULong(nCycles)) {
      throw xmdsException(nextElement,"invalid positive integer format");
    }
    i++;
  }
  else
    // default number of cycles
    nCycles=1;

  while(i<myElements->length()) {
    nextElement = dynamic_cast<const Element*> (myElements->item(i));
    if(*nextElement->nodeName() == "sequence") {
      newSegment = createxmdsSequence();
    }
    else if(*nextElement->nodeName() == "integrate") {

      // ************************************
      // find 'algorithm' so that we can create the appropriate
      // type of xmdsIntegrate class
      list<XMLString> tempXMLStringList;
      getAssignmentStrings(nextElement,"algorithm",0,1,tempXMLStringList);
      if(tempXMLStringList.size()==0) {
	// need default algorithm
	if(simulation()->parameters()->stochastic) {
	  printf("Algorithm defaulting to SIEX.\n");
	  tempXMLStringList.push_back("SIEX");
	}
	else {
	  printf("Algorithm defaulting to RK4EX.\n");
	  tempXMLStringList.push_back("RK4EX");
	}
      }
      if(*tempXMLStringList.begin()=="RK4EX") {
	newSegment = createxmdsIntegrateRK4EX();
      }
      else if(*tempXMLStringList.begin()=="RK4IP") {
	newSegment = createxmdsIntegrateRK4IP();
      }
      else if(*tempXMLStringList.begin()=="ARK45IP") {
	newSegment = createxmdsIntegrateARK45IP();
      }
      else if(*tempXMLStringList.begin()=="ARK45EX") {
	newSegment = createxmdsIntegrateARK45EX();
      }
      else if(*tempXMLStringList.begin()=="SIEX") {
	newSegment = createxmdsIntegrateSIEX();
      }
      else if(*tempXMLStringList.begin()=="SIIP") {
	newSegment = createxmdsIntegrateSIIP();
      }
      else {
	sprintf(errorMessage(),"Integrate algorithm '%s' unknown.\nCurrent algorithms are SIEX, SIIP, ARK45IP, ARK45EX, RK4EX, and RK4IP.",
		tempXMLStringList.begin()->c_str());
	throw xmdsException(yourElement,errorMessage());
      }
    }
    else if(*nextElement->nodeName() == "filter") {
      newSegment = createxmdsFilter();
    }
    else {
      sprintf(errorMessage(),"segment <%s> unknown",nextElement->nodeName()->c_str());
      throw xmdsException(yourElement,errorMessage());
    }
    newSegment->processElement(nextElement);
    i++;
  }
};

// ******************************************************************************
void xmdsSequence::outputSampleCount() const {
  if(debugFlag) {
    printf("xmdsSequence::outputSampleCount\n");
  }

  for(unsigned long i=0;i<nCycles;i++) {
    for(list<const xmdsSegment*>:: const_iterator ppxmdsSegment
	  = mySegmentList.begin(); ppxmdsSegment !=
	  mySegmentList.end(); ppxmdsSegment++) {
      (*ppxmdsSegment)->outputSampleCount();
    }
  }
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsSequence private
// ******************************************************************************
// ******************************************************************************

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

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

  if(segmentNumber==0) {
    fprintf(outfile,
	    "// ********************************************************\n"
	    "// segment defines\n"
	    "\n");
  }

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

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

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

  if(segmentNumber==0) {
    fprintf(outfile,
	    "// ********************************************************\n"
	    "// segment globals\n"
	    "\n");
  }

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

// ******************************************************************************
void xmdsSequence::writePrototypes(
								   FILE *const outfile) const {
	if(debugFlag) {
		printf("xmdsSequence::writePrototypes\n");
	}
	
	if(verbose()) {
		printf("Writing sequence prototypes ...\n");
	}
	
	if(segmentNumber==0) {
		fprintf(outfile,
				"// ********************************************************\n"
				"// segment prototypes\n"
				"\n");
		if(!simulation()->parameters()->usempi)
			fprintf(outfile,"void _segment0(); // top level sequence\n");
		else if(simulation()->parameters()->mpiMethod=="Scheduling")
			fprintf(outfile,"void _segment0(const int);                        /*top level sequence*/\n");
		else 
			fprintf(outfile,"void _segment0(const int& size); // top level sequence\n");
		
		fprintf(outfile,"\n");
	}
	else
		fprintf(outfile,"void _segment%li(unsigned long cycle); // sequence\n",segmentNumber);        
	
	xmdsElement::writePrototypes(outfile);
	fprintf(outfile,"\n");
};

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

  if(segmentNumber==0) {
    // I am the main routine

    if(verbose()) {
      printf("Writing main routine ...\n");
    }

    fprintf(outfile,
	    "// ********************************************************\n"
	    "//		main routine\n"
	    "// ********************************************************\n"
	    "\n"
	    "// *************************\n");
    if (verbose()) {
      printf("argStruct()->nameList.size() = %i\n",simulation()->argStruct()->nameList.size());
      printf("argStruct()->typeList.size() = %i\n",simulation()->argStruct()->typeList.size()); 
      printf("argStruct()->defaultValueList.size() = %i\n",simulation()->argStruct()->defaultValueList.size());
      printf("argStruct()->nameList.size() = %i\n",simulation()->argStruct()->shortOptionList.size());
      printf("argStruct()->typeConversionList.size() = %i\n",simulation()->argStruct()->typeConversionList.size());

      list<string>::const_iterator inameList = simulation()->argStruct()->nameList.begin();
      list<string>::const_iterator itypeList = simulation()->argStruct()->typeList.begin();
      list<string>::const_iterator idefaultValueList = simulation()->argStruct()->defaultValueList.begin();
      list<string>::const_iterator ishortOptionList = simulation()->argStruct()->shortOptionList.begin();
      list<string>::const_iterator itypeConversionList = simulation()->argStruct()->typeConversionList.begin();

      printf("argvStruct()->nameList = \n");
      for (unsigned long int i=0; i<simulation()->argStruct()->nameList.size(); i++) {
	printf("\t%s\n",inameList->c_str());
	inameList++;
      }
      printf("argvStruct()->typeList = \n");
      for (unsigned long int i=0; i<simulation()->argStruct()->typeList.size(); i++) {
	printf("\t%s\n",itypeList->c_str());
	itypeList++;
      }
      printf("argvStruct()->defaultValueList = \n");
      for (unsigned long int i=0; i<simulation()->argStruct()->defaultValueList.size(); i++) {
	printf("\t%s\n",idefaultValueList->c_str());
	idefaultValueList++;
      }
      printf("argvStruct()->shortOptionList = \n");
      for (unsigned long int i=0; i<simulation()->argStruct()->shortOptionList.size(); i++) {
 	printf("\t%s\n",ishortOptionList->c_str());
 	ishortOptionList++;
      }
      printf("argvStruct()->typeConversionList = \n");
      for (unsigned long int i=0; i<simulation()->argStruct()->typeConversionList.size(); i++) {
	printf("\t%s\n",itypeConversionList->c_str());
	itypeConversionList++;
      }
    }
    // if(simulation()->parameters()->usempi) {
    if(simulation()->argStruct()->nameList.size() == 0) {
      fprintf(outfile,"int main(int argc, char *argv[]) {\n"); 
    }
    else if (simulation()->argStruct()->nameList.size() != 0) {
      fprintf(outfile,
	      "int main(int argc, char **argv) {\n"
	      "int resp;\n"
	      "while (1) {\n"
	      "\tstatic struct option long_options[] = \n"
	      "\t\t{\n");
      list<string>::const_iterator inameList = simulation()->argStruct()->nameList.begin();
      //     list<string>::const_iterator itypeList = simulation()->argStruct()->typeList.begin();
      //     list<string>::const_iterator idefaultValueList = simulation()->argStruct()->defaultValueList.begin();
      list<string>::const_iterator ishortOptionList = simulation()->argStruct()->shortOptionList.begin();
      list<string>::const_iterator itypeConversionList = simulation()->argStruct()->typeConversionList.begin();

      fprintf(outfile,"\t\t\t{\"help\", no_argument, 0, 'h'},\n");
      for (unsigned long int i=0; i<simulation()->argStruct()->nameList.size(); i++) {
	fprintf(outfile,"\t\t\t{\"%s\", required_argument, 0,'%s'},\n",inameList->c_str(),ishortOptionList->c_str());
	inameList++;
	ishortOptionList++;
      }
      fprintf(outfile,
	      "\t\t\t{0, 0, 0, 0}\n"
	      "\t\t};\n"
	      "\tint option_index = 0;\n"
	      "\tresp = getopt_xmds_long(argc, argv, \"");
      ishortOptionList = simulation()->argStruct()->shortOptionList.begin();
      for (unsigned long int i=0; i<simulation()->argStruct()->shortOptionList.size(); i++) {
	fprintf(outfile,"%s:",ishortOptionList->c_str());
	ishortOptionList++;
      }
      // also need to output the short option version of "help" so that 
      // '-h' works
      fprintf(outfile, "h");
      fprintf(outfile,
	      "\", long_options, &option_index);\n"
	      "\tif (resp == -1) {\n"
	      "\t\tbreak;\n"
	      "\t}\n"
	      "\tswitch (resp) {\n"
	      "\t\tcase 'h':\n"
	      "\t\t\tprintf(\"Usage: %s",
	      simulation()->parameters()->simulationName.c_str());
      ishortOptionList = simulation()->argStruct()->shortOptionList.begin();
      list<string>::const_iterator itypeList = simulation()->argStruct()->typeList.begin();
      for (unsigned long int i=0; i<simulation()->argStruct()->nameList.size(); i++) {
	fprintf(outfile," -%s < %s>",ishortOptionList->c_str(),itypeList->c_str());
	ishortOptionList++;
	itypeList++;
      }
      fprintf(outfile,
	      "\\n\\n\");\n"
	      "\t\t\tprintf(\"Details:\\n\");\n"
	      "\t\t\tprintf(\"Option\\t\\tType\\t\\tDefault value\\n\");\n");
      inameList = simulation()->argStruct()->nameList.begin();
      ishortOptionList = simulation()->argStruct()->shortOptionList.begin();
      itypeList = simulation()->argStruct()->typeList.begin();
      list<string>::const_iterator idefaultValueList = simulation()->argStruct()->defaultValueList.begin();
      for (unsigned long int i=0; i<simulation()->argStruct()->nameList.size(); i++) {
	fprintf(outfile,"\t\t\tprintf(\"-%s, --%s\\t%s\\t\\t%s\\n\");\n",ishortOptionList->c_str(),inameList->c_str(),itypeList->c_str(),idefaultValueList->c_str());
	inameList++;
	ishortOptionList++;
	itypeList++;
	idefaultValueList++;
      }
      fprintf(outfile,"\t\t\texit(0);\n");

      inameList = simulation()->argStruct()->nameList.begin();
      ishortOptionList = simulation()->argStruct()->shortOptionList.begin();
      itypeConversionList = simulation()->argStruct()->typeConversionList.begin();
      for (unsigned long int i=0; i<simulation()->argStruct()->nameList.size(); i++) {
	fprintf(outfile,"\t\tcase '%s':\n",ishortOptionList->c_str());
	if (*itypeConversionList == "atof") {
	  fprintf(outfile,"\t\t\t%s = atof(optarg);\n",inameList->c_str());
	  fprintf(outfile,"\t\t\tbreak;\n");
	}
	else if (*itypeConversionList == "atoi") {
	  fprintf(outfile,"\t\t\t%s = atoi(optarg);\n",inameList->c_str());
	  fprintf(outfile,"\t\t\tbreak;\n");
	}
	else if (*itypeConversionList == "") {
	  fprintf(outfile,"\t\t\t%s = optarg;\n",inameList->c_str());
	  fprintf(outfile,"\t\t\tbreak;\n");
	}
	inameList++;
	ishortOptionList++;
	itypeConversionList++;
      }
      fprintf(outfile,
	      "\t\tdefault:\n"
	      "\t\t\texit(10);\n"
	      "\t\t}\n"
	      "\t}\n"
	      "\tif (optind_xmds < argc) {\n"
	      "\t\tprintf(\"Sorry, I didn't understand some (or all) of your arguments\\n\");\n"
	      "\t\t\tprintf(\"Usage: %s",simulation()->parameters()->simulationName.c_str());
      ishortOptionList = simulation()->argStruct()->shortOptionList.begin();
      itypeList = simulation()->argStruct()->typeList.begin();
      for (unsigned long int i=0; i<simulation()->argStruct()->nameList.size(); i++) {
	fprintf(outfile," -%s < %s>",ishortOptionList->c_str(),itypeList->c_str());
	ishortOptionList++;
	itypeList++;
      }
      fprintf(outfile,"\\n\\n\");\n"
	      "\t\t\tprintf(\"Details:\\n\");\n"
	      "\t\t\tprintf(\"Option\\t\\tType\\t\\tDefault value\\n\");\n");
      inameList = simulation()->argStruct()->nameList.begin();
      ishortOptionList = simulation()->argStruct()->shortOptionList.begin();
      itypeList = simulation()->argStruct()->typeList.begin();
      idefaultValueList = simulation()->argStruct()->defaultValueList.begin();
      for (unsigned long int i=0; i<simulation()->argStruct()->nameList.size(); i++) {
	fprintf(outfile,"\t\t\tprintf(\"-%s, --%s\\t%s\\t\\t%s\\n\");\n",ishortOptionList->c_str(),inameList->c_str(),itypeList->c_str(),idefaultValueList->c_str());
	inameList++;
	ishortOptionList++;
	itypeList++;
	idefaultValueList++;
      }
      fprintf(outfile,
	      "\t\t\texit(0);\n"
	      "\t}\n\n");
    }      
    else {
      fprintf(outfile,"int main() {\n");
    }
    
    fprintf(outfile,"\n");

    if (simulation()->parameters()->binaryOutput) {
      fprintf(outfile,
	      "#ifndef CPU_IS_BIG_ENDIAN\n"
	      "   #ifndef CPU_IS_LITTLE_ENDIAN\n"
	      "   cout << \"Byte ordering isn't defined for this system!\\n\";\n"
	      "   cout << \"Please define CPU_IS_BIG_ENDIAN and CPU_IS_LITTLE_ENDIAN in config.h\\n\";\n"
	      "   exit(255);\n"
	      "   #endif\n"
	      "#endif\n");
    }

    if(simulation()->parameters()->usempi) {
      fprintf(outfile,
	      "MPI_Init(&argc, &argv);\n"
	      "MPI_Comm_size(MPI_COMM_WORLD, &size);\n"
	      "MPI_Comm_rank(MPI_COMM_WORLD, &rank);\n"
	      "\n");
    }

    // if(simulation()->field()->needsFFTWPlans()) {
    //  It's possible to need FFTW plans in output but not the main field through post-propagation
    fprintf(outfile,"/* set up fftw plans */\n");
    if(simulation()->parameters()->nThreads>1) {
      fprintf(outfile,"int fftw_threads_init(void);\n");
        }
    fprintf(outfile,"int _fftw_lattice[64];\n");
    fprintf(outfile,"\n");
    simulation()->field()->writePlanCreationCalls(outfile,1,simulation()->parameters()->useWisdom);
    simulation()->output()->writePlanCreationCalls(outfile,0,simulation()->parameters()->useWisdom);
    //	}

    if(simulation()->parameters()->benchmark) {
		if(simulation()->parameters()->usempi){
			fprintf(outfile,"double startTime = MPI_Wtime();\n");
		}
		else{
			fprintf(outfile,"time_t startTime = time(NULL);\n");
		}
    }

    if((simulation()->parameters()->usempi)&!(simulation()->parameters()->stochastic)) {
		fprintf(outfile,"fftwnd_mpi_local_sizes(_main_forward_plan, &local_nx, &local_x_start, &local_ny_after_transpose, &local_y_start_after_transpose, &total_local_size);\n\n");
		list<XMLString> vectorNamesList;
		const xmdsVector* tempVector;
		const char *typeName;
		simulation()->field()->vectorNames(vectorNamesList);
		for(list<XMLString>::const_iterator pXMLString = vectorNamesList.begin(); pXMLString != vectorNamesList.end(); pXMLString++) {
			if(!simulation()->field()->getVector(*pXMLString,tempVector))  {
				throw xmdsException("Lost one of the vectors in xmdssequence::writeRoutines.  Bad!  Bad dog!\n\n");
			}
			if(tempVector->vectorType()==DOUBLE) {
				typeName="double";
			}
			else {
				typeName="complex";
			}
			fprintf(outfile,"_%s_%s = new %s[total_local_size*_%s_%s_ncomponents];\n",
					simulation()->field()->name()->c_str(),tempVector->name()->c_str(),typeName,
					simulation()->field()->name()->c_str(),tempVector->name()->c_str());
			fprintf(outfile,"_%s_%s_work = new %s[total_local_size*_%s_%s_ncomponents];\n",
					simulation()->field()->name()->c_str(),tempVector->name()->c_str(),typeName,
					simulation()->field()->name()->c_str(),tempVector->name()->c_str());
		}
		fprintf(outfile,"\n");
		
		for(unsigned long i=0; i<simulation()->output()->nMomentGroups(); i++) {
			if(!simulation()->output()->momentGroup(i)->needscomplexRawVector()) {
				typeName="double";
			}
			else {
				typeName="complex";
			}
			fprintf(outfile,"_mg%li_raw = new %s[_mg%li_size*_mg%li_raw_ncomponents];\n",i,typeName,i,i);
			fprintf(outfile,"_mg%li_fullstep = new double[_mg%li_size*_mg%li_fullstep_ncomponents];\n",i,i,i);
			if(simulation()->parameters()->errorCheck) {
				fprintf(outfile,"_mg%li_halfstep = new double[_mg%li_size*_mg%li_halfstep_ncomponents];\n",i,i,i);
			}
		}
		fprintf(outfile,"\n");
		
    }

    simulation()->field()->assignActiveVectorPointers(outfile,"");
    simulation()->output()->assignActiveVectorPointers(outfile,"");

    if(simulation()->parameters()->stochastic) {
      for(unsigned long i=0;i<simulation()->output()->nMomentGroups();i++) {
	fprintf(outfile,"_mg%li_fullstep_initialise();\n",i);
	fprintf(outfile,"_mg%li_fullstep_sd_initialise();\n",i);
	if(simulation()->parameters()->errorCheck) {
	  fprintf(outfile,"_mg%li_halfstep_initialise();\n",i);
	  fprintf(outfile,"_mg%li_halfstep_sd_initialise();\n",i);
	}
	fprintf(outfile,"\n");
      }
    }

	if((simulation()->parameters()->mpiMethod=="Scheduling")&&(simulation()->parameters()->usempi)&&(simulation()->parameters()->stochastic)) {
		if(simulation()->parameters()->errorCheck) {
			if(simulation()->parameters()->nPaths>1) {
				fprintf(outfile,"printf(\"Rank[%%i] beginning full step paths\\n\",rank);\n");
			}
			else {
				fprintf(outfile,"printf(\"Beginning full step integration ...\\n\");\n");
			}
			fprintf(outfile,"\n");
			fprintf(outfile,"_half_step=0;\n");
			fprintf(outfile,"\n");
		}
		else {
			if(simulation()->parameters()->nPaths>1) {
				fprintf(outfile,"printf(\"Rank[%%i] beginning paths\\n\",rank);\n");
			}
			else {
				fprintf(outfile,"printf(\"Beginning integration ...\\n\");\n");
			}
			fprintf(outfile,"\n");
		}
		if(simulation()->parameters()->errorCheck) {
			fprintf(outfile,"_gen1[0] = %li+rank;\n",simulation()->parameters()->seed[0]);
			fprintf(outfile,"_gen2[0] = %li+rank;\n",simulation()->parameters()->seed[1]);
			fprintf(outfile,"_gen1[1] = 0;\n");
			fprintf(outfile,"_gen2[1] = 0;\n");
			fprintf(outfile,"_gen1[2] = 0;\n");
			fprintf(outfile,"_gen2[2] = 0;\n");
			fprintf(outfile,"erand48(_gen1);\n");
			fprintf(outfile,"erand48(_gen2);\n");
		}
		else {
			fprintf(outfile,"_gen[0] = %li+rank;\n",simulation()->parameters()->seed[0]);
			fprintf(outfile,"_gen[1] = 0;\n");
			fprintf(outfile,"_gen[2] = 0;\n");
			fprintf(outfile,"erand48(_gen);\n");
		}
		
		if(simulation()->parameters()->nPaths>1) {
			fprintf(outfile,"\ntotal_paths = _n_paths;\n\n");
		}
		else {
			fprintf(outfile,"\ntotal_paths = 1;\n\n");
		}
		
		fprintf(outfile,"//Rank 0 becomes master, other processes become slaves\n"
				"if (size > 1){\n"
				"	if (rank==0)\n"
				"		master();\n"
				"	else {\n"
				"		slave((void *) NULL);\n"
				"	}\n"
				"}\n"
				"//Run the simulation using the following if only one processor used\n"
				"else {\n");
		if(simulation()->parameters()->errorCheck)
			fprintf(outfile,"	_half_step = 0;\n");
		fprintf(outfile,
				"	_segment0(total_paths);\n\n");
		
		if(simulation()->parameters()->errorCheck) {
			fprintf(outfile,"   _gen1[0] = %li+rank;\n",simulation()->parameters()->seed[0]);
			fprintf(outfile,"   _gen2[0] = %li+rank;\n",simulation()->parameters()->seed[1]);
			fprintf(outfile,"   _gen1[1] = 0;\n");
			fprintf(outfile,"   _gen2[1] = 0;\n");
			fprintf(outfile,"   _gen1[2] = 0;\n");
			fprintf(outfile,"   _gen2[2] = 0;\n");
			fprintf(outfile,"   erand48(_gen1);\n");
			fprintf(outfile,"   erand48(_gen2);\n\n"
		
					"   _half_step = 1;\n");
			if(simulation()->parameters()->nPaths>1) {
				fprintf(outfile,"    printf(\"Rank[%%i] beginning half step paths\\n\",rank);\n");
			}
			else {
				fprintf(outfile,"    printf(\"Beginning half step integration ...\\n\");\n");
			}
			fprintf(outfile,"   _segment0(total_paths);\n\n");
		}
		fprintf(outfile,"goto write_output;\n"
				"}\n\n");
		
	}
	else {
		if(simulation()->parameters()->errorCheck) {
			if(simulation()->parameters()->nPaths>1) {
				if(simulation()->parameters()->usempi) {
					fprintf(outfile,"printf(\"Rank[%%i] beginning full step paths\\n\",rank);\n");
				}
				else {
					fprintf(outfile,"printf(\"Beginning full step paths\\n\");\n");
				}
			}
			else {
				if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
					fprintf(outfile,"printf(\"Rank[%%i] beginning full step integration ...\\n\",rank);\n");
				}
				else {
					fprintf(outfile,"printf(\"Beginning full step integration ...\\n\");\n");
				}
			}
			fprintf(outfile,"\n");
			fprintf(outfile,"_half_step=0;\n");
			fprintf(outfile,"\n");
		}
		else {
			if(simulation()->parameters()->nPaths>1) {
				if(simulation()->parameters()->usempi) {
					fprintf(outfile,"printf(\"Rank[%%i] beginning paths\\n\",rank);\n");
				}
				else {
					fprintf(outfile,"printf(\"Beginning paths\\n\");\n");
				}
			}
			else {
				if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
					fprintf(outfile,"printf(\"Rank[%%i] beginning integration ...\\n\",rank);\n");
				}
				else {
					fprintf(outfile,"printf(\"Beginning integration ...\\n\");\n");
				}
			}
			fprintf(outfile,"\n");
		}
		
		if(simulation()->parameters()->usempi) {
			fprintf(outfile,"_segment0(size);\n");
		}
		else {
			fprintf(outfile,"_segment0();\n");
		}
		fprintf(outfile,"\n");
		
		if(simulation()->parameters()->errorCheck) {
			if(simulation()->parameters()->nPaths>1) {
				if(simulation()->parameters()->usempi) {
					fprintf(outfile,"printf(\"Rank[%%i] beginning half step paths\\n\",rank);\n");
				}
				else {
					fprintf(outfile,"printf(\"Beginning half step paths\\n\");\n");
				}
			}
			else {
				if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
					fprintf(outfile,"printf(\"Rank[%%i] beginning half step integration ...\\n\",rank);\n");
				}
				else {
					fprintf(outfile,"printf(\"Beginning half step integration ...\\n\");\n");
				}
			}
			fprintf(outfile,"\n");
			fprintf(outfile,"_half_step=1;\n");
			fprintf(outfile,"\n");
			
			if(simulation()->parameters()->usempi) {
				fprintf(outfile,"_segment0(size);\n");
			}
			else {
				fprintf(outfile,"_segment0();\n");
			}
			fprintf(outfile,"\n");
		}
	}
	
    if(simulation()->parameters()->usempi&simulation()->parameters()->stochastic) {
		
		fprintf(outfile,"double* _temp_vector;\n");
		fprintf(outfile,"unsigned long _length;\n");
		fprintf(outfile,"\n");
		
		for(unsigned long i=0;i<simulation()->output()->nMomentGroups();i++) {
			fprintf(outfile,"_length = _mg%li_size*_mg%li_fullstep_ncomponents;\n",i,i);
			fprintf(outfile,"_temp_vector = new double[_length];\n");
			fprintf(outfile,"\n");
			fprintf(outfile,"MPI_Reduce(_mg%li_fullstep,_temp_vector,_length,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n",i);
			fprintf(outfile,"\n");
			fprintf(outfile,"for(unsigned long _i0=0; _i0<_length; _i0++)\n");
			fprintf(outfile,"	_mg%li_fullstep[_i0]=_temp_vector[_i0];\n",i);
			fprintf(outfile,"\n");
			if(simulation()->parameters()->nPaths>1) {
				fprintf(outfile,"MPI_Reduce(_mg%li_fullstep_sd,_temp_vector,_length,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n",i);
				fprintf(outfile,"\n");
				fprintf(outfile,"for(unsigned long _i0=0; _i0<_length; _i0++)\n");
				fprintf(outfile,"	_mg%li_fullstep_sd[_i0]=_temp_vector[_i0];\n",i);
				fprintf(outfile,"\n");
			}
			
			if(simulation()->parameters()->errorCheck) {
				fprintf(outfile,"MPI_Reduce(_mg%li_halfstep,_temp_vector,_length,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n",i);
				fprintf(outfile,"\n");
				fprintf(outfile,"for(unsigned long _i0=0; _i0<_length; _i0++)\n");
				fprintf(outfile,"	_mg%li_halfstep[_i0]=_temp_vector[_i0];\n",i);
				fprintf(outfile,"\n");
				if(simulation()->parameters()->nPaths>1) {
					fprintf(outfile,"MPI_Reduce(_mg%li_halfstep_sd,_temp_vector,_length,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n",i);
					fprintf(outfile,"\n");
					fprintf(outfile,"for(unsigned long _i0=0; _i0<_length; _i0++)\n");
					fprintf(outfile,"	_mg%li_halfstep_sd[_i0]=_temp_vector[_i0];\n",i);
					fprintf(outfile,"\n");
				}
			}
			
			fprintf(outfile,"delete[] _temp_vector;\n");
			fprintf(outfile,"\n");
		}
		if(simulation()->parameters()->mpiMethod=="Scheduling") {
			fprintf(outfile,"	/* Adaptive Scheduler single processor break point */\n"
					"write_output:\n\n");
		}
		
		fprintf(outfile,"if(rank==0)\n");
		fprintf(outfile,"	_write_output();\n");
		fprintf(outfile,"\n");
    }
    else if(simulation()->parameters()->usempi){
		fprintf(outfile,"if(rank==0)\n");
		fprintf(outfile,"	_write_output();\n");
		fprintf(outfile,"\n");
    }
    else {
      fprintf(outfile,"_write_output();\n");
      fprintf(outfile,"\n");
    }

    if(simulation()->parameters()->benchmark) {
		if(simulation()->parameters()->usempi){
			fprintf(outfile,"double endTime = MPI_Wtime();\n");
			fprintf(outfile,"if(rank==0)\n   ");
		}
		else{
			fprintf(outfile,"time_t endTime = time(NULL);\n");
		}
		fprintf(outfile,"printf(\"Time elapsed for simulation is: %%5.5f seconds\\n\",endTime-startTime);\n");
    }

    simulation()->field()->writePlanDeletionCalls(outfile);
    simulation()->output()->writePlanDeletionCalls(outfile);

    if(simulation()->parameters()->usempi) {
      fprintf(outfile,"MPI_Finalize();\n");
      fprintf(outfile,"\n");
    }

    if (simulation()->parameters()->bing) {
      // this should be the last statement in the main() routine
      fprintf(outfile,"printf(\"\\a\");\n");
    }
    
    fprintf(outfile,"return 0 ;\n");
    fprintf(outfile,"}\n");
    fprintf(outfile,"\n");

	if((simulation()->parameters()->mpiMethod=="Scheduling")&&(simulation()->parameters()->usempi)&&(simulation()->parameters()->stochastic)) {
		fprintf(outfile,
				"int master(){                                                                                                             \n"
				"                                                                                                                          \n"
				"  int outstanding = 0;       /*number of slaves that are still processing tasks*/                                         \n"
				"  int schedule[size];        /*Batch size scheduled for each slave [resetted every iteration]*/                           \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  int sendtag;               /*MPI Tag value, fullstep -> 1, halfstep -> 2*/                                              \n");
		else
			fprintf(outfile,
					"  int sendtag = 1;                                                          \n");
		fprintf(outfile,
				"                                                                                                                          \n"
				"  double timing[size];       /*Timing function to determine computation to communication ratio*/                          \n"
				"  int partitions[size];      /*Batch size scheduled for each slave [resetted after completion]*/                          \n"
				"                                                                                                                          \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  double *active_throughput; /*Active pointers, used for clarity*/                                                        \n"
					"  double *active_commave;                                                                                                 \n"
					"  int *active_assigned;                                                                                                   \n"
					"                                                                                                                          \n");
		fprintf(outfile,	
				"  int i, j, k;               /*indexes, and size of MPI_COMM_WORLD*/                                                      \n"
				"  double bufs[size];         /*MPI Input buffer*/                                                                         \n"
				"  int idle = size-1;         /*Number of idle processors not including master (size - outstanding -1)*/                   \n"
				"                                                                                                                          \n"
				"  MPI_Status stats[size];    /*MPI Structures*/                                                                           \n"
				"  MPI_Request reqs[size];                                                                                                 \n"
				"  int indices[size];                                                                                                      \n"
				"  int ndone;                                                                                                              \n"
				"                                                                                                                          \n"
				"  double schedtime=0.0;      /*time spent deciding and dispatching schedules*/                                            \n"
				"  double commtime=0.0;       /*index for communication latency*/                                                          \n"
				"  double totaltime=0.0;      /*index for seconds per schedule*/                                                           \n"
				"  double totalcommtime=0.0;  /*total communication latency*/                                                              \n"
				"  double paratime=0.0;       /*total parallel walltime for slaves excluding mslave*/                                      \n"
				"                                                                                                                          \n"
				"  /************* Scheduling Parameters **************/                                                                    \n"
				"  double calpha = 0.1;        /*weighting for communication average*/                                                     \n"
				"  double talpha = 0.1;        /*weighting for throughput average*/                                                        \n"
				"  double rcalpha = 1- calpha;                                                                                             \n"
				"  double rtalpha = 1- talpha;                                                                                             \n"
				"                                                                                                                          \n"
				"  double epsilon = 0.005;     /*maximum tolerated communication overhead*/                                                \n"
				"  double lower = 2.0;         /*minimum tolerated resolution in seconds*/                                                 \n"
				"  double upper = 10.0;        /*maximum tolerated resolution seconds*/                                                    \n"
				"  /***************************************************/                                                                   \n"
				"                                                                                                                          \n"
				"  pthread_t helper;           /*handle for slave thread*/                                                                 \n"
				"                                                                                                                          \n"
				"  double tp1, tp2;                                                                                                        \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  int master_half_step;       /*This indicates whether the master is sending*/                                            \n"
					"                              /*out half step paths or not - redundant to the index k*/                                   \n"
					"                                                                                                                          \n");
		fprintf(outfile,
				"  //Allocate memories for global data structures                                                                          \n"
				"  slave_stat = new int[size]; //(int *)malloc(sizeof(int)*size);                                                          \n"
				"  fschedcount = new int[size]; //(int *)malloc(sizeof(int)*size);                                                         \n"
				"  fthroughput = new double[size];                                                                                         \n"
				"  fcommave = new double[size];                                                                                            \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  hschedcount = new int[size];                                                                                            \n"
					"  hthroughput = new double[size];                                                                                         \n"
					"  hcommave = new double[size];                                                                                            \n");
		fprintf(outfile,
				"                                                                                                                          \n"
				"  //Initialise slave status arrays                                                                                        \n"
				"  for (i=0;i<size;i++){                                                                                                   \n"
				"	  slave_stat[i]=0;                                                                                                     \n"
				"	  partitions[i]=0;                                                                                                     \n"
				"	                                                                                                                       \n"
				"	  fschedcount[i]=0;                                                                                                    \n"
				"	  fcommave[i]=0.0;                                                                                                     \n"
				"	  fthroughput[i]=0.0;                                                                                                  \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"	                                                                                                                       \n"
					"	  hschedcount[i]=0;                                                                                                    \n"
					"	  hcommave[i]=0.0;                                                                                                     \n"
					"	  hthroughput[i]=0.0;                                                                                                  \n");
		fprintf(outfile,
				"  }                                                                                                                       \n"
				"                                                                                                                          \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  //Initially set to full step paths                                                                                      \n"
					"  _half_step = 0;                                                                                                         \n"
					"  master_half_step = _half_step;                                                                                          \n"
					"  sendtag = 1;                                                                                                            \n");
		fprintf(outfile,
				"                                                                                                                          \n"
				"  //Pthreads are always busy doing something                                                                              \n"
				"  slave_stat[0] = 1;                                                                                                      \n"
				"                                                                                                                          \n"
				"  //Initialise mutual exclusion mechanism                                                                                 \n"
				"  pthread_mutex_init(&tasklock, NULL);                                                                                    \n"
				"  pthread_mutex_init(&finlock, NULL);                                                                                     \n"
				"  pthread_mutex_lock(&finlock);                                                                                           \n"
				"                                                                                                                          \n"
				"  //Create a thread to act as a slave                                                                                     \n"
				"  if (pthread_create(&helper, NULL, mslave, NULL)!=0)                                                                     \n"
				"    printf(\"Thread creation failed\\n\");                                                                               \n"
				"                                                                                                                          \n"
				"  //Listen for messages from all slaves, including the master for some reason                                             \n"
				"  for (i=0;i<size;i++)                                                                                                    \n"
				"    MPI_Irecv(&bufs[i], 1, MPI_DOUBLE, i, MPI_ANY_TAG,MPI_COMM_WORLD,&reqs[i]);                                           \n"
				"                                                                                                                          \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  //Loop until all paths are finished, and all results recieved                                                           \n"
					"  //      -repeat                                                                                                         \n"
					"  //      -Test for messages from slaves if                                                                               \n"
					"  //      -Determine path schedule for idle slaves                                                                        \n"
					"  //      -Send schedule to slaves                                                                                        \n"
					"                                                                                                                          \n"
					"  /* The first iteration of this for loop allocates full paths                                                            \n"
					 "     The second iteration allocates half paths - master_half_step is a redundant                                          \n"
					 "     indicator for this too. We could use k interchangeably with it.*/                                                    \n"
					"                                                                                                                          \n"
					"  for (k=0;k<2;k++){                                                                                                      \n"
					"                                                                                                                          \n"
					"    if (master_half_step == 0){                                                                                           \n"
					"      active_throughput = fthroughput;                                                                                    \n"
					"      active_commave = fcommave;                                                                                          \n"
					"      active_assigned = &fullpaths_assigned;                                                                              \n"
					"    }                                                                                                                     \n"
					"    else if (master_half_step == 1){                                                                                      \n"
					"      active_throughput = hthroughput;                                                                                    \n"
					"      active_commave = hcommave;                                                                                          \n"
					"      active_assigned = &halfpaths_assigned;                                                                              \n"
					"    }                                                                                                                     \n"
					"                                                                                                                          \n"
					"                                                                                                                          \n"
					"    /**************LISTEN FOR \"TASKS COMPLETED\" MESSAGES FROM SLAVES*******************/                              \n"
					"                                                                                                                          \n"
					"    /*Loop until all fullpaths or,  if during the second iteration, halfpaths are assigned                                \n"
					 "      the outstanding > 0 condition makes the master poll for \"complete\" messages when                                  \n"
					 "      it has actually allocated all its full paths and half paths already.*/                                              \n"
					"                                                                                                                          \n"
					"    while ((master_half_step == 0 && fullpaths_assigned < total_paths)||                                                  \n"
					"	   (master_half_step == 1 && (halfpaths_assigned < total_paths || outstanding > 0))){                                  \n"
					"                                                                                                                          \n"
					"      /*Test whether there is an incoming message from slaves                                                             \n"
					 "        the idle == 0 condition is for the situation when SOME slaves                                                     \n"
					 "        have not finished their full paths whilst the master is ready to                                                  \n"
					 "        allocate half paths - this allows the master to allocate half paths                                               \n"
					 "        in the meantime, instead of waiting for completion of all full paths*/                                            \n"
					"                                                                                                                          \n"
					"      if (outstanding > 0 && (idle == 0 || (fullpaths_assigned + halfpaths_assigned >= 2*total_paths))){                  \n");
		else fprintf(outfile,
					 "  //Loop until all tasks are finished, and all results recieved                    \n"
					 "  //      -Check for requests from slaves for more tasks                           \n"
					 "  //      -Cyclically allocate a task to idle slaves to thread                     \n"
					 "                                                                                   \n"
					 "  while (fullpaths_assigned < total_paths || outstanding > 0){                     \n"
					 "                                                                                   \n"
					 "    if (outstanding > 0){                                                          \n");
		fprintf(outfile,
				"	//Wait for messages from slaves                                                                                        \n"
				"	                                                                                                                       \n"
				"      Testsome:                                                                                                           \n"
				"	MPI_Testsome(size,reqs, &ndone, indices, stats);                                                                       \n"
				"	                                                                                                                       \n"
				"	if (ndone == 0){                                                                                                       \n"
				"	  /*If no messages, sleep for x secs before checking again*/                                                           \n"
				"	  usleep(10000);                                                                                                       \n"
				"	  goto Testsome;                                                                                                       \n"
				"	}                                                                                                                      \n"
				"	                                                                                                                       \n"
				"	                                                                                                                       \n"
				"	for (i=0; i<ndone; i++){                                                                                               \n"
				"	  //Deal with incoming messages                                                                                        \n"
				"	  j = indices[i];                                                                                                      \n"
				"	                                                                                                                       \n"
				"	  //Dynamically determine bandwidth and throughput                                                                     \n"
				"	  totaltime = MPI_Wtime() - timing[j];                                                                                 \n"
				"	  commtime = totaltime - bufs[j];                                                                                      \n"
				"	                                                                                                                       \n"
				"	  //Calculate average communication time and average throughput                                                        \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"	  //MPI_TAG == 1 means fullpaths completed                                                                             \n"
					"	  //MPI_TAG == 2 means halfpaths completed                                                                             \n"
					"                                                                                                                          \n"
					"	  if (stats[i].MPI_TAG == 1){                                                                                          \n");
		fprintf(outfile,
				"	    if (fcommave[j] == 0.0){                                                                                           \n"
				"	      fcommave[j] = commtime;                                                                                          \n"
				"	      fthroughput[j] = partitions[j]/totaltime;                                                                        \n"
				"	    }	                                                                                                               \n"
				"	    else {                                                                                                             \n"
				"	      fcommave[j] = fcommave[j]*rcalpha + (commtime * calpha);                                                         \n"
				"	      fthroughput[j] = fthroughput[j]*rtalpha + (partitions[j]/totaltime)*talpha;                                      \n"
				"	    }                                                                                                                  \n"
				"                                                                                                                          \n"
				"#ifdef DEBUG_MPI_Scheduling                                                                                               \n"
				"	    fschedcount[j]+=partitions[j];                                                                                     \n"
				"	    printf(\"%%d: Rank %%d, Completed %%d tasks in %%5.5f seconds ::\", stats[i].MPI_TAG,j,partitions[j],commtime+bufs[j]);\n"
				"	    printf(\" T %%5.5f : C %%5.5f : R %%5.5f %%%%\\n\",fthroughput[j],fcommave[j],commtime/(bufs[j]+commtime)*100.0);        \n"
				"#endif                                                                                                                    \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"	  }                                                                                                                    \n"
					"	  else if (stats[i].MPI_TAG == 2){                                                                                     \n"
					"	    if (hcommave[j] == 0.0){                                                                                           \n"
					"	      hcommave[j] = commtime;                                                                                          \n"
					"	      hthroughput[j] = partitions[j]/totaltime;                                                                        \n"
					"	    }	                                                                                                               \n"
					"	    else {                                                                                                             \n"
					"	      hcommave[j] = hcommave[j]*rcalpha + (commtime * calpha);                                                         \n"
					"	      hthroughput[j] = hthroughput[j]*rtalpha + (partitions[j]/totaltime)*talpha;                                      \n"
					"	    }                                                                                                                  \n"
					"#ifdef DEBUG_MPI_Scheduling                                                                                               \n"
					"	    hschedcount[j]+=partitions[j];                                                                                     \n"
					"	    printf(\"%%d: Rank %%d, Completed %%d tasks in %%5.5f seconds ::\", stats[i].MPI_TAG,j,partitions[j],commtime+bufs[j]);\n"
					"	    printf(\" T %%5.5f : C %%5.5f : R %%5.5f %%%%\\n\",hthroughput[j],hcommave[j],commtime/(bufs[j]+commtime)*100.0);        \n"
					"#endif	                                                                                                                   \n"
					"	  }                                                                                                                    \n");
		fprintf(outfile,
				"	                                                                                                                       \n"
				"	  totalcommtime += commtime;                                                                                           \n"
				"	  paratime += totaltime;                                                                                               \n"
				"	  slave_stat[j] = 0;                                                                                                   \n"
				"	  idle++;                                                                                                              \n"
				"	  outstanding--;                                                                                                       \n"
				"	  MPI_Irecv(&bufs[j], 1, MPI_DOUBLE, j, MPI_ANY_TAG,MPI_COMM_WORLD,&reqs[j]);                                          \n"
				"	}                                                                                                                      \n"
				"      }                                                                                                                   \n"
				"                                                                                                                          \n"
				"      //Continue to listen for \"completed\" messages from slaves, if no more tasks need to be assigned                   \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"      if (halfpaths_assigned + fullpaths_assigned >= 2*total_paths)                                                       \n");
		else
			fprintf(outfile,
					"      if (fullpaths_assigned >= total_paths)																			   \n");
		fprintf(outfile,
				"			continue;                                                                                                      \n"
				"                                                                                                                          \n"
				"      /********************SCHEDULE MORE TASKS FOR IDLE SLAVES***********************/                                    \n"
				"                                                                                                                          \n"
				"      for (i=0;i<size;i++)                                                                                                \n"
				"	schedule[i]=0;                                                                                                         \n"
				"                                                                                                                          \n"
				"      tp1 = MPI_Wtime();                                                                                                  \n"
				"                                                                                                                          \n"
				"      pthread_mutex_lock(&tasklock);                                                                                      \n"
				"                                                                                                                          \n"
				"      //allocate tasks to free processors                                                                                 \n"
				"      //other scheduling algorithms may be inserted here if needed                                                        \n"
				"      //this one is a FIFO algorithm                                                                                      \n"
				"      //scheduling must be mutually exclusive as the slave thread                                                         \n"
				"      //also modifies the global variables below for self-scheduling                                                      \n"
				"                                                                                                                          \n"
				"      for (i=1;i<size;i++){                                                                                               \n"
				"	                                                                                                                       \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"	if (*active_assigned >= total_paths)                                                                                   \n");
		else
			fprintf(outfile,
					"	if (fullpaths_assigned >= total_paths)                                                                                 \n");
		fprintf(outfile,
				"	  break;                                                                                                               \n"
				"	                                                                                                                       \n"
				"	//only allocate more tasks to slaves that are idle                                                                     \n"
				"	if (slave_stat[i] == 0){                                                                                               \n"
				"	  //printf(\"allocating tasks to rank %%d\\n\",i);                                                                       \n"
				"	  slave_stat[i]=1;                                                                                                     \n"
				"	  idle--;                                                                                                              \n"
				"                                                                                                                          \n"
				"	  //Determine partition size based on slave throughput and                                                             \n"
				"	  //communication overhead. Preferable estimated computing times                                                       \n"
				"	  //for each schedule is high enough to reduce comm overhead                                                           \n"
				"	  //and between upper and lower                                                                                        \n"
				"                                                                                                                          \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"	  partitions[i] = (int) (MAX((active_commave[i]*active_throughput[i])/epsilon, active_throughput[i]*lower));           \n"
					"	  partitions[i] = (int) (MIN(partitions[i], active_throughput[i]*upper));                                              \n"
					"	  partitions[i] = (int) (MAX(partitions[i], 1));                                                                       \n"
					"	                                                                                                                       \n"
					"	  if (*active_assigned + partitions[i] <= total_paths){                                                                \n"
					"	    schedule[i] = partitions[i];                                                                                       \n"
					"	    *active_assigned += partitions[i];                                                                                 \n"
					"	  }                                                                                                                    \n"
					"	  else {                                                                                                               \n"
					"	    schedule[i] = total_paths-*active_assigned;                                                                        \n"
					"	    partitions[i] = schedule[i];                                                                                       \n"
					"	    *active_assigned += schedule[i];                                                                                   \n");
		else 
			fprintf(outfile,
					"	  	partitions[i] = (int) (MAX((fcommave[i]*fthroughput[i])/epsilon, fthroughput[i]*lower));							\n"
					"       partitions[i] = (int) (MIN(partitions[i],fthroughput[i]*upper));													\n"
					"		partitions[i] = (int) (MAX(partitions[i], 1));																		\n"
					"											\n"
					"		if (fullpaths_assigned + partitions[i] <= total_paths){																\n"
					"			schedule[i] = partitions[i];																					\n"
					"			fullpaths_assigned += partitions[i];																			\n"
					"		}																													\n"
					"		else {																												\n"
					"			schedule[i] = total_paths-fullpaths_assigned;																	\n"
					"			partitions[i] = schedule[i];																					\n"
					"			fullpaths_assigned += schedule[i];                                                                              \n");
		fprintf(outfile,
				"	  }                                                                                                                    \n"
				"	}	                                                                                                                   \n"
				"      }                                                                                                                   \n"
				"                                                                                                                          \n"
				"      pthread_mutex_unlock(&tasklock);                                                                                    \n"
				"                                                                                                                          \n"
				"      /**************************SEND SCHEDULE TO SLAVE(S)********************/                                           \n"
				"      //send instructions via MPI to slaves that have been shcheduled                                                     \n"
				"      for (i=1;i<size;i++){                                                                                               \n"
				"	if (schedule[i]>0){                                                                                                    \n"
				"	  timing[i] = MPI_Wtime();                                                                                             \n"
				"	  MPI_Send(&schedule[i],1,MPI_INT,i,sendtag,MPI_COMM_WORLD);                                                           \n"
				"	  outstanding++;                                                                                                       \n"
				"	}                                                                                                                      \n"
				"      }                                                                                                                   \n"
				"      tp2 = MPI_Wtime() - tp1;                                                                                            \n"
				"      schedtime += tp2;                                                                                                   \n"
				"    }                                                                                                                     \n"
				"                                                                                                                          \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"    //Initialise variables to do half steps                                                                               \n"
					"    master_half_step = 1;                                                                                                 \n"
					"    sendtag = 2;                                                                                                          \n"
					"  }                                                                                                                       \n"
					"                                                                                                                          \n");
		fprintf(outfile,
				"  //Block until the thread slave has completed processing                                                                 \n"
				"  pthread_mutex_lock(&finlock);                                                                                           \n"
				"                                                                                                                          \n"
				"  //Tell slave processes to Reduce then exit                                                                              \n"
				"  for (i=0;i<size;i++){                                                                                                   \n"
				"    MPI_Send(NULL, 0, MPI_INT,i,0,MPI_COMM_WORLD);                                                                        \n"
				"  }                                                                                                                       \n"
				"                                                                                                                          \n"
				"  //Kill slave thread                                                                                                     \n"
				"  pthread_cancel(helper);                                                                                                 \n"
				"                                                                                                                          \n"
				"#ifdef DEBUG_MPI_Scheduling                                                                                               \n"
				"  printf(\"\\n**************************************\");                                                                   \n"
				"  printf(\"\\nFullpaths returned by each process:\\n\");                                                                    \n"
				"  print_array2(fschedcount, size);                                                                                        \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  printf(\"\\nHalfpaths returned by each process:\\n\");                                                                    \n"
					"  print_array2(hschedcount, size);                                                                                        \n");
		fprintf(outfile,
				"                                                                                                                          \n"
				"  printf(\"\\nEstimated fullpath throughput and commtime\\n\");                                                             \n"
				"  print_array(fthroughput, size);printf(\"\\n\");                                                                          \n"
				"  print_array(fcommave, size);                                                                                            \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  printf(\"\\nEstimated halfpath throughput and commtime\\n\");                                                             \n"
					"  print_array(hthroughput, size);printf(\"\\n\");                                                                          \n"
					"  print_array(hcommave, size);                                                                                            \n");
		fprintf(outfile,
				"                                                                                                                          \n"
				"  printf(\"\\n\");                                                                                                         \n"
				"  printf(\"Scheduling time = %%5.5f sec\\n\",schedtime);                                                                    \n"
				"  printf(\"Communication time = %%5.5f sec\\n\", totalcommtime);                                                            \n"
				"  printf(\"Total parallel time (excl thread) = %%5.5f sec\\n\",paratime);                                                   \n"
				"  printf(\"**************************************\\n\\n\");                                                                 \n"
				"  fflush(stdout);                                                                                                         \n"
				"#endif                                                                                                                    \n"
				"                                                                                                                          \n"
				"  /*Free memory*/	                                                                                                       \n"
				"  delete[] slave_stat;                                                                                                    \n"
				"  delete[] fschedcount;                                                                                                   \n"
				"  delete[] fthroughput;                                                                                                   \n"
				"  delete[] fcommave;                                                                                                      \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  delete[] hschedcount;                                                                                                   \n"
					"  delete[] hthroughput;                                                                                                   \n"
					"  delete[] hcommave;                                                                                                      \n"
					"                                                                                                                          \n");
		fprintf(outfile,
				"  return 0;                                                                                                               \n"
				"}                                                                                                                         \n"
				"                                                                                                                          \n"
				"                                                                                                                          \n"
				"                                                                                                                          \n"
				"void *mslave(void *ptr){                                                                                                  \n"
				"  //This is the threaded slave, which is spawned by the master to do tasks in parallel                                    \n"
				"  //it self schedules itself more task batches. It takes in one argument, which is                                        \n"
				"  //the size of the solution array.                                                                                       \n"
				"                                                                                                                          \n"
				"  int thr_batchsize;                                                                                                      \n"
				"  double thr_throughput=0.0;                                                                                              \n"
				"  double thr_time_per_batch=2.0;                                                                                          \n"
				"  double thr_talpha=0.1;                                                                                                  \n"
				"                                                                                                                          \n"
				"  int i;                                                                                                                  \n"
				"  struct timeval *tp1,*tp2;                                                                                               \n"
				"  double tp3;                                                                                                             \n"
				"                                                                                                                          \n"
				"                                                                                                                          \n"
				"  //Initialise timing structure                                                                                           \n"
				"  tp1 = new struct timeval; //(struct timeval *)malloc(sizeof(struct timeval));                                           \n"
				"  tp2 = new struct timeval; //(struct timeval *)malloc(sizeof(struct timeval));                                           \n"
				"                                                                                                                          \n"
				"  thr_batchsize = 1;                                                                                                      \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  _half_step = 0;                                                                                                         \n"
					"                                                                                                                          \n"
					"  while(fullpaths_assigned + halfpaths_assigned < 2*total_paths){                                                         \n");
		else fprintf(outfile,
					 "  while(fullpaths_assigned < total_paths){																				   \n");
		fprintf(outfile,
				"    //Self schedule more tasks to process                                                                                 \n"
				"    gettimeofday(tp1,NULL);                                                                                               \n"
				"                                                                                                                          \n"
				"    /********************SCHEDULE MORE TASKS*************************/                                                    \n"
				"    pthread_mutex_lock(&tasklock);                                                                                        \n"
				"                                                                                                                          \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"    if (halfpaths_assigned + fullpaths_assigned >= 2*total_paths)                                                         \n");
		else fprintf(outfile,
					 "    if (fullpaths_assigned >= total_paths)																				   \n");
		fprintf(outfile,
				"      break;                                                                                                              \n"
				"                                                                                                                          \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"    //Set to half step (once), if all full step paths are assigned                                                        \n"
					"    if (fullpaths_assigned >= total_paths && _half_step == 0){                                                            \n"
					"      printf(\"Rank[%%i] beginning half step paths\\n\",rank);                                                              \n"
					"      fthroughput[0]=thr_throughput;                                                                                      \n"
					"      _half_step = 1;                                                                                                     \n"
					"      thr_batchsize = 1;                                                                                                  \n"
					"      thr_throughput = 0.0;                                                                                               \n"
					"                                                                                                                          \n"
					"      _gen1[0] = 346+rank;                                                                                                \n"
					"      _gen2[0] = 983+rank;                                                                                                \n"
					"      _gen1[1] = 0;                                                                                                       \n"
					"      _gen2[1] = 0;                                                                                                       \n"
					"      _gen1[2] = 0;                                                                                                       \n"
					"      _gen2[2] = 0;                                                                                                       \n"
					"      erand48(_gen1);                                                                                                     \n"
					"      erand48(_gen2);                                                                                                     \n"
					"    }                                                                                                                     \n"
					"                                                                                                                          \n"
					"    if (_half_step==0){                                                                                                   \n");
		fprintf(outfile,
				"      if (fullpaths_assigned + thr_batchsize <= total_paths)                                                              \n"
				"	fullpaths_assigned += thr_batchsize;                                                                                   \n"
				"      else {                                                                                                              \n"
				"	thr_batchsize = (total_paths - fullpaths_assigned);                                                                    \n"
				"       fullpaths_assigned += thr_batchsize;                                                                               \n"
				"      }                                                                                                                   \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"    }                                                                                                                     \n"
					"    else if (_half_step == 1){                                                                                            \n"
					"     if (halfpaths_assigned+thr_batchsize <= total_paths)                                                                 \n"
					"       halfpaths_assigned += thr_batchsize;                                                                               \n"
					"     else {                                                                                                               \n"
					"       thr_batchsize = (total_paths - halfpaths_assigned);                                                                \n"
					"       halfpaths_assigned += thr_batchsize;                                                                               \n"
					"     }                                                                                                                    \n"
					"    }                                                                                                                     \n");
		fprintf(outfile,
				"                                                                                                                          \n"
				"   pthread_mutex_unlock(&tasklock);                                                                                       \n"
				"                                                                                                                          \n"
				"   /*****************************************************************/                                                    \n"
				"                                                                                                                          \n"
				"   _segment0(thr_batchsize);                                                                                              \n"
				"                                                                                                                          \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"#ifdef DEBUG_MPI_Scheduling                                                                                               \n"
					"   if (_half_step == 0)                                                                                                   \n");
		fprintf(outfile,
				"     fschedcount[0]+=thr_batchsize;                                                                                       \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"   else                                                                                                                   \n"
					"     hschedcount[0]+=thr_batchsize;                                                                                       \n"
					"#endif                                                                                                                    \n");
		fprintf(outfile,
				"                                                                                                                          \n"
				"   gettimeofday(tp2, NULL);                                                                                               \n"
				"                                                                                                                          \n"
				"   /********************CALCULATE NEW BATCH SIZE*********************/                                                    \n"
				"                                                                                                                          \n"
				"   tp3 = ((tp2->tv_sec - tp1->tv_sec)+((tp2->tv_usec-tp1->tv_usec)/1000000.0));                                           \n"
				"                                                                                                                          \n"
				"   if (thr_throughput == 0.0)                                                                                             \n"
				"     thr_throughput = thr_batchsize/tp3;                                                                                  \n"
				"   else                                                                                                                   \n"
				"     thr_throughput = (1-thr_talpha)*thr_throughput + thr_talpha * (thr_batchsize/tp3);                                   \n"
				"                                                                                                                          \n"
				"#ifdef DEBUG_MPI_Scheduling                                                                                               \n"
				"   printf(\"Rank 0, Completed %%d tasks in %%5.5f secs, T %%5.5f\\n\",thr_batchsize, tp3,thr_throughput);                     \n"
				"#endif                                                                                                                    \n"
				"                                                                                                                          \n"
				"   thr_batchsize = MAX(1,(int) (thr_throughput * thr_time_per_batch));                                                    \n"
				"   /*****************************************************************/                                                    \n"
				"  }                                                                                                                       \n"
				"                                                                                                                          \n"
				"  //Unlocking indicates that the thread slave has finished processing                                                     \n"
				"  //and allows the master to gather results it has produced                                                               \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"  hthroughput[0] = thr_throughput;                                                                                        \n"
					"  delete tp1;                                                                                                             \n"
					"  delete tp2;	                                                                                                           \n");
		else
			fprintf(outfile,
					"  fthroughput[0] = thr_throughput;                                                                                        \n");
		fprintf(outfile,
				"  pthread_mutex_unlock(&finlock);                                                                                         \n"
				"  return 0;																													\n"
				"}                                                                                                                         \n"
				"                                                                                                                          \n"
				"void *slave(void *ptr){                                                                                                   \n"
				"  //This is a slave that runs on another process, it requests more tasks from the master                                  \n"
				"  //when it finishes its current task. It takes the length of the solution array as                                       \n"
				"  //argument.                                                                                                             \n"
				"                                                                                                                          \n"
				"  MPI_Status stat;                                                                                                        \n"
				"  int count;                    /*batch size*/                                                                            \n"
				"  int i;                                                                                                                  \n"
				"  double tp1, tp2;                                                                                                        \n"
				"  unsigned long _length;                                                                                                  \n"
				"                                                                                                                          \n"
				"  while(1){                                                                                                               \n"
				"    //Wait for task from master                                                                                           \n"
				"    MPI_Recv(&count, 1, MPI_INT,0,MPI_ANY_TAG,MPI_COMM_WORLD,&stat);                                                      \n"
				"                                                                                                                          \n"
				"    tp1 = MPI_Wtime();                                                                                                    \n"
				"                                                                                                                          \n"
				"    //Check message type                                                                                                  \n"
				"    if (stat.MPI_TAG == 0) {                                                                                              \n"
				"      break;                                                                                                              \n"
				"    }                                                                                                                     \n");
		if(simulation()->parameters()->errorCheck) 
			fprintf(outfile,
					"    else if (stat.MPI_TAG == 1){                                                                                          \n"
					"      _half_step = 0;                                                                                                     \n"
					"    }                                                                                                                     \n"
					"    else if (stat.MPI_TAG == 2 && _half_step == 0){                                                                       \n"
					"      printf(\"Rank[%%i] beginning half step paths\\n\",rank);                                                              \n"
					"      _half_step = 1;                                                                                                     \n"
					"      _gen1[0] = 346+rank;                                                                                                \n"
					"      _gen2[0] = 983+rank;                                                                                                \n"
					"      _gen1[1] = 0;                                                                                                       \n"
					"      _gen2[1] = 0;                                                                                                       \n"
					"      _gen1[2] = 0;                                                                                                       \n"
					"      _gen2[2] = 0;                                                                                                       \n"
					"      erand48(_gen1);                                                                                                     \n"
					"      erand48(_gen2);                                                                                                     \n"
					"    }                                                                                                                     \n");
		fprintf(outfile,
				"                                                                                                                          \n"
				"    _segment0(count);                                                                                                     \n"
				"                                                                                                                          \n"
				"    tp2 = MPI_Wtime() - tp1;                                                                                              \n"
				"                                                                                                                          \n"
				"    //Send completion notice to the master                                                                                \n"
				"    MPI_Send(&tp2,1,MPI_DOUBLE,0,stat.MPI_TAG,MPI_COMM_WORLD);                                                            \n"
				"  }                                                                                                                       \n"
				"  //Tell master that slave is done                                                                                        \n"
				"  MPI_Send(NULL, 0, MPI_DOUBLE,0,0,MPI_COMM_WORLD);                                                                       \n"
				"  return 0;																													\n"
				"}                                                                                                                         \n"
				"	                                                                                                                       \n"
				"//Debugging functions                                                                                                     \n"
				"                                                                                                                          \n"
				"int print_array(double *array, int size){                                                                                 \n"
				"  int i, j;                                                                                                               \n"
				"  for (i=0, j=0 ; i<size ; i++){                                                                                          \n"
				"    printf(\"%%9.2f \", array[i]);                                                                                         \n"
				"    j++;                                                                                                                  \n"
				"    if (j==10){                                                                                                           \n"
				"      printf(\"\\n\");                                                                                                     \n"
				"      j = 0;                                                                                                              \n"
				"    }                                                                                                                     \n"
				"  }                                                                                                                       \n"
				"  return 0;                                                                                                               \n"
				"}                                                                                                                         \n"
				"                                                                                                                          \n"
				"int print_array2(int *array, int size){                                                                                   \n"
				" int i, j;                                                                                                                \n"
				"  for (i=0, j=0 ; i<size ; i++){                                                                                          \n"
				"    printf(\"%%d: %%d\\n\", i, array[i]);                                                                                    \n"
				"    j++;                                                                                                                  \n"
				"    if (j==10){                                                                                                           \n"
				"      printf(\"\\n\");                                                                                                     \n"
				"      j = 0;                                                                                                              \n"
				"    }                                                                                                                     \n"
				"  }                                                                                                                       \n"
				"  return 0;                                                                                                               \n"
				"}                                                                                                                         \n\n");
	}


    fprintf(outfile,"// *************************\n");

	if(!simulation()->parameters()->usempi)
			fprintf(outfile,"void _segment0() {\n");
	  else if(simulation()->parameters()->mpiMethod=="Scheduling")
		  fprintf(outfile,"void _segment0(int count) {\n");
	  else 
		  fprintf(outfile,"void _segment0(const int& size) {\n");
	  
	  
    fprintf(outfile,"\n");

    for(list<const xmdsSegment*>:: const_iterator ppxmdsSegment
	  = mySegmentList.begin(); ppxmdsSegment !=
	  mySegmentList.end(); ppxmdsSegment++) {
      (*ppxmdsSegment)->writeInitialisationCalls(outfile);
    }

    if((simulation()->parameters()->stochastic)&&(!(simulation()->parameters()->mpiMethod=="Scheduling")||(!simulation()->parameters()->usempi))) {
		if(simulation()->parameters()->errorCheck) {
			if(simulation()->parameters()->usempi) {
				fprintf(outfile,"_gen1[0] = %li+rank;\n",simulation()->parameters()->seed[0]);
				fprintf(outfile,"_gen2[0] = %li+rank;\n",simulation()->parameters()->seed[1]);
			}
			else {
				fprintf(outfile,"_gen1[0] = %li;\n",simulation()->parameters()->seed[0]);
				fprintf(outfile,"_gen2[0] = %li;\n",simulation()->parameters()->seed[1]);
			}
			fprintf(outfile,"_gen1[1] = 0;\n");
			fprintf(outfile,"_gen2[1] = 0;\n");
			fprintf(outfile,"_gen1[2] = 0;\n");
			fprintf(outfile,"_gen2[2] = 0;\n");
			fprintf(outfile,"erand48(_gen1);\n");
			fprintf(outfile,"erand48(_gen2);\n");
		}
		else {
			if(simulation()->parameters()->usempi) {
				fprintf(outfile,"_gen[0] = %li+rank;\n",simulation()->parameters()->seed[0]);
			}
			else  {
				fprintf(outfile,"_gen[0] = %li;\n",simulation()->parameters()->seed[0]);
			}
			
			fprintf(outfile,"_gen[1] = 0;\n");
			fprintf(outfile,"_gen[2] = 0;\n");
			fprintf(outfile,"erand48(_gen);\n");
		}
    }

    const char* indent;

	if((simulation()->parameters()->mpiMethod=="Scheduling")&&(simulation()->parameters()->stochastic)&&(simulation()->parameters()->usempi)) {
		fprintf(outfile,"  for(long _i0=0;_i0<count;_i0++)");
		indent="	";
		fprintf(outfile,"   {\n");
	}
	else {
		if(simulation()->parameters()->nPaths>1) {
			if(simulation()->parameters()->usempi) {
				fprintf(outfile,"  for(long _i0=rank;_i0<_n_paths;_i0+=size)");
			}
			else {
				fprintf(outfile,"  for(long _i0=0;_i0<_n_paths;_i0++)");
			}
			indent="	";
			fprintf(outfile,"   {\n");
		}
		else
			indent="";
		fprintf(outfile,"\n");
	}
	
    fprintf(outfile,"%s%s=0;\n",indent,simulation()->parameters()->propDimName.c_str());
    fprintf(outfile,"\n");

    simulation()->field()->writeInitialisationCalls(outfile,indent);

    for(unsigned long i=0;i<simulation()->output()->nMomentGroups();i++) {
      fprintf(outfile,"%s_mg%li_sample_pointer=0;\n",indent,i);
      if(simulation()->output()->momentGroup(i)->requiresIntegrations()|(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic)) {
	fprintf(outfile,"%s_mg%li_raw_initialise();\n",indent,i);
      }
    }
    fprintf(outfile,"\n");

    simulation()->field()->writeSampleCalls(outfile,indent);

	if((simulation()->parameters()->mpiMethod=="Scheduling")&&(simulation()->parameters()->stochastic)&&(simulation()->parameters()->usempi)) {
		//Write nothing
	} 
	else if(simulation()->parameters()->nPaths>1) {
		fprintf(outfile,"%sprintf(\"Starting path %%li\\n\",_i0+1);\n",indent);
    }
	
    for(list<const xmdsSegment*>::const_iterator ppxmdsSegment =
	  mySegmentList.begin(); ppxmdsSegment !=
	  mySegmentList.end(); ppxmdsSegment++) {
	fprintf(outfile,"%s_segment%li(1);\n",indent,(*ppxmdsSegment)->segmentNumber);
      }
    fprintf(outfile,"\n");

    for(unsigned long i=0;i<simulation()->output()->nMomentGroups();i++) {
      fprintf(outfile,"%s_mg%li_process();\n",indent,i);
    }

	if((simulation()->parameters()->mpiMethod=="Scheduling")&&(simulation()->parameters()->stochastic)&&(simulation()->parameters()->usempi)) {
		fprintf(outfile,"	}\n");
	}
	else if(simulation()->parameters()->nPaths>1) {
		fprintf(outfile,"	}\n");
    }
	
    fprintf(outfile,"}\n");
    fprintf(outfile,"\n");
  }
  else {
    // I am a sub-sequence

    if(verbose()) {
      printf("Writing sequence routine ...\n");
    }

    fprintf(outfile,"// ********************************************************\n");
    fprintf(outfile,"// segment %li (sequence) routines\n",segmentNumber);
    fprintf(outfile,"\n");

    fprintf(outfile,"// *************************\n");
    fprintf(outfile,"void _segment%li(unsigned long cycle) {\n",segmentNumber);
        
    for(list<const xmdsSegment*>:: const_iterator ppxmdsSegment
	  = mySegmentList.begin(); ppxmdsSegment !=
	  mySegmentList.end(); ppxmdsSegment++) {
      (*ppxmdsSegment)->writeInitialisationCalls(outfile);
    }

    fprintf(outfile,"for(unsigned long i=0; i<%li; i++) {\n",nCycles);
    for(list<const xmdsSegment*>::const_iterator ppxmdsSegment =
	  mySegmentList.begin(); ppxmdsSegment !=
	  mySegmentList.end(); ppxmdsSegment++) {
      fprintf(outfile,"	_segment%li(i+1);\n",
	      (*ppxmdsSegment)->segmentNumber);
    }
    fprintf(outfile,"	}\n");
    fprintf(outfile,"}\n");
    fprintf(outfile,"\n");
  }

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

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

  xmdsSegment* newxmdsSegment = new xmdsSequence(simulation(),verbose());
  addChild((xmdsElement*) newxmdsSegment);
  mySegmentList.push_back(newxmdsSegment);
  return newxmdsSegment;
}

// ******************************************************************************
xmdsSegment* xmdsSequence::createxmdsIntegrateRK4EX() {
  if(debugFlag) {
    printf("xmdsSequence::createxmdsIntegrateRK4EX\n");
  }

  xmdsSegment* newxmdsSegment = new xmdsIntegrateRK4EX(simulation(),verbose());
  addChild((xmdsElement*) newxmdsSegment);
  mySegmentList.push_back(newxmdsSegment);
  return newxmdsSegment;
}

// ******************************************************************************
xmdsSegment* xmdsSequence::createxmdsIntegrateRK4IP() {
  if(debugFlag) {
    printf("xmdsSequence::createxmdsIntegrateRK4IP\n");
  }

  xmdsSegment* newxmdsSegment = new xmdsIntegrateRK4IP(simulation(),verbose());
  addChild((xmdsElement*) newxmdsSegment);
  mySegmentList.push_back(newxmdsSegment);
  return newxmdsSegment;
}

// ******************************************************************************
xmdsSegment* xmdsSequence::createxmdsIntegrateARK45EX() {
  if(debugFlag) {
    printf("xmdsSequence::createxmdsIntegrateARK45EX\n");
  }

  xmdsSegment* newxmdsSegment = new xmdsIntegrateARK45EX(simulation(),verbose());
  addChild((xmdsElement*) newxmdsSegment);
  mySegmentList.push_back(newxmdsSegment);
  return newxmdsSegment;
}

// ******************************************************************************
xmdsSegment* xmdsSequence::createxmdsIntegrateARK45IP() {
  if(debugFlag) {
    printf("xmdsSequence::createxmdsIntegrateARK45IP\n");
  }

  xmdsSegment* newxmdsSegment = new xmdsIntegrateARK45IP(simulation(),verbose());
  addChild((xmdsElement*) newxmdsSegment);
  mySegmentList.push_back(newxmdsSegment);
  return newxmdsSegment;
}

// ******************************************************************************
xmdsSegment* xmdsSequence::createxmdsIntegrateSIEX() {
  if(debugFlag) {
    printf("xmdsSequence::createxmdsIntegrateSIEX\n");
  }

  xmdsSegment* newxmdsSegment = new xmdsIntegrateSIEX(simulation(),verbose());
  addChild((xmdsElement*) newxmdsSegment);
  mySegmentList.push_back(newxmdsSegment);
  return newxmdsSegment;
}

// ******************************************************************************
xmdsSegment* xmdsSequence::createxmdsIntegrateSIIP() {
  if(debugFlag) {
    printf("xmdsSequence::createxmdsIntegrateSIIP\n");
  }

  xmdsSegment* newxmdsSegment = new xmdsIntegrateSIIP(simulation(),verbose());
  addChild((xmdsElement*) newxmdsSegment);
  mySegmentList.push_back(newxmdsSegment);
  return newxmdsSegment;
}

// ******************************************************************************
xmdsSegment* xmdsSequence::createxmdsFilter() {
  if(debugFlag) {
    printf("xmdsSequence::createxmdsIntegrateFilter\n");
  }

  xmdsSegment* newxmdsSegment = new xmdsFilter(simulation(),verbose());
  addChild((xmdsElement*) newxmdsSegment);
  mySegmentList.push_back(newxmdsSegment);
  return newxmdsSegment;
}