Free cookie consent management tool by TermsFeed Policy Generator

source: branches/Parameter-less Population Pyramid/HeuristicLab.Algorithms.ParameterlessPopulationPyramid/3.3/LinkageTree.cs @ 11672

Last change on this file since 11672 was 11672, checked in by bgoldman, 9 years ago

#2282 Commenting, added some basic unit tests for P3.

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