simulation_failed_multihive.cu

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
#include <iostream>
#include <thrust/transform.h>
#include <thrust/fill.h>
#include <math.h> 
#include <thrust/shuffle.h>
#include <thrust/random.h>
#include <curand.h>
#include <curand_kernel.h>
#include <math.h>

#define PI 3.14159265
#define beeMax 10000000
#define threadsPerBlock 256

using namespace std;
using namespace thrust;

struct plot
{
    unsigned long long total_flowers;
    unsigned long long pollinated_flowers = 0;
    unsigned long long pollinated_last_year = 0;
    unsigned long long pollinated_today = 0;
    unsigned long long flowers_available = 0;
};

struct hive
{
    int total_nutrition;
    int should_bee_die = 0;
    int should_brood_die = 0;
    int total_bees = 0;
    int total_nurse_bees = 0;
    int total_worker_bees = 0;
    int total_broods = 0;
    int unattained_broods = 0;
    int foraged_increment = 0;
    int nutrition_deficit = 0;
    int eggLaid = 0;
    bool shouldSwarm = false;
    int swarm_bees;
    int swam_nutrition;
    int nutrition_last_year = 0;
};

struct bee
{
    hive* hive_ptr;
    int age;
    float nutrition; //decay by 0.01 per day, dies after negative -1
    int type; //0 = worker, 1 = drone

    __device__ __host__ bee() {}

    __device__ __host__ bee(hive *h, int a, float n, int t)
    {
        hive_ptr = h;
        age = a;
        nutrition = n;
        type = t;
    }
};

__host__ __device__ bool operator<(const bee &a, const bee &b) { return (a.age < b.age); };

struct brood
{
    int age;
    hive* hive_ptr;

    __device__ __host__ brood() {}

    __device__ __host__ brood(hive *h, int a) 
    {
        hive_ptr = h;
        age = a;
    }
};

__host__ __device__ bool operator<(const brood &a, const brood &b) { return (a.age < b.age); };

int bee_cnt = 0;
device_vector<plot> d_plots(18*18);//18 km x 18 km
device_vector<bee> d_bees(beeMax);
device_vector<brood> d_broods(beeMax);

int hive_cnt = 1;
device_vector<hive> d_hives(100);
host_vector<hive> h_hives(100);

//statics
__device__ int broodConsumption;
__device__ int beeConsumption;
double brood_mature_death_rate = 0.02;
__device__ double egg_laying_rate_multiplier = 1;
__device__ int egg_laying_rate_base = 100;
__device__ double starvation_death_winter_multiplier = 1.5;
__device__ double starvation_death_winter_divisor = 10;
__device__ double starvation_deficit_power = 0.5;
__device__ double starvation_deficit_multiplier = 0;
//
int bc = 2;
int pc = 6;

__device__ float total = 0;

__device__ double germination_rate = 0.9;

__global__ void plotPollinationUpdate(plot *p, int plot_cnt, hive *h, double efficiency, int day)
{
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    if(index < plot_cnt)
    {
        if(day >= 50 && day <= 100)
        {
            unsigned long long blossomedFlowers = (p[index].pollinated_last_year*35) * germination_rate;
            p[index].total_flowers += blossomedFlowers/(101-day);
        }

        if(p[index].pollinated_today > p[index].flowers_available)
        {
            //this will be replaced with something that addresses the number of pollination per plot by any bee at a given time soon
            //most likely another set of arrays per plot where it stores the amount of pollination done by any bee given i from a hive.
            h[index].total_nutrition -= (int)efficiency * (p[index].pollinated_flowers - p[index].total_flowers)/26;
            //I see
            p[index].pollinated_today = p[index].flowers_available;
        }

        p[index].pollinated_flowers += p[index].pollinated_today;

        if(day >= 320 && day <= 364)
        {
            p[index].total_flowers -= p[index].total_flowers/(365 - day); // p[index].backup_total * (1.0-((day-240.0)/30.0));
            p[index].pollinated_last_year = p[index].pollinated_flowers;
        }
        if(day == 364)
        {
            p[index].pollinated_flowers = 0;
            p[index].total_flowers = 0;
        }
        
        p[index].flowers_available = p[index].total_flowers - p[index].pollinated_flowers;
        if(p[index].total_flowers < p[index].pollinated_flowers) p[index].flowers_available = 0;
    }
}

__global__ void broodUpdate(brood *b, int brood_cnt)
{
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    if(index < brood_cnt)
    {
        atomicAdd(&(*(b[index].hive_ptr)).total_broods, 1);//hmm consumes 1 only, wtf man
        b[index].age += 1;
        if(b[index].age >= 3 && b[index].age <= 8)
        {
            //atomicAdd(&(*(b[index].hive_ptr)).total_nutrition, -6);
            atomicAdd(&(*(b[index].hive_ptr)).total_nutrition, -broodConsumption);
        }
    }
}

__device__ double eggLayingRate(double day)
{
    double rate = 0.0000000001 * pow(day, 5) - 0.0000000823  * pow(day, 4) + 0.0000143494  * pow(day, 3) + 0.0001276996  * pow(day, 2) - 0.0670571925 * day + 0.5824552258;
    if(rate < 0 || day < 50) rate = 0;
    //return rate * 30 + 100;
    return (rate * 200 + 100) * egg_laying_rate_multiplier;
    //return (rate * 200 + 100)*egg_laying_rate_multiplier; //with a multiplier of 200, you will have 2000-70000 population
    //return (cos((day/365 - 0.5) * PI * 2 ) + 1) * 1000 + 500;
}

__device__ double eggLayingRate2(double day)
{
    double rate = - 0.0000000174 * pow(day, 5) + 0.0000160050 * pow(day, 4) - 0.0049974512 * pow(day, 3) + 0.5490013102 * pow(day, 2) - 5.2106195939 * day + 54.5018666934;
    if(rate <= 0)
    {
        rate = 100;
    }
    return rate;
}

__global__ void hiveUpdate(hive *h, int hive_cnt, int day, plot *p, curandState *randstate)
{
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    if(index < hive_cnt)
    {
        int remainingDays = 0;

        if(day < 150)
        {
            remainingDays = 150 - day;
        }
        else if(day > 280)
        {
            remainingDays = 150 + (365 - day);
        }

        int totalNutritionRequirements = (int)((h[index].total_bees*beeConsumption + h[index].total_broods*broodConsumption));//brood 6 death
        h[index].nutrition_deficit = h[index].foraged_increment - totalNutritionRequirements;;

        int insufficientNutrition = 0;

        if(remainingDays > 0) 
        {
            //h[index].should_die += (0.01 * h[index].total_bees);
            insufficientNutrition += max((int)(totalNutritionRequirements * starvation_death_winter_multiplier) - h[index].total_nutrition/remainingDays, 0)/(remainingDays/starvation_death_winter_divisor);

            //h[index].should_die += max(totalNutritionRequirements - h[index].total_nutrition/remainingDays, 0)/remainingDays;
        }
        else
        {
            insufficientNutrition += (int)pow((-min(h[index].nutrition_deficit, 0)), starvation_deficit_power) * starvation_deficit_multiplier;
        }

        //h[index].eggLaid = eggLayingRate2(day);
        h[index].eggLaid = ((p[0].flowers_available / 20000000000000)+1) * eggLayingRate(day);// / hive_cnt;
        //h[index].eggLaid = eggLayingRate(day);

        h[index].should_bee_die = insufficientNutrition/(2*beeConsumption);
        h[index].should_brood_die = insufficientNutrition/(2*broodConsumption);

        h[index].nutrition_deficit += insufficientNutrition;

        if(h[index].total_nutrition < 0)
        {
            h[index].should_bee_die = h[index].total_bees;
            h[index].should_brood_die = h[index].total_broods;
        }

        //h[index].shouldSwarm = true;
       // h[index].shouldSwarm = true;

        if(day == 190)
        {
            h[index].shouldSwarm = true;
            if(h[index].eggLaid > 1500)
            {
                if(abs(h[index].nutrition_last_year - h[index].total_nutrition) < 2500)
                {
                    //stabilized enough for swarming
                    //if random that is...
                    h[index].shouldSwarm = true;
                    if(curand_uniform(&(randstate[index])) <= 0.4)
                    {
                        //swarming
                        h[index].shouldSwarm = true;
                    }
                }
            }
            h[index].nutrition_last_year = h[index].total_nutrition;
        }

        h[index].unattained_broods = 0;
        //h[index].unattained_broods = max(h[index].total_broods - h[index].total_nurse_bees*4, 0);
        h[index].total_broods = 0;
        h[index].total_bees = 0;
        h[index].total_nurse_bees = 0;
        h[index].total_worker_bees = 0;
        h[index].foraged_increment = 0;
    }
}

__global__ void beeGlobalUpdate(bee *b, int bee_cnt, int lifespan)
{
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    if(index < bee_cnt)
    {
        b[index].age += 1;
        //(*(b[blockIdx.x].hive)).total_nutrition -= 0.01;
        //if((*(b[index].hive_ptr)).total_nutrition < 10)
        atomicAdd(&(*(b[index].hive_ptr)).total_nutrition, -beeConsumption);
        atomicAdd(&(*(b[index].hive_ptr)).total_bees, 1);

        if(b[index].age >= 0 && b[index].age <= lifespan*0.7)
        {
            atomicAdd(&(*(b[index].hive_ptr)).total_nurse_bees, 1);
        }
    }
};

__global__ void beeForageUpdate(bee *b, int bee_cnt, plot *p, double efficiency, double lifespan)
{
    int index = threadIdx.x + blockIdx.x * blockDim.x;
    if(index < bee_cnt)
    {
        if(b[index].type == 0 && b[index].age >= 0.5*lifespan)
        {
            //p[0].pollinated_flowers += 1;
            atomicAdd(&(*p).pollinated_today, 111*efficiency/8.57087);
            atomicAdd(&(*(b[index].hive_ptr)).total_nutrition, (int)efficiency);
            atomicAdd(&(*(b[index].hive_ptr)).total_worker_bees, 1);
            atomicAdd(&(*(b[index].hive_ptr)).foraged_increment, (int)efficiency);
        }
    }
}

struct bee_death_functor : thrust::unary_function<bee, bool>
{
    const float lifespan;

    bee_death_functor(float _a) : lifespan(_a) {}

    __host__ __device__
    bool operator()(const bee &b)
    {
        return b.age <= lifespan;
    }
};

struct bee_hive_filter_functor : thrust::unary_function<bee, bool>
{
    const hive* hive_ptr;

    bee_hive_filter_functor(hive* _a) : hive_ptr(_a) {}

    __host__ __device__
    bool operator()(const bee &b)
    {
        return b.hive_ptr == hive_ptr;
    }
};

struct bee_hive_not_filter_functor : thrust::unary_function<bee, bool>
{
    const hive* hive_ptr;

    bee_hive_not_filter_functor(hive* _a) : hive_ptr(_a) {}

    __host__ __device__
    bool operator()(const bee &b)
    {
        return b.hive_ptr != hive_ptr;
    }
};


struct brood_hive_filter_functor : thrust::unary_function<brood, bool>
{
    const hive* hive_ptr;

    brood_hive_filter_functor(hive* _a) : hive_ptr(_a) {}

    __host__ __device__
    bool operator()(const brood &b)
    {
        return b.hive_ptr == hive_ptr;
    }
};

struct brood_hive_not_filter_functor : thrust::unary_function<brood, bool>
{
    const hive* hive_ptr;

    brood_hive_not_filter_functor(hive* _a) : hive_ptr(_a) {}

    __host__ __device__
    bool operator()(const brood &b)
    {
        return b.hive_ptr != hive_ptr;
    }
};

struct bee_hive_change_functor
{
    hive* hive_ptr;

    bee_hive_change_functor(hive* _a) : hive_ptr(_a) {}

    __host__ __device__
        bee operator()(bee& b) const 
        {
            b.hive_ptr = hive_ptr;
            return b;
        }
};

struct brood_mature_functor : thrust::unary_function<brood, bool>
{
    hive* hive_ptr;

    brood_mature_functor(hive* _a) : hive_ptr(_a) {}

    __host__ __device__
    bool operator()(const brood &x)
    {
        return x.age < 21 && x.hive_ptr != hive_ptr;
    }
};

int getPlotNum(int& x, int& y)
{
    return y*18+x;
}

void addBee(bee& add_type, const int& count)
{
    thrust::fill(d_bees.begin()+bee_cnt, d_bees.begin()+bee_cnt+count, add_type);
    bee_cnt += count;
}

void addRatioBee(const bee& b, const int& count)
{
    if(count <= 0) return;
    int drone = count/101;
    thrust::fill(d_bees.begin()+bee_cnt, d_bees.begin()+bee_cnt+drone, bee(b.hive_ptr, b.age, b.nutrition, 1));
    thrust::fill(d_bees.begin()+bee_cnt+drone, d_bees.begin()+bee_cnt+count, bee(b.hive_ptr, b.age, b.nutrition, 0));
    bee_cnt += count;
}

int simulationDays = 50;
bool verbose = true;

void outValue(string name, int val)
{
    if(verbose) cout << name << ": " << val << " | ";
    else cout << val << "|";
}

void outValue(string name, float val)
{
    if(verbose) cout << name << ": " << val << " | ";
    else cout << val << "|";
}

void outValue(string name, unsigned long long val)
{
    if(verbose) cout << name << ": " << val << " | ";
    else cout << val << "|";
}

void outValue(string name, double val)
{
    if(verbose) cout << name << ": " << val << " | ";
    else cout << val << "|";
}


double forageEfficiency(double day)
{
    if(day < 60 || day > 340) return 0;
    double efficiency = 0.0000000247 * pow(day, 4) - 0.0000193509 * pow(day, 3) + 0.0044706630 * pow(day, 2) - 0.2590425605 * day + 1.6599060149;
    if(efficiency < 0) efficiency = 0;
    return efficiency * 0.64;
}

int lifeSpan(double day)
{
    return (int)(0.0067211910 * pow(day, 2) - 2.5117894242 * day + 252.9914557348);
}

int main(int argc, char** argv)
{

    hive h;
    h.total_nutrition = 5000000;

    for(int i = 1; i < argc; i++)
    {
        if(!strcmp(argv[i], "-s"))
        {
            verbose = false;
        }

        if(!strcmp(argv[i], "-d"))
        {
            if(i+1 < argc)
            {
                simulationDays = atoi(argv[i+1]);
            }
        }

        if(!strcmp(argv[i], "-bc"))
        {
            if(i+1 < argc)
            {
                bc = atoi(argv[i+1]);
            }
        }
        cudaMemcpyToSymbol(beeConsumption, &bc, sizeof(int), 0, cudaMemcpyHostToDevice);

        if(!strcmp(argv[i], "-pc"))
        {
            if(i+1 < argc)
            {
                pc = atoi(argv[i+1]);
            }
        }
        cudaMemcpyToSymbol(broodConsumption, &pc, sizeof(int), 0, cudaMemcpyHostToDevice);

        if(!strcmp(argv[i], "-e"))
        {
            if(i+1 < argc)
            {
                double eg = atof(argv[i+1]);
                cudaMemcpyToSymbol(egg_laying_rate_multiplier, &eg, sizeof(double), 0, cudaMemcpyHostToDevice);
            }
        }

        if(!strcmp(argv[i], "-eb"))
        {
            if(i+1 < argc)
            {
                int eb = atoi(argv[i+1]);
                cudaMemcpyToSymbol(egg_laying_rate_base, &eb, sizeof(int), 0, cudaMemcpyHostToDevice);
            }
        }

        if(!strcmp(argv[i], "-sw"))
        {
            if(i+1 < argc)
            {
                double stvw = atof(argv[i+1]);
                cudaMemcpyToSymbol(starvation_death_winter_multiplier, &stvw, sizeof(double), 0, cudaMemcpyHostToDevice);
            }
        }

        if(!strcmp(argv[i], "-swd"))
        {
            if(i+1 < argc)
            {
                double stvw = atof(argv[i+1]);
                cudaMemcpyToSymbol(starvation_death_winter_divisor, &stvw, sizeof(double), 0, cudaMemcpyHostToDevice);
            }
        }

        if(!strcmp(argv[i], "-sdm")) { if(i+1 < argc) {
            double sdm = atof(argv[i+1]);
            cudaMemcpyToSymbol(starvation_deficit_multiplier, &sdm, sizeof(double), 0, cudaMemcpyHostToDevice);
        } }

        if(!strcmp(argv[i], "-sdp")) { if(i+1 < argc) {
            double sdp = atof(argv[i+1]);
            cudaMemcpyToSymbol(starvation_deficit_power, &sdp, sizeof(double), 0, cudaMemcpyHostToDevice);
        } }

        if(!strcmp(argv[i], "-in"))
        {
            if(i+1 < argc)
            {
                h.total_nutrition = atoi(argv[i+1]);
            }
        }
    }

    plot h_plot;
    //h_plot.total_flowers = 20000000000;
    h_plot.total_flowers = 200000000000;
    h_plot.pollinated_flowers = 0;

    thrust::fill(d_plots.begin(), d_plots.end(), h_plot);
    plot* d_plot_array = thrust::raw_pointer_cast(d_plots.data());
    
    thrust::fill(d_hives.begin(), d_hives.begin()+1, h);
    hive* d_hive_array = thrust::raw_pointer_cast(d_hives.data());

    host_vector<bee> h_bees(beeMax);

    //initialize bees
    //thrust::fill(d_bees.begin(), d_bees.begin()+bee_cnt, new_bee);

    for(int age = 0; age <= 30; age++)
    {
        addRatioBee(bee(&d_hive_array[0], age, 1, 0), h.total_nutrition/24000);
    }

    //50000 + 9900*14 - 11000*10
    //5000 + 9900*14 - 15000*10 = -6400
    //6400/10 = 640

    //brood new_brood = brood(&d_hive_array[0], 0);

    //initialize broods
    int brood_cnt = 1500;
    int matured_brood = 0;
    thrust::fill(d_broods.begin(), d_broods.begin() + brood_cnt, brood(&d_hive_array[0], 0));
    brood* d_brood_array = thrust::raw_pointer_cast(d_broods.data());

    bee* d_bee_array = thrust::raw_pointer_cast(d_bees.data());

    cudaDeviceSynchronize();

    curandState *d_state;
    cudaMalloc(&d_state, sizeof(curandState));

    device_vector<bee> hive_bees(beeMax);
    device_vector<brood> hive_broods(beeMax);

    for(int day = 0; day < simulationDays; day++)
    {
        outValue("Day", day);
        outValue("Brood", brood_cnt);
        outValue("Bee", bee_cnt);

        //thrust::sort(d_bees.begin(), d_bees.begin() + bee_cnt);
        thrust::default_random_engine g;
        thrust::shuffle(d_bees.begin(), d_bees.begin() + bee_cnt, g);

        /*h_bees = d_bees;
        outValue("FIRST", h_bees[0].age);
        outValue("LAST", h_bees[bee_cnt-1].age);*/
        double efficiency = forageEfficiency(day%365);
        outValue("EF", efficiency);
        double lspan = lifeSpan(day%365);
        beeForageUpdate<<<(bee_cnt + threadsPerBlock-1)/threadsPerBlock, threadsPerBlock>>>(d_bee_array, bee_cnt, d_plot_array, efficiency, lspan);
        cudaDeviceSynchronize();

        plotPollinationUpdate<<<(1 + threadsPerBlock-1)/threadsPerBlock, threadsPerBlock>>>(d_plot_array, 1, d_hive_array, efficiency, day%365);
        cudaDeviceSynchronize();

        cudaMemcpy(&h_plot, &d_plot_array[0], sizeof(plot), cudaMemcpyDeviceToHost);
        outValue("UF", h_plot.flowers_available);
        outValue("PF", h_plot.pollinated_flowers);
        outValue("TF", h_plot.total_flowers);

        beeGlobalUpdate<<<(bee_cnt + threadsPerBlock-1)/threadsPerBlock, threadsPerBlock>>>(d_bee_array, bee_cnt, lspan);
        cudaDeviceSynchronize();

        outValue("LSP", lspan);

        detail::normal_iterator<device_ptr<bee>> alive_bee_end = thrust::copy_if(d_bees.begin(), d_bees.begin() + bee_cnt, d_bees.begin(), bee_death_functor(lifeSpan(day%365)));
        bee_cnt = alive_bee_end - d_bees.begin();

        //thrust::transform(d_broods.begin(), d_broods.begin() + brood_cnt, d_broods.begin(), brood_functor());
        broodUpdate<<<(brood_cnt + threadsPerBlock-1)/threadsPerBlock, threadsPerBlock>>>(d_brood_array, brood_cnt);
        cudaDeviceSynchronize();

        hive o_hive;

        cudaMemcpy(&o_hive, &d_hive_array[0], sizeof(hive), cudaMemcpyDeviceToHost);

        outValue("Nurse", o_hive.total_nurse_bees);
        outValue("Worker", o_hive.total_worker_bees);

        hiveUpdate<<<(hive_cnt + threadsPerBlock-1)/threadsPerBlock, threadsPerBlock>>>(d_hive_array, hive_cnt, day%365, d_plot_array, d_state);
        cudaDeviceSynchronize();

        h_hives = d_hives;

        bool swarmModified = false;

        for(int i = 0; i < hive_cnt; i++)
        {
            if(h_hives[i].shouldSwarm)
            {

                cout << "SWARMED" << endl;
                //change half of the bees
                //add it to the list
                swarmModified = true;
                detail::normal_iterator<device_ptr<bee>> hive_bees_end = thrust::copy_if(d_bees.begin(), d_bees.begin() + bee_cnt, hive_bees.begin(), bee_hive_filter_functor(&d_hive_array[i]));
                int hive_bees_cnt = hive_bees_end - hive_bees.begin();

                int newHiveIndex = hive_cnt;
                hive_cnt++;
                h_hives[newHiveIndex].total_bees = hive_bees_cnt/2;
                h_hives[newHiveIndex].total_nutrition = h_hives[i].total_nutrition * 4;
                h_hives[i].total_nutrition *= 4;//-= h_hives[newHiveIndex].total_nutrition;
                h_hives[i].total_bees -= h_hives[newHiveIndex].total_bees;

                thrust::sort(hive_bees.begin(), hive_bees_end);

                thrust::transform(hive_bees.begin(), hive_bees.begin() + h_hives[newHiveIndex].total_bees, hive_bees.begin(), bee_hive_change_functor(&d_hive_array[newHiveIndex]));

                detail::normal_iterator<device_ptr<bee>> master_bee_end = thrust::copy_if(d_bees.begin(), d_bees.begin() + bee_cnt, d_bees.begin(), bee_hive_not_filter_functor(&d_hive_array[i]));
                
                thrust::merge(hive_bees.begin(), hive_bees_end, d_bees.begin(), master_bee_end,  d_bees.begin());

                h_hives[i].shouldSwarm = false;
            }
        }

        if(swarmModified)
        {
            d_hives = h_hives;
            swarmModified = false;
        }

        for(int i = 0; i < hive_cnt; i++)
        {
            if(i == 0) outValue("HC", hive_cnt);

            cudaMemcpy(&o_hive, &d_hive_array[i], sizeof(hive), cudaMemcpyDeviceToHost);

            if(i == 0) outValue("Net", o_hive.nutrition_deficit);

            //int brood_die = (int)o_hive.should_die/2;

            //if(i == 0) outValue("SD", o_hive.should_bee_die);

            //starvation death
            //alive_bee_end -= min((int)o_hive.should_die, (int)(alive_bee_end - d_bees.begin()));


            //Fix it here
            detail::normal_iterator<device_ptr<bee>> hive_bees_end = thrust::copy_if(d_bees.begin(), d_bees.begin() + bee_cnt, hive_bees.begin(), bee_hive_filter_functor(&d_hive_array[i]));
            int originalCount = hive_bees_end - hive_bees.begin();
            hive_bees_end -= min((int)o_hive.should_bee_die, (int)(hive_bees_end - hive_bees.begin()));

            if(i == 0) outValue("D", (int)(originalCount - (hive_bees_end - hive_bees.begin())));

            detail::normal_iterator<device_ptr<bee>> master_bee_end = thrust::copy_if(d_bees.begin(), d_bees.begin() + bee_cnt, d_bees.begin(), bee_hive_not_filter_functor(&d_hive_array[i]));
            thrust::merge(hive_bees.begin(), hive_bees_end, d_bees.begin(), master_bee_end,  d_bees.begin());

            bee_cnt = (master_bee_end - d_bees.begin()) + (hive_bees_end - hive_bees.begin());


            //666
            //if(hive_cnt > 1) return;

            if(i == 0) outValue("S", o_hive.should_bee_die + o_hive.should_brood_die);

            detail::normal_iterator<device_ptr<brood>> hive_broods_end = thrust::copy_if(d_broods.begin(), d_broods.begin() + brood_cnt, hive_broods.begin(), brood_hive_filter_functor(&d_hive_array[i]));
            int originalBroodCnt = hive_broods_end - hive_broods.begin();

            hive_broods_end -= min(o_hive.should_brood_die + o_hive.unattained_broods, (int)(hive_broods_end - hive_broods.begin()));

            detail::normal_iterator<device_ptr<brood>> master_brood_end = thrust::copy_if(d_broods.begin(), d_broods.begin() + brood_cnt, d_broods.begin(), brood_hive_not_filter_functor(&d_hive_array[i]));
            
            thrust::merge(hive_broods.begin(), hive_broods_end, d_broods.begin(), master_brood_end,  d_broods.begin());

            brood_cnt = (master_brood_end - d_broods.begin()) + (hive_broods_end - hive_broods.begin());

            //brood_cnt -= min(o_hive.should_brood_die + o_hive.unattained_broods, brood_cnt);
            //brood_cnt -= min(o_hive.unattained_broods, brood_cnt);
        
            if(i == 0) outValue("U", o_hive.unattained_broods);

            detail::normal_iterator<device_ptr<brood>> unmatured_bee_end = thrust::copy_if(d_broods.begin(), d_broods.begin()+brood_cnt, d_broods.begin(), brood_mature_functor(&d_hive_array[i]));
            matured_brood =(d_broods.begin() + brood_cnt) - unmatured_bee_end;
            matured_brood = int(matured_brood * (1-brood_mature_death_rate));
            if(i == 0) outValue("M", matured_brood);

            //int obc = bee_cnt;
            
            // if(day > 673) return;
            addRatioBee(bee(&d_hive_array[i], 0, 1, 0), matured_brood);

            //if(i == 0) outValue("NC", bee_cnt - obc);

            brood_cnt = unmatured_bee_end - d_broods.begin();

            if(i == 0)  outValue("N", o_hive.total_nutrition);

            //lay eggs
            int eggLaid = o_hive.eggLaid;
            if(i == 0) outValue("Egg", eggLaid);

            thrust::fill(unmatured_bee_end, unmatured_bee_end+eggLaid, brood(&d_hive_array[i], 0));
            brood_cnt+=eggLaid;

            cudaDeviceSynchronize();
        }

        cout << endl;
    }

    cudaDeviceSynchronize();

    /*h_hives = d_hives;
    cout << o_hive.total_nutrition << endl;*/
}