#q()
#quit()


a <- 1
 a=1
1 -> a

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

#help(array)
#help(matrix)}

   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)

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

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

x<-read.csv("hendryetAl.csv")


length(x[,1])
length(x[1,]) 

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

   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))

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


   shapiro.test(x$Haldanes)
   ks.test(x$Haldanes,"pnorm", 
            mean=mean(x$Haldanes,na.rm=T), 
             sd=sd(x$Haldanes,na.rm=T))
mydata<-read.csv("fruits.csv")
   str(mydata)
#   help(aggregate)
   aggregate(mydata$Fruit, list(mydata$Manger), mean)
   aggregate(mydata$Fruit, list(Manger=mydata$Manger), mean)
   aggregate(x=mydata$Fruit, by=list(Manger=mydata$Manger), FUN=mean)
#   aggregate(mydata$Fruit, mydata$Manger, mean)

   n.fruit <- aggregate(mydata$Fruit, list(Manger=mydata$Manger), length)
   mean.fruit <- aggregate(mydata$Fruit, list(Manger=mydata$Manger), mean)
   sd.fruit <- aggregate(mydata$Fruit, list(Manger=mydata$Manger), 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"))

mydata<-read.csv("croissance.csv")
   str(mydata)
   aggregate(mydata$gain, list(mydata$supplement,mydata$diete), mean)
   mean.gain<-tapply(mydata$gain, list(mydata$supplement,mydata$diete), mean)

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

   barplot(mean.fruit,  
          xlab="Traitement",  
          ylab="Production de fruit",
          ylim=c(0,100))
   mids<-barplot(mean.fruit,  
          xlab="Traitement",  
          ylab="Production de fruit" , 
          ylim=c(0,100))
arrows(mids,mean.fruit+sd.fruit,  
              mids, mean.fruit - sd.fruit)
# help(arrows)

   mids<-barplot(mean.fruit,  
          xlab="Traitement",  
          ylab="Production de fruit"  ,
          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))

myData<-read.csv("croissance.csv")  
mean.gain<-tapply(myData$gain, 
           list(myData$supplement,myData$diete), 
           mean) 
sd.gain<-tapply(myData$gain, 
           list(myData$supplement,myData$diete), 
           sd) 
n.gain<-tapply(myData$gain, 
           list(myData$supplement,myData$diete), 
           length) 
 barplot(mean.gain) 

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

   mydata<-read.csv("fruits.csv")
   plot(mydata$Racine, mydata$Fruit)
   plot(mydata$Racine, mydata$Fruit,  
   xlab="Largeur des racines",  
   ylab="Production de fruit")  
   plot(mydata$Racine, mydata$Fruit,  
   xlab="Largeur des racines",  
   ylab="Production de fruit", 
   pch=21,bg="grey",cex=2.0) 
   clr<-ifelse(mydata$Manger == "Mange", "Green","Blue")
   plot(mydata$Racine, mydata$Fruit,  
   xlab="Largeur des racines",  
   ylab="Production de fruit" ,
   pch=21,bg=clr,cex=2.0) 
legend("topleft",  
   legend=c("Mange", "Intact"), pch=21,  
   pt.bg=c("Green", "Blue"), pt.cex=2.0)

#   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))

   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))

#   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))

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

   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]

   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]


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

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


    my.lm<-lm(x$Haldanes~ x$GLength)
    summary(lm(x$Haldanes~ x$GLength))
    summary(my.lm)
    attributes(my.lm)
    my.lm$residuals
    my.lm$coefficients

    my.summ.lm<-summary(my.lm)
    attributes(my.summ.lm)
    my.summ.lm$r.squared
    my.summ.lm$fstatistic
plot(x$Haldanes~ x$GLength)
   plot(x$GLength,x$Haldanes)

   lm(x$HaldanesAbs~ x$GLength)
   summary(lm(x$HaldanesAbs~ x$GLength))

plot(x$HaldanesAbs~ x$GLength)
   help.search("regression")
# install.packages()
library(car)
#help(reg.line)

rg<-x$HaldanesAbs~ x$GLength
plot(rg,xlab="Glength",ylab="Absolute Haldanes")
reg.line(lm(rg))

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

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

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

# ANCOVA

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

ResFull.lm<-lm(ResFull)
ResCD.lm<-lm(ResCD)
ResE.lm<-lm(ResE)
plot(ResFull.lm) 
plot(ResCD.lm)
plot(ResE.lm)


summary(ResFull.lm)
summary(ResCD.lm)
summary(ResE.lm)
anova(ResFull.lm, ResE.lm) 
anova(ResFull.lm, ResCD.lm, ResE.lm) 

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


# Fibonacci

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 
} 

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



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]) 
} 

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) 
}

chaos(1000,3.5,0.5) 
c<-chaos(1000,3.5,0.5) 
plot(c) 
plot(chaos(1000,3.5,0.5)) 

mySum <- function(N){
   return(sum(seq(1,N)))
}

mySum2 <- function(N){
   return((N*(N+1))/2)
}


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

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



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

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

#####################
M1 <- matrix(c(0,0.7,0,1.29,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)
#####################

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

M1 <- matrix(c(0,0.7,0,1.29,0,0.2,0.702,0,0), ncol=3, nrow=3)
inv(M1) # (fait partie du paquet {\it car})  
M1 %*% inv(M1)
det(M1)
diag(M1)

# install.packages("deSolve")
library(deSolve)
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)

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)) 

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="-")