package main;

import java.util.Random;
import java.util.Date;

import main.Cell;

/** This class includes a bunch of extra mathematical functions and other math
tidbits that are useful to have around. All the functions are static and all
the important ones are public, so you never actually need to instantiate this class. */
public class MoreMath extends Object {
	public static final int		UNIFORM=0, GAUSSIAN=1, BROWN_NOISE=2;
	private	static Random	rng = new Random();
	
	static
	{
		setRandomSeed((int)((new Date()).getTime() % 100000000));
	}

	/** A record class that's handy when assembling a Jacobian matrix from the Affectors'
	formulas rather than numerically. */
	public class PartialDerivRec extends Object	{
		float				dfdy;
		int					CellNum;
		int					NodeNum;
		int					side;
	}
	
	/** Returns num * num */
	static public double sqr(double num) { return num * num; }
	static public float sqr(float num) { return num * num; }
	
	/** Returns the standard deviation. Sum is the sum of the values, sum_sqr is the
	sum of the squares of the values, and num is the number of values. */
	static public double getStdDev(double sum, double sum_sqr, int num) {
		return Math.sqrt((sum_sqr - (sum * sum / num)) / (num - 1.0));	
	}
	
	static public float getStdDev(float sum, float sum_sqr, int num) {
		return (float)Math.sqrt((sum_sqr - (sum * sum / num)) / (num - 1.0));	
	}
	
	static public void setRandomSeed(int seed) { rng.setSeed(seed); }

	static public float random() { return random(UNIFORM); }

	static public float random(int mode) { 
		if(mode == UNIFORM) return rng.nextFloat(); 
		else if(mode == GAUSSIAN) return (float)rng.nextGaussian();
		else return rng.nextFloat();
	}
	
	static public float randNeg(int mode)	{
		if(mode == UNIFORM)	{
			float f = rng.nextFloat();
			if(f < 0.5f)	{
				return -(rng.nextFloat());
			}
			else {
				return rng.nextFloat();
			}
		}
		else if(mode == GAUSSIAN)	{
			return (float)rng.nextGaussian();
		}
		else return rng.nextFloat();
	}

	/** Returns an integer between min and max - 1. */
	static public int randomInt(int min, int max)
	{
		int num = (int)(min + (max - min) * random(UNIFORM));
		if(num == max) num = max - 1;
		return num;
	}

	static public float pickRandomLinear(float min, float max) {
		return(min + (max - min) * random(UNIFORM));
	}
	
	static public float pickRandomLogarithmic(float min, float max) {
		double expRange = Math.log((double)max) -  Math.log((double)min);
		double exp = (random() * expRange) + Math.log((double)min);
		return (float)Math.exp(exp);
	}


	/** Randomizes the elements in the array. */
	static public void randomize(int [] array)
	{
		for(int i = 0; i < array.length; i++)
		{
			int rnd = (int)(random() * (array.length - i));
			int temp = array[i + rnd];
			array[i + rnd] = array[i];
			array[i] = temp;
		}
	}
	

	/** Calculates the jacobian matrix for the equations in a network around
	the Node values of the input cell. epsilon determines the amount to change
	each Node when finding the derivative. */
	public static void calcJacobian(float [] matrix, Cell cell, int num_nodes, float epsilon) {
		int		n1, n2;
		float	val, valplus, valminus;

		for(n1 = 0; n1 < num_nodes; n1++) {
			for(n2 = 0; n2 < num_nodes; n2++) {
				val = cell.getValue(n1);
				cell.setValue(n1, val + epsilon);
				valplus = cell.nodes[n2].getAffectorsValue(0);
				cell.setValue(n1, val - epsilon);
				valminus = cell.nodes[n2].getAffectorsValue(0);
				matrix[n2 * num_nodes + n1] = (valplus - valminus) / (2.0f * epsilon);
				cell.setValue(n1, val);
			}
		}
	}
	


	/** Numerical recipes in C code for finding eigenvalues of a general, real, non-symmetric matrix.

	 Note - this function destroys the matrix a.
	 a is the square input matrix, n is the length of a, wr are the real parts of the eigenvalues,
	 wi are the imaginary parts.*/
	public static void eigenvalues(float [] a, int n, float [] wr, float [] wi) throws Exception {
		balanc(a, n);
		elmhes(a, n);
		hqr(a, n, wr, wi);
	}

	private static void hqr(float [] a, int n, float [] wr, float [] wi) throws Exception {
		int nn, m, l, k, j, its, i, mmin;
		float	z, y, x, w, v, u, t, s, r=0.0f, q=0.0f, p=0.0f, anorm;
		
		anorm = 0.0f;
		for(i = 0; i < n; i++) {
			j = i-1; if(j < 0) j = 0;
			for(; j < n; j++) anorm += Math.abs(a[i*n+j]);
		}
		nn = n-1;
		t = 0.0f;
		while(nn >= 0) {
			its = 0;
			do {
				for(l = nn; l >= 1; l--) {
					s = Math.abs(a[(l-1) * n + l-1]) + Math.abs(a[l*n+l]);
					if(s == 0.0f) s = anorm;
					if((float)(Math.abs(a[l*n+l-1]) + s) == s) break;
				}
				x = a[nn * n + nn];
				if(l == nn) {
					wr[nn] = x + t;
					wi[nn--] = 0.0f;
				}
				else {
					y = a[(nn - 1) * n + nn - 1];
					w = a[nn*n+nn-1] * a[(nn-1)*n+nn];
					if(l == (nn-1)) {
						p = 0.5f *(y-x);
						q = p*p+w;
						z = (float)Math.sqrt(Math.abs(q));
						x += t;
						if(q >= 0.0) {
							z = p + SIGN(z,p);
							wr[nn-1] = wr[nn] = x + z;
							if(z != 0) wr[nn] = x - w/z;
							wi[nn-1] = wi[nn] = 0.0f;
						}
						else {
							wr[nn-1]=wr[nn]=x + p;
							wi[nn-1] = -(wi[nn]=z);
						}
						nn -= 2;
					}
					else {
						if(its == 30) throw new Exception("Too many iterations in hqr");
						if(its == 10 || its == 20) {
							t += x;
							for(i = 0; i < nn; i++) a[i*n+i] -= x;
							s = Math.abs(a[nn*n+nn-1]) + Math.abs(a[(nn-1)*n+nn-2]);
							y=x=0.75f*s;
							w = -0.4375f*s*s;
						}
						++its;
						for(m=(nn-2); m >=l; m--) {
							z = a[m*n+m];
							r = x-z;
							s = y-z;
							p = (r*s - w) / a[(m+1)*n+m] + a[m*n+m+1];
							q = a[(m+1)*n + m+1] - z - r - s;
							r = a[(m+2)*n+m+1];
							s=Math.abs(p)+Math.abs(q)+Math.abs(r);
							p /= s;
							q /= s;
							r /= s;
							if(m == l) break;
							u = Math.abs(a[m*n+m-1]) * (Math.abs(q) + Math.abs(r));
							v = Math.abs(p) * (Math.abs(a[(m-1) * n + m - 1]) + Math.abs(z) + Math.abs(a[(m+1)*n+m+1]));
							if((float)(u + v) == v) break;
						}
						for(i = m+2; i <= nn; i++) {
							a[i*n+i-2] = 0.0f;
							if(i != m+2) a[i*n+i-3] = 0.0f;
						}
						for(k = m; k <= nn-1; k++) {
							if(k != m) {
								p = a[k*n+k-1];
								q=a[(k+1)*n+k-1];
								r = 0.0f;
								if(k != nn-1) r = a[(k+2)*n+k-1];
								if((x=Math.abs(p) + Math.abs(q) + Math.abs(r)) != 0.0f) {
									p /= x;
									q /= x;
									r /= x;
								}
							}
							if((s = SIGN((float)Math.sqrt(p*p+q*q+r*r), p)) != 0.0f) {
								if(k == m) {
									if(l != m) a[k*n+k-1] = -a[k*n+k-1];
								}
								else a[k*n+k-1] = -s*x;
								p += s;
								x = p/s;
								y = q/s;
								z = r/s;
								q /= p;
								r /= p;
								for(j = k; j <= nn; j++) {
									p = a[k*n+j] + q*a[(k+1)*n + j];
									if(k != nn-1) {
										p += r*a[(k+2) * n+j];
										a[(k+2)*n+j] -= p*z;
									}
									a[(k+1)*n+j] -= p*y;
									a[k*n+j] -= p*x;
								}
								mmin = nn < k+3 ? nn : k+3;
								for(i = l; i <= mmin; i++) {
									p = x * a[i*n+k] + y * a[i*n+k+1];
									if(k != nn-1) {
										p += z * a[i*n+k+2];
										a[i*n+k+2] -= p*r;
									}
									a[i*n+k+1] -= p*q;
									a[i*n+k] -= p;
								}
							}
						}
					}
				}
			} while(l < nn-1);
		}
	}
	
	private static void elmhes(float [] a, int n) {
		int		m, j, i;
		float	y, x, temp;
		
		for(m = 1; m < n-1; m++) {
			x = 0.0f;
			i = m;
			for(j = m; j < n; j++) {
				if(Math.abs(a[j*n+m-1]) > Math.abs(x)) {
					x = a[j*n+m-1];
					i=j;
				}
			}
			if(i != m) {
				for(j = m-1; j < n; j++) { temp = a[i*n+j]; a[i*n+j] = a[m*n+j]; a[m*n+j] = temp; }
				for(j = 0; j < n; j++) { temp = a[j*n+i]; a[j*n+i] = a[j*n+m]; a[j*n+m] = temp; }
			}
			if(x != 0) {
				for(i = m+1; i < n; i++) {
					if((y = a[i*n+m-1]) != 0.0f) {
						y /= x;
						a[i*n+m-1] = y;
						for(j = m; j < n; j++) a[i*n+j] -= y*a[m*n+j];
						for(j = 0; j < n; j++) a[j*n+m] += y * a[j*n+i];
					}
				}
			}
		}
	}
	
	private static void balanc(float [] a, int n) {
		int	last, j, i;
		float	s, r, g, f, c, sqrdx;
		final int		RADIX = 2;
		
		sqrdx = RADIX * RADIX;
		last = 0;
		while(last == 0) {
			last = 1;
			for(i = 0; i < n; i++) {
				r = c = 0.0f;
				for(j = 0; j < n; j++) {
					if(j != i) {
						c += Math.abs(a[j*n+i]);
						r += Math.abs(a[i*n+j]);
					}
				}
				if(c != 0 && r != 0) {
					g = r / RADIX;
					f = 1.0f;
					s = c + r;
					while(c < g) {
						f *= RADIX;
						c *= sqrdx;
					}
					g = r * RADIX;
					while(c > g) {
						f /= RADIX;
						c /= sqrdx;
					}
					if((c + r) / f < 0.95f * s) {
						last = 0;
						g = 1.0f / f;
						for(j = 0; j < n; j++) a[i*n+j] *= g;
						for(j = 0; j < n; j++) a[j*n+i] *= f;
					}
				}
			}
		}				
	}
	
	private static float SIGN(float a, float b) {
		if(b < 0) return -Math.abs(a);
		else return Math.abs(a);
	}
	
	
	public static void sort(Object [] a, Comparer comparer) { sort(a, null, 0, a.length - 1, true, comparer); }

	/** Sort a and b, using a as the reference - stolen from somewhere on internet */
	public static void sort(Object [] a, Object [] b, int from, int to, boolean ascending, Comparer comparer) {
		 // No sort 
		 if (a == null || a.length < 2) return;

		 // sort using Quicksort
		 int i = from, j = to;
		 Object center = a[ (from + to) / 2 ];
		 do {
		   if (ascending) {
		     while( (i < to) && (comparer.compare(  center, a[i]) > 0) ) i++;
		     while( (j > from) && (comparer.compare(center, a[j]) < 0) ) j--;
		     } 
		   else {
		     // Decending sort
		     while( (i < to) && (comparer.compare(  center, a[i]) < 0) ) i++;
		     while( (j > from) && (comparer.compare(center, a[j]) > 0) ) j--;
		     }
		   if (i < j) {
		     // Swap elements
		     Object temp = a[i]; a[i] = a[j]; a[j] = temp;
		     // Swap in b array if needed
		     if (b != null) {
		       temp = b[i]; b[i] = b[j]; b[j] = temp;
		     }
		   }
		   if (i <= j) { i++; j--; }
		   } while(i <= j);
		 // Sort the rest
		 if (from < j) sort(a, b, from, j, ascending, comparer);
		 if (i < to) sort(a, b, i, to, ascending, comparer);
	 }

	public static interface Comparer {
	   /**
	    * The interface implementation should compare the two
	    * objects and return an int using these rules:
	    * if (a > b)  return > 0;
	    * if (a == b) return 0;
	    * if (a < b)  return < 0;
	    */
	    public int compare(Object a, Object b);
    }
}