ValueLayerND.h

/*
 *   Repast for High Performance Computing (Repast HPC)
 *
 *   Copyright (c) 2010 Argonne National Laboratory
 *   All rights reserved.
 *
 *   Redistribution and use in source and binary forms, with
 *   or without modification, are permitted provided that the following
 *   conditions are met:
 *
 *       Redistributions of source code must retain the above copyright notice,
 *       this list of conditions and the following disclaimer.
 *
 *       Redistributions in binary form must reproduce the above copyright notice,
 *       this list of conditions and the following disclaimer in the documentation
 *       and/or other materials provided with the distribution.
 *
 *       Neither the name of the Argonne National Laboratory nor the names of its
 *     contributors may be used to endorse or promote products derived from
 *     this software without specific prior written permission.
 *
 *   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 *   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 *   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
 *   PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE TRUSTEES OR
 *   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 *   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 *   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 *   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 *   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 *   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
 *   EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 *
 *  ValueLayerND.h
 *
 *  Created on: July 18, 2016
 *      Author: jtm
 */
 
#ifndef VALUELAYERND_H_
#define VALUELAYERND_H_
 
#include <fstream>
 
#include "mpi.h"
 
#include "Point.h"
#include "GridDimensions.h"
#include "RepastProcess.h"
 
 
using namespace std;
 
namespace repast {
 
struct RankDatum{
  int            rank;
  MPI_Datatype   datatype;
  int            sendPtrOffset;
  int            receivePtrOffset;
  int            sendDir;  // Integer representing the direction a send will be sent, in N-space
  int            recvDir;  // Integer representing the direction a receive will have been sent, in N-space
};
 
template<typename T>
class DimensionDatum{
public:
 
  DimensionDatum(int indx, GridDimensions globalBoundaries, GridDimensions localBoundaries, int buffer, bool isPeriodic);
 
  virtual ~DimensionDatum(){}
 
  int  globalCoordinateMin, globalCoordinateMax;            // Coordinates of the global simulation boundaries
  int  localBoundariesMin, localBoundariesMax;              // Coordinates of the local process boundaries in this layer
  int  simplifiedBoundariesMin, simplifiedBoundariesMax;    // Coordinates of the 'simplified' coordinates (assuming local boundary origins and NO wrapping)
  int  leftBufferSize, rightBufferSize;                     // Buffer size
  int  matchingCoordinateMin, matchingCoordinateMax;        // Region within simplified boundaries within which coordinates match Global Simulation Coordinate system
  bool periodic;                                            // True if the global simulation space is periodic on this dimension
  bool atLeftBound, atRightBound;                           // True if the local boundaries abut the global boundaries (whether periodic or not)
  bool spaceContinuesLeft, spaceContinuesRight;             // False if the local boundaries abut the (non-periodic) global boundaries
  int  globalWidth;                                         // Global width of the simulation boundaries in simulation units
  int  localWidth;                                          // Width of the local boundaries
  int  width;                                               // Total width (units = sizeof T) on this dimension (local extents + buffer sizes if not against global nonperiodic bounds)
  int  widthInBytes;
 
  int  getSendReceiveSize(int relativeLocation);
 
  int getTransformedCoord(int originalCoord);
 
  int getIndexedCoord(int originalCoord, bool isSimplified = false);
 
  bool isInLocalBounds(int originalCoord);
 
 
  void report(int dimensionNumber){
    std::cout << repast::RepastProcess::instance()->rank() << " " << dimensionNumber << " " <<
                 "global (" << globalCoordinateMin     << ", " << globalCoordinateMax     << ") " <<
                 "local  (" << localBoundariesMin      << ", " << localBoundariesMax      << ") " <<
                 "simple (" << simplifiedBoundariesMin << ", " << simplifiedBoundariesMax << ") " <<
                 "match  (" << matchingCoordinateMin   << ", " << matchingCoordinateMax   << ") " <<
                 "globalWidth = " << globalWidth << " localWidth = " << localWidth << " width = " << width << " bytes = " << widthInBytes << std::endl;
  }
};
 
template<typename T>
DimensionDatum<T>::DimensionDatum(int indx, GridDimensions globalBoundaries, GridDimensions localBoundaries, int buffer, bool isPeriodic):
    leftBufferSize(buffer), rightBufferSize(buffer), periodic(isPeriodic){
  globalCoordinateMin = globalBoundaries.origin(indx);
  globalCoordinateMax = globalBoundaries.origin(indx) + globalBoundaries.extents(indx);
  localBoundariesMin  = localBoundaries.origin(indx);
  localBoundariesMax  = localBoundaries.origin(indx) + localBoundaries.extents(indx);
 
  atLeftBound  = localBoundariesMin == globalCoordinateMin;
  atRightBound = localBoundariesMax == globalCoordinateMax;
 
  spaceContinuesLeft  = !atLeftBound  || periodic;
  spaceContinuesRight = !atRightBound || periodic;
 
  simplifiedBoundariesMin  = localBoundariesMin - leftBufferSize;
  simplifiedBoundariesMax  = localBoundariesMax + rightBufferSize;
 
  matchingCoordinateMin    = localBoundariesMin;
  if(spaceContinuesLeft  && !atLeftBound ) matchingCoordinateMin -= leftBufferSize;
 
  matchingCoordinateMax    = localBoundariesMax;
  if(spaceContinuesRight && !atRightBound) matchingCoordinateMax += rightBufferSize;
 
  globalWidth = globalCoordinateMax - globalCoordinateMin;
  localWidth = localBoundariesMax - localBoundariesMin;
  width = leftBufferSize + localWidth + rightBufferSize;
  widthInBytes = width * (sizeof(T));
}
 
 
template<typename T>
int DimensionDatum<T>::getSendReceiveSize(int relativeLocation){
  switch(relativeLocation){
    case -1:  return leftBufferSize;
    case  1:  return rightBufferSize;
    case  0:
    default:
      return localWidth;
  }
}
 
template<typename T>
int DimensionDatum<T>::getTransformedCoord(int originalCoord){
  if(originalCoord < matchingCoordinateMin){        // Assume (!) original is on right (!) side of periodic boundary, starting at some value
    return matchingCoordinateMax + (originalCoord - globalCoordinateMin);
  }
  else if(originalCoord > matchingCoordinateMax){
    return matchingCoordinateMin - (globalCoordinateMax - originalCoord);
  }
  else return originalCoord; // Within matching boundaries; no need to transform
 
}
 
template<typename T>
int DimensionDatum<T>::getIndexedCoord(int originalCoord, bool isSimplified){
  return (isSimplified ? originalCoord : getTransformedCoord(originalCoord)) - simplifiedBoundariesMin;
}
 
template<typename T>
bool DimensionDatum<T>::isInLocalBounds(int originalCoord){
  return originalCoord >= localBoundariesMin && originalCoord < localBoundariesMax;
}
 
/*******************************************************************/
 
template<typename T>
class AbstractValueLayerND{
 
private:
  bool dummy; // Used for cases when error flag is not requested
 
protected:
  CartesianTopology*         cartTopology;
  GridDimensions             localBoundaries;
  int                        length;                 // Total length of the entire array (one data space)
 
  int                        numDims;                // Number of dimensions
  bool                       globalSpaceIsPeriodic;  // True if the global space is periodic
 
  vector<int>                places;                 // Multipliers to calculate index, for each dimension
  vector<int>                strides;                // Sizes of each dimensions, in bytes
  vector<DimensionDatum<T> > dimensionData;          // List of data for each dimension
  RankDatum*                 neighborData;           // List of data for each adjacent rank
  int                        neighborCount;          // Count of adjacent ranks
  MPI_Request*               requests;               // Pointer to MPI requests (for wait operations)
 
  int                        instanceID;             // Unique ID for managing MPI requests without mix-ups
  int                        syncCount;
 
  AbstractValueLayerND(vector<int> processesPerDim, GridDimensions globalBoundaries, int bufferSize, bool periodic);
  virtual ~AbstractValueLayerND();
 
 
public:
 
  static int instanceCount;
 
  virtual bool isInLocalBounds(vector<int> coords);
 
  /*
   * Returns true only if the coordinates given are within the local boundaries
   *
   * @param coords Coordinates to be tested
   * @return true if the coordinates are within the local boundaries
   */
  virtual bool isInLocalBounds(Point<int> location);
 
  const GridDimensions& getLocalBoundaries(){
    return localBoundaries;
  }
 
protected:
  // Methods implemented in this class but visible only to child classes:
 
  vector<int> getIndexes(vector<int> location, bool isSimplified = false);
 
  int getIndex(vector<int> location, bool isSimplified = false);
 
 
  int getIndex(Point<int> location);
 
 
  // Virtual methods (implemented by child classes
 
  virtual void initialize(T initialValue, bool fillBufferZone = false, bool fillLocal = true) = 0;
 
  virtual void initialize(T initialLocalValue, T initialBufferZoneValue) = 0;
 
  virtual T addValueAt(T val, Point<int> location, bool& errFlag) = 0;
 
  virtual T addValueAt(T val, vector<int> location, bool& errFlag) = 0;
 
  virtual T setValueAt(T val, Point<int> location, bool& errFlag) = 0;
 
  virtual T setValueAt(T val, vector<int> location, bool& errFlag) = 0;
 
  virtual T getValueAt(Point<int> location, bool& errFlag) = 0;
 
  virtual T getValueAt(vector<int> location, bool& errFlag) = 0;
 
  virtual void synchronize() = 0;
 
private:
 
  void getMPIDataType(RelativeLocation relLoc, MPI_Datatype &datatype);
 
  void getMPIDataType(int radius, MPI_Datatype &datatype);
 
  void getMPIDataType(vector<int> sideLengths, MPI_Datatype &datatype, int dimensionIndex);
 
  MPI_Datatype getRawMPIDataType();
 
 
  int getSendPointerOffset(RelativeLocation relLoc);
 
  int getReceivePointerOffset(RelativeLocation relLoc);
 
 
};
 
 
 
template<typename T>
int AbstractValueLayerND<T>::instanceCount = 0;
 
template<typename T>
AbstractValueLayerND<T>::AbstractValueLayerND(vector<int> processesPerDim, GridDimensions globalBoundaries,int bufferSize, bool periodic): globalSpaceIsPeriodic(periodic), syncCount(0){
  instanceID = AbstractValueLayerND<T>::instanceCount;
  AbstractValueLayerND<T>::instanceCount++;
  cartTopology = RepastProcess::instance()->getCartesianTopology(processesPerDim, periodic);
  // Calculate the size to be used for the buffers
  numDims = processesPerDim.size();
 
  int rank = RepastProcess::instance()->rank();
  localBoundaries = cartTopology->getDimensions(rank, globalBoundaries);
 
  // First create the basic coordinate data per dimension
  length = 1;
  int val = 1;
  for(int i = 0; i < numDims; i++){
    DimensionDatum<T> datum(i, globalBoundaries, localBoundaries, bufferSize, periodic);
    length *= datum.width;
    dimensionData.push_back(datum);
    places.push_back(val);
    strides.push_back(val * sizeof(T));
    val *= dimensionData[i].width;
  }
 
  // Now create the rank-based data per neighbor
  RelativeLocation relLoc(numDims);
  RelativeLocation relLocTrimmed = cartTopology->trim(rank, relLoc); // Initialized to minima
 
  vector<int> myCoordinates;
  cartTopology->getCoordinates(rank, myCoordinates);
 
  neighborData = new RankDatum[relLoc.getMaxIndex()];
  neighborCount = 0;
  do{
    if(relLoc.validNonCenter()){ // Skip 0,0,0,0,0
      RankDatum* datum;
      datum = &neighborData[neighborCount];
      // Collect the information about this rank here
      getMPIDataType(relLoc, datum->datatype);
      datum->sendPtrOffset    = getSendPointerOffset(relLoc);
      datum->receivePtrOffset = getReceivePointerOffset(relLoc);
      vector<int> current = relLoc.getCurrentValue();
      datum->rank = cartTopology->getRank(myCoordinates, current);
      datum->sendDir = RelativeLocation::getDirectionIndex(current);
      datum->recvDir = RelativeLocation::getReverseDirectionIndex(current);
 
      neighborCount++;
    }
  }while(relLoc.increment());
 
  // Create arrays for MPI requests and results (statuses)
  requests = new MPI_Request[neighborCount * 2];
}
 
template<typename T>
AbstractValueLayerND<T>::~AbstractValueLayerND(){
  delete[] neighborData; // Should Free MPI Datatypes first...
  delete[] requests;
}
 
template<typename T>
bool AbstractValueLayerND<T>::isInLocalBounds(vector<int> coords){
  for(int i = 0; i < numDims; i++){
    DimensionDatum<T>* datum = &dimensionData[i];
    if(!datum->isInLocalBounds(coords[i])) return false;
  }
  return true;
}
 
template<typename T>
bool AbstractValueLayerND<T>::isInLocalBounds(Point<int> location){
  return isInLocalBounds(location.coords());
}
 
template<typename T>
vector<int> AbstractValueLayerND<T>::getIndexes(vector<int> location, bool isSimplified){
  vector<int> ret;
  ret.assign(numDims, 0); // Make the right amount of space
  for(int i = 0; i < numDims; i++) ret[i] = dimensionData[i].getIndexedCoord(location[i], isSimplified);
  return ret;
}
 
template<typename T>
int AbstractValueLayerND<T>::getIndex(vector<int> location, bool isSimplified){
  vector<int> indexed = getIndexes(location, isSimplified);
  int val = 0;
  for(int i = numDims - 1; i >= 0; i--) val += indexed[i] * places[i];
  if(val < 0 || val > length) val = -1;
  return val;
}
 
template<typename T>
int AbstractValueLayerND<T>::getIndex(Point<int> location){
  return getIndex(location.coords());
}
 
 
template<typename T>
void AbstractValueLayerND<T>::getMPIDataType(RelativeLocation relLoc, MPI_Datatype &datatype){
  vector<int> sideLengths;
  for(int i = 0; i < numDims; i++) sideLengths.push_back(dimensionData[i].getSendReceiveSize(relLoc[i]));
  getMPIDataType(sideLengths, datatype, numDims - 1);
}
 
template<typename T>
void AbstractValueLayerND<T>::getMPIDataType(int radius, MPI_Datatype &datatype){
  vector<int> sideLengths;
  sideLengths.assign(numDims, 2 * radius + 1);
  getMPIDataType(sideLengths, datatype, numDims - 1);
}
 
template<typename T>
void AbstractValueLayerND<T>::getMPIDataType(vector<int> sideLengths, MPI_Datatype &datatype, int dimensionIndex){
  if(dimensionIndex == 0){
    MPI_Type_contiguous(sideLengths[dimensionIndex], getRawMPIDataType(), &datatype);
  }
  else{
    MPI_Datatype innerType;
    getMPIDataType(sideLengths, innerType, dimensionIndex - 1);
    MPI_Type_create_hvector(sideLengths[dimensionIndex], // Count
                     1,                                                                        // BlockLength: just one of the inner data type
                     strides[dimensionIndex],                                                  // Stride, in bytes
                     innerType,                                                                // Inner Datatype
                     &datatype);
  }
  // Commit?
  MPI_Type_commit(&datatype);
}
 
 
template<typename T>
int AbstractValueLayerND<T>::getSendPointerOffset(RelativeLocation relLoc){
  int ret = 0;
  for(int i = 0; i < numDims; i++){
    DimensionDatum<T>* datum = &dimensionData[i];
    ret += (relLoc[i] <= 0 ? datum->leftBufferSize : datum->width - (2 * datum->rightBufferSize)) * places[i];
  }
  return ret;
}
 
template<typename T>
int AbstractValueLayerND<T>::getReceivePointerOffset(RelativeLocation relLoc){
  int ret = 0;
  for(int i = 0; i < numDims; i++){
    DimensionDatum<T>* datum = &dimensionData[i];
    ret += (relLoc[i] < 0 ? 0 : (relLoc[i] == 0 ? datum->leftBufferSize : datum->width - datum->rightBufferSize)) * places[i];
  }
  return ret;
}
 
 
template<typename T>
class ValueLayerND: public AbstractValueLayerND<T>{
 
private:
  T* dataSpace;              // Pointer to the data space
 
public:
 
  ValueLayerND(vector<int> processesPerDim, GridDimensions globalBoundaries, int bufferSize,
      bool periodic, T initialValue = 0, T initialBufferZoneValue = 0);
  virtual ~ValueLayerND();
 
  virtual void initialize(T initialValue, bool fillBufferZone = false, bool fillLocal = true);
 
  virtual void initialize(T initialLocalValue, T initialBufferZoneValue);
 
  virtual T addValueAt(T val, Point<int> location, bool& errFlag);
 
  virtual T addValueAt(T val, vector<int> location, bool& errFlag);
 
  virtual T setValueAt(T val, Point<int> location, bool& errFlag);
 
  virtual T setValueAt(T val, vector<int> location, bool& errFlag);
 
  virtual T getValueAt(Point<int> location, bool& errFlag);
 
  virtual T getValueAt(vector<int> location, bool& errFlag);
 
  virtual void synchronize();
 
  void write(string fileLocation, string filetag, bool writeSharedBoundaryAreas = false);
 
 
private:
 
 
  void fillDimension(T localValue, T bufferZoneValue, bool doBufferZone, bool doLocal, T* dataSpacePointer, int dimIndex);
 
  /*
   * Writes one dimension's information to the specified csv file.
   *
   * @param outfile output file
   * @param dataSpacePointer pointer to the data space to be written
   * @param currentPosition position currently being written (for recursive calls)
   * @param dimIndex dimension currently being written (for recursive calls)
   * @param writeSharedBoundaryAreas if true, write the areas that are non-local to this process
   */
  void writeDimension(std::ofstream& outfile, T* dataSpacePointer, int* currentPosition, int dimIndex, bool writeSharedBoundaryAreas = false);
 
};
 
 
 
 
 
template<typename T>
class ValueLayerNDSU: public AbstractValueLayerND<T>{
 
protected:
 
  T*                dataSpace1;             // Permanent pointer to bank 1 of the data space
  T*                dataSpace2;             // Permanent pointer to bank 2 of the data space
  T*                currentDataSpace;       // Temporary pointer to the active data space
  T*                otherDataSpace;         // Temporary pointer to the inactive data space
 
public:
 
  ValueLayerNDSU(vector<int> processesPerDim, GridDimensions globalBoundaries, int bufferSize, bool periodic, T initialValue = 0, T initialBufferZoneValue = 0);
  virtual ~ValueLayerNDSU();
 
  virtual void initialize(T initialValue, bool fillBufferZone = false, bool fillLocal = true);
 
  virtual void initialize(T initialLocalValue, T initialBufferZoneValue);
 
  virtual T addValueAt(T val, Point<int> location, bool& errFlag);
 
  virtual T addValueAt(T val, vector<int> location, bool& errFlag);
 
  virtual T setValueAt(T val, Point<int> location, bool& errFlag);
 
  virtual T setValueAt(T val, vector<int> location, bool& errFlag);
 
  virtual T getValueAt(Point<int> location, bool& errFlag);
 
  virtual T getValueAt(vector<int> location, bool& errFlag);
 
  virtual void synchronize();
 
  virtual void write(string fileLocation, string filetag, bool writeSharedBoundaryAreas = false);
 
  void switchValueLayer();
 
  virtual T addSecondaryValueAt(T val, Point<int> location, bool& errFlag);
 
  virtual T addSecondaryValueAt(T val, vector<int> location, bool& errFlag);
 
  virtual T setSecondaryValueAt(T val, Point<int> location, bool& errFlag);
 
  virtual T setSecondaryValueAt(T val, vector<int> location, bool& errFlag);
 
  virtual T getSecondaryValueAt(Point<int> location, bool& errFlag);
 
  virtual T getSecondaryValueAt(vector<int> location, bool& errFlag);
 
  virtual void copyCurrentToSecondary();
 
  virtual void copySecondaryToCurrent();
 
  void flowback();
 
private:
 
  void fillDimension(T localValue, T bufferZoneValue, bool doBufferZone, bool doLocal, T* dataSpace1Pointer, T* dataSpace2Pointer, int dimIndex);
 
  /*
   * Writes one dimension's information to the specified csv file.
   *
   * @param outfile output file
   * @param dataSpacePointer pointer to the data space to be written
   * @param currentPosition position currently being written (for recursive calls)
   * @param dimIndex dimension currently being written (for recursive calls)
   * @param writeSharedBoundaryAreas if true, write the areas that are non-local to this process
   */
  void writeDimension(std::ofstream& outfile, T* dataSpacePointer, int* currentPosition, int dimIndex, bool writeSharedBoundaryAreas = false);
 
  void sumInto(T* dataSpace1Pointer, T* dataSpace2Pointer, int dimIndex);
 
};
 
 
 
template<typename T>
ValueLayerND<T>::ValueLayerND(vector<int> processesPerDim, GridDimensions globalBoundaries, int bufferSize, bool periodic,
    T initialValue, T initialBufferZoneValue): AbstractValueLayerND<T>(processesPerDim, globalBoundaries, bufferSize, periodic){
 
  // Create the actual arrays for the data
  dataSpace = new T[AbstractValueLayerND<T>::length];
 
  // Finally, fill the data with the initial values
  initialize(initialValue, initialBufferZoneValue);
 
  // And synchronize
  synchronize();
 
}
 
template<typename T>
ValueLayerND<T>::~ValueLayerND(){
  delete[] dataSpace;
}
 
template<typename T>
void ValueLayerND<T>::initialize(T initialValue, bool fillBufferZone, bool fillLocal){
  fillDimension(initialValue, initialValue, fillBufferZone, fillLocal, dataSpace, AbstractValueLayerND<T>::numDims - 1);
}
 
template<typename T>
void ValueLayerND<T>::initialize(T initialLocalValue, T initialBufferZoneValue){
  fillDimension(initialLocalValue, initialBufferZoneValue, true, true, dataSpace, AbstractValueLayerND<T>::numDims - 1);
}
 
template<typename T>
T ValueLayerND<T>::addValueAt(T val, Point<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return val;
  }
  T* pt = &dataSpace[indx];
  return (*pt = *pt + val);
}
 
template<typename T>
T ValueLayerND<T>::addValueAt(T val, vector<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return val;
  }
 
  T* pt = &dataSpace[indx];
  return (*pt = *pt + val);
}
 
template<typename T>
T ValueLayerND<T>::setValueAt(T val, Point<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return val;
  }
  T* pt = &dataSpace[indx];
  return (*pt = val);
}
 
template<typename T>
T ValueLayerND<T>::setValueAt(T val, vector<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return val;
  }
  T* pt = &dataSpace[indx];
  return (*pt = val);
}
 
template<typename T>
T ValueLayerND<T>::getValueAt(vector<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return 0;
  }
  return dataSpace[indx];
}
 
template<typename T>
T ValueLayerND<T>::getValueAt(Point<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return 0;
  }
  return dataSpace[indx];
}
 
template<typename T>
void ValueLayerND<T>::synchronize(){
  AbstractValueLayerND<T>::syncCount++;
  if(AbstractValueLayerND<T>::syncCount > 9) AbstractValueLayerND<T>::syncCount = 0;
  int mpiTag = AbstractValueLayerND<T>::instanceID * 10 + AbstractValueLayerND<T>::syncCount;
  // Note: the syncCount and send/recv directions are used to create a unique tag value for the
  // mpi sends and receives. The tag value must be unique in two ways: first, successive calls to this
  // function must be different enough that they can't be confused. The 'syncCount' value is used to
  // achieve this, and it will loop from 0-9 and then repeat. The second, the tag must sometimes
  // differentiate between sends and receives that are going to the same rank. If a dimension
  // has only 2 processes but wrap-around borders, then one process may be sending to the other
  // process twice (once left and once right). The 'sendDir' and 'recvDir' values trap this
 
  // For each entry in neighbors:
  MPI_Status statuses[AbstractValueLayerND<T>::neighborCount * 2];
  for(int i = 0; i < AbstractValueLayerND<T>::neighborCount; i++){
    MPI_Isend(&dataSpace[AbstractValueLayerND<T>::neighborData[i].sendPtrOffset], 1, AbstractValueLayerND<T>::neighborData[i].datatype,
        AbstractValueLayerND<T>::neighborData[i].rank, 10 * (AbstractValueLayerND<T>::neighborData[i].sendDir + 1) + mpiTag, AbstractValueLayerND<T>::cartTopology->topologyComm, &AbstractValueLayerND<T>::requests[i]);
    MPI_Irecv(&dataSpace[AbstractValueLayerND<T>::neighborData[i].receivePtrOffset], 1, AbstractValueLayerND<T>::neighborData[i].datatype,
        AbstractValueLayerND<T>::neighborData[i].rank, 10 * (AbstractValueLayerND<T>::neighborData[i].recvDir + 1) + mpiTag, AbstractValueLayerND<T>::cartTopology->topologyComm, &AbstractValueLayerND<T>::requests[AbstractValueLayerND<T>::neighborCount + i]);
  }
  MPI_Waitall(AbstractValueLayerND<T>::neighborCount * 2, AbstractValueLayerND<T>::requests, statuses);
}
 
 
template<typename T>
void ValueLayerND<T>::write(string fileLocation, string fileTag, bool writeSharedBoundaryAreas){
  std::ofstream outfile;
  std::ostringstream stream;
  int rank = repast::RepastProcess::instance()->rank();
  stream << fileLocation << "ValueLayer_" << fileTag << "_" << rank << ".csv";
  std::string filename = stream.str();
 
  const char * c = filename.c_str();
  outfile.open(c, std::ios_base::trunc | std::ios_base::out); // it will not delete the content of file, will add a new line
 
  // Write headers
  for(int i = 0; i < AbstractValueLayerND<T>::numDims; i++) outfile << "DIM_" << i << ",";
  outfile << "VALUE" << endl;
 
  int* positions = new int[AbstractValueLayerND<T>::numDims];
  for(int i = 0; i < AbstractValueLayerND<T>::numDims; i++) positions[i] = 0;
 
  writeDimension(outfile, dataSpace, positions, AbstractValueLayerND<T>::numDims - 1, writeSharedBoundaryAreas);
 
  outfile.close();
}
 
 
template<typename T>
void ValueLayerND<T>::fillDimension(T localValue, T bufferValue, bool doBufferZone, bool doLocal, T* dataSpacePointer, int dimIndex){
  if(!doBufferZone && !doLocal) return;
  int bufferEdge = AbstractValueLayerND<T>::dimensionData[dimIndex].leftBufferSize;
  int localEdge  = bufferEdge + AbstractValueLayerND<T>::dimensionData[dimIndex].localWidth;
  int upperBound = localEdge + AbstractValueLayerND<T>::dimensionData[dimIndex].rightBufferSize;
 
  int pointerIncrement = AbstractValueLayerND<T>::places[dimIndex];
 
 
  int i = 0;
  for(; i < bufferEdge; i++){
    if(doBufferZone){
      if(dimIndex == 0){
        *dataSpacePointer = bufferValue;
      }
      else{
        fillDimension(bufferValue, bufferValue, doBufferZone, doLocal, dataSpacePointer, dimIndex - 1);
      }
    }
    // Increment the pointers
    dataSpacePointer += pointerIncrement;
  }
  for(; i < localEdge; i++){
    if(doLocal){
      if(dimIndex == 0){
        *dataSpacePointer = localValue;
      }
      else{
        fillDimension(localValue, bufferValue, doBufferZone, doLocal, dataSpacePointer, dimIndex - 1);
      }
    }
    // Increment the pointers
    dataSpacePointer += pointerIncrement;
  }
  if(doBufferZone){ // Note: we don't need to finish this at all if not doing buffer zone
    for(; i < upperBound; i++){
      if(dimIndex == 0){
        *dataSpacePointer = bufferValue;
      }
      else{
        fillDimension(bufferValue, bufferValue, doBufferZone, doLocal, dataSpacePointer, dimIndex - 1);
      }
    }
    dataSpacePointer += pointerIncrement;
  }
 
}
 
template<typename T>
void ValueLayerND<T>::writeDimension(std::ofstream& outfile, T* dataSpacePointer, int* currentPosition, int dimIndex, bool writeSharedBoundaryAreas){
  int bufferEdge = AbstractValueLayerND<T>::dimensionData[dimIndex].leftBufferSize;
  int localEdge  = bufferEdge + AbstractValueLayerND<T>::dimensionData[dimIndex].localWidth;
  int upperBound = localEdge + AbstractValueLayerND<T>::dimensionData[dimIndex].rightBufferSize;
 
  int pointerIncrement = AbstractValueLayerND<T>::places[dimIndex];
  int i = 0;
  for(; i < bufferEdge; i++){
    currentPosition[dimIndex] = i;
    if(writeSharedBoundaryAreas){
      if(dimIndex == 0){
        T val = *dataSpacePointer;
        if(val != 0){
          for(int j = 0; j < AbstractValueLayerND<T>::numDims; j++) outfile << (currentPosition[j] - AbstractValueLayerND<T>::dimensionData[j].leftBufferSize + AbstractValueLayerND<T>::dimensionData[j].localBoundariesMin) << ",";
          outfile << val << endl;
        }
      }
      else{
        writeDimension(outfile, dataSpacePointer, currentPosition, dimIndex - 1, writeSharedBoundaryAreas);
      }
    }
    // Increment the pointers
    dataSpacePointer += pointerIncrement;
  }
  for(; i < localEdge; i++){
    currentPosition[dimIndex] = i;
    if(dimIndex == 0){
        T val = *dataSpacePointer;
        if(val != 0){
          for(int j = 0; j < AbstractValueLayerND<T>::numDims; j++) outfile << (currentPosition[j] - AbstractValueLayerND<T>::dimensionData[j].leftBufferSize + AbstractValueLayerND<T>::dimensionData[j].localBoundariesMin) << ",";
          outfile << val << endl;
        }
    }
    else{
      writeDimension(outfile, dataSpacePointer, currentPosition, dimIndex - 1, writeSharedBoundaryAreas);
    }
    // Increment the pointers
    dataSpacePointer += pointerIncrement;
  }
  if(writeSharedBoundaryAreas){ // Note: we don't need to finish this at all if not doing buffer zone
    for(; i < upperBound; i++){
      currentPosition[dimIndex] = i;
      if(dimIndex == 0){
        T val = *dataSpacePointer;
        if(val != 0){
          for(int j = 0; j < AbstractValueLayerND<T>::numDims; j++) outfile << (currentPosition[j] - AbstractValueLayerND<T>::dimensionData[j].leftBufferSize + AbstractValueLayerND<T>::dimensionData[j].localBoundariesMin) << ",";
          outfile << *dataSpacePointer << endl;
        }
      }
      else{
        writeDimension(outfile, dataSpacePointer, currentPosition, dimIndex - 1, writeSharedBoundaryAreas);
      }
    }
    dataSpacePointer += pointerIncrement;
  }
 
}
 
 
 
 
template<typename T>
ValueLayerNDSU<T>::ValueLayerNDSU(vector<int> processesPerDim, GridDimensions globalBoundaries, int bufferSize, bool periodic,
    T initialValue, T initialBufferZoneValue): AbstractValueLayerND<T>(processesPerDim, globalBoundaries, bufferSize, periodic){
 
  // Create the actual arrays for the data
  dataSpace1 = new T[AbstractValueLayerND<T>::length];
  dataSpace2 = new T[AbstractValueLayerND<T>::length];
  currentDataSpace = dataSpace1;
  otherDataSpace   = dataSpace2;
 
  // Finally, fill the data with the initial values
  initialize(initialValue, initialBufferZoneValue);
 
  // And synchronize
  synchronize();
 
}
 
template<typename T>
ValueLayerNDSU<T>::~ValueLayerNDSU(){
  delete[] currentDataSpace;
  delete[] otherDataSpace;
}
 
template<typename T>
void ValueLayerNDSU<T>::initialize(T initialValue, bool fillBufferZone, bool fillLocal){
  fillDimension(initialValue, initialValue, fillBufferZone, fillLocal, dataSpace1, dataSpace2, AbstractValueLayerND<T>::numDims - 1);
}
 
template<typename T>
void ValueLayerNDSU<T>::initialize(T initialLocalValue, T initialBufferZoneValue){
  fillDimension(initialLocalValue, initialBufferZoneValue, true, true, dataSpace1, dataSpace2, AbstractValueLayerND<T>::numDims - 1);
}
 
template<typename T>
T ValueLayerNDSU<T>::addValueAt(T val, Point<int> location, bool& errFlag){
  int indx = this->getIndex(location);
  if(indx == -1) return nan("");
  T* pt = &currentDataSpace[indx];
  return (*pt = *pt + val);
}
 
template<typename T>
T ValueLayerNDSU<T>::addValueAt(T val, vector<int> location, bool& errFlag){
  int indx = this->getIndex(location);
  if(indx == -1) return nan("");
  T* pt = &currentDataSpace[indx];
  return (*pt = *pt + val);
}
 
template<typename T>
T ValueLayerNDSU<T>::setValueAt(T val, Point<int> location, bool& errFlag){
  int indx = this->getIndex(location);
  if(indx == -1) return nan("");
  T* pt = &currentDataSpace[indx];
  return (*pt = val);
}
 
template<typename T>
T ValueLayerNDSU<T>::setValueAt(T val, vector<int> location, bool& errFlag){
  int indx = this->getIndex(location);
  if(indx == -1) return nan("");
  T* pt = &currentDataSpace[indx];
  return (*pt = val);
}
 
template<typename T>
T ValueLayerNDSU<T>::getValueAt(Point<int> location, bool& errFlag){
  int indx = this->getIndex(location);
  if(indx == -1) return nan("");
  return currentDataSpace[indx];
}
 
template<typename T>
T ValueLayerNDSU<T>::getValueAt(vector<int> location, bool& errFlag){
  int indx = this->getIndex(location);
  if(indx == -1) return nan("");
  return currentDataSpace[indx];
}
 
 
template<typename T>
void ValueLayerNDSU<T>::synchronize(){
  AbstractValueLayerND<T>::syncCount++;
  if(AbstractValueLayerND<T>::syncCount > 9) AbstractValueLayerND<T>::syncCount = 0;
  int mpiTag = AbstractValueLayerND<T>::instanceID * 10 + AbstractValueLayerND<T>::syncCount;
  // Note: the syncCount and send/recv directions are used to create a unique tag value for the
  // mpi sends and receives. The tag value must be unique in two ways: first, successive calls to this
  // function must be different enough that they can't be confused. The 'syncCount' value is used to
  // achieve this, and it will loop from 0-9 and then repeat. The second, the tag must sometimes
  // differentiate between sends and receives that are going to the same rank. If a dimension
  // has only 2 processes but wrap-around borders, then one process may be sending to the other
  // process twice (once left and once right). The 'sendDir' and 'recvDir' values trap this
 
  // For each entry in neighbors:
  MPI_Status statuses[AbstractValueLayerND<T>::neighborCount * 2];
  for(int i = 0; i < AbstractValueLayerND<T>::neighborCount; i++){
    MPI_Isend(&currentDataSpace[AbstractValueLayerND<T>::neighborData[i].sendPtrOffset], 1, AbstractValueLayerND<T>::neighborData[i].datatype,
        AbstractValueLayerND<T>::neighborData[i].rank, 10 * (AbstractValueLayerND<T>::neighborData[i].sendDir + 1) + mpiTag, AbstractValueLayerND<T>::cartTopology->topologyComm, &AbstractValueLayerND<T>::requests[i]);
    MPI_Irecv(&currentDataSpace[AbstractValueLayerND<T>::neighborData[i].receivePtrOffset], 1, AbstractValueLayerND<T>::neighborData[i].datatype,
        AbstractValueLayerND<T>::neighborData[i].rank, 10 * (AbstractValueLayerND<T>::neighborData[i].recvDir + 1) + mpiTag, AbstractValueLayerND<T>::cartTopology->topologyComm, &AbstractValueLayerND<T>::requests[AbstractValueLayerND<T>::neighborCount + i]);
  }
  MPI_Waitall(AbstractValueLayerND<T>::neighborCount * 2, AbstractValueLayerND<T>::requests, statuses);
}
 
template<typename T>
void ValueLayerNDSU<T>::write(string fileLocation, string fileTag, bool writeSharedBoundaryAreas){
  std::ofstream outfile;
  std::ostringstream stream;
  int rank = repast::RepastProcess::instance()->rank();
  stream << fileLocation << "ValueLayer_" << fileTag << "_" << rank << ".csv";
  std::string filename = stream.str();
 
  const char * c = filename.c_str();
  outfile.open(c, std::ios_base::trunc | std::ios_base::out); // it will not delete the content of file, will add a new line
 
  // Write headers
  for(int i = 0; i < AbstractValueLayerND<T>::numDims; i++) outfile << "DIM_" << i << ",";
  outfile << "VALUE" << endl;
 
  int* positions = new int[AbstractValueLayerND<T>::numDims];
  for(int i = 0; i < AbstractValueLayerND<T>::numDims; i++) positions[i] = 0;
 
  writeDimension(outfile, currentDataSpace, positions, AbstractValueLayerND<T>::numDims - 1, writeSharedBoundaryAreas);
 
  outfile.close();
}
 
template<typename T>
void ValueLayerNDSU<T>::switchValueLayer(){
  // Switch the data banks
  T* tempDataSpace = currentDataSpace;
  currentDataSpace      = otherDataSpace;
  otherDataSpace        = tempDataSpace;
}
 
template<typename T>
T ValueLayerNDSU<T>::addSecondaryValueAt(T val, Point<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return val;
  }
  T* pt = &otherDataSpace[indx];
  return (*pt = *pt + val);
}
 
template<typename T>
T ValueLayerNDSU<T>::addSecondaryValueAt(T val, vector<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return val;
  }
  T* pt = &otherDataSpace[indx];
  return (*pt = *pt + val);
}
 
template<typename T>
T ValueLayerNDSU<T>::setSecondaryValueAt(T val, Point<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return val;
  }
  T* pt = &otherDataSpace[indx];
  return (*pt = val);
}
 
template<typename T>
T ValueLayerNDSU<T>::setSecondaryValueAt(T val, vector<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return val;
  }
  T* pt = &otherDataSpace[indx];
  return (*pt = val);
}
 
template<typename T>
T ValueLayerNDSU<T>::getSecondaryValueAt(Point<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return 0;
  }
  return otherDataSpace[indx];
}
 
template<typename T>
T ValueLayerNDSU<T>::getSecondaryValueAt(vector<int> location, bool& errFlag){
  errFlag = false;
  int indx = this->getIndex(location);
  if(indx == -1){
    errFlag = true;
    return 0;
  }
  return otherDataSpace[indx];
}
 
template<typename T>
void ValueLayerNDSU<T>::copyCurrentToSecondary(){
  T d = 0;
  memcpy(otherDataSpace, currentDataSpace, AbstractValueLayerND<T>::length * sizeof d);
}
 
template<typename T>
void ValueLayerNDSU<T>::copySecondaryToCurrent(){
  T d = 0;
  memcpy(currentDataSpace, otherDataSpace, AbstractValueLayerND<T>::length * sizeof d);
}
 
template<typename T>
void ValueLayerNDSU<T>::flowback(){
  std::cout << " valueLayer is executing flowback on rank " << repast::RepastProcess::instance()->rank() << std::endl;
  // Fill the other data space with zeros
  T d = 0;
  memset(otherDataSpace, 0, AbstractValueLayerND<T>::length * sizeof d);
 
  // Loop through the (3^N - 1)/2 directions
 
  // Note: TODO: Do all the sends as 'waitall', then process all the receives, with the final
  // sends waitall being the block before returning to simulation operation
 
  // For each direction:
  RelativeLocation relLoc(AbstractValueLayerND<T>::numDims);
  int listSize = AbstractValueLayerND<T>::neighborCount;
  std::cout << " valueLayer is executing flowback on rank " << repast::RepastProcess::instance()->rank() << " LIST SIZE = " << listSize << std::endl;
  MPI_Request* flowbackRequests = new MPI_Request[4]; // No more than four at a time (may be only 2)
  do{
    int LDir = relLoc.getIndex();
    int RDir = RelativeLocation::getReverseDirectionIndex(relLoc.getCurrentValue());
    // Need to create a send to the L and a send to the R, then matching receives
    // We create these in pairs because we can take a small advantage of non-blocking
    // sends but be assured that the L and R data are not going to interfere
    RankDatum* datum;
 
 
    int countOfRequests = 0;
    for(int i = 0; i < listSize; i++){ // To do: Shouldn't loop through this, but instead should index once
      datum = & AbstractValueLayerND<T>::neighborData[i];
      // Note: (!!!) Yes, we are using the 'send' pointer for the _receive_ and the 'receive' pointer for the send...
      if(datum->sendDir == LDir){
        MPI_Isend(&currentDataSpace[datum->receivePtrOffset], 1, datum->datatype,
                  datum->rank, 1001, AbstractValueLayerND<T>::cartTopology->topologyComm, &flowbackRequests[countOfRequests]);
        countOfRequests++;
      }
      else if(datum->sendDir == RDir){
        MPI_Isend(&currentDataSpace[datum->receivePtrOffset], 1, datum->datatype,
                  datum->rank, 2002, AbstractValueLayerND<T>::cartTopology->topologyComm, &flowbackRequests[countOfRequests]);
        countOfRequests++;
      }
      if(datum->recvDir == LDir){
        MPI_Irecv(&otherDataSpace[datum->sendPtrOffset], 1, datum->datatype,
                          datum->rank, 1001, AbstractValueLayerND<T>::cartTopology->topologyComm, &flowbackRequests[countOfRequests]);
        countOfRequests++;
      }
      else if(datum->recvDir == RDir){
        MPI_Irecv(&otherDataSpace[datum->sendPtrOffset], 1, datum->datatype,
                          datum->rank, 2002, AbstractValueLayerND<T>::cartTopology->topologyComm, &flowbackRequests[countOfRequests]);
        countOfRequests++;
      }
    }
    MPI_Status statuses[countOfRequests];
    // Perform the sends and receives
    std::cout << " valueLayer is executing flowback on rank " << repast::RepastProcess::instance()->rank() << " DOING WAITALL WITH LDIRECTION " << LDir << std::endl;
    MPI_Waitall(countOfRequests, flowbackRequests, statuses);
    std::cout << " valueLayer is executing flowback on rank " << repast::RepastProcess::instance()->rank() << " DONE WITH WAITALL... " << std::endl;
    // Copy the received data from the other data space into
    // the current data space, summing the values and clearing
    // the other data space to zeros
    sumInto(currentDataSpace, otherDataSpace, AbstractValueLayerND<T>::numDims - 1);
 
    relLoc.increment();
  }while(relLoc.validNonCenter()); // Note: STOP at 0,0,0,0...: because each location creates a send and a receive in both L and R, only need to do half
 
  delete[] flowbackRequests; // Cleanup
 
}
 
template<typename T>
void ValueLayerNDSU<T>::fillDimension(T localValue, T bufferValue, bool doBufferZone, bool doLocal, T* dataSpace1Pointer, T* dataSpace2Pointer, int dimIndex){
  if(!doBufferZone && !doLocal) return;
  int bufferEdge = AbstractValueLayerND<T>::dimensionData[dimIndex].leftBufferSize;
  int localEdge  = bufferEdge + AbstractValueLayerND<T>::dimensionData[dimIndex].localWidth;
  int upperBound = localEdge + AbstractValueLayerND<T>::dimensionData[dimIndex].rightBufferSize;
 
  int pointerIncrement = AbstractValueLayerND<T>::places[dimIndex];
 
 
  int i = 0;
  for(; i < bufferEdge; i++){
    if(doBufferZone){
      if(dimIndex == 0){
        *dataSpace1Pointer = bufferValue;
        *dataSpace2Pointer = bufferValue;
      }
      else{
        fillDimension(bufferValue, bufferValue, doBufferZone, doLocal, dataSpace1Pointer, dataSpace2Pointer, dimIndex - 1);
      }
    }
    // Increment the pointers
    dataSpace1Pointer += pointerIncrement;
    dataSpace2Pointer += pointerIncrement;
  }
  for(; i < localEdge; i++){
    if(doLocal){
      if(dimIndex == 0){
        *dataSpace1Pointer = localValue;
        *dataSpace2Pointer = localValue;
      }
      else{
        fillDimension(localValue, bufferValue, doBufferZone, doLocal, dataSpace1Pointer, dataSpace2Pointer, dimIndex - 1);
      }
    }
    // Increment the pointers
    dataSpace1Pointer += pointerIncrement;
    dataSpace2Pointer += pointerIncrement;
  }
  if(doBufferZone){ // Note: we don't need to finish this at all if not doing buffer zone
    for(; i < upperBound; i++){
      if(dimIndex == 0){
        *dataSpace1Pointer = bufferValue;
        *dataSpace2Pointer = bufferValue;
      }
      else{
        fillDimension(bufferValue, bufferValue, doBufferZone, doLocal, dataSpace1Pointer, dataSpace2Pointer, dimIndex - 1);
      }
    }
    dataSpace1Pointer += pointerIncrement;
    dataSpace2Pointer += pointerIncrement;
  }
 
}
 
template<typename T>
void ValueLayerNDSU<T>::writeDimension(std::ofstream& outfile, T* dataSpacePointer, int* currentPosition, int dimIndex, bool writeSharedBoundaryAreas){
  int bufferEdge = AbstractValueLayerND<T>::dimensionData[dimIndex].leftBufferSize;
  int localEdge  = bufferEdge + AbstractValueLayerND<T>::dimensionData[dimIndex].localWidth;
  int upperBound = localEdge + AbstractValueLayerND<T>::dimensionData[dimIndex].rightBufferSize;
 
  int pointerIncrement = AbstractValueLayerND<T>::places[dimIndex];
  int i = 0;
  for(; i < bufferEdge; i++){
    currentPosition[dimIndex] = i;
    if(writeSharedBoundaryAreas){
      if(dimIndex == 0){
        T val = *dataSpacePointer;
        if(val != 0){
          for(int j = 0; j < AbstractValueLayerND<T>::numDims; j++) outfile << (currentPosition[j] - AbstractValueLayerND<T>::dimensionData[j].leftBufferSize + AbstractValueLayerND<T>::dimensionData[j].localBoundariesMin) << ",";
          outfile << val << endl;
        }
      }
      else{
        writeDimension(outfile, dataSpacePointer, currentPosition, dimIndex - 1, writeSharedBoundaryAreas);
      }
    }
    // Increment the pointers
    dataSpacePointer += pointerIncrement;
  }
  for(; i < localEdge; i++){
    currentPosition[dimIndex] = i;
    if(dimIndex == 0){
        T val = *dataSpacePointer;
        if(val != 0){
          for(int j = 0; j < AbstractValueLayerND<T>::numDims; j++) outfile << (currentPosition[j] - AbstractValueLayerND<T>::dimensionData[j].leftBufferSize + AbstractValueLayerND<T>::dimensionData[j].localBoundariesMin) << ",";
          outfile << val << endl;
        }
    }
    else{
      writeDimension(outfile, dataSpacePointer, currentPosition, dimIndex - 1, writeSharedBoundaryAreas);
    }
    // Increment the pointers
    dataSpacePointer += pointerIncrement;
  }
  if(writeSharedBoundaryAreas){ // Note: we don't need to finish this at all if not doing buffer zone
    for(; i < upperBound; i++){
      currentPosition[dimIndex] = i;
      if(dimIndex == 0){
        T val = *dataSpacePointer;
        if(val != 0){
          for(int j = 0; j < AbstractValueLayerND<T>::numDims; j++) outfile << (currentPosition[j] - AbstractValueLayerND<T>::dimensionData[j].leftBufferSize + AbstractValueLayerND<T>::dimensionData[j].localBoundariesMin) << ",";
          outfile << *dataSpacePointer << endl;
        }
      }
      else{
        writeDimension(outfile, dataSpacePointer, currentPosition, dimIndex - 1, writeSharedBoundaryAreas);
      }
    }
    dataSpacePointer += pointerIncrement;
  }
 
}
 
 
template<typename T>
void ValueLayerNDSU<T>::sumInto(T* dataSpace1Pointer, T* dataSpace2Pointer, int dimIndex){
  //std::cout << " ON RANK " << repast::RepastProcess::instance()->rank() << " SUMINTO RUNNING WITH DIM INDEX OF " << dimIndex << std::endl;
  int bufferEdge = AbstractValueLayerND<T>::dimensionData[dimIndex].leftBufferSize;
  int localEdge  = bufferEdge + AbstractValueLayerND<T>::dimensionData[dimIndex].localWidth;
  //int upperBound = localEdge + AbstractValueLayerND<T>::dimensionData[dimIndex].rightBufferSize;
 
  int pointerIncrement = AbstractValueLayerND<T>::places[dimIndex];
  int leftPointerSkip  = pointerIncrement * AbstractValueLayerND<T>::dimensionData[dimIndex].leftBufferSize;
  //int rightPointerSkip = pointerIncrement * AbstractValueLayerND<T>::dimensionData[dimIndex].rightBufferSize;
 
  // Skip the left buffer zone
  dataSpace1Pointer += leftPointerSkip;
  dataSpace2Pointer += leftPointerSkip;
 
  // Loop local edge to local edge
  for(int i = bufferEdge; i < localEdge; i++){
    if(dimIndex == 0){
      //if(*dataSpace2Pointer > 0) std::cout << " ON RANK " << repast::RepastProcess::instance()->rank() << " SUMINTO FOUND VALUE OF " << *dataSpace2Pointer << " ADDING TO " << *dataSpace1Pointer << " at " << i << " dimIndex " << dimIndex << std::endl;
      *dataSpace1Pointer += *dataSpace2Pointer; // Add space 2's value into space 1
      *dataSpace2Pointer = 0;                   // Zero out space 2
    }
    else sumInto(dataSpace1Pointer, dataSpace2Pointer, dimIndex - 1); // Recursive call
    dataSpace1Pointer += pointerIncrement;
    dataSpace2Pointer += pointerIncrement;
  }
 
}
 
 
 
 
 
 
 
}
 
#endif /* VALUELAYERND_H_ */