/**
   * 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
   */
  
  /*
  #define CAPACITY(store) store[0]
  #define NTH_LIST(store, index) (store + 1 + index * (CAPACITY(store) + 1))
  #define LIST_SIZE(list) list[0]
  #define LIST_GET(list, index) list[1 + index]
  #define LIST_SET(list, index, val) list[1 + index] = val
  #define LIST_APPEND(list, val) LIST_SET(list, atomic_inc(list), val)
  */
  
  /*
   * CSR representation (for 15 lists of max 10 items)
   * 10 2 1 5 0 0 0 0 0 0 0 0
   *    1 5 0 0 0 0 0 0 0 0 0
   *    5 8 7 3 1 4 0 0 0 0 0
   *    0 0 0 0 0 0 0 0 0 0 0
   */
   
   /* CSR Representation
    */
    
  /*
  #define NTH_CSR(store, iSim)
  #define NTH_LIST(store, index)
  */
  
  #define ACCESSIBLE(topo, iMM, iOM) (topo[iMM * world->nbOM + iOM] != -1)
  
  /**
   * KERNELS
   */
  
  /*
   * Dimensions expectations
   *
   * Global : nbSim * nbMM, nbOM
   */
  kernel void topology(
      global const MM mmList[],
      global const OM omList[],
      global const MiorWorld * world,
      global int topo[])
  {
      const int iSim = get_group_id(0) / world->nbMM;
      const int iMM = get_global_id(0) % world->nbMM;
      const int iOM = get_global_id(1);
      const int nbMM = world->nbMM;
      const int nbOM = world->nbOM;
      
      const int mmIndex = iSim * nbMM + iMM;
      const int omIndex = iSim * nbOM + iOM;
      
      const float dist = hypot(
          omList[omIndex].x - mmList[mmIndex].x,
          omList[omIndex].y - mmList[mmIndex].y
      );
      
      if (dist <= world->RA) {
          topo[mmIndex * nbOM + iOM] = 0;
  	}
  }
  
  /**
   * Dimensions expectations
   *
   * Parameters MUST already include the offset for this simulation
   * Local  : nbOM
   */
  kernel void scatter(
      global OM omList[],
      global const MiorWorld * world,
      global const int topo[],
      global int parts[])
  {
      // Retrieve current OM index
      const int iOM  = get_local_id(0);
      
      // Number of parts distributed
      int nbParts = 0;
      
      // Compute number of parts parts
      for (int iMM = 0; iMM < world->nbMM; iMM++) {
          if (ACCESSIBLE(topo, iMM, iOM)) {
              nbParts ++;
          }
      }
      
      // Compute the carbon part to allocate to each MM
      const int carbonPart = (world->K * omList[iOM].carbon) / nbParts;
      
      // Allocate parts
      for (int iMM = 0; iMM < world->nbMM; iMM++) {
          if (ACCESSIBLE(topo, iMM, iOM)) {
              parts[iMM * world->nbOM + iOM] = carbonPart;
          }
      }
      
      omList[iOM].carbon -= nbParts * carbonPart;
  }
  
  /**
   * Dimensions expectations
   *
   * Parameters MUST already include the offset for this simulation
   * Local  : nbMM
   */
  kernel void live(
      global MM mmList[],
      global MiorWorld * world,
      global int topo[],
      global int parts[])
  {
      const int iMM = get_local_id(0);
      global MM * currentMM = mmList + iMM;
      
      // Compute needs
      const int breathNeed = currentMM->carbon * world->RR;
      const int growthNeed = currentMM->carbon * world->GR;
      
      // BREATHING CHECK
      
      int iOM = 0, nbOM = world->nbOM;
      int remainingNeed = breathNeed;
      
      while (remainingNeed > 0 && iOM < nbOM) {
          if (ACCESSIBLE(topo, iMM, iOM)) {
              remainingNeed -= parts[iMM * world->nbOM + iOM];
          }
          iOM++;
      }
      
      // DORMANCY CHECK
      
      if (remainingNeed > 0) {
          currentMM->dormancy = 1;
          return;
      }
      
      // ACTUAL BREATHING
      
      iOM = 0;
      remainingNeed = breathNeed;
      currentMM->dormancy = 0;
      
      while (iOM < nbOM) {
          if (ACCESSIBLE(topo, iMM, iOM)) {
              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 iOM !
              }
          }
          
          iOM++;
      }
      
      // ENVIRONMENT UPDATE
      
      atomic_add(&(world->CO2), breathNeed);
      
      // GROWTH
      
      remainingNeed = growthNeed;
      
      while (remainingNeed > 0 && iOM < nbOM) {
          if (ACCESSIBLE(topo, iMM, iOM)) {
              const int offer = parts[iMM * world->nbOM + iOM];
              const int consum = min(offer, remainingNeed);
              
              remainingNeed -= consum;
              parts[iMM * world->nbOM + iOM] = offer - consum;
          }
          
          iOM++;
      }
          
      // MM CARBON UPDATE
      
      currentMM->carbon += (growthNeed - remainingNeed);
  }
      
      
  /**
   * Dimensions expectations
   *
   * Local  : nbOM
   */
  kernel void gather(
      global OM omList[],
      global const MiorWorld * world,
      global const int topo[],
      global const int parts[])
  {
      // Retrieve current OM index
      const int iOM = get_local_id(0);
      
      // Update current OM carbon value
      int carbon = omList[iOM].carbon;
      
      for (int iMM = 0; iMM < world->nbMM; iMM++) {
          if (ACCESSIBLE(topo, iMM, iOM)) {
              carbon += parts[iMM * world->nbOM + iOM];
          }
      }
      
      omList[iOM].carbon = carbon;
  }
  
  /**
   * Dimensions expectations
   *
   * Global : nbSim * max(nbOM, nbMM)
   * Local  : max(nbOM, nbOM)
   */
  kernel void autolive(
      global MM mmLists[],
      global OM omLists[],
      global MiorWorld * worldList,
      global const int topoList[],
      global int partsList[])
  {
      const int iSim = get_group_id(0);
      const int iAgent = get_local_id(0);
      const int nbMM = worldList->nbMM;
      const int nbOM = worldList->nbOM;
      
      // Compute parameters offset of this simulation
      global MM * mmList = mmLists + (iSim * nbMM);
      global OM * omList = omLists + (iSim * nbOM);
      global MiorWorld * world = worldList + iSim;
      global const int * topo = topoList + (iSim * nbMM * nbOM);
      global int * parts = partsList + (iSim * nbMM * nbOM);
      
      int CO2total = -1;
      
      //for (int i = 0; i < 1; i++) {
      while (CO2total == -1 || CO2total != world->CO2) {
          CO2total = world->CO2;
          
          if (iAgent < world->nbOM) {
              scatter(omList, world, topo, parts);
          }
          
          barrier(CLK_GLOBAL_MEM_FENCE);
          
          if (iAgent < world->nbMM) {
              live(mmList, world, topo, parts);
          }
          
          barrier(CLK_GLOBAL_MEM_FENCE);
          
          if (iAgent < world->nbOM) {
              gather(omList, world, topo, parts);
          }
          
          barrier(CLK_GLOBAL_MEM_FENCE);
      }
  
  }