/*
 *  Draco Econometrics
 *  Copyright (C) 2006-2007 Approximatrix, LLC
 *  http://www.approximatrix.com
 *  
 *  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 3 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, see <http://www.gnu.org/licenses/>.
 */

package com.approximatrix.draco.math;
import org.netlib.blas.*;

/*
 * Created on Feb 7, 2005
 *
 * TODO To change the template for this generated file go to
 * Window - Preferences - Java - Code Style - Code Templates
 */

/**
 * @author jeff
 *
 * TODO To change the template for this generated type comment go to
 * Window - Preferences - Java - Code Style - Code Templates
 */
public class VectorSupport {

	public static double VectorElementSum(double[] vector) {
		double holding = 0.0;
		for(int i=0;i<vector.length;i++) {
			holding = holding + vector[i];
		}
		return holding;
	}
	
	public static double[] VectorOfZeros(int length) {
		double[] holding = new double[length];
		for(int i=0;i<length;i++) {
			holding[i] = 0.0;
		}
		return holding;
	}
	
	public static double[] VectorOfOnes(int length) {
		double[] holding = new double[length];
		for(int i=0;i<length;i++) {
			holding[i] = 1.0;
		}
		return holding;
	}
	
	public static double MaximumValue(double[] vector) {
		double holding = -1.0*Math.abs(vector[IDAMAX.IDAMAX(vector.length,vector,1)-1]);
		for(int i=0;i<vector.length;i++) {
			if(vector[i] > holding) {
				holding = vector[i];
			}
		}
		return holding;		
	}
	
	public static double MinimumValue(double[] vector) {
		double holding = Math.abs(vector[IDAMAX.IDAMAX(vector.length,vector,1)-1]);
		for(int i=0;i<vector.length;i++) {
			if(vector[i] < holding) {
				holding = vector[i];
			}
		}
		return holding;		
	}
	
	public static double[] DivideElements(double[] vector1, double[] vector2) {
		double[] return_val = (double[])vector1.clone();
		if(vector1.length != vector2.length) {
			return null;
		}
		for(int i=0;i<return_val.length;i++) {
			return_val[i] = return_val[i]/vector2[i];
		}
		return return_val;
	}
	
	public static double[] MultiplyElements(double[] vector1, double[] vector2) {
		double[] return_val = (double[])vector1.clone();
		if(vector1.length != vector2.length) {
			return null;
		}
		for(int i=0;i<return_val.length;i++) {
			return_val[i] = return_val[i]*vector2[i];
		}
		return return_val;
	}
	
	public static double[] SqrtElements(double[] vector) {
		double[] return_val = (double[])vector.clone();
		for(int i=0;i<return_val.length;i++) {
			if(return_val[i] > 0.0) {
				return_val[i] = Math.sqrt(return_val[i]);
			} else {
				return_val[i] = -1.0;
				System.err.println("Negative Sqrt in VectorSupport.SqrtElements");
				throw(new java.lang.IllegalArgumentException("Negative Sqrt in VectorSupport.SqrtElements"));
			}
		}
		return return_val;
	}
	
	public static double Magnitude(double[] vector) {
		//double value = DDOT.DDOT(vector.length,vector,1,vector,1);
		//if(value>=0.0)
		//	value = Math.sqrt(value);
		//return value;
		return DNRM2.DNRM2(vector.length,vector,1);
	}
	
	public static double Mean(double[] vector) {
		if(vector == null) return Double.NaN;
		return VectorElementSum(vector)/((double)vector.length);
	}
	
	public static double StandardDeviation(double[] vector) {
		
		if(vector == null) return Double.NaN;
		
		double[] result_vector = (double[])vector.clone();
		double[] mean_vector = VectorOfOnes(vector.length);
		double stddev_denom;
		
		DAXPY.DAXPY(vector.length,-1.0*Mean(vector),
				    mean_vector,1,result_vector,1);
		
		if(result_vector.length == 1) {
			stddev_denom = 1.0;
		} else {
			stddev_denom = Math.sqrt((double)(result_vector.length-1)/1.0);
		}
		
		return new Double(DNRM2.DNRM2(result_vector.length,result_vector,1)/stddev_denom);
	}
	
}
