Author: scooter
Date: 2010-11-18 17:48:12 -0800 (Thu, 18 Nov 2010)
New Revision: 22920
Added:
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/KCluster.java
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/RunSCPS.java
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/SCPSCluster.java
Modified:
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/ClusterMaker.java
Log:
Add first pass at SCPS clustering
Modified:
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/ClusterMaker.java
===================================================================
---
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/ClusterMaker.java
2010-11-19 01:24:27 UTC (rev 22919)
+++
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/ClusterMaker.java
2010-11-19 01:48:12 UTC (rev 22920)
@@ -66,6 +66,7 @@
import clusterMaker.algorithms.MCODE.MCODECluster;
import clusterMaker.algorithms.glay.GLayCluster;
import clusterMaker.algorithms.ConnectedComponents.ConnectedComponentsCluster;
+import clusterMaker.algorithms.SCPS.SCPSCluster;
// import clusterMaker.algorithms.QT.QTCluster;
// import clusterMaker.algorithms.Spectral.SpectralCluster;
// import clusterMaker.algorithms.CP.CPCluster;
@@ -109,15 +110,17 @@
addClusterAlgorithm(menu, new HierarchicalCluster());
addClusterAlgorithm(menu, new KMeansCluster());
// addClusterAlgorithm(menu, new QTCluster());
+ menu.addSeparator();
+ addClusterAlgorithm(menu, new APCluster());
+ addClusterAlgorithm(menu, new ConnectedComponentsCluster());
+ addClusterAlgorithm(menu, new GLayCluster());
addClusterAlgorithm(menu, new MCODECluster());
addClusterAlgorithm(menu, new MCLCluster());
// addClusterAlgorithm(menu, new SpectralCluster());
// addClusterAlgorithm(menu, new CPCluster());
- addClusterAlgorithm(menu, new APCluster());
- addClusterAlgorithm(menu, new GLayCluster());
// addClusterAlgorithm(menu, new FORCECluster());
+ addClusterAlgorithm(menu, new SCPSCluster());
addClusterAlgorithm(menu, new TransClustCluster());
- addClusterAlgorithm(menu, new ConnectedComponentsCluster());
// addClusterAlgorithm(new HOPACHCluster());
menu.addSeparator();
Added:
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/KCluster.java
===================================================================
---
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/KCluster.java
(rev 0)
+++
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/KCluster.java
2010-11-19 01:48:12 UTC (rev 22920)
@@ -0,0 +1,526 @@
+/* vim: set ts=2: */
+/**
+ * Copyright (c) 2008 The Regents of the University of California.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions, and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions, and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * 3. Redistributions must acknowledge that this software was
+ * originally developed by the UCSF Computer Graphics Laboratory
+ * under support by the NIH National Center for Research Resources,
+ * grant P41-RR01081.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER "AS IS" AND ANY
+ * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
+ * OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
+ * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
+ * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
+ * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ */
+package clusterMaker.algorithms.SCPS;
+
+import java.awt.GridLayout;
+import java.lang.Math;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Date;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Random;
+import javax.swing.JPanel;
+
+
+
+// clusterMaker imports
+import cern.colt.matrix.doublealgo.Statistic;
+
+import cern.colt.matrix.DoubleFactory2D;
+import cern.colt.matrix.DoubleMatrix1D;
+import cern.colt.matrix.DoubleMatrix2D;
+import cern.colt.matrix.impl.SparseDoubleMatrix2D;
+
+import cern.colt.list.IntArrayList;
+import cern.colt.list.DoubleArrayList;
+
+
+public class KCluster {
+
+
+
+
+ // Random seeds for uniform()
+ static int seed1 = 0;
+ static int seed2 = 0;
+ static boolean debug = false;
+
+
+
+
+ //perform Kmeans clustering on columns of matrix, the size of whose
row is K
+ public static int kmeans(int nClusters, int nIterations, DoubleMatrix2D
matrix,
+ int[] clusterID) {
+
+ int nelements = matrix.columns();
+ int ifound = 1;
+
+ int[] tclusterid = new int[nelements];
+
+ int[] saved = new int[nelements];
+
+ int[] mapping = new int[nClusters];
+ int[] counts = new int[nClusters];
+
+ double error = Double.MAX_VALUE;
+
+ // This matrix will store the centroid data. KxK matrix
+ SparseDoubleMatrix2D cData = new
SparseDoubleMatrix2D(nClusters, nClusters);
+
+ // Outer initialization
+ if (nIterations <= 1) {
+ for (int i=0; i < clusterID.length; i++) {
+ tclusterid[i] = clusterID[i];
+ }
+ nIterations = 1;
+ } else {
+ for (int i = 0; i < nelements; i++)
+ clusterID[i] = 0;
+ }
+
+ int iteration = 0;
+ do {
+ double total = Double.MAX_VALUE;
+ int counter = 0;
+ int period = 10;
+
+ // Randomly assign elements to clusters
+ if (nIterations != 0) randomAssign(nClusters,
nelements, tclusterid);
+
+
+ // Initialize
+ for (int i = 0; i < nClusters; i++) counts[i] = 0;
+ for (int i = 0; i < nelements; i++)
counts[tclusterid[i]]++;
+
+ while (true) {
+ double previous = total;
+ total = 0.0;
+ if (counter % period == 0) // Save the current
cluster assignments
+ {
+ for (int i = 0; i < nelements; i++)
+ saved[i] = tclusterid[i];
+ if (period < Integer.MAX_VALUE / 2)
+ period *= 2;
+ }
+ counter++;
+
+ //assign centroids as to maximize orthogonality
on first iteration of kmeans
+ if(iteration == 0)
+ selectCentroidsOrthogonally(nClusters,
matrix, cData);
+
+ // Find the centeroids based on cluster means
on all other iterations
+ else
+ getClusterMeans(nClusters, matrix, cData,
tclusterid);
+
+
+ for (int i = 0; i < nelements; i++) {
+ // Calculate the distances
+ double distance;
+ int k = tclusterid[i];
+ if (counts[k]==1) continue;
+
+ // Get the distance
+
+ distance =
getDistance(i,k,matrix,cData);
+
+ for (int j = 0; j < nClusters; j++) {
+ double tdistance;
+ if (j==k) continue;
+
+ tdistance =
getDistance(i,k,matrix,cData);
+ if (tdistance < distance)
+ {
+ distance = tdistance;
+ counts[tclusterid[i]]--;
+ tclusterid[i] = j;
+ counts[j]++;
+ }
+ }
+ total += distance;
+ }
+
+ if (total>=previous) break;
+ /* total>=previous is FALSE on some machines even if total and previous
+ * are bitwise identical. */
+ int i;
+ for (i = 0; i < nelements; i++)
+ if (saved[i]!=tclusterid[i]) break;
+ if (i==nelements)
+ break; /* Identical solution found; break out of this loop */
+ }
+
+ if (nIterations<=1)
+ { error = total;
+ break;
+ }
+
+ for (int i = 0; i < nClusters; i++) mapping[i] = -1;
+
+ int element = 0;
+ for (element = 0; element < nelements; element++)
+ {
+ int j = tclusterid[element];
+ int k = clusterID[element];
+ if (mapping[k] == -1)
+ mapping[k] = j;
+ else if (mapping[k] != j)
+ {
+ if (total < error)
+ {
+ ifound = 1;
+ error = total;
+ // System.out.println("Mapping
tclusterid to clusterid");
+ for (int i = 0; i < nelements; i++) clusterID[i] =
tclusterid[i];
+ }
+ break;
+ }
+ }
+ if (element==nelements) ifound++; /* break statement not encountered */
+ } while (++iteration < nIterations);
+
+
+ return ifound;
+ }
+
+
+
+ //Assign centroids to rows with maximal orthogonality
+ private static void selectCentroidsOrthogonally(int nClusters,
DoubleMatrix2D data, DoubleMatrix2D cdata){
+
+ //number of centroids allready set
+ int centroid_count = 0;
+
+ //centroid assigned element == 1 if centroid assigned, zero
otherwise
+ int[] centroid_assigned = new int[nClusters];
+ for (int i = 0; i < nClusters; i++){centroid_assigned[i] = 0;}
+
+ //randomly select first centroid
+ int centroid_id = binomial(nClusters-1,1.0/(double)(nClusters-1));
+
+ for(int i = 0; i < nClusters; i++){
+
+ //update centroid assinged array
+ centroid_assigned[centroid_id] = 1;
+
+ DoubleMatrix1D centroid = data.viewRow(centroid_id);
+
+ //copy newly selected centroid from data matrix
+ for(int j = 0; j < centroid.size(); j++)
+ cdata.set(centroid_count,j,centroid.get(j));
+
+ centroid_count++;
+
+ if(centroid_count == nClusters)
+ break;
+
+ int min_cosine = 10000;
+ int new_centroid_id = -1;
+
+ //loop through rows of data matrix, seach for next centroid,
which will minimize the cosine angle (dot product) with current centroid
+ for(int j = 0; j < data.rows(); i++){
+
+ //ignore centroids that allready have been set
+ if(centroid_assigned[j] == 1)
+ continue;
+
+ double cosine = centroid.zDotProduct(data.viewRow(j));
+
+ //if new cosine value < min cosine value, update min_cosine
and new_centroid_id to reflect taht
+ if(min_cosine > cosine){
+
+ cosine = min_cosine;
+ new_centroid_id = j;
+ }
+ }
+
+
+ centroid_id = new_centroid_id;
+
+ }
+ }
+
+
+
+
+ private static void getClusterMeans(int nClusters, DoubleMatrix2D data,
DoubleMatrix2D cdata, int[] clusterid) {
+
+ double[][]cmask = new double[nClusters][nClusters];
+
+ for (int i = 0; i < nClusters; i++) {
+ for (int j = 0; j < data.columns(); j++) {
+ //cdata.setValue(i, j, null);
+ cmask[i][j] = 0.0;
+ }
+ }
+
+ for (int k = 0; k < data.rows(); k++) {
+
+ //cetroid id
+ int i = clusterid[k];
+
+
+
+
+
+ for (int j = 0; j < data.columns(); j++) {
+
+
+ if (data.get(k,j) != 0) {
+ double cValue = 0.0;
+ double dataValue = data.get(k,j);
+ if (cdata.get(i,j) != 0) {
+ cValue = cdata.get(i,j);
+ }
+ cdata.set(i,j, cValue+dataValue);
+ cmask[i][j] = cmask[i][j] + 1.0;
+ }
+ }
+ }
+ for (int i = 0; i < nClusters; i++) {
+ for (int j = 0; j < data.columns(); j++) {
+ if (cmask[i][j] > 0.0) {
+ double cData = cdata.get(i,j) / cmask[i][j];
+ cdata.set(i,j,cData);
+ }
+ }
+ }
+ }
+
+
+ //get Euclidian Distance between two rows of matrix. Use colt Euclidian
Distance function
+ private static double getDistance(int row1_id, int row2_id, DoubleMatrix2D
matrix, DoubleMatrix2D cdata){
+
+
+ DoubleMatrix1D row1 = matrix.viewRow(row1_id);
+ DoubleMatrix1D row2 = cdata.viewRow(row2_id);
+
+ return Statistic.EUCLID.apply(row1,row2);
+ }
+
+
+ private static void randomAssign (int nClusters, int nElements, int[]
clusterID) {
+ int n = nElements - nClusters;
+ int k = 0;
+ int i = 0;
+ for (i = 0; i < nClusters-1; i++) {
+ double p = 1.0/(nClusters-1);
+ int j = binomial(n, p);
+ n -= j;
+ j += k+1; // Assign at least one element to cluster i
+ for (;k<j; k++) clusterID[k] = i;
+ }
+ // Assign the remaining elements to the last cluster
+ for (; k < nElements; k++) clusterID[k] = i;
+
+ // Create a random permutation of the cluster assignments
+ for (i = 0; i < nElements; i++) {
+ int j = (int) (i + (nElements-i)*uniform());
+ k = clusterID[j];
+ clusterID[j] = clusterID[i];
+ clusterID[i] = k;
+ }
+ }
+
+ // Debug version of "randomAssign" that isn't random
+ private static void debugAssign (int nClusters, int nElements, int[]
clusterID) {
+ for (int element = 0; element < nElements; element++) {
+ clusterID[element] = element%nClusters;
+ }
+ }
+
+ /**
+ * This routine generates a random number between 0 and n inclusive,
following
+ * the binomial distribution with probability p and n trials. The
routine is
+ * based on the BTPE algorithm, described in:
+ *
+ * Voratas Kachitvichyanukul and Bruce W. Schmeiser:
+ * Binomial Random Variate Generation
+ * Communications of the ACM, Volume 31, Number 2, February 1988, pages
216-222.
+ *
+ * @param p The probability of a single event. This should be less
than or equal to 0.5.
+ * @param n The number of trials
+ * @return An integer drawn from a binomial distribution with
parameters (p, n).
+ */
+
+ private static int binomial (int n, double p) {
+ double q = 1 - p;
+ if (n*p < 30.0) /* Algorithm BINV */
+ {
+ double s = p/q;
+ double a = (n+1)*s;
+ double r = Math.exp(n*Math.log(q)); /* pow() causes a
crash on AIX */
+ int x = 0;
+ double u = uniform();
+ while(true)
+ {
+ if (u < r) return x;
+ u-=r;
+ x++;
+ r *= (a/x)-s;
+ }
+ }
+ else /* Algorithm BTPE */
+ { /* Step 0 */
+ double fm = n*p + p;
+ int m = (int) fm;
+ double p1 = Math.floor(2.195*Math.sqrt(n*p*q) -4.6*q) +
0.5;
+ double xm = m + 0.5;
+ double xl = xm - p1;
+ double xr = xm + p1;
+ double c = 0.134 + 20.5/(15.3+m);
+ double a = (fm-xl)/(fm-xl*p);
+ double b = (xr-fm)/(xr*q);
+ double lambdal = a*(1.0+0.5*a);
+ double lambdar = b*(1.0+0.5*b);
+ double p2 = p1*(1+2*c);
+ double p3 = p2 + c/lambdal;
+ double p4 = p3 + c/lambdar;
+ while (true)
+ { /* Step 1 */
+ int y;
+ int k;
+ double u = uniform();
+ double v = uniform();
+ u *= p4;
+ if (u <= p1) return (int)(xm-p1*v+u);
+ /* Step 2 */
+ if (u > p2)
+ { /* Step 3 */
+ if (u > p3)
+ { /* Step 4 */
+ y =
(int)(xr-Math.log(v)/lambdar);
+ if (y > n) continue;
+ /* Go to step 5 */
+ v = v*(u-p3)*lambdar;
+ }
+ else
+ {
+ y =
(int)(xl+Math.log(v)/lambdal);
+ if (y < 0) continue;
+ /* Go to step 5 */
+ v = v*(u-p2)*lambdal;
+ }
+ }
+ else
+ {
+ double x = xl + (u-p1)/c;
+ v = v*c + 1.0 - Math.abs(m-x+0.5)/p1;
+ if (v > 1) continue;
+ /* Go to step 5 */
+ y = (int)x;
+ }
+ /* Step 5 */
+ /* Step 5.0 */
+ k = Math.abs(y-m);
+ if (k > 20 && k < 0.5*n*p*q-1.0)
+ { /* Step 5.2 */
+ double rho = (k/(n*p*q))*((k*(k/3.0 +
0.625) + 0.1666666666666)/(n*p*q)+0.5);
+ double t = -k*k/(2*n*p*q);
+ double A = Math.log(v);
+ if (A < t-rho) return y;
+ else if (A > t+rho) continue;
+ else
+ { /* Step 5.3 */
+ double x1 = y+1;
+ double f1 = m+1;
+ double z = n+1-m;
+ double w = n-y+1;
+ double x2 = x1*x1;
+ double f2 = f1*f1;
+ double z2 = z*z;
+ double w2 = w*w;
+ if (A > xm * Math.log(f1/x1) +
(n-m+0.5)*Math.log(z/w)
+ +
(y-m)*Math.log(w*p/(x1*q))
+ +
(13860.-(462.-(132.-(99.-140./f2)/f2)/f2)/f2)/f1/166320.
+ +
(13860.-(462.-(132.-(99.-140./z2)/z2)/z2)/z2)/z/166320.
+ +
(13860.-(462.-(132.-(99.-140./x2)/x2)/x2)/x2)/x1/166320.
+ +
(13860.-(462.-(132.-(99.-140./w2)/w2)/w2)/w2)/w/166320.)
+ continue;
+ return y;
+ }
+ }
+ else
+ { /* Step 5.1 */
+ int i;
+ double s = p/q;
+ double aa = s*(n+1);
+ double f = 1.0;
+ for (i = m; i < y; f *= (aa/(++i)-s));
+ for (i = y; i < m; f /= (aa/(++i)-s));
+ if (v > f) continue;
+ return y;
+ }
+ }
+ }
+ }
+
+ /**
+ * This routine returns a uniform random number between 0.0 and 1.0.
Both 0.0
+ * and 1.0 are excluded. This random number generator is described in:
+ *
+ * Pierre l'Ecuyer
+ * Efficient and Portable Combined Random Number Generators
+ * Communications of the ACM, Volume 31, Number 6, June 1988, pages
742-749,774.
+ *
+ * The first time this routine is called, it initializes the random
number
+ * generator using the current time. First, the current epoch time in
seconds is
+ * used as a seed for the random number generator in the C library. The
first two
+ * random numbers generated by this generator are used to initialize
the random
+ * number generator implemented in this routine.
+ *
+ * NOTE: how different is this from Java's Math.random() or
Random.nextDouble()?
+ *
+ * @return A double-precison number between 0.0 and 1.0.
+ */
+ private static double uniform() {
+ int z;
+ int m1 = 2147483563;
+ int m2 = 2147483399;
+ double scale = 1.0/m1;
+
+ if (seed1==0 || seed2==0) /* initialize */
+ {
+ Date date = new Date();
+ int initseed = (int) date.getTime();
+ Random r = new Random(initseed);
+ seed1 = r.nextInt();
+ seed2 = r.nextInt();
+ }
+
+ do
+ {
+ int k;
+ k = seed1/53668;
+ seed1 = 40014*(seed1-k*53668)-k*12211;
+ if (seed1 < 0) seed1+=m1;
+ k = seed2/52774;
+ seed2 = 40692*(seed2-k*52774)-k*3791;
+ if(seed2 < 0) seed2+=m2;
+ z = seed1-seed2;
+ if(z < 1) z+=(m1-1);
+ } while (z==m1); /* To avoid returning 1.0 */
+
+ return z*scale;
+ }
+}
Added:
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/RunSCPS.java
===================================================================
---
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/RunSCPS.java
(rev 0)
+++
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/RunSCPS.java
2010-11-19 01:48:12 UTC (rev 22920)
@@ -0,0 +1,418 @@
+/**
+ * Copyright (c) 2010 The Regents of the University of California.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions, and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions, and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * 3. Redistributions must acknowledge that this software was
+ * originally developed by the UCSF Computer Graphics Laboratory
+ * under support by the NIH National Center for Research Resources,
+ * grant P41-RR01081.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER "AS IS" AND ANY
+ * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
+ * OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
+ * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
+ * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
+ * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ */
+package clusterMaker.algorithms.SCPS;
+
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Collection;
+import java.util.Comparator;
+import java.util.HashMap;
+import java.util.Iterator;
+import java.util.List;
+import java.util.Map;
+import java.util.Set;
+import java.lang.Math;
+
+import cytoscape.Cytoscape;
+import cytoscape.CyNetwork;
+import cytoscape.CyEdge;
+import cytoscape.CyNode;
+import cytoscape.data.CyAttributes;
+import cytoscape.groups.CyGroup;
+import cytoscape.groups.CyGroupManager;
+import cytoscape.logger.CyLogger;
+import cytoscape.task.TaskMonitor;
+
+import clusterMaker.algorithms.NodeCluster;
+import clusterMaker.algorithms.DistanceMatrix;
+import clusterMaker.algorithms.ClusterResults;
+
+import cern.colt.function.IntIntDoubleFunction;
+import cern.colt.matrix.DoubleFactory2D;
+import cern.colt.matrix.DoubleMatrix1D;
+import cern.colt.matrix.DoubleMatrix2D;
+import cern.colt.matrix.impl.SparseDoubleMatrix2D;
+import cern.colt.matrix.linalg.EigenvalueDecomposition;
+import cern.colt.matrix.linalg.SingularValueDecomposition;
+
+
+
+import cern.colt.matrix.linalg.Algebra;
+import cern.colt.list.IntArrayList;
+import cern.colt.list.DoubleArrayList;
+
+import clusterMaker.algorithms.SCPS.KCluster;
+
+public class RunSCPS {
+
+
+
+ private List<CyNode> nodes;
+ private List<CyEdge> edges;
+ private boolean canceled = false;
+ private CyLogger logger;
+ public final static String GROUP_ATTRIBUTE = "__SCPSGroups";
+ protected int clusterCount = 0;
+ private boolean createMetaNodes = false;
+ private DistanceMatrix distanceMatrix = null;
+ private DoubleMatrix2D matrix = null;
+ private boolean debug = false;
+
+ private double epsilon;
+ private int kvalue;
+ private int rnumber;
+ private DoubleMatrix2D LMat;
+
+
+ private HashMap<Integer, NodeCluster> clusterMap;
+
+ private HashMap<Integer,Integer> index2indexMap;
+
+
+
+ public RunSCPS(DistanceMatrix dMat, double epsilon, int kvalue, int
rnumber, CyLogger logger )
+ {
+ this.distanceMatrix = dMat;
+ this.epsilon = epsilon;
+ this.kvalue = kvalue;
+ this.rnumber = rnumber;
+
+ this.logger = logger;
+ this.clusterMap = new HashMap();
+ this.clusterCount = 0;
+ nodes = distanceMatrix.getNodes();
+ edges = distanceMatrix.getEdges();
+ this.matrix = distanceMatrix.getDistanceMatrix();
+
+
+ //index2index maps indices of filtered nodes in new uMatrix to
the the original in the complete, unfiltered matrix
+ this.index2indexMap = new HashMap();
+ }
+
+
+ public void halt () { canceled = true; }
+
+ public List<NodeCluster> run(TaskMonitor monitor)
+ {
+
+ monitor.setStatus("Formatting Matrix Data");
+ DoubleMatrix2D sMat = getSMat(this.distanceMatrix);
+ DoubleMatrix2D LMat = getLMat(sMat);
+ monitor.setStatus("Calculating Eigenvalues");
+ DoubleMatrix2D uMat = getUMat(LMat);
+ monitor.setStatus("Running kmeans clustering");
+ doKMeansClustering(uMat);
+
+ //clusterMap calculated in getSMat and doKMeansClustering steps.
Simply return the results
+ return new ArrayList(this.clusterMap.values());
+
+
+ }
+
+ //map index2index key-value pair
+ public void setMap(int old_index, int new_index){
+
+ this.index2indexMap.put(new Integer(new_index), new
Integer(old_index));
+ }
+
+ //return value of old_index in index2index map
+ public int getMap(int new_index){
+
+ return ((Integer)this.index2indexMap.get(new
Integer(new_index))).intValue();
+ }
+
+ //Get Connected Components, cluster all components <= |5|, and connect
the remaining components with random lowscoring edges
+ public DoubleMatrix2D getSMat(DistanceMatrix distanceMatrix){
+
+ //Matrix prior to filtration modification
+ DoubleMatrix2D unfiltered_mat = distanceMatrix.getDistanceMatrix();
+
+ //Size of newly created Umat after filtering of small components
+ int sMat_rows = 0;
+
+ HashMap<Integer, List<CyNode>> filtered_cmap = new HashMap();
+
+ //Connected Componets
+ Map<Integer, List<CyNode>> cMap =
distanceMatrix.findConnectedComponents();
+
+ IntArrayList rowList = null;
+ IntArrayList columnList = null;
+ DoubleArrayList valueList = null;
+
+ //Iterate through connected components
+ Iterator it = cMap.entrySet().iterator();
+
+ while(it.hasNext()){
+
+ Map.Entry entry = (Map.Entry)it.next();
+ List<CyNode> component = (List<CyNode>)entry.getValue();
+
+ //Size <= 5. Automatically create cluster and increment
clusterCount. Delete component from Map.
+ if(component.size() <= 5){
+
+ NodeCluster iCluster = new NodeCluster();
+ iCluster.add(component,this.clusterCount);
+ this.clusterMap.put(new Integer(clusterCount),iCluster);
+ this.clusterCount++;
+
+
+
+ }
+
+ //iterate through components and assign them index mappings in
new uMatrix
+ else{
+
+ for(int i = 0; i < component.size(); i++){
+
+ CyNode n = component.get(i);
+
+ //set mapping of new matrix index to old index
+ setMap(this.nodes.indexOf(n), sMat_rows);
+ sMat_rows++;
+ }
+
+ }
+
+ }
+
+ SparseDoubleMatrix2D sMat = new SparseDoubleMatrix2D(sMat_rows,
sMat_rows);
+
+ //set diagnols of sMat to one
+ for(int i = 0; i < sMat_rows; i++)
+ sMat.set(i,i,1);
+
+
+ //iterate through nonzero edges. If both nodes in new index map,
transfer the edge to new matrix
+ unfiltered_mat.getNonZeros(rowList,columnList,valueList);
+
+ for(int i = 0; i<rowList.size(); i++){
+
+ int row_id = rowList.get(i);
+ int column_id = columnList.get(i);
+
+ int new_row_id = getMap(row_id);
+ int new_column_id = getMap(column_id);
+ double value = valueList.get(i);
+
+ //Set symmetrically the values in new matrix
+ if(new_row_id > -1 && new_column_id > -1)
+ {
+ sMat.set(new_row_id,new_column_id,value);
+ sMat.set(new_column_id,new_row_id,value);
+ }
+
+ }
+
+ //Normalize sMat
+ sMat.forEachNonZero(new Normalize(Double.MIN_VALUE,Double.MAX_VALUE,
1));
+
+ return sMat;
+
+ }
+
+ //Calculate negative square root of matrix using singular value
decomposition
+ //http://en.wikipedia.org/wiki/Square_root_of_a_matrix
+ public DoubleMatrix2D getNegSqrRoot(DoubleMatrix2D A){
+
+ //A = USV, where S is Diagnol Matrix
+ SingularValueDecomposition decomp = new
SingularValueDecomposition(A);
+ DoubleMatrix2D U = decomp.getU();
+ DoubleMatrix2D S = decomp.getS();
+ DoubleMatrix2D V = decomp.getV();
+
+ //S^1/2 = Square root of every value in diangol matrix
+ for(int i = 0; i < S.rows(); i++)
+ S.set(i,i,Math.pow(S.get(i,i),.5));
+
+ //A^1/2 = VS^1/2U
+ Algebra alg = new Algebra();
+ DoubleMatrix2D sqrtA = alg.mult(alg.mult(V,S),U);
+
+ //return A^-1/2
+ return alg.inverse(sqrtA);
+
+
+ }
+
+
+ // L = D^-1/2 * S * D^-1/2
+ public DoubleMatrix2D getLMat(DoubleMatrix2D sMat){
+
+ Algebra alg = new Algebra();
+ DoubleMatrix2D transDMat = getNegSqrRoot(sMat);
+
+ return alg.mult(transDMat,alg.mult(sMat,transDMat));
+ }
+
+ //D is Diagnol Matrix formed of vertex degree Dii = Sum Columns j over
row Si
+
+ public DoubleMatrix2D getDMat(DoubleMatrix2D sMat){
+
+
+ DoubleMatrix2D dMat = sMat.like();
+
+ for(int i = 0; i < sMat.rows(); i++){
+
+ //set the Diagnal (i,i) to sum of columns over row i
+ dMat.set(i,i, sMat.viewRow(i).zSum());
+ }
+
+
+ return dMat;
+
+ }
+
+ //U constructed from top K Eigenvectors of L. After construction, each
row of U is normalized to unit length.
+
+ public DoubleMatrix2D getUMat(DoubleMatrix2D LMat){
+
+ int k;
+ DoubleMatrix2D uMat;
+ double prevLamb;
+ double nextLamb;
+
+ IntArrayList indexList = null;
+ DoubleArrayList valueList = null;
+
+ EigenvalueDecomposition decomp = new
EigenvalueDecomposition(LMat);
+
+ //eigenvectors
+ DoubleMatrix2D eigenVect = decomp.getV();
+
+ //eigenvalues
+ DoubleMatrix1D eigenVal = decomp.getRealEigenvalues();
+
+ //set K. Use Epsilon to calculate K is this.kvalue not set
+ if(this.kvalue > -1)
+ k = this.kvalue;
+
+ //set K to smallest integer such that LambdaK+1/LambdaK >
epsilon
+ else{
+
+ prevLamb = eigenVal.get(0);
+
+ for(k = 1; k < eigenVal.size(); k++){
+
+ nextLamb = eigenVal.get(k);
+
+ if(nextLamb/prevLamb > this.epsilon)
+ break;
+
+ prevLamb = nextLamb;
+ }
+ }
+
+
+ //construct matrix U from first K eigenvectors
+ uMat = eigenVect.viewPart(0,0,eigenVect.rows()-1,k);
+
+ //Normalize each row of matrix U
+ for(int i = 0; i < uMat.columns(); i++){
+
+ DoubleMatrix1D row = uMat.viewRow(i);
+ double sum = row.zSum();
+ row.getNonZeros(indexList,valueList);
+
+ //normalize each Nozero value in row
+
+ for(int j = 0; j < indexList.size(); j ++){
+
+ int index = indexList.get(j);
+ double value = valueList.get(j)/sum;
+
+ uMat.set(i,index,value);
+ }
+
+ }
+
+
+ return uMat;
+ }
+
+ public void doKMeansClustering(DoubleMatrix2D uMat){
+
+ int k = uMat.rows();
+
+ int[] clusters = new int[uMat.columns()];
+
+ //do kmeans clustering
+ KCluster.kmeans(k,rnumber,uMat,clusters);
+
+ //Loop through clustering results, getting the clusters by order
+
+ for(int cluster_id = 0; cluster_id < k; cluster_id++){
+
+ NodeCluster iCluster = new NodeCluster();
+ List<CyNode> node_list = new ArrayList();
+
+ for(int j = 0; j < clusters.length; j++){
+
+ //node j in uMatrix belongs to cluster k
+ if(clusters[j] == k)
+ node_list.add(this.nodes.get(getMap(j)));
+ }
+
+ iCluster.add(node_list,this.clusterCount);
+ this.clusterMap.put(new Integer(clusterCount),iCluster);
+ this.clusterCount++;
+ }
+
+
+
+ }
+
+
+ /**
+ * Normalize normalizes a cell in the matrix
+ */
+ private class Normalize implements IntIntDoubleFunction {
+ private double maxValue;
+ private double minValue;
+ private double span;
+ private double factor;
+ private double minWeight = Double.MAX_VALUE;
+
+ public Normalize(double minValue, double maxValue, double
factor) {
+ this.maxValue = maxValue;
+ this.minValue = minValue;
+ this.factor = factor;
+ span = maxValue - minValue;
+ }
+
+ public double apply(int row, int column, double value) {
+ return ((value-minWeight)/span)*factor;
+ }
+ }
+
+
+}
+
Added:
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/SCPSCluster.java
===================================================================
---
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/SCPSCluster.java
(rev 0)
+++
csplugins/trunk/ucsf/scooter/clusterMaker/src/clusterMaker/algorithms/SCPS/SCPSCluster.java
2010-11-19 01:48:12 UTC (rev 22920)
@@ -0,0 +1,203 @@
+/* vim: set ts=2: */
+/**
+ * Copyright (c) 2008 The Regents of the University of California.
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions, and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above
+ * copyright notice, this list of conditions, and the following
+ * disclaimer in the documentation and/or other materials provided
+ * with the distribution.
+ * 3. Redistributions must acknowledge that this software was
+ * originally developed by the UCSF Computer Graphics Laboratory
+ * under support by the NIH National Center for Research Resources,
+ * grant P41-RR01081.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER "AS IS" AND ANY
+ * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+ * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE REGENTS BE LIABLE
+ * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
+ * OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
+ * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
+ * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
+ * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
+ * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *
+ */
+package clusterMaker.algorithms.SCPS;
+
+import java.awt.Dimension;
+import java.awt.GridLayout;
+import java.awt.event.ActionEvent;
+import java.awt.event.ActionListener;
+import java.util.ArrayList;
+import java.util.Arrays;
+import java.util.Comparator;
+import java.util.List;
+import javax.swing.JPanel;
+
+// Cytoscape imports
+import cytoscape.CyEdge;
+import cytoscape.CyNetwork;
+import cytoscape.CyNode;
+import cytoscape.Cytoscape;
+import cytoscape.data.CyAttributes;
+import cytoscape.layout.Tunable;
+import cytoscape.layout.TunableListener;
+import cytoscape.logger.CyLogger;
+import cytoscape.task.TaskMonitor;
+
+import clusterMaker.ClusterMaker;
+import clusterMaker.algorithms.AbstractNetworkClusterer;
+import clusterMaker.algorithms.NodeCluster;
+import clusterMaker.algorithms.ClusterAlgorithm;
+import clusterMaker.algorithms.ClusterResults;
+import clusterMaker.algorithms.DistanceMatrix;
+import clusterMaker.algorithms.EdgeAttributeHandler;
+import clusterMaker.ui.ClusterViz;
+import clusterMaker.ui.NewNetworkView;
+
+// clusterMaker imports
+
+public class SCPSCluster extends AbstractNetworkClusterer {
+
+ double epsilon = 1.0;
+ int rNumber = 8;
+ int knumber = -1;
+
+ EdgeAttributeHandler edgeAttributeHandler = null;
+
+ TaskMonitor monitor = null;
+ CyLogger logger = null;
+ RunSCPS runSCPS = null;
+
+ public SCPSCluster() {
+ super();
+ clusterAttributeName =
Cytoscape.getCurrentNetwork().getIdentifier()+"_SCPS_cluster";
+ logger = CyLogger.getLogger(SCPSCluster.class);
+ initializeProperties();
+ }
+
+ public String getShortName() {return "SCPS";};
+ public String getName() {return "SCPS cluster";};
+
+ public JPanel getSettingsPanel() {
+ // Everytime we ask for the panel, we want to update our
attributes
+ edgeAttributeHandler.updateAttributeList();
+
+ return clusterProperties.getTunablePanel();
+ }
+
+ public ClusterViz getVisualizer() {
+ return new NewNetworkView(true);
+ }
+
+ public void initializeProperties() {
+ super.initializeProperties();
+
+ /**
+ * Tuning values
+ */
+ clusterProperties.add(new Tunable("tunables_panel",
+ "SCPS Tuning",
+ Tunable.GROUP, new
Integer(3)));
+
+ // Lambda Parameter
+ clusterProperties.add(new Tunable("epsilon",
+ "epsilon Parameter",
+ Tunable.DOUBLE, new
Double(1.0),
+ (Object)null, (Object)null,
0));
+
+
+ // Number of iterations
+ clusterProperties.add(new Tunable("rNumber",
+ "Number of iterations",
+ Tunable.INTEGER, new
Integer(10),
+ (Object)null, (Object)null,
0));
+
+
+ // K. Trumps epsilon if set to greater then zero
+ clusterProperties.add(new Tunable("knumber",
+ "Number of clusters",
+ Tunable.INTEGER, new
Integer(-1),
+ (Object)null, (Object)null,
0));
+
+
+ // Use the standard edge attribute handling stuff....
+ edgeAttributeHandler = new
EdgeAttributeHandler(clusterProperties, true);
+
+ super.advancedProperties();
+
+ clusterProperties.initializeProperties();
+ updateSettings(true);
+ }
+
+ public void updateSettings() {
+ updateSettings(false);
+ }
+
+ public void updateSettings(boolean force) {
+ clusterProperties.updateValues();
+ super.updateSettings(force);
+
+ Tunable t = clusterProperties.get("epsilon");
+ if ((t != null) && (t.valueChanged() || force))
+ epsilon = ((Double) t.getValue()).doubleValue();
+
+ t = clusterProperties.get("knumber");
+ if ((t != null) && (t.valueChanged() || force))
+ knumber = ((Integer) t.getValue()).intValue();
+
+
+ t = clusterProperties.get("rNumber");
+ if ((t != null) && (t.valueChanged() || force))
+ rNumber = ((Integer) t.getValue()).intValue();
+
+ edgeAttributeHandler.updateSettings(force);
+ }
+
+ public void doCluster(TaskMonitor monitor) {
+ this.monitor = monitor;
+ CyNetwork network = Cytoscape.getCurrentNetwork();
+ String networkID = network.getIdentifier();
+
+ CyAttributes netAttributes = Cytoscape.getNetworkAttributes();
+ CyAttributes nodeAttributes = Cytoscape.getNodeAttributes();
+
+ DistanceMatrix matrix = edgeAttributeHandler.getMatrix();
+
+ //Cluster the nodes
+ runSCPS = new RunSCPS(matrix, epsilon, knumber, rNumber,
logger);
+
+ List<NodeCluster> clusters = runSCPS.run(monitor);
+
+ logger.info("Removing groups");
+
+ // Remove any leftover groups from previous runs
+ removeGroups(netAttributes, networkID);
+
+ logger.info("Creating groups");
+ monitor.setStatus("Creating groups");
+
+ List<List<CyNode>> nodeClusters =
+ createGroups(netAttributes, networkID, nodeAttributes,
clusters);
+
+ ClusterResults results = new ClusterResults(network,
nodeClusters);
+ monitor.setStatus("Done. SCPS results:\n"+results);
+
+
+ // Tell any listeners that we're done
+ pcs.firePropertyChange(ClusterAlgorithm.CLUSTER_COMPUTED, null,
this);
+ }
+
+ public void halt() {
+ runSCPS.halt();
+ }
+
+}
--
You received this message because you are subscribed to the Google Groups
"cytoscape-cvs" group.
To post to this group, send email to [email protected].
To unsubscribe from this group, send email to
[email protected].
For more options, visit this group at
http://groups.google.com/group/cytoscape-cvs?hl=en.
