From a89a1ecee977fe2d26bbd15bd48aa0b5b33e455c Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Sun, 2 Feb 2025 17:25:20 +0100
Subject: [PATCH 01/14] update calibration cropland land conversion cost FAO

---
 config/default.cfg                        |  8 ++-
 scripts/calibration/landconversion_cost.R | 72 ++++++++++++++---------
 scripts/start/extra/recalibrate_default.R |  2 +-
 scripts/start_functions.R                 |  3 +-
 4 files changed, 53 insertions(+), 32 deletions(-)

diff --git a/config/default.cfg b/config/default.cfg
index eff8e6a56..3f843bfc5 100644
--- a/config/default.cfg
+++ b/config/default.cfg
@@ -95,7 +95,7 @@ cfg$recalibrate_landconversion_cost <- "ifneeded" #def "ifneeded"
 # Up to which accuracy shall be recalibrated?
 cfg$calib_accuracy_landconversion_cost <- 0.01         # def = 0.01
 # What is the maximum number of iterations if the precision goal is not reached?
-cfg$calib_maxiter_landconversion_cost <- 40           # def = 40
+cfg$calib_maxiter_landconversion_cost <- 20           # def = 20
 # Restart from existing calibration factors (TRUE or FALSE)
 cfg$restart_landconversion_cost <- FALSE           # def = FALSE
 # Number of lowpass filter iterations applied on calibration factors
@@ -104,11 +104,15 @@ cfg$lowpass_filter_landconversion_cost <- 1        # def= 1
 # Set upper limit for cropland calibration factor
 cfg$cost_calib_max_landconversion_cost <- 3         # def= 3
 # Set lower limit for cropland calibration factor
-cfg$cost_calib_min_landconversion_cost <- 0.05            # def= 0.05
+cfg$cost_calib_min_landconversion_cost <- 0.2            # def= 0.2
 # Selection type of calibration factors.
 # If FALSE, calibration factors from the last iteration are used.
 # If TRUE, calibration factors from the iteration with the lowest divergence are used.
 cfg$best_calib_landconversion_cost <- FALSE     # def = FALSE
+# Target data set that will be used for cropland calibration at regional level 
+# * "MAgPIEown": same data set as used for initialization of cropland
+# * "FAO": Data from FAOSTAT on cropland area
+cfg$cost_calib_hist_data <- "FAO"            # def = "FAO"
 
 # Settings for NPI/NDC recalculation
 # * (TRUE): NPI/NDC recalculation will be performed
diff --git a/scripts/calibration/landconversion_cost.R b/scripts/calibration/landconversion_cost.R
index 12ab1bb0e..5a103467c 100644
--- a/scripts/calibration/landconversion_cost.R
+++ b/scripts/calibration/landconversion_cost.R
@@ -35,15 +35,27 @@ calibration_run <- function(putfolder, calib_magpie_name, logoption = 3, s_use_g
 
 # get ratio between modelled area and reference area
 
-getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpass_filter = 1) {
+getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpass_filter = 1, histData = "FAO") {
   require(magclass)
   require(magpie4)
   require(gdx2)
   y <- readGDX(gdx_file,"t")
   magpie <- land(gdx_file)[, y, "crop"]
-  hist <- dimSums(readGDX(gdx_file, "f10_land")[, , "crop"], dim = 1.2)
-  data <- hist[, y, "crop"]
-  shrLost <- (setYears(hist[getRegions(magpie), 2015, ], NULL) - setYears(hist[getRegions(magpie), 1995, ], NULL)) / setYears(hist[getRegions(magpie), 1995, ], NULL)
+  if (histData == "MAgPIEown") {
+    hist <- dimSums(readGDX(gdx_file, "f10_land")[, , "crop"], dim = 1.2)
+    data <- hist[, y, "crop"]
+  } else if (histData == "FAO") {
+    if(file.exists("calib_data.rds")) {
+      val <- readRDS("calib_data.rds")
+    } else {
+      val <- read.report("input/validation.mif", as.list = FALSE)
+      val <- val[getRegions(magpie),getYears(magpie),"historical.FAO_crop_past.Resources|Land Cover|+|Cropland (million ha)"]
+      names(dimnames(val)) <- names(dimnames(magpie))
+      getNames(val) <- "crop"
+      saveRDS(val, file = "calib_data.rds")
+    }
+    data <- val
+  }
   if (nregions(magpie) != nregions(data) | !all(getRegions(magpie) %in% getRegions(data))) {
     stop("Regions in MAgPIE do not agree with regions in reference calibration area data set!")
   }
@@ -54,16 +66,19 @@ getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpa
     getNames(out) <- NULL
 
     out[which(out < 0, arr.ind = T)] <- 1
+    out <- lowpass(out,i = lowpass_filter)
   } else if (mode == "reward") {
-    out <- magpie / data - 1
+    shrLostProj <- (setYears(magpie[, 2015, ], NULL) - setYears(magpie[, 1995, ], NULL)) / setYears(magpie[, 1995, ], NULL)
+    shrLostHist <- (setYears(data[getRegions(magpie), 2015, ], NULL) - setYears(data[getRegions(magpie), 1995, ], NULL)) / setYears(data[getRegions(magpie), 1995, ], NULL)
+    out <- shrLostHist/shrLostProj
     out[is.na(out)] <- 0
+    out[is.infinite(out)] <- 0
     getNames(out) <- NULL
-
+    
     # only reward if share of cropland lost between 1995 and 2015 exceeds a certain threshold. Otherwise set to 0.
-    out[rownames(which(shrLost > -calib_accuracy, arr.ind = T)),,] <- 0
+    out[rownames(which(shrLostHist > -calib_accuracy, arr.ind = T)),,] <- 0
     out[which(out < 0, arr.ind = T)] <- 0
   }
-  out <- lowpass(out,i = lowpass_filter)
   return(magpiesort(out))
 }
 
@@ -79,19 +94,13 @@ time_series_cost <- function(calib_factor) {
 
 time_series_reward <- function(calib_factor) {
   out2 <- new.magpie(getRegions(calib_factor), years = c(seq(1995, 2015, by = 5), seq(2050, 2150, by = 5)), fill = 0)
-  out2[, getYears(calib_factor), ] <- calib_factor
+  out2[, seq(2000, 2015, by = 5), ] <- calib_factor
   out2 <- time_interpolate(out2, seq(2020, 2050, by = 5), integrate_interpolated_years = T)
   return(out2)
 }
 
-getHistCrop <- function() {
-  rep <- read.report("input/validation.mif", as.list = FALSE)
-  crop <- collapseNames(rep[, , "historical.FAO_crop_past.Resources|Land Cover|+|Cropland (million ha)"])
-  return(crop)
-}
-
 # Calculate the correction factor and save it
-update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, calib_file, cost_max = 3, cost_min = 0.05, calibration_step = "", n_maxcalib = 40, best_calib = FALSE) {
+update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, calib_file, cost_max = 3, cost_min = 0.2, calibration_step = "", n_maxcalib = 20, best_calib = FALSE, histData = "FAO") {
   require(magclass)
   require(magpie4)
   require(gdx2)
@@ -103,7 +112,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
   calib_divergence_cost <- abs(calib_correction_cost - 1)
 
   calib_correction_reward <- getCalibFactor(gdx_file, mode = "reward", calib_accuracy = calib_accuracy, lowpass_filter = lowpass_filter)
-  calib_divergence_reward <- abs(calib_correction_reward)
+  calib_divergence_reward <- abs(calib_correction_reward - 1)
   
   ### -> in case it is the first step, it forces the initial factors to be equal to 1
   if (file.exists(calib_file)) {
@@ -111,12 +120,11 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
     start_flag <- FALSE
   } else {
     old_calib <- new.magpie(cells_and_regions = getCells(calib_divergence_cost), years = y, names = c("cost", "reward"), fill = 1)
-    old_calib[,,"reward"] <- 0
     start_flag <- TRUE
   }
 
   calib_factor_cost <- setNames(old_calib[, , "cost"], NULL) * calib_correction_cost
-  calib_factor_reward <- setNames(old_calib[, , "reward"], NULL) + calib_correction_reward
+  calib_factor_reward <- setNames(old_calib[, , "reward"], NULL) * calib_correction_reward
   calib_factor_reward[calib_factor_reward < 0] <- 0
   
   if (!start_flag) {
@@ -124,9 +132,11 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
     # use stricter divergence threshold in first 5 calibration_step steps
     cost_acc_reached <- calib_divergence_cost <= ifelse(calibration_step < 6, 0.01, calib_accuracy)
     calib_factor_cost[cost_acc_reached] <- setNames(old_calib[, , "cost"], NULL)[cost_acc_reached]
-
+    calib_correction_cost[cost_acc_reached] <- 1
+    
     reward_acc_reached <- calib_divergence_reward <= ifelse(calibration_step < 6, 0.01, calib_accuracy)
     calib_factor_reward[reward_acc_reached] <- setNames(old_calib[, , "reward"], NULL)[reward_acc_reached]
+    calib_correction_reward[reward_acc_reached] <- 0
   }
 
   if (!is.null(cost_max)) {
@@ -139,9 +149,11 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
       calib_factor_cost["IND", , ] <- cost_max
     }
     
-    # set value for CAZ to cost_max for better convergence
-    if ("CAZ" %in% getRegions(calib_factor_cost)) {
-      calib_factor_cost["CAZ", , ] <- cost_max
+    # set value for CAZ to cost_max for better convergence; only needed for LUH2 due to mismatch between LUH2 and FAO historical data for cropland.
+    if (histData == "MAgPIEown") {
+      if ("CAZ" %in% getRegions(calib_factor_cost)) {
+        calib_factor_cost["CAZ", , ] <- cost_max
+      }
     }
     
     above_limit <- (calib_factor_cost >= cost_max)
@@ -201,6 +213,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
         " unit: -",
         paste0(" note: Best calibration factor from the run"),
         " origin: scripts/calibration/landconversion_cost.R (path relative to model main directory)",
+        paste(" Calibration settings:",  "calib_accuracy=", calib_accuracy, "lowpass_filter=", lowpass_filter, "cost_max=", cost_max, "cost_min=", cost_min, "n_maxcalib=",n_maxcalib, "best_calib=",best_calib, "histData=",histData),
         paste0(" creation date: ", date())
       )
       write.magpie(round(calib_best_full, 3), calib_file, comment = comment)
@@ -227,6 +240,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
       " unit: -",
       paste0(" note: Calibration step ", calibration_step),
       " origin: scripts/calibration/landconversion_cost.R (path relative to model main directory)",
+      paste(" Calibration settings:",  "calib_accuracy=", calib_accuracy, "lowpass_filter=", lowpass_filter, "cost_max=", cost_max, "cost_min=", cost_min, "n_maxcalib=",n_maxcalib, "best_calib=",best_calib, "histData=",histData),
       paste0(" creation date: ", date())
     )
 
@@ -238,11 +252,11 @@ update_calib <- function(gdx_file, calib_accuracy = 0.05, lowpass_filter = 1, ca
 
 
 
-calibrate_magpie <- function(n_maxcalib = 40,
-                             restart = TRUE,
-                             calib_accuracy = 0.05,
+calibrate_magpie <- function(n_maxcalib = 20,
+                             restart = FALSE,
+                             calib_accuracy = 0.01,
                              cost_max = 3,
-                             cost_min = 0.05,
+                             cost_min = 0.2,
                              calib_magpie_name = "magpie_calib",
                              lowpass_filter = 1,
                              calib_file = "modules/39_landconversion/input/f39_calib.csv",
@@ -250,7 +264,8 @@ calibrate_magpie <- function(n_maxcalib = 40,
                              data_workspace = NULL,
                              logoption = 3,
                              debug = FALSE,
-                             best_calib = FALSE) {
+                             best_calib = FALSE,
+                             histData = "FAO") {
   require(magclass)
 
   if (!restart) {
@@ -280,6 +295,7 @@ calibrate_magpie <- function(n_maxcalib = 40,
   # delete calib_magpie_gms in the main folder
   unlink(paste0(calib_magpie_name, ".*"))
   unlink("fulldata.gdx")
+  unlink("calib_data.rds")
 
   cat("\nLand conversion cost calibration finished\n")
 }
diff --git a/scripts/start/extra/recalibrate_default.R b/scripts/start/extra/recalibrate_default.R
index 78de9855e..b792deb4e 100644
--- a/scripts/start/extra/recalibrate_default.R
+++ b/scripts/start/extra/recalibrate_default.R
@@ -25,4 +25,4 @@ cfg$output <- c("rds_report")
 cfg$force_replace <- TRUE
 cfg$qos <- "priority"
 start_run(cfg,codeCheck=FALSE)
-magpie4::submitCalibration("H12")
+magpie4::submitCalibration("H12_FAO")
diff --git a/scripts/start_functions.R b/scripts/start_functions.R
index 24ae800f6..cfc129755 100644
--- a/scripts/start_functions.R
+++ b/scripts/start_functions.R
@@ -476,7 +476,8 @@ start_run <- function(cfg, scenario = NULL, codeCheck = TRUE, lock_model = TRUE,
                      data_workspace = cfg$val_workspace,
                      logoption = 3,
                      debug = cfg$debug,
-                     best_calib = cfg$best_calib_landconversion_cost)
+                     best_calib = cfg$best_calib_landconversion_cost,
+                     histData = cfg$cost_calib_hist_data)
     cat("Land conversion cost calibration factor calculated!\n")
   }
 

From d29c1b5b05e66afff59db7ed3208fcb67c8e5da7 Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Sun, 2 Feb 2025 18:02:38 +0100
Subject: [PATCH 02/14] bugfix

---
 scripts/calibration/landconversion_cost.R | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/scripts/calibration/landconversion_cost.R b/scripts/calibration/landconversion_cost.R
index 5a103467c..0e955b00d 100644
--- a/scripts/calibration/landconversion_cost.R
+++ b/scripts/calibration/landconversion_cost.R
@@ -94,7 +94,7 @@ time_series_cost <- function(calib_factor) {
 
 time_series_reward <- function(calib_factor) {
   out2 <- new.magpie(getRegions(calib_factor), years = c(seq(1995, 2015, by = 5), seq(2050, 2150, by = 5)), fill = 0)
-  out2[, seq(2000, 2015, by = 5), ] <- calib_factor
+  out2[, getYears(calib_factor), ] <- calib_factor
   out2 <- time_interpolate(out2, seq(2020, 2050, by = 5), integrate_interpolated_years = T)
   return(out2)
 }
@@ -125,6 +125,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
 
   calib_factor_cost <- setNames(old_calib[, , "cost"], NULL) * calib_correction_cost
   calib_factor_reward <- setNames(old_calib[, , "reward"], NULL) * calib_correction_reward
+  calib_factor_reward[,1995,] <- 0
   calib_factor_reward[calib_factor_reward < 0] <- 0
   
   if (!start_flag) {
@@ -159,6 +160,10 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
     above_limit <- (calib_factor_cost >= cost_max)
     calib_factor_cost[above_limit] <- cost_max
     calib_divergence_cost[getRegions(calib_factor_cost), , ][above_limit] <- 0
+    
+    above_limit <- (calib_factor_reward >= cost_max)
+    calib_factor_reward[above_limit] <- cost_max
+    calib_divergence_reward[getRegions(calib_factor_reward), , ][above_limit] <- 0
   }
 
   if (!is.null(cost_min)) {

From 7acf6dcf08c473ff77f1757c930a8f373c409284 Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Sun, 2 Feb 2025 19:02:25 +0100
Subject: [PATCH 03/14] bugfix

---
 scripts/calibration/landconversion_cost.R | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/scripts/calibration/landconversion_cost.R b/scripts/calibration/landconversion_cost.R
index 0e955b00d..73871ec6d 100644
--- a/scripts/calibration/landconversion_cost.R
+++ b/scripts/calibration/landconversion_cost.R
@@ -70,7 +70,8 @@ getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpa
   } else if (mode == "reward") {
     shrLostProj <- (setYears(magpie[, 2015, ], NULL) - setYears(magpie[, 1995, ], NULL)) / setYears(magpie[, 1995, ], NULL)
     shrLostHist <- (setYears(data[getRegions(magpie), 2015, ], NULL) - setYears(data[getRegions(magpie), 1995, ], NULL)) / setYears(data[getRegions(magpie), 1995, ], NULL)
-    out <- shrLostHist/shrLostProj
+    out <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)
+    out[ , seq(2000, 2015, by = 5), ] <- shrLostHist / shrLostProj
     out[is.na(out)] <- 0
     out[is.infinite(out)] <- 0
     getNames(out) <- NULL
@@ -113,6 +114,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
 
   calib_correction_reward <- getCalibFactor(gdx_file, mode = "reward", calib_accuracy = calib_accuracy, lowpass_filter = lowpass_filter)
   calib_divergence_reward <- abs(calib_correction_reward - 1)
+  calib_divergence_reward[calib_correction_reward == 0] <- 0
   
   ### -> in case it is the first step, it forces the initial factors to be equal to 1
   if (file.exists(calib_file)) {
@@ -125,9 +127,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
 
   calib_factor_cost <- setNames(old_calib[, , "cost"], NULL) * calib_correction_cost
   calib_factor_reward <- setNames(old_calib[, , "reward"], NULL) * calib_correction_reward
-  calib_factor_reward[,1995,] <- 0
-  calib_factor_reward[calib_factor_reward < 0] <- 0
-  
+
   if (!start_flag) {
     # use calibration factors where accuracy was reached
     # use stricter divergence threshold in first 5 calibration_step steps

From bf6fe1208f4c2226560527cbd44a8e2b8af06071 Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Sun, 2 Feb 2025 21:48:31 +0100
Subject: [PATCH 04/14] bugfix

---
 scripts/calibration/landconversion_cost.R | 25 +++++++++++++----------
 1 file changed, 14 insertions(+), 11 deletions(-)

diff --git a/scripts/calibration/landconversion_cost.R b/scripts/calibration/landconversion_cost.R
index 73871ec6d..44ac1fcb8 100644
--- a/scripts/calibration/landconversion_cost.R
+++ b/scripts/calibration/landconversion_cost.R
@@ -64,22 +64,25 @@ getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpa
     out[out == 0] <- 1
     out[is.na(out)] <- 1
     getNames(out) <- NULL
-
-    out[which(out < 0, arr.ind = T)] <- 1
-    out <- lowpass(out,i = lowpass_filter)
+    out[out < 0] <- 1
+    out[,1,] <- 0
   } else if (mode == "reward") {
-    shrLostProj <- (setYears(magpie[, 2015, ], NULL) - setYears(magpie[, 1995, ], NULL)) / setYears(magpie[, 1995, ], NULL)
-    shrLostHist <- (setYears(data[getRegions(magpie), 2015, ], NULL) - setYears(data[getRegions(magpie), 1995, ], NULL)) / setYears(data[getRegions(magpie), 1995, ], NULL)
-    out <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)
-    out[ , seq(2000, 2015, by = 5), ] <- shrLostHist / shrLostProj
+    shrLostProj <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)
+    shrLostHist <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)
+    for (i in 2:length(y)) {
+      shrLostProj[ , y[i], ] <- (setYears(magpie[, y[i], ], NULL) - setYears(magpie[, y[i-1], ], NULL)) / setYears(magpie[, y[i-1], ], NULL)
+      shrLostHist[ , y[i], ] <- (setYears(data[, y[i], ], NULL) - setYears(data[, y[i-1], ], NULL)) / setYears(data[, y[i-1], ], NULL)
+    }
+    out <- shrLostHist / shrLostProj
     out[is.na(out)] <- 0
     out[is.infinite(out)] <- 0
     getNames(out) <- NULL
     
     # only reward if share of cropland lost between 1995 and 2015 exceeds a certain threshold. Otherwise set to 0.
-    out[rownames(which(shrLostHist > -calib_accuracy, arr.ind = T)),,] <- 0
-    out[which(out < 0, arr.ind = T)] <- 0
+    out[shrLostHist > -calib_accuracy] <- 0
+    out[out < 0] <- 0
   }
+  out <- lowpass(out,i = lowpass_filter)
   return(magpiesort(out))
 }
 
@@ -147,13 +150,13 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
     
     # set value for India to cost_max for better convergence
     if ("IND" %in% getRegions(calib_factor_cost)) {
-      calib_factor_cost["IND", , ] <- cost_max
+      calib_factor_cost["IND", -1 , ] <- cost_max
     }
     
     # set value for CAZ to cost_max for better convergence; only needed for LUH2 due to mismatch between LUH2 and FAO historical data for cropland.
     if (histData == "MAgPIEown") {
       if ("CAZ" %in% getRegions(calib_factor_cost)) {
-        calib_factor_cost["CAZ", , ] <- cost_max
+        calib_factor_cost["CAZ", -1 , ] <- cost_max
       }
     }
     

From 56da0441ec5d4dc4b700602ae84aaffb7936d14e Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Sun, 2 Feb 2025 23:16:35 +0100
Subject: [PATCH 05/14] update

---
 scripts/calibration/landconversion_cost.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/scripts/calibration/landconversion_cost.R b/scripts/calibration/landconversion_cost.R
index 44ac1fcb8..50ed7ecee 100644
--- a/scripts/calibration/landconversion_cost.R
+++ b/scripts/calibration/landconversion_cost.R
@@ -35,7 +35,7 @@ calibration_run <- function(putfolder, calib_magpie_name, logoption = 3, s_use_g
 
 # get ratio between modelled area and reference area
 
-getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpass_filter = 1, histData = "FAO") {
+getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpass_filter = 1, histData = "FAO", cost_min = 0.2) {
   require(magclass)
   require(magpie4)
   require(gdx2)
@@ -65,7 +65,7 @@ getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpa
     out[is.na(out)] <- 1
     getNames(out) <- NULL
     out[out < 0] <- 1
-    out[,1,] <- 0
+    out[,1,] <- cost_min
   } else if (mode == "reward") {
     shrLostProj <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)
     shrLostHist <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)

From 9f3f32a7bcea1dcbed46358b1884b2738193eccf Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Mon, 3 Feb 2025 09:19:55 +0100
Subject: [PATCH 06/14] update

---
 scripts/calibration/landconversion_cost.R | 27 ++++++++++++++---------
 1 file changed, 16 insertions(+), 11 deletions(-)

diff --git a/scripts/calibration/landconversion_cost.R b/scripts/calibration/landconversion_cost.R
index 50ed7ecee..e0820b014 100644
--- a/scripts/calibration/landconversion_cost.R
+++ b/scripts/calibration/landconversion_cost.R
@@ -40,6 +40,7 @@ getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpa
   require(magpie4)
   require(gdx2)
   y <- readGDX(gdx_file,"t")
+  y <- seq(1995,2015,by=5)
   magpie <- land(gdx_file)[, y, "crop"]
   if (histData == "MAgPIEown") {
     hist <- dimSums(readGDX(gdx_file, "f10_land")[, , "crop"], dim = 1.2)
@@ -65,7 +66,7 @@ getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpa
     out[is.na(out)] <- 1
     getNames(out) <- NULL
     out[out < 0] <- 1
-    out[,1,] <- cost_min
+    out <- lowpass(out,i = lowpass_filter)
   } else if (mode == "reward") {
     shrLostProj <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)
     shrLostHist <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)
@@ -73,16 +74,18 @@ getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpa
       shrLostProj[ , y[i], ] <- (setYears(magpie[, y[i], ], NULL) - setYears(magpie[, y[i-1], ], NULL)) / setYears(magpie[, y[i-1], ], NULL)
       shrLostHist[ , y[i], ] <- (setYears(data[, y[i], ], NULL) - setYears(data[, y[i-1], ], NULL)) / setYears(data[, y[i-1], ], NULL)
     }
-    out <- shrLostHist / shrLostProj
+    
+    out <- magpie / data - 1
     out[is.na(out)] <- 0
     out[is.infinite(out)] <- 0
     getNames(out) <- NULL
     
-    # only reward if share of cropland lost between 1995 and 2015 exceeds a certain threshold. Otherwise set to 0.
-    out[shrLostHist > -calib_accuracy] <- 0
+    # set reward to 0 if no cropland was lost in historic data set
+    out[shrLostHist < -calib_accuracy] <- 0
     out[out < 0] <- 0
+    out <- lowpass(out,i = lowpass_filter)
+    out[,1,] <- 0
   }
-  out <- lowpass(out,i = lowpass_filter)
   return(magpiesort(out))
 }
 
@@ -116,7 +119,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
   calib_divergence_cost <- abs(calib_correction_cost - 1)
 
   calib_correction_reward <- getCalibFactor(gdx_file, mode = "reward", calib_accuracy = calib_accuracy, lowpass_filter = lowpass_filter)
-  calib_divergence_reward <- abs(calib_correction_reward - 1)
+  calib_divergence_reward <- abs(calib_correction_reward)
   calib_divergence_reward[calib_correction_reward == 0] <- 0
   
   ### -> in case it is the first step, it forces the initial factors to be equal to 1
@@ -125,12 +128,14 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
     start_flag <- FALSE
   } else {
     old_calib <- new.magpie(cells_and_regions = getCells(calib_divergence_cost), years = y, names = c("cost", "reward"), fill = 1)
+    old_calib[,,"reward"] <- 0
     start_flag <- TRUE
   }
 
   calib_factor_cost <- setNames(old_calib[, , "cost"], NULL) * calib_correction_cost
-  calib_factor_reward <- setNames(old_calib[, , "reward"], NULL) * calib_correction_reward
-
+  calib_factor_reward <- setNames(old_calib[, , "reward"], NULL) + calib_correction_reward
+  calib_factor_reward[calib_factor_reward < 0] <- 0
+  
   if (!start_flag) {
     # use calibration factors where accuracy was reached
     # use stricter divergence threshold in first 5 calibration_step steps
@@ -164,9 +169,9 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
     calib_factor_cost[above_limit] <- cost_max
     calib_divergence_cost[getRegions(calib_factor_cost), , ][above_limit] <- 0
     
-    above_limit <- (calib_factor_reward >= cost_max)
-    calib_factor_reward[above_limit] <- cost_max
-    calib_divergence_reward[getRegions(calib_factor_reward), , ][above_limit] <- 0
+    # above_limit <- (calib_factor_reward >= cost_max)
+    # calib_factor_reward[above_limit] <- cost_max
+    # calib_divergence_reward[getRegions(calib_factor_reward), , ][above_limit] <- 0
   }
 
   if (!is.null(cost_min)) {

From 3f139119fa7216fdd0a4d4014569ea8e0c7b94eb Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Mon, 3 Feb 2025 10:59:15 +0100
Subject: [PATCH 07/14] test

---
 scripts/calibration/landconversion_cost.R | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/scripts/calibration/landconversion_cost.R b/scripts/calibration/landconversion_cost.R
index e0820b014..99c3a063a 100644
--- a/scripts/calibration/landconversion_cost.R
+++ b/scripts/calibration/landconversion_cost.R
@@ -150,8 +150,8 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
 
   if (!is.null(cost_max)) {
     # if reward exists, set cost calibration to cost_max
-    reward_exists <- (calib_factor_reward > 0)
-    calib_factor_cost[reward_exists] <- cost_max
+    # reward_exists <- (calib_factor_reward > 0)
+    # calib_factor_cost[reward_exists] <- cost_max
     
     # set value for India to cost_max for better convergence
     if ("IND" %in% getRegions(calib_factor_cost)) {

From 2350ddbcb607b08e0480002e1b728b5a60866ba6 Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Mon, 3 Feb 2025 14:23:34 +0100
Subject: [PATCH 08/14] update default values

---
 config/default.cfg                        |  2 +-
 scripts/calibration/landconversion_cost.R | 13 +++----------
 2 files changed, 4 insertions(+), 11 deletions(-)

diff --git a/config/default.cfg b/config/default.cfg
index 175c1806d..5d9542f13 100644
--- a/config/default.cfg
+++ b/config/default.cfg
@@ -102,7 +102,7 @@ cfg$restart_landconversion_cost <- FALSE           # def = FALSE
 # for time steps 1995-2015
 cfg$lowpass_filter_landconversion_cost <- 1        # def= 1
 # Set upper limit for cropland calibration factor
-cfg$cost_calib_max_landconversion_cost <- 3         # def= 3
+cfg$cost_calib_max_landconversion_cost <- 2.5         # def= 2.5
 # Set lower limit for cropland calibration factor
 cfg$cost_calib_min_landconversion_cost <- 0.2            # def= 0.2
 # Selection type of calibration factors.
diff --git a/scripts/calibration/landconversion_cost.R b/scripts/calibration/landconversion_cost.R
index 99c3a063a..4ba2cf1e3 100644
--- a/scripts/calibration/landconversion_cost.R
+++ b/scripts/calibration/landconversion_cost.R
@@ -40,7 +40,6 @@ getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpa
   require(magpie4)
   require(gdx2)
   y <- readGDX(gdx_file,"t")
-  y <- seq(1995,2015,by=5)
   magpie <- land(gdx_file)[, y, "crop"]
   if (histData == "MAgPIEown") {
     hist <- dimSums(readGDX(gdx_file, "f10_land")[, , "crop"], dim = 1.2)
@@ -107,7 +106,7 @@ time_series_reward <- function(calib_factor) {
 }
 
 # Calculate the correction factor and save it
-update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, calib_file, cost_max = 3, cost_min = 0.2, calibration_step = "", n_maxcalib = 20, best_calib = FALSE, histData = "FAO") {
+update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, calib_file, cost_max = 2.5, cost_min = 0.2, calibration_step = "", n_maxcalib = 20, best_calib = FALSE, histData = "FAO") {
   require(magclass)
   require(magpie4)
   require(gdx2)
@@ -149,10 +148,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
   }
 
   if (!is.null(cost_max)) {
-    # if reward exists, set cost calibration to cost_max
-    # reward_exists <- (calib_factor_reward > 0)
-    # calib_factor_cost[reward_exists] <- cost_max
-    
+
     # set value for India to cost_max for better convergence
     if ("IND" %in% getRegions(calib_factor_cost)) {
       calib_factor_cost["IND", -1 , ] <- cost_max
@@ -169,9 +165,6 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
     calib_factor_cost[above_limit] <- cost_max
     calib_divergence_cost[getRegions(calib_factor_cost), , ][above_limit] <- 0
     
-    # above_limit <- (calib_factor_reward >= cost_max)
-    # calib_factor_reward[above_limit] <- cost_max
-    # calib_divergence_reward[getRegions(calib_factor_reward), , ][above_limit] <- 0
   }
 
   if (!is.null(cost_min)) {
@@ -268,7 +261,7 @@ update_calib <- function(gdx_file, calib_accuracy = 0.01, lowpass_filter = 1, ca
 calibrate_magpie <- function(n_maxcalib = 20,
                              restart = FALSE,
                              calib_accuracy = 0.01,
-                             cost_max = 3,
+                             cost_max = 2.5,
                              cost_min = 0.2,
                              calib_magpie_name = "magpie_calib",
                              lowpass_filter = 1,

From 58b3f50e86397b0f2f58478c97d93e1f7af980aa Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Mon, 3 Feb 2025 14:26:19 +0100
Subject: [PATCH 09/14] cleanup code

---
 scripts/calibration/landconversion_cost.R | 2 --
 1 file changed, 2 deletions(-)

diff --git a/scripts/calibration/landconversion_cost.R b/scripts/calibration/landconversion_cost.R
index 4ba2cf1e3..bf3356d86 100644
--- a/scripts/calibration/landconversion_cost.R
+++ b/scripts/calibration/landconversion_cost.R
@@ -67,10 +67,8 @@ getCalibFactor <- function(gdx_file, mode = "cost", calib_accuracy = 0.05, lowpa
     out[out < 0] <- 1
     out <- lowpass(out,i = lowpass_filter)
   } else if (mode == "reward") {
-    shrLostProj <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)
     shrLostHist <- new.magpie(getRegions(magpie), getYears(magpie), fill = 0)
     for (i in 2:length(y)) {
-      shrLostProj[ , y[i], ] <- (setYears(magpie[, y[i], ], NULL) - setYears(magpie[, y[i-1], ], NULL)) / setYears(magpie[, y[i-1], ], NULL)
       shrLostHist[ , y[i], ] <- (setYears(data[, y[i], ], NULL) - setYears(data[, y[i-1], ], NULL)) / setYears(data[, y[i-1], ], NULL)
     }
     

From db50da04170732d7c7f41c7ef50aee880cc10b7e Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Mon, 3 Feb 2025 14:28:36 +0100
Subject: [PATCH 10/14] changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index a854989ca..e08e6d563 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,6 +12,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - **config.cfg** default for `cfg$gms$s29_treecover_max` changed from 0.4 to 1
 - **config.cfg** default for `cfg$gms$s29_fallow_max ` changed from 0.4 to 0
 - **config.cfg** default for `cfg$gms$s35_forest_damage ` changed from 2 to 0
+- **scripts** land conversion cost calibration for cropland - FAO as target data set instead of MAgPIEown
 
 ### added
 - **scenario_config.csv** added column `NPI-revert`

From c1cd8b186fa5af92f0ede50798aba6ab936be6d5 Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Mon, 3 Feb 2025 20:27:30 +0100
Subject: [PATCH 11/14] update start scripts and default.cfg

---
 config/default.cfg                     | 4 ++--
 scripts/start/extra/recalibrateH16.R   | 2 +-
 scripts/start/extra/recalibrate_FSEC.R | 3 +--
 3 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/config/default.cfg b/config/default.cfg
index 5d9542f13..3c255d640 100644
--- a/config/default.cfg
+++ b/config/default.cfg
@@ -26,7 +26,7 @@ cfg$input <- c(regional    = "rev4.116_h12_magpie.tgz",
                cellular    = "rev4.116_h12_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz",
                validation  = "rev4.116_h12_validation.tgz",
                additional  = "additional_data_rev4.60.tgz",
-               calibration = "calibration_H12_27Sep24.tgz")
+               calibration = "calibration_H12_FAO30_03Feb25.tgz")
 
 # NOTE: It is recommended to recalibrate the model when changing cellular input data
 #       as well as for any other setting that would affect initial values in the model,
@@ -100,7 +100,7 @@ cfg$calib_maxiter_landconversion_cost <- 20           # def = 20
 cfg$restart_landconversion_cost <- FALSE           # def = FALSE
 # Number of lowpass filter iterations applied on calibration factors
 # for time steps 1995-2015
-cfg$lowpass_filter_landconversion_cost <- 1        # def= 1
+cfg$lowpass_filter_landconversion_cost <- 0        # def= 0
 # Set upper limit for cropland calibration factor
 cfg$cost_calib_max_landconversion_cost <- 2.5         # def= 2.5
 # Set lower limit for cropland calibration factor
diff --git a/scripts/start/extra/recalibrateH16.R b/scripts/start/extra/recalibrateH16.R
index 4cc1c33ff..f85889477 100644
--- a/scripts/start/extra/recalibrateH16.R
+++ b/scripts/start/extra/recalibrateH16.R
@@ -30,4 +30,4 @@ cfg$output <- c("rds_report")
 cfg$force_replace <- TRUE
 cfg$qos <- "priority"
 start_run(cfg,codeCheck=FALSE)
-magpie4::submitCalibration("H16")
+magpie4::submitCalibration("H16_FAO")
diff --git a/scripts/start/extra/recalibrate_FSEC.R b/scripts/start/extra/recalibrate_FSEC.R
index f651a5f03..b4fd187ed 100644
--- a/scripts/start/extra/recalibrate_FSEC.R
+++ b/scripts/start/extra/recalibrate_FSEC.R
@@ -25,9 +25,8 @@ cfg$results_folder                  <- "output/:title:"
 cfg$recalibrate                     <- TRUE # required when penality_apr22 activated
 cfg$best_calib                      <- TRUE
 cfg$recalibrate_landconversion_cost <- TRUE
-cfg$best_calib_landconversion_cost  <- TRUE
 cfg$output                          <- c("rds_report")
 cfg$force_replace                   <- TRUE
 cfg$qos <- "priority"
 start_run(cfg, codeCheck = FALSE)
-magpie4::submitCalibration("FSEC")
+magpie4::submitCalibration("FSEC_FAO")

From da55637c3947a58741fd073099edf7b1a5f68ad8 Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Mon, 3 Feb 2025 20:33:54 +0100
Subject: [PATCH 12/14] changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index e08e6d563..31e584de0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - **config.cfg** default for `cfg$gms$s29_fallow_max ` changed from 0.4 to 0
 - **config.cfg** default for `cfg$gms$s35_forest_damage ` changed from 2 to 0
 - **scripts** land conversion cost calibration for cropland - FAO as target data set instead of MAgPIEown
+- **default.cfg** settings for  land conversion cost calibration updated
 
 ### added
 - **scenario_config.csv** added column `NPI-revert`

From 96cead0509000cdafc6973e06e737ae711c2f921 Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Tue, 4 Feb 2025 08:35:58 +0100
Subject: [PATCH 13/14] update calib factor FSEC

---
 config/projects/scenario_config_fsec.csv | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/config/projects/scenario_config_fsec.csv b/config/projects/scenario_config_fsec.csv
index c830f7fed..636aa5fb0 100644
--- a/config/projects/scenario_config_fsec.csv
+++ b/config/projects/scenario_config_fsec.csv
@@ -79,6 +79,5 @@ gms$c73_build_demand;;;;;;;;;;;;;;;;;;;;;;;;50pc;;;;;;;;;;;;;;;;;;;;;;;;;;;
 input['cellular'];rev4.116_FSEC_3c888fa5_cellularmagpie_c200_MRI-ESM2-0-ssp460_lpjml-8e6c5eb1.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;rev4.116_FSEC_0bd54110_cellularmagpie_c200_MRI-ESM2-0-ssp119_lpjml-8e6c5eb1.tgz;rev4.116_FSEC_6819938d_cellularmagpie_c200_MRI-ESM2-0-ssp126_lpjml-8e6c5eb1.tgz;;rev4.116_FSEC_1b5c3817_cellularmagpie_c200_MRI-ESM2-0-ssp245_lpjml-8e6c5eb1.tgz;rev4.116_FSEC_3c888fa5_cellularmagpie_c200_MRI-ESM2-0-ssp460_lpjml-8e6c5eb1.tgz;rev4.116_FSEC_fd712c0b_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1.tgz;rev4.116_FSEC_09a63995_cellularmagpie_c200_MRI-ESM2-0-ssp585_lpjml-8e6c5eb1.tgz;;;
 input['regional'];rev4.116_FSEC_magpie.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 input['validation'];rev4.116_FSEC_validation.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-input['additional'];additional_data_rev4.59.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-input['calibration'];calibration_FSEC_19Nov24.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+input['calibration'];calibration_FSEC_FAO_03Feb25.tgz;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 magicc_emis_scen;bjoernAR6_C_SSP2-NDC.mif;;;bjoernAR6_C_SSP2-PkBudg900.mif;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;bjoernAR6_C_SSP1-NDC.mif;;;;;;;;;;;;bjoernAR6_C_RemSDP-900-MagSSP1.mif;;

From 0a36a6ddc98419bc96c6a8c19836d2ef8cdfb4fe Mon Sep 17 00:00:00 2001
From: florianh <humpenoeder@pik-potsdam.de>
Date: Tue, 4 Feb 2025 09:30:32 +0100
Subject: [PATCH 14/14] callib factors H16

---
 scripts/start/extra/recalibrateH16.R   | 2 +-
 scripts/start/projects/project_ABCDR.R | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/scripts/start/extra/recalibrateH16.R b/scripts/start/extra/recalibrateH16.R
index f85889477..b4c1071b6 100644
--- a/scripts/start/extra/recalibrateH16.R
+++ b/scripts/start/extra/recalibrateH16.R
@@ -19,7 +19,7 @@ source("scripts/start_functions.R")
 source("config/default.cfg")
 cfg$input['regional'] <- "rev4.116_36f73207_magpie.tgz"
 cfg$input['validation'] <- "rev4.116_36f73207_validation.tgz"
-cfg$input['calibration'] <- "calibration_H16_27Sep24.tgz"
+cfg$input['calibration'] <- "calibration_H16_FAO_03Feb25.tgz"
 cfg$input['cellular'] <- "rev4.116_36f73207_44a213b6_cellularmagpie_c400_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1_clusterweight-ba4466a8.tgz"
 
 cfg$results_folder <- "output/:title:"
diff --git a/scripts/start/projects/project_ABCDR.R b/scripts/start/projects/project_ABCDR.R
index c37eb4194..dfbd722ca 100644
--- a/scripts/start/projects/project_ABCDR.R
+++ b/scripts/start/projects/project_ABCDR.R
@@ -33,7 +33,7 @@ cfg$qos <- "standby_dayMax"
 
 cfg$input['regional'] <- "rev4.116_36f73207_magpie.tgz"
 cfg$input['validation'] <- "rev4.116_36f73207_validation.tgz"
-cfg$input['calibration'] <- "calibration_H16_27Sep24.tgz"
+cfg$input['calibration'] <- "calibration_H16_FAO_03Feb25.tgz"
 cfg$input['cellular'] <- "rev4.116_36f73207_bd86374e_cellularmagpie_c200_MRI-ESM2-0-ssp370_lpjml-8e6c5eb1_clusterweight-ba4466a8.tgz"
 
 ssp <- "SSP2"