# program to calculate response to directional mass selection
# input parameters

p=0.01    # selected proportion
varam=0.1 # genetic variancve in phenotype
varav=0.05 # genetic variance Ve
varp=1.0   # phenotypic variance
rgamav=0.0   # genetic correlation between am and av

# derived parameters
Ve=varp-varam

x<-qnorm(p,mean=0,sd=1,lower.tail=FALSE,log.p=FALSE) # truncation point normal distribution
d<-dnorm(x,mean=0,sd=1,log=FALSE) #value of the normal density function at truncation point x
i<-d/p # selection intensity
seldif_var=i*x

# b-values without correlation
b1<-varam/varp # heritability
b2<-varav/(2*(varp^2)+3*varav) #heritability Ve

# response in mean and in variance  (question 1)
delta_Am<-i*b1*sqrt(varp)
delta_Av<-seldif_var*b2*varp

delta_Am
delta_Av

# Ve and heritability in next generation (question 2)
Ve_1=Ve+delta_Av
heritability<-varam/(varam+Ve_1)	#assumption is that varam stays the same


# matrix of P and P2, if rgamav is nonzero  (question 4)

P<-array(0,c(2,2))
P[1,1]<-varp
P[1,2]<-3*rgamav*sqrt(varav*varam)
P[2,1]<-P[1,2]
P[2,2]<-2*(varp^2)+3*varav

invP<-solve(P)

G<-array(0,c(2,2))
G[1,1]<-varam
G[1,2]<-rgamav*sqrt(varav*varam)
G[2,1]<-rgamav*sqrt(varav*varam)
G[2,2]<-varav

b<-array(0,c(2,2))

b<-invP%*%G

response<-array(0,c(2,2))

x<-array(0,c(2))
x[1]<-i
x[2]<-seldif_var

response<-t(b)%*%x     # genetic response is the regression of the breeding values on the selection differentials

# question 3 varying rgamav
rgamav=-1
results_rg=data.frame()
for (ii in 1:19) {
rgamav=rgamav+0.1

P<-array(0,c(2,2))
P[1,1]<-varp
P[1,2]<-3*rgamav*sqrt(varav*varam)
P[2,1]<-P[1,2]
P[2,2]<-2*(varp^2)+3*varav

invP<-solve(P)

G<-array(0,c(2,2))
G[1,1]<-varam
G[1,2]<-rgamav*sqrt(varav*varam)
G[2,1]<-rgamav*sqrt(varav*varam)
G[2,2]<-varav

b<-array(0,c(2,2))

b<-invP%*%G

response<-array(0,c(2,2))

x<-array(0,c(2))
x[1]<-i
x[2]<-seldif_var

response<-t(b)%*%x     # genetic response is the regression of the breeding values on the selection differentials

results_rg <- rbind(results_rg,data.frame(rgamav=rgamav,response_am=response[1],response_av=response[2] ))
}
directory=getwd()
write("OUTPUT: results_rg written to file","results_rg.txt",append=FALSE)
write(c(directory,date()),"results_rg.txt",append=TRUE)
write.table(results_rg,"results_rg.txt",row.names=FALSE,quote=FALSE,append=TRUE)
print (paste("File results_rg.txt written to directory ",directory))