/**
 * 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);
    }
}