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

#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;
using HEAL.Attic;

namespace HeuristicLab.Problems.DynamicalSystemsModelling {
  [Item("OdeParameterIdentification", "TODO")]
  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 120)]
  [StorableType("93C151A2-9579-4013-9E67-7611F9378962")]
  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(StorableConstructorFlag _) : base(_) { }
    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 => 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);
    }
    #endregion
  }
}