# slide 5
#q()
#quit()

#slide 6
#setwd("~/Desktop/<yourName>")

#slide 8
a <- 1
a=1
1 -> a
b <-c(1,2,3)

#slide 9
a <- array(NA,dim=10)
a[4] <- 5
b <- matrix(NA,ncol=10,nrow=30)
b[30,3]<-1

#slide 10
#help(array)
#help(matrix)}

#slide 11
a<-array(seq(1,10,2))
   a<-array(rnorm(10,mean=15,sd=3))
   b<-matrix(seq(1,20),ncol=2,nrow=10)
   b<-matrix(runif(21),ncol=3,nrow=7)

#slide 12 
#####################################
c<-c(rnorm(5000,mean=10,sd=4), rnorm(5000,mean=15,sd=4))
#####################################

#slide 13
a+a
a+5
1+b
5*a
a*a
#a*b
b*b
a-a
a-5
a/2

#slide 18
x<-read.csv("hendryEtAl.csv")

#slide 20
length(x[,1])
length(x[1,]) 

#slide 21
x$Years
x[,18]
#   x$years

#slide 23
   x<-read.csv("hendryEtAl.csv")
   hist(x$Haldanes)
   hist(x$Haldanes,breaks=100)
#   hist(..., breaks=seq(from=-1.2,to=0.8,by=0.1))

#slide 24
#   hist(...,col=2)
#   hist(...,border=3)

#slide 25
shapiro.test(x$Haldanes)
   ks.test(x$Haldanes,"pnorm", 
            mean=mean(x$Haldanes,na.rm=T), 
             sd=sd(x$Haldanes,na.rm=T))

#slide 26
mydata<-read.csv("fruits.csv")
   str(mydata)
#   help(aggregate)
   aggregate(mydata$Fruit, list(mydata$State), mean)
   aggregate(mydata$Fruit, list(State=mydata$State), mean)
   aggregate(x=mydata$Fruit, by=list(State=mydata$State), FUN=mean)
#   aggregate(mydata$Fruit, mydata$State, mean)

#slide 27
   n.fruit <- aggregate(mydata$Fruit, list(State=mydata$State), length)
   mean.fruit <- aggregate(mydata$Fruit, list(State=mydata$State), mean)
   sd.fruit <- aggregate(mydata$Fruit, list(State=mydata$State), sd)
summary.table <- cbind(n.fruit[,2], mean.fruit[,2], sd.fruit[,2])
summary.table 
dimnames(summary.table) <- list(n.fruit[,1], c("n","mean","SD"))

#slide 28
mydata<-read.csv("gain.csv")
   str(mydata)
   aggregate(mydata$growth, list(mydata$experiment,mydata$food), mean)
   mean.growth<-tapply(mydata$growth, list(mydata$experiment,mydata$food), mean)

#slide 29
mydata<-read.csv("fruits.csv")
   n.fruit <- tapply(mydata$Fruit, list(State=mydata$State), length)
   mean.fruit <- tapply(mydata$Fruit, list(State=mydata$State), mean)
   sd.fruit <- tapply(mydata$Fruit, list(State=mydata$State), sd)
   barplot(mean.fruit)
#    help(barplot)}

#slide 30
barplot(mean.fruit,  
          xlab="Treatment",  
          ylab="Fruit production",
          ylim=c(0,100))
   mids<-barplot(mean.fruit,  
          xlab="Treatment",  
          ylab="Fruit production" , 
          ylim=c(0,100))
   arrows(mids,mean.fruit+sd.fruit,  
              mids, mean.fruit - sd.fruit)
   # help(arrows)

#slide 31
mids<-barplot(mean.fruit,  
          xlab="Treatment",  
          ylab="Fruit production"  ,
          ylim=c(0,100))
   arrows(mids,mean.fruit+sd.fruit,  
              mids, mean.fruit - sd.fruit,  
              angle=90, code=3)
   text(mids,5, paste("N = ", n.fruit))

#slide 32
myData<-read.csv("gain.csv")  
   mean.growth<-tapply(myData$growth, 
           list(myData$experiment,myData$food), 
           mean) 
   sd.growth<-tapply(myData$growth, 
           list(myData$experiment,myData$food), 
           sd) 
   n.growth<-tapply(myData$growth, 
           list(myData$experiment,myData$food), 
           length) 
   barplot(mean.growth) 

#slide 33
mids<- barplot(mean.growth,beside=T, 
       xlab="Food type", 
       ylab="Gain", 
       ylim=c(0,35), 
       col=grey(c(0,0.3,0.6,1))) 
   arrows(mids, mean.growth+sd.growth, 
       mids, mean.growth-sd.growth, 
       angle=90, code=3, length=0.1) 
   text(mids, 2, paste(n.growth), col=c("white", rep("black",3))) 
   legend("topleft", legend=rownames(mean.growth), fill=grey(c(0,0.3,0.6,1))) 

# slide 34
mydata<-read.csv("fruits.csv")
   plot(mydata$Root, mydata$Fruit)
   plot(mydata$Root, mydata$Fruit,  
   xlab="Root size",  
   ylab="Fruit production")  
   plot(mydata$Root, mydata$Fruit,  
   xlab="Root size",  
   ylab="Fruit production", 
   pch=21,bg="grey",cex=2.0) 

# slide 35
clr<-ifelse(mydata$State == "Eaten", "Green","Blue")
   plot(mydata$Root, mydata$Fruit,  
   xlab="Root size",  
   ylab="Fruit production" ,
   pch=21,bg=clr,cex=2.0) 
legend("topleft",  
   legend=c("Eaten", "Intact"), pch=21,  
   pt.bg=c("Green", "Blue"), pt.cex=2.0)

# slide 36
#   plot(..., main="Titre")
#   plot(..., xlab="nomX",ylab="nomY")
#   plot(..., cex=2.0)
#   plot(..., cex.lab=2.0)
#   plot(..., cex.axis=2.0)
#   plot(..., xlim=c(0,100),ylim=c(0,2))

# slide 37
   points(x=5,y=80)
#   points(...,pch=2)
#   points(..., col=2)
#   points(..., cex=2.0)
#   points(x=c(1,2,3),y=c(1,2,3))

# slide 38
#   lines(x=c(x1,x2),y=c(y1,y2))
   lines(x=c(5,10),y=c(100,20))
#   lines(...,lty=2)
#   lines(..., col=2)
#   lines(..., lwd=2.0)
#   lines(x=c(x1,x2,x3),y=c(y1,y2,y3))

# slide 40
plot(1:10, xlab=expression(paste( "Text, Greek ", lambda[subscript]," ",
and^{Power})), ylab="Simple text", main="Cristian you better love this")

# slide 41
#   jpeg("fileName.jpg")
#   plot(...)
#   dev.off()
#   help(jpeg)

#slide 43
plot(x$Years,x$Haldanes)
   plot(x$Years,abs(x$Haldanes),
      xlab="Years", ylab="Haldanes",
      main="Absolute values")
   plot(x$Years,abs(x$Haldanes),
      xlim=c(0,100), ylim=c(0,2))

#slide 44
   x<-read.csv("hendryEtAl.csv")
   x$Years
   x$Years>100
   x$Years==111
   x$Years!=124
   cond<-x$Years>100
   x$Years[cond]
   x$Haldanes[cond]
   x$Haldanes[x$Years>100]

#slide 46
   cond<-x$Years>100 & x$Years<113
   cond<-x$Years<100 | x$Years>113
   cond<-x$Years>100 & x$Haldanes>=0
   x$Years[cond]
   x$Haldanes[cond]

#slide 47
y<-array(1,length(x$Haldanes))
   y[x$Haldanes<0]<-2
   plot(x$Years,abs(x$Haldanes),col=y,
      xlab="Years",ylab="Haldanes")
   legend("topright",
      c("positive","negative"), col=c(1,2),
      pch=c(1,1))


#slide 48
#   save(x,y,z,file="saveXYZ.RData")
#   save.image(file="workspace.RData")

#slide 49
#   load("saveXYZ.RData")
#   load("workspace.RData")
   ls()

#slide 50
lm(x$Haldanes~x$GLength)
   my.lm<-lm(x$Haldanes~ x$GLength)
    summary(lm(x$Haldanes~ x$GLength))
    summary(my.lm)

#slide 51
attributes(my.lm)
    my.lm$residuals
    my.lm$coefficients

#slide 52
    my.summ.lm<-summary(my.lm)
    attributes(my.summ.lm)
    my.summ.lm$r.squared
    my.summ.lm$fstatistic

#slide 53
plot(x$Haldanes~ x$GLength)
   plot(x$GLength,x$Haldanes)

#slide 54
lm(x$HaldanesAbs~ x$GLength)
   summary(lm(x$HaldanesAbs~ x$GLength))

#slide 55
plot(x$HaldanesAbs~ x$GLength)
#   help.search("regression")

#slide 56
# install.packages()

#slide 57
library(car)
#  help(regLine)
   rg<-x$HaldanesAbs~ x$GLength
   plot(rg,xlab="Glength",ylab="Absolute Haldanes")
   regLine(lm(rg))

#slide 60
m1<- x$HaldanesAbs~ x$GLength+x$Years
   m2<- x$HaldanesAbs~ x$GLength:x$Years
   m3<- x$HaldanesAbs~ x$GLength*x$Years
   m3<- x$HaldanesAbs~ x$GLength + x$Years + x$GLength:x$Years

#slide 62
z<-read.table("anova.txt")
names(z) <- c("response", "category", "replicat", "coVar")
attach(z)
response
detach(z)
#response
attach(z)

# slide 63
category<-factor(category)

# slide 64
mod1<-response~ category 
   boxplot(mod1)
   mod1.lm<-lm(mod1)
   plot(mod1.lm)
   summary(mod1.lm)
   anova(mod1.lm)

# ANCOVA

#slide 65
as.factor(category) 
   plot(response ~  coVar,pch=as.numeric(category))
   ResE<-response ~  coVar 
   ResCD<-response ~  category+ coVar 
   ResFull<-response ~  category+ coVar + category:coVar

#slide 66
ResE.lm<-lm(ResE)
ResCD.lm<-lm(ResCD)
ResFull.lm<-lm(ResFull)
plot(ResE.lm)
plot(ResCD.lm)
plot(ResFull.lm) 

#slide 67
summary(ResE.lm)
summary(ResCD.lm)
summary(ResFull.lm)
anova(ResE.lm,ResFull.lm, ) 
anova(ResE.lm,ResCD.lm,ResFull.lm) 

#slide 68
tapply(response, category, var, na.rm=TRUE) 
tapply(response, category, function(x) shapiro.test(x)) 


# Fibonacci

# slide 71
x<- array(1,dim=10) 
for (i in seq(3,length(x))) {  
   x[i]<-x[i-1]+x[i-2] 
} 

x<- array(1,dim=10) 
i<-3
while (i <= length(x)) {  
   x[i]<-x[i-1]+x[i-2] 
   i<-i+1 
} 

#slide 74
x<-25
if (x>0) { 
   x<-x+1
} else { 
   x<-x-1 
} 

#slide 76
x<-array(dim=1000)
x[1]=0.5
r<-3.7
for (t in 2:length(x)){ 
   x[t]<-r*x[t-1]*(1-x[t-1]) 
} 

#slide 79
chaos <- function(time,r,popS){ 
   x<-array(dim=time)
   x[1]=popS
   for (t in 2:length(x)){ 
      x[t]<-r*x[t-1]*(1-x[t-1]) 
   } 
   return(x) 
}

#slide 80
chaos(1000,3.5,0.5) 
c<-chaos(1000,3.5,0.5) 
plot(c) 
plot(chaos(1000,3.5,0.5)) 

#slide 82
V1 <- as.vector(seq(1,10)) 
V2 <- as.vector(rnorm(10)) 
V3<- V1 * V2 
s1<- t(V1) %*% V2 
M1<- V1 %*% t(V2) 

#slide 83
M3 <- matrix(c(1,2,3,4,5,6),ncol=3,nrow=2)
   M4 <- matrix(c(1,2,3,4,5,6),ncol=2,nrow=3)

#slide 84
M1 <- matrix(c(1,2,3,4,5,6), ncol=3, nrow=2)
   V1 <- as.vector(c(7,8,9))
   M1 %*% V1
#    V1 %*% M1

#slide 85
M1 <- matrix(c(1,2,3,4,5,6), ncol=3, nrow=2)
   M2 <- matrix(c(7,8,9,10,11,12), ncol=2, nrow=3)
   M1 %*% M2
#  M2 %*% M1

#slide 86
#####################
M1 <- matrix(c(0,0.7,0,0.63,0,0.2,0.702,0,0), ncol=3, nrow=3)
V0 <- as.vector(c(1,0,0))
for (i in 1:100){
   V1<-M1 %*% V0
   V0<-V1
}
prop<-V0/sum(V0)
#####################

#slide 87
M1 <- matrix(c(0,0.7,0,0.63,0,0.2,0.702,0,0), ncol=3, nrow=3)
eigen(M1)
################
prop2<-eigen(M1)$vectors[,1]/sum(eigen(M1)$vectors[,1])
################


#slide 88
M1 <- matrix(c(0,0.7,0,0.63,0,0.2,0.702,0,0), ncol=3, nrow=3)
solve(M1)   
M1 %*% solve(M1)
det(M1)
diag(M1)

#slide 89
# install.packages("deSolve")
library(deSolve)

#slide 90
parameters<-c(r=1.5,K=100) 
state<-c(X=80) 
logistic <- function(t, state, parameters){ 
   with(as.list(c(state, parameters)), 
   { 
      dX=r*X*(K-X)/K 
      return(list(dX)) 
   }) 
}
times <- seq(1,100,by=0.1) 
out <- as.data.frame(ode(y=state, times=times, func=logistic, parms=parameters))
plot(out$X)

#slide 92
parameters<-c(b=1.0, d=0.2, beta=0.1, gamma=0.7, D=0.3) 
state<-c(X=15, Y=4, Z=0) 
SIR<- function(t,state,parameters){ 
      with(as.list(c(state,parameters)),{ 
      dX = b*(d*X+D*Y+d*Z) - beta*X*Y - d*X 
      dY = beta*X*Y - (D+gamma)*Y 
      dZ = gamma*Y -d*Z 
      return(list(c(dX,dY,dZ))) 
  }) 
} 
times<- seq(1,100,by=1) 
out<-as.data.frame(ode(y=state,times=times,  
  func=SIR,parms=parameters)) 

#slide 93
par(mfrow=c(2,2), oma=c(0,0,3,0)) 
plot(times,out$X, type="l",      main="Susceptible", xlab="time", ylab="-") 
plot(times,out$Y, type="l",      main="Infecte", xlab="time", ylab="-") 
plot(times,out$Z, type="l",      main="Gueri", xlab="time", ylab="-") 

###########################################################
# Solution slide 94
parameters<-c(r=2.0, a=0.1, p=0.01, m=0.1)
state<-c(X=200, Y=150)

lvM<- function(t,state,parameters){
   with(as.list(c(state,parameters)),{
# Equations
      dX = r*X - p*X*Y
      dY =  a*p*X*Y -m*Y
      return(list(c(dX,dY)))
   })
}

times<- seq(1,100,by=0.01)
out<-as.data.frame(ode(y=state,times=times,func=lvM,parms=parameters))

pdf("lvGraph.pdf")
plot(out$X, out$Y, type="l", xlab="H", ylab="P")

x<-with(as.list(parameters), {
peq<-r/p
heq<-m/(a*p)
return(c(peq,heq))
})

peq<-x[1]
heq<-x[2]


lines(c(heq,heq),c(-10000,10000),col=2)
lines(c(-10000,10000),c(peq,peq),col=3)
legend("topright", c("Predator","Prey"),col=c(2,3),lty=1)
###########################################################