###################################################
### chunk number 1: load STAR and data etc
###################################################
library(STAR)
data(e060817citron)
data(e060817terpi)
data(e060817mix)
deltaCitron <- 6-attr(e060817citron[[1]],"stimTimeCourse")[1]
deltaTerpi <- 6-attr(e060817terpi[[1]],"stimTimeCourse")[1]
deltaMix <- 6-attr(e060817mix[[1]],"stimTimeCourse")[1]
allCitron <- lapply(1:20,
                    function(idx) {
                      st <- sort(c(e060817citron[[1]][[idx]],
                                   e060817citron[[2]][[idx]],
                                   e060817citron[[3]][[idx]])
                                 ) + deltaCitron
                      st <- st[1 < st & st < 14]
                      unique(st)
                    }
                    )
allCitron <- as.repeatedTrain(allCitron)
allTerpi <- lapply(1:20,
                   function(idx) {
                     st <- sort(c(e060817terpi[[1]][[idx]],
                                  e060817terpi[[2]][[idx]],
                                  e060817terpi[[3]][[idx]])
                                ) + deltaTerpi
                     st <- st[1 < st & st < 14]
                     unique(st)
                   }
                   )
allTerpi <- as.repeatedTrain(allTerpi)
allMix <- lapply(1:20,
                 function(idx) {
                   st <- sort(c(e060817mix[[1]][[idx]],
                                e060817mix[[2]][[idx]],
                                e060817mix[[3]][[idx]])
                              ) + deltaMix
                   st <- st[1 < st & st < 14]
                   unique(st)
                 }
                 )
allMix <- as.repeatedTrain(allMix)


###################################################
### chunk number 2: make raster figure
###################################################
png(file="CitronellalRaster.png",width=800,heigh=800)
par(cex=2)
plot(allCitron,
     stimTimeCourse=c(6,6.5),
     xlab="Temps (s)",
     ylab="Stimulation",
     main="Réponses au citronellal",
     xlim=c(4,14))
dev.off()


###################################################
### chunk number 3: get collapsed histograms
###################################################
## 50 ms bin width
allCitronRG_50 <- hist(unlist(allCitron),breaks=seq(1,14,0.05),plot=FALSE)
allTerpiRG_50 <- hist(unlist(allTerpi),breaks=seq(1,14,0.05),plot=FALSE)
allMixRG_50 <- hist(unlist(allMix),breaks=seq(1,14,0.05),plot=FALSE)
YYcitron_50 <- allCitronRG_50$counts/20/0.05
YYterpi_50 <- allTerpiRG_50$counts/20/0.05
YYmix_50 <- allMixRG_50$counts/20/0.05
bigMax_50 <- max(c(YYcitron_50,
                   YYterpi_50,
                   YYmix_50))
yMax_50 <- round(bigMax_50/1000,digits=1)*1000
XX_50 <- allCitronRG_50$mids

## 10 ms bin width
allCitronRG_10 <- hist(unlist(allCitron),breaks=seq(1,14,0.01),plot=FALSE)
allTerpiRG_10 <- hist(unlist(allTerpi),breaks=seq(1,14,0.01),plot=FALSE)
allMixRG_10 <- hist(unlist(allMix),breaks=seq(1,14,0.01),plot=FALSE)
YYcitron_10 <- allCitronRG_10$counts/20/0.01
YYterpi_10 <- allTerpiRG_10$counts/20/0.01
YYmix_10 <- allMixRG_10$counts/20/0.01
bigMax_10 <- max(c(YYcitron_10,
                   YYterpi_10,
                   YYmix_10))
yMax_10 <- round(bigMax_10/1000,digits=1)*1000
XX_10 <- allCitronRG_10$mids

## 250 ms bin width
allCitronRG_250 <- hist(unlist(allCitron),breaks=seq(1,14,0.25),plot=FALSE)
allTerpiRG_250 <- hist(unlist(allTerpi),breaks=seq(1,14,0.25),plot=FALSE)
allMixRG_250 <- hist(unlist(allMix),breaks=seq(1,14,0.25),plot=FALSE)
YYcitron_250 <- allCitronRG_250$counts/20/0.25
YYterpi_250 <- allTerpiRG_250$counts/20/0.25
YYmix_250 <- allMixRG_250$counts/20/0.25
bigMax_250 <- max(c(YYcitron_250,
                   YYterpi_250,
                   YYmix_250))
yMax_250 <- round(bigMax_250/1000,digits=1)*1000
XX_250 <- allCitronRG_250$mids

bigMax <- max(c(bigMax_50,bigMax_10,bigMax_250))
yMax <- round(bigMax/1000,digits=1)*1000


###################################################
### chunk number 4: psths 1
###################################################
ylim <- c(0,3*yMax)
## 50 ms bin width
png(file="allPSTHs_50.png",width=1000,height=1000)
par(mar=c(2,2,1,1),cex=2)
plot(XX_50,
     YYcitron_50,
     type="n",
     axes=FALSE,
     xlab="",
     ylab="",
     xlim=c(1,14),
     ylim=ylim)
rect(6,0,6.5,ylim[2],col="grey80",border=NA)
abline(h=50,col="grey30",lty=2)
abline(h=100,col="grey30",lty=2)
abline(h=150,col="grey30",lty=2)
lines(XX_50,
      YYcitron_50,
      type="s",
      lwd=2,
      col=2)
abline(h=0,col="grey30")
segments(0.5,100,0.5,150,lwd=2)
text(1.5,125,"50 (Hz)")
text(10,125,"Citronellal")
abline(h=50+yMax,col="grey30",lty=2)
abline(h=100+yMax,col="grey30",lty=2)
abline(h=150+yMax,col="grey30",lty=2)
lines(XX_50,
      YYterpi_50+yMax,
      type="s",
      lwd=2,
      col=2)
abline(h=yMax,col="grey30")
##segments(0.5,100+yMax,0.5,150+yMax,lwd=2)
##text(1,125+yMax,"50")
text(10,125+yMax,"Terpineol")
abline(h=50+2*yMax,col="grey30",lty=2)
abline(h=100+2*yMax,col="grey30",lty=2)
abline(h=150+2*yMax,col="grey30",lty=2)
lines(XX_50,
      YYmix_50+2*yMax,
      type="s",
      lwd=2,
      col=2)
##segments(0.5,100+2*yMax,0.5,150+2*yMax,lwd=2)
abline(h=2*yMax,col="grey30")
##text(1,125+2*yMax,"50")
text(10,125+2*yMax,"Mélange")
dev.off()

## 10 ms bin width
png(file="allPSTHs_10.png",width=1000,height=1000)
par(mar=c(2,2,1,1),cex=2)
plot(XX_10,
     YYcitron_10,
     type="n",
     axes=FALSE,
     xlab="",
     ylab="",
     xlim=c(1,14),
     ylim=ylim)
rect(6,0,6.5,ylim[2],col="grey80",border=NA)
abline(h=50,col="grey30",lty=2)
abline(h=100,col="grey30",lty=2)
abline(h=150,col="grey30",lty=2)
lines(XX_10,
      YYcitron_10,
      type="s",
      lwd=2,
      col=2)
abline(h=0,col="grey30")
segments(0.5,100,0.5,150,lwd=2)
text(1.5,125,"50 (Hz)")
text(10,125,"Citronellal")
abline(h=50+yMax,col="grey30",lty=2)
abline(h=100+yMax,col="grey30",lty=2)
abline(h=150+yMax,col="grey30",lty=2)
lines(XX_10,
      YYterpi_10+yMax,
      type="s",
      lwd=2,
      col=2)
abline(h=yMax,col="grey30")
##segments(0.5,100+yMax,0.5,150+yMax,lwd=2)
##text(1,125+yMax,"50")
text(10,125+yMax,"Terpineol")
abline(h=50+2*yMax,col="grey30",lty=2)
abline(h=100+2*yMax,col="grey30",lty=2)
abline(h=150+2*yMax,col="grey30",lty=2)
lines(XX_10,
      YYmix_10+2*yMax,
      type="s",
      lwd=2,
      col=2)
##segments(0.5,100+2*yMax,0.5,150+2*yMax,lwd=2)
abline(h=2*yMax,col="grey30")
##text(1,125+2*yMax,"50")
text(10,125+2*yMax,"Mélange")
dev.off()

## 250 ms bin width
png(file="allPSTHs_250.png",width=1000,height=1000)
par(mar=c(2,2,1,1),cex=2)
plot(XX_250,
     YYcitron_250,
     type="n",
     axes=FALSE,
     xlab="",
     ylab="",
     xlim=c(1,14),
     ylim=ylim)
rect(6,0,6.5,ylim[2],col="grey80",border=NA)
abline(h=50,col="grey30",lty=2)
abline(h=100,col="grey30",lty=2)
abline(h=150,col="grey30",lty=2)
lines(XX_250,
      YYcitron_250,
      type="s",
      lwd=2,
      col=2)
abline(h=0,col="grey30")
segments(0.5,100,0.5,150,lwd=2)
text(1.5,125,"50 (Hz)")
text(10,125,"Citronellal")
abline(h=50+yMax,col="grey30",lty=2)
abline(h=100+yMax,col="grey30",lty=2)
abline(h=150+yMax,col="grey30",lty=2)
lines(XX_250,
      YYterpi_250+yMax,
      type="s",
      lwd=2,
      col=2)
abline(h=yMax,col="grey30")
##segments(0.5,100+yMax,0.5,150+yMax,lwd=2)
##text(1,125+yMax,"50")
text(10,125+yMax,"Terpineol")
abline(h=50+2*yMax,col="grey30",lty=2)
abline(h=100+2*yMax,col="grey30",lty=2)
abline(h=150+2*yMax,col="grey30",lty=2)
lines(XX_250,
      YYmix_250+2*yMax,
      type="s",
      lwd=2,
      col=2)
##segments(0.5,100+2*yMax,0.5,150+2*yMax,lwd=2)
abline(h=2*yMax,col="grey30")
##text(1,125+2*yMax,"50")
text(10,125+2*yMax,"Mélange")
dev.off()


###################################################
### chunk number 5: psth 2
###################################################
png(file="mixPSTHs_comp.png",width=1000,height=1000)
par(mar=c(2,2,3,1),cex=2)
plot(XX_250,
     YYmix_250,
     type="n",
     axes=FALSE,
     xlab="",
     ylab="",
     xlim=c(5.5,8.5),
     ylim=c(0,yMax),
     main="Réponse au mélange")
rect(6,0,6.5,yMax,col="grey80",border=NA)
abline(h=50,col="grey30",lty=2)
abline(h=100,col="grey30",lty=2)
abline(h=150,col="grey30",lty=2)
abline(h=0,col="grey30")
segments(0.5,100,0.5,150,lwd=2)
text(1.5,125,"50 (Hz)")
lines(XX_10-diff(XX_10)[1]/2,
      YYmix_10,
      type="s",
      lwd=2,
      col=2)
lines(XX_50-diff(XX_50)[1]/2,
      YYmix_50,
      type="s",
      lwd=3,
      col=4)
lines(XX_250-diff(XX_250)[1]/2,
      YYmix_250,
      type="s",
      lwd=3,
      col=1)
legend(6.75,170,c("10 ms","50 ms","250 ms"),
       col=c(2,4,1),lwd=c(2,3,3),bty="n"
       )
dev.off()


###################################################
### chunk number 6: total count at 50 ms
###################################################
citronC <- sapply(allCitron,function(l) hist(l,seq(1,14,0.05),plot=FALSE)$counts)
citronC.hat <- apply(citronC,1,sum)


###################################################
### chunk number 7: bootstrap citronC
###################################################
set.seed(20061001,"Mersenne-Twister")
citronC.boot <- sapply(1:1000,
                       function(idx) {
                         myCol <- sample(1:20,20,replace=TRUE)
                         bootD <- citronC[,myCol]
                         apply(bootD,1,sum)
                         }
                       )
citronC.hat.var <- apply(citronC.boot,1,var)


###################################################
### chunk number 8: figure citronC
###################################################
png(file="VarCitronC_50.png",width=1000,height=1000)
par(mar=c(4,4,3,1),cex=2)
plot(citronC.hat,
     citronC.hat.var,
     xlab=expression(eta(x[i])),
     ylab=expression(sigma^2),
     main="Réponse au citronellal (bin de 50 ms)")
dev.off()


###################################################
### chunk number 9: sqrt of counts
###################################################
citronSC.hat <- sqrt(apply(citronC,1,sum))
citronSC.hat.var <- apply(sqrt(citronC.boot),1,var)


###################################################
### chunk number 10: figure citronSC
###################################################
png(file="VarCitronSC_50.png",width=1000,height=1000)
par(mar=c(4,4,3,1),cex=2)
plot(citronSC.hat,
     citronSC.hat.var,
     xlab=expression(sqrt(eta(x[i]))),
     ylab=expression(sigma^2),
     main="Réponse au citronellal (bin de 50 ms)")
abline(h=0.25,col=2,lwd=2)
dev.off()


###################################################
### chunk number 11: figure gogo
###################################################
png(file="CitronSCdist_50.png",width=1000,height=1000)
layout(matrix(1:4,nrow=2,byrow=TRUE))
par(mar=c(4,4,2,1),cex=2)
qqnorm(sqrt(citronC.boot[50,]),
       xlab="Quantiles théoriques",
       ylab="Quantiles de l'échantillon au bin 50",
       main="")
qqnorm(sqrt(citronC.boot[100,]),
       xlab="Quantiles théoriques",
       ylab="Quantiles de l'échantillon au bin 100",
       main="")
qqnorm(sqrt(citronC.boot[150,]),
       xlab="Quantiles théoriques",
       ylab="Quantiles de l'échantillon au bin 150",
       main="")
qqnorm(sqrt(citronC.boot[200,]),
       xlab="Quantiles théoriques",
       ylab="Quantiles de l'échantillon au bin 200",
       main="")
dev.off()


###################################################
### chunk number 12: gss fit1
###################################################
XX <- seq(1.025,13.975,0.05)
GF1 <- ssanova0(citronSC.hat ~ XX)
estGF1 <- predict(GF1,newdata=data.frame(XX=XX),se=TRUE)


###################################################
### chunk number 13: fig gss fit1
###################################################
png(file="GF1.png",width=800,height=800)
par(mar=c(4,4,3,1),cex=2)
plot(XX,citronSC.hat,type="n",
     xlab="Temps (s)",
     ylab=expression(sqrt(y)),
     main="ssanova0(citronSC.hat ~ XX)"
     )
polygon(c(XX,rev(XX)),
        c(estGF1$fit+1.96*estGF1$se.fit,
          rev(estGF1$fit-1.96*estGF1$se.fit)
          ),
        col="grey50",
        border=NA)
lines(XX,estGF1$fit,
      col=2,lwd=2)
points(XX,citronSC.hat)
dev.off()


###################################################
### chunk number 14: gss fit2
###################################################
s2 <- mean(citronSC.hat.var)
GF2 <- ssanova0(citronSC.hat ~ XX,method="u",varht=s2)
estGF2 <- predict(GF2,newdata=data.frame(XX=XX),se=TRUE)


###################################################
### chunk number 15: fig gss fit2
###################################################
png(file="GF2.png",width=800,height=800)
par(mar=c(4,4,3,1),cex=2)
plot(XX,citronSC.hat,type="n",
     xlab="Temps (s)",
     ylab=expression(sqrt(y)),
     main="ssanova0(citronSC.hat ~ XX,method=\"u\",varht=s2)"
     )
polygon(c(XX,rev(XX)),
        c(estGF2$fit+1.96*estGF2$se.fit,
          rev(estGF2$fit-1.96*estGF2$se.fit)
          ),
        col="grey50",
        border=NA)
lines(XX,estGF2$fit,
      col=2,lwd=2)
points(XX,citronSC.hat)
dev.off()


###################################################
### chunk number 16: smooth.spline fit
###################################################
SF <- smooth.spline(XX,citronSC.hat,all.knots=TRUE)


###################################################
### chunk number 17: fig smooth.spline fit
###################################################
png(file="SF.png",width=800,height=800)
par(mar=c(4,4,3,1),cex=2)
plot(XX,citronSC.hat,type="n",
     xlab="Temps (s)",
     ylab=expression(sqrt(y)),
     main="smooth.spline(XX,citronSC.hat,all.knots=TRUE)"
     )
polygon(c(XX,rev(XX)),
        c(estGF1$fit+1.96*estGF1$se.fit,
          rev(estGF1$fit-1.96*estGF1$se.fit)
          ),
        col="grey50",
        border=NA)
lines(SF,
      col=2,lwd=2)
points(XX,citronSC.hat)
dev.off()


###################################################
### chunk number 18: other fits
###################################################
citronGF <- ssanova0(sqrt(allCitronRG_50$counts) ~ XX)
citronPred <- predict(citronGF,newdata=data.frame(XX=XX),se=TRUE)
citronU <- (citronPred$fit + 1.96*citronPred$se.fit)^2
citronL <- (citronPred$fit - 1.96*citronPred$se.fit)^2
terpiGF <- ssanova0(sqrt(allTerpiRG_50$counts) ~ XX)
terpiPred <- predict(terpiGF,newdata=data.frame(XX=XX),se=TRUE)
terpiU <- (terpiPred$fit + 1.96*terpiPred$se.fit)^2
terpiL <- (terpiPred$fit - 1.96*terpiPred$se.fit)^2
mixGF <- ssanova0(sqrt(allMixRG_50$counts) ~ XX)
mixPred <- predict(mixGF,newdata=data.frame(XX=XX),se=TRUE)
mixU <- (mixPred$fit + 1.96*mixPred$se.fit)^2
mixL <- (mixPred$fit - 1.96*mixPred$se.fit)^2
yMax <- max(c(citronU,terpiU,mixU))


###################################################
### chunk number 19: final fig
###################################################
png(file="allSPSTHs_50.png",width=1000,height=1000)
par(mar=c(2,2,1,1),cex=2)
plot(XX,
     citronU,
     type="n",
     axes=FALSE,
     xlab="Temps (s)",
     ylab="Freq. (Hz)",
     xlim=c(1,14),
     ylim=c(0,yMax))
rect(6,0,6.5,yMax,col="grey80",border=NA)
abline(h=0,col=1,lwd=2)
abline(h=50,col="grey30",lty=2)
abline(h=100,col="grey30",lty=2)
abline(h=150,col="grey30",lty=2)
polygon(c(XX,rev(XX)),
        c(citronU,rev(citronL)),
        col=rgb(1,0,0,0.5),
        border=NA)
polygon(c(XX,rev(XX)),
        c(terpiU,rev(terpiL)),
        col=rgb(0,0,1,0.5),
        border=NA)
polygon(c(XX,rev(XX)),
        c(mixU,rev(mixL)),
        col=rgb(1,0,1,0.5),
        border=NA)
legend(8,0.9*yMax,c("Terpineol","Citronellal","Mélange"),col=c(rgb(0,0,1,0.5),rgb(1,0,0,0.5),rgb(1,0,1,0.5)),lwd=7)
dev.off()