graph.splitSignature <- function(sgn) {
 L <- sgn[[1]];
 S <- sgn[[2]];
 possibilities <- 1;
 if (length(L) < 2) {
   return(1);
 }
 while ( length(which( L == 2 )) > 0 ) {
   if ( length(which( L == 2 )) > 1 ) {
     long_seq_len <- sum(S[which(L == 2)]);
     for (k in which(L == 2)) {
       possibilities <- possibilities * choose(long_seq_len, S[k]);
       long_seq_len <- long_seq_len - S[k];
     }
     long_seq_len <- sum(S[which(L == 2)]);
     for (k in which(L == 1 & S > 1)) {
       possibilities <- possibilities / choose(long_seq_len, S[k] - 1);
       long_seq_len <- long_seq_len - S[k] + 1;
     }
   }
   mins <- which.min(L);
   L <- L[-mins];
   S <- S[-mins];
   L <- L - 1;
 }
 return( possibilities );
}

# Intern: Calculates actual signature
graph.compressToSignature <- function(Expression, Model) {
 perm <- sort(Expression, index.return=TRUE, decreasing=TRUE)[[2]];
 A <- diag(rep(1, ncol(Model)));
 sum <- c(rep(0, ncol(Model)));
 M <- Model;
 for (i in perm) {
   sum <- sum + A[,i];
   A <- pmax(A, (A[,i] %*% t(M[,i])));
   adjacent <- which(M[,i] == 1);
   for (j in adjacent) {
     for (k in adjacent) {
       if (k != j) {
         M[j,k] <- 1;
       }
     }
   }
   M[,i] <- 0;
   M[i,] <- 0;
 }
 return(list(c(sum), colSums(A)));
}


# Compute probability of a particular data row
graph.pr <- function(ExpressionSeries, Model) {
  signature <-  graph.compressToSignature(ExpressionSeries,Model);
  nfactorial <- factorial(length(ExpressionSeries));
  cellsize <- graph.splitSignature(signature);
  return(cellsize/nfactorial);
}


#
cyclohedronP<-function(delta,n) {
  S<-list(c(delta[1]))
  value<-1
  for (m in 2:(n-1)) {
	um <- delta[m]+1
	if (um>n) um<-1
	dm <- delta[m]-1
	if (dm<1) dm<-n
	Left <- c()
	Right <- c()	
	for (s in S) {
		if (um %in% s) Left<-s
		if (dm %in% s) Right<-s
	}
	S <- c(S, list(union(c(delta[m]),union(Left,Right))))
	value<-value*choose(length(Left)+length(Right),length(Left))
  	}
  return(value)
  }

#Make graphs
buildGraph <- function(n, type) {
  if (type=='line') {
    model_line <- t(array(c(c(0,1,0),rep(0,n-3),c(1)), dim = c(n,n)));
    return(model_line);
  }
  if (type=='circle') {
    model_circle <- buildGraph(n,'line');
    model_circle[1, n] <- 1;
    model_circle[n, 1] <- 1;
  }
  return(model_circle);
}

# Compute preimages of a given dataset wrt G
preimagesizes <-function(dataframe,G,sortedlog=FALSE) {
  dset <-as.matrix(dataframe)
  preimagesizevector <- matrix(0,nrow(dset),1)
  for (i in 1:nrow(dset)) {
    preimagesizevector[i,1]<-graph.splitSignature(graph.compressToSignature(dset[i,],G))
  }
  if (sortedlog) {
    return(log(sort(preimagesizevector, decreasing=FALSE))) 
  }
 return(preimagesizevector)
}


# Compute preimages of a given dataset wrt G
cycleCounts <-function(dataframe,sortedlog=FALSE) {
  dset <-as.matrix(dataframe)
  n<-ncol(dset)
  preimagesizevector <- matrix(0,nrow(dset),1)
  for (i in 1:nrow(dset)) {
    preimagesizevector[i,1]<-cyclohedronP(sort(dset[i,], index.return=TRUE, decreasing=TRUE)[[2]],n)
  }
#  if (sortedlog) {
#    return(log(sort(preimagesizevector, decreasing=FALSE))) 
#  }
 return(preimagesizevector)
}



#load data in csv format with one column of probsets for row names, one header row
#returns dataframe
loaddata <-function(filename) {
return (read.table(filename, row.names=1, header=TRUE, fill=TRUE, sep=",", strip.white=TRUE))
}


#load gene names and whatnot
loadnames <- function(filename,numcols=3) {
temp<-read.table(filename, header=TRUE, fill=TRUE, sep="\t", strip.white=TRUE)
return(as.matrix(temp)[,1:numcols])
}

# namedp<-data.frame(p,row.names=rownames(mydf))
#
rankby<-function(mynames,myvalues,decreasing=FALSE) {
ord<-order(myvalues, decreasing=decreasing)
mynames<-as.matrix(mynames)
myvalues<-as.matrix(myvalues)
N<-nrow(myvalues)
M<-ncol(mynames)
newnames<-matrix(0,nrow=N,ncol=M)
for (i in 1:N) {
newnames[i,]<-mynames[ord[i],]
}
return(newnames)
}


# Make random data, uniform on permutations 
uniformperms<-function(Nvertices,Nrows) {
  randomdata<-matrix(0,Nrows,Nvertices)
  for (i in 1:Nrows) {
    randomdata[i,]<-sample(1:Nvertices,Nvertices)
  }
  return(randomdata)
}


rhoofdata<-function(datav) {
  datav<-as.matrix(datav)
  delta <- sort(datav, index.return=TRUE, decreasing=TRUE)[[2]];
  n<- length(delta)
  rhov <- vector(mode="integer",n);
  for (i in 1:n){
    rhov[delta[i]]<-n-i+1
  }
  return(rhov)
}

rhoofdelta<-function(delta) {
  delta<-as.matrix(delta)
  n<- length(delta)
  rhov <- vector(mode="integer",n);
  for (i in 1:n){
    rhov[delta[i]]<-n-i+1
  }
  return(rhov)
}

descentperm<-function(datav) {
return(sort(datav, index.return=TRUE, decreasing=TRUE)[[2]])
}

#Plot the height of the topographic map, given a data permutation or raw data and a graph G
plottubing<-function(datav,G){
  gname<-rownames(datav)
  v<-as.matrix(datav)
  n<-length(v)
  plot(1:n,graph.compressToSignature(v,G)[[1]],'l',xlab='Position',ylab='Height',main=gname)
}