Use previous demand results as initial conditions

GAMS/elasticDemand.gms

@@ -10,6 +10,9 @@ PARAMETER
   delta(i)     Parameter of MAIDADS
   tau(i)       Parameter of MAIDADS
   omega        Parameter of MAIDADS
+  previousU
+  previousSubs(i)
+  previousDisc(i)
 ;
 
 $gdxin %gdx_parameters% 
@@ -32,10 +35,10 @@ PARAMETER
 
 $gdxin %gdx_prices_and_gdp% 
 *$gdxin "../../GAMS/PricesAndGdp.gdx" 
-$load p=price, m=gdp_pc
+$load p=price, m=gdp_pc, previousU, previousSubs, previousDisc
 $gdxin
 
-display p, m , alpha, beta, kappa, delta, tau, omega;
+display p, m, alpha, beta, kappa, delta, tau, omega;
 
 w.up(j) = 1;
 
@@ -69,12 +72,14 @@ model MAIDADS_Sim "MAIDADS model for simulation" /
 MAIDADS_Sim.optfile = 1;
 option cns=path;
 
-u = 0;
-SubsistenceC(j)  = 0.05;
-DiscretionaryC(i)$p(i) = 0.05;
+u = previousU; 
+SubsistenceC(i)  = previousSubs(i);
+DiscretionaryC(i) = previousDisc(i);
 DiscretionaryC.LO(i) = 1E-11;
 SubsistenceC.LO(j) = 0;
-w(j)       = (SubsistenceC(j)+DiscretionaryC(j))*p(j)/m;
+w(j) = (SubsistenceC(j)+DiscretionaryC(j))*p(j)/m;
+
+display u.l, DiscretionaryC.l, SubsistenceC.l;
 
 solve MAIDADS_Sim using cns;
 
@@ -82,13 +87,8 @@ Scalar ms 'model status', ss 'solve status';
 ms=MAIDADS_Sim.modelstat;
 ss=MAIDADS_Sim.solvestat;
 
-display DiscretionaryC.l, alpha, beta, kappa,u.l;
-
-parameter totalDemand(i);
-totalDemand(i) = SubsistenceC(i) + DiscretionaryC(i);
-
-parameter disc(i);
-
-disc(i) = [m-sum(j,p(j)*SubsistenceC(j))];
+display DiscretionaryC.l, u.l;
 
-display disc;
+*parameter disc(i);
+*disc(i) = [m-sum(j,p(j)*SubsistenceC(j))];
+*display disc;

debug_config.properties

@@ -2,15 +2,15 @@ BASE_DIR=..
 
 YIELD_DIR=/Users/peteralexander/Documents/LURG/LPJ/LPJGPLUM_expt1.1_2006-2100_PLUM6xtra_20180412171654/rcp60
 
-DEBUG_LIMIT_COUNTRIES=false
-DEBUG_COUNTRY_NAME=China
+DEBUG_LIMIT_COUNTRIES=true
+DEBUG_COUNTRY_NAME=Japan
 
-IS_CALIBRATION_RUN = true
+IS_CALIBRATION_RUN = false
 GENERATE_NEW_YIELD_CLUSTERS = false
 NUM_YIELD_CLUSTERS=8000
 
-END_TIMESTEP=0
-TIMESTEP_SIZE=1
+END_TIMESTEP=5
+TIMESTEP_SIZE=5
 
 INTERPOLATE_OUTPUT_YEARS = false
 

src/ac/ed/lurg/ModelConfig.java

@@ -378,5 +378,4 @@ 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 double NON_FOOD_PRICE = getDoubleProperty("NON_FOOD_PRICE", 0.06);
 }

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

@@ -72,7 +72,7 @@ public abstract class AbstractCountryAgent {
 			LogWriter.println(commodity.getGamsName() + ": " + commPrice);
 		}
 		
-		prices.put(CommodityType.NONFOOD, ModelConfig.NON_FOOD_PRICE);
+		prices.put(CommodityType.NONFOOD, Double.NaN);
 		return prices;
 	}
 

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

+package ac.ed.lurg.country.gams;
+
+import ac.ed.lurg.types.CommodityType;
+
+public class GamsCommodityDemand {
+
+	private CommodityType commodity;
+	private double subsistence;
+	private double discretionary;
+
+	public GamsCommodityDemand(CommodityType commodity, double subsistence, double discretionary) {
+		this.commodity = commodity;
+		this.subsistence = subsistence;
+		this.discretionary = discretionary;
+	}
+
+	public double getSubsistence() {
+		return subsistence;
+	}
+
+	public double getDiscretionary() {
+		return discretionary;
+	}
+
+	public double getKcalTotal() {
+		return (subsistence+ discretionary) * 2000;
+	}
+
+	public double getPlumDemand() {
+		return getKcalTotal()*365/commodity.getkcalPerT();
+	}
+
+	public CommodityType getCommodity() {
+		return commodity;
+	}
+}
\ No newline at end of file

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

+package ac.ed.lurg.country.gams;
+
+import java.util.Map;
+
+import ac.ed.lurg.country.SingleCountry;
+import ac.ed.lurg.types.CommodityType;
+
+public class GamsDemandInput {
+	private SingleCountry country;
+	private double gdpPc;
+	private Map<CommodityType, Double> prices;
+	private double usaGdpPc;
+	private GamsDemandOutput previousGamsDemands;
+	
+	public GamsDemandInput(SingleCountry country, double gdpPc, Map<CommodityType, Double> prices, double usaGdpPc, GamsDemandOutput previousGamsDemands) {
+		super();
+		this.country = country;
+		this.gdpPc = gdpPc;
+		this.prices = prices;
+		this.usaGdpPc = usaGdpPc;
+		this.previousGamsDemands = previousGamsDemands;
+	}
+	
+	public SingleCountry getCountry() {
+		return country;
+	}
+	public double getGdpPc() {
+		return gdpPc;
+	}
+	public Map<CommodityType, Double> getPrices() {
+		return prices;
+	}
+	public double getUsaGdpPc() {
+		return usaGdpPc;
+	}
+	public GamsDemandOutput getPreviousDemands() {
+		return previousGamsDemands;
+	}
+}
\ No newline at end of file

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

 package ac.ed.lurg.country.gams;
 
 import java.io.File;
-import java.util.HashMap;
+import java.util.HashSet;
 import java.util.Map;
 import java.util.Vector;
 
@@ -11,29 +11,24 @@ 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.country.SingleCountry;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.utils.LogWriter;
 
 public class GamsDemandOptimiser {
 
-	private SingleCountry country;
-	private double gdpPc;
-	private Map<CommodityType, Double> prices;
-	private double usaGdpPc;
+	private GamsDemandInput inputData;
 
-	public GamsDemandOptimiser(SingleCountry country, double gdpPc, Map<CommodityType, Double> prices, double usaGdpPc) {
-		this.country = country;
-		this.gdpPc = gdpPc;
-		this.prices = prices;
-		this.usaGdpPc = usaGdpPc;
+	public GamsDemandOptimiser(GamsDemandInput inputData) {
+		this.inputData = inputData;
 	}
 	
-	public Map<CommodityType, Double> getDemandPc() {
+	public GamsDemandOutput getDemandPc() {
 		File workingDirectory = new File(ModelConfig.TEMP_DIR);
 		workingDirectory.mkdir();
 		GAMSWorkspaceInfo  wsInfo  = new GAMSWorkspaceInfo();
@@ -46,8 +41,6 @@ public class GamsDemandOptimiser {
 		GAMSJob gamsJob = ws.addJobFromFile(ModelConfig.DEMAND_GAMS_MODEL);
 		GAMSOptions opt = ws.addOptions();
 
-//		opt.setAllModelTypes("conopt");
-		
 		opt.defines("gdx_parameters", new File(ModelConfig.DEMAND_PARAM_FILE).getAbsolutePath());
 		opt.defines("gdx_prices_and_gdp", inDB.getName());
 
@@ -60,25 +53,49 @@ public class GamsDemandOptimiser {
 
 	private void setupInDB(GAMSDatabase inDB) {
 		GAMSParameter gdpPcP = inDB.addParameter("gdp_pc", 0);
-		LogWriter.println(String.format("gdp_pc:\t= %.5f", gdpPc)); 
+		double gdpPc = inputData.getGdpPc();
+		LogWriter.println(String.format("gdp_pc = %.5f", gdpPc)); 
 		setGamsParamValue(gdpPcP.addRecord(), gdpPc, 5);
+		GamsDemandOutput previousGamsDemands = inputData.getPreviousDemands();
+		
+		GAMSParameter prevUP = inDB.addParameter("previousU", 0);
+		double previousUtility = previousGamsDemands == null ? 0 : previousGamsDemands.getUtility();//-16.3829624069452 + 4.3659621833299 * Math.log10(gdpPc);
+		LogWriter.println(String.format("previousU = %.5f, estimated = %.5f", previousUtility, -16.3829624069452 + 4.3659621833299 * Math.log10(gdpPc))); 
+		setGamsParamValue(prevUP.addRecord(), previousUtility, 5);
 
 		LogWriter.println("\nPrice");	
-		GAMSParameter priceP = inDB.addParameter("price", 1);		
-				
-		for (Map.Entry<CommodityType, Double> entry : prices.entrySet()) {
+		GAMSParameter priceP = inDB.addParameter("price", 1);	
+		GAMSParameter prevSubsP = inDB.addParameter("previousSubs", 1);
+		GAMSParameter prevDiscP = inDB.addParameter("previousDisc", 1);
+
+		for (Map.Entry<CommodityType, Double> entry : inputData.getPrices().entrySet()) {
 			CommodityType comm = entry.getKey();
 			Vector<String> v = new Vector<String>();
 			v.add(comm.getGamsName());
+			double priceComm, prevSubs = 0.1, prevDisc = 0.1;
+			
+			if (comm == CommodityType.NONFOOD) {
+				prevSubs = 5;
+				prevDisc = 1;
+				priceComm = 45.376+53.036*gdpPc/inputData.getUsaGdpPc();
+			}
+			else
+				priceComm = entry.getValue() * 1000 / comm.getkcalPerT()*2000*365;  //price per calorie of commodity x per year assuming 2000 kcal a day for whole year
 			
-			double priceComm;
-			if (comm == CommodityType.NONFOOD)
-				priceComm = 45.376+53.036*gdpPc/usaGdpPc;
-			else 
-				priceComm = 	entry.getValue() * 1000 / comm.getkcalPerT()*2000*365;  //price per calorie of commodity x per year assuming 2000 kcal a day for whole year
-				
-			LogWriter.println(String.format("%14s:  price= %.4f", comm.getGamsName(), priceComm)); 
+			if (previousGamsDemands != null) {
+				GamsCommodityDemand demand = previousGamsDemands.getGamsDemands(comm);
+				if (demand == null)
+					LogWriter.printlnError("Null GamsCommodityDemand for " + comm);
+				else {
+					prevSubs = demand.getSubsistence();
+					prevDisc = demand.getDiscretionary();
+				}
+			}
+			
+			LogWriter.println(String.format("%14s:  price=%.4f, previousSubs=%.4f, previousDisc=%.4f", comm.getGamsName(), priceComm, prevSubs, prevDisc));
 			setGamsParamValue(priceP.addRecord(v), priceComm, 4);
+			setGamsParamValue(prevSubsP.addRecord(v), prevSubs, 4);
+			setGamsParamValue(prevDiscP.addRecord(v), prevDisc, 4);
 		}
 	}
 
@@ -88,36 +105,43 @@ public class GamsDemandOptimiser {
 		return dOut;
 	}
 		
-	private Map<CommodityType, Double> handleResults(GAMSDatabase outDB) {
+	private GamsDemandOutput handleResults(GAMSDatabase outDB) {
 		int modelStatus = (int) outDB.getParameter("ms").findRecord().getValue();
-		LogWriter.println(String.format("\nDemamd %s: Modelstatus %s, Solvestatus %s", country,
+		LogWriter.println(String.format("\nDemamd %s: Modelstatus %s, Solvestatus %s", inputData.getCountry(),
 				GAMSGlobals.ModelStat.lookup( modelStatus ),
 				GAMSGlobals.SolveStat.lookup((int) outDB.getParameter("ss").findRecord().getValue()) ));
-
-		GAMSParameter parmProd = outDB.getParameter("totalDemand");
 		
-		Map<CommodityType, Double> demandMap = new HashMap<CommodityType, Double>();
+		GAMSVariable varU = outDB.getVariable("u");
+		double utility = varU.getFirstRecord().getLevel();
+		LogWriter.println(String.format("u = %.4f", utility)); 
+		
+		GAMSVariable varSubs = outDB.getVariable("SubsistenceC");
+		GAMSVariable varDisc = outDB.getVariable("DiscretionaryC");
+
+		HashSet<GamsCommodityDemand> demandMap = new HashSet<GamsCommodityDemand>();
 		
-		for (GAMSParameterRecord rec : parmProd) {
+		for (GAMSVariableRecord rec : varSubs) {
 			String commodityName = rec.getKeys()[0];
 			CommodityType commodity = CommodityType.getForGamsName(commodityName, true);
+			
 			if (commodity == null)
 				LogWriter.println(String.format("Can't find commodity from Gams demand, so ignoring: " + commodityName)); 
 			else {
-				 //gams is calculating demand in number of 2000 kcal servings per day so * by 2000 and * 365 to get kcal per year
-				double demandGams = rec.getValue();
-				
+				double subsGams = rec.getLevel();
+				double discGams = varDisc.findRecord(commodityName).getLevel();
+
 				if (commodity == CommodityType.NONFOOD) { // we don't want to consider non-food, so not added to demandMap
-					LogWriter.println(String.format("%14s:  demand(gams)= %.3f", commodityName, demandGams)); 
+					LogWriter.println(String.format("%14s:  subs=%.3f, disc=%.3f", commodityName, subsGams, discGams)); 
 				}
-				else {
-					double demandPc = rec.getValue()*2000*365*1/commodity.getkcalPerT();
-					LogWriter.println(String.format("%14s:  demand(gams)= %.3f, demand(kcal)=%.1f, demand(plum)= %.3f", commodityName, demandGams, demandGams*2000, demandPc)); 
-					demandMap.put(commodity, demandPc);
+				else {	//gams is calculating demand in number of 2000 kcal servings per day so * by 2000 and * 365 to get kcal per year
+					GamsCommodityDemand demand = new GamsCommodityDemand(commodity, subsGams, discGams);
+					
+					LogWriter.println(String.format("%14s:  subs=%.3f, disc=%.3f, demand(kcal)=%.1f, demand(plum)= %.3f", commodityName, subsGams, discGams, demand.getKcalTotal(), demand.getPlumDemand())); 
+					demandMap.add(demand);
 				}
 			}
 		}
 		
-		return demandMap;
+		return new GamsDemandOutput(demandMap, utility);
 	}
 }

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

+package ac.ed.lurg.country.gams;
+
+import java.util.Collection;
+import java.util.HashMap;
+import java.util.Map;
+
+import ac.ed.lurg.types.CommodityType;
+
+public class GamsDemandOutput {
+	public Collection<GamsCommodityDemand> demands;
+	public double utility;
+
+	public GamsDemandOutput(Collection<GamsCommodityDemand> demands, double utility) {
+		this.demands = demands;
+		this.utility = utility;
+	}
+
+	public Map<CommodityType, Double> getPlumDemands() {
+		Map<CommodityType, Double> plumDemands = new HashMap<CommodityType, Double>();
+		
+		for (GamsCommodityDemand demand : demands)
+			plumDemands.put(demand.getCommodity(), demand.getPlumDemand());
+		
+		return plumDemands;
+	}
+	
+	public GamsCommodityDemand getGamsDemands(CommodityType comm) {
+		for (GamsCommodityDemand demand : demands)
+			if (demand.getCommodity() == comm)
+				return demand;
+		
+		return null;
+	}
+
+	public double getUtility() {
+		return utility;
+	}
+}

src/ac/ed/lurg/demand/ElasticDemandManager.java

@@ -6,12 +6,15 @@ import java.util.Map;
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.country.CompositeCountryManager;
 import ac.ed.lurg.country.SingleCountry;
+import ac.ed.lurg.country.gams.GamsDemandInput;
 import ac.ed.lurg.country.gams.GamsDemandOptimiser;
+import ac.ed.lurg.country.gams.GamsDemandOutput;
 import ac.ed.lurg.types.CommodityType;
 
 public class ElasticDemandManager extends AbstractSSPDemandManager {
 	
-	protected Map<SingleCountry, Map<CommodityType, Double>> baseCpcCache = new HashMap<SingleCountry, Map<CommodityType, Double>>();
+	private Map<SingleCountry, Map<CommodityType, Double>> baseCpcCache = new HashMap<SingleCountry, Map<CommodityType, Double>>();
+	private Map<SingleCountry, GamsDemandOutput> previousGamsDemands = new HashMap<SingleCountry, GamsDemandOutput>();
 	
 	public ElasticDemandManager(String ssp_scenario, BaseConsumpManager baseConsumpManager, CompositeCountryManager compositeCountryManager) {
 		super(ssp_scenario, baseConsumpManager, compositeCountryManager);
@@ -20,8 +23,12 @@ public class ElasticDemandManager extends AbstractSSPDemandManager {
 	@Override
 	Map<CommodityType, Double> getFoodDemandMap(SingleCountry c, double gdpPc, double population, Map<CommodityType, Double> prices, double usaGdpPc) {
 		
+		GamsDemandInput inputData = new GamsDemandInput(c, gdpPc, prices, usaGdpPc, previousGamsDemands.get(c));
+		
 		// Do the projection
-		Map<CommodityType, Double> foodPc = new GamsDemandOptimiser(c, gdpPc, prices, usaGdpPc).getDemandPc();
+		GamsDemandOutput gamsOutput = new GamsDemandOptimiser(inputData).getDemandPc();
+		previousGamsDemands.put(c, gamsOutput);
+		Map<CommodityType, Double> foodPc = gamsOutput.getPlumDemands();
 		
 		Map<CommodityType, Double> foodDemands = new HashMap<CommodityType, Double>();
 		SspData baseSspData = sspManager.get(ssp_scenario, ModelConfig.BASE_YEAR, c);