#functions for doing diurnal rhythm analyses
#
#
#a function to estimate diurnal phase of mood data
#the input is a data frame or matrix with 
#time of measurement (in 24 hour clock)
#and  then the mood measures (1 or many) 
#Version of October 22, 2008
#seriously revised April 12, 2009
#
#find the best fitting phase  (in hours)
"cosinor" <- 
function(angle,x=NULL,code=NULL,period=24,plot=FALSE,opti=FALSE) {
if (is.null(x)) {x <- angle[2:ncol(angle)]
                 angle<- angle[1]}   #this allows for easier input if the time is the first variable
nvar <- dim(x)
if(is.null(code)) { fit <- cosinor1(angle,x,period=period,opti=opti)  #if there is a code (i.e., subject id), then do the fits for each separate subject
m.resp <- mean(x[1])
s.resp <- sd(x[1])
 x <- NA   #just to avoid being told that x is a global
if(plot) {curve(cos((x-fit[1,1])*pi/12)*s.resp+m.resp,add=TRUE) } #this draws the first fitted variable
  }  else  { fit.list <- by(x,code,cosinor1,period=period,opti=opti)  #this is the case if code is specified
   ncases <- length(fit.list)
   fit<- matrix(unlist(fit.list),nrow=ncases,byrow=TRUE)
   colnames(fit) <- c(paste(colnames(x), "phase",sep="."),paste(colnames(x[2:(nvar+1)]), "fit",sep="."))  
 }
return(fit)
}

# cosinor1 actually does the work
# it either uses a fitting function (optimize) from core R
# or calls a linear regression fit
"cosinor1" <-
function(angle,x,period,opti) {
response <- x
n.var <- dim(x)[2] 
fit <- matrix(NA,nrow=n.var,ncol=2)
for (i in 1:n.var) {
if(opti) {fits <- optimize(f=phaser,c(0,24),time=angle,response=x[,i],period=period,maximum=TRUE)  #iterative fit 
fit[i,1] <- fits$maximum
fit[i,2] <- fits$objective
} else {fits <- cosinor.lm2 (angle,x[,i],period=24)  #simple linear regression based upon sine and cosine of time
   fit[i,1] <- 	fits[[1]]
   fit[i,2] <- fits[[2]]
   }
}
colnames(fit) <- c("phase","fit")
rownames(fit) <- colnames(x)
return(fit)
}

"phaser" <- function(phase,time,response,period) {       #this is used in the iterative fit procedure
   phaser <- cor(cos(((time-phase)*2*pi)/period),response,use="pairwise")}
   
   
#the alternative to the iterative procedure is simple linear regression of the cosine and sine    
   "cosinor.lm2" <- 
function(angle,y,period=24) {
p2 <- period/2

        cos.t <- cos(angle*pi/p2)
        sin.t <- sin(angle *pi/p2)
      
    dat.df <- data.frame(iv1=cos.t,iv2=sin.t,y)
    
    cor.dat <- cor(dat.df,use="pairwise")
    beta1 <- (cor.dat[3,1] - cor.dat[3,2] * cor.dat[1,2])/(cor.dat[1,1]-cor.dat[1,2]^2)
    beta2 <- (cor.dat[3,2] - cor.dat[3,1] * cor.dat[1,2])/(cor.dat[1,1]-cor.dat[1,2]^2)
  #  mod <- lm(y ~cos.t+ sin.t
  #  mod <- lm(y~ iv1+ iv2,data=dat.df)
  # phase <- atan(-mod$coefficients[3]/mod$coefficients[2])*12/pi
   phase <- ( sign(beta2) *acos( beta1/sqrt(beta1^2 + beta2^2)))*p2/pi 
   phase <- phase %% period
   #phase <- atan(beta2/beta1)*p2/pi + p2
   R <- sqrt(cor.dat[3,1]*beta1 + cor.dat[3,2]*beta2)
    fit <- list(phase=phase,R=R,beta1,beta2)
return(fit)
}

   
   
##
#
#find the mean phase of output from cosiner or any other circular data set
#can find the mean phase of data in radians or hours (default)
#
"circadian.mean" <- 
function(angle,hours=TRUE) {
if(hours) { angle <- angle*2*pi/24 }
x <- cos(angle)
y <- sin(angle)
if (is.vector(angle)) {
mx <- mean(x)
my <- mean(y) } else { 
mx <- colMeans(x)
my <- colMeans(y) }
mean.angle <- sign(my) * acos((mx)/sqrt(mx^2+my^2)) 
# mean.angle <- atan(my/mx)   #according to circular stats, but the other form is clearer
if (hours) {mean.angle <- mean.angle*24/(2*pi)
mean.angle[mean.angle <= 0] <-  mean.angle[mean.angle<=0] + 24}
return(mean.angle)
}


## The circular correlation matrix of phase data
#adapted from the circStats package
#with various modifications for the study of mood data 
#one change is not use atan but rather use cosine over length
#

"circadian.cor"  <- 
function(angle,hours=TRUE) {
if(hours) { angle <- angle*2*pi/24 }
nvar <- dim(angle)[2]
correl <- diag(nvar)
x <- cos(angle)
y <- sin(angle)
mx <- colMeans(x)
my <- colMeans(y) 
mean.angle <- sign(my) * acos((mx)/sqrt(mx^2+my^2)) 

for (i in 1:nvar) {#the logic is that larger deviations are weighted more, up to the sin(theta)
  for (j in 1:i) {covar <-  sum(sin(angle[,i] -mean.angle[i]) *sin(angle[,j] -mean.angle[j]))   
  correl[i,j] <- correl[j,i] <- covar}
}
sd <- diag(sqrt(1/diag(correl)))
correl <- sd %*% correl %*% sd
colnames(correl) <- rownames(correl) <- colnames(angle)
return(correl) }


#to find the correlation of a linear variable (e.g., a personality trait) with a circular one (e.g., phase)
"circadian.linear.cor" <- 
function(angle,x=NULL,hours=TRUE) {
if(hours) { angle <- angle*2*pi/24 }
if(is.null(x)) {x <- angle[2:dim(angle)[2]]
                  angle <- angle[1]}
cos.angle <- cos(angle)
sin.angle <- sin(angle)
cor.cos <- cor(cos.angle,x,use="pairwise")
cor.sin <- cor(sin.angle,x,use="pairwise")
if(!is.vector(angle)) {cor.cs <- diag(cor(cos.angle,sin.angle,use="pairwise"))} else {cor.cs <- cor(cos.angle,sin.angle,use="pairwise")}
R <- sqrt((cor.cos^2 + cor.sin^2 - 2 * cor.cos * cor.sin * cor.cs)/(1-cor.cs^2))*sign(cor.cos)
return(R) }


"circadian.plot" <-
function(angle,x=NULL,hours=TRUE,title="Polar plot") {
if(hours) { angle <- angle*2*pi/24 }
x1 <- cos(angle) * x
y1 <- sin(angle) * x
maxx <- max(x)
plot(x1,y1,axes=FALSE,xlab="",ylab="",xlim=c(-maxx,maxx),ylim=c(-maxx,maxx),asp=1) 
segments  <- 51
 angles <- (0:segments) * 2 * pi/segments
unit.circle <- cbind(cos(angles), sin(angles))
points(unit.circle*maxx,typ="l")
text(maxx,0,"24",pos=4)
text(-maxx,0,"12",pos=2)
text(0,maxx,"6",pos=3)
text(0,-maxx,"18",pos=1)
}

#deprecated
"circadian.linear.cor.2" <- 
function(angle,x,hours=TRUE) {
if(hours) { angle <- angle*2*pi/24 }
cos.angle <- cos(angle)
sin.angle <- sin(angle)
cor.cos <- cor(cos.angle,x,use="pairwise")
cor.sin <- cor(sin.angle,x,use="pairwise")
if(!is.vector(angle)) {cor.cs <- diag(cor(cos.angle,sin.angle))} else {cor.cs <- cor(cos.angle,sin.angle)}
R <- sqrt((cor.cos^2 + cor.sin^2 - 2 * cor.cos * cor.sin * cor.cs)/(1-cor.cs^2))
return(R) }