package iterators;

import java.io.*;
import main.BetterTokenizer;
import main.GeneralInput;
import java.lang.*;
import main.Network;
import main.Model;
import main.Cell;
import main.Node;
import affectors.Affector;
import parameters.ParameterSet;
import stoppers.NoChangeStop;
import stoppers.SimpleStop;
import main.MoreMath;
import main.Globals;

public class ACSCIterator extends ModelIterator implements Runnable {
	float		toleranceRel, toleranceAbs;
	float		epsilon = 0.0001f;
	int		[]	nodes = new int[4];
	int			numNodes;
	int			modelNumNodes;
	int		[]	numIncrs = new int[4];
	float	[]	startVals = new float[4];
	float	[]	incrFactors = new float[4];
	float	[]	curVals = new float[4];
	int		[]	startIdentical = new int[4];
	int			numIdentical;
	int		[]	numTimes = new int[20];
	int		[]	ststDistances = new int[30];
	int			outputType = 0;		// 0 = standard output, 1 = save st st's output

// Put this here so don't have to keep remaking it for each use
	private SimpleStop	afterStop = new SimpleStop(1000),
				afterStop2 = new SimpleStop(1000),
				hysteriaStop = new SimpleStop(300);

	private	int		num0Points;

// Put these here so I don't have to pull contortions to pass around pointers to numPts
	private float	[]	ststPts;
	private int	[]	numPtsFound = new int[20];
	private int		numPts = 0;

	int			numHysterias = 0, avg1protdown_no = 0, avg2protdown_no = 0,
				avg1rnadown_no = 0, avg1rnaup_no = 0, avg2rnadown_no = 0, avg2rnaup_no = 0;
	float		avg1rnaup = 0, avg1rnaupsqr = 0, avg1rnadown = 0, avg1rnadownsqr = 0,
				avg2rnaup = 0, avg2rnaupsqr = 0, avg2rnadown = 0, avg2rnadownsqr = 0,
				avg1protdown = 0, avg1protdownsqr = 0, avg2protdown = 0, avg2protdownsqr = 0;
	int			num2protdown=0, num1protdown=0, num2rnadown=0, num1rnadown=0, 
				num2rnaup=0, num1rnaup=0;

	public ACSCIterator() { setPrint(true); }	// Have a no argument initializer so we can do new_by_name

	public void init(Network network, Model model)  {				
		super.init(network, model);
		
		//modelNumNodes = 4; //model.net.numNodes;
		modelNumNodes = 1;
		int		i;
		for(i = 0; i < 20; i++) numTimes[i] = 0;
		for(i = 0; i < 30; i++) ststDistances[i] = 0;
	}
	
	public ModelIterator copy()  throws Exception {
		ACSCIterator newIter = new ACSCIterator();
		newIter.init(this.network, this.model);
		newIter.nParsTV = this.nParsTV;
		newIter.parsTV = this.parsTV.copy();
		newIter.theFunction = this.theFunction.copy();
		newIter.toleranceAbs = this.toleranceAbs;
		newIter.toleranceRel = this.toleranceRel;
		newIter.nodes = new int[nodes.length];
		System.arraycopy(nodes, 0, newIter.nodes, 0, nodes.length);
		newIter.startVals = new float[startVals.length];
		System.arraycopy(startVals, 0, newIter.startVals, 0, startVals.length);
		newIter.incrFactors = new float[incrFactors.length];
		System.arraycopy(incrFactors, 0, newIter.incrFactors, 0, incrFactors.length);
		newIter.numIncrs = new int[numIncrs.length];
		System.arraycopy(numIncrs, 0, newIter.numIncrs, 0, numIncrs.length);
		newIter.curVals = new float[curVals.length];
		System.arraycopy(curVals, 0, newIter.curVals, 0, curVals.length);
		newIter.startIdentical = new int[startIdentical.length];
		System.arraycopy(startIdentical, 0, newIter.startIdentical, 0, startIdentical.length);

		return newIter;
	}

	public void saveOutputTags(PrintWriter ps, String inset) {
		super.saveOutputTags(ps, inset);
		ps.println(inset + "&Proportion0\tnumber");
	}

	public void saveOutput(PrintWriter ps, String inset) {
		super.saveOutput(ps, inset);
		int	num = Cell.arrayWidth * Cell.arrayHeight;
		for(int n = 2; n < numNodes; n++) num *= numIncrs[n];
		ps.println(inset + "&Proportion0\t" + num0Points / (float)num);
	}

	public void loadParameters(BetterTokenizer tokenizer) throws Exception {
		super.loadParameters(tokenizer);	// The grunt work is done by parent class

		// This is initialization stuff we need to do after figuring out how many
		// parameters will be varied.		
	}
	
	
	// Put any parameters specific to this iterator class in here, as if clauses.
	protected void loadParameter(String info, BetterTokenizer tokenizer) throws Exception {
		if(info.equals("Node")) {
			if(numNodes == nodes.length) {	// If too many nodes, expand arrays
				int [] temp_nodes = new int[numNodes + 4];
				System.arraycopy(nodes, 0, temp_nodes, 0, numNodes);
				nodes = temp_nodes;
				float [] temp_starts = new float[numNodes + 4];
				System.arraycopy(startVals, 0, temp_starts, 0, numNodes);
				startVals = temp_starts;
				float [] temp_incrs = new float[numNodes + 4];
				System.arraycopy(incrFactors, 0, temp_incrs, 0, numNodes);
				incrFactors = temp_incrs;
				int [] temp_numincrs = new int[numNodes + 4];
				System.arraycopy(numIncrs, 0, temp_numincrs, 0, numNodes);
				numIncrs = temp_numincrs;
				curVals = new float[numNodes + 4];
			}
			GeneralInput.nextToken(tokenizer); 
			nodes[numNodes] = cells[0].getNodeNum(tokenizer.sval);
			GeneralInput.nextToken(tokenizer);
			startVals[numNodes] = (float)tokenizer.nval;
			GeneralInput.nextToken(tokenizer);
			incrFactors[numNodes] = (float)tokenizer.nval;
			GeneralInput.nextToken(tokenizer);
			numIncrs[numNodes] = (int)tokenizer.nval;
			numNodes++;
			
			modelNumNodes = numNodes + numIdentical;	// Assume that the number of important nodes for finding st st, etc.
										// is just the nodes that user input
			ststPts = new float[20 * modelNumNodes];
		}
		else if(info.equals("StartIdentical")) {
			if(numIdentical * 2 == startIdentical.length) {
				int [] temp_ident = new int[numIdentical * 2 + 8];
				System.arraycopy(startIdentical, 0, temp_ident, 0, numIdentical * 2);
				startIdentical = temp_ident;
			}
			GeneralInput.nextToken(tokenizer);
			startIdentical[numIdentical * 2] = cells[0].getNodeNum(tokenizer.sval);
			GeneralInput.nextToken(tokenizer);
			startIdentical[numIdentical * 2 + 1] = cells[0].getNodeNum(tokenizer.sval);
			numIdentical++;

			modelNumNodes = numNodes + numIdentical;	// Assume that the number of important nodes for finding st st, etc.
										// is just the nodes that user input
			ststPts = new float[20 * modelNumNodes];
		}
		else if(info.equals("RelativeTolerance")) {GeneralInput.nextToken(tokenizer); toleranceRel = (float)tokenizer.nval; }
		else if(info.equals("AbsoluteTolerance")) {GeneralInput.nextToken(tokenizer); toleranceAbs = (float)tokenizer.nval; }
		else if(info.equals("Epsilon")) {GeneralInput.nextToken(tokenizer); epsilon = (float)tokenizer.nval; }
		else if(info.equals("OutputType")) { 
			GeneralInput.nextToken(tokenizer);
			if(tokenizer.sval.equals("Standard")) outputType = 0;
			else outputType = 1;
		}
		else super.loadParameter(info, tokenizer);
	}

	public void doRun()  {
		int			i, j, x, y, n;
		Cell		cell;
		int		[]	num_incrs = new int[curVals.length];
		
		reset();
		numPts = 0;
		num0Points = 0;
		
		for(n = 0; n < numNodes; n++) {
			curVals[n] = startVals[n];
			num_incrs[n] = 0;
		}
		boolean done_cycling = false;
		while(!done_cycling) {
			// Set up all the cells. The width of the grid gets increasing values of nodes[0].
			// The height of the grid gets increasing values of nodes[1]. Other nodes are set as 
			// in curVals or else left at their default initial values. If a node is in the
			// startIdentical list then it is started the same as the node it refers to.
			curVals[1] = startVals[1];
			for(y = 0; y < Cell.arrayHeight; y++) {
				curVals[0] = startVals[0];
				for(x = 0; x < Cell.arrayWidth; x++) {
					cell = cells[y * Cell.arrayWidth + x];
					for(n = 0; n < numNodes; n++)
						cell.setInitialValue(nodes[n], curVals[n]);
					for(n = 0; n < numIdentical; n++)
						cell.setInitialValue(startIdentical[n*2], cell.getInitialValue(startIdentical[n*2+1]));						
					curVals[0] *= incrFactors[0];
				}
				curVals[1] *= incrFactors[1];
			}

		// The meat of it
			findStStForArray();
		
		// Move to the next array set
			n = 2;	// The first two nodes are always cycled completely in a run, along the x and y dimensions of the array of cells
			if(n < numNodes) {
				curVals[n] *= incrFactors[n];
				num_incrs[n] += 1;
			}
			else done_cycling = true;
			while(n < numNodes && num_incrs[n] >= numIncrs[n]) {
				curVals[n] = startVals[n];
				num_incrs[n] = 0;
				n++;
				if(n == numNodes) done_cycling = true;
				else {
					curVals[n] *= incrFactors[n];
					num_incrs[n] += 1;
				}
			}
		}

	// Figure out which st st's are really good ones
		System.out.println("stability");
		checkStStStability();
		
	// Save results to file
		if(outputType == 1) {
			print("Params");
			parsTV.getModel();
			for(i = 0; i < nParsTV; i++)
				print("\t" + parsTV.getValue(i));
			println("");
			print("Num st st pts = " + numPts + "  :");
			numTimes[ classifyStSt(ststPts, modelNumNodes, numPts) ]++;
			for(i = 0; i < 14; i++) print("  " + numTimes[i]);
			println("");
			for(i = 0; i < numPts; i++) {
				for(j = 0; j < modelNumNodes; j++) print(ststPts[i*modelNumNodes+j] + "\t");
				println("");
			}
			if(numPts == 2 && classifyStSt(ststPts, modelNumNodes, numPts) == 4) {
				// Check for histeresis
				if(ststPts[nodes[0]] > toleranceAbs)
					checkHisteresis(ststPts[nodes[0]], ststPts[nodes[1]]);
				else checkHisteresis(ststPts[modelNumNodes + nodes[0]], ststPts[modelNumNodes + nodes[1]]);
				print(num1rnaup + "\t" + avg1rnaup/num1rnaup + "\t" + MoreMath.getStdDev(avg1rnaup, avg1rnaupsqr, num1rnaup) + "\t" + avg1rnaup_no + "\t");
				print(num1rnadown + "\t" + avg1rnadown/num1rnadown + "\t" + MoreMath.getStdDev(avg1rnadown, avg1rnadownsqr, num1rnadown) + "\t" + avg1rnadown_no + "\t");
				print(num2rnaup + "\t" + avg2rnaup/num2rnaup + "\t" + MoreMath.getStdDev(avg2rnaup, avg2rnaupsqr, num2rnaup) + "\t" + avg2rnaup_no + "\t");
				print(num2rnadown + "\t" + avg2rnadown/num2rnadown + "\t" + MoreMath.getStdDev(avg2rnadown, avg2rnadownsqr, num2rnadown) + "\t" + avg2rnadown_no + "\t");
				print(num1protdown + "\t" + avg1protdown/num1protdown + "\t" + MoreMath.getStdDev(avg1protdown, avg1protdownsqr, num1protdown) + "\t" + avg1protdown_no + "\t");
				print(num2protdown + "\t" + avg2protdown/num2protdown + "\t" + MoreMath.getStdDev(avg2protdown, avg2protdownsqr, num2protdown) + "\t" + avg2protdown_no + "\t");
				println();
			}

		}
		finalScore = classifyStSt(ststPts, modelNumNodes, numPts);

		super.stopRun();
	
		System.out.println("leave");
	}


	private void findStStForArray() {
		int			i, j;
		float	[]	matrix = new float[modelNumNodes*modelNumNodes];
		float	[]	vec = new float[modelNumNodes];
		float	[]	values = new float[modelNumNodes];
		float	[]	deriv = new float[modelNumNodes];
		float		dist, step_size;
		Cell		cell;
		
		F(p);	// Don't care about value, I think, except maybe if its oscillating?
		// The stopper should have been simple stop, going for a short time to get near st st
		for(int c = 0; c < numCells; c++) {
			boolean done = false, bad = false;
			cell = cells[c];
			for(i = 0; i < modelNumNodes; i++)
				values[i] = cell.getValue(i);
			dist = 0;
			for(i = 0; i < modelNumNodes; i++) dist += sqr(values[i]);
			dist = (float)Math.sqrt(dist);
			step_size = 0.1f;
			double last_dist2 = 1000000;
			int	num_times = 0;
			while(!done && num_times < 100) {
				// Figure out the move we need to make
				MoreMath.calcJacobian(matrix, cell, modelNumNodes, epsilon);
				for(i = 0; i < modelNumNodes; i++)
					deriv[i] = -cell.nodes[i].getAffectorsValue(0);
				int	[]	index = new int[modelNumNodes];
				ludcmp(matrix, modelNumNodes, index);
				lubksb(matrix, modelNumNodes, index, deriv);
				
				double dist2 = 0;
				for(i = 0; i < modelNumNodes; i++) dist2 += sqr(deriv[i]);
				dist2 = Math.sqrt(dist2);
				if(dist2 > last_dist2) { done = true; bad = true; }
				if(dist2 < 0.000001) done = true;
				else {
					// Make sure we don't move too far
					if(dist2 > step_size * dist) {
						for(i = 0; i < modelNumNodes; i++)
							deriv[i] *= step_size * dist / dist2;
					}
					// Make the move
					for(i = 0; i < modelNumNodes; i++) {
						values[i] += deriv[i];
						cell.setValue(i, values[i]);
					}

					for(i = 0; i < modelNumNodes; i++) dist += sqr(values[i]);
					dist = (float)Math.sqrt(dist);
					if(dist < 0.01) dist = 0.01f;	// To prevent one stall mechanism
				}
				last_dist2 = dist2;
				num_times++;
			}	// while(!done)
			if(num_times > 99) { bad = true; 
					System.out.println("Tried > 99 times for st st " + last_dist2); }

			if(!bad) {
				int		found = -1;
				float	min_dist = 10000;
				i = 0;
				while(i < numPts) {
					if(isClose(values, 0, ststPts, i*modelNumNodes, modelNumNodes, toleranceRel, toleranceAbs))
					{
						dist = calcDistance(values, 0, ststPts, i*modelNumNodes, modelNumNodes);
						if(found < 0) { min_dist = dist; found = i; }
						else if(dist < min_dist) { min_dist = dist; found = i; }
					}
					i++;
				}
				if(found >= 0) {
					numPtsFound[found]++;
				}
				else if(numPts < 20) {
					for(i = 0; i < modelNumNodes; i++)
						ststPts[numPts*modelNumNodes+i] = values[i];
					numPts++;
				}
				if(classifyStSt(values, modelNumNodes, 1) == 1)  num0Points++;
			}	// if(! bad)
		}	// for c < numCells
	}

	private void checkStStStability() {
	// Check each of the steady state points to see if its stable
		boolean	good = true;
		float		dist, min_dist;
		int			min_pt, store_num_cells;
		int			i, j, k;
		float	[]	matrix = new float[modelNumNodes*modelNumNodes];
		float	[]	values = new float[modelNumNodes];
		float []	eigen_real = new float[modelNumNodes];
		float []	eigen_imag = new float[modelNumNodes];
		float [] 	stst = new float[modelNumNodes];
		boolean	[]	good_points = new boolean[numPts],
						good_points2 = new boolean[numPts];

// First flag ones with negative eigenvalues as being definitely good.
		for(i = 0; i < numPts; i++) good_points[i] = false;
		i = 0;
		while(i < numPts) {
			try {
				for(j = 0; j < modelNumNodes; j++)
					cells[0].setValue(j, ststPts[i*modelNumNodes+j]);
				MoreMath.calcJacobian(matrix, cells[0], modelNumNodes, epsilon);
				MoreMath.eigenvalues(matrix, modelNumNodes, eigen_real, eigen_imag);
				good_points[i] = true;
				for(j = 0; j < modelNumNodes; j++)
					if(eigen_real[j] > 0) good_points[i] = false;
				i++;
			} catch(Exception e) { System.out.println(e.toString()); i++; }
		}

	// Run the st st points which are left out for long time to see if they stick around or
	// go far away. Do this 4 times, so we start from different random deviations from each st st point (in case its a saddle)
		
		for(i = 0; i < numPts; i++) good_points2[i] = true;
		for(int a = 0; a < 4; a++) {
			k = 0;
			for(i = 0; i < numPts; i++) {
				if(!good_points[i] && good_points2[i]) {
					// Start cells a little distance away from each st st
					for(j = 0; j < modelNumNodes; j++)
						cells[k].setInitialValue(j, ststPts[i*modelNumNodes+j] + MoreMath.random() * 0.01f - 0.005f);
					k++;
				}
			}
	
			if(k > 0) {
				setStopper(afterStop);
				store_num_cells = model.numCells;
				model.numCells = k;
				F(p);	// Just to run for a while
				model.numCells = store_num_cells;
				resetStopper();
				
				// Now go through and check each cell's nodes to see if they are close to any of 
				// the other steady states
				int num_removed = 0;
				i = 0;
				k = 0;
				for(i = 0; i < numPts; i++) {
					if(!good_points[i] && good_points2[i]) {
						min_dist = 100000;
						min_pt = -1;
						for(j = 0; j < modelNumNodes; j++) values[j] = cells[k].getValue(j);
						for(j = 0; j < numPts; j++) {
							dist = calcDistance(values, 0, ststPts, j*modelNumNodes, modelNumNodes);
							if(dist < min_dist) { min_dist = dist; min_pt = j; }
						}
						if(min_dist < 0.01 && min_pt != i) {
							// It fell into another st st
							good_points2[i] = false;
						}
						else if(min_dist > toleranceAbs) good_points2[i] = false;
						k++;
					}
				}
			}
		}

		// Get rid of the ones that failed both tests
		for(i = numPts - 1; i >= 0; i--) {	// Go backwards so we don't have to keep track of array positions
			if(!good_points[i] && !good_points2[i]) {
				removePtFromArray(ststPts, numPts, modelNumNodes, i);
				numPts--;
				System.out.println("Removed point that was not stable");
			}
		}

// For each steady state that is within 0.2 of another st st, put a point inbetween
// the nearest neighbors and see if that point goes definitively towards one or
// another of the st st's. If not, combine them together into one. Repeat as neccesary.
/*		float []	values2 = new float[modelNumNodes];
		i = 0;
		while(i < numPts) {
			// First find the nearest st st to this one
			min_dist = 1000000;
			min_pt = -1;
			for(j = 0; j < numPts; j++) {
				if(j != i) {
					dist = calcDistance(ststPts, i*modelNumNodes, ststPts, j*modelNumNodes, modelNumNodes);
					if(dist < min_dist) {
						min_dist = dist;
						min_pt = j;
					}
				}
			}
			if(min_dist < 0.2) {
				// Put a point in the middle
				for(j = 0; j < modelNumNodes; j++) {
					values2[j] = (ststPts[i*modelNumNodes+j] + ststPts[min_pt*modelNumNodes+j]) / 2;
					cells[0].setInitialValue(j, values2[j]);
				}
				// Run this for 1000 timesteps
				setStopper(afterStop2);
				store_num_cells = model.numCells;
				model.numCells = 1;
				F(p);	// Just to run for a while
				model.numCells = store_num_cells;
				resetStopper();
				
				// Check whether cell's value went significantly towards one or the other st st
				for(j = 0; j < modelNumNodes; j++)
					values[j] = cells[0].getValue(j);
				if(calcDistance(ststPts, i*modelNumNodes, values, 0, modelNumNodes) > 0.2 * min_dist &&
						calcDistance(ststPts, min_pt*modelNumNodes, values, 0, modelNumNodes) > 0.2 * min_dist &&
						calcDistance(values, 0, values2, 0, modelNumNodes) < 0.05) {
					// Combine the points together
					removePtFromArray(ststPts, numPts, modelNumNodes, min_pt);
					numPts--;
					if(min_pt < i) i--;
				}
				else i++;
			}
			else i++;
		}
*/	}
	
	public boolean didPass() { return true; }
	

	private boolean isClose(float [] values1, int pt1, float [] values2, int pt2, int numpts, float toleranceRel, float toleranceAbs) {
		int i = 0;
		boolean close = true;
		while(i < numpts && close) {
			close &= (Math.abs((values1[pt1+i] - values2[pt2+i]) / values2[pt2+i]) < toleranceRel) 
						|| (Math.abs(values1[pt1+i] - values2[pt2+i]) < toleranceAbs);
			i++;
		}
		return close;
	}

	private final float sqr(float val) { return val * val; }

	private float calcDistance(float [] vec1, int pos1, float [] vec2, int pos2, int n) {
		float dist = 0;
		for(int i = 0; i < n; i++)
			dist += sqr(vec1[pos1+i] - vec2[pos2+i]);
		return((float)Math.sqrt(dist));
	}

	private void removePtFromArray(float [] pts, int num_pts, int num_axes, int remove_pt) {
		int		i, j;
		for(i = remove_pt; i < num_pts - 1; i++) {
			for(j = 0; j < num_axes; j++)
				pts[i*modelNumNodes+j] = pts[(i+1)*modelNumNodes+j];
		}
	}

	private void gaussj(float [] a, int n, float [] b, int m) {
		int	[] indxc = new int[n];
		int	[] indxr = new int[n];
		int	[] ipiv = new int[n];
		int	i, icol = 0, irow = 0, j, k, l, ll;
		float	big, dum, pivinv, temp;
		
		for(j = 0; j < n; j++) ipiv[j] = 0;
		for(i = 0; i < n; i++) {
			big = 0.0f;
			for(j = 0; j < n; j++) {
				if(ipiv[j] != 1) {
					for(k = 0; k < n; k++) {
						if(ipiv[k] == 0) {
							if(Math.abs(a[j * n + k]) >= big) {
								big = Math.abs(a[j * n + k]);
								irow = j;
								icol = k;
							}
						}
						else if(ipiv[k] > 1) System.out.println("Error: singular matrix");
					}
				}
			}
			ipiv[icol]++;
			
			if(irow != icol) {
				for(l = 0; l < n; l++) {
					temp = a[irow * n + l];
					a[irow * n + l] = a[icol * n + l];
					a[icol * n + l] = temp;
				}
				for(l = 0; l < m; l++) {
					temp = b[irow * m + l];
					b[irow * m + l] = b[icol * m + l];
					b[icol * m + l] = temp;
				}
			}
			indxr[i] = irow;
			indxc[i] = icol;
			if(a[icol * n + icol] == 0) System.out.println("Error: singular matrix");
			pivinv = 1.0f / a[icol * n + icol];
			a[icol * n + icol] = 1.0f;
			for(l = 0; l < n; l++)
				a[icol * n + l] *= pivinv;
			for(l = 0; l < m; l++)
				b[icol * m + l] *= pivinv;

			for(ll = 0; ll < n; ll++) {
				if(ll != icol) {
					dum = a[ll * n + icol];
					a[ll * n + icol] = 0.0f;
					for(l = 0; l < n; l++)
						a[ll * n + l] -= a[icol * n + l] * dum;
					for(l = 0; l < m; l++)
						b[ll * m + l] -= b[icol * m + l] * dum;
				}
			}
		}
		for(l = n - 1; l >= 0; l--) {
			if(indxr[l] != indxc[l]) {
				for(k = 0; k < n; k++) {
					temp = a[k * n + indxr[l]];
					a[k * n + indxr[l]] = a[k * n + indxc[l]];
					a[k * n + indxc[l]] = temp;
				}
			}
		}		
	}
	

	private void ludcmp(float [] a, int n, int [] indx) {
		int		i, imax = 0, j, k;
		float	big, dum, sum, temp;
		float	[] vv = new float[n];
		float	d;
		final float TINY = 1.0e-20f;
		
		d = 1.0f;
		for(i = 0; i < n; i++) {
			big = 0.0f;
			for(j = 1; j < n; j++)
				if((temp = Math.abs(a[i*n+j])) > big) big = temp;
			if(big == 0.0f) { System.out.println("Error: Singluar matrix"); return; }
			vv[i] = 1.0f/big;
		}
		for(j = 0; j < n; j++) {
			for(i = 0; i < n; i++) {
				sum=a[i*n+j];
				for(k = 1; k < i; k++) sum -= a[i*n+k] * a[k*n+j];
				a[i*n+j] = sum;
			}
			big = 0.0f;
			for(i = j; i < n; i++) {
				sum = a[i*n+j];
				for(k = 0; k < j; k++)
					sum -= a[i*n+k]*a[k*n+j];
				a[i*n+j] = sum;
				if((dum = vv[i] * Math.abs(sum)) >= big) {
					big = dum;
					imax = i;
				}
			}
			if(j != imax) {
				for(k = 0; k < n; k++) {
					dum = a[imax*n+k];
					a[imax*n+k] = a[j*n+k];
					a[j*n+k] = dum;
				}
				d = -d;
				vv[imax] = vv[j];
			}
			indx[j] = imax;
			if(a[j*n+j] == 0.0f) a[j*n+j] = TINY;
			if(j != n) {
				dum = 1.0f/a[j*n+j];
				for(i = j+1; i < n; i++) a[i*n+j] *= dum;
			}
		}		
	}


	private void lubksb(float [] a, int n, int [] indx, float [] b) {
		int		i, ii=0, ip, j;
		float	sum;
		
		for(i = 0; i < n; i++) {
			ip = indx[i];
			sum = b[ip];
			b[ip] = b[i];
			if(ii != 0)
				for(j = ii; j < i-1; j++) sum -= a[i*n+j]*b[j];
			else if(sum != 0) ii = i;
			b[i] = sum;
		}
		for(i = n-1; i >= 0; i--) {
			sum = b[i];
			for(j=i+1;j < n; j++) sum -= a[i*n+j]*b[j];
			b[i]=sum/a[i*n+i];
		}
	}
	
/* Notes - this classification business is kind of tough. This was a rough cut. In my first
	analysis, I'm not combining together some of the following scores. 11 and 16 both say
	that the two top ac states are equal while one of those has a lower sc concentration than
	the other, so I'm sticking those together. Similarly, 9 and 14 both say that ac is 0, 0.5, 1
	and sc is 0, 0.5, 1 (in that combination - ie. when ac is 0.5, sc is 0.5). 14 just has total ac 
	concentration less than sc concentration, which now doesn't seem relevant to anything.
	12 must go with 11 + 16, except with sc presentation high in both states rather than ac (by a
	process of elimination). It's odd, though that # of 12 != # 11 + # 16... 
	For presentation (but not for x-correlation)
	I'm combining 8 and 13, one of which is ac 0, 0, 1, sc 0, 1, 1, while the other is vis-versa.
	10 and 15 are, not surprisingly, absent from the solution set. */
	// This one is pretty specific to my Ac/Sc model, probably not useful elsewhere.
	// It assumes 4 nodes, with protein/rna right next to each other
	// Will also now work for 2 nodes
	// Assumes there will never be > 4 st st points.
	// Also assumes Nodes come in pairs, rna and protein
	private int classifyStSt(float [] stst_pts, int num_nodes, int num_pts) {
		int		zero_state = -1, zero_high_state = -1;
		int		i;
		
		if(numNodes< 4) return classifyStSt2(stst_pts, num_nodes, num_pts);
		switch(num_pts) {
			case 0:
				return 0;
			case 1:
				if(stst_pts[0] < toleranceAbs && stst_pts[3] < toleranceAbs) return 1;
				if(stst_pts[0] > toleranceAbs && stst_pts[3] > toleranceAbs) return 2;
				return 3;
			case 2:
				i = 0;
				while(i < num_pts && zero_state == -1) {
					if(stst_pts[i*num_nodes + 0] < toleranceAbs && stst_pts[i*num_nodes + 3] < toleranceAbs)
						zero_state = i;
					i++;
				}
				i = 0;
				while(i < num_pts && zero_high_state == -1) {
					if((stst_pts[i*num_nodes + 0] < toleranceAbs && stst_pts[i*num_nodes + 3] > toleranceAbs)
						|| (stst_pts[i*num_nodes + 0] > toleranceAbs && stst_pts[i*num_nodes + 3] < toleranceAbs))
						zero_high_state = i;
					i++;
				}
				if(zero_state >= 0 && zero_high_state < 0) return 4;
				if(zero_state >= 0 && zero_high_state >= 0) return 5;
				return 6;
			case 3:
				i = 0;
				while(i < num_pts && zero_state == -1) {
					if(stst_pts[i*num_nodes + 0] < toleranceAbs && stst_pts[i*num_nodes + 3] < toleranceAbs)
						zero_state = i;
					i++;
				}
				if(zero_state < 0) return 7;
				if(zero_state != 0) {
					for(i = 0; i < num_nodes; i++) {
						float temp = stst_pts[zero_state*num_nodes + i];
						stst_pts[zero_state*num_nodes + i] = stst_pts[i];
						stst_pts[i] = temp;
					}
				}
			// 8-11 are SC high, 13 - 16 are AC high (or maybe vis-versa)
				float sum0 = stst_pts[0*num_nodes+0] + stst_pts[1*num_nodes+0] + stst_pts[2*num_nodes+0];
				float sum3 = stst_pts[0*num_nodes+3] + stst_pts[1*num_nodes+3] + stst_pts[2*num_nodes+3];
				if(stst_pts[1*num_nodes+0] < toleranceAbs || stst_pts[2*num_nodes+0] < toleranceAbs)
					return 8;
				if(stst_pts[1*num_nodes+3] < toleranceAbs || stst_pts[2*num_nodes+3] < toleranceAbs)
					return 13;
				if(stst_pts[1*num_nodes+0] < stst_pts[2*num_nodes+0] * (1.0f - toleranceRel)) {
					if(stst_pts[1*num_nodes+3] < stst_pts[2*num_nodes+3] * (1.0f - toleranceRel)) {
						if(sum0 < sum3) return 9;
						else return 14;
					}
					if(stst_pts[1*num_nodes+3] > stst_pts[2*num_nodes+3] * (1.0f + toleranceRel)) {
						if(sum0 < sum3) return 10;
						else return 15;
					}
					if(sum0 < sum3) return 11;
					else return 16;
				}
				if(stst_pts[1*num_nodes+0] > stst_pts[2*num_nodes+0] * (1.0f + toleranceRel)) {
					if(stst_pts[1*num_nodes+3] > stst_pts[2*num_nodes+3] * (1.0f + toleranceRel)) {
						if(sum0 < sum3) return 9;
						else return 14;
					}
					if(stst_pts[1*num_nodes+3] < stst_pts[2*num_nodes+3] * (1.0f - toleranceRel)) {
						if(sum0 < sum3) return 10;
						else return 15;
					}
					if(sum0 < sum3) return 11;
					else return 16;
				}
				return 12;
			case 4:
				return 17;
		}
		return 18;
	}
	
	
	/** This version is for classifying steady states when there is only one gene
		being considered. It uses the same numbering scheme as above, but some numbers
		would be missing. */
	private int classifyStSt2(float [] stst_pts, int num_nodes, int num_pts) {
		int		zero_state = -1, zero_high_state = -1;
		int		i;
		
		switch(num_pts) {
			case 0:
				return 0;
			case 1:
				if(stst_pts[0] < toleranceAbs) return 1;
				if(stst_pts[0] > toleranceAbs) return 2;
			case 2:
				i = 0;
				while(i < num_pts && zero_state == -1) {
					if(stst_pts[i*num_nodes + 0] < toleranceAbs)
						zero_state = i;
					i++;
				}
				if(zero_state >= 0) return 4;
				return 6;
			case 3:
				i = 0;
				while(i < num_pts && zero_state == -1) {
					if(stst_pts[i*num_nodes + 0] < toleranceAbs)
						zero_state = i;
					i++;
				}
				if(zero_state < 0) return 7;
				if(zero_state != 0) {
					for(i = 0; i < num_nodes; i++) {
						float temp = stst_pts[zero_state*num_nodes + i];
						stst_pts[zero_state*num_nodes + i] = stst_pts[i];
						stst_pts[i] = temp;
					}
				}
				return 8;
			case 4:
				return 17;
		}
		return 18;
	}
	
	private void checkHisteresis(float ststnode0, float ststnode1) {
		int		n, param;
		float	x, y, addition;
		
		System.out.println("Checking hysteria\r\r");
		Globals.flag = 1;

		setStopper(hysteriaStop);
		int store_num_cells = model.numCells;
		model.numCells = 1;


		addition = findSwitchPointUp("alpha_WG", "", ststnode0, ststnode1);
		num1rnaup++;
		if(addition < 5) {
			avg1rnaup += addition;
			avg1rnaupsqr += addition * addition;
		}
		else if(addition < 5000) avg1rnaup_no++;
		else num1rnaup--;

		addition = findSwitchPointUp("alpha_WG", "beta_WG", ststnode0, ststnode1);
		num2rnaup++;
		if(addition < 5) {
			avg2rnaup += addition;
			avg2rnaupsqr += addition * addition;
		}
		else if(addition < 5000) avg2rnaup_no++;
		else num2rnaup--;
		
		addition = findSwitchPointDown("alpha_WG", "", ststnode0, ststnode1);
		num1rnadown++;
		if(addition < 5) {
			avg1rnadown += addition;
			avg1rnadownsqr += addition * addition;
		}
		else if(addition < 5000) avg1rnadown_no++;
		else num1rnadown--;
		
		addition = findSwitchPointDown("alpha_WG", "beta_WG", ststnode0, ststnode1);
		num2rnadown++;
		if(addition < 5) {
			avg2rnadown += addition;
			avg2rnadownsqr += addition * addition;
		}
		else if(addition < 5000) avg2rnadown_no++;
		else num2rnadown--;
		
// Now do proteins
		addition = findSwitchPointDown("alpha_EMC", "", ststnode0, ststnode1);
		num1protdown++;
		if(addition < 5) {
			avg1protdown += addition;
			avg1protdownsqr += addition * addition;
		}
		else if(addition < 5000) avg1protdown_no++;
		else num1protdown--;
		
		addition = findSwitchPointDown("alpha_EMC", "beta_EMC", ststnode0, ststnode1);
		num2protdown++;
		if(addition < 5) {
			avg2protdown += addition;
			avg2protdownsqr += addition * addition;
		}
		else if(addition < 5000) avg2protdown_no++;
		else num2protdown--;

		numHysterias++;

		model.numCells = store_num_cells;
		resetStopper();
		Globals.flag = 0;

/*		try {
			param1 = parsTV.getPosition("alpha_WG");
			param2 = parsTV.getPosition("beta_WG");
		} catch(Exception e) { return; }
		for(n = 0; n < numNodes; n++) cells[0].setInitialValue(nodes[n], 0.0001f);
		for(n = 0; n < numIdentical; n++) cells[0].setInitialValue(startIdentical[n*2], 0.0001f);
		boolean done = false;
		p[param2] = 0;
		for(x = 0.0f, y = 0.001f; x < 2.0f && !done; x += y, y *= 1.5f) {
			p[param1] = x;
			F(p);	// Just to run for a while
//			println(x + "\t" + cells[0].getValue(nodes[0]) + "\t" + cells[0].getValue(nodes[1]));
			if(cells[0].getValue(nodes[1]) > ststnode1 - 0.02)
				done = true;
		}
		if(cells[0].getValue(nodes[1]) > ststnode1 - 0.02) {
			avg1rnaup += x / y;
			avg1rnaupsqr += (x/y) * (x/y);
		}
		else avg1rnaup_no++; 
		boolean done = false;
		for(x = 0.0f, y = 0.001f; x < 2.0f && !done; x += y, y *= 1.5f) {
			p[param1] = x;
			p[param2] = x;
			F(p);	// Just to run for a while
			if(cells[0].getValue(nodes[0]) > ststnode0 - 0.02 && cells[0].getValue(nodes[1]) > ststnode1 - 0.02)
				done = true;
		}
		if(cells[0].getValue(nodes[0]) > ststnode0 - 0.02 && cells[0].getValue(nodes[1]) > ststnode1 - 0.02) {
			avg2rnaup += x / y;
			avg2rnaupsqr += (x/y) * (x/y);
		}
		else avg2rnaup_no++; 
		p[param1] = 0.00001f;
		p[param2] = 0.00001f;

// Now try going down
		try {
			param1 = parsTV.getPosition("alpha_EMC");
			param2 = parsTV.getPosition("beta_EMC");
		} catch(Exception e) { return; }
		cells[0].setInitialValue(nodes[0], ststnode0);
		cells[0].setInitialValue(nodes[1], ststnode1);
		for(n = 0; n < numIdentical; n++) cells[0].setInitialValue(startIdentical[n*2], cells[0].getInitialValue(startIdentical[n*2+1]));
		done = false;
		p[param2] = 0.00001f;
		for(x = 0.0f; x > -2.0f && !done; x -= 0.005f) {
			p[param] = x;
			F(p);	// Just to run for a while
			println(x + "\t" + cells[0].getValue(nodes[0]) + "\t" + cells[0].getValue(nodes[1]));
			if(cells[0].getValue(nodes[0]) < -0.1 || cells[0].getValue(nodes[0]) > 1.5)
				done = true;
		}
		p[param] = 0.00001f;

// Now try proteins
		try {
			param = parsTV.getPosition("alpha_EMC");
		} catch(Exception e) { return; }
		cells[0].setInitialValue(nodes[0], ststnode0);
		cells[0].setInitialValue(nodes[1], ststnode1);
		for(n = 0; n < numIdentical; n++) cells[0].setInitialValue(startIdentical[n*2], cells[0].getInitialValue(startIdentical[n*2+1]));
		done = false;
		for(x = 0.0f; x > -2.0f && !done; x -= 0.005f) {
			p[param] = x;
			F(p);	// Just to run for a while
			println(x + "\t" + cells[0].getValue(nodes[0]) + "\t" + cells[0].getValue(nodes[1]));
			if(cells[0].getValue(nodes[0]) < -0.1 || cells[0].getValue(nodes[0]) > 1.5)
				done = true;
		}
		p[param] = 0.00001f;
*/		
	}

	private float findSwitchPointUp(String param_str1, String param_str2, float ststnode0, float ststnode1) {
		float 	x, y;
		int		n, param1, param2;
		final float multiplier = 1.5f;

		try {
			param1 = parsTV.getPosition(param_str1);
			if(param_str2.length() > 0)
				param2 = parsTV.getPosition(param_str2);
			else param2 = -1;
		} catch(Exception e) { return 10000; }
		for(n = 0; n < numNodes; n++) cells[0].setInitialValue(nodes[n], 0.0001f);
		for(n = 0; n < numIdentical; n++) cells[0].setInitialValue(startIdentical[n*2], 0.0001f);
		boolean done = false, bad = false;
		for(x = 0.0f, y = 0.001f; x < 2.0f && !done; x += y, y *= multiplier) {
			p[param1] = x;
			if(param2 >= 0) p[param2] = x;
			F(p);	// Just to run for a while
			if(cells[0].getValue(nodes[0]) > ststnode0 - 0.02 && cells[0].getValue(nodes[1]) > ststnode1 - 0.02)
				done = true;
			else if(cells[0].getValue(nodes[0]) > 1.01 || cells[0].getValue(nodes[1]) > 1.01) {
				done = true; bad = true;
			}
		}

		p[param1] = 0.00001f;
		if(param2 >= 0) p[param2] = 0.00001f;
		
		float val0 = cells[0].getValue(nodes[0]);
		float val1 = cells[0].getValue(nodes[1]);
		if(!bad && val0 > ststnode0 - 0.02 && val1 > ststnode1 - 0.02) {
			return x - (y/multiplier);
		}
		else if(bad && val0 > ststnode0 - 0.02 && val0 < 1.1 && val1 > ststnode1 && val1 < 1.1)
			return x - (y/multiplier);
		else if(bad && (val0 < ststnode0 - 0.05 || val1 < ststnode1 - 0.05))
			return 1000;
		else return 10000;
	}


	private float findSwitchPointDown(String param_str1, String param_str2, float ststnode0, float ststnode1) {
		float	x, y;
		final float multiplier = 1.5f;
		int		n, param1, param2;
		
		try {
			param1 = parsTV.getPosition(param_str1);
			if(param_str2.length() > 0)
				param2 = parsTV.getPosition(param_str2);
			else param2 = -1;
		} catch(Exception e) { return 10000; }
		cells[0].setInitialValue(nodes[0], ststnode0);
		cells[0].setInitialValue(nodes[1], ststnode1);
		for(n = 0; n < numIdentical; n++) cells[0].setInitialValue(startIdentical[n*2], 0.0001f);
		boolean done = false, bad = false;
		for(x = 0.0f, y = -0.001f; x > -2.0f && !done; x += y, y *= multiplier) {
			p[param1] = x;
			if(param2 >= 0) p[param2] = x;
			F(p);	// Just to run for a while
			if(cells[0].getValue(nodes[0]) < toleranceAbs && cells[0].getValue(nodes[1]) < toleranceAbs)
				done = true;
			else if(cells[0].getValue(nodes[0]) < -0.001 || cells[0].getValue(nodes[0]) < -0.001) {
				done = true; bad = true;
			}
		}
		
		p[param1] = 0.00001f;
		if(param2 >= 0) p[param2] = 0.00001f;
		
		float val0 = cells[0].getValue(nodes[0]);
		float val1 = cells[0].getValue(nodes[1]);
		if(!bad && val0 < toleranceAbs && val1 < toleranceAbs) {
			return x - (y/multiplier);
		}
		else if(bad && val0 < toleranceAbs && val0 > -0.1 && val1 < toleranceAbs && val1 > -0.1)
			return x - (y/multiplier);
		else if(bad && (val0 > toleranceAbs || val1 > toleranceAbs))
			return 1000;
		else return 10000;
	}
	
	private float findNoiseThreshold() {
		return 0f;
		
		
	}
}