ModelMain.java

package ac.ed.lurg;

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Collection;
import java.util.HashMap;
import java.util.HashSet;
import java.util.List;
import java.util.Map;
import java.util.Map.Entry;

import ac.ed.lurg.country.CompositeCountry;
import ac.ed.lurg.country.CompositeCountryManager;
import ac.ed.lurg.country.CountryAgent;
import ac.ed.lurg.country.CountryBoundaryRaster;
import ac.ed.lurg.country.CountryBoundaryReader;
import ac.ed.lurg.country.gams.GamsRasterOutput;
import ac.ed.lurg.demand.BaseConsumpManager;
import ac.ed.lurg.demand.DemandManager;
import ac.ed.lurg.landuse.AreasItem;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.CropUsageReader;
import ac.ed.lurg.landuse.IntensitiesItem;
import ac.ed.lurg.landuse.IrrigationCostItem;
import ac.ed.lurg.landuse.IrrigiationCostReader;
import ac.ed.lurg.landuse.LandCoverItem;
import ac.ed.lurg.landuse.LandCoverReader;
import ac.ed.lurg.output.LpjgOutputer;
import ac.ed.lurg.types.CropToDoubleMap;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.types.ModelFitType;
import ac.ed.lurg.utils.LogWriter;
import ac.ed.lurg.yield.LPJYieldResponseMapReader;
import ac.ed.lurg.yield.YieldRaster;
import ac.sac.raster.IntegerRasterItem;
import ac.sac.raster.InterpolatingRasterSet;
import ac.sac.raster.RasterHeaderDetails;
import ac.sac.raster.RasterKey;
import ac.sac.raster.RasterOutputer;
import ac.sac.raster.RasterSet;

public class ModelMain {

	private Collection<CountryAgent> countryAgents;
	private CountryBoundaryRaster countryBoundaryRaster;
	private DemandManager demandManager;
	private CompositeCountryManager compositeCountryManager;
	private RasterHeaderDetails desiredProjection;

	private Map<CropType, Double> prevWorldPrices;
	private RasterSet<IntensitiesItem> prevIntensityRaster;
	private RasterSet<AreasItem> prevCropAreaRaster;
	
	
	public static void main(String[] args)  {
		ModelMain theModel = new ModelMain();
		theModel.setup();
		theModel.run();
	}

	/* setup models, reading inputs, etc. */
	private void setup() {
		desiredProjection = RasterHeaderDetails.getGlobalHeaderFromCellSize(ModelConfig.CELL_SIZE_X, ModelConfig.CELL_SIZE_Y, "999");
		
		BaseConsumpManager baseConsumpManager = new BaseConsumpManager();
		compositeCountryManager = new CompositeCountryManager(baseConsumpManager);
		demandManager = new DemandManager(ModelFitType.LOGISTIC, ModelConfig.SSP_SCENARIO, baseConsumpManager, compositeCountryManager);
		
		countryBoundaryRaster = getCountryBoundaryRaster();
		countryAgents = createCountryAgents(compositeCountryManager.getAll());

		// in first timestep we don't have this info, but ok as constrained to import/export specified amount
		prevWorldPrices = new HashMap<CropType, Double>();  
		for (CropType c : CropType.getImportedTypes()) 
			prevWorldPrices.put(c, 0.3);
		
		prevWorldPrices.put(CropType.STARCHY_ROOTS, 1.4);
	}

	/* run the model */
	private void run() {
		for (int i = ModelConfig.START_TIMESTEP; i <= ModelConfig.END_TIMESTEP; i++) {
			Timestep timestep = new Timestep(i);
			try {
				doTimestep(timestep);
			}
			catch (Exception e) {
				LpjgOutputer.writeMarkerFile(timestep.getYear(), true);
				throw new RuntimeException(e);
			}
		}
	}

	private void doTimestep(Timestep timestep) {
		LogWriter.println(timestep.toString());
		
		YieldRaster yieldSurfaces = getYieldSurfaces(timestep);  // this will wait for the marker file from LPJ if configured to do so

		//		YieldResponsesItem yresp = yieldSurfaces.getFromCoordinates(-50.0, -4.0);
		//		LogWriter.printlnError("Test key: " + yresp.getYieldMax(CropType.CEREALS)  + ", " + yresp.getYieldFertOnly(CropType.CEREALS) + ", " + yresp.getYieldIrrigOnly(CropType.CEREALS));

		CropToDoubleMap totalImportCommodities = new CropToDoubleMap();
		CropToDoubleMap totalExportCommodities = new CropToDoubleMap();

		RasterSet<IntensitiesItem> globalIntensityRaster = new RasterSet<IntensitiesItem>(desiredProjection);
		RasterSet<AreasItem> globalCropAreaRaster = new RasterSet<AreasItem>(desiredProjection);
		RasterSet<IntegerRasterItem> globalLocationIdRaster = new RasterSet<IntegerRasterItem>(desiredProjection);
				
		for (CountryAgent ca : countryAgents) {

			LogWriter.println("Country " + ca.getCountry());
			Collection<RasterKey> countryKeys = countryBoundaryRaster.getKeysFor(ca.getCountry());
			YieldRaster countryYieldSurfaces = yieldSurfaces.getSubsetRasterForKeys(countryKeys);
						
			GamsRasterOutput result = null;
			
			// do the optimization
			try {
				result = ca.determineProduction(timestep, countryYieldSurfaces, prevWorldPrices);

			}
			catch (Exception e) {
				LogWriter.printlnError("Exception processing " + ca.getCountry() + " will continue with other countries");			
				e.printStackTrace();
			}

			if (result == null) {
				LogWriter.printlnError("No results for " + ca.getCountry());
				continue;
			}

			// update global rasters
			globalIntensityRaster.putAll(result.getIntensityRaster());
			globalCropAreaRaster.putAll(result.getCropAreaRaster());
			globalLocationIdRaster.putAll(result.getLocationIdRaster());
			
			// Get values for world input costs
			Map<CropType, CropUsageData> cropUsage = result.getCropUsageData();

			for (Entry<CropType, CropUsageData> entry : cropUsage.entrySet()) {
				CropType c = entry.getKey();
				double countryNetImports = entry.getValue().getNetImports();

				if (countryNetImports > 0)
					totalImportCommodities.incrementValue(c, countryNetImports);
				else
					totalExportCommodities.incrementValue(c, -countryNetImports);
			}
		}
		
		// Look at trade balance and adjust appropriately
		for (CropType crop : CropType.getImportedTypes()) {
			double prevPrice = prevWorldPrices.get(crop);

			if (!totalImportCommodities.containsKey(crop) || !totalExportCommodities.containsKey(crop)) {
				LogWriter.printlnError(String.format("Price for %s not updated (still %.3f), as no imports or no exports for it", crop.getGamsName(), prevPrice));
				continue;
			}
			
			double imports = totalImportCommodities.get(crop);
			double exports = totalExportCommodities.get(crop) * (1.0-ModelConfig.TRANSPORT_LOSSES);
			double adjustedPrice = updateMarketPrices(prevPrice, imports, exports);
			LogWriter.println(String.format("Price for %s updated from %.3f (imports %.0f, exports %.0f) to %.3f ", crop.getGamsName(), prevPrice, imports, exports, adjustedPrice));
			prevWorldPrices.put(crop, adjustedPrice);
		}
		
		
		// output results
		outputTimestepResults(timestep, globalIntensityRaster, globalCropAreaRaster, globalLocationIdRaster, yieldSurfaces);
		
		// keep last to allow interpolation
		prevIntensityRaster = globalIntensityRaster;
		prevCropAreaRaster = globalCropAreaRaster;
	}
	
	public double updateMarketPrices(double previousPrice, double demand, double supply) {
		if (demand > 0 || supply > 0) {
			double ratio;
			
			if (demand > supply)
				ratio = (demand-supply)/demand;
			else
				ratio = (demand-supply)/supply;

			double adjustment = Math.exp(ratio * ModelConfig.MARKET_LAMBA);
			return previousPrice * adjustment;
		}
		else {
			return previousPrice;
		}
	}

	private BufferedWriter getFileWriter(Timestep timestep, String file, String columnHeadings) throws IOException {
		if (timestep.isInitialTimestep()) {
			FileWriter fstream = new FileWriter(file,false);
			BufferedWriter outputFile = new BufferedWriter(fstream);
			outputFile.write(columnHeadings);
			outputFile.newLine();
			return outputFile;
		}
		else {
			FileWriter fstream = new FileWriter(file,true);
			return new BufferedWriter(fstream);
		}
	}
	
	private double getTotalArea(LandCoverType lcType, RasterSet<AreasItem> globalCropAreaRaster) {
		double totalArea = 0;
		for (AreasItem item : globalCropAreaRaster.values()) {
			totalArea += item.getLandCoverArea(lcType);
		}
		return totalArea;
	}

	private void writeMarketFile(Timestep timestep, RasterSet<AreasItem> cropAreaRaster) {
		try {
			StringBuffer sbHeadings = new StringBuffer("Year, Cropland (Mha), Pasture (Mha), Natural (Mha)");
			for (CropType crop : CropType.getImportedTypes())
				sbHeadings.append(",Px_" + crop.getGamsName());
			
			BufferedWriter outputFile = getFileWriter(timestep, ModelConfig.TOTAL_LAND_COVER_FILE, sbHeadings.toString());

			StringBuffer sbData = new StringBuffer();
			sbData.append(String.format("%d,%.1f,%.1f,%.1f", 
					timestep.getYear(), getTotalArea(LandCoverType.CROPLAND, cropAreaRaster), getTotalArea(LandCoverType.PASTURE, cropAreaRaster), getTotalArea(LandCoverType.OTHER_NATURAL, cropAreaRaster)));
			
			for (CropType crop : CropType.getImportedTypes() )
				sbData.append(String.format(",%.3f", prevWorldPrices.get(crop)));

			outputFile.write(sbData.toString());
			outputFile.newLine();
			outputFile.close();
		}
		catch (IOException e) {
			e.printStackTrace();
		}
	}

	private void outputTimestepResults(Timestep timestep, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster, 
			RasterSet<IntegerRasterItem> locationIdRaster, YieldRaster yieldSurfaces) {
		
		writeMarketFile(timestep, cropAreaRaster);
		
		if (ModelConfig.OUTPUT_FOR_LPJG) {
			for (int outputYear : timestep.getYearsFromLast()) {
				LogWriter.printlnError("Outputing Year: " + outputYear);
				RasterSet<AreasItem> areasToOutput;
				
				if (outputYear == timestep.getYear()) {
					areasToOutput = cropAreaRaster;
				}
				else {
					InterpolatingRasterSet<AreasItem> intermediateAreas = new InterpolatingRasterSet<AreasItem>(cropAreaRaster.getHeaderDetails()) {
						private static final long serialVersionUID = 1306045141011047760L;
						protected AreasItem createRasterData() {
							return new AreasItem();
						}
					};
					intermediateAreas.setup(prevCropAreaRaster, cropAreaRaster, timestep.getPreviousTimestep().getYear(), timestep.getYear(), outputYear);
					areasToOutput = intermediateAreas;
				}
				
				LpjgOutputer lpjOutputer = new LpjgOutputer(outputYear, intensityRaster, areasToOutput, yieldSurfaces);
				lpjOutputer.writeOutput();
			}
		}
		
		outputLandCover(timestep.getYear(), cropAreaRaster, LandCoverType.CROPLAND);
		outputLandCover(timestep.getYear(), cropAreaRaster, LandCoverType.PASTURE); 
	}

	private void outputLandCover(int year, RasterSet<AreasItem> cropAreaRaster, final LandCoverType lcType) {
		new RasterOutputer<AreasItem>(cropAreaRaster, lcType.getName() + "Area" + year) {
			@Override
			public Double getValue(RasterKey location) {
				AreasItem area = results.get(location);
				if (area == null)
					return null;

				return area.getLandCoverArea(lcType);
			}
		}.writeOutput();
	}

	public CountryBoundaryRaster getCountryBoundaryRaster() {
		CountryBoundaryRaster countryBoundaries = new CountryBoundaryRaster(desiredProjection, compositeCountryManager);
		CountryBoundaryReader countryReader = new CountryBoundaryReader(countryBoundaries);
		countryReader.getRasterDataFromFile(ModelConfig.COUNTRY_BOUNDARY_FILE);
		
		return countryBoundaries;
	}

	public Collection<CountryAgent> createCountryAgents(Collection<CompositeCountry> countryGrouping) {
		Collection<CountryAgent> countryAgents = new HashSet<CountryAgent>();

		RasterSet<LandCoverItem> initLC = getInitialLandCover();
		RasterSet<IrrigationCostItem> allIrrigationCosts = getIrrigationCosts();
		Map<CompositeCountry, Map<CropType, CropUsageData>> cropUsageDataMap = new CropUsageReader(compositeCountryManager).getCommodityData();
				
		for (CompositeCountry cc : countryGrouping) {
			
			// DEBUG code
			if (ModelConfig.DEBUG_LIMIT_COUNTRIES) {
				if (!(cc.getName().equals("United States of America") || cc.getName().equals("Russian Federationxx") || cc.getName().equals("South Asia_otherxx")) ) {
					continue;
				}
			}
			
			List<RasterKey> keys = countryBoundaryRaster.getKeysFor(cc);
			RasterSet<LandCoverItem> initCountryLC = initLC.popSubsetForKeys(new RasterSet<LandCoverItem>(initLC.getHeaderDetails()), keys);
			RasterSet<IrrigationCostItem> irrigationCosts = allIrrigationCosts.popSubsetForKeys(new RasterSet<IrrigationCostItem>(allIrrigationCosts.getHeaderDetails()), keys);
			Map<CropType, CropUsageData> countryCommodityData = cropUsageDataMap.get(cc);

			if (countryCommodityData == null) {
				LogWriter.printlnError("No commodities data for " + cc + ", so skipping");
			}
			else if (initCountryLC.size() == 0) {
				LogWriter.printlnError("No initial land cover for " +cc + ", so skipping");
			}
			else {
				CountryAgent ca = new CountryAgent(demandManager, cc, initCountryLC, irrigationCosts, countryCommodityData);
				countryAgents.add(ca);
				LogWriter.println("Creating country agent for: " + cc );
			}
		}

		return countryAgents;
	}

	private RasterSet<LandCoverItem> getInitialLandCover() {
		LandCoverReader lcReader = new LandCoverReader(desiredProjection);
		RasterSet<LandCoverItem>  initLC = lcReader.getRasterDataFromFile(ModelConfig.INITAL_LAND_COVER_FILE);
		return initLC;
	}

	private YieldRaster getYieldSurfaces(Timestep timestep) {
		LPJYieldResponseMapReader yieldReader = new LPJYieldResponseMapReader(desiredProjection);
		return yieldReader.getRasterData(timestep); 
	}

	private RasterSet<IrrigationCostItem> getIrrigationCosts() {
		RasterSet<IrrigationCostItem> irigCosts = new RasterSet<IrrigationCostItem>(desiredProjection) {
			private static final long serialVersionUID = 8393130687550888654L;

			protected IrrigationCostItem createRasterData() {
				return new IrrigationCostItem();
			}
		};

		IrrigiationCostReader irigCostReader = new IrrigiationCostReader(irigCosts);
		irigCostReader.getRasterDataFromFile(ModelConfig.IRRIGATION_COST_FILE); 
		return irigCosts;
	}
}