ModelMain.java

package ac.ed.lurg;

import java.io.BufferedWriter;
import java.io.File;
import java.io.FileInputStream;
import java.io.FileOutputStream;
import java.io.IOException;
import java.io.ObjectInputStream;
import java.io.ObjectOutputStream;
import java.util.Collection;
import java.util.Map;
import java.util.Map.Entry;

import ac.ed.lurg.country.AbstractCountryAgent;
import ac.ed.lurg.country.AnimalRateManager;
import ac.ed.lurg.country.CompositeCountry;
import ac.ed.lurg.country.CompositeCountryManager;
import ac.ed.lurg.country.CountryAgentManager;
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.demand.AbstractDemandManager;
import ac.ed.lurg.demand.BaseConsumpManager;
import ac.ed.lurg.demand.CalorieManager;
import ac.ed.lurg.demand.DemandManagerFromFile;
import ac.ed.lurg.demand.DemandManagerSSP;
import ac.ed.lurg.demand.ElasticDemandManager;
import ac.ed.lurg.landuse.CropUsageData;
import ac.ed.lurg.landuse.CropUsageReader;
import ac.ed.lurg.landuse.FPUManager;
import ac.ed.lurg.landuse.IrrigationConstraintReader;
import ac.ed.lurg.landuse.IrrigationItem;
import ac.ed.lurg.landuse.IrrigationMaxAmountReader;
import ac.ed.lurg.landuse.IrrigationRasterSet;
import ac.ed.lurg.landuse.IrrigiationCostReader;
import ac.ed.lurg.landuse.LandCoverItem;
import ac.ed.lurg.landuse.LandCoverReader;
import ac.ed.lurg.landuse.LandUseItem;
import ac.ed.lurg.landuse.MaxCropAreaReader;
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.types.CommodityType;
import ac.ed.lurg.types.CropType;
import ac.ed.lurg.types.LandCoverType;
import ac.ed.lurg.utils.FileWriterHelper;
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.IntegerRasterReader;
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 CountryAgentManager countryAgents;
	private CountryBoundaryRaster countryBoundaryRaster;
	private AbstractDemandManager demandManager;
	private AnimalRateManager animalRateManager;
	private CompositeCountryManager compositeCountryManager;
	LPJYieldResponseMapReader lpjYieldReader;
	private RasterHeaderDetails desiredProjection;

	private InternationalMarket internationalMarket;
	private IrrigationRasterSet currentIrrigationData;
	private RasterSet<LandUseItem> globalLandUseRaster;
	private RasterSet<IntegerRasterItem> clusterIdRaster;

	public static void main(String[] args) {
		ModelMain theModel = new ModelMain();
	    	System.out.println("Working Directory = " + System.getProperty("user.dir"));
		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();
		CalorieManager calorieManager = new CalorieManager();
		compositeCountryManager = new CompositeCountryManager(baseConsumpManager);
		lpjYieldReader = new LPJYieldResponseMapReader(desiredProjection);
		animalRateManager = new AnimalRateManager(compositeCountryManager);
		
		if (ModelConfig.DEMAND_FROM_FILE)
			demandManager = new DemandManagerFromFile(compositeCountryManager,calorieManager);
		else if (ModelConfig.PRICE_ELASTIC_DEMAND)
			demandManager = new ElasticDemandManager(ModelConfig.SSP_SCENARIO, baseConsumpManager,calorieManager, compositeCountryManager);
		else
			demandManager = new DemandManagerSSP(ModelConfig.SSP_SCENARIO, baseConsumpManager,calorieManager, compositeCountryManager);
		
		if (ModelConfig.SHOCKS_POSSIBLE) ModelConfig.readInShocksFile();
			
		currentIrrigationData = getFixedIrrigationData();
		countryBoundaryRaster = getCountryBoundaryRaster();
		clusterIdRaster = ModelConfig.GENERATE_NEW_YIELD_CLUSTERS ? new RasterSet<IntegerRasterItem>(desiredProjection) : getClusterRaster();

		globalLandUseRaster = new RasterSet<LandUseItem>(desiredProjection);
		internationalMarket = new InternationalMarket();
		createCountryAgents(compositeCountryManager.getAll());
	}

	/* 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: " + timestep.toString());

		YieldRaster yieldSurfaces = getYieldSurfaces(timestep); // this will wait for the marker file from LPJ if configured to do so
		getUpdateIrrigationData(timestep, yieldSurfaces); // updating currentIrrigationData

		// When running half earth we can to alter protected areas data at a point in time
		if(ModelConfig.HALFEARTH && ModelConfig.FORCE_PROTECTED_AREAS_START_YEAR == timestep.getYear() && !ModelConfig.IS_CALIBRATION_RUN) {
			new ProtectedAreasReader(globalLandUseRaster).getRasterDataFromFile(ModelConfig.HALF_EARTH_FILE);
			countryAgents.updateProtectedAreasForAll(globalLandUseRaster);
		}
		
		double previousGen2EcDDemand = (timestep.isInitialTimestep() || ModelConfig.IS_CALIBRATION_RUN ) ? 0: demandManager.getSecondGenBioenergyDemand(timestep.getPreviousTimestep());
		double gen2EcDDemand = demandManager.getSecondGenBioenergyDemand(ModelConfig.IS_CALIBRATION_RUN ? new Timestep(1) : timestep);
		double gen2Increase = (gen2EcDDemand>previousGen2EcDDemand) ? gen2EcDDemand - previousGen2EcDDemand : 0.0;
		
		countryAgents.determineProductionForAll(timestep, yieldSurfaces, currentIrrigationData, gen2Increase);
		internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, timestep);
		
		// loop for iterations.  Could check within a tolerance using internationalMarket.findMaxPriceDiff, not doing so as doesn't find a solution due to inelastic consumption
		for (int i=0; i < ModelConfig.DEMAND_RECALC_MAX_ITERATIONS; i++) {
			LogWriter.println("\n++ Re-estimating prices and demand " + i);
			countryAgents.recalculateDemandForAll(); // recalculate demand from new prices
			countryAgents.updateNetImportsForAll(); // calculate imports and exports
			internationalMarket.determineInternationalTrade(countryAgents.getAll(), gen2EcDDemand, timestep); // calculate prices
		}
		internationalMarket.applyPriceShocks(timestep);
		
		// output results
		outputTimestepResults(timestep, globalLandUseRaster);
	}

	private void writeLandCoverFile(Timestep timestep, RasterSet<LandUseItem> landUseRaster) {
		try {
			StringBuffer sbHeadings = new StringBuffer("Year,Cropland,Pasture,ManForest,UnmanForest,Natural,AbPasture,Suitable,EnergyCrop,FertCrop,IrrigCrop");
			BufferedWriter outputFile = FileWriterHelper.getFileWriter(timestep, ModelConfig.LAND_COVER_OUTPUT_FILE, sbHeadings.toString());

			StringBuffer sbData = new StringBuffer();
			sbData.append(String.format("%d,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f,%.1f", timestep.getYear(),
					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.CROPLAND),
					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.PASTURE),
					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.MANAGED_FOREST),
					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.UNMANAGED_FOREST),
					LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.OTHER_NATURAL),
					LandUseItem.getAbandonedPasture(landUseRaster.values()),
					LandUseItem.getSuitableTotal(landUseRaster.values(), timestep.getYear()))
					);
						
			sbData.append(String.format(",%.1f", LandUseItem.getTotalCropArea(landUseRaster.values(), CropType.ENERGY_CROPS)));
			sbData.append(String.format(",%.1f", LandUseItem.getFertiliserTotal(landUseRaster.values(), CropType.getCropsLessPasture()) / 1000));
			sbData.append(String.format(",%.1f", LandUseItem.getIrrigationTotal(landUseRaster.values(), CropType.getCropsLessPasture())));

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

	private void writeGlobalMarketFile(Timestep timestep) {
		try {
			BufferedWriter outputFile = FileWriterHelper.getFileWriter(timestep, ModelConfig.PRICES_OUTPUT_FILE, "Year,Crop,Imports (Mt),Exports (Mt),New export price, Stock Levels (Mt)");
			internationalMarket.writeGlobalMarketFile(timestep, outputFile);
			outputFile.close();
		} catch (IOException e) {
			LogWriter.print(e);
		}
	}
	
	private void writeGlobalFoodBalanceSheet(Timestep timestep, RasterSet<LandUseItem> landUseRaster) {
		try {
			BufferedWriter outputFile = FileWriterHelper.getFileWriter(timestep, ModelConfig.FOOD_BALANCE_SHEET_FILE, "Year,Crop,Production,Imports,Export,TransportLosses,StockVar,Supply,MonogastricsFeed,RuminantsFeed,SeedAndOtherLosses,FoodAnd1stGen,ProdArea");
			Map<CropType, GlobalPrice> worldPrices = internationalMarket.getWorldPrices();
			double harvestedAndFallowArea = 0;
			
			for (CropType crop : CropType.getCropsLessPasture()) {
				GlobalPrice priceQuantity = worldPrices.get(crop); // some specific logic for import/exports and this has been aggregated already, so best to use it
				double prod=0, prodArea=0, feedMonogastrics=0, feedRuminants=0;
				double exportsBeforeTL=0, imports=0, transportloss=0, stockChange=0;
				
				if (priceQuantity != null) {
					exportsBeforeTL = priceQuantity.getExportsBeforeTransportLoss();
					imports = priceQuantity.getImportAmount();
					transportloss = priceQuantity.getTransportLosses();
					stockChange = priceQuantity.getStockChange();
				}
				
				for (AbstractCountryAgent ca : countryAgents.getAll()) {
					Map<CropType, CropUsageData> allCropUsage = ca.getCropUsageData();
					CropUsageData cropUsage = allCropUsage.get(crop);
					if (cropUsage != null) {
						prod += cropUsage.getProductionExpected();
						prodArea += cropUsage.getArea();
						feedMonogastrics += cropUsage.getMonogastricFeed();
						feedRuminants += cropUsage.getRuminantFeed();
					}
				}

				double seedAndWaste = prod * crop.getSeedAndWasteRate();
				double netSupply = prod - exportsBeforeTL + imports;
				double foodAnd1stGen = netSupply - feedMonogastrics - feedRuminants - seedAndWaste;
				
				if (!crop.equals(CropType.SETASIDE))
					prodArea *= (1-ModelConfig.UNHANDLED_CROP_RATE);  // remove unhandled crop area
				
				harvestedAndFallowArea += prodArea;
				
				StringBuffer sbData = new StringBuffer();
				sbData.append(String.format("%d,%s", timestep.getYear(), crop.getGamsName()));
				sbData.append(String.format(",%.2f", prod));
				sbData.append(String.format(",%.2f,%.2f,%.2f,%.2f", imports, exportsBeforeTL, transportloss, stockChange));
				sbData.append(String.format(",%.2f,%.2f,%.2f,%.2f,%.2f", netSupply, feedMonogastrics, feedRuminants, seedAndWaste, foodAnd1stGen));
				sbData.append(String.format(",%.2f", prodArea));
				outputFile.write(sbData.toString());
				outputFile.newLine();
			}

			double cropLandArea = LandUseItem.getTotalLandCover(landUseRaster.values(), LandCoverType.CROPLAND);
			double unhandledArea = cropLandArea - harvestedAndFallowArea;
			outputFile.write(String.format("%d,%s,,,,,,,,,,,%.2f", timestep.getYear(), "unhandled", unhandledArea));
			outputFile.newLine();

			outputFile.close();
		} catch (IOException e) {
			LogWriter.print(e);
		}
	}


	private void writeDemandFile(Timestep timestep) {
		try {
			BufferedWriter outputFile = FileWriterHelper.getFileWriter(timestep, ModelConfig.DEMAND_OUTPUT_FILE, "Year,Commodity,Amount (Mt)");

			for (CommodityType comm : CommodityType.getAllFoodItems()) {
				double demandAmount = 0;

				for (AbstractCountryAgent country : countryAgents.getAll()) {
					Map<CommodityType, Double> demands = country.getCurrentProjectedDemand();
					if (demands == null) {
						LogWriter.printlnError(country.getCountry() + " " + comm);
					}
					
					Double d = demands.get(comm);
					if (d != null) {
						demandAmount += d.doubleValue();
						LogWriter.println(String.format("%s,%s,%.4f", country.getCountry(), comm.getGamsName(), d));
					}
				}
				StringBuffer sbData = new StringBuffer();
				sbData.append(String.format("%d,%s", timestep.getYear(), comm.getGamsName()));
				sbData.append(String.format(",%.1f", demandAmount));

				LogWriter.println("Global demand " + timestep.getYear() + " " + comm.getGamsName() + " " + demandAmount + "\n");
				outputFile.write(sbData.toString());
				outputFile.newLine();
			}
			outputFile.close();
		} catch (IOException e) {
			LogWriter.print(e);
		}
	}

	private void writeDomesticProductionFile(Timestep timestep) {
		try {
			StringBuffer sbHeadings = new StringBuffer("Year, Country, Crop, Area, Production, Production_cost, Import_price, Export_price, Net_imports, Net_import_cost, Prod_shock, Rum_feed_amount, Mon_feed_amount");
			BufferedWriter outputFile = FileWriterHelper.getFileWriter(timestep, ModelConfig.DOMESTIC_OUTPUT_FILE, sbHeadings.toString());

			for (CropType crop : CropType.getAllItems()) {
				for (AbstractCountryAgent country : countryAgents.getAll()) {
					
					Map<CropType, CropUsageData> cropUsageAllCrops = country.getCropUsageData();
					CropUsageData cropUsage = cropUsageAllCrops.get(crop);
					
					if (cropUsage == null) 
						continue;
					
					double prodCosts = cropUsage.getProdCost();
					double prod = cropUsage.getProductionExpected();
					double prodShock = cropUsage.getProductionShock();
					double area = cropUsage.getArea();
					double rumFeedAmount = cropUsage.getRuminantFeed();
					double monFeedAmount = cropUsage.getMonogastricFeed();
					
					double importPrice = 0;
					double exportPrice = 0 ;
					double netImports = 0;
					double netImportCost = 0;

					if (crop.isImportedCrop()) {
						CountryPrice px = country.getCurrentCountryPrices().get(crop);
						importPrice = px.getImportPrice(); 
						exportPrice = px.getExportPrice();
						netImports = cropUsage.getNetImportsExpected();	//this isn't accounting for transport losses in exports 
						netImportCost = cropUsage.getNetImportCostExpected();
					}
	
					StringBuffer sbData = new StringBuffer();
					sbData.append(String.format("%d,%s,%s", timestep.getYear(), country.getCountry(), crop.getGamsName()));
					sbData.append(String.format(",%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f,%.3f", area, prod, prodCosts, importPrice, exportPrice, netImports, netImportCost, prodShock, rumFeedAmount, monFeedAmount));

					outputFile.write(sbData.toString());
					outputFile.newLine();
				}
			}
			outputFile.close();

		} catch (IOException e) {
			LogWriter.print(e);
		}
	}

	private void writeCountryDemandFile(Timestep timestep){

		try {
			StringBuffer sbHeadings = new StringBuffer("Year, Country, Commodity, Demand, BioenergyDemand");
			BufferedWriter outputFile = FileWriterHelper.getFileWriter(timestep, ModelConfig.COUNTRY_DEMAND_FILE, sbHeadings.toString());

			for (AbstractCountryAgent country : countryAgents.getAll()) {
				for (CommodityType commodity : CommodityType.getAllFoodItems()) {
					double bioenergyDemand = demandManager.getFirstGenBioenergyDemand(country.getCountry(), timestep.getYear(), commodity);
					double  demand =  country.getCurrentProjectedDemand().get(commodity);

					StringBuffer sbData = new StringBuffer();
					sbData.append(String.format("%d,%s,%s", timestep.getYear(), country.getCountry(), commodity.getGamsName()));
					sbData.append(String.format(",%.3f,%.3f", demand, bioenergyDemand));

					outputFile.write(sbData.toString());
					outputFile.newLine();
				}
			}
			outputFile.close();

		} catch (IOException e) {
			LogWriter.print(e);
		}
	}
	
	private void writeAnimalNumber(Timestep timestep) {
		try {
			StringBuffer sbHeadings = new StringBuffer("Year,Country,FAOItem,Heads(M)");
			BufferedWriter outputFile = FileWriterHelper.getFileWriter(timestep, ModelConfig.ANIMAL_NUMBERS_OUTPUT_FILE, sbHeadings.toString());
		
			for (AbstractCountryAgent country : countryAgents.getAll()) {
				Map<CropType, CropUsageData> cropUsageAllCrops = country.getCropUsageData();
				for (CropType crop : CropType.getMeatTypes()) {
					CropUsageData cropusage = cropUsageAllCrops.get(crop);
					if (cropusage == null)
						continue;
					double prod = cropusage.getProductionExpected();
					
					Map<String, Double> animalRates = animalRateManager.getAnimalRates(country.getCountry(), crop);
					for (Entry<String, Double> entry : animalRates.entrySet()) {
						StringBuffer sbData = new StringBuffer();
						double animalNum = prod * entry.getValue();
						sbData.append(String.format("%d,%s,%s,%.4f", timestep.getYear(), country.getCountry(), entry.getKey(), animalNum));		
						outputFile.write(sbData.toString());
						outputFile.newLine();
					}
			}
		}
		outputFile.close();
			
		} catch (IOException e) {
			LogWriter.print(e);
		}
	}
	
	private void outputTimestepResults(Timestep timestep, RasterSet<LandUseItem> landUseRaster) {

		writeLandCoverFile(timestep, landUseRaster);
		writeGlobalMarketFile(timestep);
		writeDemandFile(timestep);
		writeDomesticProductionFile(timestep);
		writeCountryDemandFile(timestep);
		writeGlobalFoodBalanceSheet(timestep, landUseRaster);
		writeAnimalNumber(timestep);

		if (ModelConfig.OUTPUT_FOR_LPJG) {
			for (int outputYear : timestep.getYearsFromLast()) {
				LogWriter.println("Outputing Year: " + outputYear);
				RasterSet<LandUseItem> landUseToOutput = null;

				if (outputYear == timestep.getYear()) {
					landUseToOutput = landUseRaster;
				} 
				else if (ModelConfig.INTERPOLATE_OUTPUT_YEARS) {
					InterpolatingRasterSet<LandUseItem> intermediateLandUse = new InterpolatingRasterSet<LandUseItem>( landUseRaster.getHeaderDetails()) {
						private static final long serialVersionUID = 1306045141011047760L;
						protected LandUseItem createRasterData() {
							return new LandUseItem();
						}
					};
					intermediateLandUse.setup(globalLandUseRaster, landUseRaster, timestep.getPreviousTimestep().getYear(), timestep.getYear(), outputYear);
					landUseToOutput = intermediateLandUse;
				}

				if (landUseToOutput != null) {
					LpjgOutputer lpjOutputer = new LpjgOutputer(outputYear, landUseToOutput);
					lpjOutputer.writeOutput();
				}
			}
			outputWaterAvailablity(timestep, currentIrrigationData);  // uses the year directory structure created above
		}

		if (ModelConfig.IS_CALIBRATION_RUN) {
			serializeLandUse(landUseRaster);
			countryAgents.serializeCropUsageForAll();
			internationalMarket.serializeGlobalPrices();
		}

		if (timestep.isInitialTimestep() && ModelConfig.GENERATE_NEW_YIELD_CLUSTERS)
			outputClusters(clusterIdRaster);

		// Output LandUses to tabular file, for analysis (perhaps)
		LogWriter.println("Outputing land uses Year: " + timestep.getYear());
		LandUseOutputer landuseOutputer = new LandUseOutputer(timestep.getYear(), landUseRaster);
		landuseOutputer.writeOutput();

		// don't really need this a LPJ outputs have same data, although in a slightly different format
		// outputLandCover(timestep.getYear(), landUseRaster, LandCoverType.CROPLAND);
		// outputLandCover(timestep.getYear(), landUseRaster, LandCoverType.PASTURE);
	}

	private void outputWaterAvailablity(Timestep timestep, IrrigationRasterSet irrigiationRS) {
		new RasterOutputer<Double, IrrigationItem>(irrigiationRS, ModelConfig.OUTPUT_DIR + File.separator + timestep.getYear() + File.separator + "IrrigConstraint.asc") {
			@Override
			public Double getValue(RasterKey location) {
				IrrigationItem item = results.get(location);
				if (item == null)
					return null;

				return item.getIrrigConstraint();
			}
		}.writeOutput();
	}

	private void outputClusters(RasterSet<IntegerRasterItem> landUseRaster) {
		new RasterOutputer<Integer, IntegerRasterItem>(landUseRaster, ModelConfig.CLUSTERED_YIELD_FILE) {
			@Override
			public Integer getValue(RasterKey location) {
				IntegerRasterItem item = results.get(location);
				if (item == null)
					return null;

				return item.getInt();
			}
		}.writeOutput();
	}

	public RasterSet<IntegerRasterItem> getClusterRaster() {
		RasterSet<IntegerRasterItem> clusters = new RasterSet<IntegerRasterItem>(desiredProjection) {
			private static final long serialVersionUID = 2467452274591854417L;

			@Override
			protected IntegerRasterItem createRasterData() {
				return new IntegerRasterItem(0);
			}
		};

		IntegerRasterReader clusterReader = new IntegerRasterReader(clusters);
		clusterReader.getRasterDataFromFile(ModelConfig.CLUSTERED_YIELD_FILE);
		return clusters;
	}

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

	public void createCountryAgents(Collection<CompositeCountry> countryGrouping) {
		countryAgents = new CountryAgentManager(compositeCountryManager, demandManager, countryBoundaryRaster, internationalMarket, clusterIdRaster, globalLandUseRaster);
		Map<CompositeCountry, Map<CropType, CropUsageData>> cropUsageDataMap = getInitialCropUsageData();
		RasterSet<LandUseItem> initLU = getInitialLandUse();

		for (CompositeCountry cc : countryGrouping) {
			countryAgents.addForCountry(cc, cropUsageDataMap, initLU);
		}
	}

	private RasterSet<LandUseItem> getInitialLandUse() {
		RasterSet<LandUseItem> initialLU;
		if (ModelConfig.IS_CALIBRATION_RUN)
			initialLU = getLandUseFromBaseline();
		else
			initialLU = deserializeLandUse();
		
		new ProtectedAreasReader(initialLU).getRasterDataFromFile(ModelConfig.PROTECTED_AREAS_FILE);
		return initialLU;
	}
	
	private Map<CompositeCountry, Map<CropType, CropUsageData>> getInitialCropUsageData() {
		Map<CompositeCountry, Map<CropType, CropUsageData>> cropUsageDataMap;
		if (ModelConfig.IS_CALIBRATION_RUN)
			cropUsageDataMap = new CropUsageReader(compositeCountryManager).getCommodityData();
		else
			cropUsageDataMap = deserializeCropUsage();

		return cropUsageDataMap;
	}

	private void serializeLandUse(RasterSet<LandUseItem> landUseRaster) {
		try {
			LogWriter.println("Starting serializing LandUse to " + ModelConfig.SERIALIZED_LAND_USE_FILE);
			FileOutputStream fileOut = new FileOutputStream(ModelConfig.SERIALIZED_LAND_USE_FILE);
			ObjectOutputStream out = new ObjectOutputStream(fileOut);
			out.writeObject(landUseRaster);
			out.close();
			fileOut.close();
			LogWriter.println("Serialized data is saved");
		} catch (IOException i) {
			i.printStackTrace();
		}
	}

	@SuppressWarnings("unchecked")
	private RasterSet<LandUseItem> deserializeLandUse() {
		try {
			RasterSet<LandUseItem> initLU;
			FileInputStream fileIn = new FileInputStream(ModelConfig.SERIALIZED_LAND_USE_FILE);
			ObjectInputStream in = new ObjectInputStream(fileIn);
			initLU = (RasterSet<LandUseItem>) in.readObject();
			in.close();
			fileIn.close();
			LogWriter.println("Deserialized " + ModelConfig.SERIALIZED_LAND_USE_FILE);
			return initLU;
		} catch (IOException i) {
			LogWriter.printlnError("Problem deserializing " + ModelConfig.SERIALIZED_LAND_USE_FILE);
			LogWriter.print(i);
			return null;
		} catch (ClassNotFoundException c) {
			LogWriter.printlnError("RasterSet<LandUseItem> class not found");
			c.printStackTrace();
			return null;
		}
	}
	
	@SuppressWarnings("unchecked")
	private Map<CompositeCountry, Map<CropType, CropUsageData>> deserializeCropUsage() {
		try {
			Map<CompositeCountry, Map<CropType, CropUsageData>> initCropUsage;
			FileInputStream fileIn = new FileInputStream(ModelConfig.SERIALIZED_CROP_USAGE_FILE);
			ObjectInputStream in = new ObjectInputStream(fileIn);
			initCropUsage = (Map<CompositeCountry, Map<CropType, CropUsageData>>) in.readObject();
			in.close();
			fileIn.close();
			LogWriter.println("Deserialized " + ModelConfig.SERIALIZED_CROP_USAGE_FILE);
			return initCropUsage;
		} catch (IOException i) {
			LogWriter.printlnError("Problem deserializing " + ModelConfig.SERIALIZED_CROP_USAGE_FILE);
			LogWriter.print(i);
			return null;
		} catch (ClassNotFoundException c) {
			LogWriter.printlnError("Map<CompositeCountry, Map<CropType, CropUsageData>> not found");
			c.printStackTrace();
			return null;
		}
	}
	

	/** this is if we are starting from Hurtt of other initial land covers (so we don't have land uses and intensity data) */
	private RasterSet<LandUseItem> getLandUseFromBaseline() {
		RasterSet<LandCoverItem> initialLC = new RasterSet<LandCoverItem>(desiredProjection) {
			private static final long serialVersionUID = 4642550777741425501L;

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

		new MaxCropAreaReader(initialLC).getRasterDataFromFile(ModelConfig.HIGH_SLOPE_AREAS_FILE);
		new LandCoverReader(initialLC).getRasterDataFromFile(ModelConfig.INITAL_LAND_COVER_FILE);

		RasterSet<LandUseItem> landUseRaster = new RasterSet<LandUseItem>(initialLC.getHeaderDetails());

		for (Map.Entry<RasterKey, LandCoverItem> entry : initialLC.entrySet())
			landUseRaster.put(entry.getKey(), new LandUseItem(entry.getValue()));

		return landUseRaster;
	}

	private YieldRaster getYieldSurfaces(Timestep timestep) {
		return lpjYieldReader.getRasterData(timestep);
	}

	/** Get irrigation data that does not change with time, should only be called once */
	private IrrigationRasterSet getFixedIrrigationData() {
		IrrigationRasterSet fixedIrrigData = new IrrigationRasterSet(desiredProjection, new FPUManager(desiredProjection));
		new IrrigiationCostReader(fixedIrrigData).getRasterDataFromFile(ModelConfig.IRRIGATION_COST_FILE);
		new IrrigationConstraintReader(fixedIrrigData).getRasterDataFromFile(ModelConfig.IRRIGATION_CONSTRAINT_FILE);
		
		String baseTimestepRootDir = Timestep.getYearSubDir(ModelConfig.YIELD_DIR, ModelConfig.ELLIOTT_BASEYEAR); // needs to be Elliott base timestep
		new RunOffReader(fixedIrrigData, true).getRasterDataFromFile(baseTimestepRootDir + File.separator + ModelConfig.IRRIG_RUNOFF_FILE);
		
		fixedIrrigData.calcIrrigConstraintOffsets();  // should have everything we need to calc offset between Elliott and LPJ data
		return fixedIrrigData;
	}

	/** Ugly in situ update of currentIrrigationData, better if IrrigationRasterSets were handled more immutably */
	private void getUpdateIrrigationData(Timestep timestep, YieldRaster yieldSurfaces) {
		String rootDir = timestep.getYearSubDir(ModelConfig.YIELD_DIR);
		
		new IrrigationMaxAmountReader(currentIrrigationData).getRasterDataFromFile(rootDir + File.separator + ModelConfig.IRRIG_MAX_WATER_FILENAME);

		if (!ModelConfig.USE_BLUE_WATER_FILE_IRRIG_CONSTRAINT) {
			new RunOffReader(currentIrrigationData, false).getRasterDataFromFile(rootDir + File.separator + ModelConfig.IRRIG_RUNOFF_FILE);
			currentIrrigationData.updateIrrigConstraints(timestep);
		}
	}	
}