library(stockassessment)
setwd("run")
load("data.RData")
conf<-loadConf(dat,"../conf/model.cfg", patch=TRUE)
par<-defpar(dat,conf)
fit<-sam.fit(dat,conf,par)
if(fit$opt$convergence!=0) stop("Model did not converge.")


runwithnewm<-function(fit, M){
  dat<-fit$data
  dat$natMor<-M
  fit<-sam.fit(dat,conf,par)  
  fit
}
numberEaten<-function(a,y, fit, pcc){
    Ntable <- ntable(fit)
    Ftab <- faytable(fit)
    CC <- pcc[[as.character(y)]][-c(1:5),-c(1:5)]    
    C1 <- CC[,((1:ncol(CC))%%2)==1]
    C2 <- CC[,((1:ncol(CC))%%2)==0]
    idx <- c(1:ncol(C1), rep(ncol(C1),4))
    C1 <- C1[,idx]
    C2 <- C2[,idx]
    thisN <- Ntable[rownames(Ntable)==y,]
    thisF <- Ftab[rownames(Ftab)==y,]
    mo <- fit$data$propMat
    thisO <-mo[rownames(mo)==y,]
    thisM <- fit$data$natMor[rownames(fit$data$natMor)==y,]
    ret <- sum(C1[a-2,]*thisN*exp(-(thisF+thisM)*.25)*(1-thisO))+sum(C2[a-2,]*thisN*exp(-(thisF+thisM)*.75))
    return(ret)  
}



result<-list()
result[[1]]<-fit

maxyear<-max(fit$data$years)

for(run in 2:5){
  yrrange<-1984:maxyear
  fit<-result[[run-1]]  
  Ntable <- ntable(fit)
  Ftab <-faytable(fit)
  yridx<-rownames(Ntable)%in%yrrange  
  oldN7<-Ntable[yridx,1:5]
  checkN<-oldN7
  oldN7[1:(nrow(oldN7)-1),1:4]<-NA
  oldM<-fit$data$natMor[yridx,1:5]
  oldF<-Ftab[yridx,1:5]
  Z<-fit$data$natMor+Ftab
  modelC<-Ftab/Z*(1-exp(-Z))*Ntable
   modelC<-modelC[rownames(modelC)%in%yrrange,1:4] 
  #oldC<-cn[rownames(cn)%in%yrrange,1:4] 
  con<-modelC*0
  con<-con[-nrow(con),]
  for(y in 1:nrow(con)){
    for(a in 1:ncol(con)){
      cat(y," ",a," \n")  
      con[y,a]<-numberEaten(as.integer(colnames(con))[a],as.integer(rownames(con))[y], fit, pcc)   
    }
  }
  N7<-oldN7
  minAgem1<-3-1
  minYearm1<-1984-1
  for(y in (maxyear-1):1984){
    for(a in 6:3){  
      N7[y-minYearm1,a-minAgem1]<-N7[y+1-minYearm1,a+1-minAgem1]*exp(oldM[y-minYearm1,a-minAgem1])+
                                  modelC[y-minYearm1,a-minAgem1]*exp(0.5*oldM[y-minYearm1,a-minAgem1])+
                                  con[y-minYearm1,a-minAgem1]*exp(0.5*oldM[y-minYearm1,a-minAgem1])
    }
  }

  MM<-oldF*0
  for(y in (maxyear-1):1984){
    for(a in 6:3){  
      MM[y-minYearm1,a-minAgem1]<-log(N7[y-minYearm1,a-minAgem1]/N7[y+1-minYearm1,a+1-minAgem1])-oldF[y-minYearm1,a-minAgem1]-oldM[y-minYearm1,a-minAgem1]
    }
  }

  MM[nrow(MM),] <- MM[nrow(MM)-1,] 
  MM[MM<0]<-0

  newM<-nmOldAdj

  ridx<-rownames(newM)%in%rownames(MM)
  cidx<-colnames(newM)%in%colnames(MM)
  newM[ridx,cidx]<- MM+.2
  result[[run]]<-runwithnewm(fit, newM)
}


Rsum<-lapply(result, function(x){xx<-rectable(x); sum(xx[rownames(xx)%in%(1984:maxyear),1])})


fit<-result[[5]]


save(fit, con, file="model.RData")