/**
   * DATA STRUCTURES
   */
  
  typedef struct MM {
      float   x;
      float   y;
      int     carbone;
      int     dormance;
  } MM;
  
  typedef struct OM {
      float   x;
      float   y;
      int     carbone;
      int     lock;
  } OM;
  
  typedef struct MiorWorld {
      int     nbMM;
      int     nbOM;
      int     RA;
      float   RR;
      float   GR;
      float   K;
      int     width;
      int     minSize;
      int     CO2;
      int     lock;
  } MiorWorld;
  
  
  typedef struct RandomGen {
  	ulong a;
  	ulong b;
  	ulong c;
  } RandomGen;
  
  /**
   * HELPER FUNCTIONS
   */
  
  // Compatibility macros for OpenCL <= 1.1
  
  #if __OPENCL_VERSION__ <= CL_VERSION_1_0
      #define atomic_inc(p) atom_inc(p)
      #define atomic_add(p, v) atom_add(p, v)
      //#define atomic_cmpxchg(p, o, v) atom_cmpxchg(p, o, v)
  #endif
  
  // End of compatibility macros
  
  //#define LOCK_P(lock, id) while (atomic_cmpxchg(&(lock), -1, id) != id) {}
  //#define LOCK_V(lock, id) while (atomic_cmpxchg(&(lock), id, -1) != -1) {}
  #define store_CO2(world, val) atomic_add(&(world->CO2), val)
  #define val(arr, iMM, iOM) arr[iMM * world->nbOM + iOM]
  
  // RNG
  
  void rng_init(RandomGen *r, long seed)
  {
      r->a = seed;
      r->b = 0;
      r->c = 362436;
  }
  
  unsigned long rng_rand(RandomGen *r)
  {
      const long old = r->b;
      r->b = r->a * 1103515245 + 12345;
      r->a = (~old ^ (r->b >> 3)) - r->c++;
      return r->b;
  }
  
  float rng_rand_01(RandomGen *r)
  {
      return (rng_rand(r) & 4294967295) / 4294967295.0f;
  }
  
  void shuffle_indexes(global int * array, int size, RandomGen *r)
  {   
      for (int i = size - 1; i >= 1; i--) {
  		const int j = rng_rand_01(r) * (i + 1);
  		
  		const int t = array[j];
  		array[j] = array[i];
  		array[i] = t;
  	}
  }
  
  /**
   * KERNELS
   */
  
  kernel void topology(
      global MM *mmList,
      global OM *omList,
      global int *topo,
      global MiorWorld *world,
      global int *totals)
  {
      const int iMM = get_global_id(0);
      const int iOM = get_global_id(1);
      
  	const float dx = omList[iOM].x - mmList[iMM].x;
  	const float dy = omList[iOM].y - mmList[iMM].y;
  	const float dist = hypot(dx, dy);
  	
  	if (dist <= world->RA) {
  	    val(topo, iMM, iOM) = 0;
  	    atomic_inc(totals + iOM);
  	} else {
  	    val(topo, iMM, iOM) = -1;
  	}
  }
  
  kernel void carbon_scatter(
      global MM * mmList,
      global OM * omList,
      global int * topo,
      global MiorWorld * world,
      global int * totals)
  {
      const int iOM = get_global_id(0);
      const int total = totals[iOM];
      
      if (total == 0) {
          return;
      }
      
      const int part = (world->K * omList[iOM].carbone) / total; 
      int totalCarbon = 0;
      
      for (int iMM = 0; iMM < world->nbMM; iMM++) {
          const int offer = val(topo, iMM, iOM);
          
          if (offer >= 0) {
              val(topo, iMM, iOM) = part;
              totalCarbon += part;
          }
      }
      
      omList[iOM].carbone -= totalCarbon;
  }
  
  kernel void carbon_reduce(
      global MM * mmList,
      global OM * omList,
      global int * topo,
      global MiorWorld * world)
  {
      const int iOM = get_global_id(0);
      int totalCarbon = 0;
      
      for (int iMM = 0; iMM < world->nbMM; iMM++) {
          const int part = val(topo, iMM, iOM);
          
          if (part >= 0) {
              val(topo, iMM, iOM) = 0;
              totalCarbon += part;
          }
      }
      
      omList[iOM].carbone += totalCarbon;
  }
  
  kernel void live(
      global MM *mmList,
      global OM *omList,
      global int *topo,
      global MiorWorld *world,
      global int *omIndexes,
      long   seed)
  {   
  
      int iMM = get_global_id(0);
      global int * targets = topo + iMM * world->nbOM;
      global int * indexes = omIndexes + iMM * world->nbOM;
      RandomGen rng;
      
      for (int i = 0; i < world->nbOM; i++) {
          indexes[i] = i;
      }
      
      if (seed != 0) {
          rng_init(&rng, seed + iMM);
          shuffle_indexes(indexes, world->nbOM, &rng);
      }
      
      
      const int breathNeed = world->RR * mmList[iMM].carbone;
      int remainingBreathNeed = breathNeed;
      int iOM;
      
      // Check breath requirement
      for (int i = 0; i < world->nbOM; i++) {
          iOM = indexes[i];
          
          if (val(topo, iMM, iOM) >= 0) {
              remainingBreathNeed -= val(topo, iMM, iOM);
              
              if (remainingBreathNeed <= 0) {
                  break;
              }
          }
      }
      
      // Go into dormancy if breath requirement is not met
      if (remainingBreathNeed > 0) {
          mmList[iMM].dormance = 1;
          return;
      } else {
          mmList[iMM].dormance = 0;
      }
      
      // Actual breathing and growth
      remainingBreathNeed = breathNeed;
      const int growthNeed  = world->GR * mmList[iMM].carbone + breathNeed;
      int remainingGrowthNeed = growthNeed;
      int consum = 0, growth = 0;
      
      for (int i = 0; i < world->nbOM; i++) {
          iOM = indexes[i];
          
          int offer = val(topo, iMM, iOM);
          
          if (offer >= 0) {
              // If breathing phase is not finished
              if (remainingBreathNeed > 0) {
                  consum = min(offer, remainingBreathNeed);
                  atomic_add(&(world->CO2), consum);
                  remainingBreathNeed -= consum;
                  offer -= consum;
              }
              
              consum = min(offer, remainingGrowthNeed);
              remainingGrowthNeed -= consum;
              growth += consum;
              offer -= consum;
              
              val(topo, iMM, iOM) = offer;
          }
          
          if (remainingGrowthNeed <= 0) {
              break;
          }
      }
      
      mmList[iMM].carbone += growth;
  }