/**
   * 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)
  
  /**
   * KERNELS
   */
  
  /*
   * Dimensions expectations
   *
   * 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) / 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) {
          LIST_APPEND(NTH_LIST(omCSR, mmIndex), omIndex);
          LIST_APPEND(NTH_LIST(mmCSR, omIndex), mmIndex);
          topo[mmIndex * nbOM + iOM] = 0;
  	}
  }
  
  /**
   * Dimensions expectations
   *
   * Global : nbSim, nbOM
   */
  kernel void scatter(
      global OM omList[],
      global const MiorWorld * world,
      global const int mmCSR[],
      global int parts[])
  {
      // Retrieve current OM index
      // const int iSim = get_global_id(0);
      // const int iLocal = get_global_id(1);
      const int iOM = get_global_id(0) * world->nbOM + get_global_id(1);
      
      // Retrieve list and number of accessible MM (CSR)
      global const int * mmIndexes = NTH_LIST(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
   */
  kernel void live(
      global MM mmList[],
      global MiorWorld * world,
      global int omCSR[],
      global int parts[])
  {
      // const int iSim = get_global_id(0);
      // const int iLocal = get_global_id(1);
      const int iMM  = get_global_id(0) * world->nbMM + get_global_id(1);
      
      global MM * currentMM = mmList + iMM;
      //int localParts[NBOM];
      
      /*
      for (int i = 0; i < 310; i++) {
          localParts[i] = parts[iMM * world->nbOM + i];
      }
      */
      // 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 int * omIndexes = NTH_LIST(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];
          //remainingNeed -= localParts[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 = localParts[iOM];
          const int offer = parts[iMM * world->nbOM + iOM];
          const int consum = min(offer, remainingNeed);
          
          remainingNeed -= consum;
          parts[iMM * world->nbOM + iOM] = offer - consum;
          //localParts[iOM] = offer - consum;
          
          if (remainingNeed == 0) {
              //localParts[iOM] = offer - consum;
              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 offer = localParts[iOM];
          const int consum = min(offer, remainingNeed);
          
          remainingNeed -= consum;
          parts[iMM * world->nbOM + iOM] = offer - consum;
          //localParts[iOM] = offer - consum;
          
          i++;
      }
          
      // MM CARBON UPDATE
      
      currentMM->carbon += (growthNeed - remainingNeed);
      
      /*for (int i = 0; i < 310; i++) {
          parts[iMM * world->nbOM + i] = localParts[i];
      }*/
  }
      
      
  /**
   * Dimensions expectations
   *
   * Global : nbSim, nbOM
   */
  kernel void gather(
      global OM omList[],
      global const MiorWorld * world,
      global const int mmCSR[],
      global const int parts[])
  {
      // Retrieve current OM index
      // const int iSim = get_global_id(0);
      // const int iLocal = get_global_id(1);
      const int iOM = get_global_id(0) * world->nbOM + get_global_id(1);
      
      // Retrieve list and number of accessible MM (CSR)
      global const int * mmIndexes = NTH_LIST(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)
   */
  kernel void autolive(
      global MM mmList[],
      global OM omList[],
      global MiorWorld * world,
      global const int mmCSR[],
      global const int omCSR[],
      global int parts[])
  {
      const int i = get_global_id(1);
      int CO2total = -1;
      
      while (CO2total == -1 || CO2total != world->CO2) {
          if (i < world->nbOM) {
              scatter(omList, world, mmCSR, parts);
          }
          
          barrier(CLK_GLOBAL_MEM_FENCE);
          
          if (i < world->nbMM) {
              live(mmList, world, omCSR, parts);
          }
          
          if (i < world->nbOM) {
              gather(omList, world, mmCSR, parts);
          }
          
          barrier(CLK_GLOBAL_MEM_FENCE);
          
          //if (CO2total == world->CO2) break;
          CO2total = world->CO2;
      }
  
  }