TrakEM2_/src/main/java/ini/trakem2/analysis/Centrality.java

   1 package ini.trakem2.analysis;
   2
   3 import java.util.ArrayList;
   4 import java.util.Collection;
   5 import java.util.HashMap;
   6 import java.util.HashSet;
   7 import java.util.Iterator;
   8 import java.util.LinkedList;
   9 import java.util.List;
  10 import java.util.Map;
  11 import java.util.Set;
  12
  13 /** Following pseudo-code from:
  14  *  Ulrik Brandes. 2001. A faster algorithm for betweenness centrality.
  15  *  Journal of Mathematical Sociology 25(2):163-177, (2001). */
  16 public class Centrality {
  17
  18         /** Like @see compute, but operates on a copy of all Vertex instances,
  19          * and return a Collection with the same order as @param vs. */
  20         static public final<T> ArrayList<Vertex<T>> safeCompute(final Collection<Vertex<T>> vs) {
  21                 final ArrayList<Vertex<T>> copies = Vertex.clone(vs);
  22                 compute(copies);
  23                 return copies;
  24         }
  25
  26         /** Computes betweenness centrality of each vertex in the collection of vertices @param vs
  27          *  where all vertices are part of the same graph.
  28          *  Assumes the internal variables of each Vertex are reset to its initialization values,
  29          *  and that the neighbors array is proper as well.
  30          *  When done, the centrality of each Vertex is set in the homonymous Vertex field. */
  31         static public final<T> void compute(final Collection<Vertex<T>> vs) {
  32                 if (1 == vs.size()) {
  33                         return;
  34                 }
  35                 final LinkedList<Vertex<T>> stack = new LinkedList<Vertex<T>>();
  36                 final LinkedList<Vertex<T>> queue = new LinkedList<Vertex<T>>();
  37
  38                 for (final Vertex<T> s : vs) {
  39                         // Reset all:
  40                         for (final Vertex<T> t : vs) {
  41                                 t.predecessors.clear();
  42                                 t.sigma = 0;
  43                                 t.d = -1;
  44                         }
  45                         // Prepare s:
  46                         s.sigma = 1;
  47                         s.d = 0;
  48                         queue.add(s);
  49                         //
  50                         while (!queue.isEmpty()) {
  51                                 final Vertex<T> v = queue.removeFirst();
  52                                 stack.add(v);
  53                                 for (final Vertex<T> w : v.neighbors) {
  54                                         // w found for the first time?
  55                                         if (-1 == w.d) {
  56                                                 queue.add(w);
  57                                                 w.d = v.d + 1;
  58                                         }
  59                                         // shortest path to w via v?
  60                                         if (w.d == v.d + 1) {
  61                                                 w.sigma += v.sigma; // sigma is always 1 in a tree
  62                                                 w.predecessors.add(v);
  63                                         }
  64                                 }
  65                         }
  66                         for (final Vertex<T> v : vs) {
  67                                 v.delta = 0;
  68                         }
  69                         while (!stack.isEmpty()) {
  70                                 final Vertex<T> w = stack.removeLast();
  71                                 for (final Vertex<T> v : w.predecessors) {
  72                                         v.delta += (v.sigma / (float)w.sigma) * (1 + w.delta); // sigma is always 1 in a tree
  73                                 }
  74                                 if (w != s) {
  75                                         w.centrality += w.delta;
  76                                 }
  77                         }
  78                 }
  79         }
  80
  81         /** Find the chain of nodes, over branch points if necessary, that have not yet been processed. */
  82         static private final <T> List<Vertex<T>> findChain(
  83                         final Vertex<T> origin,
  84                         final Vertex<T> parent,
  85                         final Set<Vertex<T>> processed) {
  86                 final ArrayList<Vertex<T>> chain = new ArrayList<Vertex<T>>();
  87                 chain.add(origin);
  88
  89                 Vertex<T> o = origin;
  90                 Vertex<T> p = parent;
  91
  92                 A: while (true) {
  93                         if (1 == o.neighbors.size()) { break; }
  94                         Vertex<T> o2 = o;
  95                         for (final Vertex<T> v : o.neighbors) {
  96                                 if (v == p || processed.contains(v) || v.centrality < p.centrality) {
  97                                         continue;
  98                                 }
  99                                 chain.add(v);
 100                                 p = o;
 101                                 o = v;
 102                                 continue A; // there can only be one unexplored path
 103                         }
 104                         if (o2 == o) break;
 105                 }
 106                 return chain;
 107         }
 108
 109         static public final<T> ArrayList<EtchingStep<T>> branchWise(final Collection<Vertex<T>> vs_, final int etching_multiplier) {
 110
 111                 // Copy
 112                 final Set<Vertex<T>> vs = new HashSet<Vertex<T>>(Vertex.clone(vs_));
 113
 114                 // Map of original vertex instances
 115                 final HashMap<T,Vertex<T>> m = new HashMap<T,Vertex<T>>();
 116                 for (final Vertex<T> v : vs_) {
 117                         m.put(v.data, v);
 118                 }
 119
 120                 // A set of the remaining branch vertices
 121                 final Set<Vertex<T>> branch_vertices = new HashSet<Vertex<T>>();
 122                 for (final Vertex<T> v : vs) {
 123                         if (v.getNeighborCount() > 2) branch_vertices.add(v);
 124                 }
 125
 126                 // Map of the count of remaining branch vertices and vertices removed at that step
 127                 final ArrayList<EtchingStep<T>> steps = new ArrayList<EtchingStep<T>>();
 128
 129                 final HashSet<Vertex<T>> processed = new HashSet<Vertex<T>>();
 130
 131                 // TODO the node centrality needs to be computed only once!
 132                 // TODO the centrality value is thrown out!
 133
 134                 while (vs.size() > 0) {
 135                         // Reset all internal vars related to computing centrality
 136                         for (final Vertex<T> v : vs) {
 137                                 v.reset();
 138                         }
 139                         // Recompute centrality for the now smaller graph
 140                         Centrality.compute(vs);
 141
 142                         // Remove all vertices whose centrality falls below a certain threshold
 143                         final HashSet<Vertex<T>> removed = new HashSet<Vertex<T>>();
 144                         final int vs_size = vs.size();
 145                         for (final Iterator<Vertex<T>> it = vs.iterator(); it.hasNext(); ) {
 146                                 final Vertex<T> v = it.next();
 147                                 if (v.centrality < etching_multiplier * vs_size) {
 148                                         it.remove();
 149                                         removed.add(v);
 150                                 }
 151                         }
 152
 153                         // Determine which branches have been removed at this etching step
 154                         final Set<Collection<Vertex<T>>> etched_branches = new HashSet<Collection<Vertex<T>>>();
 155                         for (final Vertex<T> r : removed) {
 156                                 for (final Vertex<T> w : r.neighbors) {
 157                                         if (w.isBranching() && w.centrality >= r.centrality) {
 158                                                 final List<Vertex<T>> branch = findChain(m.get(r.data), m.get(w.data), processed);
 159                                                 processed.addAll(branch);
 160                                                 etched_branches.add(branch);
 161                                         }
 162                                 }
 163                         }
 164
 165                         if (etched_branches.size() > 0) {
 166                                 steps.add(new EtchingStep<T>(branch_vertices.size(), etched_branches));
 167                         }
 168
 169                         // Fix neighbors of remaining vertices
 170                         for (final Vertex<T> v : removed) {
 171                                 for (final Iterator<Vertex<T>> it = v.neighbors.iterator(); it.hasNext(); ) {
 172                                         final Vertex<T> w = it.next();
 173                                         if (removed.contains(w)) {
 174                                                 it.remove();
 175                                                 continue;
 176                                         }
 177                                         // to the vertex that remains, remove the vertex v from its neighbors:
 178                                         w.neighbors.remove(v);
 179                                 }
 180                         }
 181                         // Remove branch vertices which aren't branch vertices anymore
 182                         for (final Iterator<Vertex<T>> it = branch_vertices.iterator(); it.hasNext(); ) {
 183                                 final Vertex<T> v = it.next();
 184                                 if (v.neighbors.size() < 3 || removed.contains(v)) {
 185                                         it.remove();
 186                                 }
 187                         }
 188                         if (0 == branch_vertices.size()) {
 189                                 break;
 190                         }
 191                 }
 192
 193                 for (final EtchingStep<T> es : steps) {
 194                         for (final Collection<Vertex<T>> b : es.branches) {
 195                                 List<Vertex<T>> b2 = new ArrayList<Vertex<T>>(b);
 196                                 b.clear();
 197                                 for (final Vertex<T> v : b2) {
 198                                         b.add(m.get(v.data));
 199                                 }
 200                         }
 201                 }
 202
 203                 // The last, central branch
 204                 Set<T> done = new HashSet<T>();
 205                 for (final Vertex<T> v : processed) {
 206                         done.add(v.data);
 207                 }
 208                 HashMap<T,Vertex<T>> last = new HashMap<T,Vertex<T>>(m);
 209                 last.keySet().removeAll(done);
 210                 HashSet<Collection<Vertex<T>>> bs = new HashSet<Collection<Vertex<T>>>();
 211                 bs.add(last.values());
 212                 steps.add(new EtchingStep<T>(0, bs));
 213
 214                 return steps;
 215         }
 216
 217         /** An entry of the number of remaining branch vertices versus
 218          *  the set of branches (each as a List<Vertex<T>>) removed. */
 219         static public class EtchingStep<T> implements Map.Entry<Integer,Set<Collection<Vertex<T>>>> {
 220                 public int remaining_branch_vertices = 0;
 221                 public Set<Collection<Vertex<T>>> branches = null;
 222                 public EtchingStep(final int remaining_branch_vertices, final Set<Collection<Vertex<T>>> branches) {
 223                         this.remaining_branch_vertices = remaining_branch_vertices;
 224                         this.branches = branches;
 225                 }
 226                 @Override
 227                 public Integer getKey() {
 228                         return remaining_branch_vertices;
 229                 }
 230                 @Override
 231                 public Set<Collection<Vertex<T>>> getValue() {
 232                         return branches;
 233                 }
 234                 @Override
 235                 public Set<Collection<Vertex<T>>> setValue(final Set<Collection<Vertex<T>>> value) {
 236                         throw new UnsupportedOperationException();
 237                 }
 238         }
 239 }