package iterators;

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

public class ACSC_SwitchIterator extends ModelIterator implements Runnable {
	final static int	REPRESSORS = 0, ACTIVATORS = 1;
	float		toleranceAbs = 0.001f;
	float		epsilon = 0.0001f;
	float		externalNodeValue = 0.5f;
	int			repressionTime = 240;
	int			numCycles = 1000;
	int			mode = REPRESSORS;
	int		[]	nodes = new int[4];
	int		[]	repressors = new int[4];
	int		[]	numSwitches = new int[4];
	int			numNodes = 0, numRepressors = 0;
	ParameterSet	params = new ParameterSet();
	private MultiNodeThresholdStop	switchStop = new MultiNodeThresholdStop(300);
	private SimpleStop	simpleStop = new SimpleStop(300);
	private SimpleStop	ststStop = new SimpleStop(100);

	public ACSC_SwitchIterator() { setPrint(true); }	// Have a no argument initializer so we can do new_by_name
	
	public ModelIterator copy()  throws Exception {
		ACSC_SwitchIterator newIter = new ACSC_SwitchIterator();
		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.nodes = new int[nodes.length];
		System.arraycopy(nodes, 0, newIter.nodes, 0, nodes.length);
		newIter.repressors = new int[repressors.length];
		System.arraycopy(repressors, 0, newIter.repressors, 0, repressors.length);

		newIter.params = params.copy();

		return newIter;
	}

	public void saveOutputTags(PrintWriter ps, String inset) {
		super.saveOutputTags(ps, inset);
		for(int i = 0; i < numRepressors; i++) {
			ps.println(inset + "&Prop" + network.getNode(repressors[i]).getName() + "Switched\tnumber");
		}
	}

	public void saveOutput(PrintWriter ps, String inset) {
		super.saveOutput(ps, inset);
		for(int i = 0; i < numRepressors; i++) {
			ps.println(inset + "&Prop" + network.getNode(repressors[i]).getName() + "Switched\t" + numSwitches[i] / (float)numCycles);
		}
	}
	
	public void loadParameters(BetterTokenizer tokenizer) throws Exception {
		super.loadParameters(tokenizer);
		
		switchStop.init(model);
		if(mode == ACTIVATORS)
			switchStop.setMode(MultiNodeThresholdStop.HIGH);
		else switchStop.setMode(MultiNodeThresholdStop.LOW);
		switchStop.setEpsilon(0.01f);	// Note - I could set this to toleranceAbs, but
										// I choose to set it to a small value so it will be
										// rare that its not in attractive basin but still
										// within the tolerance. Since the activator/repressor
										// will get taken away before the end of the run, 
										// the system should escape back towards the initial state
										// before we check against toleranceAbs at the end of the run
		switchStop.setInterval(1f);
	}
	
	// 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;
			}
			GeneralInput.nextToken(tokenizer); 
			nodes[numNodes] = cells[0].getNodeNum(tokenizer.sval);
			numNodes++;
		}
		else if(info.equals("Repressor")) {
			if(numRepressors == repressors.length) {
				int [] temp_repress = new int[numRepressors + 4];
				System.arraycopy(repressors, 0, temp_repress, 0, numRepressors);
				repressors = temp_repress;
				numSwitches = new int[numRepressors + 4];
			}
			GeneralInput.nextToken(tokenizer);
			repressors[numRepressors] = cells[0].getNodeNum(tokenizer.sval);
			numRepressors++;
		}
		else if(info.equals("RepressorParam")) {
			GeneralInput.nextToken(tokenizer);
			Parameter p = Affector.getParameter(tokenizer.sval);
			if(p == null) throw new Exception("Couldn't find repressor parameter " + tokenizer.sval);
			params.addParameter(p);
		}
		else if(info.equals("RepressionTime")) {
			GeneralInput.nextToken(tokenizer); 
			repressionTime = (int)tokenizer.nval; 
			switchStop.setStopTime(repressionTime + 10);	// Note - assumption here that 10 minutes is long enough for repression to decay away
		}
		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("NumCycles")) {GeneralInput.nextToken(tokenizer); numCycles = (int)tokenizer.nval; }
		else if(info.equals("ExternalNodeValue")) {GeneralInput.nextToken(tokenizer); externalNodeValue = (float)tokenizer.nval; }
		else if(info.equals("Mode")) { 
			GeneralInput.nextToken(tokenizer);
			if(tokenizer.sval.equals("Repression")) mode = REPRESSORS;
			else mode = ACTIVATORS;
		}
		else super.loadParameter(info, tokenizer);
	}

	public void doRun()  {
		int			i, n, c;
		float		dist, step_size;
		boolean		is0, isHigh;
		Cell		cell;
		float	[]	values = new float[numNodes];
		float	[]	node_values = new float[numNodes];
		
		reset();
		for(i = 0; i < numRepressors; i++) numSwitches[i] = 0;
		
		// Note : in this code I am assuming that AC, ac, SC, and sc are the first four node values,
		// with the protein and rna of each gene in consecutive slots (order of AC and SC doesn't matter).
		
		// First find the upper steady state point. Do this by setting nodes to 1,
		// running model for a little while, then zooming in with newton's algorithm
		// Throughout here, I'm assuming there is just a single cell
		cell = cells[0];
		for(n = 0; n < numNodes; n++) {
			cell.setInitialValue(nodes[i], 1.0f);
		}

		// Get rid of any repression
		for(n = 0; n < numRepressors; n++) {
			cells[n].setInitialValue(repressors[n], 0.0f);
		}

		// Find out what the upper steady state is
		SimpleStop stopper = getStopper();
		// The stopper should be simple stop, going for a short time to get near st st
		setStopper(ststStop);
		F(p);	// Don't care about value
		setStopper(stopper);
		if( ! moveToStSt(cell, numNodes, values)) {
			System.out.println("Couldn't find upper steady state for this parameter set");
			finalScore = 0;
			return;
		}

		// Turn each external node on in one of the cells
		// Assume that there are at least enough cells to
		// use one cell per external node so can do all nodes in each run
		
		for(n = 0; n < numRepressors; n++) {
			cells[n].setInitialValue(repressors[n], externalNodeValue);
		}

		// if checking repression, set each cell to the high steady state
		// otherwise set each cell to 0
		for(n = 0; n < numNodes; n++) {
			for(c = 0; c < numRepressors; c++) {
				if(mode == REPRESSORS)
					cells[c].setInitialValue(nodes[n], values[n]);
				else cells[c].setInitialValue(nodes[n], 0f);
			}
		}
		
		// Set up the switchStop
		try {
			switchStop.clearNodes();
			for(n = 0; n < numNodes; n++) {
				if(mode == REPRESSORS) switchStop.addNode(network.getNode(nodes[n]).getName(), 0);
				else switchStop.addNode(network.getNode(nodes[n]).getName(), values[n]);
			}
		} catch(Exception e) {
			System.out.println("problem in iterator: " + e.toString());
		}
		
		//setStopper(simpleStop);
		setStopper(switchStop);
		// Now pick random parameter values for each of the repressor parameters and see
		// whether each repressor can switch its cell to the off state
		for(i = 0; i < numCycles; i++) {
			
			// Pick a random set of parameters and set the model to them
			params.makeNewPoint(MoreMath.UNIFORM);
			params.setModel();
			
			// Run the model
			F(p);

			for(c = 0; c < numRepressors; c++) {
				// Get to the closest stable steady state
				boolean at_st_st = false;
				float stop_time = switchStop.getStopTime();
				while(!at_st_st && switchStop.getStopTime() < stop_time + 2000) {	// If can't get there in 2000 minutes, give up
					switchStop.setStopTime(switchStop.getStopTime() + 100);
					continueF(p);
					moveToStSt(cells[c], numNodes, node_values);
				
					// Check to make sure that we went to one of the stable steady states.
					at_st_st = true;
					n = 0;
					while(n < numNodes && at_st_st) {
						if(node_values[n] > toleranceAbs && node_values[n] < values[n] - toleranceAbs)
							at_st_st = false;
						n++;
					}
				}
				switchStop.setStopTime(stop_time);
			
				// If these are repressors, see if any of the cells have switched down near 0
				if(mode == REPRESSORS) {
					is0 = true;
					for(n = 0; n < numNodes; n++) {
						if(node_values[n] > toleranceAbs) {
//							if(is0)
//								System.out.println(values[0] + "\t" + values[2] + "\t" + node_values[0] + "\t" + node_values[1] +
//										"\t" + node_values[2] + "\t" + node_values[3]);
							
							is0 = false;
						}
					}
					if(is0) numSwitches[c]++;
				}
				else {	// Activators, so see if any cells have switched up above high steady state
					isHigh = true;
					for(n = 0; n < numNodes; n++) {
						if(node_values[n] < values[n] - toleranceAbs) {
//							if(isHigh) 	
//								System.out.println(values[0] + "\t" + values[2] + "\t" + cells[c].getValue(nodes[0]) + "\t" + cells[c].getValue(nodes[1]) + "\t" +
//												cells[c].getValue(nodes[2]) + "\t" + cells[c].getValue(nodes[3]));
							isHigh = false;
						}
					}
					if(isHigh) numSwitches[c]++;
				}
			}
		}

		resetStopper();
		
		finalScore = 1;
		for(c = 0; c < numRepressors; c++) finalScore *= numSwitches[c] / (float)numCycles;

		super.stopRun();
	}
	
	public boolean didPass() { return true; }
	
	
	private final boolean moveToStSt(Cell cell, int num_nodes, float [] ststvalues) {
		int			i;
		float		dist, step_size;
		boolean		done, bad;
		float	[]	deriv = new float[numNodes];
		int		[]	index = new int[numNodes];
		float	[]	matrix = new float[numNodes*numNodes];
		float	[]	orig_values = new float[num_nodes];

		for(i = 0; i < num_nodes; i++) {
			ststvalues[i] = cell.getValue(i);
			orig_values[i] = cell.getValue(i);
		}
		dist = 0;
		for(i = 0; i < num_nodes; i++) dist += sqr(ststvalues[i]);
		dist = (float)Math.sqrt(dist);
		step_size = 0.1f;
		float last_dist2 = 1000000f;
		int	num_times = 0;
		done = false; bad = false;
		while(!done && num_times < 100) {
			// Figure out the move we need to make
			MoreMath.calcJacobian(matrix, cell, num_nodes, epsilon);
			for(i = 0; i < num_nodes; i++)
				deriv[i] = -cell.nodes[i].getAffectorsValue(0);
			ludcmp(matrix, num_nodes, index);
			lubksb(matrix, num_nodes, index, deriv);
			
			float dist2 = 0;
			for(i = 0; i < num_nodes; i++) dist2 += sqr(deriv[i]);
			dist2 = (float)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 < num_nodes; i++)
						deriv[i] *= step_size * dist / dist2;
				}
				// Make the move
				for(i = 0; i < num_nodes; i++) {
					ststvalues[i] += deriv[i];
					cell.setValue(i, ststvalues[i]);
				}
			}
			last_dist2 = dist2;
			num_times++;
		}	// while(!done)

		// Set the values in the cells back to where they were to start with
		for(i = 0; i < num_nodes; i++) {
			cell.setValue(i, orig_values[i]);
		}

		if( bad || num_times > 99 ) return false;
		return true;
	}

	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 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];
		}
	}
	
}