source: branches/NCA/HeuristicLab.Algorithms.NCA/3.3/NeighborhoodComponentsAnalysis.cs @ 8441

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2012 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * This file is part of HeuristicLab.
 *
 * HeuristicLab 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.
 *
 * HeuristicLab 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 HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.Collections.Generic;
using System.Linq;
using HeuristicLab.Algorithms.DataAnalysis;
using HeuristicLab.Analysis;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.PluginInfrastructure;
using HeuristicLab.Problems.DataAnalysis;

namespace HeuristicLab.Algorithms.NCA {
  public delegate void Reporter(double quality, double[] coefficients);
  /// <summary>
  /// Neighborhood Components Analysis
  /// </summary>
  [Item("Neighborhood Components Analysis", "NCA is described in J. Goldberger, S. Roweis, G. Hinton, R. Salakhutdinov. 2005. Neighbourhood Component Analysis. Advances in Neural Information Processing Systems, 17. pp. 513-520.")]
  [Creatable("Data Analysis")]
  [StorableClass]
  public sealed class NeighborhoodComponentsAnalysis : FixedDataAnalysisAlgorithm<IClassificationProblem> {
    #region Parameter Properties
    public IValueLookupParameter<IntValue> KParameter {
      get { return (IValueLookupParameter<IntValue>)Parameters["k"]; }
    }
    public IValueLookupParameter<IntValue> ReduceDimensionsParameter {
      get { return (IValueLookupParameter<IntValue>)Parameters["ReduceDimensions"]; }
    }
    private IConstrainedValueParameter<INCAInitializer> InitializationParameter {
      get { return (IConstrainedValueParameter<INCAInitializer>)Parameters["Initialization"]; }
    }
    #endregion

    #region Properties
    public IntValue K {
      get { return KParameter.Value; }
    }
    public IntValue ReduceDimensions {
      get { return ReduceDimensionsParameter.Value; }
    }
    #endregion

    [StorableConstructor]
    private NeighborhoodComponentsAnalysis(bool deserializing) : base(deserializing) { }
    private NeighborhoodComponentsAnalysis(NeighborhoodComponentsAnalysis original, Cloner cloner) : base(original, cloner) { }
    public NeighborhoodComponentsAnalysis()
      : base() {
      Parameters.Add(new ValueLookupParameter<IntValue>("k", "The k for the nearest neighbor.", new IntValue(1)));
      Parameters.Add(new ValueLookupParameter<IntValue>("ReduceDimensions", "The number of dimensions that NCA should reduce the data to.", new IntValue(2)));
      Parameters.Add(new ConstrainedValueParameter<INCAInitializer>("Initialization", "Which method should be used to initialize the matrix. Typically LDA (linear discriminant analysis) should provide a good estimate."));

      INCAInitializer defaultInitializer = null;
      foreach (var initializer in ApplicationManager.Manager.GetInstances<INCAInitializer>().OrderBy(x => x.ItemName)) {
        if (initializer is LDAInitializer) defaultInitializer = initializer;
        InitializationParameter.ValidValues.Add(initializer);
      }
      if (defaultInitializer != null) InitializationParameter.Value = defaultInitializer;

      Problem = new ClassificationProblem();
    }

    public override IDeepCloneable Clone(Cloner cloner) {
      return new NeighborhoodComponentsAnalysis(this, cloner);
    }

    public override void Prepare() {
      if (Problem != null) base.Prepare();
    }

    protected override void Run() {
      var k = K.Value;
      var dimensions = ReduceDimensions.Value;
      var initializer = InitializationParameter.Value;

      var clonedProblem = (IClassificationProblemData)Problem.ProblemData.Clone();
      var model = Train(clonedProblem, k, dimensions, initializer, ReportQuality);
      var classification = new NCAClassificationSolution(clonedProblem, model);
      Results.Add(new Result("ClassificationSolution", "The classification solution.", classification));
    }

    private void ReportQuality(double func, double[] coefficients) {
      var instances = Problem.ProblemData.TrainingIndices.Count();
      DataTable qualities;
      if (!Results.ContainsKey("Optimization")) {
        qualities = new DataTable("Optimization");
        qualities.Rows.Add(new DataRow("Quality", string.Empty));
        Results.Add(new Result("Optimization", qualities));
      } else qualities = (DataTable)Results["Optimization"].Value;
      qualities.Rows["Quality"].Values.Add(-func / instances);

      if (!Results.ContainsKey("Quality")) {
        Results.Add(new Result("Quality", new DoubleValue(-func / instances)));
      } else ((DoubleValue)Results["Quality"].Value).Value = -func / instances;
    }

    public static INCAModel Train(IClassificationProblemData data, int k, int dimensions, INCAInitializer initializer, Reporter reporter = null) {
      var instances = data.TrainingIndices.Count();
      var attributes = data.AllowedInputVariables.Count();

      double[] matrix = initializer.Initialize(data, dimensions);

      var info = new OptimizationInfo(data, dimensions, reporter);
      alglib.mincgstate state;
      alglib.mincgreport rep;

      alglib.mincgcreate(matrix, out state);
      alglib.mincgsetcond(state, 0, 1e-05, 0, 20);
      alglib.mincgsetxrep(state, true);
      alglib.mincgoptimize(state, Gradient, Report, info);
      alglib.mincgresults(state, out matrix, out rep);

      var transformationMatrix = new double[attributes, dimensions];
      var counter = 0;
      for (var i = 0; i < attributes; i++)
        for (var j = 0; j < dimensions; j++)
          transformationMatrix[i, j] = matrix[counter++];

      var transformedTrainingset = new double[instances, dimensions];
      var rowCount = 0;
      foreach (var r in data.TrainingIndices) {
        var i = 0;
        foreach (var v in data.AllowedInputVariables) {
          var val = data.Dataset.GetDoubleValue(v, r);
          for (var j = 0; j < dimensions; j++)
            transformedTrainingset[rowCount, j] += val * transformationMatrix[i, j];
          i++;
        }
        rowCount++;
      }

      var ds = data.Dataset;
      var targetVariable = data.TargetVariable;
      return new NCAModel(transformedTrainingset, info.Scaling, transformationMatrix, k, data.TargetVariable, data.AllowedInputVariables,
        data.TrainingIndices.Select(i => ds.GetDoubleValue(targetVariable, i)).ToArray());
    }

    private static void Report(double[] A, double func, object obj) {
      var info = (OptimizationInfo)obj;
      if (info.Reporter != null) info.Reporter(func, A);
    }

    private static void Gradient(double[] A, ref double func, double[] grad, object obj) {
      var info = (OptimizationInfo)obj;
      var data = info.Data;
      var classes = info.TargetValues;
      var instances = info.Instances;
      var attributes = info.Attributes;

      var AMatrix = new Matrix(A, A.Length / info.ReduceDimensions, info.ReduceDimensions);

      alglib.sparsematrix probabilities;
      alglib.sparsecreate(instances, instances, out probabilities);
      var transformedDistances = new double[instances];
      for (int i = 0; i < instances - 1; i++) {
        var iVector = new Matrix(GetRow(data, i), data.GetLength(1));
        var denom = 0.0;
        for (int k = 0; k < instances; k++) {
          if (k == i) continue;
          var kVector = new Matrix(GetRow(data, k));
          transformedDistances[k] = Math.Exp(-iVector.Multiply(AMatrix).Subtract(kVector.Multiply(AMatrix)).SquaredVectorLength());
          denom += transformedDistances[k];
        }
        if (denom > 1e-05) {
          for (int j = i + 1; j < instances; j++) {
            if (i == j) continue;
            var v = transformedDistances[j] / denom;
            alglib.sparseset(probabilities, i, j, v);
          }
        }
      }
      alglib.sparseconverttocrs(probabilities); // needed to enumerate in order (top-down and left-right)

      int t0 = 0, t1 = 0, r, c;
      double val;
      var pi = new double[instances];
      while (alglib.sparseenumerate(probabilities, ref t0, ref t1, out r, out c, out val)) {
        if (classes[r].IsAlmost(classes[c])) {
          pi[r] += val;
        }
      }

      var innerSum = new double[attributes, attributes];
      while (alglib.sparseenumerate(probabilities, ref t0, ref t1, out r, out c, out val)) {
        var vector = new Matrix(GetRow(data, r)).Subtract(new Matrix(GetRow(data, c)));
        vector.OuterProduct(vector).Multiply(2.0 * val * pi[r]).AddTo(innerSum);

        if (classes[r].IsAlmost(classes[c])) {
          vector.OuterProduct(vector).Multiply(-2.0 * val).AddTo(innerSum);
        }
      }

      func = -2.0 * pi.Sum();

      r = 0;
      var newGrad = AMatrix.Multiply(-2.0).Transpose().Multiply(new Matrix(innerSum)).Transpose();
      foreach (var g in newGrad) {
        grad[r++] = g;
      }
    }

    #region Helpers
    private static IEnumerable<double> GetRow(double[,] data, int row) {
      for (int i = 0; i < data.GetLength(1); i++)
        yield return data[row, i];
    }

    private class OptimizationInfo {
      public Scaling Scaling { get; private set; }
      public double[,] Data { get; private set; }
      public double[] TargetValues { get; private set; }
      public int ReduceDimensions { get; private set; }
      public int Instances { get; private set; }
      public int Attributes { get; private set; }
      public Reporter Reporter { get; private set; }

      public OptimizationInfo(IClassificationProblemData data, int reduceDimensions, Reporter reporter) {
        this.Scaling = new Scaling(data.Dataset, data.AllowedInputVariables, data.TrainingIndices);
        this.Data = AlglibUtil.PrepareAndScaleInputMatrix(data.Dataset, data.AllowedInputVariables, data.TrainingIndices, Scaling);
        this.TargetValues = data.Dataset.GetDoubleValues(data.TargetVariable, data.TrainingIndices).ToArray();
        this.ReduceDimensions = reduceDimensions;
        this.Instances = data.TrainingIndices.Count();
        this.Attributes = data.AllowedInputVariables.Count();
        this.Reporter = reporter;
      }
    }
    #endregion
  }
}