#simulations to check model works if all assumptions are met
#show the estimates for varying values of p,s in NBS and BS
#and varying intercepts for b 


#set 1: high p (0.5), high s (0.7), high sB (0.8), high pB (0.6), strong effect on b(slope = -1)

#set 2: high p (0.5), high s (0.7), low sB (0.4), low pB (0.3),  strong effect on b(slope = -1)

#set 3: low p (0.2), low s (0.3), high sB (0.8), high pB (0.6), strong effect on b(slope = -1)

#set 4: low p (0.2), low s (0.3), low sB (0.4), low pB (0.3),  strong effect on b(slope = -1)


#set 5: high p (0.5), high s (0.7), high sB (0.8), high pB (0.6), mild effect on b(slope = 1)

#set 6: high p (0.5), high s (0.7), low sB (0.4), low pB (0.3),  mild effect on b(slope = 1)

#set 7: low p (0.2), low s (0.3), high sB (0.8), high pB (0.6), mild effect on b(slope = 1)

#set 8: low p (0.2), low s (0.3), low sB (0.4), low pB (0.3),  mild effect on b(slope = 1)


source("AllFunctions.R")


#------------------------------------------------------------------------------------------------#
# beta(t)phi(c)p(c)s(c)psi(lb)xi(lb)eta(lb)


mLL.betat.phic.pc.sc.psilb.xilb.etalb <- function(invect)
   {
   	this.NM    <- exp(invect[1])
   
   	this.NU    <- exp(invect[2])

   	this.beta <- c(lgp(invect[3:(K_NB + 1)]))

	this.phi <- matrix(expo(invect[(K_NB + 2)]), K_NB-1, K_NB-1, byrow = TRUE);
	this.phi[lower.tri(this.phi)] <- 0

    this.p <- matrix(expo(invect[K_NB + 3]), nrow = K_NB, ncol = K_NB, byrow = TRUE);
   	this.p[ ,what[1:K_NB] == 2] <- 0
    this.p[lower.tri(this.p)] <- 0

    this.p1 <- matrix(expo(invect[K_NB + 3]), nrow = K_NB, ncol = K_NB, byrow = TRUE);
   	this.p1[ ,what[1:K_NB] == 2] <- 0
    this.p1[lower.tri(this.p1)] <- 0

    this.s <- matrix(expo(invect[K_NB + 4]), nrow = K_NB, ncol = K_NB, byrow = TRUE);
   	this.s[ ,what[1:K_NB] == 1] <- 0
    this.s[lower.tri(this.s)] <- 0

    this.s1 <- matrix(expo(invect[K_NB + 4]), nrow = K_NB, ncol = K_NB, byrow = TRUE);
   	this.s1[ ,what[1:K_NB] == 1] <- 0
    this.s1[lower.tri(this.s1)] <- 0

	this.psi <- expo(invect[K_NB + 5] + invect[K_NB + 6]*scale(1:K_NB)[,1])

	this.xi <- expo(invect[K_NB + 7] + invect[K_NB + 8]*scale(1:K_NB)[,1])

	this.pb <- expo(invect[K_NB + 9])
	
	this.sb <- expo(invect[K_NB + 10])
	
	this.c <- expo(invect[K_NB + 11] + invect[K_NB + 12]*scale(1:K_NB)[,1])


LL.WithBreeding.c(this.NM, this.NU, this.beta, this.phi, this.p, this.p1, this.s, this.s1,
this.psi, this.xi, this.pb, this.sb, this.c)
   }

fit.betat.phic.pc.sc.psilb.xilb.etalb <- function(pars.start){

    this.fit <-
            try(optim(par = pars.start, fn = mLL.betat.phic.pc.sc.psilb.xilb.etalb, 
            lower = c(minLogNM, minLogNU, rep(-Inf, npar - 2)),
            method = "L-BFGS-B", hessian = TRUE, control = list(maxit = 10000)), TRUE)

minvalue  <- this.fit$value 

hess <- this.fit$hessian

npar <- length(this.fit$par) 

AIC <- 2*minvalue + 2*npar

pars <- this.fit$par

 	NM.est    <- exp(pars[1])
 
   	NU.est    <- exp(pars[2])

   	beta.est <- c(lgp(pars[3:(K_NB + 1)]))

	phi.est <- matrix(expo(pars[(K_NB + 2)]), K_NB - 1, K_NB - 1, byrow = TRUE)
	phi.est[lower.tri(phi.est)] <- 0

    p.est <- matrix(expo(pars[K_NB + 3]), nrow = K_NB, ncol = K_NB, byrow = TRUE)
   	p.est[ ,what[1:K_NB] == 2] <- 0
    p.est[lower.tri(p.est)] <- 0

    p1.est <- matrix(expo(pars[K_NB + 3]), nrow = K_NB, ncol = K_NB, byrow = TRUE)
   	p1.est[ ,what[1:K_NB] == 2] <- 0
    p1.est[lower.tri(p1.est)] <- 0

    s.est <- matrix(expo(pars[K_NB + 4]), nrow = K_NB, ncol = K_NB, byrow = TRUE)
   	s.est[ ,what[1:K_NB] == 1] <- 0
    s.est[lower.tri(s.est)] <- 0
                      
   	s1.est <- matrix(expo(pars[K_NB + 4]), nrow = K_NB, ncol = K_NB, byrow = TRUE)
   	s1.est[ ,what[1:K_NB] == 1] <- 0
    s1.est[lower.tri(s1.est)] <- 0

	psi.est <- expo(pars[K_NB + 5] + pars[K_NB + 6]*scale(1:K_NB)[,1])

	xi.est <- expo(pars[K_NB + 7] + pars[K_NB + 8]*scale(1:K_NB)[,1])

	pb.est <- expo(pars[K_NB + 9])
	
	sb.est <- expo(pars[K_NB + 10])
	
	c.est <- expo(pars[K_NB + 11] + pars[K_NB + 12]*scale(1:K_NB)[,1])


list("minvalue" = minvalue, "hessian" = hess, "AIC" = AIC, "npar" = npar,
"NM.est" = round(NM.est, 3), "NU.est" = round(NU.est, 3), "beta.est" = round(beta.est,3),
"phi.est" = round(phi.est, 3), "p.est" = round(p.est, 3), "p1.est" = round(p1.est, 3),
"s.est" = round(s.est, 3), "s1.est" = round(s1.est, 3), "psi.est" = round(psi.est, 3),
"xi.est" = round(xi.est, 3), "pb.est" = round(pb.est, 3), "sb.est" = round(sb.est, 3),
"c.est" = round(c.est, 3), "pars" = pars) }


K_NB <- 10
K_B <- 4
K <- K_NB + K_B
N_Mtr <- 800
N_Utr <- 1000

what <- c(1, 1, 2, 2, 3, 3, 2, 2, 1, 1, 0, 0, 0, 0)

beta.tr <- c(0.1, 0.04, 0.18, 0.12, 0.01, 0.2, 0.3, 0.04, 0.01, 0, 0, 0, 0, 0)
length(beta.tr)
sum(beta.tr);

plot(beta.tr, type ="l", ylim=c(0,1))

phi.tr <- matrix(c(rep(0.8, K_NB-1), rep(1, 4)), K_NB-1 + 4, K_NB-1 + 4, byrow = T)

p1.tr <- matrix(0.2, K_NB, K_NB)

p2.tr <- matrix(0.2, K_NB, K_NB)

p1.tr[,what[1:K_NB]==2] <- p2.tr[,what[1:K_NB]==2] <- 0

s.tr <- matrix(0.3, K_NB, K_NB)

s1.tr <- matrix(0.3, K_NB, K_NB)

s.tr[,what[1:K_NB]==1] <- s1.tr[,what[1:K_NB]==1] <- 0

psi.tr <- c(expo(-1 - 4*scale(1:K_NB)[,1]))

xi.tr <- c(expo(-1 - 2*scale(1:K_NB)[,1]))

sb.tr <- 0.8

pb.tr <- 0.6

c.tr <- c(expo(-1 - 3*scale(1:K_NB)[,1]))

NMest3 <- c()
NUest3 <- c()
slpsiest3 <- c()
slksiest3 <- c()
sletaest3 <- c()

psiest3 <- matrix(NA, nrow=100, ncol = K_NB)
xiest3 <- matrix(NA, nrow=100, ncol = K_NB)
etaest3 <- matrix(NA, nrow=100, ncol = K_NB)

for(k in 1:100){

print(k)

	DATA<- data.simulation(K_NB,N_Mtr,N_Utr,beta.tr,phi.tr,p1.tr,p2.tr,s.tr,s1.tr,psi.tr,xi.tr,pb.tr,sb.tr,c.tr,what)

	x.matM <- DATA$UCH_M[,1:(K_NB+4)]
	y.vectM <- DATA$UCH_M[,K_NB+5]
	DDM <- sum(y.vectM)

	x.matU <- DATA$UCH_U[,1:(K_NB+4)]
	y.vectU <- DATA$UCH_U[,K_NB+5]
	DDU <- sum(y.vectU)


pars.start <- c(log(600), log(600), plg(rep(1/K_NB, K_NB)), logit(runif(3, 0.01, 0.99)), rep(0, 8))
npar <- length(pars.start)
minLogNM <- log(DDM); minLogNU <- log(DDU); 

fit <- try(fit.betat.phic.pc.sc.psilb.xilb.etalb(pars.start), T)

if(is.list(fit)){
NMest3 <- c(NMest3, exp(fit$pars[1]))
NUest3 <- c(NUest3, exp(fit$pars[2]))
slpsiest3 <- c(slpsiest3, fit$pars[16])
slksiest3 <- c(slksiest3, fit$pars[18])
sletaest3 <- c(sletaest3, fit$pars[22])


psiest3[k,] <- fit$psi.est
xiest3[k,] <- fit$xi.est
etaest3[k,] <- fit$c.est

					}
}




par(mar = c(5.1, 5.1, 1, 1))
boxplot(psiest1, xlab = "b", ylab = expression(psi[b]),
cex.lab=1.4, cex.axis=1.1,axes=F)
points(c(expo(-1 - 4*scale(1:K_NB)[,1])),pch=19)
axis(2,at=seq(0,1,0.1),tcl=0.5)
axis(1,at=seq(1,10,1),tcl=0.5)

par(mar = c(5.1, 5.1, 1, 1))
boxplot(xiest1, xlab = "b", ylab = expression(xi[b]),
cex.lab=1.4, cex.axis=1.1,axes=F)
points(c(expo(-1 - 2*scale(1:K_NB)[,1])),pch=19)
axis(2,at=seq(0,1,0.1),tcl=0.5)
axis(1,at=seq(1,10,1),tcl=0.5)

par(mar = c(5.1, 5.1, 1, 1))
boxplot(etaest1, xlab = "b", ylab = expression(eta[b]),
cex.lab=1.4, cex.axis=1.1,axes=F)
points(c(expo(-1 - 3*scale(1:K_NB)[,1])),pch=19)
axis(2,at=seq(0,1,0.1),tcl=0.5)
axis(1,at=seq(1,10,1),tcl=0.5)

par(mar = c(5.1, 5.1, 1, 1))
boxplot(psiest4, xlab = "b", ylab = expression(psi[b]),
cex.lab=1.4, cex.axis=1.1,axes=F)
points(c(expo(-1 - 4*scale(1:K_NB)[,1])),pch=19)
axis(2,at=seq(0,1,0.1),tcl=0.5)
axis(1,at=seq(1,10,1),tcl=0.5)

boxplot(xiest4, xlab = "b", ylab = expression(xi[b]),
cex.lab=1.4, cex.axis=1.1,axes=F)
points(c(expo(-1 - 2*scale(1:K_NB)[,1])),pch=19)
axis(2,at=seq(0,1,0.1),tcl=0.5)
axis(1,at=seq(1,10,1),tcl=0.5)

par(mar = c(5.1, 5.1, 1, 1))
boxplot(etaest4, xlab = "b", ylab = expression(eta[b]),
cex.lab=1.4, cex.axis=1.1,axes=F)
points(c(expo(-1 - 3*scale(1:K_NB)[,1])),pch=19)
axis(2,at=seq(0,1,0.1),tcl=0.5)
axis(1,at=seq(1,10,1),tcl=0.5)