source: branches/3087_Ceres_Integration/HeuristicLab.Tests/HeuristicLab.Problems.DataAnalysis.Symbolic-3.4/SymbolicDataAnalysisExpressionTreeInterpreterTest.cs @ 18006

Last change on this file since 18006 was 18006, checked in by gkronber, 3 months ago

#3087: merged r17784:18004 from trunk to branch to prepare for trunk reintegration (fixed a conflict in CrossValidation.cs)

File size: 37.7 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using System.Collections.Generic;
24using System.Globalization;
25using System.Linq;
26using HeuristicLab.Encodings.SymbolicExpressionTreeEncoding;
27using HeuristicLab.Random;
28using Microsoft.VisualStudio.TestTools.UnitTesting;
29namespace HeuristicLab.Problems.DataAnalysis.Symbolic.Tests {
30
31
32  [TestClass]
33  public class SymbolicDataAnalysisExpressionTreeInterpreterTest {
34    private const int N = 1000;
35    private const int Rows = 1000;
36    private const int Columns = 50;
37
38    private static Dataset ds = new Dataset(new string[] { "Y", "A", "B" }, new double[,] {
39        { 1.0, 1.0, 1.0 },
40        { 2.0, 2.0, 2.0 },
41        { 3.0, 1.0, 2.0 },
42        { 4.0, 1.0, 1.0 },
43        { 5.0, 2.0, 2.0 },
44        { 6.0, 1.0, 2.0 },
45        { 7.0, 1.0, 1.0 },
46        { 8.0, 2.0, 2.0 },
47        { 9.0, 1.0, 2.0 },
48        { 10.0, 1.0, 1.0 },
49        { 11.0, 2.0, 2.0 },
50        { 12.0, 1.0, 2.0 }
51      });
52
53    [TestMethod]
54    [TestCategory("Problems.DataAnalysis.Symbolic")]
55    [TestProperty("Time", "long")]
56    public void StandardInterpreterTestTypeCoherentGrammarPerformance() {
57      TestTypeCoherentGrammarPerformance(new SymbolicDataAnalysisExpressionTreeInterpreter(), 12.5e6);
58    }
59    [TestMethod]
60    [TestCategory("Problems.DataAnalysis.Symbolic")]
61    [TestProperty("Time", "long")]
62    public void StandardInterpreterTestFullGrammarPerformance() {
63      TestFullGrammarPerformance(new SymbolicDataAnalysisExpressionTreeInterpreter(), 12.5e6);
64    }
65    [TestMethod]
66    [TestCategory("Problems.DataAnalysis.Symbolic")]
67    [TestProperty("Time", "long")]
68    public void StandardInterpreterTestArithmeticGrammarPerformance() {
69      TestArithmeticGrammarPerformance(new SymbolicDataAnalysisExpressionTreeInterpreter(), 12.5e6);
70    }
71
72    [TestMethod]
73    [TestCategory("Problems.DataAnalysis.Symbolic")]
74    [TestProperty("Time", "long")]
75    public void CompiledInterpreterTestTypeCoherentGrammarPerformance() {
76      TestTypeCoherentGrammarPerformance(new SymbolicDataAnalysisExpressionCompiledTreeInterpreter(), 12.5e6);
77    }
78    [TestMethod]
79    [TestCategory("Problems.DataAnalysis.Symbolic")]
80    [TestProperty("Time", "long")]
81    public void CompiledInterpreterTestFullGrammarPerformance() {
82      TestFullGrammarPerformance(new SymbolicDataAnalysisExpressionCompiledTreeInterpreter(), 12.5e6);
83    }
84    [TestMethod]
85    [TestCategory("Problems.DataAnalysis.Symbolic")]
86    [TestProperty("Time", "long")]
87    public void CompiledInterpreterTestArithmeticGrammarPerformance() {
88      TestArithmeticGrammarPerformance(new SymbolicDataAnalysisExpressionCompiledTreeInterpreter(), 12.5e6);
89    }
90
91    [TestMethod]
92    [TestCategory("Problems.DataAnalysis.Symbolic")]
93    [TestProperty("Time", "long")]
94    public void ILEmittingInterpreterTestTypeCoherentGrammarPerformance() {
95      TestTypeCoherentGrammarPerformance(new SymbolicDataAnalysisExpressionTreeILEmittingInterpreter(), 7.5e6);
96    }
97    [TestMethod]
98    [TestCategory("Problems.DataAnalysis.Symbolic")]
99    [TestProperty("Time", "long")]
100    public void ILEmittingInterpreterTestArithmeticGrammarPerformance() {
101      TestArithmeticGrammarPerformance(new SymbolicDataAnalysisExpressionTreeILEmittingInterpreter(), 7.5e6);
102    }
103
104    [TestMethod]
105    [TestCategory("Problems.DataAnalysis.Symbolic")]
106    [TestProperty("Time", "long")]
107    public void LinearInterpreterTestTypeCoherentGrammarPerformance() {
108      TestTypeCoherentGrammarPerformance(new SymbolicDataAnalysisExpressionTreeLinearInterpreter(), 12.5e6);
109    }
110    [TestMethod]
111    [TestCategory("Problems.DataAnalysis.Symbolic")]
112    [TestProperty("Time", "long")]
113    public void LinearInterpreterTestFullGrammarPerformance() {
114      TestFullGrammarPerformance(new SymbolicDataAnalysisExpressionTreeLinearInterpreter(), 12.5e6);
115    }
116    [TestMethod]
117    [TestCategory("Problems.DataAnalysis.Symbolic")]
118    [TestProperty("Time", "long")]
119    public void LinearInterpreterTestArithmeticGrammarPerformance() {
120      TestArithmeticGrammarPerformance(new SymbolicDataAnalysisExpressionTreeLinearInterpreter(), 12.5e6);
121    }
122
123    [TestMethod]
124    [TestCategory("Problems.DataAnalysis.Symbolic")]
125    [TestProperty("Time", "long")]
126    public void NativeInterpreterTestTypeCoherentGrammarPerformance() {
127      TestTypeCoherentGrammarPerformance(new NativeInterpreter(), 12.5e6);
128    }
129    [TestMethod]
130    [TestCategory("Problems.DataAnalysis.Symbolic")]
131    [TestProperty("Time", "long")]
132    public void NativeInterpreterTestFullGrammarPerformance() {
133      TestFullGrammarPerformance(new NativeInterpreter(), 12.5e6);
134    }
135    [TestMethod]
136    [TestCategory("Problems.DataAnalysis.Symbolic")]
137    [TestProperty("Time", "long")]
138    public void NativeInterpreterTestArithmeticGrammarPerformance() {
139      TestArithmeticGrammarPerformance(new NativeInterpreter(), 12.5e6);
140    }
141
142    [TestMethod]
143    [TestCategory("Problems.DataAnalysis.Symbolic")]
144    [TestProperty("Time", "long")]
145    public void NativeInterpreterTestCeres() {
146      var parser = new InfixExpressionParser();
147      var random = new FastRandom(1234);
148      const int nRows = 20;
149
150      var x1 = Enumerable.Range(0, nRows).Select(_ => UniformDistributedRandom.NextDouble(random, -1, 1)).ToArray();
151      var x2 = Enumerable.Range(0, nRows).Select(_ => UniformDistributedRandom.NextDouble(random, -1, 1)).ToArray();
152      var x3 = Enumerable.Range(0, nRows).Select(_ => UniformDistributedRandom.NextDouble(random, -1, 1)).ToArray();
153
154      var optimalAlpha = new double[] { -2, -3, -5 };
155      var y = Enumerable.Range(0, nRows).Select(i =>
156          Math.Exp(x1[i] * optimalAlpha[0]) +
157          Math.Exp(x2[i] * optimalAlpha[1]) +
158          Math.Exp(x3[i] * optimalAlpha[2])).ToArray();
159
160      var initialAlpha = Enumerable.Range(0, 3).Select(_ => UniformDistributedRandom.NextDouble(random, -1, 1)).ToArray();
161      var ds = new Dataset(new[] { "x1", "x2", "x3", "y" }, new[] { x1, x2, x3, y });
162
163      var expr = "EXP(x1) + EXP(x2) + EXP(x3)";
164      var tree = parser.Parse(expr);
165      var rows = Enumerable.Range(0, nRows).ToArray();
166      var options = new SolverOptions {
167        Minimizer = (int)MinimizerType.TRUST_REGION,
168        Iterations = 20,
169        TrustRegionStrategy = (int)TrustRegionStrategyType.LEVENBERG_MARQUARDT,
170        LinearSolver = (int)LinearSolverType.DENSE_QR
171      };
172
173      var nodesToOptimize = new HashSet<ISymbolicExpressionTreeNode>(tree.IterateNodesPrefix().Where(x => x is VariableTreeNode));
174      int idx = 0;
175      foreach(var node in nodesToOptimize) {
176        (node as VariableTreeNode).Weight = initialAlpha[idx++];
177        Console.WriteLine((node as VariableTreeNode).Weight);
178
179      }
180
181      var summary = new OptimizationSummary();
182      var parameters = ParameterOptimizer.OptimizeTree(tree, ds, rows, "y", nodesToOptimize, options, ref summary);
183
184      Console.Write("Optimized parameters: ");
185      foreach (var t in parameters) {
186        Console.Write(t.Value + " ");
187      }
188      Console.WriteLine();
189
190      Console.WriteLine("Optimization summary:");
191      Console.WriteLine("Initial cost:         " + summary.InitialCost);
192      Console.WriteLine("Final cost:           " + summary.FinalCost);
193      Console.WriteLine("Successful steps:     " + summary.SuccessfulSteps);
194      Console.WriteLine("Unsuccessful steps:   " + summary.UnsuccessfulSteps);
195      Console.WriteLine("Residual evaluations: " + summary.ResidualEvaluations);
196      Console.WriteLine("Jacobian evaluations: " + summary.JacobianEvaluations);
197    }
198
199    [TestMethod]
200    [TestCategory("Problems.DataAnalysis.Symbolic")]
201    [TestProperty("Time", "long")]
202    public void NativeInterpreterTestCeresVariableProjection() {
203      var parser = new InfixExpressionParser();
204      var random = new FastRandom(1234);
205      const int nRows = 20;
206
207      var x1 = Enumerable.Range(0, nRows).Select(_ => UniformDistributedRandom.NextDouble(random, -1, 1)).ToArray();
208      var x2 = Enumerable.Range(0, nRows).Select(_ => UniformDistributedRandom.NextDouble(random, -1, 1)).ToArray();
209      var x3 = Enumerable.Range(0, nRows).Select(_ => UniformDistributedRandom.NextDouble(random, -1, 1)).ToArray();
210
211      var optimalAlpha = new double[] { -2, -3, -5 };
212      var y = Enumerable.Range(0, nRows).Select(i =>
213        Math.Exp(x1[i] * optimalAlpha[0]) +
214        Math.Exp(x2[i] * optimalAlpha[1]) +
215        Math.Exp(x3[i] * optimalAlpha[2])).ToArray();
216
217      var initialAlpha = Enumerable.Range(0, 3).Select(_ => UniformDistributedRandom.NextDouble(random, -1, 1)).ToArray();
218      var ds = new Dataset(new[] { "x1", "x2", "x3", "y" }, new[] { x1, x2, x3, y });
219
220      var expr = new[] { "EXP(x1)", "EXP(x2)", "EXP(x3)" };
221      var trees = expr.Select(x => parser.Parse(x)).ToArray();
222      var rows = Enumerable.Range(0, nRows).ToArray();
223      var options = new SolverOptions {
224        Minimizer = (int)MinimizerType.TRUST_REGION,
225        Iterations = 100,
226        TrustRegionStrategy = (int)TrustRegionStrategyType.LEVENBERG_MARQUARDT,
227        LinearSolver = (int)LinearSolverType.DENSE_QR
228      };
229
230      var summary = new OptimizationSummary();
231
232      var nodesToOptimize = new HashSet<ISymbolicExpressionTreeNode>(trees.SelectMany(t => t.IterateNodesPrefix().Where(x => x is VariableTreeNode)));
233      int idx = 0;
234      Console.Write("Initial parameters: ");
235      foreach (var node in nodesToOptimize) {
236        (node as VariableTreeNode).Weight = initialAlpha[idx++];
237        Console.Write((node as VariableTreeNode).Weight + " ");
238      }
239      Console.WriteLine();
240
241      var coeff = new double[trees.Length + 1];
242      var parameters = ParameterOptimizer.OptimizeTree(trees, ds, rows, "y", nodesToOptimize, options, coeff, ref summary);
243      Console.Write("Optimized parameters: ");
244      foreach (var t in parameters) {
245        Console.Write(t.Value + " ");
246      }
247      Console.WriteLine();
248
249      Console.Write("Coefficients: ");
250      foreach (var v in coeff) Console.Write(v + " ");
251      Console.WriteLine();
252
253      Console.WriteLine("Optimization summary:");
254      Console.WriteLine("Initial cost:         " + summary.InitialCost);
255      Console.WriteLine("Final cost:           " + summary.FinalCost);
256      Console.WriteLine("Successful steps:     " + summary.SuccessfulSteps);
257      Console.WriteLine("Unsuccessful steps:   " + summary.UnsuccessfulSteps);
258      Console.WriteLine("Residual evaluations: " + summary.ResidualEvaluations);
259      Console.WriteLine("Jacobian evaluations: " + summary.JacobianEvaluations);
260    }
261
262    [TestMethod]
263    [TestCategory("Problems.DataAnalysis.Symbolic")]
264    [TestProperty("Time", "long")]
265    public void BatchInterpreterTestTypeCoherentGrammarPerformance() {
266      TestTypeCoherentGrammarPerformance(new SymbolicDataAnalysisExpressionTreeBatchInterpreter(), 12.5e6);
267    }
268    [TestMethod]
269    [TestCategory("Problems.DataAnalysis.Symbolic")]
270    [TestProperty("Time", "long")]
271    public void BatchInterpreterTestArithmeticGrammarPerformance() {
272      TestArithmeticGrammarPerformance(new SymbolicDataAnalysisExpressionTreeBatchInterpreter(), 12.5e6);
273    }
274
275    private void TestTypeCoherentGrammarPerformance(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter, double nodesPerSecThreshold) {
276      var twister = new MersenneTwister(31415);
277      var dataset = Util.CreateRandomDataset(twister, Rows, Columns);
278
279      var grammar = new TypeCoherentExpressionGrammar();
280      grammar.ConfigureAsDefaultRegressionGrammar();
281
282      var randomTrees = Util.CreateRandomTrees(twister, dataset, grammar, N, 1, 100, 0, 0);
283      foreach (ISymbolicExpressionTree tree in randomTrees) {
284        Util.InitTree(tree, twister, new List<string>(dataset.VariableNames));
285      }
286      double nodesPerSec = Util.CalculateEvaluatedNodesPerSec(randomTrees, interpreter, dataset, 3);
287      //mkommend: commented due to performance issues on the builder
288      // Assert.IsTrue(nodesPerSec > nodesPerSecThreshold); // evaluated nodes per seconds must be larger than 15mNodes/sec
289    }
290
291    private void TestFullGrammarPerformance(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter, double nodesPerSecThreshold) {
292      var twister = new MersenneTwister(31415);
293      var dataset = Util.CreateRandomDataset(twister, Rows, Columns);
294
295      var grammar = new FullFunctionalExpressionGrammar();
296      var randomTrees = Util.CreateRandomTrees(twister, dataset, grammar, N, 1, 100, 0, 0);
297      foreach (ISymbolicExpressionTree tree in randomTrees) {
298        Util.InitTree(tree, twister, new List<string>(dataset.VariableNames));
299      }
300      double nodesPerSec = Util.CalculateEvaluatedNodesPerSec(randomTrees, interpreter, dataset, 3);
301      //mkommend: commented due to performance issues on the builder
302      //Assert.IsTrue(nodesPerSec > nodesPerSecThreshold); // evaluated nodes per seconds must be larger than 15mNodes/sec
303    }
304
305    private void TestArithmeticGrammarPerformance(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter, double nodesPerSecThreshold) {
306      var twister = new MersenneTwister(31415);
307      var dataset = Util.CreateRandomDataset(twister, Rows, Columns);
308
309      var grammar = new ArithmeticExpressionGrammar();
310      var randomTrees = Util.CreateRandomTrees(twister, dataset, grammar, N, 1, 100, 0, 0);
311      foreach (SymbolicExpressionTree tree in randomTrees) {
312        Util.InitTree(tree, twister, new List<string>(dataset.VariableNames));
313      }
314
315      double nodesPerSec = Util.CalculateEvaluatedNodesPerSec(randomTrees, interpreter, dataset, 3);
316      //mkommend: commented due to performance issues on the builder
317      //Assert.IsTrue(nodesPerSec > nodesPerSecThreshold); // evaluated nodes per seconds must be larger than 15mNodes/sec
318    }
319
320
321    /// <summary>
322    ///A test for Evaluate
323    ///</summary>
324    [TestMethod]
325    [TestCategory("Problems.DataAnalysis.Symbolic")]
326    [TestProperty("Time", "short")]
327    public void StandardInterpreterTestEvaluation() {
328      var interpreter = new SymbolicDataAnalysisExpressionTreeInterpreter();
329      EvaluateTerminals(interpreter, ds);
330      EvaluateOperations(interpreter, ds);
331      EvaluateLaggedOperations(interpreter, ds);
332      EvaluateSpecialFunctions(interpreter, ds);
333      EvaluateAdf(interpreter, ds);
334    }
335
336    /// <summary>
337    ///A test for Evaluate
338    ///</summary>
339    [TestMethod]
340    [TestCategory("Problems.DataAnalysis.Symbolic")]
341    [TestProperty("Time", "short")]
342    public void ILEmittingInterpreterTestEvaluation() {
343      var interpreter = new SymbolicDataAnalysisExpressionTreeILEmittingInterpreter();
344      EvaluateTerminals(interpreter, ds);
345      EvaluateOperations(interpreter, ds);
346      EvaluateLaggedOperations(interpreter, ds);
347      EvaluateSpecialFunctions(interpreter, ds);
348    }
349
350    [TestMethod]
351    [TestCategory("Problems.DataAnalysis.Symbolic")]
352    [TestProperty("Time", "short")]
353    public void CompiledInterpreterTestEvaluation() {
354      var interpreter = new SymbolicDataAnalysisExpressionCompiledTreeInterpreter();
355      EvaluateTerminals(interpreter, ds);
356      EvaluateOperations(interpreter, ds);
357      EvaluateSpecialFunctions(interpreter, ds);
358    }
359
360    [TestMethod]
361    [TestCategory("Problems.DataAnalysis.Symbolic")]
362    [TestProperty("Time", "short")]
363    public void LinearInterpreterTestEvaluation() {
364      var interpreter = new SymbolicDataAnalysisExpressionTreeLinearInterpreter();
365      //ADFs are not supported by the linear interpreter
366      EvaluateTerminals(interpreter, ds);
367      EvaluateOperations(interpreter, ds);
368      EvaluateLaggedOperations(interpreter, ds);
369      EvaluateSpecialFunctions(interpreter, ds);
370    }
371
372    [TestMethod]
373    [TestCategory("Problems.DataAnalysis.Symbolic")]
374    [TestProperty("Time", "long")]
375    public void TestInterpretersEstimatedValuesConsistency() {
376      var twister = new MersenneTwister();
377      int seed = twister.Next(0, int.MaxValue);
378      twister.Seed((uint)seed);
379      const int numRows = 100;
380      var dataset = Util.CreateRandomDataset(twister, numRows, Columns);
381
382      var grammar = new TypeCoherentExpressionGrammar();
383
384      var interpreters = new ISymbolicDataAnalysisExpressionTreeInterpreter[] {
385        new SymbolicDataAnalysisExpressionTreeLinearInterpreter(),
386        new SymbolicDataAnalysisExpressionTreeInterpreter(),
387      };
388
389      var rows = Enumerable.Range(0, numRows).ToList();
390      var randomTrees = Util.CreateRandomTrees(twister, dataset, grammar, N, 1, 10, 0, 0);
391      foreach (ISymbolicExpressionTree tree in randomTrees) {
392        Util.InitTree(tree, twister, new List<string>(dataset.VariableNames));
393      }
394
395      for (int i = 0; i < randomTrees.Length; ++i) {
396        var tree = randomTrees[i];
397        var valuesMatrix = interpreters.Select(x => x.GetSymbolicExpressionTreeValues(tree, dataset, rows)).ToList();
398        for (int m = 0; m < interpreters.Length - 1; ++m) {
399          var sum = valuesMatrix[m].Sum();
400          for (int n = m + 1; n < interpreters.Length; ++n) {
401            var s = valuesMatrix[n].Sum();
402            if (double.IsNaN(sum) && double.IsNaN(s)) continue;
403
404            string errorMessage = string.Format("Interpreters {0} and {1} do not agree on tree {2} (seed = {3}).", interpreters[m].Name, interpreters[n].Name, i, seed);
405            Assert.AreEqual(sum, s, 1e-12, errorMessage);
406          }
407        }
408      }
409    }
410
411    [TestMethod]
412    [TestCategory("Problems.DataAnalysis.Symbolic")]
413    [TestProperty("Time", "long")]
414    public void TestCompiledInterpreterEstimatedValuesConsistency() {
415      const double delta = 1e-8;
416
417      var twister = new MersenneTwister();
418      int seed = twister.Next(0, int.MaxValue);
419      twister.Seed((uint)seed);
420
421      Console.WriteLine(seed);
422
423      const int numRows = 100;
424      var dataset = Util.CreateRandomDataset(twister, numRows, Columns);
425
426      var grammar = new TypeCoherentExpressionGrammar();
427      grammar.ConfigureAsDefaultRegressionGrammar();
428      grammar.Symbols.First(x => x.Name == "Power Functions").Enabled = true;
429      grammar.Symbols.First(x => x is Cube).Enabled = true;
430      grammar.Symbols.First(x => x is CubeRoot).Enabled = true;
431      grammar.Symbols.First(x => x is Square).Enabled = true;
432      grammar.Symbols.First(x => x is SquareRoot).Enabled = true;
433      grammar.Symbols.First(x => x is Absolute).Enabled = true;
434      grammar.Symbols.First(x => x is Sine).Enabled = true;
435      grammar.Symbols.First(x => x is Cosine).Enabled = true;
436      grammar.Symbols.First(x => x is Tangent).Enabled = true;
437      grammar.Symbols.First(x => x is Root).Enabled = false;
438      grammar.Symbols.First(x => x is Power).Enabled = false;
439
440      var randomTrees = Util.CreateRandomTrees(twister, dataset, grammar, N, 1, 10, 0, 0);
441      foreach (ISymbolicExpressionTree tree in randomTrees) {
442        Util.InitTree(tree, twister, new List<string>(dataset.VariableNames));
443      }
444
445      var interpreters = new ISymbolicDataAnalysisExpressionTreeInterpreter[] {
446        new SymbolicDataAnalysisExpressionCompiledTreeInterpreter(),
447        new SymbolicDataAnalysisExpressionTreeInterpreter(),
448        new SymbolicDataAnalysisExpressionTreeLinearInterpreter(),
449      };
450      var rows = Enumerable.Range(0, numRows).ToList();
451      var formatter = new SymbolicExpressionTreeHierarchicalFormatter();
452
453      for (int i = 0; i < randomTrees.Length; ++i) {
454        var tree = randomTrees[i];
455        var valuesMatrix = interpreters.Select(x => x.GetSymbolicExpressionTreeValues(tree, dataset, rows).ToList()).ToList();
456        for (int m = 0; m < interpreters.Length - 1; ++m) {
457          for (int n = m + 1; n < interpreters.Length; ++n) {
458            for (int row = 0; row < numRows; ++row) {
459              var v1 = valuesMatrix[m][row];
460              var v2 = valuesMatrix[n][row];
461              if (double.IsNaN(v1) && double.IsNaN(v2)) continue;
462              if (v1 != v2 && Math.Abs(1.0 - v1 / v2) >= delta) {
463                Console.WriteLine(formatter.Format(tree));
464                foreach (var node in tree.Root.GetSubtree(0).GetSubtree(0).IterateNodesPrefix().ToList()) {
465                  var rootNode = (SymbolicExpressionTreeTopLevelNode)grammar.ProgramRootSymbol.CreateTreeNode();
466                  if (rootNode.HasLocalParameters) rootNode.ResetLocalParameters(twister);
467                  rootNode.SetGrammar(grammar.CreateExpressionTreeGrammar());
468
469                  var startNode = (SymbolicExpressionTreeTopLevelNode)grammar.StartSymbol.CreateTreeNode();
470                  if (startNode.HasLocalParameters) startNode.ResetLocalParameters(twister);
471                  startNode.SetGrammar(grammar.CreateExpressionTreeGrammar());
472
473                  rootNode.AddSubtree(startNode);
474                  var t = new SymbolicExpressionTree(rootNode);
475                  var start = t.Root.GetSubtree(0);
476                  var p = node.Parent;
477                  start.AddSubtree(node);
478                  Console.WriteLine(node);
479
480                  var y1 = interpreters[m].GetSymbolicExpressionTreeValues(t, dataset, new[] { row }).First();
481                  var y2 = interpreters[n].GetSymbolicExpressionTreeValues(t, dataset, new[] { row }).First();
482
483                  if (double.IsNaN(y1) && double.IsNaN(y2)) continue;
484                  string prefix = Math.Abs(y1 - y2) > delta ? "++" : "==";
485                  Console.WriteLine("\t{0} Row {1}: {2} {3}, Deviation = {4}", prefix, row, y1, y2, Math.Abs(y1 - y2));
486                  node.Parent = p;
487                }
488              }
489              string errorMessage = string.Format("Interpreters {0} and {1} do not agree on tree {2} and row {3} (seed = {4}).", interpreters[m].Name, interpreters[n].Name, i, row, seed);
490              Assert.IsTrue(double.IsNaN(v1) && double.IsNaN(v2) ||
491                            v1 == v2 || // in particular 0 = 0
492                            Math.Abs(1.0 - v1 / v2) < delta, errorMessage);
493            }
494          }
495        }
496      }
497    }
498
499    private void EvaluateTerminals(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter, IDataset ds) {
500      // constants
501      Evaluate(interpreter, ds, "(+ 1.5 3.5)", 0, 5.0);
502
503      // variables
504      Evaluate(interpreter, ds, "(variable 2.0 a)", 0, 2.0);
505      Evaluate(interpreter, ds, "(variable 2.0 a)", 1, 4.0);
506    }
507
508    private void EvaluateAdf(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter, IDataset ds) {
509
510      // ADF     
511      Evaluate(interpreter, ds, @"(PROG
512                                    (MAIN
513                                      (CALL ADF0))
514                                    (defun ADF0 1.0))", 1, 1.0);
515      Evaluate(interpreter, ds, @"(PROG
516                                    (MAIN
517                                      (* (CALL ADF0) (CALL ADF0)))
518                                    (defun ADF0 2.0))", 1, 4.0);
519      Evaluate(interpreter, ds, @"(PROG
520                                    (MAIN
521                                      (CALL ADF0 2.0 3.0))
522                                    (defun ADF0
523                                      (+ (ARG 0) (ARG 1))))", 1, 5.0);
524      Evaluate(interpreter, ds, @"(PROG
525                                    (MAIN (CALL ADF1 2.0 3.0))
526                                    (defun ADF0
527                                      (- (ARG 1) (ARG 0)))
528                                    (defun ADF1
529                                      (+ (CALL ADF0 (ARG 1) (ARG 0))
530                                         (CALL ADF0 (ARG 0) (ARG 1)))))", 1, 0.0);
531      Evaluate(interpreter, ds, @"(PROG
532                                    (MAIN (CALL ADF1 (variable 2.0 a) 3.0))
533                                    (defun ADF0
534                                      (- (ARG 1) (ARG 0)))
535                                    (defun ADF1                                                                             
536                                      (CALL ADF0 (ARG 1) (ARG 0))))", 1, 1.0);
537      Evaluate(interpreter, ds,
538               @"(PROG
539                                    (MAIN (CALL ADF1 (variable 2.0 a) 3.0))
540                                    (defun ADF0
541                                      (- (ARG 1) (ARG 0)))
542                                    (defun ADF1                                                                             
543                                      (+ (CALL ADF0 (ARG 1) (ARG 0))
544                                         (CALL ADF0 (ARG 0) (ARG 1)))))", 1, 0.0);
545    }
546
547    private void EvaluateSpecialFunctions(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter, IDataset ds) {
548      // special functions
549      Action<double> checkAiry = (x) => {
550        double ai, aip, bi, bip;
551        alglib.airy(x, out ai, out aip, out bi, out bip);
552        Evaluate(interpreter, ds, "(airya " + x + ")", 0, ai);
553        Evaluate(interpreter, ds, "(airyb " + x + ")", 0, bi);
554      };
555
556      Action<double> checkBessel = (x) => {
557        Evaluate(interpreter, ds, "(bessel " + x + ")", 0, alglib.besseli0(x));
558      };
559
560      Action<double> checkSinCosIntegrals = (x) => {
561        double si, ci;
562        alglib.sinecosineintegrals(x, out si, out ci);
563        Evaluate(interpreter, ds, "(cosint " + x + ")", 0, ci);
564        Evaluate(interpreter, ds, "(sinint " + x + ")", 0, si);
565      };
566      Action<double> checkHypSinCosIntegrals = (x) => {
567        double shi, chi;
568        alglib.hyperbolicsinecosineintegrals(x, out shi, out chi);
569        Evaluate(interpreter, ds, "(hypcosint " + x + ")", 0, chi);
570        Evaluate(interpreter, ds, "(hypsinint " + x + ")", 0, shi);
571      };
572      Action<double> checkFresnelSinCosIntegrals = (x) => {
573        double c = 0, s = 0;
574        alglib.fresnelintegral(x, ref c, ref s);
575        Evaluate(interpreter, ds, "(fresnelcosint " + x + ")", 0, c);
576        Evaluate(interpreter, ds, "(fresnelsinint " + x + ")", 0, s);
577      };
578      Action<double> checkNormErf = (x) => {
579        Evaluate(interpreter, ds, "(norm " + x + ")", 0, alglib.normaldistribution(x));
580        Evaluate(interpreter, ds, "(erf " + x + ")", 0, alglib.errorfunction(x));
581      };
582
583      Action<double> checkGamma = (x) => {
584        Evaluate(interpreter, ds, "(gamma " + x + ")", 0, alglib.gammafunction(x));
585      };
586      Action<double> checkPsi = (x) => {
587        try {
588          Evaluate(interpreter, ds, "(psi " + x + ")", 0, alglib.psi(x));
589        } catch (alglib.alglibexception) { // ignore cases where alglib throws an exception
590        }
591      };
592      Action<double> checkDawson = (x) => {
593        Evaluate(interpreter, ds, "(dawson " + x + ")", 0, alglib.dawsonintegral(x));
594      };
595      Action<double> checkExpInt = (x) => {
596        Evaluate(interpreter, ds, "(expint " + x + ")", 0, alglib.exponentialintegralei(x));
597      };
598
599      foreach (var e in new[] { -2.0, -1.0, 0.0, 1.0, 2.0 }) {
600        checkAiry(e);
601        checkBessel(e);
602        checkSinCosIntegrals(e);
603        checkGamma(e);
604        checkExpInt(e);
605        checkDawson(e);
606        checkPsi(e);
607        checkNormErf(e);
608        checkFresnelSinCosIntegrals(e);
609        checkHypSinCosIntegrals(e);
610      }
611    }
612
613    private void EvaluateLaggedOperations(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter, IDataset ds) {
614      // lag
615      Evaluate(interpreter, ds, "(lagVariable 1.0 a -1) ", 1, ds.GetDoubleValue("A", 0));
616      Evaluate(interpreter, ds, "(lagVariable 1.0 a -1) ", 2, ds.GetDoubleValue("A", 1));
617      Evaluate(interpreter, ds, "(lagVariable 1.0 a 0) ", 2, ds.GetDoubleValue("A", 2));
618      Evaluate(interpreter, ds, "(lagVariable 1.0 a 1) ", 0, ds.GetDoubleValue("A", 1));
619
620      // integral
621      Evaluate(interpreter, ds, "(integral -1.0 (variable 1.0 a)) ", 1, ds.GetDoubleValue("A", 0) + ds.GetDoubleValue("A", 1));
622      Evaluate(interpreter, ds, "(integral -1.0 (lagVariable 1.0 a 1)) ", 1, ds.GetDoubleValue("A", 1) + ds.GetDoubleValue("A", 2));
623      Evaluate(interpreter, ds, "(integral -2.0 (variable 1.0 a)) ", 2, ds.GetDoubleValue("A", 0) + ds.GetDoubleValue("A", 1) + ds.GetDoubleValue("A", 2));
624      Evaluate(interpreter, ds, "(integral -1.0 (* (variable 1.0 a) (variable 1.0 b)))", 1, ds.GetDoubleValue("A", 0) * ds.GetDoubleValue("B", 0) + ds.GetDoubleValue("A", 1) * ds.GetDoubleValue("B", 1));
625      Evaluate(interpreter, ds, "(integral -2.0 3.0)", 1, 9.0);
626
627      // derivative
628      // (f_0 + 2 * f_1 - 2 * f_3 - f_4) / 8; // h = 1
629      Evaluate(interpreter, ds, "(diff (variable 1.0 a)) ", 5, (ds.GetDoubleValue("A", 5) + 2 * ds.GetDoubleValue("A", 4) - 2 * ds.GetDoubleValue("A", 2) - ds.GetDoubleValue("A", 1)) / 8.0);
630      Evaluate(interpreter, ds, "(diff (variable 1.0 b)) ", 5, (ds.GetDoubleValue("B", 5) + 2 * ds.GetDoubleValue("B", 4) - 2 * ds.GetDoubleValue("B", 2) - ds.GetDoubleValue("B", 1)) / 8.0);
631      Evaluate(interpreter, ds, "(diff (* (variable 1.0 a) (variable 1.0 b)))", 5, +
632        (ds.GetDoubleValue("A", 5) * ds.GetDoubleValue("B", 5) +
633        2 * ds.GetDoubleValue("A", 4) * ds.GetDoubleValue("B", 4) -
634        2 * ds.GetDoubleValue("A", 2) * ds.GetDoubleValue("B", 2) -
635        ds.GetDoubleValue("A", 1) * ds.GetDoubleValue("B", 1)) / 8.0);
636      Evaluate(interpreter, ds, "(diff -2.0 3.0)", 5, 0.0);
637
638      // timelag
639      Evaluate(interpreter, ds, "(lag -1.0 (lagVariable 1.0 a 2)) ", 1, ds.GetDoubleValue("A", 2));
640      Evaluate(interpreter, ds, "(lag -2.0 (lagVariable 1.0 a 2)) ", 2, ds.GetDoubleValue("A", 2));
641      Evaluate(interpreter, ds, "(lag -1.0 (* (lagVariable 1.0 a 1) (lagVariable 1.0 b 2)))", 1, ds.GetDoubleValue("A", 1) * ds.GetDoubleValue("B", 2));
642      Evaluate(interpreter, ds, "(lag -2.0 3.0)", 1, 3.0);
643    }
644
645    private void EvaluateOperations(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter, IDataset ds) {
646      // addition
647      Evaluate(interpreter, ds, "(+ (variable 2.0 a ))", 1, 4.0);
648      Evaluate(interpreter, ds, "(+ (variable 2.0 a ) (variable 3.0 b ))", 0, 5.0);
649      Evaluate(interpreter, ds, "(+ (variable 2.0 a ) (variable 3.0 b ))", 1, 10.0);
650      Evaluate(interpreter, ds, "(+ (variable 2.0 a) (variable 3.0 b ))", 2, 8.0);
651      Evaluate(interpreter, ds, "(+ 8.0 2.0 2.0)", 0, 12.0);
652
653      // subtraction
654      Evaluate(interpreter, ds, "(- (variable 2.0 a ))", 1, -4.0);
655      Evaluate(interpreter, ds, "(- (variable 2.0 a ) (variable 3.0 b))", 0, -1.0);
656      Evaluate(interpreter, ds, "(- (variable 2.0 a ) (variable 3.0 b ))", 1, -2.0);
657      Evaluate(interpreter, ds, "(- (variable 2.0 a ) (variable 3.0 b ))", 2, -4.0);
658      Evaluate(interpreter, ds, "(- 8.0 2.0 2.0)", 0, 4.0);
659
660      // multiplication
661      Evaluate(interpreter, ds, "(* (variable 2.0 a ))", 0, 2.0);
662      Evaluate(interpreter, ds, "(* (variable 2.0 a ) (variable 3.0 b ))", 0, 6.0);
663      Evaluate(interpreter, ds, "(* (variable 2.0 a ) (variable 3.0 b ))", 1, 24.0);
664      Evaluate(interpreter, ds, "(* (variable 2.0 a ) (variable 3.0 b ))", 2, 12.0);
665      Evaluate(interpreter, ds, "(* 8.0 2.0 2.0)", 0, 32.0);
666
667      // division
668      Evaluate(interpreter, ds, "(/ (variable 2.0 a ))", 1, 1.0 / 4.0);
669      Evaluate(interpreter, ds, "(/ (variable 2.0 a ) 2.0)", 0, 1.0);
670      Evaluate(interpreter, ds, "(/ (variable 2.0 a ) 2.0)", 1, 2.0);
671      Evaluate(interpreter, ds, "(/ (variable 3.0 b ) 2.0)", 2, 3.0);
672      Evaluate(interpreter, ds, "(/ 8.0 2.0 2.0)", 0, 2.0);
673
674      // gt
675      Evaluate(interpreter, ds, "(> (variable 2.0 a) 2.0)", 0, -1.0);
676      Evaluate(interpreter, ds, "(> 2.0 (variable 2.0 a))", 0, -1.0);
677      Evaluate(interpreter, ds, "(> (variable 2.0 a) 1.9)", 0, 1.0);
678      Evaluate(interpreter, ds, "(> 1.9 (variable 2.0 a))", 0, -1.0);
679      Evaluate(interpreter, ds, "(> (log -1.0) (log -1.0))", 0, -1.0); // (> nan nan) should be false
680
681      // lt
682      Evaluate(interpreter, ds, "(< (variable 2.0 a) 2.0)", 0, -1.0);
683      Evaluate(interpreter, ds, "(< 2.0 (variable 2.0 a))", 0, -1.0);
684      Evaluate(interpreter, ds, "(< (variable 2.0 a) 1.9)", 0, -1.0);
685      Evaluate(interpreter, ds, "(< 1.9 (variable 2.0 a))", 0, 1.0);
686      Evaluate(interpreter, ds, "(< (log -1.0) (log -1.0))", 0, -1.0); // (< nan nan) should be false
687
688      // If
689      Evaluate(interpreter, ds, "(if -10.0 2.0 3.0)", 0, 3.0);
690      Evaluate(interpreter, ds, "(if -1.0 2.0 3.0)", 0, 3.0);
691      Evaluate(interpreter, ds, "(if 0.0 2.0 3.0)", 0, 3.0);
692      Evaluate(interpreter, ds, "(if 1.0 2.0 3.0)", 0, 2.0);
693      Evaluate(interpreter, ds, "(if 10.0 2.0 3.0)", 0, 2.0);
694      Evaluate(interpreter, ds, "(if (log -1.0) 2.0 3.0)", 0, 3.0); // if(nan) should return the else branch
695
696      // NOT
697      Evaluate(interpreter, ds, "(not -1.0)", 0, 1.0);
698      Evaluate(interpreter, ds, "(not -2.0)", 0, 1.0);
699      Evaluate(interpreter, ds, "(not 1.0)", 0, -1.0);
700      Evaluate(interpreter, ds, "(not 2.0)", 0, -1.0);
701      Evaluate(interpreter, ds, "(not 0.0)", 0, 1.0);
702      Evaluate(interpreter, ds, "(not (log -1.0))", 0, 1.0);
703
704      // AND
705      Evaluate(interpreter, ds, "(and -1.0 -2.0)", 0, -1.0);
706      Evaluate(interpreter, ds, "(and -1.0 2.0)", 0, -1.0);
707      Evaluate(interpreter, ds, "(and 1.0 -2.0)", 0, -1.0);
708      Evaluate(interpreter, ds, "(and 1.0 0.0)", 0, -1.0);
709      Evaluate(interpreter, ds, "(and 0.0 0.0)", 0, -1.0);
710      Evaluate(interpreter, ds, "(and 1.0 2.0)", 0, 1.0);
711      Evaluate(interpreter, ds, "(and 1.0 2.0 3.0)", 0, 1.0);
712      Evaluate(interpreter, ds, "(and 1.0 -2.0 3.0)", 0, -1.0);
713      Evaluate(interpreter, ds, "(and (log -1.0))", 0, -1.0); // (and NaN)
714      Evaluate(interpreter, ds, "(and (log -1.0)  1.0)", 0, -1.0); // (and NaN 1.0)
715
716      // OR
717      Evaluate(interpreter, ds, "(or -1.0 -2.0)", 0, -1.0);
718      Evaluate(interpreter, ds, "(or -1.0 2.0)", 0, 1.0);
719      Evaluate(interpreter, ds, "(or 1.0 -2.0)", 0, 1.0);
720      Evaluate(interpreter, ds, "(or 1.0 2.0)", 0, 1.0);
721      Evaluate(interpreter, ds, "(or 0.0 0.0)", 0, -1.0);
722      Evaluate(interpreter, ds, "(or -1.0 -2.0 -3.0)", 0, -1.0);
723      Evaluate(interpreter, ds, "(or -1.0 -2.0 3.0)", 0, 1.0);
724      Evaluate(interpreter, ds, "(or (log -1.0))", 0, -1.0); // (or NaN)
725      Evaluate(interpreter, ds, "(or (log -1.0)  1.0)", 0, -1.0); // (or NaN 1.0)
726
727      // XOR
728      Evaluate(interpreter, ds, "(xor -1.0 -2.0)", 0, -1.0);
729      Evaluate(interpreter, ds, "(xor -1.0 2.0)", 0, 1.0);
730      Evaluate(interpreter, ds, "(xor 1.0 -2.0)", 0, 1.0);
731      Evaluate(interpreter, ds, "(xor 1.0 2.0)", 0, -1.0);
732      Evaluate(interpreter, ds, "(xor 0.0 0.0)", 0, -1.0);
733      Evaluate(interpreter, ds, "(xor -1.0 -2.0 -3.0)", 0, -1.0);
734      Evaluate(interpreter, ds, "(xor -1.0 -2.0 3.0)", 0, 1.0);
735      Evaluate(interpreter, ds, "(xor -1.0 2.0 3.0)", 0, -1.0);
736      Evaluate(interpreter, ds, "(xor 1.0 2.0 3.0)", 0, 1.0);
737      Evaluate(interpreter, ds, "(xor (log -1.0))", 0, -1.0);
738      Evaluate(interpreter, ds, "(xor (log -1.0)  1.0)", 0, 1.0);
739
740      // sin, cos, tan
741      Evaluate(interpreter, ds, "(sin " + Math.PI.ToString(NumberFormatInfo.InvariantInfo) + ")", 0, 0.0);
742      Evaluate(interpreter, ds, "(sin 0.0)", 0, 0.0);
743      Evaluate(interpreter, ds, "(cos " + Math.PI.ToString(NumberFormatInfo.InvariantInfo) + ")", 0, -1.0);
744      Evaluate(interpreter, ds, "(cos 0.0)", 0, 1.0);
745      Evaluate(interpreter, ds, "(tan " + Math.PI.ToString(NumberFormatInfo.InvariantInfo) + ")", 0, Math.Tan(Math.PI));
746      Evaluate(interpreter, ds, "(tan 0.0)", 0, Math.Tan(Math.PI));
747
748      // exp, log
749      Evaluate(interpreter, ds, "(log (exp 7.0))", 0, Math.Log(Math.Exp(7)));
750      Evaluate(interpreter, ds, "(exp (log 7.0))", 0, Math.Exp(Math.Log(7)));
751      Evaluate(interpreter, ds, "(log -3.0)", 0, Math.Log(-3));
752
753      // power
754      Evaluate(interpreter, ds, "(pow 2.0 3.0)", 0, 8.0);
755      Evaluate(interpreter, ds, "(pow 4.0 0.5)", 0, 1.0); // interpreter should round to the nearest integer value value (.5 is rounded to the even number)
756      Evaluate(interpreter, ds, "(pow 4.0 2.5)", 0, 16.0); // interpreter should round to the nearest integer value value (.5 is rounded to the even number)
757      Evaluate(interpreter, ds, "(pow -2.0 3.0)", 0, -8.0);
758      Evaluate(interpreter, ds, "(pow 2.0 -3.0)", 0, 1.0 / 8.0);
759      Evaluate(interpreter, ds, "(pow -2.0 -3.0)", 0, -1.0 / 8.0);
760
761      // root
762      Evaluate(interpreter, ds, "(root 9.0 2.0)", 0, 3.0);
763      Evaluate(interpreter, ds, "(root 27.0 3.0)", 0, 3.0);
764      Evaluate(interpreter, ds, "(root 2.0 -3.0)", 0, Math.Pow(2.0, -1.0 / 3.0));
765
766      // mean
767      Evaluate(interpreter, ds, "(mean -1.0 1.0 -1.0)", 0, -1.0 / 3.0);
768    }
769
770    private void Evaluate(ISymbolicDataAnalysisExpressionTreeInterpreter interpreter, IDataset ds, string expr, int index, double expected) {
771      var importer = new SymbolicExpressionImporter();
772      ISymbolicExpressionTree tree = importer.Import(expr);
773
774      double actual = interpreter.GetSymbolicExpressionTreeValues(tree, ds, Enumerable.Range(index, 1)).First();
775
776      Assert.IsFalse(double.IsNaN(actual) && !double.IsNaN(expected));
777      Assert.IsFalse(!double.IsNaN(actual) && double.IsNaN(expected));
778      if (!double.IsNaN(actual) && !double.IsNaN(expected))
779        Assert.AreEqual(expected, actual, 1.0E-12, expr);
780    }
781  }
782}
Note: See TracBrowser for help on using the repository browser.