nephOM <- function(Stock, SR, Sel, recmod = 'srr', nits = 1, trec = 1, tspwn = 1, omega = 0.5) {
 
	# checks
 
	if(!is.FLStock(Stock)) 
		stop("Stock is not an FLStock object")
 
	if(!is.FLSR(SR))
		stop("SR not an FLSR object")
 
	if(recmod != 'srr' && recmod != 'fixed' && recmod != 'bootstrap')
		stop("No recruitment dynamics specified")
 
	if(!all(dim(Stock@stock.n) == dim(Sel)))
		stop("Dimensions of stock and selctivity quants not the same")
 
	# take what we need and strip out the dimnames
 
	dms.n <- c(dim(stock.n(Stock))[1:5],nits)
	dms <- c(dim(catch(Stock))[1:5],nits)
	stk.n <- array(dim=dms.n)
	stk.n[,,,,,1] <- stock.n(Stock)@.Data[,,,,,1]
	stk.c <- array(dim=dms)
	stk.c[] <- catch(Stock)@.Data[]
	sel <- array(dim=c(dim(Sel)[1:5],nits))
	sel[] <- Sel@.Data[]
	m <- array(dim=dms.n)
	m[] <- m(Stock)@.Data[]
	m.spwn <- array(dim=dms.n)
	m.spwn[] <- m.spwn(Stock)@.Data[]
	stk.e <- array(dim=dms)
	stk.ssb <- stk.e
	h <- stk.n
	h[] <- switch(units(harvest(Stock)),'f' = 1-exp(-harvest(Stock)@.Data[]),'hr' = harvest(Stock)@.Data[])
	stk.cn <- stk.n
	mat <- array(dim=dms.n)
	mat[] <- mat(Stock)@.Data[]
	stk.wt <- array(dim=dms.n)
	stk.wt[] <- stock.wt(Stock)@.Data[]
 
	# sr params and checks
 
	rec.age <- dims(stock.n(Stock))[['min']]
	if(recmod == 'srr') {
 
		sr.mod <- model(SR)
		sr.par <- params(SR)
		sr.var <- var(SR)
		sr.vara <- varacorr(SR)
		sr.rho <- sqrt(1-sr.vara/sr.var)
 
		sr.func <- function(.ssb,sr.mod,sr.par) {
 
			if(sr.mod == 'bevholt') 
				.rec <- sr.par[,'alpha'] * .ssb[] /(sr.par[,'beta'] + .ssb[])
			if(sr.mod == 'ricker')
				.rec <- sr.par[,'alpha'] * .ssb[] * exp(-sr.par[,'beta'] * .ssb[])
			if(sr.mod == 'segreg')
				.rec <- apply(.ssb,c(1:6),function(x,sr.par){x <- ifelse(x < sr.par[,'beta'],sr.par[,'alpha']*x,sr.par[,'alpha']*sr.par[,'beta'])},sr.par)
 
			return(.rec)
		} 
	}
 
	if(recmod == 'fixed') {
 
		if(is.na(var(SR)) | is.na(varacorr(SR))) 
			stop("Error in nephOM: variance and auto-correlated variance are missing")
		R0m <- as.double(stock.n(Stock)[1,1,1,trec,,])
		R0f <- as.double(stock.n(Stock)[1,1,1,trec,,])
		sr.var <- var(SR)
		sr.vara <- varacorr(SR)
		sr.rho <- ifelse((sr.var && sr.vara) > 0,sqrt(1-sr.vara/sr.var),0) 
	}
 
	if(recmod == 'bootstrap') {
 
		if(!all(dim(rec(SR)) == dim(stock.n(Stock)[1,,,,,])))
			stop("Error in nephOM: bootstrap recruitment series not same as recruit dimensions of stock object")
	}
 
	# now we perform the stock dynamics going forward
	# ::
	# eqm (exploited or unexploited) initial dynamics first
 
	amax <- dms.n[1]
	ymax <- dms.n[2]
	ns <- dms.n[4]
 
	if(trec > 1) 
		stk.n[1,1,,1:(trec-1),,] <- 0
 
	if(recmod != 'bootstrap') {
		eps <- exp(rnorm(nits,0,sqrt(sr.vara)))
		stk.n[1,1,1,trec,,] <- R0m * eps[]
		stk.n[1,1,2,trec,,] <- R0f * eps[]
	} else {
		.rec <- boot(rec(SR),~year,nits)
		stk.n[1,1,1,trec,,] <- .rec[,1,1,trec,,]
		stk.n[1,1,2,trec,,] <- .rec[,1,2,trec,,] 
	}
 
	# switch for seasonal/non-seasonal model
 
	if(ns > 1) {
 
		if(trec < ns) {
			for(s in (trec+1):ns) 
				stk.n[1,1,,s,,] <- stk.n[1,1,,s-1,,] * exp(-m[1,1,,s-1,,]) * (1-h[1,1,,s-1,,])
		}
 
		for(a in 2:amax) {
 
			stk.n[a,1,,1,,] <- stk.n[a-1,1,,ns,,] * exp(-m[a-1,1,,ns,,]) * (1-h[a-1,1,,ns,,])
			if(a == amax) {
 
				tmp1 <- 1
				for(s in 1:ns)
					tmp1 <- tmp1 * exp(-m[a,1,,s,,]) * (1-h[a,1,,s,,])
				stk.n[a,1,,1,,] <- stk.n[a-1,1,,ns,,] * exp(-m[a-1,1,,ns,,]) * (1-h[a-1,1,,ns,,]) / (1 - tmp1)
			}
 
			for(s in 2:ns) 
				stk.n[a,1,,s,,] <- stk.n[a,1,,s-1,,] * exp(-m[a,1,,s-1,,]) * (1-h[a,1,,s-1,,])
		}
 
	} else {
 
		# non seasonal
 
		for(a in 2:amax)
			stk.n[a,1,,,,] <- stk.n[a-1,1,,,,] * exp(-m[a-1,1,,,,]) * (1-h[a-1,1,,,,])
 
		# plus group
 
		stk.n[amax,1,,,,] <- stk.n[amax,1,,,,]/(1-exp(-m[amax,1,,,,]) * (1-h[amax,1,,,,]))
	}
 
	# calculate the initial harvest rates
	# :: 
	# first exploitable biomass
 
 
	stk.e[,1,,,,] <- quantSums(FLQuant(stk.n * stk.wt * sel))@.Data[,1,,,,]
	h.tot <- stk.e
	h.tot[,1,,,,] <- stk.c[,1,,,,]/stk.e[,1,,,,]
	h.tot <- apply(h.tot,c(1:6),function(x) {x <- min(1,x)})
	for(a in 1:amax)
		h[a,1,,,,] <- h.tot[,1,,,,] * sel[a,1,,,,]
 
	# initial ssb
 
	stk.ssb[,1,,,,] <- quantSums(FLQuant(stk.n * stk.wt * mat * exp(-m * m.spwn)))@.Data[,1,,,,]
 
	# catches
 
	stk.cn[,1,,,,] <- stk.n[,1,,,,] * h[,1,,,,]
 
	# loop through the years
 
	for(y in 2:ymax) {
 
		# no recruits before recruitment season
 
		stk.n[1,y,,1:(trec-1),,] <- 0
 
		# season 1 stuff
 
		stk.n[-c(1),y,,1,,] <- stk.n[-c(amax),y-1,,ns,,] * exp(-m[-c(amax),y-1,,ns,,]) * (1-h[-c(amax),y-1,,ns,,])
 
		# plus group
 
		stk.n[amax,y,,1,,] <- stk.n[amax,y,,1,,] + stk.n[amax,y-1,,ns,,] * exp(-m[amax,y-1,,ns,,]) * (1-h[amax,y-1,,ns,,])
 
		# ssb
 
		stk.ssb[,y,,1,,] <- quantSums(FLQuant(stk.n * stk.wt * mat * exp(-m * m.spwn)))@.Data[,y,,1,,]
 
		if(trec == 1) {
 
			# recruitment
 
			if(recmod == 'srr') {
 
				# first create the autocorrelated noise terms
 
				eps <- exp(sr.rho * log(eps) + rnorm(nits,0,sqrt(sr.vara)))
				rec.mf <- sr.func(stk.ssb[,y-rec.age,2,tspwn,,],sr.mod,sr.par) * eps
 
				# now use omega to split them male and female
 
				stk.n[1,y,1,trec,,] <- omega * rec.mf[]
				stk.n[1,y,2,trec,,] <- (1-omega) * rec.mf[]
 
			} 
 
			if(recmod == 'fixed') {
 
				# first create the autocorrelated noise terms
 
				eps <- exp(sr.rho * log(eps) + rnorm(nits,0,sqrt(sr.vara)))
 
				stk.n[1,y,1,trec,,] <- R0m * eps[]
				stk.n[1,y,2,trec,,] <- R0f * eps[] 
			}
 
			if(recmod == 'bootstrap') {
 
				.rec <- boot(rec(SR),~year,nits)
				stk.n[1,y,1,trec,,] <- .rec[,y,1,trec,,]
				stk.n[1,y,2,trec,,] <- .rec[,y,2,trec,,] 
			}
		}
 
		# harvest rates
 
		stk.e[,y,,1,,] <- quantSums(FLQuant(stk.n * stk.wt * sel))@.Data[,y,,1,,]
		h.tot[,y,,1,,] <- stk.c[,y,,1,,]/stk.e[,y,,1,,]
		h.tot <- apply(h.tot,c(1:6),function(x) {x <- min(1,x)})
		for(a in 1:amax)
			h[a,y,,1,,] <- h.tot[,y,,1,,] * sel[a,y,,1,,] 
 
		if(ns > 1) {	
 
			# seasons 2 to ns now if needed
 
			for(s in 2:ns) {
 
				# recruits after spawning
 
				if(s > trec)
					stk.n[1,y,,s,,] <- stk.n[1,y,,s-1,,] * exp(-m[1,y,,s-1,,]) * (1-h[1,y,,s-1,,])
 
				# rest of the ages
 
				stk.n[-c(1),y,,s,,] <- stk.n[-c(1),y,,s-1,,] * exp(-m[-c(1),y,,s-1,,]) * (1-h[-c(1),y,,s-1,,])
 
				# ssb
 
				stk.ssb[,y,,s,,] <- quantSums(FLQuant(stk.n * stk.wt * mat * exp(-m.spwn * m)))@.Data[,y,,s,,]
 
				if(s == trec) {
 
					# recruitment
 
					if(recmod == 'srr') {
 
						# first create the autocorrelated noise terms
 
						eps <- exp(sr.rho * log(eps) + rnorm(nits,0,sqrt(sr.vara)))
						rec.mf <- sr.func(stk.ssb[,y-rec.age,2,tspwn,,],sr.mod,sr.par) 
 
						# now use omega to split them male and female
 
						stk.n[1,y,1,trec,,,] <- omega * rec.mf[] * eps[]
						stk.n[1,y,2,trec,,,] <- (1-omega) * rec.mf[] * eps[]
					} 
 
					if(recmod == 'fixed') {
 
						# first create the autocorrelated noise terms
 
						eps <- exp(sr.rho * log(eps) + rnorm(nits,0,sqrt(sr.vara)))
 
						stk.n[1,y,1,trec,,] <- R0m * eps[]
						stk.n[1,y,2,trec,,] <- R0f * eps[] 
					}	
 
					if(recmod == 'bootstrap') {
 
						.rec <- boot(rec(SR),~year,nits)
						stk.n[1,y,1,trec,,] <- .rec[,y,1,trec,,]
						stk.n[1,y,2,trec,,] <- .rec[,y,2,trec,,] 
					} 
				}
 
				# harvest rates
 
				stk.e[,y,,s,,] <- quantSums(FLQuant(stk.n * stk.wt * sel))@.Data[,y,,s,,]
				h.tot[,y,,s,,] <- stk.c[,y,,s,,]/stk.e[,y,,s,,]
				h.tot <- apply(h.tot,c(1:6),function(x) {x <- min(1,x)})
				for(a in 1:amax)
					h[a,y,,s,,] <- h.tot[,y,,s,,] * sel[a,y,,s,,]
			}
		}
 
		# catches
 
		stk.cn[,y,,,,] <- stk.n[,y,,,,] * h[,y,,,,]
	}
 
	# now we create new FLQuants to return back in the res stock object
 
	res <- Stock
	dimn.n <- dimnames(stock.n(res))
	dimn <- dimnames(stock(res))
	dimn.n$iter <- seq(1,nits,1)
	dimn$iter <- seq(1,nits,1)
	stock.n(res) <- FLQuant(quant='age',stk.n,dimnames=dimn.n)
	catch.n(res) <- FLQuant(quant='age',stk.cn,dimnames=dimn.n) 
	harvest(res) <- FLQuant(quant='age',h,dimnames=dimn.n) 
	units(harvest(res)) <- 'hr'
	stock(res) <- FLQuant(quant='all',stk.e,dimnames=dimn) 
 
	return(res)
}
 
efimas1/wp4/cs4/main/bom.txt · Last modified: 2007/04/02 18:32 by rh
 
Except where otherwise noted, content on this wiki is licensed under the following license:CC Attribution-Noncommercial-Share Alike 3.0 Unported
Recent changes RSS feed Donate Powered by PHP Valid XHTML 1.0 Valid CSS Driven by DokuWiki