mior_model.cl 7.01 KB
/**
 * 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;


/**
 * CONSTANTS
 */

#define OFFER_INVALID -1
#define OFFER_UNINIT  -2

/**
 * HELPER FUNCTIONS
 */

// Compatibility macros for OpenCL <= 1.1

#if __OPENCL_VERSION__ <= CL_VERSION_1_0
    #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) {}

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

void shuffle_help(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;
	}
}


kernel void shuffle_test(global int * array, int size)
{
	RandomGen rng;
	rng_init(&rng, 42L);
	
	for (int i = 0; i < size; i++) {
		array[i] = i;
	}
	
	shuffle_help(array, size, &rng);
}
	
// Cleanup

void topo_clean(global int * targets, int size)
{
    for (int i = 0; i < size; i++) {
        if (targets[i] != OFFER_INVALID) {
            targets[i] = OFFER_UNINIT;
        }
    }
}

/**
 * KERNELS
 */

kernel void topology(
    global MM *mmList,
    global OM *omList,
    global int *topo,
    global MiorWorld *world)
{
    const int iMM = get_global_id(0);
    const int iOM = get_global_id(1);
    
	//float dx = omList[iOM].x - mmList[iMM].x;
	//float dy = omList[iOM].y - mmList[iMM].y;
	//float distance = sqrt(dx * dx + dy * dy);
	
    //if (distance <= world->RA) {
	//    topo[iMM * world->nbOM + iOM] = OFFER_UNINIT;
	//} else {
	//    topo[iMM * world->nbOM + iOM] = OFFER_INVALID;
	//}
	
	const float dx = omList[iOM].x - mmList[iMM].x;
	const float dy = omList[iOM].y - mmList[iMM].y;
	const float distance = hypot(dx, dy);
	
	if (distance <= world->RA) {
	    topo[iMM * world->nbOM + iOM] = OFFER_UNINIT;
	} else {
	    topo[iMM * world->nbOM + iOM] = OFFER_INVALID;
	}
}

int live_breath(
    int iMM,
    global MM * mmList,
    global OM * omList,
    global int * targets,
    global MiorWorld * world,
    global int *omIndexes)
{
    global MM * currentMM = mmList + iMM;
    global OM * currentOM;
    
    const int totalNeed = currentMM->carbone * world->RR;
    int remainingNeed = totalNeed;
    const int nbOM = world->nbOM;
    int i = 0, iOM = 0;
    int offer, consum;
    
    while (i < nbOM && remainingNeed > 0) {
        iOM = omIndexes[i];
        
   		if (targets[iOM] != OFFER_INVALID) {
    		currentOM = omList + iOM;
            
            offer = world->K * currentOM->carbone;
            consum = min(offer, remainingNeed);
            remainingNeed -= consum;
    	}
    	
    	i++;
    }
    
    if (remainingNeed > 0) {
    	currentMM->dormance = 1;
    } else {
    	currentMM->dormance = 0;
    	remainingNeed = totalNeed;
  
		i = 0;
		while (i < nbOM && remainingNeed > 0) {
		    iOM = omIndexes[i];
		    
		    if (targets[iOM] != OFFER_INVALID) {
		        currentOM = omList + iOM;
		         
		        LOCK_P(currentOM->lock, iMM);
		        
		        if (targets[iOM] != OFFER_UNINIT) {
		            offer = targets[iOM];
		        } else {
		            offer = world->K * currentOM->carbone;
		        }
		        
		        consum = min(offer, remainingNeed);
		        
		        remainingNeed -= consum;
		        currentOM->carbone -= consum;
		        
		        LOCK_V(currentOM->lock, iMM);
		        
		        LOCK_P(world->lock, iMM);
		        world->CO2 += consum;
		        LOCK_V(world->lock, iMM);
		        
		        targets[iOM] = offer - consum;
		    }
		    
		    i++;
		}
	}
    
    return remainingNeed;
}

int live_grow(
    int iMM,
    global MM * mmList,
    global OM * omList,
    global int * targets,
    global MiorWorld * world,
    global int *omIndexes)
{
    global MM * currentMM = mmList + iMM;
    global OM * currentOM;
    
    const int totalNeed = currentMM->carbone * world->GR;
    int remainingNeed = totalNeed;
    const int nbOM = world->nbOM;
    int i = 0, iOM = 0;
    int offer, consum;
    
    while (i < nbOM && remainingNeed > 0) {
        iOM = omIndexes[i];
        
        if (targets[iOM] != OFFER_INVALID) {
            currentOM = omList + iOM;
            
            LOCK_P(currentOM->lock, iMM);
            
            if (targets[iOM] != OFFER_UNINIT) {
                offer = targets[iOM];
            } else {
                offer = world->K * currentOM->carbone;
            }
            
            consum = min(offer, remainingNeed);
            
            remainingNeed -= consum;
            currentOM->carbone -= consum;
            
            LOCK_V(currentOM->lock, iMM);
            
            currentMM->carbone += consum;
        }
        
        i++;
    }
    
    return remainingNeed;
}

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);
    }
    
    // BREATHING
    live_breath(iMM, mmList, omList, targets, world, indexes);
    
    if (mmList[iMM].dormance == 0) {
    
    	// GROWTH
    	live_grow(iMM, mmList, omList, targets, world, indexes);
    }
    
    // CLEANUP
    topo_clean(targets, world->nbOM);
}

kernel void autolive(
    global MM *mmList,
    global OM *omList,
    global int *topo,
    global MiorWorld *world,
    global int *omIndexes,
    long   seed)
{
    LOCK_P(world->lock, get_global_id(0));
    int CO2 = world->CO2;
    LOCK_V(world->lock, get_global_id(0));
    
    int oldCO2 = -1;
    
    while (CO2 != oldCO2) {
        live(mmList, omList, topo, world, omIndexes, seed);
        
        oldCO2 = CO2;
        LOCK_P(world->lock, get_global_id(0));
        CO2 = world->CO2;
        LOCK_V(world->lock, get_global_id(0));
    }   
}