# First of all we need to load the squangle polynomials.

squang1 <- scan("sqg1.txt")
squang2 <- scan("sqg2.txt")
squang3 <- scan("sqg3.txt")

# Now we load the packages we will require.

library(tensor) 

# "polyeval" is a function which evaluates homogeneous polynomials represented 
# as a vector.

polyeval <- function(poly,values,degree){

	len <- length(poly)
	
	seqs <- list()
	for(i in 1:(degree+1)) seqs[[i]] <- seq(i,len,by=(degree+1))

	coeffs <- poly[seqs[[1]]]
	
	variables <- list()
	for(i in 1:degree) variables[[i]] <- poly[seqs[[i+1]]]
	
	evalvar <- list()
	for(i in 1:degree) evalvar[[i]] <- values[variables[[i]]]

	for(i in 1:degree) coeffs <- coeffs*evalvar[[i]]
	
	final <- sum(coeffs)
	final
}

# To improve the speed of working with simulations of many runs we load
# the full sequence file using scan.

read.seq <- function(seq.file){scan(seq.file,what="character")}	

# "quartpatt" computes the pattern frequencies of the nth quartet of sequences
# from the value of read.seq: "allseq".

quartpatt <- function(allseq,vec=c(1,2,3,4)){
	
	seqn <- 2*vec+2
	quartet <- as.list(allseq[seqn])
	quartet2 <- list()
	for(i in 1:length(quartet)) quartet2[[i]] <- c(strsplit(quartet[[i]],split=vector()),recursive=TRUE)
	names(quartet) <- c(allseq[2*vec[1]+1],allseq[2*vec[2]+1],allseq[2*vec[3]+1],allseq[2*vec[4]+1])
	almost <- table(quartet2)
	final <- almost[1:4,1:4,1:4,1:4]
	final
}

quartpatt.boot <- function(allseq,vec=c(1,2,3,4)){
	
	seqn <- 2*vec+2
	quartet <- as.list(allseq[seqn])
	quartet2 <- list()
	for(i in 1:length(quartet)) quartet2[[i]] <- c(strsplit(quartet[[i]],split=vector()),recursive=TRUE)
	len <- length(quartet2[[1]])
	samp <- sample(1:len,replace=TRUE)
	for(i in 1:length(quartet)) quartet2[[i]] <- quartet2[[i]][samp]
	names(quartet) <- c(allseq[2*vec[1]+1],allseq[2*vec[2]+1],allseq[2*vec[3]+1],allseq[2*vec[4]+1])
	almost <- table(quartet2)
	final <- almost[1:4,1:4,1:4,1:4]
	final
}

            
# "squangle" is a function which computes the squangles from "allseq". Included
# transformation of "patterns" data to the basis used in the squangle polynomials. 
# This is executed as a tensor multiplication with the operator "trans" using the 
# R-package "tensor".

squangle <- function(allseq,vec=c(1,2,3,4)){

	patt <- quartpatt(allseq,vec)

	trans <- diag(4)
	trans[1,] <- 1

	pattern1 <- tensor(trans,patt,c(2),c(1))
	pattern2 <- tensor(trans,pattern1,c(2),c(2))
	pattern3 <- tensor(trans,pattern2,c(2),c(3))
	pattern4 <- tensor(trans,pattern3,c(2),c(4))

	finalpatt <- vector()
	for(i in 1:4) for(j in 1:4) for(k in 1:4) for(l in 1:4) finalpatt[(l-1)*4^0+(k-1)*4^1+(j-1)*4^2+(i-1)*4^3+1] <- pattern4[i,j,k,l]

  len <- sum(patt)

	ttt1 <- polyeval(squang1,finalpatt,5)/len^5
	ttt2 <- polyeval(squang2,finalpatt,5)/len^5
	ttt3 <- polyeval(squang3,finalpatt,5)/len^5

	c(-(3/2)*ttt1 + ttt2 + 2*ttt3,-(3/2)*ttt1 + 2*ttt2 + ttt3,-ttt2+ttt3)
}

squangle.boot <- function(allseq,vec=c(1,2,3,4)){

	patt<-quartpatt.boot(allseq,vec)

	trans <- diag(4)
	trans[1,] <- 1

	pattern1 <- tensor(trans,patt,c(2),c(1))
	pattern2 <- tensor(trans,pattern1,c(2),c(2))
	pattern3 <- tensor(trans,pattern2,c(2),c(3))
	pattern4 <- tensor(trans,pattern3,c(2),c(4))

	finalpatt<- vector()
	for(i in 1:4) for(j in 1:4) for(k in 1:4) for(l in 1:4) finalpatt[(l-1)*4^0+(k-1)*4^1+(j-1)*4^2+(i-1)*4^3+1]<-pattern4[i,j,k,l]

  len <- sum(patt)

	ttt1 <- polyeval(squang1,finalpatt,5)/len^5
	ttt2 <- polyeval(squang2,finalpatt,5)/len^5
	ttt3 <- polyeval(squang3,finalpatt,5)/len^5

	c(-(3/2)*ttt1 + ttt2 + 2*ttt3,-(3/2)*ttt1 + 2*ttt2 + ttt3,-ttt2+ttt3)
}


# The following defines the likelihood routine. We assume
# the squangles are stochastically independent and normally distributed. 

# "tree" finds the maximum likelihood tree 
# given squangles "x".
# tree1=(12,34), tree2=(13,24), tree3=(14,23)

tree <- function(x) {

      u1 <- max(c(0,(x[2]-x[3])/2))
      u2 <- max(c(0,(x[1]+x[3])/2))
      u3 <- max(c(0,-(x[1]+x[2])/2))
        
      if (u1==0) e1 <- Inf else e1 <- x[1]^2+(x[2]-u1)^2+(x[3]+u1)^2
      if (u2==0) e2 <- Inf else e2 <- x[2]^2+(x[1]-u2)^2+(x[3]-u2)^2
      if (u3==0) e3 <- Inf else e3 <- x[3]^2+(x[1]+u3)^2+(x[2]+u3)^2
    
      
      ee <-  c(e1,e2,e3)
      
      ans <- which(ee==min(ee))
      paste("tree",ans,sep="")

}



# Here we define "write.squang" and "read.squang" for convenience. 
# "squangs" is a list of n squangle vectors.

write.squang<-function(squangs,file.name){

	squangs2<-c(squangs,recursive=TRUE)
	write(squangs2,file.name)
	print("squangles written to file")

}

read.squang<-function(file.name){

	squangs1<-scan(file.name)
	squangs2<-list()
	n<-length(squangs1)/3
	for(i in 1:n) squangs2[[i]]<-squangs1[(3*i-2):(3*i)]
	squangs2

}

# AIC computes the Aikake weights. We normalize to the interval [0,.5].

aicW<-function(squangs){

  x<-squangs
  if(sum(abs(x))==0) (1/3)*c(1,1,1)  else
  {likes<-tree.loglike(x)
  aic<- -likes
  mi<-min(aic)
  del<-aic-mi
  weight<-exp(-del)/(sum(exp(-del)))
  weight}
  
}

# Now we print some comments.

print("Squangle package loaded.",quote=FALSE)