Pixel.cpp

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2015 by Laboratoire d'Economie Forestière               *
 *   http://ffsm-project.org                                               *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 3 of the License, or     *
 *   (at your option) any later version, given the compliance with the     *
 *   exceptions listed in the file COPYING that is distribued together     *
 *   with this file.                                                       *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
#include "Pixel.h"
#include "ThreadManager.h"
#include "Scheduler.h"
#include "Init.h"
 
Pixel::Pixel(double ID_h, ThreadManager* MTHREAD_h): ID(ID_h)
{
  MTHREAD=MTHREAD_h;
  int nft = MTHREAD->MD->getForTypeIds().size();
  vector<double> temp(nft,1);
  //vector<double> temp2(nft,0);
  spMods = temp;
  avalCoef = temp;
  //vMort  = temp2;
  //std::fill(v.begin(), v.end(), 0);
}
 
Pixel::~Pixel()
{
}
 
/**
The function return a vector of pointers to Pixels at the gived distance from the caller pixel.\\
The list start with those on the Top, then add those on the right, those on the bottom and those on the left. Finally it had the corner pixels (that are more far).\\
It takes into consideration borders correctly.
 
Fully tested on internal points as well semi-border cases, border cases and corner cases. ALL OK.
 
@param distLevel_h Distance in number of adiacent pixels. It has to be at least 1 (the function return an error if it is 0).
*/
vector <Pixel *>
Pixel::getPixelsAtDistLevel(int distLevel_h) const{
 
  if (distLevel_h<1) {
    msgOut(MSG_CRITICAL_ERROR, "getPixelsAtDistLevel() is defined for distances of at least 1 !");
  }
 
  vector <Pixel *> toReturn;
  int xNPixels = MTHREAD->GIS->getXNPixels();
  int yNPixels = MTHREAD->GIS->getYNPixels();
  int thisX = this->getX();
  int thisY = this->getY();
  int minX = max(0        , (thisX - distLevel_h)+1);
  int maxX = min(xNPixels , thisX + distLevel_h);
  int minY = max(0        , (thisY - distLevel_h)+1);
  int maxY = min(yNPixels , thisY + distLevel_h);
 
  // getting the top pixels (corner exluded)...
  if (thisY-distLevel_h >=0){
    for(int i=minX;i<maxX;i++){
      toReturn.push_back(MTHREAD->GIS->getPixel(i,thisY-distLevel_h));
    }
  }
  // getting the right pixels (corner exluded)...
  if (thisX+distLevel_h < xNPixels){
    for(int i=minY;i<maxY;i++){
      toReturn.push_back(MTHREAD->GIS->getPixel(thisX+distLevel_h,i));
    }
  }
  // getting the bottom pixels (corner exluded)...
  if (thisY+distLevel_h < yNPixels){
    for(int i=minX;i<maxX;i++){
      toReturn.push_back(MTHREAD->GIS->getPixel(i,thisY+distLevel_h));
    }
  }
  // getting the left pixels (corner exluded)...
  if (thisX-distLevel_h >= 0){
    for(int i=minY;i<maxY;i++){
      toReturn.push_back(MTHREAD->GIS->getPixel(thisX-distLevel_h,i));
    }
  }
 
  // getting the corners (correctly at the end, after already retrieved the other pixels..)...
  // top-left..
  if (thisX-distLevel_h >= 0 && thisY-distLevel_h >=0){
    toReturn.push_back(MTHREAD->GIS->getPixel(thisX-distLevel_h,thisY-distLevel_h));
  }
  // top-right..
  if (thisX+distLevel_h < xNPixels && thisY-distLevel_h >=0){
    toReturn.push_back(MTHREAD->GIS->getPixel(thisX+distLevel_h,thisY-distLevel_h));
  }
  // bottom-right..
  if (thisX+distLevel_h < xNPixels && thisY+distLevel_h <yNPixels){ // bug discovered 20070719
    toReturn.push_back(MTHREAD->GIS->getPixel(thisX+distLevel_h,thisY+distLevel_h));
  }
  // bottom-left..
  if (thisX-distLevel_h >= 0 && thisY+distLevel_h <yNPixels){
    toReturn.push_back(MTHREAD->GIS->getPixel(thisX-distLevel_h,thisY+distLevel_h));
  }
  return toReturn;
} 
 
/*void //moved as inline function
Pixel::setValue(const string & layerName_h, const double & value_h){
 
  //tempValuePair.first = layerName_h;     // type of first is string
  //tempValuePair.second = value_h;            // type of second is double
  //tempValuePair = make_pair (layerName_h,value_h);
  values.insert(pair<string, double>(layerName_h, value_h));
  //values.insert(tempValuePair);
 
 
}*/
 
/*
inline void
Pixel::setValue (const string& parName, const string& forName, const string& dClass, const int& year, const double& value_h){
  values.insert(pair<string, double>(MTHREAD->GIS->pack(parName, forName, dClass, year), value_h));
}
*/
 
 
void
Pixel::changeValue(const string &layerName_h, const double &value_h, const bool &setNoValueForZero){
  map<string, double>::iterator p;
  p=values.find(layerName_h);
  if(p != values.end()){
    if(setNoValueForZero && value_h == 0){
      p->second = MTHREAD->GIS->getNoValue();
    } else {
      p->second = value_h;
    }
  } else {
    msgOut(MSG_ERROR, "Coud not change pixel value for layer "+layerName_h+". Layer don't found.");
  }
  return;
}
 
/*
void
Pixel::changeValue (const double &value_h, const string& parName, const string& forName, const string &dClass, const int &year, const bool &setNoValueForZero){
  changeValue(MTHREAD->GIS->pack(parName, forName, dClass, year), value_h, setNoValueForZero);
}
*/
 
double
Pixel::getDoubleValue(const string &layerName_h, const bool &returnZeroForNoValue) const{
  vIter=values.find(layerName_h);
  if(vIter != values.end()) {
    if(returnZeroForNoValue && vIter->second==MTHREAD->GIS->getNoValue()){
      return 0.0;
    } else {
      return vIter->second;
    }
  } else {
    msgOut(MSG_WARNING, "No layer \""+layerName_h+"\" found on pixel ("+i2s(getX())+","+i2s(getY())+"). Sure you didn't mispelled it?");
    if(returnZeroForNoValue){
      return 0.0;
    } else {
      return MTHREAD->GIS->getNoValue();
    }
  }
}
 
/**
getMultiplier() returns the value of the multiplier as memorized in the spatialDataSubfolder layers or in the forData table.
It will looks for exact match or for lower years if available.
If no layers are available or the usePixelData option is not used, it will return 1.
If the tp_multiplier is asked for, it will adjusts for spatial variance coefficient.
If the mortCoef_multiplier is used and we are in the table settings it will adjust it by mortCoef_link.
*/
double
Pixel::getMultiplier(const string& multiplierName, const string& forName, int year){
 
 
  if(year==DATA_NOW){year = MTHREAD->SCD->getYear();}
 
 
    double multiplierSpVar = (multiplierName == "tp_multiplier")?getSpModifier(forName):1.0;
 
    vector <string> modifiersFromTable = MTHREAD->MD->getStringVectorSetting("modifiersFromTable");
 
    if(std::find(modifiersFromTable.begin(), modifiersFromTable.end(), multiplierName) != modifiersFromTable.end()) {
      // load multiplier from forData table..
      int regId = getMyRegion()->getRegId();
      double multiplier = MTHREAD->MD->getForData(multiplierName, regId, forName, "",year);
      if (multiplierName == "mortCoef_multiplier"){
        return pow(multiplier,MTHREAD->MD->getDoubleSetting("mortMultiplier_link",DATA_NOW))*multiplierSpVar; //Added to account that our multipliers are based on probability of presence and not on planted/managed forests, where mortality is somhow reduced
      }
      return multiplier*multiplierSpVar;
 
    } else {
      // load multiplier from layer file..
      return getSTData(multiplierName, forName, year, "", 1.0);
 
//      // return 1 if not using pixel mode
//      if(!MTHREAD->MD->getBoolSetting("usePixelData")) return 1.0;
//      string search_for = multiplierName+"#"+forName+"##"+i2s(year);
//      map <string,double>::const_iterator i = values.upper_bound(search_for); //return the position always upper to the found one, even if it's an equal match.
//      if(i!= values.begin())  i--; // this rewind the position to the one just before or equal
//      const string& key = i->first;
//      string search_base = search_for.substr(0,search_for.size()-4);
//      if (key.compare(0, search_base.size(), search_base) == 0){
//          //cout << "MATCH: " << search_for <<", "<< i->first << ", " << i->second << endl;
//          //if(i->second != 1){
//          //  cout << "NOT ONE: " << search_for <<", "<< i->first << ", " << i->second << endl;
//          //  exit(0);
//          //}
//          return i->second*multiplierSpVar;
//      } else {
//          //cout << "NOTM:  " << search_for <<", "<< i->first << endl;
//          return 1.0*multiplierSpVar;
//      }
 
    }
}
 
/**
getSTData() (standing for "get spatial temporal data") returns the value memorized in a spatial layer whose name follows the above convention
 
parName#fType#dimension2#year
 
It queries the layer framework with the following logic:
 
@param parName    Parameter name (exact match)
@param forName    Forest type (exact match or look for "")
@param year       Year (exact match or look for the highest lower value). Accept enum DATA_NOW
@param d2         Optional dimension (exact match or look for "", string. Default = "")
@param naToReturn Optional value to return if no match (default to secnario-specific novalue)
 
 
*/
 
double
Pixel::getSTData(const string& parName, const string& forName, int year, const string& d2, double naToReturn){
 
  double defaultNoValue = MTHREAD->MD->getDoubleSetting("noValue");
  if (naToReturn == RETNA) naToReturn = defaultNoValue;
 
  // return NA if not using pixel mode
  if(!MTHREAD->MD->getBoolSetting("usePixelData")) return naToReturn;
 
  // If year is DATA_NOW get current year
  if (year == DATA_NOW) year = MTHREAD->SCD->getYear();
 
  // Looking for an exact match in forName
  string search_for = parName+"#"+forName+"#"+d2+"#"+i2s(year);
  map <string,double>::const_iterator i = values.upper_bound(search_for); //return the position always upper to the found one, even if it's an equal match.
  if(i!= values.begin())  i--; // this rewind the position to the one just before or equal
  const string& key = i->first;
  string search_base = search_for.substr(0,search_for.size()-4);
  if (key.compare(0, search_base.size(), search_base) == 0){
      //cout << "MATCH: " << search_for <<", "<< i->first << ", " << i->second << endl;
      //if(i->second != 1){
      //  cout << "NOT ONE: " << search_for <<", "<< i->first << ", " << i->second << endl;
      //  exit(0);
      //}
      return i->second == defaultNoValue ? naToReturn : i->second;
  } else {
      // An exact match in forName has not being found, look for no forName:
      string search_for2 = parName+"#"+"#"+d2+"#"+i2s(year);
      map <string,double>::const_iterator i2 = values.upper_bound(search_for2); //return the position always upper to the found one, even if it's an equal match.
      if(i2!= values.begin())  i2--; // this rewind the position to the one just before or equal
      const string& key2 = i2->first;
      string search_base2 = search_for2.substr(0,search_for2.size()-4);
      if (key2.compare(0, search_base2.size(), search_base2) == 0){
        //cout << "MATCH: " << search_for <<", "<< i->first << ", " << i->second << endl;
        //if(i->second != 1){
        //  cout << "NOT ONE: " << search_for <<", "<< i->first << ", " << i->second << endl;
        //  exit(0);
        //}
        return i2->second == defaultNoValue ? naToReturn : i2->second;
      } else {
          // An exact match in forName or noforName has not being found, look for no dc:
          string search_for2 = parName+"#"+forName+"#"+"#"+i2s(year);
          map <string,double>::const_iterator i2 = values.upper_bound(search_for2); //return the position always upper to the found one, even if it's an equal match.
          if(i2!= values.begin())  i2--; // this rewind the position to the one just before or equal
          const string& key2 = i2->first;
          string search_base2 = search_for2.substr(0,search_for2.size()-4);
          if (key2.compare(0, search_base2.size(), search_base2) == 0){
            //cout << "MATCH: " << search_for <<", "<< i->first << ", " << i->second << endl;
            //if(i->second != 1){
            //  cout << "NOT ONE: " << search_for <<", "<< i->first << ", " << i->second << endl;
            //  exit(0);
            //}
            return i2->second == defaultNoValue ? naToReturn : i2->second;
          } else {
              // An exact match, or partial match in forName or dc alones has not being found, look for no forName AND no dc:
              string search_for2 = parName+"#"+"#"+"#"+i2s(year);
              map <string,double>::const_iterator i2 = values.upper_bound(search_for2); //return the position always upper to the found one, even if it's an equal match.
              if(i2!= values.begin())  i2--; // this rewind the position to the one just before or equal
              const string& key2 = i2->first;
              string search_base2 = search_for2.substr(0,search_for2.size()-4);
              if (key2.compare(0, search_base2.size(), search_base2) == 0){
                //cout << "MATCH: " << search_for <<", "<< i->first << ", " << i->second << endl;
                //if(i->second != 1){
                //  cout << "NOT ONE: " << search_for <<", "<< i->first << ", " << i->second << endl;
                //  exit(0);
                //}
                return i2->second == defaultNoValue ? naToReturn : i2->second;
              } else {
                 //cout << "NOTM:  " << search_for <<", "<< i->first << endl;
                 return naToReturn;
              }
           }
      }
  }
}
 
 
/**
The mortality returned is the increased yearly mortality due to any affecting pathogenes.
The function load the relevant pathogen mortality rule(s), for each of them check for how many years the phatogen is present with concentrations
above the threshold and returns the relavant increase in mortality (summing them in case of multiple pathogens).
 
*/
double
Pixel::getPathMortality(const string& forType, const string& dC, int year){
  if(!MTHREAD->MD->getBoolSetting("usePathogenModule")) return 0.0;
 
  string debug=forType;
  int initialOptYear = MTHREAD->MD->getIntSetting("initialOptYear");
  int simulationYears = MTHREAD->MD->getIntSetting("simulationYears");
 
  int maxYear = initialOptYear + simulationYears;
 
  vector<pathRule*> pathRules = MTHREAD->MD->getPathMortalityRule(forType,dC);
 
  double pathMort = 0.0;
  if(year==DATA_NOW){year = MTHREAD->SCD->getYear();}
 
  for(uint r=0;r<pathRules.size();r++){
    string pathId=pathRules[r]->pathId;
    double pres_min=pathRules[r]->pres_min;
    vector<double> mortCoefficients=pathRules[r]->mortCoefficents;
    double pathMort_thispath = 0.0;
    for(uint y=year;y>(year-mortCoefficients.size());y--){
      int i =year-y;
      int y2 = y;
      if(y>=maxYear){
        y2=maxYear-1;
      }
 
      string layerName="pathogen_pp#"+pathId+"#"+i2s(y2);
      if(MTHREAD->GIS->layerExist(layerName)){
        if (this->getDoubleValue(layerName,true)>= pres_min){
          pathMort_thispath = mortCoefficients[i];
        }
      }
    }
    pathMort += pathMort_thispath;
  }
  return pathMort;
 
}
 
void
Pixel::correctInputMultiplier (const string& multiplierName, const string& forName, double coefficient){
  string search_for = multiplierName+"#"+forName+"#";
  for (std::map<string,double>::iterator it=values.lower_bound(search_for); it!=values.end(); ++it){
    if (it->first.compare(0, search_for.size(), search_for) == 0){
            //cout << ID << ";" << forName << ";" << coefficient << endl;
      it->second = it->second * coefficient;
    }
  }
}
 
double
Pixel::getDoubleValue (const string& parName, const string& forName, const string& dClass, const int& year, const  bool& returnZeroForNoValue){
  return getDoubleValue(MTHREAD->GIS->pack(parName, forName, dClass, year), returnZeroForNoValue);
}
 
void
Pixel::newYear(){
 
}
 
double
Pixel::getPastRegArea(const int& ft_idx, const int& year){
  map <int,vector<double> >::const_iterator i=regArea.find(year);
  if(i != regArea.end()) {
    return i->second.at(ft_idx);
  } else {
    msgOut(MSG_ERROR, "Asking for a pastRegArea of a not-registered year. I don't have year "+i2s(year)+"!");
  }
}
 
void
Pixel::setPastRegArea(const double& value, const int& ft_idx, const int& year){
  msgOut(MSG_CRITICAL_ERROR,"TODO");
  /*map <int,vector<double> >::const_iterator i=regArea.find(year);
  if(i != regArea.end()) {
    // we already have this year, let's see if the vector is bif enough
    int currsize = i->second.size();
    for(j=0;j<ft_idx-currside;j++){
 
    }
    return i->second.at(ft_idx);
  } else {
    // new year
  }
 
 
  pair<int,vector<double> newRegArea;
  */
 
 
}
 
void
Pixel::swap(const int& swap_what){
  switch (swap_what){
    case VAR_VOL:
      vol_l = vol;
      break;
    case VAR_AREA:
      area_l = area;
      break;
    default:
      msgOut(MSG_CRITICAL_ERROR,"Don't know how to swap "+swap_what);
  }
}
 
 
double
Pixel::getSpModifier(const string& ft){
  vector<string>ftypes = MTHREAD->MD->getForTypeIds();
  for (int i=0;i<ftypes.size();i++){
      if (ftypes[i] == ft){
          return spMods.at(i);
      }
  }
  msgOut(MSG_CRITICAL_ERROR,"Asked spatial modifier for a forest type that doesn't exist");
 
}
 
ModelRegion*
Pixel::getMyRegion(const int& rLevel){
    if(rLevel==2){
        return l2region;
    } else if (rLevel==1) {
        return l2region->getParent();
    } else {
        msgOut(MSG_ERROR, "Requested a unknown level region code in getMyRegion().");
    }
}