Free cookie consent management tool by TermsFeed Policy Generator

source: trunk/HeuristicLab.Algorithms.ParameterlessPopulationPyramid/3.3/LinkageTree.cs @ 17277

Last change on this file since 17277 was 17180, checked in by swagner, 5 years ago

#2875: Removed years in copyrights

File size: 10.6 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 * and the BEACON Center for the Study of Evolution in Action.
5 *
6 * This file is part of HeuristicLab.
7 *
8 * HeuristicLab is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * HeuristicLab is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
20 */
21#endregion
22
23using System;
24using System.Collections.Generic;
25using System.Linq;
26using HeuristicLab.Common;
27using HeuristicLab.Core;
28using HeuristicLab.Encodings.BinaryVectorEncoding;
29using HEAL.Attic;
30using HeuristicLab.Random;
31
32namespace HeuristicLab.Algorithms.ParameterlessPopulationPyramid {
33  // This code is based off the publication
34  // B. W. Goldman and W. F. Punch, "Parameter-less Population Pyramid," GECCO, pp. 785–792, 2014
35  // and the original source code in C++11 available from: https://github.com/brianwgoldman/Parameter-less_Population_Pyramid
36  [StorableType("40E75AC2-ABD0-43A9-AC91-7478FDD6A399")]
37  public class LinkageTree : DeepCloneable {
38    [Storable]
39    private readonly int[][][] occurances;
40    [Storable]
41    private readonly List<int>[] clusters;
42    [Storable]
43    private List<int> clusterOrdering;
44    [Storable]
45    private readonly int length;
46    [Storable]
47    private readonly IRandom rand;
48    [Storable]
49    private bool rebuildRequired = false;
50
51
52    [StorableConstructor]
53    protected LinkageTree(StorableConstructorFlag _) { }
54
55
56    protected LinkageTree(LinkageTree original, Cloner cloner) : base(original, cloner) {
57      occurances = new int[original.occurances.Length][][];
58      //mkommend: first entry is not used, cf. ctor line 83
59      for (int i = 1; i < original.occurances.Length; i++) {
60        occurances[i] = new int[original.occurances[i].Length][];
61        for (int j = 0; j < original.occurances[i].Length; j++)
62          occurances[i][j] = original.occurances[i][j].ToArray();
63      }
64
65      clusters = original.clusters.Select(c => c.ToList()).ToArray();
66      clusterOrdering = new List<int>(original.clusterOrdering);
67      length = original.length;
68      rand = cloner.Clone(original.rand);
69      rebuildRequired = original.rebuildRequired;
70    }
71
72
73    public override IDeepCloneable Clone(Cloner cloner) {
74      return new LinkageTree(this, cloner);
75    }
76
77    public LinkageTree(int length, IRandom rand) {
78      this.length = length;
79      this.rand = rand;
80      occurances = new int[length][][];
81
82      // Create a lower triangular matrix without the diagonal
83      for (int i = 1; i < length; i++) {
84        occurances[i] = new int[i][];
85        for (int j = 0; j < i; j++) {
86          occurances[i][j] = new int[4];
87        }
88      }
89      clusters = new List<int>[2 * length - 1];
90      for (int i = 0; i < clusters.Length; i++) {
91        clusters[i] = new List<int>();
92      }
93      clusterOrdering = new List<int>();
94
95      // first "length" clusters just contain a single gene
96      for (int i = 0; i < length; i++) {
97        clusters[i].Add(i);
98      }
99    }
100
101    public void Add(BinaryVector solution) {
102      if (solution.Length != length) throw new ArgumentException("The individual has not the correct length.");
103      for (int i = 1; i < solution.Length; i++) {
104        for (int j = 0; j < i; j++) {
105          // Updates the entry of the 4 long array based on the two bits
106
107          var pattern = (Convert.ToByte(solution[j]) << 1) + Convert.ToByte(solution[i]);
108          occurances[i][j][pattern]++;
109        }
110      }
111      rebuildRequired = true;
112    }
113
114    // While "total" always has an integer value, it is a double to reduce
115    // how often type casts are needed to prevent integer divison
116    // In the GECCO paper, calculates Equation 2
117    private static double NegativeEntropy(int[] counts, double total) {
118      double sum = 0;
119      for (int i = 0; i < counts.Length; i++) {
120        if (counts[i] != 0) {
121          sum += ((counts[i] / total) * Math.Log(counts[i] / total));
122        }
123      }
124      return sum;
125    }
126
127    // Uses the frequency table to calcuate the entropy distance between two indices.
128    // In the GECCO paper, calculates Equation 1
129    private double EntropyDistance(int i, int j) {
130      int[] bits = new int[4];
131      // This ensures you are using the lower triangular part of "occurances"
132      if (i < j) {
133        int temp = i;
134        i = j;
135        j = temp;
136      }
137      var entry = occurances[i][j];
138      // extracts the occurrences of the individual bits
139      bits[0] = entry[0] + entry[2];  // i zero
140      bits[1] = entry[1] + entry[3];  // i one
141      bits[2] = entry[0] + entry[1];  // j zero
142      bits[3] = entry[2] + entry[3];  // j one
143      double total = bits[0] + bits[1];
144      // entropy of the two bits on their own
145      double separate = NegativeEntropy(bits, total);
146      // entropy of the two bits as a single unit
147      double together = NegativeEntropy(entry, total);
148      // If together there is 0 entropy, the distance is zero
149      if (together.IsAlmost(0)) {
150        return 0.0;
151      }
152      return 2 - (separate / together);
153    }
154
155
156
157    // Performs O(N^2) clustering based on the method described in:
158    // "Optimal implementations of UPGMA and other common clustering algorithms"
159    // by I. Gronau and S. Moran
160    // In the GECCO paper, Figure 2 is a simplified version of this algorithm.
161    private void Rebuild() {
162      double[][] distances = null;
163      if (distances == null) {
164        distances = new double[clusters.Length * 2 - 1][];
165        for (int i = 0; i < distances.Length; i++)
166          distances[i] = new double[clusters.Length * 2 - 1];
167      }
168
169
170      // Keep track of which clusters have not been merged
171      var topLevel = new List<int>(length);
172      for (int i = 0; i < length; i++)
173        topLevel.Add(i);
174
175      bool[] useful = new bool[clusters.Length];
176      for (int i = 0; i < useful.Length; i++)
177        useful[i] = true;
178
179      // Store the distances between all clusters
180      for (int i = 1; i < length; i++) {
181        for (int j = 0; j < i; j++) {
182          distances[i][j] = EntropyDistance(clusters[i][0], clusters[j][0]);
183          // make it symmetric
184          distances[j][i] = distances[i][j];
185        }
186      }
187      // Each iteration we add some amount to the path, and remove the last
188      // two elements.  This keeps track of how much of usable is in the path.
189      int end_of_path = 0;
190
191      // build all clusters of size greater than 1
192      for (int index = length; index < clusters.Length; index++) {
193        // Shuffle everything not yet in the path
194        topLevel.ShuffleInPlace(rand, end_of_path, topLevel.Count - 1);
195
196        // if nothing in the path, just add a random usable node
197        if (end_of_path == 0) {
198          end_of_path = 1;
199        }
200        while (end_of_path < topLevel.Count) {
201          // last node in the path
202          int final = topLevel[end_of_path - 1];
203
204          // best_index stores the location of the best thing in the top level
205          int best_index = end_of_path;
206          double min_dist = distances[final][topLevel[best_index]];
207          // check all options which might be closer to "final" than "topLevel[best_index]"
208          for (int option = end_of_path + 1; option < topLevel.Count; option++) {
209            if (distances[final][topLevel[option]] < min_dist) {
210              min_dist = distances[final][topLevel[option]];
211              best_index = option;
212            }
213          }
214          // If the current last two in the path are minimally distant
215          if (end_of_path > 1 && min_dist >= distances[final][topLevel[end_of_path - 2]]) {
216            break;
217          }
218
219          // move the best to the end of the path
220          topLevel.Swap(end_of_path, best_index);
221          end_of_path++;
222        }
223        // Last two elements in the path are the clusters to join
224        int first = topLevel[end_of_path - 2];
225        int second = topLevel[end_of_path - 1];
226
227        // Only keep a cluster if the distance between the joining clusters is > zero
228        bool keep = !distances[first][second].IsAlmost(0.0);
229        useful[first] = keep;
230        useful[second] = keep;
231
232        // create the new cluster
233        clusters[index] = clusters[first].Concat(clusters[second]).ToList();
234        // Calculate distances from all clusters to the newly created cluster
235        int i = 0;
236        int end = topLevel.Count - 1;
237        while (i <= end) {
238          int x = topLevel[i];
239          // Moves 'first' and 'second' to after "end" in topLevel
240          if (x == first || x == second) {
241            topLevel.Swap(i, end);
242            end--;
243            continue;
244          }
245          // Use the previous distances to calculate the joined distance
246          double first_distance = distances[first][x];
247          first_distance *= clusters[first].Count;
248          double second_distance = distances[second][x];
249          second_distance *= clusters[second].Count;
250          distances[x][index] = ((first_distance + second_distance)
251              / (clusters[first].Count + clusters[second].Count));
252          // make it symmetric
253          distances[index][x] = distances[x][index];
254          i++;
255        }
256
257        // Remove first and second from the path
258        end_of_path -= 2;
259        topLevel.RemoveAt(topLevel.Count - 1);
260        topLevel[topLevel.Count - 1] = index;
261      }
262      // Extract the useful clusters
263      clusterOrdering.Clear();
264      // Add all useful clusters. The last one is never useful.
265      for (int i = 0; i < useful.Length - 1; i++) {
266        if (useful[i]) clusterOrdering.Add(i);
267      }
268
269      // Shuffle before sort to ensure ties are broken randomly
270      clusterOrdering.ShuffleInPlace(rand);
271      clusterOrdering = clusterOrdering.OrderBy(i => clusters[i].Count).ToList();
272    }
273
274    public IEnumerable<List<int>> Clusters {
275      get {
276        // Just in time rebuilding
277        if (rebuildRequired) Rebuild();
278        foreach (var index in clusterOrdering) {
279          // Send out the clusters in the desired order
280          yield return clusters[index];
281        }
282      }
283    }
284  }
285}
Note: See TracBrowser for help on using the repository browser.