source: branches/3106_AnalyticContinuedFractionsRegression/HeuristicLab.Algorithms.DataAnalysis/3.4/ContinuedFractionRegression/Algorithm.cs @ 17972

Last change on this file since 17972 was 17972, checked in by gkronber, 6 weeks ago

#3106: minor changes to alg but no big improvement

File size: 16.1 KB
Line 
1using System;
2using System.Collections.Generic;
3using System.Linq;
4using System.Threading;
5using HEAL.Attic;
6using HeuristicLab.Common;
7using HeuristicLab.Core;
8using HeuristicLab.Data;
9using HeuristicLab.Problems.DataAnalysis;
10using HeuristicLab.Random;
11
12namespace HeuristicLab.Algorithms.DataAnalysis.ContinuedFractionRegression {
13  [Item("Continuous Fraction Regression (CFR)", "TODO")]
14  [Creatable(CreatableAttribute.Categories.DataAnalysisRegression, Priority = 999)]
15  [StorableType("7A375270-EAAF-4AD1-82FF-132318D20E09")]
16  public class Algorithm : FixedDataAnalysisAlgorithm<IRegressionProblem> {
17    public override IDeepCloneable Clone(Cloner cloner) {
18      throw new NotImplementedException();
19    }
20
21    protected override void Run(CancellationToken cancellationToken) {
22      var problemData = Problem.ProblemData;
23
24      var x = problemData.Dataset.ToArray(problemData.AllowedInputVariables.Concat(new[] { problemData.TargetVariable }),
25        problemData.TrainingIndices);
26      var nVars = x.GetLength(1) - 1;
27      var rand = new MersenneTwister(31415);
28      CFRAlgorithm(nVars, depth: 6, 0.10, x, out var best, out var bestObj, rand, numGen: 200, stagnatingGens: 5, cancellationToken);
29    }
30
31    private void CFRAlgorithm(int nVars, int depth, double mutationRate, double[,] trainingData,
32      out ContinuedFraction best, out double bestObj,
33      IRandom rand, int numGen, int stagnatingGens,
34      CancellationToken cancellationToken) {
35      /* Algorithm 1 */
36      /* Generate initial population by a randomized algorithm */
37      var pop = InitialPopulation(nVars, depth, rand, trainingData);
38      best = pop.pocket;
39      bestObj = pop.pocketObjValue;
40      var bestObjGen = 0;
41      for (int gen = 1; gen <= numGen && !cancellationToken.IsCancellationRequested; gen++) {
42        /* mutate each current solution in the population */
43        var pop_mu = Mutate(pop, mutationRate, rand);
44        /* generate new population by recombination mechanism */
45        var pop_r = RecombinePopulation(pop_mu, rand, nVars);
46
47        /* local search optimization of current solutions */
48        foreach (var agent in pop_r.IterateLevels()) {
49          LocalSearchSimplex(agent.current, ref agent.currentObjValue, trainingData, rand);
50        }
51
52        foreach (var agent in pop_r.IteratePostOrder()) agent.MaintainInvariant(); // Deviates from Alg1 in paper
53
54        /* replace old population with evolved population */
55        pop = pop_r;
56
57        /* keep track of the best solution */
58        if (bestObj > pop.pocketObjValue) {
59          best = pop.pocket;
60          bestObj = pop.pocketObjValue;
61          bestObjGen = gen;
62          Results.AddOrUpdateResult("MSE (best)", new DoubleValue(bestObj));
63        }
64
65
66        if (gen > bestObjGen + stagnatingGens) {
67          bestObjGen = gen; // wait at least stagnatingGens until resetting again
68          // Reset(pop, nVars, depth, rand, trainingData);
69          InitialPopulation(nVars, depth, rand, trainingData);
70        }
71      }
72    }
73
74    private Agent InitialPopulation(int nVars, int depth, IRandom rand, double[,] trainingData) {
75      /* instantiate 13 agents in the population */
76      var pop = new Agent();
77      // see Figure 2
78      for (int i = 0; i < 3; i++) {
79        pop.children.Add(new Agent());
80        for (int j = 0; j < 3; j++) {
81          pop.children[i].children.Add(new Agent());
82        }
83      }
84
85      foreach (var agent in pop.IteratePostOrder()) {
86        agent.current = new ContinuedFraction(nVars, depth, rand);
87        agent.pocket = new ContinuedFraction(nVars, depth, rand);
88
89        agent.currentObjValue = Evaluate(agent.current, trainingData);
90        agent.pocketObjValue = Evaluate(agent.pocket, trainingData);
91
92        /* within each agent, the pocket solution always holds the better value of guiding
93         * function than its current solution
94         */
95        agent.MaintainInvariant();
96      }
97      return pop;
98    }
99
100    // TODO: reset is not described in the paper
101    private void Reset(Agent root, int nVars, int depth, IRandom rand, double[,] trainingData) {
102      root.pocket = new ContinuedFraction(nVars, depth, rand);
103      root.current = new ContinuedFraction(nVars, depth, rand);
104
105      root.currentObjValue = Evaluate(root.current, trainingData);
106      root.pocketObjValue = Evaluate(root.pocket, trainingData);
107
108      /* within each agent, the pocket solution always holds the better value of guiding
109       * function than its current solution
110       */
111      root.MaintainInvariant();
112    }
113
114
115
116    private Agent RecombinePopulation(Agent pop, IRandom rand, int nVars) {
117      var l = pop;
118      if (pop.children.Count > 0) {
119        var s1 = pop.children[0];
120        var s2 = pop.children[1];
121        var s3 = pop.children[2];
122        l.current = Recombine(l.pocket, s1.current, SelectRandomOp(rand), rand, nVars);
123        s3.current = Recombine(s3.pocket, l.current, SelectRandomOp(rand), rand, nVars);
124        s1.current = Recombine(s1.pocket, s2.current, SelectRandomOp(rand), rand, nVars);
125        s2.current = Recombine(s2.pocket, s3.current, SelectRandomOp(rand), rand, nVars);
126      }
127
128      foreach (var child in pop.children) {
129        RecombinePopulation(child, rand, nVars);
130      }
131      return pop;
132    }
133
134    private Func<bool[], bool[], bool[]> SelectRandomOp(IRandom rand) {
135      bool[] union(bool[] a, bool[] b) {
136        var res = new bool[a.Length];
137        for (int i = 0; i < a.Length; i++) res[i] = a[i] || b[i];
138        return res;
139      }
140      bool[] intersect(bool[] a, bool[] b) {
141        var res = new bool[a.Length];
142        for (int i = 0; i < a.Length; i++) res[i] = a[i] && b[i];
143        return res;
144      }
145      bool[] symmetricDifference(bool[] a, bool[] b) {
146        var res = new bool[a.Length];
147        for (int i = 0; i < a.Length; i++) res[i] = a[i] ^ b[i];
148        return res;
149      }
150      switch (rand.Next(3)) {
151        case 0: return union;
152        case 1: return intersect;
153        case 2: return symmetricDifference;
154        default: throw new ArgumentException();
155      }
156    }
157
158    private static ContinuedFraction Recombine(ContinuedFraction p1, ContinuedFraction p2, Func<bool[], bool[], bool[]> op, IRandom rand, int nVars) {
159      ContinuedFraction ch = new ContinuedFraction() { h = new Term[p1.h.Length] };
160      /* apply a recombination operator chosen uniformly at random on variable sof two parents into offspring */
161      ch.vars = op(p1.vars, p2.vars);
162
163      /* recombine the coefficients for each term (h) of the continued fraction */
164      for (int i = 0; i < p1.h.Length; i++) {
165        var coefa = p1.h[i].coef; var varsa = p1.h[i].vars;
166        var coefb = p2.h[i].coef; var varsb = p2.h[i].vars;
167
168        /* recombine coefficient values for variables */
169        var coefx = new double[nVars];
170        var varsx = new bool[nVars]; // TODO: deviates from paper -> check
171        for (int vi = 1; vi < nVars; vi++) {
172          if (ch.vars[vi]) {
173            if (varsa[vi] && varsb[vi]) {
174              coefx[vi] = coefa[vi] + (rand.NextDouble() * 5 - 1) * (coefb[vi] - coefa[vi]) / 3.0;
175              varsx[vi] = true;
176            } else if (varsa[vi]) {
177              coefx[vi] = coefa[vi];
178              varsx[vi] = true;
179            } else if (varsb[vi]) {
180              coefx[vi] = coefb[vi];
181              varsx[vi] = true;
182            }
183          }
184        }
185        /* update new coefficients of the term in offspring */
186        ch.h[i] = new Term() { coef = coefx, vars = varsx };
187        /* compute new value of constant (beta) for term hi in the offspring solution ch using
188         * beta of p1.hi and p2.hi */
189        ch.h[i].beta = p1.h[i].beta + (rand.NextDouble() * 5 - 1) * (p2.h[i].beta - p1.h[i].beta) / 3.0;
190      }
191      /* update current solution and apply local search */
192      // return LocalSearchSimplex(ch, trainingData); // Deviates from paper because Alg1 also has LocalSearch after Recombination
193      return ch;
194    }
195
196    private Agent Mutate(Agent pop, double mutationRate, IRandom rand) {
197      foreach (var agent in pop.IterateLevels()) {
198        if (rand.NextDouble() < mutationRate) {
199          if (agent.currentObjValue < 1.2 * agent.pocketObjValue ||
200             agent.currentObjValue > 2 * agent.pocketObjValue)
201            ToggleVariables(agent.current, rand); // major mutation
202          else
203            ModifyVariable(agent.current, rand); // soft mutation
204        }
205      }
206      return pop;
207    }
208
209    private void ToggleVariables(ContinuedFraction cfrac, IRandom rand) {
210      double coinToss(double a, double b) {
211        return rand.NextDouble() < 0.5 ? a : b;
212      }
213
214      /* select a variable index uniformly at random */
215      int N = cfrac.vars.Length;
216      var vIdx = rand.Next(N);
217
218      /* for each depth of continued fraction, toggle the selection of variables of the term (h) */
219      foreach (var h in cfrac.h) {
220        /* Case 1: cfrac variable is turned ON: Turn OFF the variable, and either 'Remove' or
221         * 'Remember' the coefficient value at random */
222        if (cfrac.vars[vIdx]) {
223          h.vars[vIdx] = false;
224          h.coef[vIdx] = coinToss(0, h.coef[vIdx]);
225        } else {
226          /* Case 2: term variable is turned OFF: Turn ON the variable, and either 'Remove'
227           * or 'Replace' the coefficient with a random value between -3 and 3 at random */
228          if (!h.vars[vIdx]) {
229            h.vars[vIdx] = true;
230            h.coef[vIdx] = coinToss(0, rand.NextDouble() * 6 - 3);
231          }
232        }
233      }
234      /* toggle the randomly selected variable */
235      cfrac.vars[vIdx] = !cfrac.vars[vIdx];
236    }
237
238    private void ModifyVariable(ContinuedFraction cfrac, IRandom rand) {
239      /* randomly select a variable which is turned ON */
240      var candVars = cfrac.vars.Count(vi => vi);
241      if (candVars == 0) return; // no variable active
242      var vIdx = rand.Next(candVars);
243
244      /* randomly select a term (h) of continued fraction */
245      var h = cfrac.h[rand.Next(cfrac.h.Length)];
246
247      /* modify the coefficient value*/
248      if (h.vars[vIdx]) {
249        h.coef[vIdx] = 0.0;
250      } else {
251        h.coef[vIdx] = rand.NextDouble() * 6 - 3;
252      }
253      /* Toggle the randomly selected variable */
254      h.vars[vIdx] = !h.vars[vIdx];
255    }
256
257    private static double Evaluate(ContinuedFraction cfrac, double[,] trainingData) {
258      var dataPoint = new double[trainingData.GetLength(1) - 1];
259      var yIdx = trainingData.GetLength(1) - 1;
260      double sum = 0.0;
261      for (int r = 0; r < trainingData.GetLength(0); r++) {
262        for (int c = 0; c < dataPoint.Length; c++) {
263          dataPoint[c] = trainingData[r, c];
264        }
265        var y = trainingData[r, yIdx];
266        var pred = Evaluate(cfrac, dataPoint);
267        var res = y - pred;
268        sum += res * res;
269      }
270      var delta = 0.1; // TODO
271      return sum / trainingData.GetLength(0) * (1 + delta * cfrac.vars.Count(vi => vi));
272    }
273
274    private static double Evaluate(ContinuedFraction cfrac, double[] dataPoint) {
275      var res = 0.0;
276      for (int i = cfrac.h.Length - 1; i > 1; i -= 2) {
277        var hi = cfrac.h[i];
278        var hi1 = cfrac.h[i - 1];
279        var denom = hi.beta + dot(hi.vars, hi.coef, dataPoint) + res;
280        var numerator = hi1.beta + dot(hi1.vars, hi1.coef, dataPoint);
281        res = numerator / denom;
282      }
283      return res;
284    }
285
286    private static double dot(bool[] filter, double[] x, double[] y) {
287      var s = 0.0;
288      for (int i = 0; i < x.Length; i++)
289        if (filter[i])
290          s += x[i] * y[i];
291      return s;
292    }
293
294
295    private static void LocalSearchSimplex(ContinuedFraction ch, ref double quality, double[,] trainingData, IRandom rand) {
296      double uniformPeturbation = 1.0;
297      double tolerance = 1e-3;
298      int maxEvals = 250;
299      int numSearches = 4;
300      var numRows = trainingData.GetLength(0);
301      int numSelectedRows = numRows / 5; // 20% of the training samples
302
303      quality = Evaluate(ch, trainingData); // get quality with origial coefficients
304
305      double[] origCoeff = ExtractCoeff(ch);
306      if (origCoeff.Length == 0) return; // no parameters to optimize
307
308      var bestQuality = quality;
309      var bestCoeff = origCoeff;
310
311      var fittingData = SelectRandomRows(trainingData, numSelectedRows, rand);
312
313      double objFunc(double[] curCoeff) {
314        SetCoeff(ch, curCoeff);
315        return Evaluate(ch, fittingData);
316      }
317
318      for (int count = 0; count < numSearches; count++) {
319
320        SimplexConstant[] constants = new SimplexConstant[origCoeff.Length];
321        for (int i = 0; i < origCoeff.Length; i++) {
322          constants[i] = new SimplexConstant(origCoeff[i], uniformPeturbation);
323        }
324
325        RegressionResult result = NelderMeadSimplex.Regress(constants, tolerance, maxEvals, objFunc);
326
327        var optimizedCoeff = result.Constants;
328        SetCoeff(ch, optimizedCoeff);
329
330        var newQuality = Evaluate(ch, trainingData);
331
332        // TODO: optionally use regularization (ridge / LASSO)
333
334        if (newQuality < bestQuality) {
335          bestCoeff = optimizedCoeff;
336          bestQuality = newQuality;
337        }
338      } // reps
339
340      SetCoeff(ch, bestCoeff);
341      quality = bestQuality;
342    }
343
344    private static double[,] SelectRandomRows(double[,] trainingData, int numSelectedRows, IRandom rand) {
345      var numRows = trainingData.GetLength(0);
346      var numCols = trainingData.GetLength(1);
347      var selectedRows = Enumerable.Range(0, numRows).Shuffle(rand).Take(numSelectedRows).ToArray();
348      var subset = new double[numSelectedRows, numCols];
349      var i = 0;
350      foreach (var r in selectedRows) {
351        for (int c = 0; c < numCols; c++) {
352          subset[i, c] = trainingData[r, c];
353        }
354        i++;
355      }
356      return subset;
357    }
358
359    private static double[] ExtractCoeff(ContinuedFraction ch) {
360      var coeff = new List<double>();
361      foreach (var hi in ch.h) {
362        coeff.Add(hi.beta);
363        for (int vIdx = 0; vIdx < hi.vars.Length; vIdx++) {
364          if (hi.vars[vIdx]) coeff.Add(hi.coef[vIdx]);
365        }
366      }
367      return coeff.ToArray();
368    }
369
370    private static void SetCoeff(ContinuedFraction ch, double[] curCoeff) {
371      int k = 0;
372      foreach (var hi in ch.h) {
373        hi.beta = curCoeff[k++];
374        for (int vIdx = 0; vIdx < hi.vars.Length; vIdx++) {
375          if (hi.vars[vIdx]) hi.coef[vIdx] = curCoeff[k++];
376        }
377      }
378    }
379  }
380
381  public class Agent {
382    public ContinuedFraction pocket;
383    public double pocketObjValue;
384    public ContinuedFraction current;
385    public double currentObjValue;
386
387    public IList<Agent> children = new List<Agent>();
388
389    public IEnumerable<Agent> IterateLevels() {
390      var agents = new List<Agent>() { this };
391      IterateLevelsRec(this, agents);
392      return agents;
393    }
394    public IEnumerable<Agent> IteratePostOrder() {
395      var agents = new List<Agent>();
396      IteratePostOrderRec(this, agents);
397      return agents;
398    }
399
400    internal void MaintainInvariant() {
401      foreach (var child in children) {
402        MaintainInvariant(parent: this, child);
403      }
404      if (currentObjValue < pocketObjValue) {
405        Swap(ref pocket, ref current);
406        Swap(ref pocketObjValue, ref currentObjValue);
407      }
408    }
409
410
411    private static void MaintainInvariant(Agent parent, Agent child) {
412      if (child.pocketObjValue < parent.pocketObjValue) {
413        Swap(ref child.pocket, ref parent.pocket);
414        Swap(ref child.pocketObjValue, ref parent.pocketObjValue);
415      }
416    }
417
418    private void IterateLevelsRec(Agent agent, List<Agent> agents) {
419      foreach (var child in agent.children) {
420        agents.Add(child);
421      }
422      foreach (var child in agent.children) {
423        IterateLevelsRec(child, agents);
424      }
425    }
426
427    private void IteratePostOrderRec(Agent agent, List<Agent> agents) {
428      foreach (var child in agent.children) {
429        IteratePostOrderRec(child, agents);
430      }
431      agents.Add(agent);
432    }
433
434
435    private static void Swap<T>(ref T a, ref T b) {
436      var temp = a;
437      a = b;
438      b = temp;
439    }
440  }
441}
Note: See TracBrowser for help on using the repository browser.