source: branches/2925_AutoDiffForDynamicalModels/HeuristicLab.Problems.DynamicalSystemsModelling/3.3/OdeParameterIdentification.cs @ 16597

#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.Linq;
using System.Threading;
using HeuristicLab.Algorithms.DataAnalysis;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.DataAnalysis;
using HeuristicLab.Problems.DataAnalysis.Symbolic;
using HeuristicLab.Random;
using System.Collections.Generic;

namespace HeuristicLab.Problems.DynamicalSystemsModelling {
  [Item("OdeParameterIdentification", "TODO")]
  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 120)]
  [StorableClass]
  public sealed class OdeParameterIdentification : FixedDataAnalysisAlgorithm<Problem> {
    private const string RegressionSolutionResultName = "Regression solution";
    private const string ModelStructureParameterName = "Model structure";
    private const string IterationsParameterName = "Iterations";
    private const string RestartsParameterName = "Restarts";
    private const string SetSeedRandomlyParameterName = "SetSeedRandomly";
    private const string SeedParameterName = "Seed";
    private const string InitParamsRandomlyParameterName = "InitializeParametersRandomly";

    public IValueParameter<StringArray> ModelStructureParameter {
      get { return (IValueParameter<StringArray>)Parameters[ModelStructureParameterName]; }
    }
    public IFixedValueParameter<IntValue> IterationsParameter {
      get { return (IFixedValueParameter<IntValue>)Parameters[IterationsParameterName]; }
    }

    public IFixedValueParameter<BoolValue> SetSeedRandomlyParameter {
      get { return (IFixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyParameterName]; }
    }

    public IFixedValueParameter<IntValue> SeedParameter {
      get { return (IFixedValueParameter<IntValue>)Parameters[SeedParameterName]; }
    }

    public IFixedValueParameter<IntValue> RestartsParameter {
      get { return (IFixedValueParameter<IntValue>)Parameters[RestartsParameterName]; }
    }

    public IFixedValueParameter<BoolValue> InitParametersRandomlyParameter {
      get { return (IFixedValueParameter<BoolValue>)Parameters[InitParamsRandomlyParameterName]; }
    }

    public StringArray ModelStructure {
      get { return ModelStructureParameter.Value; }
      set { ModelStructureParameter.Value = value; }
    }

    public int Iterations {
      get { return IterationsParameter.Value.Value; }
      set { IterationsParameter.Value.Value = value; }
    }

    public int Restarts {
      get { return RestartsParameter.Value.Value; }
      set { RestartsParameter.Value.Value = value; }
    }

    public int Seed {
      get { return SeedParameter.Value.Value; }
      set { SeedParameter.Value.Value = value; }
    }

    public bool SetSeedRandomly {
      get { return SetSeedRandomlyParameter.Value.Value; }
      set { SetSeedRandomlyParameter.Value.Value = value; }
    }

    public bool InitializeParametersRandomly {
      get { return InitParametersRandomlyParameter.Value.Value; }
      set { InitParametersRandomlyParameter.Value.Value = value; }
    }

    [StorableConstructor]
    private OdeParameterIdentification(bool deserializing) : base(deserializing) { }
    private OdeParameterIdentification(OdeParameterIdentification original, Cloner cloner)
      : base(original, cloner) {
    }
    public OdeParameterIdentification()
      : base() {
      Problem = new Problem();
      Parameters.Add(new ValueParameter<StringArray>(ModelStructureParameterName, "The function for which the parameters must be fit (only numeric constants are tuned).", new StringArray(new string[] { "1.0 * x*x + 0.0" })));
      Parameters.Add(new FixedValueParameter<IntValue>(IterationsParameterName, "The maximum number of iterations for constants optimization.", new IntValue(200)));
      Parameters.Add(new FixedValueParameter<IntValue>(RestartsParameterName, "The number of independent random restarts (>0)", new IntValue(10)));
      Parameters.Add(new FixedValueParameter<IntValue>(SeedParameterName, "The PRNG seed value.", new IntValue()));
      Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyParameterName, "Switch to determine if the random number seed should be initialized randomly.", new BoolValue(true)));
      Parameters.Add(new FixedValueParameter<BoolValue>(InitParamsRandomlyParameterName, "Switch to determine if the real-valued model parameters should be initialized randomly in each restart.", new BoolValue(false)));

      SetParameterHiddenState();

      InitParametersRandomlyParameter.Value.ValueChanged += (sender, args) => {
        SetParameterHiddenState();
      };
    }

    private void SetParameterHiddenState() {
      var hide = !InitializeParametersRandomly;
      RestartsParameter.Hidden = hide;
      SeedParameter.Hidden = hide;
      SetSeedRandomlyParameter.Hidden = hide;
    }

    [StorableHook(HookType.AfterDeserialization)]
    private void AfterDeserialization() {
    }

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

    #region nonlinear regression
    protected override void Run(CancellationToken cancellationToken) {
      if (SetSeedRandomly) Seed = (new System.Random()).Next();
      var rand = new MersenneTwister((uint)Seed);
      if (InitializeParametersRandomly) {
        throw new NotImplementedException();
        // var qualityTable = new DataTable("RMSE table");
        // qualityTable.VisualProperties.YAxisLogScale = true;
        // var trainRMSERow = new DataRow("RMSE (train)");
        // trainRMSERow.VisualProperties.ChartType = DataRowVisualProperties.DataRowChartType.Points;
        // var testRMSERow = new DataRow("RMSE test");
        // testRMSERow.VisualProperties.ChartType = DataRowVisualProperties.DataRowChartType.Points;
        // 
        // qualityTable.Rows.Add(trainRMSERow);
        // qualityTable.Rows.Add(testRMSERow);
        // Results.Add(new Result(qualityTable.Name, qualityTable.Name + " for all restarts", qualityTable));

        // CreateSolution(Problem, ModelStructure.ToArray(), Iterations, rand);
        // 
        // for (int r = 0; r < Restarts; r++) {
        //   CreateSolution(Problem, ModelStructure.ToArray(), Iterations, rand);
        //   trainRMSERow.Values.Add(solution.TrainingRootMeanSquaredError);
        //   testRMSERow.Values.Add(solution.TestRootMeanSquaredError);
        //   if (solution.TrainingRootMeanSquaredError < bestSolution.TrainingRootMeanSquaredError) {
        //     bestSolution = solution;
        //   }
        // }
      } else {
        CreateSolution(Problem, ModelStructure.ToArray(), Iterations, rand);
      }

      // Results.Add(new Result(RegressionSolutionResultName, "The nonlinear regression solution.", bestSolution));
      // Results.Add(new Result("Root mean square error (train)", "The root of the mean of squared errors of the regression solution on the training set.", new DoubleValue(bestSolution.TrainingRootMeanSquaredError)));
      // Results.Add(new Result("Root mean square error (test)", "The root of the mean of squared errors of the regression solution on the test set.", new DoubleValue(bestSolution.TestRootMeanSquaredError)));

    }

    public void CreateSolution(Problem problem, string[] modelStructure, int maxIterations, IRandom rand) {
      var parser = new InfixExpressionParser();
      var trees = modelStructure.Select(expr => Convert(parser.Parse(expr))).ToArray();
      var names = problem.Encoding.Encodings.Select(enc => enc.Name).ToArray();
      if (trees.Length != names.Length) throw new ArgumentException("The number of expressions must match the number of target variables exactly");

      var scope = new Scope();
      for (int i = 0; i < names.Length; i++) {
        scope.Variables.Add(new Core.Variable(names[i], trees[i]));
      }
      var ind = problem.Encoding.GetIndividual(scope);
      var quality = problem.Evaluate(ind, rand);
      problem.Analyze(new[] { ind }, new[] { quality }, Results, rand);
    }

    private ISymbolicExpressionTree Convert(ISymbolicExpressionTree tree) {
      return new SymbolicExpressionTree(Convert(tree.Root));
    }


    // for translation from symbolic expressions to simple symbols
    private static Dictionary<Type, string> sym2str = new Dictionary<Type, string>() {
      {typeof(Addition), "+" },
      {typeof(Subtraction), "-" },
      {typeof(Multiplication), "*" },
      {typeof(Sine), "sin" },
      {typeof(Cosine), "cos" },
      {typeof(Square), "sqr" },
    };

    private ISymbolicExpressionTreeNode Convert(ISymbolicExpressionTreeNode node) {
      if (sym2str.ContainsKey(node.Symbol.GetType())) {
        var children = node.Subtrees.Select(st => Convert(st)).ToArray();
        return Make(sym2str[node.Symbol.GetType()], children);
      } else if (node.Symbol is ProgramRootSymbol) {
        var child = Convert(node.GetSubtree(0));
        node.RemoveSubtree(0);
        node.AddSubtree(child);
        return node;
      } else if (node.Symbol is StartSymbol) {
        var child = Convert(node.GetSubtree(0));
        node.RemoveSubtree(0);
        node.AddSubtree(child);
        return node;
      } else if (node.Symbol is Division) {
        var children = node.Subtrees.Select(st => Convert(st)).ToArray();
        if (children.Length == 1) {
          return Make("%", new[] { new SimpleSymbol("θ", 0).CreateTreeNode(), children[0] });
        } else if (children.Length != 2) throw new ArgumentException("Division is not supported for multiple arguments");
        else return Make("%", children);
      } else if (node.Symbol is Constant) {
        return new SimpleSymbol("θ", 0).CreateTreeNode();
      } else if (node.Symbol is DataAnalysis.Symbolic.Variable) {
        var varNode = node as VariableTreeNode;
        if (!varNode.Weight.IsAlmost(1.0)) throw new ArgumentException("Variable weights are not supported");
        return new SimpleSymbol(varNode.VariableName, 0).CreateTreeNode();
      } else throw new ArgumentException("Unsupported symbol: " + node.Symbol.Name);
    }

    private ISymbolicExpressionTreeNode Make(string op, ISymbolicExpressionTreeNode[] children) {
      if (children.Length == 1) {
        var s = new SimpleSymbol(op, 1).CreateTreeNode();
        s.AddSubtree(children.First());
        return s;
      } else {
        var s = new SimpleSymbol(op, 2).CreateTreeNode();
        var c0 = children[0];
        var c1 = children[1];
        s.AddSubtree(c0);
        s.AddSubtree(c1);
        for (int i = 2; i < children.Length; i++) {
          var sn = new SimpleSymbol(op, 2).CreateTreeNode();
          sn.AddSubtree(s);
          sn.AddSubtree(children[i]);
          s = sn;
        }
        return s;
      }
    }
    #endregion
  }
}