#include <iostream>
#include <vector>
#include <ctime>
#include <sys/stat.h>
#include "netcdf.h"
#include "dictionary.h"
#include "iniparser.h"

#include "config.h"
#include "grid_info.h"
#include "gridded_input.h"
#include "zone.h"

#include "controller.h"

using namespace std;
using namespace glm;

void Controller::sanity_checks() {
   if (RES == 0.25) {
      if (NCCODE != 199) {
         cerr << " Incorrect number of countries " << NCCODE << " for given resolution " << RES << "\n";
         exit(0);
      }
   }
   else if (RES == 0.5) {
      if (NCCODE != 192) {
         cerr << " Incorrect number of countries " << NCCODE << " for given resolution " << RES << "\n";
         exit(0);
      }     
   }
}

/******************************************************************************/
void Controller::initialize (char* name, const char* confFileName) {

   cout << "\n\nInitializing GLM\n\n";
   sanity_checks();
      
   now = 0;
   next = 1;

   config.load(name, confFileName);
   cout << "loaded config\n";

   if (mkdir(config.outputDir_, 0755) != 0 && errno != EEXIST) {
      cerr << "Failed to create output directory\n";
      exit(1);
   }

   countries.resize(NCCODE);

   grid = new GridInfo (countries, config);

   optionSettings(*grid);

   cout << "opening inputs\n";
   inputs = new InputFiles();
   inputs->init(grid, config);

   // Load CFT fraction data for first year, so that we can read in CFT fractions
   // in disaggregateInitialLandUse
   inputs->openCFTFractionFile(config.cft_frac_file);
   if (! config.extensionRun) {
   inputs->loadCFTFractionData(config.start_year, countries);
   }
   else {
   inputs->loadCFTFractionData(2015, countries);
   }
   // count valid grid cells
   // while we are at it, initialize zdis
   size_t count = 0;
   for (size_t y=0; y<NY; y++) {
      for (size_t x=0; x<NX; x++) {
        if (grid->ice_->get(y, x) < 1 | grid->countryCode[y][x] > 0) {
            count++;
         }
         zdis[y][x] = NZ - 1;
         Ldis[y][x] = NL - 1;
         
      }
   }
   // make the cells vector the right size to hold all the cells we need
   cells.resize(count);

   // Initialize all the newly created cells in the vector
   size_t i = 0;
   for (size_t y=0; y<NY; y++) {
      for (size_t x=0; x<NX; x++) {
        if (grid->ice_->get(y, x) < 1 | grid->countryCode[y][x] > 0) {
            cells[i].init(now, y, x, &countries[grid->countryIndex[y][x]],
                          *grid, *inputs, config);
            i++;
	  
         }
      }
   }
 
   for (CellIterator c=cells.begin(); c!= cells.end(); ++c) {
     c->frac[now].RSForestLossRemaining = (grid->rs_forest_loss_2012[c->y][c->x]);
     c->frac[now].landsat = grid->rs_forest_loss_2012[c->y][c->x];
   }
   
   
   string outFile = config.outputDir_;
   outFile += "/";
   outputter.init(outFile, config.start_year, grid->lat, grid->lon, config);

}

/******************************************************************************/
void Controller::step (int year, GridInfo& grid) {

   clock_t smartFlowTime = 0;
   clock_t adjSmartFlowTime = 0; 

   cout << "starting step: curr_year=" << year 
        << " start_year=" << config.start_year
        << " stop_year=" << config.stop_year << "\n";

   if (year >= config.start_year && year <= config.stop_year) {

for (ZoneIterator z=countries.begin(); z!=countries.end(); ++z) {
     if (year <1850){
	      z->priority = VRGN_PRIORITY;
          z->force_priority = VRGN_PRIORITY;
	      z->includeClearingInHarvest = 1;
     }
     else if (year>=1850 & year < 1920){
          z->includeClearingInHarvest = 1- (double(year) - 1850) / (1920 - 1850);
          if (z->includeClearingInHarvest > 1){
	          z->includeClearingInHarvest = 1;
          }
          if (z->includeClearingInHarvest < 0){
	          z->includeClearingInHarvest = 0;
          }
          if (config.BEST_CASE) {
             if ( (z->continentCode == 3) || (z->continentCode == 4) ) {
	            z->priority = VRGN_PRIORITY;
             } else {
	            z->priority = VRGN_PRIORITY;
             }
             if (config.smart_flow_option == 1) {
	            z->force_priority = VRGN_PRIORITY;
             } else {
	           z->force_priority = VRGN_PRIORITY;
             }
         } else {
             if (config.smart_flow_option == 1) {
	           z->priority = VRGN_PRIORITY;
	           z->force_priority = VRGN_PRIORITY;
             } else {
	           z->priority = VRGN_PRIORITY;
               z->force_priority = VRGN_PRIORITY;
             }
        }
     }
     else{
       if (config.BEST_CASE) {
         if ( (z->continentCode == 3) || (z->continentCode == 4) ) {
	        z->priority = VRGN_PRIORITY;
         } else {
	        z->priority = VRGN_PRIORITY;
         }
         if (config.smart_flow_option == 1) {
	        z->force_priority = VRGN_PRIORITY;
         } else {
	        z->force_priority = VRGN_PRIORITY;
         }
       } else {
          if (config.smart_flow_option == 1) {
	        z->priority = VRGN_PRIORITY;
	        z->force_priority = VRGN_PRIORITY;
         } else {
	        z->priority = VRGN_PRIORITY;
            z->force_priority = VRGN_PRIORITY;
         }
      }
      z->includeClearingInHarvest = 0;
         }
 }

      updateZDistance();
      updateLDistance(year);
      combineZLDistances(year);
      
      for (ZoneIterator z=countries.begin(); z!=countries.end(); ++z) {
         z->resetTotals();
      }
      inputs->loadNextStates (year);
      double global_crop_now = 0.0;
      double global_crop_next = 0.0;
      
      for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
         global_crop_now += c->frac[now].lu[CROP]*c->cellArea;
         global_crop_next += c->frac[next].lu[CROP]*c->cellArea;
      }
      
      for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
         c->transitions (now, next, *inputs, config,
                         smartFlowTime, adjSmartFlowTime, year, grid);
      }

      // Read in CFT fraction data
      if ( (! config.futureRun) & (! config.extensionRun)) {
         if (year==2016){
            inputs->loadCFTFractionData(2015, countries);
         }
         else{
            inputs->loadCFTFractionData(year, countries);
         }
      }
      else {
            inputs->loadCFTFractionData(2015, countries);
      }
         
      if (config.TOTAL_HARVEST_SWITCH) {
#if 1
	inputs->loadWoodHarvestData (year - 1, countries);
      
	europe.unmet = 0;
	europe.availBiomass = 0;

for (int nr = 0; nr<NREG; nr++){
    futureReg[nr].unmet = 0;
    futureReg[nr].availBiomass = 0;
}

	for (ZoneIterator z=countries.begin(); z!=countries.end(); ++z) {
	      	   z->totalAvailBiomass = 0;
	   for (int znum=0; znum<NZ; znum++) {
	  
	     z->totalAvailBiomass += z->vrgnForestBiomass[znum];
	     z->totalAvailBiomass += z->scndZdisForestBiomass[znum];
	   }
	  z->totalAvailBiomass += z->vrgnNonForestBiomass + z->scndNonForestBiomass;
	  z-> harvestRemaining = z->harvestDemand;     // units = MgC
	  z->calcConvertedForest();
	  
	  if ((z->continentCode == 3 | z->countryCode == 901)&(!config.futureRun)&(!config.extensionRun)) {
	    if (z->harvestRemaining  <= z->totalAvailBiomass){
	      z->totalAvailBiomass -= z->harvestRemaining;
	      europe.availBiomass += z->totalAvailBiomass;
	    }
	    else {
	      europe.unmet += (z->harvestRemaining - z->totalAvailBiomass);
	      z->harvestRemaining = z->totalAvailBiomass;
	      z->totalAvailBiomass = 0;
	    }
	  }
	  
	  if ((config.futureRun==1) | (config.extensionRun==1)) {
	    if (z->harvestRemaining  <= z->totalAvailBiomass){
	      z->totalAvailBiomass -= z->harvestRemaining;
	      futureReg[z->regionCode-1].availBiomass += z->totalAvailBiomass;
	    }
	    else {
	      futureReg[z->regionCode-1].unmet += (z->harvestRemaining - z->totalAvailBiomass);
	      z->harvestRemaining = z->totalAvailBiomass;
	      z->totalAvailBiomass = 0;
	    }
	  }
	  
	}

	 for (ZoneIterator z=countries.begin(); z!=countries.end(); ++z) {
	   z->newHarvestDemand = 0;
	   if ((z->continentCode == 3 | z->countryCode == 901)&(!config.futureRun)&(!config.extensionRun)){
	     if (europe.unmet < europe.availBiomass) {
	       z->harvestRemaining += z->totalAvailBiomass / europe.availBiomass * europe.unmet;
	     }
	     else {
	       if (europe.availBiomass > 0){
		 z->harvestRemaining += z->totalAvailBiomass;
	 
	       }
	     }
	   }
	   
	   if ((config.futureRun ==1) | (config.extensionRun==1)){
	     if (futureReg[z->regionCode-1].unmet < futureReg[z->regionCode-1].availBiomass) {
	       z->harvestRemaining += z->totalAvailBiomass / futureReg[z->regionCode-1].availBiomass * futureReg[z->regionCode-1].unmet;
	     }
	     else {
	       if (futureReg[z->regionCode-1].availBiomass > 0){
		      z->harvestRemaining += z->totalAvailBiomass;
	       }
	     }
	   }
	   
	   z->precomputeHarvestRatios (config);
	   
	   z->totalHarvested = 0;
	   z->newHarvestDemand = z->harvestRemaining;
	 }
	 
#endif
	 
     if (config.run_management == 1){
	    if ( (! config.futureRun) & (!config.extensionRun)) {
	       if (year==2016){
	          inputs->loadTillageData(2015, countries); 
	          inputs->loadFertData (2015, countries, config);
	          inputs->loadFuelWoodData (2015, countries);
              inputs->loadCropBiofuelData (2015, countries);
              inputs->loadCropRotationsData(2015, countries);
	       }
	       else{
	       
	       inputs->loadTillageData(year, countries); 
	       
	          inputs->loadFertData (year, countries, config);
	          inputs->loadFuelWoodData (year, countries);
              inputs->loadCropBiofuelData (year, countries);
              inputs->loadCropRotationsData(year, countries);
           }
	    }
	    else {
	    if (config.extensionRun){
	       inputs->loadTillageData(2015, countries);
	    }
	    else {
           inputs->loadTillageData(2015, countries);
           }
           if (year == (config.start_year+1)){
              inputs->loadFertData (config.start_year, countries, config);
           }
	       inputs->loadFertData (year, countries, config);
	       if (year == (config.start_year+1)){
              inputs->loadIrrigData (config.start_year, countries, config);
           }
	       inputs->loadIrrigData (year, countries, config);
	       inputs->loadFuelWoodData (year, countries);
	       inputs->loadNationalFloodedData (year, countries);
           inputs->loadCropBiofuelData (2015, countries);
           if (config.extensionRun){
           inputs->loadCropRotationsData(2015, countries);
           }
           else {
           inputs->loadCropRotationsData(2015, countries);
           }
	    }
    }

         for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
#if 1
	   c->harvest (now, *inputs, config, year);
	   
	   }
	   
#endif
	    
            c->doSubTransitions(now, next, config);
            if (config.futureRun) {
               if (config.c4p_flag==1){
               c->expandC4Perennials(now, next);
               }
               if (config.c3p_flag==1 ){
               c->expandC3Perennials(now, next);
               }
               if (config.c4a_flag==1 ){
               c->expandC4Annuals(now, next);
               }
               if (config.c3n_flag == 1 ) {
               c->expandC3NFixing(now, next);
               }
            }
            if (config.crop_rotations) {
               c->cropRotations(now, next);
            }
            
            c->updateStates (now, next, config);
            if (config.futureRun) {
               c->updateCropTypesRatios (next, config);
            }
            c->updateCropTypeTallies (next);
            c->updateFertTallies (now);
            c->updateIrrigTallies (now, next, config);
	    
            c->updateSecondaryLands(now, next);
} // end cell iterator
 double temp;
 
// harmonize fertilizer data
if (config.futureRun | config.extensionRun){
 for (ZoneIterator z=countries.begin(); z!=countries.end(); ++z) {

 temp = z->irrigData;
 z->irrigData = ((z->irrigData - z->irrigDataPrev) + z->irrigCurrent);
 z->irrigDataPrev = temp;

z->irrigChange = z->irrigData - z->irrigCurrent;
z->irrigDefaultChange = z->irrigCurrentNew - z->irrigCurrent;
z->irrigNeededChange = z->irrigChange - z->irrigDefaultChange;

if (z->irrigNeededChange <= 0){
   if (z->irrigCurrent>0.000001){
      z->irrigDecFrac = z->irrigNeededChange / z->irrigCurrentNew;
      if (z->irrigDecFrac < -1) {
         z->irrigDecFrac = -1;
      }
   }
   else {
      z->irrigDecFrac = 0;
   }
}
else {
   if (z->irrigAvail>0.000001){
      if (z->irrigNeededChange <= z->irrigAvail){
         z->irrigIncFrac1 = z->irrigNeededChange/ z->irrigAvail;
      }
      else {
         z->irrigIncFrac1 = 1;
         z->irrigUnmet = z->irrigNeededChange - z->irrigAvail;
         
         if (z->irrigNew >0.000001){
            z->irrigIncFrac2 = z->irrigUnmet / z->irrigNew;
            if (z->irrigIncFrac2 > 1){
               z->irrigIncFrac2 = 1;
            }
         }
         else {
            z->irrigIncFrac2 = 0;
         }    
      }
   }
   else {
      z->irrigIncFrac1 = 0;
   }
   
}
 
 }
 }
 
for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
            if (config.run_management == 1) {
               if (!config.futureRun & !config.extensionRun){
               c->computeIrrigatedArea (next);
               }
               else if (config.extensionRun){
               c->computeIrrigatedAreaExtension (next);
               }
               else {
               c->computeIrrigatedAreaFuture (now, next);
               }
               if (!config.futureRun & !config.extensionRun){
               c->computeFertilizerAmount (next);
               }
               else if (config.extensionRun){
               c->computeFertilizerAmountExtension (next);
               }
               else {
               c->computeFertilizerAmountFuture (now,next,config);
               }
               if (config.futureRun==1){
                c->computeNationalFloodedFraction(now,next);
               }
               if (config.futureRun==0 & config.extensionRun==0){
                  c->computeFloodedFraction (next);
               }
               else if (config.extensionRun==1){
               c->computeFloodedFractionExtension (next);
               }
               c->computeTillageRate (next);
               c->computeWHBiofuelFractions (next);
               c->computeCroplandBiofuelFractions (next);
               c->computeFutureCroplandBiofuelFractions (next);
               c->computeBiomassHarvestFraction (next);
            }
 
            outputter.updateOutputs(*c, now);
	    
         }
    
      }
      
      outputter.flush();

      next = next ? 0 : 1;
      now = now ? 0 : 1;

   }
}


/******************************************************************************/
void Controller::runModel () {
   for (int year=config.start_year+1; year<=config.stop_year; year++) {
     step (year, *grid);
   }
}


/******************************************************************************/
void Controller::finalize () {
   // write out the final states and mngmnt
   for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
      outputter.updateOutputs(*c, now);
   }
   outputter.flush(false);
}


/******************************************************************************/
void Controller::optionSettings (GridInfo& grid) {


   for (ZoneIterator z=countries.begin(); z!=countries.end(); ++z) {
      
      if (config.BEST_CASE) {
         if ( (z->continentCode == 3) || (z->continentCode == 4) ) {
            z->priority = SCND_PRIORITY;
         } else {
            z->priority = VRGN_PRIORITY;
         }
         if (config.smart_flow_option == 1) {
            z->force_priority = VRGN_PRIORITY;
         } else {
            z->force_priority = SCND_PRIORITY;
         }
      } else {
         if (config.smart_flow_option == 1) {
            z->priority = VRGN_PRIORITY;
            z->force_priority = VRGN_PRIORITY;
         } else {
            z->priority = SCND_PRIORITY;
            z->force_priority = SCND_PRIORITY;
         }
      }

      z->useZDistance = true;
      
      z->includeClearingInHarvest = 0;

   }
}

/***************************************************************/
void Controller::updateLDistance (int year) {
   // Ldis represents distance from a Landsat cell
   // The graphic below helps understand how zdis is calculated for each cell
   // Here y,x represent cell location i.e. lat,lon
   // For the focal cell, zdis is 0, all others have been initialized to NL
   // earlier in the code
   //+-------+-------+-------+
   //|       |       |       |
   //| x - 1 |   x   | x + 1 |
   //| y + 1 | y + 1 | y + 1 |
   //|       |       |       |
   //|       |       |       |
   //+-----------------------+
   //|       |       |       |
   //| x - 1 |   x   | x + 1 |
   //|   y   |   y   |   y   |
   //|       | Focal |       |
   //|       | Cell  |       |
   //|       |       |       |
   //+-----------------------+
   //|       |       |       |
   //| x - 1 |   x   | x + 1 |
   //| y - 1 | y - 1 | y - 1 |
   //|       |       |       |
   //|       |       |       |
   //|       |       |       |
   //+-------+-------+-------+

for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
         Ldis[c->y][c->x] = NL - 1;
   }
   
   // determine focal cell for zdis, initialize all others to NZ
   for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
      if (c->hasLandsat(now)) {
      if (year>=1950 & year<=2050 ){
         Ldis[c->y][c->x] = 1;
         }
         else {
         Ldis[c->y][c->x] = 0;
         }
         
      }
   }

int min_L = 0;
if (year>=1950 & year<=2050 ){
min_L = 1;
}
   // check all other cells relative to focal cell
   for (int L=min_L; L<NL; L++){
      for (int y=0; y<NY; y++) {
         for (int x=0; x<NX; x++) {
            if (Ldis[y][x] == L) {

               // Check for y and x to be greater than 0 so that we do not compute
               // zdis for cells outside of our boundary
               if ((y > 0) && (x > 0) && (Ldis[y-1][x-1] > L + 1)) {
                  Ldis[y-1][x-1] = L + 1;
               }

               if ((y > 0) && (Ldis[y-1][x] > L + 1)) {
                  Ldis[y-1][x] = L + 1;
               }

               if ((y > 0) &&  (x + 1 < NX) && (Ldis[y-1][x+1] > L + 1)) {
                  Ldis[y-1][x+1] = L + 1;
               }

               if ((x > 0) && (Ldis[y][x-1] > L + 1)) {
                  Ldis[y][x-1] = L + 1;
               }

               if ((x + 1 < NX) && (Ldis[y][x+1] > L + 1)) {
                  Ldis[y][x+1] = L + 1;
               }

               if ((y + 1 < NY) && (x > 0) && (Ldis[y+1][x-1] > L + 1)) {
                  Ldis[y+1][x-1] = L + 1;
               }

               if ((y + 1 < NY) && (Ldis[y+1][x] > L + 1)) {
                  Ldis[y+1][x] = L + 1;
               }

               if ((y + 1 < NY) && (x + 1 < NX) && (Ldis[y+1][x+1] > L + 1)) {
                  Ldis[y+1][x+1] = L + 1;
               }
            }
         } // end x loop
      } // end y loop
   } // end z loop


   // map zdis back to cells
   for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
    
    if ((c->frac[now].landsat * c->cellArea /1000/1000 > 125)&(year>=1950)&(year<=2050)){
         Ldis[c->y][c->x] = 0;
    }
         
      c->Ldis = Ldis[c->y][c->x];
   }
}
/******************************************************************************/
void Controller::updateZDistance () {
   // zdis represents distance from a focal cell
   // The graphic below helps understand how zdis is calculated for each cell
   // Here y,x represent cell location i.e. lat,lon
   // For the focal cell, zdis is 0, all others have been initialized to NZ
   // earlier in the code
   //+-------+-------+-------+
   //|       |       |       |
   //| x - 1 |   x   | x + 1 |
   //| y + 1 | y + 1 | y + 1 |
   //|       |       |       |
   //|       |       |       |
   //+-----------------------+
   //|       |       |       |
   //| x - 1 |   x   | x + 1 |
   //|   y   |   y   |   y   |
   //|       | Focal |       |
   //|       | Cell  |       |
   //|       |       |       |
   //+-----------------------+
   //|       |       |       |
   //| x - 1 |   x   | x + 1 |
   //| y - 1 | y - 1 | y - 1 |
   //|       |       |       |
   //|       |       |       |
   //|       |       |       |
   //+-------+-------+-------+

 for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
         zdis[c->y][c->x] = NZ - 1;
   }
   
   // determine focal cell for zdis, initialize all others to NZ
   for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
      if (c->hasLandUse(now)) {
         zdis[c->y][c->x] = 0;
         
      }
   }

   // check all other cells relative to focal cell
   for (int z=0; z<NZ; z++){
      for (int y=0; y<NY; y++) {
         for (int x=0; x<NX; x++) {
            if (zdis[y][x] == z) {

               // Check for y and x to be greater than 0 so that we do not compute
               // zdis for cells outside of our boundary
               if ((y > 0) && (x > 0) && (zdis[y-1][x-1] > z + 1)) {
                  zdis[y-1][x-1] = z + 1;
               }

               if ((y > 0) && (zdis[y-1][x] > z + 1)) {
                  zdis[y-1][x] = z + 1;
               }

               if ((y > 0) &&  (x + 1 < NX) && (zdis[y-1][x+1] > z + 1)) {
                  zdis[y-1][x+1] = z + 1;
               }

               if ((x > 0) && (zdis[y][x-1] > z + 1)) {
                  zdis[y][x-1] = z + 1;
               }

               if ((x + 1 < NX) && (zdis[y][x+1] > z + 1)) {
                  zdis[y][x+1] = z + 1;
               }

               if ((y + 1 < NY) && (x > 0) && (zdis[y+1][x-1] > z + 1)) {
                  zdis[y+1][x-1] = z + 1;
               }

               if ((y + 1 < NY) && (zdis[y+1][x] > z + 1)) {
                  zdis[y+1][x] = z + 1;
               }

               if ((y + 1 < NY) && (x + 1 < NX) && (zdis[y+1][x+1] > z + 1)) {
                  zdis[y+1][x+1] = z + 1;
               }
            }
         } // end x loop
      } // end y loop
   } // end z loop

   // map zdis back to cells
   for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
      c->zdis = zdis[c->y][c->x];
   }
}

/******************************************************************************/
void Controller::combineZLDistances (int year) {

for (CellIterator c=cells.begin(); c!=cells.end(); ++c) {
if (year<1901){
c->zdis = c->zdis;
}
else if (year >=1901 & year <2001){
if (year<1950){
c->zdis = floor((1-0.5*(year-1901)/100)*(c->zdis) + 0.5*(c->Ldis)*(year-1901)/100);
}
else{
c->zdis = floor((1-0.5*(year-1901)/100)*(c->zdis+1) + 0.5*(c->Ldis)*(year-1901)/100);
}
      
  }
  else if (year >=2001 & year<2016 ){
    c->zdis = floor(0.5*(c->zdis+1) + 0.5*(c->Ldis) );
  }
  else if (year >= 2016 & year<2100){
   c->zdis = floor( (1-0.5*(2100-year)/84)*(c->zdis+1) + (0.5*(2100-year)/84)*(c->Ldis) );
  }
  else if (year>=2100) {
      c->zdis = floor((c->zdis+1));
  }
  else {
   cout << "year out of range \n";
   }
}


}