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

#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.Common;
using HeuristicLab.Problems.DataAnalysis;

namespace HeuristicLab.Algorithms.NCA {
  public class NeighborhoodComponentsAnalysis {

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

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

      alglib.mincgstate state;
      alglib.mincgreport rep;

      // first run
      alglib.mincgcreate(matrix, out state);
      alglib.mincgsetcond(state, 0.0000000001, 0, 0, 0);
      alglib.mincgoptimize(state, Gradient, null, new OptimizationInfo(data, reduceDimensions));
      alglib.mincgresults(state, out matrix, out rep);

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

      var transformedTrainingset = new double[instances, reduceDimensions];
      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 < reduceDimensions; j++)
            transformedTrainingset[rowCount, j] += val * transformationMatrix[i, j];
          i++;
        }
        rowCount++;
      }

      return new NCAModel(transformedTrainingset, transformationMatrix, k, data.TargetVariable, data.AllowedInputVariables,
        data.Dataset.GetDoubleValues(data.TargetVariable)
          .Select((v, i) => new { I = i, V = v })
          .Where(x => x.I >= data.TrainingPartition.Start && x.I < data.TrainingPartition.End
              && !(x.I >= data.TestPartition.Start && x.I < data.TestPartition.End))
          .Select(x => x.V).ToArray());
    }

    private static void Gradient(double[] A, ref double func, double[] grad, object obj) {
      var info = (OptimizationInfo)obj;
      var instances = info.ProblemData.TrainingIndices.ToArray();
      var attributes = info.ProblemData.AllowedInputVariables.Count();
      var AMatrix = new Matrix(A, A.Length / info.ReduceDimensions, info.ReduceDimensions);

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

      int t0 = 0, t1 = 0, r, c;
      double val;
      var classes = info.ProblemData.Dataset.GetDoubleValues(info.ProblemData.TargetVariable, instances).ToArray();
      var pi = new double[instances.Length];
      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(info.ProblemData, instances[r])).Subtract(new Matrix(GetRow(info.ProblemData, instances[c]))).Apply();
        vector.OuterProduct(vector).Multiply(val * pi[r]).AddTo(innerSum);

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

      func = -pi.Sum();

      grad = AMatrix.Multiply(-2.0).Transpose().Multiply(new Matrix(innerSum)).Transpose().ToArray();
    }

    private static IEnumerable<double> GetRow(IClassificationProblemData data, int row) {
      return data.AllowedInputVariables.Select(v => data.Dataset.GetDoubleValue(v, row));
    }

    public static NCAClassificationSolution CreateNCASolution(IClassificationProblemData problemData, int k, int reduceDimensions, INCAInitializer initializer) {
      return new NCAClassificationSolution(problemData, Train(problemData, k, reduceDimensions, initializer));
    }

    private class OptimizationInfo {
      public IClassificationProblemData ProblemData { get; set; }
      public int ReduceDimensions { get; set; }
      public OptimizationInfo(IClassificationProblemData problem, int reduceDimensions) {
        this.ProblemData = problem;
        this.ReduceDimensions = reduceDimensions;
      }
    }
  }
}