From 9d92de5621f7ae3ec517d409e1012b411a94c655 Mon Sep 17 00:00:00 2001
From: Jean-Daniel Fekete <Jean-Daniel.Fekete@inria.fr>
Date: Sat, 18 Apr 2015 08:29:24 -0400
Subject: [PATCH] barycenter returns the configuration with lowest crossings

---
 Makefile                     |   3 +-
 reorder.v1.js                |  36 ++++++---
 reorder.v1.min.js            |   2 +-
 src/barycenter.js            |  36 ++++++---
 src/count_crossings.js       |  56 ++++++++++++++
 src/package.js               |   3 +
 test/barycenter-test.js      |   2 +
 test/count_crossings-test.js | 142 +++++++++++++++++++++++++++++++++++
 8 files changed, 255 insertions(+), 25 deletions(-)
 create mode 100644 src/count_crossings.js
 create mode 100644 test/count_crossings-test.js

diff --git a/Makefile b/Makefile
index 69a948d..555cb83 100644
--- a/Makefile
+++ b/Makefile
@@ -1,6 +1,5 @@
 NODE_PATH = ./node_modules
 JS_COMPILER = $(NODE_PATH)/uglify-js/bin/uglifyjs
-JS_TESTER = $(NODE_PATH)/vows/bin/vows --nocolor -v
 
 all: \
 	reorder.v1.js \
@@ -40,7 +39,7 @@ reorder.v1.js: \
 	src/ca.js
 
 test: all
-	@$(JS_TESTER)
+	@npm test
 
 %.min.js: %.js Makefile
 	@rm -f $@
diff --git a/reorder.v1.js b/reorder.v1.js
index 0af19ef..c74be39 100644
--- a/reorder.v1.js
+++ b/reorder.v1.js
@@ -726,19 +726,20 @@ function median(neighbors) {
 	return (neighbors[lm]*rspan + neighbors[rm]*lspan) / (lspan+rspan);
 }
 
-reorder.barycenter1 = function(graph, comp, iter) {
+reorder.barycenter1 = function(graph, comp, max_iter) {
     var nodes = graph.nodes(),
-	layer1, layer2,
-	layer,
+	layer1, layer2, crossings, iter,
+	best_layer1, best_layer2, best_crossings, best_iter,
+	layer, 
 	i, v, neighbors;
 
     if (comp.length < 3)
 	return comp;
 
-    if (! iter)
-	iter = 24;
-    else if ((iter%2)==1)
-	iter++; // want even number of iterations
+    if (! max_iter)
+	max_iter = 24;
+    else if ((max_iter%2)==1)
+	max_iter++; // want even number of iterations
 
     layer1 = comp.filter(function(n) {
 	return graph.outEdges(n).length!=0;
@@ -750,9 +751,14 @@ reorder.barycenter1 = function(graph, comp, iter) {
     for (i = 0; i < layer2.length; i++)
 	nodes[layer2[i]].pos = i;
 
-    for (layer = layer1;
-	 iter--;
-	 layer = (layer == layer1) ? layer2 : layer1) {
+    best_crossings = count_crossings(graph, layer1, layer2);
+    best_layer1 = layer1;
+    best_layer2 = layer2;
+    best_iter = 0;
+
+    for (layer = layer1, iter = 0;
+	 iter < max_iter;
+	 iter++, layer = (layer == layer1) ? layer2 : layer1) {
 	for (i = 0; i < layer.length; i++) {
 	    // Compute the median/barycenter for this node and set
 	    // its (real) value into node.mval
@@ -782,8 +788,16 @@ reorder.barycenter1 = function(graph, comp, iter) {
 	});
 	for (i = 0; i < layer2.length; i++)
 	    nodes[layer[i]].pos = i;
+	crossings = count_crossings(graph, layer1, layer2);
+	if (crossings < best_crossings) {
+	    best_crossings = crossings;
+	    best_layer1 = layer1;
+	    best_layer2 = layer2;
+	    best_iter = iter;
+	}
     }
-    return [layer1, layer2];
+    console.log('Best iter: '+best_iter);
+    return [best_layer1, best_layer2, best_crossings];
 };
 reorder.dijkstra = function(graph) {
     var g = graph, dijkstra = {};
diff --git a/reorder.v1.min.js b/reorder.v1.min.js
index c3992ef..71c200b 100644
--- a/reorder.v1.min.js
+++ b/reorder.v1.min.js
@@ -1 +1 @@
-(function(exports){reorder={version:"0.0.2"};reorder.dot=science.lin.dot;reorder.length=science.lin.length;reorder.normalize=science.lin.normalize;reorder.printmat=function(m){var i,j,row,line;for(i=0;i<m.length;i++){row=m[i];line="";for(j=0;j<row.length;j++){if(line.length!=0)line+=", ";line+=row[j].toFixed(4)}console.log(i.toPrecision(3)+": "+line)}};reorder.assert=function(v,msg){if(!v){console.log(msg);throw msg||"Assertion failed"}};reorder.printhcluster=function(cluster,indent){if(cluster.left==null)return Array(indent+1).join(" ")+"id: "+cluster.id;return Array(indent+1).join(" ")+"id: "+cluster.id+", dist: "+cluster.dist+"\n"+reorder.printhcluster(cluster.left,indent+1)+"\n"+reorder.printhcluster(cluster.right,indent+1)};reorder.mean=science.stats.mean;reorder.meantranspose=function(v,j){var n=v.length;if(n==0)return NaN;var o=v[0].length,m=0,i=-1,row;while(++i<n)m+=(v[i][j]-m)/(i+1);return m};reorder.meancolumns=function(v){var n=v.length;if(n==0)return NaN;var o=v[0].length,m=v[0].slice(0),i=0,j,row;while(++i<n){row=v[i];for(j=0;j<o;j++)m[j]+=(row[j]-m[j])/(i+1)}return m};reorder.sum=function(v){var i=v.length,s=v[0];while(i-->1)s+=v[i];return s};function isNum(a,b){return!(isNaN(a)||isNaN(b)||a==Infinity||b==Infinity)}reorder.distance={euclidean:function(a,b){var i=a.length,s=0,x;while(i-->0){if(isNum(a[i],b[i])){x=a[i]-b[i];s+=x*x}}return Math.sqrt(s)},manhattan:function(a,b){var i=a.length,s=0;while(i-->0){if(isNum(a[i],b[i])){s+=Math.abs(a[i]-b[i])}}return s},minkowski:function(p){return function(a,b){var i=a.length,s=0;while(i-->0){if(isNum(a[i],b[i])){s+=Math.pow(Math.abs(a[i]-b[i]),p)}}return Math.pow(s,1/p)}},chebyshev:function(a,b){var i=a.length,max=0,x;while(i-->0){if(isNum(a[i],b[i])){x=Math.abs(a[i]-b[i]);if(x>max)max=x}}return max},hamming:function(a,b){var i=a.length,d=0;while(i-->0){if(isNum(a[i],b[i])){if(a[i]!==b[i])d++}}return d},jaccard:function(a,b){var n=0,i=a.length,s=0;while(i-->0){if(isNum(a[i],b[i])){if(a[i]===b[i])s++;n++}}if(n==0)return 0;return s/n},braycurtis:function(a,b){var i=a.length,s0=0,s1=0,ai,bi;while(i-->0){ai=a[i];bi=b[i];if(isNum(ai,bi)){s0+=Math.abs(ai-bi);s1+=Math.abs(ai+bi)}}if(s1==0)return 0;return s0/s1}};reorder.range=function(start,stop,step){if(arguments.length<3){step=1;if(arguments.length<2){stop=start;start=0}}var range=[],i=start;if(step<0)for(;i>stop;i+=step)range.push(i);else for(;i<stop;i+=step)range.push(i);return range};reorder.transpose=science.lin.transpose;reorder.transposeSlice=function(a,start,end){if(arguments.length<3){end=a[0].length;if(arguments.length<2){start=0}}var m=a.length,n=end,i=start-1,j,b=new Array(end-start);while(++i<n){b[i]=new Array(m);j=-1;while(++j<m)b[i-start][j]=a[j][i]}return b};reorder.correlation={pearson:function(a,b){var ma=science.stats.mean(a),mb=science.stats.mean(b),s1=0,s2=0,s3=0,i,dx,dy,n=Math.min(a.length,b.length);if(n===0)return NaN;for(i=0;i<n;i++){dx=a[i]-ma;dy=b[i]-mb;s1+=dx*dy;s2+=dx*dx;s3+=dy*dy}return s1/Math.sqrt(s2*s3)},pearsonMatrix:function(matrix){var a,ma,i,j,dx,n=matrix.length,ret,mx,sx,sx2;if(n===0)return NaN;mx=Array(n);sx=science.zeroes(n);sx2=science.zeroes(n);for(i=0;i<n;i++){mx[i]=science.stats.mean(matrix[i])}for(i=0;i<n;i++){a=matrix[i];ma=mx[i];for(j=0;j<n;j++){dx=a[j]-ma;sx[j]+=dx;sx2[j]+=dx*dx}}ret=Array(n);for(i=0;i<n;i++){ret[i]=Array(n);for(j=0;j<n;j++){ret[i][j]=sx[i]*sx[j]/Math.sqrt(sx2[i]*sx2[j])}}return ret}};reorder.heap=function(comp){var heap=[],index={},Heap={};Heap.length=function(){return heap.length};Heap.insert=function(o){var i=heap.length;heap.push(null);percolateUp(i,o)};Heap.isEmpty=function(){return heap.length==0};Heap.peek=function(){if(heap.length==0)throw{error:"Empty heap"};return heap[0]};Heap.pop=function(){if(heap.length==0)throw{error:"Empty heap"};var top=heap[0];if(top==null)return top;var boto=heap[heap.length-1];heap[0]=boto;index[boto]=0;heap.pop();if(heap.length>1)percolateDown(0);delete index[top];return top};function lChild(i){return i*2+1}function rChild(i){return i*2+2}function parent(i){return Math.floor((i-1)/2)}function swap(i,j){var io=heap[i],jo=heap[j];heap[i]=jo;index[jo]=i;heap[j]=io;index[io]=j}function update(o){var i=percolateUp(index[o],o);percolateDown(i)}function percolateDown(cur){while(true){var left=lChild(cur),right=rChild(cur);var smallest;if(left<heap.length&&comp(heap[left],heap[cur])<0)smallest=left;else smallest=cur;if(right<heap.length&&comp(heap[right],heap[smallest])<0)smallest=right;if(cur==smallest)return;swap(cur,smallest);cur=smallest}}function percolateUp(cur,o){var i=cur;for(par=parent(i);i>0&&comp(heap[par],o)>0;par=parent(par)){var p=heap[par];heap[i]=p;index[p]=i;i=par}heap[i]=o;index[o]=i;return i}return Heap};reorder.permutation=reorder.range;function inverse_permutation(perm){var inv={};for(var i=0;i<perm.length;i++){inv[perm[i]]=i}return inv}reorder.inverse_permutation=inverse_permutation;reorder.graph=function(nodes,links,directed){var graph={},linkDistance=1,edges,inEdges,outEdges,distances,components;graph.nodes=function(x){if(!arguments.length)return nodes;nodes=x;return graph};graph.links=function(x){if(!arguments.length)return links;links=x;return graph};graph.linkDistance=function(x){if(!arguments.length)return linkDistance;linkDistance=typeof x==="function"?x:+x;return graph};graph.directed=function(x){if(!arguments.length)return directed;directed=x;return graph};function init(){var i,o,n=nodes.length,m=links.length;for(i=0;i<n;++i){(o=nodes[i]).index=i;o.weight=0}for(i=0;i<m;++i){(o=links[i]).index=i;if(typeof o.source=="number")o.source=nodes[o.source];if(typeof o.target=="number")o.target=nodes[o.target];++o.source.weight;++o.target.weight}distances=[];if(typeof linkDistance==="function")for(i=0;i<m;++i)distances[i]=+linkDistance.call(this,links[i],i);else for(i=0;i<m;++i)distances[i]=linkDistance;edges=Array(nodes.length);for(i=0;i<nodes.length;++i){edges[i]=[]}if(directed){inEdges=Array(nodes.length);outEdges=Array(nodes.length);for(i=0;i<nodes.length;++i){inEdges[i]=[];outEdges[i]=[]}}for(i=0;i<links.length;++i){o=links[i];edges[o.source.index].push(o);if(directed)inEdges[o.source.index].push(o);if(o.source.index!=o.target.index)edges[o.target.index].push(o);if(directed)outEdges[o.target.index].push(o)}return graph}graph.init=init;graph.edges=function(node){return edges[node]};function inEdges(node){if(directed)return inEdges[node];return edges[node]}graph.inEdges=inEdges;function outEdges(node){if(directed)return outEdges[node];return edges[node]}graph.outEdges=outEdges;graph.sinks=function(){var sinks=[],i;if(!directed)return reorder.range(nodes.length);for(i=0;i<nodes.length;i++){if(outEdges(i).length==0)sinks.push(i)}return sinks};graph.sources=function(){var sources=[],i;if(!directed)return reorder.range(nodes.length);for(i=0;i<nodes.length;i++){if(inEdges(i).length==0)sources.push(i)}return sources};function distance(i){return distances[i]}graph.distance=distance;function neighbors(node){var e=edges[node],ret=[];for(var i=0;i<e.length;++i){var o=e[i];if(o.source.index==node)ret.push(o.target);else ret.push(o.source)}return ret}graph.neighbors=neighbors;graph.other=function(o,node){if(typeof o=="number")o=links[o];if(o.source.index==node)return o.target;else return o.source};function compute_components(){var stack=[],comp=0,comps=[],ccomp,n=nodes.length,i,j,v,l,o,e;for(i=0;i<n;i++)nodes[i].comp=0;for(j=0;j<n;j++){if(nodes[j].comp!=0)continue;comp=comp+1;nodes[j].comp=comp;stack.push(j);ccomp=[j];while(stack.length){v=stack.shift();l=edges[v];for(i=0;i<l.length;i++){e=l[i];o=e.source;if(o.index==v)o=e.target;if(o.index==v)continue;if(o.comp==0){o.comp=comp;ccomp.push(o.index);stack.push(o.index)}}}if(ccomp.length){ccomp.sort(function(a,b){return a-b});comps.push(ccomp)}}return comps}graph.components=function(){if(!components)components=compute_components();return components};return graph};reorder.printmat=function(mat,rowperm,colperm){var i,j,row,col,str;console.log("Matrix:");for(i=0;i<mat.length;i++){row=rowperm?mat[rowperm[i]]:mat[i];str="";for(j=0;j<row.length;j++){col=colperm?row[colperm[j]]:row[j];str+=col?"*":" "}console.log(str)}};reorder.mat2graph=function(mat,directed){var n=mat.length,nodes=[],links=[],i,j,v,m;for(i=0;i<n;i++)nodes.push({id:i});for(i=0;i<n;i++){v=mat[i];m=directed?v.length:i+1;m=Math.min(m,v.length);for(j=0;j<m;j++){if(j==nodes.length)nodes.push({id:j});if(v[j]!=0){links.push({source:i,target:j,value:v[j]})}}}return reorder.graph(nodes,links,directed).init()};reorder.graph2mat=function(graph,directed){var nodes=graph.nodes(),links=graph.links(),n=nodes.length,i,l,mat;if(!directed)directed=graph.directed();if(directed){var rows=n,cols=n;for(i=n-1;i>=0;i--){if(graph.inEdges(i).length!=0)break;else rows--}for(i=n-1;i>=0;i--){if(graph.outEdges(i).length!=0)break;else cols--}mat=science.zeroes(rows,cols);for(i=0;i<links.length;i++){l=links[i];mat[l.source.index][l.target.index]=l.value?l.value:1}}else{mat=science.zeroes(n,n);for(i=0;i<links.length;i++){l=links[i];mat[l.source.index][l.target.index]=l.value?l.value:1;mat[l.target.index][l.source.index]=l.value?l.value:1}}return mat};function count_crossings(graph,north,south){var i,j,n,firstIndex,treeSize,tree,index,weightSum,invert=false,crosscount;if(north.length<south.length){var tmp=north;north=south;south=tmp;invert=true}var south_inv=inverse_permutation(south),southsequence=[];for(i=0;i<north.length;i++){if(invert){n=graph.inEdges(north[i]).map(function(e){return south_inv[e.target.index]})}else{n=graph.outEdges(north[i]).map(function(e){return south_inv[e.source.index]})}n.sort(function(a,b){return a-b});southsequence=southsequence.concat(n)}firstIndex=1;while(firstIndex<south.length)firstIndex<<=1;treeSize=2*firstIndex-1;firstIndex-=1;tree=science.zeroes(treeSize);crosscount=0;for(i=0;i<southsequence.length;i++){index=southsequence[i]+firstIndex;tree[index]++;while(index>0){if(index%2)crosscount+=tree[index+1];index=index-1>>1;tree[index]++}}return crosscount}reorder.count_crossings=count_crossings;reorder.barycenter=function(graph,iter,comps){var perm=[];if(!comps){comps=graph.components()}for(var i=0;i<comps.length;i++)perm=perm.concat(reorder.barycenter1(graph,comps[i],iter));return perm};function median(neighbors){if(neighbors.length==0)return-1;if(neighbors.length==1)return neighbors[0];if(neighbors.length==2)return(neighbors[0]+neighbors[1])/2;neighbors.sort(function(a,b){return a-b});if(neighbors.length%2)return neighbors[(neighbors.length-1)/2];var rm=neighbors.length/2,lm=rm-1,rspan=neighbors[neighbors.length-1]-neighbors[rm],lspan=neighbors[lm]-neighbors[0];if(lspan==rspan)return(neighbors[lm]+neighbors[rm])/2;else return(neighbors[lm]*rspan+neighbors[rm]*lspan)/(lspan+rspan)}reorder.barycenter1=function(graph,comp,iter){var nodes=graph.nodes(),layer1,layer2,layer,i,v,neighbors;if(comp.length<3)return comp;if(!iter)iter=24;else if(iter%2==1)iter++;layer1=comp.filter(function(n){return graph.outEdges(n).length!=0});layer2=comp.filter(function(n){return graph.inEdges(n).length!=0});for(i=0;i<layer2.length;i++)nodes[layer2[i]].pos=i;for(layer=layer1;iter--;layer=layer==layer1?layer2:layer1){for(i=0;i<layer.length;i++){v=nodes[layer[i]];if(layer==layer1)neighbors=graph.outEdges(v.index);else neighbors=graph.inEdges(v.index);neighbors=neighbors.map(function(e){var n=e.source==v?e.target:e.source;return nodes[n.index].pos});v.median=median(neighbors)}layer.sort(function(a,b){var d=nodes[a].median-nodes[b].median;if(d==0){d=graph.edges(b).length%2-graph.edges(a).length%2}if(d<0)return-1;else if(d>0)return 1;return 0});for(i=0;i<layer2.length;i++)nodes[layer[i]].pos=i}return[layer1,layer2]};reorder.dijkstra=function(graph){var g=graph,dijkstra={};function allShortestPaths(from,queued){var preds={},edges,v,v2,e,d,D,queue=reorder.heap(function(i,j){return preds[i].weight-preds[j].weight}),p={edge:-1,vertex:from,weight:0},p2;if(!queued)queued={};queued[from]=p;queue.insert(from);while(!queue.isEmpty()){p=queued[queue.pop()];v=p.vertex;edges=graph.edges(v);for(var i=0;i<edges.length;i++){e=edges[i].index;d=p.weight+graph.distance(e);v2=graph.other(e,v).index;D=queued[v2];if(!D){p2={edge:e,vertex:v2,weight:d};queue.insert(v2);queued[v2]=p2}else if(D.weight>d){D.weight=d;D.edge=e;queue.update(D.vertex)}}}return queued}dijkstra.shortestPath=function(from,to){var map=allShortestPaths(from),path,v;v=map[to];path=[v];while(v.edge!=-1){v=map[graph.other(v.edge,v.vertex).index];path.unshift(v)}return path};return dijkstra};reorder.distmax=function(distMatrix){var max=0,n=distMatrix.length,i,j,row;for(i=0;i<n;i++){row=distMatrix[i];for(j=i+1;j<n;j++)if(row[j]>max)max=row[j]}return max};reorder.distmin=function(distMatrix){var min=Infinity,n=distMatrix.length,i,j,row;for(i=0;i<n;i++){row=distMatrix[i];for(j=i+1;j<n;j++)if(row[j]<min)min=row[j]}return min};reorder.dist=function(){var distance=reorder.distance.euclidean;function dist(vectors){var n=vectors.length,distMatrix=[];for(var i=0;i<n;i++){var d=[];distMatrix[i]=d;for(var j=0;j<n;j++){if(j<i){d.push(distMatrix[j][i])}else if(i===j){d.push(0)}else{d.push(distance(vectors[i],vectors[j]))}}}return distMatrix}dist.distance=function(x){if(!arguments.length)return distance;distance=x;return dist};return dist};reorder.dist_remove=function(dist,n,m){if(arguments.length<3)m=n+1;var i;dist.splice(n,m-n);for(i=dist.length;i-->0;)dist[i].splice(n,m-n);return dist};reorder.randomPermute=function(array,i,j){if(arguments.length<3){j=array.length;if(arguments.length<2){i=0}}var m=j-i,t,k;while(m>0){k=i+Math.floor(Math.random()*m--);t=array[i+m];array[i+m]=array[k];array[k]=t}return array};reorder.randomPermutation=function(n){return reorder.randomPermute(reorder.permutation(n))};reorder.randomMatrix=function(p,n,m,sym){if(!m)m=n;if(n!=m)sym=false;else if(!sym)sym=true;var mat=science.zeroes(n,m),i,j,cnt;if(sym){for(i=0;i<n;i++){cnt=0;for(j=0;j<i+1;j++){if(Math.random()<p){mat[i][j]=mat[j][i]=1;cnt++}}if(cnt==0){j=Math.floor(Math.random()*n/2);mat[i][j]=mat[j][i]=1}}}else{for(i=0;i<n;i++){cnt=0;for(j=0;j<m;j++){if(Math.random()<p){mat[i][j]=1;cnt++}}if(cnt==0)mat[i][Math.floor(Math.random()*m)]=1}}return mat};reorder.permute=function(list,indexes){var m=indexes.length;var copy=list.slice(0);while(m--)copy[m]=list[indexes[m]];return copy};reorder.permutetranspose=function(array,indexes){var m=array.length;while(m-->0)array[m]=reorder.permute(array[m],indexes);return array};reorder.stablepermute=function(list,indexes){var p=reorder.permute(list,indexes);if(p[0]>p[p.length-1])p.reverse();return p};if(typeof science=="undefined"){science={version:"1.9.1"};science.stats={}}science.stats.hcluster=function(){var distance=reorder.distance.euclidean,linkage="single",distMatrix=null;function hcluster(vectors){var n=vectors.length,dMin=[],cSize=[],clusters=[],c1,c2,c1Cluster,c2Cluster,p,root,i,j,id=0;if(distMatrix==null){distMatrix=[];i=-1;while(++i<n){dMin[i]=0;distMatrix[i]=[];j=-1;while(++j<n){distMatrix[i][j]=i===j?Infinity:distance(vectors[i],vectors[j]);if(distMatrix[i][dMin[i]]>distMatrix[i][j])dMin[i]=j}}}else{if(distMatrix.length<n||distMatrix[0].length<n)throw{error:"Provided distance matrix length "+distMatrix.length+" instead of "+n};i=-1;while(++i<n){dMin[i]=0;j=-1;while(++j<n){if(i===j)distMatrix[i][j]=Infinity;if(distMatrix[i][dMin[i]]>distMatrix[i][j])dMin[i]=j}}}i=-1;while(++i<n){clusters[i]=[];clusters[i][0]={left:null,right:null,dist:0,centroid:vectors[i],id:id++,size:1,depth:0};cSize[i]=1}for(p=0;p<n-1;p++){c1=0;for(i=0;i<n;i++){if(distMatrix[i][dMin[i]]<distMatrix[c1][dMin[c1]])c1=i}c2=dMin[c1];c1Cluster=clusters[c1][0];c2Cluster=clusters[c2][0];var newCluster={left:c1Cluster,right:c2Cluster,dist:distMatrix[c1][c2],centroid:calculateCentroid(c1Cluster.size,c1Cluster.centroid,c2Cluster.size,c2Cluster.centroid),id:id++,size:c1Cluster.size+c2Cluster.size,depth:1+Math.max(c1Cluster.depth,c2Cluster.depth)};clusters[c1].splice(0,0,newCluster);cSize[c1]+=cSize[c2];for(j=0;j<n;j++){switch(linkage){case"single":if(distMatrix[c1][j]>distMatrix[c2][j])distMatrix[j][c1]=distMatrix[c1][j]=distMatrix[c2][j];break;case"complete":if(distMatrix[c1][j]<distMatrix[c2][j])distMatrix[j][c1]=distMatrix[c1][j]=distMatrix[c2][j];break;case"average":distMatrix[j][c1]=distMatrix[c1][j]=(cSize[c1]*distMatrix[c1][j]+cSize[c2]*distMatrix[c2][j])/(cSize[c1]+cSize[j]);break}}distMatrix[c1][c1]=Infinity;for(i=0;i<n;i++)distMatrix[i][c2]=distMatrix[c2][i]=Infinity;for(j=0;j<n;j++){if(dMin[j]==c2)dMin[j]=c1;if(distMatrix[c1][j]<distMatrix[c1][dMin[c1]])dMin[c1]=j}root=newCluster}return root}hcluster.linkage=function(x){if(!arguments.length)return linkage;linkage=x;return hcluster};hcluster.distance=function(x){if(!arguments.length)return distance;distance=x;return hcluster};hcluster.distanceMatrix=function(x){if(!arguments.length)return distMatrix;distMatrix=x.map(function(y){return y.slice(0)});return hcluster};return hcluster};function calculateCentroid(c1Size,c1Centroid,c2Size,c2Centroid){var newCentroid=[],newSize=c1Size+c2Size,n=c1Centroid.length,i=-1;while(++i<n){newCentroid[i]=(c1Size*c1Centroid[i]+c2Size*c2Centroid[i])/newSize}return newCentroid}reorder.leafOrder=function(){var distanceMatrix=null,distance=reorder.distance.euclidean,linkage="complete",debug=0,leavesMap={},orderMap={};function isLeaf(n){return n.depth==0}function leaves(n){if(n==null)return[];if(n.id in leavesMap)return leavesMap[n.id];return leavesMap[n.id]=_leaves(n)}function _leaves(n){if(n==null)return[];if(n.depth==0)return[n.id];return leaves(n.left).concat(leaves(n.right))}function order(v,i,j){var key="k"+v.id+"-"+i+"-"+j;if(key in orderMap)return orderMap[key];return orderMap[key]=_order(v,i,j)}function _order(v,i,j){if(v.depth==0)return[0,[v.id]];var l=v.left,r=v.right;var L=leaves(l),R=leaves(r);var w,x;if(L.indexOf(i)!=-1&&R.indexOf(j)!=-1){w=l;x=r}else if(R.indexOf(i)!=-1&&L.indexOf(j)!=-1){w=r;x=l}else throw{error:"Node is not common ancestor of "+i+", "+j};var Wl=leaves(w.left),Wr=leaves(w.right);var Ks=Wr.indexOf(i)!=-1?Wl:Wr;if(Ks.length==0)Ks=[i];var Xl=leaves(x.left),Xr=leaves(x.right);var Ls=Xr.indexOf(j)!=-1?Xl:Xr;if(Ls.length==0)Ls=[j];var max=Infinity,optimal_order=[];for(var k=0;k<Ks.length;k++){var w_max=order(w,i,Ks[k]);for(var m=0;m<Ls.length;m++){var x_max=order(x,Ls[m],j);var sim=w_max[0]+distanceMatrix[Ks[k]][Ls[m]]+x_max[0];if(sim<max){max=sim;optimal_order=w_max[1].concat(x_max[1])}}}return[max,optimal_order]}function orderFull(v){leavesMap={};orderMap={};var max=Infinity,optimal_order=[],left=leaves(v.left),right=leaves(v.right);if(debug)console.log(reorder.printhcluster(v,0));for(var i=0;i<left.length;i++){for(var j=0;j<right.length;j++){var so=order(v,left[i],right[j]);if(so[0]<max){max=so[0];optimal_order=so[1]}}}distanceMatrix=null;return optimal_order}function leafOrder(vector){if(distanceMatrix==null)distanceMatrix=reorder.dist().distance(distance)(vector);var hcluster=science.stats.hcluster().linkage(linkage).distanceMatrix(distanceMatrix);return orderFull(hcluster(vector))}leafOrder.debug=function(x){if(!arguments.length)return debug;debug=x;return leafOrder};leafOrder.distance=function(x){if(!arguments.length)return distance;distance=x;distanceMatrix=null;return leafOrder};leafOrder.linkage=function(x){if(!arguments.length)return linkage;linkage=x;return leafOrder};leafOrder.distanceMatrix=function(x){if(!arguments.length)return distanceMatrix;distanceMatrix=x.map(function(y){return y.slice(0)});return leafOrder};leafOrder.orderFull=orderFull;return leafOrder};reorder.order=function(){var distance=reorder.distance.euclidean,ordering=reorder.leafOrder,linkage="complete",distanceMatrix,vector,except=[],debug=0,i=0,j=Infinity;function _reset(){distance=reorder.distance.euclidean;ordering=reorder.leafOrder;linkage="complete";distanceMatrix=null;vector=null;except=[];debug=0;i=0;j=Infinity}function order(v){vector=v;j=Math.min(j,v.length);var i0=i>0?i-1:0,j0=j<vector.length?j+1:j,k,low,high;for(k=except.length-1;k>0;k-=2){low=except[k-1];high=except[k];if(high>=j0){if(j0>j){j0=Math.min(j0,low+1);except.splice(k-1,2)}else{high=j0}}else if(low<=i0){if(i0<i){i0=Math.max(i0,high-1);except.splice(k-1,2)}else{low=i0}}else if(high-low<3)except.splice(k-1,2)}try{return _order_limits(i0,j0)}finally{_reset()}}function _order_limits(i0,j0){var orig=vector,perm,row,k,l;vector=vector.slice(i0,j0);if(i==0&&j==vector.length)return _order_except();if(debug)console.log("i0="+i0+" j0="+j0);if(distanceMatrix!=null){if(j0!=vector.length)reorder.dist_remove(distanceMatrix,j0,vector.length);if(i0>0)reorder.dist_remove(distanceMatrix,0,i0)}else{_compute_dist()}var max=reorder.distmax(distanceMatrix);if(i0<i){row=distanceMatrix[0];for(k=row.length;k-->1;)row[k]+=max;for(k=distanceMatrix.length;k-->1;)distanceMatrix[k][0]+=max;max+=max;if(i0!=0){for(k=0;k<except.length;k++)except[k]-=i0}}if(j0>j){l=distanceMatrix.length-1;row=distanceMatrix[l];for(k=l;k-->0;){row[k]+=max;distanceMatrix[k][l]+=max}}perm=_order_except();if(i0<i){if(perm[0]!=0)perm.reverse();if(j0>j){reorder.assert(perm[0]==0&&perm[perm.length-1]==perm.length-1,"Invalid constrained permutation endpoints")}else{reorder.assert(perm[0]==0,"Invalid constrained permutation start")}}else if(j0>j){if(perm[perm.length-1]!=perm.length-1)perm=perm.reverse();reorder.assert(perm[perm.length-1]==perm.length-1,"Invalid constrained permutation end")}if(i0!=0){perm=reorder.permutation(i0).concat(perm.map(function(v){return v+i0}))}if(orig.length>j0){perm=perm.concat(reorder.range(j0,orig.length))}return perm}function _order_except(){var perm,k,l,low,high,pos;if(except.length==0)return _order_equiv();_compute_dist();for(k=except.length-1;k>0;k-=2){low=except[k-1];high=except[k];distanceMatrix=reorder.dist_remove(distanceMatrix,low+1,high-1);vector.splice(low+1,high-low-2);if(debug)console.log("Except["+low+", "+high+"]");if(distanceMatrix[low][low+1]!=0){distanceMatrix[low][low+1]=distanceMatrix[low+1][low]=-1}}perm=_order_equiv();for(k=0;k<except.length;k+=2){low=except[k];high=except[k+1];for(l=0;l<perm.length;l++){if(perm[l]>low)perm[l]+=high-low-2;else if(perm[l]==low)pos=l}if(pos>0&&perm[pos-1]==high-1){Array.prototype.splice.apply(perm,[pos,0].concat(reorder.range(high-2,low,-1)))}else if(perm[pos+1]==high-1){Array.prototype.splice.apply(perm,[pos+1,0].concat(reorder.range(low+1,high-1)))}else{throw"Range not respected"}}return perm}function _order_equiv(){var perm,row,e,j,k,l,m,n,has_1=false,equiv=[],fix_except={};_compute_dist();for(k=0;k<distanceMatrix.length-1;k++){row=distanceMatrix[k];e=[];j=row.indexOf(-1);if(j!=-1){fix_except[k]=[k,j];has_1=true}for(l=row.length;--l>k;){if(row[l]==0){j=distanceMatrix[l].indexOf(-1);if(j!=-1){fix_except[k]=[l,j];distanceMatrix[j][k]=row[j]=-1;has_1=true}e.unshift(l);distanceMatrix=reorder.dist_remove(distanceMatrix,l);vector.splice(l,1)}else if(row[l]<0)has_1=true}if(e.length!=0){e.unshift(k);equiv.push(e)}}if(has_1){for(k=0;k<distanceMatrix.length-1;k++){row=distanceMatrix[k];for(l=k+1;l<row.length-1;l++){if(distanceMatrix[l][l+1]==-1){distanceMatrix[l+1][l]=distanceMatrix[l][l+1]=0}}}}perm=_order();for(k=equiv.length;k-->0;){e=equiv[k];l=perm.indexOf(e[0]);m=fix_except[e[0]];if(m&&m[0]==e[0]){l=_fix_exception(perm,l,m[0],m[1],0);m=undefined}for(n=1;n<e.length;n++){perm=_perm_insert(perm,l,e[n]);if(m&&m[0]==e[n]){l=_fix_exception(perm,l,m[0],m[1],n);m=undefined}}}return perm}function _fix_exception(perm,l,m,next,len){var i,j,k;if(l>0&&perm[l-1]==next){_swap(perm,l,perm.indexOf(m));return l+1}else if(perm[l+len+1]==next){_swap(perm,l+len,perm.indexOf(m));return l}else throw"Index not found"}function _swap(perm,a,b){if(a==b)return;var c=perm[a];perm[a]=perm[b];perm[b]=c}function _order(){if(debug>1)reorder.printmat(distanceMatrix);if(debug>2)reorder.printmat(vector);var perm=ordering().debug(debug).linkage(linkage).distanceMatrix(distanceMatrix)(vector);if(debug)console.log("Permutation: "+perm);return perm}function _perm_insert(perm,i,nv){perm=perm.map(function(v){return v<nv?v:v+1});perm.splice(i,0,nv);return perm}order.debug=function(x){if(!arguments.length)return debug;debug=x;return order};function _compute_dist(){if(distanceMatrix==null)distanceMatrix=reorder.dist().distance(distance)(vector);return distanceMatrix}order.distance=function(x){if(!arguments.length)return distance;distance=x;return order};order.linkage=function(x){if(!arguments.length)return linkage;linkage=x;return order};order.limits=function(x,y){if(!arguments.length)return[i,j];i=x;j=y;return order};order.except=function(list){if(!arguments.length)return except.slice(0);except=list.sort(function(a,b){if(a>=b)throw"Invalid list, indices not sorted";return a-b});return order};function _orderExcept(vector,i,j){var distanceMatrix=reorder.dist().distance(distance)(vector);var row,k,l,rev=false,args,pos=-1;distanceMatrix[i][i+1]=0;distanceMatrix[i+1][i]=0;var perm=ordering().distanceMatrix(distanceMatrix)(vector);pos=perm.indexOf(i);for(k=0;k<perm.length;k++){l=perm[k];if(l>i)perm[k]+=j-i-2}if(pos!=0&&perm[pos-1]==j-1)rev=true;if(rev){perm.reverse();pos=perm.length-pos-1}args=[pos+1,0].concat(reorder.range(i+1,j-1));Array.prototype.splice.apply(perm,args);return perm}order.orderrowsexcept=order.orderexcept;return order};reorder.covariance=science.lin.dot;reorder.covariancetranspose=function(v,a,b){var n=v.length,cov=0,i;for(i=0;i<n;i++){cov+=v[i][a]*v[i][b]}return cov};reorder.variancecovariance=function(v){var o=v[0].length,cov=Array(o),i,j;for(i=0;i<o;i++){cov[i]=Array(o)}for(i=0;i<o;i++){for(j=i;j<o;j++)cov[i][j]=cov[j][i]=reorder.covariancetranspose(v,i,j)}return cov};function normalize(v){var norm=science.lin.length(v),i=v.length;while(i-->0)v[i]/=norm;return v}reorder.poweriteration=function(v,eps){if(arguments.length<2)eps=1e-4;var n=v.length,b=Array(n),i,j,tmp=Array(n),norm,s=10;reorder.assert(n==v[0].length,"poweriteration needs a square matrix");for(i=0;i<n;i++)b[i]=Math.random();b=normalize(b);while(s-->0){for(i=0;i<n;i++){tmp[i]=0;for(j=0;j<n;j++)tmp[i]+=v[i][j]*b[j]}tmp=normalize(tmp);var t=tmp;tmp=b;b=t}return tmp};reorder.sortorder=function(v){return reorder.range(0,v.length).sort(function(a,b){return v[a]-v[b]})};reorder.center=function(v){var n=v.length;if(n==0)return null;var mean=reorder.meancolumns(v),o=mean.length,v1=Array(n),i,j,row;for(i=0;i<n;i++){row=v[i].slice(0);for(j=0;j<o;j++){row[j]-=mean[j]}v1[i]=row}return v1};reorder.pca1d=function(v,eps){if(arguments.length<2)eps=1e-4;var n=v.length;if(v.length==0)return null;v=reorder.center(v);var cov=reorder.variancecovariance(v);return reorder.poweriteration(cov,eps)};reorder.pca1dorder=function(v,eps){return reorder.sortorder(pca1d(v,eps))};reorder.sumlines=function(v){var n=v.length,o=v[0].length,sumline=Array(n),i,j,row,s;for(i=0;i<n;i++){row=v[i];s=0;for(j=0;j<o;j++)s+=row[j];sumline[i]=s}return sumline};reorder.sumcols=function(v){var n=v.length,o=v[0].length,sumcol=science.zeroes(o),i,j,row;for(i=0;i<n;i++){row=v[i];for(j=0;j<o;j++)sumcol[j]+=row[j]}return sumcol};reorder.ca=function(v,eps){if(arguments.length<2)eps=1e-4;var n=v.length,o=v[0].length,sumline=reorder.sumlines(v),sumcol=reorder.sumcols(v),s=reorder.sum(sumcol),i,j,row;for(i=0;i<n;i++){v[i]=v[i].map(function(a){return a/s})}sumline=reorder.sumlines(v);sumcol=reorder.sumcols(v);var v2=Array(n),ep;for(i=0;i<n;i++){v2[i]=Array(o);for(j=0;j<o;j++){ep=sumline[i]*sumcol[j];v2[i][j]=(v[i][j]-ep)/Math.sqrt(ep)}}for(i=0;i<n;i++){for(j=0;j<o;j++){ep=sumline[i]*sumcol[j];v[i][j]=(v[i][j]-ep)/(sumcol[j]*Math.sqrt(sumline[i]))}}var cov=Array(n);for(i=0;i<n;i++)cov[i]=Array(n);for(i=0;i<n;i++)for(j=i;j<n;j++)cov[i][j]=cov[j][i]=reorder.covariance(v2[i],v2[j]);var eigenvector1=reorder.poweriteration(cov,eps),eigenvalue1=0;for(i=0;i<n;i++)eigenvalue1+=eigenvector1[i]*cov[i][0];eigenvalue1/=eigenvector1[0];return eigenvector1}})(this);
\ No newline at end of file
+(function(exports){reorder={version:"0.0.2"};reorder.dot=science.lin.dot;reorder.length=science.lin.length;reorder.normalize=science.lin.normalize;reorder.printmat=function(m){var i,j,row,line;for(i=0;i<m.length;i++){row=m[i];line="";for(j=0;j<row.length;j++){if(line.length!=0)line+=", ";line+=row[j].toFixed(4)}console.log(i.toPrecision(3)+": "+line)}};reorder.assert=function(v,msg){if(!v){console.log(msg);throw msg||"Assertion failed"}};reorder.printhcluster=function(cluster,indent){if(cluster.left==null)return Array(indent+1).join(" ")+"id: "+cluster.id;return Array(indent+1).join(" ")+"id: "+cluster.id+", dist: "+cluster.dist+"\n"+reorder.printhcluster(cluster.left,indent+1)+"\n"+reorder.printhcluster(cluster.right,indent+1)};reorder.mean=science.stats.mean;reorder.meantranspose=function(v,j){var n=v.length;if(n==0)return NaN;var o=v[0].length,m=0,i=-1,row;while(++i<n)m+=(v[i][j]-m)/(i+1);return m};reorder.meancolumns=function(v){var n=v.length;if(n==0)return NaN;var o=v[0].length,m=v[0].slice(0),i=0,j,row;while(++i<n){row=v[i];for(j=0;j<o;j++)m[j]+=(row[j]-m[j])/(i+1)}return m};reorder.sum=function(v){var i=v.length,s=v[0];while(i-->1)s+=v[i];return s};function isNum(a,b){return!(isNaN(a)||isNaN(b)||a==Infinity||b==Infinity)}reorder.distance={euclidean:function(a,b){var i=a.length,s=0,x;while(i-->0){if(isNum(a[i],b[i])){x=a[i]-b[i];s+=x*x}}return Math.sqrt(s)},manhattan:function(a,b){var i=a.length,s=0;while(i-->0){if(isNum(a[i],b[i])){s+=Math.abs(a[i]-b[i])}}return s},minkowski:function(p){return function(a,b){var i=a.length,s=0;while(i-->0){if(isNum(a[i],b[i])){s+=Math.pow(Math.abs(a[i]-b[i]),p)}}return Math.pow(s,1/p)}},chebyshev:function(a,b){var i=a.length,max=0,x;while(i-->0){if(isNum(a[i],b[i])){x=Math.abs(a[i]-b[i]);if(x>max)max=x}}return max},hamming:function(a,b){var i=a.length,d=0;while(i-->0){if(isNum(a[i],b[i])){if(a[i]!==b[i])d++}}return d},jaccard:function(a,b){var n=0,i=a.length,s=0;while(i-->0){if(isNum(a[i],b[i])){if(a[i]===b[i])s++;n++}}if(n==0)return 0;return s/n},braycurtis:function(a,b){var i=a.length,s0=0,s1=0,ai,bi;while(i-->0){ai=a[i];bi=b[i];if(isNum(ai,bi)){s0+=Math.abs(ai-bi);s1+=Math.abs(ai+bi)}}if(s1==0)return 0;return s0/s1}};reorder.range=function(start,stop,step){if(arguments.length<3){step=1;if(arguments.length<2){stop=start;start=0}}var range=[],i=start;if(step<0)for(;i>stop;i+=step)range.push(i);else for(;i<stop;i+=step)range.push(i);return range};reorder.transpose=science.lin.transpose;reorder.transposeSlice=function(a,start,end){if(arguments.length<3){end=a[0].length;if(arguments.length<2){start=0}}var m=a.length,n=end,i=start-1,j,b=new Array(end-start);while(++i<n){b[i]=new Array(m);j=-1;while(++j<m)b[i-start][j]=a[j][i]}return b};reorder.correlation={pearson:function(a,b){var ma=science.stats.mean(a),mb=science.stats.mean(b),s1=0,s2=0,s3=0,i,dx,dy,n=Math.min(a.length,b.length);if(n===0)return NaN;for(i=0;i<n;i++){dx=a[i]-ma;dy=b[i]-mb;s1+=dx*dy;s2+=dx*dx;s3+=dy*dy}return s1/Math.sqrt(s2*s3)},pearsonMatrix:function(matrix){var a,ma,i,j,dx,n=matrix.length,ret,mx,sx,sx2;if(n===0)return NaN;mx=Array(n);sx=science.zeroes(n);sx2=science.zeroes(n);for(i=0;i<n;i++){mx[i]=science.stats.mean(matrix[i])}for(i=0;i<n;i++){a=matrix[i];ma=mx[i];for(j=0;j<n;j++){dx=a[j]-ma;sx[j]+=dx;sx2[j]+=dx*dx}}ret=Array(n);for(i=0;i<n;i++){ret[i]=Array(n);for(j=0;j<n;j++){ret[i][j]=sx[i]*sx[j]/Math.sqrt(sx2[i]*sx2[j])}}return ret}};reorder.heap=function(comp){var heap=[],index={},Heap={};Heap.length=function(){return heap.length};Heap.insert=function(o){var i=heap.length;heap.push(null);percolateUp(i,o)};Heap.isEmpty=function(){return heap.length==0};Heap.peek=function(){if(heap.length==0)throw{error:"Empty heap"};return heap[0]};Heap.pop=function(){if(heap.length==0)throw{error:"Empty heap"};var top=heap[0];if(top==null)return top;var boto=heap[heap.length-1];heap[0]=boto;index[boto]=0;heap.pop();if(heap.length>1)percolateDown(0);delete index[top];return top};function lChild(i){return i*2+1}function rChild(i){return i*2+2}function parent(i){return Math.floor((i-1)/2)}function swap(i,j){var io=heap[i],jo=heap[j];heap[i]=jo;index[jo]=i;heap[j]=io;index[io]=j}function update(o){var i=percolateUp(index[o],o);percolateDown(i)}function percolateDown(cur){while(true){var left=lChild(cur),right=rChild(cur);var smallest;if(left<heap.length&&comp(heap[left],heap[cur])<0)smallest=left;else smallest=cur;if(right<heap.length&&comp(heap[right],heap[smallest])<0)smallest=right;if(cur==smallest)return;swap(cur,smallest);cur=smallest}}function percolateUp(cur,o){var i=cur;for(par=parent(i);i>0&&comp(heap[par],o)>0;par=parent(par)){var p=heap[par];heap[i]=p;index[p]=i;i=par}heap[i]=o;index[o]=i;return i}return Heap};reorder.permutation=reorder.range;function inverse_permutation(perm){var inv={};for(var i=0;i<perm.length;i++){inv[perm[i]]=i}return inv}reorder.inverse_permutation=inverse_permutation;reorder.graph=function(nodes,links,directed){var graph={},linkDistance=1,edges,inEdges,outEdges,distances,components;graph.nodes=function(x){if(!arguments.length)return nodes;nodes=x;return graph};graph.links=function(x){if(!arguments.length)return links;links=x;return graph};graph.linkDistance=function(x){if(!arguments.length)return linkDistance;linkDistance=typeof x==="function"?x:+x;return graph};graph.directed=function(x){if(!arguments.length)return directed;directed=x;return graph};function init(){var i,o,n=nodes.length,m=links.length;for(i=0;i<n;++i){(o=nodes[i]).index=i;o.weight=0}for(i=0;i<m;++i){(o=links[i]).index=i;if(typeof o.source=="number")o.source=nodes[o.source];if(typeof o.target=="number")o.target=nodes[o.target];++o.source.weight;++o.target.weight}distances=[];if(typeof linkDistance==="function")for(i=0;i<m;++i)distances[i]=+linkDistance.call(this,links[i],i);else for(i=0;i<m;++i)distances[i]=linkDistance;edges=Array(nodes.length);for(i=0;i<nodes.length;++i){edges[i]=[]}if(directed){inEdges=Array(nodes.length);outEdges=Array(nodes.length);for(i=0;i<nodes.length;++i){inEdges[i]=[];outEdges[i]=[]}}for(i=0;i<links.length;++i){o=links[i];edges[o.source.index].push(o);if(directed)inEdges[o.source.index].push(o);if(o.source.index!=o.target.index)edges[o.target.index].push(o);if(directed)outEdges[o.target.index].push(o)}return graph}graph.init=init;graph.edges=function(node){return edges[node]};function inEdges(node){if(directed)return inEdges[node];return edges[node]}graph.inEdges=inEdges;function outEdges(node){if(directed)return outEdges[node];return edges[node]}graph.outEdges=outEdges;graph.sinks=function(){var sinks=[],i;if(!directed)return reorder.range(nodes.length);for(i=0;i<nodes.length;i++){if(outEdges(i).length==0)sinks.push(i)}return sinks};graph.sources=function(){var sources=[],i;if(!directed)return reorder.range(nodes.length);for(i=0;i<nodes.length;i++){if(inEdges(i).length==0)sources.push(i)}return sources};function distance(i){return distances[i]}graph.distance=distance;function neighbors(node){var e=edges[node],ret=[];for(var i=0;i<e.length;++i){var o=e[i];if(o.source.index==node)ret.push(o.target);else ret.push(o.source)}return ret}graph.neighbors=neighbors;graph.other=function(o,node){if(typeof o=="number")o=links[o];if(o.source.index==node)return o.target;else return o.source};function compute_components(){var stack=[],comp=0,comps=[],ccomp,n=nodes.length,i,j,v,l,o,e;for(i=0;i<n;i++)nodes[i].comp=0;for(j=0;j<n;j++){if(nodes[j].comp!=0)continue;comp=comp+1;nodes[j].comp=comp;stack.push(j);ccomp=[j];while(stack.length){v=stack.shift();l=edges[v];for(i=0;i<l.length;i++){e=l[i];o=e.source;if(o.index==v)o=e.target;if(o.index==v)continue;if(o.comp==0){o.comp=comp;ccomp.push(o.index);stack.push(o.index)}}}if(ccomp.length){ccomp.sort(function(a,b){return a-b});comps.push(ccomp)}}return comps}graph.components=function(){if(!components)components=compute_components();return components};return graph};reorder.printmat=function(mat,rowperm,colperm){var i,j,row,col,str;console.log("Matrix:");for(i=0;i<mat.length;i++){row=rowperm?mat[rowperm[i]]:mat[i];str="";for(j=0;j<row.length;j++){col=colperm?row[colperm[j]]:row[j];str+=col?"*":" "}console.log(str)}};reorder.mat2graph=function(mat,directed){var n=mat.length,nodes=[],links=[],i,j,v,m;for(i=0;i<n;i++)nodes.push({id:i});for(i=0;i<n;i++){v=mat[i];m=directed?v.length:i+1;m=Math.min(m,v.length);for(j=0;j<m;j++){if(j==nodes.length)nodes.push({id:j});if(v[j]!=0){links.push({source:i,target:j,value:v[j]})}}}return reorder.graph(nodes,links,directed).init()};reorder.graph2mat=function(graph,directed){var nodes=graph.nodes(),links=graph.links(),n=nodes.length,i,l,mat;if(!directed)directed=graph.directed();if(directed){var rows=n,cols=n;for(i=n-1;i>=0;i--){if(graph.inEdges(i).length!=0)break;else rows--}for(i=n-1;i>=0;i--){if(graph.outEdges(i).length!=0)break;else cols--}mat=science.zeroes(rows,cols);for(i=0;i<links.length;i++){l=links[i];mat[l.source.index][l.target.index]=l.value?l.value:1}}else{mat=science.zeroes(n,n);for(i=0;i<links.length;i++){l=links[i];mat[l.source.index][l.target.index]=l.value?l.value:1;mat[l.target.index][l.source.index]=l.value?l.value:1}}return mat};function count_crossings(graph,north,south){var i,j,n,firstIndex,treeSize,tree,index,weightSum,invert=false,crosscount;if(north.length<south.length){var tmp=north;north=south;south=tmp;invert=true}var south_inv=inverse_permutation(south),southsequence=[];for(i=0;i<north.length;i++){if(invert){n=graph.inEdges(north[i]).map(function(e){return south_inv[e.target.index]})}else{n=graph.outEdges(north[i]).map(function(e){return south_inv[e.source.index]})}n.sort(function(a,b){return a-b});southsequence=southsequence.concat(n)}firstIndex=1;while(firstIndex<south.length)firstIndex<<=1;treeSize=2*firstIndex-1;firstIndex-=1;tree=science.zeroes(treeSize);crosscount=0;for(i=0;i<southsequence.length;i++){index=southsequence[i]+firstIndex;tree[index]++;while(index>0){if(index%2)crosscount+=tree[index+1];index=index-1>>1;tree[index]++}}return crosscount}reorder.count_crossings=count_crossings;reorder.barycenter=function(graph,iter,comps){var perm=[];if(!comps){comps=graph.components()}for(var i=0;i<comps.length;i++)perm=perm.concat(reorder.barycenter1(graph,comps[i],iter));return perm};function median(neighbors){if(neighbors.length==0)return-1;if(neighbors.length==1)return neighbors[0];if(neighbors.length==2)return(neighbors[0]+neighbors[1])/2;neighbors.sort(function(a,b){return a-b});if(neighbors.length%2)return neighbors[(neighbors.length-1)/2];var rm=neighbors.length/2,lm=rm-1,rspan=neighbors[neighbors.length-1]-neighbors[rm],lspan=neighbors[lm]-neighbors[0];if(lspan==rspan)return(neighbors[lm]+neighbors[rm])/2;else return(neighbors[lm]*rspan+neighbors[rm]*lspan)/(lspan+rspan)}reorder.barycenter1=function(graph,comp,max_iter){var nodes=graph.nodes(),layer1,layer2,crossings,iter,best_layer1,best_layer2,best_crossings,best_iter,layer,i,v,neighbors;if(comp.length<3)return comp;if(!max_iter)max_iter=24;else if(max_iter%2==1)max_iter++;layer1=comp.filter(function(n){return graph.outEdges(n).length!=0});layer2=comp.filter(function(n){return graph.inEdges(n).length!=0});for(i=0;i<layer2.length;i++)nodes[layer2[i]].pos=i;best_crossings=count_crossings(graph,layer1,layer2);best_layer1=layer1;best_layer2=layer2;best_iter=0;for(layer=layer1,iter=0;iter<max_iter;iter++,layer=layer==layer1?layer2:layer1){for(i=0;i<layer.length;i++){v=nodes[layer[i]];if(layer==layer1)neighbors=graph.outEdges(v.index);else neighbors=graph.inEdges(v.index);neighbors=neighbors.map(function(e){var n=e.source==v?e.target:e.source;return nodes[n.index].pos});v.median=median(neighbors)}layer.sort(function(a,b){var d=nodes[a].median-nodes[b].median;if(d==0){d=graph.edges(b).length%2-graph.edges(a).length%2}if(d<0)return-1;else if(d>0)return 1;return 0});for(i=0;i<layer2.length;i++)nodes[layer[i]].pos=i;crossings=count_crossings(graph,layer1,layer2);if(crossings<best_crossings){best_crossings=crossings;best_layer1=layer1;best_layer2=layer2;best_iter=iter}}console.log("Best iter: "+best_iter);return[best_layer1,best_layer2,best_crossings]};reorder.dijkstra=function(graph){var g=graph,dijkstra={};function allShortestPaths(from,queued){var preds={},edges,v,v2,e,d,D,queue=reorder.heap(function(i,j){return preds[i].weight-preds[j].weight}),p={edge:-1,vertex:from,weight:0},p2;if(!queued)queued={};queued[from]=p;queue.insert(from);while(!queue.isEmpty()){p=queued[queue.pop()];v=p.vertex;edges=graph.edges(v);for(var i=0;i<edges.length;i++){e=edges[i].index;d=p.weight+graph.distance(e);v2=graph.other(e,v).index;D=queued[v2];if(!D){p2={edge:e,vertex:v2,weight:d};queue.insert(v2);queued[v2]=p2}else if(D.weight>d){D.weight=d;D.edge=e;queue.update(D.vertex)}}}return queued}dijkstra.shortestPath=function(from,to){var map=allShortestPaths(from),path,v;v=map[to];path=[v];while(v.edge!=-1){v=map[graph.other(v.edge,v.vertex).index];path.unshift(v)}return path};return dijkstra};reorder.distmax=function(distMatrix){var max=0,n=distMatrix.length,i,j,row;for(i=0;i<n;i++){row=distMatrix[i];for(j=i+1;j<n;j++)if(row[j]>max)max=row[j]}return max};reorder.distmin=function(distMatrix){var min=Infinity,n=distMatrix.length,i,j,row;for(i=0;i<n;i++){row=distMatrix[i];for(j=i+1;j<n;j++)if(row[j]<min)min=row[j]}return min};reorder.dist=function(){var distance=reorder.distance.euclidean;function dist(vectors){var n=vectors.length,distMatrix=[];for(var i=0;i<n;i++){var d=[];distMatrix[i]=d;for(var j=0;j<n;j++){if(j<i){d.push(distMatrix[j][i])}else if(i===j){d.push(0)}else{d.push(distance(vectors[i],vectors[j]))}}}return distMatrix}dist.distance=function(x){if(!arguments.length)return distance;distance=x;return dist};return dist};reorder.dist_remove=function(dist,n,m){if(arguments.length<3)m=n+1;var i;dist.splice(n,m-n);for(i=dist.length;i-->0;)dist[i].splice(n,m-n);return dist};reorder.randomPermute=function(array,i,j){if(arguments.length<3){j=array.length;if(arguments.length<2){i=0}}var m=j-i,t,k;while(m>0){k=i+Math.floor(Math.random()*m--);t=array[i+m];array[i+m]=array[k];array[k]=t}return array};reorder.randomPermutation=function(n){return reorder.randomPermute(reorder.permutation(n))};reorder.randomMatrix=function(p,n,m,sym){if(!m)m=n;if(n!=m)sym=false;else if(!sym)sym=true;var mat=science.zeroes(n,m),i,j,cnt;if(sym){for(i=0;i<n;i++){cnt=0;for(j=0;j<i+1;j++){if(Math.random()<p){mat[i][j]=mat[j][i]=1;cnt++}}if(cnt==0){j=Math.floor(Math.random()*n/2);mat[i][j]=mat[j][i]=1}}}else{for(i=0;i<n;i++){cnt=0;for(j=0;j<m;j++){if(Math.random()<p){mat[i][j]=1;cnt++}}if(cnt==0)mat[i][Math.floor(Math.random()*m)]=1}}return mat};reorder.permute=function(list,indexes){var m=indexes.length;var copy=list.slice(0);while(m--)copy[m]=list[indexes[m]];return copy};reorder.permutetranspose=function(array,indexes){var m=array.length;while(m-->0)array[m]=reorder.permute(array[m],indexes);return array};reorder.stablepermute=function(list,indexes){var p=reorder.permute(list,indexes);if(p[0]>p[p.length-1])p.reverse();return p};if(typeof science=="undefined"){science={version:"1.9.1"};science.stats={}}science.stats.hcluster=function(){var distance=reorder.distance.euclidean,linkage="single",distMatrix=null;function hcluster(vectors){var n=vectors.length,dMin=[],cSize=[],clusters=[],c1,c2,c1Cluster,c2Cluster,p,root,i,j,id=0;if(distMatrix==null){distMatrix=[];i=-1;while(++i<n){dMin[i]=0;distMatrix[i]=[];j=-1;while(++j<n){distMatrix[i][j]=i===j?Infinity:distance(vectors[i],vectors[j]);if(distMatrix[i][dMin[i]]>distMatrix[i][j])dMin[i]=j}}}else{if(distMatrix.length<n||distMatrix[0].length<n)throw{error:"Provided distance matrix length "+distMatrix.length+" instead of "+n};i=-1;while(++i<n){dMin[i]=0;j=-1;while(++j<n){if(i===j)distMatrix[i][j]=Infinity;if(distMatrix[i][dMin[i]]>distMatrix[i][j])dMin[i]=j}}}i=-1;while(++i<n){clusters[i]=[];clusters[i][0]={left:null,right:null,dist:0,centroid:vectors[i],id:id++,size:1,depth:0};cSize[i]=1}for(p=0;p<n-1;p++){c1=0;for(i=0;i<n;i++){if(distMatrix[i][dMin[i]]<distMatrix[c1][dMin[c1]])c1=i}c2=dMin[c1];c1Cluster=clusters[c1][0];c2Cluster=clusters[c2][0];var newCluster={left:c1Cluster,right:c2Cluster,dist:distMatrix[c1][c2],centroid:calculateCentroid(c1Cluster.size,c1Cluster.centroid,c2Cluster.size,c2Cluster.centroid),id:id++,size:c1Cluster.size+c2Cluster.size,depth:1+Math.max(c1Cluster.depth,c2Cluster.depth)};clusters[c1].splice(0,0,newCluster);cSize[c1]+=cSize[c2];for(j=0;j<n;j++){switch(linkage){case"single":if(distMatrix[c1][j]>distMatrix[c2][j])distMatrix[j][c1]=distMatrix[c1][j]=distMatrix[c2][j];break;case"complete":if(distMatrix[c1][j]<distMatrix[c2][j])distMatrix[j][c1]=distMatrix[c1][j]=distMatrix[c2][j];break;case"average":distMatrix[j][c1]=distMatrix[c1][j]=(cSize[c1]*distMatrix[c1][j]+cSize[c2]*distMatrix[c2][j])/(cSize[c1]+cSize[j]);break}}distMatrix[c1][c1]=Infinity;for(i=0;i<n;i++)distMatrix[i][c2]=distMatrix[c2][i]=Infinity;for(j=0;j<n;j++){if(dMin[j]==c2)dMin[j]=c1;if(distMatrix[c1][j]<distMatrix[c1][dMin[c1]])dMin[c1]=j}root=newCluster}return root}hcluster.linkage=function(x){if(!arguments.length)return linkage;linkage=x;return hcluster};hcluster.distance=function(x){if(!arguments.length)return distance;distance=x;return hcluster};hcluster.distanceMatrix=function(x){if(!arguments.length)return distMatrix;distMatrix=x.map(function(y){return y.slice(0)});return hcluster};return hcluster};function calculateCentroid(c1Size,c1Centroid,c2Size,c2Centroid){var newCentroid=[],newSize=c1Size+c2Size,n=c1Centroid.length,i=-1;while(++i<n){newCentroid[i]=(c1Size*c1Centroid[i]+c2Size*c2Centroid[i])/newSize}return newCentroid}reorder.leafOrder=function(){var distanceMatrix=null,distance=reorder.distance.euclidean,linkage="complete",debug=0,leavesMap={},orderMap={};function isLeaf(n){return n.depth==0}function leaves(n){if(n==null)return[];if(n.id in leavesMap)return leavesMap[n.id];return leavesMap[n.id]=_leaves(n)}function _leaves(n){if(n==null)return[];if(n.depth==0)return[n.id];return leaves(n.left).concat(leaves(n.right))}function order(v,i,j){var key="k"+v.id+"-"+i+"-"+j;if(key in orderMap)return orderMap[key];return orderMap[key]=_order(v,i,j)}function _order(v,i,j){if(v.depth==0)return[0,[v.id]];var l=v.left,r=v.right;var L=leaves(l),R=leaves(r);var w,x;if(L.indexOf(i)!=-1&&R.indexOf(j)!=-1){w=l;x=r}else if(R.indexOf(i)!=-1&&L.indexOf(j)!=-1){w=r;x=l}else throw{error:"Node is not common ancestor of "+i+", "+j};var Wl=leaves(w.left),Wr=leaves(w.right);var Ks=Wr.indexOf(i)!=-1?Wl:Wr;if(Ks.length==0)Ks=[i];var Xl=leaves(x.left),Xr=leaves(x.right);var Ls=Xr.indexOf(j)!=-1?Xl:Xr;if(Ls.length==0)Ls=[j];var max=Infinity,optimal_order=[];for(var k=0;k<Ks.length;k++){var w_max=order(w,i,Ks[k]);for(var m=0;m<Ls.length;m++){var x_max=order(x,Ls[m],j);var sim=w_max[0]+distanceMatrix[Ks[k]][Ls[m]]+x_max[0];if(sim<max){max=sim;optimal_order=w_max[1].concat(x_max[1])}}}return[max,optimal_order]}function orderFull(v){leavesMap={};orderMap={};var max=Infinity,optimal_order=[],left=leaves(v.left),right=leaves(v.right);if(debug)console.log(reorder.printhcluster(v,0));for(var i=0;i<left.length;i++){for(var j=0;j<right.length;j++){var so=order(v,left[i],right[j]);if(so[0]<max){max=so[0];optimal_order=so[1]}}}distanceMatrix=null;return optimal_order}function leafOrder(vector){if(distanceMatrix==null)distanceMatrix=reorder.dist().distance(distance)(vector);var hcluster=science.stats.hcluster().linkage(linkage).distanceMatrix(distanceMatrix);return orderFull(hcluster(vector))}leafOrder.debug=function(x){if(!arguments.length)return debug;debug=x;return leafOrder};leafOrder.distance=function(x){if(!arguments.length)return distance;distance=x;distanceMatrix=null;return leafOrder};leafOrder.linkage=function(x){if(!arguments.length)return linkage;linkage=x;return leafOrder};leafOrder.distanceMatrix=function(x){if(!arguments.length)return distanceMatrix;distanceMatrix=x.map(function(y){return y.slice(0)});return leafOrder};leafOrder.orderFull=orderFull;return leafOrder};reorder.order=function(){var distance=reorder.distance.euclidean,ordering=reorder.leafOrder,linkage="complete",distanceMatrix,vector,except=[],debug=0,i=0,j=Infinity;function _reset(){distance=reorder.distance.euclidean;ordering=reorder.leafOrder;linkage="complete";distanceMatrix=null;vector=null;except=[];debug=0;i=0;j=Infinity}function order(v){vector=v;j=Math.min(j,v.length);var i0=i>0?i-1:0,j0=j<vector.length?j+1:j,k,low,high;for(k=except.length-1;k>0;k-=2){low=except[k-1];high=except[k];if(high>=j0){if(j0>j){j0=Math.min(j0,low+1);except.splice(k-1,2)}else{high=j0}}else if(low<=i0){if(i0<i){i0=Math.max(i0,high-1);except.splice(k-1,2)}else{low=i0}}else if(high-low<3)except.splice(k-1,2)}try{return _order_limits(i0,j0)}finally{_reset()}}function _order_limits(i0,j0){var orig=vector,perm,row,k,l;vector=vector.slice(i0,j0);if(i==0&&j==vector.length)return _order_except();if(debug)console.log("i0="+i0+" j0="+j0);if(distanceMatrix!=null){if(j0!=vector.length)reorder.dist_remove(distanceMatrix,j0,vector.length);if(i0>0)reorder.dist_remove(distanceMatrix,0,i0)}else{_compute_dist()}var max=reorder.distmax(distanceMatrix);if(i0<i){row=distanceMatrix[0];for(k=row.length;k-->1;)row[k]+=max;for(k=distanceMatrix.length;k-->1;)distanceMatrix[k][0]+=max;max+=max;if(i0!=0){for(k=0;k<except.length;k++)except[k]-=i0}}if(j0>j){l=distanceMatrix.length-1;row=distanceMatrix[l];for(k=l;k-->0;){row[k]+=max;distanceMatrix[k][l]+=max}}perm=_order_except();if(i0<i){if(perm[0]!=0)perm.reverse();if(j0>j){reorder.assert(perm[0]==0&&perm[perm.length-1]==perm.length-1,"Invalid constrained permutation endpoints")}else{reorder.assert(perm[0]==0,"Invalid constrained permutation start")}}else if(j0>j){if(perm[perm.length-1]!=perm.length-1)perm=perm.reverse();reorder.assert(perm[perm.length-1]==perm.length-1,"Invalid constrained permutation end")}if(i0!=0){perm=reorder.permutation(i0).concat(perm.map(function(v){return v+i0}))}if(orig.length>j0){perm=perm.concat(reorder.range(j0,orig.length))}return perm}function _order_except(){var perm,k,l,low,high,pos;if(except.length==0)return _order_equiv();_compute_dist();for(k=except.length-1;k>0;k-=2){low=except[k-1];high=except[k];distanceMatrix=reorder.dist_remove(distanceMatrix,low+1,high-1);vector.splice(low+1,high-low-2);if(debug)console.log("Except["+low+", "+high+"]");if(distanceMatrix[low][low+1]!=0){distanceMatrix[low][low+1]=distanceMatrix[low+1][low]=-1}}perm=_order_equiv();for(k=0;k<except.length;k+=2){low=except[k];high=except[k+1];for(l=0;l<perm.length;l++){if(perm[l]>low)perm[l]+=high-low-2;else if(perm[l]==low)pos=l}if(pos>0&&perm[pos-1]==high-1){Array.prototype.splice.apply(perm,[pos,0].concat(reorder.range(high-2,low,-1)))}else if(perm[pos+1]==high-1){Array.prototype.splice.apply(perm,[pos+1,0].concat(reorder.range(low+1,high-1)))}else{throw"Range not respected"}}return perm}function _order_equiv(){var perm,row,e,j,k,l,m,n,has_1=false,equiv=[],fix_except={};_compute_dist();for(k=0;k<distanceMatrix.length-1;k++){row=distanceMatrix[k];e=[];j=row.indexOf(-1);if(j!=-1){fix_except[k]=[k,j];has_1=true}for(l=row.length;--l>k;){if(row[l]==0){j=distanceMatrix[l].indexOf(-1);if(j!=-1){fix_except[k]=[l,j];distanceMatrix[j][k]=row[j]=-1;has_1=true}e.unshift(l);distanceMatrix=reorder.dist_remove(distanceMatrix,l);vector.splice(l,1)}else if(row[l]<0)has_1=true}if(e.length!=0){e.unshift(k);equiv.push(e)}}if(has_1){for(k=0;k<distanceMatrix.length-1;k++){row=distanceMatrix[k];for(l=k+1;l<row.length-1;l++){if(distanceMatrix[l][l+1]==-1){distanceMatrix[l+1][l]=distanceMatrix[l][l+1]=0}}}}perm=_order();for(k=equiv.length;k-->0;){e=equiv[k];l=perm.indexOf(e[0]);m=fix_except[e[0]];if(m&&m[0]==e[0]){l=_fix_exception(perm,l,m[0],m[1],0);m=undefined}for(n=1;n<e.length;n++){perm=_perm_insert(perm,l,e[n]);if(m&&m[0]==e[n]){l=_fix_exception(perm,l,m[0],m[1],n);m=undefined}}}return perm}function _fix_exception(perm,l,m,next,len){var i,j,k;if(l>0&&perm[l-1]==next){_swap(perm,l,perm.indexOf(m));return l+1}else if(perm[l+len+1]==next){_swap(perm,l+len,perm.indexOf(m));return l}else throw"Index not found"}function _swap(perm,a,b){if(a==b)return;var c=perm[a];perm[a]=perm[b];perm[b]=c}function _order(){if(debug>1)reorder.printmat(distanceMatrix);if(debug>2)reorder.printmat(vector);var perm=ordering().debug(debug).linkage(linkage).distanceMatrix(distanceMatrix)(vector);if(debug)console.log("Permutation: "+perm);return perm}function _perm_insert(perm,i,nv){perm=perm.map(function(v){return v<nv?v:v+1});perm.splice(i,0,nv);return perm}order.debug=function(x){if(!arguments.length)return debug;debug=x;return order};function _compute_dist(){if(distanceMatrix==null)distanceMatrix=reorder.dist().distance(distance)(vector);return distanceMatrix}order.distance=function(x){if(!arguments.length)return distance;distance=x;return order};order.linkage=function(x){if(!arguments.length)return linkage;linkage=x;return order};order.limits=function(x,y){if(!arguments.length)return[i,j];i=x;j=y;return order};order.except=function(list){if(!arguments.length)return except.slice(0);except=list.sort(function(a,b){if(a>=b)throw"Invalid list, indices not sorted";return a-b});return order};function _orderExcept(vector,i,j){var distanceMatrix=reorder.dist().distance(distance)(vector);var row,k,l,rev=false,args,pos=-1;distanceMatrix[i][i+1]=0;distanceMatrix[i+1][i]=0;var perm=ordering().distanceMatrix(distanceMatrix)(vector);pos=perm.indexOf(i);for(k=0;k<perm.length;k++){l=perm[k];if(l>i)perm[k]+=j-i-2}if(pos!=0&&perm[pos-1]==j-1)rev=true;if(rev){perm.reverse();pos=perm.length-pos-1}args=[pos+1,0].concat(reorder.range(i+1,j-1));Array.prototype.splice.apply(perm,args);return perm}order.orderrowsexcept=order.orderexcept;return order};reorder.covariance=science.lin.dot;reorder.covariancetranspose=function(v,a,b){var n=v.length,cov=0,i;for(i=0;i<n;i++){cov+=v[i][a]*v[i][b]}return cov};reorder.variancecovariance=function(v){var o=v[0].length,cov=Array(o),i,j;for(i=0;i<o;i++){cov[i]=Array(o)}for(i=0;i<o;i++){for(j=i;j<o;j++)cov[i][j]=cov[j][i]=reorder.covariancetranspose(v,i,j)}return cov};function normalize(v){var norm=science.lin.length(v),i=v.length;while(i-->0)v[i]/=norm;return v}reorder.poweriteration=function(v,eps){if(arguments.length<2)eps=1e-4;var n=v.length,b=Array(n),i,j,tmp=Array(n),norm,s=10;reorder.assert(n==v[0].length,"poweriteration needs a square matrix");for(i=0;i<n;i++)b[i]=Math.random();b=normalize(b);while(s-->0){for(i=0;i<n;i++){tmp[i]=0;for(j=0;j<n;j++)tmp[i]+=v[i][j]*b[j]}tmp=normalize(tmp);var t=tmp;tmp=b;b=t}return tmp};reorder.sortorder=function(v){return reorder.range(0,v.length).sort(function(a,b){return v[a]-v[b]})};reorder.center=function(v){var n=v.length;if(n==0)return null;var mean=reorder.meancolumns(v),o=mean.length,v1=Array(n),i,j,row;for(i=0;i<n;i++){row=v[i].slice(0);for(j=0;j<o;j++){row[j]-=mean[j]}v1[i]=row}return v1};reorder.pca1d=function(v,eps){if(arguments.length<2)eps=1e-4;var n=v.length;if(v.length==0)return null;v=reorder.center(v);var cov=reorder.variancecovariance(v);return reorder.poweriteration(cov,eps)};reorder.pca1dorder=function(v,eps){return reorder.sortorder(pca1d(v,eps))};reorder.sumlines=function(v){var n=v.length,o=v[0].length,sumline=Array(n),i,j,row,s;for(i=0;i<n;i++){row=v[i];s=0;for(j=0;j<o;j++)s+=row[j];sumline[i]=s}return sumline};reorder.sumcols=function(v){var n=v.length,o=v[0].length,sumcol=science.zeroes(o),i,j,row;for(i=0;i<n;i++){row=v[i];for(j=0;j<o;j++)sumcol[j]+=row[j]}return sumcol};reorder.ca=function(v,eps){if(arguments.length<2)eps=1e-4;var n=v.length,o=v[0].length,sumline=reorder.sumlines(v),sumcol=reorder.sumcols(v),s=reorder.sum(sumcol),i,j,row;for(i=0;i<n;i++){v[i]=v[i].map(function(a){return a/s})}sumline=reorder.sumlines(v);sumcol=reorder.sumcols(v);var v2=Array(n),ep;for(i=0;i<n;i++){v2[i]=Array(o);for(j=0;j<o;j++){ep=sumline[i]*sumcol[j];v2[i][j]=(v[i][j]-ep)/Math.sqrt(ep)}}for(i=0;i<n;i++){for(j=0;j<o;j++){ep=sumline[i]*sumcol[j];v[i][j]=(v[i][j]-ep)/(sumcol[j]*Math.sqrt(sumline[i]))}}var cov=Array(n);for(i=0;i<n;i++)cov[i]=Array(n);for(i=0;i<n;i++)for(j=i;j<n;j++)cov[i][j]=cov[j][i]=reorder.covariance(v2[i],v2[j]);var eigenvector1=reorder.poweriteration(cov,eps),eigenvalue1=0;for(i=0;i<n;i++)eigenvalue1+=eigenvector1[i]*cov[i][0];eigenvalue1/=eigenvector1[0];return eigenvector1}})(this);
\ No newline at end of file
diff --git a/src/barycenter.js b/src/barycenter.js
index c53d2ae..4963544 100644
--- a/src/barycenter.js
+++ b/src/barycenter.js
@@ -32,19 +32,20 @@ function median(neighbors) {
 	return (neighbors[lm]*rspan + neighbors[rm]*lspan) / (lspan+rspan);
 }
 
-reorder.barycenter1 = function(graph, comp, iter) {
+reorder.barycenter1 = function(graph, comp, max_iter) {
     var nodes = graph.nodes(),
-	layer1, layer2,
-	layer,
+	layer1, layer2, crossings, iter,
+	best_layer1, best_layer2, best_crossings, best_iter,
+	layer, 
 	i, v, neighbors;
 
     if (comp.length < 3)
 	return comp;
 
-    if (! iter)
-	iter = 24;
-    else if ((iter%2)==1)
-	iter++; // want even number of iterations
+    if (! max_iter)
+	max_iter = 24;
+    else if ((max_iter%2)==1)
+	max_iter++; // want even number of iterations
 
     layer1 = comp.filter(function(n) {
 	return graph.outEdges(n).length!=0;
@@ -56,9 +57,14 @@ reorder.barycenter1 = function(graph, comp, iter) {
     for (i = 0; i < layer2.length; i++)
 	nodes[layer2[i]].pos = i;
 
-    for (layer = layer1;
-	 iter--;
-	 layer = (layer == layer1) ? layer2 : layer1) {
+    best_crossings = count_crossings(graph, layer1, layer2);
+    best_layer1 = layer1;
+    best_layer2 = layer2;
+    best_iter = 0;
+
+    for (layer = layer1, iter = 0;
+	 iter < max_iter;
+	 iter++, layer = (layer == layer1) ? layer2 : layer1) {
 	for (i = 0; i < layer.length; i++) {
 	    // Compute the median/barycenter for this node and set
 	    // its (real) value into node.mval
@@ -88,6 +94,14 @@ reorder.barycenter1 = function(graph, comp, iter) {
 	});
 	for (i = 0; i < layer2.length; i++)
 	    nodes[layer[i]].pos = i;
+	crossings = count_crossings(graph, layer1, layer2);
+	if (crossings < best_crossings) {
+	    best_crossings = crossings;
+	    best_layer1 = layer1;
+	    best_layer2 = layer2;
+	    best_iter = iter;
+	}
     }
-    return [layer1, layer2];
+    console.log('Best iter: '+best_iter);
+    return [best_layer1, best_layer2, best_crossings];
 };
diff --git a/src/count_crossings.js b/src/count_crossings.js
new file mode 100644
index 0000000..6d3cd86
--- /dev/null
+++ b/src/count_crossings.js
@@ -0,0 +1,56 @@
+// Wilhelm Barth, Petra Mutzel, Michael Jünger: 
+// Simple and Efficient Bilayer Cross Counting.
+// J. Graph Algorithms Appl. 8(2): 179-194 (2004)
+function count_crossings(graph, north, south) {
+    var i, j, n,
+	firstIndex, treeSize, tree, index, weightSum,
+	invert = false, crosscount;
+
+    // Choose the smaller axis
+    if (north.length < south.length) {
+	var tmp = north;
+	north = south;
+	south = tmp;
+	invert = true;
+    }
+
+    var south_inv = inverse_permutation(south),
+	southsequence = [];
+
+    for (i = 0; i < north.length; i++) {
+	if (invert) {
+	    n = graph.inEdges(north[i])
+		.map(function(e) {
+		    return south_inv[e.target.index];
+		});
+	}
+	else {
+	    n = graph.outEdges(north[i])
+		.map(function(e) {
+		    return south_inv[e.source.index];
+		});
+	}
+	n.sort(function(a,b){ return a-b; });
+	southsequence = southsequence.concat(n);
+    }
+    
+    firstIndex = 1;
+    while (firstIndex < south.length)
+	firstIndex <<= 1;
+    treeSize = 2 * firstIndex - 1;
+    firstIndex -= 1;
+    tree = science.zeroes(treeSize);
+
+    crosscount = 0;
+    for (i = 0; i < southsequence.length; i++) {
+	index = southsequence[i] + firstIndex;
+	tree[index]++;
+	while (index > 0) {
+	    if (index%2) crosscount += tree[index+1];
+	    index = (index - 1) >> 1;
+	    tree[index]++;
+	}
+    }
+    return crosscount;
+}
+reorder.count_crossings = count_crossings;
diff --git a/src/package.js b/src/package.js
index 705916f..4fe9b8d 100644
--- a/src/package.js
+++ b/src/package.js
@@ -23,6 +23,9 @@ require("util").puts(JSON.stringify({
 	"vows": "0.8.1",
 	"jsonfile": "2.0.0"
     },
+    "scripts": {
+	"test": "vows; echo"
+    },
     "licenses": [
 	{
 	    "type": "BSD-3",
diff --git a/test/barycenter-test.js b/test/barycenter-test.js
index 394a185..65c001a 100644
--- a/test/barycenter-test.js
+++ b/test/barycenter-test.js
@@ -19,6 +19,8 @@ suite.addBatch({
 		expect = [0, 1, 3, 2, 4];
 	    reorder.printmat(mat);
 	    perm = reorder.barycenter(graph);
+	    console.log('VOrder: %j, HOrder: %j, Crossings: %d',
+			perm[1], perm[0], perm[2]);
 	    reorder.printmat(mat, perm[1], perm[0]);
 	    //assert.deepEqual(perm, expect);
 	}
diff --git a/test/count_crossings-test.js b/test/count_crossings-test.js
new file mode 100644
index 0000000..a38666b
--- /dev/null
+++ b/test/count_crossings-test.js
@@ -0,0 +1,142 @@
+require("science");
+require("../reorder.v1");
+require("../reorder.v1");
+
+var vows = require("vows"),
+    assert = require("assert"),
+    seedrandom = require('seedrandom');
+//var    jf = require('jsonfile');
+
+Math.seedrandom('reorder');
+
+var suite = vows.describe("reorder.count_crossings");
+
+function naive_count_crossings(graph, north, south) {
+    var i, j, e1, e2, v, count = 0,
+	inv_north = reorder.inverse_permutation(north),
+	inv_south = reorder.inverse_permutation(south),
+	links = [];
+    for (i = 0; i < north.length; i++) {
+	v = north[i];
+	links = links.concat(graph.outEdges(v).map(function(e) {
+	    return [ inv_north[e.target.index], inv_south[e.source.index] ];
+	}));
+    }
+    for (i = 0; i < links.length; i++) {
+	e1 = links[i];
+	for (j = i+1; j < links.length; j++) {
+	    e2 = links[j];
+	    if ((e1[0] < e2[0] && e1[1] > e2[1]) 
+		|| (e1[0] > e2[0] && e1[1] < e2[1])) 
+		count++;
+	}
+    }
+    return count;
+}
+
+function dohard(mat) {
+    var graph = reorder.mat2graph(mat, true),
+	comps = graph.components(),
+	comp = comps.reduce(function(a, b) {
+	    return  (a.length > b.length) ? a : b;
+	});
+    comp.sort(function(a,b){return a-b; });
+    var layer1 = comp.filter(function(n) {
+	return graph.outEdges(n).length!=0;
+    }),
+	layer2 = comp.filter(function(n) {
+	    return graph.inEdges(n).length!=0;
+	});
+    //console.time('fast_crossings');
+    var c1 = reorder.count_crossings(graph, layer1, layer2);
+    //console.timeEnd('fast_crossings');
+    //console.time('naive_crossings');
+    var c2 = naive_count_crossings(graph, layer1, layer2);
+    //console.timeEnd('naive_crossings');
+    // if (c1 != c2) {
+    // 	var file = 'error_count_crossings.json';
+    // 	jf.writeFile(file, mat, function(err) {
+    // 	    console.log(err);
+    // 	});
+    // }
+    assert.equal(c1, c2);
+}
+
+suite.addBatch({
+    "count_crossings": {
+	"simple": function() {
+	    var graph = reorder.graph()
+		    .nodes([{id: 0}, {id: 1}, {id: 2}, {id: 3}, {id: 4}, {id: 5}])
+		    .links([{source: 0, target: 0},
+			    {source: 1, target: 1},
+			    {source: 1, target: 2},
+			    {source: 2, target: 0},
+			    {source: 2, target: 3},
+			    {source: 2, target: 4},
+			    {source: 3, target: 0},
+			    {source: 3, target: 2},
+			    {source: 4, target: 3},
+			    {source: 5, target: 2},
+			    {source: 5, target: 4}])
+		    .directed(true)
+		    .init(),
+		comp = graph.components()[0],
+		layer1 = comp.filter(function(n) {
+		    return graph.outEdges(n).length!=0;
+		}),
+		layer2 = comp.filter(function(n) {
+		    return graph.inEdges(n).length!=0;
+		});
+	    
+	    
+	    assert.equal(naive_count_crossings(graph, layer1, layer2),
+			 12);
+	    assert.equal(reorder.count_crossings(graph, layer1, layer2),
+			 12);
+	},
+	"bug": function() {
+	    dohard([
+    [0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0],
+    [0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+    [1, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
+    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0],
+    [1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    [0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    [1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0],
+    [0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
+    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1],
+    [0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1],
+    [0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0],
+    [0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    [0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 1],
+    [0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
+    [0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+    [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1],
+    [1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    [0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    [0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0],
+    [0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    [0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+    [1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0],
+    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+    [0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1],
+    [0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+    [1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0, 0],
+    [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0],
+    [0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0]
+	    ]);
+	},
+	"hard": function() {
+	    for (var i = 10; i < 100; i += 20) {
+		for (var j = 10; j < 100; j += 20) {
+		    var mat = reorder.randomMatrix(0.2, i, j, false);
+		    dohard(mat);
+		}
+	    }
+	}
+    }
+});
+
+suite.export(module);