source: branches/DataAnalysis/HeuristicLab.Problems.DataAnalysis.MultiVariate/3.3/Evaluators/OnlineMeanMahalanobisDistanceEvaluator.cs @ 8543

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Problems.DataAnalysis.Evaluators;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Parameters;

namespace HeuristicLab.Problems.DataAnalysis.MultiVariate.Evaluators {
  public class OnlineMeanMahalanobisDistanceEvaluator : IMultiVariateOnlineEvaluator {
    private int n;
    private double distance;
    private double[,] covMatrix;
    private double[] diff;
    private double[] target;
    public double MeanMahalanobisDistance {
      get {
        if (n == 0) throw new InvalidOperationException("no elements");
        else
          return Math.Sqrt(distance) / n;
      }
    }

    public double MeanGeneralizedSquaredInterpointDistance {
      get {
        if (n == 0) throw new InvalidOperationException("no elements");
        else
          return distance / n;
      }
    }

    public OnlineMeanMahalanobisDistanceEvaluator() {
      Reset();
    }

    #region IMultiVariateOnlineEvaluator Members
    public double Value {
      get { return MeanMahalanobisDistance; }
    }

    public void Add(IEnumerable<double> original, IEnumerable<double> estimated) {
      if (covMatrix == null) throw new InvalidOperationException("Covariance matrix must be initialized before values can be added.");

      {
        // calculate difference vector
        var originalEnumerator = original.GetEnumerator();
        var estimatedEnumerator = estimated.GetEnumerator();
        int i = 0;
        while (originalEnumerator.MoveNext() & estimatedEnumerator.MoveNext() && i < diff.Length) {
          diff[i++] = originalEnumerator.Current - estimatedEnumerator.Current;
        }
        if (originalEnumerator.MoveNext() | estimatedEnumerator.MoveNext() || i < diff.Length) {
          throw new ArgumentException("Number of elements of original and estimated doesn't match or is not compatible with covariance matrix.");
        }
      }

      {
        // calculate mahalanobis distance using covariance matrix
        // covMatrix^(-1) * diff => target
        alglib.ablas.rmatrixmv(covMatrix.GetLength(0), covMatrix.GetLength(1), covMatrix, 0, 0, 0, diff, 0, ref target, 0);

        // diff^T * (covMatrix^(-1) * diff) => sum
        double sum = 0.0;
        for (int i = 0; i < diff.Length; i++) {
          sum += diff[i] * target[i];
        }
        distance += sum;
        n++;
      }
    }

    public void Reset() {
      n = 0;
      distance = 0.0;
      covMatrix = null;
      diff = null;
      target = null;
    }

    #endregion

    public void InitializeCovarianceMatrixFromSamples(params IEnumerable<double>[] samples) {
      covMatrix = new double[samples.Length, samples.Length];
      OnlineCovarianceEvaluator covEvaluator = new OnlineCovarianceEvaluator();
      for (int i = 0; i < samples.Length; i++) {
        for (int j = i; j < samples.Length; j++) {
          var xEnumerator = samples[i].GetEnumerator();
          var yEnumerator = samples[j].GetEnumerator();
          covEvaluator.Reset();
          while (xEnumerator.MoveNext() & yEnumerator.MoveNext()) {
            covEvaluator.Add(xEnumerator.Current, yEnumerator.Current);
          }
          if (xEnumerator.MoveNext() | yEnumerator.MoveNext()) {
            throw new ArgumentException("Number of elements must be the same in all enumerations.");
          }
          covMatrix[i, j] = covEvaluator.Covariance;
          covMatrix[j, i] = covEvaluator.Covariance;
        }
      }
      int info = 0;
      alglib.matinv.matinvreport report = new alglib.matinv.matinvreport();
      alglib.matinv.rmatrixinverse(ref covMatrix, covMatrix.GetLength(0), ref info, report);
      if (info != 1) throw new InvalidOperationException("Can't invert covariance matrix.");
      diff = new double[samples.Length];
      target = new double[samples.Length];
    }
  }
}