/**
   * pre-OpenCL 1.1 compatibility macros
  */
  
  #if __OPENCL_VERSION__ <= CL_VERSION_1_0
      #define atomic_add(p, v) atom_add(p, v)
      #define atomic_inc(p) atom_inc(p)
  #endif
  
  /**
   * DATA STRUCTURES
   */
   
  // Microbial colony
  typedef struct MM {
      float   x;
      float   y;
      int     carbon;
      int     dormancy;
  } MM;
  
  // Carbon deposit
  typedef struct OM {
      float   x;
      float   y;
      int     carbon;
      int     lock;
  } OM;
  
  // Environment
  typedef struct MiorWorld {
      int     nbMM;
      int     nbOM;
      int     RA;
      float   RR;
      float   GR;
      float   K;
      int     width;
      int     minSize;
      int     CO2;
      int     lock;
  } MiorWorld;
  
  /**
   * HELPERS
   */
  
  /**
   * BUFFERS FORMAT EXPECTIONS
   * 
   * mmList    : static array of nbMM MM structs
   * mmList(n) : n x mmList
   *
   * omList    : static array of nbOM OM structs
   * omList(n) : n x omList
   * 
   * world(n)  : static array of n MiorWorld structs
   *
   * mmCSR     : static array of nbOM x (nbMM + 1) integers => items are MM indexes
   * mmCSR(n)  : n x mmCSR
   *
   * omCSR     : static array of nbMM x (nbOM + 1) integers => items are OM indexes
   * omCSR(n)  : n x omCSR
   *
   * topo      : static array of nbOM * nbMM integers
   * topo(n)   : n x topo 
   */
  
  /* STRUCTURE SIZE DEFINITIONS */
  
  #define PARTS_SIZE   (NB_MM * NB_OM)
  #define TOPO_SIZE    (NB_MM * NB_OM)
  #define MMLIST_SIZE  (NB_MM + 1)
  #define OMLIST_SIZE  (NB_OM + 1)
  #define MMCSR_SIZE   (NB_OM * MMLIST_SIZE)
  #define OMCSR_SIZE   (NB_MM * OMLIST_SIZE)
  
  /* HELPER FUNCTIONS */
  
  #define NTH_OMLIST(omCSR, index) (omCSR + index * OMLIST_SIZE)
  #define NTH_MMLIST(mmCSR, index) (mmCSR + index * MMLIST_SIZE)
  
  #define LIST_SIZE(list)         list[0]
  #define LIST_GET(list, index)   list[1 + index]
  #define LIST_APPEND(list, item) list[1 + atomic_inc(list)] = item
  
  /**
   * KERNELS
   */
  
  /*
   * Dimensions expectations
   *
   * //Global : nbMM, nbSim * nbOM
   * //Local  : 1, nbOM
   * Global : nbSim * nbMM, nbOM
   */
  kernel void topology(
      global const MM mmList[],
      global const OM omList[],
      global const MiorWorld * world,
      global int mmCSR[],
      global int omCSR[]/*,
      global int topo[]*/)
  {
      const int iSim = get_global_id(0) / NB_MM;
      const int iMM  = get_global_id(0) % NB_MM;
      const int iOM  = get_global_id(1);
      //const int iSim = get_group_id(1);
      //const int iMM  = get_global_id(0);
      //const int iOM  = get_local_id(1);
      const int nbMM = world->nbMM;
      const int nbOM = world->nbOM;
      
      if (iMM > nbMM || iOM > nbOM) {
      	return;
      }
      
      const int mmIndex = iSim * NB_MM + iMM;
      const int omIndex = iSim * NB_OM + iOM;
      
      const float dist = hypot(
          omList[omIndex].x - mmList[mmIndex].x,
          omList[omIndex].y - mmList[mmIndex].y
      );
      
      if (dist <= world->RA) {
          LIST_APPEND(NTH_OMLIST(omCSR, mmIndex), iOM);
          LIST_APPEND(NTH_MMLIST(mmCSR, omIndex), iMM);
          //topo[mmIndex * nbOM + iOM] = 0;
  	}
  }
  
  /**
   * Dimensions expectations
   *
   * Global : nbSim * nbOM
   * Local  : nbOM
   */
  void scatter(
      global OM omList[],
      global const MiorWorld * world,
      global const int mmCSR[],
      global int parts[])
  {
      // Retrieve current OM index
      const int iOM = get_local_id(0);
      
      // Retrieve list and number of accessible MM (CSR)
      global const int * const mmIndexes = NTH_MMLIST(mmCSR, iOM);
      const int nbIndexes = LIST_SIZE(mmIndexes);
      
      // Compute the carbon part to allocate to each MM
      const int carbonPart = (world->K * omList[iOM].carbon) / nbIndexes;
      
      // Allocate parts
      for (int i = 0; i < nbIndexes; i++) {
          const int iMM = LIST_GET(mmIndexes, i);       
          parts[iMM * world->nbOM + iOM] = carbonPart;
      }
      
      omList[iOM].carbon -= nbIndexes * carbonPart;
  }
  
  /**
   * Dimensions expectations
   *
   * Global : nbSim * nbMM
   * Local  : nbMM
   */
  void live(
      global MM mmList[],
      global MiorWorld * world,
      global const int omCSR[],
      global int parts[])
  {
      // Retrieve current MM index
      const int iMM  = get_local_id(0);
      global MM * const currentMM = mmList + iMM;
      
      // Compute needs
      const int breathNeed = currentMM->carbon * world->RR;
      const int growthNeed = currentMM->carbon * world->GR;
      
      // Retrieve list and number of accessible OM (CSR)
      global const int * const omIndexes = NTH_OMLIST(omCSR, iMM);
      const int nbIndexes = LIST_SIZE(omIndexes);
      
      // BREATHING CHECK
      
      int i = 0;
      int remainingNeed = breathNeed;
      
      while (remainingNeed > 0 && i < nbIndexes) {
          const int iOM = LIST_GET(omIndexes, i);
          remainingNeed -= parts[iMM * world->nbOM + iOM];
          i++;
      }
      
      // DORMANCY CHECK
      
      if (remainingNeed > 0) {
          currentMM->dormancy = 1;
          return;
      }
      
      // ACTUAL BREATHING
      
      i = 0;
      remainingNeed = breathNeed;
      currentMM->dormancy = 0;
      
      while (i < nbIndexes) {
          const int iOM = LIST_GET(omIndexes, i);
          const int offer = parts[iMM * world->nbOM + iOM];
          const int consum = min(offer, remainingNeed);
          
          remainingNeed -= consum;
          parts[iMM * world->nbOM + iOM] = offer - consum;
          
          if (remainingNeed == 0) {
              break; // don't increment i !
          }
          
          i++;
      }
      
      // ENVIRONMENT UPDATE
      
      atomic_add(&(world->CO2), breathNeed);
      
      // GROWTH
      
      remainingNeed = growthNeed;
      
      while (remainingNeed > 0 && i < nbIndexes) {
          const int iOM = LIST_GET(omIndexes, i);
          
          const int offer = parts[iMM * world->nbOM + iOM];
          const int consum = min(offer, remainingNeed);
          
          remainingNeed -= consum;
          parts[iMM * world->nbOM + iOM] = offer - consum;
          
          i++;
      }
          
      // MM CARBON UPDATE
      
      currentMM->carbon += (growthNeed - remainingNeed);
  }
      
      
  /**
   * Dimensions expectations
   *
   * Global : nbSim * nbOM
   * Local  : nbOM
   */
  void gather(
      global OM omList[],
      global const MiorWorld * world,
      global const int mmCSR[],
      global const int parts[])
  {
      // Retrieve current OM index
      const int iOM = get_local_id(0);
      
      // Retrieve list and number of accessible MM (CSR)
      global const int * const mmIndexes = NTH_MMLIST(mmCSR, iOM);
      const int nbIndexes = LIST_SIZE(mmIndexes);
      
      // Update current OM carbon value
      int carbon = omList[iOM].carbon;
      
      for (int i = 0; i < nbIndexes; i++) {
          const int iMM = LIST_GET(mmIndexes, i);
          carbon += parts[iMM * world->nbOM + iOM];
      }
      
      omList[iOM].carbon = carbon;
  }
  
  /*
   * Dimensions expectations
   * 
   * Global : nbSim * max(nbOM, nbMM)
   * Local  : max(nbOM, nbOM)
   */
  
  kernel void simulate(
      global MM           allMMList[],
      global OM           allOMList[],
      global MiorWorld    allWorlds[],
      global const int    allMMCSR[],
      global const int    allOMCSR[],
      global int          allParts[])
  {
      const int iSim = get_group_id(0);
      const int i    = get_local_id(0);
      
      global MiorWorld * world = allWorlds + iSim;
      global MM * mmList       = allMMList + iSim * NB_MM;
      global OM * omList       = allOMList + iSim * NB_OM;
      global const int * mmCSR = allMMCSR  + iSim * MMCSR_SIZE;
      global const int * omCSR = allOMCSR  + iSim * OMCSR_SIZE;
      global int * parts       = allParts  + iSim * PARTS_SIZE;
      
      //barrier(CLK_GLOBAL_MEM_FENCE);
      
      if (i < world->nbOM) {
          scatter(omList, world, mmCSR, parts);
      }
      
      barrier(CLK_GLOBAL_MEM_FENCE);
      
      if (i < world->nbMM) {
          live(mmList, world, omCSR, parts);
      }
      
      barrier(CLK_GLOBAL_MEM_FENCE);
      
      if (i < world->nbOM) {
          gather(omList, world, mmCSR, parts);
      }
      
      barrier(CLK_GLOBAL_MEM_FENCE);
  }
      
      
  /**
   * Dimensions expectations
   *
   * Global : nbSim * max(nbOM, nbMM)
   * Local  : max(nbOM, nbMM)
   */
  kernel void autolive(
      global MM allMMList[],
      global OM allOMList[],
      global MiorWorld allWorlds[],
      global const int allMMCSR[],
      global const int allOMCSR[],
      global int allParts[])
  {
      const int iSim = get_group_id(0);
      const int i    = get_local_id(0);
      int CO2total   = -1;
      
      global MiorWorld * world = allWorlds + iSim;
       
      while (CO2total == -1 || CO2total != world->CO2) {
          CO2total = world->CO2;
          simulate(allMMList, allOMList, allWorlds, allMMCSR, allOMCSR, allParts);
      }
  }