package medusa.georgios.enhanced_mcl;
import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Iterator;
/*
* Copyright (c) 2006 Gregor Heinrich. 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.
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESSED 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 AUTHOR 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.
*/
/*
* Created on 03.08.2006
*/
/**
* My_MarkovClustering implements the Markov clustering (MCL) algorithm for graphs,
* using a HashMap-based sparse representation of a Markov matrix, i.e., an
* adjacency matrix m that is normalised to one. Elements in a column / node can
* be interpreted as decision probabilities of a random walker being at that
* node. Note: whereas we explain the algorithms with columns, the actual
* implementation works with rows because the used sparse matrix has row-major
* format.
*
* The basic idea underlying the MCL algorithm is that dense regions in sparse
* graphs correspond with regions in which the number of k-length paths is
* relatively large, for small k in N, which corresponds to multiplying
* probabilities in the matrix appropriately. Random walks of length k have
* higher probability (product) for paths with beginning and ending in the same
* dense region than for other paths.
*
* The algorithm starts by creating a Markov matrix from the graph, for which
* first the adjacency matrix is added diagonal elements to include self-loops
* for all nodes, i.e., probabilities that the random walker stays at a
* particular node. After this initialisation, the algorithm works by
* alternating two operations, expansion and inflation, iteratively recomputing
* the set of transition probabilities. The expansion step corresponds to matrix
* multiplication (on stochastic matrices), the inflation step corresponds with
* a parametrized inflation operator Gamma_r, which acts column-wise on (column)
* stochastic matrices (here, we use row-wise operation, which is analogous).
*
* The inflation operator transforms a stochastic matrix into another one by
* raising each element to a positive power p and re-normalising columns to keep
* the matrix stochastic. The effect is that larger probabilities in each column
* are emphasised and smaller ones deemphasised. On the other side, the matrix
* multiplication in the expansion step creates new non-zero elements, i.e.,
* edges. The algorithm converges very fast, and the result is an idempotent
* Markov matrix, M = M * M, which represents a hard clustering of the graph
* into components.
*
* Expansion and inflation have two opposing effects: While expansion flattens
* the stochastic distributions in the columns and thus causes paths of a random
* walker to become more evenly spread, inflation contracts them to favoured
* paths.
*
* Description is based on the introduction of Stijn van Dongen's thesis Graph
* Clustering by Flow Simulation (2000); for a mathematical treatment of the
* algorithm and the associated MCL process, see there.
*
* @author gregor :: arbylon . net
*/
public class My_MarkovClustering {
public static SparseMatrix aa;
public static double[][] a;
public static List2Graph g;
int clusters = 0;
String report;
/**
* test the MCL algorithm with the matrix loaded from the file in the
* argument. If no argument is given, the default file for the matrix T(G3 +
* I) from van Dongen (2000), page 50.
*
* @param args
*/
public static void main(String[] args) {
String filename;
int elements;
System.out.println("\nrun like java -jar far_filename input_file inflation\n\n");
// filename=args[0];
filename = "example.txt";
// filename="gavin_best_MCL.txt";
// filename="tong2004_MCL.txt";
// filename="MIPS_MCL.txt";
// filename="DIP_MCL.txt";
// filename="BIOGRID_MCL.txt";
//filename="2.txt";
g = new List2Graph(filename);
a = g.return_double_matrix();
aa = new SparseMatrix(a);
// we use the transpose because our sparse matrices are row-major
aa = aa.transpose();
double maxResidual = 0.001;
// double gammaExp = 2;
double gammaExp = Double.parseDouble(args[1]);
double loopGain = 0.05;//einai i main
double zeroMax = 0.001;
System.out.println("filename:" + filename);
System.out.println("gamma:" + gammaExp);
My_MarkovClustering mcl = new My_MarkovClustering();
aa = mcl.run(aa, maxResidual, gammaExp, loopGain, zeroMax);
mcl.show_clusters();
mcl.Save_To_File("MCL_results.txt", true);
}
/**
* run the MCL process.
*
* @param a matrix
* @param maxResidual maximum difference between row elements and row square
* sum (measure of idempotence)
* @param pGamma inflation exponent for Gamma operator
* @param loopGain values for cycles
* @param maxZero maximum value considered zero for pruning operations
* @return the resulting matrix
*/
public void show_clusters() {
StringBuffer str = new StringBuffer();
ArrayList used = new ArrayList();
for (int i = 0; i < a.length; i++) {
if(!used.contains(i))
if (aa.getColum(i).getLength() != 0) {
{
clusters++;
str.append("----------------- Cluster: " + clusters + " -----------------\n");
}
Iterator it = aa.getColum(i).keySet().iterator();
while (it.hasNext()) {
int index = (Integer) it.next();
{
used.add(index);
str.append(g.labels.get(index) + "\n");
}
}
}
}
System.out.println(str.toString());
}
public int getNumberOfClusters() {
return this.clusters;
}
public String Save_To_File(String filename, boolean storeFile) {
try {
clusters=0;
StringBuffer str = new StringBuffer();
ArrayList used = new ArrayList();
for (int i = 0; i < a.length; i++) {
if(!used.contains(i))
if (aa.getColum(i).getLength() != 0) {
{
clusters++;
str.append("----------------- Cluster: " + clusters + " -----------------\n");
}
Iterator it = aa.getColum(i).keySet().iterator();
while (it.hasNext()) {
int index = (Integer) it.next();
{
used.add(index);
str.append(g.labels.get(index) + "\n");
}
}
}
}
// System.out.println(str.toString());
if(storeFile==true){
BufferedWriter outp = new BufferedWriter(new FileWriter(filename));
outp.append(str.toString());
outp.close();}
this.report = str.toString();
return str.toString();
} catch (IOException e) {
}
return null;
}
public String getReport() {
return this.report;
}
public void setReport(String str) {
this.report = str;
}
public SparseMatrix run(SparseMatrix a, double maxResidual, double pGamma,
double loopGain, double maxZero) {
System.out.println("Inflation: " + pGamma);
// Vectors.setFormat(8, 3);
// print(a, "original matrix");
// add cycles
addLoops(a, loopGain);
if (loopGain > 0) {
print(a, "... with added loops");
}
// make stochastic
a.normaliseRows();
// print(a, "normalised");
double residual = 1.;
int i = 0;
// main iteration
while (residual > maxResidual) {
i++;
a = expand(a);
//print(a, "iteration " + i + ", after expansion");
residual = inflate(a, pGamma, maxZero);
// print(a, "iteration " + i + ", after inflation");
// System.out.println("residual energy = " + residual);
}
return a;
}
private static void print(SparseMatrix a, String label) {
//System.out.println(label + ":");
// use transpose to compare with van Dongen's thesis
//System.out.println(a.transpose().toStringDense());
}
public My_MarkovClustering() {
// nothing yet
}
public My_MarkovClustering(String filename) {
g = new List2Graph(filename);
a = g.return_double_matrix();
aa = new SparseMatrix(a);
// we use the transpose because our sparse matrices are row-major
aa = aa.transpose();
double maxResidual = 0.001;
double gammaExp = 2.0;
double loopGain = 0.001;
double zeroMax = 0.001;
My_MarkovClustering mcl = new My_MarkovClustering();
aa = mcl.run(aa, maxResidual, gammaExp, loopGain, zeroMax);
}
public My_MarkovClustering(ArrayList list) {
g = new List2Graph(list);
a = g.return_double_matrix();
aa = new SparseMatrix(a);
// we use the transpose because our sparse matrices are row-major
aa = aa.transpose();
double maxResidual = 0.001;
double gammaExp = 2.0;
double loopGain = 0.05;//edw---------------
double zeroMax = 0.001;
My_MarkovClustering mcl = new My_MarkovClustering();
aa = mcl.run(aa, maxResidual, gammaExp, loopGain, zeroMax);
}
/**
* inflate stochastic matrix by Hadamard (elementwise) exponentiation,
* pruning and normalisation :
*
* result = Gamma ( m, p ) = normalise ( prune ( m .^ p ) ).
*
* By convention, normalisation is done along rows (SparseMatrix has
* row-major representation)
*
* @param m matrix (mutable)
* @param p exponent as a double
* @param zeromax below which elements are pruned from the sparse matrix
* @return residuum value, m is modified.
*/
public double inflate(SparseMatrix m, double p, double zeromax) {
double res = 0.;
// matlab: m = m .^ p
m.hadamardPower(p);
// matlab: m(find(m < threshold)) = 0
m.prune(zeromax);
// matlab [for cols]: dinv = diag(1./sum(m)); m = m * dinv; return
// sum(m)
SparseVector rowsums = m.normalise(1.);
// check if done: if the maximum element
for (int i : rowsums.keySet()) {
SparseVector row = m.get(i);
double max = row.max();
double sumsq = row.sum(2.);
res = Math.max(res, max - sumsq);
}
return res;
}
/**
* expand stochastic quadratic matrix by sqaring it with itself: result = m *
* m. Here normalisation is rowwise.
*
* @param matrix
* @return new matrix (pointer != argument)
*/
public SparseMatrix expand(SparseMatrix m) {
m = m.times(m);
return m;
}
/**
* add loops with specific energy, which corresponds to adding loopGain to
* the diagonal elements.
*
* @param a
* @param loopGain
*/
private void addLoops(SparseMatrix a, double loopGain) {
if (loopGain <= 0) {
return;
}
for (int i = 0; i < a.size(); i++) {
a.add(i, i, loopGain);
}
}
}
Ads help maintain this website.