Changes for new LPJ data version. Also, some demand work.

GAMS/IntExtOpt.gms

@@ -119,7 +119,8 @@ $gdxin
        CROP_INCREASE_CALC(location)                     
        CROP_DECREASE_CALC(location)                     
        PASTURE_INCREASE_CONV_CALC(location)             
-       PASTURE_DECREASE_CONV_CALC(location)             
+       PASTURE_DECREASE_CONV_CALC(location)          
+       PASTURE_TOTAL_CHANGE_CONSTRAINT(location)   
        COST_EQ                                        total cost objective function;
  
  UNIT_COST_EQ(crop, location) .. unitCost(crop, location) =E=  (     baseCost(crop) * 0.5 +           
@@ -177,6 +178,7 @@ $gdxin
  CROP_DECREASE_CALC(location) .. cropDecrease(location) =G= -sum(crop_less_pasture, area(crop_less_pasture, location) - previousArea(crop_less_pasture, location));
  PASTURE_INCREASE_CONV_CALC(location) .. pastureIncrease(location) =G= area('pasture', location) - previousArea('pasture', location);
  PASTURE_DECREASE_CONV_CALC(location) .. pastureDecrease(location) =G= -(area('pasture', location) - previousArea('pasture', location));
+ PASTURE_TOTAL_CHANGE_CONSTRAINT(location) .. pastureIncrease(location) -pastureDecrease(location) =G= area('pasture', location) - previousArea('pasture', location);
   
  COST_EQ .. total_cost =E= 
          ( 

GAMS/demand.gms

@@ -3,49 +3,53 @@
  ALIAS (com, com2);
  ALIAS (com, com3);
       
- SCALAR income / 100 /;
+      
+ SCALAR income;
+      
+$gdxin %demandgdx% 
+$load income
+$gdxin    
       
  TABLE subElasticities(com, com2)
                      ruminants  wheat   maize    rice  oilcrops  pulses  starchyRoots
      monogastrics          0.7   0.15     0.1     0.1      0.1      0.1          0.1
      ruminants                    0.1     0.1     0.1      0.1      0.1          0.1
-     wheat                                0.6     0.6      0.4      0.4          0.4
-     maize                                        0.6      0.4      0.4          0.4
-     rice                                                  0.4      0.4          0.4
+     wheat                                0.6     0.6      0.3      0.3          0.3
+     maize                                        0.6      0.3      0.3          0.3
+     rice                                                  0.3      0.3          0.3
      oilcrops                                                       0.5          0.5
      pulses                                                                      0.5 ;
      
  subElasticities(com, com2)$subElasticities(com2, com) = subElasticities(com2, com);     
  
  PARAMETER calories(com)
-   / monogastrics          3,
-     ruminants             5,
-     wheat                 1,
-     maize                 1,
+   / monogastrics          1.5,
+     ruminants             1.7,
+     wheat                 1.02,
+     maize                 0.95,
      rice                  1,
-     oilcrops              2,
-     pulses                1,
-     starchyRoots          1  / ;
- calories(com) = 1;
+     oilcrops              1.3,
+     pulses                1.05,
+     starchyRoots          0.9  / ;
  
  PARAMETER idealDiet(com)
-   / monogastrics          1,
+   / monogastrics          1.1,
      ruminants             1,
-     wheat                 1,
-     maize                 1,
-     rice                  1,
-     oilcrops              1,
-     pulses                1,
-     starchyRoots          1  / ;
+     wheat                 0.9,
+     maize                 0.4,
+     rice                  0.5,
+     oilcrops              0.6,
+     pulses                0.5,
+     starchyRoots          0.8  / ;
       
  PARAMETER price(com)
-   / monogastrics          4,
-     ruminants             5,
-     wheat                 2,
-     maize                 2,
-     rice                  2,
-     oilcrops              2,
-     pulses                1,
+   / monogastrics          5,
+     ruminants             7.7,
+     wheat                 1.4,
+     maize                 1.25,
+     rice                  1.5,
+     oilcrops              3,
+     pulses                1.5,
      starchyRoots          1  / ;
        
  PARAMETER comUtility(com);
@@ -105,7 +109,7 @@
  
  total_cost_calc ..  totalCost =e= sum(com, consumption(com) * price(com));
  
- net_utility_calc ..  netUtility =e=   ((totalCost/income)) * 25 + foodUtility;
+ net_utility_calc ..  netUtility =e=   ((totalCost/income)) * 30 + foodUtility;
 
  MODEL DEMAND /ALL/ ;
  
@@ -121,4 +125,8 @@
   
  PARAMETER check(com);
  check(com) = 0;
-*display check;
\ No newline at end of file
+*display check;
+
+ Scalar ms 'model status', ss 'solve status'; 
+ ms=DEMAND.modelstat;
+ ss=DEMAND.solvestat;
\ No newline at end of file

debug_config.properties

 BASE_DIR=/Users/peteralexander/Documents/R_Workspace/UNPLUM
 #YIELD_DIR=/Users/peteralexander/Documents/LURG/LPJ/yield_anpp_rcps-2016-01-29/rcp8p5_past_fert/output-2016-01-29
 
-YIELD_DIR=/Users/peteralexander/Documents/LURG/LPJ/sam-03Nov2016
+YIELD_DIR=/Users/peteralexander/Documents/LURG/LPJ/10Jan2017/combined
 SPATIAL_DIR_NAME=halfdeg
 CELL_SIZE_X=0.5
-ANPP_FILENAME=anpp.out.2011-2015
-YIELD_FILENAME=yield.out.2011-2015
-IRRIG_MAX_WATER_FILENAME=gsirrigation.out.2011-2015
+#ANPP_FILENAME=anpp.out
+#YIELD_FILENAME=yield.out
+#IRRIG_MAX_WATER_FILENAME=gsirrigation.out
 MIN_FERT_AMOUNT=0
 MID_FERT_AMOUNT=200
 MAX_FERT_AMOUNT=1000
@@ -16,7 +16,7 @@ PASTURE_FERT_RESPONSE_FROM_LPJ=false
 CLEANUP_GAMS_DIR=false
 
 # Properties for testing
-CHANGE_YIELD_DATA_YEAR=false
+CHANGE_YIELD_DATA_YEAR=true
 CHANGE_DEMAND_YEAR=false
 DEBUG_LIMIT_COUNTRIES=false
 DEBUG_COUNTRY_NAME=Spain
@@ -25,7 +25,7 @@ SSP_SCENARIO=SSP2_v9_130325
 MAX_IMPORT_CHANGE=0.0
 MARKET_ADJ_PRICE=false
 
-IS_CALIBRATION_RUN = true
+IS_CALIBRATION_RUN = false
 
 #LAND_CHANGE_COST=.6
 
@@ -36,7 +36,7 @@ IS_CALIBRATION_RUN = true
 #NUM_CEREAL_CATEGORIES = 10
 #NUM_PASTURE_CATEGORIES = 3
 
-END_TIMESTEP=28
+END_TIMESTEP=18
 TIMESTEP_SIZE=5
 
 INTERPOLATE_OUTPUT_YEARS = false
\ No newline at end of file

src/ac/ed/lurg/ModelConfig.java

@@ -102,8 +102,8 @@ public class ModelConfig {
 	// yield data
 	public static final String YIELD_DIR = getProperty("YIELD_DIR");
 	public static final int LPJG_MONITOR_TIMEOUT_SEC = getIntProperty("LPJG_MONITOR_TIMEOUT", 60*60*2);
-	public static final String ANPP_FILENAME = getProperty("ANPP_FILENAME", "anpp_plum.out");
-	public static final String YIELD_FILENAME = getProperty("YIELD_FILENAME", "yield_plum.out");
+	public static final String ANPP_FILENAME = getProperty("ANPP_FILENAME", "anpp.out");
+	public static final String YIELD_FILENAME = getProperty("YIELD_FILENAME", "yield.out");
 	public static final boolean PASTURE_FERT_RESPONSE_FROM_LPJ = getBooleanProperty("PASTURE_FERT_RESPONSE_FROM_LPJ", true);;
 
 	public static final double CALIB_FACTOR_WHEAT = getDoubleProperty("CALIB_FACTOR_WHEAT", 0.966);
@@ -122,7 +122,7 @@ public class ModelConfig {
 	public static final String COUNTRY_BOUNDARY_FILE = SPATIAL_DATA_DIR + File.separator + "country_boundaries.asc";
 	public static final String IRRIGATION_COST_FILE = SPATIAL_DATA_DIR + File.separator + "irrigation_cost.asc";
 	public static final String IRRIGATION_CONSTRAINT_FILE = SPATIAL_DATA_DIR + File.separator + "blue_water_available_pseudoCRU_rcp8p5_2004_2013_grid_allhdyro_mm.txt";
-	public static final String IRRIG_MAX_WATER_FILENAME = getProperty("IRRIG_MAX_WATER_FILENAME", "gsirrigation_plum.out");
+	public static final String IRRIG_MAX_WATER_FILENAME = getProperty("IRRIG_MAX_WATER_FILENAME", "gsirrigation.out");
 	public static final String PROTECTED_AREAS_FILE = SPATIAL_DATA_DIR + File.separator + "protected_areas.txt";
 	
 	// Output

src/ac/ed/lurg/Timestep.java

@@ -88,6 +88,6 @@ public class Timestep {
 		if (ModelConfig.CHANGE_YIELD_DATA_YEAR)
 			return rootDir + File.separator + getYieldYear();
 		else
-			return rootDir;
+			return rootDir + File.separator + ModelConfig.BASE_YEAR;
 	}
 }

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

+package ac.ed.lurg.country.gams;
+
+import java.io.BufferedWriter;
+import java.io.File;
+import java.io.FileWriter;
+import java.io.IOException;
+import java.util.HashMap;
+import java.util.Map;
+
+import com.gams.api.GAMSDatabase;
+import com.gams.api.GAMSGlobals;
+import com.gams.api.GAMSJob;
+import com.gams.api.GAMSOptions;
+import com.gams.api.GAMSParameter;
+import com.gams.api.GAMSParameterRecord;
+import com.gams.api.GAMSVariable;
+import com.gams.api.GAMSVariableRecord;
+import com.gams.api.GAMSWorkspace;
+import com.gams.api.GAMSWorkspaceInfo;
+
+import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.utils.LogWriter;
+
+public class GamsDemandOptimiser {
+	private static final boolean DEBUG = true;
+
+	private double income;
+
+
+	public static void main(String[] args)  {
+		BufferedWriter demandWriter = null;
+
+		try {
+			demandWriter = new BufferedWriter(new FileWriter(ModelConfig.OUTPUT_DIR + "/demand.out", false));
+
+			StringBuffer sbHeader = new StringBuffer("income");
+			for (CropType crop : CropType.getImportedTypes())
+				sbHeader.append("," + crop.getGamsName());
+
+			demandWriter.write(sbHeader.toString());
+			demandWriter.newLine();
+
+			for (int i=0; i<50; i++) {
+				double income = 20.0 + i * 2.0;
+				GamsDemandOptimiser theModel = new GamsDemandOptimiser(income);
+				Map<CropType, Double> results = theModel.run();
+
+				StringBuffer sbData = new StringBuffer(String.format("%.6f", income));
+				for (CropType crop : CropType.getImportedTypes()) {
+					Double d = results.get(crop);
+					sbData.append(String.format(",%.6f", d));
+				}
+
+				demandWriter.write(sbData.toString());
+				demandWriter.newLine();
+			}
+		}
+		catch (IOException e) {
+			LogWriter.print(e);
+		}
+		finally {
+			if (demandWriter != null) {
+				try {
+					demandWriter.close();
+				} 
+				catch (IOException e) {
+					LogWriter.print(e);
+				}
+			}
+
+		}
+	}
+
+
+	public GamsDemandOptimiser(double income) {
+		this.income = income;
+	}
+
+	public Map<CropType, Double> run() {
+		File workingDirectory = new File(ModelConfig.TEMP_DIR, "GamsTmp");
+		workingDirectory.mkdir();
+
+		GAMSWorkspaceInfo  wsInfo  = new GAMSWorkspaceInfo();
+		wsInfo.setWorkingDirectory(workingDirectory.getAbsolutePath());
+		//	wsInfo.setDebugLevel(DebugLevel.VERBOSE);
+
+		GAMSWorkspace ws = new GAMSWorkspace(wsInfo);
+		GAMSDatabase inDB = ws.addDatabase();
+
+		setupInDB(inDB);
+
+		GAMSJob gamsJob = ws.addJobFromFile(ModelConfig.BASE_DIR + File.separator + "GAMS/demand.gms");
+		GAMSOptions opt = ws.addOptions();
+		//		opt.setAllModelTypes("conopt");
+		opt.defines("demandgdx", inDB.getName());
+
+		long startTime = System.currentTimeMillis();
+		gamsJob.run(opt, inDB);		
+		LogWriter.println("Took " + (System.currentTimeMillis() - startTime) + " ms to run");
+
+		return handleResults(gamsJob.OutDB());
+	}
+
+	private void setupInDB(GAMSDatabase inDB) {
+		addScalar(inDB, "income", income, 3);
+	}
+
+	private void addScalar(GAMSDatabase gamsDb, String recordName, double val, int places) {
+		GAMSParameter param = gamsDb.addParameter(recordName, 0);
+		setGamsParamValue(param.addRecord(), val, places);
+		if (DEBUG) LogWriter.println(recordName + ": " + val);
+	}
+
+	private void setGamsParamValue(GAMSParameterRecord param, double val, int places) {
+		double dOut = Math.round(val * Math.pow(10,places)) / Math.pow(10,places);
+		param.setValue(dOut);
+	}
+
+	private Map<CropType, Double> handleResults(GAMSDatabase outDB) {
+		int modelStatus = (int) outDB.getParameter("ms").findRecord().getValue();
+		LogWriter.println(String.format("\nModelstatus %s, Solvestatus %s", 
+				GAMSGlobals.ModelStat.lookup( modelStatus ),
+				GAMSGlobals.SolveStat.lookup((int) outDB.getParameter("ss").findRecord().getValue()) ));
+
+		GAMSVariable varConsump = outDB.getVariable("consumption");
+
+		double consumpValue;
+		Map<CropType, Double> consumptionData = new HashMap<CropType, Double>();
+
+		for (GAMSVariableRecord rec : varConsump) {
+			String itemName = rec.getKeys()[0];
+			consumpValue = rec.getLevel();
+
+			CropType cropType = CropType.getForGamsName(itemName);
+
+			consumptionData.put(cropType, consumpValue);
+			if (DEBUG) LogWriter.println(String.format("\n%s:\tconsumption= %.3f", itemName, consumpValue)); 
+		}
+
+		return consumptionData;
+	}
+}