library(stockassessment)
load("run/model.RData")


pearsonresiduals<-function(object){
  o <- object$data$logobs
  p <- object$rep$predObs
  oo <- o;
  oo[is.na(o)]<-p[is.na(o)]
  r <- oo-p
  cholP <- lapply(object$rep$obsCov, function(S)chol(solve(S)))
  for(i in 1:nrow(object$data$idx1)){
    for(j in 1:ncol(object$data$idx1)){
      if(!is.na(object$data$idx1[i,j])){
        i1<-object$data$idx1[i,j]+1
        i2<-object$data$idx2[i,j]+1
        r[i1:i2] <- cholP[[i]]%*%r[i1:i2]
      }
    }
  }
  ret <- as.data.frame(cbind(object$data$aux, residual=r))
  attr(ret, "fleetNames") <- attr(object$data, "fleetNames")
  class(ret) <- "samres"
  ret
}


RES<-pearsonresiduals(fit)
RESP<-procres(fit)
save(RES, RESP, file="run/residuals.RData")