(* Mathematica Package *)

BeginPackage["LamDelosmeSA`", { "MultivariateStatistics`"}]
(* Exported symbols added here with SymbolName::usage *)  

LamDelosmeSA::usage = "LamDelosmeSA[metric, variables] performs Lam-Delosme simulated annealing on a function given by metric. Format for variables is {name, low_end_of_estimated_range, high_end_of_estimated_range}."
ldsa::init = 
  "Failed to find valid intitialization in `1` tries.";
  
Begin["`Private`"] (* Begin Private Context *) 

SetAttributes[LamDelosmeSA, HoldFirst];

DefaultNumForgetMoves      = 10000;
DefaultNumInitMoves        = 10000;
DefaultInitInverseTemp     = 0;
DefaultLambda              = .01;
DefaultMeanMemoryProduce   = 600;
DefaultStdDevMemoryProduce = 30000;
DefaultEstimateUpdateFreq  = 100;
DefaultTemperUpdateFreq    = 10;    
DefaultPercentOfRange      = .03;
DefaultReportingFunc       = Report;
DefaultMoveGenerator	   = ReinitzUpdate;
DefaultStoppingCriterion   = "DeltaMeanEnergy";
DefaultStoppingDeltaE      = 10^-3;
DefaultMaxFrozenIters      = 3;
DefaultSweepsThetaUpdate   = 100;
DefaultIntializeAttempts   = 10^6;

Options[LamDelosmeSA] = {
		InitialInverseTemp -> DefaultInitInverseTemp,
		NumForgetMoves -> DefaultNumForgetMoves,
		NumInitMoves -> DefaultNumInitMoves,
		Lambda -> DefaultLambda, 
		MeanMemProd -> DefaultMeanMemoryProduce, 
		StdDevMemProd -> DefaultStdDevMemoryProduce, 
		EstimateUpdateFreq -> DefaultEstimateUpdateFreq, 
		TemperatureUpdateFreq -> DefaultTemperUpdateFreq,
		MoveGenerator -> DefaultMoveGenerator,
		InitState -> {},
		MaximumFrozenIterations -> DefaultMaxFrozenIters,
		InitialThetas -> {},
		ReportFunc -> DoNothing,
		TalkativeInit -> False,
		TalkativeErrors -> False,
		UseSlaves -> False,
		MixingFrequency -> 10,
		StoppingCriterion -> DefaultStoppingCriterion,
		StoppingDeltaEnergy -> DefaultStoppingDeltaE,
		Debugging -> False,
		SweepsBetweenThetaUpdates -> DefaultSweepsThetaUpdate, 
		InitializeTries -> DefaultIntializeAttempts
	};

LamDelosmeSA[metric_, variables_List,  opts:OptionsPattern[]] := 
	LamDelosmeSA[metric, variables, {True}, opts]; 
	
LamDelosmeSA[metric_, variables_List, hardConstraints_:{True}, OptionsPattern[]] := 		
	Module[{variableNames, estVariableRanges, initState, userInit, needsInit, initRest, 
		state, thetas, energy, sum, sqrdSum, numInitMoves, totalNumMoves, estimateMean, estimateVrnc,
		estimateStDv, acceptRatio, lambda, estimateUpdateFreq, meanEstimator, stdvEstimator,
		s, deltaS, tempUpdateFreq, alpha, numTemp, numAccepted, keepRunning, maxFrozen, timesFrozen,
		bestEnergyEver, bestStateEver, initThetas,report, talkativeInit,
		talkativeErrs, constraint, moveGenerator, mixingIters, mean, deltaMean, oldMean,
		numParallelProcessors, estimateUpdateFreqEach, mmp, sdmp,variance,
		maxParallelProcessors,parSuccess,slaves, userSlaves, mixingFreq,stopCrit, stopDeltaEnergy, debug, allstats,
		acceptStats, emptyStats, thetaMin, thetaMax, updateTheta, gain, numForgetMoves, thetaUpdateSweeps, movesPerThetaUpdt,
		completedMoves, movesPerProcessor, movesPerUpdate, movesLeft, movesThisIteration, tsum, tsqrdSum, tnumAccepted},
		
		If[OptionValue[UseSlaves],
			numParallelProcessors = Length[Parallel`Parallel`$Slaves];
			Parallel`Parallel`ExportEnvironment[UpdateEstimator, Estimate, ReinitzUpdate, InitialRandomValue, 
				              					InitializeVariables, InitializeThetas, parameterA, parameterB, Report, 
				              					NewEstimator, initializationPeriod, sweep, DefaultNumInitMoves, DefaultLambda, 
				              					DefaultMeanMemoryProduce, DefaultStdDevMemoryProduce, DefaultEstimateUpdateFreq, 
				              					DefaultTemperUpdateFreq, DefaultPercentOfRange, DefaultReportingFunc, DefaultMaxFrozenIters,
				              					DefaultMoveGenerator, DefaultStoppingCriterion, DefaultStoppingDeltaE, 
				              					forgettingPeriod, AddTrial, AddSuccess, DefaultInitInverseTemp,
				              					aggregateParallelStats, AccumulateStats, DefaultSweepsThetaUpdate, DefaultIntializeAttempts, 
				              					ValidlyInitializeVariables];,
			(* else *)
			numParallelProcessors = 0
		];

		(* True if we're in debugging mode. False otherwise. *)		
		debug = If[OptionValue[Debugging], Print["debugging mode on."]; True, False, False];

		(* List of the names of the variables we're optimizing. *)		
		variableNames     = Map[#[[1]]&, variables];
		(* Approximate range over which to optimize each variable. *) 
		estVariableRanges = Map[#[[1]] -> {#[[2]],#[[3]]}&, variables];

		(* Single conjunction that encodes all of the hard constraints *)		
		constraint = Fold[(#1 && #2)&, True, hardConstraints];

		(* At the moment, the idea of a user-supplied initialization exists, but is poorly thought-out. 
		   Right here, what we have is that if we're using more than one processor we initialize all of the 
		   processors' states independently, but if we are using only one processor we allow the user to set 
		   some subset of the initial state. *)		
		If[numParallelProcessors > 0,
			state = Table[ValidlyInitializeVariables[variableNames, estVariableRanges, constraint, OptionValue[InitializeTries]], {numParallelProcessors}];,
			(* else *)
			initState = OptionValue[InitState];
			userInit = Cases[initState, (x_ -> y_) /; MemberQ[variableNames, x] ];
			needsInit = Cases[variableNames, x_ /; ! MemberQ[Map[#[[1]] &, userInit], x]];
			initRest = ValidlyInitializeVariables[needsInit, estVariableRanges, constraint, OptionValue[InitializeTries]];			
			state = Union[userInit, initRest];
		];
		acceptStats = Map[(#1 -> Stats[0, 0])&, variableNames];
		emptyStats = acceptStats; (* so we don't have to keep redoing the initialization *)

		(* Initialization for the theta values *)				
		initThetas = OptionValue[InitialThetas];
		userInit = Cases[initThetas, (x_ -> y_) /; MemberQ[variableNames, x]];
		needsInit = Cases[variableNames, x_ /; ! MemberQ[Map[#[[1]] &, userInit], x]];
		initRest = InitializeThetas[needsInit, estVariableRanges];					
		thetas = Union[userInit, initRest];		
		initThetas = thetas;
		
		thetaMin = 10^-5; (* Maybe this should be made a parameter? *)
		thetaMax = Map[#[[1]]->#[[2]]*(1/DefaultPercentOfRange)*10&,thetas];
		gain = 3; (* Maybe this should be made a parameter? *)
		updateTheta[name_ -> theta_, stats_] := 
			name->Min[updateTheta[theta, name/.stats], name /.thetaMax];
		updateTheta[theta_?NumberQ, Stats[accepted_, hits_]] := 
			Max[E^(Log[theta] + gain*((accepted/hits)-.44)), thetaMin];		
		
		stopCrit          = OptionValue[StoppingCriterion];
		stopDeltaEnergy   = OptionValue[StoppingDeltaEnergy];
		talkativeInit     = OptionValue[TalkativeInit];
		talkativeErrs     = OptionValue[TalkativeErrors];		
		moveGenerator     = OptionValue[MoveGenerator];		
		report            = OptionValue[ReportFunc];		
		numInitMoves      = OptionValue[NumInitMoves];
(*		numInitMoves      = (numInitMoves*Length[variableNames]);*)
		numForgetMoves    = OptionValue[NumForgetMoves];
		mixingFreq    	  = OptionValue[MixingFrequency];
		thetaUpdateSweeps = OptionValue[SweepsBetweenThetaUpdates];
		movesPerThetaUpdt = Ceiling[(thetaUpdateSweeps*Length[variableNames]) / (* it's important for this to be a multiple of the
																				  number of processors. *)
								Max[1,numParallelProcessors]]*Max[1,numParallelProcessors]; 
		
		If[debug,
			Print["reportFunction:<",report, "> talkativeInit:<", talkativeInit, ">"];
			Print["testing your report function, with empty state, energy 1, 0 iterations."];
			report[{}, 1, 0];
		];
		
		If[!MatchQ[stopCrit, "DeltaMeanEnergy"] && !MatchQ[stopCrit, "AcceptRatio"],
			Print["Stopping criterion given is <", stopCrit, ">, which is neither DeltaMeanEnergy nor AcceptRatio.\n",
				  "Defaulting to <" , DefaultStoppingCriterion, ">"];
			stopCrit = DefaultStoppingCriterion;
		];

		s = OptionValue[InitialInverseTemp];
		(* Perform the "forgetting" decorrelation steps *)
		completedMoves = 0;
		While[completedMoves < numForgetMoves,
			movesLeft = numForgetMoves - completedMoves;
			If[numParallelProcessors > 0,
				movesPerProcessor = Ceiling[Min[movesLeft,movesPerThetaUpdt]/numParallelProcessors];		
				movesPerUpdate = movesPerProcessor * numParallelProcessors;				
				Parallel`Parallel`ExportEnvironment[moveGenerator, Evaluate[OptionValue[MoveGenerator]]];
				Parallel`Parallel`ExportEnvironment[report, Evaluate[OptionValue[ReportFunc]]];			
				Module[{results,movesHolder,namesHolder,movegenHolder,thetasHolder,metricHolder,constraintHolder,
					    talkErrHolder,talkInitHolder,invTempHolder,statsHolder},
					(* OMFG, HAXX! 
					   I can't figure out a better way to get it to evaluate these variables & pass their *values* to the remote kernels.
					   If you're thinking, 'just ExportEnv them,' then you clearly haven't looked at previous svn check-ins.
					*)
					results = ReleaseHold[ReplaceAll[Hold[Parallel`Evaluate`ParallelMap[forgettingPeriod[#,movesHolder,namesHolder,movegenHolder,
																										 thetasHolder,metricHolder,constraintHolder,
																										 talkErrHolder,talkInitHolder,invTempHolder,statsHolder]&, 
						                                    							state]],
						                            {movesHolder->movesPerProcessor, namesHolder->variableNames,movegenHolder->moveGenerator, 
						                             thetasHolder->thetas, metricHolder->metric, constraintHolder->constraint, 
						                             talkErrHolder->talkativeErrs, talkInitHolder->talkativeInit,invTempHolder->s,
						                             statsHolder->emptyStats}]];
					state = Map[#[[1]]&, results];
					energy = Map[#[[2]]&, results];			
					allstats = Map[#[[3]]&, results];
					acceptStats = Map[aggregateParallelStats[#, allstats]&, variableNames];
				],
				(* else *)
				movesPerUpdate = Min[movesLeft,movesPerThetaUpdt];
				{state, energy, acceptStats} = forgettingPeriod[state, movesPerUpdate, variableNames, moveGenerator, 
						                                        thetas, metric, constraint, talkativeErrs, talkativeInit, s, emptyStats];
			];
			thetas = Map[updateTheta[#, acceptStats]&, thetas];
			completedMoves += movesPerUpdate;
		];	
		
		totalNumMoves = completedMoves;

		numAccepted = sum = sqrdSum = 0;
		acceptStats = emptyStats;
		(* Perform the intialization steps *)		
		completedMoves = 0;
		While[completedMoves < numInitMoves,
			movesLeft = numInitMoves - completedMoves;
			If[numParallelProcessors > 0,
				movesPerProcessor = Ceiling[Min[movesPerThetaUpdt, movesLeft]/numParallelProcessors];		
				movesPerUpdate = movesPerProcessor * numParallelProcessors;	
				Module[{results,a,b,c,d,e,f,g,h,invTempH, statsH},
					(* OMFG, HAXX! 
					   I can't figure out a better way to get it to evaluate these variables & pass their *values* to the remote kernels.
					   If you're thinking, 'just ExportEnv them,' then you clearly haven't looked at previous svn check-ins.
					*)
					results = ReleaseHold[ReplaceAll[Hold[Parallel`Evaluate`ParallelMap[initializationPeriod[#,a,b,c,d,e,f,g,h,
																											 invTempH, statsH]&, 
						                                    							state]],
						                            {a->movesPerProcessor, b->variableNames, c->moveGenerator, 
						                             d->thetas, e->metric, f->constraint, g->talkativeErrs, h->talkativeInit,
						                             invTempH->s, statsH->emptyStats} ]];
					state = Map[#[[1]]&, results];
					sum += Total[Map[#[[2]]&, results]];
					sqrdSum += Total[Map[#[[3]]&, results]];
					energy = Map[#[[4]]&, results];
					allstats = Map[#[[5]]&, results];
					acceptStats = Map[aggregateParallelStats[#, allstats]&, variableNames];
					numAccepted += Total[Map[#[[6]]&, results]];					
				],
				(* else *)
				movesPerUpdate = Min[movesPerThetaUpdt, movesLeft];
				{state, tsum, tsqrdSum, energy, acceptStats, tnumAccepted} = initializationPeriod[state, movesPerUpdate, variableNames, moveGenerator, 
						                                                			              thetas, metric, constraint, talkativeErrs, talkativeInit, s, emptyStats];			
				sum += tsum;
				sqrdSum += tsqrdSum;              
				numAccepted += tnumAccepted;
			];
			thetas = Map[updateTheta[#, acceptStats]&, thetas];			
			completedMoves += movesPerUpdate;
		];		
		totalNumMoves += completedMoves;
		
		estimateMean = sum/completedMoves;
		estimateVrnc = sqrdSum/completedMoves - estimateMean^2;
		(* floating point weirdness sometimes results in variances that should be 0 coming out negative *)
		If[estimateVrnc < 0, estimateVrnc = 0];
		estimateStDv = Sqrt[estimateVrnc];
		
		acceptRatio = numAccepted/completedMoves;
		
		If[debug,
			Print["sum         : ", sum];
			Print["sqrdSum     : ", sqrdSum];
			Print["completedMvs: ", completedMoves];
			Print["estimateMean: ", estimateMean];
			Print["estimateVrnc: ", estimateVrnc];
			Print["estimateStdv: ", estimateStDv];
			Print["acceptRatio : ", acceptRatio];
		];

		lambda = OptionValue[Lambda];
		estimateUpdateFreq = OptionValue[EstimateUpdateFreq];
(*		estimateUpdateFreq = (estimateUpdateFreq*Length[variableNames]);*)
		If[numParallelProcessors > 0,
			estimateUpdateFreqEach = Ceiling[estimateUpdateFreq/numParallelProcessors];
			estimateUpdateFreq = numParallelProcessors*estimateUpdateFreqEach;,
			(*else*)
			estimateUpdateFreqEach = estimateUpdateFreq;
		];						
		
		mmp = OptionValue[MeanMemProd];
		meanEstimator = NewEstimator[mmp, 
									 estimateUpdateFreq,
									 lambda,
									 estimateMean, estimateVrnc/(estimateMean^2)];

		sdmp = OptionValue[StdDevMemProd];
		stdvEstimator = NewEstimator[sdmp, 
									 estimateUpdateFreq,
									 lambda,
									 estimateStDv, estimateStDv/estimateMean];		
			
		alpha = 4 * acceptRatio * ((1-acceptRatio)/(2-acceptRatio))^2;						 
		
		deltaS = .5/estimateStDv;
		
		If[debug,
			Print["deltaS: ", deltaS];
		];
		
		tempUpdateFreq = OptionValue[TemperatureUpdateFreq];

		numTemp = 0;					
		keepRunning = True;
		maxFrozen = OptionValue[MaximumFrozenIterations];
		timesFrozen = 0;
		bestEnergyEver = Infinity;
		bestStateEver = {};
		mixingIters = 0;
		mean = estimateMean;
		variance = estimateVrnc;
		completedMoves = 0;
		deltaMean = 0;

		numAccepted = sum = sqrdSum = 0;					

		While[keepRunning,
			movesThisIteration = Min[estimateUpdateFreq-Mod[completedMoves,estimateUpdateFreq], 
									 movesPerThetaUpdt -Mod[completedMoves,movesPerThetaUpdt ]];
			Module[{ratio, a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,statsH,numProcH,printablestats},
				If[numParallelProcessors > 0,		
					If[debug,					
						Print["s<", s, "> deltaS<", deltaS, ">"];
						Print["eststddv<", Estimate[stdvEstimator, s], "> estmean<", Estimate[meanEstimator, s], ">"];
					];
					Module[{results},
						results = ReleaseHold[ReplaceAll[Hold[Parallel`Evaluate`ParallelMap[sweep[#[[1]], #[[2]], o, 
																								  a,b,c,d,e,f,g,h,i,j,k,l,m,n,p,statsH,numProcH
																								 ]&, 
																						    Thread[{state, energy}]]],
											             {a->meanEstimator, b->stdvEstimator, c->variableNames, d->moveGenerator, 
											              e->thetas, f->constraint, g->metric, h->talkativeErrs, i->deltaS, j->numTemp,
											              k->tempUpdateFreq, l->lambda, m->alpha, n->movesThisIteration/numParallelProcessors, o->s,
											              p->estimateMean,statsH->emptyStats,numProcH->numParallelProcessors}
														]];
						state = Map[#[[1]]&, results];
						sum += Total[Map[#[[2]]&, results]];
						sqrdSum += Total[Map[#[[3]]&, results]];
						energy = Map[#[[4]]&, results];
						numAccepted += Total[Map[#[[5]]&, results]];
						(* computation of these is deterministic, so we'll just take the one from the first processor. *)
						s = results[[1,6]];
						deltaS = results[[1,7]];
						allstats = Append[Map[#[[8]]&, results],acceptStats];
						acceptStats = Map[aggregateParallelStats[#, allstats]&, variableNames];						
						numTemp = tempUpdateFreq - Mod[(movesThisIteration/numParallelProcessors)-numTemp, tempUpdateFreq];						
					],

					(* else *)
					{state, tsum, tsqrdSum, energy, tnumAccepted, s, deltaS, acceptStats} = 
						sweep[state, energy, s, meanEstimator, stdvEstimator, variableNames, 
						   	  moveGenerator, thetas, constraint, metric, talkativeErrs, 
						      deltaS, numTemp, tempUpdateFreq, lambda, alpha, movesThisIteration, estimateMean, acceptStats, 1];							
				    sum += tsum;
				    sqrdSum += tsqrdSum;
				    numAccepted += tnumAccepted;
					numTemp = tempUpdateFreq - Mod[movesThisIteration-numTemp, tempUpdateFreq];												      
				];	
				completedMoves += movesThisIteration;
				totalNumMoves += movesThisIteration;
							
				(* do things that happen some multiple of every TAU updates *)
				If[Mod[completedMoves, estimateUpdateFreq]==0,
					If[numParallelProcessors > 0,
						mixingIters++;
						If[mixingIters >= mixingFreq,
							Module[{boltzmannFactors, total, probs,  sepairs, matters, lowestEnergy, lowestEposit,
								th, keep, energies, states},
								matters = Map[-#*s > Log[$MinNumber]&, energy];
								If[matters == Table[False, {numParallelProcessors}],
									lowestEnergy = Min[energy];
									lowestEposit = Position[energy, lowestEnergy][[1,1]];
									state = Table[Part[state,lowestEposit], {numParallelProcessors}];
									energy = Table[lowestEnergy, {numParallelProcessors}];, 
									(*else*)
									th = Thread[{matters, state, energy}];
									keep = Cases[th, {mm_,ss_,ee_}  /; mm];
									energies = Map[#[[3]]&, keep];
									states = Map[#[[2]]&, keep];
									boltzmannFactors = Map[E^(Max[-#*s,Log[$MinNumber]+1])&, energies];
									total = Total[boltzmannFactors];
									probs = Map[#/total&, boltzmannFactors];
									sepairs = Table[RandomChoice[probs->Thread[{states, energies}]], {numParallelProcessors}];
									state = Map[#[[1]]&, sepairs];
									energy = Map[#[[2]]&, sepairs];								
								];							
							];
							mixingIters = 0;
						];
					];
					oldMean = mean;				
					mean = sum / estimateUpdateFreq;
					variance = sqrdSum / estimateUpdateFreq;
					acceptRatio = numAccepted / estimateUpdateFreq;
				
					numAccepted = sqrdSum = sum = 0;
									
					deltaMean = Abs[mean - oldMean];
				
					If[MatchQ[stopCrit, "AcceptRatio"] && acceptRatio == 0,
						timesFrozen++,
						(* else *)
						If[MatchQ[stopCrit, "DeltaMeanEnergy"] && deltaMean < stopDeltaEnergy,
							timesFrozen++,
							timesFrozen = 0
						];,
						(* fail! *)
						Print["This shouldn't happen: neither stopping criterion?"];
					];
				
					If[timesFrozen == maxFrozen,
						keepRunning = False
					];
					
					If[debug,
						Print["mean: ", mean];
						Print["vrnc: ", variance];
						Print["s   : ", s];
					];
				
					meanEstimator = UpdateEstimator[meanEstimator, s, mean];

					If[variance > 0, 
						stdvEstimator = UpdateEstimator[stdvEstimator, s, Sqrt[variance]];
					];
					estimateMean = Estimate[meanEstimator, s];
					If[MatchQ[estimateMean,Indeterminate],
						Print[meanEstimator];
						Abort[];
					];
					estimateStDv = Estimate[stdvEstimator, s];
					
					ratio = (1-acceptRatio)/(2-acceptRatio);
					alpha = 4 * acceptRatio * ratio^2;

				];				  	
				printablestats = acceptStats;				
				(* do things that happen some multiple of every theta-update-cycle updates *)
				If[Mod[completedMoves, movesPerThetaUpdt]==0,
					thetas = Map[updateTheta[#, acceptStats]&, thetas];							
					acceptStats = emptyStats;
				];
				
				(* might as well just always do this stuff *)
				Module[{lowestEnergy, lowestEposit, lowestEState},				 
					If[MatchQ[energy, _List],				
						lowestEnergy = Min[energy];
						lowestEposit = Position[energy, lowestEnergy][[1,1]];
						lowestEState = Part[state,lowestEposit];,
						(* else *)
						lowestEnergy = energy;
						lowestEState = state;
					];
					If[lowestEnergy < bestEnergyEver,
						bestEnergyEver = lowestEnergy;
						bestStateEver = lowestEState;
						report[bestStateEver, bestEnergyEver, (totalNumMoves)];
					];
				];
				
				If[debug,
					Print["iter<",completedMoves,"> mean<", mean, "> dM<", deltaMean, ">" ];
					Print["S<", s, "> aR<", acceptRatio, ">"];
					Module[{ma, la},
						ma = Fold[If[#1[[2, 1]]/Max[1,#1[[2, 2]]] > #2[[2, 1]]/Max[1,#2[[2, 2]]], #1, #2] &, First[printablestats], Rest[printablestats]];
						la = Fold[If[#1[[2, 1]]/Max[1,#1[[2, 2]]] < #2[[2, 1]]/Max[1,#2[[2, 2]]], #1, #2] &, First[printablestats], Rest[printablestats]];
						Print["most <", ma, ">\nleast<", la, ">"];
					];
				];	


			];
		];
		Module[{lowestEnergy, lowestEposit, lowestEState},				 
			If[MatchQ[energy, _List],				
				lowestEnergy = Min[energy];
				lowestEposit = Position[energy, lowestEnergy][[1,1]];
				lowestEState = Part[state, lowestEposit];,
				(* else *)
				lowestEnergy = energy;
				lowestEState = state;
			];
			{{lowestEnergy, lowestEState}, {bestEnergyEver, bestStateEver}, totalNumMoves}			
		]
	];
	
Report[state_, energy_, iter_] :=
	Module[{},
		Print["New best energy ever: <", energy, "> at iter <",iter,">\n<", state, ">"];
	];

NewEstimator[prod_, freq_, lambda_, invY_, pA_] := 
	Module[{y,z},
		y = 1/invY;
		z = (freq/(prod/lambda));
		Estimator[Max[10^-10, 1-z], 1 , 0 , 0, 0, y, y^2 , pA, y]
	]; 

parameterA[e_Estimator] := e[[8]];
parameterB[e_Estimator] := e[[9]];

InitializeThetas[names_, ranges_] := 
	MapThread[#1->(DefaultPercentOfRange * (Max[#2]-Min[#2]))&, {names, names /. ranges}];

ValidlyInitializeVariables[variableNames_, variableRanges_, constraint_, tries_] := 
	Module[{state, tried=0},
		state = InitializeVariables[variableNames, variableRanges];
		While[tried < tries && !(constraint /. state), 
			state = InitializeVariables[variableNames, variableRanges];
		];
		If[tried == tries && !(constraint /. state),
			Message[ldsa::init, tries]
		];
		state		
	];	
	 
InitializeVariables[variableNames_, variableRanges_] := 
	Map[InitialRandomValue[#, (# /. variableRanges)]&, variableNames];	

InitialRandomValue[name_, {min_, max_}] := (name -> RandomReal[UniformDistribution[{min,max}]]);

ReinitzUpdate[name_, current_, theta_] :=
	current + (RandomChoice[{1,-1}]*theta*Log[RandomReal[]]);
	
Estimate[estimator_Estimator, temp_] := 
	Module[{denom},
		denom =  ((parameterA[estimator] * temp) + parameterB[estimator]);
		If[denom == 0,
			Print["got denominator 0. temperature is <", temp ,"> here is the estimator"];
			Print[FullForm[estimator]];
			10^30 (* i dunno, some big value? *), 

			(* else *)
			1/denom,

			(* if the comparison fails, then ... yeah, i've got no idea then. *)
			Print["could not compare denominator <", denom, "> to 0. temperature is <", temp ,"> here is the estimator"];
			Print[FullForm[estimator]];			
			-1
		]		
	];

	
UpdateEstimator[estimator_Estimator, temp_, invY_] := 
	Module[{y, weight, sum, sumXX, sumXY, sumX, sumY, sumYY, pa, pb},
		y = 1/invY;
		weight = estimator[[1]];
		sum    = estimator[[2]] * weight + 1 ;
		sumXX  = estimator[[3]] * weight + temp^2;
		sumXY  = estimator[[4]] * weight + temp*y;		
		sumX   = estimator[[5]] * weight + temp;
		sumY   = estimator[[6]] * weight + y;				
		sumYY  = estimator[[7]] * weight + y^2;
		If[MemberQ[Map[realQ[#]&, {sum, sumXX, sumXY, sumX, sumY, sumYY}], False], 
			Print["boo: ", {sum, sumXX, sumXY, sumX, sumY, sumYY, temp, invY}];
			Abort[];
		];
		pa = (sum * sumXY - sumX * sumY)/Max[(sum * sumXX - sumX^2),10^-30];
		pb = (sumY - pa * sumX)/sum;
		Estimator[weight, sum, sumXX, sumXY, sumX, sumY, sumYY, pa,	pb]
	];

	
forgettingPeriod[startLocalState_, numForgottenMovesEach_, variableNames_, moveGenerator_, thetas_, 
				 metric_, constraint_, talkativeErrs_, talkativeInit_, 
				 forgettingInvTemp_, stats_] :=
    Module[ {iter, whoseTurn, oldVal, newVal, var, 
    		 localstats, varStat, nVarStat, 
    	     plocalState, localState, plocalEnergy, localEnergy, deltaEnergy},
    	localstats = stats;    	     
        localState = startLocalState;
        localEnergy = metric /. localState;
        For[iter = 0, iter < numForgottenMovesEach, iter++, 
            whoseTurn = Mod[iter, Length[variableNames]] + 1;
            var = variableNames[[whoseTurn]];
            varStat     = var /. localstats;
            nVarStat    = AddTrial[varStat];            
            oldVal =  var /. localState;            
            newVal = moveGenerator[var, oldVal, var /. thetas];
            plocalState = ReplaceAll[localState, (var -> oldVal)  -> (var -> newVal)];
            localstats  = ReplaceAll[localstats, (var -> varStat) -> (var -> nVarStat)];
            If[ constraint /. plocalState,
                plocalEnergy = metric /. plocalState;
                If[ (talkativeErrs && !NumberQ[localEnergy]),
                    Print["Got an energy that wasn't a number (during forgetting)."];
                    Print["Tried to change <", variableNames[[whoseTurn]], "> from <", oldVal, "> to <", 
                    	  newVal,">. Got <",plocalEnergy,">"];
                ];
                deltaEnergy = plocalEnergy - localEnergy;
                If[ deltaEnergy <= 0 || (deltaEnergy*-forgettingInvTemp > Log[$MinNumber] && RandomReal[] <= E^(deltaEnergy*-forgettingInvTemp)),                
                    localState = plocalState;
                    localEnergy = plocalEnergy;
		            localstats  = ReplaceAll[localstats, (var -> nVarStat) -> 
		            	                                 (var -> AddSuccess[nVarStat])];    
                ];
            ];
(*            If[ talkativeInit && Mod[iter, Floor[numForgottenMovesEach/10]]==0,
                Print["Working hard. Forgetting iter: ", iter];
            ];*)
        ];
        If[ talkativeInit,
            Print["Working hard. Forgetting iter: ", iter];
            If[numForgottenMovesEach == 0,
            	Print["Dude, check your damn math."];
            ];
        ];
        {localState, localEnergy, localstats}
    ];

initializationPeriod[startLocalState_, numInitMovesEach_, variableNames_, moveGenerator_, thetas_, 
					 metric_, constraint_, talkativeErrs_, talkativeInit_, 
					 initializationInvTemp_, stats_] :=
    Module[ {iter, whoseTurn, oldVal, newVal, lSum, lSqrdSum, localstats, 
    	     var, varStat, nVarStat, deltaEnergy, numAccepted,
    	     plocalState, localState, plocalEnergy, localEnergy},
    	localstats = stats;
        lSum = lSqrdSum = 0;
        localState = startLocalState;
        localEnergy = metric /. localState;        
        numAccepted = 0;
        For[iter = 0, iter < numInitMovesEach, iter++,
            whoseTurn   = Mod[iter, Length[variableNames]] + 1;
            var         = variableNames[[whoseTurn]];
            varStat     = var /. localstats;
            nVarStat    = AddTrial[varStat];
            oldVal      = var/. localState;            
            newVal      = moveGenerator[var, oldVal, var /. thetas];
            plocalState = ReplaceAll[localState, (var -> oldVal)  -> (var -> newVal)  ];
            localstats  = ReplaceAll[localstats, (var -> varStat) -> (var -> nVarStat)];
            If[ constraint /. plocalState,
                plocalEnergy = metric /. plocalState;
                If[ (talkativeErrs && !NumberQ[localEnergy]),
                    Print["Got an energy that wasn't a number (during initialization)."];
                    Print["Tried to change <", variableNames[[whoseTurn]], "> from <", oldVal, 
                    	  "> to <", newVal,">. Got<",plocalEnergy,">"];
                ];
                deltaEnergy = plocalEnergy - localEnergy;

                If[ deltaEnergy <= 0 || (deltaEnergy*-initializationInvTemp > Log[$MinNumber] && RandomReal[] <= E^(deltaEnergy*-initializationInvTemp)),
                    localState  = plocalState;
                    localEnergy = plocalEnergy;
		            localstats  = ReplaceAll[localstats, (var -> nVarStat) -> 
		            	                                 (var -> AddSuccess[nVarStat])];                    
		            numAccepted++;
                ];
            ];
            lSum += localEnergy;
            lSqrdSum += localEnergy^2;
(*            If[ talkativeInit && Mod[iter, Floor[numInitMovesEach/10]]==0,
                Print["Working hard. Iter: ", iter];
            ];*)
        ];
        If[ talkativeInit,
            Print["Working hard. Iter: ", iter];
        ];
        {localState, lSum, lSqrdSum, localEnergy, localstats, numAccepted}
    ];	
			
sweep[startLocalState_,  startLocalEnergy_, startTemp_, meanEstimator_, stdvEstimator_, 
	  variableNames_, moveGenerator_, thetas_, constraint_, metric_, talkativeErrs_, 
	  startDS_, startNumTemp_, tempUpdateFreq_, lambda_, alpha_, estimateUpdateFreqEach_, emn_,
	  stats_, numProcessors_] :=
	Module[{estmtCounter, propState, newEnergy, energyChange, mean, variance, ratio, 
		    oldVal, newVal, whoseTurn, lSum, lSqrdSum, 
		    localState, localEnergy, lNumAccepted, s, estimateMean, estimateStDv, 
		    deltaS, numTemp,
		    localstats, var, varStat, nVarStat},	     	

    	localstats = stats;		    
	    s = startTemp;
	    deltaS = startDS;
	    numTemp = startNumTemp;	    	    
		estimateMean = emn;
		localState = startLocalState;
		localEnergy = startLocalEnergy;
		lSum = lSqrdSum = lNumAccepted = 0;
		
		For[estmtCounter=0, estmtCounter < estimateUpdateFreqEach , estmtCounter++,
			whoseTurn = Mod[estmtCounter, Length[variableNames]] + 1;					

            var         = variableNames[[whoseTurn]];
            varStat     = var /. localstats;
            nVarStat    = AddTrial[varStat];
            oldVal      = var/. localState;            
            newVal      = moveGenerator[var, oldVal, var /. thetas];
            propState   = ReplaceAll[localState, (var -> oldVal)  -> (var -> newVal)  ];
            localstats  = ReplaceAll[localstats, (var -> varStat) -> (var -> nVarStat)];			
			
			If[ constraint /. propState,
				newEnergy = metric /. propState;
				If[!NumberQ[newEnergy],
					If[talkativeErrs,
						Print["Got an energy that wasn't a number (during annealing).\nTried to change <", 
							  var, "> from <", oldVal, "> to <", newVal,">."];
					], 
					(* else *)
					energyChange = newEnergy - localEnergy;
					If[energyChange < 0 || (-s*energyChange > Log[$MinNumber]  && E^(-s*energyChange) > RandomReal[]),
						localState = propState;
						localEnergy = newEnergy;
						localstats  = ReplaceAll[localstats, (var -> nVarStat) -> 
		            	                                     (var -> AddSuccess[nVarStat])];         						
						lNumAccepted++;	
					];	
				];						
			];
			lSum += localEnergy;
			lSqrdSum += (localEnergy - estimateMean)^2;
			If[MatchQ[lSqrdSum,Indeterminate],
				Print["booooooooooooooooooo"];
				Print[lSum, " ", localEnergy, " ", estimateMean];
				Abort[];
			];
			If[--numTemp <= 0,
				s += (deltaS * numProcessors);
				estimateMean = Estimate[meanEstimator, s];
				estimateStDv = Estimate[stdvEstimator, s];
				deltaS = lambda * (alpha / (s^2 * estimateStDv^3)) * tempUpdateFreq;
				numTemp = tempUpdateFreq;
			];			
		];
		{localState, lSum, lSqrdSum, localEnergy, lNumAccepted, s, deltaS, localstats}		
	];		

AddTrial[Stats[successes_, trials_]] := Stats[successes, trials+1];
AddSuccess[Stats[successes_, trials_]] := Stats[successes+1, trials];
AccumulateStats[stl_] := 
	Fold[(Stats[#1[[1]] + #2[[1]], #1[[2]] + #2[[2]]]) &, Stats[0, 0], stl];
aggregateParallelStats[vrblNm_, statsList_] := 
	Module[{allstats}, 
		allstats = Map[(vrblNm /. #)&, statsList];
		(vrblNm -> AccumulateStats[allstats])
	];	
	
realQ[x_] := NumberQ[x] && ! MatchQ[x, _Complex];	
	
End[]; (* End Private Context *)

EndPackage[];