source: branches/2974_Constants_Optimization/HeuristicLab.Problems.DataAnalysis.Symbolic/3.4/ConstantsOptimization/LMConstantsOptimizer.cs @ 16514

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2018 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.Encodings.SymbolicExpressionTreeEncoding;

namespace HeuristicLab.Problems.DataAnalysis.Symbolic.ConstantsOptimization {
  public class LMConstantsOptimizer {

    private LMConstantsOptimizer() { }

    /// <summary>
    /// Method to determine whether the numeric constants of the tree can be optimized. This depends primarily on the symbols occuring in the tree.
    /// </summary>
    /// <param name="tree">The tree that should be analyzed</param>
    /// <returns>A flag indicating whether the numeric constants of the tree can be optimized</returns>
    public static bool CanOptimizeConstants(ISymbolicExpressionTree tree) {
      return AutoDiffConverter.IsCompatible(tree);
    }

    /// <summary>
    /// Optimizes the numeric constants in a symbolic expression tree in place.
    /// </summary>
    /// <param name="tree">The tree for which the constants should be optimized</param>
    /// <param name="dataset">The dataset containing the data.</param>
    /// <param name="targetVariable">The target variable name.</param>
    /// <param name="rows">The rows for which the data should be extracted.</param>
    /// <param name="applyLinearScaling">A flag to determine whether linear scaling should be applied during the optimization</param>
    /// <param name="maxIterations">The maximum number of iterations of the Levenberg-Marquard algorithm.</param>
    /// <returns></returns>
    public static double OptimizeConstants(ISymbolicExpressionTree tree,
      IDataset dataset, string targetVariable, IEnumerable<int> rows,
      bool applyLinearScaling, int maxIterations = 10) {
      if (tree == null) throw new ArgumentNullException("tree");
      if (dataset == null) throw new ArgumentNullException("dataset");
      if (!dataset.ContainsVariable(targetVariable)) throw new ArgumentException("The dataset does not contain the provided target variable.");

      var allVariables = Util.ExtractVariables(tree);
      var numericNodes = Util.ExtractNumericNodes(tree);

      AutoDiff.IParametricCompiledTerm term;
      if (!AutoDiffConverter.TryConvertToAutoDiff(tree, applyLinearScaling, numericNodes, allVariables, out term))
        throw new NotSupportedException("Could not convert symbolic expression tree to an AutoDiff term due to not supported symbols used in the tree.");

      //Variables of the symbolic expression tree correspond to parameters in the term
      //Hence if no parameters are present no variables occur in the tree and the R² = 0
      if (term.Parameters.Count == 0) return 0.0;

      var initialConstants = Util.ExtractConstants(numericNodes, applyLinearScaling);
      double[] constants;
      double[,] x = Util.ExtractData(dataset, allVariables, rows);
      double[] y = dataset.GetDoubleValues(targetVariable, rows).ToArray();

      var result = OptimizeConstants(term, initialConstants, x, y, maxIterations, out constants);
      if (result > 0.0 && constants.Length != 0)
        Util.UpdateConstants(numericNodes, constants);

      return result;
    }

    /// <summary>
    /// Optimizes the numeric coefficents of an AutoDiff Term using the Levenberg-Marquard algorithm.
    /// </summary>
    /// <param name="term">The AutoDiff term for which the numeric coefficients should be optimized.</param>
    /// <param name="initialConstants">The starting values for the numeric coefficients.</param>
    /// <param name="x">The input data for the optimization.</param>
    /// <param name="y">The target values for the optimization.</param>
    /// <param name="maxIterations">The maximum number of iterations of the Levenberg-Marquard</param>
    /// <param name="constants">The opitmized constants.</param>
    /// <param name="LM_IterationCallback">An optional callback for detailed analysis that is called in each algorithm iteration.</param>
    /// <returns>The R² of the term evaluated on the input data x and the target data y using the optimized constants</returns>
    public static double OptimizeConstants(AutoDiff.IParametricCompiledTerm term, double[] initialConstants, double[,] x, double[] y,
      int maxIterations, out double[] constants, Action<double[], double, object> LM_IterationCallback = null) {

      if (term.Parameters.Count == 0) {
        constants = new double[0];
        return 0.0;
      }

      var optimizedConstants = (double[])initialConstants.Clone();
      int numberOfRows = x.GetLength(0);
      int numberOfColumns = x.GetLength(1);
      int numberOfConstants = optimizedConstants.Length;

      alglib.lsfitstate state;
      alglib.lsfitreport rep;
      alglib.ndimensional_rep xrep = (p, f, obj) => LM_IterationCallback(p, f, obj);
      int retVal;

      try {
        alglib.lsfitcreatefg(x, y, optimizedConstants, numberOfRows, numberOfColumns, numberOfConstants, cheapfg: false, state: out state);
        alglib.lsfitsetcond(state, 0.0, 0.0, maxIterations);
        alglib.lsfitsetxrep(state, LM_IterationCallback != null);
        alglib.lsfitfit(state, Evaluate, EvaluateGradient, xrep, term);
        alglib.lsfitresults(state, out retVal, out optimizedConstants, out rep);
      } catch (ArithmeticException) {
        constants = new double[0];
        return double.NaN;
      } catch (alglib.alglibexception) {
        constants = new double[0];
        return double.NaN;
      }

      constants = optimizedConstants;
      return rep.r2;
    }


    private static void Evaluate(double[] c, double[] x, ref double fx, object o) {
      AutoDiff.IParametricCompiledTerm term = (AutoDiff.IParametricCompiledTerm)o;
      fx = term.Evaluate(c, x);
    }

    private static void EvaluateGradient(double[] c, double[] x, ref double fx, double[] grad, object o) {
      AutoDiff.IParametricCompiledTerm term = (AutoDiff.IParametricCompiledTerm)o;
      Tuple<double[], double> result = term.Differentiate(c, x);
      fx = result.Item2;
      Array.Copy(result.Item1, grad, grad.Length);
    }
  }
}