/*
  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: xmdsintegrate.cc,v 1.23 2005/05/02 02:35:58 sebwuester Exp $
*/

/*! @file xmdsintegrate.cc
  @brief Integrate element parsing classes and methods

  More detailed explanation...
*/

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

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrate public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsIntegrates=0;  //!< Number of xmds integrate objects

// ******************************************************************************
xmdsIntegrate::xmdsIntegrate(
			     const xmdsSimulation *const yourSimulation,
			     const bool& yourVerboseMode, const bool& adaptiveIP) :
  xmdsSegment(yourSimulation,yourVerboseMode) {
  myAdaptiveIP = adaptiveIP;
  if(debugFlag) {
    nxmdsIntegrates++;
    printf("xmdsIntegrate::xmdsIntegrate\n");
    printf("nxmdsIntegrates=%li\n",nxmdsIntegrates);
  }
  myNextcoKey=0;
};

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

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


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

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

  // ************************************
  // find 'interval'

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

  myInterval=*myXMLStringList.begin();

  if(verbose()) {
    printf("integrate interval = %s\n",myInterval.c_str());
  }

  // ************************************
  // find 'lattice'

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

  myLattice=*myULongList.begin();

  if(myLattice==0) {
    throw xmdsException(yourElement,"Integration lattice must be >= 1 !");
  }

  if(verbose()) {
    printf("integrate lattice has %li points\n",myLattice);
  }

  // ************************************
  // find 'samples'

  getAssignmentULongs(yourElement,"samples",1,simulation()->output()->nMomentGroups(),mySamplesList);

  for(long unsigned int i=0;i<mySamplesList.size();i++) {
    if(samples(i) != 0) {
      if(myLattice != samples(i) * (myLattice/samples(i))) {
	throw xmdsException(yourElement,"moment group samples must be either 0 or else factors of the integration lattice");
      }
      if(verbose()) {
	printf("moment group #%li is sampled every %li points\n",i+1,myLattice/samples(i));
      }
    }
    else {
      if(verbose()) {
	printf("moment group #%li is not sampled\n",i+1);
      }
    }
  }
  
 // ************************************
  // find 'smallmemory'
  
  if(myAdaptiveIP){

  list<bool> myBoolList;

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

   if(myBoolList.size()==1) {
     
    mySmallmemory=*myBoolList.begin();
    if(verbose()) 
      if(mySmallmemory) {
        printf("\"save memory\" in k_propagate is 'yes'\n");
      }
      else {
        printf("\"save memory\" in k _propagate 'no'\n");
      }
  }
  else {
      printf("\"save memory\" in k_propagate defaulting to no \n");
      mySmallmemory=false;
       }
    }

  // ************************************
  // find nonoises flag

  list<bool> tempNonoise;
  myNonoises=false;
  
  getAssignmentBools(yourElement,"no_noise",0,1,tempNonoise);
  
  if(tempNonoise.size()==1)
    {
    if(*tempNonoise.begin()) {
        myNonoises = true;
        if(verbose())
            printf("No noises will be definied in the filter.\n");
        }
  }

  if(simulation()->field()->geometry()->nDims()>0) {

    // ************************************
    // find k_operators element

    const NodeList *const candidateElements = yourElement->getElementsByTagName("k_operators",0);

    if(candidateElements->length()>1) {
      throw xmdsException(yourElement,"Multiple <k_operators> elements defined");
    }

    if(candidateElements->length()==1) {

      // we have a k_operators element

      const Element *const kPropagationElement = dynamic_cast<const Element*> (candidateElements->item(0));

      if(verbose()) {
	printf("Processing <k_operators> element ...\n");
      }

      // ************************************
      // find constantK

      list<bool> aBoolList;

      getAssignmentBools(kPropagationElement,"constant",0,1,aBoolList);

      if(aBoolList.size()>0) {
	myConstantK=*aBoolList.begin();
	if(verbose()) {
	  if(myConstantK) {
	    printf("k-operators are propagation independent\n");
	  }
	  else {
	    printf("k-propagation are propagation dependent\n");
	  }
	}
      }
      else {
	myConstantK=1;
	printf("k-operators defaulting to propagation independent\n");
      }
      
    if(!myConstantK && myAdaptiveIP){  
      throw xmdsException("ARK45IP does not work with time-dependent k-operators. Use ARK45EX instead.");
    }

      // ************************************
      // find vectors

      getAssignmentStrings(kPropagationElement,"vectors",0,0,myKVectorNamesList);

      if(myKVectorNamesList.size()==0) {
	// no vectors specified, therefore assume using only main vector
	myKVectorNamesList.push_back("main");
      }

      simulation()->field()->processVectors(myKVectorNamesList,simulation()->field()->geometry()->fullSpace());

      // ************************************
      // find operators

      getAssignmentStrings(kPropagationElement,"operator_names",1,0,myKOperatorNamesList);

      // ************************************
      // find code

      myKOperatorsCode=*kPropagationElement->textContent(0);

      // check for non-white space charaters in code:

      if(myKOperatorsCode.isAllWhiteSpace()) {
	throw xmdsException(kPropagationElement,"No <k_operators> code defined!");
      }

      if(verbose()) {
	printf("<k_operators> code loaded\n");
      }
    }
  }

  // ************************************
  // find vectors

  getAssignmentStrings(yourElement,"vectors",1,0,myVectorNamesList);

  // Vectors are now required by law in order to place the code element
  /*  if(myVectorNamesList.size()==0) {
  // no vectors specified, therefore assume using only main vector
  myVectorNamesList.push_back("main");
  }*/

  simulation()->field()->processVectors(myVectorNamesList,0);

  // ************************************
  // find code

  myPropagationCode=*yourElement->textContent(0);

  // check for non-white space charaters in code:

  if(myPropagationCode.isAllWhiteSpace()) {
    throw xmdsException(yourElement,"No propagation code defined!");
  }

  if(verbose()) {
    printf("processing propagation code ...\n");
  }


  if(simulation()->field()->geometry()->nDims()>0) {

    // ************************************
    // find cross_propagation

    const NodeList *const candidateElements = yourElement->getElementsByTagName("cross_propagation",0);

    if(candidateElements->length() > 1) {
      throw xmdsException(yourElement,"only one <cross_propagation> allowed");
    }

    if(candidateElements->length() == 1) {

      // we have a cross_propagation element
                
      const Element *const crossPropagationElement = dynamic_cast<const Element*> (candidateElements->item(0));

      // ...but not if we're distributing the vector over many processors!
      if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	throw xmdsException(crossPropagationElement,"I really can't handle using MPI when there are cross-propagation elements.");
      }

      // ************************************
      // find cross propagating vectors

      getAssignmentStrings(crossPropagationElement,"vectors",1,0,myCrossVectorNamesList);

      for(list<XMLString>::const_iterator pXMLString = myCrossVectorNamesList.begin(); pXMLString != myCrossVectorNamesList.end(); pXMLString++) {
	if(*pXMLString=="main") {
	  throw xmdsException(crossPropagationElement,"Cannot cross-propagate 'main' vector!");
	}
      }

      simulation()->field()->processVectors(myCrossVectorNamesList,0);

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

      list<XMLString> anXMLStringList;

      getAssignmentStrings(crossPropagationElement,"prop_dim",1,1,anXMLStringList);

      // now find which member of the main field dimensions it is
      if(!simulation()->field()->geometry()->getDimNumber(*anXMLStringList.begin(),myCrossDimNumber)) {
	sprintf(errorMessage(),"prop_dim '%s' not a dimension of main field",anXMLStringList.begin()->c_str());
	throw xmdsException(crossPropagationElement,errorMessage());
      }

      if(verbose()) {
	printf("cross prop_dim is field dimension number %li\n",myCrossDimNumber);
      }

      // ************************************
      // find cross_propagation code

      myCrossPropagationCode=*crossPropagationElement->textContent(0);

      // check for non-white space charaters in code:

      if(myCrossPropagationCode.isAllWhiteSpace()) {
	throw xmdsException(crossPropagationElement,"No <cross_propagation> code defined!");
      }

      if(verbose()) {
	printf("cross_propagation code loaded\n");
      }
    }
  }
 

  // ************************************
  // process elements in order to process and place the code elements

  myNumNonLoopPropagation=0;
  myNumIntegrateMomentGroups=0;
  
  long unsigned int i;
  const NodeList* myElements;
  const Element* nextElement;

  myElements = yourElement->getElementsByTagName("*",0);
  i=0;

  while(i<myElements->length()) {
    nextElement = dynamic_cast<const Element*> (myElements->item(i));
    if(*nextElement->nodeName() == "functions") {
      XMLString someLoopPropagationCode = *nextElement->textContent(0);
      // check for non-white space charaters in code:
      if(someLoopPropagationCode.isAllWhiteSpace()) {
        throw xmdsException(nextElement,"No <functions> code defined!");
      }
      if(verbose()) {
        printf("    <functions> code loaded\n");
      }
      myNumNonLoopPropagation++;
      myNonLoopPropagationCodeList.push_back(someLoopPropagationCode);
      myCodeElementList.push_back(*nextElement->nodeName());
    }
    else if(*nextElement->nodeName() == "moment_group") {
      integrateMomentGroup tempIntegrateMG;
      
      XMLString someCode = *nextElement->textContent(0);
      // check for non-white space charaters in code:
      if(someCode.isAllWhiteSpace()) {
        throw xmdsException(nextElement,"No <moment_group> code defined!");
      }
      if(verbose()) {
        printf("   <moment_group> code loaded\n");
      }
      tempIntegrateMG.integrateMomentGroupCode = someCode;

      getAssignmentStrings(nextElement,"moments",1,0,tempIntegrateMG.momentNameList);      
      
      getAssignmentBools(nextElement,"integrate_dimension",1,simulation()->field()->geometry()->nDims(),tempIntegrateMG.integrateDimensionList);
      
      myNumIntegrateMomentGroups++;
      myIntegrateMomentGroupList.push_back(tempIntegrateMG);
      myCodeElementList.push_back(*nextElement->nodeName());
    }
    else if(*nextElement->nodeName() == "vectors") {
      myCodeElementList.push_back(*nextElement->nodeName());
    }

    i++;
  }
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrate protected
// ******************************************************************************
// ******************************************************************************
// ******************************************************************************

const bool xmdsIntegrate::Smallmemory() const {
  if(debugFlag) {
    printf("xmdsIntegrate::Smallmemory\n");
  }

  return mySmallmemory;
};

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

  if(usesKOperators()) {
    if((myAdaptiveIP && !Smallmemory()) || (!myAdaptiveIP && constantK())) {
      fprintf(outfile,"// integrate defines\n");
      fprintf(outfile,"\n");
      fprintf(outfile,"#define _segment%li_nkoperators %li\n",segmentNumber,(long)myKOperatorNamesList.size());
      fprintf(outfile,"\n");
    }
  }
};

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

if(usesKOperators()) {
 if(myAdaptiveIP){
   if(usesKOperators()) 
    if(!Smallmemory()) {
      fprintf(outfile,"// integrate globals\n");
      fprintf(outfile,"\n");
      
      for(int n=1;n<6;n++)
      if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	fprintf(outfile,"complex *_segment%li_k_operator_field_%i;\n",segmentNumber,n);
      }
      else {
	fprintf(outfile,"complex *_segment%li_k_operator_field_%i = new complex[_%s_size*_segment%li_nkoperators];\n",
		segmentNumber,n,simulation()->field()->name()->c_str(),segmentNumber);
      }
      fprintf(outfile,"\n");
  }
       
}else{
    if(constantK()) {
      fprintf(outfile,"// integrate globals\n");
      fprintf(outfile,"\n");
      if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	fprintf(outfile,"complex *_segment%li_k_operator_field;\n",segmentNumber);
      }
      else {
	fprintf(outfile,"complex *_segment%li_k_operator_field = new complex[_%s_size*_segment%li_nkoperators];\n",
		segmentNumber,simulation()->field()->name()->c_str(),segmentNumber);
      }
      fprintf(outfile,"\n");
    }
   } 
  }
};

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

  writexSpacePrototype(outfile);
};

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

  writexSpaceRoutine(outfile);
};

// ******************************************************************************
const XMLString* xmdsIntegrate::interval() const {
  if(debugFlag) {
    printf("xmdsIntegrate::interval\n");
  }

  return &myInterval;
};

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

  return myLattice;
};

// ******************************************************************************
unsigned long xmdsIntegrate::samples(
				     const unsigned long& index) const {
  if(debugFlag) {
    printf("xmdsIntegrate::KVectorNamesList\n");
  }

  if(index>=mySamplesList.size()) {
    throw xmdsException("Internal range error in xmdsIntegrate::samples");
  }

  list<unsigned long>::const_iterator pULong = mySamplesList.begin();
  for(unsigned long i=0; i<index; i++) {
    pULong++;
  }
  return *pULong;
};

// ******************************************************************************
const list<XMLString>* xmdsIntegrate::KVectorNamesList() const {
  if(debugFlag) {
    printf("xmdsIntegrate::KVectorNamesList\n");
  }

  return &myKVectorNamesList;
};

// ******************************************************************************
bool xmdsIntegrate::constantK() const {
  if(debugFlag) {
    printf("xmdsIntegrate::constantK\n");
  }

  return myConstantK;
};

// ******************************************************************************
bool xmdsIntegrate::noNoises() const {
  if(debugFlag) {
    printf("xmdsIntegrate::noNoises\n");
  }

  return myNonoises;
};

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

  return myKOperatorNamesList.size();
};

// ******************************************************************************
const XMLString* xmdsIntegrate::KOperator(
					  const unsigned long& index) const {
  if(debugFlag) {
    printf("xmdsIntegrate::KOperator\n");
  }

  if(index>=myKOperatorNamesList.size()) {
    throw xmdsException("Internal range error in xmdsIntegrate::KOperator");
  }

  list<XMLString>::const_iterator pXMLString = myKOperatorNamesList.begin();
  for(unsigned long i=0; i<index; i++) {
    pXMLString++;
  }

  return &*pXMLString;
};

// ******************************************************************************
bool xmdsIntegrate::getKOperator(
				 const XMLString& operatorName,
				 unsigned long& index) const {
  if(debugFlag) {
    printf("xmdsField::getVector\n");
  }

  index=0;
  for(list<XMLString>::const_iterator pXMLString = myKOperatorNamesList.begin(); pXMLString != myKOperatorNamesList.end(); pXMLString++) {
    if(*pXMLString==operatorName) {
      return 1;
    }
    index++;
  }

  return 0;
};

// ******************************************************************************
const XMLString* xmdsIntegrate::KOperatorsCode() const {
  if(debugFlag) {
    printf("xmdsIntegrate::KOperatorsCode\n");
  }

  return &myKOperatorsCode;
};

// ******************************************************************************
list<XMLString>* xmdsIntegrate::vectorNamesList() {
  if(debugFlag) {
    printf("xmdsIntegrate::vectorNamesList\n");
  }

  return &myVectorNamesList;
};

// ******************************************************************************
const list<XMLString>* xmdsIntegrate::vectorNamesList() const {
  if(debugFlag) {
    printf("xmdsIntegrate::vectorNamesList\n");
  }

  return &myVectorNamesList;
};

// ******************************************************************************
XMLString* xmdsIntegrate::propagationCode() {
  if(debugFlag) {
    printf("xmdsIntegrate::propagationCode\n");
  }

  return &myPropagationCode;
};

// ******************************************************************************
const XMLString* xmdsIntegrate::propagationCode() const {
  if(debugFlag) {
    printf("xmdsIntegrate::propagationCode\n");
  }

  return &myPropagationCode;
};

// ******************************************************************************
const list<XMLString>* xmdsIntegrate::crossVectorNamesList() const {
  if(debugFlag) {
    printf("xmdsIntegrate::crossVectorNamesList\n");
  }

  return &myCrossVectorNamesList;
};

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

  return myCrossDimNumber;
};

// ******************************************************************************
const XMLString* xmdsIntegrate::crossPropagationCode() const {
  if(debugFlag) {
    printf("xmdsIntegrate::crossPropagationCode\n");
  }

  return &myCrossPropagationCode;
};

// ******************************************************************************
list<XMLString>* xmdsIntegrate::nonLoopPropagationCodeList() {
  if(debugFlag) {
    printf("xmdsIntegrate::nonLoopPropagationCodeList\n");
  }

  return &myNonLoopPropagationCodeList;
};

// ******************************************************************************
const list<XMLString>* xmdsIntegrate::nonLoopPropagationCodeList() const {
  if(debugFlag) {
    printf("xmdsIntegrate::nonLoopPropagationCodeList\n");
  }

  return &myNonLoopPropagationCodeList;
};

// ******************************************************************************
list<XMLString>* xmdsIntegrate::codeElementList() {
  if(debugFlag) {
    printf("xmdsIntegrate::codeElementList\n");
  }

  return &myCodeElementList;
};

// ******************************************************************************
const list<XMLString>* xmdsIntegrate::codeElementList() const {
  if(debugFlag) {
    printf("xmdsIntegrate::codeElementList\n");
  }

  return &myCodeElementList;
};

// ******************************************************************************
bool xmdsIntegrate::usesKOperators() const {
  if(debugFlag) {
    printf("xmdsIntegrate::usesKOperators\n");
  }

  return (mycoStructList.size()>0);
};

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

const bool xmdsIntegrate::AdaptiveIP() const {
  if(debugFlag) {
    printf("xmdsIntegrate::AdaptiveIP\n");
  }

  return myAdaptiveIP;
};

// ******************************************************************************
long unsigned int xmdsIntegrate::numNonLoopPropagation() const {
  if(debugFlag) {
    printf("xmdsIntegrate::numNonLoopPropagation\n");
  }

  return (myNumNonLoopPropagation);
};

// ******************************************************************************
long unsigned int xmdsIntegrate::numIntegrateMomentGroups() const {
  if(debugFlag) {
    printf("xmdsIntegrate::numIntegrateMomentGroups\n");
  }

  return (myNumIntegrateMomentGroups);
};

// ******************************************************************************
list<integrateMomentGroup>* xmdsIntegrate::integrateMomentGroupList() {
  if(debugFlag) {
    printf("xmdsIntegrate::integrateMomentGroupList\n");
  }

  return &myIntegrateMomentGroupList;
};

// ******************************************************************************
const list<integrateMomentGroup>* xmdsIntegrate::integrateMomentGroupList() const {
  if(debugFlag) {
    printf("xmdsIntegrate::integrateMomentGroupList\n");
  }

  return &myIntegrateMomentGroupList;
};

// ******************************************************************************
bool xmdsIntegrate::getcoStruct(
				const unsigned long& componentNumber,
				const coStruct*& thecoStruct) const {
  if(debugFlag) {
    printf("xmdsIntegrate::getcoStruct\n");
  }

  for(list<coStruct>::const_iterator pcoStruct = mycoStructList.begin(); pcoStruct != mycoStructList.end(); pcoStruct++) {
    if(pcoStruct->componentNumber == componentNumber) {
      thecoStruct = &*pcoStruct;
      return 1;
    }
  }
  return 0;
};

// ******************************************************************************
bool xmdsIntegrate::getcoKey(
			     const unsigned long& componentNumber,
			     const unsigned long& operatorNumber,
			     unsigned long& coKey) const {
  if(debugFlag) {
    printf("xmdsIntegrate::getcoPairKey\n");
  }

  const coStruct* thecoStruct;

  if(!getcoStruct(componentNumber,thecoStruct)) {
    return 0;
  }

  list<unsigned long>::const_iterator pULong1 = thecoStruct->operatorNumbersList.begin();
  list<unsigned long>::const_iterator pULong2 = thecoStruct->coKeysList.begin();

  while(pULong1 != thecoStruct->operatorNumbersList.end()) {
    if(*pULong1==operatorNumber) {
      coKey = *pULong2;
      return 1;
    }
    pULong1++;
    pULong2++;
  }

  return 0;
};

// ******************************************************************************
unsigned long xmdsIntegrate::addcoPair(
				       const unsigned long& componentNumber,
				       const unsigned long& operatorNumber) {
  if(debugFlag) {
    printf("xmdsIntegrate::addcoPair\n");
  }

  coStruct* existingcoStruct = 0;

  for(list<coStruct>::iterator pcoStruct = mycoStructList.begin(); pcoStruct != mycoStructList.end(); pcoStruct++) {
    if(pcoStruct->componentNumber == componentNumber) {
      existingcoStruct = &*pcoStruct;
    }
  }

  if(existingcoStruct != 0) {
    existingcoStruct->operatorNumbersList.push_back(operatorNumber);
    existingcoStruct->coKeysList.push_back(myNextcoKey);
  }
  else {
    coStruct nextcoStruct;
    nextcoStruct.componentNumber = componentNumber;
    nextcoStruct.operatorNumbersList.push_back(operatorNumber);
    nextcoStruct.coKeysList.push_back(myNextcoKey);
    mycoStructList.push_back(nextcoStruct);
  }

  myNextcoKey++;
  return myNextcoKey-1;
};

// ******************************************************************************
bool xmdsIntegrate::findNextcoPair(
				   XMLString& operatorName,
				   XMLString& componentName,
				   unsigned long& start,
				   unsigned long& end) const {
  if(debugFlag) {
    printf("xmdsIntegrate::findNextcoPair\n");
  }

  // here we are looking for the next **** operatorName S? '[' S? componentName S? ']' **** production
  // starting from the offset provided

  const char* code=myPropagationCode.c_str();
  end = start;
  unsigned long openBrace;
  unsigned long nameStart;
  unsigned long nameEnd;

  // find next [

  while((code[end] != '[') && (code[end] != 0)) {
    end++;
  }

  if(code[end] == 0) {
    return 0;
  }

  // ok, we found the first [
  openBrace = end;

  // now skip white space
  end++;
  while(XMLChar::isWhiteSpace(code[end]) && (code[end] != 0)) {
    end++;
  }

  if(code[end] == 0) {
    return 0;
  }

  if(!(XMLChar::isLatinLetter(code[end]) | (code[end]=='_'))) {
    return 0;
  }

  nameStart = end;

  end++;
  while((XMLChar::isLatinLetter(code[end]) | XMLChar::isLatinDigit(code[end]) | (code[end]=='_'))&(code[end] != 0)) {
    end++;
  }

  if(code[end] == 0) {
    return 0;
  }

  nameEnd = end;

  //check that we have a ] on the end

  // skip white space
  while(XMLChar::isWhiteSpace(code[end]) && (code[end] != 0)) {
    end++;
  }

  if(code[end] != ']') {
    return 0;
  }

  // ok, we definititely had a nice clean [name]
  end++;

  myPropagationCode.subString(componentName,nameStart,nameEnd);

  // now work backwards to find the operator name
  start=openBrace-1;

  // skip any white space
  while(XMLChar::isWhiteSpace(code[start])&(start>0)) {
    start--;
  }

  if(!(XMLChar::isLatinLetter(code[start]) | XMLChar::isLatinDigit(code[start]) | (code[start]=='_'))) {
    return 0;
  }

  nameEnd=start+1;

  while((XMLChar::isLatinLetter(code[start]) | XMLChar::isLatinDigit(code[start]) | (code[start]=='_'))&(start>0)) {
    start--;
  }

  if(!(XMLChar::isLatinLetter(code[start]) | XMLChar::isLatinDigit(code[start]) | (code[start]=='_'))) {
    start++;
  }

  if(!(XMLChar::isLatinLetter(code[start]) | (code[start]=='_'))) {
    return 0;
  }

  nameStart=start;

  myPropagationCode.subString(operatorName,nameStart,nameEnd);

  return 1;
};

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

  fprintf(outfile,"void _segment%li_calculate_delta_a(\n",segmentNumber);
  fprintf(outfile,"	const double& _step");
  fprintf(outfile,",\n	const unsigned long cycle");
  if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
    fprintf(outfile,",\n	const double *const _noise_vector");
  }
  fprintf(outfile,");\n");
  fprintf(outfile,"\n");
};

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

  const unsigned long nDims = simulation()->field()->geometry()->nDims();
  const char *const fieldName = simulation()->field()->name()->c_str();

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li_calculate_delta_a(\n",segmentNumber);
  fprintf(outfile,"	const double& _step");
  fprintf(outfile,",\n	const unsigned long cycle");
  if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
    fprintf(outfile,",\n	const double *const _noise_vector");
  }
  fprintf(outfile,") {\n");
  fprintf(outfile,"\n");

  const xmdsVector* mainVector;

  if(!simulation()->field()->getVector("main",mainVector)) {
    throw xmdsException("Internal error in xmdsIntegrate::writexSpaceRoutine: cannot find 'main' vector");
  }

  const char* typeName="";
  if(mainVector->vectorType()==COMPLEX) {
    typeName="complex";
  }
  else if(mainVector->vectorType()==DOUBLE) {
    typeName="double";
  }

  for(unsigned long i=0;i<mainVector->nComponents();i++) {
    fprintf(outfile,"%s d%s_d%s;\n",typeName,mainVector->componentName(i)->c_str(),simulation()->parameters()->propDimName.c_str());
  }
  fprintf(outfile,"\n");

  // add cross vectors to total vectors to use, but no duplicating names!

  list<XMLString> myTotalVectorsList = *vectorNamesList();

  for(list<XMLString>::const_iterator pXMLString = crossVectorNamesList()->begin(); pXMLString != crossVectorNamesList()->end(); pXMLString++) {
    list<XMLString>::const_iterator pXMLString2 = myTotalVectorsList.begin();
    while((pXMLString2 != myTotalVectorsList.end()) && (*pXMLString2 != *pXMLString)) {
      pXMLString2++;
    }
    if(*pXMLString2 != *pXMLString) {
      myTotalVectorsList.push_back(*pXMLString);
    }
  }

  simulation()->field()->vectors2space(outfile,0,myTotalVectorsList,"");

  if(crossVectorNamesList()->size() > 0) {
    fprintf(outfile,"_segment%li_calculate_cross_field();\n",segmentNumber);
    fprintf(outfile,"\n");
  }

  if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
    fprintf(outfile,"const double* _noises = _noise_vector;\n");
    fprintf(outfile,"\n");
  }


  list<XMLString>::const_iterator nextNLPElement = myNonLoopPropagationCodeList.begin();
  list<integrateMomentGroup>::const_iterator nextMGElement = myIntegrateMomentGroupList.begin();

  list<XMLString> theCodeList = *codeElementList();
  if(!(theCodeList.size()==1+numIntegrateMomentGroups()+numNonLoopPropagation())) {
    throw xmdsException("The list of integrate code elements is the wrong length!");
  }
    
  long unsigned int whichMG = 0;

  list<long> deleteMGArrayList;

  for(list<XMLString>::const_iterator codeElement = theCodeList.begin(); codeElement != theCodeList.end(); codeElement++) {
    if(!strcmp(codeElement->c_str(),"vectors")) {
      simulation()->field()->openLoops(outfile,0,myTotalVectorsList);
        
      fprintf(outfile,"\n");
        
      char indent[64];
      for(unsigned long i=0;i<nDims;i++) {
	indent[i]=0x09;
      }
      indent[nDims]=0;

      // integrate moment group pointer definition

      long unsigned int tempWhichMG = 0;
      for(list<integrateMomentGroup>::const_iterator tempNextMGElement = myIntegrateMomentGroupList.begin(); 
	  tempNextMGElement != myIntegrateMomentGroupList.end(); tempNextMGElement++) {

	long unsigned int nDimsRemaining = simulation()->field()->geometry()->nDims();
	tempWhichMG++;
			  
	for(list<bool>::const_iterator nextIntegrateDimension = tempNextMGElement->integrateDimensionList.begin(); 
	    nextIntegrateDimension != tempNextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
	  if(*nextIntegrateDimension) {
	    nDimsRemaining--;
	  }
	}
				
	if(nDimsRemaining>0) {
	  fprintf(outfile,"        long _img%li_pointer = ",tempWhichMG);
	  for(long unsigned int j=0; j<nDimsRemaining-1; j++) {
	    fprintf(outfile,"(");
	  }
	  long unsigned int whichMoment = 0;
	  list<bool>::const_iterator nextIntegrateDimension = tempNextMGElement->integrateDimensionList.begin();

	  while(*nextIntegrateDimension) {
	    nextIntegrateDimension++;
	    whichMoment++;
	  }
	  fprintf(outfile,"_i%li",whichMoment);  
                  
	  nextIntegrateDimension++;
	  whichMoment++;
	  for(long unsigned int k=0; k<nDimsRemaining-1; k++) {
	    while(*nextIntegrateDimension) {
	      nextIntegrateDimension++;
	      whichMoment++;
	    }
	    fprintf(outfile,"*_%s_lattice%li+_i%li)",fieldName,whichMoment,whichMoment);
	    nextIntegrateDimension++;
	    whichMoment++;
	  }

	  fprintf(outfile,"*%i;\n",tempNextMGElement->momentNameList.size());
	}
		
      }
        
      fprintf(outfile,"// ************** propagation code **************\n");
      fprintf(outfile,"%s\n",propagationCode()->c_str());
      fprintf(outfile,"// **********************************************\n");
      fprintf(outfile,"\n");
      for(long unsigned int i=0;i<mainVector->nComponents();i++) {
	fprintf(outfile,"%s_active_%s_main[_%s_main_index_pointer + %li] = d%s_d%s*_step;\n",
		indent,fieldName,fieldName,i,mainVector->componentName(i)->c_str(),simulation()->parameters()->propDimName.c_str());
      }
      fprintf(outfile,"\n");
        
      if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
	fprintf(outfile,"%s_noises += _n_noises;\n",indent);
	fprintf(outfile,"\n");
      }
        
      simulation()->field()->closeLoops(outfile,0,myTotalVectorsList);
      fprintf(outfile,"\n");
    }
    else if(!strcmp(codeElement->c_str(),"functions")) {       
      fprintf(outfile,"// ************** propagation code **************\n");
      fprintf(outfile,"%s\n",nextNLPElement->c_str());
      fprintf(outfile,"// **********************************************\n");
      fprintf(outfile,"\n");
      nextNLPElement++;
    }
    else if(!strcmp(codeElement->c_str(),"moment_group")) {  
      whichMG++;  
              
      long unsigned int nDimsRemaining = simulation()->field()->geometry()->nDims();

      for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	  nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
	if(*nextIntegrateDimension) {
	  nDimsRemaining--;
	}
      }

      if(nDimsRemaining>0) {
	deleteMGArrayList.push_back(whichMG);
      }

      fprintf(outfile,"// ************** integrate moment group code **************\n");
                
      //setup defines for the integrate_moment_groups if they are arrays
      long unsigned int whichMoment=0;
      if(nDimsRemaining>0) {
	for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin(); 
	    nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++){
	  fprintf(outfile,"#define %s _img%li_array[_img%li_pointer+%li]\n",nextMomentName->c_str(),whichMG,whichMG,whichMoment);
	  whichMoment++;
	}
      }
      fprintf(outfile,"\n");

      //setup actual arrays for the integrate_moment_groups (MPI and non-MPI)
      if(nDimsRemaining>0) {
	fprintf(outfile,"complex *_img%li_array = new complex[%i",whichMG,nextMGElement->momentNameList.size());
	if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	  whichMoment=0;
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	    if(!*nextIntegrateDimension) {
	      if(whichMoment==0){
		fprintf(outfile,"*local_nx");
	      }
	      else {
		fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	      }
	    }
	    whichMoment++;
	  }
	}
	else {
	  whichMoment=0;
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	    if(!*nextIntegrateDimension) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	    }
	    whichMoment++;
	  }
	}
	fprintf(outfile,"];\n");
      }
      else {
	for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin(); 
	    nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++) {
	  fprintf(outfile,"complex %s = rcomplex(0,0);\n",nextMomentName->c_str());
	}
      }
      fprintf(outfile,"\n");
              
      //zero the arrays and setup a local environment for the sum
      if(nDimsRemaining>0) {
	fprintf(outfile,"for(long _i0=0; _i0<%i",nextMGElement->momentNameList.size());
	if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	  whichMoment=0;
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	    if(!*nextIntegrateDimension) {
	      if(whichMoment==0){
		fprintf(outfile,"*local_nx");
	      }
	      else {
		fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	      }
	    }
	    whichMoment++;
	  }
	}
	else {
	  whichMoment=0;
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	    if(!*nextIntegrateDimension) {
	      fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	    }
	    whichMoment++;
	  }
	}
	fprintf(outfile,"; _i0++) {\n");
	fprintf(outfile,"  _img%li_array[_i0] = rcomplex(0,0);\n",whichMG);
	fprintf(outfile,"}\n");
      }

      fprintf(outfile,"{\n");
      // open loops
      simulation()->field()->openLoops(outfile,0,myTotalVectorsList);
      // pointer definition
      if(nDimsRemaining>0) {
	fprintf(outfile,"        long _img%li_pointer = ",whichMG);
	for(long unsigned int j=0; j<nDimsRemaining-1; j++) {
	  fprintf(outfile,"(");
	}
	whichMoment = 0;
	list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin();

	while(*nextIntegrateDimension) {
	  nextIntegrateDimension++;
	  whichMoment++;
	}
	fprintf(outfile,"_i%li",whichMoment);  
                  
	nextIntegrateDimension++;
	whichMoment++;
	for(long unsigned int k=0; k<nDimsRemaining-1; k++) {
	  while(*nextIntegrateDimension) {
	    nextIntegrateDimension++;
	    whichMoment++;
	  }
	  fprintf(outfile,"*_%s_lattice%li+_i%li)",fieldName,whichMoment,whichMoment);
	  nextIntegrateDimension++;
	  whichMoment++;
	}

	fprintf(outfile,"*%i;\n",nextMGElement->momentNameList.size());
      }
      // code
      fprintf(outfile,"%s\n",nextMGElement->integrateMomentGroupCode.c_str());

      // close loops
      simulation()->field()->closeLoops(outfile,0,myTotalVectorsList);
			  
      fprintf(outfile,"}\n");

      // if there was any integration, multiply by the relevant volume element
      if(nDimsRemaining<simulation()->field()->geometry()->nDims()) {
	if(nDimsRemaining>0) {
	  fprintf(outfile,"for(long _i0=0; _i0<%i",nextMGElement->momentNameList.size());
	  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	    whichMoment=0;
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	      if(!*nextIntegrateDimension) {
		if(whichMoment==0){
		  fprintf(outfile,"*local_nx");
		}
		else {
		  fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
		}
	      }
	      whichMoment++;
	    }
	  }
	  else {
	    whichMoment=0;
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	      if(!*nextIntegrateDimension) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	      }
	      whichMoment++;
	    }
	  }
	  whichMoment=0;
	  fprintf(outfile,"; _i0++) {\n");
	  fprintf(outfile,"   _img%li_array[_i0] = _img%li_array[_i0]",whichMG,whichMG);
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
	    if(*nextIntegrateDimension) {
	      fprintf(outfile,"*_%s_dx%li",fieldName,whichMoment);
	    }
	    whichMoment++;
	  }
	  fprintf(outfile,";\n");
	  fprintf(outfile," }\n");
	}
	else {
	  for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin(); 
	      nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++) {
	    whichMoment = 0;
	    fprintf(outfile,"  %s = %s",nextMomentName->c_str(),nextMomentName->c_str());
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
	      if(*nextIntegrateDimension) {
		fprintf(outfile,"*_%s_dx%li",fieldName,whichMoment);
	      }
	      whichMoment++;
	    }
	    fprintf(outfile,";\n");
	  }
	}
      }

      // if mpi and first variable integrated then mpireduce
      if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	if(* nextMGElement->integrateDimensionList.begin()) {
	  if(nDimsRemaining>0) {
	    fprintf(outfile,"\n  complex *_temp_img%li_array = new complex[%i",whichMG,nextMGElement->momentNameList.size());
	    whichMoment=0;
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	      if(!*nextIntegrateDimension) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	      }
	      whichMoment++;
	    }
	    fprintf(outfile,"];\n");
	    fprintf(outfile,"  MPI_Reduce(_img%li_array,_temp_img%li_array,2*%i",whichMG,whichMG,nextMGElement->momentNameList.size());
	    whichMoment=0;
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++) {
	      if(!*nextIntegrateDimension) {
		fprintf(outfile,"*_%s_lattice%li",fieldName,whichMoment);
	      }
	      whichMoment++;
	    }
	    fprintf(outfile,",MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n");
	    fprintf(outfile,"\n  delete[] _temp_img%li_array;\n",whichMG);					
	  }
	  else {
	    fprintf(outfile,"\n  complex *_temp%li_mpireduce = new complex;\n",whichMG);
	    for(list<XMLString>::const_iterator nextMomentName = nextMGElement->momentNameList.begin(); 
		nextMomentName != nextMGElement->momentNameList.end(); nextMomentName++) {
	      fprintf(outfile,"  MPI_Reduce(&%s,_temp%li_mpireduce,2,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);\n",nextMomentName->c_str(),whichMG);
	      fprintf(outfile,"\n  delete _temp%li_mpireduce;\n",whichMG);
	    }
	  }
	}
      }

      fprintf(outfile,"\n");
      fprintf(outfile,"// *********************************************************\n");
      fprintf(outfile,"\n");
      nextMGElement++;           
    }
    else {
      throw xmdsException("Unknown code element in the integrate block!");
    }
  }

  // delete any arrays that were generated
  for(list<long>::const_iterator deleteMGArray = deleteMGArrayList.begin(); 
      deleteMGArray != deleteMGArrayList.end(); deleteMGArray++) {
    fprintf(outfile,"delete[] _img%li_array;\n",*deleteMGArray);					
  }
  
  fprintf(outfile,"}\n");
  fprintf(outfile,"\n");
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsIntegrate private
// ******************************************************************************
// ******************************************************************************

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

  simulation()->output()->addSamples(mySamplesList);
};

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

  if(usesKOperators()&&!Smallmemory()) {
   if(myAdaptiveIP){   
      for(int i=1;i<6;i++)
       if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
         fprintf(outfile,"_segment%li_k_operator_field_%i = new complex[total_local_size*_segment%li_nkoperators];\n",
	      segmentNumber,i,segmentNumber);
         fprintf(outfile,"\n");
       }    
       }else{
     if(constantK()) 
       fprintf(outfile,"_segment%li_calculate_k_operator_field();\n",segmentNumber);
    }
  }
};