Allow setting of subsidy rates

GAMS/IntExtOpt.gms

@@ -34,6 +34,7 @@
  PARAMETER irrigConstraint(location)         max water available for irrigation in litre per m2;
  PARAMETER minDemandPerCereal(cereal_crop)   min demand for each cereal crop as factor of all cereals;
  PARAMETER seedAndWasteRate(all_types)       rate of use for seed and waste combined;
+ PARAMETER subsidyRate(crop)                 rates of subsidy compared to costs;
  
  PARAMETER agriExpansionCost(location)       price for increase agricultural area varies based on amount of managed or unmanaged forest;
  SCALAR pastureIncCost                       price for increasing pasture area;
@@ -56,7 +57,7 @@ $gdxin %gdxincname%
 $load location, suitableLandArea, demand, agriExpansionCost, cropIncCost, pastureIncCost, cropDecCost, pastureDecCost
 $load previousArea, previousFertIntensity, previousIrrigIntensity, previousOtherIntensity, previousRuminantFeed, previousMonogastricFeed, previousImportAmount, previousExportAmount
 $load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth
-$load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImport, minNetImport, unhandledCropRate, setAsideRate, maxLandExpansionRate
+$load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxNetImport, minNetImport, unhandledCropRate, setAsideRate, maxLandExpansionRate, subsidyRate
 $load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandPerCereal, seedAndWasteRate, animalFeedFromOtherSources
 $gdxin    
  
@@ -203,7 +204,7 @@ $gdxin
 
  COST_EQ .. total_cost =E= 
          ( 
-              (   SUM((crop, location), area(crop, location) * unitCost(crop, location)) +
+              (   SUM((crop, location), area(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) +
                   
                   SUM(location, 
                      agriExpansionCost(location) * agriLandExpansion(location) +

GAMS/LUOpt.gms

@@ -35,6 +35,7 @@
  PARAMETER irrigConstraint(location)         max water available for irrigation in litre per m2;
  PARAMETER minDemandPerCereal(cereal_crop)   min demand for each cereal crop as factor of all cereals;
  PARAMETER seedAndWasteRate(all_types)       rate of use for seed and waste combined;
+ PARAMETER subsidyRate(crop)                 rates of subsidy compared to costs;
  
  PARAMETER agriExpansionCost(location)       price for increase agricultural area varies based on amount of managed or unmanaged forest;
  SCALAR pastureIncCost                       price for increasing pasture area;
@@ -57,7 +58,7 @@ $gdxin %gdxincname%
 $load location, suitableLandArea, demand, agriExpansionCost, cropIncCost, pastureIncCost, cropDecCost, pastureDecCost
 $load previousArea, previousFertIntensity, previousIrrigIntensity, previousOtherIntensity, previousRuminantFeed, previousMonogastricFeed, previousImportAmount
 $load yieldNone, yieldFertOnly, yieldIrrigOnly, yieldBoth
-$load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxProduction, minProduction, unhandledCropRate, setAsideRate, maxLandExpansionRate
+$load fertParam, irrigParam, otherIParam, exportPrices, importPrices, maxProduction, minProduction, unhandledCropRate, setAsideRate, maxLandExpansionRate, subsidyRate
 $load meatEfficency, otherICost, irrigCost, irrigMaxRate, irrigConstraint, fertiliserUnitCost, domesticPriceMarkup, minDemandPerCereal, seedAndWasteRate, animalFeedFromOtherSources
 
 $gdxin    
@@ -212,7 +213,7 @@ $gdxin
          ( 
               SUM(traded_commodity, domesticallyMetDemand(traded_commodity) * importPrices(traded_commodity) + export(traded_commodity) * exportPrices(traded_commodity)) - 
               
-              (   SUM((crop, location), area(crop, location) * unitCost(crop, location)) +
+              (   SUM((crop, location), area(crop, location) * unitCost(crop, location) * (1-subsidyRate(crop))) +
                   
                   SUM(location, 
                      agriExpansionCost(location) * agriLandExpansion(location) +

data/subsidyrates.csv

+country,crop,rate

src/ac/ed/lurg/ModelConfig.java

@@ -166,6 +166,7 @@ public class ModelConfig {
 	public static final String OTHER_WATER_USES_FILE = DATA_DIR + File.separator + "other_water_uses.csv";
 	public static final String BASE_CEREAL_DEMAND_FILE = DATA_DIR + File.separator + "base_cereal_demand.csv";
 	public static final String SHOCKS_PARAMETER_FILE = OUTPUT_DIR + File.separator+ "shocks.csv";
+	public static final String SUBSIDY_RATE_FILE = DATA_DIR + File.separator + "subsidyrates.csv";
 
 	// yield data
 	public static final String YIELD_DIR_BASE = getProperty("YIELD_DIR_BASE");
@@ -190,7 +191,6 @@ public class ModelConfig {
 	public static final double INITAL_PRICE_MAIZE = getDoubleProperty("INITAL_PRICE_MAIZE", 0.152 * ModelConfig.INITIAL_PRICE_SHIFT);
 	public static final double INITAL_PRICE_RICE = getDoubleProperty("INITAL_PRICE_RICE", 0.182 * ModelConfig.INITIAL_PRICE_SHIFT);
 	public static final double INITAL_PRICE_OILCROPS = getDoubleProperty("INITAL_PRICE_OILCROPS", (0.820 * .4 + 0.314 * .6) * 0.3 * ModelConfig.INITIAL_PRICE_SHIFT);
-
 	public static final double INITAL_PRICE_PULSES = getDoubleProperty("INITAL_PRICE_PULSES", 0.4 * ModelConfig.INITIAL_PRICE_SHIFT);
 	public static final double INITAL_PRICE_STARCHYROOTS = getDoubleProperty("INITAL_PRICE_STARCHYROOTS", 0.1 * ModelConfig.INITIAL_PRICE_SHIFT);
 	public static final double INITAL_PRICE_MONOGASTRICS = getDoubleProperty("INITAL_PRICE_MONOGASTRICS", 0.4 * 0.5 * ModelConfig.INITIAL_PRICE_SHIFT); // quantities is in feed equivalent term (0.4 is weighted average price per feed, and 0.5 accounts for mark-up for additional processing)
@@ -360,4 +360,6 @@ public class ModelConfig {
 	public static final boolean FORCE_PROTECTED_AREAS = IS_CALIBRATION_RUN ? false : getBooleanProperty("FORCE_PROTECTED_AREAS", false);
 	public static final int FORCE_PROTECTED_AREAS_START_YEAR  = getIntProperty("FORCE_PROTECTED_AREAS_START_YEAR", 2020);
 	public static final int FORCE_PROTECTED_AREAS_END_YEAR  = getIntProperty("FORCE_PROTECTED_AREAS_END_YEAR", 2050);
+	public static final String COUNTRY_FOR_SUBSIDY = getProperty("COUNTRY_FOR_SUBSIDY", "United States of America");
+	public static final double SUBSIDY_RATE = getDoubleProperty("SUBSIDY_RATE", 0.39);
 }

src/ac/ed/lurg/ModelMain.java

@@ -23,6 +23,7 @@ import ac.ed.lurg.country.CountryBoundaryRaster;
 import ac.ed.lurg.country.CountryBoundaryReader;
 import ac.ed.lurg.country.CountryPrice;
 import ac.ed.lurg.country.GlobalPrice;
+import ac.ed.lurg.country.SubsidyRateManager;
 import ac.ed.lurg.country.TradeManager;
 import ac.ed.lurg.demand.AbstractDemandManager;
 import ac.ed.lurg.demand.BaseConsumpManager;
@@ -44,8 +45,8 @@ import ac.ed.lurg.landuse.ProtectedAreasReader;
 import ac.ed.lurg.landuse.RunOffReader;
 import ac.ed.lurg.output.LandUseOutputer;
 import ac.ed.lurg.output.LpjgOutputer;
-import ac.ed.lurg.shock.YieldShocks;
 import ac.ed.lurg.shock.YieldShockReader;
+import ac.ed.lurg.shock.YieldShocks;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.types.CropType;
 import ac.ed.lurg.types.LandCoverType;
@@ -66,6 +67,7 @@ public class ModelMain {
 	private Collection<CountryAgent> countryAgents;
 	private CountryBoundaryRaster countryBoundaryRaster;
 	private AbstractDemandManager demandManager;
+	private SubsidyRateManager subsidyRateManager;
 	private TradeManager tradeManager;
 	private CompositeCountryManager compositeCountryManager;
 	private RasterHeaderDetails desiredProjection;
@@ -88,7 +90,8 @@ public class ModelMain {
 
 		BaseConsumpManager baseConsumpManager = new BaseConsumpManager();
 		compositeCountryManager = new CompositeCountryManager(baseConsumpManager);
-
+		subsidyRateManager = new SubsidyRateManager(compositeCountryManager);
+		
 		if (ModelConfig.DEMAND_FROM_FILE)
 			demandManager = new DemandManagerFromFile(compositeCountryManager);
 		else
@@ -487,8 +490,10 @@ private void writeDomesticProductionFile(Timestep timestep) {
 					LogWriter.printlnError("No initial land use for " + cc + ", so skipping");
 					continue;
 				}
+				
+				Map<CropType, Double> subsidyRates = subsidyRateManager.getSubsidyRates(cc);
 
-				CountryAgent ca = new CountryAgent(demandManager, cc, initCountryLandUse, countryCommodityData, countryTradeBarriers, yieldClusters);
+				CountryAgent ca = new CountryAgent(demandManager, cc, initCountryLandUse, countryCommodityData, countryTradeBarriers, yieldClusters, subsidyRates);
 				countryAgents.add(ca);
 				LogWriter.println("Creating country agent for: " + cc);
 			}

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

@@ -42,14 +42,17 @@ public class CountryAgent {
 	private Map<CropType, Double> tradeBarriers;
 	private RasterSet<IntegerRasterItem> yieldClusters;
 	private double animalFeedFromOtherSources;
+	private Map<CropType, Double> subsidyRates;
 
 	public CountryAgent(AbstractDemandManager demandManager,CompositeCountry country, RasterSet<LandUseItem> cropAreaRaster,
-			Map<CropType, CropUsageData> cropUsageData, Map<CropType, Double> tradeBarriers, RasterSet<IntegerRasterItem> yieldClusters) {
+			Map<CropType, CropUsageData> cropUsageData, Map<CropType, Double> tradeBarriers, RasterSet<IntegerRasterItem> yieldClusters,
+			Map<CropType, Double> subsidyRates) {
 
 		this.demandManager = demandManager;
 		this.country = country;
 		this.tradeBarriers = tradeBarriers;
 		this.yieldClusters = yieldClusters;
+		this.subsidyRates = subsidyRates;
 
 		GamsRasterOutput initialData = new GamsRasterOutput(cropAreaRaster, cropUsageData);
 		previousGamsRasterOutput = initialData;
@@ -266,7 +269,7 @@ public class CountryAgent {
 		Map<CropType, Double> minCerealFract = getMinCerealFraction(currentTimestep);
 
 		GamsCountryInput countryLevelInputs = new GamsCountryInput(country, currentProjectedDemand, currentCountryPrices, importConstraints, 
-				previousGamsRasterOutput.getCropUsageData(), minCerealFract, animalFeedFromOtherSources);	
+				previousGamsRasterOutput.getCropUsageData(), minCerealFract, animalFeedFromOtherSources, subsidyRates);	
 		GamsRasterInput input = new GamsRasterInput(currentTimestep, countryYieldSurfaces, previousGamsRasterOutput.getLandUses(), irrigData, countryLevelInputs);
 
 		return input;

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

+package ac.ed.lurg.country;
+
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+
+import ac.ed.lurg.ModelConfig;
+import ac.ed.lurg.types.CropType;
+import ac.ed.lurg.utils.LogWriter;
+import ac.ed.lurg.utils.StringTabularReader;
+
+public class SubsidyRateManager {
+	protected CompositeCountryManager compositeCountryManager;
+	private StringTabularReader tabularReader;
+
+	public SubsidyRateManager(CompositeCountryManager compositeCountryManager) {
+		this.compositeCountryManager = compositeCountryManager;
+		tabularReader = new StringTabularReader(",", new String[]{"country", "crop", "rate"});
+		tabularReader.read(ModelConfig.SUBSIDY_RATE_FILE);
+	}
+
+	public Map<CropType, Double> getSubsidyRates(CompositeCountry cc) {
+
+		Map<CropType, Double> subsidyRates = new HashMap<CropType, Double>();
+		for (SingleCountry c : compositeCountryManager.getAllForCompositeCountry(cc)) {
+			Map<CropType, Double> singleCMap = getSubsidyRates(c);
+
+			for (CropType crop : CropType.getAllItems()) {
+				Double newRateD = singleCMap.get(crop);
+				Double prevRateD = subsidyRates.get(crop);
+
+				// just takes max which isn't right. Should instead weight by production or something else
+				double newRate = Math.max(getDoubleOrZero(newRateD), getDoubleOrZero(prevRateD));
+				subsidyRates.put(crop, newRate);
+			}
+		}
+		return subsidyRates;
+	}
+
+	private double getDoubleOrZero(Double d) {
+		return (d==null ? 0.0 : d.doubleValue());
+	}
+
+	private Map<CropType, Double> getSubsidyRates(SingleCountry c) {
+		Map<CropType, Double> rates = new HashMap<CropType, Double>();
+
+		for (CropType crop : CropType.getAllItems()) {
+			Map<String, String> queryMap = new HashMap<String, String>();
+			queryMap.put("country", c.getCountryName());
+			queryMap.put("crop", crop.getGamsName());
+			try {
+				double subsidy = 0 ;
+				List<Map<String, String>> rows = tabularReader.query(queryMap);
+				if (rows.size() == 1) {
+					String subsidyS = rows.get(0).get("rate");
+					subsidy = Double.valueOf(subsidyS);
+				}
+				rates.put(crop, subsidy);
+			}
+			catch (Exception e) {
+					LogWriter.println("Problem finding cereal demand: " + crop.getFaoName() + ", " + c.getCountryName());
+			}
+		}
+		return rates;
+	}
+}

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

@@ -20,10 +20,11 @@ public class GamsCountryInput {
 	private Map<CropType, CropUsageData> previousCropUsageData;
 	private Map<CropType, Double> minCerealFraction;
 	private double animalFeedFromOtherSources;
+	private Map<CropType, Double> subsidyRates;
 
 	public GamsCountryInput(CompositeCountry country, Map<CommodityType, Double> projectedDemand, Map<CropType, CountryPrice> countryPrices, 
 			Map<CropType, TradeOrProductionConstraint> importConstraints, Map<CropType, CropUsageData> previousCropUsageData, Map<CropType, Double> minCerealFraction,
-			double animalFeedFromOtherSources) {
+			double animalFeedFromOtherSources, Map<CropType, Double> subsidyRates) {
 		super();
 		this.country = country;
 		this.projectedDemand = projectedDemand;
@@ -32,6 +33,7 @@ public class GamsCountryInput {
 		this.previousCropUsageData = previousCropUsageData;
 		this.minCerealFraction = minCerealFraction;
 		this.animalFeedFromOtherSources = animalFeedFromOtherSources;
+		this.subsidyRates = subsidyRates;
 	}
 		
 	public CompositeCountry getCountry() { 
@@ -83,4 +85,8 @@ public class GamsCountryInput {
 	public double getAnimalFeedFromOtherSources() {
 		return animalFeedFromOtherSources;
 	}
+	
+	public Map<CropType, Double> getSubsidyRates() {
+		return subsidyRates;
+	}
 }
\ No newline at end of file

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

@@ -232,6 +232,14 @@ public class GamsLocationOptimiser {
 			}
 		}
 
+		if (DEBUG) LogWriter.println("\nCrop, subsidy rate");
+		GAMSParameter subsideRateP = inDB.addParameter("subsidyRate", 1);
+		for (CropType crop : CropType.getNonMeatTypes()) {
+			double subsidyRate = countryInput.getSubsidyRates().get(crop);
+			if (DEBUG) LogWriter.println(String.format("%15s,\t %.4f", crop.getGamsName(), subsidyRate));
+			setGamsParamValue(subsideRateP.addRecord(crop.getGamsName()), subsidyRate, 4);
+		}
+		
 		if (DEBUG) LogWriter.println("\nImport-export, min trade/prod, max trade/prod, global import price, global export price, imports,   exports, ruminantFeed, monogastricFeed, seedAndWasteRate");
 		
 		GAMSParameter minTradeOrProdP = null;

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

@@ -3,7 +3,6 @@ package ac.ed.lurg.demand;
 import ac.ed.lurg.ModelConfig;
 import ac.ed.lurg.types.CommodityType;
 import ac.ed.lurg.types.ModelFitType;
-import ac.ed.lurg.utils.LogWriter;
 
 
 public class DemandCurve {