merged juliette/documentation with master

src/ac/ed/lurg/country/CountryAgent.java

@@ -38,7 +38,6 @@ public class CountryAgent extends AbstractCountryAgent {
 
 	private GamsRasterOutput previousGamsRasterOutput;
 	private RasterSet<IntegerRasterItem> yieldClusters;
-	private double animalFeedFromOtherSources;
 	private Map<CropType, Double> subsidyRates;
 	private boolean saveGamsGdxFiles;
 
@@ -127,11 +126,6 @@ public class CountryAgent extends AbstractCountryAgent {
 		if (saveGamsGdxFiles && ModelConfig.PRICE_ELASTIC_DEMAND)
 			saveGDXFile("demand");
 		
-		if (currentTimestep.isInitialTimestep()) {
-			double totalAnimalProductDemand = currentProjectedDemand.get(CommodityType.MONOGASTRICS) + currentProjectedDemand.get(CommodityType.RUMINANTS);
-			animalFeedFromOtherSources = totalAnimalProductDemand * ModelConfig.ANIMAL_FEED_FROM_OTHER_SOURCES_RATE;
-		}
-						
 		if (currentProjectedDemand.size() == 0) {
 			LogWriter.printlnError("No demand for country " + country + " so skipping it");
 		}
@@ -221,7 +215,7 @@ public class CountryAgent extends AbstractCountryAgent {
 		}
 
 		GamsCountryInput countryLevelInputs = new GamsCountryInput(country, currentProjectedDemand, currentCountryPrices, importConstraints, 
-				previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, animalFeedFromOtherSources, subsidyRates);	
+				previousGamsRasterOutput.getCropUsageData(), currentMinDemandFract, subsidyRates);	
 		GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData, countryLevelInputs);
 
 		return input;

src/ac/ed/lurg/country/GlobalPrice.java

@@ -73,7 +73,12 @@ public class GlobalPrice implements Serializable {
 			double transportLossRate = (!ModelConfig.SHOCKS_POSSIBLE) ? ModelConfig.TRANSPORT_LOSSES : ModelConfig.getParameter(Parameter.TRANSPORT_LOSSES, timestep.getYear());
 			double newTransportLosses = newExportAmountBeforeLoss * transportLossRate;
 			double stockChange = newExportAmountBeforeLoss - newTransportLosses - newImports - oldDiff;
-			double updatedStock = stockLevel + stockChange;
+			
+			double updatedStock;
+			if (ModelConfig.IS_CALIBRATION_RUN  && timestep.getTimestep() <= ModelConfig.END_FIRST_STAGE_CALIBRATION)
+				updatedStock = stockLevel;  // don't update stock in inital stage of calibration
+			else
+				updatedStock = stockLevel + stockChange;
 				
 			LogWriter.println(String.format("     imports %.2f, exports %.2f", newImports, newExportAmountBeforeLoss - newTransportLosses));
 			LogWriter.println(String.format("     updatedStock %.2f, previous %.2f (last timestep %.2f), stockChange %.2f, oldDiff %.2f", updatedStock, stockLevel, stockLevel-oldDiff, stockChange, oldDiff));

src/ac/ed/lurg/country/gams/GamsCountryInput.java

@@ -19,12 +19,11 @@ public class GamsCountryInput {
 	private Map<CropType, CountryPrice> countryPrices;
 	private Map<CropType, CropUsageData> previousCropUsageData;
 	private Map<CommodityType, Map<CropType, Double>> minDemandFractions;
-	private double animalFeedFromOtherSources;
 	private Map<CropType, Double> subsidyRates;
 
 	public GamsCountryInput(CompositeCountry country, Map<CommodityType, Double> projectedDemand, Map<CropType, CountryPrice> countryPrices, 
 			Map<CropType, TradeConstraint> importConstraints, Map<CropType, CropUsageData> previousCropUsageData, 
-			Map<CommodityType, Map<CropType, Double>> minDemandFracts, double animalFeedFromOtherSources, Map<CropType, Double> subsidyRates) {
+			Map<CommodityType, Map<CropType, Double>> minDemandFracts, Map<CropType, Double> subsidyRates) {
 		super();
 		this.country = country;
 		this.projectedDemand = projectedDemand;
@@ -32,7 +31,6 @@ public class GamsCountryInput {
 		this.countryPrices = countryPrices;
 		this.previousCropUsageData = previousCropUsageData;
 		this.minDemandFractions = minDemandFracts;
-		this.animalFeedFromOtherSources = animalFeedFromOtherSources;
 		this.subsidyRates = subsidyRates;
 	}
 		
@@ -82,10 +80,6 @@ public class GamsCountryInput {
 		return minDemandFractions;
 	}
 	
-	public double getAnimalFeedFromOtherSources() {
-		return animalFeedFromOtherSources;
-	}
-	
 	public Map<CropType, Double> getSubsidyRates() {
 		return subsidyRates;
 	}

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);
@@ -310,7 +315,6 @@ public class GamsLocationOptimiser {
 		addScalar(inDB, "unhandledCropRate", ModelConfig.UNHANDLED_CROP_RATE, 3);
 		addScalar(inDB, "setAsideRate", ModelConfig.SETASIDE_RATE, 5);
 		addScalar(inDB, "domesticPriceMarkup", ModelConfig.DOMESTIC_PRICE_MARKUP, 3);
-		addScalar(inDB, "animalFeedFromOtherSources", countryInput.getAnimalFeedFromOtherSources(), 2);		
 		double maxExpansion = 1.0;
 		if (!ModelConfig.IS_CALIBRATION_RUN && countryInput.getCountry().getName().equals("China")) {
 			maxExpansion = ModelConfig.MAX_CHINA_LAND_EXPANSION_RATE;

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);
 			}
 		}
 		

src/ac/ed/lurg/utils/cluster/Cluster.java

@@ -55,7 +55,9 @@ public class Cluster<K, P extends ClusteringPoint<K>> {
     protected double distanceFromCentroid(ClusteringPoint<K> p) {
     	double squaredTotal=0;
     	for (K key : centroid.getAllClusteringKeys()) {
-       		squaredTotal += Math.pow(centroid.getClusteringValue(key)-p.getClusteringValue(key), 2);
+    		double a = centroid.getClusteringValue(key);
+    		double b = p.getClusteringValue(key);
+    		squaredTotal += Math.pow(a-b, 2);
     	}
     	
         return Math.sqrt(squaredTotal);

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;
 	}
@@ -50,7 +51,11 @@ public class YieldClusterPoint implements ClusteringPoint<String> {
 
 	@Override
 	public double getClusteringValue(String key) {
-
+		double v = getClusteringValueNaN(key);		
+		return Double.isNaN(v) ? 0.0 : v;  // LPJ with emulator provides NaN yields.  How to cluster these is tricky, but assuming they are 0 
+	}
+	
+	private double getClusteringValueNaN(String key) {
 		switch (key) {
 			case PASTURE:  return pasture;
 			case WHEAT_MIN:  return wheatMin;
@@ -60,9 +65,9 @@ public class YieldClusterPoint implements ClusteringPoint<String> {
 			case MAIZE_MAX:  return maizeMax;
 			case IRRIG_WATER_AVAL:  return irrigWaterAval;
 			case PROTECTED_AREA:  return protectedArea;
+			default:
+				throw new RuntimeException("YieldClusterPoint.getClusteringValue: got unknown value " + key);
 		}
-		LogWriter.printlnError("YieldClusterPoint.getClusteringValue: got unknown value " + key);
-		return Double.NaN;
 	}
 
 	@Override

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) {
+			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);
+			}
 		}
+		
+		// old style, if not wanting extrapolate or data doesn't allow
+		if (yieldsExtrapolated == null) {
+			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 (yM < 0 || yH < 0 || fM * yH > yM) {
+			LogWriter.printlnError(String.format("yieldAsymptoteSolve: less than zero yield or 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

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();
 	}