--- a/create_heatmap.R	Fri Nov 17 11:37:25 2017 -0500
+++ b/create_heatmap.R	Fri Nov 17 11:37:33 2017 -0500
@@ -11,25 +11,52 @@
 args <- parse_args(parser, positional_arguments=TRUE)
 opt <- args$options
 
-create_heatmap<-function(data_matrix, statecolor=NULL, markcolor=NULL, cols=c("white","dark blue")) {
-    k = dim(data_matrix)[2];
-    l = dim(data_matrix)[1];
+state_color <- function(statemean, markcolor=NULL)
+{    
+    if(length(markcolor) == 0) {
+        markcolor = rep("", dim(statemean)[2]);
+        markcolor[order(apply(statemean, 2, sd), decreasing=T)] = hsv((1:dim(statemean)[2]-1) / dim(statemean)[2], 1, 1);
+        markcolor = t(col2rgb(markcolor));
+    }
+
+    rg = apply(statemean, 1, range);
+    mm = NULL;
+    for(i in 1:dim(statemean)[1]) {
+        mm = rbind(mm, (statemean[i,] - rg[1,i] + 1e-10) / (rg[2,i] - rg[1,i] + 1e-10));
+    }
+    mm = mm^6;
+    if(dim(mm)[2] > 1) {
+        mm = mm / (apply(mm, 1, sum) + 1e-10);
+    }
+    mycol = mm%*%markcolor;
+    s = apply(statemean, 1, max);
+    s = (s - min(s)) / (max(s) - min(s) + 1e-10);
+    h = t(apply(mycol, 1, function(x){rgb2hsv(x[1], x[2], x[3])}));
+    h[,2] = h[,2] * s;
+    h = apply(h, 1, function(x){hsv(x[1], x[2], x[3])});
+    rt = cbind(apply(t(col2rgb(h)), 1, function(x){paste(x, collapse=",")}) ,h);
+    return(rt);
+}
+
+create_heatmap<-function(data_frame, statecolor=NULL, markcolor=NULL, cols=c("white","dark blue")) {
+    k = dim(data_frame)[2];
+    l = dim(data_frame)[1];
     p = (sqrt(9 + 8 * (k - 1)) - 3) / 2;
-    m = as.matrix(data_matrix[,1+1:p]/data_matrix[,1]);
-    colnames(m) = colnames(data_matrix)[1+1:p];
-    marks = colnames(m);
-    rownames(m) = paste(1:l-1," (", round(data_matrix[,1] / sum(data_matrix[,1]) * 10000) / 100, "%)", sep="");
+    data_matrix = as.matrix(data_frame[,1+1:p] / data_frame[,1]);
+    colnames(data_matrix) = colnames(data_frame)[1+1:p];
+    marks = colnames(data_matrix);
+    rownames(data_matrix) = paste(1:l," (", round(data_frame[,1] / sum(data_frame[,1]) * 10000) / 100, "%)", sep="");
     pdf(file=opt$output);
     par(mar=c(6,1,1,6));
-    rg = range(m);
+    rg = range(data_matrix);
     colors = 0:100/100*(rg[2]-rg[1])+rg[1];
     my_palette = colorRampPalette(cols)(n=100);
     defpalette = palette(my_palette);
 
     plot(NA, NA, xlim=c(0,p+0.7), ylim=c(0,l), xaxt="n", yaxt="n", xlab=NA, ylab=NA, frame.plot=F);
-    axis(1, at=1:p-0.5, labels=colnames(m), las=2);
-    axis(4, at=1:l-0.5, labels=rownames(m), las=2);
-    rect(rep(1:p-1,l), rep(1:l-1,each=p), rep(1:p,l), rep(1:l,each=p), col=round((t(m)-rg[1])/(rg[2]-rg[1])*100));
+    axis(1, at=1:p-0.5, labels=colnames(data_matrix), las=2);
+    axis(4, at=1:l-0.5, labels=rownames(data_matrix), las=2);
+    rect(rep(1:p-1,l), rep(1:l-1,each=p), rep(1:p,l), rep(1:l,each=p), col=round((t(data_matrix)-rg[1])/(rg[2]-rg[1])*100));
 
     if(length(statecolor)==0) {
         if(length(markcolor)==0) {
@@ -76,15 +103,13 @@
                 }
             }
         }
-        statecolor = stateColor(m, markcolor)[,2];
+        sc = state_color(data_matrix, markcolor)[,2];
     }
-    rect(rep(p+0.2,l), 1:l-0.8, rep(p+0.8,l), 1:l-0.2, col=statecolor);
-
+    rect(rep(p+0.2,l), 1:l-0.8, rep(p+0.8,l), 1:l-0.2, col=sc);
     palette(defpalette);
     dev.off();
-    #return(statecolor);
 }
 
 # Read the input.
-data_matrix <- read.table(opt$input, header=T);
-create_heatmap(data_matrix);
+data_frame <- read.table(opt$input, comment="!", header=T);
+create_heatmap(data_frame);