From bc47f15efc69620cbd45c7ef1ee920ced4b72033 Mon Sep 17 00:00:00 2001
From: Peter Alexander <peter@blackhillock.co.uk>
Date: Thu, 10 Nov 2016 14:54:10 +0000
Subject: [PATCH] Some fixes. To try to stop pasture disappearing.

---
 GAMS/IntExtOpt.gms                            | 13 +++-------
 src/ac/ed/lurg/ModelConfig.java               |  2 +-
 .../country/gams/GamsRasterOptimiser.java     | 25 +++++++++++--------
 3 files changed, 19 insertions(+), 21 deletions(-)

diff --git a/GAMS/IntExtOpt.gms b/GAMS/IntExtOpt.gms
index 81cb9de9..bdda8f68 100644
--- a/GAMS/IntExtOpt.gms
+++ b/GAMS/IntExtOpt.gms
@@ -36,7 +36,7 @@
  SCALAR unhandledCropArea                    includes fruit veg forage crops set aside and failed crop;
  SCALAR domesticPriceMarkup                  factor price increased from cost of production;
 
-*$gdxin "/Users/peteralexander/Documents/R_Workspace/temp/GamsTest/_gams_java_1091455539.gdx"
+*$gdxin "/Users/peteralexander/Documents/R_Workspace/UNPLUM/temp/GamsTmp/_gams_java_gdb1.gdx"
 $gdxin %gdxincname% 
 $load location, suitableLandArea, previousArea, demand, landChangeCost
 $load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth
@@ -183,15 +183,10 @@ $gdxin
                ) * domesticPriceMarkup + 
               
                SUM(import_crop, importAmount(import_crop) * worldImportPrices(import_crop) - exportAmount(import_crop) * worldExportPrices(import_crop)) 
-         ) / 1000000;
+         );
  
  MODEL LAND_USE /ALL/ ;
  
-   
- fertI.L(crop, location) = 0.5;
- irrigI.L(crop, location) = 0.5;
- otherIntensity.L(crop, location) = 0.5;
-
  importAmount.L(import_crop)$((maxNetImport(import_crop) + minNetImport(import_crop)) gt 0) = (maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
  exportAmount.L(import_crop)$((maxNetImport(import_crop) + minNetImport(import_crop)) lt 0) = -(maxNetImport(import_crop) + minNetImport(import_crop)) / 2;
        
@@ -199,7 +194,7 @@ $gdxin
  SOLVE LAND_USE USING NLP MINIMIZING total_cost;
       
 * display agriLandExpansion.l, previousArea, area.l, net_supply.l, demand, feedAmount.l, yield.l, importAmount.l, exportAmount.l;
-
+* display unitCost.l, area.l, pastureDecrease.l, cropIncrease.l
 
 * Calculate summary information used in Java process 
  parameter totalProd(all_types);
@@ -225,7 +220,7 @@ $gdxin
  feedCost(feed_crop)$(netImportAmount(feed_crop) le 0 and totalProd(feed_crop) gt 0) = totalProdCost(feed_crop) * feedAmount.l(feed_crop) / totalProd(feed_crop);
  totalProdCost('meat') = sum(feed_crop, feedCost(feed_crop));
  
-* display totalProdCost, totalProd, totalImportCost, feedCost, totalProdCost
+* display totalProdCost, totalProd, netImportCost
          
  Scalar ms 'model status', ss 'solve status'; 
  ms=LAND_USE.modelstat;
diff --git a/src/ac/ed/lurg/ModelConfig.java b/src/ac/ed/lurg/ModelConfig.java
index 57c3eecc..45385279 100644
--- a/src/ac/ed/lurg/ModelConfig.java
+++ b/src/ac/ed/lurg/ModelConfig.java
@@ -139,7 +139,7 @@ public class ModelConfig {
 	// Fertiliser application rates in kg/ha
 	public static final double MIN_FERT_AMOUNT = getDoubleProperty("MIN_FERT_AMOUNT", 5.0);
 	public static final double MID_FERT_AMOUNT = getDoubleProperty("MID_FERT_AMOUNT", 20.0);
-	public static final double MAX_FERT_AMOUNT = getDoubleProperty("MID_FERT_AMOUNT", 200.0);
+	public static final double MAX_FERT_AMOUNT = getDoubleProperty("MAX_FERT_AMOUNT", 200.0);
 	public static final int FERT_AMOUNT_PADDING = getIntProperty("FERT_AMOUNT_PADDING", 3);;
 
 	// Other model parameters
diff --git a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
index 18c12dc2..e232e78e 100644
--- a/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
+++ b/src/ac/ed/lurg/country/gams/GamsRasterOptimiser.java
@@ -2,7 +2,6 @@ package ac.ed.lurg.country.gams;
 
 import java.util.ArrayList;
 import java.util.Collections;
-import java.util.HashMap;
 import java.util.HashSet;
 import java.util.List;
 import java.util.Map;
@@ -53,7 +52,7 @@ public class GamsRasterOptimiser {
 
 		GamsLocationOutput result = null;
 
-		// this code does similar calibration to that in the GAMS script
+		/*  This code used to try to do some calibration, but was causing more distortions than it fixed.
 		if (gamsInput.getCountryInput().isCalibrateToObserved()) {
 			Map<CropType, Double> cropAdjustments = new HashMap<CropType, Double>();  // recreate, but might have been null is this case anyway
 			for (CropType c : CropType.getNonMeatTypes()) {
@@ -63,8 +62,6 @@ public class GamsRasterOptimiser {
 			double baseCropland = LandUseItem.getTotalLandCover(gamsInput.getPreviousLandUse().values(), LandCoverType.CROPLAND);
 			double basePasture = LandUseItem.getTotalLandCover(gamsInput.getPreviousLandUse().values(), LandCoverType.PASTURE);
 
-			LogWriter.println("baseCropland=" + baseCropland + ", basePasture=" + basePasture);
-
 			for (int i=0; i<ModelConfig.NUM_CALIBRATION_ITERATIONS; i++) {
 
 				GamsLocationInput newInput = GamsLocationInput.createWithNewAdjustments(gamsInput, cropAdjustments);
@@ -87,10 +84,10 @@ public class GamsRasterOptimiser {
 				}
 			}
 		}
-		else {
+		else { */
 			GamsLocationOptimiser opti = new GamsLocationOptimiser(gamsInput);		
 			result = opti.run();
-		}
+	//	}
 		return result;
 	}
 
@@ -360,17 +357,19 @@ public class GamsRasterOptimiser {
 							aggIrig.setMaxIrrigAmount(crop, aggregateMean(aggIrig.getMaxIrrigAmount(crop), suitableAreaSoFar, irrigMax, suitableAreaThisTime));
 					}
 					
-					double areaSoFar = aggLandUse.getCropArea(crop);
-					double areaThisTime = landUseItem.getCropArea(crop);
-
 					for (YieldType yieldType : YieldType.values()) {
 						double y = yresp.getYield(yieldType, crop);
 						if (Double.isNaN(y))
 							logErrorWithCoord("Problem getting yield for type:" + yieldType + ", ", key, yieldRaster, crop);
-						else
-							aggYResp.setYield(yieldType, crop, aggregateMean(aggYResp.getYield(yieldType, crop), areaSoFar, y, areaThisTime));
+						else {
+							double yieldSoFar = aggYResp.getYield(yieldType, crop);
+							double updatedYield = aggregateMean(yieldSoFar, suitableAreaSoFar, y, suitableAreaThisTime);
+							aggYResp.setYield(yieldType, crop, updatedYield);
+						}
 					}
 				
+					double areaSoFar = aggLandUse.getCropArea(crop);
+					double areaThisTime = landUseItem.getCropArea(crop);
 					double updateCropland = (landUseItem.getLandCoverArea(LandCoverType.CROPLAND) + aggLandUse.getLandCoverArea(LandCoverType.CROPLAND));
 
 					if (updateCropland > 0)
@@ -402,6 +401,10 @@ public class GamsRasterOptimiser {
 				LogWriter.println("");
 			} */
 		}
+		
+		double baseCropland = LandUseItem.getTotalLandCover(aggregatedAreas.values(), LandCoverType.CROPLAND);
+		double basePasture = LandUseItem.getTotalLandCover(aggregatedAreas.values(), LandCoverType.PASTURE);
+		LogWriter.println("\nbaseCropland=" + baseCropland + ", basePasture=" + basePasture);
 
 		return new GamsLocationInput(rasterInputData.getTimestep(), aggregatedYields, aggregatedAreas, aggregatedIrrigCosts, rasterInputData.getCountryInput());
 	}
-- 
GitLab