package iterators;

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

/** <p> The VPCDoseResponseIterator is a special iterator class used to look for Parameter sets for
 * which the dose response curve for activated Map Kinase in response to increasing EGF input
 * exhibits an initial plateau, an intermediate level plateau, and highest response.  The iterator
 * computes the dose response curve and then finds inflection points along it, i.e. points where the
 * second derivative of the response changes sign (thus indicating local maxima and minima for the
 * first derivative function.
 * <p>	Should be run with a stopping condition that detects stable solutions.
 * <p> The iterator rejects this parameter set if:
 *
 * 	<blockquote>  The stopper rejects any one of the solutions.
 * 	<blockquote>  	the difference between the response to high EGF and the response to low EGF
 * 		falls below HighLowThreshold. </blockquote>
 * 	<blockquote>  	The number of inflection points is less than 4. </blockquote>
 * 	<blockquote>  	the slope at the minima (i.e.the bottom of the troughs) exceeds the minFlatness. </blockquote>
 *
 * <p>  If these criteria are met, the iterator writes the does response curve to its local output file and
 * calculates a score equal to the sum of the inverses of the trough widths (the fraction of the total
 * parametter value interval for which the derivative is less than minflatness, and the differences between
 * medium and low values, and high and medium values, scaled by the difference between high and low values.
*/
public class VPCDoseResponseIterator extends ModelIterator implements Runnable {
	
	static int MAX_INFLECTIONS = 10;
		

	/** The number of steps to break the parameter interval into. */
	int						numSteps=100;
	
	/** The threshold stopper value above which the iterator considers trajectory for each
	 control par value okay.
	 */
	float 					stabilityThreshold = 0.1f;
	
	/** The threshold value that determines whether or not the iterator passes the current parameter set. */
	float 					scoreThreshold = 4.0f;
	
	float 					minFlatness;
	Node 					targetNode;
	String 					nodeName;
	PrintWriter				localOut=null;
	String					outfileName;
	float					highLowThreshold;
	
	public VPCDoseResponseIterator() {	}	// Have a no argument initializer so we can do new_by_name

	public ModelIterator copy()  throws Exception {
		VPCDoseResponseIterator newIter = new VPCDoseResponseIterator();
		newIter.init(this.network, this.model);
		newIter.nParsTV = this.nParsTV;
		newIter.parsTV = this.parsTV.copy();
		newIter.theFunction = this.theFunction.copy();
		return newIter;
	}


	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.
		try {
			targetNode = model.cells[0].getNode(nodeName);
		}catch(Exception e) {
			throw new Exception("PlateauIterator failed to find a Node corresponding to " + nodeName);
		}
			
		
		Calendar	calendar = new GregorianCalendar();
		outfileName = outfileName + "_" + (calendar.get(Calendar.MONTH) + 1) + "_" + calendar.get(Calendar.DAY_OF_MONTH) + "_" + calendar.get(Calendar.HOUR_OF_DAY) + "_" + calendar.get(Calendar.MINUTE) + "_" + calendar.get(Calendar.SECOND);
		try {
			localOut = new PrintWriter(new FileOutputStream(outFileName), true);
		} catch(IOException e) { throw new Exception ("Error creating TransectIterator output file: " + e.toString()); }
	}
	
	
	// 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("HighLowThreshold")) {GeneralInput.nextToken(tokenizer); highLowThreshold = (float)tokenizer.nval; }
		else if(info.equals("MinFlatness")) {GeneralInput.nextToken(tokenizer); minFlatness = (float)tokenizer.nval; }
		else if(info.equals("NumSteps")) {GeneralInput.nextToken(tokenizer); numSteps = (int)tokenizer.nval; }
		else if(info.equals("TargetNode")) {GeneralInput.nextToken(tokenizer); nodeName = tokenizer.sval; }
		else if(info.equals("PlateauFile")) {GeneralInput.nextToken(tokenizer); outfileName = tokenizer.sval; }
		else if(info.equals("scoreThreshold")) {GeneralInput.nextToken(tokenizer); scoreThreshold = (float)tokenizer.nval; }
		else super.loadParameter(info, tokenizer);
	}

	/** Run the model at a incrementing set of values of one particular
	 * parameter value. Only the first parameter in the parameter
	 * array is used.
	 */
	public void doRun()  {
		float	low, high, incr, closest_dist, val;
		int		i, closest;
		
		//	Values of the parameter for which we evaluate model performance.
		float []	parVals = new float[numSteps+1];;
		
		// Model performance scores for each of the sampled parameter values.
		float 	[] 	nodeVals = new float[numSteps+1];
	
		reset();
		low = parsTV.getLowerBound(0);
		high = parsTV.getUpperBound(0);
		
		if(parsTV.getVariationMode(0) == Parameter.LOGARITHMIC) {
			low = (float)Math.log(low);
			high = (float)Math.log(high);
		}
		incr = (high - low) / numSteps;

		for(i = 0, val=low; i <= numSteps; i++, val+= incr) {
			if(parsTV.getVariationMode(0) == Parameter.LOGARITHMIC)
				p[0] = (float)Math.exp(val);
			else p[0] = val;
			parVals[i] = p[0];
			
			// If not stable, we don't care about the rest of the parameter range
			if(F(p) > stabilityThreshold) {
				finalScore = 10000;
				super.stopRun();
				return;
			}
			
			
			// Store the value of the node
			nodeVals[i] = targetNode.getValue(0);
		}
		// Now that we have stable values of the node for each value of the parameter,
		// Check whether the difference in the highest and lowest values of
		// the node are larger than a threshold
//System.out.println("high = " + nodeVals[numSteps] + " low = " + nodeVals[0] + " thr = " + highLowThreshold);
		if(nodeVals[numSteps] - nodeVals[0] < highLowThreshold) {
			finalScore = 10000;
//			super.stopRun();
			return;
		}
		
		// compute a scaled min flatness relative to the average slope over the interval.
		float scaledMinFlatness = minFlatness*(nodeVals[numSteps] - nodeVals[0])/(high - low);
		// Find the derivative curve
		float[]	derivatives = new float[numSteps];
		
		for(i = 0; i < numSteps; i++) {
			derivatives[i] = (nodeVals[i+1] - nodeVals[i])/incr;
		}
//System.out.println("firstDerivative = " + derivatives[0] + " diff = " + (derivatives[1] - derivatives[0]) + " SMF = " + scaledMinFlatness);
		
		// we expect a sufficiently flat minimum at 0.  Otherwise, done;
		if(derivatives[0] > scaledMinFlatness || derivatives[1] - derivatives[0] < 0) {
			finalScore = 10000;
//			super.stopRun();
			return;
		}
			
		
		int [] inflectionPoints = new int[MAX_INFLECTIONS];
		
		// count the first value as an inflection (a minimum)
		inflectionPoints[0] = 0;
		int numInflections = 1;
		boolean rising = true;
		float max = derivatives[0];
		float min = derivatives[0];
		
		// now search for and log consecutive maxima and minima.
		for(i = 1; i < numSteps; i++) {
			if(rising)  {
				if(max - derivatives[i] > scaledMinFlatness) {  // found a maximum
					inflectionPoints[numInflections] = i;
					numInflections++;
					rising = false;
					min = derivatives[i];
				}
				else {
					max = Math.max(max,derivatives[i]);
				}
			}
			else {
				if(derivatives[i] - min > scaledMinFlatness) {	// found a minimum
					inflectionPoints[numInflections] = i;
					numInflections++;
					rising = true;
					max = derivatives[i];
				}
				else {
					min = Math.min(min,derivatives[i]);
				}
			}
		}
		if(rising) { // we call the last point a maximum
			inflectionPoints[numInflections] = i-1;
			numInflections++;
		}
//System.out.println("numInflections = " + numInflections);
		// Reject if we don't get two clearly defined sufficiently flat minima
		if(numInflections != 4 || derivatives[inflectionPoints[2]] > scaledMinFlatness)  {
//if(derivatives[inflectionPoints[2]] > scaledMinFlatness)
//	System.out.println("not flat enough");
			finalScore = 10000;
//			super.stopRun();
			return;
		}
		
		// now go back and calculate a specific score
		
		float binWidth = 1f/numSteps, width1 = 0, width2 = 0;
		
		
		i = 0;
		while(i < numSteps && derivatives[i] < scaledMinFlatness) {
			width1 += binWidth;
			++i;
		}

		i = inflectionPoints[2];
		while(i >= 0 && derivatives[i] < scaledMinFlatness) {
			width2 += binWidth;
			--i;
		}
		i = inflectionPoints[2] + 1;
		while(i < numSteps && derivatives[i] < scaledMinFlatness) {
			width2 += binWidth;
			++i;
		}
		
		// height of the jump between plateaus scaled by total jump from minVal to maxVal
		float height1 = (nodeVals[inflectionPoints[2]] - nodeVals[0])/(nodeVals[numSteps] - nodeVals[0]);
		float height2 = (nodeVals[numSteps] - nodeVals[inflectionPoints[2]])/(nodeVals[numSteps] - nodeVals[0]);
		
		
		// making any one of these too small kills the score
		finalScore = 0.3f/width1 + 0.3f/width2 + 0.5f/height1 + 0.5f/height2;
//System.out.println("finalScore = " + finalScore);
		// Save to special output file
		localOut.print(parsTV.getName(0) + "\t");
		for(i = 0; i <= numSteps; i++) {
			localOut.print(parVals[i] + "\t");
		}
		localOut.println();
		localOut.print(parsTV.getName(0) + "\t");
		for(i = 0; i <= numSteps; i++) {
			localOut.print(nodeVals[i] + "\t");
		}
		localOut.println();
			
		
		// Leave model in its original state - p will be set with orig values now
		parsTV.setFrom(p);
		parsTV.setModel();
//		super.stopRun();
	}
	
	/** Override's ModelIterator's didPass method so that it can compare an internally computed score to
	 * it's scoreThreshold.
	 */
	public boolean didPass()  {
		return finalScore < scoreThreshold;
	}
}