Free cookie consent management tool by TermsFeed Policy Generator

source: branches/3106_AnalyticContinuedFractionsRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/ContinuedFractionRegression/Problem.cs @ 17849

Last change on this file since 17849 was 17849, checked in by gkronber, 3 years ago

#3106 register event handlers in storable hook

File size: 42.1 KB
Line 
1using System;
2using System.Collections.Generic;
3using System.Linq;
4using HEAL.Attic;
5using HeuristicLab.Analysis;
6using HeuristicLab.Common;
7using HeuristicLab.Core;
8using HeuristicLab.Data;
9using HeuristicLab.Encodings.BinaryVectorEncoding;
10using HeuristicLab.Encodings.RealVectorEncoding;
11using HeuristicLab.Optimization;
12using HeuristicLab.Parameters;
13using HeuristicLab.Problems.DataAnalysis;
14using HeuristicLab.Problems.Instances;
15
16namespace HeuristicLab.Algorithms.DataAnalysis.ContinuedFractionRegression {
17  [Item("Continued Fraction Regression (CFR)", "Attempts to find a continued fraction minimizing the regularized MSE / R^2 to given data.")]
18  [StorableType("CAC0F743-8524-436C-B81E-7C628A302DF8")]
19  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 999)]
20  public class Problem : HeuristicLab.Optimization.SingleObjectiveBasicProblem<MultiEncoding>
21    /*, IProblemInstanceConsumer<IRegressionProblemData>, IProblemInstanceExporter<IRegressionProblemData>*/  // only if we can change the code to work with a single dataset
22    {
23    private const double epsilon = 1e-6;
24    enum Fitting { MSE, RegMSE, ElasticNet }
25    private readonly Fitting fit = Fitting.ElasticNet;
26
27    #region parameter properties
28    public IFixedValueParameter<DoubleValue> L1L2MixingParameter {
29      get { return (IFixedValueParameter<DoubleValue>)Parameters["L1L2_mixing"]; }
30    }
31    public IFixedValueParameter<DoubleValue> L1L2WeightParameter {
32      get { return (IFixedValueParameter<DoubleValue>)Parameters["L1L2_weight"]; }
33    }
34    public IFixedValueParameter<DoubleValue> PenaltyParameter {
35      get { return (IFixedValueParameter<DoubleValue>)Parameters["penalty"]; }
36    }
37    public IFixedValueParameter<IntValue> FractionDepthParameter {
38      get { return (IFixedValueParameter<IntValue>)Parameters["fraction_depth"]; }
39    }
40    public IValueParameter<Dataset> DatasetParameter {
41      get { return (IValueParameter<Dataset>)Parameters["dataset"]; }
42    }
43    public IValueParameter<Dataset> TrainingTrainingDatasetParameter {
44      get { return (IValueParameter<Dataset>)Parameters["datasetTrainingTraining"]; }
45    }
46    public IValueParameter<Dataset> TrainingTestDatasetParameter {
47      get { return (IValueParameter<Dataset>)Parameters["datasetTrainingTest"]; }
48    }
49    public IValueParameter<Dataset> Phi0_08DatasetParameter {
50      get { return (IValueParameter<Dataset>)Parameters["data_phi_0_08"]; }
51    }
52    #endregion
53
54    #region properties
55    public double L1L2Mixing {
56      get { return L1L2MixingParameter.Value.Value; }
57      set { L1L2MixingParameter.Value.Value = value; }
58    }
59    public double L1L2Weight {
60      get { return L1L2WeightParameter.Value.Value; }
61      set { L1L2WeightParameter.Value.Value = value; }
62    }
63    public double Penalty {
64      get { return PenaltyParameter.Value.Value; }
65      set { PenaltyParameter.Value.Value = value; }
66    }
67    public int FractionDepth {
68      get { return FractionDepthParameter.Value.Value; }
69      set { FractionDepthParameter.Value.Value = value; }
70    }
71    public Dataset Dataset {
72      get { return DatasetParameter.Value; }
73      set { DatasetParameter.Value = value; }
74    }
75    public Dataset TrainingTrainingDataset {
76      get { return TrainingTrainingDatasetParameter.Value; }
77      set { TrainingTrainingDatasetParameter.Value = value; }
78    }
79    public Dataset TrainingTestDataset {
80      get { return TrainingTestDatasetParameter.Value; }
81      set { TrainingTestDatasetParameter.Value = value; }
82    }
83    public Dataset Phi0_08Dataset {
84      get { return Phi0_08DatasetParameter.Value; }
85      set { Phi0_08DatasetParameter.Value = value; }
86    }
87    #endregion
88
89
90    public override bool Maximization { get { return false; } }
91
92    public IRegressionProblemData ProblemData { get; set; }
93
94    // BEWARE local fields
95    // DO NOT store or clone
96    // these variables are initialized on the first call of Evaluate.
97    // whenever a dataset parameter is changed (which has an effect on these variables)
98    // all variables must be cleared to ensure that they are recalculated on the next call of Evaluate
99    private DoubleMatrix dataMatrix; // data for training
100    private DoubleMatrix dataNelderMead; // data for optimizing coefficients
101    private DoubleMatrix dataMatrixTest; // data for testing
102    private Transformation transformation;
103    private double arithmeticMean;
104    private double SQT; //Sum of Squares Total, short: SQT
105    private void ClearLocalVariables() {
106      dataMatrix = null;
107      dataNelderMead = null;
108      dataMatrixTest = null;
109      transformation = null;
110      arithmeticMean = 0.0;
111      SQT = 0.0;
112    }
113    // END local fields
114
115
116    // cloning ctor
117    public Problem(Problem orig, Cloner cloner) : base(orig, cloner) {
118    }
119
120    public Problem() : base() {
121      Parameters.Add(new FixedValueParameter<DoubleValue>("L1L2_mixing", "TODO Description", new DoubleValue(0.2)));
122      Parameters.Add(new FixedValueParameter<DoubleValue>("L1L2_weight", "TODO Description", new DoubleValue(1)));
123      Parameters.Add(new FixedValueParameter<DoubleValue>("penality", "TODO Description", new DoubleValue(0.01)));
124      Parameters.Add(new FixedValueParameter<IntValue>("fraction_depth", "TODO Description", new IntValue(4)));
125      Parameters.Add(new ValueParameter<Dataset>("dataset", "TODO Description", new Dataset()));
126      Parameters.Add(new ValueParameter<Dataset>("datasetTrainingTraining", "TODO Description", new Dataset()));
127      Parameters.Add(new ValueParameter<Dataset>("data_phi_0_08", "TODO Description", new Dataset()));
128      Parameters.Add(new ValueParameter<Dataset>("datasetTrainingTest", "TODO Description", new Dataset()));
129
130      foreach (var temperature in new[] { 350, 375, 400, 425, 450, 475, 500 }) {
131        foreach (var phip in new[] { "0_001", "0_01", "0_1", "1" }) {
132          // TODO: refactor
133          Parameters.Add(new ValueParameter<Dataset>("data_" + temperature + "_" + phip, "TODO Description", new Dataset()));
134        }
135      }
136
137      ClearLocalVariables();
138
139      DatasetParameter.ValueChanged += DatasetParameter_ValueChanged;
140      TrainingTrainingDatasetParameter.ValueChanged += DatasetParameter_ValueChanged;
141      TrainingTestDatasetParameter.ValueChanged += DatasetParameter_ValueChanged;
142      Phi0_08DatasetParameter.ValueChanged += DatasetParameter_ValueChanged;
143
144      UpdateEncoding();
145    }
146
147    [StorableHook(HookType.AfterDeserialization)]
148    private void AfterDeserialization() {
149      DatasetParameter.ValueChanged += DatasetParameter_ValueChanged;
150      TrainingTrainingDatasetParameter.ValueChanged += DatasetParameter_ValueChanged;
151      TrainingTestDatasetParameter.ValueChanged += DatasetParameter_ValueChanged;
152      Phi0_08DatasetParameter.ValueChanged += DatasetParameter_ValueChanged;
153    }
154
155    #region event handlers for clearing local variables whenever a dataset is changed
156    private void DatasetParameter_ValueChanged(object sender, EventArgs e) {
157      ClearLocalVariables();
158      UpdateEncoding();
159    }
160
161    private void UpdateEncoding() {
162      int fraction_depth = FractionDepth;
163      int number_variables = Dataset.Columns;
164      int vector_length = (2 * fraction_depth - 1) * number_variables; //number of variables n because n-1 coefficients + 1 constant
165
166      Encoding = new MultiEncoding()
167        .Add(new RealVectorEncoding("coeff", length: vector_length, min: -3, max: 3))
168        .Add(new BinaryVectorEncoding("binar", length: vector_length))
169      ;
170    }
171    #endregion
172
173    public override IDeepCloneable Clone(Cloner cloner) {
174      return new Problem(this, cloner);
175    }
176
177    public override double Evaluate(Individual individual, IRandom random) {
178      if(dataMatrix == null) {
179        InitializeTransformedDatasets();
180      }
181
182      double LS = 0.0; //sum of squared residuals -> SQR = LS
183      double quality = 0.0;
184
185      double[] coeff = individual.RealVector("coeff").ToArray();
186      bool[] binar = individual.BinaryVector("binar").ToArray();
187
188      LS = calculateLS(coeff, binar, dataMatrix);
189      quality = calculateQuality(coeff, binar, LS);
190
191      #region Nelder-Mead-Methode
192      double uniformPeturbation = 1.0;
193      double tolerance = 1e-3;
194      int maxEvals = 250;
195      int numSearches = 4;
196      int numPoints = dataMatrix.Rows / 5; // 20% of the training samples
197
198      double[] optimizedCoeffShort = new double[numVariablesUsed(binar)];
199      double[] optimizedCoeffLong = new double[coeff.Length];
200
201      for (int count = 0; count < numSearches; count++) {
202        int indexShort = 0;
203        for (int indexLong = 0; indexLong < coeff.Length; indexLong++) {
204          if (binar[indexLong] == true) {
205            optimizedCoeffShort[indexShort] = coeff[indexLong];
206            indexShort++;
207          }
208        }
209
210        SimplexConstant[] constants = new SimplexConstant[numVariablesUsed(binar)];
211        for (int i = 0; i < numVariablesUsed(binar); i++) {
212          constants[i] = new SimplexConstant(optimizedCoeffShort[i], uniformPeturbation);
213        }
214
215        //ObjFunctionNelderMead NelderMead = new ObjFunctionNelderMead(this, dataNelderMead, coeff, binar);
216        ObjFunctionNelderMead NelderMead = new ObjFunctionNelderMead(this, dataMatrix, numPoints, coeff, binar);
217
218        ObjectiveFunctionDelegate objFunction = new ObjectiveFunctionDelegate(NelderMead.objFunctionNelderMead);
219        RegressionResult result = NelderMeadSimplex.Regress(constants, tolerance, maxEvals, objFunction);
220
221        optimizedCoeffShort = result.Constants;
222        int shortIndex = 0;
223        for (int i = 0; i < coeff.Length; i++) {
224          if (binar[i] == true) {
225            optimizedCoeffLong[i] = optimizedCoeffShort[shortIndex];
226            shortIndex++;
227          } else {
228            optimizedCoeffLong[i] = coeff[i];
229          }
230        }
231
232        double newLS = calculateLS(optimizedCoeffLong, binar, dataMatrix);
233        double newQuality = calculateQuality(optimizedCoeffLong, binar, newLS);
234
235        if (newQuality < quality) {
236          // dont set new coeff!
237          individual["coeff"] = new RealVector(optimizedCoeffLong);
238          LS = newLS;
239          quality = newQuality;
240        }
241      }
242      #endregion
243
244      //doubleThrowInvOrNaN(quality, "quality");
245      individual["LS"] = new DoubleValue(LS); //Least Squares - Sum of Squares Residual, short: SQR
246      individual["Depth"] = new DoubleValue(calculateDepth(coeff, binar));
247      individual["CoeffNumber"] = new DoubleValue((double)numVariablesUsed(binar));
248      individual["MSE"] = new DoubleValue(LS / dataMatrix.Rows);
249      individual["MSETest"] = new DoubleValue(0.0); // only calculated for the best solution in a generation (see Analyze())
250      individual["R2"] = new DoubleValue(1.0 - LS / SQT); //R^2 = 1 - SQR/SQT
251
252      return quality;
253    }
254
255    private void InitializeTransformedDatasets() {
256      //dataMatrix = new DoubleMatrix(this.createMatrix("dataset", 1000));
257      dataMatrix = new DoubleMatrix(TrainingTrainingDataset.ToArray(TrainingTrainingDataset.VariableNames, Enumerable.Range(0, TrainingTrainingDataset.Rows)));
258      dataNelderMead = new DoubleMatrix(Phi0_08Dataset.ToArray(Phi0_08Dataset.VariableNames, Enumerable.Range(0, Phi0_08Dataset.Rows)));
259      dataMatrixTest = new DoubleMatrix(TrainingTestDataset.ToArray(TrainingTestDataset.VariableNames, Enumerable.Range(0, TrainingTestDataset.Rows)));
260      //minMaxTransformation(0);
261      //minMaxTransformation(1);
262      //log10Transformation(2); // no negativ values!
263      //minMaxTransformation(2); // Min and Max have to be different!
264
265      transformation = new Transformation(dataMatrix);
266      // vars["dataMatrix"] = dataMatrix;
267
268      transformation.useSameTransformation(dataNelderMead);
269      // vars["dataNelderMead"] = dataNelderMead;
270
271      transformation.useSameTransformation(dataMatrixTest);
272      // vars["dataMatrixTest"] = dataMatrixTest;
273
274      arithmeticMean = calculateArithmeticMean();
275      SQT = calculateSQT(arithmeticMean);
276    }
277
278    public override void Analyze(Individual[] individuals, double[] qualities, ResultCollection results, IRandom random) {
279      // Use vars.yourVariable to access variables in the variable store i.e. yourVariable
280      // Write or update results given the range of vectors and resulting qualities
281
282      Individual bestIndividual = null;
283      double bestQuality = double.MaxValue;
284      int theBest = 0;
285
286      for (int i = 0; i < qualities.Count(); i++) {
287        if (qualities[i] < bestQuality) {
288          bestQuality = qualities[i];
289          bestIndividual = individuals[i].Copy();
290          theBest = i;
291        }
292      }
293
294      /*
295      // Uncomment the following lines if you want to retrieve the best individual
296      var orderedIndividuals = individuals.Zip(qualities, (i, q) => new { Individual = i, Quality = q }).OrderBy(z => z.Quality);
297      var best = Maximization ? orderedIndividuals.Last().Individual : orderedIndividuals.First().Individual;
298      if (!results.ContainsKey("Best Solution")) {
299        results.Add(new Result("Best Solution", typeof(RealVector)));
300       
301      }
302      results["Best Solution"].Value = (IItem)best.RealVector("coeff").Clone();
303      //
304      */
305
306      bool[] binar = bestIndividual.BinaryVector("binar").ToArray();
307      double[] coeff = bestIndividual.RealVector("coeff").ToArray();
308
309      // calculate and set test-MSE (only) for the individual with best training-MSE
310      double bestLSTest = calculateLS(coeff, binar, dataMatrixTest);
311      double bestQualityTest = calculateQuality(coeff, binar, bestLSTest);
312      individuals[theBest]["MSETest"] = new DoubleValue(bestQualityTest);
313
314      InitDataTables(results);
315      InitScatterPlots(results);
316
317      #region set dataTables     
318      var datatable = (DataTable)results["datatable"].Value;
319      var coefficients = (DataTable)results["coefficients"].Value;
320      var pearsonR2 = (DataTable)results["personR2"].Value;
321
322      datatable.Rows["MSE"].Values.Add(calculateRegMSE(((DoubleValue)bestIndividual["LS"]).Value, binar, 0.0));
323      datatable.Rows["RegMSE"].Values.Add(calculateRegMSE(((DoubleValue)bestIndividual["LS"]).Value, binar, Penalty));
324      datatable.Rows["ElasticNet"].Values.Add(calculateElasticNet(((DoubleValue)bestIndividual["LS"]).Value, coeff, binar));
325
326      coefficients.Rows["Depth"].Values.Add(((DoubleValue)bestIndividual["Depth"]).Value);
327      coefficients.Rows["Number"].Values.Add(((DoubleValue)bestIndividual["CoeffNumber"]).Value);
328
329      pearsonR2.Rows["R2"].Values.Add(((DoubleValue)bestIndividual["R2"]).Value);
330      #endregion
331
332      #region set scatterPlots
333      var curves0_001 = (ScatterPlot)results["Curves0_001"].Value;
334      var curves0_01 = (ScatterPlot)results["Curves0_01"].Value;
335      var curves0_1 = (ScatterPlot)results["Curves0_1"].Value;
336      var curves1 = (ScatterPlot)results["Curves1"].Value;
337      #region curves 0.001
338      curves0_001.Rows["Temp350Convergence"].Points.Clear();
339      curves0_001.Rows["Temp375Convergence"].Points.Clear();
340      curves0_001.Rows["Temp400Convergence"].Points.Clear();
341      curves0_001.Rows["Temp425Convergence"].Points.Clear();
342      curves0_001.Rows["Temp450Convergence"].Points.Clear();
343      curves0_001.Rows["Temp475Convergence"].Points.Clear();
344      curves0_001.Rows["Temp500Convergence"].Points.Clear();
345      sampleCurve(curves0_001, "Temp350Convergence", coeff, binar, 350.0, 0.001, 100, 50.0);
346      sampleCurve(curves0_001, "Temp375Convergence", coeff, binar, 375.0, 0.001, 100, 50.0);
347      sampleCurve(curves0_001, "Temp400Convergence", coeff, binar, 400.0, 0.001, 100, 50.0);
348      sampleCurve(curves0_001, "Temp425Convergence", coeff, binar, 425.0, 0.001, 100, 50.0);
349      sampleCurve(curves0_001, "Temp450Convergence", coeff, binar, 450.0, 0.001, 100, 50.0);
350      sampleCurve(curves0_001, "Temp475Convergence", coeff, binar, 475.0, 0.001, 100, 50.0);
351      sampleCurve(curves0_001, "Temp500Convergence", coeff, binar, 500.0, 0.001, 100, 50.0);
352
353      results.AddOrUpdateResult("Curves Phi' = 0.001", curves0_001);
354      #endregion
355      #region curves 0.01
356      curves0_01.Rows["Temp350Convergence"].Points.Clear();
357      curves0_01.Rows["Temp375Convergence"].Points.Clear();
358      curves0_01.Rows["Temp400Convergence"].Points.Clear();
359      curves0_01.Rows["Temp425Convergence"].Points.Clear();
360      curves0_01.Rows["Temp450Convergence"].Points.Clear();
361      curves0_01.Rows["Temp475Convergence"].Points.Clear();
362      curves0_01.Rows["Temp500Convergence"].Points.Clear();
363      sampleCurve(curves0_01, "Temp350Convergence", coeff, binar, 350.0, 0.01, 100, 60.0);
364      sampleCurve(curves0_01, "Temp375Convergence", coeff, binar, 375.0, 0.01, 100, 60.0);
365      sampleCurve(curves0_01, "Temp400Convergence", coeff, binar, 400.0, 0.01, 100, 60.0);
366      sampleCurve(curves0_01, "Temp425Convergence", coeff, binar, 425.0, 0.01, 100, 60.0);
367      sampleCurve(curves0_01, "Temp450Convergence", coeff, binar, 450.0, 0.01, 100, 60.0);
368      sampleCurve(curves0_01, "Temp475Convergence", coeff, binar, 475.0, 0.01, 100, 60.0);
369      sampleCurve(curves0_01, "Temp500Convergence", coeff, binar, 500.0, 0.01, 100, 60.0);
370
371      results.AddOrUpdateResult("Curves Phi' = 0.01", curves0_01);
372      #endregion
373      #region curves 0.1
374      curves0_1.Rows["Temp350Convergence"].Points.Clear();
375      curves0_1.Rows["Temp375Convergence"].Points.Clear();
376      curves0_1.Rows["Temp400Convergence"].Points.Clear();
377      curves0_1.Rows["Temp425Convergence"].Points.Clear();
378      curves0_1.Rows["Temp450Convergence"].Points.Clear();
379      curves0_1.Rows["Temp475Convergence"].Points.Clear();
380      curves0_1.Rows["Temp500Convergence"].Points.Clear();
381      sampleCurve(curves0_1, "Temp350Convergence", coeff, binar, 350.0, 0.1, 100, 70.0);
382      sampleCurve(curves0_1, "Temp375Convergence", coeff, binar, 375.0, 0.1, 100, 70.0);
383      sampleCurve(curves0_1, "Temp400Convergence", coeff, binar, 400.0, 0.1, 100, 70.0);
384      sampleCurve(curves0_1, "Temp425Convergence", coeff, binar, 425.0, 0.1, 100, 70.0);
385      sampleCurve(curves0_1, "Temp450Convergence", coeff, binar, 450.0, 0.1, 100, 70.0);
386      sampleCurve(curves0_1, "Temp475Convergence", coeff, binar, 475.0, 0.1, 100, 70.0);
387      sampleCurve(curves0_1, "Temp500Convergence", coeff, binar, 500.0, 0.1, 100, 70.0);
388
389      results.AddOrUpdateResult("Curves Phi' = 0.1", curves0_1);
390      #endregion
391      #region curves 1
392      curves1.Rows["Temp350Convergence"].Points.Clear();
393      curves1.Rows["Temp375Convergence"].Points.Clear();
394      curves1.Rows["Temp400Convergence"].Points.Clear();
395      curves1.Rows["Temp425Convergence"].Points.Clear();
396      curves1.Rows["Temp450Convergence"].Points.Clear();
397      curves1.Rows["Temp475Convergence"].Points.Clear();
398      curves1.Rows["Temp500Convergence"].Points.Clear();
399      sampleCurve(curves1, "Temp350Convergence", coeff, binar, 350.0, 1, 100, 80.0);
400      sampleCurve(curves1, "Temp375Convergence", coeff, binar, 375.0, 1, 100, 80.0);
401      sampleCurve(curves1, "Temp400Convergence", coeff, binar, 400.0, 1, 100, 80.0);
402      sampleCurve(curves1, "Temp425Convergence", coeff, binar, 425.0, 1, 100, 80.0);
403      sampleCurve(curves1, "Temp450Convergence", coeff, binar, 450.0, 1, 100, 80.0);
404      sampleCurve(curves1, "Temp475Convergence", coeff, binar, 475.0, 1, 100, 80.0);
405      sampleCurve(curves1, "Temp500Convergence", coeff, binar, 500.0, 1, 100, 80.0);
406
407      results.AddOrUpdateResult("Curves Phi' = 1", curves1);
408      #endregion
409      #endregion
410    }
411
412    private void InitDataTables(ResultCollection results) {
413      DataTable datatable, coefficients, pearsonR2;
414      if (!results.ContainsKey("Datatable")) {
415        datatable = new DataTable("Fitness");
416        DataRow MSE = new DataRow("MSE");
417        DataRow RegMSE = new DataRow("RegMSE");
418        DataRow ElasticNet = new DataRow("ElasticNet");
419        datatable.Rows.Add(MSE);
420        datatable.Rows.Add(RegMSE);
421        datatable.Rows.Add(ElasticNet);
422        results.Add(new Result("Datatable", datatable));
423      }
424
425
426      if (!results.ContainsKey("Coefficients")) {
427        coefficients = new DataTable("Coefficients");
428        DataRow Depth = new DataRow("Depth");
429        DataRow Number = new DataRow("Number");
430        coefficients.Rows.Add(Depth);
431        coefficients.Rows.Add(Number);
432        results.Add(new Result("Coefficients", coefficients));
433      }
434
435      if (!results.ContainsKey("Pearson R2")) {
436        pearsonR2 = new DataTable("PearsonR2");
437        DataRow R2 = new DataRow("R2");
438        pearsonR2.Rows.Add(R2);
439        results.Add(new Result("Pearson R2", pearsonR2));
440      }
441    }
442
443    private void InitScatterPlots(ResultCollection results) {
444      #region curves 0.001
445      ScatterPlot curves0_001 = new ScatterPlot("Curves0_001", "Kurven mit Phi'=0.001");
446      readCurve(curves0_001, "data6082_350_0_001", "Temp350", System.Drawing.Color.Blue);
447      readCurve(curves0_001, "data6082_375_0_001", "Temp375", System.Drawing.Color.Orange);
448      readCurve(curves0_001, "data6082_400_0_001", "Temp400", System.Drawing.Color.Red);
449      readCurve(curves0_001, "data6082_425_0_001", "Temp425", System.Drawing.Color.Green);
450      readCurve(curves0_001, "data6082_450_0_001", "Temp450", System.Drawing.Color.Gray);
451      readCurve(curves0_001, "data6082_475_0_001", "Temp475", System.Drawing.Color.Olive);
452      readCurve(curves0_001, "data6082_500_0_001", "Temp500", System.Drawing.Color.Gold);
453
454      var empty0_001 = new Point2D<double>[0];
455      ScatterPlotDataRow Temp350Convergence0_001 = new ScatterPlotDataRow("Temp350Convergence", "Temp350Convergence", empty0_001);
456      Temp350Convergence0_001.VisualProperties.Color = System.Drawing.Color.Blue;
457      Temp350Convergence0_001.VisualProperties.PointSize = 2;
458      curves0_001.Rows.Add(Temp350Convergence0_001);
459
460      ScatterPlotDataRow Temp375Convergence0_001 = new ScatterPlotDataRow("Temp375Convergence", "Temp375Convergence", empty0_001);
461      Temp375Convergence0_001.VisualProperties.Color = System.Drawing.Color.Orange;
462      Temp375Convergence0_001.VisualProperties.PointSize = 2;
463      curves0_001.Rows.Add(Temp375Convergence0_001);
464
465      ScatterPlotDataRow Temp400Convergence0_001 = new ScatterPlotDataRow("Temp400Convergence", "Temp400Convergence", empty0_001);
466      Temp400Convergence0_001.VisualProperties.Color = System.Drawing.Color.Red;
467      Temp400Convergence0_001.VisualProperties.PointSize = 2;
468      curves0_001.Rows.Add(Temp400Convergence0_001);
469
470      ScatterPlotDataRow Temp425Convergence0_001 = new ScatterPlotDataRow("Temp425Convergence", "Temp425Convergence", empty0_001);
471      Temp425Convergence0_001.VisualProperties.Color = System.Drawing.Color.Green;
472      Temp425Convergence0_001.VisualProperties.PointSize = 2;
473      curves0_001.Rows.Add(Temp425Convergence0_001);
474
475      ScatterPlotDataRow Temp450Convergence0_001 = new ScatterPlotDataRow("Temp450Convergence", "Temp450Convergence", empty0_001);
476      Temp450Convergence0_001.VisualProperties.Color = System.Drawing.Color.Gray;
477      Temp450Convergence0_001.VisualProperties.PointSize = 2;
478      curves0_001.Rows.Add(Temp450Convergence0_001);
479
480      ScatterPlotDataRow Temp475Convergence0_001 = new ScatterPlotDataRow("Temp475Convergence", "Temp475Convergence", empty0_001);
481      Temp475Convergence0_001.VisualProperties.Color = System.Drawing.Color.Olive;
482      Temp475Convergence0_001.VisualProperties.PointSize = 2;
483      curves0_001.Rows.Add(Temp475Convergence0_001);
484
485      ScatterPlotDataRow Temp500Convergence0_001 = new ScatterPlotDataRow("Temp500Convergence", "Temp500Convergence", empty0_001);
486      Temp500Convergence0_001.VisualProperties.Color = System.Drawing.Color.Gold;
487      Temp500Convergence0_001.VisualProperties.PointSize = 2;
488      curves0_001.Rows.Add(Temp500Convergence0_001);
489      results.Add(new Result("Curves0_001", curves0_001));
490
491      #endregion
492      #region curves 0.01
493      ScatterPlot curves0_01 = new ScatterPlot("Curves0_01", "Kurven mit Phi'=0.01");
494      readCurve(curves0_01, "data6082_350_0_01", "Temp350", System.Drawing.Color.Blue);
495      readCurve(curves0_01, "data6082_375_0_01", "Temp375", System.Drawing.Color.Orange);
496      readCurve(curves0_01, "data6082_400_0_01", "Temp400", System.Drawing.Color.Red);
497      readCurve(curves0_01, "data6082_425_0_01", "Temp425", System.Drawing.Color.Green);
498      readCurve(curves0_01, "data6082_450_0_01", "Temp450", System.Drawing.Color.Gray);
499      readCurve(curves0_01, "data6082_475_0_01", "Temp475", System.Drawing.Color.Olive);
500      readCurve(curves0_01, "data6082_500_0_01", "Temp500", System.Drawing.Color.Gold);
501
502      var empty0_01 = new Point2D<double>[0];
503      ScatterPlotDataRow Temp350Convergence0_01 = new ScatterPlotDataRow("Temp350Convergence", "Temp350Convergence", empty0_01);
504      Temp350Convergence0_01.VisualProperties.Color = System.Drawing.Color.Blue;
505      Temp350Convergence0_01.VisualProperties.PointSize = 2;
506      curves0_01.Rows.Add(Temp350Convergence0_01);
507
508      ScatterPlotDataRow Temp375Convergence0_01 = new ScatterPlotDataRow("Temp375Convergence", "Temp375Convergence", empty0_01);
509      Temp375Convergence0_01.VisualProperties.Color = System.Drawing.Color.Orange;
510      Temp375Convergence0_01.VisualProperties.PointSize = 2;
511      curves0_01.Rows.Add(Temp375Convergence0_01);
512
513      ScatterPlotDataRow Temp400Convergence0_01 = new ScatterPlotDataRow("Temp400Convergence", "Temp400Convergence", empty0_01);
514      Temp400Convergence0_01.VisualProperties.Color = System.Drawing.Color.Red;
515      Temp400Convergence0_01.VisualProperties.PointSize = 2;
516      curves0_01.Rows.Add(Temp400Convergence0_01);
517
518      ScatterPlotDataRow Temp425Convergence0_01 = new ScatterPlotDataRow("Temp425Convergence", "Temp425Convergence", empty0_01);
519      Temp425Convergence0_01.VisualProperties.Color = System.Drawing.Color.Green;
520      Temp425Convergence0_01.VisualProperties.PointSize = 2;
521      curves0_01.Rows.Add(Temp425Convergence0_01);
522
523      ScatterPlotDataRow Temp450Convergence0_01 = new ScatterPlotDataRow("Temp450Convergence", "Temp450Convergence", empty0_01);
524      Temp450Convergence0_01.VisualProperties.Color = System.Drawing.Color.Gray;
525      Temp450Convergence0_01.VisualProperties.PointSize = 2;
526      curves0_01.Rows.Add(Temp450Convergence0_01);
527
528      ScatterPlotDataRow Temp475Convergence0_01 = new ScatterPlotDataRow("Temp475Convergence", "Temp475Convergence", empty0_01);
529      Temp475Convergence0_01.VisualProperties.Color = System.Drawing.Color.Olive;
530      Temp475Convergence0_01.VisualProperties.PointSize = 2;
531      curves0_01.Rows.Add(Temp475Convergence0_01);
532
533      ScatterPlotDataRow Temp500Convergence0_01 = new ScatterPlotDataRow("Temp500Convergence", "Temp500Convergence", empty0_01);
534      Temp500Convergence0_01.VisualProperties.Color = System.Drawing.Color.Gold;
535      Temp500Convergence0_01.VisualProperties.PointSize = 2;
536      curves0_01.Rows.Add(Temp500Convergence0_01);
537      results.Add(new Result("Curves0_01", curves0_01));
538
539      #endregion
540      #region curves 0.1
541      ScatterPlot curves0_1 = new ScatterPlot("Curves0_1", "Kurven mit Phi'=0.1");
542      readCurve(curves0_1, "data6082_350_0_1", "Temp350", System.Drawing.Color.Blue);
543      readCurve(curves0_1, "data6082_375_0_1", "Temp375", System.Drawing.Color.Orange);
544      readCurve(curves0_1, "data6082_400_0_1", "Temp400", System.Drawing.Color.Red);
545      readCurve(curves0_1, "data6082_425_0_1", "Temp425", System.Drawing.Color.Green);
546      readCurve(curves0_1, "data6082_450_0_1", "Temp450", System.Drawing.Color.Gray);
547      readCurve(curves0_1, "data6082_475_0_1", "Temp475", System.Drawing.Color.Olive);
548      readCurve(curves0_1, "data6082_500_0_1", "Temp500", System.Drawing.Color.Gold);
549
550      var empty0_1 = new Point2D<double>[0];
551      ScatterPlotDataRow Temp350Convergence0_1 = new ScatterPlotDataRow("Temp350Convergence", "Temp350Convergence", empty0_1);
552      Temp350Convergence0_1.VisualProperties.Color = System.Drawing.Color.Blue;
553      Temp350Convergence0_1.VisualProperties.PointSize = 2;
554      curves0_1.Rows.Add(Temp350Convergence0_1);
555
556      ScatterPlotDataRow Temp375Convergence0_1 = new ScatterPlotDataRow("Temp375Convergence", "Temp375Convergence", empty0_1);
557      Temp375Convergence0_1.VisualProperties.Color = System.Drawing.Color.Orange;
558      Temp375Convergence0_1.VisualProperties.PointSize = 2;
559      curves0_1.Rows.Add(Temp375Convergence0_1);
560
561      ScatterPlotDataRow Temp400Convergence0_1 = new ScatterPlotDataRow("Temp400Convergence", "Temp400Convergence", empty0_1);
562      Temp400Convergence0_1.VisualProperties.Color = System.Drawing.Color.Red;
563      Temp400Convergence0_1.VisualProperties.PointSize = 2;
564      curves0_1.Rows.Add(Temp400Convergence0_1);
565
566      ScatterPlotDataRow Temp425Convergence0_1 = new ScatterPlotDataRow("Temp425Convergence", "Temp425Convergence", empty0_1);
567      Temp425Convergence0_1.VisualProperties.Color = System.Drawing.Color.Green;
568      Temp425Convergence0_1.VisualProperties.PointSize = 2;
569      curves0_1.Rows.Add(Temp425Convergence0_1);
570
571      ScatterPlotDataRow Temp450Convergence0_1 = new ScatterPlotDataRow("Temp450Convergence", "Temp450Convergence", empty0_1);
572      Temp450Convergence0_1.VisualProperties.Color = System.Drawing.Color.Gray;
573      Temp450Convergence0_1.VisualProperties.PointSize = 2;
574      curves0_1.Rows.Add(Temp450Convergence0_1);
575
576      ScatterPlotDataRow Temp475Convergence0_1 = new ScatterPlotDataRow("Temp475Convergence", "Temp475Convergence", empty0_1);
577      Temp475Convergence0_1.VisualProperties.Color = System.Drawing.Color.Olive;
578      Temp475Convergence0_1.VisualProperties.PointSize = 2;
579      curves0_1.Rows.Add(Temp475Convergence0_1);
580
581      ScatterPlotDataRow Temp500Convergence0_1 = new ScatterPlotDataRow("Temp500Convergence", "Temp500Convergence", empty0_1);
582      Temp500Convergence0_1.VisualProperties.Color = System.Drawing.Color.Gold;
583      Temp500Convergence0_1.VisualProperties.PointSize = 2;
584      curves0_1.Rows.Add(Temp500Convergence0_1);
585      results.Add(new Result("Curves0_1", curves0_1));
586      #endregion
587      #region curves 1
588      ScatterPlot curves1 = new ScatterPlot("Curves1", "Kurven mit Phi'=1");
589      readCurve(curves1, "data6082_350_1", "Temp350", System.Drawing.Color.Blue);
590      readCurve(curves1, "data6082_375_1", "Temp375", System.Drawing.Color.Orange);
591      readCurve(curves1, "data6082_400_1", "Temp400", System.Drawing.Color.Red);
592      readCurve(curves1, "data6082_425_1", "Temp425", System.Drawing.Color.Green);
593      readCurve(curves1, "data6082_450_1", "Temp450", System.Drawing.Color.Gray);
594      readCurve(curves1, "data6082_475_1", "Temp475", System.Drawing.Color.Olive);
595      readCurve(curves1, "data6082_500_1", "Temp500", System.Drawing.Color.Gold);
596
597      var empty1 = new Point2D<double>[0];
598      ScatterPlotDataRow Temp350Convergence1 = new ScatterPlotDataRow("Temp350Convergence", "Temp350Convergence", empty1);
599      Temp350Convergence1.VisualProperties.Color = System.Drawing.Color.Blue;
600      Temp350Convergence1.VisualProperties.PointSize = 2;
601      curves1.Rows.Add(Temp350Convergence1);
602
603      ScatterPlotDataRow Temp375Convergence1 = new ScatterPlotDataRow("Temp375Convergence", "Temp375Convergence", empty1);
604      Temp375Convergence1.VisualProperties.Color = System.Drawing.Color.Orange;
605      Temp375Convergence1.VisualProperties.PointSize = 2;
606      curves1.Rows.Add(Temp375Convergence1);
607
608      ScatterPlotDataRow Temp400Convergence1 = new ScatterPlotDataRow("Temp400Convergence", "Temp400Convergence", empty1);
609      Temp400Convergence1.VisualProperties.Color = System.Drawing.Color.Red;
610      Temp400Convergence1.VisualProperties.PointSize = 2;
611      curves1.Rows.Add(Temp400Convergence1);
612
613      ScatterPlotDataRow Temp425Convergence1 = new ScatterPlotDataRow("Temp425Convergence", "Temp425Convergence", empty1);
614      Temp425Convergence1.VisualProperties.Color = System.Drawing.Color.Green;
615      Temp425Convergence1.VisualProperties.PointSize = 2;
616      curves1.Rows.Add(Temp425Convergence1);
617
618      ScatterPlotDataRow Temp450Convergence1 = new ScatterPlotDataRow("Temp450Convergence", "Temp450Convergence", empty1);
619      Temp450Convergence1.VisualProperties.Color = System.Drawing.Color.Gray;
620      Temp450Convergence1.VisualProperties.PointSize = 2;
621      curves1.Rows.Add(Temp450Convergence1);
622
623      ScatterPlotDataRow Temp475Convergence1 = new ScatterPlotDataRow("Temp475Convergence", "Temp475Convergence", empty1);
624      Temp475Convergence1.VisualProperties.Color = System.Drawing.Color.Olive;
625      Temp475Convergence1.VisualProperties.PointSize = 2;
626      curves1.Rows.Add(Temp475Convergence1);
627
628      ScatterPlotDataRow Temp500Convergence1 = new ScatterPlotDataRow("Temp500Convergence", "Temp500Convergence", empty1);
629      Temp500Convergence1.VisualProperties.Color = System.Drawing.Color.Gold;
630      Temp500Convergence1.VisualProperties.PointSize = 2;
631      curves1.Rows.Add(Temp500Convergence1);
632      results.Add(new Result("Curves1", curves1));
633      #endregion
634    }
635
636    public override IEnumerable<Individual> GetNeighbors(Individual individual, IRandom random) {
637      // Use vars.yourVariable to access variables in the variable store i.e. yourVariable
638      // Create new vectors, based on the given one that represent small changes
639      // This method is only called from move-based algorithms (Local Search, Simulated Annealing, etc.)
640      while (true) {
641        // Algorithm will draw only a finite amount of samples
642        // Change to a for-loop to return a concrete amount of neighbors
643        var neighbor = individual.Copy();
644        // For instance, perform a single bit-flip in a binary parameter
645        //var bIndex = random.Next(neighbor.BinaryVector("b").Length);
646        //neighbor.BinaryVector("b")[bIndex] = !neighbor.BinaryVector("b")[bIndex];
647        yield return neighbor;
648      }
649    }
650    #region Import & Export
651    public void Load(IRegressionProblemData data) {
652      ProblemData = data;
653    }
654
655    public IRegressionProblemData Export() {
656      return ProblemData;
657    }
658    #endregion
659
660    // Implement further classes and methods
661
662    private int numVariablesUsed(bool[] binar) {
663      int num = 0;
664      for (int i = 0; i < binar.Count(); i++) {
665        if (binar[i] == true) // eventually we need coeff < valuenearzero too
666          num++;
667      }
668      return num;
669    }
670
671
672    /* old code
673    private double[,] createMatrix(string name, int numRows) {
674      var dataset = (HeuristicLab.Problems.DataAnalysis.Dataset)vars[name];
675
676      IEnumerable<string> variables = dataset.VariableNames;
677      // var rows = dataset.Rows;
678      var rows = Enumerable.Range(0, numRows).ToArray();
679
680      double[,] dataMatrix = dataset.ToArray(variables, rows);
681
682      return dataMatrix;
683    }
684    */
685
686
687    private double evaluateContinuedFraction(double[] coeff, bool[] binar, double[] dataPoint) {
688      double value = 0.0;
689      bool firstCoeffUnequalZero = false;
690
691      for (var i = coeff.Count() - 1; i > 0; i = i - 2 * dataPoint.Count()) {
692
693        if ((linearFunctionEqualZero(coeff, binar, i - dataPoint.Count() + 1, i) == false)) {
694          firstCoeffUnequalZero = true;
695          value = this.evaluateLinearFunction(coeff, binar, i - dataPoint.Count() + 1, i, dataPoint) + value;
696        }
697
698        if (firstCoeffUnequalZero == true && i > 2 * dataPoint.Count()) { // don't take the first coeffVecPart and don't take the last coeffVecPart (both are only to add)
699          if (linearFunctionEqualZero(coeff, binar, i - 2 * dataPoint.Count() + 1, i - dataPoint.Count()) == false) {
700            if (valueNearZero(value) == true)  // no division by zero
701              return double.MaxValue;
702            else
703              value = this.evaluateLinearFunction(coeff, binar, i - 2 * dataPoint.Count() + 1, i - dataPoint.Count(), dataPoint) / value;
704          } else return double.MaxValue; // don't allow coeffVecParts in middle to be equal zero
705        }
706      }
707
708      //doubleThrowInvOrNaN(value, "evaluateContinuedFractionValue");
709      return value;
710    }
711
712    private bool linearFunctionEqualZero(double[] coeff, bool[] binar, int start, int end) {
713      for (int i = start; i <= end; i++) {
714        if (binar[i] == true && valueNearZero(coeff[i]) == false)
715          return false;
716      }
717      return true;
718    }
719
720    private double evaluateLinearFunction(double[] coeff, bool[] binar, int start, int end, double[] dataPoint) {
721      double value = 0.0;
722
723      for (int i = start; i < end; i++) {
724        if (binar[i] == true)
725          value += dataPoint[i - start] * coeff[i];
726      }
727      if (binar[end] == true)
728        value += coeff[end];
729
730      //doubleThrowInvOrNaN(value, "evaluateLinearFunctionValue");
731      return value;
732    }
733
734    private bool valueNearZero(double value) {
735      if (Math.Abs(value) < epsilon)
736        return true;
737      else
738        return false;
739    }
740
741    private double calculateDepth(double[] coeff, bool[] binar) {
742      int coeffUnequalZero = -1;
743      double depth = 0.0;
744
745      for (int i = 0; i < coeff.Count(); i++) {
746        if (valueNearZero(coeff[i]) == false && binar[i] == true)
747          coeffUnequalZero = i;
748      }
749      if (coeffUnequalZero == -1) return 0.0; // equal zero
750      if (coeffUnequalZero < dataMatrix.Columns) return 1.0; // linear function
751      depth = (double)(coeffUnequalZero - dataMatrix.Columns + 1);
752      depth = Math.Ceiling(depth / (2 * dataMatrix.Columns)) + 1.0;
753      return depth;
754    }
755
756    private double calculateRegMSE(double quality, bool[] binar, double penalty) {
757      return (quality / dataMatrix.Rows) * (1 + penalty * (double)numVariablesUsed(binar));
758    }
759
760    private double calculateElasticNet(double quality, double[] coeff, bool[] binar) {
761      double valueL1 = 0.0;
762      double valueL2 = 0.0;
763      double elasticNet = 0.0;
764      double L1L2weight = L1L2Weight;
765      double L1L2mixing = L1L2Mixing;
766
767      for (int i = 0; i < coeff.Count(); i++) {
768        if (binar[i] == true) {
769          valueL1 += Math.Abs(coeff[i]);
770          valueL2 += Math.Pow(coeff[i], 2);
771        }
772      }
773      elasticNet = quality + (L1L2weight * L1L2mixing * valueL1) + (L1L2weight * (1 - L1L2mixing) * valueL2);
774      return elasticNet / dataMatrix.Rows;
775    }
776
777    public double calculateLS(double[] coeff, bool[] binar, DoubleMatrix dataMatrix) {
778      double LS = 0.0;
779      double[] dataPoint = new double[dataMatrix.Columns];
780      for (var i = 0; i < dataMatrix.Rows; i++) {
781        for (var j = 0; j < dataMatrix.Columns; j++) {
782          dataPoint[j] = dataMatrix[i, j];
783        }
784        double continuedFractionValue = evaluateContinuedFraction(coeff, binar, dataPoint);
785        if (continuedFractionValue == double.MaxValue)
786          return double.MaxValue;
787        else
788          LS += Math.Pow((continuedFractionValue - dataMatrix[i, dataMatrix.Columns - 1]), 2);
789      }
790      return LS;
791    }
792
793    public double calculateQuality(double[] coeff, bool[] binar, double LS) {
794      switch (fit) {
795        case Fitting.MSE:
796          return calculateRegMSE(LS, binar, 0.0); //MSE
797        case Fitting.RegMSE:
798          return calculateRegMSE(LS, binar, Penalty); //RegMSE
799        case Fitting.ElasticNet:
800          return calculateElasticNet(LS, coeff, binar); // Elastic Net
801        default:
802          return calculateElasticNet(LS, coeff, binar);
803      }
804    }
805
806    private double calculateArithmeticMean() {
807      double arithmeticMean = 0.0;
808      for (int i = 0; i < dataMatrix.Rows; i++) {
809        arithmeticMean = arithmeticMean + dataMatrix[i, dataMatrix.Columns - 1];
810      }
811      return arithmeticMean / dataMatrix.Rows;
812    }
813
814    private double calculateSQT(double arithmeticMean) {
815      double SQT = 0.0;
816      for (int i = 0; i < dataMatrix.Rows; i++) {
817        SQT += Math.Pow((dataMatrix[i, dataMatrix.Columns - 1] - arithmeticMean), 2);
818      }
819      return SQT;
820    }
821
822    private static void doubleThrowInvOrNaN(double value, string valueName) {
823      if (double.IsInfinity(value) || double.IsNaN(value))
824        throw new InvalidProgramException(valueName + " is Infinity or NaN");
825    }
826
827
828    private void readCurve(ScatterPlot scatterPlot, string matrixName, string RowName, System.Drawing.Color color) {
829      var ds = ((IValueParameter<Dataset>)Parameters[matrixName]).Value;
830      DoubleMatrix dataMatrix = new DoubleMatrix(ds.ToArray(ds.VariableNames, Enumerable.Range(0, ds.Rows))); // TODO: performance / refactoring
831      var points = new Point2D<double>[dataMatrix.Rows];
832      for (int i = 0; i < dataMatrix.Rows; i++) {
833        points[i] = new Point2D<double>(dataMatrix[i, 1], dataMatrix[i, 3]);
834      }
835      ScatterPlotDataRow Curve = new ScatterPlotDataRow(RowName, RowName, points);
836      Curve.VisualProperties.Color = color;
837      scatterPlot.Rows.Add(Curve);
838    }
839
840    private void sampleCurve(ScatterPlot scatterPlot, string scatterPlotDataRow, double[] coeff, bool[] binar, double temp, double phip, int subdivisions, double maxStress) {
841      var points = new Point2D<double>[subdivisions + 1];
842      double phiMin = 0.0;
843      double phiMax = 0.70;
844      double step = (phiMax - phiMin) / subdivisions; // subdivisions > 0 !
845      double minStress = 0.0;
846
847      double[] dataPoint = new double[3] { temp, 0.0, phip };
848      //vars["sample0"] = new DoubleArray(dataPoint);
849      double[] dataPointNew = new double[4];
850      double x = 0.0;
851      //scatterPlot.Rows[scatterPlotDataRow].Points.Clear();
852      for (int i = 0; i <= subdivisions; i++) {
853        x = phiMin + (double)i * step;
854        dataPoint[0] = temp;
855        dataPoint[1] = x;
856        dataPoint[2] = phip;
857        //vars["sample1"] = new DoubleArray(dataPoint);
858        dataPoint = transformation.transform0(dataPoint);
859        //vars["sample2"] = new DoubleArray(dataPoint);
860        dataPointNew[0] = dataPoint[0];
861        dataPointNew[1] = dataPoint[1];
862        dataPointNew[2] = dataPoint[2];
863        dataPointNew[3] = 0.0;
864        points[i] = new Point2D<double>(x, evaluateContinuedFraction(coeff, binar, dataPointNew));
865        if (points[i].Y >= minStress && points[i].Y <= maxStress) {
866          scatterPlot.Rows[scatterPlotDataRow].Points.Add(points[i]);
867        }
868      }
869      //scatterPlot.Rows[ScatterPlotDataRow].Point2D.Add(points);
870    }
871  }
872}
873
Note: See TracBrowser for help on using the repository browser.