#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 . */ #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 { 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 ModelStructureParameter { get { return (IValueParameter)Parameters[ModelStructureParameterName]; } } public IFixedValueParameter IterationsParameter { get { return (IFixedValueParameter)Parameters[IterationsParameterName]; } } public IFixedValueParameter SetSeedRandomlyParameter { get { return (IFixedValueParameter)Parameters[SetSeedRandomlyParameterName]; } } public IFixedValueParameter SeedParameter { get { return (IFixedValueParameter)Parameters[SeedParameterName]; } } public IFixedValueParameter RestartsParameter { get { return (IFixedValueParameter)Parameters[RestartsParameterName]; } } public IFixedValueParameter InitParametersRandomlyParameter { get { return (IFixedValueParameter)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(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(IterationsParameterName, "The maximum number of iterations for constants optimization.", new IntValue(200))); Parameters.Add(new FixedValueParameter(RestartsParameterName, "The number of independent random restarts (>0)", new IntValue(10))); Parameters.Add(new FixedValueParameter(SeedParameterName, "The PRNG seed value.", new IntValue())); Parameters.Add(new FixedValueParameter(SetSeedRandomlyParameterName, "Switch to determine if the random number seed should be initialized randomly.", new BoolValue(true))); Parameters.Add(new FixedValueParameter(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) { IRegressionSolution bestSolution = null; 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 sym2str = new Dictionary() { {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 } }