/*  
 *  src/org/metaqtl/algo/SVDAlgorithm.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.algo;

import org.metaqtl.MultiFitResult;
import org.metaqtl.numrec.NumericalUtilities;
/**
 * This class defines methods to solve usual linear algebra
 * problems by using Singular Value Decomposition (SVD) strategy.
 * 
 * The methods of this class have been adapted from the
 * C implementation available in : 
 * 
 * Numerical Recipes in C, 
 * The Art of Scientific Computing by W. H. Press, B. P. Flannery,
 *  S. A. Teukolsky and W. T. Vetterling.
 * 1990, Cambridge University Press, Cambridge, England
 * ISBN 0-521-35465-X  (the book)
 */
public final strictfp class SVDAlgorithm {
	
	/**
	 * The tolerance for singular value balancing.
	 */
	private static double TOL = 1.0e-5;
	/**
	 * The maximum number of iteration for SVD algorithm.
	 */
	private static final int ITMAX = 30; 
		
	/**
	 * Given a set of data points y[1..m] with indiviudal weights
	 * w[1..m], use chi2 minimization to determine the coefficients a[1..n] of the fitting
	 * function y = X*a. Here we solve the fitting equations using singular value 
	 * decomposition of the X matrix.
	 * @param X the design matrix.
	 * @param y the observations.
	 * @param w the weigths, i.e the invsers of the standard deviations.
	 * @param m the number of observations.
	 * @param n the number of parameters.
	 * @return the parameter estimates.
	 * 
	 * @see MultiFitResult
	 * @see SVDAlgorithm#SVDBackSub(double[][], double[], double[][], double[], int, int)
	 * @see SVDAlgorithm#SVDecomposition(double[][], double[][], double[], int, int)
	 */
	public static MultiFitResult SVDMFit(double[][] X, double[] y, double[] w, int m, int n)
	{
		
		int k,j,i;
		double tmax, tmp, thresh, sum;
		MultiFitResult result = new MultiFitResult(m, n);
		double[][] U  = new double[m][n];
		double[][] V  = new double[n][n];
		double[] t  = new double[n];
		double[] b  = new double[m];
		double[] wti  = new double[m];
		/* Accumulate coefficients of the fitting matrix */
		for(i=0; i< m; i++) {
			for(j=0; j<n; j++) {
				U[i][j]=X[i][j]*w[i];
			}
			b[i]=y[i]*w[i];			
		}
		
		SVDecomposition(U, V, t, m, n);
		
		tmax = 0.0;
		for(j=0; j<n; j++) {
			if (t[j] > tmax) tmax=t[j];
		}
			
		thresh = TOL*tmax;
		
		for(j=0; j<n; j++) 
			if (t[j] < thresh) t[j] = 0.0;
			
		result.theta = SVDBackSub(U, t, V, b, m, n);
		
		/* Evaluate Chi2 */
		result.chi2 = 0.0;
		for(i=0; i<m; i++) {
			for(sum=0.0, j=0; j<n; j++) sum += result.theta[j]*X[i][j];
			tmp=(y[i]-sum)*w[i];
			result.rsd[i] = tmp;
			tmp*=tmp;
			result.chi2 += tmp;
		}
		/* DOF and P-Value */ 
		result.dof    = m-n;
		result.pvalue = NumericalUtilities.gammp(0.5*result.dof, 0.5*result.chi2);
		/* Evaluate cov */
		for(i=0; i<n; i++) {
			wti[i] = 0.0;
			if (t[i] != 0.0) wti[i] = 1.0/(t[i]*t[i]);
		}
		for(i=0; i<n; i++) {
			for(j=0; j<=i; j++) {
				for(sum=0.0, k=0; k<n; k++) sum += V[i][k]*V[j][k]*wti[k];
				result.cov[j][i]=result.cov[i][j]=sum;
			}
		}
		
		result.m = n;
		result.n = m;
		
		return result;
	}
	/**
	 * Solves A.X=B for a vector X, where <code>A</code> is specified by the
	 * matrix <code>U</code>,<code>V</code> and the vector <code>w</code>
	 * as returned by the SVDecomposition() method. No input are destroyed, so
	 * the method may be called sequentially with different <code>b</code> 's.
	 * Note that this routine presumes that you have already zeroed the small
	 * wj's. It does not do this for you. If you haven't zeroed the small wj's,
	 * the this routine is just as ill-conditioned as any direct method, and you
	 * are misusing SVD.
	 * 
	 * @param A the design matrix.
	 * @param w the singular values.
	 * @param V the right singular vectors.
	 * @param b the vector.
	 * @return X the vector solution of A.X = b. 
	 */
	 public static double[] SVDBackSub(										
	 							double[][] U,
								double[]   w,
								double[][] V,
								double[]   b,
								int n,
								int m
							   )
		{
			int jj, j, i;
			double s;
			double[] x   = new double[n];
			double[] tmp = new double[n];
			
			for(j=0; j<m; j++) {
				s = 0.0;
				if (w[j] != 0.0) {
					for(i=0; i< n; i++) s += U[i][j] * b[i];
					s /= w[j];
				}
				tmp[j] = s;
			}
			
			for(j=0; j<m; j++) {
				s = 0.;
				for(jj=0; jj<m;jj++) 
					s += V[j][jj]*tmp[jj];
				x[j]=s;		
			}
			tmp = null;
			
			return x;
		}
		/**
		 * Given a matrix <code>A</code>, this routine computes its singular
		 * value decomposition A = U.W.V' . The matrix U replaces A on output.
		 * The diagonal matrix of singular values W is output as a vector
		 * <code>w</code>. The matrix V (not the transpose V') is output as
		 * the matrix <code>V</code>. The original routine comes from :
		 * Forsythe, G.E., Malcolm, M.A. and Moler, C.B. 1977, Computer Methods
		 * for Mathematical Computations (Englewood Cliffs, NJ: Prentice-Hall),
		 * Chapter 9. which is in turned based on the routine from : Golub, G.H.
		 * and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore:
		 * Johns Hopkins University Press), p8.3 and chapter 12.
		 * 
		 * @param A the original matrix, after the left singular vectors.
		 * @param V the right singular vectors.
		 * @param w the singular values.
		 */
		public static void SVDecomposition(
																	double[][] A,
																	double[][] V,
																	double[]   w,
																	int m,
																	int n
															   )
		{
					int flag, i, its, j, jj, k, l, nm;
					double anorm, c, f, g, h, s, scale, x, y, z;
					double[] rv1 = new double[m];
					
					g=scale=anorm=0.0;
					l = nm = 0;
					
					for(i=0; i<n; i++) {
						l=i+1;
						rv1[i] = scale*g;
						g=s=scale=0.0;
						if (i<m) {
							for(k=i;k<m;k++) scale += Math.abs(A[k][i]);
							if (scale != 0.) {
								for(k=i; k<m; k++) {
									A[k][i] /= scale;
									s += A[k][i]*A[k][i];
								}
								f=A[i][i];
								g = -NumericalUtilities.SIGN(Math.sqrt(s), f);
								h=f*g-s;
								A[i][i] =  f-g;
								for(j=l; j<n; j++) {
									for(s=0., k=i; k<m; k++) {
										s += A[k][i]*A[k][j];
										}									 
									f=s/h;
									for(k=i; k<m; k++) A[k][j] += f*A[k][i];
								}
								for(k=i; k<m; k++) A[k][i] = A[k][i]*scale;
							}
						}
						w[i] = scale*g;
						g=s=scale=0.0;
						if (i < m && i != n-1) {
							for(k=l; k<n; k++) scale += Math.abs(A[i][k]);
							if (scale != 0.) {
								for(k=l; k<n; k++) {
									A[i][k] = A[i][k]/scale;
									s += A[i][k]*A[i][k]; 
								}
								f=A[i][l];
								g = -NumericalUtilities.SIGN(Math.sqrt(s), f);
								h=f*g-s;
								A[i][l] =  f-g;
								for(k=l; k<n; k++) rv1[k] = A[i][k]/h;
								for(j=l; j<m; j++) {
									for(s=0., k=l; k<n; k++) s+= A[j][k]*A[i][k];
									for(k=l; k<n; k++) A[j][k] = A[j][k] + s*rv1[k];
								}
								for(k=l; k<n; k++) A[i][k] = A[i][k]*scale;
							}
						}
						anorm = NumericalUtilities.DMAX(anorm, (Math.abs(w[i]) + Math.abs(rv1[i])));
					}
					for(i=n-1;i>=0; i--) {
						/* Accumulation of right-hand transformations. */
						if (i < n-1) {
							if (g != 0.) {
								for(j=l;j<n; j++)
									V[j][i] = (A[i][j]/A[i][l])/g;
								for(j=l; j<n; j++) {
									for(s=0., k=l; k<n; k++) s += A[i][k]*V[k][j];
									for(k=l; k<n; k++) V[k][j] += V[k][i]*s;
								}
							}
							for(j=l; j<n; j++)	V[i][j] = V[j][i] = 0.0;
						}					
						V[i][i] = 1.0;
						g=rv1[i];
						l=i;
					}
					for(i=NumericalUtilities.IMIN(m-1, n-1); i>=0; i--) {
						l=i+1;
						g=w[i];
						for(j=l; j<n; j++)	A[i][j] = 0.0;
						if (g != 0.0) {
							g = 1.0/g;
							for(j=l; j<n; j++) {
								for(s=0.0, k=l; k<m; k++) s += A[k][i]*A[k][j];
								f = (s/A[i][i])*g;
								for(k=i; k<m; k++) A[k][j] += f*A[k][i];
							}
							for(j=i; j<m; j++) A[j][i] *= g;
						} else for(j=i; j<m; j++) A[j][i] = 0.0;
						++A[i][i];
					}
					/*
					 * Diagonalization of the bidiagonal form : Loop over
					 * singular values and over allowed values
					 */
					 for(k=n-1;k>=0; k--) {
					 	for(its=1; its<=ITMAX; its++) {
					 		flag = 1;
					 		for(l=k; l>=0; l--) {
					 			nm = l-1;
					 			if (Math.abs(rv1[l])+anorm == anorm) {
					 				flag = 0;
					 				break;
					 			}
					 			if (Math.abs(w[nm])+anorm == anorm) break;				 								 	 
					 		}
					 		if (flag != 0) {
					 			/* Cancelation of rv1[l], if l > 1 */
					 			c=0.0;
					 			s=1.0;
					 			for(i=l; i<=k; i++) {
					 				f = s*rv1[i];
					 				rv1[i] = c*rv1[i];
					 				if (Math.abs(f)+anorm == anorm) break;
					 				g=w[i];
					 				h=NumericalUtilities.PYTHAG(f,g);
					 				w[i] = h;
					 				h = 1.0/h;
					 				c=g*h;
					 				s=-f*h;
					 				for(j=0;j<m;j++) {
					 					y=A[j][nm];
					 					z=A[j][i];
					 					A[j][nm]=y*c+z*s;
										A[j][i]=z*c-y*s;
					 				}
					 			}				 			
					 		}
					 		z=w[k];
					 		if (l == k) {
					 			if (z < 0.0) {
					 				w[k] = -z;
					 				for(j=0; j<n; j++) V[j][k]=-V[j][k];
					 			}
					 			break;
					 		}
					 		if (its == ITMAX) {
					 			System.err.println("SVD : Too many iterations\n");
					 			return;
					 		}
					 		x=w[l];
					 		nm=k-1;
					 		y=w[nm];
					 		g=rv1[nm];
					 		h=rv1[k];
					 		f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
					 		g=NumericalUtilities.PYTHAG(f,1.0);
					 		f=((x-z)*(x+z)+h*((y/(f+NumericalUtilities.SIGN(g,f)))-h))/x;
					 		c=s=1.0;
					 		/* Next QR transformation */
					 		for(j=l; j<=nm; j++) {
					 			i=j+1;
					 			g=rv1[i];
					 			y=w[i];
					 			h=s*g;
					 			g=c*g;
					 			z=NumericalUtilities.PYTHAG(f,h);
					 			rv1[j]=z;
					 			c=f/z;
					 			s=h/z;
					 			f=x*c+g*s;
					 			g=g*c-x*s;
					 			h=y*s;
					 			y *= c;
					 			for(jj=0; jj<n; jj++) {
					 				x=V[jj][j];
					 				z=V[jj][i];
									V[jj][j] = x*c + z*s;
									V[jj][i] = z*c - x*s;
					 			}
					 			z = NumericalUtilities.PYTHAG(f,h);
					 			w[j]=z;
					 			if (z != 0.0) {
					 				z=1.0/z;
					 				c=f*z;
					 				s=h*z;
					 			}
					 			f=c*g+s*y;
					 			x=c*y-s*g;
					 			for(jj=0; jj<m; jj++) {
					 				y=A[jj][j];
					 				z=A[jj][i];
									A[jj][j] = y*c+z*s;
									A[jj][i] = z*c-y*s;
					 			}
					 		}
					 		rv1[l]=0.0;
					 		rv1[k]=f;
					 		w[k]=x;
					 	}
					 }
					rv1=null;
				}
}
