Pches statesandconstraints

NAMESPACE

@@ -76,4 +76,5 @@ importFrom(tibble,tibble)
 importFrom(tidyr,gather)
 importFrom(tidyr,separate)
 importFrom(utils,capture.output)
+importFrom(utils,read.csv)
 importFrom(utils,sessionInfo)

R/ag_production_technology.R

@@ -9,8 +9,10 @@
 #' @param aLandLeaf Land leaf
 #' @param aPeriod Model time period.
 #' @param aScenarioInfo Scenario-related information, including names, logits, expectations
+#' @param aLandConstraintCost Additional cost used to constrain land
+#' @param aLandConstraintString String identifying which leafs are included in the constraint
 #' @author KVC September 2017
-AgProductionTechnology_initCalc <- function(aLandLeaf, aPeriod, aScenarioInfo) {
+AgProductionTechnology_initCalc <- function(aLandLeaf, aPeriod, aScenarioInfo, aLandConstraintCost, aLandConstraintString) {
   # First, set the nonLandVariableCost for this technology in this period.
   # If no nonLandVariableCost is read in, get the previous period cost.
   # Note: you can never overwrite a positive cost with a zero cost. If the model sees a
@@ -66,6 +68,13 @@ AgProductionTechnology_initCalc <- function(aLandLeaf, aPeriod, aScenarioInfo) {
     aLandLeaf$mYield[aPeriod] <- 0.0
   }
 
+  # Set constraint cost for this period
+  if( aLandConstraintString != "" & grepl(aLandConstraintString, aLandLeaf$mName[1]) ) {
+    aLandLeaf$mConstraintCost[aPeriod] <- aLandConstraintCost
+  } else if (length(aLandLeaf$mConstraintCost) < aPeriod) {
+    aLandLeaf$mConstraintCost[aPeriod] <- 0.0
+  }
+
   # Calculate profit rate
   AgProductionTechnology_calcProfitRate(aLandLeaf, aPeriod, aScenarioInfo)
 }
@@ -108,4 +117,7 @@ AgProductionTechnology_calcProfitRate <- function(aLandLeaf, aPeriod, aScenarioI
     # Subsidies are assumed to br in $/billion m2
     aLandLeaf$mProfitRate[aPeriod] <- aLandLeaf$mProfitRate[[aPeriod]] + aLandLeaf$mSubsidy[[aPeriod]]
   }
+
+  # Subtract constraint cost from profit
+  aLandLeaf$mProfitRate[aPeriod] <- aLandLeaf$mProfitRate[[aPeriod]] - aLandLeaf$mConstraintCost[[aPeriod]]
 }

R/generate_price_data.R

@@ -27,18 +27,19 @@ get_prices <- function(aScenType) {
         prices
     }
   } else {
-    file <- paste("./scenario-data/AgPrices_", aScenType, ".csv", sep="")
-    prices <- suppressMessages(read_csv(system.file("extdata", file, package = "gcamland"), skip = 1))
 
     # Tidy data
     if(aScenType == "PCHES") {
+      prices <- read_csv("./inst/extdata/scenario-data/AgPrices_PCHES.csv", skip=1)
       # PCHES scenarios have subregional price information
       prices %>%
         select(-scenario, -Units) %>%
         gather(year, price, -region, -sector, -subregion) %>%
-        mutate(year = as.integer(year)) ->
+        mutate(year = as.integer(substr(year, 2, 5))) ->
         prices
     } else {
+      file <- paste("./scenario-data/AgPrices_", aScenType, ".csv", sep="")
+      prices <- suppressMessages(read_csv(system.file("extdata", file, package = "gcamland"), skip = 1))
       prices %>%
         select(-scenario, -Units) %>%
         gather(year, price, -region, -sector) %>%

R/generate_scenario_info.R

@@ -39,6 +39,8 @@ DEFAULT.SCENARIO.TYPE <- "Reference"
 #' @param aRegion Region to use in the calculation.  Right now we only run a
 #' single region at a time.
 #' @param aSubRegion Subregion name. Note we can only run a full region or a subregion not both
+#' @param aIncludeConstraint Boolean flag indicating that we are constraining some sort of land
+#' @param aConstraintFile String with the name/path to the file with constraints
 #' @return New ScenarioInfo object
 #' @export
 #' @author KVC November 2017
@@ -64,7 +66,9 @@ ScenarioInfo <- function(# Currently only "Perfect", "Linear", "Adaptive", "Hybr
                          aOutputDir = "./outputs",
                          aSerialNum = NA,
                          aRegion = DEFAULT.REGION,
-                         aSubRegion = NULL) {
+                         aSubRegion = NULL,
+                         aIncludeConstraint = FALSE,
+                         aConstraintFile = "") {
 
   self <- new.env(parent=emptyenv())
   class(self) <- c("ScenarioInfo", class(self))
@@ -95,6 +99,8 @@ ScenarioInfo <- function(# Currently only "Perfect", "Linear", "Adaptive", "Hybr
   self$mPointwiseLikelihood <- data.frame() # actually log-likelihood, tabulated
                                         # by data point.
   self$mLogPost <- data.frame()
+  self$mIncludeConstraint <- aIncludeConstraint
+  self$mConstraintFile <- aConstraintFile
 
   self
 }
@@ -184,7 +190,9 @@ PCHES.SCENARIO.INFO <- ScenarioInfo(aScenarioType = "PCHES",
                               aUseZeroCost = FALSE,
                               aCalibrateShareWt = TRUE,
                               aScenarioName = paste0(DEFAULT.SCENARIO.TYPE, "_", "Perfect", "_", "PCHES"),
-                              aFileName = paste0(DEFAULT.SCENARIO.TYPE, "_", "Perfect", "_", "PCHES"))
+                              aFileName = paste0(DEFAULT.SCENARIO.TYPE, "_", "Perfect", "_", "PCHES"),
+                              aIncludeConstraint = FALSE,
+                              aConstraintFile = "./scenario-data/PCHES_constraint.csv")
 
 #' SRB.SCENARIO.INFO
 #'

R/helpers.R

@@ -112,3 +112,73 @@ sumx <- function(x)
 {
     sum(x[order(abs(x))])
 }
+
+#' LandAllocator_getLandAreaByCriteria
+#'
+#' @details Calculate total land area that meets a particular criteria.
+#'          Criteria is currently determined by a string appearing in the
+#'          name of the land type
+#' @param aLandAllocator Land allocator
+#' @param aString String to group land area by
+#' @param aPeriod Model period to get land area for
+#'
+#' @return Land area
+LandAllocator_getLandAreaByCriteria <- function(aLandAllocator, aString, aPeriod) {
+  # Set total land to 0
+  totalLand <- 0
+
+  for(child in aLandAllocator$mChildren) {
+    if(inherits(child, "LandNode")) {
+      totalLand <- totalLand + LandNode_getLandAreaByCriteria(child, aString, aPeriod)
+    } else {
+      totalLand <- totalLand + LandLeaf_getLandAreaByCriteria(child, aString, aPeriod)
+    }
+  }
+
+  return(totalLand)
+}
+
+#' LandNode_getLandAreaByCriteria
+#'
+#' @details Calculate total land area that meets a particular criteria.
+#'          Criteria is currently determined by a string appearing in the
+#'          name of the land type
+#' @param aLandNode Land allocator
+#' @param aString String to group land area by
+#' @param aPeriod Model period to get land area for
+#'
+#' @return Land area
+LandNode_getLandAreaByCriteria <- function(aLandNode, aString, aPeriod) {
+  # Set land to zero for this node
+  land <- 0
+
+  # Loop over children, calling their `getLandAreaByCriteria` function
+  for(child in aLandNode$mChildren) {
+    if(inherits(child, "LandNode")) {
+      land <- land + LandNode_getLandAreaByCriteria(child, aString, aPeriod)
+    } else {
+      land <- land + LandLeaf_getLandAreaByCriteria(child, aString, aPeriod)
+    }
+  }
+
+  return(land)
+}
+
+
+#' LandLeaf_getLandAreaByCriteria
+#'
+#' @details Calculate total land area that meets a particular criteria.
+#'          Criteria is currently determined by a string appearing in the
+#'          name of the land type
+#' @param aLandLeaf Land allocator
+#' @param aString String to group land area by
+#' @param aPeriod Model period to get land area for
+#'
+#' @return Land area for this leaf
+LandLeaf_getLandAreaByCriteria <- function(aLandLeaf, aString, aPeriod) {
+  if(grepl(aString, aLandLeaf$mName[1])) {
+    return(aLandLeaf$mLandAllocation[[aPeriod]])
+  } else {
+    return(0)
+  }
+}

R/land_leaf.R

@@ -26,6 +26,7 @@
 #' @field mCalOutput calibrated-output for this leaf
 #' @field mNonLandCostTechChange Technical change on the cost of this leaf
 #' @field mAgProdChange Technical change on yield for this leaf
+#' @field mConstraintCost Cost added to ensure constraints are met
 #'
 #' @return New, initialized LandLeaf
 #' @author KVC September 2017
@@ -52,6 +53,7 @@ LandLeaf <- function(aName, aFinalCalPeriod, aFinalPeriod) {
   self$mCalOutput = list() # Note: we read this in to calculate yield for consistency with C++
   self$mNonLandCostTechChange = list()
   self$mAgProdChange = list()
+  self$mConstraintCost = list()
   self$mSubsidy = list()
 
   class(self) <- c("LandLeaf", class(self))

R/main.R

@@ -517,13 +517,16 @@ gen_ensemble_member <- function(agFor, agForNonPast, crop, share1, share2, share
 #' @return Table of model results.
 #' @author KVC
 #' @importFrom assertthat assert_that has_attr
+#' @importFrom utils read.csv
 #' @export
 #' @examples
 #' \dontrun{
 #' run_model(SCENARIO.INFO, aVerbose=TRUE)
 #' run_model(SCENARIO.INFO, aPeriods = 1:5)
 #' }
 run_model <- function(aScenarioInfo, aPeriods=NULL, aVerbose=FALSE, agData=NULL) {
+  # Silence package checks
+  read.csv <- period <- NULL
 
   #### Step 1: Setup
   # Ensure that output directories exist
@@ -558,6 +561,47 @@ run_model <- function(aScenarioInfo, aPeriods=NULL, aVerbose=FALSE, agData=NULL)
     #### Step 3: Final calculation
     # Next, call calcFinalLandAllocation for LandAllocator
     LandAllocator_calcFinalLandAllocation(mLandAllocator, per, aScenarioInfo)
+
+    ### Step 4: Check that constraints are met
+    if ( aScenarioInfo$mIncludeConstraint & per > TIME.PARAMS[[aScenarioInfo$mScenarioType]]$FINAL_CALIBRATION_PERIOD) {
+      read.csv(system.file("extdata", aScenarioInfo$mConstraintFile, package = "gcamland")) %>%
+        filter(period == per) ->
+        constraints
+
+      if(nrow(constraints) > 0) {
+        # Loop through all constraints for this period
+        for(i in 1:nrow(constraints)) {
+          message("Constraining: ", constraints$string[i])
+          landConstraintCost <- 0
+          while( LandAllocator_getLandAreaByCriteria(mLandAllocator, constraints$string[i], per) > constraints$value[i] ) {
+            currValue <- LandAllocator_getLandAreaByCriteria(mLandAllocator, constraints$string[i], per)
+            message("Current value: ", currValue, ", Target value:", constraints$value[i])
+            if ( currValue < 2*constraints$value[i] ) {
+              # Use smaller increments when we get close to solution
+              landConstraintCost <- landConstraintCost + 1e8
+            } else if ( currValue > 20*constraints$value[i] ) {
+              # Use smaller increments when we get close to solution
+              landConstraintCost <- landConstraintCost + 1e10
+            } else {
+              landConstraintCost <- landConstraintCost + 1e9
+            }
+
+            # Constraint not met. Adjust cost and recalculate
+            Sector_initCalc(mLandAllocator, per, aScenarioInfo, landConstraintCost, constraints$string[i])
+            LandAllocator_initCalc(mLandAllocator, per, aScenarioInfo)
+            LandAllocator_calcFinalLandAllocation(mLandAllocator, per, aScenarioInfo)
+
+            # Check that we are making progress toward meeting constraint
+            if( LandAllocator_getLandAreaByCriteria(mLandAllocator, constraints$string[i], per) > 0.999999*currValue) {
+              message("Progress is not being made. New value (",
+                      LandAllocator_getLandAreaByCriteria(mLandAllocator, constraints$string[i], per),
+                      ") is less than 0.0001% smaller than old value (", currValue, ")")
+              break
+            }
+          } # End while loop to meet constraint
+        } # End for loop over constraints
+      }
+    }
   }
 
   #### Step 4: Reporting

R/perfect_expectation.R

@@ -39,8 +39,8 @@ PerfectExpectation_calcExpectedPrice <- function(aLandLeaf, aPeriod, aScenarioIn
       # PCHES uses subregional prices. First, figure out the subregion by decomposing the land leaf name
       # Note: can't use separate on this because some subregions have "_" in their name.
       subRegion <- sub(paste(aLandLeaf$mProductName, "_", sep=""), "", aLandLeaf$mName)
-      subRegion <- sub("_IRR", "", subRegion)
-      subRegion <- sub("_RFD", "", subRegion)
+      subRegion <- sub("_Irrigated", "", subRegion)
+      subRegion <- sub("_Rainfed", "", subRegion)
       expectedPrice <- price_table$price[price_table$year == y & price_table$sector == aLandLeaf$mProductName[1]
                                          & price_table$subregion == subRegion]
     } else {

R/sector.R

@@ -9,11 +9,13 @@
 #' @param aLandAllocator Land allocator to perform initializations on
 #' @param aPeriod Current model time period
 #' @param aScenarioInfo Scenario-related information, including names, logits, expectations
+#' @param aLandConstraintCost Additional cost used to constrain land
+#' @param aLandConstraintString String identifying which leafs are included in the constraint
 #' @author KVC September 2017
-Sector_initCalc <- function(aLandAllocator, aPeriod, aScenarioInfo){
+Sector_initCalc <- function(aLandAllocator, aPeriod, aScenarioInfo, aLandConstraintCost = 0.0, aLandConstraintString = ""){
   # Loop through the land allocator
   for(child in aLandAllocator$mChildren) {
-    Subsector_initCalc(child, aPeriod, aScenarioInfo)
+    Subsector_initCalc(child, aPeriod, aScenarioInfo, aLandConstraintCost, aLandConstraintString)
   }
 }
 
@@ -26,15 +28,17 @@ Sector_initCalc <- function(aLandAllocator, aPeriod, aScenarioInfo){
 #' @param aLandNode Land node to perform initializations on
 #' @param aPeriod Current model time period
 #' @param aScenarioInfo Scenario-related information, including names, logits, expectations
+#' @param aLandConstraintCost Additional cost used to constrain land
+#' @param aLandConstraintString String identifying which leafs are included in the constraint
 #' @author KVC September 2017
-Subsector_initCalc <- function(aLandNode, aPeriod, aScenarioInfo){
+Subsector_initCalc <- function(aLandNode, aPeriod, aScenarioInfo, aLandConstraintCost, aLandConstraintString){
   # loop through children
   for(child in aLandNode$mChildren) {
     if(inherits(child, "LandNode")) {
-      Subsector_initCalc(child, aPeriod, aScenarioInfo)
+      Subsector_initCalc(child, aPeriod, aScenarioInfo, aLandConstraintCost, aLandConstraintString)
     } else if (inherits(child, "LandLeaf")) {
       # If this is a leaf, call the AgProductionTechnology method
-      AgProductionTechnology_initCalc(child, aPeriod, aScenarioInfo)
+      AgProductionTechnology_initCalc(child, aPeriod, aScenarioInfo, aLandConstraintCost, aLandConstraintString)
     }
   }
 }