<commit_before>#include "imageio.hpp"
#include "io.h"

#include <assert.h>
#include <string.h>

#ifdef PV_USE_GDAL
#  include <gdal_priv.h>
#  include <ogr_spatialref.h>
#else
#  define GDAL_CONFIG_ERR_STR "PetaVision must be compiled with GDAL to use this file type\n"
#endif // PV_USE_GDAL

#undef DEBUG_OUTPUT

static void copyToLocBuffer(int buf[], PVLayerLoc * loc)
{
   buf[0] = loc->nx;
   buf[1] = loc->ny;
   buf[2] = loc->nxGlobal;
   buf[3] = loc->nyGlobal;
   buf[4] = loc->kx0;
   buf[5] = loc->ky0;
   buf[6] = loc->nb;
   buf[7] = loc->nf;
   buf[8] = loc->halo.lt;
   buf[9] = loc->halo.rt;
   buf[10] = loc->halo.dn;
   buf[11] = loc->halo.up;
}

static void copyFromLocBuffer(int buf[], PVLayerLoc * loc)
{
   loc->nx       = buf[0];
   loc->ny       = buf[1];
   loc->nxGlobal = buf[2];
   loc->nyGlobal = buf[3];
   loc->kx0      = buf[4];
   loc->ky0      = buf[5];
   loc->nb       = buf[6];
   loc->nf       = buf[7];
   loc->halo.lt  = buf[8];
   loc->halo.rt  = buf[9];
   loc->halo.dn  = buf[10];
   loc->halo.up  = buf[11];
}

int getFileType(const char * filename)
{
   const char * ext = strrchr(filename, '.');
   if (strcmp(ext, ".pvp") == 0) {
      return PVP_FILE_TYPE;
   }
   return 0;
}

/**
 * Calculates location information given processor distribution and the
 * size of the image.
 *
 * @filename the name of the image file (in)
 * @ic the inter-column communicator (in)
 * @loc location information (inout) (loc->nx and loc->ny are out)
 */
int getImageInfo(const char * filename, PV::Communicator * comm, PVLayerLoc * loc)
{
   if (getFileType(filename) == PVP_FILE_TYPE) {
      return getImageInfoPVP(filename, comm, loc);
   }
   return getImageInfoGDAL(filename, comm, loc);
}

int getImageInfoPVP(const char * filename, PV::Communicator * comm, PVLayerLoc * loc)
{
   const int locSize = sizeof(PVLayerLoc) / sizeof(int);
   int locBuf[locSize];
   int status = 0;

   // LayerLoc should contain 8 ints
   assert(locSize == 8);

   const int icCol = comm->commColumn();
   const int icRow = comm->commRow();

#ifdef DEBUG_OUTPUT
   const int nxProcs = comm->numCommColumns();
   const int nyProcs = comm->numCommRows();
   fprintf(stderr, "[%2d]: nxProcs==%d nyProcs==%d icRow==%d icCol==%d\n",
           comm->commRank(), nxProcs, nyProcs, icRow, icCol);
#endif

   if (comm->commRank() == 0) {
      int numParams, params[NUM_PAR_BYTE_PARAMS];

      FILE * fp = fopen(filename, "rb");
      assert(fp != NULL);

      numParams = pv_read_binary_params(fp, NUM_PAR_BYTE_PARAMS, params);
      fclose(fp);

      assert(numParams == NUM_PAR_BYTE_PARAMS);
      assert(params[INDEX_FILE_TYPE] == PVP_FILE_TYPE);

      const int dataSize = params[INDEX_DATA_SIZE];
      const int dataType = params[INDEX_DATA_TYPE];
//      const int nxProcs  = params[INDEX_NX_PROCS];
//      const int nyProcs  = params[INDEX_NY_PROCS];

      loc->nx       = params[INDEX_NX];
      loc->ny       = params[INDEX_NY];
      loc->nxGlobal = params[INDEX_NX_GLOBAL];
      loc->nyGlobal = params[INDEX_NY_GLOBAL];
      loc->kx0      = params[INDEX_KX0];
      loc->ky0      = params[INDEX_KY0];
      loc->nb       = params[INDEX_NB];
      loc->nf       = params[INDEX_NF];

      assert(dataSize == 1);
      assert(dataType == PV_BYTE_TYPE);

      copyToLocBuffer(locBuf, loc);
   }

#ifdef PV_USE_MPI
   // broadcast location information
   // TODO - IMPORTANT - WHICH OF THESE IS CORRECT
   //MPI_Bcast(locBuf, 1+locSize, MPI_INT, 0, comm->communicator());
   MPI_Bcast(locBuf, locSize, MPI_INT, 0, comm->communicator());
#endif

   copyFromLocBuffer(locBuf, loc);

   // fix up layer indices
   loc->kx0 = loc->nx * icCol;
   loc->ky0 = loc->ny * icRow;

   return status;
}

int getImageInfoGDAL(const char * filename, PV::Communicator * comm, PVLayerLoc * loc)
{
   int status = 0;

#ifdef PV_USE_GDAL
   const int locSize = sizeof(PVLayerLoc) / sizeof(int);
   int locBuf[locSize];

   // LayerLoc should contain 12 ints
   assert(locSize == 12);

   const int nxProcs = comm->numCommColumns();
   const int nyProcs = comm->numCommRows();

   const int icCol = comm->commColumn();
   const int icRow = comm->commRow();

#ifdef DEBUG_OUTPUT
   fprintf(stderr, "[%2d]: nxProcs==%d nyProcs==%d icRow==%d icCol==%d\n",
           comm->commRank(), nxProcs, nyProcs, icRow, icCol);
#endif

   if (comm->commRank() == 0) {
      GDALAllRegister();

      GDALDataset * dataset = (GDALDataset *) GDALOpen(filename, GA_ReadOnly);
      if (dataset == NULL) return 1;

      int xImageSize = dataset->GetRasterXSize();
      int yImageSize = dataset->GetRasterYSize();

      loc->nf = dataset->GetRasterCount();

      // calculate local layer size

      int nx = xImageSize / nxProcs;
      int ny = yImageSize / nyProcs;

      loc->nx = nx;
      loc->ny = ny;

      loc->nxGlobal = nxProcs * nx;
      loc->nyGlobal = nyProcs * ny;

      copyToLocBuffer(locBuf, loc);

      GDALClose(dataset);
   }

#ifdef PV_USE_MPI
   // broadcast location information
   MPI_Bcast(locBuf, 1+locSize, MPI_INT, 0, comm->communicator());
#endif

   copyFromLocBuffer(locBuf, loc);

   // fix up layer indices
   loc->kx0 = loc->nx * icCol;
   loc->ky0 = loc->ny * icRow;
#else
   fprintf(stderr, GDAL_CONFIG_ERR_STR);
   exit(1);
#endif // PV_USE_GDAL

   return status;
}


int gatherImageFile(const char * filename,
                    PV::Communicator * comm, const PVLayerLoc * loc, pvdata_t * pvdata_buf){
   unsigned char * char_buf;
   const int numItems = loc->nx * loc->ny * loc->nf;
   char_buf = (unsigned char *) calloc(numItems, sizeof(unsigned char));
   assert( char_buf != NULL );
   pvdata_t max_buf = -1.0e20;
   pvdata_t min_buf = 1.0e20;
   for (int i = 0; i < numItems; i++) {
      max_buf = pvdata_buf[i] > max_buf ? pvdata_buf[i] : max_buf;
      min_buf = pvdata_buf[i] < min_buf ? pvdata_buf[i] : min_buf;
   }
   pvdata_t range_buf = max_buf - min_buf;  // all char_buf == 0
   if (range_buf == 0) {
      range_buf = 1.0;
   }
   for (int i = 0; i < numItems; i++) {
      char_buf[i] = 255 * ( pvdata_buf[i] - min_buf ) / range_buf;
   }
   int status = gatherImageFile(filename, comm, loc, char_buf);
   free(char_buf);
   return status;
}

int gatherImageFile(const char * filename,
                    PV::Communicator * comm, const PVLayerLoc * loc, unsigned char * buf)
{
   if (getFileType(filename) == PVP_FILE_TYPE) {
      return gatherImageFilePVP(filename, comm, loc, buf);
   }
   return gatherImageFileGDAL(filename, comm, loc, buf);
}

int gatherImageFilePVP(const char * filename,
                       PV::Communicator * comm, const PVLayerLoc * loc, unsigned char * buf)
{
   int status = 0;
   // const int maxBands = 3;

   const int nxProcs = comm->numCommColumns();
   const int nyProcs = comm->numCommRows();

   const int icRank = comm->commRank();

   const int nx = loc->nx;
   const int ny = loc->ny;

   const int numBands = loc->nf;
   // assert(numBands <= maxBands);

   // const int nxny     = nx * ny;
   const int numItems = nx * ny * numBands;

#ifdef PV_USE_MPI
   const int tag = PVP_FILE_TYPE;
   const MPI_Comm mpi_comm = comm->communicator();
#endif // PV_USE_MPI

   if (icRank > 0) {
#ifdef PV_USE_MPI
      const int dest = 0;
      MPI_Send(buf, numItems, MPI_BYTE, dest, tag, mpi_comm);
#ifdef DEBUG_OUTPUT
      fprintf(stderr, "[%2d]: gather: sent to 0, nx==%d ny==%d size==%d\n",
              comm->commRank(), nx, ny, nx*ny);
#endif
#endif // PV_USE_MPI
   }
   else {
      int params[NUM_PAR_BYTE_PARAMS];

      const int numParams  = NUM_PAR_BYTE_PARAMS;
      const int headerSize = numParams * sizeof(int);
      const int recordSize = numItems * sizeof(unsigned char);

      FILE * fp = fopen(filename, "wb");

      params[INDEX_HEADER_SIZE] = headerSize;
      params[INDEX_NUM_PARAMS]  = numParams;
      params[INDEX_FILE_TYPE]   = PVP_FILE_TYPE;
      params[INDEX_NX]          = loc->nx;
      params[INDEX_NY]          = loc->ny;
      params[INDEX_NF]          = loc->nf;
      params[INDEX_NUM_RECORDS] = nxProcs * nyProcs;
      params[INDEX_RECORD_SIZE] = recordSize;
      params[INDEX_DATA_SIZE]   = sizeof(unsigned char);
      params[INDEX_DATA_TYPE]   = PV_BYTE_TYPE;
      params[INDEX_NX_PROCS]    = nxProcs;
      params[INDEX_NY_PROCS]    = nyProcs;
      params[INDEX_NX_GLOBAL]   = loc->nxGlobal;
      params[INDEX_NY_GLOBAL]   = loc->nyGlobal;
      params[INDEX_KX0]         = loc->kx0;
      params[INDEX_KY0]         = loc->ky0;
      params[INDEX_NB]          = loc->nb;
      params[INDEX_NBANDS]      = loc->nf;

      int numWrite = fwrite(params, sizeof(int), numParams, fp);
      assert(numWrite == numParams);

      // write local image portion
      fseek(fp, (long) headerSize, SEEK_SET);
      numWrite = fwrite(buf, sizeof(unsigned char), numItems, fp);
      assert(numWrite == numItems);

#ifdef PV_USE_MPI
      int src = -1;
      unsigned char * icBuf = (unsigned char *) malloc(numItems * sizeof(unsigned char));
      assert(icBuf != NULL);
      for (int py = 0; py < nyProcs; py++) {
         for (int px = 0; px < nxProcs; px++) {
            if (++src == 0) continue;
#ifdef DEBUG_OUTPUT
            fprintf(stderr, "[%2d]: gather: receiving from %d nx==%d ny==%d nf==%d, numItems==%d\n",
                    comm->commRank(), src, nx, ny, numBands, numItems);
#endif
            MPI_Recv(icBuf, numItems, MPI_BYTE, src, tag, mpi_comm, MPI_STATUS_IGNORE);

            long offset = headerSize + src * recordSize;
            fseek(fp, offset, SEEK_SET);
            numWrite = fwrite(icBuf, sizeof(unsigned char), numItems, fp);
            assert(numWrite == numItems);
         }
      }
      free(icBuf);
      // TODO write to a file with params[INDEX_NX_PROCS] = params[INDEX_NY_PROCS] = 1
      // so that reading the file doesn't depend on having the same number of processes.
#endif // PV_USE_MPI

      status = fclose(fp);
   }

   return status;
}

int gatherImageFileGDAL(const char * filename,
                        PV::Communicator * comm, const PVLayerLoc * loc, unsigned char * buf)
{
   int status = 0;

#ifdef PV_USE_GDAL
   // const int maxBands = 3;

   const int nxProcs = comm->numCommColumns();
   const int nyProcs = comm->numCommRows();

   const int icRank = comm->commRank();

   const int nx = loc->nx;
   const int ny = loc->ny;

   const int numBands = loc->nf;
   // assert(numBands <= maxBands);

   const int nxnynf = nx * ny * numBands;

#ifdef PV_USE_MPI
   const int tag = 14;
   const MPI_Comm mpi_comm = comm->communicator();
#endif // PV_USE_MPI

   if (icRank > 0) {
#ifdef PV_USE_MPI
      const int dest = 0;

      MPI_Send(buf, nx*ny*numBands, MPI_BYTE, dest, tag, mpi_comm);
      // for (int b = 0; b < numBands; b++) {
      //   MPI_Send(&buf[b*nxny], nxny, MPI_BYTE, dest, tag, mpi_comm);
      // }
#ifdef DEBUG_OUTPUT
      fprintf(stderr, "[%2d]: gather: sent to 0, nx==%d ny==%d nf==%d size==%d\n",
              comm->commRank(), nx, ny, numBands, nxnynf);
#endif
#endif // PV_USE_MPI
   }
   else {
      GDALAllRegister();

      GDALDriver * driver = GetGDALDriverManager()->GetDriverByName("GTiff");
      if (driver == NULL) {
         exit(1);
      }

      int xImageSize = nx * nxProcs;
      int yImageSize = ny * nyProcs;

      GDALDataset * dataset = driver->Create(filename, xImageSize, yImageSize, numBands,
                                             GDT_Byte, NULL);

      if (dataset == NULL) {
          fprintf(stderr, "[%2d]: gather: failed to open file %s\n", comm->commRank(), filename);
      }
      else {
#ifdef DEBUG_OUTPUT
          fprintf(stderr, "[%2d]: gather: opened file %s\n", comm->commRank(), filename);
#endif
      }

      // GDALRasterBand * band[maxBands];

      assert(numBands <= dataset->GetRasterCount()); // Since dataset was created using numBands, what is this assert testing for?

      // GDALRasterBand ** band = (GDALRasterBand **) malloc( numBands * sizeof(GDALRasterBand *) );

      // for (int b = 0; b < numBands; b++) {
      //    band[b] = dataset->GetRasterBand(b+1);
      // }

      // write local image portion
      dataset->RasterIO(GF_Write, 0, 0, nx, ny, buf, nx, ny, GDT_Byte,
                        numBands, NULL,numBands, numBands*nx, 1);
      // for (int b = 0; b < numBands; b++) {
      //    band[b]->RasterIO(GF_Write, 0, 0, nx, ny,
      //                      &buf[b*nxny], nx, ny, GDT_Byte, 0, 0);
      // }

#ifdef PV_USE_MPI
      int src = -1;
      unsigned char * icBuf = (unsigned char *) malloc(nxnynf * sizeof(unsigned char));
      assert(icBuf != NULL);
      for (int py = 0; py < nyProcs; py++) {
         for (int px = 0; px < nxProcs; px++) {
            if (++src == 0) continue;
            int kx = nx * px;
            int ky = ny * py;
#ifdef DEBUG_OUTPUT
            fprintf(stderr, "[%2d]: gather: receiving from %d xSize==%d"
                    " ySize==%d size==%d total==%d\n",
                    comm->commRank(), src, nx, ny, nxny*numBands,
                    numTotal*comm->commSize());
#endif
            dataset->RasterIO(GF_Write, kx, ky, nx, ny, icBuf, nx, ny, GDT_Byte,
                              numBands, NULL, numBands, numBands*nx, 1);
            // for (int b = 0; b < numBands; b++) {
            //    MPI_Recv(icBuf, nxny, MPI_BYTE, src, tag, mpi_comm, MPI_STATUS_IGNORE);
            //    band[b]->RasterIO(GF_Write, kx, ky, nx, ny,
            //                      tmp, nx, ny, GDT_Byte, 0, 0);
            // }
         }
      }
      free(icBuf);
#endif // PV_USE_MPI
      GDALClose(dataset);
   }
#else
   fprintf(stderr, GDAL_CONFIG_ERR_STR);
   exit(1);
#endif // PV_USE_GDAL

   return status;
}

int scatterImageFile(const char * filename, int xOffset, int yOffset,
                     PV::Communicator * comm, PVLayerLoc * loc, unsigned char * buf)
{
   if (getFileType(filename) == PVP_FILE_TYPE) {
      // TODO - add xOffset and yOffset parameters
      return scatterImageFilePVP(filename, comm, loc, buf);
   }
   return scatterImageFileGDAL(filename, xOffset, yOffset, comm, loc, buf);
}

int scatterImageFilePVP(const char * filename,
                        PV::Communicator * comm, PVLayerLoc * loc, unsigned char * buf)
{
   int status = 0;
   // const int maxBands = 3;

   const int icRank = comm->commRank();

   const int nx = loc->nx;
   const int ny = loc->ny;

   const int numBands = loc->nf;
   // assert(numBands <= maxBands);

   // const int nxny     = nx * ny;
   const int numItems = nx * ny * numBands;

   if (icRank > 0) {
#ifdef PV_USE_MPI
      const int src = 0;
      const int tag = PVP_FILE_TYPE;
      const MPI_Comm mpi_comm = comm->communicator();

      MPI_Recv(buf, numItems, MPI_BYTE, src, tag, mpi_comm, MPI_STATUS_IGNORE);

#ifdef DEBUG_OUTPUT
      fprintf(stderr, "[%2d]: scatter: received from 0, nx==%d ny==%d nf==%d numItems==%d\n",
              comm->commRank(), nx, ny, numBands, numItems);
#endif
#endif // PV_USE_MPI
   }
   else {
      int params[NUM_PAR_BYTE_PARAMS];
      int numParams, numRead, type, nxIn, nyIn, nfIn;

      FILE * fp = pv_open_binary(filename, &numParams, &type, &nxIn, &nyIn, &nfIn);
      assert(fp != NULL);
      assert(numParams == NUM_PAR_BYTE_PARAMS);
      assert(type      == PVP_FILE_TYPE);

      status = pv_read_binary_params(fp, numParams, params);
      assert(status == numParams);

      const size_t headerSize = (size_t) params[INDEX_HEADER_SIZE];

      const int numRecords = params[INDEX_NUM_RECORDS];
      const int dataSize = params[INDEX_DATA_SIZE];
      const int dataType = params[INDEX_DATA_TYPE];
      const int nxProcs  = params[INDEX_NX_PROCS];
      const int nyProcs  = params[INDEX_NY_PROCS];

      loc->nx       = params[INDEX_NX];
      loc->ny       = params[INDEX_NY];
      loc->nxGlobal = params[INDEX_NX_GLOBAL];
      loc->nyGlobal = params[INDEX_NY_GLOBAL];
      loc->kx0      = params[INDEX_KX0];
      loc->ky0      = params[INDEX_KY0];
      loc->nb       = params[INDEX_NB];
      loc->nf       = params[INDEX_NF];

      assert(dataSize == 1);
      assert(dataType == PV_BYTE_TYPE);
      assert(nxProcs == comm->numCommColumns());
      assert(nyProcs == comm->numCommRows());
      assert(numRecords == comm->commSize());
      // TODO when in MPI, reading a file should not depend
      // on the number of processors being the same between
      // reading and writing

#ifdef PV_USE_MPI
      int dest = -1;
      const int tag = PVP_FILE_TYPE;
      const MPI_Comm mpi_comm = comm->communicator();
      const int recordSize = numItems * sizeof(unsigned char);

      for (int py = 0; py < nyProcs; py++) {
         for (int px = 0; px < nxProcs; px++) {
            if (++dest == 0) continue;

#ifdef DEBUG_OUTPUT
            fprintf(stderr, "[%2d]: scatter: sending to %d xSize==%d"
                    " ySize==%d size==%d total==%d\n",
                    comm->commRank(), dest, nx, ny, nx*ny,
                    nx*ny*comm->commSize());
#endif
            long offset = headerSize + dest * recordSize;
            fseek(fp, offset, SEEK_SET);
            numRead = fread(buf, sizeof(unsigned char), numItems, fp);
            assert(numRead == numItems);
            MPI_Send(buf, numItems, MPI_BYTE, dest, tag, mpi_comm);
         }
      }
#endif // PV_USE_MPI

      // get local image portion
      fseek(fp, (long) headerSize, SEEK_SET);
      numRead = fread(buf, sizeof(unsigned char), numItems, fp);
      assert(numRead == numItems);

      status = pv_close_binary(fp);
   }

   return status;
}

int scatterImageFileGDAL(const char * filename, int xOffset, int yOffset,
                         PV::Communicator * comm, PVLayerLoc * loc, unsigned char * buf)
{
   int status = 0;

#ifdef PV_USE_GDAL
   // const int maxBands = 3;

   const int nxProcs = comm->numCommColumns();
   const int nyProcs = comm->numCommRows();

   const int icRank = comm->commRank();

   const int nx = loc->nx;
   const int ny = loc->ny;

   const int numBands = loc->nf;

   // assert(numBands <= maxBands);

   // const int nxny = nx * ny;
   const int numTotal = nx * ny * numBands;

   if (icRank > 0) {
#ifdef PV_USE_MPI

      const int src = 0;
      const int tag = 13;
      const MPI_Comm mpi_comm = comm->communicator();

      MPI_Recv(buf, numTotal, MPI_BYTE, src, tag, mpi_comm, MPI_STATUS_IGNORE);
      // for (int b = 0; b < numBands; b++) {
      //    MPI_Recv(&buf[b*nxny], numTotal, MPI_BYTE, src, tag, mpi_comm, MPI_STATUS_IGNORE);
      // }
#ifdef DEBUG_OUTPUT
      fprintf(stderr, "[%2d]: scatter: received from 0, nx==%d ny==%d size==%d\n",
              comm->commRank(), nx, ny, numTotal);
#endif
#endif // PV_USE_MPI
   }
   else {<commit_msg>Added lines to ifdef PV_USE_MPI to eliminate compiler warnings<commit_after>#include "imageio.hpp"
#include "io.h"

#include <assert.h>
#include <string.h>

#ifdef PV_USE_GDAL
#  include <gdal_priv.h>
#  include <ogr_spatialref.h>
#else
#  define GDAL_CONFIG_ERR_STR "PetaVision must be compiled with GDAL to use this file type\n"
#endif // PV_USE_GDAL

#undef DEBUG_OUTPUT

static void copyToLocBuffer(int buf[], PVLayerLoc * loc)
{
   buf[0] = loc->nx;
   buf[1] = loc->ny;
   buf[2] = loc->nxGlobal;
   buf[3] = loc->nyGlobal;
   buf[4] = loc->kx0;
   buf[5] = loc->ky0;
   buf[6] = loc->nb;
   buf[7] = loc->nf;
   buf[8] = loc->halo.lt;
   buf[9] = loc->halo.rt;
   buf[10] = loc->halo.dn;
   buf[11] = loc->halo.up;
}

static void copyFromLocBuffer(int buf[], PVLayerLoc * loc)
{
   loc->nx       = buf[0];
   loc->ny       = buf[1];
   loc->nxGlobal = buf[2];
   loc->nyGlobal = buf[3];
   loc->kx0      = buf[4];
   loc->ky0      = buf[5];
   loc->nb       = buf[6];
   loc->nf       = buf[7];
   loc->halo.lt  = buf[8];
   loc->halo.rt  = buf[9];
   loc->halo.dn  = buf[10];
   loc->halo.up  = buf[11];
}

int getFileType(const char * filename)
{
   const char * ext = strrchr(filename, '.');
   if (strcmp(ext, ".pvp") == 0) {
      return PVP_FILE_TYPE;
   }
   return 0;
}

/**
 * Calculates location information given processor distribution and the
 * size of the image.
 *
 * @filename the name of the image file (in)
 * @ic the inter-column communicator (in)
 * @loc location information (inout) (loc->nx and loc->ny are out)
 */
int getImageInfo(const char * filename, PV::Communicator * comm, PVLayerLoc * loc)
{
   if (getFileType(filename) == PVP_FILE_TYPE) {
      return getImageInfoPVP(filename, comm, loc);
   }
   return getImageInfoGDAL(filename, comm, loc);
}

int getImageInfoPVP(const char * filename, PV::Communicator * comm, PVLayerLoc * loc)
{
   const int locSize = sizeof(PVLayerLoc) / sizeof(int);
   int locBuf[locSize];
   int status = 0;

   // LayerLoc should contain 8 ints
   assert(locSize == 8);

   const int icCol = comm->commColumn();
   const int icRow = comm->commRow();

#ifdef DEBUG_OUTPUT
   const int nxProcs = comm->numCommColumns();
   const int nyProcs = comm->numCommRows();
   fprintf(stderr, "[%2d]: nxProcs==%d nyProcs==%d icRow==%d icCol==%d\n",
           comm->commRank(), nxProcs, nyProcs, icRow, icCol);
#endif

   if (comm->commRank() == 0) {
      int numParams, params[NUM_PAR_BYTE_PARAMS];

      FILE * fp = fopen(filename, "rb");
      assert(fp != NULL);

      numParams = pv_read_binary_params(fp, NUM_PAR_BYTE_PARAMS, params);
      fclose(fp);

      assert(numParams == NUM_PAR_BYTE_PARAMS);
      assert(params[INDEX_FILE_TYPE] == PVP_FILE_TYPE);

      const int dataSize = params[INDEX_DATA_SIZE];
      const int dataType = params[INDEX_DATA_TYPE];
//      const int nxProcs  = params[INDEX_NX_PROCS];
//      const int nyProcs  = params[INDEX_NY_PROCS];

      loc->nx       = params[INDEX_NX];
      loc->ny       = params[INDEX_NY];
      loc->nxGlobal = params[INDEX_NX_GLOBAL];
      loc->nyGlobal = params[INDEX_NY_GLOBAL];
      loc->kx0      = params[INDEX_KX0];
      loc->ky0      = params[INDEX_KY0];
      loc->nb       = params[INDEX_NB];
      loc->nf       = params[INDEX_NF];

      assert(dataSize == 1);
      assert(dataType == PV_BYTE_TYPE);

      copyToLocBuffer(locBuf, loc);
   }

#ifdef PV_USE_MPI
   // broadcast location information
   // TODO - IMPORTANT - WHICH OF THESE IS CORRECT
   //MPI_Bcast(locBuf, 1+locSize, MPI_INT, 0, comm->communicator());
   MPI_Bcast(locBuf, locSize, MPI_INT, 0, comm->communicator());
#endif

   copyFromLocBuffer(locBuf, loc);

   // fix up layer indices
   loc->kx0 = loc->nx * icCol;
   loc->ky0 = loc->ny * icRow;

   return status;
}

int getImageInfoGDAL(const char * filename, PV::Communicator * comm, PVLayerLoc * loc)
{
   int status = 0;

#ifdef PV_USE_GDAL
   const int locSize = sizeof(PVLayerLoc) / sizeof(int);
   int locBuf[locSize];

   // LayerLoc should contain 12 ints
   assert(locSize == 12);

   const int nxProcs = comm->numCommColumns();
   const int nyProcs = comm->numCommRows();

   const int icCol = comm->commColumn();
   const int icRow = comm->commRow();

#ifdef DEBUG_OUTPUT
   fprintf(stderr, "[%2d]: nxProcs==%d nyProcs==%d icRow==%d icCol==%d\n",
           comm->commRank(), nxProcs, nyProcs, icRow, icCol);
#endif

   if (comm->commRank() == 0) {
      GDALAllRegister();

      GDALDataset * dataset = (GDALDataset *) GDALOpen(filename, GA_ReadOnly);
      if (dataset == NULL) return 1;

      int xImageSize = dataset->GetRasterXSize();
      int yImageSize = dataset->GetRasterYSize();

      loc->nf = dataset->GetRasterCount();

      // calculate local layer size

      int nx = xImageSize / nxProcs;
      int ny = yImageSize / nyProcs;

      loc->nx = nx;
      loc->ny = ny;

      loc->nxGlobal = nxProcs * nx;
      loc->nyGlobal = nyProcs * ny;

      copyToLocBuffer(locBuf, loc);

      GDALClose(dataset);
   }

#ifdef PV_USE_MPI
   // broadcast location information
   MPI_Bcast(locBuf, 1+locSize, MPI_INT, 0, comm->communicator());
#endif

   copyFromLocBuffer(locBuf, loc);

   // fix up layer indices
   loc->kx0 = loc->nx * icCol;
   loc->ky0 = loc->ny * icRow;
#else
   fprintf(stderr, GDAL_CONFIG_ERR_STR);
   exit(1);
#endif // PV_USE_GDAL

   return status;
}


int gatherImageFile(const char * filename,
                    PV::Communicator * comm, const PVLayerLoc * loc, pvdata_t * pvdata_buf){
   unsigned char * char_buf;
   const int numItems = loc->nx * loc->ny * loc->nf;
   char_buf = (unsigned char *) calloc(numItems, sizeof(unsigned char));
   assert( char_buf != NULL );
   pvdata_t max_buf = -1.0e20;
   pvdata_t min_buf = 1.0e20;
   for (int i = 0; i < numItems; i++) {
      max_buf = pvdata_buf[i] > max_buf ? pvdata_buf[i] : max_buf;
      min_buf = pvdata_buf[i] < min_buf ? pvdata_buf[i] : min_buf;
   }
   pvdata_t range_buf = max_buf - min_buf;  // all char_buf == 0
   if (range_buf == 0) {
      range_buf = 1.0;
   }
   for (int i = 0; i < numItems; i++) {
      char_buf[i] = 255 * ( pvdata_buf[i] - min_buf ) / range_buf;
   }
   int status = gatherImageFile(filename, comm, loc, char_buf);
   free(char_buf);
   return status;
}

int gatherImageFile(const char * filename,
                    PV::Communicator * comm, const PVLayerLoc * loc, unsigned char * buf)
{
   if (getFileType(filename) == PVP_FILE_TYPE) {
      return gatherImageFilePVP(filename, comm, loc, buf);
   }
   return gatherImageFileGDAL(filename, comm, loc, buf);
}

int gatherImageFilePVP(const char * filename,
                       PV::Communicator * comm, const PVLayerLoc * loc, unsigned char * buf)
{
   int status = 0;
   // const int maxBands = 3;

   const int nxProcs = comm->numCommColumns();
   const int nyProcs = comm->numCommRows();

   const int icRank = comm->commRank();

   const int nx = loc->nx;
   const int ny = loc->ny;

   const int numBands = loc->nf;
   // assert(numBands <= maxBands);

   // const int nxny     = nx * ny;
   const int numItems = nx * ny * numBands;

#ifdef PV_USE_MPI
   const int tag = PVP_FILE_TYPE;
   const MPI_Comm mpi_comm = comm->communicator();
#endif // PV_USE_MPI

   if (icRank > 0) {
#ifdef PV_USE_MPI
      const int dest = 0;
      MPI_Send(buf, numItems, MPI_BYTE, dest, tag, mpi_comm);
#ifdef DEBUG_OUTPUT
      fprintf(stderr, "[%2d]: gather: sent to 0, nx==%d ny==%d size==%d\n",
              comm->commRank(), nx, ny, nx*ny);
#endif
#endif // PV_USE_MPI
   }
   else {
      int params[NUM_PAR_BYTE_PARAMS];

      const int numParams  = NUM_PAR_BYTE_PARAMS;
      const int headerSize = numParams * sizeof(int);
      const int recordSize = numItems * sizeof(unsigned char);

      FILE * fp = fopen(filename, "wb");

      params[INDEX_HEADER_SIZE] = headerSize;
      params[INDEX_NUM_PARAMS]  = numParams;
      params[INDEX_FILE_TYPE]   = PVP_FILE_TYPE;
      params[INDEX_NX]          = loc->nx;
      params[INDEX_NY]          = loc->ny;
      params[INDEX_NF]          = loc->nf;
      params[INDEX_NUM_RECORDS] = nxProcs * nyProcs;
      params[INDEX_RECORD_SIZE] = recordSize;
      params[INDEX_DATA_SIZE]   = sizeof(unsigned char);
      params[INDEX_DATA_TYPE]   = PV_BYTE_TYPE;
      params[INDEX_NX_PROCS]    = nxProcs;
      params[INDEX_NY_PROCS]    = nyProcs;
      params[INDEX_NX_GLOBAL]   = loc->nxGlobal;
      params[INDEX_NY_GLOBAL]   = loc->nyGlobal;
      params[INDEX_KX0]         = loc->kx0;
      params[INDEX_KY0]         = loc->ky0;
      params[INDEX_NB]          = loc->nb;
      params[INDEX_NBANDS]      = loc->nf;

      int numWrite = fwrite(params, sizeof(int), numParams, fp);
      assert(numWrite == numParams);

      // write local image portion
      fseek(fp, (long) headerSize, SEEK_SET);
      numWrite = fwrite(buf, sizeof(unsigned char), numItems, fp);
      assert(numWrite == numItems);

#ifdef PV_USE_MPI
      int src = -1;
      unsigned char * icBuf = (unsigned char *) malloc(numItems * sizeof(unsigned char));
      assert(icBuf != NULL);
      for (int py = 0; py < nyProcs; py++) {
         for (int px = 0; px < nxProcs; px++) {
            if (++src == 0) continue;
#ifdef DEBUG_OUTPUT
            fprintf(stderr, "[%2d]: gather: receiving from %d nx==%d ny==%d nf==%d, numItems==%d\n",
                    comm->commRank(), src, nx, ny, numBands, numItems);
#endif
            MPI_Recv(icBuf, numItems, MPI_BYTE, src, tag, mpi_comm, MPI_STATUS_IGNORE);

            long offset = headerSize + src * recordSize;
            fseek(fp, offset, SEEK_SET);
            numWrite = fwrite(icBuf, sizeof(unsigned char), numItems, fp);
            assert(numWrite == numItems);
         }
      }
      free(icBuf);
      // TODO write to a file with params[INDEX_NX_PROCS] = params[INDEX_NY_PROCS] = 1
      // so that reading the file doesn't depend on having the same number of processes.
#endif // PV_USE_MPI

      status = fclose(fp);
   }

   return status;
}

int gatherImageFileGDAL(const char * filename,
                        PV::Communicator * comm, const PVLayerLoc * loc, unsigned char * buf)
{
   int status = 0;

#ifdef PV_USE_GDAL
   // const int maxBands = 3;

   const int nxProcs = comm->numCommColumns();
   const int nyProcs = comm->numCommRows();

   const int icRank = comm->commRank();

   const int nx = loc->nx;
   const int ny = loc->ny;

   const int numBands = loc->nf;
   // assert(numBands <= maxBands);

#ifdef PV_USE_MPI
   const int nxnynf = nx * ny * numBands;

   const int tag = 14;
   const MPI_Comm mpi_comm = comm->communicator();
#endif // PV_USE_MPI

   if (icRank > 0) {
#ifdef PV_USE_MPI
      const int dest = 0;

      MPI_Send(buf, nx*ny*numBands, MPI_BYTE, dest, tag, mpi_comm);
      // for (int b = 0; b < numBands; b++) {
      //   MPI_Send(&buf[b*nxny], nxny, MPI_BYTE, dest, tag, mpi_comm);
      // }
#ifdef DEBUG_OUTPUT
      fprintf(stderr, "[%2d]: gather: sent to 0, nx==%d ny==%d nf==%d size==%d\n",
              comm->commRank(), nx, ny, numBands, nxnynf);
#endif
#endif // PV_USE_MPI
   }
   else {
      GDALAllRegister();

      GDALDriver * driver = GetGDALDriverManager()->GetDriverByName("GTiff");
      if (driver == NULL) {
         exit(1);
      }

      int xImageSize = nx * nxProcs;
      int yImageSize = ny * nyProcs;

      GDALDataset * dataset = driver->Create(filename, xImageSize, yImageSize, numBands,
                                             GDT_Byte, NULL);

      if (dataset == NULL) {
          fprintf(stderr, "[%2d]: gather: failed to open file %s\n", comm->commRank(), filename);
      }
      else {
#ifdef DEBUG_OUTPUT
          fprintf(stderr, "[%2d]: gather: opened file %s\n", comm->commRank(), filename);
#endif
      }

      // GDALRasterBand * band[maxBands];

      assert(numBands <= dataset->GetRasterCount()); // Since dataset was created using numBands, what is this assert testing for?

      // GDALRasterBand ** band = (GDALRasterBand **) malloc( numBands * sizeof(GDALRasterBand *) );

      // for (int b = 0; b < numBands; b++) {
      //    band[b] = dataset->GetRasterBand(b+1);
      // }

      // write local image portion
      dataset->RasterIO(GF_Write, 0, 0, nx, ny, buf, nx, ny, GDT_Byte,
                        numBands, NULL,numBands, numBands*nx, 1);
      // for (int b = 0; b < numBands; b++) {
      //    band[b]->RasterIO(GF_Write, 0, 0, nx, ny,
      //                      &buf[b*nxny], nx, ny, GDT_Byte, 0, 0);
      // }

#ifdef PV_USE_MPI
      int src = -1;
      unsigned char * icBuf = (unsigned char *) malloc(nxnynf * sizeof(unsigned char));
      assert(icBuf != NULL);
      for (int py = 0; py < nyProcs; py++) {
         for (int px = 0; px < nxProcs; px++) {
            if (++src == 0) continue;
            int kx = nx * px;
            int ky = ny * py;
#ifdef DEBUG_OUTPUT
            fprintf(stderr, "[%2d]: gather: receiving from %d xSize==%d"
                    " ySize==%d size==%d total==%d\n",
                    comm->commRank(), src, nx, ny, nxny*numBands,
                    numTotal*comm->commSize());
#endif
            dataset->RasterIO(GF_Write, kx, ky, nx, ny, icBuf, nx, ny, GDT_Byte,
                              numBands, NULL, numBands, numBands*nx, 1);
            // for (int b = 0; b < numBands; b++) {
            //    MPI_Recv(icBuf, nxny, MPI_BYTE, src, tag, mpi_comm, MPI_STATUS_IGNORE);
            //    band[b]->RasterIO(GF_Write, kx, ky, nx, ny,
            //                      tmp, nx, ny, GDT_Byte, 0, 0);
            // }
         }
      }
      free(icBuf);
#endif // PV_USE_MPI
      GDALClose(dataset);
   }
#else
   fprintf(stderr, GDAL_CONFIG_ERR_STR);
   exit(1);
#endif // PV_USE_GDAL

   return status;
}

int scatterImageFile(const char * filename, int xOffset, int yOffset,
                     PV::Communicator * comm, PVLayerLoc * loc, unsigned char * buf)
{
   if (getFileType(filename) == PVP_FILE_TYPE) {
      // TODO - add xOffset and yOffset parameters
      return scatterImageFilePVP(filename, comm, loc, buf);
   }
   return scatterImageFileGDAL(filename, xOffset, yOffset, comm, loc, buf);
}

int scatterImageFilePVP(const char * filename,
                        PV::Communicator * comm, PVLayerLoc * loc, unsigned char * buf)
{
   int status = 0;
   // const int maxBands = 3;

   const int icRank = comm->commRank();

   const int nx = loc->nx;
   const int ny = loc->ny;

   const int numBands = loc->nf;
   // assert(numBands <= maxBands);

   // const int nxny     = nx * ny;
   const int numItems = nx * ny * numBands;

   if (icRank > 0) {
#ifdef PV_USE_MPI
      const int src = 0;
      const int tag = PVP_FILE_TYPE;
      const MPI_Comm mpi_comm = comm->communicator();

      MPI_Recv(buf, numItems, MPI_BYTE, src, tag, mpi_comm, MPI_STATUS_IGNORE);

#ifdef DEBUG_OUTPUT
      fprintf(stderr, "[%2d]: scatter: received from 0, nx==%d ny==%d nf==%d numItems==%d\n",
              comm->commRank(), nx, ny, numBands, numItems);
#endif
#endif // PV_USE_MPI
   }
   else {
      int params[NUM_PAR_BYTE_PARAMS];
      int numParams, numRead, type, nxIn, nyIn, nfIn;

      FILE * fp = pv_open_binary(filename, &numParams, &type, &nxIn, &nyIn, &nfIn);
      assert(fp != NULL);
      assert(numParams == NUM_PAR_BYTE_PARAMS);
      assert(type      == PVP_FILE_TYPE);

      status = pv_read_binary_params(fp, numParams, params);
      assert(status == numParams);

      const size_t headerSize = (size_t) params[INDEX_HEADER_SIZE];

      const int numRecords = params[INDEX_NUM_RECORDS];
      const int dataSize = params[INDEX_DATA_SIZE];
      const int dataType = params[INDEX_DATA_TYPE];
      const int nxProcs  = params[INDEX_NX_PROCS];
      const int nyProcs  = params[INDEX_NY_PROCS];

      loc->nx       = params[INDEX_NX];
      loc->ny       = params[INDEX_NY];
      loc->nxGlobal = params[INDEX_NX_GLOBAL];
      loc->nyGlobal = params[INDEX_NY_GLOBAL];
      loc->kx0      = params[INDEX_KX0];
      loc->ky0      = params[INDEX_KY0];
      loc->nb       = params[INDEX_NB];
      loc->nf       = params[INDEX_NF];

      assert(dataSize == 1);
      assert(dataType == PV_BYTE_TYPE);
      assert(nxProcs == comm->numCommColumns());
      assert(nyProcs == comm->numCommRows());
      assert(numRecords == comm->commSize());
      // TODO when in MPI, reading a file should not depend
      // on the number of processors being the same between
      // reading and writing

#ifdef PV_USE_MPI
      int dest = -1;
      const int tag = PVP_FILE_TYPE;
      const MPI_Comm mpi_comm = comm->communicator();
      const int recordSize = numItems * sizeof(unsigned char);

      for (int py = 0; py < nyProcs; py++) {
         for (int px = 0; px < nxProcs; px++) {
            if (++dest == 0) continue;

#ifdef DEBUG_OUTPUT
            fprintf(stderr, "[%2d]: scatter: sending to %d xSize==%d"
                    " ySize==%d size==%d total==%d\n",
                    comm->commRank(), dest, nx, ny, nx*ny,
                    nx*ny*comm->commSize());
#endif
            long offset = headerSize + dest * recordSize;
            fseek(fp, offset, SEEK_SET);
            numRead = fread(buf, sizeof(unsigned char), numItems, fp);
            assert(numRead == numItems);
            MPI_Send(buf, numItems, MPI_BYTE, dest, tag, mpi_comm);
         }
      }
#endif // PV_USE_MPI

      // get local image portion
      fseek(fp, (long) headerSize, SEEK_SET);
      numRead = fread(buf, sizeof(unsigned char), numItems, fp);
      assert(numRead == numItems);

      status = pv_close_binary(fp);
   }

   return status;
}

int scatterImageFileGDAL(const char * filename, int xOffset, int yOffset,
                         PV::Communicator * comm, PVLayerLoc * loc, unsigned char * buf)
{
   int status = 0;

#ifdef PV_USE_GDAL
   // const int maxBands = 3;

   const int nxProcs = comm->numCommColumns();
   const int nyProcs = comm->numCommRows();

   const int icRank = comm->commRank();

   const int nx = loc->nx;
   const int ny = loc->ny;

   const int numBands = loc->nf;

   // assert(numBands <= maxBands);

#ifdef PV_USE_MPI
   // const int nxny = nx * ny;
   const int numTotal = nx * ny * numBands;
#endif

   if (icRank > 0) {
#ifdef PV_USE_MPI

      const int src = 0;
      const int tag = 13;
      const MPI_Comm mpi_comm = comm->communicator();

      MPI_Recv(buf, numTotal, MPI_BYTE, src, tag, mpi_comm, MPI_STATUS_IGNORE);
      // for (int b = 0; b < numBands; b++) {
      //    MPI_Recv(&buf[b*nxny], numTotal, MPI_BYTE, src, tag, mpi_comm, MPI_STATUS_IGNORE);
      // }
#ifdef DEBUG_OUTPUT
      fprintf(stderr, "[%2d]: scatter: received from 0, nx==%d ny==%d size==%d\n",
              comm->commRank(), nx, ny, numTotal);
#endif
#endif // PV_USE_MPI
   }
   else {
{'fraction_non_alphanumeric': 0.13613544114505446, 'fraction_numerical': 0.008117923520615253, 'mean_word_length': 2.5591142368371034, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 289, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3349005', 'n_tokens_mistral': 13669, 'n_tokens_neox': 13095, 'n_words': 3281}