source: branches/MemPRAlgorithm/HeuristicLab.Algorithms.MemPR/3.3/MemPRAlgorithm.cs @ 14496

Last change on this file since 14496 was 14496, checked in by abeham, 3 years ago

#2701:

  • Reusing similiarty calculator in BinaryMemPR
  • Fixing distance calculation for linear linkage and LinearLinkageMemPR
  • Small changes to base algorithm
  • Added biased model trainer for permutation (rank and fitness)
  • Fixing best known quality calculation for GCP
File size: 34.0 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2016 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.ComponentModel;
25using System.Linq;
26using System.Runtime.CompilerServices;
27using System.Threading;
28using HeuristicLab.Algorithms.MemPR.Interfaces;
29using HeuristicLab.Algorithms.MemPR.Util;
30using HeuristicLab.Analysis;
31using HeuristicLab.Common;
32using HeuristicLab.Core;
33using HeuristicLab.Data;
34using HeuristicLab.Optimization;
35using HeuristicLab.Parameters;
36using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
37
38namespace HeuristicLab.Algorithms.MemPR {
39  [Item("MemPR Algorithm", "Base class for MemPR algorithms")]
40  [StorableClass]
41  public abstract class MemPRAlgorithm<TProblem, TSolution, TPopulationContext, TSolutionContext> : BasicAlgorithm, INotifyPropertyChanged
42      where TProblem : class, IItem, ISingleObjectiveHeuristicOptimizationProblem, ISingleObjectiveProblemDefinition
43      where TSolution : class, IItem
44      where TPopulationContext : MemPRPopulationContext<TProblem, TSolution, TPopulationContext, TSolutionContext>, new()
45      where TSolutionContext : MemPRSolutionContext<TProblem, TSolution, TPopulationContext, TSolutionContext> {
46    private const double MutationProbabilityMagicConst = 0.1;
47
48    public override Type ProblemType {
49      get { return typeof(TProblem); }
50    }
51
52    public new TProblem Problem {
53      get { return (TProblem)base.Problem; }
54      set { base.Problem = value; }
55    }
56
57    protected string QualityName {
58      get { return Problem != null && Problem.Evaluator != null ? Problem.Evaluator.QualityParameter.ActualName : null; }
59    }
60
61    public int? MaximumEvaluations {
62      get {
63        var val = ((OptionalValueParameter<IntValue>)Parameters["MaximumEvaluations"]).Value;
64        return val != null ? val.Value : (int?)null;
65      }
66      set {
67        var param = (OptionalValueParameter<IntValue>)Parameters["MaximumEvaluations"];
68        param.Value = value.HasValue ? new IntValue(value.Value) : null;
69      }
70    }
71
72    public TimeSpan? MaximumExecutionTime {
73      get {
74        var val = ((OptionalValueParameter<TimeSpanValue>)Parameters["MaximumExecutionTime"]).Value;
75        return val != null ? val.Value : (TimeSpan?)null;
76      }
77      set {
78        var param = (OptionalValueParameter<TimeSpanValue>)Parameters["MaximumExecutionTime"];
79        param.Value = value.HasValue ? new TimeSpanValue(value.Value) : null;
80      }
81    }
82
83    public double? TargetQuality {
84      get {
85        var val = ((OptionalValueParameter<DoubleValue>)Parameters["TargetQuality"]).Value;
86        return val != null ? val.Value : (double?)null;
87      }
88      set {
89        var param = (OptionalValueParameter<DoubleValue>)Parameters["TargetQuality"];
90        param.Value = value.HasValue ? new DoubleValue(value.Value) : null;
91      }
92    }
93
94    protected FixedValueParameter<IntValue> MaximumPopulationSizeParameter {
95      get { return ((FixedValueParameter<IntValue>)Parameters["MaximumPopulationSize"]); }
96    }
97    public int MaximumPopulationSize {
98      get { return MaximumPopulationSizeParameter.Value.Value; }
99      set { MaximumPopulationSizeParameter.Value.Value = value; }
100    }
101
102    public bool SetSeedRandomly {
103      get { return ((FixedValueParameter<BoolValue>)Parameters["SetSeedRandomly"]).Value.Value; }
104      set { ((FixedValueParameter<BoolValue>)Parameters["SetSeedRandomly"]).Value.Value = value; }
105    }
106
107    public int Seed {
108      get { return ((FixedValueParameter<IntValue>)Parameters["Seed"]).Value.Value; }
109      set { ((FixedValueParameter<IntValue>)Parameters["Seed"]).Value.Value = value; }
110    }
111
112    public IAnalyzer Analyzer {
113      get { return ((ValueParameter<IAnalyzer>)Parameters["Analyzer"]).Value; }
114      set { ((ValueParameter<IAnalyzer>)Parameters["Analyzer"]).Value = value; }
115    }
116
117    public IConstrainedValueParameter<ISolutionModelTrainer<TPopulationContext>> SolutionModelTrainerParameter {
118      get { return (IConstrainedValueParameter<ISolutionModelTrainer<TPopulationContext>>)Parameters["SolutionModelTrainer"]; }
119    }
120
121    public IConstrainedValueParameter<ILocalSearch<TSolutionContext>> LocalSearchParameter {
122      get { return (IConstrainedValueParameter<ILocalSearch<TSolutionContext>>)Parameters["LocalSearch"]; }
123    }
124
125    [Storable]
126    private TPopulationContext context;
127    public TPopulationContext Context {
128      get { return context; }
129      protected set {
130        if (context == value) return;
131        context = value;
132        OnPropertyChanged("State");
133      }
134    }
135
136    [Storable]
137    private BestAverageWorstQualityAnalyzer qualityAnalyzer;
138
139    [StorableConstructor]
140    protected MemPRAlgorithm(bool deserializing) : base(deserializing) { }
141    protected MemPRAlgorithm(MemPRAlgorithm<TProblem, TSolution, TPopulationContext, TSolutionContext> original, Cloner cloner) : base(original, cloner) {
142      context = cloner.Clone(original.context);
143      qualityAnalyzer = cloner.Clone(original.qualityAnalyzer);
144      RegisterEventHandlers();
145    }
146    protected MemPRAlgorithm() {
147      Parameters.Add(new ValueParameter<IAnalyzer>("Analyzer", "The analyzer to apply to the population.", new MultiAnalyzer()));
148      Parameters.Add(new FixedValueParameter<IntValue>("MaximumPopulationSize", "The maximum size of the population that is evolved.", new IntValue(20)));
149      Parameters.Add(new OptionalValueParameter<IntValue>("MaximumEvaluations", "The maximum number of solution evaluations."));
150      Parameters.Add(new OptionalValueParameter<TimeSpanValue>("MaximumExecutionTime", "The maximum runtime.", new TimeSpanValue(TimeSpan.FromMinutes(1))));
151      Parameters.Add(new OptionalValueParameter<DoubleValue>("TargetQuality", "The target quality at which the algorithm terminates."));
152      Parameters.Add(new FixedValueParameter<BoolValue>("SetSeedRandomly", "Whether each run of the algorithm should be conducted with a new random seed.", new BoolValue(true)));
153      Parameters.Add(new FixedValueParameter<IntValue>("Seed", "The random number seed that is used in case SetSeedRandomly is false.", new IntValue(0)));
154      Parameters.Add(new ConstrainedValueParameter<ISolutionModelTrainer<TPopulationContext>>("SolutionModelTrainer", "The object that creates a solution model that can be sampled."));
155      Parameters.Add(new ConstrainedValueParameter<ILocalSearch<TSolutionContext>>("LocalSearch", "The local search operator to use."));
156
157      qualityAnalyzer = new BestAverageWorstQualityAnalyzer();
158      RegisterEventHandlers();
159    }
160
161    [StorableHook(HookType.AfterDeserialization)]
162    private void AfterDeserialization() {
163      RegisterEventHandlers();
164    }
165
166    private void RegisterEventHandlers() {
167      MaximumPopulationSizeParameter.Value.ValueChanged += MaximumPopulationSizeOnChanged;
168    }
169
170    private void MaximumPopulationSizeOnChanged(object sender, EventArgs eventArgs) {
171      if (ExecutionState == ExecutionState.Started || ExecutionState == ExecutionState.Paused)
172        throw new InvalidOperationException("Cannot change maximum population size before algorithm finishes.");
173      Prepare();
174    }
175
176    protected override void OnProblemChanged() {
177      base.OnProblemChanged();
178      qualityAnalyzer.MaximizationParameter.ActualName = Problem.MaximizationParameter.Name;
179      qualityAnalyzer.MaximizationParameter.Hidden = true;
180      qualityAnalyzer.QualityParameter.ActualName = Problem.Evaluator.QualityParameter.ActualName;
181      qualityAnalyzer.QualityParameter.Depth = 1;
182      qualityAnalyzer.QualityParameter.Hidden = true;
183      qualityAnalyzer.BestKnownQualityParameter.ActualName = Problem.BestKnownQualityParameter.Name;
184      qualityAnalyzer.BestKnownQualityParameter.Hidden = true;
185
186      var multiAnalyzer = Analyzer as MultiAnalyzer;
187      if (multiAnalyzer != null) {
188        multiAnalyzer.Operators.Clear();
189        if (Problem != null) {
190          foreach (var analyzer in Problem.Operators.OfType<IAnalyzer>()) {
191            foreach (var param in analyzer.Parameters.OfType<IScopeTreeLookupParameter>())
192              param.Depth = 1;
193            multiAnalyzer.Operators.Add(analyzer, analyzer.EnabledByDefault);
194          }
195        }
196        multiAnalyzer.Operators.Add(qualityAnalyzer, qualityAnalyzer.EnabledByDefault);
197      }
198    }
199
200    public override void Prepare() {
201      base.Prepare();
202      Results.Clear();
203      Context = null;
204    }
205
206    protected virtual TPopulationContext CreateContext() {
207      return new TPopulationContext();
208    }
209
210    protected sealed override void Run(CancellationToken token) {
211      if (Context == null) {
212        Context = CreateContext();
213        if (SetSeedRandomly) Seed = new System.Random().Next();
214        Context.Random.Reset(Seed);
215        Context.Scope.Variables.Add(new Variable("Results", Results));
216        Context.Problem = Problem;
217      }
218
219      if (MaximumExecutionTime.HasValue)
220        CancellationTokenSource.CancelAfter(MaximumExecutionTime.Value);
221
222      IExecutionContext context = null;
223      foreach (var item in Problem.ExecutionContextItems)
224        context = new Core.ExecutionContext(context, item, Context.Scope);
225      context = new Core.ExecutionContext(context, this, Context.Scope);
226      Context.Parent = context;
227
228      if (!Context.Initialized) {
229        // We initialize the population with two local optima
230        while (Context.PopulationCount < 2) {
231          var child = Create(token);
232          Context.LocalSearchEvaluations += HillClimb(child, token);
233          if (Replace(child, token) >= 0)
234            Analyze(token);
235          token.ThrowIfCancellationRequested();
236          if (Terminate()) return;
237        }
238        Context.LocalSearchEvaluations /= 2;
239        Context.Initialized = true;
240      }
241
242      while (!Terminate()) {
243        Iterate(token);
244        Analyze(token);
245        token.ThrowIfCancellationRequested();
246      }
247    }
248
249    private void Iterate(CancellationToken token) {
250      var replaced = false;
251
252      var i1 = Context.Random.Next(Context.PopulationCount);
253      var i2 = Context.Random.Next(Context.PopulationCount);
254      while (i1 == i2) i2 = Context.Random.Next(Context.PopulationCount);
255
256      var p1 = Context.AtPopulation(i1);
257      var p2 = Context.AtPopulation(i2);
258
259      var parentDist = Dist(p1, p2);
260
261      ISingleObjectiveSolutionScope<TSolution> offspring = null;
262      int replPos = -1;
263
264      if (Context.Random.NextDouble() > parentDist * parentDist) {
265        offspring = BreedAndImprove(p1, p2, token);
266        replPos = Replace(offspring, token);
267        if (replPos >= 0) {
268          replaced = true;
269          Context.ByBreeding++;
270        }
271      }
272
273      if (Context.Random.NextDouble() < Math.Sqrt(parentDist)) {
274        offspring = RelinkAndImprove(p1, p2, token);
275        replPos = Replace(offspring, token);
276        if (replPos >= 0) {
277          replaced = true;
278          Context.ByRelinking++;
279        }
280      }
281
282      offspring = PerformSampling(token);
283      replPos = Replace(offspring, token);
284      if (replPos >= 0) {
285        replaced = true;
286        Context.BySampling++;
287      }
288
289      if (!replaced) {
290        offspring = Create(token);
291        if (HillclimbingSuited(offspring)) {
292          HillClimb(offspring, token);
293          replPos = Replace(offspring, token);
294          if (replPos >= 0) {
295            Context.ByHillclimbing++;
296            replaced = true;
297          }
298        } else {
299          offspring = (ISingleObjectiveSolutionScope<TSolution>)Context.AtPopulation(Context.Random.Next(Context.PopulationCount)).Clone();
300          Mutate(offspring, token);
301          PerformTabuWalk(offspring, Context.LocalSearchEvaluations, token);
302          replPos = Replace(offspring, token);
303          if (replPos >= 0) {
304            Context.ByTabuwalking++;
305            replaced = true;
306          }
307        }
308      }
309      Context.Iterations++;
310    }
311
312    protected void Analyze(CancellationToken token) {
313      IResult res;
314      if (!Results.TryGetValue("EvaluatedSolutions", out res))
315        Results.Add(new Result("EvaluatedSolutions", new IntValue(Context.EvaluatedSolutions)));
316      else ((IntValue)res.Value).Value = Context.EvaluatedSolutions;
317      if (!Results.TryGetValue("Iterations", out res))
318        Results.Add(new Result("Iterations", new IntValue(Context.Iterations)));
319      else ((IntValue)res.Value).Value = Context.Iterations;
320      if (!Results.TryGetValue("LocalSearch Evaluations", out res))
321        Results.Add(new Result("LocalSearch Evaluations", new IntValue(Context.LocalSearchEvaluations)));
322      else ((IntValue)res.Value).Value = Context.LocalSearchEvaluations;
323      if (!Results.TryGetValue("ByBreeding", out res))
324        Results.Add(new Result("ByBreeding", new IntValue(Context.ByBreeding)));
325      else ((IntValue)res.Value).Value = Context.ByBreeding;
326      if (!Results.TryGetValue("ByRelinking", out res))
327        Results.Add(new Result("ByRelinking", new IntValue(Context.ByRelinking)));
328      else ((IntValue)res.Value).Value = Context.ByRelinking;
329      if (!Results.TryGetValue("BySampling", out res))
330        Results.Add(new Result("BySampling", new IntValue(Context.BySampling)));
331      else ((IntValue)res.Value).Value = Context.BySampling;
332      if (!Results.TryGetValue("ByHillclimbing", out res))
333        Results.Add(new Result("ByHillclimbing", new IntValue(Context.ByHillclimbing)));
334      else ((IntValue)res.Value).Value = Context.ByHillclimbing;
335      if (!Results.TryGetValue("ByTabuwalking", out res))
336        Results.Add(new Result("ByTabuwalking", new IntValue(Context.ByTabuwalking)));
337      else ((IntValue)res.Value).Value = Context.ByTabuwalking;
338
339      var sp = new ScatterPlot("Parent1 vs Offspring", "");
340      sp.Rows.Add(new ScatterPlotDataRow("corr", "", Context.BreedingStat.Select(x => new Point2D<double>(x.Item1, x.Item3))) { VisualProperties = { PointSize = 6 }});
341      if (!Results.TryGetValue("BreedingStat1", out res)) {
342        Results.Add(new Result("BreedingStat1", sp));
343      } else res.Value = sp;
344
345      sp = new ScatterPlot("Parent2 vs Offspring", "");
346      sp.Rows.Add(new ScatterPlotDataRow("corr", "", Context.BreedingStat.Select(x => new Point2D<double>(x.Item2, x.Item3))) { VisualProperties = { PointSize = 6 } });
347      if (!Results.TryGetValue("BreedingStat2", out res)) {
348        Results.Add(new Result("BreedingStat2", sp));
349      } else res.Value = sp;
350
351      sp = new ScatterPlot("Solution vs Local Optimum", "");
352      sp.Rows.Add(new ScatterPlotDataRow("corr", "", Context.HillclimbingStat.Select(x => new Point2D<double>(x.Item1, x.Item2))) { VisualProperties = { PointSize = 6 } });
353      if (!Results.TryGetValue("HillclimbingStat", out res)) {
354        Results.Add(new Result("HillclimbingStat", sp));
355      } else res.Value = sp;
356
357      sp = new ScatterPlot("Solution vs Tabu Walk", "");
358      sp.Rows.Add(new ScatterPlotDataRow("corr", "", Context.TabuwalkingStat.Select(x => new Point2D<double>(x.Item1, x.Item2))) { VisualProperties = { PointSize = 6 } });
359      if (!Results.TryGetValue("TabuwalkingStat", out res)) {
360        Results.Add(new Result("TabuwalkingStat", sp));
361      } else res.Value = sp;
362
363      RunOperator(Analyzer, Context.Scope, token);
364    }
365
366    protected int Replace(ISingleObjectiveSolutionScope<TSolution> child, CancellationToken token) {
367      if (double.IsNaN(child.Fitness)) {
368        Evaluate(child, token);
369        Context.IncrementEvaluatedSolutions(1);
370      }
371      if (IsBetter(child.Fitness, Context.BestQuality)) {
372        Context.BestQuality = child.Fitness;
373        Context.BestSolution = (TSolution)child.Solution.Clone();
374      }
375
376      var popSize = MaximumPopulationSize;
377      if (Context.Population.All(p => !Eq(p, child))) {
378
379        if (Context.PopulationCount < popSize) {
380          Context.AddToPopulation(child);
381          return Context.PopulationCount - 1;
382        }
383       
384        // The set of replacement candidates consists of all solutions at least as good as the new one
385        var candidates = Context.Population.Select((p, i) => new { Index = i, Individual = p })
386                                         .Where(x => x.Individual.Fitness == child.Fitness
387                                           || IsBetter(child, x.Individual)).ToList();
388        if (candidates.Count == 0) return -1;
389
390        var repCand = -1;
391        var avgChildDist = 0.0;
392        var minChildDist = double.MaxValue;
393        var plateau = new List<int>();
394        var worstPlateau = -1;
395        var minAvgPlateauDist = double.MaxValue;
396        var minPlateauDist = double.MaxValue;
397        // If there are equally good solutions it is first tried to replace one of those
398        // The criteria for replacement is that the new solution has better average distance
399        // to all other solutions at this "plateau"
400        foreach (var c in candidates.Where(x => x.Individual.Fitness == child.Fitness)) {
401          var dist = Dist(c.Individual, child);
402          avgChildDist += dist;
403          if (dist < minChildDist) minChildDist = dist;
404          plateau.Add(c.Index);
405        }
406        if (plateau.Count > 2) {
407          avgChildDist /= plateau.Count;
408          foreach (var p in plateau) {
409            var avgDist = 0.0;
410            var minDist = double.MaxValue;
411            foreach (var q in plateau) {
412              if (p == q) continue;
413              var dist = Dist(Context.AtPopulation(p), Context.AtPopulation(q));
414              avgDist += dist;
415              if (dist < minDist) minDist = dist;
416            }
417
418            var d = Dist(Context.AtPopulation(p), child);
419            avgDist += d;
420            avgDist /= plateau.Count;
421            if (d < minDist) minDist = d;
422
423            if (minDist < minPlateauDist || (minDist == minPlateauDist && avgDist < avgChildDist)) {
424              minAvgPlateauDist = avgDist;
425              minPlateauDist = minDist;
426              worstPlateau = p;
427            }
428          }
429          if (minPlateauDist < minChildDist || (minPlateauDist == minChildDist && minAvgPlateauDist < avgChildDist))
430            repCand = worstPlateau;
431        }
432
433        if (repCand < 0) {
434          // If no solution at the same plateau were identified for replacement
435          // a worse solution with smallest distance is chosen
436          var minDist = double.MaxValue;
437          foreach (var c in candidates.Where(x => IsBetter(child, x.Individual))) {
438            var d = Dist(c.Individual, child);
439            if (d < minDist) {
440              minDist = d;
441              repCand = c.Index;
442            }
443          }
444        }
445
446        // If no replacement was identified, this can only mean that there are
447        // no worse solutions and those on the same plateau are all better
448        // stretched out than the new one
449        if (repCand < 0) return -1;
450       
451        Context.ReplaceAtPopulation(repCand, child);
452        return repCand;
453      }
454      return -1;
455    }
456
457    [MethodImpl(MethodImplOptions.AggressiveInlining)]
458    protected bool IsBetter(ISingleObjectiveSolutionScope<TSolution> a, ISingleObjectiveSolutionScope<TSolution> b) {
459      return IsBetter(a.Fitness, b.Fitness);
460    }
461    [MethodImpl(MethodImplOptions.AggressiveInlining)]
462    protected bool IsBetter(double a, double b) {
463      return double.IsNaN(b) && !double.IsNaN(a)
464        || Problem.Maximization && a > b
465        || !Problem.Maximization && a < b;
466    }
467   
468    protected abstract bool Eq(ISingleObjectiveSolutionScope<TSolution> a, ISingleObjectiveSolutionScope<TSolution> b);
469    protected abstract double Dist(ISingleObjectiveSolutionScope<TSolution> a, ISingleObjectiveSolutionScope<TSolution> b);
470    protected abstract ISingleObjectiveSolutionScope<TSolution> ToScope(TSolution code, double fitness = double.NaN);
471    protected abstract ISolutionSubspace<TSolution> CalculateSubspace(IEnumerable<TSolution> solutions, bool inverse = false);
472    protected virtual void Evaluate(ISingleObjectiveSolutionScope<TSolution> scope, CancellationToken token) {
473      var prob = Problem as ISingleObjectiveProblemDefinition;
474      if (prob != null) {
475        var ind = new SingleEncodingIndividual(prob.Encoding, scope);
476        scope.Fitness = prob.Evaluate(ind, Context.Random);
477      } else RunOperator(Problem.Evaluator, scope, token);
478    }
479
480    #region Create
481    protected virtual ISingleObjectiveSolutionScope<TSolution> Create(CancellationToken token) {
482      var child = ToScope(null);
483      RunOperator(Problem.SolutionCreator, child, token);
484      return child;
485    }
486    #endregion
487
488    #region Improve
489    protected virtual int HillClimb(ISingleObjectiveSolutionScope<TSolution> scope, CancellationToken token, ISolutionSubspace<TSolution> subspace = null) {
490      if (double.IsNaN(scope.Fitness)) {
491        Evaluate(scope, token);
492        Context.IncrementEvaluatedSolutions(1);
493      }
494      var before = scope.Fitness;
495      var lscontext = Context.CreateSingleSolutionContext(scope);
496      LocalSearchParameter.Value.Optimize(lscontext);
497      var after = scope.Fitness;
498      Context.HillclimbingStat.Add(Tuple.Create(before, after));
499      Context.IncrementEvaluatedSolutions(lscontext.EvaluatedSolutions);
500      return lscontext.EvaluatedSolutions;
501    }
502
503    protected virtual void PerformTabuWalk(ISingleObjectiveSolutionScope<TSolution> scope, int steps, CancellationToken token, ISolutionSubspace<TSolution> subspace = null) {
504      if (double.IsNaN(scope.Fitness)) {
505        Evaluate(scope, token);
506        Context.IncrementEvaluatedSolutions(1);
507      }
508      var before = scope.Fitness;
509      var newScope = (ISingleObjectiveSolutionScope<TSolution>)scope.Clone();
510      var newSteps = TabuWalk(newScope, steps, token, subspace);
511      Context.TabuwalkingStat.Add(Tuple.Create(before, newScope.Fitness));
512      //Context.HcSteps = (int)Math.Ceiling(Context.HcSteps * (1.0 + Context.TabuwalkingStat.Count) / (2.0 + Context.TabuwalkingStat.Count) + newSteps / (2.0 + Context.TabuwalkingStat.Count));
513      if (IsBetter(newScope, scope) || (newScope.Fitness == scope.Fitness && Dist(newScope, scope) > 0))
514        scope.Adopt(newScope);
515    }
516    protected abstract int TabuWalk(ISingleObjectiveSolutionScope<TSolution> scope, int maxEvals, CancellationToken token, ISolutionSubspace<TSolution> subspace = null);
517    protected virtual void TabuClimb(ISingleObjectiveSolutionScope<TSolution> scope, int steps, CancellationToken token, ISolutionSubspace<TSolution> subspace = null) {
518      if (double.IsNaN(scope.Fitness)) {
519        Evaluate(scope, token);
520        Context.IncrementEvaluatedSolutions(1);
521      }
522      var before = scope.Fitness;
523      var newScope = (ISingleObjectiveSolutionScope<TSolution>)scope.Clone();
524      var newSteps = TabuWalk(newScope, steps, token, subspace);
525      Context.TabuwalkingStat.Add(Tuple.Create(before, newScope.Fitness));
526      //Context.HcSteps = (int)Math.Ceiling(Context.HcSteps * (1.0 + Context.TabuwalkingStat.Count) / (2.0 + Context.TabuwalkingStat.Count) + newSteps / (2.0 + Context.TabuwalkingStat.Count));
527      if (IsBetter(newScope, scope) || (newScope.Fitness == scope.Fitness && Dist(newScope, scope) > 0))
528        scope.Adopt(newScope);
529    }
530    #endregion
531   
532    #region Breed
533    protected virtual ISingleObjectiveSolutionScope<TSolution> PerformBreeding(CancellationToken token) {
534      if (Context.PopulationCount < 2) throw new InvalidOperationException("Cannot breed from population with less than 2 individuals.");
535      var i1 = Context.Random.Next(Context.PopulationCount);
536      var i2 = Context.Random.Next(Context.PopulationCount);
537      while (i1 == i2) i2 = Context.Random.Next(Context.PopulationCount);
538
539      var p1 = Context.AtPopulation(i1);
540      var p2 = Context.AtPopulation(i2);
541
542      if (double.IsNaN(p1.Fitness)) {
543        Evaluate(p1, token);
544        Context.IncrementEvaluatedSolutions(1);
545      }
546      if (double.IsNaN(p2.Fitness)) {
547        Evaluate(p2, token);
548        Context.IncrementEvaluatedSolutions(1);
549      }
550
551      return BreedAndImprove(p1, p2, token);
552    }
553
554    protected virtual ISingleObjectiveSolutionScope<TSolution> BreedAndImprove(ISingleObjectiveSolutionScope<TSolution> p1, ISingleObjectiveSolutionScope<TSolution> p2, CancellationToken token) {
555      var offspring = Cross(p1, p2, token);
556      var subspace = CalculateSubspace(new[] { p1.Solution, p2.Solution });
557      if (Context.Random.NextDouble() < MutationProbabilityMagicConst) {
558        Mutate(offspring, token, subspace); // mutate the solutions, especially to widen the sub-space
559      }
560      if (double.IsNaN(offspring.Fitness)) {
561        Evaluate(offspring, token);
562        Context.IncrementEvaluatedSolutions(1);
563      }
564      Context.BreedingStat.Add(Tuple.Create(p1.Fitness, p2.Fitness, offspring.Fitness));
565      if ((IsBetter(offspring, p1) && IsBetter(offspring, p2))
566        || Context.Population.Any(p => IsBetter(offspring, p))) return offspring;
567
568      if (HillclimbingSuited(offspring))
569        HillClimb(offspring, token, subspace); // perform hillclimb in the solution sub-space
570      return offspring;
571    }
572
573    protected abstract ISingleObjectiveSolutionScope<TSolution> Cross(ISingleObjectiveSolutionScope<TSolution> p1, ISingleObjectiveSolutionScope<TSolution> p2, CancellationToken token);
574    protected abstract void Mutate(ISingleObjectiveSolutionScope<TSolution> offspring, CancellationToken token, ISolutionSubspace<TSolution> subspace = null);
575    #endregion
576
577    #region Relink
578    protected virtual ISingleObjectiveSolutionScope<TSolution> PerformRelinking(CancellationToken token) {
579      if (Context.PopulationCount < 2) throw new InvalidOperationException("Cannot breed from population with less than 2 individuals.");
580      var i1 = Context.Random.Next(Context.PopulationCount);
581      var i2 = Context.Random.Next(Context.PopulationCount);
582      while (i1 == i2) i2 = Context.Random.Next(Context.PopulationCount);
583
584      var p1 = Context.AtPopulation(i1);
585      var p2 = Context.AtPopulation(i2);
586
587      return RelinkAndImprove(p1, p2, token);
588    }
589
590    protected virtual ISingleObjectiveSolutionScope<TSolution> RelinkAndImprove(ISingleObjectiveSolutionScope<TSolution> a, ISingleObjectiveSolutionScope<TSolution> b, CancellationToken token) {
591      var child = Relink(a, b, token);
592      if (IsBetter(child, a) && IsBetter(child, b)) return child;
593
594      var dist1 = Dist(child, a);
595      var dist2 = Dist(child, b);
596      if (dist1 > 0 && dist2 > 0) {
597        var subspace = CalculateSubspace(new[] { a.Solution, b.Solution }, inverse: true);
598        if (HillclimbingSuited(child)) {
599          HillClimb(child, token, subspace); // perform hillclimb in solution sub-space
600        }
601      }
602      return child;
603    }
604
605    protected abstract ISingleObjectiveSolutionScope<TSolution> Relink(ISingleObjectiveSolutionScope<TSolution> a, ISingleObjectiveSolutionScope<TSolution> b, CancellationToken token);
606    #endregion
607
608    #region Sample
609    protected virtual ISingleObjectiveSolutionScope<TSolution> PerformSampling(CancellationToken token) {
610      SolutionModelTrainerParameter.Value.TrainModel(Context);
611      var sample = ToScope(Context.Model.Sample());
612      Evaluate(sample, token);
613      Context.IncrementEvaluatedSolutions(1);
614      if (Context.Population.Any(p => IsBetter(sample, p) || sample.Fitness == p.Fitness)) return sample;
615
616      if (HillclimbingSuited(sample)) {
617        var subspace = CalculateSubspace(Context.Population.Select(x => x.Solution));
618        HillClimb(sample, token, subspace);
619      }
620      return sample;
621    }
622    #endregion
623
624    protected bool HillclimbingSuited(ISingleObjectiveSolutionScope<TSolution> scope) {
625      return Context.Random.NextDouble() < ProbabilityAccept(scope, Context.HillclimbingStat);
626    }
627    protected bool HillclimbingSuited(double startingFitness) {
628      return Context.Random.NextDouble() < ProbabilityAccept(startingFitness, Context.HillclimbingStat);
629    }
630    protected bool TabuwalkingSuited(ISingleObjectiveSolutionScope<TSolution> scope) {
631      return Context.Random.NextDouble() < ProbabilityAccept(scope, Context.TabuwalkingStat);
632    }
633    protected bool TabuwalkingSuited(double startingFitness) {
634      return Context.Random.NextDouble() < ProbabilityAccept(startingFitness, Context.TabuwalkingStat);
635    }
636
637    protected double ProbabilityAccept(ISingleObjectiveSolutionScope<TSolution> scope, IList<Tuple<double, double>> data) {
638      if (double.IsNaN(scope.Fitness)) {
639        Evaluate(scope, CancellationToken.None);
640        Context.IncrementEvaluatedSolutions(1);
641      }
642      return ProbabilityAccept(scope.Fitness, data);
643    }
644    protected double ProbabilityAccept(double startingFitness, IList<Tuple<double, double>> data) {
645      if (data.Count < 10) return 1.0;
646      int[] clusterValues;
647      var centroids = CkMeans1D.Cluster(data.Select(x => x.Item1).ToArray(), 2, out clusterValues);
648      var cluster = Math.Abs(startingFitness - centroids.First().Key) < Math.Abs(startingFitness - centroids.Last().Key) ? centroids.First().Value : centroids.Last().Value;
649
650      var samples = 0;
651      double meanStart = 0, meanStartOld = 0, meanEnd = 0, meanEndOld = 0;
652      double varStart = 0, varStartOld = 0, varEnd = 0, varEndOld = 0;
653      for (var i = 0; i < data.Count; i++) {
654        if (clusterValues[i] != cluster) continue;
655
656        samples++;
657        var x = data[i].Item1;
658        var y = data[i].Item2;
659
660        if (samples == 1) {
661          meanStartOld = x;
662          meanEndOld = y;
663        } else {
664          meanStart = meanStartOld + (x - meanStartOld) / samples;
665          meanEnd = meanEndOld + (x - meanEndOld) / samples;
666          varStart = varStartOld + (x - meanStartOld) * (x - meanStart) / (samples - 1);
667          varEnd = varEndOld + (x - meanEndOld) * (x - meanEnd) / (samples - 1);
668
669          meanStartOld = meanStart;
670          meanEndOld = meanEnd;
671          varStartOld = varStart;
672          varEndOld = varEnd;
673        }
674      }
675      if (samples < 5) return 1.0;
676      var cov = data.Select((v, i) => new { Index = i, Value = v }).Where(x => clusterValues[x.Index] == cluster).Select(x => x.Value).Sum(x => (x.Item1 - meanStart) * (x.Item2 - meanEnd)) / data.Count;
677
678      var biasedMean = meanEnd + cov / varStart * (startingFitness - meanStart);
679      var biasedStdev = Math.Sqrt(varEnd - (cov * cov) / varStart);
680
681      if (Problem.Maximization) {
682        var goal = Context.Population.Min(x => x.Fitness);
683        var z = (goal - biasedMean) / biasedStdev;
684        return 1.0 - Phi(z); // P(X >= z)
685      } else {
686        var goal = Context.Population.Max(x => x.Fitness);
687        var z = (goal - biasedMean) / biasedStdev;
688        return Phi(z); // P(X <= z)
689      }
690    }
691
692    protected virtual bool Terminate() {
693      return MaximumEvaluations.HasValue && Context.EvaluatedSolutions >= MaximumEvaluations.Value
694        || MaximumExecutionTime.HasValue && ExecutionTime >= MaximumExecutionTime.Value
695        || TargetQuality.HasValue && (Problem.Maximization && Context.BestQuality >= TargetQuality.Value
696                                  || !Problem.Maximization && Context.BestQuality <= TargetQuality.Value);
697    }
698
699    public event PropertyChangedEventHandler PropertyChanged;
700    protected void OnPropertyChanged(string property) {
701      var handler = PropertyChanged;
702      if (handler != null) handler(this, new PropertyChangedEventArgs(property));
703    }
704
705    #region Engine Helper
706    protected void RunOperator(IOperator op, IScope scope, CancellationToken cancellationToken) {
707      var stack = new Stack<IOperation>();
708      stack.Push(Context.CreateChildOperation(op, scope));
709
710      while (stack.Count > 0) {
711        cancellationToken.ThrowIfCancellationRequested();
712
713        var next = stack.Pop();
714        if (next is OperationCollection) {
715          var coll = (OperationCollection)next;
716          for (int i = coll.Count - 1; i >= 0; i--)
717            if (coll[i] != null) stack.Push(coll[i]);
718        } else if (next is IAtomicOperation) {
719          var operation = (IAtomicOperation)next;
720          try {
721            next = operation.Operator.Execute((IExecutionContext)operation, cancellationToken);
722          } catch (Exception ex) {
723            stack.Push(operation);
724            if (ex is OperationCanceledException) throw ex;
725            else throw new OperatorExecutionException(operation.Operator, ex);
726          }
727          if (next != null) stack.Push(next);
728        }
729      }
730    }
731    #endregion
732
733    #region Math Helper
734    // normal distribution CDF (left of x) for N(0;1) standard normal distribution
735    // from http://www.johndcook.com/blog/csharp_phi/
736    // license: "This code is in the public domain. Do whatever you want with it, no strings attached."
737    // added: 2016-11-19 21:46 CET
738    protected static double Phi(double x) {
739      // constants
740      double a1 = 0.254829592;
741      double a2 = -0.284496736;
742      double a3 = 1.421413741;
743      double a4 = -1.453152027;
744      double a5 = 1.061405429;
745      double p = 0.3275911;
746
747      // Save the sign of x
748      int sign = 1;
749      if (x < 0)
750        sign = -1;
751      x = Math.Abs(x) / Math.Sqrt(2.0);
752
753      // A&S formula 7.1.26
754      double t = 1.0 / (1.0 + p * x);
755      double y = 1.0 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.Exp(-x * x);
756
757      return 0.5 * (1.0 + sign * y);
758    }
759    #endregion
760  }
761}
Note: See TracBrowser for help on using the repository browser.