Data/WestIce2012Old/ADMButils.R

In NorskRegnesentral/pupR: Harp and hooded seal pup production estimation

# R interface to AD Model Builder files
# Author: Hans J. Skaug
# Version: Nov 07.

# Write to .dat file
# Example: dat_write("epil5sim.dat",list(n=n,p=6,q=q,X=as.matrix(d)))
dat_write = function(name,L)
{
  n = nchar(name)
  if(substring(name,n-3,n)==".dat")
    file_name = name
  else
    file_name = paste(name,".dat",sep="")

  cat("# \"",name,".dat\" produced by dat_write() from ADMButils; ",date(),"\n", file=file_name,sep="")
  for(i in 1:length(L))
  {
    x = L[[i]]

    if(data.class(x)=="numeric")
      cat("#",names(L)[i],"\n",x,"\n\n",file=file_name,append=T)

    if(data.class(x)=="matrix")
    {
      cat("#",names(L)[i],"\n",file=file_name,append=T)
      write.table(x,,col=F,row=F,quote=F,file=file_name,append=T)
      cat("\n",file=file_name,append=T)
    }

    if(data.class(x)=="list")
    {
      cat("#",names(L)[i],"\n",file=file_name,append=T)
      for(j in 1:length(x))
        if(is.numeric(x[[j]]))
          cat(x[[j]],"\n",file=file_name,append=T)
	else
	  stop("List with non-numeric elements not yet implemented")
      cat("\n",file=file_name,append=T)
    }

  }
}

# Write to file
# Example: file_write("epil5sim.dat",list(n=n,p=6,q=q,X=as.matrix(d)))
file_write = function(name,L)
{
  n = nchar(name)
  cat("#Produced by file_write() from ADMButils; ",date(),"\n", file=name,sep="")
  for(i in 1:length(L))
  {
    x = L[[i]]

    if(data.class(x)=="numeric")
      cat("#",names(L)[i],"\n",x,"\n\n",file=name,append=T)

    if(data.class(x)=="matrix")
    {
      cat("#",names(L)[i],"\n",file=name,append=T)
      write.table(x,file =name,append = T,row.names = FALSE, col.names=FALSE)
      cat("\n",file=name,append=T)
    }

    if(data.class(x)=="list")
    {
      cat("#",names(L)[i],"\n",file=name,append=T)
      for(j in 1:length(x))
        if(is.numeric(x[[j]]))
          cat(x[[j]],"\n",file=name,append=T)
	else
	  stop("List with non-numeric elements not yet implemented")
      cat("\n",file=name,append=T)
    }

  }
}


# Write to pin-file
# Example: pin_write("kalman_ar1.pin",list(log_sigma=c(0,0),a=0,p=rep(1,10)))
pin_write = function(name,L)
{
  n = nchar(name)
  if(substring(name,n-3,n)==".pin")
    file_name = name
  else
    file_name = paste(name,".pin",sep="")

  cat("# \"",name,".pin\" produced by pin_write() from ADMButils; ",date(),"\n", file=file_name,sep="")
  for(i in 1:length(L))
  {
    x = L[[i]]

    if(data.class(x)=="numeric")
      cat("#",names(L)[i],"\n",L[[i]],"\n\n",file=file_name,append=T)

    if(data.class(x)=="matrix")
    {
      cat("#",names(L)[i],"\n",file=file_name,append=T)
      write.table(L[[i]],,col=F,row=F,quote=F,file=file_name,append=T)
      cat("\n",file=file_name,append=T)
    }
  }
}

# Read par-file (or files with the same format)
#Note: matrices must be handeled by the "ncol" argument (se example below)

# Examples:
#  par_read("sea16.par")
#  par_read("sea16.rep",ncols=list(N=3,ogives=2,age_dist=50)) # N,ogives,age_dist are matrices

par_read = function(name,ncols=list()) # matrices must be specified by name and ncol
{
  n = nchar(name)
  endelse = substring(name,max(1,n-3),n)
  har_endlese = (substring(endelse,1,1) == ".")
  if(har_endlese)
    file_name = name
  else
    file_name = paste(name,".par",sep="")

  tmp = scan(file_name,what="",quiet=T)
  tmp2 = split(tmp,cumsum(tmp=="#"))
  x = tmp2

  if(endelse ==".par")
    x = x[-1]

  for(i in 1:length(x))
  {
    y = x[[i]]
    n = nchar(y[2])
    x[[i]] = as.numeric(y[-(1:2)])
    names(x)[i] = substring(y[2],1,n-1)
  }

  # Convert to matrix for those arguments relevant
  if(length(ncols)>0)
    for(i in 1:length(ncols))
    {
      NN = names(ncols)[i]
      x[[NN]] <- matrix(x[[NN]],ncol=ncols[[i]],byrow=T)
    }

  if(endelse == ".par")
  {
    x$loglik = -as.numeric(tmp2[[1]][11])
    x$gradient = -as.numeric(tmp2[[1]][16])
  }

  x
} 


# Reads std-file
std_read = function(name)
{
  n = nchar(name)
  if(substring(name,n-3,n)==".std")
    file_name = name
  else
    file_name = paste(name,".std",sep="")

  tmp = read.table(file_name,skip=1)
  est = tmp[,3]
  names(est) = tmp[,2]
  std = tmp[,4]
  names(std) = tmp[,2]

  L1 = list()
  L2 = list()
  
  for(i in unique(names(std)))
  {
    L1[[i]] = est[names(est)==i]
    L2[[i]] = std[names(std)==i]
  }

  list(est=L1,std=L2)
}


# HJS utility functions; actually part of R it turns out
cov2corr <- function(m) diag(1/sqrt(diag(m))) %*% m %*% diag(1/sqrt(diag(m)))
member <- function(x,y) !is.na(match(x,y))     
below <- function(n,strictly=F)
{
  M <- matrix(T,n,n)
  M[rep(1:n,n)<rep(1:n,rep(n,n))] <- F
  if(strictly)
    diag(M) = F
  M
}


# read Hessian of dimension n from .cor file
readH <- function(file,n,cor=F)
{
  N = n*(n+1)/2+4*n
  tmp = scan(file,what="",skip=2,quiet=T)
  if(length(tmp)<N) stop("n is too large")
  tmp = tmp[1:N]

  stdtab = numeric(n)
  H = diag(n)

  for(i in 1:n)
  {
    stdtab[i] = as.numeric(tmp[4])
    tmp = tmp[-(1:4)]
    H[i,1:i] = as.numeric(tmp[1:i])
    tmp = tmp[-(1:i)]
  }

  if(length(tmp)!=0)
  {
    print(length(tmp))
    stop("Det er noe galt")
  }
 
  H = H+t(H)	# Fill in upper diagonal
  diag(H) = 1
  

  if(!cor)
    H = diag(stdtab) %*% H %*% diag(stdtab)

  H
}



myplot <- function(Year,STD,i2=1,Areaname = "",Pupest=c(1),showtotals = FALSE,showleg = FALSE,showscale = FALSE,ylim){
	#plot(Year[1:i2],N1[1:i2],type="l",xlim = c(1978,max(Year+2)),ylim = c(0,max((N1+N0))+10),col = "blue",xlab = "Year",main = Areaname,ylab = "Number of animals",lwd = 2,bty = "l")
	#plot(Year[1:i2],N1[1:i2],type="l",xlim = c(1978,2013),ylim = c(0,max((N1[1:i2]))+10),col = "blue",xlab = "Year",main = Areaname,ylab = "Number of animals",lwd = 2,cex.lab = 1.5,cex.main = 1.5,bty = "l")
	plot(Year[1:i2],STD$N0[1:i2],type="l",xlim = c(1978,2010),ylim = c(0,ylim),col = "blue",xlab = "Year",main = Areaname,ylab = "Number of animals",lwd = 2,cex.lab = 1.5,cex.main = 1.5,bty = "l")
	#lines(REP$Year[(i2+1):length(REP$Year)],REP$N1[(i2+1):length(REP$Year)],type="l",lty = 2, col = "blue")
	#lines(Year,N1,type="l",lty = 2,lwd = 2, col = "blue")
	#lines(REP$Year[1:i2],REP$N1[1:i2]-1.96*STD$sigmaN1[1:i2],type = "l",lty=2,col="blue")
	#lines(REP$Year[1:i2],REP$N1[1:i2]+1.96*STD$sigmaN1[1:i2],type = "l",lty=2,col="blue")
	#lines(Year[1:i2],N0[1:i2],type="l",col="darkgreen",lwd = 2)
	##lines(Year,N0,type="l",col="darkgreen",lty = 2,lwd = 2)
	
	if(sum(dim(Pupest))>1) {
	lines(Pupest[,1],Pupest[,2],type = "p",col = "red",pch=16)
	#lines(1999*array(1,100),seq(1,(max(REP$N1)+10),length.out=100),lty = 2, col = "pink")
	}
	
	if(showtotals){
		lines(Year[1:i2],N1[1:i2]+N0[1:i2],type="l",lty=1,lwd =2, col = "orange2")
		lines(Year,N1+N0,type="l",lty=2,lwd =2, col = "orange2")

		}
	
	
	if(showleg){
		legend(1979,max(N1)+15,legend=c("Population trajectory N1+","Population trajectory pups"),col=c("blue","darkgreen"),lty = c(1,1),lwd = c(2,2),bty="n",cex = 1.0)
		mtext(paste("A =", A,", M = ",M,", M0 = ", M0),3)

		}
	
	if(showscale){
		lines(Year[1:i2],N0[1:i2]*3.5,type = "l",lty = 1,col = "lightblue",lwd=2)
		lines(Year[1:i2],N0[1:i2]*4.7,type = "l",lty = 1,col = "midnightblue",lwd=2)
		lines(Year,N0*3.5,type = "l",lty = 2,col = "lightblue",lwd=2)
		lines(Year,N0*4.7,type = "l",lty = 2,col = "midnightblue",lwd=2)
		if(showtotals){
			legend(1979,max(N1)+15,legend=c("Modelestimate, Total","Model estimate, N1+","Model estimate, pups","Scaled estimate, c = 3.5","Scaled estimate, c = 4.7"),col=c("orange2","blue","darkgreen","lightblue","midnightblue"),lty = c(1,1,1,1,1),lwd = c(2,2,2,2,2),bty="n",cex = 1.0)
			} else {
		legend(1979,max(N1)+15,legend=c("Model estimate N1+","Model estimate pups","Scaled estimate, c = 3.5","Scaled estimate, c = 4.7"),col=c("blue","darkgreen","lightblue","midnightblue"),lty = c(1,1,2,2),lwd = c(2,2,2,2),bty="n",cex = 1.0)}
		}
	
	#mtext(paste("M = ",STD$M,", M0 = ", STD$M0),3)
	#mtext(paste("A =", A,", M = ",M,", M0 = ", M0,", D = ",STD$D),3)
	#legend(1989,max(REP$N1)+15,legend=c("Population trajectory N1+","Population trajectory pups","Pup production estimates"),col=c("blue","darkgreen","red"),lty = c(1,1,NA),lwd = c(2,2,NA),pch=c(NA,NA,16),bty="n",cex = 1.0)
	}