ModelMain.java

package ac.ed.lurg;

import java.io.File;
import java.util.ArrayList;
import java.util.Arrays;
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.Country;
import ac.ed.lurg.country.CountryAgent;
import ac.ed.lurg.country.CountryBoundaryItem;
import ac.ed.lurg.country.CountryBoundaryReader;
import ac.ed.lurg.country.gams.GamsRasterOutput;
import ac.ed.lurg.demand.DemandManager;
import ac.ed.lurg.landuse.AreasItem;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.IntensitiesItem;
import ac.ed.lurg.landuse.Intensity;
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.RasterOutputer;
import ac.ed.lurg.types.CropToDoubleMap;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.FertiliserRate;
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.RasterHeaderDetails;
import ac.sac.raster.RasterKey;
import ac.sac.raster.RasterSet;

public class ModelMain {

	private Collection<CountryAgent> countryAgents;
	private Map<Country, List<RasterKey>> countryToKeysMap;
	private DemandManager demandManager;
	private Map<CropType, Double> prevWorldInputCost;
	private RasterHeaderDetails desiredProjection;

	public static void main(String[] args)  {
		ModelMain theModel = new ModelMain();
		theModel.setup();
		theModel.run();
	}

	/* setup models, reading inputs, etc. */
	private void setup() {
		desiredProjection = new RasterHeaderDetails(720, 360, -180, -90, 0.5, "999");
		countryToKeysMap = createCountryToKeysMap();
		demandManager = new DemandManager(ModelFitType.LOGISTIC, ModelConfig.SSP_SCENARIO);
		countryAgents = createCountryAgents();

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

	/* run the model */
	private void run() {
		for (int i = ModelConfig.START_TIMESTEP; i <= ModelConfig.END_TIMESTEP; i++)
			doTimestep(i);
	}

	private void doTimestep(int timestep) {
		YieldRaster yieldSurfaces = getYieldSurfaces(timestep);

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

		LogWriter.println("Timestep " + timestep);

		CropToDoubleMap totalQuantity = new CropToDoubleMap();
		CropToDoubleMap totalWorldInputCost = new CropToDoubleMap();
		//	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 = countryToKeysMap.get(ca.getCountry());
			YieldRaster countryYieldSurfaces = yieldSurfaces.getSubsetRasterForKeys(countryKeys);
						
			// do the optimization
			GamsRasterOutput result = ca.determineProduction(timestep, countryYieldSurfaces, prevWorldInputCost);

			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();
				CropUsageData d = entry.getValue();

				totalQuantity.incrementValue(c, d.getProduction());
				totalWorldInputCost.incrementValue(c, d.getProdCost());
			}

			/* after first timestep should look at trade balance and adjust appropriately
			for (CropType crop : CropType.values()) {
				CommodityData cd = result.getCommoditiesData().get(crop);
				double netImport = cd.getNetImports();

				if (netImport > 0)
					totalImportCommodities.incrementValue(crop, netImport);
				else
					totalExportCommodities.incrementValue(crop, -netImport);
			} */
		}

		prevWorldInputCost = totalWorldInputCost.divideBy(totalQuantity);
		
		// output results
		outputTimestepResults(timestep, globalIntensityRaster, globalCropAreaRaster, globalLocationIdRaster);
	}

	private void outputTimestepResults(int timestep, RasterSet<IntensitiesItem> intensityRaster, RasterSet<AreasItem> cropAreaRaster, RasterSet<IntegerRasterItem> locationIdRaster) {
		
		new RasterOutputer<IntegerRasterItem>(locationIdRaster, "locId" + timestep) {
			@Override
			public Double getValue(RasterKey location) {
				IntegerRasterItem locItem = results.get(location);
				if (locItem == null)
					return null;

				return (double)locItem.getInt();
			}

			@Override
			public int convertToPixelValue(double value) {
				return (int) (value * 3 + 1);
			}	
		}.writeOutput(ModelConfig.WRITE_JPEG_IMAGES);


		new RasterOutputer<IntensitiesItem>(intensityRaster, "wheatIntensity" + timestep) {
			@Override
			public Double getValue(RasterKey location) {
				IntensitiesItem intensityItem = results.get(location);
				if (intensityItem == null)
					return null;

				Intensity cropIntensity = intensityItem.getIntensity(CropType.WHEAT);

				return cropIntensity == null ? 0 : cropIntensity.getFertiliserIntensity();
			}

			@Override
			public int convertToPixelValue(double value) {
				return (int) (value * 20) + 1;
			}	
		}.writeOutput(ModelConfig.WRITE_JPEG_IMAGES);

		outputAreas(timestep, cropAreaRaster, CropType.WHEAT);
		outputAreas(timestep, cropAreaRaster, CropType.MAIZE);
		outputAreas(timestep, cropAreaRaster, CropType.PASTURE);
	}

	private void outputAreas(int timestep, RasterSet<AreasItem> cropAreaRaster, final CropType crop) {
		new RasterOutputer<AreasItem>(cropAreaRaster, crop.getGamsName() + "Area" + timestep) {
			@Override
			public Double getValue(RasterKey location) {
				AreasItem area = results.get(location);
				if (area == null)
					return null;

				return area.getCropArea(crop);
			}

			@Override
			public int convertToPixelValue(double value) {
				return value > 0 ? 10 : 1;
			}	
		}.writeOutput(ModelConfig.WRITE_JPEG_IMAGES);
	}

	public Map<Country, List<RasterKey>> createCountryToKeysMap() {

		RasterSet<CountryBoundaryItem> countryBoundaries = new RasterSet<CountryBoundaryItem>(desiredProjection) {
			private static final long serialVersionUID = -8449000692429399253L;
			protected CountryBoundaryItem createRasterData() {
				return new CountryBoundaryItem();
			}
		};

		CountryBoundaryReader countryReader = new CountryBoundaryReader(countryBoundaries);
		countryReader.getRasterDataFromFile(ModelConfig.COUNTRY_BOUNDARY_FILE);

		Map<Country, List<RasterKey>> countryMap = new HashMap<Country, List<RasterKey>>();

		for (Map.Entry<RasterKey, CountryBoundaryItem> entry : countryBoundaries.entrySet()) {
			Country c = entry.getValue().getCountry();
			if (c == null)
				continue;

			List<RasterKey> keys = countryMap.get(c);

			if (keys == null) {
				keys = new ArrayList<RasterKey>();
				countryMap.put(c, keys);
			}
			keys.add(entry.getKey());
		}
		return countryMap;
	}

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

		RasterSet<LandCoverItem> initLC = getInitialLandCover();
		RasterSet<IrrigationCostItem> allIrrigationCosts = getIrrigationCosts();
		Map<Country, Map<CropType, CropUsageData>> cropUsageDataMap = CropUsageData.readCommodityData();

		HashSet<String> countryExclusionList = new HashSet<String>(Arrays.asList("Sierra Leone", "Togo", "United Arab Emirates")); //"French Polynesia", "Cabo Verde", "Samoa", "Saint Vincent and the Grenadines"));

		for (Map.Entry<Country, List<RasterKey>> entry : countryToKeysMap.entrySet()) {
			Country country = entry.getKey();
			List<RasterKey> keys = entry.getValue();


			// DEBUG code
			if (!(country.getCountryName().equals("United States of Americaxx") || country.getCountryName().equals("Russian Federation"))) { //|| country.getCountryName().equals("China")
				continue;
			}

			if (demandManager.getPopulation(country, 2010) < 10 || countryExclusionList.contains(country.getCountryName())) {
				LogWriter.printlnError("Skipping " + country);
				continue;
			}

			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(country);

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

		return countryAgents;
	}

	private RasterSet<LandCoverItem> getInitialLandCover() {

		RasterSet<LandCoverItem> initLC = new RasterSet<LandCoverItem> (desiredProjection) {
			private static final long serialVersionUID = 4642550777741425501L;
			@Override
			protected LandCoverItem createRasterData() {
				return new LandCoverItem();
			}
		};

		String rootDir = ModelConfig.INITAL_LAND_COVER_DIR + File.separator;

		for (LandCoverType landType : LandCoverType.values()) {
			LandCoverReader lcReader = new LandCoverReader(initLC, landType);
			initLC = lcReader.getRasterDataFromFile(rootDir + landType.getFileName());
		}

		return initLC;
	}

	private YieldRaster getYieldSurfaces(int timestep) {
		/*	String rootDir = ModelConfig.YIELD_DIR + File.separator + (timestep + ModelConfig.BASE_YEAR) + File.separator;
		YieldRaster yieldSurfaces = null;

		for (CropType cropType : CropType.values()) {
			String cropDir = rootDir + cropType.getGamsName() + File.separator;

			for (YieldType yieldType : YieldType.values()) {
				YieldResponseMapReader yieldReader = new YieldResponseMapReader(yieldSurfaces, yieldType, cropType);
				yieldSurfaces = (YieldRaster)yieldReader.getRasterDataFromFile(cropDir + yieldType.getFileName()); 
			}
		} */

		LPJYieldResponseMapReader yieldReader = new LPJYieldResponseMapReader(ModelConfig.YIELD_DIR);

		for (FertiliserRate fr : FertiliserRate.values()) {
			yieldReader.getRasterDataFromFile(fr); 
		}

		return yieldReader.getYieldRaster();
	}

	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;
	}
}