1 | using System;
|
---|
2 | using System.Collections.Generic;
|
---|
3 | using System.Linq;
|
---|
4 | using System.Threading;
|
---|
5 | using HEAL.Attic;
|
---|
6 | using HeuristicLab.Common;
|
---|
7 | using HeuristicLab.Core;
|
---|
8 | using HeuristicLab.Data;
|
---|
9 | using HeuristicLab.Problems.DataAnalysis;
|
---|
10 | using HeuristicLab.Random;
|
---|
11 |
|
---|
12 | namespace 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 | }
|
---|