/*
 *  Growth under dictatorship analysis.
 *
 */
 
/* Kenneth L. Simons, 2000.  This program was written March 2000, based 
on an earlier version I created in May 1999.  It uses much-modified versions 
of two  earlier programs, arkina and arkpathp.  Notations at the top of those
programs indicate that they were transferred from Jody Overland to
Madelyn Hicks at WINPO.  Arkina determined optimal consumption policy
over time for a country ruled by a dictator who is able to control
consumption but who fears being overthrown.  Arkpathp simulated the economy 
given optimal policy.  This program carries out sensitivity analyses in 
which model parameters are varied and for each set of parameters, optimal 
policy is computed and then checks are carried out to see (a) whether the 
initial value of human capital has a bifurcation point above which human 
capital grows and below which it falls, (b) whether consumption policy as a 
function of human capital has one of several possible shapes, and (c) whether 
the dictator's value (summed discounted utility stream) as a funtion of human 
capital has one of several possible shapes.

If future users need more speed, the program might be sped up -- at the cost
of reprogramming -- using the alternating sweep method (see for example Judd's
book Numerical Methods in Economics).
*/

/* Specify the target machine below as either TARGETMACHINEISMACPPC or TARGETMACHINEISWIN . */
#define TARGETMACHINEISMACPPC 1

/* Libraries used. */
#include <stdio.h>
#include <iostream>
#include <iomanip>
#include <stdlib.h>
#include <time.h>
//#include <extmath>
//#include <valarray>
#include <cmath>
#include <new>
#include <console.h>
#include <sioux.h>
#include <fstream>

using namespace std;  //introduces namespace std

/* Functions below, except those that pertain to specific classes. */
void varyParametersAndTestPolicies(void);
void maybeTerminateByUser(void);
void newFailed();
void testWritingToDiskFile(int n);

/* Constants. */
const maxKGridSize = 6002; //5010
const maxHGridSize = 21;
const dictator = 0;
const socialPlanner = 1;
const maxAllowableTimePeriodsForKPath = 200;

/* Global variable declarations. */
char *gSpareBlockPtr = 0;
long int parameterSetNumber;

/* Global variables used only by findOptimalHs(). */
// These are declared as globals to keep down the amount of stack memory needed for local variables in that routine; otherwise it would exceed the max. allowable.
long double lowHRangeToTry[maxKGridSize], highHRangeToTry[maxKGridSize], stepOfHRangeToTry[maxKGridSize];
long double discountedUtilityByH[maxHGridSize][maxKGridSize], highestUtilitySoFar[maxKGridSize];
long int highestUtilitySoFarHGridPoint[maxKGridSize];

/* Global variables used only by computePresentAndFutureUtilityIfInPower(). */
// This is declared here for the same reason as given above.
long double capitalNext[maxKGridSize];


/* Class modelParameters */
class modelParameters
{
	protected:
	long double initK, alpha, cpform, eta, theta, delta, gamma, sigma, lambda, uMin, beta, rho, nu;
	long double logLambda, oneMinusDelta, oneMinusSigma, wageOfK, utilityConstantMultiplierTerm;
	
	public:
	modelParameters(long double pInitK, long double pAlpha, long double pCpform, long double pEta, long double pTheta, long double pDelta, long double pGamma, long double pSigma, long double pLambda, long double pUMin, long double pBeta, long double pRho, long double pNu);
};

modelParameters::modelParameters(long double pInitK = 20.0, long double pAlpha = 0.33333333, long double pCpform = 1, long double pEta = 0.002, long double pTheta = -0.02, long double pDelta = 0.07, long double pGamma = 1.0, long double pSigma = 1.0, long double pLambda = 1.0, long double pUMin = 0.0, long double pBeta = 0.95, long double pRho = 0.02, long double pNu = .005)
{
	initK = pInitK;
	alpha = pAlpha;
	cpform = pCpform;
	eta = pEta;
	theta = pTheta;
	delta = pDelta;
	gamma = pGamma;
	sigma = pSigma;
	lambda = pLambda;
	uMin = pUMin;
	beta = pBeta;
	rho = pRho;
	nu = pNu;
	logLambda = logl(lambda);
	oneMinusDelta = 1.0 - delta;
	oneMinusSigma = 1.0 - sigma;
	wageOfK = nu * alpha * powl( nu*(1 - alpha)/rho, (1 - alpha)/alpha );
	utilityConstantMultiplierTerm = powl(lambda, oneMinusSigma) / oneMinusSigma;
	
	if (gamma != 1.0) cout << "Need to modify program for gamma other than 1!  Using gamma = 1 and proceeding....\n\n";
}

/* Class optimalPolicy */
class optimalPolicy : public modelParameters
{
	protected:
	long int maxT, kGridPoints, kGridSteps, gridStepsPerOrderOfMagnitudeInK, kGridStepsMinusOne;
	long double searchHTolerance, finalHTolerance, finalVTolerance, minUtilityDiffForConvergenceAsFractionOfMaxUtility, kGridMin, kGridMax;
	bool kBoundedBelow, kBoundedAbove;
	long int hGridPoints, hGridSteps;
	long int numberOfExtraChecksOnConvergence, whatToDoIfNotConverging;
	bool useMoreExactVComputation;
	long double k[maxKGridSize], vwSum[2][maxKGridSize], prevVWSum[2][maxKGridSize];
	long double kForInterpolation[maxKGridSize], kGridMinForInterpolation, kGridMaxForInterpolation, kDiffFirstGridStepForInterpolation, kDiffLastGridStepForInterpolation;
	long int kGridStepsForInterpolation, kGridStepsMinusOneForInterpolation;
	long double newHDict[maxKGridSize], prevHDict[maxKGridSize], newHSocialPlanner[maxKGridSize], prevHSocialPlanner[maxKGridSize];
	long double productionInCountry[maxKGridSize], h[maxKGridSize], consumption[maxKGridSize], qNext[maxKGridSize], utilityIfInPower[maxKGridSize], oneMinusDeltaTimesK[maxKGridSize], presentAndFutureUtilityIfInPower[maxKGridSize];
	long double uMinTimesBetaPowerSummation;
	long double savedK[maxKGridSize], savedProductionInCountry[maxKGridSize], savedOneMinusDeltaTimesK[maxKGridSize], savedH[maxKGridSize], probInPower[maxKGridSize], prevVOrWSum[maxKGridSize];
	long double kPath[maxAllowableTimePeriodsForKPath];
	long double timePeriodsForKPath, bifurcationSensitivity, bifurcationAnomalyStepFrac;
	long double bifurcationPoint, bifPtCheckLow, bifPtCheckHigh, dictatorCutThroughSPPolicyPoint, spMinFractionalGrowthPerAnnum, lowestHSP, lowestHDict, highestHSPMinusHDict;
	long double kBiggestDropH1, kBiggestDropH10, kBiggestDropH50, kBiggestDropH250, biggestDropH1, biggestDropH10, biggestDropH50, biggestDropH250, qtplus1AtStartBiggestDropH1, qtplus1AtStartBiggestDropH10, qtplus1AtStartBiggestDropH50, qtplus1AtStartBiggestDropH250, qtplus1AtEndBiggestDropH1, qtplus1AtEndBiggestDropH10, qtplus1AtEndBiggestDropH50, qtplus1AtEndBiggestDropH250;
	int dictValFnMonotoniclyIncreasing, spValFnMonotoniclyIncreasing, dictValFnConcave, spValFnConcave, dictPolFnShape, spPolFnShape;
	long double fracDiscreteGrowthPerAnnumDictLow, fracDiscreteGrowthPerAnnumDictMedLow, fracDiscreteGrowthPerAnnumDictMed, fracDiscreteGrowthPerAnnumDictMedHigh, fracDiscreteGrowthPerAnnumDictHigh, fracDiscreteGrowthPerAnnumSPLow, fracDiscreteGrowthPerAnnumSPMedLow, fracDiscreteGrowthPerAnnumSPMed, fracDiscreteGrowthPerAnnumSPMedHigh, fracDiscreteGrowthPerAnnumSPHigh;
	long double fracDiscreteGrowthPerAnnumDict[maxKGridSize], fracDiscreteGrowthPerAnnumSP[maxKGridSize];
	long int minGridPointForGrowthRates, maxGridPointForGrowthRates;
	bool readInPastAnalyses, writeHsToFile;
	long double readInK[maxKGridSize], readInHDict[maxKGridSize], readInHSocialPlanner[maxKGridSize], readInVWSum[2][maxKGridSize];
	long double sumHDict[maxKGridSize];
	long int sumHDictNPeriodsSummed;
	long int useMeanPerturbationAtTime, useMinMaxPerturbationCode;
	
	public:
	optimalPolicy(long int pUseMeanPerturbationAtTime, long int pUseMinMaxPerturbationCode, bool pUseMoreExactVComputation, bool pReadInPastAnalyses, long int pWhatToDoIfNotConverging, bool pWriteConvergenceInfo, long int pMaxT, long double pKGridMin, long double pKGridMax, long double pInitKMin, long double pInitKMax, long double pAlpha, long double pCpform, long double pEta, long double pTheta, long double pDelta, long double pGamma, long double pSigma, long double pLambda, long double pUMin, long double pBeta, long double pRho, long double pNu);
	// ~optimalPolicy();
	bool findOptimalPolicy(void);
	void analyzeResults(void);
	void computeMoreExactPresentAndFutureUtilityIfInPower(char dictatorOrPlanner);
	void writeResultsToDiskFile(char resultsIdentifier);
	void writeFullGrowthResultsToDiskFile(char resultsIdentifier);
	void writeConvergedHFunctionsToDiskFile(bool successfulConvergence);
	void writeVWSumFunctionsToDiskFile(bool successfulConvergence);
	short int readHFunctionsFromDiskFile(long int parameterSetNumberToLoad, bool errorIfCantRead = true, bool loadValFunctions = false);
	
	private:
	long int finalPeriodPolicy(long int normalStartT);
	inline long int kGridPointForKLessThanOrEqualToValue(const long double &kValue);
	long double interpolateFnOnKGrid(long double &kValue, long double *f);
	void computeProduction(void);
	void computeQNext(void);
	void computeConsumption(void);
	void computeUtilityIfInPower(void);
	void findOptimalHs(char dictatorOrPlanner);
	void computePresentAndFutureUtilityIfInPower(char dictatorOrPlanner, long double *presentAndFutureUtilityIfInPower);
	void displayProgress(long int t, long double largestHDiffDict, long double largestHDiffSP, long double meanHDiffDict, long double meanHDiffSP, long double largestVDiffDict, long double largestVDiffSP, long double meanVDiffDict, long double meanVDiffSP);
	void writeConvergenceInfoToDiskFile(bool closeInsteadOfWriting, long int t);
	void computeKPath(long double initK);
	int kDirection(void);
	int computeKPathAndGetDirection(long double initK);
	void findBifurcationPoint(long double lowRange, long double highRange);
	void checkBifurcation(long double lowRange, long double highRange);
	void findDictatorCutThroughSPPolicyPoint(long double lowRange, long double highRange);
	void findSPMinFractionalGrowthPerAnnum(long double lowRange, long double highRange);
	void findLowestHInfo(long double lowRange, long double highRange);
	void findDictFastestHDropInfo(long double lowRange, long double highRange);
	void assessValueFunctionShape(long double lowRange, long double highRange);
	void assessPolicyFunctionShape(long double lowRange, long double highRange);
	void getGrowthRates(long double lowRange, long double highRange);
};

optimalPolicy::optimalPolicy(long int pUseMeanPerturbationAtTime, long int pUseMinMaxPerturbationCode, bool pUseMoreExactVComputation, bool pReadInPastAnalyses, long int pWhatToDoIfNotConverging, bool pWriteConvergenceInfo = false, long int pMaxT = 20, long double pKGridMin = -1.0, long double pKGridMax = -1.0, long double pInitKMin = 1.5, long double pInitKMax = 250.0, long double pAlpha = 0.33333333, long double pCpform = 1, long double pEta = 0.002, long double pTheta = -0.02, long double pDelta = 0.07, long double pGamma = 1.0, long double pSigma = 1.0, long double pLambda = 1.0, long double pUMin = 0.0, long double pBeta = 0.95, long double pRho = 0.02, long double pNu = .005) 
: modelParameters(pInitKMin, pAlpha, pCpform, pEta, pTheta, pDelta, pGamma, pSigma, pLambda, pUMin, pBeta, pRho, pNu)
{
	long int i;
	long double ordersOfMagnitudeDiffInK, klnmin, klnmax, klndiff, klndiffOverKGridSteps;
	
	// Determine the number of time periods T over which to optimize (0 or less implies T=inifinity).
	maxT = pMaxT;
	
	// The following parameters affect computational accuracy.
	searchHTolerance = 0.000001; // 0.00001
	finalHTolerance = 0.00001; // 0.0001
	finalVTolerance = 0.00001;
	gridStepsPerOrderOfMagnitudeInK = 500; // 200
	hGridPoints = 21;
	useMoreExactVComputation = pUseMoreExactVComputation;
	minUtilityDiffForConvergenceAsFractionOfMaxUtility = 0.001; // Only relevant when useMoreExactVComputation == true, for the "more exact" computation of present and future V (W).
	numberOfExtraChecksOnConvergence = 50; // 20 // Satisfy H convergence criterion this many extra times in a row before declaring convergence, just to be safe.
	whatToDoIfNotConverging = pWhatToDoIfNotConverging;
	writeHsToFile = pWriteConvergenceInfo; // Set this to true to save the current h for each k at each time step of the optimal policy search.
	
	// Determine range and number of grid points to use for the k-grid.
	kBoundedBelow = pKGridMin != -1.0;
	if (kBoundedBelow) kGridMin = pKGridMin;
	else
	{
		if (maxT > 0) kGridMin = pInitKMin * powl(oneMinusDelta, maxT - 1) * 0.99;
		else kGridMin = 0.0001;
	}
	kBoundedAbove = pKGridMax != -1.0;
	if (kBoundedAbove) kGridMax = pKGridMax;
	else
	{
		if (maxT > 0) kGridMax = pInitKMax * powl( oneMinusDelta + wageOfK, maxT - 1 ) * 1.01;
		else kGridMax = 1000000.0;
	}
	ordersOfMagnitudeDiffInK = log10l(kGridMax) - log10l(kGridMin);
	kGridPoints = 1 + ceil( ordersOfMagnitudeDiffInK * (long double) gridStepsPerOrderOfMagnitudeInK );
	if (kGridPoints > maxKGridSize) {
		cout << "Maximal k-grid size is too low, " << maxKGridSize << " when " << kGridPoints << " is needed.  Exiting....\n\n";
		exit(EXIT_FAILURE);
	}
	kGridSteps = kGridPoints - 1;
	kGridStepsMinusOne = kGridSteps - 1;
	
	// Prepare for use of the h-grid.
	hGridSteps = hGridPoints - 1;
	if (hGridPoints > maxHGridSize) {
		cout << "Maximal h-grid size is too low, " << maxHGridSize << " when " << hGridPoints << " is needed.  Exiting....\n\n";
		exit(EXIT_FAILURE);
	}
	
	// Prepare variables.
	klnmin = logl(kGridMin);
	klnmax = logl(kGridMax);
	klndiff = klnmax - klnmin;
	klndiffOverKGridSteps = klndiff / kGridSteps;
	for (i = 0; i < kGridPoints; i++) k[i] = expl(klnmin + klndiffOverKGridSteps * i);
	for (i = 0; i < kGridPoints; i++) kForInterpolation[i] = k[i]; // The kForInterpolation array is used by the interpolation routine.  Note that k will differ from kForInterpolation during the routine computeMoreExactPresentAndFutureUtilityIfInPower, hence the need for the two separate variables.
	kDiffFirstGridStepForInterpolation = k[1] - k[0];
	kDiffLastGridStepForInterpolation = k[kGridSteps] - k[kGridStepsMinusOne];
	kGridMinForInterpolation = kGridMin;
	kGridMaxForInterpolation = kGridMax;
	kGridStepsForInterpolation = kGridSteps;
	kGridStepsMinusOneForInterpolation = kGridStepsMinusOne;
	
	// The following parameters affect analysis of results.
	timePeriodsForKPath = 20;
	if (timePeriodsForKPath > maxAllowableTimePeriodsForKPath) {
		cout << "maxAllowableTimePeriodsForKPath is too low, " << maxAllowableTimePeriodsForKPath << " when " << timePeriodsForKPath << " is needed.  Exiting....\n\n";
		exit(EXIT_FAILURE);
	}
	bifurcationSensitivity = 0.005;
	bifurcationAnomalyStepFrac = 0.1;
	
	readInPastAnalyses = pReadInPastAnalyses;
	useMeanPerturbationAtTime = pUseMeanPerturbationAtTime;
	useMinMaxPerturbationCode = pUseMinMaxPerturbationCode;
}

//optimalPolicy::~optimalPolicy()
//{
//}

bool optimalPolicy::findOptimalPolicy(void)
{
	long int i, t, tDiffForSummation;
	bool notDoneYet = true;
	long double betaPowerSummation, betaToThePowerI;
	bool allHConvergedNow, allVConvergedNow, successfulConvergence;
	long int nTimesConvergedPrev = 0;
	long double largestHDiffDict, largestHDiffSP, meanHDiffDict, meanHDiffSP, largestVDiffDict, largestVDiffSP, meanVDiffDict, meanVDiffSP;
	
	// Prepare to search for the optimal policy functions.
	computeProduction();  // This always comes out the same, so only need to compute it once.
	for (i=0; i<kGridPoints; i++) oneMinusDeltaTimesK[i] = oneMinusDelta * k[i]; // This always comes out the same, so only need to compute it once.
	if (cpform == 1.0L || cpform == 3.0L) computeQNext();  // If cpform == 1 or 3, this always comes out the same, so only need to compute it once.
	if (useMinMaxPerturbationCode != 1L) for (i=0; i<kGridPoints; i++) sumHDict[i] = 0.0; else for (i=0; i<kGridPoints; i++) sumHDict[i] = 1.0; sumHDictNPeriodsSummed = 0L;
	
	// Search for the optimal policy functions for the dictator and the social planner.
	successfulConvergence = true;  // Assume successful convergence on finishing, unless proven otherwise.
	t = finalPeriodPolicy(maxT);  // IF READ IN FROM FILE, NEED TO SET T APPROPRIATELY!!!
	while (notDoneYet)
	{
		// First create uMin times the sum of 1.0 + beta + beta^2 + beta^3 + ... with terms corresponding to t+1, ..., T.
		// This is the present and future utility if not in power, evaluated at time t+1.
		// For example if T = 4 and we are now in year 3 (t=3), this is uMin * 1.0; if we are now in year 2, this is uMin * (1.0 + beta).
		// It is used in computePresentAndFutureUtilityIfInPower().
		betaPowerSummation = 0.0L;
		betaToThePowerI = 1.0L;
		tDiffForSummation = maxT-(t+1);
		for (i=0; i<=tDiffForSummation; i++)
		{
			betaPowerSummation += betaToThePowerI;
			betaToThePowerI *= beta;
		}
		uMinTimesBetaPowerSummation = uMin * betaPowerSummation;
	
		// Now find the optimal policy by dictator and social planner for this time period.
		// Note: computeProduction would also need to happen here were it not possible to do it earlier.
		findOptimalHs(dictator);
		findOptimalHs(socialPlanner);
	
		// Has the analysis finished?
		if ( (whatToDoIfNotConverging > 0) && (t < (-whatToDoIfNotConverging)) ) {notDoneYet = false; successfulConvergence = false;}
		else if (maxT > 0) {notDoneYet = t > 1; largestHDiffDict = -1.0L; largestHDiffSP = -1.0L; meanHDiffDict = -1.0L; meanHDiffSP = -1.0L; largestVDiffDict = -1.0L; largestVDiffSP = -1.0L; meanVDiffDict = -1.0L; meanVDiffSP = -1.0L;}
		else /* maxT <= 0 */
		{
			// For the infinite time horizon case, determine whether the optimal policy function has converged.
			allHConvergedNow = true;
			/* for (i=0; i<kGridPoints; i++)
				if ( (abs(newHDict[i] - prevHDict[i]) > finalHTolerance) || (abs(newHSocialPlanner[i] - prevHSocialPlanner[i]) > finalHTolerance) )
					{allHConvergedNow = false; break;} */
			largestHDiffDict = 0.0L;
			for (i=0; i<kGridPoints; i++)
				if ( largestHDiffDict < (abs(newHDict[i] - prevHDict[i])) )
					largestHDiffDict = (abs(newHDict[i] - prevHDict[i]));
			largestHDiffSP = 0.0L;
			for (i=0; i<kGridPoints; i++)
				if ( largestHDiffSP < (abs(newHSocialPlanner[i] - prevHSocialPlanner[i])) )
					largestHDiffSP = (abs(newHSocialPlanner[i] - prevHSocialPlanner[i]));
			allHConvergedNow = (largestHDiffDict <= finalHTolerance) && (largestHDiffSP <= finalHTolerance);
			
			// For the infinite time horizon case, determine whether the value function has converged.
			allVConvergedNow = true;
			/* for (i=0; i<kGridPoints; i++)
				if ( (abs(newHDict[i] - prevHDict[i]) > finalHTolerance) || (abs(newHSocialPlanner[i] - prevHSocialPlanner[i]) > finalHTolerance) )
					{allHConvergedNow = false; break;} */
			largestVDiffDict = 0.0L;
			for (i=0; i<kGridPoints; i++)
				if ( largestVDiffDict < (abs(vwSum[dictator][i] - prevVWSum[dictator][i])) )
					largestVDiffDict = (abs(vwSum[dictator][i] - prevVWSum[dictator][i]));
			largestVDiffSP = 0.0L;
			for (i=0; i<kGridPoints; i++)
				if ( largestVDiffSP < (abs(vwSum[socialPlanner][i] - prevVWSum[socialPlanner][i])) )
					largestVDiffSP = (abs(vwSum[socialPlanner][i] - prevVWSum[socialPlanner][i]));
			allVConvergedNow = (largestVDiffDict <= finalVTolerance) && (largestVDiffSP <= finalVTolerance);
			
			notDoneYet = ( nTimesConvergedPrev < numberOfExtraChecksOnConvergence ) || ( !allHConvergedNow ) || ( !allVConvergedNow );
			
			// Compute mean changes for reporting purposes.
			meanHDiffDict = 0.0L;
			for (i=0; i<kGridPoints; i++)
				meanHDiffDict += abs(newHDict[i] - prevHDict[i]);
			meanHDiffDict /= (long double) kGridPoints;
			meanHDiffSP = 0.0L;
			for (i=0; i<kGridPoints; i++)
				meanHDiffSP += abs(newHSocialPlanner[i] - prevHSocialPlanner[i]);
			meanHDiffSP /= (long double) kGridPoints;
			meanVDiffDict = 0.0L;
			for (i=0; i<kGridPoints; i++)
				meanVDiffDict += abs(vwSum[dictator][i] - prevVWSum[dictator][i]);
			meanVDiffDict /= (long double) kGridPoints;
			meanVDiffSP = 0.0L;
			for (i=0; i<kGridPoints; i++)
				meanVDiffSP += abs(vwSum[socialPlanner][i] - prevVWSum[socialPlanner][i]);
			meanVDiffSP /= (long double) kGridPoints;
			
			// If trying a perturbation to the dictator's solution based on the mean (or min or max) of the h's over 50 periods, compute the mean (or min or max) and use it in place of the h's.
			if (useMeanPerturbationAtTime != 0L && (t - useMeanPerturbationAtTime) < 50L && (t - useMeanPerturbationAtTime) >= 0L)
			{
				if (useMinMaxPerturbationCode == 0)
				{
					for (i=0; i<kGridPoints; i++) sumHDict[i] += newHDict[i];
					sumHDictNPeriodsSummed++;
				}
				else if (useMinMaxPerturbationCode == 1) for (i=0; i<kGridPoints; i++) sumHDict[i] = min(sumHDict[i], newHDict[i]);
				else if (useMinMaxPerturbationCode == 2) for (i=0; i<kGridPoints; i++) sumHDict[i] = max(sumHDict[i], newHDict[i]);
				if (t == useMeanPerturbationAtTime)
				{
					if (useMinMaxPerturbationCode == 0) for (i=0; i<kGridPoints; i++) newHDict[i] = sumHDict[i] / (long double) sumHDictNPeriodsSummed;
					else for (i=0; i<kGridPoints; i++) newHDict[i] = sumHDict[i]; // This actually contains a min or max in this case, not a sum.
					computeMoreExactPresentAndFutureUtilityIfInPower(dictator);
				}
			}
			
			// For the infinite time horizon case, save information for use in the next iteration.
			if (allHConvergedNow) nTimesConvergedPrev++; else nTimesConvergedPrev = 0;
			for (i=0; i<kGridPoints; i++) prevHDict[i] = newHDict[i];
			for (i=0; i<kGridPoints; i++) prevHSocialPlanner[i] = newHSocialPlanner[i];
			for (i=0; i<kGridPoints; i++) prevVWSum[dictator][i] = vwSum[dictator][i];
			for (i=0; i<kGridPoints; i++) prevVWSum[socialPlanner][i] = vwSum[socialPlanner][i];
			
		}
		
		// Write results to a file if requested.
		if (writeHsToFile) writeConvergenceInfoToDiskFile(false, t);
		
		// Check if the user wants to terminate the program.  Display progress.
		displayProgress(t, largestHDiffDict, largestHDiffSP, meanHDiffDict, meanHDiffSP, largestVDiffDict, largestVDiffSP, meanVDiffDict, meanVDiffSP);
		maybeTerminateByUser();
		
		// Decrease t by 1 in preparation for the next step of the loop.
		t -= 1;
	}
	
	// Now newHDict[] and newHSocialPlanner[] tell the optimal values of h as functions of k for a dictator and a social planner.
	
	// Stop write to a file if requested.
	if (writeHsToFile) writeConvergenceInfoToDiskFile(true, t);
	
	return successfulConvergence;
}

long int optimalPolicy::finalPeriodPolicy(long int normalStartT)
{
	bool readValueFunctionsIfReading;
	short int readSuccessIndicator;
	long int i;
	long int startTToUse;
	
	if (readInPastAnalyses)
	{
		readValueFunctionsIfReading = false; // DEBUG
		readSuccessIndicator = readHFunctionsFromDiskFile(parameterSetNumber, false, readValueFunctionsIfReading);
		
		if ((readValueFunctionsIfReading == false) && (readSuccessIndicator < 99))
		{
			// This is only relevant if the value function files are not read.
			computeMoreExactPresentAndFutureUtilityIfInPower(dictator);
			computeMoreExactPresentAndFutureUtilityIfInPower(socialPlanner);
		}
		if (readSuccessIndicator < 99) startTToUse = -4001;  // ASSUMING PREVIOUS RUN WAS FOR 4000 PERIODS.
		else startTToUse = normalStartT;
	}
	else
	{
		startTToUse = normalStartT;
	}
	
	if (!readInPastAnalyses || (readSuccessIndicator==99))
	{
		for (i=0; i<kGridPoints; i++) h[i] = 1.0L;
		for (i=0; i<kGridPoints; i++) newHDict[i] = h[i];
		for (i=0; i<kGridPoints; i++) newHSocialPlanner[i] = h[i];
		if (!useMoreExactVComputation)
		{
			// computeProduction would happen here (or before choosing h[i]) if it were not possible to do it earlier.
			computeConsumption();
			computeUtilityIfInPower();
			for (i=0; i<kGridPoints; i++) vwSum[dictator][i] = utilityIfInPower[i]; // vwSum[0] pertains to the dictator.
			for (i=0; i<kGridPoints; i++) vwSum[socialPlanner][i] = utilityIfInPower[i]; // vwSum[1] pertains to the social planner.
		}
		else /* useMoreExactVComputation */
		{
			computeMoreExactPresentAndFutureUtilityIfInPower(dictator);
			computeMoreExactPresentAndFutureUtilityIfInPower(socialPlanner);
		}
	}
	
	for (i=0; i<kGridPoints; i++) prevHDict[i] = newHDict[i];
	for (i=0; i<kGridPoints; i++) prevHSocialPlanner[i] = newHSocialPlanner[i];
	for (i=0; i<kGridPoints; i++) prevVWSum[dictator][i] = vwSum[dictator][i];
	for (i=0; i<kGridPoints; i++) prevVWSum[socialPlanner][i] = vwSum[socialPlanner][i];
	
	return startTToUse;
}

/* If within the k-grid, return an array index indicating the last point in the k-grid <= kValue.  (Similar to upper_bound STL algorithm.) */
/* If the value exactly equals the k-grid upper bound, the answer returned is correct, but interpolation should not be used based on the next point in the grid since there is no next point in the grid. */
/* If outside the k-grid, you must ignore the answer returned. */
inline long int optimalPolicy::kGridPointForKLessThanOrEqualToValue(const long double &kValue)
{
	long int first; /* an iterator */
	long int len; /* a distance: distance(first, last) */
	long int half; /* a distance */
	long int middle; /* a forward iterator */
	first = 0;
	len = kGridStepsForInterpolation;
	while (len > 0)
	{
		half = len / 2;
		middle = first + half + 1;
		if (kValue < kForInterpolation[middle])
			len = half;
		else
		{
			first = middle;
			len = len - half - 1;
		}
	}
	return first;
}

/* optimalPolicy::interpolateFnOnKGridForAnArray(const long double &kValues, const long double &f[], long double &results[])
{
	long int i;
	for (i=0; i<kGridPoints; i++) results[i] = interpolateFnOnKGrid(kValues[i], f);
} */

long double optimalPolicy::interpolateFnOnKGrid(long double &kValue, long double *f)
{
	long int gridPointBelow, gridPointAbove;
	long double kBelow;
	if (kValue < kGridMinForInterpolation)
	{
		// Below the grid.
		return f[0] + ((f[1]-f[0])/kDiffFirstGridStepForInterpolation) * (kValue - kGridMinForInterpolation);
	}
	else
	{
		if (kValue < kGridMaxForInterpolation)
		{
			// kValue is on [kGridMinForInterpolation, kGridMaxForInterpolation).
			gridPointBelow = kGridPointForKLessThanOrEqualToValue(kValue);
			gridPointAbove = gridPointBelow + 1;
			kBelow = kForInterpolation[gridPointBelow];
			return f[gridPointBelow] + ((f[gridPointAbove] - f[gridPointBelow]) / (kForInterpolation[gridPointAbove] - kBelow)) * (kValue - kBelow);
		}
		else
		{
			if (kValue == kGridMaxForInterpolation)
			{
				// Last point in the grid.
				return f[kGridStepsForInterpolation];
			}
			else
			{
				// Above the grid.
				return f[kGridStepsForInterpolation] + ((f[kGridStepsForInterpolation]-f[kGridStepsMinusOneForInterpolation])/kDiffLastGridStepForInterpolation) * (kValue - kGridMaxForInterpolation);
			}
		}
	}
}

void optimalPolicy::computeProduction(void)
{
	long int i;
	for (i=0; i<kGridPoints; i++) productionInCountry[i] = k[i] * wageOfK;
}

void optimalPolicy::computeQNext(void)
{
	long int i;
	if (cpform == 1.0L)
		for (i=0; i<kGridPoints; i++) qNext[i] = min( expl(eta + theta * k[i]), 1.0L );
	else if (cpform == 2.0L)
		for (i=0; i<kGridPoints; i++) qNext[i] = min( expl(eta + theta * consumption[i]), 1.0L );
	else if (cpform == 3.0L)
		for (i=0; i<kGridPoints; i++) qNext[i] = min( expl(eta + theta * logl(k[i])), 1.0L );
	else /* cpform == 4.0L */
		for (i=0; i<kGridPoints; i++) qNext[i] = min( expl(eta + theta * logl(consumption[i])), 1.0L );
}

void optimalPolicy::computeConsumption(void)
{
	long int i;
	for (i=0; i<kGridPoints; i++) consumption[i] = productionInCountry[i] * h[i];
}

void optimalPolicy::computeUtilityIfInPower(void)
{
	long int i;
	if (sigma != 1.0L)
		for (i=0; i<kGridPoints; i++) utilityIfInPower[i] = utilityConstantMultiplierTerm * powl(consumption[i], oneMinusSigma);
	else /* sigma == 1.0 */
		for (i=0; i<kGridPoints; i++) utilityIfInPower[i] = logLambda + logl(consumption[i]);
}

void optimalPolicy::findOptimalHs(char dictatorOrPlanner)
{
	long int i, hCount;
	long double initialLowH;
	bool convergedInOptimalHSearch = false;
	
	// Initialize variables to begin search trying h from 0 to 1 inclusive for each k.
	initialLowH = searchHTolerance / 2.0L;
	for (i=0; i<kGridPoints; i++) lowHRangeToTry[i] = initialLowH;
	for (i=0; i<kGridPoints; i++) highHRangeToTry[i] = 1.0L;
	
	// Search for optimal values of h in successive iterations, getting closer and closer to the solution until converged.
	while (!convergedInOptimalHSearch)
	{
		// For each k, for each h in the h-grid, compute present and future utility using the past estimates of VSum / WSum rather than the actual V (W) that would result in the infinite-horizon or rolling fixed horizon cases.
		for (i=0; i<kGridPoints; i++) stepOfHRangeToTry[i] = (highHRangeToTry[i] - lowHRangeToTry[i]) / hGridSteps;
		for (hCount = 0; hCount<hGridPoints; hCount++)
		{
			for (i=0; i<kGridPoints; i++) h[i] = lowHRangeToTry[i] + stepOfHRangeToTry[i] * hCount;
			computePresentAndFutureUtilityIfInPower(dictatorOrPlanner, discountedUtilityByH[hCount]);  // Note: The computed utility is also called VSumNew or WSumNew in the documentation.
		}
		// Determine which point(s) in the h-grid have the highest present-and-future utility for each k.
		for (i=0; i<kGridPoints; i++) highestUtilitySoFar[i] = discountedUtilityByH[0][i];
		for (i=0; i<kGridPoints; i++) highestUtilitySoFarHGridPoint[i] = 0;
		for (hCount = 1; hCount<hGridPoints; hCount++)
			for (i=0; i<kGridPoints; i++)
				if (discountedUtilityByH[hCount][i] > highestUtilitySoFar[i])
					{highestUtilitySoFar[i] = discountedUtilityByH[hCount][i]; highestUtilitySoFarHGridPoint[i] = hCount;}
		// Determine the min & max values of h to try at the next iteration.
		for (i=0; i<kGridPoints; i++) highHRangeToTry[i] = min(lowHRangeToTry[i] + stepOfHRangeToTry[i] * (long double) min(hGridPoints, highestUtilitySoFarHGridPoint[i] + 1), 1.0L);
		for (i=0; i<kGridPoints; i++) lowHRangeToTry[i] = lowHRangeToTry[i] + stepOfHRangeToTry[i] * (long double) max(0L, highestUtilitySoFarHGridPoint[i] - 1);	// Note that the new lowHRangeToTry must be computed after highHRangeToTry, because proper computation of highHRangeToTry requires the old values of lowHRangeToTry.
		// Has a solution been found to within a plus or minus searchHTolerance, for all k?  If so, the search for the optimal h's has converged.
		// Note that it does not suffice to check this just for one k, because the different k's can have different (highHRangeToTry - lowHRangeToTry), because of the max and min functions above.
		convergedInOptimalHSearch = true;
		for (i=0; i<kGridPoints; i++) if ( (highHRangeToTry[i] - lowHRangeToTry[i]) > searchHTolerance ) {convergedInOptimalHSearch = false; break;}
		
		// Check if the user wants to terminate the program.
		maybeTerminateByUser();
	}
	
	// Save the results.
	for (i=0; i<kGridPoints; i++) h[i] = (lowHRangeToTry[i] + highHRangeToTry[i])/2;
	if (dictatorOrPlanner == dictator)
	{
		for (i=0; i<kGridPoints; i++) newHDict[i] = h[i];
		if (!useMoreExactVComputation) computePresentAndFutureUtilityIfInPower(dictator, vwSum[dictator]);
		else /* useMoreExactVComputation */ computeMoreExactPresentAndFutureUtilityIfInPower(dictator);
	}
	else /* dictatorOrPlanner = socialPlanner */
	{
		for (i=0; i<kGridPoints; i++) newHSocialPlanner[i] = h[i];
		if (!useMoreExactVComputation) computePresentAndFutureUtilityIfInPower(socialPlanner, vwSum[socialPlanner]);
		else /* useMoreExactVComputation */ computeMoreExactPresentAndFutureUtilityIfInPower(socialPlanner);
	}
}

// The following routine computes the actual discounted utility stream given a particular choice of h(k).
// It is used in order to implement a value function iteration method to speed up computation.
// This routine is called only by findOptimalHs.  NOTE that this routine has not been fully debugged,
// and should be carefully checked before using any results that depend upon it.
void optimalPolicy::computeMoreExactPresentAndFutureUtilityIfInPower(char dictatorOrPlanner)
{
	long int i, tCounter;
	long double discountMultiplier, minUtilityDiffForConvergence, maxUtil;
	bool exitFromLoop;
	
	// Save some of the arrays that were computed as a function of the basic values of k,
	// because below -- in contrast to the rest of the computations in the optimalPolicy 
	// class -- k will be altered.  The saved results will be returned to their places later.
	for (i=0; i<kGridPoints; i++) savedK[i] = k[i];
	for (i=0; i<kGridPoints; i++) savedProductionInCountry[i] = productionInCountry[i];
	for (i=0; i<kGridPoints; i++) savedOneMinusDeltaTimesK[i] = oneMinusDeltaTimesK[i];
	if (dictatorOrPlanner == dictator) for (i=0; i<kGridPoints; i++) savedH[i] = newHDict[i];
	else /* dictatorOrPlanner == socialPlanner */ for (i=0; i<kGridPoints; i++) savedH[i] = newHSocialPlanner[i];
	
	// Replace the VSum or WSum array with zeros, before adding the values from each point in time.
	for (i=0; i<kGridPoints; i++) vwSum[dictatorOrPlanner][i] = 0.0L;
	
	// Compute utility and growth now and at each point in the future, and sum the discounted utility, 
	// until the results converge for each initial k or a rolling fixed time horizon is reached.
	tCounter = 1;
	discountMultiplier = 1.0L; // Start out not having discounted yet.
	for (i=0; i<kGridPoints; i++) probInPower[i] = 1.0L; // Start out being in power for all k.
	exitFromLoop = false;
	while (!exitFromLoop)
	{
		computeProduction();
		for (i=0; i<kGridPoints; i++) oneMinusDeltaTimesK[i] = oneMinusDelta * k[i];
		if (tCounter == 1)
		{
			/* No need to interpolate the value of h to use now -- just use the value in newHDict or newHSocialPlanner for each k. */
			if (dictatorOrPlanner == dictator) for (i=0; i<kGridPoints; i++) h[i] = newHDict[i];
			else /* dictatorOrPlanner == socialPlanner */ for (i=0; i<kGridPoints; i++) h[i] = newHSocialPlanner[i];
		}
		else
		{
			// From each initial k, k has now changed.  Use the function h(k), held in the array savedH, to find the correct h for each new k.
			for (i=0; i<kGridPoints; i++) h[i] = interpolateFnOnKGrid(k[i], savedH);
		}
		computeConsumption();
		computeUtilityIfInPower();
		computeQNext();
		if (dictatorOrPlanner == dictator)
			for (i=0; i<kGridPoints; i++) vwSum[dictator][i] += discountMultiplier * (probInPower[i] * utilityIfInPower[i] + (1 - probInPower[i]) * uMin);
		else /* dictatorOrPlanner == socialPlanner */
			for (i=0; i<kGridPoints; i++) vwSum[socialPlanner][i] += discountMultiplier * utilityIfInPower[i];

		// Determine whether enough time periods have been computed.
		if (maxT > 0) exitFromLoop = tCounter >= maxT;
		else
		{
			if (tCounter > 1)
			{
				exitFromLoop = true;
				for (i=0; i<kGridPoints; i++) if (abs(prevVOrWSum[i] - vwSum[dictatorOrPlanner][i]) > minUtilityDiffForConvergence) {exitFromLoop = false; break;}
			}
			else {
				maxUtil = 0.0L;
				for (i=0; i<kGridPoints; i++) if (maxUtil < abs(utilityIfInPower[i])) maxUtil = utilityIfInPower[i];
				if (maxUtil == 0.0) maxUtil = 0.000001;
				minUtilityDiffForConvergence = maxUtil * minUtilityDiffForConvergenceAsFractionOfMaxUtility;
			}
			for (i=0; i<kGridPoints; i++) prevVOrWSum[i] = vwSum[dictatorOrPlanner][i];
		}
	
		// Compute capital in next time period and place it in the array k.
		for (i=0; i<kGridPoints; i++) k[i] = oneMinusDeltaTimesK[i] + (1.0L - h[i]) * productionInCountry[i];
		if (kBoundedBelow) for (i=0; i<kGridPoints; i++) k[i] = max(k[i], kGridMin);
		if (kBoundedAbove) for (i=0; i<kGridPoints; i++) k[i] = min(k[i], kGridMax);
			
		// Compute the discount multiplier for next time period.
		discountMultiplier *= beta;
		for (i=0; i<kGridPoints; i++) probInPower[i] *= (1 - qNext[i]);
		tCounter++;
		
		// Give an opportunity for the user to stop program execution.
		maybeTerminateByUser();
	}
	
	// Restore the saved arrays.
	for (i=0; i<kGridPoints; i++) k[i] = savedK[i];
	for (i=0; i<kGridPoints; i++) productionInCountry[i] = savedProductionInCountry[i];
	for (i=0; i<kGridPoints; i++) oneMinusDeltaTimesK[i] = savedOneMinusDeltaTimesK[i];
	if (dictatorOrPlanner == dictator) for (i=0; i<kGridPoints; i++) newHDict[i] = savedH[i];
	else /* dictatorOrPlanner == socialPlanner */ for (i=0; i<kGridPoints; i++) newHSocialPlanner[i] = savedH[i];
}

// The following routine assumes only that production, and qNext if cpform == 1 or 3, has been computed; and it uses the h array.
// This routine is called only by findOptimalHs.
void optimalPolicy::computePresentAndFutureUtilityIfInPower(char dictatorOrPlanner, long double *presentAndFutureUtilityIfInPower)
{
	long int i;
	
	computeConsumption();
	computeUtilityIfInPower();
	for (i=0; i<kGridPoints; i++) capitalNext[i] = oneMinusDeltaTimesK[i] + (1.0L - h[i]) * productionInCountry[i];
	if (kBoundedBelow) for (i=0; i<kGridPoints; i++) capitalNext[i] = max(capitalNext[i], kGridMin);
	if (kBoundedAbove) for (i=0; i<kGridPoints; i++) capitalNext[i] = min(capitalNext[i], kGridMax);
	if (cpform == 2 || cpform == 4) computeQNext();
	if (dictatorOrPlanner == dictator)
		for (i=0; i<kGridPoints; i++) presentAndFutureUtilityIfInPower[i] = utilityIfInPower[i] + beta * ( (1-qNext[i]) * interpolateFnOnKGrid(capitalNext[i], vwSum[dictator]) + qNext[i] * uMinTimesBetaPowerSummation );
	else /* dictatorOrPlanner == socialPlanner */
		for (i=0; i<kGridPoints; i++) presentAndFutureUtilityIfInPower[i] = utilityIfInPower[i] + beta * interpolateFnOnKGrid(capitalNext[i], vwSum[socialPlanner]);
}

void optimalPolicy::displayProgress(long int t, long double largestHDiffDict, long double largestHDiffSP, long double meanHDiffDict, long double meanHDiffSP, long double largestVDiffDict, long double largestVDiffSP, long double meanVDiffDict, long double meanVDiffSP)
{
	#ifdef TARGETMACHINEISMACPPC
	SIOUXSetTitle("\pDictators & Growth Computation");  // FOR MAC ONLY
	#endif
	clrscr();
	cout << "Growth under Dictatorship\n\n";
	#ifdef TARGETMACHINEISMACPPC
	cout << "To stop the program, press any key, then choose Quit from the File menu to quit.\n\n";
	#else
	cout << "To stop the program, choose Exit from the File menu.\n\n";
	#endif
	cout << "Parameter set: " << parameterSetNumber << "\n\n";
	cout << "Time: " << t << "\n";
	if ( (largestHDiffDict >= 0) && (t != maxT) )
	{
		cout << "\n";
		cout << "Largest h-difference for dictator:       " << largestHDiffDict << "\n";
		cout << "Largest h-difference for social planner: " << largestHDiffSP << "\n";
		cout << "Mean h-difference for dictator:          " << meanHDiffDict << "\n";
		cout << "Mean h-difference for social planner:    " << meanHDiffSP << "\n";
		cout << "\n";
		cout << "Largest v-difference for dictator:       " << largestVDiffDict << "\n";
		cout << "Largest v-difference for social planner: " << largestVDiffSP << "\n";
		cout << "Mean v-difference for dictator:          " << meanVDiffDict << "\n";
		cout << "Mean v-difference for social planner:    " << meanVDiffSP << "\n";
	}
}

// Create a disk file that, when transposed, has the data for dictator and social planner optimal policies at each t.
void optimalPolicy::writeConvergenceInfoToDiskFile(bool closeInsteadOfWriting, long int t)
{
	static ofstream fileWritingStream;
	static bool fileOpened = false;
	long int i;
	long int outputVal;
	char fileName[] = "GDout.txt";
	
	if (closeInsteadOfWriting)
	{
		// Close the file if it was previously opened.
		if (fileOpened) fileWritingStream.close();
	}
	else
	{
		if (!fileOpened)
		{
			// Open the file to overwrite it.
			fileWritingStream.open(fileName, ios::trunc | ios::binary);
			if (!fileWritingStream.is_open())
				{cout << "Unable to open file GDOut.txt for writing!  Exiting....\n\n"; exit(EXIT_FAILURE);}
			fileOpened = true;
			
			// Write the first line of the file.
			fileWritingStream << "k\t";
			for (i=0; i<kGridPoints; i++) fileWritingStream << k[i] << "\t";
			for (i=0; i<kGridPoints; i++) fileWritingStream << k[i] << "\t";
			fileWritingStream << "\n";
		}
		
		// Write the dictator's and social planner's policies for this time period.
		fileWritingStream << "d";
		if (t < 0) {fileWritingStream << "n"; t *= -1;}
		fileWritingStream << t << "\t";
		for (i=0; i<kGridPoints; i++) {outputVal = lround(newHDict[i] * 10000.0); fileWritingStream << outputVal << "\t";}
		for (i=0; i<kGridPoints; i++) {outputVal = lround(newHSocialPlanner[i] * 10000.0); fileWritingStream << outputVal << "\t";}
		fileWritingStream << "\n";
	}
}

// Compute the capital path for the dictator, starting from some initial capital.
void optimalPolicy::computeKPath(long double initK)
{
	long int simTime, timePeriodsForKPathMinusOne;
	
	kPath[0] = initK;
	timePeriodsForKPathMinusOne = timePeriodsForKPath - 1;
	for (simTime = 0; simTime < timePeriodsForKPathMinusOne; simTime++)
		kPath[simTime+1] = ( oneMinusDelta + (1.0 - interpolateFnOnKGrid(kPath[simTime], newHDict)) * wageOfK ) * kPath[simTime];
}

// Determine whether the capital path for the dictator is monotonically increasing (1), monotonically decreasing (-1), or otherwise (0).
int optimalPolicy::kDirection(void)
{
	int kUp = 1;  // Assume a monotonic increase until proven wrong.
	int kDown = -1;  // Assume a monotonic decrease until proven wrong.
	long int simTime, timePeriodsForKPathMinusOne;
	
	timePeriodsForKPathMinusOne = timePeriodsForKPath - 1;
	for (simTime = 0; simTime < timePeriodsForKPathMinusOne; simTime++)
	{
		if (kPath[simTime+1] >= kPath[simTime]) kDown = 0;
		if (kPath[simTime+1] <= kPath[simTime]) kUp = 0;
	}
	return (kDown + kUp);
}

int optimalPolicy::computeKPathAndGetDirection(long double initK)
{
	computeKPath(initK);
	return kDirection();
}

// The following procedure searches for a bifurcation point in the economy's capital, if it is managed by the dictator, between k = lowRange and k = highRange.
// This sets the value of bifurcationPoint.  Negative values have special meanings:
// -998 = no bifurcation point found due to apparent nonmonotonicities
// -991 = capital was always increasing over the range
// -990 = capital was always decreasing over the range
void optimalPolicy::findBifurcationPoint(long double lowRange, long double highRange)
{
	long double kMinBound, kMaxBound, kDiff, kTry, apparentNonmonotonicity, kDirection;
	apparentNonmonotonicity = 0.0L;
	kMinBound = lowRange;
	kMaxBound = highRange;
	
	kDiff = kMaxBound - kMinBound;
	kTry = kMinBound + kDiff / 2.0L;
	while ((kDiff > bifurcationSensitivity) & (apparentNonmonotonicity == 0.0L))
	{
		kDirection = computeKPathAndGetDirection(kTry);
		if (kDirection == 1) kMaxBound = kTry; // A monotonic increase in k occurred.
		else if (kDirection == -1) kMinBound = kTry; // A monotonic decrease in k occurred.
		else {apparentNonmonotonicity = -998.0L; break;} //* K was constant or nonmonotonic.  In theory we should never happen to land exactly on the bifurcation point, so assume this must be a nonmonotonicity.
		kDiff = kMaxBound - kMinBound;
		kTry = kMinBound + kDiff / 2.0L;
	}
	if (apparentNonmonotonicity != 0.0L) bifurcationPoint = apparentNonmonotonicity;
	else if (kTry <= (lowRange + bifurcationSensitivity)) bifurcationPoint = -991.0; // Capital was always increasing over the range.
	else if (kTry >= (highRange - bifurcationSensitivity)) bifurcationPoint = -990.0; // Capital was always decreasing over the range.
	else bifurcationPoint = kTry;
}

// After computing a bifurcation point, check that it really is a bifurcation point.
// This sets the values of bifPtCheckLow and bifPtCheckHigh.  Values of 0.0 will 
// result if everything checks out as anticipated.  Values of -998 mean
// that no bifurcation point had been found, due to apparent nonmonotonicities.
// Other values indicate an anomoly (capital was non-decreasing over time below 
// the bifurcation point, or vice-versa above), and give the percentage of the 
// bifurcation point at which the anomoly occurred (or if bifurcationPoint == 
// -991 the % of 0.999 times the low range, or if bifurcationPoint == -990, the 
// % of 1.001 times the high range).  A value of -991 will occur for the low
// check if capital was found to be always-increasing, and a value of -990 will 
// occur for the high check if capital was found to be always-decreasing.
void optimalPolicy::checkBifurcation(long double lowRange, long double highRange)
{
	long double anomDLow, anomDHigh, anomDiff, bifPtToUse;
	int kDirection;
	
	if (bifurcationPoint == -998.0)
	{
		// Because of apparent nonmonotonicities in growth, no bifurcation point was found.
		bifPtCheckLow = -998.0;
		bifPtCheckHigh = -998.0;
	}
	else
	{
		bifPtToUse = bifurcationPoint;
		if (bifurcationPoint == -991.0) bifPtToUse = lowRange * 0.999; // Capital should be always increasing over the range.
		if (bifurcationPoint == -990.0) bifPtToUse = highRange * 1.001; // Capital should be always decreasing over the range.
		
		// Check below the bifurcation point.
		if (bifurcationPoint == -991.0) bifPtCheckLow = -991.0; // We are apparently above any bifurcation point.
		else
		{
			anomDLow = lowRange;
			anomDiff = bifPtToUse - anomDLow;
			while (anomDiff > bifurcationSensitivity)
			{
				kDirection = computeKPathAndGetDirection(anomDLow);
				if (kDirection >= 0) break; // Capital grew instead of falling over time!!  This is an anomaly.
				else
				{
					// Capital fell over time as expected.
					anomDLow = anomDLow + anomDiff * bifurcationAnomalyStepFrac;
					anomDiff = bifPtToUse - anomDLow;
				}
			}
			if (anomDiff <= bifurcationSensitivity) bifPtCheckLow = 0.0; else bifPtCheckLow = 100.0 * anomDiff / bifPtToUse;
		}
		
		// Check above the bifurcation point.
		if (bifurcationPoint == -990.0) bifPtCheckHigh = -990.0; // We are apparently below any bifurcation point.
		else
		{
			anomDHigh = highRange;
			anomDiff = anomDHigh - bifPtToUse;
			while (anomDiff > bifurcationSensitivity)
			{
				kDirection = computeKPathAndGetDirection(anomDHigh);
				if (kDirection <= 0) break; // Capital fell instead of growing over time!!  This is an anomaly.
				else
				{
					// Capital grew over time as expected.
					anomDHigh = anomDHigh - anomDiff * bifurcationAnomalyStepFrac;
					anomDiff = anomDHigh - bifPtToUse;
				}
			}
			if (anomDiff <= bifurcationSensitivity) bifPtCheckHigh = 0.0; else bifPtCheckHigh = 100.0 * anomDiff / bifPtToUse;
		}
	}
}

// The following procedure checks whether and when the dictator's optimal policy function
// cuts through the social planner's optimal policy function (on the specified range).
// The result is placed into dictatorCutThroughSPPolicyPoint.
// -99 = invalid range
// -2 = dictator's optimal policy was initially less than or equal to the social planner's
// -1 = does not cut through
// other number = point at which cut-through occurs
void optimalPolicy::findDictatorCutThroughSPPolicyPoint(long double lowRange, long double highRange)
{
	long int i;
	bool stillNeedToCheckAtFirstPoint = true;

	if ((lowRange < kGridMin) || (highRange > kGridMax) || (lowRange >= highRange))
		dictatorCutThroughSPPolicyPoint = -99.0L;
	else for (i = 0; i < kGridPoints; i++)
		{
			if (k[i] > highRange) {dictatorCutThroughSPPolicyPoint = -1.0L; break;}
			if (k[i] >= lowRange)
			{
				if (stillNeedToCheckAtFirstPoint)
				{
					if (newHSocialPlanner[i] >= newHDict[i]) {dictatorCutThroughSPPolicyPoint = -2.0L; break;}
					else stillNeedToCheckAtFirstPoint = false;
				}
				else if (newHSocialPlanner[i] >= newHDict[i]) {dictatorCutThroughSPPolicyPoint = k[i]; break;}
			}
			if (i == kGridSteps) dictatorCutThroughSPPolicyPoint = -1.0L;
		}
}

// The following procedure checks whether the social planner economy ever declines, by finding the minimum fraction per annum growth on all K in the specified range.
// The result is placed into spMinFractionalGrowthPerAnnum.
// -99 = invalid range
// other number = min. fractional per annum growth
void optimalPolicy::findSPMinFractionalGrowthPerAnnum(long double lowRange, long double highRange)
{
	long int i;
	long double kFractionalGrowthPerAnnum;
	
	spMinFractionalGrowthPerAnnum = 9999999999.0L;
	if ((lowRange < kGridMin) || (highRange > kGridMax) || (lowRange >= highRange))
		spMinFractionalGrowthPerAnnum = -99.0L;
	else for (i = 0; i < kGridPoints; i++)
		{
			if (k[i] > highRange) break;
			if (k[i] >= lowRange)
			{
				kFractionalGrowthPerAnnum = wageOfK * (1.0L - newHSocialPlanner[i]) - delta;
				if (kFractionalGrowthPerAnnum < spMinFractionalGrowthPerAnnum) spMinFractionalGrowthPerAnnum = kFractionalGrowthPerAnnum;
			}
		}
}

// The following procedure finds the lowest values of newHSocialPlanner and newHDictator, and the highest value of (newHSocialPlanner - newHDictator).
// The results are placed into lowestHSP, lowestHDict, and highestHSPMinusHDict.
// -99 = invalid range
// other number = min. value
void optimalPolicy::findLowestHInfo(long double lowRange, long double highRange)
{
	long int i;
	long double diff;
	
	lowestHSP = 1.0L;
	lowestHDict = 1.0L;
	highestHSPMinusHDict = -1.0L;
	if ((lowRange < kGridMin) || (highRange > kGridMax) || (lowRange >= highRange))
	{
		lowestHSP = -99.0L;
		lowestHDict = -99.0L;
		highestHSPMinusHDict = -99.0L;
	}
	else for (i = 0; i < kGridPoints; i++)
		{
			if (k[i] > highRange) break;
			if (k[i] >= lowRange)
			{
				if (newHSocialPlanner[i] < lowestHSP) lowestHSP = newHSocialPlanner[i];
				if (newHDict[i] < lowestHDict) lowestHDict = newHDict[i];
				diff = newHSocialPlanner[i] - newHDict[i];
				if (diff > highestHSPMinusHDict) highestHSPMinusHDict = diff;
			}
		}
}

// The following procedure finds the value of K at which the biggest drop in the dictator's h*(K) occurs.
// Results are placed into kBiggestDropH1, kBiggestDropH10, kBiggestDropH50, kBiggestDropH250,
// biggestDropH1, biggestDropH10, biggestDropH50, biggestDropH250,
// qtAtStartBiggestDropH1, qtAtStartBiggestDropH10, qtAtStartBiggestDropH50, qtAtStartBiggestDropH250,
// qtAtEndBiggestDropH1, qtAtEndBiggestDropH10, qtAtEndBiggestDropH50, qtAtEndBiggestDropH250,
// qtplus1AtStartBiggestDropH1, qtplus1AtStartBiggestDropH10, qtplus1AtStartBiggestDropH50, qtplus1AtStartBiggestDropH250,
// qtplus1AtEndBiggestDropH1, qtplus1AtEndBiggestDropH10, qtplus1AtEndBiggestDropH50, qtplus1AtEndBiggestDropH250.
// ACTUALLY, AT THE MOMENT I AM ONLY COMPUTING QTPLUS1 AT THE START AND END OF THE DROP, BECAUSE QT IS HARDER TO COMPUTE: MUST KNOW K(t-1).
// The numbers refer to the number of steps on the grid over which the drop is assessed: 1, 10, 50, or 250.
// -99 = invalid range
// other number = value of K (if multiple values of K yield identical greatest negative slopes, this returns the first of the values),
//   of the drop in h (positive for a drop), and of the values of q(t) and q(t+1) at the start and end of the drop.
void optimalPolicy::findDictFastestHDropInfo(long double lowRange, long double highRange)
{
	long int i;
	
	kBiggestDropH1 = -99.0L;
	kBiggestDropH10 = -99.0L;
	kBiggestDropH50 = -99.0L;
	kBiggestDropH250 = -99.0L;
	biggestDropH1 = -99999.0L;
	biggestDropH10 = -99999.0L;
	biggestDropH50 = -99999.0L;
	biggestDropH250 = -99999.0L;
	/* qtAtStartBiggestDropH1 = -99.0L;
	qtAtStartBiggestDropH10 = -99.0L;
	qtAtStartBiggestDropH50 = -99.0L;
	qtAtStartBiggestDropH250 = -99.0L;
	qtAtEndBiggestDropH1 = -99.0L;
	qtAtEndBiggestDropH10 = -99.0L;
	qtAtEndBiggestDropH50 = -99.0L;
	qtAtEndBiggestDropH250 = -99.0L; */
	qtplus1AtStartBiggestDropH1 = -99.0L;
	qtplus1AtStartBiggestDropH10 = -99.0L;
	qtplus1AtStartBiggestDropH50 = -99.0L;
	qtplus1AtStartBiggestDropH250 = -99.0L;
	qtplus1AtEndBiggestDropH1 = -99.0L;
	qtplus1AtEndBiggestDropH10 = -99.0L;
	qtplus1AtEndBiggestDropH50 = -99.0L;
	qtplus1AtEndBiggestDropH250 = -99.0L;
	
	// Determine the values of qNext for each K given the optimal hDict.
	computeProduction();
	for (i = 0; i < kGridSteps; i++) h[i] = newHDict[i];
	computeConsumption();
	computeQNext();
	
	if ((lowRange < kGridMin) || (highRange > kGridMax) || (lowRange >= highRange))
	{
		biggestDropH1 = -99.0L;
		biggestDropH10 = -99.0L;
		biggestDropH50 = -99.0L;
		biggestDropH250 = -99.0L;
	}
	else for (i = 0; i < kGridSteps; i++)  // i < kGridSteps, not kGridPoints, is used to ensure that point i+1 exists for each point i tried.
		{
			if (k[i] > highRange) break;
			if (k[i] >= lowRange)
			{
				if ( (newHDict[i] - newHDict[i+1]) > biggestDropH1) {biggestDropH1 = (newHDict[i] - newHDict[i+1]); kBiggestDropH1 = expl((logl(k[i]) + logl(k[i+1]))/2); qtplus1AtStartBiggestDropH1 = qNext[i]; qtplus1AtEndBiggestDropH1 = qNext[i+1];}
				if ( (i-5)>=0 && (i+5)<kGridPoints ) if ( (newHDict[i-5] - newHDict[i+5]) > biggestDropH10) {biggestDropH10 = (newHDict[i-5] - newHDict[i+5]); kBiggestDropH10 = k[i]; qtplus1AtStartBiggestDropH10 = qNext[i-5]; qtplus1AtEndBiggestDropH10 = qNext[i+5];}
				if ( (i-25)>=0 && (i+25)<kGridPoints ) if ( (newHDict[i-25] - newHDict[i+25]) > biggestDropH50) {biggestDropH50 = (newHDict[i-25] - newHDict[i+25]); kBiggestDropH50 = k[i]; qtplus1AtStartBiggestDropH50 = qNext[i-25]; qtplus1AtEndBiggestDropH50 = qNext[i+25];}
				if ( (i-125)>=0 && (i+125)<kGridPoints ) if ( (newHDict[i-125] - newHDict[i+125]) > biggestDropH250) {biggestDropH250 = (newHDict[i-125] - newHDict[i+125]); kBiggestDropH250 = k[i]; qtplus1AtStartBiggestDropH250 = qNext[i-125]; qtplus1AtEndBiggestDropH250 = qNext[i+125];}
			}
		}
}

// The following procedure assesses the shape of the value function, V(K) (not V(h,K), but rather V(h*(K),K)).
// The results are placed into dictValFnMonotoniclyIncreasing, spValFnMonotoniclyIncreasing, dictValFnConcave, spValFnConcave.
// -9 = invalid range
// other number = 2 for strictly true, 1 for non-strictly true, 0 for false
void optimalPolicy::assessValueFunctionShape(long double lowRange, long double highRange)
{
	long int i;
	long double slopeDict, slopeSP, prevSlopeDict, prevSlopeSP;
	
	dictValFnMonotoniclyIncreasing = 2;
	spValFnMonotoniclyIncreasing = 2;
	dictValFnConcave = 2;
	spValFnConcave = 2;
	if ((lowRange < kGridMin) || (highRange > kGridMax) || (lowRange >= highRange))
	{
		dictValFnMonotoniclyIncreasing = -9;
		spValFnMonotoniclyIncreasing = -9;
		dictValFnConcave = -9;
		spValFnConcave = -9;
	}
	else for (i = 1; i < kGridPoints; i++)
		{
			slopeDict = (vwSum[dictator][i] - vwSum[dictator][i-1]) / (k[i] - k[i-1]);
			slopeSP = (vwSum[socialPlanner][i] - vwSum[socialPlanner][i-1]) / (k[i] - k[i-1]);
			if (k[i] > highRange) break;
			if (k[i] >= lowRange)
			{
				if (slopeDict < 0) dictValFnMonotoniclyIncreasing = 0; else if (slopeDict == 0) dictValFnMonotoniclyIncreasing = min(1, dictValFnMonotoniclyIncreasing);
				if (slopeSP < 0) spValFnMonotoniclyIncreasing = 0; else if (slopeSP == 0) spValFnMonotoniclyIncreasing = min(1, spValFnMonotoniclyIncreasing);
				if (i > 1)
				{
					if (slopeDict > prevSlopeDict) dictValFnConcave = 0; else if (slopeDict == prevSlopeDict) dictValFnConcave = min(1, dictValFnConcave);
					if (slopeSP > prevSlopeSP) spValFnConcave = 0; else if (slopeSP == prevSlopeSP) spValFnConcave = min(1, spValFnConcave);	
				}
			}
			prevSlopeDict = slopeDict;
			prevSlopeSP = slopeSP;
		}
}

// The following procedure assesses the shape of the policy function h*(K).
// The results are placed into dictPolFnShape, spPolFnShape.
// -9 = invalid range
// other number = 2 for decrease (non-strict) and rise (non-strict), 1 for decrease (non-strict) but not rise (non-strict), 0 if a rise occurs before a decrease or if there is a decrease then a rise then a decrease.
void optimalPolicy::assessPolicyFunctionShape(long double lowRange, long double highRange)
{
	long int i;
	long double deltaHDict, deltaHSP;
	
	dictPolFnShape = 3;
	spPolFnShape = 3;
	if ((lowRange < kGridMin) || (highRange > kGridMax) || (lowRange >= highRange))
	{
		dictPolFnShape = -9;
		spPolFnShape = -9;
	}
	else for (i = 1; i < kGridPoints; i++)
		{
			if (k[i] > highRange) break;
			if (k[i] >= lowRange)
			{
				deltaHDict = newHDict[i] - newHDict[i-1];
				if (dictPolFnShape == 3) {if (deltaHDict <= 0) dictPolFnShape = 1; else dictPolFnShape = 0;}  // Transition from 3 (initial) to 1 (decrease) or 0 (not a standard form).
				if (dictPolFnShape == 1) {if (deltaHDict > 0) dictPolFnShape = 2;}  // Transition from 1 (decrease) to 2 (decrease-rise) or stay at 1.
				if (dictPolFnShape == 2) {if (deltaHDict < 0) dictPolFnShape = 0;}  // Transition from 2 (decrease-rise) to 0 (not a standard form) or stay at 2.
				
				deltaHSP = newHSocialPlanner[i] - newHSocialPlanner[i-1];
				if (spPolFnShape == 3) {if (deltaHSP <= 0) spPolFnShape = 1; else spPolFnShape = 0;}  // Transition from 3 (initial) to 1 (decrease) or 0 (not a standard form).
				if (spPolFnShape == 1) {if (deltaHSP > 0) spPolFnShape = 2;}  // Transition from 1 (decrease) to 2 (decrease-rise) or stay at 1.
				if (spPolFnShape == 2) {if (deltaHSP < 0) spPolFnShape = 0;}  // Transition from 2 (decrease-rise) to 0 (not a standard form) or stay at 2.
			}
		}
}

// The following procedure computes the growth rates under dictator and social planner optimal policies.
// The results are placed into two sets of variables: (a) 5 growth rates corresponding to the low, med.-low, medium, med.-high, and high point in the range specified (evenly spaced in log-K space except rounded to the nearest K-grid point),
// and (b) an array of values for all K along with point numbers on the K-grid that tell the lowest and highest points within the range specified.
// The variables into which results are stored are:
// (a) fracDiscreteGrowthPerAnnumDictLow, fracDiscreteGrowthPerAnnumDictMedLow, fracDiscreteGrowthPerAnnumDictMed, fracDiscreteGrowthPerAnnumDictMedHigh, fracDiscreteGrowthPerAnnumDictHigh
//     fracDiscreteGrowthPerAnnumSPLow, fracDiscreteGrowthPerAnnumSPMedLow, fracDiscreteGrowthPerAnnumSPMed, fracDiscreteGrowthPerAnnumSPMedHigh, fracDiscreteGrowthPerAnnumSPHigh
// (b) fracDiscreteGrowthPerAnnumDict[i], fracDiscreteGrowthPerAnnumSP[i], minGridPointForGrowthRates, maxGridPointForGrowthRates
// For the low, med.-low, medium, med.-high, and high points:
// -999.0 = invalid range
// other number = fractional discrete growth per annum
// For the point numbers on the K-grid: -1 = invalid range, other value = point number.
void optimalPolicy::getGrowthRates(long double lowRange, long double highRange)
{
	long int i, medLowGridPointForGrowthRates, medGridPointForGrowthRates, medHighGridPointForGrowthRates;
	
	// Compute growth rates for all K on the grid.
	for (i = 0; i < kGridSteps; i++) fracDiscreteGrowthPerAnnumDict[i] = (1.0L - newHDict[i]) * wageOfK - delta;
	for (i = 0; i < kGridSteps; i++) fracDiscreteGrowthPerAnnumSP[i] = (1.0L - newHSocialPlanner[i]) * wageOfK - delta;
	
	// Determine the range of grid points corresponding to the specified range of K.
	if ((lowRange < kGridMin) || (highRange > kGridMax) || (lowRange >= highRange))
	{
		minGridPointForGrowthRates = -1;
		maxGridPointForGrowthRates = -1;
	}
	else
	{
		minGridPointForGrowthRates = kGridPoints; // Initialize to a non-existent grid point greater than the highest grid point number.
		maxGridPointForGrowthRates = -1; // Initialize to -1.
		for (i = 1; i < kGridPoints; i++)
		{
			if (k[i] > highRange) break;
			if (k[i] >= lowRange)
			{
				minGridPointForGrowthRates = min(minGridPointForGrowthRates, i);
				maxGridPointForGrowthRates = i;
			}
		}
	}
	
	// Determine the growth rates at the low, med.-low, medium, med.-high, and high points.
	if (maxGridPointForGrowthRates == -1)
	{
		fracDiscreteGrowthPerAnnumDictLow = -999.0L;
		fracDiscreteGrowthPerAnnumDictMedLow = -999.0L;
		fracDiscreteGrowthPerAnnumDictMed = -999.0L;
		fracDiscreteGrowthPerAnnumDictMedHigh = -999.0L;
		fracDiscreteGrowthPerAnnumDictHigh = -999.0L;
		fracDiscreteGrowthPerAnnumSPLow = -999.0L;
		fracDiscreteGrowthPerAnnumSPMedLow = -999.0L;
		fracDiscreteGrowthPerAnnumSPMed = -999.0L;
		fracDiscreteGrowthPerAnnumSPMedHigh = -999.0L;
		fracDiscreteGrowthPerAnnumSPHigh = -999.0L;
	}
	else
	{
		medLowGridPointForGrowthRates = round( 0.75 * (double) minGridPointForGrowthRates + 0.25 * (double) maxGridPointForGrowthRates );
		medGridPointForGrowthRates = round( ((double) minGridPointForGrowthRates + (double) maxGridPointForGrowthRates)/2.0 );
		medHighGridPointForGrowthRates = round( 0.25 * (double) minGridPointForGrowthRates + 0.75 * (double) maxGridPointForGrowthRates );
		fracDiscreteGrowthPerAnnumDictLow = fracDiscreteGrowthPerAnnumDict[minGridPointForGrowthRates];
		fracDiscreteGrowthPerAnnumDictMedLow = fracDiscreteGrowthPerAnnumDict[medLowGridPointForGrowthRates];
		fracDiscreteGrowthPerAnnumDictMed = fracDiscreteGrowthPerAnnumDict[medGridPointForGrowthRates];
		fracDiscreteGrowthPerAnnumDictMedHigh = fracDiscreteGrowthPerAnnumDict[medHighGridPointForGrowthRates];
		fracDiscreteGrowthPerAnnumDictHigh = fracDiscreteGrowthPerAnnumDict[maxGridPointForGrowthRates];
		fracDiscreteGrowthPerAnnumSPLow = fracDiscreteGrowthPerAnnumSP[minGridPointForGrowthRates];
		fracDiscreteGrowthPerAnnumSPMedLow = fracDiscreteGrowthPerAnnumSP[medLowGridPointForGrowthRates];
		fracDiscreteGrowthPerAnnumSPMed = fracDiscreteGrowthPerAnnumSP[medGridPointForGrowthRates];
		fracDiscreteGrowthPerAnnumSPMedHigh = fracDiscreteGrowthPerAnnumSP[medHighGridPointForGrowthRates];
		fracDiscreteGrowthPerAnnumSPHigh = fracDiscreteGrowthPerAnnumSP[maxGridPointForGrowthRates];
	}
}


void optimalPolicy::analyzeResults(void)
{
	int i;
	long double minOfRange, maxOfRange;
	char resultsIdentifierChar;
	
	for (i=0; i<=1; i++)
	{
		if (i==0) {minOfRange = 1.0L; maxOfRange = 100.0L; resultsIdentifierChar = 'A';}
		else {minOfRange = 0.1L; maxOfRange = 1000.0L; resultsIdentifierChar = 'B';}
		findBifurcationPoint(minOfRange, maxOfRange);
		checkBifurcation(minOfRange, maxOfRange);
		findDictatorCutThroughSPPolicyPoint(minOfRange, maxOfRange);
		findSPMinFractionalGrowthPerAnnum(minOfRange, maxOfRange);
		findLowestHInfo(minOfRange, maxOfRange);
		findDictFastestHDropInfo(minOfRange, maxOfRange);
		assessValueFunctionShape(minOfRange, maxOfRange);
		assessPolicyFunctionShape(minOfRange, maxOfRange);
		getGrowthRates(minOfRange, maxOfRange);
		writeResultsToDiskFile(resultsIdentifierChar);
		writeFullGrowthResultsToDiskFile(resultsIdentifierChar);
	}
}

void optimalPolicy::writeResultsToDiskFile(char resultsIdentifier)
{
	static ofstream fileStream;
	char fileName[] = "GDResX.txt";
	
	fileName[5] = resultsIdentifier;  // Replace the X in "GDResX.txt" with a character, such as A or B, that identifies which set of results this is.
	
	if (parameterSetNumber==0)
	{
		// Open the file with truncation.
		fileStream.open(fileName, ios::binary | ios::trunc);
		if (!fileStream.is_open()) {cout << "Unable to open file GDResX.txt for writing!  Exiting....\n\n"; exit(EXIT_FAILURE);}
		// Write the first line of the file.
		fileStream << "n\tcpform\talpha\ttheta\teta\tdelta\tsigma\tuMin\tbeta\trho\tnu\tbifPtDic\tbpChkLow\tbpChkHi\tdspCutPt\tspMinGro\tlowHSP\tlowHDict\thiHSPDD\t";
		fileStream << "dhdk1\tdhdk10\tdhdk50\tdhdk250\tdhdh1\tdhdh10\tdhdh50\tdhdh250\tdhdqa1\tdhdqa10\tdhdqa50\tdhdqa250\tdhdqb1\tdhdqb10\tdhdqb50\tdhdqb250\t";
		fileStream << "dvMonInc\tsvMonInc\tdvConcav\tsvConcav\tdHShape\tsHShape\t";
		fileStream << "dgr1\tdgr2\tdgr3\tdgr4\tdgr5\tsgr1\tsgr2\tsgr3\tsgr4\tsgr5\n";
	}
	else
	{
		// Open the file to append to it.
		fileStream.open(fileName, ios::out | ios::binary | ios::app);
		if (!fileStream.is_open()) {cout << "Unable to open file GDResX.txt for writing!  Exiting....\n\n"; exit(EXIT_FAILURE);}
	}
	
	// Write the results for the current parameter set.
	fileStream << parameterSetNumber << "\t";
	fileStream << cpform << "\t";
	fileStream << alpha << "\t";
	fileStream << theta << "\t";
	fileStream << eta << "\t";
	fileStream << delta << "\t";
	fileStream << sigma << "\t";
	fileStream << uMin << "\t";
	fileStream << beta << "\t";
	fileStream << rho << "\t";
	fileStream << nu << "\t";
	fileStream << bifurcationPoint << "\t";
	fileStream << bifPtCheckLow << "\t";
	fileStream << bifPtCheckHigh << "\t";
	fileStream << dictatorCutThroughSPPolicyPoint << "\t";
	fileStream << spMinFractionalGrowthPerAnnum << "\t";
	fileStream << lowestHSP << "\t";
	fileStream << lowestHDict << "\t";
	fileStream << highestHSPMinusHDict << "\t";
	fileStream << kBiggestDropH1 << "\t" << kBiggestDropH10 << "\t" << kBiggestDropH50 << "\t" << kBiggestDropH250 << "\t";
	fileStream << biggestDropH1 << "\t" << biggestDropH10 << "\t" << biggestDropH50 << "\t" << biggestDropH250 << "\t";
	fileStream << qtplus1AtStartBiggestDropH1 << "\t" << qtplus1AtStartBiggestDropH10 << "\t" << qtplus1AtStartBiggestDropH50 << "\t" << qtplus1AtStartBiggestDropH250 << "\t";
	fileStream << qtplus1AtEndBiggestDropH1 << "\t" << qtplus1AtEndBiggestDropH10 << "\t" << qtplus1AtEndBiggestDropH50 << "\t" << qtplus1AtEndBiggestDropH250 << "\t";
	fileStream << dictValFnMonotoniclyIncreasing << "\t" << spValFnMonotoniclyIncreasing << "\t" << dictValFnConcave << "\t" << spValFnConcave << "\t";
	fileStream << dictPolFnShape << "\t" << spPolFnShape << "\t";
	fileStream << fracDiscreteGrowthPerAnnumDictLow << "\t" << fracDiscreteGrowthPerAnnumDictMedLow << "\t" << fracDiscreteGrowthPerAnnumDictMed << "\t" << fracDiscreteGrowthPerAnnumDictMedHigh << "\t" << fracDiscreteGrowthPerAnnumDictHigh << "\t";
	fileStream << fracDiscreteGrowthPerAnnumSPLow << "\t" << fracDiscreteGrowthPerAnnumSPMedLow << "\t" << fracDiscreteGrowthPerAnnumSPMed << "\t" << fracDiscreteGrowthPerAnnumSPMedHigh << "\t" << fracDiscreteGrowthPerAnnumSPHigh;
	
	fileStream << "\n";
	
	fileStream.close();
}

void optimalPolicy::writeFullGrowthResultsToDiskFile(char resultsIdentifier)
{
	static ofstream fileStream;
	char fileName[] = "GDGrResX.txt";
	long int i;
	
	fileName[7] = resultsIdentifier;  // Replace the X in "GDGrResX.txt" with a character, such as A or B, that identifies which set of results this is.
	
	if (parameterSetNumber==0)
	{
		// Open the file with truncation.
		fileStream.open(fileName, ios::binary | ios::trunc);
		if (!fileStream.is_open()) {cout << "Unable to open file GDGrResX.txt for writing!  Exiting....\n\n"; exit(EXIT_FAILURE);}
		// Write the first line of the file.
		fileStream << "n\tk\tgrowDict\tgrowSP\n";
	}
	else
	{
		// Open the file to append to it.
		fileStream.open(fileName, ios::out | ios::binary | ios::app);
		if (!fileStream.is_open()) {cout << "Unable to open file GDGrResX.txt for writing!  Exiting....\n\n"; exit(EXIT_FAILURE);}
	}
	
	// Write the results for the current parameter set.
	if (maxGridPointForGrowthRates > -1)
		for (i = minGridPointForGrowthRates; i <= maxGridPointForGrowthRates; i++)
			fileStream << parameterSetNumber << "\t" << k[i] << "\t" << fracDiscreteGrowthPerAnnumDict[i] << "\t" << fracDiscreteGrowthPerAnnumSP[i] << "\n";
	
	fileStream.close();
}


// Create a disk file with the data for dictator and social planner optimal policies (as close as we've come to convergence).
void optimalPolicy::writeConvergedHFunctionsToDiskFile(bool successfulConvergence)
{
	ofstream fileWritingStream;
	long int i;
	//long int outputValDict, outputValSP;
	char fileName[101];
	
	if (successfulConvergence) sprintf( fileName, "Fn%d.raw", parameterSetNumber);
	else sprintf( fileName, "UnFn%d.raw", parameterSetNumber);
	
	// Open the file to overwrite it.
	fileWritingStream.open(fileName, ios::trunc | ios::binary);
	if (!fileWritingStream.is_open())
		{cout << "Unable to open file " << fileName << " for writing!  Exiting....\n\n"; exit(EXIT_FAILURE);}
	// Write the file.
	fileWritingStream << "k\thDict\thSP\n";
	fileWritingStream.precision(15);
	for (i=0; i<kGridPoints; i++)
	{
		//outputValDict = lround(newHDict[i] * 1000000.0); // 10000.0
		//outputValSP = lround(newHSocialPlanner[i] * 1000000.0);  // 10000.0
		fileWritingStream << k[i] << "\t" << newHDict[i] << "\t" << newHSocialPlanner[i] << "\n";
		//fileWritingStream << k[i] << "\t" << outputValDict << "\t" << outputValSP << "\n";
	}
	fileWritingStream.close();
}

// Create a disk file with the data for dictator and social planner optimal policies (as close as we've come to convergence).
void optimalPolicy::writeVWSumFunctionsToDiskFile(bool successfulConvergence)
{
	ofstream fileWritingStream;
	long int i;
	char fileName[101];
	
	if (successfulConvergence) sprintf( fileName, "Val%d.raw", parameterSetNumber);
	else sprintf( fileName, "UnVa%d.raw", parameterSetNumber);
	
	// Open the file to overwrite it.
	fileWritingStream.open(fileName, ios::trunc | ios::binary);
	if (!fileWritingStream.is_open())
		{cout << "Unable to open file " << fileName << " for writing!  Exiting....\n\n"; exit(EXIT_FAILURE);}
	// Write the file.
	fileWritingStream << "k\tvDict\tvSP\n";
	fileWritingStream.precision(15);
	for (i=0; i<kGridPoints; i++)
		fileWritingStream << k[i] << "\t" << vwSum[dictator][i] << "\t" << vwSum[socialPlanner][i] << "\n";
	fileWritingStream.close();
}

// Read a disk file with the data for dictator and social planner optimal policies.
// Return a boolean, 0 if successfully read (1 if the k-grid has changed), 2 if successfully read non-converged file but not converged file (3 if the k-grid has changed), 99 if didn't successfully read.
// If loading value functions, returns the min. of this value from the two reads.
short int optimalPolicy::readHFunctionsFromDiskFile(long int parameterSetNumberToLoad, bool errorIfCantRead, bool loadValFunctions)
{
	long int i;
	ifstream fileReadingStream;
	char fileName[101];
	char strBuffer[31];
	//char junkChar;
	char firstLineOfFile [] = "k\thDict\thSP";
	char firstLineOfFile2 [] = "k\tvDict\tvSP";
	short int fileAcceptable, fileAcceptable2;
	long int firstLineStringLength, bufferStringLength, dataPointsReadIn, dataPointsReadIn2;
	long double savedKDiffFirstGridStepForInterpolation, savedKDiffLastGridStepForInterpolation;
	long double accuracyAfterDiskFileMin, accuracyAfterDiskFileMax;
	long double readInKValue;
	bool keepReading;
	int c;
	
	// Determine the file name.
	sprintf( fileName, "Fn%d.raw", parameterSetNumberToLoad);
	
	// Open the file to read it.
	fileReadingStream.open(fileName, ios::in | ios::binary);
	if (!fileReadingStream.is_open())
	{
		//if (errorIfCantRead) {cout << "Unable to open file " << fileName << " for reading!  Exiting....\n\n"; exit(EXIT_FAILURE);}
		//else 
		fileAcceptable = 99;
		fileReadingStream.close();
		fileReadingStream.clear();
	}
	else
	{
		fileAcceptable = 0;
	}
	
	// If that didn't work, try to open the "UnFn#" file.
	if (fileAcceptable == 99)
	{
		sprintf( fileName, "UnFn%d.raw", parameterSetNumberToLoad);
		//fileReadingStream.open(fileName, ios::in | ios::binary);
		fileReadingStream.open(fileName, ios::in);
		if (!fileReadingStream.is_open())
		{
			if (errorIfCantRead) {cout << "Unable to open file " << fileName << " for reading!  Exiting....\n\n"; exit(EXIT_FAILURE);}
			else fileAcceptable = 99;
		}
		else
		{
			fileAcceptable = 2;
		}
	}
		
	// Read the first line of the file, and check it's okay.
	if (fileAcceptable < 99)
	{
		i = 0; keepReading = true;
		while (keepReading)
		{
			c = fileReadingStream.get();
			if ( (c == -1) || (c == '\n') || (c == '\r') )
			{
				bufferStringLength = i;
				keepReading = false;
			}
			else strBuffer[i] = c;
			if (i >= 29)
			{
				bufferStringLength = 30;
				keepReading = false;
			}
			i++;
		}
		strBuffer[bufferStringLength] = '\0';
		
		//fileReadingStream.get(strBuffer, 30, "\n\r"); // The "\n\r" was added to deal with the \r problem.
		firstLineStringLength = strlen(firstLineOfFile);
		//bufferStringLength = strlen(strBuffer);
		
		/*
		// This new version of checking okay was introduced because some old files have \r in place of \n.
		if (firstLineStringLength > bufferStringLength)
		{
			// First line has wrong length.
			if (errorIfCantRead) {cout << "File " << fileName << " does not seem to be a void h functions file, because the first line is the wrong length!  Exiting....\n\n"; exit(EXIT_FAILURE);}
			else fileAcceptable = 99;
		}
		else
		{
			for (i = 0; i < firstLineStringLength; i++) if (strBuffer[i] != firstLineOfFile[i])
				{
					// First line is not the correct start of the file.
					if (errorIfCantRead) {cout << "File " << fileName << " does not seem to be a void h functions file, because the first line is wrong!  Exiting....\n\n"; exit(EXIT_FAILURE);}
					else fileAcceptable = 99;
				}
			if (fileAcceptable < 99) fileReadingStream.get(junkChar); // Go past the "\n" in the stream.
			if ( (firstLineStringLength < bufferStringLength) && (strBuffer[firstLineStringLength] != '\r') )
			{
				// First line has wrong length.
				if (errorIfCantRead) {cout << "File " << fileName << " does not seem to be a void h functions file, because the first line is the wrong length!  Exiting....\n\n"; exit(EXIT_FAILURE);}
				else fileAcceptable = 99;
			}
		}
		*/
		
		if (firstLineStringLength != bufferStringLength)
		{
			// First line has wrong length.
			if (errorIfCantRead) {cout << "File " << fileName << " does not seem to be a valid h functions file, because the first line is the wrong length!  Exiting....\n\n"; exit(EXIT_FAILURE);}
			else fileAcceptable = 99;
		}
		else
		{
			for (i = 0; i < firstLineStringLength; i++) if (strBuffer[i] != firstLineOfFile[i])
				{
					// First line is not the correct start of the file.
					if (errorIfCantRead) {cout << "File " << fileName << " does not seem to be a valid h functions file, because the first line is wrong!  Exiting....\n\n"; exit(EXIT_FAILURE);}
					else fileAcceptable = 99;
				}
			//if (fileAcceptable < 99) fileReadingStream.get(junkChar); // Go past the "\n" in the stream.
		}
	}
	
	// Read the file.  The data will be placed in relevant variables on the k-grid.  This will work fine even if the number of kGridPoints has changed, although the policies will no longer necessarily be optimal.
	if (fileAcceptable < 99)
	{
		i = 0;
		while (fileReadingStream.eof()==0 && fileReadingStream.good()==1)
		{
			if (i >= maxKGridSize) {cout << "File " << fileName << " contains more than " << maxKGridSize << " points, exceeding the maxKGridSize!  Exiting....\n\n"; exit(EXIT_FAILURE);}
			fileReadingStream >> readInK[i] >> readInHDict[i] >> readInHSocialPlanner[i];
			//fileReadingStream >> readInK[i] >> junkChar >> readInHDict[i] >> junkChar >> readInHSocialPlanner[i] >> junkChar;
			i++;
		}
		dataPointsReadIn = i - 1; // Note the use of i-1 here to make the answer come out right.  The last line "read from" the file is rubbish; it is not really read because in trying to read in the program figures out it has reached the end of the file.  This works properly on my Mac, but might need to be re-checked on other machines.
		
		// In an earlier version the saved values of h used to be multiplied by 10,000 and rounded off to save on disk space, so I had to divide by 10,000 after reading, but I'm no longer doing that.
		//for (i = 0; i < dataPointsReadIn; i++) readInHDict[i] = readInHDict[i] / 10000.0
		//for (i = 0; i < dataPointsReadIn; i++) readInHSocialPlanner[i] = readInHSocialPlanner[i] / 10000.0
		
		if ( ( (readInK[0] - kGridMin) / kGridMin < 0.000001 ) && ( (readInK[dataPointsReadIn-1] - kGridMax) / kGridMax < 0.000001 ) && dataPointsReadIn == kGridPoints )
		{
			// Just copy the data into the grids.  The data were computed using the same grid as the program is currently using.
			for (i = 0; i < kGridPoints; i++)
			{
				newHDict[i] = readInHDict[i];
				newHSocialPlanner[i] = readInHSocialPlanner[i];
			}
		}
		else
		{
			// The data were computed using a different grid than is currently being used.  Interpolate onto the current grid.
			
			// Save and set the k-grid function and parameters for interpolation.
			for (i = 0; i < kGridPoints; i++) savedK[i] = kForInterpolation[i];
			savedKDiffFirstGridStepForInterpolation = kDiffFirstGridStepForInterpolation;
			savedKDiffLastGridStepForInterpolation = kDiffLastGridStepForInterpolation;
			for (i = 0; i < dataPointsReadIn; i++) kForInterpolation[i] = readInK[i];
			kGridMinForInterpolation = readInK[0];
			kGridMaxForInterpolation = readInK[dataPointsReadIn-1];
			kGridStepsForInterpolation = dataPointsReadIn - 1;
			kGridStepsMinusOneForInterpolation = kGridStepsForInterpolation - 1;
			kDiffFirstGridStepForInterpolation = readInK[1] - readInK[0];
			kDiffLastGridStepForInterpolation = readInK[kGridStepsForInterpolation] - readInK[kGridStepsMinusOneForInterpolation];
			
			// Load the read-in optimal policy functions into the appropriate arrays on the k-grid.
			// For points on the k-grid that are below or above those read in, use 1.0 for the optimal policy values.
			for (i = 0; i < kGridPoints; i++)
			{
				if ( (k[i] < kGridMinForInterpolation) || (k[i] > kGridMaxForInterpolation) )
				{
					newHDict[i] = 1.0;
					newHSocialPlanner[i] = 1.0;
				}
				else
				{
					newHDict[i] = interpolateFnOnKGrid(k[i], readInHDict);
					newHSocialPlanner[i] = interpolateFnOnKGrid(k[i], readInHSocialPlanner);
				}
			}
			
			// If the number of grid points has changed, or the min or max has changed, note so.
			accuracyAfterDiskFileMin = kGridMin * 0.000001;
			accuracyAfterDiskFileMax = kGridMax * 0.000001;
			if (  (kGridPoints != dataPointsReadIn) || ( (kGridMin < (readInK[0] - accuracyAfterDiskFileMin )) || (kGridMin > (readInK[0] + accuracyAfterDiskFileMin)) ) ||  ( (kGridMax < (readInK[dataPointsReadIn-1] - accuracyAfterDiskFileMax)) || (kGridMax > (readInK[dataPointsReadIn-1] + accuracyAfterDiskFileMax)) )  )
				fileAcceptable++;
			
			// Restore the k-grid function and parameters for interpolation.
			for (i = 0; i < kGridPoints; i++) kForInterpolation[i] = savedK[i];
			kGridMinForInterpolation = kGridMin;
			kGridMaxForInterpolation = kGridMax;
			kGridStepsForInterpolation = kGridSteps;
			kGridStepsMinusOneForInterpolation = kGridStepsMinusOne;
			kDiffFirstGridStepForInterpolation = savedKDiffFirstGridStepForInterpolation;
			kDiffLastGridStepForInterpolation = savedKDiffLastGridStepForInterpolation;
		}
			
		// Could set the parameterSetNumber appropriately, but that would not set the parameter values.  So, the parameterSetNumber
		// and all the parameters should already have been correct before calling this function, if the parameters are to be used 
		// in computations.
	}
	
	// Close the file.
	fileReadingStream.close();
	fileReadingStream.clear();
	
	
	// NOW READ VALUE FUNCTIONS IF DESIRED.
	if ((fileAcceptable < 99) && loadValFunctions)
	{
		
		// Determine the file name.
		sprintf( fileName, "Val%d.raw", parameterSetNumberToLoad);
		
		// Open the file to read it.
		fileReadingStream.open(fileName, ios::in | ios::binary);
		if (!fileReadingStream.is_open())
		{
			//if (errorIfCantRead) {cout << "Unable to open file " << fileName << " for reading!  Exiting....\n\n"; exit(EXIT_FAILURE);}
			//else
			fileAcceptable2 = 99;
			fileReadingStream.close();
			fileReadingStream.clear();
		}
		else
		{
			fileAcceptable2 = 0;
		}
		
		// If that didn't work, try to open the "UnVa#" file.
		if (fileAcceptable2 == 99)
		{
			sprintf( fileName, "UnVa%d.raw", parameterSetNumberToLoad);
			fileReadingStream.open(fileName, ios::in | ios::binary);
			if (!fileReadingStream.is_open())
			{
				if (errorIfCantRead) {cout << "Unable to open file " << fileName << " for reading!  Exiting....\n\n"; exit(EXIT_FAILURE);}
				else fileAcceptable2 = 99;
			}
			else
			{
				fileAcceptable2 = 2;
			}
		}
			
		// Read the first line of the file, and check it's okay.
		if (fileAcceptable2 < 99)
		{
			i = 0; keepReading = true;
			while (keepReading)
			{
				c = fileReadingStream.get();
				if ( (c == -1) || (c == '\n') || (c == '\r') )
				{
					bufferStringLength = i;
					keepReading = false;
				}
				else strBuffer[i] = c;
				if (i >= 29)
				{
					bufferStringLength = 30;
					keepReading = false;
				}
				i++;
			}
			strBuffer[bufferStringLength] = '\0';
			firstLineStringLength = strlen(firstLineOfFile2);
			if (firstLineStringLength != bufferStringLength)
			{
				// First line has wrong length.
				if (errorIfCantRead) {cout << "File " << fileName << " does not seem to be a valid value functions file, because the first line is the wrong length!  Exiting....\n\n"; exit(EXIT_FAILURE);}
				else fileAcceptable2 = 99;
			}
			else
			{
				for (i = 0; i < firstLineStringLength; i++) if (strBuffer[i] != firstLineOfFile2[i])
					{
						// First line is not the correct start of the file.
						if (errorIfCantRead) {cout << "File " << fileName << " does not seem to be a valid value functions file, because the first line is wrong!  Exiting....\n\n"; exit(EXIT_FAILURE);}
						else fileAcceptable2 = 99;
					}
				//if (fileAcceptable < 99) fileReadingStream.get(junkChar); // Go past the "\n" in the stream.
			}
		}
		
		// Read the file.  The data will be placed in relevant variables on the k-grid.  This will work fine even if the number of kGridPoints has changed, although the policies will no longer necessarily be optimal.
		if (fileAcceptable2 < 99)
		{
			i = 0; keepReading = true;
			while (fileReadingStream.eof()==0 && fileReadingStream.good()==1 && keepReading)
			{
				if (i >= maxKGridSize) {cout << "File " << fileName << " contains more than " << maxKGridSize << " points, exceeding the maxKGridSize!  Exiting....\n\n"; exit(EXIT_FAILURE);}
				fileReadingStream >> readInKValue >> readInVWSum[dictator][i] >> readInVWSum[socialPlanner][i];
				if (readInKValue != readInK[i]) {cout << "File " << fileName << " has values of k (at " << i << ") that don't match with the policy function file!  The policy function file said " << readInK[i] << " while the value function file said " << readInKValue << ".  Exiting....\n\n"; exit(EXIT_FAILURE);}
				if (i >= (dataPointsReadIn - 1)) keepReading = false;
				i++;
			}
			dataPointsReadIn2 = i; // Note the use of i, not i-1, here to make the answer come out right.  This loop exits one time earlier than the comparable loop for the policy function, because of the line that sets keepReading to false.  This works properly on my Mac, but might need to be re-checked on other machines.
			if (dataPointsReadIn2 != dataPointsReadIn) {cout << "File " << fileName << " has a different number of data lines than the policy function file!  Exiting....\n\n"; exit(EXIT_FAILURE);}
			
			if ( ( (readInK[0] - kGridMin) / kGridMin < 0.000001 ) && ( (readInK[dataPointsReadIn-1] - kGridMax) / kGridMax < 0.000001 ) && dataPointsReadIn2 == kGridPoints )
			{
				// Just copy the data into the grids.  The data were computed using the same grid as the program is currently using.
				for (i = 0; i < kGridPoints; i++)
				{
					vwSum[dictator][i] = readInVWSum[dictator][i];
					vwSum[socialPlanner][i] = readInVWSum[socialPlanner][i];
				}
			}
			else
			{
				// The data were computed using a different grid than is currently being used.  Interpolate onto the current grid.
			
				// Save and set the k-grid function and parameters for interpolation.
				for (i = 0; i < kGridPoints; i++) savedK[i] = kForInterpolation[i];
				savedKDiffFirstGridStepForInterpolation = kDiffFirstGridStepForInterpolation;
				savedKDiffLastGridStepForInterpolation = kDiffLastGridStepForInterpolation;
				for (i = 0; i < dataPointsReadIn; i++) kForInterpolation[i] = readInK[i];
				kGridMinForInterpolation = readInK[0];
				kGridMaxForInterpolation = readInK[dataPointsReadIn-1];
				kGridStepsForInterpolation = dataPointsReadIn - 1;
				kGridStepsMinusOneForInterpolation = kGridStepsForInterpolation - 1;
				kDiffFirstGridStepForInterpolation = readInK[1] - readInK[0];
				kDiffLastGridStepForInterpolation = readInK[kGridStepsForInterpolation] - readInK[kGridStepsMinusOneForInterpolation];
				
				// Load the read-in value functions into the appropriate arrays on the k-grid.
				// For points on the k-grid that are below or above those read in, load 99.9 for the values.
				for (i = 0; i < kGridPoints; i++)
				{
					if ( (k[i] < kGridMinForInterpolation) || (k[i] > kGridMaxForInterpolation) )
					{
						vwSum[dictator][i] = 99.9;
						vwSum[socialPlanner][i] = 99.9;
					}
					else
					{
						vwSum[dictator][i] = interpolateFnOnKGrid(k[i], readInVWSum[dictator]);
						vwSum[socialPlanner][i] = interpolateFnOnKGrid(k[i], readInVWSum[socialPlanner]);
					}
				}
				
				// If the number of grid points has changed, or the min or max has changed, note so.
				accuracyAfterDiskFileMin = kGridMin * 0.000001;
				accuracyAfterDiskFileMax = kGridMax * 0.000001;
				if (  (kGridPoints != dataPointsReadIn) || ( (kGridMin < (readInK[0] - accuracyAfterDiskFileMin )) || (kGridMin > (readInK[0] + accuracyAfterDiskFileMin)) ) ||  ( (kGridMax < (readInK[dataPointsReadIn-1] - accuracyAfterDiskFileMax)) || (kGridMax > (readInK[dataPointsReadIn-1] + accuracyAfterDiskFileMax)) )  )
					fileAcceptable2++;
			
				// Restore the k-grid function and parameters for interpolation.
				for (i = 0; i < kGridPoints; i++) kForInterpolation[i] = savedK[i];
				kGridMinForInterpolation = kGridMin;
				kGridMaxForInterpolation = kGridMax;
				kGridStepsForInterpolation = kGridSteps;
				kGridStepsMinusOneForInterpolation = kGridStepsMinusOne;
				kDiffFirstGridStepForInterpolation = savedKDiffFirstGridStepForInterpolation;
				kDiffLastGridStepForInterpolation = savedKDiffLastGridStepForInterpolation;
			}
				
			// Could set the parameterSetNumber appropriately, but that would not set the parameter values.  So, the parameterSetNumber
			// and all the parameters should already have been correct before calling this function, if the parameters are to be used 
			// in computations.
		}
		
		// Close the file.
		fileReadingStream.close();
			
		fileAcceptable = min(fileAcceptable, fileAcceptable2);	
	}
	
	return fileAcceptable ;
}

void varyParametersAndTestPolicies(void)
{
	optimalPolicy	*optimalPolicyStructPtr;
	long int maxT = 0; // 250; // 20
	long double initKMin = 1.5;
	long double initKMax = 250.0;
	long double kGridMin = 0.001;
	long double kGridMax = 100000.0;
	long double alpha, cpform, eta, theta, nu, delta, gamma, sigma, lambda, uMin, beta, rho;
	long int deltaChoice, thetaChoice, etaChoice, sigmaChoice, uMinChoice, betaChoice;
	long double thetaChoiceVal, etaChoiceVal, uMinChoiceVal;
	long double FofK0;
	long int startingParameterSetNumber, analyzePoliciesInputNum, writeConvergenceInfoInputNum, whatToDoIfNotConverging, readInPastAnalysesNum, useMoreExactVComputationInputNum, saveValueFunctionsInputNum, restrictedParameterSetsNum, useMeanPerturbationAtTime, useMinMaxPerturbationCode;
	bool onlyAnalyzePolicies, writeConvergenceInfo, readInPastAnalyses, useMoreExactVComputation, saveValueFunctions, restrictedParameterSets;
	bool successfulConvergence;
	short int successfullyRead;
	int kludgeToStepThroughCpform;  /* Later addition (3 Nov. 2000) to first do all cases with cpform==1, then cpform==2. */
	
	cout << "Enter the parameter set number at which to start (0 to 1727): ";
	cin >> startingParameterSetNumber;
	if ( (startingParameterSetNumber < 0) || (startingParameterSetNumber > 1727 ) )
		{cout << "Invalid parameter set number!  Exiting....\n\n"; exit(EXIT_FAILURE);}
		
	cout << "\nEnter 0 to find optimal policies, or 1 to analyze optimal policies based on past analyses: ";
	cin >> analyzePoliciesInputNum;
	if (analyzePoliciesInputNum == 1) onlyAnalyzePolicies = true; else onlyAnalyzePolicies = false;
	
	if (!onlyAnalyzePolicies)
	{
	
		cout << "\nEnter 0, or 1 to save detailed convergence info: ";
		cin >> writeConvergenceInfoInputNum;
		if (writeConvergenceInfoInputNum == 1) writeConvergenceInfo = true; else writeConvergenceInfo = false;
	
		if (maxT == 0)
		{
			cout << "\nIf convergence is not achieved after N+1 time steps, do what?\n";
			cout << "     0 Keep trying and hopefully converge sometime...\n";
			cout << "     N Give up after N+1 and switch to the next parameter set\n";
			cout << "Your choice: ";
			cin >> whatToDoIfNotConverging;
			if (whatToDoIfNotConverging < 0)
				{cout << "Invalid number!  Exiting....\n\n"; exit(EXIT_FAILURE);}
		}
		else whatToDoIfNotConverging = 0;
		
		if (maxT == 0)
		{
			cout << "\nEnter 0, or 1 to read in past policy Fn (or UnFn) files where available for each policy: ";
			cin >> readInPastAnalysesNum;
			if (readInPastAnalysesNum == 1) readInPastAnalyses = true;
			
			cout << "\nEnter 0, or 1 to only do p. sets 406, 694, 695, 1486, 1558, and 1559: "; // 407, 
			cin >> restrictedParameterSetsNum;
			if (restrictedParameterSetsNum == 1) restrictedParameterSets = true;
		}
		
		cout << "\nEnter 0, or 1 to use more exact V computation: ";
		cin >> useMoreExactVComputationInputNum;
		if (useMoreExactVComputationInputNum == 1) useMoreExactVComputation = true; else useMoreExactVComputation = false;
		//useMoreExactVComputationInputNum = 0;
		//useMoreExactVComputation = false;
		
		cout << "\nEnter 0, or another number for a time when the mean of the previous 50 periods of h (or as many as are available) will be substituted for the dictator's h values: ";
		cin >> useMeanPerturbationAtTime;
		
		if (useMeanPerturbationAtTime != 0)
		{
			cout << "\nEnter 0, or 1 to use the min over the 50 periods, or 2 to use the max, instead of the mean: ";
			cin >> useMinMaxPerturbationCode;
			if (useMinMaxPerturbationCode < 0 || useMinMaxPerturbationCode > 2)
				{cout << "Invalid number!  Exiting....\n\n"; exit(EXIT_FAILURE);}
		}
	}
	else
	{
		writeConvergenceInfo = false;
		whatToDoIfNotConverging = 0;
		readInPastAnalyses = false;
		useMoreExactVComputation = true;
	}
	
	cout << "\nEnter 0, or 1 to save value functions: ";
	cin >> saveValueFunctionsInputNum;
	if (saveValueFunctionsInputNum == 1) saveValueFunctions = true; else saveValueFunctions = false;

	for (kludgeToStepThroughCpform = 1; kludgeToStepThroughCpform <= 2; kludgeToStepThroughCpform++) /* Kludge of 3 Nov. 2000, to first do all cases with cpform==1, then cpform==2. */
	for (parameterSetNumber = startingParameterSetNumber; parameterSetNumber < 1728; parameterSetNumber++)
	{
		if (restrictedParameterSets && !(parameterSetNumber == 406 || parameterSetNumber == 694 || parameterSetNumber == 695 || parameterSetNumber == 1486 || parameterSetNumber == 1558 || parameterSetNumber == 1559))
			continue;   //  || parameterSetNumber == 407
		
		// Determine parameter values.
		lambda = 1.0; /* Change of 3 Nov. 2000.  Previously used 1.0. */
		gamma = 1.0;
		alpha = 1.0/3.0;
		rho = 0.02;
		deltaChoice = (parameterSetNumber % 2) + 1;  // 1,2 for the different values of delta
		if (deltaChoice == 1) delta = 0.05; else delta = 0.1;
		nu = powl( 3.0*(0.03+delta) / alpha , alpha ) * powl( rho / (1.0-alpha) , 1.0-alpha );
		cpform = (long double) ((parameterSetNumber % 8) / 2 + 1); // 1,2,3,4 for the different values of cpform
		if (cpform != kludgeToStepThroughCpform) continue; /* Kludge of 3 Nov. 2000, to first do all cases with cpform==1, then cpform==2. */
		thetaChoice = (parameterSetNumber % 24) / 8 + 1;  // 1,2,3 for the different rates of fall of catastrophe probability
		if (thetaChoice == 1) thetaChoiceVal = 0.75; else if (thetaChoice == 2) thetaChoiceVal = 0.9; else thetaChoiceVal = 0.97;
		if (cpform <= 2.0) theta = logl(thetaChoiceVal);
		else /* cpform == 3 or 4 */ theta = logl(thetaChoiceVal)/logl(1.1);  // Note that log(n)/log(1.1) = log base 1.1 of n.
		etaChoice = (parameterSetNumber % 96) / 24 + 1;  // 1,2,3,4 for the different catastrophe probabilities at the end of the first year in power
		if (etaChoice == 1) etaChoiceVal = 0.05; else if (etaChoice == 2) etaChoiceVal = 0.15; else if (etaChoice == 3) etaChoiceVal = 0.3; else etaChoiceVal = 0.5;
		FofK0 = 10.0 * alpha * powl( nu*(1.0-alpha)/rho , (1.0-alpha)/alpha );
		if (cpform == 1.0) eta = logl(etaChoiceVal) - theta * 10.0;
		else if (cpform == 3.0) eta = logl(etaChoiceVal) - theta * logl(10.0);
		else
		{
			if (cpform == 2.0) eta = logl(etaChoiceVal) - theta * 0.5 * FofK0;
			else /* cpform == 4.0 */ eta = logl(etaChoiceVal) - theta * logl(0.5 * FofK0);
		}
		sigmaChoice = (parameterSetNumber % 288) / 96 + 1;  // 1,2,3 for the different values of sigma
		if (sigmaChoice == 1) sigma = 0.5; else if (sigmaChoice == 2) sigma = 1.0; else sigma = 1.5;
		uMinChoice = (parameterSetNumber % 864) / 288 + 1;  // 1,2,3 for the different values of a1a2 used in determining uMin
		if (uMinChoice == 1) uMinChoiceVal = 0.001; else if (uMinChoice == 2) uMinChoiceVal = 0.01; else uMinChoiceVal = 0.05;
		if (sigma==1.0) uMin = logl(lambda * uMinChoiceVal * FofK0); else uMin = powl( lambda * uMinChoiceVal * FofK0 , 1.0 - sigma ) / (1.0 - sigma);
		betaChoice = parameterSetNumber / 864 + 1;  // 1,2 for the different values of beta
		if (betaChoice == 1) beta = 0.9; else beta = 0.95;

		// Determine optimal policy.
		optimalPolicyStructPtr = new optimalPolicy(useMeanPerturbationAtTime, useMinMaxPerturbationCode, useMoreExactVComputation, readInPastAnalyses, whatToDoIfNotConverging, writeConvergenceInfo, maxT, kGridMin, kGridMax, initKMin, initKMax, alpha, cpform, eta, theta, delta, gamma, sigma, lambda, uMin, beta, rho, nu);
		if (!onlyAnalyzePolicies)
		{
			// If readInPastAnalyses is true, then past analyses will be read in at the step where the final period (the first computed) is computed.  Otherwise we would have to do: if (readInPastAnalyses) successfullyRead = optimalPolicyStructPtr->readHFunctionsFromDiskFile(parameterSetNumber, true, true);
			successfulConvergence = optimalPolicyStructPtr->findOptimalPolicy();
		}
		else
		{
			successfullyRead = optimalPolicyStructPtr->readHFunctionsFromDiskFile(parameterSetNumber, true, true);
			successfulConvergence = (successfullyRead < 2);
		}
		
		// Analyze optimal policy.
		if (successfulConvergence)
		{
			optimalPolicyStructPtr->analyzeResults();
		}
		if (!onlyAnalyzePolicies) optimalPolicyStructPtr->writeConvergedHFunctionsToDiskFile(successfulConvergence);
		if (saveValueFunctions)
		{
			if (onlyAnalyzePolicies) {optimalPolicyStructPtr->computeMoreExactPresentAndFutureUtilityIfInPower(dictator); optimalPolicyStructPtr->computeMoreExactPresentAndFutureUtilityIfInPower(socialPlanner);}
			optimalPolicyStructPtr->writeVWSumFunctionsToDiskFile(successfulConvergence);
		}
		delete optimalPolicyStructPtr;
	}
}

void testWritingToDiskFile(int n)
{
	static ofstream fileStream;
	char fileName[] = "TestFile.txt";
	
	if (n==0)
	{
		// Open the file with truncation.
		fileStream.open(fileName, ios::out | ios::binary | ios::trunc);
		if (!fileStream.is_open()) {cout << "Unable to open file TestFile.txt for writing!  Exiting....\n\n"; exit(EXIT_FAILURE);}
		// Write the first line of the file.
		fileStream << "First\tline.\n";
	}
	else
	{
		// Open the file to append to it.
		fileStream.open(fileName, ios::out | ios::binary | ios::app);
		if (!fileStream.is_open()) {cout << "Unable to open file TestFile.txt for writing!  Exiting....\n\n"; exit(EXIT_FAILURE);}
	}
	
	// Write the results for the current parameter set.
	fileStream << "hi" << "\t" << "ho.";
	fileStream << "\n";
	
	fileStream.close();
	
	cout << "Finished writing to test file.\n";
}


/* Main. */
int main(void)
{	

	set_new_handler( newFailed );
	gSpareBlockPtr = new char[20480];

	srand(time(NULL));
	
	varyParametersAndTestPolicies();
	
	cout << "Done with Growth under Dictatorship\n\nPress command-Q to quit.";
	
	
	delete [] gSpareBlockPtr;

	return 0;
}

void maybeTerminateByUser(void)
{
	//int charsInBuffer;
	//int nextChar;
	
	//nextChar = cin.peek();
	//cout << "\nnextChar = " << nextChar << "\n";
	//if ( nextChar == 113 ) exit(EXIT_SUCCESS);
	//charsInBuffer = cin.in_avail();
	#ifdef TARGETMACHINEISMACPPC
	if ( kbhit() ) exit(EXIT_FAILURE);
	#endif
}

void newFailed(void)
{
	if ( gSpareBlockPtr != 0 )
	{
		delete [] gSpareBlockPtr;
		gSpareBlockPtr = 0;
	}
	cout << "Out of memory!!!  Exiting....\n\n";
	exit(EXIT_FAILURE);
}