GamsLocationOptimiser.java

package ac.ed.lurg.country.gams;

import java.io.File;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;
import java.util.Vector;

import com.gams.api.GAMSDatabase;
import com.gams.api.GAMSException;
import com.gams.api.GAMSGlobals;
import com.gams.api.GAMSGlobals.ModelStat;
import com.gams.api.GAMSJob;
import com.gams.api.GAMSOptions;
import com.gams.api.GAMSParameter;
import com.gams.api.GAMSParameterRecord;
import com.gams.api.GAMSSet;
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.CountryPrice;
import ac.ed.lurg.country.ForestRecord;
import ac.ed.lurg.country.TradeConstraint;
import ac.ed.lurg.landuse.CarbonFluxItem;
import ac.ed.lurg.landuse.ConversionCostReader;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.Intensity;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.WoodYieldItem;
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.DoubleMap;
import ac.ed.lurg.utils.LazyHashMap;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.utils.TripleMap;
import ac.ed.lurg.yield.YieldResponsesItem;

public class GamsLocationOptimiser {

	private static final boolean DEBUG = false;

	private GamsLocationInput inputData;

	public GamsLocationOptimiser(GamsLocationInput inputData) {
		this.inputData = inputData;
	}

	public GamsLocationOutput run() {

		File workingDirectory = new File(ModelConfig.TEMP_DIR);
		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.GAMS_MODEL);
		GAMSOptions opt = ws.addOptions();

		opt.setAllModelTypes("conopt");

		opt.defines("gdxincname", inDB.getName());

		long startTime = System.currentTimeMillis();
		gamsJob.run(opt, inDB);	
		
		/*
		if (inputData.getCountryInput().getCountry().getName().equals("Peru & Ecuador"))
			LogWriter.println("" + inputData.getPreviousLandUse().get(2).getUnprotectedLandCoverArea(LandCoverType.CROPLAND));
		*/
		
		if (ModelConfig.CLEANUP_GAMS_DIR) 
			cleanup(ws.workingDirectory());
		
		LogWriter.println("Took " + (System.currentTimeMillis() - startTime) + " ms to run");
		

		return handleResults(gamsJob.OutDB());
		
	}

	private void setupInDB(GAMSDatabase inDB) {
		//if (DEBUG) LogWriter.println("\nLocation set");
		GAMSSet locationSet = inDB.addSet("location", 1);
		for (Integer locId : inputData.getPreviousLandUse().keySet()) {
			//if (DEBUG) LogWriter.println("     " + locId);
			locationSet.addRecord(locId.toString());
		}

		if (DEBUG) LogWriter.println("\nPrevious crop and land areas");	
		GAMSParameter prevCropP = inDB.addParameter("previousCropArea", 2);
		GAMSParameter prevFertIP = inDB.addParameter("previousFertIntensity", 2);
		GAMSParameter prevIrrigIP = inDB.addParameter("previousIrrigIntensity", 2);
		GAMSParameter prevOtherIP = inDB.addParameter("previousOtherIntensity", 2);
		
		GAMSParameter previousRuminantFeedP = inDB.addParameter("previousRuminantFeed", 1);
		GAMSParameter previousMonogastricFeedP = inDB.addParameter("previousMonogastricFeed", 1);
		GAMSParameter previousImportAmountP = inDB.addParameter("previousImportAmount", 1);
		GAMSParameter previousExportAmountP = inDB.addParameter("previousExportAmount", 1);
	
		GAMSParameter landP = inDB.addParameter("suitableLandArea", 1);
		GAMSParameter seedAndWasteRateP = inDB.addParameter("seedAndWasteRate", 1);	
		
		GAMSParameter prevLandCoverP = inDB.addParameter("previousLandCoverArea", 2);
		
		double totalAgriLand = 0;
		double totalSuitable = 0;

		for (Map.Entry<Integer, ? extends LandUseItem> entry : inputData.getPreviousLandUse().entrySet()) {
			Integer locationId = entry.getKey();
			String locString = Integer.toString(locationId);
			LandUseItem landUseItem = entry.getValue();
			
			double suitableLand = landUseItem.getSuitableArea();
			totalSuitable += suitableLand;
			if (DEBUG) LogWriter.println(String.format("  %d   %15s,\t %.3f", locationId, "suitableLand", suitableLand));
			setGamsParamValueTruncate(landP.addRecord(locString), suitableLand, 5);	
			
			for (CropType cropType : CropType.getNonMeatTypes()) {
				Vector<String> v = new Vector<String>();
				v.add(cropType.getGamsName());
				v.add(locString);
				
				double area;
				if (CropType.PASTURE == cropType)
					area = landUseItem.getLandCoverArea(LandCoverType.PASTURE);
				else
					area = landUseItem.getCropArea(cropType);
				
				double prevFertI, prevIrrigI, prevOtherI;
				Intensity intensity = landUseItem.getIntensity(cropType);

				if (intensity==null) { // could be first time through or this crop not previously grown in this location, so give it some default values
					prevFertI = CropType.PASTURE.equals(cropType) ? 0.0 : 0.5;
					prevIrrigI = CropType.PASTURE.equals(cropType) ? 0.0 : 0.5;
					prevOtherI = CropType.PASTURE.equals(cropType) ? 0.1 : 0.5;
				}
				else {
					prevOtherI = intensity.getOtherIntensity();
					if (prevOtherI == 0 & !CropType.PASTURE.equals(cropType)) {  // this is needed or optimizer gets a bit confused if some areas are kept as cropland but are not productive, due to zero other intensity (i.e. set aside)
						prevFertI = 0.5;
						prevIrrigI = 0.5;
					}
					else {
						prevFertI = intensity.getFertiliserIntensity();
						prevIrrigI = intensity.getIrrigationIntensity();
					}
				}
				
				if (DEBUG) LogWriter.println(String.format("  %d   %15s,\t %.2f,\t %.3f,\t %.3f,\t %.3f", locationId, cropType.getGamsName(), area, prevFertI, prevIrrigI, prevOtherI));
				
				setGamsParamValue(prevCropP.addRecord(v), area, 3);
				setGamsParamValue(prevFertIP.addRecord(v), prevFertI, 4);
				setGamsParamValue(prevIrrigIP.addRecord(v), prevIrrigI, 4);
				setGamsParamValue(prevOtherIP.addRecord(v), prevOtherI, 4);
				
				totalAgriLand += area;
				
			}
			
			// Previous land covers
			for (LandCoverType lc : LandCoverType.getConvertibleTypes()) {
				Vector<String> v = new Vector<String>();
				v.add(lc.getName());
				v.add(locString);
				setGamsParamValueTruncate(prevLandCoverP.addRecord(v), landUseItem.getUnprotectedLandCoverArea(lc), 5);
			}
				
		}
		
		if (DEBUG) LogWriter.println(String.format("  Total agricultural %.1f,\t suitable %.1f", totalAgriLand, totalSuitable));

		if (DEBUG) LogWriter.println("\nIrrigation data (cost, constraint)");
		GAMSParameter irrigCostP = inDB.addParameter("irrigCost", 1);		
		GAMSParameter irrigConstraintP = inDB.addParameter("irrigConstraint", 1);
		Map<Integer, ? extends IrrigationItem> irrigationData = inputData.getIrrigationCosts();
		
		double irrigCostMultiplier = (!ModelConfig.SHOCKS_POSSIBLE) ? 1.0 : ModelConfig.getParameter(Parameter.IRRIG_COST_MULTIPLIER, inputData.getTimestep().getYear());
		
		for (Entry<Integer, ? extends IrrigationItem> entry : irrigationData.entrySet()) {
			Integer locationId = entry.getKey();
			IrrigationItem irrigCostItem = entry.getValue();
			double irrigCost = irrigCostItem.getIrrigCost()*irrigCostMultiplier;
			double irrigConstraint = irrigCostItem.getIrrigConstraint();
			if (DEBUG) LogWriter.println(String.format("  %d  \t %.5f,\t %.1f", locationId, irrigCost, irrigConstraint));
			setGamsParamValue(irrigCostP.addRecord(Integer.toString(locationId)), irrigCost, 5);
			setGamsParamValue(irrigConstraintP.addRecord(Integer.toString(locationId)), irrigConstraint, 3);
		}

		if (DEBUG) LogWriter.println("\nDemand: " + inputData.getCountryInput().getCountry() + " " + inputData.getTimestep().getYear());
		GamsCountryInput countryInput = inputData.getCountryInput();
		addCommodityMapParm(inDB.addParameter("demand", 1), countryInput.getProjectedDemand(), -1);

		GAMSParameter minCerealFracP = inDB.addParameter("minDemandPerCereal", 1);
		GAMSParameter minOilcropsFracP = inDB.addParameter("minDemandPerOilcrop", 1);
		if (DEBUG) LogWriter.println("\nMinDemandFractions");

		for (Entry<CommodityType, Map<CropType, Double>> entry : countryInput.getMinDemandFractions().entrySet()) {
			CommodityType comm = entry.getKey();
			for (Map.Entry<CropType, Double> entry2 : entry.getValue().entrySet()) {
				CropType crop = entry2.getKey();
				double minCommFract = ModelConfig.LIMIT_DEMAND_FRACTION ? entry2.getValue() : 0.0;			
				GAMSParameter minCropsFracP = (comm == CommodityType.CEREALS ? minCerealFracP : minOilcropsFracP);
				setGamsParamValueTruncate(minCropsFracP.addRecord(crop.getGamsName()), minCommFract, 4);
				if (DEBUG) LogWriter.println(String.format("  %15s,  %10s,  %.4f", comm.getGamsName(), crop.getGamsName(), minCommFract));
			}
		}

		if (DEBUG) LogWriter.println("\nYield (fert/irrig) None/None, Max/None, None/Max, Max/Max, Shock,\t [fert p],\t [irrig p],\t {max irrig}");
		GAMSParameter yNoneP = inDB.addParameter("yieldNone", 2);
		GAMSParameter y_fert = inDB.addParameter("yieldFertOnly", 2);
		GAMSParameter y_irrig = inDB.addParameter("yieldIrrigOnly", 2);
		GAMSParameter y_both = inDB.addParameter("yieldBoth", 2);
		GAMSParameter y_shock = inDB.addParameter("yieldShock", 2);
		GAMSParameter fert_p = inDB.addParameter("fertParam", 2);
		GAMSParameter irrig_p = inDB.addParameter("irrigParam", 2);
		GAMSParameter irrigMaxP = inDB.addParameter("irrigMaxRate", 2);

		for (Entry<Integer, ? extends YieldResponsesItem> entry : inputData.getYields().entrySet()) {
			Integer locationId = entry.getKey();
			String locString = Integer.toString(locationId);
			YieldResponsesItem yresp = entry.getValue();
			IrrigationItem irrigationItem = irrigationData.get(locationId);

			for (CropType crop : CropType.getNonMeatTypes()) {
				Vector<String> v = new Vector<String>();
				v.add(crop.getGamsName());
				v.add(locString);
				
				if (crop.equals(CropType.SETASIDE)) {
					setGamsParamValue(irrigMaxP.addRecord(v), 1000, 3);  // need to set this to any positive value to give an incentive not to irrigate setaside
					continue;
				}
				
				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.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);
				setGamsParamValue(irrigMaxP.addRecord(v), maxIrrig, 3);
			}
		}

		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 minTradeP = null;
		GAMSParameter maxTradeP = null;
		minTradeP = inDB.addParameter("minNetImport", 1);
		maxTradeP = inDB.addParameter("maxNetImport", 1);
		
		GAMSParameter importPriceP = inDB.addParameter("importPrices", 1);
		GAMSParameter exportPriceP = inDB.addParameter("exportPrices", 1);
		
		for (CropType crop : CropType.getImportedTypes()) {		
			
			TradeConstraint iec = countryInput.getTradeConstraints().get(crop);
			CountryPrice gp = countryInput.getCountryPrices().get(crop);			
			double minTrade = iec.getMinConstraint();
			double maxTrade = iec.getMaxConstraint();
			double importPrice = gp.getImportPrice();
			double exportPrice = gp.getExportPrice();
			
			CropUsageData cu = countryInput.getPreviousCropUsageData().get(crop);
			double netImports = cu.getNetImportsExpected();
			double imports = netImports>0 ? netImports : 0;
			double exports = netImports<0 ? -netImports : 0;

			double ruminantFeed = cu.getRuminantFeed();
			double monogastricFeed = cu.getMonogastricFeed();
			double seedAndWasteRate = crop.getSeedAndWasteRate();
					
			if (DEBUG) LogWriter.println(String.format("     %15s, \t %5.1f, \t %5.1f, \t %5.3f, \t %5.3f,     \t %5.1f, \t %5.1f, \t %5.1f, \t %5.1f, \t %5.3f", 
					crop.getGamsName(), minTrade, maxTrade, importPrice, exportPrice, imports, exports, ruminantFeed, monogastricFeed, seedAndWasteRate));
			
			setGamsParamValue(minTradeP.addRecord(crop.getGamsName()), minTrade, 3);
			setGamsParamValue(maxTradeP.addRecord(crop.getGamsName()), maxTrade, 3);
			setGamsParamValue(importPriceP.addRecord(crop.getGamsName()), importPrice, 3);
			setGamsParamValue(exportPriceP.addRecord(crop.getGamsName()), exportPrice, 3);
			setGamsParamValue(previousImportAmountP.addRecord(crop.getGamsName()), imports, 3);
			setGamsParamValue(previousExportAmountP.addRecord(crop.getGamsName()), exports, 3);
			setGamsParamValue(previousRuminantFeedP.addRecord(crop.getGamsName()), ruminantFeed, 3);
			setGamsParamValue(previousMonogastricFeedP.addRecord(crop.getGamsName()), monogastricFeed, 3);
			setGamsParamValue(seedAndWasteRateP.addRecord(crop.getGamsName()), seedAndWasteRate, 3);
		}
		

		//below in case running without shocks to use original values in model config 
		
		double meatEff = (!ModelConfig.SHOCKS_POSSIBLE) ? ModelConfig.MEAT_EFFICIENCY: ModelConfig.getParameter(Parameter.MEAT_EFFICIENCY, inputData.getTimestep().getYear());
		double fertCost = (!ModelConfig.SHOCKS_POSSIBLE) ? ModelConfig.FERTILISER_MAX_COST : ModelConfig.getParameter(Parameter.FERTILISER_COST_PER_T, inputData.getTimestep().getYear()) * ModelConfig.MAX_FERT_AMOUNT/1000;
		double otherIntCost = (!ModelConfig.SHOCKS_POSSIBLE) ? ModelConfig.OTHER_INTENSITY_COST : ModelConfig.getParameter(Parameter.OTHER_INTENSITY_COST, inputData.getTimestep().getYear());

		
		LogWriter.print("\n");	
		addScalar(inDB, "meatEfficency", meatEff, 3);		
		addScalar(inDB, "fertiliserUnitCost", fertCost, 3);
		addScalar(inDB, "otherICost",otherIntCost, 3);
		addScalar(inDB, "otherIParam", ModelConfig.OTHER_INTENSITY_PARAM, 3);
		addScalar(inDB, "unhandledCropRate", ModelConfig.UNHANDLED_CROP_RATE, 3);
		addScalar(inDB, "setAsideRate", ModelConfig.SETASIDE_RATE, 5);
		addScalar(inDB, "domesticPriceMarkup", ModelConfig.DOMESTIC_PRICE_MARKUP, 3);
		double maxExpansion = 1.0;
		if (!ModelConfig.IS_CALIBRATION_RUN && countryInput.getCountry().getName().equals("China")) {
			maxExpansion = ModelConfig.MAX_CHINA_LAND_EXPANSION_RATE;
		}
		addScalar(inDB, "maxLandExpansionRate", maxExpansion, 3);
		
		addScalar(inDB, "carbonPrice", countryInput.getCarbonPrice().getExportPrice(), 3);
		addScalar(inDB, "woodPrice", countryInput.getTimberPrice().getExportPrice(), 3);
		
		// Carbon fluxes
		GAMSParameter cFluxRateP = inDB.addParameter("carbonFluxRate", 3);
		
		for (Entry<Integer, ? extends CarbonFluxItem> entry : inputData.getCarbonFluxes().entrySet()) {
			Integer locationId = entry.getKey();
			String locString = Integer.toString(locationId);
			CarbonFluxItem cFlux = entry.getValue();
			

			for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
				for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
					Vector<String> v = new Vector<String>();
					v.add(prevLC.getName());
					v.add(newLC.getName());
					v.add(locString);
					// TODO Need a better way to deal with missing values
					if (cFlux.checkForKeys(prevLC, newLC))
						setGamsParamValue(cFluxRateP.addRecord(v), cFlux.getCarbonFlux(prevLC, newLC), 3);
				}
			}		
		}
		
		// Wood yields
		Map<Integer, ? extends WoodYieldItem> woodYieldData = inputData.getWoodYields();
		Map<Integer, Double> naturalWoodYieldFactor = inputData.getNaturalWoodYieldFactor();
		
		// Scale yields for natural areas
		for (Entry<Integer, Double> entry : naturalWoodYieldFactor.entrySet()) {
			int locId = entry.getKey();
			double factor = entry.getValue();
			for (LandCoverType toLC : LandCoverType.getConvertibleTypes()) {
				WoodYieldItem item = woodYieldData.get(locId);
				double maxYield = item.getWoodYield(LandCoverType.NATURAL, toLC);
				item.setWoodYield(LandCoverType.NATURAL, toLC, maxYield * factor);
			}
		}
		
		GAMSParameter woodYieldP = inDB.addParameter("woodYield", 3);
		
		for (Entry<Integer, ? extends WoodYieldItem> entry : woodYieldData.entrySet()) {
			Integer locationId = entry.getKey();
			String locString = Integer.toString(locationId);
			WoodYieldItem wYield = entry.getValue();

			for (LandCoverType prevLC : LandCoverType.getConvertibleTypes()) {
				for (LandCoverType newLC : LandCoverType.getConvertibleTypes()) {
					Vector<String> v = new Vector<String>();
					v.add(prevLC.getName());
					v.add(newLC.getName());
					v.add(locString);
					if (wYield.checkForKeys(prevLC, newLC))
						setGamsParamValue(woodYieldP.addRecord(v), wYield.getWoodYield(prevLC, newLC), 3);
				}
			}			
		}
		
		// Timber yields
		GAMSParameter timberForestYieldP = inDB.addParameter("timberForestYield", 1);
		
		Map<Integer, Map<LandCoverType, Double>> timberYield = inputData.getTimberYield().getMap();
		
		if (!timberYield.isEmpty()) {
			for (Entry<Integer, Map<LandCoverType, Double>> locMap : timberYield.entrySet()) {
				Integer locationId = locMap.getKey();
				String locString = Integer.toString(locationId);
				
				for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) {
					Vector<String> v = new Vector<String>();
					LandCoverType landCover = coverMap.getKey();
					if (!landCover.equals(LandCoverType.TIMBER_FOREST))
						continue;
					v.add(locString);
					setGamsParamValueTruncate(timberForestYieldP.addRecord(v), coverMap.getValue(), 5);
				}
			}
		}
		
		// Minimum land covers
		GAMSParameter minimumLandCoverP = inDB.addParameter("minimumLandCoverArea", 2);
				
		for (Entry<Integer, Map<LandCoverType, Double>> locMap : inputData.getMinimumLandCover().getMap().entrySet()) {
			Integer locationId = locMap.getKey();
			String locString = Integer.toString(locationId);
			
			for (Entry<LandCoverType, Double> coverMap : locMap.getValue().entrySet()) {
				Vector<String> v = new Vector<String>();
				LandCoverType landCover = coverMap.getKey();
				v.add(landCover.getName());
				v.add(locString);
				double minCover = coverMap.getValue();
				double previousCover = inputData.getPreviousLandUse().get(locationId).getUnprotectedLandCoverArea(landCover);
				// Correction for floating point errors
				double minCoverCorrected = Math.min(minCover, previousCover);
				setGamsParamValueTruncate(minimumLandCoverP.addRecord(v), minCoverCorrected, 5);
			}
		}
		
		// Land cover conversion cost
		GAMSParameter conversionCostP = inDB.addParameter("conversionCost", 2);
		DoubleMap<LandCoverType, LandCoverType, Double> conversionCosts = inputData.getConversionCosts();
		
		for (Entry<LandCoverType, Map<LandCoverType, Double>> fromMap : conversionCosts.getMap().entrySet()) {
			String fromName = fromMap.getKey().getName();
			for (Entry<LandCoverType, Double> toMap : fromMap.getValue().entrySet()) {
				String toName = toMap.getKey().getName();
				double cost = toMap.getValue();
				Vector<String> v = new Vector<String>();
				v.add(fromName);
				v.add(toName);
				setGamsParamValue(conversionCostP.addRecord(v), cost, 5);
			}
		}
	}

	private void addScalar(GAMSDatabase gamsDb, String recordName, double val, int places) {
		GAMSParameter param = gamsDb.addParameter(recordName, 0);
		double dOut = setGamsParamValue(param.addRecord(), val, places);
		if (DEBUG) LogWriter.println(recordName + ": " + dOut);
	}

	private double setGamsParamValue(GAMSParameterRecord param, double val, int places) {
		double dOut = places >= 0 ? Math.round(val * Math.pow(10,places)) / Math.pow(10,places) : val;
		param.setValue(dOut);
		return dOut;
	}
	
	private double setGamsParamValueTruncate(GAMSParameterRecord param, double val, int places) {
		double dOut = ((int)(val * Math.pow(10,places))) / Math.pow(10,places);
		param.setValue(dOut);
		return dOut;
	}
	
	@SuppressWarnings("serial")
	private GamsLocationOutput handleResults(GAMSDatabase outDB) {
		int modelStatusInt = (int) outDB.getParameter("ms").findRecord().getValue();
		ModelStat modelStatus = GAMSGlobals.ModelStat.lookup(modelStatusInt);
		String contextString = String.format("%s %s: Modelstatus %s, Solvestatus %s", 
    			inputData.getCountryInput().getCountry(),
			inputData.getTimestep().getYear(),
			modelStatus.toString(),
			GAMSGlobals.SolveStat.lookup((int) outDB.getParameter("ss").findRecord().getValue()));
        LogWriter.println("\n" + contextString);

		if (modelStatus != ModelStat.OPTIMAL_LOCAL) {
			LogWriter.printlnError("Critical!!! Land use incorrectly solved. " + contextString);
		}
		
		GAMSVariable varAreas = outDB.getVariable("cropArea");
		GAMSVariable varFertIntensities = outDB.getVariable("fertI");
		GAMSVariable varIrrigIntensities = outDB.getVariable("irrigI");
		GAMSVariable varOtherIntensities = outDB.getVariable("otherIntensity");
		GAMSVariable varRuminantFeed = outDB.getVariable("ruminantFeed");
		GAMSVariable varMonogastricFeed = outDB.getVariable("monogastricFeed");
		GAMSParameter parmNetImports = outDB.getParameter("netImportAmount");
		GAMSParameter parmNetImportCost = outDB.getParameter("netImportCost");
		GAMSVariable varYields = outDB.getVariable("yield");
		GAMSVariable varUnitEnergies = outDB.getVariable("unitCost");
		GAMSParameter parmProd = outDB.getParameter("totalProd");
		GAMSParameter parmProdCost = outDB.getParameter("totalProdCost");
		GAMSParameter parmCroplandArea = outDB.getParameter("totalCropland");
		GAMSParameter parmTotalArea = outDB.getParameter("totalArea");
		GAMSParameter parmProdShock = outDB.getParameter("productionShock");		
		
		double totalCropArea = 0;
		double totalPastureArea = 0;
		double area, fertIntensity, irrigIntensity, otherIntensity = Double.NaN, ruminantFeed, monogastricFeed, netImport, netImportCost, yield, unitCost, prod, prodCost;	
		
		final LazyHashMap<Integer, LandUseItem> landUses = new LazyHashMap<Integer, LandUseItem>() { 
			protected LandUseItem createValue() { return new LandUseItem(); }
		};		
		
		Map<Integer, ? extends IrrigationItem> allIrrigationRefData = inputData.getIrrigationCosts();

		Map<CropType, CropUsageData> cropUsageData = new HashMap<CropType, CropUsageData>();

		for (GAMSVariableRecord rec : varAreas) {
			String itemName = rec.getKeys()[0];
			String locationName = rec.getKeys()[1];
			area = rec.getLevel();
			fertIntensity = varFertIntensities.findRecord(itemName, locationName).getLevel();
			irrigIntensity = varIrrigIntensities.findRecord(itemName, locationName).getLevel();
			otherIntensity = varOtherIntensities.findRecord(itemName, locationName).getLevel();
			yield = varYields.findRecord(itemName, locationName).getLevel();
			unitCost = varUnitEnergies.findRecord(itemName, locationName).getLevel();

			int locId = Integer.parseInt(locationName);
			CropType cropType = CropType.getForGamsName(itemName);

			if (!cropUsageData.containsKey(cropType)) {  // then we must not have seen this crop type before, so need to do all non location specific stuff
				ruminantFeed = varRuminantFeed.findRecord(itemName).getLevel();
				monogastricFeed = varMonogastricFeed.findRecord(itemName).getLevel();
				netImport = cropType.isImportedCrop() ? getParmValue(parmNetImports, itemName) : 0;
				netImportCost = cropType.isImportedCrop() ? getParmValue(parmNetImportCost, itemName) : 0;
				prod =  getParmValue(parmProd, itemName);
				prodCost = getParmValue(parmProdCost, itemName);
				double totalArea = getParmValue(parmTotalArea, itemName);
				double prodShock = getParmValue(parmProdShock, itemName);

				cropUsageData.put(cropType, new CropUsageData(ruminantFeed, monogastricFeed, netImport, netImportCost, prod, prodCost, totalArea, prodShock));
				if (DEBUG) LogWriter.println(String.format("\n%s:\tarea= %.1f,\tmonogastricFeed= %.1f,\truminantFeed= %.1f,\tnetImports= %.3f,\tnetImportCost= %.3f,\tprod= %.3f, \tprodCost= %.3f,\tprodCostRate= %.3f,\tprodShock= %.3f", itemName, totalArea, monogastricFeed, ruminantFeed, netImport, netImportCost, prod, prodCost, prodCost/prod, prodShock)); 
			
			}

			LandUseItem landUseItem = landUses.lazyGet(locId); // returns landUseItem for location. If does not exist, creates new one

			if (area > 0) { 

				if (DEBUG) LogWriter.println(String.format("\t location %s, %s:\tarea= %.1f,\tfert= %.3f,\tirrg= %.3f,\tintensity= %.3f", locationName, itemName, area, fertIntensity, irrigIntensity, otherIntensity));				
				IrrigationItem irrigRefData = allIrrigationRefData.get(locId);
				landUseItem.setIntensity(cropType, new Intensity(fertIntensity, irrigIntensity, otherIntensity, yield, unitCost, irrigRefData.getMaxIrrigAmount(cropType)));
			}
			
			
			double croplandArea = getParmValue(parmCroplandArea, locationName);
			if (cropType.equals(CropType.PASTURE)) {
				landUseItem.setLandCoverArea(LandCoverType.PASTURE, area);
				totalPastureArea += area;
			}
			else {
				landUseItem.setLandCoverArea(LandCoverType.CROPLAND, croplandArea);  // will set this multiple times, once for each arable crop, but doesn't really matter
				totalCropArea += area;
			}
			
			landUseItem.setCropFraction(cropType, croplandArea > 0 ? area/croplandArea : 0);
			
		}

		for (CropType meatTypes : CropType.getMeatTypes()) {
			netImport = getParmValue(parmNetImports, meatTypes.getGamsName());
			netImportCost= getParmValue(parmNetImportCost, meatTypes.getGamsName());
			prod = getParmValue(parmProd, meatTypes.getGamsName());
			prodCost = getParmValue(parmProdCost, meatTypes.getGamsName());

			cropUsageData.put(meatTypes, new CropUsageData(0.0, 0.0, netImport, netImportCost, prod, prodCost, Double.NaN, 0));
			if (DEBUG) LogWriter.println(String.format("\n%s:\t\t\t\t\tnetImports= %.3f,\tnetImportCost= %.3f,\tprod= %.3f,\tprodCost= %.3f", meatTypes.getGamsName(), netImport, netImportCost, prod, prodCost)); 
		}
		LogWriter.println(String.format("\n%s %s: Total area= %.1f (crop=%.1f, pasture %.1f)", 
				inputData.getCountryInput().getCountry(), inputData.getTimestep().getYear(), 
				totalCropArea+totalPastureArea, totalCropArea, totalPastureArea));
			
		// Land cover change
		GAMSVariable varLandCoverChange = outDB.getVariable("landCoverChange");
		
		TripleMap<Integer, LandCoverType, LandCoverType, Double> landCoverChanges = new TripleMap<Integer, LandCoverType, LandCoverType, Double>();
		
		for (GAMSVariableRecord rec : varLandCoverChange) {
			String fromLC = rec.getKeys()[0];
			String toLC = rec.getKeys()[1];
			String locationName = rec.getKeys()[2];
			int locId = Integer.parseInt(locationName);
			
			double change = rec.getLevel();
			
			if (!fromLC.equals(toLC)) { //Important as don't want to include unchanged land cover
				landCoverChanges.put(locId, LandCoverType.getForName(fromLC), LandCoverType.getForName(toLC), change);
			}					
		}
		
		// Minimum land cover additions. Need to keep track of new forest areas to restrict conversion
		// Also keep track of natural areas
		List<ForestRecord> newForest = new ArrayList<ForestRecord>();
		DoubleMap<Integer, LandCoverType, Double> totalArea = new DoubleMap<Integer, LandCoverType, Double>();
		DoubleMap<Integer, LandCoverType, Double> areaTimesYield = new DoubleMap<Integer, LandCoverType, Double>();

		for (Entry<Integer, DoubleMap<LandCoverType, LandCoverType, Double>> locMap : landCoverChanges.getMap().entrySet()) {
			int locId = locMap.getKey();
			DoubleMap<LandCoverType, LandCoverType, Double> changeMap = locMap.getValue();
			for (Entry<LandCoverType, Map<LandCoverType, Double>> fromMap : changeMap.getMap().entrySet()) {
				LandCoverType fromLC = fromMap.getKey();
				for (Entry<LandCoverType, Double> toMap : fromMap.getValue().entrySet()) {
					LandCoverType toLC = toMap.getKey();
					double newArea = toMap.getValue();
					if (!toLC.isNatural() || newArea == 0.0)
						continue;
					Double minArea = inputData.getMinimumLandCover().get(locId, fromLC);
					minArea = (minArea == null) ? 0.0 : minArea;
					double newForestArea = newArea - minArea; // subtract unchanged forest area
					double woodYield = (toLC.equals(LandCoverType.TIMBER_FOREST)) ? inputData.getWoodYields().get(locId).getWoodYield(fromLC, toLC) : 0.0;					
					totalArea.addTo(locId, toLC, newForestArea);
					areaTimesYield.addTo(locId, toLC, woodYield * newForestArea);
				}
			}
		}
		
		// Calculate weighted average yield
		for (Entry<Integer, Map<LandCoverType, Double>> locMap : totalArea.getMap().entrySet()) {
			Integer locId = locMap.getKey();
			for (Entry<LandCoverType, Double> lcMap : locMap.getValue().entrySet()) {
				LandCoverType lcType = lcMap.getKey();
				double a = lcMap.getValue();
				double ya = areaTimesYield.get(locId, lcType);
				double y = ya / a;
				newForest.add(new ForestRecord(locId, lcType, a, y));
			}
		}
		
			
		// reserved Forest which was deforested
		GAMSVariable varReservedDeforested = outDB.getVariable("reservedAreaDeforested");
		DoubleMap<Integer, LandCoverType, Double> reservedDeforested = new DoubleMap<Integer, LandCoverType, Double>();
		
		for (GAMSVariableRecord rec : varReservedDeforested) {
			String landCover = rec.getKeys()[0];
			String locationName = rec.getKeys()[1];
			int locId = Integer.parseInt(locationName);
			
			double deforestedArea = rec.getLevel();
			
			if (deforestedArea > 0.00001) {
				reservedDeforested.put(locId, LandCoverType.getForName(landCover), deforestedArea);
			}							
		}
		
		// Carbon flux
		double netCarbonFlux = outDB.getParameter("netCarbonFlux").getFirstRecord().getValue();
		
		// Timber harvest
		double totalTimberProduction = outDB.getParameter("totalTimberProduction").getFirstRecord().getValue();

		GamsLocationOutput results = new GamsLocationOutput(modelStatus, landUses, cropUsageData, landCoverChanges, newForest, reservedDeforested, netCarbonFlux, totalTimberProduction);
		return results;
	}
	
	/*
	private double getPrevUnmanagedForestProp(Integer locId) {
		double prevOtherNaturalArea = inputData.getPreviousLandUse().get(locId).getUnprotectedLandCoverArea(LandCoverType.OTHER_NATURAL);
		double prevUnmanagedForestArea = inputData.getPreviousLandUse().get(locId).getUnprotectedLandCoverArea(LandCoverType.UNMANAGED_FOREST);
		double prevUnmanagedForestProp = (prevOtherNaturalArea + prevUnmanagedForestArea) != 0 ? prevUnmanagedForestArea / (prevUnmanagedForestArea + prevOtherNaturalArea) : 0;
		return prevUnmanagedForestProp;
	}
	*/

	private double getParmValue(GAMSParameter aParm, String itemName) {
		try {
			GAMSParameterRecord record = aParm.findRecord(itemName);
			double d = record.getValue();
			return d;
		}
		catch (GAMSException gamsEx) {
			//LogWriter.println("GAMSException thrown for " + itemName);
			return 0;
		}
	}

	private void addCommodityMapParm(GAMSParameter parm, Map<CommodityType, Double> itemMap, int places) {
		for (Map.Entry<CommodityType, Double> entry : itemMap.entrySet()) {
			double d = entry.getValue();
			if (DEBUG) LogWriter.println(String.format("     %15s,\t %.1f", entry.getKey().getGamsName(), d));
			if (!Double.isNaN(d))
				setGamsParamValue(parm.addRecord(entry.getKey().getGamsName()), d, places);
		}
	}


	private void cleanup(String directory)  {
		File directoryToDelete = new File(directory);
		String files[] = directoryToDelete.list();
		for (String file : files) {
			File fileToDelete = new File(directoryToDelete, file);
			try {
				fileToDelete.delete();
			} catch(Exception e){
				LogWriter.print(e);
			}
		}
		try {
			directoryToDelete.delete();
		} catch(Exception e) {
			LogWriter.print(e);
		}
	}
}