using System;
using System.Threading;
using System.Collections;
using System.Linq;
using System.Collections.Generic;
using MarketData.MarketDataModel;
using MarketData.Utils;

// Filename: Numerics.cs
// Author:Sean Kessler / Werner Stanzl
// Date:07/2003

namespace MarketData.Numerical
{
  public class Numerics
  {
    public enum CorrelationType{PerfectDownhill,StrongDownhill,ModerateDownhill,WeakDownhill,NoLinearRelationship,WeakUphill,ModerateUphill,StrongUphill,PerfectUphill};

    /// <summary>
    /// Discount - Discounts the given number by a percentage.  For example to discount $600.00 by 20% call Discount(600.00,.20)
    /// </summary>
    /// <param name="number"></param>
    /// <param name="discountPercentDecimal"></param>
    /// <returns></returns>
    public static double Discount(double number, double discountPercentDecimal)
    {
      return number * (1.00 - discountPercentDecimal);
    }

// (i.e.) 2.789  will round to 2.79
    public static double Round(double x)
		{
		  return Math.Round(x, 2, MidpointRounding.AwayFromZero);
		}
    public static void ZeroForNaNOrInfinity(ref float[] floatArray)
    {
      for(int index=0;index<floatArray.Length;index++)
			{
			  if(float.IsNaN(floatArray[index]))floatArray[index]=0.0f;
        else if(float.IsInfinity(floatArray[index]))floatArray[index]=0.0f;
			}
    }
    public static double[] ToDouble(float[] floatArray)
    {
      double[] doubleArray=new double[floatArray.Length];
      for(int index=0;index<floatArray.Length;index++)doubleArray[index]=floatArray[index];
      return doubleArray;
    }
    public static float[] ToFloat(double[] doubleArray)
    {
      float[] floatArray=new float[doubleArray.Length];
      for(int index=0;index<doubleArray.Length;index++)floatArray[index]=(float)doubleArray[index];
      return floatArray;
    }
    public static float[] Reverse(ref float[] x)
    {
      float[] y=new float[x.Length];
      for(int index=0;index<x.Length;index++)
      {
        y[index]=x[x.Length-index-1];
      }
      return y;
    }
    public static double[] Reverse(ref double[] x)
    {
      double[] y=new double[x.Length];
      for(int index=0;index<x.Length;index++)
      {
        y[index]=x[x.Length-index-1];
      }
      return y;
    }
// Removes the outliers from the observations
    public static double[] RemoveOutliers(double[] observations,int deviations=2)
    {
      double stddev=StdDev(ref observations);
      deviations=Math.Abs(deviations);
      if(deviations<1)deviations=1;
      double[] newObservations=(from double value in observations where value<deviations*stddev && value > (deviations*-1)*stddev select value).ToArray<double>();
      return newObservations;
    }
    public static LeastSquaresResult LeastSquares(double[] observations)
    {
      return LeastSquaresHelper.CalculateLeastSquares(observations);
    }
    public static LeastSquaresResultWithR2 LeastSquaresWithR2(double[] observations)
    {
      return LeastSquaresHelper.CalculateLeastSquaresWithR2(observations);
    }
    public static double RSquared(double[] observations)
    {
      try
      {
        double[] estimated = new double[observations.Length];
        double[] estimatedLessMeanSquared = new double[observations.Length];
        double sumOfEstimatedLessMeanSquared = 0.00;
        double rSquared = 0.00;
        LeastSquaresResult leastSquaresResult = null;
        double mean = Numerics.Mean(ref observations);
        double sumofsquares = 0.00;
        foreach (double observation in observations) sumofsquares += Math.Pow(observation - mean, 2);
        leastSquaresResult = LeastSquares(observations);
        double slope = leastSquaresResult.Slope * -1.00;  // LeastSquares generator reverses the slope.  Because prices fed to it are in reverse order
        double yIntercept = leastSquaresResult.YIntercept;
        for (int index = 0; index < observations.Length; index++) estimated[index] = yIntercept + (slope * (index + 1));
        for (int index = 0; index < observations.Length; index++) estimatedLessMeanSquared[index] = Math.Pow(estimated[index] - mean, 2);
        foreach (double value in estimatedLessMeanSquared) sumOfEstimatedLessMeanSquared += value;
        rSquared = sumOfEstimatedLessMeanSquared / sumofsquares;
        return rSquared;
      }
      catch (Exception)
      {
        return double.NaN;
      }
    }
/// <summary>YIntercept - Calculates slope of values.</summary>
/// <returns>YIntercept</returns>
    public static double YIntercept(double xMean,double yMean,double slope)
    {
      return yMean-(slope*xMean);
    }
    public static double Slope(ref double[] x,ref double[] y)
    {
      double sumxy=0.00;
      double sumx=0.00;
      double sumy=0.00;
      double sumx2=0.00;

      if(x.Length!=y.Length)return double.NaN;
      for(int index=0;index<x.Length;index++)
      {
        sumxy+=(x[index]*y[index]);
        sumx+=x[index];
        sumy+=y[index];
        sumx2+=(x[index]*x[index]);
      }
      return ((sumxy-sumx*sumy/x.Length)/(sumx2-sumx*sumx/x.Length));
    }
/// <summary>Slope - Calculates slope of values.</summary>
/// <param name="y">y-values x-values are assumed to be sequential</param>
/// <returns>Slope</returns>
    public static float Slope(float[] y)
    {
      float sum1 = 0.00F;
      float sum2 = 0.00F;
      float[] x=new float[y.Length];
      for (int index = 0; index < y.Length; index++) x[index] = index + 1;
      float averageY = Mean(ref y);  
      float averageX = Mean(ref x);
      for (int index = 0; index < y.Length; index++)
      {
        float currX = x[index];
        sum1+=(currX-averageX)*(y[index]-averageY);
        sum2 += (float)Math.Pow(currX-averageX,2); 
      }
      if (0 == sum2) return float.NaN;
      return sum1 / sum2;
    }
/// <summary>Slope - Calculates slope of values.</summary>
/// <param name="y">y-values x-values are assumed to be sequential</param>
/// <returns>Slope</returns>
    public static double Slope(ref double[] y)
    {
      double sum1 = 0.00F;
      double sum2 = 0.00F;
      double[] x=new double[y.Length];
      for (int index = 0; index < y.Length; index++) x[index] = index + 1;
      double averageY = Mean(ref y);  
      double averageX = Mean(ref x);
      for (int index = 0; index < y.Length; index++)
      {
        double currX = x[index];
        sum1+=(currX-averageX)*(y[index]-averageY);
        sum2 += (float)Math.Pow(currX-averageX,2); 
      }
      if (0 == sum2) return float.NaN;
      return sum1 / sum2;
    }
    /// <summary>Slope - Calculates slope of values.</summary>
    /// <param name="y">y-values x-values are assumed to be sequential</param>
    /// <returns>Slope</returns>
    public static double Slope(double[] y)
    {
      double sum1 = 0.00F;
      double sum2 = 0.00F;
      double[] x = new double[y.Length];
      for (int index = 0; index < y.Length; index++) x[index] = index + 1;
      double averageY = Mean(ref y);
      double averageX = Mean(ref x);
      for (int index = 0; index < y.Length; index++)
      {
        double currX = x[index];
        sum1 += (currX - averageX) * (y[index] - averageY);
        sum2 += (float)Math.Pow(currX - averageX, 2);
      }
      if (0 == sum2) return float.NaN;
      return sum1 / sum2;
    }
/// <summary>Median - Calculates median of values.</summary>
/// <param name="values">values</param>
/// <returns>Median</returns>
  	public static double Median(ref double[] values)
  	{
  	  int length=values.Length;
  	  if (length == 0)throw new Exception("Numerics::cannot calculate median of an empty array");
  	  if (length == 1 )return values[0];
      double [] sorted = new double[length];
   	  Array.Copy(values, 0, sorted, 0, length);
  	  Array.Sort(sorted);
  	  if( (length%2) == 0 ) return (sorted[length/2]+sorted[length/2-1])/2;
  	  return sorted[(length-1)/2];
  	}
/// <summary>Floor - Calculates floor of values.</summary>
/// <param name="values">values</param>
/// <returns>floor</returns>
    public static void Floor(ref double[] values)
    {
	    for (int i=0; i< values.Length; i++)
	    {
	      values[i] = Math.Floor(values[i]);
      }
  	}
/// <summary>Percentile - Calculates percentile of values.</summary>
/// <param name="values">values</param>
/// <returns>Percentile</returns>
    public static float Percentile(int perc,ref float[] values)
    {
      if(values.Length<5)return float.NaN;
      int index;
      Array.Sort(values);
      index=(int)Math.Floor(perc*(values.Length)/100.00F)-1;
      if(index<0)index=0;
      else if(index>=values.Length)index=values.Length-1;
      return values[index];
    }
/// <summary>Min - Calculates minimum of values.</summary>
/// <param name="values">values</param>
/// <returns>minimum</returns>
    public static float Min(ref float[] values)
    {
      if(0==values.Length)return 0;
      float min=values[0];
      for(int index=1;index<values.Length;index++)
      {
        if(values[index]<min)min=values[index];
      }
      return min;
    }
/// <summary>Min - Calculates minimum of values.</summary>
/// <param name="values">values</param>
/// <returns>minimum</returns>
    public static double Min(ref double[] values)
    {
      if(0==values.Length)return 0;
      double min=values[0];
      for(int index=1;index<values.Length;index++)
      {
        if(values[index]<min)min=values[index];
      }
      return min;
    }
/// <summary>Max - Calculates maximum of values.</summary>
/// <param name="values">values</param>
/// <returns>maximum</returns>
    public static float Max(ref float[] values)
    {
      if(0==values.Length)return 0;
      float max=values[0];
      for(int index=1;index<values.Length;index++)
      {
        if(values[index]>max)max=values[index];
      }
      return max;
    }
/// <summary>Max - Calculates maximum of values.</summary>
/// <param name="values">values</param>
/// <returns>maximum</returns>
    public static double Max(ref double[] values)
    {
      double max=values[0];
      for(int index=1;index<values.Length;index++)
      {
        if(values[index]>max)max=values[index];
      }
      return max;
    }
/// <summary>Mean - Calculates mean of values.</summary>
/// <param name="values">values</param>
/// <returns>Mean</returns>
    public static float Mean(ref float[] values)
    {
      if (null == values) return float.NaN;
      int length=values.Length;
      float mean=0;
      if(0==length)return float.NaN;
      for(int index=0;index<length;index++)
      {
        mean+=values[index];
      }
      return mean/length;
    }
/// <summary>Mean - Calculates mean of values.</summary>
/// <param name="values">values</param>
/// <returns>Mean</returns>
    public static double Mean(ref double[] values)
    {
      if (null == values) return float.NaN;
      int length=values.Length;
      double mean=0;
      if(0==values.Length)return double.NaN;
      for(int index=0;index<length;index++)
      {
        mean+=values[index];
      }
      return mean/length;
    }
/// <summary>MeanAbs - Calculates mean of abs(values).</summary>
/// <param name="values">values</param>
/// <returns>Mean</returns>
    public static float MeanAbs(ref float[] values)
    {
      int length=values.Length;
      float mean=0;
      if(0==length)return mean;
      for(int index=0;index<length;index++)
      {
        mean+=Math.Abs(values[index]);
      }
      return mean/length;
    }
/// <summary>GeometricMean - Calculates mean of values.</summary>
/// <param name="values">values</param>
/// <returns>Mean</returns>
    public static float GeometricMean(ref float[] values)
    {
      float geometricMean=1.00F;
      int length=values.Length;
      if(0==length)return geometricMean;
      for(int index=0;index<length;index++)
      {
        if(0.00==values[index])geometricMean*=1.00F;
        else geometricMean*=values[index];
      }
      bool isNegative=geometricMean<0.00?true:false;
      double result=Math.Pow(isNegative?geometricMean*-1.00:geometricMean,1.00/(double)length);
      if(isNegative)result*=-1.00;
      return (float)result;
    }
/// <summary>Sum - Calculates Sum.</summary>
/// <param name="x">float array x.</param>
/// <returns>float Sum</returns>
    public static float Sum(ref float[] x)
    {
      float sum=0;
      int length=x.Length;
      if(0==length)return sum;
      for(int index=0;index<length;index++)sum+=x[index];
      return sum;
    }
	  /// <summary>Sum - Calculates Sum.</summary>
	  /// <param name="x">double array x.</param>
	  /// <returns>double Sum</returns>
	  public static double Sum(ref double[] x)
	  {
		  double sum=0;
		  int length=x.Length;
		  if(0==length)return sum;
		  for(int index=0;index<length;index++)sum+=x[index];
		  return sum;
	  }

/// <summary>AnnualizedVolatility - Calculates annualized volatility of values.  can be any number of period</summary>
/// <param name="values">values</param>
/// <returns>volatility</returns>
    public static float AnnualizedVolatility(ref float[] values)
    {
      return (float)(Math.Sqrt(252.00)*StdDev(ref values));  // 252 is average number of trading days in a year
    }
/// <summary>AnnualizedVolatility - Calculates annualized volatility of values.  can be any number of periods</summary>
/// <param name="values">values</param>
/// <returns>volatility</returns>
    public static double AnnualizedVolatility(ref double[] values)
    {
      return Math.Sqrt(252)*StdDev(ref values);  // 252 is average number of trading days in a year
    }
/// <summary>Volatility - Calculates volatility of values.</summary>
/// <param name="values">values</param>
/// <returns>volatility</returns>
    public static float Volatility(ref float[] values)
    {
      return StdDev(ref values);
    }
/// <summary>Volatility - Calculates volatility of values.</summary>
/// <param name="values">values</param>
/// <returns>volatility</returns>
    public static double Volatility(ref double[] values)
    {
      return StdDev(ref values);
    }
/// <summary>Standardize - Standardizes the input.</summary>
/// <param name="values">values</param>
/// <returns>volatility</returns>
    public static float[] Standardize(ref float[] values)
    {
      float[] returnData=new float[values.Length];
      float valuesMean=Mean(ref values);
      float valuesStd=Volatility(ref values);

      if(0==valuesStd)return values;
      for(int index=0;index<values.Length;index++)
      {
        returnData[index]=(values[index]-valuesMean)/valuesStd;
      }
      return returnData;
    }
/// <summary>DownStd - Returns downside standard deviation of values.</summary>
/// <param name="values">values</param>
/// <returns>downside standard deviation</returns>
    public static float DownStd(ref float[] values)
    {
      float average;
      float downStd=0;
      if(1>=values.Length)return downStd;
      average=Mean(ref values);
      foreach(float value in values)
      {
        if(value<average)downStd+=(float)Math.Pow(value-average,2);
      }
      return (float)Math.Sqrt(downStd/(values.Length-1));
    }
/// <summary>BattingAverage - Returns batting average of values.</summary>
/// <param name="values">values</param>
/// <returns>batting average</returns>
    public static float BattingAverage(ref float[] values)
    {
      float battingAverage=0;
      if(0==values.Length)return float.NaN;
      foreach(float value in values)
      {
        if(value>0)battingAverage++;
      }
      return battingAverage/values.Length;
    }
/// <summary>AverageReturnWithSpline -This is a collection of prices ordered highest date first.</summary>
/// <param name="values">values</param>
/// <returns>AverageReturnWithSpline</returns>
    public static double AverageReturnWithSpline(TimeSeriesCollection timeSeriesCollection)
    {
      List<Element> sourceElements = new List<Element>();
      List<Element> destElements = null;
      List<float> values = new List<float>();

      if (null == timeSeriesCollection || 0 == timeSeriesCollection.Count) return double.NaN;
      for (int index = timeSeriesCollection.Count-1; index >=0; index--)
      {
        TimeSeriesElement timeSeriesElement=timeSeriesCollection[index];
        sourceElements.Add(new Element(Utility.DateToLong(timeSeriesElement.AsOf), timeSeriesElement.Value));
      }
      DateTime startDate = timeSeriesCollection[timeSeriesCollection.Count - 1].AsOf;
      DateTime endDate = timeSeriesCollection[0].AsOf;
      DateTime valueDate = startDate;

      while (valueDate <= endDate)
      {
        destElements = new List<Element>();
        destElements.Add(new Element(Utility.DateToLong(valueDate), 0));
        CatmullRom.PerformSpline((Element[])sourceElements.ToArray(), (Element[])destElements.ToArray());
        double value = destElements[0].Row;
        values.Add((float)value);
        valueDate = new DateTime(valueDate.Year+1,valueDate.Month,valueDate.Day);
      }
      List<float> inverseValues = new List<float>();
      for (int index = values.Count - 1; index >= 0; index--)
      {
        inverseValues.Add(values[index]);
      }
      float[] fvalues = inverseValues.ToArray();
      double averageReturn = Numerics.AverageReturn(ref fvalues);
      return averageReturn;
    }
    /// <summary>AverageReturnWithSpline -This is a collection of prices ordered highest date first.</summary>
    /// <param name="values">values</param>
    /// <returns>AverageReturnWithSpline</returns>
    public static double AverageReturnWithSpline(TimeSeriesCollection timeSeriesCollection,ReturnItems returnItems)
    {
      List<Element> sourceElements = new List<Element>();
      List<Element> destElements = null;

      returnItems.Clear();
      if (null == timeSeriesCollection || 0 == timeSeriesCollection.Count) return double.NaN;
      for (int index = timeSeriesCollection.Count - 1; index >= 0; index--)
      {
        TimeSeriesElement timeSeriesElement = timeSeriesCollection[index];
        sourceElements.Add(new Element(Utility.DateToLong(timeSeriesElement.AsOf), timeSeriesElement.Value));
      }
      DateTime startDate = timeSeriesCollection[timeSeriesCollection.Count - 1].AsOf;
      DateTime endDate = timeSeriesCollection[0].AsOf;
      DateTime valueDate = startDate;

      while (valueDate <= endDate)
      {
        destElements = new List<Element>();
        destElements.Add(new Element(Utility.DateToLong(valueDate), 0));
        CatmullRom.PerformSpline((Element[])sourceElements.ToArray(), (Element[])destElements.ToArray());
        double value = destElements[0].Row;
        returnItems.Add(new ReturnItem(valueDate,value));
        valueDate = new DateTime(valueDate.Year + 1, valueDate.Month, valueDate.Day);
      }
      returnItems.CalculateReturns();
      float[] values = null;
      values = returnItems.ToFloat();
      double averageReturn = Numerics.Mean(ref values);
      return averageReturn;
    }
    /// <summary>AverageReturn -Assumes this is a collection of prices ordered highest date first.</summary>
    /// <param name="timeSeriesColletion">timeSeriesCollection,if exclusionThreshold is provided then return values>than threshold will be thrown away</param>
    /// <returns>AverageReturn</returns>
    public static double AverageReturn(TimeSeriesCollection timeSeriesCollection)
    {
      if(null==timeSeriesCollection || 0==timeSeriesCollection.Count)return double.NaN;
      float[] values = timeSeriesCollection.ToFloat();
      return AverageReturn(ref values);
    }
    /// <summary>AverageReturn -Assumes this is a collection of prices ordered highest date first.</summary>
    /// <param name="values">values,if exclusionThreshold is provided then return values>than threshold will be thrown away</param>
    /// <returns>AverageReturn</returns>
    public static double AverageReturn(ref float[] values)
    {
      try
      {
        if (0 == values.Length) return double.NaN;
        List<double> returns = new List<double>();
        for (int index = 0; index < values.Length - 1; index++)
        {
          float current = values[index];
          float prev = values[index + 1];
          if (0 == prev) returns.Add(0.00);
          else
          {
            double currentReturn=double.NaN;
            if (prev == 0.00) currentReturn = 0;
            else currentReturn = (float)((current - prev) / Math.Abs(prev));
            returns.Add(currentReturn);   
          }
        }
        double[] returnsArray = returns.ToArray();
        return Mean(ref returnsArray);
      }
      catch (Exception)
      {
        return double.NaN;
      }
    }
    /// <summary>AverageReturnTop -Assumes this is a collection of prices ordered highest date first.</summary>
    /// <param name="values">values</param>
    /// <returns>AverageReturnTop - Calculates return of top n elements, count must be greater than 1</returns>
    public static double AverageReturnTop(ref float[] values,int top)
    {
      try
      {
        if (0 == values.Length || 1==values.Length || values.Length<top) return double.NaN;
        double[] returns = new double[top-1];
        for (int index = 0; index < top-1; index++)
        {
          float current = values[index];
          float prev = values[index + 1];
          if (0 == prev) returns[index] = 0;
          else returns[index] = (float)((current - prev) / Math.Abs(prev));
        }
        return Mean(ref returns);
      }
      catch (Exception)
      {
        return double.NaN;
      }
    }
/// <summary>PeriodicReturn -Assumes this is a collection of daily prices ordered highest date first.</summary>
/// <param name="values">values</param>
/// <returns>PeriodicReturn - Calculates return of elements</returns>
    public static double AnnualReturn(ref float[] values)
    {
      double periodicReturn = double.NaN;
      double p0 = values[values.Length - 1];  
      double p1 = values[0];
      double returnValue=(p1-p0)/Math.Abs(p0);
      double n = values.Length / 365;
      periodicReturn = Math.Pow(1 + returnValue, 1 / n) - 1;
      return periodicReturn;
    }
/// <summary>PowerHit - Returns power hitting of vector x.</summary>
/// <param name="values">values</param>
/// <returns>power hit</returns>
    public static float PowerHitting(ref float[] values)
    {
      float win=0;
      long winCount=0;
      float lose=0;
      long loseCount=0;

      if(0==values.Length)return float.NaN;
      foreach(float value in values)
      {
        if(value>0)
        {
          win+=value;
          winCount++;
        }
        else if(value<0)
        {
          lose+=(-value);
          loseCount++;
        }
      }
      return (loseCount*win/(winCount*lose));
    }
/// <summary>Product - Returns sum of products.</summary>
/// <param name="values">values</param>
/// <returns>product</returns>
    public static float Product(ref float[] values)
    {
      float product=1;
      if(0==values.Length)return 0;
      foreach(float value in values)
      {
        product*=value;
      }
      return product;
    }
/// <summary>CumProduct - Returns cumulative product of values.</summary>
/// <param name="values">values</param>
/// <returns>cumulative product</returns>
    public static float[] CumProduct(ref float[] values)
    {
      if(0==values.Length)return new float[0];
      float[] cumProduct=new float[values.Length];
      cumProduct[0]=values[0];
      for(int index=1;index<values.Length;index++)
      {
        cumProduct[index]=cumProduct[index-1]*values[index];
      }
      return cumProduct;
    }
/// <summary>StdDev - Calculates standard deviation of values.</summary>
/// <param name="values">values</param>
/// <returns>standard deviation</returns>
    public static float StdDev(ref float[] values)
    {
      long length = 0;// values.Length;
      float value;
      float sum=0;
      float sqr=0;

      if(null==values || 1>=(length=values.Length))return float.NaN;
      for(int index=0;index<length;index++)
      {
        value=values[index];
        sum+=value;
        sqr+=(value*value);
      }
      return (float)Math.Sqrt(Math.Abs((length*sqr-(sum*sum))/(length*(length-1))));
    }
/// <summary>StdDev - Calculates standard deviation of values.</summary>
/// <param name="values">values</param>
/// <returns>standard deviation</returns>
    public static double StdDev(ref double[] values)
    {
      long length = 0; ;
      double value;
      double sum=0;
      double sqr=0;

      if(null==values || 1>=(length=values.Length))return float.NaN;
      for(int index=0;index<length;index++)
      {
        value=values[index];
        sum+=value;
        sqr+=(value*value);
      }
      return (double)Math.Sqrt(Math.Abs((length*sqr-(sum*sum))/(length*(length-1))));
    }
/// <summary>StDevWithDecay - Calculates StDev with a Decay.</summary>
/// <param name="x">float array x.</param>
/// <returns>float StDev</returns>
    public static float StdDevWithDecay(ref float[] x)
    {
      if(1>=x.Length)return 0;
      return ApplyDecay(ref x);
    }
/// <summary>Covariance - Calculates covariance of x and y.</summary>
/// <param name="x">vector x</param>
/// <param name="y">vector y.</param>
/// <returns>covariance</returns>
///In probability theory and statistics, covariance is a measure of the joint variability of two random variables.
///[1] If the greater values of one variable mainly correspond with the greater values of the other variable, 
///and the same holds for the lesser values, (i.e., the variables tend to show similar behavior), 
///the covariance is positive.[2] In the opposite case, when the greater values of one variable mainly 
///correspond to the lesser values of the other, (i.e., the variables tend to show opposite behavior), 
///the covariance is negative. The sign of the covariance therefore shows the tendency in the linear 
///relationship between the variables. The magnitude of the covariance is not easy to interpret because 
///it is not normalized and hence depends on the magnitudes of the variables. The normalized version of the 
///covariance, the correlation coefficient, however, shows by its magnitude the strength of the linear relation.  
    public static float Covariance(ref float[] x,ref float[] y)
    {
      float meanx;
      float meany;
      int length;
      float product;

      if(0==x.Length && 0==y.Length)return 0;
      if((length=x.Length)!=y.Length)return 0;
      product=0;
      meanx=Mean(ref x);
      meany=Mean(ref y);
      for(int index=0;index<length;index++)
      {
        product+=(x[index]-meanx)*(y[index]-meany);
      }
      return product/length;
    }
/// <summary>Covariance - Calculates covariance of x and y.</summary>
/// <param name="x">vector x</param>
/// <param name="y">vector y.</param>
/// <returns>covariance</returns>
    public static double Covariance(ref double[] x,ref double[] y)
    {
      double meanx;
      double meany;
      int length;
      double product;

      if(0==x.Length && 0==y.Length)return 0;
      if((length=x.Length)!=y.Length)return 0;
      product=0;
      meanx=Mean(ref x);
      meany=Mean(ref y);
      for(int index=0;index<length;index++)
      {
        product+=(x[index]-meanx)*(y[index]-meany);
      }
      return product/length;
    }

/// <summary>Variance - Calculates variance of x (population).</summary>
/// <param name="x">vector x</param>
/// <returns>variance</returns>
	  public static float Variance(ref float[] xarray)
	  {
     return Covariance(ref xarray, ref xarray);
	  }
/// <summary>Variance - Calculates variance of x (population).</summary>
/// <param name="x">vector x</param>
/// <returns>variance</returns>
	  public static double Variance(ref double[] xarray)
	  {
     return Covariance(ref xarray, ref xarray);
	  }
/// <summary>Skewness - Calculates skewness of x (population).</summary>
/// <param name="x">vector x</param>
/// <returns>skewness</returns>
	  public static float Skewness(ref float[] x)
	  {
		  float meanx;
		  float sd_x;
		  int length;
		  float product;

		  if(0==x.Length)return float.NaN;
		  product=0;
		  meanx=Mean(ref x);
		  sd_x=StdDev(ref x);
		  if(float.NaN==sd_x)return float.NaN;
		  length=x.Length;
		  for(int index=0;index<length;index++)
		  {
			  product+=(x[index]-meanx)*(x[index]-meanx)*(x[index]-meanx);
		  }
		  return product/(length*sd_x*sd_x*sd_x);
	  }
/// <summary>Kurtosis - Calculates kurtosis of x.</summary>
/// <param name="x">vector x</param>
/// <returns>kurtosis</returns>
	  public static float Kurtosis(ref float[] x)
	  {
		  float meanx;
		  float sd_x;
		  int length;
		  float product;

		  if(0==x.Length)return float.NaN;
		  product=0;
		  meanx=Mean(ref x);
		  sd_x=StdDev(ref x);
		  if(float.NaN==sd_x) return float.NaN;
		  length=x.Length;
		  for(int index=0;index<length;index++)
		  {
			  product+=(x[index]-meanx)*(x[index]-meanx)*(x[index]-meanx)*(x[index]-meanx);
		  }
		  return product/(length*sd_x*sd_x*sd_x*sd_x);
	  }

/// <summary>RemoveDuplicates
/// Removes duplicates from x and y and stores the rest of the values in result.
/// This is a helper function for KStest.</summary>
/// <param name="x">vector x</param>
/// <param name="x">vector y</param>
/// <param name="x">ArrayList result</param>
/// <returns>x without duplicates</returns>
	  private static float[] RemoveDuplicates(float[] x,float[] y,ArrayList result)
	  {
		  ArrayList list= new ArrayList();
		  int newLength=0;
		  for(int i=0; i<x.Length; i++)
		  {
			  if(!list.Contains(x[i]))
			  {
				  list.Add(x[i]);
				  result.Add(y[i]);
				  newLength++;
			  }
		  }
		  float[] ret = new float[newLength];
		  list.CopyTo(ret);

		  return ret;
	  }


/* KSTEST Single sample Kolmogorov-Smirnov goodness-of-fit hypothesis test.
   KStest(X) performs a Kolmogorov-Smirnov (K-S) test
   to determine if a random sample X could have the hypothesized, continuous
   cumulative distribution function CDF (standard normal) at the default = 0.05
   significance level; this is a 2-tailed test. The following info is returned
   as part of KSTest:
   (1) H, indicating the result of the hypothesis test:
       H = 0 => Do not reject the null hypothesis at 0.05 significance level.
       H = 1 => Reject the null hypothesis at 0.05 significance level.
   (2) pValue: the asymptotic P-value P.
   (3) KSstat: the K-S test statistic, T = max|S(x) - CDF|, where S(x) is the
       empirical c.d.f. estimated from the sample vector X and CDF is the
	   standard normal.
   (4) criticalValue: the critical value of the test, on which the decision
       to reject the null hypothesis is based.

   NaN's are removed and a blank is returned if the length of the input array is 0.

   The implementation of this function is based on the Matlab function kstest.m
*/
	  public static KSTest KSTest(ref float[]x)
	  {
		  float alpha=0.05f;
		   KSTest blank=new KSTest(0,Double.NaN,Double.NaN,Double.NaN);

		  ArrayList xNoNaNs=new ArrayList();
		  for(int i=0; i<x.Length; i++)
			  if(!Single.IsNaN(x[i]))
				  xNoNaNs.Add(x[i]);
		  float[] xNoNaNsarr=new float[xNoNaNs.Count];
		  xNoNaNs.CopyTo(xNoNaNsarr);
		  x=xNoNaNsarr;
		  if(0==x.Length)
		  {
			  return blank;
		  }

		  int n=x.Length;
		  for(int i=0; i<n; i++)
			  if(Single.IsNaN(x[i]))
				  return blank;


		  float [] yCDFcdfcalc=new float[n+1];


		  Array.Sort(x);
		  yCDFcdfcalc[0]=0;
		  for(int i=1; i<=n;i++)
		  {
			  yCDFcdfcalc[i]=(float)i/n;
		  }
		  ArrayList sampleCDFarr=new ArrayList();
		  float[] xCDFcdfcalc=RemoveDuplicates(x,yCDFcdfcalc,sampleCDFarr);
		  sampleCDFarr.Add(yCDFcdfcalc[yCDFcdfcalc.Length-1]);
		  float [] sampleCDF=new float[sampleCDFarr.Count];

		  sampleCDFarr.CopyTo(sampleCDF);


		  float[] xCDF=new float[xCDFcdfcalc.Length];
		  Array.Copy(xCDFcdfcalc,xCDF,xCDFcdfcalc.Length);

		  double[] minusZsqrt = new double[xCDFcdfcalc.Length];
		  for(int i=0; i<xCDFcdfcalc.Length; i++)
			  minusZsqrt[i]=(xCDFcdfcalc[i]*-1.0)/Math.Sqrt(2);

		  double[] yCDF=ERFCore(minusZsqrt,1);
		  for(int i=0; i<yCDF.Length; i++)
		  {
			  yCDF[i]=yCDF[i]*.5;
		  }

		  double[] nullCDF= new double[yCDF.Length];
		  Array.Copy(yCDF,nullCDF,yCDF.Length);

		  double[] delta1= new double[sampleCDF.Length-1];
		  double[] delta2= new double[sampleCDF.Length-1];
		  for(int i=0; i<delta1.Length; i++)
		  {
			  delta1[i]=sampleCDF[i]-nullCDF[i];
			  delta2[i]=sampleCDF[i+1]-nullCDF[i];
		  }
		  double[] deltaCDF= new double[delta1.Length+delta2.Length];
		  Array.Copy(delta1,0,deltaCDF,0,delta1.Length);
		  Array.Copy(delta2,0,deltaCDF,delta1.Length,delta2.Length);
		  for(int i=0; i<deltaCDF.Length; i++)
		  {
			  deltaCDF[i]=Math.Abs(deltaCDF[i]);
		  }
		  float alpha1=alpha/2;

		  Array.Sort(deltaCDF);
		  double KSstatistic=deltaCDF[deltaCDF.Length-1];

		  double criticalValue;
		  if(n<=20)
		  {
			  double[] a1=new double[] {0.00500,  0.01000,  0.02500,  0.05000,  0.10000};
			  double[][] exact= new double[][]{new double[]{0.99500,  0.99000,  0.97500,  0.95000,  0.90000},new double[]{0.92929,  0.90000,  0.84189,  0.77639,  0.68377},new double[]{0.82900,  0.78456,  0.70760,  0.63604,  0.56481},new double[]{0.73424,  0.68887,  0.62394,  0.56522,  0.49265},new double[]{0.66853,  0.62718,  0.56328,  0.50945,  0.44698},new double[]{0.61661,  0.57741,  0.51926,  0.46799,  0.41037},new double[]{0.57581,  0.53844,  0.48342,  0.43607,  0.38148},new double[]{0.54179,  0.50654,  0.45427,  0.40962,  0.35831},new double[]{0.51332,  0.47960,  0.43001,  0.38746,  0.33910},new double[]{0.48893,  0.45662,  0.40925,  0.36866,  0.32260},new double[]{0.46770,  0.43670,  0.39122,  0.35242,  0.30829},new double[]{0.44905,  0.41918,  0.37543,  0.33815,  0.29577},new double[]{0.43247,  0.40362,  0.36143,  0.32549,  0.28470},new double[]{0.41762,  0.38970,  0.34890,  0.31417,  0.27481},new double[]{0.40420,  0.37713,  0.33760,  0.30397,  0.26588},new double[]{0.39201,  0.36571,  0.32733,  0.29472,  0.25778},new double[]{0.38086,  0.35528,  0.31796,  0.28627,  0.25039},new double[]{0.37062,  0.34569,  0.30936,  0.27851,  0.24360},new double[]{0.36117,  0.33685,  0.30143,  0.27136,  0.23735},new double[]{0.35241,  0.32866,  0.29408,  0.26473,  0.23156}};
			  criticalValue=exact[n-1][2];
		  }
		  else
		  {
			  double A=0.09037*Math.Pow((-1*Math.Log10(alpha1)),1.5) + 0.01515*Math.Pow((Math.Log10(alpha1)),2) - 0.08467*alpha1 - 0.11143;
			  double asymptoticStat=Math.Sqrt(-0.5*Math.Log(alpha1)/n);
			  double criticalValue1=asymptoticStat - 0.16693/n - A/Math.Pow(n,1.5);
			  criticalValue=Math.Min(criticalValue1,1-alpha1);
		  }
		  int H=KSstatistic>criticalValue?1:0;
		  double lambda=Math.Max((Math.Sqrt(n)+0.12+0.11/Math.Sqrt(n))*KSstatistic,0);
		  int[] j=new int[101];
		  double sumP=0;
		  for(int i=0; i<j.Length; i++)
		  {
			  j[i]=i+1;
			  sumP+=(Math.Pow(-1,j[i]-1)*Math.Exp(-2*lambda*lambda*Math.Pow(j[i],2)));
		  }
		  double pValue=2*sumP;
		  if(pValue<0)pValue=0;
		  if(pValue>1)pValue=1;

		  KSTest ret=new KSTest(H,pValue,KSstatistic,criticalValue);
		  return ret;
	  }

/* ERFC Complementary error function for each element of x.
 * The function is defined as:
 * erfc(x) = 2/sqrt(pi) * integral from x to inf of exp(-t^2) dt.
 *         = 1 - erf(x).
 *
 * This is a helper function to KStest and is based on the Matlab function erfc.m
 * /
 */
	  private static double[] ERFCore(double[] x, int jint)
	  {
		  double[] result = new double[x.Length];
		  for(int i=0; i<x.Length;i++)
			  result[i]=Single.NaN;

		  double xbreak=0.46875;

		  ArrayList kList = new ArrayList();
		  int newLength=0;
		  for(int i=0; i<x.Length; i++)
			  if(Math.Abs(x[i])<=xbreak)
			  {
				  kList.Add(i);
				  newLength++;
			  }

		  int[] k=new int[newLength];
		  kList.CopyTo(k);

		  if(k.Length>0)
		  {
			  double[] a=new double[] {3.16112374387056560e00, 1.13864154151050156e02, 3.77485237685302021e02, 3.20937758913846947e03,1.85777706184603153e-1};
			  double[] b=new double[] {2.36012909523441209e01, 2.44024637934444173e02, 1.28261652607737228e03, 2.84423683343917062e03};

			  double[] y=new double[k.Length];
			  double[] z=new double[k.Length];
			  double[] xnum=new double[k.Length];
			  double[] xden=new double[k.Length];
			  for(int i=0; i<y.Length; i++)
			  {
				  y[i]=Math.Abs(x[k[i]]);
				  z[i]=y[i]*y[i];
				  xnum[i]=a[4]*z[i];
				  xden[i]=z[i];
			  }
			  for(int i=0; i<3; i++)
			  {
				  for(int j=0; j<y.Length; j++)
				  {
					  xnum[j]=(xnum[j]+a[i])*z[j];
					  xden[j]=(xden[j]+b[i])*z[j];
				  }
			  }
			  for(int i=0; i<k.Length; i++)
			  {
				  result[k[i]]=x[k[i]]*(xnum[i]+a[3])/(xden[i]+b[3]);
				  result[k[i]]=1-result[k[i]];
			  }

		  }


		  ArrayList kList1 = new ArrayList();
		  int newLength1=0;
		  for(int i=0; i<x.Length; i++)
			  if(Math.Abs(x[i])>xbreak && Math.Abs(x[i])<=4.0)
			  {
				  kList1.Add(i);
				  newLength1++;
			  }

		  int[] k1=new int[newLength1];
		  kList1.CopyTo(k1);


		  if(k1.Length>0)
		  {
			  double[] c=new double[] {5.64188496988670089e-1, 8.88314979438837594e00, 6.61191906371416295e01, 2.98635138197400131e02, 8.81952221241769090e02, 1.71204761263407058e03, 2.05107837782607147e03, 1.23033935479799725e03, 2.15311535474403846e-8};
			  double[] d=new double[] {1.57449261107098347e01, 1.17693950891312499e02, 5.37181101862009858e02, 1.62138957456669019e03, 3.29079923573345963e03, 4.36261909014324716e03, 3.43936767414372164e03, 1.23033935480374942e03};

			  double[] y=new double[k1.Length];
			  double[] z=new double[k1.Length];
			  double[] xnum=new double[k1.Length];
			  double[] xden=new double[k1.Length];
			  double[] del=new double[k1.Length];
			  for(int i=0; i<y.Length; i++)
			  {
				  y[i]=Math.Abs(x[k1[i]]);
				  xnum[i]=c[8]*y[i];
				  xden[i]=y[i];
			  }
			  for(int i=0; i<7; i++)
			  {
				  for(int j=0; j<y.Length; j++)
				  {
					  xnum[j]=(xnum[j]+c[i])*y[j];
					  xden[j]=(xden[j]+d[i])*y[j];
				  }
			  }
			  for(int i=0; i<k1.Length; i++)
			  {
				  result[k1[i]]=(xnum[i]+c[7])/(xden[i]+d[7]);
			  }

			  for(int i=0; i<y.Length; i++)
			  {
				  z[i]=(y[i]*16)>0?Math.Floor(y[i]*16)/16:Math.Ceiling(y[i]*16)/16;
				  del[i]=(y[i]-z[i])*(y[i]+z[i]);
				  result[k1[i]]=Math.Exp(z[i]*z[i]*-1)*Math.Exp(del[i]*-1)*result[k1[i]];
			  }
		  }

		  ArrayList kList2 = new ArrayList();
		  int newLength2=0;
		  for(int i=0; i<x.Length; i++)
			  if(Math.Abs(x[i])>4.0)
			  {
				  kList2.Add(i);
				  newLength2++;
			  }

		  int[] k2=new int[newLength2];
		  kList2.CopyTo(k2);

		  if(k2.Length>0)
		  {
			  double[] p=new double[] {3.05326634961232344e-1, 3.60344899949804439e-1, 1.25781726111229246e-1, 1.60837851487422766e-2, 6.58749161529837803e-4, 1.63153871373020978e-2};
			  double[] q=new double[] {2.56852019228982242e00, 1.87295284992346047e00, 5.27905102951428412e-1, 6.05183413124413191e-2, 2.33520497626869185e-3};

			  double[] y=new double[k2.Length];
			  double[] z=new double[k2.Length];
			  double[] xnum=new double[k2.Length];
			  double[] xden=new double[k2.Length];
			  double[] del=new double[k2.Length];
			  for(int i=0; i<y.Length; i++)
			  {
				  y[i]=Math.Abs(x[k2[i]]);
				  z[i]=1/(y[i]*y[i]);
				  xnum[i]=p[5]*z[i];
				  xden[i]=z[i];
			  }
			  for(int i=0; i<4; i++)
			  {
				  for(int j=0; j<y.Length; j++)
				  {
					  xnum[j]=(xnum[j]+p[i])*z[j];
					  xden[j]=(xden[j]+q[i])*z[j];
				  }
			  }
			  for(int i=0; i<k2.Length; i++)
			  {
				  result[k2[i]]=z[i]*(xnum[i]+p[4])/(xden[i]+q[4]);
				  result[k2[i]]=(1/Math.Sqrt(Math.PI) - result[k2[i]])/y[i];
			  }

			  for(int i=0; i<y.Length; i++)
			  {
				  z[i]=(y[i]*16)>0?Math.Floor(y[i]*16)/16:Math.Ceiling(y[i]*16)/16;
				  del[i]=(y[i]-z[i])*(y[i]+z[i]);
				  result[k2[i]]=Math.Exp(z[i]*z[i]*-1)*Math.Exp(del[i]*-1)*result[k2[i]];
			  }


			  ArrayList kList3 = new ArrayList();
			  int newLength3=0;
			  for(int i=0; i<result.Length; i++)
				  if(Double.IsInfinity(result[i]) || Double.IsNaN(result[i]) || Double.IsNegativeInfinity(result[i]) || Double.IsPositiveInfinity(result[i]))
				  {
					  kList3.Add(i);
					  newLength3++;
				  }

			  int[] k3=new int[newLength3];
			  kList3.CopyTo(k3);

			  for(int i=0; i<k3.Length; i++)
			  {
				  result[k3[i]]=0*k3[i];
			  }
		  }

		  ArrayList kList4 = new ArrayList();
		  int newLength4=0;
		  for(int i=0; i<x.Length; i++)
			  if(x[i]<(xbreak*-1))
			  {
				  kList4.Add(i);
				  newLength4++;
			  }

		  int[] k4=new int[newLength4];
		  kList4.CopyTo(k4);

		  for(int i=0; i<k4.Length; i++)
		  {
			  result[k4[i]]=2.0-result[k4[i]];
		  }

		  double[] resultRet= new double[result.Length];
		  Array.Copy(result,resultRet,result.Length);
		  return resultRet;
	  }
/// <summary>Correlation - Calculates correlation coefficient from input vectors.</summary>
/// <param name="x">vector x</param>
/// <param name="y">vector y.</param>
/// <returns>correlation coefficient</returns>
/// Interpretation of results:
/// -1 :  A perfect downhill (negative) linear relationship
/// -0.70 : A strong downhill (negative) linear relationship
/// �0.50. A moderate downhill (negative) relationship
/// �0.30. A weak downhill (negative) linear relationship
/// 0. No linear relationship
/// +0.30. A weak uphill (positive) linear relationship
/// +0.50. A moderate uphill (positive) relationship
/// +0.70. A strong uphill (positive) linear relationship
/// Exactly +1. A perfect uphill (positive) linear relationship
///In probability theory and statistics, covariance is a measure of the joint variability of two random variables.
///[1] If the greater values of one variable mainly correspond with the greater values of the other variable, 
///and the same holds for the lesser values, (i.e., the variables tend to show similar behavior), 
///the covariance is positive.[2] In the opposite case, when the greater values of one variable mainly 
///correspond to the lesser values of the other, (i.e., the variables tend to show opposite behavior), 
///the covariance is negative. The sign of the covariance therefore shows the tendency in the linear 
///relationship between the variables. The magnitude of the covariance is not easy to interpret because 
///it is not normalized and hence depends on the magnitudes of the variables. The normalized version of the 
///covariance, the correlation coefficient, however, shows by its magnitude the strength of the linear relation.  
    public static float Correlation(ref float[] x, ref float[] y)
    {
			float sd_x = StdDev(ref x);
			float sd_y = StdDev(ref y);
			return Covariance(ref x, ref y)*x.Length/((sd_x * sd_y)*(x.Length-1));
    }
    public static Numerics.CorrelationType GetCorrelationType(float correlationCovarianeResult)
    {
      if(correlationCovarianeResult>=1)return CorrelationType.PerfectUphill;
      if(correlationCovarianeResult>=.70)return CorrelationType.StrongUphill;
      if(correlationCovarianeResult>=.50)return CorrelationType.ModerateUphill;
      if(correlationCovarianeResult>=.50)return CorrelationType.ModerateUphill;
      if(correlationCovarianeResult>=.30)return CorrelationType.WeakUphill;
      if(correlationCovarianeResult<=-1)return CorrelationType.PerfectDownhill;
      if(correlationCovarianeResult<=-.70)return CorrelationType.StrongDownhill;
      if(correlationCovarianeResult<=-.50)return CorrelationType.ModerateDownhill;
      if(correlationCovarianeResult<=-.30)return CorrelationType.WeakDownhill;
      return CorrelationType.NoLinearRelationship;
    }

/// <summary>Correlation - Calculates correlation matrix from input matrix.</summary>
/// <param name="matrix">source matrix</param>
/// <param name="rows">number of columns in source.</param>
/// <param name="cols">number of columns in source.</param>
/// <returns>tranposed array</returns>
    public static float[,] Correlation(ref float[,] matrix,int rows,int cols)
    {
      float[,] result=new float[cols,cols];
      float[] x=new float[rows];
      float[] y=new float[rows];
      int col=0;
      int row;

      for(int srcIndex=0;srcIndex<cols;srcIndex++)
      {
        row=0;
        for(int dstIndex=0;dstIndex<cols;dstIndex++)
        {
          GetXY(srcIndex,dstIndex,ref matrix,ref x,ref y,rows);
          result[row,col]=Correlation(ref x,ref y);
          row++;
        }
        col++;
      }
      return result;
    }
/// <summary>CreateCorrelationMatrix - Creates correlation matrix from input matrix.</summary>
/// <param name="values">the input matrix.</param>
/// <param name="transpose">indicates whether to transpose the input matrix.</param>
/// <returns>float[,] correlation matrix</returns>
    public static float[,] CreateCorrelationMatrix(ref float[,] values,bool transpose)
    {
      int columns=values.GetLength(1);
      int rows=values.GetLength(0);
      if(!transpose)return Correlation(ref values,rows,columns);
      float[,] transposeMatrix=Transpose(ref values,ref rows,ref columns);
      return Correlation(ref transposeMatrix,rows,columns);
    }
/// <summary>Transpose - Convert horizontal range of cells to vertical.</summary>
/// <param name="srcRows">number of rows in source</param>
/// <param name="srcCols">number of columns in source.</param>
/// <param name="source">the source matrix.</param>
/// <returns>tranposed array</returns>
    public static float[] Transpose(ref int srcRows,ref int srcCols,ref float[] source)
    {
      float[] result=new float[srcRows*srcCols];
      long dstRow=0;
      long dstCol=0;

      for(int srcRow=0;srcRow<srcRows;srcRow++)
      {
        dstRow=0;
        for(int srcCol=0;srcCol<srcCols;srcCol++)
        {
          result[(dstRow*srcRows)+dstCol]=source[(srcRow*srcCols)+srcCol];
          dstRow++;
        }
        dstCol++;
      }
      int tmpRows=srcRows;
      srcRows=srcCols;
      srcCols=tmpRows;
      return result;
    }

/// <summary>Transpose - Convert horizontal range of cells to vertical.</summary>
/// <param name="srcRows">number of rows in source</param>
/// <param name="srcCols">number of columns in source.</param>
/// <param name="source">the source matrix.</param>
/// <returns>tranposed array</returns>
    public static float[,] Transpose(ref float[,] source,ref int srcRows,ref int srcCols)
    {
      float[,] result=new float[srcCols,srcRows];
      long dstRow=0;
      long dstCol=0;

      for(int srcRow=0;srcRow<srcRows;srcRow++)
      {
        dstRow=0;
        for(int srcCol=0;srcCol<srcCols;srcCol++)
        {
          result[dstRow,dstCol]=source[srcRow,srcCol];
          dstRow++;
        }
        dstCol++;
      }
      int tmpRows=srcRows;
      srcRows=srcCols;
      srcCols=tmpRows;
      return result;
    }
/// <summary>TransposeXAY - Transpose(x)*A*y.</summary>
/// <param name="x">row vector</param>
/// <param name="a">matrix.</param>
/// <param name="y">column vector.</param>
/// <param name="elements">number of elements in each dimension.</param>
/// <returns>tranposed array</returns>
    public static float TransposeXAY(ref float[] x,ref float[,] a,ref float[] y,int elements)
    {
      float[] acol=new float[elements];
      float[] intermediate=new float[elements];

      for(int col=0;col<elements;col++)
      {
        GetX(col,ref a,ref acol,elements);
	    intermediate[col]=MMult(ref x,ref acol);
      }
	  return MMult(ref intermediate,ref y);
    }
/// <summary>MMult - Performs matrix multiplication of vectors.</summary>
/// <param name="rowVector">the row vector</param>
/// <param name="colVector">the column vector</param>
/// <returns>Product</returns>
    public static float MMult(ref float[] rowVector,ref float[] colVector)
    {
      float product=0;
      int length;

      if((length=rowVector.Length)!=colVector.Length)return product;
      for(int index=0;index<length;index++)product+=rowVector[index]*colVector[index];
      return product;
    }
/// <summary>GetXY - Get vertical columns from matrix[,] into x and y.</summary>
/// <param name="xColumn">column from matrix to copy to x array</param>
/// <param name="yColumn">column from matrix to copy to y array</param>
/// <param name="matrix">source matrix</param>
/// <param name="x">destination array x</param>
/// <param name="y">destination array y</param>
/// <returns>None</returns>
    private static void GetXY(int xColumn,int yColumn,ref float[,] matrix,ref float[] x,ref float[] y,int rows)
    {
      for(int row=0;row<rows;row++)
      {
        x[row]=matrix[row,xColumn];
        y[row]=matrix[row,yColumn];
      }
    }
/// <summary>GetX - Get vertical column from matrix[,] into x.</summary>
/// <param name="xColumn">column from matrix to copy</param>
/// <param name="matrix">the source matrix.</param>
/// <param name="x">the destination vector.</param>
/// <param name="rows">number of rows in matrix.</param>
/// <returns>None</returns>
    private static void GetX(int xColumn,ref float[,] matrix,ref float[] x,int rows)
    {
      for(int row=0;row<rows;row++)
      {
        x[row]=matrix[row,xColumn];
      }
    }
/// <summary>Beta - Calculates Beta.</summary>
/// <param name="x">ticker returns.</param>
/// <param name="y">sector returns.</param>
/// <returns>Beta value</returns>
/// A beta of 1.00 indicates that the security's price will move with the market
/// A beta of less than 1.00 indicates that the security will be less volatile than the market
/// A beta of more than 1.00 indicates that the security will be more volatile than the market
/// The revised code (i.e.) Beta=Covariance(x,y)/Variance(x) was taken from Investopedia 
    public static double Beta(ref double[] assetReturns, ref double[] benchmarkReturns)
    {
      double variance = Variance(ref benchmarkReturns);
      if(0.00 == variance)return double.NaN;
      return Covariance(ref assetReturns, ref benchmarkReturns) / variance;
    }
/// <summary>ApplyDecay - Apply exponential weighting.</summary>
/// <param name="samples">the samples to weight.</param>
/// <returns>The weighted standard deviation</returns>
    private static float ApplyDecay(ref float[] samples)
    {
      int numSamples=samples.Length;
      Decay decay=new Decay(numSamples);
      float average=0;
      float weightedVariance=0;
      float weightedError=0;
      float stddev=0;
      float error=0;

      for(int item=0;item<numSamples;item++)average+=(float)samples[item];
      average/=numSamples;
      for(int item=0;item<numSamples;item++)
      {
        error=(float)Math.Pow(samples[item]-average,2);
	      weightedError=error*decay[item]*(numSamples);
        weightedVariance+=weightedError;
      }
      weightedVariance/=(numSamples-1);
      stddev=(float)Math.Sqrt(Math.Abs(weightedVariance));
      return stddev;
    }
/// <summary>TStatistic - Calculate TStatistic.</summary>
/// <param name="average">The average for the samples.</param>
/// <param name="stddev">Standard deviation of the samples.</param>
/// <param name="samples">The number of samples.</param>
/// <returns>TStatistic</returns>
    public static float TStatistic(float average,float stddev,int samples)
    {
      if(0==samples||0==stddev)return float.NaN;
      return (float)((average/stddev)*Math.Sqrt(samples));
    }
/// <summary>TStatistic - Calculate TStatistic.</summary>
/// <param name="leastSquaresEstimates">Least squares estimates.</param>
/// <param name="standardErrors">Standard errors.</param>
/// <returns>TStatistic</returns>
     public static double[] TStatistic(ref double[] leastSquaresEstimates,ref double[] standardErrors)
	   {
	    double[] tStats;
	    int length;
	    if((length=leastSquaresEstimates.Length)!=standardErrors.Length)return new double[0];
	    tStats=new double[length];
	    for(int index=0;index<length;index++)
	    {
	      if(0==standardErrors[index])tStats[index]=0;
		  else tStats[index]=leastSquaresEstimates[index]/standardErrors[index];
	    }
	    return tStats;
	   }
// Calculate max drawdown of series
     public static double MaxDrawdown(double[] prices) 
     {
       if (prices.Length <= 1) return 0;
       double maxPrice = prices[0];
       double maxDd = 0;
       for (int i = 1; i < prices.Length; i++)
       {
         if (prices[i] > maxPrice) maxPrice = prices[i];
         else if (prices[i] < maxPrice) maxDd = Math.Min(maxDd, prices[i] / maxPrice - 1);
       }
       return maxDd;
     }
// a float version
     public static float MaxDrawdown(float[] prices) 
     {
       if (prices.Length <= 1) return 0;
       float maxPrice = prices[0];
       float maxDd = 0;
       for (int i = 1; i < prices.Length; i++)
       {
         if (prices[i] > maxPrice) maxPrice = prices[i];
         else if (prices[i] < maxPrice) maxDd = Math.Min(maxDd, prices[i] / maxPrice - 1);
       }
       return maxDd;
     }
// float MaxUpside
     public static float MaxUpside(float[] prices) 
     {
       if (prices.Length <= 1) return 0;
       float minPrice = prices[0];
       float maxUpside = 0;
       for (int i = 1; i < prices.Length; i++)
       {
         if (prices[i] < minPrice) minPrice = prices[i];
         else if (prices[i] > minPrice) maxUpside = Math.Max(maxUpside, prices[i] / minPrice - 1);
       }
       return maxUpside;
     }
  }
}