diff --git a/debug_config.properties b/debug_config.properties
index 717dffaff034bd4bac6453e69eb9d95c661faabf..6f97253c832cd0d8951ab8fb1e968b6b2fcf3d2f 100644
--- a/debug_config.properties
+++ b/debug_config.properties
@@ -1,30 +1,29 @@
BASE_DIR=..
-YIELD_DIR=/Users/peteralexander/Documents/LURG/LPJ/LPJGPLUM_remap6p7_20190225/rcp26
-YIELD_FILENAME=yield.out
+YIELD_DIR=/Users/palexand/Documents/LURG/LPJ/CMIP_emulated_test_forPLUM_20200821180759/GFDL-ESM4_ssp126_v20200821_rmol/sim_LPJ-GUESS
-DEBUG_LIMIT_COUNTRIES=false
-DEBUG_COUNTRY_NAME=Mongolia
+END_TIMESTEP=0
-IS_CALIBRATION_RUN = true
-GENERATE_NEW_YIELD_CLUSTERS = false
-NUM_YIELD_CLUSTERS=8000
+END_FIRST_STAGE_CALIBRATION=10
+TIMESTEP_SIZE=1
-END_TIMESTEP=0
-TIMESTEP_SIZE=10
+IS_CALIBRATION_RUN=true
-INTERPOLATE_OUTPUT_YEARS = false
+GENERATE_NEW_YIELD_CLUSTERS=false
-CHANGE_YIELD_DATA_YEAR=false
+DEMAND_PRICE_IMPORT_AND_PROD_COST=true
-ORIG_LEAST_COST_MIN=true
+YIELD_FILENAME=yield_intpinfs_rmol.out
+IRRIG_MAX_WATER_FILENAME=gsirrigation_intpinfs_rmol.out
-DEBUG_JUST_DEMAND_OUTPUT=false
-PRICE_ELASTIC_DEMAND=true
+MIN_FERT_AMOUNT=10
+MID_FERT_AMOUNT=60
+MAX_FERT_AMOUNT=200
+FERT_AMOUNT_PADDING=3
-GAMS_COUNTRY_TO_SAVE=Brazil
+LPJG_TIMESTEP_SIZE=10
-SSP_SCENARIO=SSP2_v9_130325
-DONT_REBASE_DEMAND=false
+EXTRAPOLATE_YIELD_FERT_RESPONSE=false
-DEMAND_RECALC_MAX_ITERATIONS=0
\ No newline at end of file
+DEBUG_LIMIT_COUNTRIES=true
+DEBUG_COUNTRY_NAME=Italy
\ No newline at end of file
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 1e7a00fe0665a84f3c4f317b8d825ec13fb139b6..ab610f85d8818fc9b931a14039b5cd4290565f9d 100755
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -155,7 +155,7 @@ public class ModelConfig {
public static final boolean DONT_REBASE_DEMAND = getBooleanProperty("DONT_REBASE_DEMAND", false);;
public static final String DEMAND_CURVES_FILE = getProperty("DEMAND_CURVES_FILE", DATA_DIR + File.separator + "com_curves.csv"); // either DEMAND_CURVES_FILE or DEMAND_CONSUMPTION_FILE is used, but not both
public static final String DEMAND_CONSUMPTION_FILE = getProperty("DEMAND_CONSUMPTION_FILE", DATA_DIR + File.separator + "hist_comsump.csv");
- public static final boolean DEMAND_PRICE_IMPORT_AND_PROD_COST = getBooleanProperty("DEMAND_PRICE_IMPORT_AND_PROD_COST", false);
+ public static final boolean DEMAND_PRICE_IMPORT_AND_PROD_COST = getBooleanProperty("DEMAND_PRICE_IMPORT_AND_PROD_COST", true);
public static final String SSP_FILE = DATA_DIR + File.separator + "ssp.csv";
public static final String BASELINE_CONSUMP_FILE = DATA_DIR + File.separator + "base_consump.csv";
public static final String CALORIE_PER_T_FILE = DATA_DIR + File.separator + "calories_per_t.csv";
@@ -304,7 +304,6 @@ public class ModelConfig {
public static final double MAX_INCOME_PROPORTION_FOOD_SPEND = getDoubleProperty("MAX_INCOME_PROPORTION_FOOD_SPEND", 0.7);
public static final double PASTURE_HARVEST_FRACTION = getDoubleProperty("PASTURE_HARVEST_FRACTION", 0.5);
- public static final double ANIMAL_FEED_FROM_OTHER_SOURCES_RATE = getDoubleProperty("ANIMAL_FEED_FROM_OTHER_SOURCES_RATE", 0.127); // animal nutrition coming from sources other crops modelled
public static final double MEAT_EFFICIENCY = getDoubleProperty("MEAT_EFFICIENCY", 1.0); // 'meat' is includes feed conversion ratio already, this is tech. change or similar
public static final double IRRIGIATION_EFFICIENCY = getDoubleProperty("IRRIGIATION_EFFICIENCY", 0.5);
public static final int ELLIOTT_BASEYEAR = 2010;
@@ -408,4 +407,6 @@ public class ModelConfig {
public static final boolean USE_CRAFTY_COUNTRIES = getBooleanProperty("USE_CRAFTY_COUNTRIES", false);
public static final String CRAFTY_PRODUCTION_DIR = getProperty("CRAFTY_PRODUCTION_DIR", OUTPUT_DIR + File.separator + "crafty");
+
+ public static final boolean EXTRAPOLATE_YIELD_FERT_RESPONSE = getBooleanProperty("EXTRAPOLATE_YIELD_FERT_RESPONSE", false);
}
diff --git a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
index 5f7d2f0dd3e92bc12fc3fcede25aa59f29dcf91a..af4e20c6a75a76683bf519b623c553fb65697916 100644
--- a/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsLocationOptimiser.java
@@ -31,6 +31,7 @@ import ac.ed.lurg.types.CommodityType;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.Parameter;
+import ac.ed.lurg.types.YieldType;
import ac.ed.lurg.utils.LazyHashMap;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.yield.YieldResponsesItem;
@@ -227,13 +228,17 @@ public class GamsLocationOptimiser {
double maxIrrig = irrigationItem.getMaxIrrigAmount(crop);
if (DEBUG) LogWriter.println(String.format("%d %15s,\t %.1f,\t %.1f, \t %.1f,\t %.1f,\t %.2f,\t\t [%.2f],\t [%.2f],\t {%.2f}",
- locationId, crop.getGamsName(), yresp.getYieldNone(crop), yresp.getYieldFertOnly(crop), yresp.getYieldIrrigOnly(crop), yresp.getYieldMax(crop), yresp.getShockRate(crop),
- yresp.getFertParam(crop), yresp.getIrrigParam(crop), maxIrrig));
-
- setGamsParamValue(yNoneP.addRecord(v), yresp.getYieldNone(crop), 4);
- setGamsParamValue(y_fert.addRecord(v), yresp.getYieldFertOnly(crop), 4);
- setGamsParamValue(y_irrig.addRecord(v), yresp.getYieldIrrigOnly(crop), 4);
- setGamsParamValue(y_both.addRecord(v), yresp.getYieldMax(crop), 4);
+ locationId, crop.getGamsName(),
+ yresp.getExtrapolatedYield(YieldType.NO_FERT_NO_IRRIG, crop),
+ yresp.getExtrapolatedYield(YieldType.FERT_MAX_NO_IRRIG, crop),
+ yresp.getExtrapolatedYield(YieldType.NO_FERT_IRRIG_MAX, crop),
+ yresp.getExtrapolatedYield(YieldType.FERT_MAX_IRRIG_MAX, crop),
+ yresp.getShockRate(crop), yresp.getFertParam(crop), yresp.getIrrigParam(crop), maxIrrig));
+
+ setGamsParamValue(yNoneP.addRecord(v), yresp.getExtrapolatedYield(YieldType.NO_FERT_NO_IRRIG, crop), 4);
+ setGamsParamValue(y_fert.addRecord(v), yresp.getExtrapolatedYield(YieldType.FERT_MAX_NO_IRRIG, crop), 4);
+ setGamsParamValue(y_irrig.addRecord(v), yresp.getExtrapolatedYield(YieldType.NO_FERT_IRRIG_MAX, crop), 4);
+ setGamsParamValue(y_both.addRecord(v), yresp.getExtrapolatedYield(YieldType.FERT_MAX_IRRIG_MAX, crop), 4);
setGamsParamValue(y_shock.addRecord(v), yresp.getShockRate(crop), 4);
setGamsParamValue(fert_p.addRecord(v), yresp.getFertParam(crop), 4);
setGamsParamValue(irrig_p.addRecord(v), yresp.getIrrigParam(crop), 4);
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index b6602c587c6356b96f3672708019d982c47bd8a2..75ff7437b552ea2461cb28b7485c0842684bc48e 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -290,8 +290,8 @@ public class GamsRasterOptimiser {
RasterKey key = entry.getKey();
CropType crop = CropType.WHEAT;
- if (yresp != null && (yresp.getYieldMax(crop) < yresp.getYieldFertOnly(crop) || yresp.getYieldMax(crop) < yresp.getYieldIrrigOnly(crop))) {
- logWarningWithCoord("Inconsistency F only:" + yresp.getYieldFertOnly(crop) + ", I only" + yresp.getYieldIrrigOnly(crop) + ", max " + yresp.getYieldMax(crop) + " at ", key, yieldRaster);
+ if (yresp != null && (yresp.getYield(YieldType.FERT_MAX_IRRIG_MAX, crop) < yresp.getYield(YieldType.FERT_MAX_NO_IRRIG, crop) || yresp.getYield(YieldType.FERT_MAX_IRRIG_MAX, crop) < yresp.getYield(YieldType.NO_FERT_IRRIG_MAX, crop))) {
+ logWarningWithCoord("Inconsistency F only:" + yresp.getYield(YieldType.FERT_MAX_NO_IRRIG, crop) + ", I only" + yresp.getYield(YieldType.NO_FERT_IRRIG_MAX, crop) + ", max " + yresp.getYield(YieldType.FERT_MAX_IRRIG_MAX, crop) + " at ", key, yieldRaster);
}
}
diff --git a/src/ac/ed/lurg/yield/YieldClusterPoint.java b/src/ac/ed/lurg/yield/YieldClusterPoint.java
index f8f5cfc3f07cfa0a30945401b82c0a2b3545362d..91f9eb806ce9ddf2cc0013026b19b985a130f115 100644
--- a/src/ac/ed/lurg/yield/YieldClusterPoint.java
+++ b/src/ac/ed/lurg/yield/YieldClusterPoint.java
@@ -5,6 +5,7 @@ import java.util.Collection;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.types.YieldType;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.utils.cluster.ClusteringPoint;
import ac.sac.raster.RasterKey;
@@ -34,12 +35,12 @@ public class YieldClusterPoint implements ClusteringPoint<String> {
this.rasterKey = rasterKey;
// not sure if we be better to get a reference to the YieldResponsesItem, rather than caching these values?
- this.wheatMin = yields.getYieldNone(CropType.WHEAT);
- this.riceMax = yields.getYieldMax(CropType.RICE);
- this.maizeMin = yields.getYieldNone(CropType.MAIZE);
- this.pasture = yields.getYieldNone(CropType.PASTURE);
- this.wheatMax = yields.getYieldMax(CropType.WHEAT);
- this.maizeMax = yields.getYieldMax(CropType.MAIZE);
+ this.wheatMin = yields.getYield(YieldType.NO_FERT_NO_IRRIG, CropType.WHEAT);
+ this.riceMax = yields.getYield(YieldType.FERT_MAX_IRRIG_MAX, CropType.RICE);
+ this.maizeMin = yields.getYield(YieldType.NO_FERT_NO_IRRIG, CropType.MAIZE);
+ this.pasture = yields.getYield(YieldType.NO_FERT_NO_IRRIG, CropType.PASTURE);
+ this.wheatMax = yields.getYield(YieldType.FERT_MAX_IRRIG_MAX, CropType.WHEAT);
+ this.maizeMax = yields.getYield(YieldType.FERT_MAX_IRRIG_MAX, CropType.MAIZE);
this.irrigWaterAval = irrigItem.getIrrigConstraint();
this.protectedArea = protectedArea;
}
diff --git a/src/ac/ed/lurg/yield/YieldResponse.java b/src/ac/ed/lurg/yield/YieldResponse.java
index 86e1e9b2a66110c20f554b4cfd9a50422a2d7230..fb4b7c4a58a68f3f3a574c00040e3c1323fee796 100644
--- a/src/ac/ed/lurg/yield/YieldResponse.java
+++ b/src/ac/ed/lurg/yield/YieldResponse.java
@@ -5,9 +5,11 @@ import java.util.Map;
import ac.ed.lurg.ModelConfig;
import ac.ed.lurg.types.YieldType;
+import ac.ed.lurg.utils.LogWriter;
public class YieldResponse {
private Map<YieldType, Double> yields = new HashMap<YieldType, Double>();
+ private Map<YieldType, Double> yieldsExtrapolated;
private double yieldShock;
private double fertParm;
private double irrigParm;
@@ -15,7 +17,7 @@ public class YieldResponse {
private static double defaultParam;
static {
- defaultParam = calcParam(0, 0.8, 1, 0, 0.5, 1.0); // default assumes 80% at mid point
+ defaultParam = calcParam(0, 0.8, 1, 0.5); // default assumes 80% at mid point
}
@@ -36,14 +38,59 @@ public class YieldResponse {
return d.doubleValue();
}
+ public double getExtrapolatedYield(YieldType type) {
+ if (yieldsExtrapolated == null)
+ calcParams();
+
+ Double d = yieldsExtrapolated.get(type);
+
+ if (d == null)
+ return Double.NaN;
+
+ return d.doubleValue();
+ }
+
public double getFertParam() {
- if (fertParm == 0) {
- // fertParm = calcParam (0, 0.5, 1, 5.0/200, 50.0/200, 200.0/200); // we do have MID fert data, but looks wrong
- fertParm = calcParam(getYield(YieldType.NO_FERT_NO_IRRIG), getYield(YieldType.FERT_MID_NO_IRRIG), getYield(YieldType.FERT_MAX_NO_IRRIG),
- ModelConfig.MIN_FERT_AMOUNT/ModelConfig.MAX_FERT_AMOUNT, ModelConfig.MID_FERT_AMOUNT/ModelConfig.MAX_FERT_AMOUNT, ModelConfig.MAX_FERT_AMOUNT/ModelConfig.MAX_FERT_AMOUNT);
+ if (fertParm == 0)
+ calcParams();
+
+ return fertParm;
+ }
+
+ private double calcParams() {
+ double yMinNoI = getYield(YieldType.NO_FERT_NO_IRRIG);
+ double yMidNoI = getYield(YieldType.FERT_MID_NO_IRRIG);
+ double yMaxNoI = getYield(YieldType.FERT_MAX_NO_IRRIG);
+
+ double fMid = (ModelConfig.MID_FERT_AMOUNT-ModelConfig.MIN_FERT_AMOUNT)/(ModelConfig.MAX_FERT_AMOUNT-ModelConfig.MIN_FERT_AMOUNT);
+
+ if (ModelConfig.EXTRAPOLATE_YIELD_FERT_RESPONSE && !(yMidNoI <= yMinNoI || yMaxNoI <= yMidNoI)) {
+ double yMinI = getYield(YieldType.NO_FERT_IRRIG_MAX);
+ double yMidI = getYield(YieldType.FERT_MID_IRRIG_MAX);
+ double yMaxI = getYield(YieldType.FERT_MAX_IRRIG_MAX);
+
+ double asymptoteYieldIncNoI = yieldAsymptoteSolve(yMidNoI-yMinNoI, yMaxNoI-yMinNoI, fMid);
+ fertParm = -1/fMid * Math.log(1-(yMidI-yMinI)/asymptoteYieldIncNoI);
+ double asymptoteYieldIncI = (yMaxI-yMinI) / (1 - Math.exp(-fertParm));
+ //LogWriter.println(String.format("fertParm from solver: asymptoteYieldIncNoI: %.2f, asymptoteYieldIncI: %.2f, fertParm %.4f", asymptoteYieldIncNoI, asymptoteYieldIncI, fertParm));
- //LogWriter.println(String.format("%s %s %s", getYield(YieldType.NO_FERT_NO_IRRIG), getYield(YieldType.FERT_MID_NO_IRRIG), getYield(YieldType.FERT_MAX_NO_IRRIG)));
+ if (Double.isNaN(asymptoteYieldIncI) || Double.isNaN(fertParm)) {
+ LogWriter.println("Not finding a suitable extrapolating solution. Defaulting to old style");
+ }
+ else {
+ yieldsExtrapolated = new HashMap<YieldType, Double>();
+ yieldsExtrapolated.put(YieldType.NO_FERT_NO_IRRIG, yMinNoI);
+ yieldsExtrapolated.put(YieldType.FERT_MAX_NO_IRRIG, asymptoteYieldIncNoI+yMinI);
+ yieldsExtrapolated.put(YieldType.NO_FERT_IRRIG_MAX, yMinI);
+ yieldsExtrapolated.put(YieldType.FERT_MAX_IRRIG_MAX, asymptoteYieldIncI+yMinI);
+ }
+ }
+
+ if (yieldsExtrapolated == null) { // Default behaviour if EXTRAPOLATE_YIELD_FERT_RESPONSE is false or a problem with the other approach
+ fertParm = calcParam(yMinNoI, yMidNoI, yMaxNoI, fMid);
+ yieldsExtrapolated = yields;
}
+
return fertParm;
}
@@ -55,11 +102,11 @@ public class YieldResponse {
return irrigParm;
}
- private static double calcParam(double yMin, double yMid, double yMax, double xMin, double xMid, double xMax) {
+ private static double calcParam(double yMin, double yMid, double yMax, double xMid) {
if (yMid <= yMin || yMax <= yMid)
return defaultParam;
-
- return -Math.log(1 - ((yMid-yMin)/(yMax-yMin))) * (xMax-xMin) / (xMid-xMin); // 1-exp(-x) function
+
+ return -Math.log(1 - ((yMid-yMin)/(yMax-yMin))) / xMid; // 1-exp(-x) function
// return Math.log((yMid - yMin)/(yMax - yMin)) / Math.log((xMid - xMin)/(xMax - xMin)) / xMax; // Cobb-Douglas
}
@@ -70,4 +117,48 @@ public class YieldResponse {
public double getShock() {
return yieldShock;
}
+
+ private static double f(double x, double yM, double yH, double fM) {
+ return 1-yM/x - Math.pow(1-yH/x, fM);
+ }
+
+ private double yieldAsymptoteSolve (double yM, double yH, double fM) {
+
+ double a = yH;
+ double b = yH * 20;
+ double epsilon = 0.0001;
+
+ if (fM * yH > yM) {
+ LogWriter.printlnError(String.format("yieldAsymptoteSolve: not diminishing yield increases to N: yM=%.3f, yH=%.3f, fM=%.3f", yM, yH, fM));
+ return Double.NaN;
+ }
+
+ double m = Double.NaN;
+ double fm;
+ double fa = f(a, yM, yH, fM);
+ double fb = f(b, yM, yH, fM);
+
+ if (fa * fb > 0) {
+ LogWriter.printlnError(String.format("yieldAsymptoteSolve: not bracketing root: yM=%.3f, yH=%.3f, fM=%.3f", yM, yH, fM));
+ return Double.NaN;
+ }
+
+ while ((b-a) > epsilon)
+ {
+ m = (a+b)/2; // Mid point
+ fm = f(m, yM, yH, fM);
+
+ if ( fm * fa > 0 )
+ { // f(a) and f(m) have same signs: move a
+ a = m;
+ fa = fm;
+ }
+ else
+ { // f(a) and f(m) have different signs: move b
+ b = m;
+ fb = fm;
+ }
+ }
+ return m;
+ }
}
\ No newline at end of file
diff --git a/src/ac/ed/lurg/yield/YieldResponsesItem.java b/src/ac/ed/lurg/yield/YieldResponsesItem.java
index e78fa2ff79316084cd243949d99aa0d0a9265720..e6723643f8a33dee46399e90f8eb44eaeef77d1b 100644
--- a/src/ac/ed/lurg/yield/YieldResponsesItem.java
+++ b/src/ac/ed/lurg/yield/YieldResponsesItem.java
@@ -29,26 +29,11 @@ public class YieldResponsesItem implements RasterItem {
return d;
}
- public double getYieldNone(CropType crop) {
- return getYield(YieldType.NO_FERT_NO_IRRIG, crop);
+ public double getExtrapolatedYield(YieldType yieldType, CropType crop) {
+ YieldResponse yr = getYieldResponseForCrop(crop);
+ return yr.getExtrapolatedYield(yieldType);
}
- public double getYieldMidFertOnly(CropType crop) {
- return getYield(YieldType.FERT_MID_NO_IRRIG, crop);
- }
-
- public double getYieldFertOnly(CropType crop) {
- return getYield(YieldType.FERT_MAX_NO_IRRIG, crop);
- }
-
- public double getYieldIrrigOnly(CropType crop) {
- return getYield(YieldType.NO_FERT_IRRIG_MAX, crop);
- }
-
- public double getYieldMax(CropType crop) {
- return getYield(YieldType.FERT_MAX_IRRIG_MAX, crop);
- }
-
public double getFertParam(CropType crop) {
return getYieldResponseForCrop(crop).getFertParam();
}