load("fsa.RData") # gets "dat"

par <- list(
  logN1Y=rep(0,nrow(dat$M)),
  logN1A=rep(0,ncol(dat$M)-1),
  logFY=rep(0,ncol(dat$M)),
  logFApar=rep(0,4),
  logSdCatch=0, 
  logQ=rep(0,length(unique(dat$age[dat$fleet==2]))),
  logSdSurvey=0
)

getAll<-function(...) ## helper function borrowed from RTMB  
{
    fr <- parent.frame()
    x <- c(...)
    nm <- names(x)
    for (i in seq_along(x)) {
        fr[[nm[i]]] <- x[[i]]
    }
    invisible(NULL)
}

nll<-function(parvec){
  par<-relist(parvec, skeleton=par)  
  getAll(par, dat)

  na <- max(age)-min(age)+1
  ny <- max(year)-min(year)+1

  ## setup F
  logFA <- c(logFApar,0,0,0)
  F <- outer(exp(logFA),exp(logFY))

  ## setup N
  logN <- matrix(0, nrow=na, ncol=ny)
  logN[,1] <- logN1Y  
  for(y in 2:ny){
    logN[1,y] <- logN1A[y-1]
    for(a in 2:na){
      logN[a,y] <- logN[a-1,y-1]-F[a-1,y-1]-M[a-1,y-1]
    }
  } 

  # Match to observations
  logObs <- log(obs)
  logPred <- numeric(length(logObs))
  sdvec <- numeric(length(logObs))
  for(i in 1:length(logObs)){
    a <- age[i]-min(age)+1
    y <- year[i]-min(year)+1
    if(fleet[i]==1){
      logPred[i] <- log(F[a,y])-log(F[a,y]+M[a,y])+log(1-exp(-F[a,y]-M[a,y]))+logN[a,y]
      sdvec[i] <- exp(logSdCatch)
    }else{
      logPred[i] <- logQ[a]-(F[a,y]+M[a,y])*surveyTime+logN[a,y]
      sdvec[i] <- exp(logSdSurvey)
    }    
  }
  
  ans <- -sum(dnorm(logObs,logPred,sdvec,TRUE))

  logssb <<- log(apply(exp(logN)*stockMeanWeight*propMature,2,sum))
  return(ans)
}

opt <- nlminb(unlist(par), nll, control=list(iter.max=1000,eval.max=1000))

library(numDeriv)
getLogSSB<-function(theta){nll(theta);logssb}
J<-jacobian(getLogSSB, unlist(opt$par))
H<-hessian(nll, opt$par)
sd<-sqrt(diag(J%*%solve(H)%*%t(J)))

yr<-sort(unique(dat$year))
plot(yr, exp(logssb), type="l", lwd=5, col="red", ylim=c(0,550000), xlab="Year", ylab="SSB")
lines(yr, exp(logssb-2*sd), type="l", lwd=1, col="red")
lines(yr, exp(logssb+2*sd), type="l", lwd=1, col="red")