/*  
 *  src/org/metaqtl/numrec/NumericalUtilities.java
 * 
 *  $Author: veyrieras $, $Date: 2006-10-17 15:41:55 $, $Version$
 *
 *  MetaQTL : a Java package for QTL Meta-Analysis
 * 
 *  Copyright (C) 2003-2006  Jean-Baptiste Veyrieras, INRA, France.
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA * 
 */
package org.metaqtl.numrec;

/**
 * 
 * Class Description Here
 * 
 * @author Jean-Baptiste Veyrieras
 *  
 */
public final strictfp class NumericalUtilities {
	
	public static double[] GAMCOF = {76.18009173, -86.50532033, 24.01409822,
			  -1.231739516, 0.120858003e-2, -0.536382e-5};
	
	private static final int ITMAX = 100;
	
	private static final double EPS = 3.e-8;

	private static final int NTAB=10;

	private static final double BIG = 1.e+30;

	private static final double CON = 1.4;
	
	private static final double CON2 = CON*CON;
	
	private static final double SAFE = 2.0;

	/**
	 * Returns abs(a) if b > 0, -abs(a) otherwise.
	 * 
	 * @param a
	 * @param b
	 * @return +/- abs(a)
	 */	
	public static double SIGN(double a, double b) {
		return (b > 0.0) ? java.lang.Math.abs(a) : -java.lang.Math.abs(a);		
	}

	/**
	 * Returns the max between <code>a</code> and <code>b</code>
	 * 
	 * @param a
	 * @param b
	 * @return max(a,b)
	 */
	public static double DMAX(double a, double b) {
		return (a > b) ? a : b;
	}
	/**
	 * Returns the min between <code>a</code> and <code>b</code>
	 * 
	 * @param a
	 * @param b
	 * @return min(a,b)
	 */
	public static int IMIN(int a, int b) {
		return (a < b) ? a : b;
	}
	/**
	 * Computes sqrt(a^2+b^2)
	 * 
	 * @param a
	 * @param b
	 * @return sqrt(a^2+b^2) the hypothenus
	 */
	public static double PYTHAG(double a, double b) {
		double absa, absb;
		absa = Math.abs(a);
		absb = Math.abs(b);
		if (absa > absb) return absa*Math.sqrt(1.0 + SQR(absb/absa));
		else return (absb == 0.0) ? 0.0 : absb*Math.sqrt(1.0 + SQR(absa/absb));
	}
	/**
	 * Returns a^2
	 * 
	 * @param a
	 * @return a^2
	 */
	public static double SQR(double a) {
		return (a == 0.0) ? 0.0 : a*a;
	}
	/**
	 * This routine returns the value at x of the normal
	 * centered gaussian density.
	 */
	public static double gauss(double x) {
		return Math.exp(-0.5*x*x)/Math.sqrt(2.0*Math.PI);
	}
	/**
	 * @param sd
	 * @return
	 */
	public static double sd2ci(double sd) {
		return sd*2*1.96;
	}
	/**
	 * 
	 * @param ci
	 * @return
	 */
	public static double ci2sd(double ci) {
		return ci/(2*1.96);
	}
	
	
	public static double gammln(double xx)
	{
	  int j;
	  double x, tmp, ser;
	  
	  x = xx - 1.0;
	  tmp = x + 5.5;
	  tmp -= (x + 0.5) * Math.log(tmp);
	  ser = 1.0;
	  for (j = 0; j <= 5; j++) {
	    x += 1.0;
	    ser += GAMCOF[j] / x;
	  }
	  return -tmp + Math.log(2.50662827465 * ser);
	}
	
	public static double gammp(double a, double  x)
	{
	  double y[];

	  if (x < 0.0 || a <= 0.0)
	    return 0.0;
	  if (x < (a + 1.0)) {
	    y = gser(a, x);
	    return y[0];
	  }
	  else {
	    y = gcf(a, x);
	    return 1.0 - y[0];
	  }
	  
	}

	public static  double[] gcf(double  a,double  x)
	{
	  int n;
	  double gold = 0.0, g, fac = 1.0, b1 = 1.0;
	  double b0 = 0.0, anf, ana, an, a1, a0 = 1.0;
	  double[] ret = new double[2];

	  ret[1] = gammln(a);
	  
	  a1 = x;
	  for (n = 1; n <= ITMAX; n++) {
	    an = (double) n;
	    ana = an - a;
	    a0 = (a1 + a0 * ana) * fac;
	    b0 = (b1 + b0 * ana) * fac;
	    anf = an * fac;
	    a1 = x * a0 + anf * a1;
	    b1 = x * b0 + anf * b1;
	    if (a1!=0.0) {
	      fac = 1.0 / a1;
	      g = b1 * fac;
	      if (Math.abs((g - gold) / g) < EPS) {
		ret[0] = Math.exp(-x + a * Math.log(x) - ret[1]) * g;
		return ret;
	      }
	      gold = g;
	    }
	  }
	  return ret;
	}



	public static  double[] gser(double a,double x)
	{
	  int n;
	  double sum, del, ap;
	  double[] ret = new double[2];

	  ret[1] = gammln(a);
	  
	  if (x <= 0.0) {
	    if (x < 0.0)
	     ret[0] = 0.0;
	    return ret;
	  }
	  else {
	    ap = a;
	    del = sum = 1.0 / a;
	    for (n = 1; n <= ITMAX; n++) {
	      ap += 1.0;
	      del *= x / ap;
	      sum += del;
	      if (Math.abs(del) < Math.abs(sum) * EPS) {
			ret[0] = sum * Math.exp(-x + a * Math.log(x) - ret[1]);
			return ret;
	      }
	    }
	    return ret;
	  }
	}
	
    /**
	 * Used by betai : Evaluates continued fraction for incomplete beta function
	 * by modified Lentz's method.
	 * 
	 * @param a
	 * @param b
	 * @param x
	 * @return
	 */
	public static double betacf(double a,double  b,double  x)
	{
	  double qap, qam, qab, em, tem, d;
	  double bz, bm = 1.0, bp, bpp;
	  double az = 1.0, am = 1.0, ap, app, aold;
	  int m;

	  qab = a + b;
	  qap = a + 1.0;
	  qam = a - 1.0;
	  bz = 1.0 - qab * x / qap;
	  for (m = 1; m <= ITMAX; m++) {
	    em = (float) m;
	    tem = em + em;
	    d = em * (b - em) * x / ((qam + tem) * (a + tem));
	    ap = az + d * am;
	    bp = bz + d * bm;
	    d = -(a + em) * (qab + em) * x / ((qap + tem) * (a + tem));
	    app = ap + d * az;
	    bpp = bp + d * bz;
	    aold = az;
	    am = ap / bpp;
	    bm = bp / bpp;
	    az = app / bpp;
	    bz = 1.0;
	    if (Math.abs(az - aold) < (EPS * Math.abs(az)))
	      return(az);
	  }
	  return(0.0);
	}
	/**
	 * Returns the incomplete beta function Ix(a,b)
	 * 
	 * @param a
	 * @param b
	 * @param x
	 * @return
	 */
	public static  double betai(double a,double  b,double  x)
	{
	  double bt;

	  if (x < 0.0 || x > 1.0)
	    return 0.0; 
	  if (x == 0.0 || x == 1.0)
	    bt = 0.0;
	  else
	    bt = Math.exp(gammln(a + b) - gammln(a) - gammln(b) + a * Math.log(x) + b * Math.log(1.0 - x));
	  if (x < (a + 1.0) / (a + b + 2.0))
	    return bt * betacf(a, b, x) / a;
	  else
	    return 1.0 - bt * betacf(b, a, 1.0 - x) / b;
	}
	/**
	 * Returns the beta function I(a,b)
	 * 
	 * @param z
	 * @param w
	 * @return
	 */
	public static double beta(double z, double  w)
	{
	  return  Math.exp(gammln(z) + gammln(w) - gammln(z + w)) ;
	}
	
	/**
	 * Sort the given array of double arrin[1..n] without destroying it and
	 * returns the sorted indices of the array elements into indx[1..n].
	 * 
	 * @param n
	 *            the array size
	 * @param arrin
	 * @param indx
	 */
	public static void indexx(int n, double[] arrin, int[] indx)
	{
		  int l, j, ir, indxt, i;
		  double q;

		  for (j=1;j<=n;j++) indx[j] = j;
		  l = (n >> 1) + 1;
		  
		  ir = n;
		  
		  for (;;) {
		    if (l > 1)
		      q = arrin[(indxt = indx[--l])];
		    else {
		      q = arrin[(indxt = indx[ir])];
		      indx[ir] = indx[1];
		      if (--ir == 1) {
			indx[1] = indxt;
			return;
		      }
		    }
		    i = l;
		    j = l << 1;
		    while (j <= ir) {
		      if (j < ir && arrin[indx[j]] < arrin[indx[j + 1]])
			j++;
		      if (q < arrin[indx[j]]) {
			indx[i] = indx[j];
			j += (i = j);
		      }
		      else
			j = ir + 1;
		    }
		    indx[i] = indxt;
		  }
		  		  
	}
	/**
	 * 
	 * @param function
	 * @param x1
	 * @param x2
	 * @param xacc
	 * @return
	 */
	public static double rtsafe(NRFunction funcd, double x1, double x2, double xacc){
		int j;
		double df,dx,dxold,f,fh,f1;
		double temp,xh,xl,rts;
		
		funcd.funcd(x1);
		f1 = funcd.getF();
		funcd.funcd(x2);
		fh = funcd.getDF();
		
		if ((f1 > 0.0 && fh > 0.0) || (f1 < 0.0 && fh < 0.0)) {
			System.err.println("Root must be bracketed in rtsafe");
			return 0.0;
		}
		
		if (f1 == 0.0) return x1;
		if (fh == 0.0) return x2;
		
		if (f1 < 0.0) {
			xl=x1;
			xh=x2;
		} else {
			xh=x1;
			xl=x2;
		}
		
		rts = 0.5*(x1+x2);
		dxold=Math.abs(x2-x1);
		dx=dxold;
		
		funcd.funcd(rts);
		f  = funcd.getF();
		df = funcd.getDF();
		
		for(j=1;j<=ITMAX;j++) {
			if ((((rts-xh)*df-f)*((rts-x1)*df-f) > 0.0) || (Math.abs(2.0*f) > Math.abs(dxold*df))) {
				dxold = dx;
				dx = 0.5*(xh-xl);
				rts=xl+dx;
				if (xl == rts) return rts;
			} else {
				dxold=dx;
				dx=f/df;
				temp=rts;
				rts -= dx;
				if (temp == rts) return rts;
			}
			if (Math.abs(dx) < xacc) return rts;
			funcd.funcd(rts);
			f  = funcd.getF();
			df = funcd.getDF();
			if (f < 0.0) xl = rts;
			else xh = rts;
		} 
		
		System.err.println("Maximum number of iterations exceeded in rtsafe");
		
		return 0.0;
	}
	/**
	 * Returns the derivative of a function func at a point x by Ridders' method of
	 * polynomial extrapolation. The value h is input as an extimated initial stepsize;
	 * it need not to be small, but rather should an increment in x over which func changes
	 * substantially. The error is returned and the derivative is set in function f.
	 * 
	 */
	public static double dfridr(Function f, double x, double h) {
		int i,j;
		double err,errt,fac,hh,ans;
		double[][] a;
		
		if (h == 0.0) return 0.0;
		
		a = new double[NTAB+1][NTAB+1];
		hh = h;
		ans = a[1][1] = (f.func(x+hh)-f.func(x-hh))/(2.0*hh);
		err     = BIG;
		
		for(i=2;i<=NTAB;i++) {
			hh /= CON;
			a[1][i] = (f.func(x+hh)-f.func(x-hh))/(2.0*hh);
			fac=CON2;
			for(j=2;j<=i;j++) {
				a[j][i]=(a[j-1][i]*fac-a[j-1][i-1])/(fac-1.0);
				fac=CON2*fac;
				errt=Math.max(Math.abs(a[j][i]-a[j-1][i]), Math.abs(a[j][i]-a[j-1][i-1]));
				if (errt <= err) {
					err = errt;
					ans = a[j][i];
				}
			}
			if (Math.abs(a[i][i]-a[i-1][i-1]) >= SAFE*err) break;
		}
		
		f.setDF(ans);
		
		return err;
	}
	
	public interface Function {
		/**
		 * Returns the value of the function in x.
		 * @param x
		 * @return f(x)
		 */
		public double func(double x);
		/**
		 * Returns the derivative of the function.
		 * in x.
		 * @return
		 */
		public double getDF();
		/**
		 * Sets the derivative of the function at x.
		 * @return
		 */
		public void setDF(double x);
		
	}
	
	/**
	 * 
	 * @author jveyrier
	 *
	 * One dimensional function definition for Newton Raphson root
	 * finding algorithm.
	 */
	public interface NRFunction {
		/**
		 * Compute both the value and the derivate of the
		 * function in x.
		 * @param x
		 */
		public void funcd(double x);
		/**
		 * Returns the value of the function in x.
		 * @return
		 */
		public double getF();
		/**
		 * Returns the value of the derivate in x;
		 * @return
		 */
		public double getDF();
	}

	/**
	 * Swaps elements i and j in a.
	 * @param a
	 * @param i
	 * @param j
	 */
	public static void ISWAP(int[] a, int i, int j) {
		int tmp = a[i];
		a[i]=a[j];
		a[j]=tmp;
	}
	/**
	 * Swaps elements i and j in a.
	 * @param a
	 * @param i
	 * @param j
	 */
	public static void DSWAP(double[] a, int i, int j) {
		double tmp = a[i];
		a[i]=a[j];
		a[j]=tmp;
	}
}
