package iterators;

import java.io.*;
import main.BetterTokenizer;
import main.Model;
import main.GeneralInput;
import affectors.Affector;
import parameters.ParameterSet;
import java.util.Random;

import java.util.Calendar;
import java.util.GregorianCalendar;

public class EllipIterator extends ModelIterator {
	final static float hhuge = 1.7758E30f;
	final static float tiny = 2.25E-12f;
	static String whystop[] = {
    	new String("k > Kmax with feasible point found "),
    	new String("k > Kmax with no feasible point found"),
		new String("Q matrix is non Positive-Definite, but feasible point found"),
		new String("Q matrix is non Positive-Definite -- no feasible point found "),
		new String("Gradient vector = 0, but feasible point found "),
		new String("Gradient vector = 0, -- no feasible point found "),
		new String(""),
	};
	float [] delta;	// Used in gradients
	
	float GEPS = 0.001f;	// Eli, this is an important number. For reasons I do not understand
   										// the performance of this ea3() method depends stronglyh on GEPS for 
   										// the test problemincluded in this file
	
	int		MAXITT = 150;		// assert the maximum number of iterations you will allow the ea3() algorithm to make
	int		numCycles = 1;		// The number of times to call ea3 per call to this optimizer
								// supposedly, calling ea3 several times for a few iterations is more efficient than calling it once for many iterations
	float	finishThreshold = 0.0f;	// If value of function gets below this, finish
	float	metaBestScore = 9999999;
	int		metaNumFuncEvals = 0;

// I use a Random object so that I can set the seed for debugging purposes
	Calendar calendar = new GregorianCalendar();
	Random					rng = new Random(calendar.get(Calendar.SECOND) * calendar.get(Calendar.MINUTE));
	
	public EllipIterator() {}
		
	public void reset()  {

// ELI 5/28 - using super.reset() instead of some of this
/*		if(verbose)  {
			try {
				ps = new PrintStream(new FileOutputStream(outFileName));
			} catch(IOException e) { ps = null; }
		}		
		nParsTV = parsTV.getNumParams();	// set local values of nParsTV 
		p = parsTV.getParVals();  			// and p from parsTV
		bombed = false;
		numFunctionCalls = 0;
*/		super.reset();
		metaNumFuncEvals = 0;
		delta = new float[nParsTV];
		metaBestScore = 9999999;
	}
	

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

		
		
	
	// 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("GEPS")) {GeneralInput.nextToken(tokenizer); GEPS = (float)tokenizer.nval; }
		else if(info.equals("MaxIterations")) {GeneralInput.nextToken(tokenizer); MAXITT = (int)tokenizer.nval; }
		else if(info.equals("NumCycles")) {GeneralInput.nextToken(tokenizer); numCycles = (int)tokenizer.nval; }
		else if(info.equals("RandomSeed")) {GeneralInput.nextToken(tokenizer); rng.setSeed((int)tokenizer.nval); }
		else if(info.equals("FinishThreshold")) {GeneralInput.nextToken(tokenizer); finishThreshold = (float)tokenizer.nval; }
		else super.loadParameter(info, tokenizer);
	}


	
	/*****************************************************************
	// ******* This version of theFunction for diagnostic purposes only
	public float F(float [] thePoint) {
		double score,x3eight,sinx3eight,fac0,fac1,fac2,fac3;

		numFunctionCalls++;

		x3eight= 8.0 * thePoint[3];
		sinx3eight = Math.sin(x3eight);
		sinx3eight = thePoint[3] * sinx3eight*sinx3eight;
		fac0 = 5.0 * (thePoint[0] - sinx3eight);
		score = fac0*fac0; 																
		fac1 = 5.0 * (thePoint[1] - sinx3eight);
		score += fac1*fac1; 
		fac2 = Math.cos(x3eight);
		fac2 = 5.0 * (thePoint[2] - fac2*fac2);
		score += fac2*fac2;
		fac3 = 0.5 * Math.abs(thePoint[3]);
		score += fac3;  
		return (float)score;
	}
	// ******* This version of theFunction for diagnostic purposes only
	//*****************************************************************/
	

//ELI 3/23 Changed from run to doRun
	public void doRun() {
		float		f_best;
		int 		n_evals;		// a counter of the number of function evaluations made
		int 		n_pars;			// the number of parameters to be optomized
		int 		n_constraints;	// the number of constraints to be imposed
		int 		return_code;	// index into the whystop[] array telling the reason the ellipsoid algorithm quit

		int 		i,kpass;		// junk storage for this demo main() function
		float	 	initEllipSize = 2.0f;
																

	// The initial guess values of the parameters are in p	
	
		reset();  // reload p, nParsTV, delta because these may have changed
		float [] x_guess = new float[nParsTV];	
		float [] x_best = new float[nParsTV];														
		for(kpass=0; kpass < numCycles && metaBestScore >= finishThreshold; kpass++) {
			for(i = 0; i < nParsTV; i++) x_best[i] = p[i]; // transfer best pars found so far into initial guess variable for ea3();
							  	
			metaNumFuncEvals=0; // reset meta function eval num counter.
			//fprintf(stderr,"\n\n\n\n Fresh start (pass %d) for ea3()...",kpass);
		 
			// the next line calls the ea3() function
			return_code = ea3(initEllipSize, MAXITT,x_best);
			initEllipSize = .5f*initEllipSize; // on next pass, use smaller ellipse
		
	    // by code in my funcs() function, the bes parameter set so far visited will
		// be in metaBestPars[]. The ea3() function does not do this.. The next kpass
		// iteration will reload these best parameter values into the x_guess array

	    }
		finalScore = metaBestScore;		// finalScore now holds the final score
		System.out.println("Minimum score found = " + metaBestScore + " at:");
		for(i = 0; i < nParsTV; i++)
			System.out.println(parsTV.getName(i) + " = " + p[i]);
		super.stopRun();
 	}


	/************************ ea3() ******************************/
		/* 
		initEllipSize specifies the "size" of the initial ellipse surrounding the 
				initial guess point.
		kmax is the maximum number of iterations allowed (supplied by user).
		xc[1],...,xc[n] is the user--supplies initial guess, and is the 
		         array that returns the best answer found.
		*k is the number of iterations the algorithm makes. To start algorithm,
			user supplies *k=0. The iteration count is returned by ea3() in *k.
		*fr is the least vbalue found for the function to be minimized = value
			of fcn() at xr[] (returned by ea3()).
								*/
	/* Translation of Kupferschmid/Ecker ellipsoid algorithm... */
	private int ea3(float initEllipSize, int kmax, float [] xc) {
		float fr = hhuge;	// These two variables used to be passed in as pointers (in C code)
		int k = 0;		// Since I couldn't see that we'd use them for anything, it was easier to just
						// make them local for now.
	
		float		fn,mu,e,c,f,gmax,gqc,dum;
		int			rc = 1,i,ii,j,jj,l,jjj;
		int			someConstraintViolated;
			// Note - q's indexing starts at 1, not at 0, like fortran arrays - 
			// the logic that the original author used was too complicated for
			// me to decipher, so I've left it and hope it works.
		float []	q = new float[(nParsTV+1)*(nParsTV+1)];
		float []	g = new float[nParsTV],
					d = new float[nParsTV], 
					xr = new float[nParsTV];

	// check for sensible input
		if(kmax <= 0 ) return (rc=8);

	// xr[0],...,xr[n-1] is the vector of the best set of parameters found (returned by ea3()).

		initEllip(xc,initEllipSize, q); // set up initial ellipse
		for(j = 0; j < nParsTV; j++) xr[j] = xc[j];
		
	// Find constants used repeatedly in the updates
		fn = nParsTV;
		mu = 1.0f/(1.0f+fn);
		e = 0.25f;
		c = 0.0f;
		if(nParsTV > 1) {
		    e = fn*fn/(fn*fn - 1.0f);
		    c = 2.0f*mu;
		}
	
	// prepare to examine the constraints in cyclical order
		i=0;
	
		while((k++ <= kmax) && (metaNumFuncEvals <= kmax) && metaBestScore > finishThreshold) { // New version uses count of func evals
									// generated both by gradients() function as well as by funcs() function
			
		    // begin the next iteration
	
			someConstraintViolated=0;
			// Find a violated inequality if there is one
			ii=0;
			while((ii < nParsTV) && someConstraintViolated==0) {
			    i++;
			    if(i > nParsTV) i = 1;
			    if(ea3Function(xc, i) >= 0) someConstraintViolated = i;
				++ii;
			}
								 
			if(someConstraintViolated==0) { // all constraints are currently met.
			    rc=0;        
			    i = 0; // this picks objective function to be minimized
			    if( (f=ea3Function(xc,i)) < fr ) { // Update record value and the record point
			        fr=f; // new best record point found
			        for(j=0; j < nParsTV; j++) xr[j] = xc[j];
				}
				System.out.print("Eval " + metaNumFuncEvals + " @ (" + xc[0]);
				for(int m = 1; m < nParsTV; m++)
					System.out.print(", " + xc[m]);
				System.out.println(") = " + f);
			}
			else i = someConstraintViolated;
		
		// Construct the next ellipsoid ... find the gradient and its absolutely largest element
			ea3Gradients(xc, i, g);
		
			gmax = 0;
			for(j=0; j < nParsTV; j++)
			    if( (dum=Math.abs(g[j])) > gmax) gmax = dum;
		// check that the gradient is not (effectively) the zero vector
			if(gmax <= .0000001) {
			      for(j = 0; j < nParsTV; j++) xc[j] = xr[j];    	
			      return ( (rc+=4));
			}
		
		// Normalize the gradient vector so largest element has absolute value equal to 1.0
			for (j = 0; j < nParsTV; j++) g[j] = g[j] / gmax;
		
		// find the direction qg and the scale factor gqc
			gqc = 0.0f;
			for(j = 1; j <= nParsTV; j++) {
				d[j-1] = q[1 + ( (j-1)*j)/2]*g[0];
				if(nParsTV > 1) {
					for(jj = 2; jj <= nParsTV; jj++) {
						if(jj<j) jjj = jj + ( (j-1)*j)/2;
						else jjj = j + ( (jj-1)*jj)/2;
						d[j-1] = d[j-1] + q[jjj]*g[jj-1];
					}
				}
				gqc += g[j-1]*d[j-1];
			}

	
		// check to see whether Q has become non-pd
			if(gqc < tiny) {
			    for(j = 0; j < nParsTV; j++) xc[j] = xr[j];    	
			 	return( (rc += 2));
			}        
			    
		// scale the direction vector
			gqc=(float)(1./Math.sqrt(gqc));
			for(j = 0; j < nParsTV; j++) d[j] = d[j] * gqc;
		
		// update xc and q
			l = 1;
			for(j = 0; j < nParsTV; j++) {
			    xc[j]=xc[j]-mu*d[j];
			    for(jj=0; jj <= j; jj++)   {
					q[l] = e*(q[l]-c*d[j]*d[jj]);
					l++;
			    }
			}
		}	// end of while k < kmax
	
		for(j = 0; j < nParsTV; j++) xc[j] = xr[j];
		return (rc);
	}
	
	private float ea3Function(float [] pt, int funcNum) {
		if(funcNum == 0 || funcNum > nParsTV) {
			metaNumFuncEvals++;
			float score = F(pt);

			if(score < metaBestScore) { // Yes! only if it's a feasible point
				for(int jcon = 1; jcon <= nParsTV; jcon++)
					if(ea3Function(pt, jcon) > 0)  // Note this function calls itself
						return score;
				 // if fall through, all the constraints are satisfied
				for(int i = 0; i < nParsTV; i++) p[i] = pt[i];
				metaBestScore = score;				 
			}
			return score;
		}
		
		funcNum -= 1;
		if(pt[funcNum] < parsTV.getLowerBound(funcNum))
			return parsTV.getLowerBound(funcNum) - pt[funcNum];
		else if(pt[funcNum] > parsTV.getUpperBound(funcNum))
			return pt[funcNum] - parsTV.getUpperBound(funcNum);
		else return -tiny;	// Just return a small negative number, so there's not too much of a discontinuity
	}


/******************************  gradients() **********************/
/* sample function that returns gradientsof the objective and constraint functions.
ea3() calls this function repeatedly, supplying defined values of:
	x[1],x[2],....,x[n] is the parameter vector. 
	funcNum = the number of the function ea3() wants the gradient of at the point x[].
   	for funcNum=1,...n_constraints, this routine must return in the g[] array
		the gradient of the funcNum-th constraint function.  
	for funcNum == n_constraints+1, or funcNum==0 this routine must return in the g[] array
   		the gradient of the objective function.  
 */
	private void ea3Gradients(float [] x, int funcNum, float [] g) {
		float []	xdiff = new float[nParsTV];
		float		fdiff,dum,fctr;
		int 		numGrads = 0;

		int j;
	
// this demo routine illustrates how to compute numerically the gradient of the objective function.
    
    	if(funcNum == (nParsTV+1) || funcNum == 0) {	// return gradients of objective function in the g[] array
			for(j = 0; j < nParsTV;j++) {
				// randomize the perturbations used in computing gradients each pass through
				dum = (float)((1.1 + 2.*(rng.nextFloat()-.5))*GEPS*(.1 + Math.abs(x[j])));
						 /* so the perturbation away from x[j] could be as small as
						 .01*GEPS and could be as large as 2.1*EPS*(.1 + |x[j]|)    
						 which, for GEPS==.001,  is .21% be as large as the parameter's absolute value. */
				if(rng.nextFloat() > 0.5) delta[j] = dum; /*flip the sign at random */
				else delta[j] = (-dum); 
			}
		     
			// this initializes the array of infinitesimals...  Using same deltaX for each parameter.  Probably not generally a good idea...
        	++numGrads;
			//fprintf(stderr,"G{%d}",numGrads);  /* announce computation of gradient */

			for(j = 0; j < nParsTV; j++) xdiff[j] = x[j];
			// get function value at base point
			fctr = ea3Function(x, 0);
			for(j = 0; j < nParsTV; j++) {
				xdiff[j] += delta[j];
				// Check the constraints for that parameter to make sure our random addition didn't violate
				if(ea3Function(xdiff, j + 1) > 0) {
					xdiff[j] -= delta[j];
					delta[j] = -delta[j];	// Just go the other way, and assume that if we're at one boundary we can't be near the other boundary
					xdiff[j] += delta[j];
				}
				fdiff = ea3Function(xdiff, 0);
				xdiff[j] -= delta[j];
				g[j] = (fdiff-fctr)/delta[j];				 
			}
		}

		else   {    // gradients of the constraint functions
			funcNum--;
			for(j = 0; j < nParsTV;j++) g[j] = 0.0f;  // start with all zero components of the gradient vector
			
			if(x[funcNum] < parsTV.getLowerBound(funcNum)) g[funcNum] = -1;
			else if(x[funcNum] > parsTV.getUpperBound(funcNum)) g[funcNum] = 1;
		}

	}


	/************************ initEllip() ******************************/
	/* this defines an initial ellipsoidal region that ea3() 
	iteratively subdivides. You must call this function before calling ea3() */
	private void initEllip(float [] x_guess, float size, float [] qqq) {
		int		mi,i,j;							  
		float	dum = nParsTV;

	// Initialize qqq to 0 - not in original code, but seems neccesary to me ?? anyway, can't hurt
		for(i = 0; i < (dum+1)*(dum+1); i++) qqq[i] = 0f;
	

	// Note - as mentioned above, qqq index starts at 1, like fortran.  Other arrays start at 0
		qqq[1]=dum*(0.1f + size*x_guess[0]*x_guess[0] );
		j = 1;
		for(i = 2; i <= nParsTV; i++) {
			for(mi = 1; mi < i; mi++) {
				++j;
				qqq[j]=0.0f;
			}
			++j;
			qqq[j]=dum*(0.1f+ size*x_guess[i-1]*x_guess[i-1]);
		}
	}
}