#include "main.h"


// TODO: check if no bias in selfing: can indiv pick itself?
// TODO: while we spawn idividuals: new inds have been added to the population
//       Can they be used as parents? MAJOR ISSSUE !
//       => Not as long as index & demogrphy are not updated?


/* * * * * * * * * * * REMINDER * * * * * * * * * * * * * * * * * * * * * * * *
 *
 * The program is supposed to go as the following:
 * - Load lua script
 * - Load species, trait and individuals
 *
 * Then starts the main loop,
 *    Iterate on cells,
 *        Iterate on species (pops in cell)
 *            Apply events to individuals accordig to the scheduler
 *
 *
 * Clock object: iterates tics
 * Calculates events at each tick (or at the beginning of the program?)
 * Set a priority list in the script to determine the order of events occurring a the same tick
 *
 * If nothing happens, no need to update each individual (age), fitness
 * Age should be updatable since the last event (use birth date and not age?=
 * Should I divide populations in age classes?
 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */


// Default seeds
double S1 = 10;
double S2 = 11;

// ENVT PARAMETERS
float STORM_PROB = 0.1;


//bool verbose = true;
bool checkRandom = false;
bool saveTime = true;
bool saveByLocus = false;
bool saveGeneflow = false;
bool geneOut = true;

bool verbose = false;

indId individual::counter_created(0);
indId individual::counter_alive(0);


void randomValues()
{
    // check random randomness
    cout << "r: " ;
    for (int i=0;i<5;i++)
    {
        //cout << random() << " ; ";
    }
    cout << endl;
    cout << "g: " ;
    for (int i=0;i<5;i++)
    {
        cout << gasdev() << " ; ";
    }
    cout << endl;

    cout << "R: " ;
    for (int i=0;i<5;i++)
    {
        cout << Random() << " ; ";
    }
    cout << endl;
}

int main()
{

    ofstream eventFile;
    ofstream mortalityFile;

    vector<int> secondsByTick;
    time_t tInit, t0, t1;
    tInit = 0;
    t0 = 0;
    t1 = 0;

    if (saveTime)
    {
        tInit = time(NULL);
    }

    // Lua main script name declaration
    string file = "luaSrc/main.lua";

    // New lua state (for now the only state?)
    cout << "\n\nStrating Lua scripts.."  << endl;
    lua_wrapper state(file);
    cout << "End of Lua initialization, Starting C++\n" << endl;


    //ProfilerStart("/tmp/profile");

    cout << "Loading parameters from config.lua" << endl;
    // Open config.lua table
    S1 = 0;
    S2 = 0;
    state.open_global_table("configSeed"); // <- ERROR: Open global table, getglobal got wrong type or non existant value for seed
    //Get seeds for random generator
    S1 = state.get_field<int>("seed1");
    S2 = state.get_field<int>("seed2");
    state.close_table();

    //cout << "GET SEEDS" << S1 << " " << S2 << endl;

    state.open_global_table("config");
    //Eref (when plasticity only)
    //int Eref = state.get_field<int>("Eref");
    //cout << "ER" << Eref << endl;
    // Init random function
    Random(S1);

    if (checkRandom)
    {
        randomValues();
    }


    // Get map size parameters
    int mw = state.get_field<int>("mapLength");
    int ml = state.get_field<int>("mapWidth");


    // Clock parameters
    double simLength = state.get_field<double>("simLength");
    double timestep = state.get_field<double>("timeStep");

    double savestep = state.get_field<double>("saveStep");
    bool isRegression = state.get_field<bool>("isRegression");
    unsigned int cellRef = state.get_field<int>("cellRef");

    cout << "CELL REF FOR GENEFLUX OUTPUTS:" << cellRef << endl;


    // Output management
    string outputDir = state.get_field<string>("outputDir");
    cout << "OUT:" << outputDir << endl;


    eventFile.open(outputDir+"eventsRecap.txt"); // flush the event file
    eventFile << "Step Type Name CutValue (Cell)" << endl;
    eventFile.close();
    mortalityFile.open(outputDir+"mortalityRecap.txt"); // flush the event file
    mortalityFile << "Step Name Cell CutValue Impacted" << endl;
    mortalityFile.close();

    // Close table
    state.close_table();
    cout << "End of loading." << endl;


    /*********************** Map Creation (with an lenght) ********************/
    cout << "Generating Map"  << endl;
    landscape myMap(mw, ml);


    /*********************** Loading events:***********************************/
    cout << "Loading Events"  << endl;
    vector<Event> allEvents;
    load_events(state, "eventsTable", allEvents);
    /*Event eInsect("insect", false, "T1", 0.1, -3.0, 0.7);
    Event eLatefrost("late frost", true, "T1", 0.1, -3.5, 0.8);
    allEvents.push_back(eInsect);
    allEvents.push_back(eLatefrost);*/


    /************************  Loading species: *******************************/
    cout << "Loading Species"  << endl;
    vector<species> speciesVect;
    load_species(state, "speciesTable", speciesVect);

    // visualisation of chromosom map ?

    species * currSpe = NULL;
    for (unsigned int i = 0; i < speciesVect.size(); ++i)
    {
            currSpe = &(speciesVect[i]);
            currSpe->checkChromoStruct(true);

            //initiate empty pops
            myMap.addIsEmptyForSpe(currSpe);

            currSpe->update_SelInt_matrix(myMap);
    }

    /****** Loading events: ******/

    vector<Event> globalEvents;
    vector<Event> localEvents;

    for (int i=0;i<allEvents.size();i++)
    {
        if (allEvents[i].get_isGlobal())
        {
            globalEvents.push_back(allEvents[i]);
        }
        else
        {
            localEvents.push_back(allEvents[i]);
        }
    }


    /********************* Loading intitial populations ***********************/
    cout << "Loading Individuals"  << endl;
    //TODO OVERLOAD LOAD_POPULATION FOR DIFFERENT INITIAL AGE
    load_populations(state, "populations", speciesVect, myMap);

    for (unsigned int i = 0; i < speciesVect.size(); ++i)
    {
            currSpe = &(speciesVect[i]);
            myMap.refreshIsEmptyForSpe(currSpe);
    }

    /************************* Init clock *************************************/
    cout << "Init Clock"  << endl;
    mClock myClock(timestep, simLength);

    cout << "End of initialization\n\n" << endl;


    cell * currCell = NULL;
    population * currPop  = NULL;

    valarray<double> sfit; //temporary fitnesse storage
    vector<trait *> vTraits;
    //individual_index indInd;

    //cout << "Checkpoint Alpha"  << endl;
    //Output Useful infos that are not changed through simulation
    output_ploidy(speciesVect, outputDir);
    //cout << "Checkpoint Beta"  << endl;


    if (saveTime)
    {
        t1 = time(NULL);
        secondsByTick.push_back(difftime(t1, tInit));
    }
    cout << "INITIALIZATION TIME:" << difftime(t1, tInit) << endl;
    cout << "Main loop" << endl;

    /**************************************************************************/
    /****************************    MAIN LOOP    *****************************/
    /**************************************************************************/
    while (myClock.inc_step())
    {
        if (saveTime)
        {
            t0 = t1;
        }
        cout << "\n\n### Clock " << myClock.get_time() << " ###\n" << endl;
        cout << "Alive ind: " << individual::counter_alive << endl;
        cout << "Created ind: " << individual::counter_created << endl;

        /*********************   Genetic outputs    ***************************/
         if (( myClock.get_time()==1)
            || (!isRegression && myClock.get_time()%(int)savestep == 0)
            || (isRegression && (myClock.get_time()==250 || myClock.get_time()==1000))
           )
        {
            if (geneOut)
            {
              cout << "Genet out" << endl;
              for (unsigned int i = 0; i < speciesVect.size(); i++)
              {
                  currSpe = &(speciesVect[i]);
                  output_genetic(myMap, currSpe, myClock, outputDir);
                  if (saveGeneflow)
                  {
                    output_geneflow(myMap, currSpe, myClock, outputDir, cellRef);
                  }
              }
            }
        }

        /**************** Update ages / traits / fitness **********************/
        cout << "Update Ages/trait/fitness" << endl;
        for (unsigned int i = 0; i < speciesVect.size(); ++i)
        {
            currSpe = &(speciesVect[i]);
            currSpe->get_age_classes(myMap);
        }

        for (unsigned int i = 0; i < speciesVect.size(); ++i)
        {
            currSpe = &(speciesVect[i]);
            myMap.refreshIsEmptyForSpe(currSpe);
            for(unsigned int cellId = 0; cellId < myMap.max_size(); ++cellId)
            {
                currCell = myMap.get_cell(cellId);
                currPop = currCell->get_pop(currSpe);
                currPop->make_ages_index(myClock);
                currPop->calc_traits();
                currPop->refresh_handicap();

            }
        }
        for (unsigned int i = 0; i < speciesVect.size(); ++i)
        {
            currSpe = &(speciesVect[i]);
            myMap.refresh_pop_ageClass(currSpe);

            vector<int> getC = myMap.get_pop_ageClass(currSpe);//USELESS?
        }


        /*****************  Making group index     ****************************/
        cout << "Make groupe index" << endl;
        for (unsigned int i = 0; i < speciesVect.size(); ++i)
        {
            currSpe = &(speciesVect[i]);

            if (currSpe->is_assort_strategy("unmatchable"))
            {
                // clean unmatchable vector
                currSpe->clear_unmatchable();
            }

            currSpe->calc_demo(myMap);

            vTraits = currSpe->get_traits();
            for(trait * currTrait : vTraits)
            {
               currSpe->make_gIndex(myMap, currTrait->get_name() + "-pheno", 1); // SEGFAULT WHEN POP IS EMPTY
            }
            currSpe->make_gIndex(myMap, "fitness", 1);
            currSpe->make_gIndex(myMap, "age", 1); // WHy?
            //currSpe->make_gIndex(myMap, "age", 1); // WHy? DUPLICATE?

            //Update assortMating threshold
            if(currSpe->get_assortMating() != "")
            {
                for(trait * currTrait : vTraits)
                {
                   currSpe->make_gIndex(myMap, currTrait->get_name() + "-phenoWithE", 1); // SEGFAULT WHEN POP IS EMPTY
                }
                currSpe->update_assort_threshold(myMap);
            }
        }

        bool isOutputTime = (( myClock.get_time()==1)
                || (!isRegression && myClock.get_time()%(int)savestep == 0)
                || (isRegression && (myClock.get_time()==250 || myClock.get_time()==1000))
	        );
        bool isEventOutputTime = true; //isOutputTime;

        /*************************   Outputs values    ************************/
        if (isOutputTime)
        {
            cout << "Outputs" << endl;
            for (species& s : speciesVect)
            {
                currSpe = &s;
                //Print age and fitness
                output_value(myMap, currSpe, myClock, "age", outputDir);
                output_value(myMap, currSpe, myClock, "fitness", outputDir);


                vTraits = currSpe->get_traits();
                for(trait* t : vTraits)
                {
                   output_value(myMap,  currSpe, myClock, t->get_name() + "-add", outputDir);
                   output_value(myMap,  currSpe, myClock, t->get_name() + "-pheno", outputDir);
                   output_value(myMap,  currSpe, myClock, t->get_name() + "-dz", outputDir);
                   if(currSpe->get_assortMating() != "")
                   {
                       output_value(myMap,  currSpe, myClock, t->get_name() + "-phenoWithE", outputDir);
                   }
                   output_value(myMap,  currSpe, myClock, t->get_name() + "-epsi", outputDir);

		               if (saveByLocus)
                   {
                        output_addByLocus(myMap,  currSpe, myClock, t, outputDir);
                   }

                   // check if several ltypes, generates separated ltypes add output files if more than on
                   if (t->get_nLtype() > 1)
                   {
                        output_addByLtype(myMap,  currSpe, myClock, t, outputDir);
                   }
                }
            }
        }

        /*******************  Refresh environment *****************************/
        int actualT = myClock.get_time();

        eventFile.open(outputDir+"eventsRecap.txt", std::ios_base::app);
        // OCCURENCE OF ENVIRONMENTAL EVENTS
      	for(unsigned int cellId = 0; cellId < myMap.max_size(); ++cellId)
        {
            currCell = myMap.get_cell(cellId);
            currCell->reset_events();
      	}
	cout << "Events" << endl;
        // 1) Global events: affect the whole map
        for (Event& e : globalEvents)
        {
            Event* currEvent = &e;
            float cutValue = currEvent->pick_cutValue();


            EventOccur eOcc(e, cutValue);
            //cout << "add event " << eOcc.get_name() << " to all cells with cutValue " << cutValue << endl;
            if (isEventOutputTime)
            {
              eventFile << myClock.get_time() << " GE " << eOcc.get_name() << " " << cutValue << " All" << endl;
            }

            for(unsigned int cellId = 0; cellId < myMap.max_size(); ++cellId)
            {
                currCell = myMap.get_cell(cellId);
                currCell->add_event(eOcc);
            }
        }

	//cout << "Add local events: ";
        for(unsigned int cellId = 0; cellId < myMap.max_size(); ++cellId)
        {
            currCell = myMap.get_cell(cellId);
            //currCell->reset_events();
            // 2) Local events: specific on patch (insect invasion)
            for (Event& e : localEvents)
            {
                Event* currEvent = &e;
                float cutValue = currEvent->pick_cutValue();
                EventOccur eOcc(e, cutValue);
                currCell->add_event(eOcc);
                //cout <<  eOcc.get_name() << " [ " << cellId << ": " << cutValue << "] ; ";
                if (isEventOutputTime)
                {
                  eventFile << myClock.get_time() << " LE " << eOcc.get_name() << " " << cutValue << " " << cellId << endl;
                }
            }
        }
    	//cout << endl;

        // STORM CASIS // TODO FIND LAST STEP BEFORE LAST CUT?
        // memory for storm events impact
        // stormTouched[speciesId][patchId] = nbKilledByStorm
        bool stormTouched = false;
        if (Random() < STORM_PROB)
        {
            stormTouched = true;
            cout << "Storm !" << endl;
            if (isEventOutputTime)
            {
                eventFile << myClock.get_time() << " GE Storm " << 0 << " All" << endl;
            }
        }

        eventFile.close();

        // RAW CHANGES OF SPECIES PARAMETERS
        for (species& s : speciesVect)
        {
            currSpe = &s;
            // REFRESH FLUXMAT
            currSpe->adapt_fluxMat(actualT);

            vTraits = currSpe->get_traits();

            bool isChangeSelint = false;
            // REFRESH E AND ZOPT
            for(trait* t : vTraits)
            {
                if (t->isKeyInTemporalE(actualT))
                {
                    if (actualT > 1)
                    {
                        cout << "CHANGE ENVIRONMENT FOR TRAIT " << t-> get_name() << endl;
                    }
                    t->set_E(actualT);
                }

                if (t->isKeyInTemporalZopt(actualT))
                {
                    if (actualT > 1)
                    {
                        cout << "CHANGE ZOPT FOR TRAIT " << t-> get_name() << endl;
                    }
                    t->set_Zopt(actualT);
                }

                if (t->isKeyInTemporalSelInt(actualT))
                {
                    if (actualT > 1)
                    {
                        cout << "CHANGE SELINT FOR TRAIT " << t-> get_name() << endl;
                    }
                    t->set_SelInt(actualT);
                    isChangeSelint = true;
                }
            }
            if (isChangeSelint)
            {
                currSpe->update_SelInt_matrix(myMap);
            }
        }


        mortalityFile.open(outputDir+"mortalityRecap.txt", std::ios_base::app);

        for(unsigned int cellId = 0; cellId < myMap.max_size(); ++cellId)
        {

            //Iteration on species
            for (unsigned int i = 0; i < speciesVect.size(); ++i)
            {
                currCell = myMap.get_cell(cellId);
                currSpe = &(speciesVect[i]);
                currPop = currCell->get_pop(currSpe);


                vector<EventOccur> eOccs = currCell->get_events();
                for (EventOccur& eOcc : eOccs)
                {
                    EventOccur* currEventOcc = &eOcc;
                    if (verbose) cout << "EVENT " << currEventOcc->get_name() << " in cell " << cellId << endl;
                    // TODO INJECT INTO OUTPUT
                    if (currPop->get_age() <= currSpe->getAgeEnds()[0]-1 && isEventOutputTime) // last step as young
                    {
                        int nselected = currPop->abs_select(currEventOcc->get_relatedTrait(), currEventOcc->get_cutValue(), currEventOcc->get_lethality());
                        mortalityFile << myClock.get_time() << " " << eOcc.get_name() << " " << cellId << " " << currEventOcc->get_cutValue() <<  " " << nselected << endl;
                    }
                    currPop->impact_height(currEventOcc->get_relatedTrait(), currEventOcc->get_cutValue(), currEventOcc->get_handicap());
                }

		            int kbSize = currPop->purge_killBuffer();
                if (currPop->get_age() <= currSpe->getAgeEnds()[0]-1) // last step as young
                {
                    mortalityFile << myClock.get_time() << " " << "All" << " " << cellId << " " << "NaN" <<  " " << kbSize << endl;
                }

            		// REFRESH HEIGHT
            		currPop->calc_height();
	          }
        }

        if (stormTouched)
        {
            for(unsigned int cellId = 0; cellId < myMap.max_size(); ++cellId)
            {
                currCell = myMap.get_cell(cellId);
                //Iteration on species
                for (unsigned int i = 0; i < speciesVect.size(); ++i)
                {
                    currSpe = &(speciesVect[i]);
                    vector<int> currAges = currSpe->getAgeEnds();
                    vector<int> currSizes = currSpe->getAgeSizes();
                    currPop = currCell->get_pop(currSpe);
                    int impacted = currPop->storm_impact();
                    if (impacted > 0)
                    {
                      mortalityFile << myClock.get_time() << " " << "Storm" << " " << cellId << " " << currPop->get_age() << " " << impacted << endl;
                    }
                }
            }
        }


        /*******************  Demography & Repro  *****************************/
        cout << "Demography/cut" << endl;
        //vector<vector<population *>> newSpePops;
        vector<vector<int>> speNewPatches(speciesVect.size());
        vector<vector<int>> speCutValue(speciesVect.size());
        //TODO HANDLE DIFFERENT AGE CLASSES - EXTERNALIZE IN SPE?
        //vector<vector<int>> curPops(speciesVect.size(), vector<int>(myMap.max_size()));

        for(unsigned int cellId = 0; cellId < myMap.max_size(); ++cellId)
        {
            currCell = myMap.get_cell(cellId);

            //Iteration on species
            for (unsigned int i = 0; i < speciesVect.size(); ++i)
            {
                //vector<population *> newPops(myMap.max_size());

                currSpe = &(speciesVect[i]);
                vector<int> currAges = currSpe->getAgeEnds();
                vector<int> currSizes = currSpe->getAgeSizes();
                currPop = currCell->get_pop(currSpe);

                currPop->update_ages(myClock); // SHOULD BE DONE FOR ALL BEFORE ANY ELDERING?
            }
        }

        for(unsigned int cellId = 0; cellId < myMap.max_size(); ++cellId)
        {
            currCell = myMap.get_cell(cellId);

            //Iteration on species
            for (unsigned int i = 0; i < speciesVect.size(); ++i)
            {
                currSpe = &(speciesVect[i]);
                vector<int> currAges = currSpe->getAgeEnds();
                vector<int> currSizes = currSpe->getAgeSizes();
                currPop = currCell->get_pop(currSpe);

                currPop->update_ages(myClock); // SHOULD BE DONE FOR ALL BEFORE ANY ELDERING? DOUBLON?

                /* CHECK IF CHANGE OF AGE CLASS */
                int agePatch = currPop->get_age()+1;//get_rInd(currPop->pick_random()).get_age();
                if (agePatch != currPop->get_rInd(currPop->pick_random()).get_age())
                {
                  cout << "Compare ages:" << agePatch << " V S " << currPop->get_rInd(currPop->pick_random()).get_age() << endl;
                }
                //int id = currPop->get_rInd(currPop->pick_random()).get_id();
                //cout << "age?" << cellId << " " << id << " " << agePatch << endl;
                //int maxAge = currAges[currAges.size()-1];

                bool is_new_pop = false; // will become true only if end of last age class
                int cutValue = 0;
                //bool is_select_action = false; // will become true if end of class != last

                int c=0;
                bool found = false;
                while (c<currAges.size() && !found)
                {
                    if (agePatch == currAges[c])
                    {
                        found = true;
                        if (c == currAges.size()-1)
                        {
                            cout << "new pop in cell " << cellId << endl;
                            is_new_pop = true;
                            speNewPatches[i].push_back(cellId);
                        }
                        else
                        {
                            //cout << "is transition:" << agePatch << endl;
                            //is_cut = true;
                            cutValue = currSizes[c+1];
                            //cout << "cut" << cellId << ";" << cutValue << endl;
                        }
                    }
                    c++;
                }
                speCutValue[i].push_back(cutValue);
                if (is_new_pop)
                {
                    //currPop->kill_range("age", maxAge, maxAge+1); // PUT KILL HERE ?
                    int newP = currPop->demography_new(); //Calculate demography

                    if (newP == 0)
                    {
                        myMap.setIsEmpty(currSpe, cellId, true);
                    }
                    else
                    {
                        myMap.setIsEmpty(currSpe, cellId, false);
                    }
                    currPop->adjust(myMap, myClock);
                }

                //cout << cellId << "cut ? " << cutValue;
                if (cutValue > 0)
                {
                    //1) if young: clean impacted by event
                    /*if (currPop->get_age() == currSpe->getAgeEnds()[0]-1) // last step as young
                    {
                       vector<EventOccur> eOccs = currCell->get_events();
                       for (EventOccur& eOcc : eOccs)
                       {
                            EventOccur* currEventOcc = &eOcc;
                            cout << "EVENT " << currEventOcc->get_name() << " in cell " << cellId << endl;
                            currPop->abs_select(currEventOcc->get_relatedTrait(), currEventOcc->get_cutValue(), currEventOcc->get_lethality());
                       }

                    }*/
                    //2) Adjust remaning pop
                    currPop->force_targetPop(cutValue);
                    currPop->adjust_generationmode(myMap, myClock); //Resolves demography
                }
            }
        }

        if (isOutputTime)
        {
            for (unsigned int i = 0; i < speciesVect.size(); ++i)
            {
                currSpe = &(speciesVect[i]);
                output_value(myMap,  currSpe, myClock, "height", outputDir);
            }
        }

        // need refresh of demo ?
        for (unsigned int i = 0; i < speciesVect.size(); ++i)
        {
            currSpe = &(speciesVect[i]);
            currSpe->calc_demo(myMap);
        }

        /*********************** Kill above threshold *************************/
        cout << "Kill" << endl;
        for(unsigned int cellId = 0; cellId < myMap.max_size(); ++cellId)
        {
            currCell = myMap.get_cell(cellId);

            //Iteration on species
            for (unsigned int i = 0; i < speciesVect.size(); ++i)
            {
                currSpe = &(speciesVect[i]);
                int maxAge = currSpe->get_max_age();
               //cout << currSpe << "/" << cellId << endl;
                currPop = currCell->get_pop(currSpe);
                currPop->kill_range("age", maxAge-1, maxAge);


                currPop->purge_killBuffer();

                // INJECT NEW POPULATION
                if (speNewPatches[i].size() > 0)
                {
                    for (int p=0;p<speNewPatches[i].size();p++)
                    {
                        if (speNewPatches[i][p] == cellId)
                        {
                            //cout << "Will change pop" << endl;
                            currPop->set_new_population();
                            //cout << "NEW POP SIZE:" << currPop->get_nInd() << endl;
                        }
                    }
                }
            }
        }
        mortalityFile.close();

        if (saveTime)
        {
            t1 = time(NULL);
            secondsByTick.push_back(difftime(t1, t0) );
        }
    }

    //End time loop
    cout << "End of simulation." << endl;
    if (saveTime)
    {
        output_time(secondsByTick, outputDir);
    }

    //ProfilerStop();
    // DELETE ALL OBJECTS
    // pops
    for(unsigned int cellId = 0; cellId < myMap.max_size(); ++cellId)
    {
        currCell = myMap.get_cell(cellId);

        //Iteration on species
        for (unsigned int i = 0; i < speciesVect.size(); ++i)
        {
            currSpe = &(speciesVect[i]);
            currPop = currCell->get_pop(currSpe);
            delete currPop;
        }
    }

    return 0;
}
