/*
  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: xmdsfilter.cc,v 1.18 2004/10/21 09:42:03 paultcochrane Exp $
*/

/*! @file xmdsfilter.cc
  @brief Filter element parsing classes and methods

  More detailed explanation...
*/

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

// ******************************************************************************
// ******************************************************************************
//                              xmdsFilter public
// ******************************************************************************
// ******************************************************************************

extern bool debugFlag;

long nxmdsFilters=0;  //!< Number of xmds filter objects

// ******************************************************************************
xmdsFilter::xmdsFilter(
		       const xmdsSimulation *const yourSimulation,
		       const bool& yourVerboseMode) :
  xmdsSegment(yourSimulation,yourVerboseMode) {
  if(debugFlag) {
    nxmdsFilters++;
    printf("xmdsFilter::xmdsFilter\n");
    printf("nxmdsFilters=%li\n",nxmdsFilters);
  }
};

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

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

  const long unsigned int nDims = simulation()->field()->geometry()->nDims();

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

  if(nDims>0) {

    // ************************************
    // find space

    list<bool> aSpaceList;

    getAssignmentBools(yourElement,"fourier_space",0,simulation()->field()->geometry()->nDims(),aSpaceList);

    if(aSpaceList.size() == 0) {
      printf("Filter space for defaulting to x-space.\n");
      mySpace = 0;
    }
    else {
      list<bool>::const_iterator pBool = aSpaceList.begin();
      for(long unsigned int i=0;i<simulation()->field()->geometry()->nDims();i++) {
	if(verbose()) {
	  if(*pBool) {
	    printf("Filter will be performed with dimension #%li in fourier space\n",i+1);
	  }
	  else {
	    printf("Filter will be performed with dimension #%li in normal space\n",i+1);
	  }
	}
	pBool++;
      }
      mySpace = spaceList2ULong(aSpaceList);
    }
  }

  // ************************************
  // 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");
        }
  }

  // ************************************
  // 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,mySpace);

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

  myCode=*yourElement->textContent(0);

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

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

  if(verbose()) {
    printf("filter 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++;
  }
};

// ******************************************************************************
// ******************************************************************************
//                              xmdsFilter protected
// ******************************************************************************
// ******************************************************************************

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

  return &myCodeElementList;
};

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

  return &myCodeElementList;
};
// ******************************************************************************
list<XMLString>* xmdsFilter::nonLoopPropagationCodeList() {
  if(debugFlag) {
    printf("xmdsFilter::nonLoopPropagationCodeList\n");
  }

  return &myNonLoopPropagationCodeList;
};

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

  return &myNonLoopPropagationCodeList;
};
// ******************************************************************************
long unsigned int xmdsFilter::numNonLoopPropagation() const {
  if(debugFlag) {
    printf("xmdsFilter::numNonLoopPropagation\n");
  }

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

  return (myNumIntegrateMomentGroups);
};

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

  return &myIntegrateMomentGroupList;
};

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

  return &myIntegrateMomentGroupList;
};



// ******************************************************************************
// ******************************************************************************
//                              xmdsFilter private
// ******************************************************************************
// ******************************************************************************

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

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

  fprintf(outfile,"void _segment%li(unsigned long cycle); // filter\n",segmentNumber);
};

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

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

  long unsigned int i;

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

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

  fprintf(outfile,"// *************************\n");
  fprintf(outfile,"void _segment%li(unsigned long cycle) {\n",segmentNumber);
  fprintf(outfile,"\n");

  if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
    fprintf(outfile,"const double _var = 1.0");
    for(i=0;i<nDims;i++) {
      if(space(i)) {
	fprintf(outfile,"/_%s_dk%li",fieldName,i);
      }
      else {
	fprintf(outfile,"/_%s_dx%li",fieldName,i);
      }
    }
    fprintf(outfile,";\n");
    fprintf(outfile,"double _noises[_n_noises];\n");
    if(simulation()->parameters()->errorCheck) {
      fprintf(outfile,"double _noises2[_n_noises];\n");
    }
    fprintf(outfile,"\n");
  }

  simulation()->field()->vectors2space(outfile,mySpace,myVectorNamesList,"");
  
  bool swapped=false;
  if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
    if((mySpace&1)&((mySpace>>1)&1)) {
      swapped=true;
    }
  }

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

  list<XMLString> theCodeList = *codeElementList();

  long unsigned int filterCodeElemLen = 1+numIntegrateMomentGroups()+numNonLoopPropagation();
  if(!(theCodeList.size()==filterCodeElemLen)) {
    throw xmdsException("The list of filter 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,mySpace,myVectorNamesList);
		
      char indent[64];
      for(i=0;i<nDims;i++) {
	indent[i]=0x09;
      }
      indent[nDims]=0;
		
      if((simulation()->parameters()->stochastic)&&(!myNonoises)) {
	if(simulation()->parameters()->errorCheck) {
	  fprintf(outfile,"%s_make_noises(_gen1,_var/2,_noises,_n_noises);\n",indent);
	  fprintf(outfile,"%s_make_noises(_gen2,_var/2,_noises2,_n_noises);\n",indent);
	  fprintf(outfile,"%sfor(unsigned long _s0=0;_s0<_n_noises;_s0++)\n",indent);
	  fprintf(outfile,"%s	_noises[_s0] += _noises2[_s0];\n",indent);
	  fprintf(outfile,"\n");
	}
	else {
	  fprintf(outfile,"%s_make_noises(_gen,_var,_noises,_n_noises);\n",indent);
	  fprintf(outfile,"\n");
	}
      }
      // 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++;
	    }
	    if(swapped&&(whichMoment==1)) {
	      fprintf(outfile,"*_%s_lattice0+_i1)",fieldName);
	    }
	    else {
	      fprintf(outfile,"*_%s_lattice%li+_i%li)",fieldName,whichMoment,whichMoment);
	    }
	    nextIntegrateDimension++;
	    whichMoment++;
	  }
				
	  fprintf(outfile,"*%i;\n",tempNextMGElement->momentNameList.size());
	}
			
      }
		
      fprintf(outfile,"// *************** Filter code ******************\n");
      fprintf(outfile,"%s\n",myCode.c_str());
      fprintf(outfile,"// **********************************************\n");
				
      simulation()->field()->closeLoops(outfile,mySpace,myVectorNamesList);
				
      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) {
		if(swapped) {
		  fprintf(outfile,"*_%s_lattice0",fieldName);
		}
		else {
		  fprintf(outfile,"*local_nx");
		}
	      }
	      else if(whichMoment==1) {
		if(!swapped) {
		  fprintf(outfile,"*_%s_lattice1",fieldName);
		}
		else {
		  fprintf(outfile,"*local_ny_after_transpose");
		}
	      }
	      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) {
		if(swapped) {
		  fprintf(outfile,"*_%s_lattice0",fieldName);
		}
		else {
		  fprintf(outfile,"*local_nx");
		}
	      }
	      else if(whichMoment==1) {
		if(!swapped) {
		  fprintf(outfile,"*_%s_lattice1",fieldName);
		}
		else {
		  fprintf(outfile,"*local_ny_after_transpose");
		}
	      }
	      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,mySpace,myVectorNamesList);
      // 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++;
	  }
	  if(swapped&&(whichMoment==1)) {
	    fprintf(outfile,"*_%s_lattice0+_i1)",fieldName);
	  }
	  else {
	    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,mySpace,myVectorNamesList);
			  
      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,"\nfor(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) {
		  if(swapped) {
		    fprintf(outfile,"*_%s_lattice0",fieldName);
		  }
		  else {
		    fprintf(outfile,"*local_nx");
		  }
		}
		else if(whichMoment==1) {
		  if(!swapped) {
		    fprintf(outfile,"*_%s_lattice1",fieldName);
		  }
		  else {
		    fprintf(outfile,"*local_ny_after_transpose");
		  }
		}
		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);
	  long unsigned int tempSpace = mySpace;
	  for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
	      nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
	    if(*nextIntegrateDimension) {
	      if(tempSpace&1) {
		fprintf(outfile,"*_%s_dk%li",fieldName,whichMoment);
	      }
	      else {
		fprintf(outfile,"*_%s_dx%li",fieldName,whichMoment);
	      }
	    }
	    tempSpace >>= 1;
	    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());
	    long unsigned int tempSpace = mySpace;
	    for(list<bool>::const_iterator nextIntegrateDimension = nextMGElement->integrateDimensionList.begin(); 
		nextIntegrateDimension != nextMGElement->integrateDimensionList.end(); nextIntegrateDimension++){
	      if(*nextIntegrateDimension) {
		if(tempSpace&1) {
		  fprintf(outfile,"*_%s_dk%li",fieldName,whichMoment);
		}
		else {
		  fprintf(outfile,"*_%s_dx%li",fieldName,whichMoment);
		}
	      }
	      tempSpace >>= 1;
	      whichMoment++;
	    }
	    fprintf(outfile,";\n");
	  }
	}
      }

      // if mpi and first variable integrated then mpireduce
      if(simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
	if(swapped) {
	  // have to predicate on the second field at least
	  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);
	      }
	    }
	  }
	}
	else {
	  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");
};

// ******************************************************************************
//! Returns space at index.  Don't ask me why this returns boolean
bool xmdsFilter::space(
		       const long unsigned int& index) const {
  if(debugFlag) {
    printf("xmdsFilter::space\n");
  }

  if(index>=simulation()->field()->geometry()->nDims()) {
    throw xmdsException("Internal range error in xmdsFilter::space()");
  }

  return (mySpace >> index)&1;
};


syntax highlighted by Code2HTML, v. 0.9.1