/*
  Copyright (C) 2000-2008

  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: xmds_vector_element.cc 1885 2008-03-18 15:24:56Z paultcochrane $
*/

/*! @file xmds_vector_element.cc
  @brief Vector element parsing classes and methods

  More detailed explanation...
*/

#include <xmds_common.h>
#include <xmds_vector_element.h>
#include <xmds_simulation.h>

// **************************************************************************
// **************************************************************************
//                              xmdsVectorElement public
// **************************************************************************
// **************************************************************************

extern bool debugFlag;

long nxmdsVectorElements=0;   //!< The number of xmds vector elements

enum {
  ASCII = 1,
  BINARY = 2,
  XSIL = 3
};

// **************************************************************************
xmdsVectorElement::xmdsVectorElement(
                                     const xmdsSimulation *const yourSimulation,
                                     const bool& yourVerboseMode,
                                     const xmdsField *const yourField) :
  xmdsVector(yourField),
  xmdsElement(yourSimulation, yourVerboseMode) {
  if (debugFlag) {
    nxmdsVectorElements++;
    printf("xmdsVectorElement::xmdsVectorElement\n");
    printf("nxmdsVectorElements=%li\n", nxmdsVectorElements);
  }
}

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

// **************************************************************************
void xmdsVectorElement::processElement(
                                       const Element *const yourElement) {
  if (debugFlag) {
    printf("xmdsVectorElement::processElement\n");
  }
  
  list<XMLString> anXMLStringList;
  list<long> aLongList;
  XMLString anXMLString;
  
  if (verbose()) {
    printf("Processing vector element ...\n");
  }
  
  // ************************************
  // find name
  
  getAssignmentStrings(yourElement, "name", 1, 1, anXMLStringList);
  
  const xmdsVector* possibleTwin;
  
  if (field()->getVector(*anXMLStringList.begin(), possibleTwin)) {
    sprintf(errorMessage(), "a vector of name '%s' already exists in this field", anXMLStringList.begin()->c_str());
    throw xmdsException(yourElement, errorMessage());
  }
  
  setName(*anXMLStringList.begin());
  if (verbose()) {
    printf("vector name is '%s'\n", name()->c_str());
  }
  
  // ************************************
  // find type
  
  getAssignmentStrings(yourElement, "type", 0, 1, anXMLStringList);
  
  if (anXMLStringList.size()==1) {
    
    if (*anXMLStringList.begin()=="complex") {
      setVectorType(COMPLEX);
      if (verbose()) {
        printf("vector type is 'complex'\n");
      }
    }
    else if (*anXMLStringList.begin()=="double") {
      setVectorType(DOUBLE);
      if (verbose()) {
        printf("vector type is 'double'\n");
      }
    }
    else {
    }
  }
  else {
    setVectorType(COMPLEX);
    printf("vector type defaulting to 'complex'\n");
  }
  
  // ************************************
  // find components
  
  getAssignmentVectorStrings(yourElement, "components", 1, 0, anXMLStringList, aLongList);
  
  if (anXMLStringList.size()==0) {
    throw xmdsException(yourElement, "No vector components specified!");
  }
  
  setComponents(anXMLStringList);
  setLengths(aLongList);
  
  if (verbose()) {
    for (unsigned long i=0; i<nComponents(); i++) {
      printf("adding vector component '%s'\n", componentName(i)->c_str());
    }
  }
  
  // ************************************
  // find 'filename'
  
  getAssignmentStrings(yourElement, "filename", NOT_REQD, 1, anXMLStringList);
  
  if (anXMLStringList.size()==1) {
    // we have a filename assignent
    
    getAttributeStrings(yourElement, "filename", "format", NOT_REQD, anXMLString);
    myInputFileFormat = ASCII;  // the default value
    
    setInitialSpace(0);
    
    if (anXMLString != EMPTY_STRING) {
      // a format attribute has been set
      if (anXMLString == "binary") {
        myInputFileFormat = BINARY;
      }
      else if (anXMLString == "bin") {
        myInputFileFormat = BINARY;
      }
      else if (anXMLString == "ascii") {
        myInputFileFormat = ASCII;
      }
      else if (anXMLString == "text") {
        myInputFileFormat = ASCII;
      }
      else if (anXMLString == "txt") {
        myInputFileFormat = ASCII;
      }
      else if (anXMLString == "xsil") {
        myInputFileFormat = XSIL;
        
        // find the moment group that we are initialising from
        XMLString momentGroupName;
        getAttributeStrings(yourElement, "filename", "moment_group", NOT_REQD, momentGroupName);
        if (momentGroupName == EMPTY_STRING)
          myInitialisationMomentGroupName = "NULL";
        else
          myInitialisationMomentGroupName = string("moment_group_") + momentGroupName.c_str();
        
        
        // find the geometry matching mode
        XMLString geometryMatchMode;
        getAttributeStrings(yourElement, "filename", "geometry_matching_mode", NOT_REQD, geometryMatchMode);
        if (geometryMatchMode != EMPTY_STRING) {
          if (geometryMatchMode == "strict")
            isGeometryMatchingModeStrict = true;
          else if (geometryMatchMode == "loose")
            isGeometryMatchingModeStrict = false;
          else
            throw xmdsException(yourElement, "The geometry matching mode must be either 'strict' or 'loose'");
        }
        else
          isGeometryMatchingModeStrict = true;
        
        if (field()->geometry()->nDims()>0) {
          
          // ************************************
          // find space
          
          list<bool> aSpaceList;
          
          getAssignmentBools(yourElement, "fourier_space", 0, field()->geometry()->nDims(), aSpaceList);
          
          if (aSpaceList.size() == 0) {
            printf("Initialisation space for vector '%s' defaulting to x-space.\n", name()->c_str());
          }
          else {
            list<bool>::const_iterator pBool = aSpaceList.begin();
            for (unsigned long i=0; i<field()->geometry()->nDims(); i++) {
              if (verbose()) {
                if (*pBool) {
                  printf("initialisation will be performed with dimension #%li in fourier space\n", i+1);
                }
                else {
                  printf("initialisation will be performed with dimension #%li in normal space\n", i+1);
                }
              }
              pBool++;
            }
            setInitialSpace(spaceList2ULong(aSpaceList));
          }
        }
        // ************************************
        // find vectors
        
        getAssignmentStrings(yourElement, "vectors", 0, 0, myVectorNamesList);
        
        field()->processVectors(myVectorNamesList, initialSpace());
        
        // ************************************
        // find code
        
        myCode=*yourElement->textContent(0);
        
      }
      else {
        throw xmdsException(yourElement,
                  "Warning: Unknown file format attribute value in 'filename' tag\nI expected either 'binary', 'ascii' or 'xsil'.");
      }
    }
    
    myFileName=*anXMLStringList.begin();
    if (verbose()) {
      printf("vector initialisation from file '%s'\n", myFileName.c_str());
    }
    
  }
  else {
    // initialisation from code
    if (verbose()) {
      printf("initialisation from user code\n");
    }
    
    setInitialSpace(0);
    
    if (field()->geometry()->nDims()>0) {
      
      // ************************************
      // find space
      
      list<bool> aSpaceList;
      
      getAssignmentBools(yourElement, "fourier_space", 0, field()->geometry()->nDims(), aSpaceList);
      
      if (aSpaceList.size() == 0) {
        printf("Initialisation space for vector '%s' defaulting to x-space.\n", name()->c_str());
      }
      else {
        list<bool>::const_iterator pBool = aSpaceList.begin();
        for (unsigned long i=0; i<field()->geometry()->nDims(); i++) {
          if (verbose()) {
            if (*pBool) {
              printf("initialisation will be performed with dimension #%li in fourier space\n", i+1);
            }
            else {
              printf("initialisation will be performed with dimension #%li in normal space\n", i+1);
            }
          }
          pBool++;
        }
        setInitialSpace(spaceList2ULong(aSpaceList));
      }
    }
    
    // ************************************
    // find vectors
    
    getAssignmentStrings(yourElement, "vectors", 0, 0, myVectorNamesList);
    
    field()->processVectors(myVectorNamesList, initialSpace());
    
    // ************************************
    // find code
    
    myCode=*yourElement->textContent(0);
    
    // check for non-white space charaters in code:
    
    if (myCode.isAllWhiteSpace()) {
      throw xmdsException(yourElement, "No initialisation code defined!");
    }
    
    if (verbose()) {
      printf("initialisation code loaded\n");
    }
  }
}

// **************************************************************************
// **************************************************************************
//                              xmdsVectorElement private
// **************************************************************************
// **************************************************************************

// **************************************************************************
void xmdsVectorElement::writeInitialiseRoutine(FILE *const outfile) const
{
  if (debugFlag) {
    printf("xmdsVectorElement::writeInitialise\n");
  }

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

  fprintf(outfile,
          "// *************************\n"
          "void _%s_%s_initialise() {\n\n",
          fieldName,
          vectorName);

  if (simulation()->parameters()->useOpenMP) {
    // This may seem a little bizarre, but this is important for supercomputers that use a NUMA
    // (non-uniform memory addressing) architecture (e.g. ANU's ac supercomputer)
    // Memory is assigned to the process that first accesses it, so if we initialise the whole field to zero
    // in a parallel fashion, then memory should be assigned to the thread that will use it in future. If we just
    // did this with a single thread, then all the memory would be 'local' on the first thread, and any time any other
    // thread tries to access this memory, it will need to talk to the node on which this memory is stored.

    fprintf(outfile, "#ifdef _OPENMP\n" // this code is unnecessary if OpenMP isn't actually available
            "{\n");  // ensure that loop variables are all local

    // Begin loop over the grid (for this vector)
    list<XMLString> vectorNameList;
    vectorNameList.push_back(*name());
    field()->openLoops(outfile, initialSpace(), vectorNameList, PARALLELISE_LOOP);
    for (long unsigned int i=0; i<nComponents(); i++) {
      fprintf(outfile, "\t %s = 0.0;\n", componentName(i)->c_str());
    }
    field()->closeLoops(outfile, initialSpace(), vectorNameList);
    fprintf(outfile, "}\n"
            "#endif // _OPENMP\n\n");
  }

  if (myFileName.length()>0) {
    // initialisation from file

    if (myInputFileFormat == ASCII) {

      fprintf(outfile,
              "FILE* infile = fopen(\"%s\", \"r\");\n"
              "\n"
              "if (infile==0) {\n"
              "\t printf(\"Error opening file '%s' for initialisation of vector '%s'\\n\");\n"
              "\t exit(255);\n"
              "\t }\n"
              "\n"
              "unsigned long _i0=0;\n"
              "\n",
              myFileName.c_str(),
              myFileName.c_str(),
              vectorName);
      if (simulation()->parameters()->usempi&!simulation()->parameters()->stochastic) {
        //forwarding to "ranks' point in file. Slow but safe
        fprintf(outfile, "double _garbage=0.0;\n");
        fprintf(outfile, "for (unsigned long _i0=0; _i0<(_%s_size/_%s_lattice0)*_%s_%s_ncomponents*local_x_start; _i0++){\n", fieldName, fieldName, fieldName, vectorName);

        if (vectorType()==COMPLEX) {
          fprintf(outfile, "\t if (fscanf(infile, \"%%lf %%lf\", &_garbage, &_garbage) != 2) \n");
        }
        else if (vectorType()==DOUBLE) {
          fprintf(outfile, "\t if (fscanf(infile, \"%%lf\", &_garbage) != 1) \n");
        }
        fprintf(outfile,
                "\t\t printf(\"Rank %%i Error forwarding '%s' from file '%s': either bad float format or insufficient data\\n\", rank);\n", vectorName,
                myFileName.c_str());

        fprintf(outfile, "}\n");
        fprintf(outfile, "\n");
        fprintf(outfile, "while(_i0<total_local_size*_%s_%s_ncomponents) {\n",
                fieldName,
                fieldName);
      }
      else{
        fprintf(outfile, "while(_i0<_%s_size*_%s_%s_ncomponents) {\n",
                fieldName,
                fieldName,
                vectorName);
      }

      if (vectorType()==COMPLEX) {
        fprintf(outfile, "\t if (fscanf(infile, \"%%lf %%lf\", &_%s_%s[_i0].re, &_%s_%s[_i0].im) != 2) {\n",
                fieldName, vectorName, fieldName, vectorName);
      }
      else if (vectorType()==DOUBLE) {
        fprintf(outfile, "\t if (fscanf(infile, \"%%lf\", &_%s_%s[_i0]) != 1) {\n", fieldName, vectorName);
      }
      fprintf(outfile,
              "\t\t printf(\"Error loading '%s' from file '%s': either bad float format or insufficient data\\n\");\n", vectorName,
              myFileName.c_str());

      if (simulation()->parameters()->usempi&!simulation()->parameters()->stochastic)
        fprintf(outfile, "\t\t _i0=total_local_size*_%s_%s_ncomponents;\n",
                fieldName,
                vectorName);
      else
        fprintf(outfile, "\t\t _i0=_%s_size*_%s_%s_ncomponents;\n",
                fieldName,
                fieldName,
                vectorName);
      /*  This next line kills the program at this point.
          We may or may not want to do this, although I just had a user
          request that the program not continue past this error. (Wasted time)
          Reverting is as simple as deleting the next line.  (JJH)  */
      fprintf(outfile, "\t\t exit(255);\n"
              "\t }\n"
              "\t _i0++;\n"
              "}\n"
              "\n"
              "fclose(infile);\n");

    }
    else if (myInputFileFormat == BINARY) {
      fprintf(outfile,
              "FILE* infile = fopen(\"%s\", \"rb\");\n"
              "\n"
              "if (infile==0) {\n"
              "\t printf(\"Error opening file '%s' for initialisation of vector '%s'\\n\");\n"
              "\t return;\n"
              "\t }\n"
              "\n",
              myFileName.c_str(),
              myFileName.c_str(),
              vectorName);

      // Begin loop over the grid (for this vector)
      list<XMLString> vectorNameList;
      vectorNameList.push_back(*name());
      // This loop isn't parallelisable due to the use of fseek/fread
      field()->openLoops(outfile, initialSpace(), vectorNameList, DO_NOT_PARALLELISE_LOOP);
      // Index pointer into the input array
      fprintf(outfile, "unsigned long _inputfield_index_pointer=0;\n");

      for (long unsigned int i=0; i<field()->geometry()->nDims(); i++) {
        string dimName;
        if (!((initialSpace()>>i)&1)) {
          // x space
          dimName = field()->geometry()->dimension(i)->name.c_str();
          fprintf(outfile, "long _inputfield_i%1$li = lround((%2$s - _%3$s_xmin%1$li)/_%3$s_dx%1$li);\n\n",
                  i, dimName.c_str(), fieldName);
        }
        else {
          // k space
          dimName = string("k") + field()->geometry()->dimension(i)->name.c_str();
          fprintf(outfile, "long _inputfield_i%1$li = lround(%2$s/_%3$s_dk%1$li) + (_%3$s_lattice%1$li/2);\n"
                  , i, dimName.c_str(), fieldName);
        }

        // Calculate the input field index pointer
        fprintf(outfile, "_inputfield_index_pointer += _inputfield_i%li", i);
        for (long unsigned int j=i+1; j<field()->geometry()->nDims(); j++) {
          fprintf(outfile, "*_%s_lattice%li", fieldName, j);
        }
        fprintf(outfile, ";\n\n");

      }


      fprintf(outfile, "fseek(infile, _inputfield_index_pointer*sizeof(%s)*_%s_%s_ncomponents, SEEK_SET);\n", vectorType()==COMPLEX ? "complex" : "double", fieldName, vectorName);
      fprintf(outfile, "fread(_%1$s_%2$s + _%1$s_%2$s_index_pointer, sizeof(%3$s), _%1$s_%2$s_ncomponents, infile);\n", fieldName, vectorName, vectorType()==COMPLEX ? "complex" : "double");

      field()->closeLoops(outfile, initialSpace(), vectorNameList);
      /*
        if (vectorType()==COMPLEX) {
        //  We are cutting the temporary variable due to stack overflow.  Why was this done this way originally?
        // fprintf(outfile, "\t complex _tempInput[_%s_size*_%s_%s_ncomponents];\n",
        //   fieldName, fieldName, vectorName);
        // fprintf(outfile, "\t fread(&_tempInput, sizeof(complex), _%s_size*_%s_%s_ncomponents, infile);\n",
        //   fieldName, fieldName, vectorName);
        fprintf(outfile, "\t fread(_%s_%s, sizeof(complex), _%s_size*_%s_%s_ncomponents, infile);\n", fieldName, vectorName,
        fieldName, fieldName, vectorName);
        // fprintf(outfile, "\t _%s_%s = _tempInput;\n", fieldName, vectorName);
        }
        else if (vectorType()==DOUBLE) {
        // fprintf(outfile, "\t double _tempInput[_%s_size*_%s_%s_ncomponents];\n",
        //   fieldName, fieldName, vectorName);
        // fprintf(outfile, "\t fread(&_tempInput, sizeof(double), _%s_size*_%s_%s_ncomponents, infile);\n",
        //   fieldName, fieldName, vectorName);
        fprintf(outfile, "\t fread(_%s_%s, sizeof(double), _%s_size*_%s_%s_ncomponents, infile);\n", fieldName, vectorName,
        fieldName, fieldName, vectorName);
        // fprintf(outfile, "\t _%s_%s = _tempInput;\n", fieldName, vectorName);
        }*/
      fprintf(outfile,
              "\n"
              "fclose(infile);\n");
    }
    else if (myInputFileFormat == XSIL) {
      // Create noises if needed
      if (simulation()->parameters()->stochastic) {
        fprintf(outfile, "const double _var = 1.0");
        for (unsigned long 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 = new double[_n_noises];\n");
        if (simulation()->parameters()->errorCheck) {
          fprintf(outfile, "double *_noises2 = new double[_n_noises];\n");
        }
        fprintf(outfile, "\n");
      }
      // Write variables for this routine
      fprintf(outfile, "\t char **dimNames  = new char* [_%1$s_ndims];\n"
              "\t double *dimDelta = new double[_%1$s_ndims];\n" // use the first argument. Saves writing it out half a dozen times
              "\t double *dimMin   = new double[_%1$s_ndims];\n"
              "\t unsigned long *dimLattice = new unsigned long[_%1$s_ndims];\n\n",
              fieldName);
      for (long unsigned int i=0; i<field()->geometry()->nDims(); i++) {
        fprintf(outfile, "\t dimNames[%1$li]   = \"%2$s%3$s\";\n"
                "\t dimDelta[%1$li]   = _%4$s_d%5$s%1$li;\n",
                i,
                (initialSpace()>>i)&1 ? "k" : "",
                field()->geometry()->dimension(i)->name.c_str(),
                fieldName,
                (initialSpace()>>i)&1 ? "k" : "x");
        if ((initialSpace()>>i)&1)
          fprintf(outfile, "\t dimMin[%1$li]     = (-(_%2$s_lattice%1$li)/2)*_%2$s_dk%1$li;\n",
                  i, fieldName);
        else
          fprintf(outfile, "\t dimMin[%1$li]     = _%2$s_xmin%1$li;\n",
                  i, fieldName);
        fprintf(outfile, "\t dimLattice[%1$li] = _main_lattice%1$li;\n\n", i);
      }

      int vectorComponents = nComponents();

      fprintf(outfile, "\t char **componentNames = new char*[%1$i+1];\n"
              "\t int *componentFieldIndices = new int[%1$i];\n\n",
              vectorType()==DOUBLE ? vectorComponents : 2*vectorComponents);

      int index=0;
      for (int naturalIndex=0; naturalIndex<vectorComponents; naturalIndex++) {
        fprintf(outfile, "\t componentNames[%i] = \"%sR\";\n",
                index, componentName(naturalIndex)->c_str());
        index++;

        if (vectorType()==COMPLEX) {
          fprintf(outfile, "\t componentNames[%i] = \"%sI\";\n", index, componentName(naturalIndex)->c_str());
          index++;
        }
        fprintf(outfile, "\n");
      }
      fprintf(outfile, "\t componentNames[%i] = NULL;\n\n", index);

      // This are the output variables for initialiseFieldFromXSILFile
      fprintf(outfile, "\t char *binaryDataFilename;\n"
              "\t int unsignedLongSize;\n"
              "\t bool dataEncodingIsNative;\n"
              "\t bool isPrecisionDouble;\n"
              "\t unsigned long nDataComponents;\n"
              "\t unsigned long *inputLattice;\n"
              "\t int *componentInputIndices;\n\n"
              "\t if (!initialiseFieldFromXSILFile(\"%s\", \"%s\", _%s_ndims, dimNames, componentNames, \n" //file, moment group and field
              "// output variables\n"
              "\t\t\t &binaryDataFilename, &unsignedLongSize, &dataEncodingIsNative, &isPrecisionDouble, &nDataComponents, &inputLattice, &componentInputIndices)) {\n",
              myFileName.c_str(), myInitialisationMomentGroupName.c_str(), fieldName);
      if (simulation()->parameters()->usempi && !simulation()->parameters()->stochastic) {
        fprintf(outfile, "\t\t printf(\"Rank [%%i]: Unable to load data from XSIL file.\\nExiting.\\n\", rank);\n"
                "\t\t MPI_Abort(MPI_COMM_WORLD, 1);\n"
                "\t\t exit(1);\n");
      }
      else {
        fprintf(outfile, "\t\t printf(\"Unable to load data from XSIL file.\\nExiting.\\n\");\n"
                "\t\t exit(1);\n");
      }
      fprintf(outfile, "\t }\n\n");

      // Open the binary data file
      fprintf(outfile, "\t FILE *inputFile = fopen(binaryDataFilename, \"rb\");\n"
              "\t if (inputFile == NULL) {\n"
              "\t\t printf(\"Unable to open binary input file %%s\\n\", binaryDataFilename);\n"
              "\t\t exit(1);\n"
              "\t }\n"
              "// I would like to point out that the binary data described by the XSIL file does *NOT* follow the file format described in the spec.\n\n"
              "typedef union {\n"
              "\t float *floatPtr;\n"
              "\t double *doublePtr;\n"
              "} GenericRealNumberPtr;\n\n"
              "typedef union {\n"
              "\t unsigned long ulong;\n"
              "\t uint32_t uint32;\n"
              "\t uint64_t uint64;\n"
              "} GenericUnsignedLong;\n\n"
              "\t GenericUnsignedLong size;\n"
              "\t GenericRealNumberPtr *inputData = new GenericRealNumberPtr[_%s_ndims];\n"
              "\t for (int __i=0; __i<_%s_ndims; __i++) {\n" // loop over the dimension elements in the binary file
              "\t\t uint64_t sizeValue;\n"
              "\t\t switch(unsignedLongSize) {\n", fieldName, fieldName); // deal with the different possible sizes for the unsigned long size
      const char *string1[] = {"0", "4", "8"};
      const char *string2[] = {"ulong", "uint32", "uint64"};
      const char *string3[] = {"unsigned long", "uint32_t", "uint64_t"};
      const char *string4[] = {"XMDSSwapInt32", "XMDSSwapInt32", "XMDSSwapInt64"};
      for (int i=0; i<3; i++) {
        fprintf(outfile, "\t\t\t case %s:\n"
                "\t\t\t\t fread(&size.%s, sizeof(%s), 1, inputFile);\n"
                "\t\t\t\t if (dataEncodingIsNative)\n"
                "\t\t\t\t\t sizeValue = size.%s;\n"
                "\t\t\t\t else\n"
                "\t\t\t\t\t sizeValue = %s(size.%s);\n"
                "\t\t\t\t unsignedLongSize = sizeof(%s);\n"
                "\t\t\t\t break;\n\n",
                string1[i], string2[i], string3[i], string2[i], string4[i], string2[i], string3[i]);
      }
      fprintf(outfile, "\t\t }\n");         // Read in the data (in either double or single precision)
      fprintf(outfile, "\t\t if (isPrecisionDouble) {\n"
              "\t\t\t inputData[__i].doublePtr = new double[sizeValue];\n"
              "\t\t\t fread(inputData[__i].doublePtr, sizeof(double), sizeValue, inputFile);\n"
              "\t\t\t if (!dataEncodingIsNative) {\n"
              "\t\t\t\t for (int j=0; j<sizeValue; j++)\n"
              "\t\t\t\t\t inputData[__i].doublePtr[j] = XMDSSwapDouble(inputData[__i].doublePtr[j]);\n"
              "\t\t\t }\n"
              "\t\t }\n"
              "\t\t else {\n"
              "\t\t\t inputData[__i].floatPtr = new float[sizeValue];\n"
              "\t\t\t fread(inputData[__i].floatPtr, sizeof(float), sizeValue, inputFile);\n"
              "\t\t\t if (!dataEncodingIsNative) {\n"
              "\t\t\t\t for (int j=0; j<sizeValue; j++)\n"
              "\t\t\t\t\t inputData[__i].floatPtr[j] = XMDSSwapFloat(inputData[__i].floatPtr[j]);\n"
              "\t\t\t }\n"
              "\t\t }\n"
              //               "\t\t printf(\"We loaded dimension %%i with %%lu points\\n\", i, sizeValue);\n"
              "\t }\n"
              "// Assert that the deltas are the same up to 1%%\n"    // One restriction on the XSIL input is that the grid spacing must be the same
              "\t for (int __i=0; __i<_%s_ndims; __i++) {\n" // field name
              "\t\t double difference;\n"
              "\t\t if (isPrecisionDouble)\n"
              "\t\t\t difference = inputData[__i].doublePtr[1] - inputData[__i].doublePtr[0];\n"
              "\t\t else\n"
              "\t\t\t difference = inputData[__i].floatPtr[1] - inputData[__i].floatPtr[0];\n"
              "\t\t if (fabs(dimDelta[__i] - difference) > 0.01*dimDelta[__i]) {\n"
              "\t\t\t printf(\"The step size in the '%%s' dimension of the input data and the simulation grid do not match\\n\", dimNames[__i]);\n"
              "\t\t\t printf(\"The step size in the '%%s' dimension was %%e, while the input data had a step size of %%e\\n\", dimNames[__i], dimDelta[__i], difference);\n"
              "\t\t\t exit(1);\n"
              "\t\t }\n"
              "\t\t else {\n"
              //               "\t\t\t printf(\"Step requirement succeeded. Field delta: %%e, input delta: %%e\\n\", dimDelta[__i], difference);\n"
              "\t\t }\n"
              "\t }\n", fieldName);
      // In strict mode we require that the grid be exactly the same
      if (isGeometryMatchingModeStrict) {
        fprintf(outfile, "// STRICT MODE: assert that the start point is the same to within 10%% of one delta, and that the number of lattice points is the same\n"
                "\t for (int __i=0; __i<_%s_ndims; __i++) {\n" // field name
                "\t\t double start;\n"
                "\t\t if (isPrecisionDouble)\n"
                "\t\t\t start = inputData[__i].doublePtr[0];\n"
                "\t\t else\n"
                "\t\t\t start = inputData[__i].floatPtr[0];\n\n"
                "\t\t if (fabs(dimMin[__i] - start) > dimDelta[__i]/10) {\n"
                "\t\t\t printf(\"Geometry matching mode is strict, so the starting coordinate of each dimension in the field must be the same as in the input grid\\n\");\n"
                "\t\t\t printf(\"The problem is with dimension '%%s'\\n\", dimNames[__i]);\n"
                "\t\t\t exit(1);\n"
                "\t\t }\n"
                "\t\t if (dimLattice[__i] != inputLattice[__i]) {\n"
                "\t\t\t printf(\"Geometry matching mode is strict, so the number of lattice points in each dimension of the field must be the same as in the input grid\\n\");\n"
                "\t\t\t printf(\"The problem is with dimension '%%s'\\n\", dimNames[__i]);\n"
                "\t\t\t exit(1);\n"
                "\t\t }\n"
                "\t }\n\n", fieldName);
      }

      fprintf(outfile, "long binaryFileBaseOffset = ftell(inputFile);\n");
      fprintf(outfile, "long realNumberSize = isPrecisionDouble ? sizeof(double) : sizeof(float);\n");
      fprintf(outfile, "off_t vectorFieldSize = unsignedLongSize + realNumberSize");
      for (long unsigned int i=0; i<field()->geometry()->nDims(); i++) {
        fprintf(outfile, " * inputLattice[%li]", i);
      }
      fprintf(outfile, ";\n");

      // Transform vectors to appropriate space
      field()->vectors2space(outfile, initialSpace(), myVectorNamesList, "");

      // I'm honestly not sure why this code is here
      list<XMLString> vectorNameList = myVectorNamesList;
      vectorNameList.push_back(*name());

      // Create minimum variables for x-space dimensions
      for (long unsigned int i=0; i<field()->geometry()->nDims(); i++) {
        if ((initialSpace()>>i)&1)
          continue;

        fprintf(outfile, "double _inputfield_min%1$li;\n"
                "if (isPrecisionDouble)\n"
                "\t _inputfield_min%1$li = inputData[%1$li].doublePtr[0];\n"
                "else\n"
                "\t _inputfield_min%1$li = inputData[%1$li].floatPtr[0];\n\n",
                i);
      }

      fprintf(outfile, "{ // Put the loop in a block\n");
      field()->openLoops(outfile, initialSpace(), vectorNameList, PARALLELISE_LOOP);
      if (simulation()->parameters()->stochastic) {
        if (simulation()->parameters()->errorCheck) {
          if (simulation()->parameters()->noiseKind == "poissonian") {
            fprintf(outfile, "_make_noises(_gen1, _var/2, 0, _noises, _n_noises);\n");
            fprintf(outfile, "_make_noises(_gen2, _var/2, 0, _noises2, _n_noises);\n");
          }
          else {
            fprintf(outfile, "_make_noises(_gen1, _var/2, _noises, _n_noises);\n");
            fprintf(outfile, "_make_noises(_gen2, _var/2, _noises2, _n_noises);\n");
          }
          fprintf(outfile, "for (unsigned long _s0=0; _s0<_n_noises; _s0++)\n");
          fprintf(outfile, " _noises[_s0] += _noises2[_s0];\n");
          fprintf(outfile, "\n");
        }
        else {
          if (simulation()->parameters()->noiseKind == "poissonian") {
            fprintf(outfile, "_make_noises(_gen, _var, 0, _noises, _n_noises);\n");
          }
          else {
            fprintf(outfile, "_make_noises(_gen, _var, _noises, _n_noises);\n");
          }
          fprintf(outfile, "\n");
        }
      }

      // Set the input field index pointer from the grid coordinates

      // Initialise the components to zero
      fprintf(outfile, "for (unsigned int _component=0; _component<_%1$s_%2$s_ncomponents; _component++)\n"
              "\t _%1$s_%2$s[_%1$s_%2$s_index_pointer+_component] = 0.0;\n\n", fieldName, vectorName);
      // Include any initialisation code, if it exists
      if (!myCode.isAllWhiteSpace()) {
        fprintf(outfile, "// ********** Code from vector element***********\n");
        fprintf(outfile, "%s\n", myCode.c_str());
        fprintf(outfile, "// **********************************************\n");
        fprintf(outfile, "\n");
      }
      // Of course, initialisation code is only used if the input grid doesn't specify a value there.
      // This is so, if for some bizarre reason, a vector needs to be set to '1' (or something else) outside the input grid...
      // If you want to do something where you override the input data, then that can be done easily enough with a filter element
      // at the start of the simulation.

      field()->closeLoops(outfile, initialSpace(), vectorNameList);
      fprintf(outfile, "} // End block for loop\n");


      // Index pointer into the input array (new and old)
      fprintf(outfile, "off_t _inputfield_index_pointer, _inputfield_old_index_pointer;\n");

      fprintf(outfile, "for (unsigned int _component=0; _component<%i; _component++) {\n"
              "\tif (componentInputIndices[_component] == -1)\n"
              "\t\tcontinue;\n\n"
              "_inputfield_index_pointer = -42; // Just so that we always seek the first time\n",
              vectorType() == DOUBLE ? vectorComponents : 2*vectorComponents);

      field()->openLoops(outfile, initialSpace(), vectorNameList, DO_NOT_PARALLELISE_LOOP);

      // Save the old index pointer, and clear the new one
      fprintf(outfile, "_inputfield_old_index_pointer = _inputfield_index_pointer;\n"
              "_inputfield_index_pointer = 0;\n\n");


      for (long unsigned int i=0; i<field()->geometry()->nDims(); i++)
        fprintf(outfile, "long _inputfield_i%li;\n", i);

      for (long unsigned int i=0; i<field()->geometry()->nDims(); i++) {
        string dimName;
        if (!((initialSpace()>>i)&1)) {
          // x space
          dimName = field()->geometry()->dimension(i)->name.c_str();
          fprintf(outfile, "_inputfield_i%1$li = lround((%2$s - _inputfield_min%1$li)/dimDelta[%1$li]);\n\n",
                  i, dimName.c_str());
        }
        else {
          // k space
          dimName = string("k") + field()->geometry()->dimension(i)->name.c_str();
          fprintf(outfile, "_inputfield_i%1$li = lround(%2$s/dimDelta[%1$li]) + (inputLattice[%1$li]/2);\n",
                  i, dimName.c_str());
        }

        // check to see if the input grid has this coordinate value
        fprintf(outfile, "if (_inputfield_i%1$li < 0 || _inputfield_i%1$li >= inputLattice[%1$li])\n"
                "// The input grid doesn't cover this point\n"
                "\t goto ENDLOOP;\n\n", i);

        // check that the input coordinate matches up with the field coordinate
        fprintf(outfile, "{\n\tdouble _inputfield_%1$s;\n"
                "\tif (isPrecisionDouble)\n"
                "\t\t _inputfield_%1$s = inputData[%2$li].doublePtr[_inputfield_i%2$li];\n"
                "\telse\n"
                "\t\t _inputfield_%1$s = inputData[%2$li].floatPtr[_inputfield_i%2$li];\n\n"
                "\tif (fabs(_inputfield_%1$s - %1$s) > dimDelta[%2$li]/10) {\n"
                "\t// This only happens if your input field grid does not exactly match up with the simulation grid\n"
                "\t\t printf(\"The input field coordinate in the '%1$s' dimension does not match up with the field coordinate\\n\");\n"
                "\t\t printf(\"_i%2$li: %%li, %1$s: %%e, _inputfield_i%2$li: %%li, _inputfield_%1$s: %%e, d%1$s: %%e, diff/Delta: %%e\\n\", \n"
                "\t\t\t\t\t\t  _i%2$li, %1$s, _inputfield_i%2$li, _inputfield_%1$s, dimDelta[%2$li], fabs(_inputfield_%1$s - %1$s)/dimDelta[%2$li]);\n"
                "\t\t exit(1);\n"
                "\t}\n}\n\n",
                dimName.c_str(), i);

        // Calculate the input field index pointer
        fprintf(outfile, "_inputfield_index_pointer += _inputfield_i%li", i);
        for (long unsigned int j=i+1; j<field()->geometry()->nDims(); j++) {
          fprintf(outfile, "*inputLattice[%li]", j);
        }
        fprintf(outfile, ";\n\n");

      }

      fprintf(outfile, "double value;\n");
      fprintf(outfile, "if (_inputfield_index_pointer != _inputfield_old_index_pointer + 1) \n"
              "\t fseeko(inputFile, binaryFileBaseOffset + (componentInputIndices[_component]-_%s_ndims)*vectorFieldSize + \n"
              "\t\t unsignedLongSize + _inputfield_index_pointer*realNumberSize, SEEK_SET);\n", fieldName);
      fprintf(outfile, "if (isPrecisionDouble) {\n"
              "\t fread(&value, sizeof(double), 1, inputFile);\n"
              "\t if (!dataEncodingIsNative)\n"
              "\t\t value = XMDSSwapDouble(value);\n"
              " }\n"
              "\t else {\n"
              "\t float temp;\n"
              "\t fread(&temp, sizeof(float), 1, inputFile);\n"
              "\t if (!dataEncodingIsNative)\n"
              "\t\t temp = XMDSSwapFloat(temp);\n"
              "\t value = (double)temp;\n"
              " }\n");
      if (vectorType() == DOUBLE)
        fprintf(outfile, "_%1$s_%2$s[_%1$s_%2$s_index_pointer+_component] = value;\n", fieldName, vectorName);
      else // COMPLEX
        fprintf(outfile, "if (_component & 1) \n"
                "\t _%1$s_%2$s[_%1$s_%2$s_index_pointer+_component/2].im = value;\n"
                "else \n"
                "\t _%1$s_%2$s[_%1$s_%2$s_index_pointer+_component/2].re = value;\n\n", fieldName, vectorName);

      fprintf(outfile, "ENDLOOP:\n\n");
      field()->closeLoops(outfile, initialSpace(), vectorNameList);
      fprintf(outfile, "} // end loop over components \n");
      fprintf(outfile, "\t fclose(inputFile);\n"
              "\t delete [] dimNames;\n"
              "\t delete [] componentNames;\n"
              "\t free(binaryDataFilename);\n"
              "\t delete [] inputLattice;\n"
              "\t delete [] componentInputIndices;\n"
              "\t for (int i =0; i<_%s_ndims; i++) {\n"
              "\t\t if (isPrecisionDouble)\n"
              "\t\t\t delete [] inputData[i].doublePtr;\n"
              "\t\t else\n"
              "\t\t\t delete [] inputData[i].floatPtr;\n"
              "\t }\n"
              "\t delete [] inputData;\n", fieldName);
      // Delete noises (if needed)
      if (simulation()->parameters()->stochastic) {
        fprintf(outfile, "    delete[] _noises;\n");

        if (simulation()->parameters()->errorCheck) {
          fprintf(outfile, "    delete[] _noises2;\n");
        }
      }
    }
    else {
      // something must have really screwed up to get here...
      throw xmdsException("For some reason the input file is neither ascii, binary or xsil...\nHow did we get here??");
    }
  }
  else {
    // initialisation from code
    if (simulation()->parameters()->stochastic) {
      fprintf(outfile, "const double _var = 1.0");
      for (unsigned long 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 = new double[_n_noises];\n");
      if (simulation()->parameters()->errorCheck) {
        fprintf(outfile, "double *_noises2 = new double[_n_noises];\n");
      }
      fprintf(outfile, "\n");
    }

    field()->vectors2space(outfile, initialSpace(), myVectorNamesList, "");

    list<XMLString> vectorNameList = myVectorNamesList;
    vectorNameList.push_back(*name());

    // The code could be doing anything, let's be safe
    field()->openLoops(outfile, initialSpace(), vectorNameList, DO_NOT_PARALLELISE_LOOP);

    char indent[64];
    for (unsigned long i=0; i<nDims; i++) {
      indent[i]=0x09;
    }
    indent[nDims]=0;

    if (simulation()->parameters()->stochastic) {
      if (simulation()->parameters()->errorCheck) {
        if (simulation()->parameters()->noiseKind == "poissonian") {
          fprintf(outfile, "%s_make_noises(_gen1, _var/2, 0, _noises, _n_noises);\n", indent);
          fprintf(outfile, "%s_make_noises(_gen2, _var/2, 0, _noises2, _n_noises);\n", indent);
        }
        else {
          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 {
        if (simulation()->parameters()->noiseKind == "poissonian") {
          fprintf(outfile, "%s_make_noises(_gen, _var, 0, _noises, _n_noises);\n", indent);
        }
        else {
          fprintf(outfile, "%s_make_noises(_gen, _var, _noises, _n_noises);\n", indent);
        }
        fprintf(outfile, "\n");
      }
    }

    fprintf(outfile, "// ********** Code from vector element***********\n");
    fprintf(outfile, "%s\n", myCode.c_str());
    fprintf(outfile, "// **********************************************\n");
    fprintf(outfile, "\n");

    field()->closeLoops(outfile, initialSpace(), vectorNameList);

    if (simulation()->parameters()->stochastic) {
      fprintf(outfile, "    delete[] _noises;\n");

      if (simulation()->parameters()->errorCheck) {
        fprintf(outfile, "    delete[] _noises2;\n");
      }
    }

  }

  if (needsFFTWRoutines()) {
    fprintf(outfile, "\n");
    fprintf(outfile, "_%s_%s_space=%li;\n", field()->name()->c_str(), name()->c_str(), initialSpace());
  }

  fprintf(outfile, "}\n");
  fprintf(outfile, "\n");
}

/*
 * Local variables:
 * c-indentation-style: bsd
 * c-basic-offset: 2
 * indent-tabs-mode: nil
 * End:
 *
 * vim: tabstop=2 expandtab shiftwidth=2:
 */