Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2708_ScopedAlgorithms/HeuristicLab.Algorithms.ParameterlessPopulationPyramid/3.3/LinkageTree.cs @ 16124

Last change on this file since 16124 was 12012, checked in by ascheibe, 10 years ago

#2212 merged r12008, r12009, r12010 back into trunk

File size: 9.5 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2015 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 HeuristicLab.Random;
30
31namespace HeuristicLab.Algorithms.ParameterlessPopulationPyramid {
32  // This code is based off the publication
33  // B. W. Goldman and W. F. Punch, "Parameter-less Population Pyramid," GECCO, pp. 785–792, 2014
34  // and the original source code in C++11 available from: https://github.com/brianwgoldman/Parameter-less_Population_Pyramid
35  public class LinkageTree {
36    private readonly int[][][] occurances;
37    private readonly List<int>[] clusters;
38    private List<int> clusterOrdering;
39    private readonly int length;
40    private readonly IRandom rand;
41    private bool rebuildRequired = false;
42
43    public LinkageTree(int length, IRandom rand) {
44      this.length = length;
45      this.rand = rand;
46      occurances = new int[length][][];
47
48      // Create a lower triangular matrix without the diagonal
49      for (int i = 1; i < length; i++) {
50        occurances[i] = new int[i][];
51        for (int j = 0; j < i; j++) {
52          occurances[i][j] = new int[4];
53        }
54      }
55      clusters = new List<int>[2 * length - 1];
56      for (int i = 0; i < clusters.Length; i++) {
57        clusters[i] = new List<int>();
58      }
59      clusterOrdering = new List<int>();
60
61      // first "length" clusters just contain a single gene
62      for (int i = 0; i < length; i++) {
63        clusters[i].Add(i);
64      }
65    }
66
67    public void Add(BinaryVector solution) {
68      if (solution.Length != length) throw new ArgumentException("The individual has not the correct length.");
69      for (int i = 1; i < solution.Length; i++) {
70        for (int j = 0; j < i; j++) {
71          // Updates the entry of the 4 long array based on the two bits
72
73          var pattern = (Convert.ToByte(solution[j]) << 1) + Convert.ToByte(solution[i]);
74          occurances[i][j][pattern]++;
75        }
76      }
77      rebuildRequired = true;
78    }
79
80    // While "total" always has an integer value, it is a double to reduce
81    // how often type casts are needed to prevent integer divison
82    // In the GECCO paper, calculates Equation 2
83    private static double NegativeEntropy(int[] counts, double total) {
84      double sum = 0;
85      for (int i = 0; i < counts.Length; i++) {
86        if (counts[i] != 0) {
87          sum += ((counts[i] / total) * Math.Log(counts[i] / total));
88        }
89      }
90      return sum;
91    }
92
93    // Uses the frequency table to calcuate the entropy distance between two indices.
94    // In the GECCO paper, calculates Equation 1
95    private int[] bits = new int[4];
96    private double EntropyDistance(int i, int j) {
97      // This ensures you are using the lower triangular part of "occurances"
98      if (i < j) {
99        int temp = i;
100        i = j;
101        j = temp;
102      }
103      var entry = occurances[i][j];
104      // extracts the occurrences of the individual bits
105      bits[0] = entry[0] + entry[2];  // i zero
106      bits[1] = entry[1] + entry[3];  // i one
107      bits[2] = entry[0] + entry[1];  // j zero
108      bits[3] = entry[2] + entry[3];  // j one
109      double total = bits[0] + bits[1];
110      // entropy of the two bits on their own
111      double separate = NegativeEntropy(bits, total);
112      // entropy of the two bits as a single unit
113      double together = NegativeEntropy(entry, total);
114      // If together there is 0 entropy, the distance is zero
115      if (together.IsAlmost(0)) {
116        return 0.0;
117      }
118      return 2 - (separate / together);
119    }
120
121
122
123    // Performs O(N^2) clustering based on the method described in:
124    // "Optimal implementations of UPGMA and other common clustering algorithms"
125    // by I. Gronau and S. Moran
126    // In the GECCO paper, Figure 2 is a simplified version of this algorithm.
127    private double[][] distances;
128    private void Rebuild() {
129      if (distances == null) {
130        distances = new double[clusters.Length * 2 - 1][];
131        for (int i = 0; i < distances.Length; i++)
132          distances[i] = new double[clusters.Length * 2 - 1];
133      }
134
135
136      // Keep track of which clusters have not been merged
137      var topLevel = new List<int>(length);
138      for (int i = 0; i < length; i++)
139        topLevel.Add(i);
140
141      bool[] useful = new bool[clusters.Length];
142      for (int i = 0; i < useful.Length; i++)
143        useful[i] = true;
144
145      // Store the distances between all clusters
146      for (int i = 1; i < length; i++) {
147        for (int j = 0; j < i; j++) {
148          distances[i][j] = EntropyDistance(clusters[i][0], clusters[j][0]);
149          // make it symmetric
150          distances[j][i] = distances[i][j];
151        }
152      }
153      // Each iteration we add some amount to the path, and remove the last
154      // two elements.  This keeps track of how much of usable is in the path.
155      int end_of_path = 0;
156
157      // build all clusters of size greater than 1
158      for (int index = length; index < clusters.Length; index++) {
159        // Shuffle everything not yet in the path
160        topLevel.ShuffleInPlace(rand, end_of_path, topLevel.Count - 1);
161
162        // if nothing in the path, just add a random usable node
163        if (end_of_path == 0) {
164          end_of_path = 1;
165        }
166        while (end_of_path < topLevel.Count) {
167          // last node in the path
168          int final = topLevel[end_of_path - 1];
169
170          // best_index stores the location of the best thing in the top level
171          int best_index = end_of_path;
172          double min_dist = distances[final][topLevel[best_index]];
173          // check all options which might be closer to "final" than "topLevel[best_index]"
174          for (int option = end_of_path + 1; option < topLevel.Count; option++) {
175            if (distances[final][topLevel[option]] < min_dist) {
176              min_dist = distances[final][topLevel[option]];
177              best_index = option;
178            }
179          }
180          // If the current last two in the path are minimally distant
181          if (end_of_path > 1 && min_dist >= distances[final][topLevel[end_of_path - 2]]) {
182            break;
183          }
184
185          // move the best to the end of the path
186          topLevel.Swap(end_of_path, best_index);
187          end_of_path++;
188        }
189        // Last two elements in the path are the clusters to join
190        int first = topLevel[end_of_path - 2];
191        int second = topLevel[end_of_path - 1];
192
193        // Only keep a cluster if the distance between the joining clusters is > zero
194        bool keep = !distances[first][second].IsAlmost(0.0);
195        useful[first] = keep;
196        useful[second] = keep;
197
198        // create the new cluster
199        clusters[index] = clusters[first].Concat(clusters[second]).ToList();
200        // Calculate distances from all clusters to the newly created cluster
201        int i = 0;
202        int end = topLevel.Count - 1;
203        while (i <= end) {
204          int x = topLevel[i];
205          // Moves 'first' and 'second' to after "end" in topLevel
206          if (x == first || x == second) {
207            topLevel.Swap(i, end);
208            end--;
209            continue;
210          }
211          // Use the previous distances to calculate the joined distance
212          double first_distance = distances[first][x];
213          first_distance *= clusters[first].Count;
214          double second_distance = distances[second][x];
215          second_distance *= clusters[second].Count;
216          distances[x][index] = ((first_distance + second_distance)
217              / (clusters[first].Count + clusters[second].Count));
218          // make it symmetric
219          distances[index][x] = distances[x][index];
220          i++;
221        }
222
223        // Remove first and second from the path
224        end_of_path -= 2;
225        topLevel.RemoveAt(topLevel.Count - 1);
226        topLevel[topLevel.Count - 1] = index;
227      }
228      // Extract the useful clusters
229      clusterOrdering.Clear();
230      // Add all useful clusters. The last one is never useful.
231      for (int i = 0; i < useful.Length - 1; i++) {
232        if (useful[i]) clusterOrdering.Add(i);
233      }
234
235      // Shuffle before sort to ensure ties are broken randomly
236      clusterOrdering.ShuffleInPlace(rand);
237      clusterOrdering = clusterOrdering.OrderBy(i => clusters[i].Count).ToList();
238    }
239
240    public IEnumerable<List<int>> Clusters {
241      get {
242        // Just in time rebuilding
243        if (rebuildRequired) Rebuild();
244        foreach (var index in clusterOrdering) {
245          // Send out the clusters in the desired order
246          yield return clusters[index];
247        }
248      }
249    }
250  }
251}
Note: See TracBrowser for help on using the repository browser.