// ***************** COPYRIGHT (c) 2010 STEFAN WANER ******************
// *********************** ALL RIGHTS RESERVED ************************

// ***TODO when SSE isNaN we need a message to increase grad steps etc.

//***HERE TODO IN CANCAVITY TRY TAKING DOUBLE STEPS ONCE k IS SELECTED
// *** perhaps mutliply by factor of 1.5 each time...
// ** arithmetic pregression here...

// 1. convert gradient function to balanced dq



// ********** Gradient Descent Utilities *****************
// ***** temporary testing input data:
var $1, $2, $3, $4, $5, $6, $7, $8;// must be global...
var xValsCF = [0], yValsCF = [0];
// xValsCF = [0,0, 1, 2, 3, 4, 5, 6, 7, 8, 9];
// yValsCF = [0,24, 28, 32, 37, 42, 48, 54, 59, 63, 65];


var numParamsCF = 0;
var CFnumPoints = 10;
var theVariablesCF = []; // names of the variables
var theModelCF = ""; // will be the curve model
var thePhaseCF = 1;
var fitSuccessCF = false;
var numCyclesCF = 0; // number of fit attempts

// ** some monitored stats vars

var numAtitkenCD = 0;
var numJumpsCF = 0;
var numJumpTriesCF = 0;
var lastJumpAtCF = 0; // step of last jump
var numGD = 0;
var jumpFracCF = .5 // will jump with a prob of .5 each step
var foundStartCF = false;
var cutsCounter = 0;
var littleJumpsCF = 0;


var locationCF = []; // location
var locPrevCF = [], locPrev2CF = []; // previous and one before
var locAitken = []; // to try for Aitken
var bestLocationCF = []; // for best of worst in multiple attempts


var currentFnCF = ""; // the function with the variables set to numbers
var hCF = .00001 // for computing derivatives
var kCF = 1; //for testing concavity
var minSSE = infinity;
var prevSSE = 0;
var bestPointCF = [];
var okToRollCF = true;
var abendCF = false;
var displayAccuracyCF = 6; // for showing params etc
var accCF = 8; // number sig digits for phase 1 SSE convergence accuracy
var acc2CF = 8; // sig digits for end result convergence
var isCloseEnuffCF = false;
var furtherFactor = .75; // extending GD step factor gtrial

// various maxima
var numAttemptsCF = 200; // number of trials for starting point
var maxNumGradSteps = 120; // max # successive grad descent steps
var maxNumPhase1 = 20; // max # successive grad descent steps
var maxNumIteratesCF = 20;
var mustGuessCF = true;

var maxNumCutsCF = 15;

var testFunctions = ["$1 + $2*x", "$1 + $2*x + $3*x^2", "($1-$3)*exp(-$2*x)+$3", "$1/(1+$2*$3^(-x))", "$1*$2^x", "$1 + $2*$3^x" , "$1*ln(x) + $2", "$1*sin(x-$3)+$4"]
// newton, logistic
// logistic: best SSE = 7.590133 params = 78.690723,2.4716967,1.3190381
// another: Best SSE = 5.90095; params = 87.0824,2.71888,1.27283

var dataSamples = [];
dataSamples[0] = ' 0, 210\n 1, 198\n 2, 187\n 5, 155\n 10, 121\n 15, 103\n 30, 77\n 45, 71\n 60, 70';
dataSamples[1] = ' 0, 24\n 1, 28\n 2, 32\n 3, 37\n 4, 42\n 5, 48\n 6, 54\n 7, 59\n 8, 63\n 9, 65';
dataSamples[2] = '5, 199\n6, 210\n7, 222\n8, 235\n9, 250\n10, 268\n11, 271\n12, 265\n13, 272\n14, 275\n15, 289\n16, 299\n17, 300\n18, 304\n19, 309';
dataSamples[3] = '0, 125\n1, 153\n2, 170\n3, 180\n4, 167\n5, 140\n6, 130\n7, 137\n8, 158\n9, 180\n10, 180\n11, 151\n12, 131\n13, 126\n14, 134\n15, 168';

// *** alarms
var dataAlarmCF = false; // if no data to work with
var modelAlarmCF = false; // if no data to work with

// messages
var tryAgainMsg = "Sorry. I could not fit this data to this model to the precision required. Either lower the precision, supply another starting point, or start again by pressing 'Fit Curve.'"
var pressImproveMsg = "The parameters don't seem to be converging well. Press 'Improve fit' a few times to remedy this, or else just start again by pressing 'Fit Curve.'"
var pressWeakImproveMsg = "The parameters may not have converged completely. Press 'Improve fit' a few times to remedy this."
// *** end of globals
var noDataMsg = "You must enter data below to fit.";
var badDataMsg = "Something is wrong with the format of the data points you entered. Format should be  x, y on each line (see the examples)).";
var noModelMsg = 'You must first enter a model. Press "Examples" to see some.'

if (theLanguage == "es") {

noDataMsg = "Debe primero ingresar datos para ajustar.";
badDataMsg = "Hay algo equivocado conj el formato de los puntos datos. El formato debe ser x, y en cada fila (vea los ejemplos))."
tryAgainMsg = "Disculpe. No puedo ajustar estos dataos a este modelo con la precisión requisita. Baje la precisión, ingrese una nueva adivina, o pruebe de nuevo por pulsar 'Ajustar curva.'";
noModelMsg = 'Debe primero ingresar un modelo. Pulse "Ejemplos" para ver algunos modelos.'
}

// *** output 
var theResultsCF = "<span style = 'width: 50%; border-style:solid; border-width:1px; border-color: goldenrod'>";





// *** globals for Pezzulo's code
var stErCF = []; // st error of coeff
var thePInputCF = "";
var CR = unescape("%0D");
var LF = unescape("%0A");
var Tb = unescape("%09");
var NL = CR + LF;
function createArray() {
	this.length = createArray.arguments.length
	for (var i = 0; i < this.length; i++) { 
		this[i] = createArray.arguments[i] 
		}
	}
var SEP = new createArray(1,1,1,1,1,1,1,1);
var Der = new createArray(0,0,0,0,0,0,0,0,0);
var Arr = new createArray(0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0);
var Cov = new createArray(0,0,0,0,0,0,0,0,   0,0,0,0,0,0,0,0,   0,0,0,0,0,0,0,0,   0,0,0,0,0,0,0,0,   0,0,0,0,0,0,0,0,   0,0,0,0,0,0,0,0,   0,0,0,0,0,0,0,0,   0,0,0,0,0,0,0,0  );
var xArr = new createArray(0,0,0,0,0,0,0,0,0);

function roundCF(inArr) {
var outArr = []
for (var i = 0; i <= inArr.length-1; i++) outArr [i] = (inArr[i]).toPrecision(displayAccuracyCF);
return(outArr);
	}

function doOneGD() {
var maxNumStepssave = maxNumGradSteps;
maxNumGradSteps = 1;
gradientDescent ();
maxNumGradSteps = maxNumStepssave;
}


// *** main routine here
function multipleFits() {
var bestSSE = infinity;
var maxFitTriesCF = 10;
fitSuccessCF = false;
numCyclesCF = 0;
var theBestResultsCF = "";
theResultsCF = '';
dataAlarmCF = false; modelAlarmCF = false; 
// turn off alarms - they get turned on again in initiate if necesary
while ((numCyclesCF < maxFitTriesCF) && (!fitSuccessCF) && (!dataAlarmCF) && (!modelAlarmCF)) {
	numCyclesCF ++;

	mustGuessCF = false;
	thePhaseCF = 1;
	if (stripSpaces(document.theFormP.params.value) == '') { mustGuessCF = true;}
	try{
		fitTheCurve();
		}
	catch(error) {
		stripSpaces(document.theFormP.params.value) == '';
//testDisplay("*****ERRRROORRR*****")
		}
	if ((minSSE < bestSSE)&&(isGood(locationCF)) ) {
		theBestResultsCF = theResultsCF;
		bestSSE = minSSE;
		for (var i = 0; i <= numParamsCF-1; i++) bestLocationCF[i] = locationCF[i];
			} // if better
	} // while
if (numCyclesCF >= maxFitTriesCF) {
	//show the best one
	document.getElementById("fitResults").innerHTML = theBestResultsCF;
	document.theFormP.sse.value = bestSSE.toPrecision(acc2CF);
	showParams(bestLocationCF);
	
	}
else {
	document.getElementById("fitResults").innerHTML = theBestResultsCF;
	}
} // multipleFits


function fitTheCurve() {

initiateCF();
if (!okToRollCF) return(false);
	
gradientDescent ();

for (var i = 1; i <= maxNumIteratesCF; i++) {
	if (okToRollCF) {
		iterateCF();
		moveTo(locationCF);
		prevSSE = minSSE;
		minSSE = SSECF(locationCF);
		document.theFormP.sse.value = minSSE.toPrecision(acc2CF);;
		if ((closeEnough()) || (!okToRollCF)) i = maxNumIteratesCF + 1;
		} // if OK
	} // i

if (!isNaN(minSSE) && (document.theFormP.sse.value != "infinity")) {
	fitSuccessCF = true;
	// do the next step at the end only to get the actual function
	substituteParamsCF(displayAccuracyCF);
	// now display results
	showParams(locationCF);
	//search for a "yi" to put the equation
	// this should be in the form (default y1);
	var jCF = 1;
	theInstruction = "document.theFormA.y"+ jCF +".value = currentFnCF;";
	doIt = eval(theInstruction);

	plotCF();
//testDisplay("********** Found it on attempt # " + numCyclesCF)
	} // if you got a good SSE
else {
//testDisplay("****FAILED TRYING AGAIN")
	document.theFormP.params.value = "" // must start fresh


}


if (!fitSuccessCF) alert(tryAgainMsg);



} // fitTheCurve

function closeEnough() {
var hasConverged = true;
if (prevSSE.toPrecision(acc2CF) != minSSE.toPrecision(acc2CF) ) hasConverged = false;
else {
	for (var i = 0; i <= numParamsCF-1; i++)  if(locPrevCF [i].toPrecision(acc2CF)!= locationCF[i].toPrecision(acc2CF)) hasConverged = false;
	}
return(hasConverged)
} // closeEnough



function abandonShip () {
// move to prev location

abendCF = true;
//testDisplay("\n !!!Interate abend!!!")
document.theFormP.sse.value = "infinity";
document.theFormP.params.value = "";
okToRollCF = false;
//moveTo(locPrevCF);
// document.theFormP.sse.value = SSECF(locationCF).toPrecision(acc2CF);
//substituteParamsCF(displayAccuracyCF);
//plotCF();
} // 

function initiateCF() {
isBadSSE = false;
okToRollCF = true;
abendCF = false;
locationCF = []; locPrevCF = []; locPrev2CF = [];locAitken = [];bestLocationCF = [];
maxNumGradSteps = parseInt(document.theFormP.gradsteps.value);
// alert(maxNumGradSteps);
//document.resultsForm.output.value = '';
dataAlarmCF = false;
if (stripSpaces(document.theFormP.thePoints.value) == '') {
	okToRollCF = false; 
	alert(noDataMsg); 
	dataAlarmCF = true;
	return(false);
	}

okToRollCF = true;
abendCF = false;
maxNumGradSteps = parseInt(document.theFormP.gradsteps.value);
// alert(maxNumGradSteps);
// document.resultsForm.output.value = '';
readInTextarea();
if (okToRollCF == false) {alert(badDataMsg); return(false)};

// now grab the globals the above function gives:
CFnumPoints = numPoints;
xValsCF = []; yValsCF = [];
for (var i = 1; i <= numPoints; i++) {
	xValsCF[i] = thePlottedPoints[i][1];
	yValsCF[i] = thePlottedPoints[i][2];
	} // i
xValsCF[0] = xValsCF[1] 
yValsCF[0] = yValsCF[1]
// only for computing .min later
// the zeroth coordinate is not used in the curve fit
// initialize the logging counts:
//testDisplay(xValsCF);
//testDisplay(yValsCF);
//alert("here")

numJumpsCF = 0;
numJumpTriesCF = 0;
lastJumpAtCF = 0;
cutsCounter = 0;
isCloseEnuffCF = false;
// thePhaseCF = 1; 
littleJumpsCF = 0;
littleJumpsTryCF = 0;
improvingCF = false;
minSSE = infinity;
numGD = 0; // number of grad descent steps
// get the curve equation from the input field
modelAlarmCF = false;
if (stripSpaces(document.theFormP.model.value) == '') {
	okToRollCF = false; 
	alert(noModelMsg); 
	modelAlarmCF = true;
	return(false);
	}
theModelCF = document.theFormP.model.value;
// get the number of parameters

var len = theModelCF.length;
theVariablesCF = []; // reinitialize this;
var j = 0;
var theIndex = 0;
for (var i = 0; i <= len; i++) {
	if (theModelCF.charAt(i) == "$") {
		if (isNumberChar(theModelCF.charAt(i+2) )) theIndex  = theModelCF.substring(i+1,i+2);
		else theIndex  = theModelCF.charAt(i+1);
		var test = (theVariablesCF.join("")).indexOf(theIndex);
		if (test == -1) { theVariablesCF[j] = "$"+ theIndex; j++}
		i++;
		}
	} // i

numParamsCF = theVariablesCF.length;
// now supply an initial guess for the location

var startingAt = findStartPt();

// prepare points data for pelluzo's algorithm
thePInputCF = "";
for (var i = 1; i <= CFnumPoints; i++) {
	thePInputCF += xValsCF[i] + "," + yValsCF[i]
	if (i < CFnumPoints) thePInputCF += CR;
	} // i


} // initiateCF


function gradientDescent() {
// with random tests between each step

while ((numGD <= maxNumGradSteps) && (!isCloseEnuffCF)) {
if (numGD > maxNumPhase1) thePhaseCF = 2;
gradientDescentStep ();
if (Math.random() < jumpFracCF) var jt = jumpTry(1);
numGD += 1;
	} // while

displayGDREsults()
} // gradientDescent

function displayGDREsults() {

// *** Display Results
showParams(locationCF);


//testDisplay("***DID " + numGD + " **GRAD DESCENTS/n Number of jumps = " + numJumpsCF + " out of " + numJumpTriesCF + "tries.\n Last jump at step " + lastJumpAtCF + "\n cut grad stepsize " + cutsCounter + " times." + "\n Aitkensworked in concave down regions " + numAtitkenCD + " times." + "\n Number of little jumps = " + littleJumpsCF + " in " + littleJumpsTryCF + " tries." + "\n Best SSE was " + minSSE);
}

function gradientDescentStep() {
if (okToRollCF) {
// test concavity and move to low point if concave up
// if not, takes bigger and bigger steps downward to a local min
// testDisplay("SSE CURRENTLY " + f0);
	var theGrad = negGrad(locationCF);

		// first take a step in the direction of the neg gradient:
		// but step size should be inverse prop to the gradient here
		var kNew = 0.5/matrixNormSq (theGrad);
		var newPositionCF = matrixSum (locationCF, scalarMult (kNew, theGrad));
		// now adjust by more steps if necessary
		var testSSE = SSECF(newPositionCF);
		if (testSSE < minSSE) {
			// go there or further
			minSSE = testSSE;
			var thePosition = jiggleK(kNew,theGrad, newPositionCF);
			moveTo(thePosition);
			
			}
		else {
			var thePosition = jiggleK(hCF*kNew,theGrad, locationCF);
			moveTo(thePosition);
			testSSE = SSECF(thePosition);
			if (testSSE < minSSE) {
				
				}
			else {
				littleJumpsTryCF ++;
				for (var i = 100; i >=1; i--) {
					var jt = jumpTry (i*hCF);
					if (jt) { i = 1; littleJumpsCF ++}
					}
				if (!jt) {
//testDisplay("**IN TROUBLE HERE FINDING A STEP**");
				isCloseEnuffCF = true;
				// get out of here and hope for the best...
					}
				}

			}

// * try Aitken to see if that is perhaps even better
			
			var AtestSSE = SSECF(locAitken);
			if (AtestSSE < minSSE) {
				moveTo(locAitken);
				minSSE = AtestSSE;
				numAtitkenCD++;
				} //
	


		updateSSE();
//testDisplay("\n ***CONCAVE DOWN SSE =" + minSSE);
// document.theFormP.params.value = roundCF(locationCF);
// document.theFormP.sse.value = minSSE;
	
		
	// newPositionCF = jiggleK(kNew,theGrad,newPositionCF);


		


// document.theFormP.params.value = roundCF(locationCF);
updateSSE(minSSE)


} //if okToRollCF

//testDisplay("\n SSE =" + minSSE);
} // gradientDescentStep


function jumpTry (theJumpFactor) {
numJumpTriesCF++;
var presentLoc = [];
var xPick = Math.random();
var jumpSize = theJumpFactor *2* xPick/(1-Math.abs(xPick));
for (var i = 0; i <= numParamsCF-1; i++)  presentLoc[i] = locationCF[i];
jiggleCoords(presentLoc, jumpSize);
var testSSE = SSECF(presentLoc);
if (testSSE < minSSE) {
	moveTo(presentLoc);
	prevSSE = minSSE;
	minSSE = testSSE;
	numJumpsCF ++;
	lastJumpAtCF = numGD;
	return(true);
// testDisplay("\n Jumped to better SSE: " + minSSE);
	}
else {
// testDisplay("\n attempted to jump but SSE was: " + testSSE);
	return(false);
}
} //jumpTry

function negGrad(inArr) {
// computes the *unit* negative gradient at inArr
// returns an array 
// if error returns 666 in the numParamsCF-position
var outArr = [];

var f0 = 0, f1 = 0;
// old one-ord diff q approx
// f0 = SSECF(inArr);
// for (var i = 0; i <= numParamsCF-1; i++) {
//	inArr [i] += hCF;
//	if (i > 0) inArr [i-1] -= hCF;
//	f1 = SSECF(inArr);
//	outArr[i] = (f0-f1)/hCF; // note negative grad
//	} // i

// balanced dq here:
for (var i = 0; i <= numParamsCF-1; i++) {
inArr [i] += hCF;
f1 = SSECF(inArr);
inArr [i] -= 2*hCF;
f0 = SSECF(inArr);
inArr [i] += hCF; // back to normal
outArr[i] = (f0-f1)/(2*hCF); // note negative grad
} // i


// normalize -- dont do that
//var norm = 0;
//for (var i = 0; i <= numParamsCF-1; i++) {
//	norm += outArr[i]* outArr[i];
//	} // i
//norm = Math.sqrt(norm);
//for (var i = 0; i <= numParamsCF-1; i++) {
//	outArr[i]= outArr[i]/norm;
//	} // i


return(outArr);

//	}
// catch(error)
//	{
// alert("***ERROR**")
//	outArr[numParamsCF] = 666;
//	}

} // grad


function initializeFunction(inArr){
// sets each variable equal to its value in inArr
// currentFnCF is the global string that results

//alert(inArr);
var theCommand = "";
for (var i = 0; i <= inArr.length-1; i++) {
	theCommand += theVariablesCF [i] + " = parseFloat(" + inArr [i] + ");";
		} // i
var doIt = eval(theCommand); 
// now have function with numbers


//alert(currentFnCF);
} // initializeFunction


function SSECF(thePlace) {
// compute SSE at thePlace
// if thePhaseCF = 1 it actually Sums 4th powers
initializeFunction(thePlace);
var theResult = 0;
var theResidual = 0;
for (var i = 1; i <= CFnumPoints; i++) {
	x = parseFloat(xValsCF[i]);
	theResidual = myEval(theModelCF)-yValsCF[i];
	var theIncr = theResidual*theResidual;
	if (thePhaseCF == 1) theIncr *= theIncr; // sum of 4th powers
	theResult += theIncr;
	

	} // i
return(theResult);
} // SSE




function findStartPt() {
// start at 0 unless have somewhere already
// best not to do a random walk...
var numTried = 0; var f0 = 0;
var atPt = [], bestPt = [];

// alert(mustGuessCF);
if (!mustGuessCF) {
		var theString = document.theFormP.params.value;
		theString = replaceChar(theString,comma,tab);
		var thePt = parser(theString, tab);
		for (var i = 0; i <= numParamsCF-1; i++) {
			
			bestPt[i] = parseFloat(thePt [i+1]);
			atPt[i] = bestPt[i];
//testDisplay( "We are the following point " + bestPt[i]);
			} // i	 
	f0 = SSECF(bestPt);
// testDisplay("started SSE being " + f0);
	if ((isNaN(f0)) || (f0 == infinity) || (f0 == -infinity)) {
		for (var i = 0; i <= numParamsCF-1; i++) { 	
			atPt[i] = 0; bestPt[i] = 0;
			}// i
		} // bad start pt
	} // if a pt already there
else {
thePhaseCF = 1;
	for (var i = 0; i <= numParamsCF-1; i++) { 	
		atPt[i] = 0; bestPt[i] = 0;
		}// i
	}
while ((okToRollCF) && (numTried < numAttemptsCF)) {
	numTried += 1;
	try {



	f0 = SSECF(atPt);
// testDisplay(numTried + ":" + atPt + "SSE = " + f0);

	if (isNaN(f0)) {
//testDisplay("hit NAN");
		// this means that we are at a bad pt
		// move a little more agressively to get out of there:	
		var trying = true, theNum = 0;
		while ((trying) || (theNum < numAttemptsCF)) {
			theNum += 1;
			var xPick = Math.random();
			var jumpSize = 2* xPick/(1-Math.abs(xPick));
			atPt = jiggleCoords(atPt, jumpSize);
			f0 = SSECF(atPt);
			if (!isNaN(f0)) {
				trying = false;
// testDisplay("FOUND GOOD POINT AT " + atPt);
				}
			if ((theNum == numAttemptsCF) && (trying)) {
//testDisplay("FAILED TO FIND START PT.");
 
				okToRollCF = false;
				} // failed
			} // while		
		} // hit a bad point

	if (f0 < minSSE) {
			for (var i = 0; i <= numParamsCF-1; i++) bestPt[i] = atPt[i];
			minSSE = f0;
			} 
		} // try
	catch(error) {

//testDisplay("ERROR")

		}
	// now jiggle coords

	pickRandomPt(atPt);

	} // while
if (minSSE == infinity) okToRollCF = false;
document.theFormP.params.value = roundCF(bestPt);
document.theFormP.sse.value = minSSE.toPrecision(acc2CF);;
// ** testing
//testDisplay("Start Point is :\n" + roundCF(bestPt) + "\n Starting SSE = " + minSSE);
// ** end testing

// set all three locations to current loc
for (var i = 0; i <= numParamsCF-1; i++) locPrevCF [i] = bestPt [i];
for (var i = 0; i <= numParamsCF-1; i++) locPrev2CF [i] = bestPt [i];
for (var i = 0; i <= numParamsCF-1; i++) locationCF [i] = bestPt [i];


foundStartCF = true;
return(bestPt);
} // findStartPt


function moveCF(inArr, jumpSize) {
// move to location by a random jump of jumpSize in *every* coordinate
for (var i = 0; i <= numParamsCF-1; i++) inArr [i] += jumpSize*pickOne(-1,1);
return(inArr);
}

function pickRandomPt(inArr) {
var xPick = 0;
for (var i = 0; i <= numParamsCF-1; i++) { 
	xPick = Math.random();
	inArr[i] =  xPick/(1-Math.abs(xPick));
	}
} // pickRandomPt()


function jiggleCoords(inArr, jumpSize) {
// change each coord by 0 or +- jumpSize
for (var i = 0; i <= numParamsCF-1; i++) inArr [i] += jumpSize*pickOne(-1,0,1);
return(inArr);
}


function moveTo(inArr) {
for (var i = 0; i <= numParamsCF-1; i++) {
	locPrev2CF [i] = locPrevCF [i];
	locPrevCF [i] = locationCF [i];
	locationCF[i] = inArr[i];
	locAitken [i] = (locPrev2CF [i]* locationCF [i] - locPrevCF[i]* locPrevCF[i])/( locationCF[i]-2* locPrevCF[i]+ locPrev2CF[i]);
	} // i
} // moveTo



function reportResults() {


}

function matrixSum(A,B) {
var outArr = [];
for (var i = 0; i <= A.length-1; i++) {
	if ((typeof A[i] == 'number')&& (typeof B[i] == 'number')) outArr[i] = A[i]+B[i];
	} //i
return(outArr)
}

function scalarMult(k,A) {
var outArr = [];
for (var i = 0; i <= A.length-1; i++) {
	if (typeof A[i] == 'number') outArr[i] = k*A[i];
	} //i
return(outArr)
}

function matrixNormSq(A) {
var theSum = 0;
for (var i = 0; i <= A.length-1; i++) {
	if (typeof A[i] == 'number') theSum += A[i]*A[i];
	} //i
return(theSum)
}

function testDisplay(theString) {
// for testing
document.resultsForm.output.value += cr + theString;
}


// ********* plotting

function plotCF() {
var xMinCF = xValsCF.min();
var xMaxCF = xValsCF.max();
var dlx = .2*(xMaxCF- xMinCF)
//xMinCF -= dlx;
//xMaxCF += dlx;
var yMinCF = yValsCF.min();
var yMaxCF = yValsCF.max();

var dly = .2*(yMaxCF- yMinCF)
yMinCF -= dly;
yMaxCF += dly;
// now round to decent values if at least 10 units apart
var theLog = Math.floor( Math.log(xMaxCF - xMinCF)/Math.log(10) );
//alert(theLog)
var roundTo = Math.pow(10,theLog);
xMinCF = roundToNearest(xMinCF - .51*roundTo, roundTo);
xMaxCF = roundToNearest(xMaxCF + .51*roundTo, roundTo);
document.theFormB.a.value = xMinCF.toPrecision(3);
document.theFormB.b.value = xMaxCF.toPrecision(3);
document.theFormB.c.value = yMinCF.toPrecision(3);
document.theFormB.d.value = yMaxCF.toPrecision(3);
document.theFormB.autoYButton.checked = false;
// now graph the curve(s)
setUp();
readBasics(myGraph);
var L = readAndPlotCurves (myGraph);
cleanUp(myGraph) ;

} // plotCF


function substituteParamsCF(acc) {
// sets currentFnCF = actual function
currentFnCF  = theModelCF;
for (var i = 0; i <= numParamsCF-1; i++) {
	currentFnCF = replaceSubstring(currentFnCF, theVariablesCF[i], locationCF [i].toPrecision(acc));
	} // i
} // substituteParamsCF


function jiggleK(kIn,theGrad,theLocation) {
// used in concave down regions
// tried to extend stepsize by factors of 2
// up to a maximum of maxMoves times
var theJump = kIn; // start small
var maxMoves = 40, numMoves = 0;
var gettingBetter = true;
var goingPositive = false;
var currSSE = minSSE;
var thePosition1 = matrixSum (theLocation, scalarMult (theJump, theGrad));
var testSSE = SSECF(thePosition1);
	if (testSSE < currSSE) goingPositive = true;
	else {
//testDisplay("delta SSE is " + (testSSE-currSSE).toString() + "\n");
		thePosition1 = matrixSum (theLocation, scalarMult (-theJump, theGrad));
		testSSE = SSECF(thePosition1);
		if (testSSE < currSSE) goingPositive = false;
		else return(theLocation); // no better on either sides
		}
	
if (goingPositive == false) {
	theJump *= -1;
//testDisplay("Reversed Direction");
	}
while ((gettingBetter) && (numMoves < maxMoves)) {
	numMoves ++;
	theJump = theJump*2;
	thePosition1 = matrixSum (theLocation, scalarMult (theJump, theGrad));
	var cSSE = SSECF(thePosition1)
	if (cSSE < testSSE) {
		testSSE = cSSE;
		}
	else {
		// move back
		theJump = theJump/2;
		thePosition1 = matrixSum (theLocation, scalarMult (theJump, theGrad));
		gettingBetter = false; // get out
		}
	}
//testDisplay("Jiggled k " + (numMoves-1) + "times and k is " + theJump);
minSSE = testSSE;
return(thePosition1);

} // jiggleK()



// *** Pezzulo's code follows


function StudT(t,n) {
    t=Math.abs(t); var w=t/Math.sqrt(n); var th=Math.atan(w)
    if(n==1) { return 1-th/PiD2 }
    var sth=Math.sin(th); var cth=Math.cos(th)
    if((n%2)==1)
        { return 1-(th+sth*cth*StatCom(cth*cth,2,n-3,-1))/PiD2 }
        else
        { return 1-sth*StatCom(cth*cth,1,n-3,-1) }
    }


function AStudT(p,n) { var v=0.5; var dv=0.5; var t=0
    while(dv>1e-6) { t=1/v-1; dv=dv/2; if(StudT(t,n)>p) { v=v-dv } else { v=v+dv } }
    return t
    }
function AFishF(p,n1,n2) { var v=0.5; var dv=0.5; var f=0
    while(dv>1e-6) { f=1/v-1; dv=dv/2; if(FishF(f,n1,n2)>p) { v=v-dv } else { v=v+dv } }
    return f
    }


function Max(x1,x2) { return (x1>x2)?x1:x2 }
function Min(x1,x2) { return (x1<x2)?x1:x2 }

function Fmt(x) { var v;
	if(Math.abs(x)<0.00005) { x=0 }
	if(x>=0) { v='          '+(x+0.00005) } else { v='          '+(x-0.00005) }
	v = v.substring(0,v.indexOf('.')+5)
	return v.substring(v.length-10,v.length)
	}

function vFmt(x) { var v;
	if(Math.abs(x)<0.0000005) { x=0 }
	if(x>=0) { v='                   '+(x+0.0000005) } else { v='               '+(x-0.0000005) }
	//v = v.substring(0,v.indexOf('.')+7)
	//return v.substring(v.length-18,v.length)
// my amended for debuf
// uncomment the above two later
return(x.toPrecision(4))
	}

function Xlate(s,from,to) { var v = s;
	var l=v.indexOf(from);
	while(l>-1) {
		v = v.substring(0,l) + to + v.substring(l+1,v.length);
		l=v.indexOf(from)
		}
	return v
    }

function createArray() {
	this.length = createArray.arguments.length
	for (var i = 0; i < this.length; i++) { this[i] = createArray.arguments[i] }
	}

function ix(j,k) { return j*( numParamsCF +1)+k }






var i = 0; var j = 0; var k = 0; var l = 0; var m = 0;

var nVarCF = 1; // number of independent variables


var ccSW = 1;
var ccAv = 0;
var ccSD = 1;


function iterateCF() {

theResultsCF = "<div style = 'width: 50%; border-style:solid; border-width:1px; border-color: goldenrod'>";

var cccSW = 0;
var cccAv = 0;
var cccSD = 0;

dgfr = CFnumPoints - numParamsCF;
var St95 = AStudT( 0.05 , dgfr );
var Pctile = .5; // not used, I think



var da = Xlate(thePInputCF,Tb,",");
thePInputCF = da;
if( da.indexOf(NL)==-1 ) { if( da.indexOf(CR)>-1 ) { NL = CR } else { NL = LF } }

var vCentered = 1; // centering approx for derivs
var SSq = 0;
for (j = 0; j<numParamsCF*(numParamsCF+1); j++) { Arr[j] = 0; }
for (i = 1; i<=CFnumPoints; i++) {
	l = da.indexOf(NL); if( l==-1 ) { l = da.length };
	var v = da.substring(0,l);
	da = da.substring(l+NL.length,da.length);
	for (j = 0; j<nVarCF; j++) {
		l = v.indexOf(","); if( l==-1 ) { l = v.length };
		xArr[j] = eval(v.substring(0,l));
		
		v = v.substring(l+1,v.length);
		}
	x  = xArr[0]; 
	initializeFunction(locationCF);
	var yc = myEval(theModelCF);

	l = v.indexOf(","); if( l==-1 ) { l = v.length };
	y = eval(v.substring(0,l));
	v = v.substring(l+1,v.length);

	var vSEy = 1;
	yo = y;
	w = 1;
	var yT = y;

	y = yT;

	cccSW = cccSW + 1 / ( w * w );
	cccAv = cccAv + yc / ( w * w );
	cccSD = cccSD + ( y - ccAv ) * ( y - ccAv ) / ( w * w );

	for (var j=0; j<numParamsCF; j++) {
		var Save = locationCF[j]; if(Save==0) {var Del = 0.0001} else {var Del = Save/1000}
		locationCF[j] = Save + Del;
//if (isNaN(locationCF[j]))  testDisplay("/n" + "Step1 FAILED");
	initializeFunction(locationCF);
	ycIncr = myEval(theModelCF);

		Der[j] = ( ycIncr - yc ) / ( Del * w );
		if(vCentered=="1") {
			locationCF[j] = Save - Del;
			if (isNaN(locationCF[j]))  {
			okToRollCF = false;
			abandonShip ();
//testDisplay("/n" + "Step2 FAILED");
				} // if failed
			initializeFunction(locationCF);
			ycDecr = myEval(theModelCF);
			Der[j] = ( ycIncr - ycDecr ) / ( 2 * Del * w );
			}
		locationCF[j] = Save;
if (isNaN(locationCF[j]))  testDisplay("/n" + "Step3 FAILED");
		initializeFunction(locationCF);
		}
	Der[numParamsCF] = (y - yc) / w;
	SSq = SSq + Der[numParamsCF]*Der[numParamsCF];
	for (j=0; j<numParamsCF; j++) {
		for (k=0; k<=numParamsCF; k++) {
			Arr[ix(j,k)] = Arr[ix(j,k)] + Der[j] * Der[k]
			}
		}
	var SEest = 0;
	for (j=0; j<numParamsCF; j++) {
		SEest = SEest + Cov[ix(j,j)] * Der[j] * Der[j];
		for (k=j+1; k<numParamsCF; k++) {
			SEest = SEest + 2 * Cov[ix(j,k)] * Der[j] * Der[k];
			}
		}
	SEest=w*Math.sqrt(SEest);
	var yco=yc; var ycl=yc-St95*SEest; var ych=yc+St95*SEest;



	}

ccSW = cccSW;
ccAv = cccAv / ccSW;
ccSD = cccSD / ccSW;
var GenR2 = (ccSD-(SSq/ccSW))/ccSD;
var GenR = Math.sqrt(GenR2);

RMS = Math.sqrt(SSq/Max(1,dgfr));
if (theLanguage == "en") theResultsCF += "<p>RMS Error = " + vFmt(RMS) + "; &nbsp; &nbsp; <span class = math>df</span> = " + dgfr + "; &nbsp; &nbsp; SSE = " + vFmt(SSq);
else theResultsCF += "<p>Error RMS = " + vFmt(RMS) + "; &nbsp; &nbsp; <span class = math>df</span> = " + dgfr + "; &nbsp; &nbsp; SSE = " + vFmt(SSq);

theResultsCF += "<br><span class = math>r</span> = " + vFmt(GenR) + "; &nbsp; &nbsp; <span class = math>r</span><sup>2</sup> = " + vFmt(GenR2);

AIC = CFnumPoints * Math.log( SSq / CFnumPoints ) + 2 * numParamsCF;
AICc = AIC + ( 2 * numParamsCF * (numParamsCF + 1) ) / ( CFnumPoints - numParamsCF - 1 );

for (i=0; i<numParamsCF; i++) { var s = Arr[ix(i,i)]; Arr[ix(i,i)] = 1;
	for (k=0; k<=numParamsCF; k++) { Arr[ix(i,k)] = Arr[ix(i,k)] / s; }
	for (j=0; j<numParamsCF; j++) {
		if (i!=j) { s = Arr[ix(j,i)]; Arr[ix(j,i)] = 0;
			for (k=0; k<=numParamsCF; k++) {
				Arr[ix(j,k)] = Arr[ix(j,k)] - s * Arr[ix(i,k)];
				}
			}
		}
	}
var FRelax = 1.0; // rate of conv
if (theLanguage == "en") theResultsCF += "<p><b>Parameter Estimates</b><table>";
else theResultsCF += "<p><b>Estimaci&oacute;n de par&aacute;metros</b><table>";
for( i=0; i<numParamsCF; i++) {

	locationCF[i] = locationCF[i] + FRelax * Arr[ix(i,numParamsCF)];
	if (isNaN(locationCF[i]))   {
//testDisplay("/n" + "Step4 FAILED" + "\n" + ix(i,numParamsCF) + " arr val = " + Arr[ix(i,numParamsCF)]);
		abandonShip ();
		i = numParamsCF;
		}
	
	
	SEP[i] = RMS * Math.sqrt(Arr[ix(i,i)]);
	theResultsCF += "<tr><td>" + theVariablesCF[i] +" = "+  locationCF[i].toPrecision(displayAccuracyCF) +" &plusmn; " + SEP[i].toPrecision(displayAccuracyCF) + ";  &nbsp; &nbsp; 	<span class = math>p</span> = "+ (StudT(locationCF[i]/SEP[i] ,dgfr)).toPrecision(displayAccuracyCF) + "</td></tr>";
	}

theResultsCF += "</table></div>"

// Now clear it out for later



for (var i = 0; i <= numParamsCF-1; i++) {
	stErCF[i] = SEP[i];
	pValCF = vFmt(StudT(locationCF[i]/SEP[i],dgfr));
	} // i


}

function isGood(theLocation) {
for (var i = 0; i <= numParamsCF-1; i++) if (isNaN(theLocation[i])) return (false);
return(true);

} // isGood

// *** 
function updateSSE(val) {
// if no value there it updates and displays SSE
// otherwise just updates
if (!(val > 0)) {
	prevSSE = minSSE;
	minSSE = SSECF(locationCF);
	}
if (isNaN(minSSE)) abandonShip();
document.theFormP.sse.value = minSSE.toPrecision(acc2CF);
}

function showParams(inArr) {
document.theFormP.params.value = replaceSubstring(roundCF(inArr).toString(), ",", ", ");
}