From 5c972e1d10f376a1e2a7a0352af5b09fffbc29a3 Mon Sep 17 00:00:00 2001 From: Philipp Boersch-Supan Date: Mon, 18 Dec 2017 16:08:14 -0500 Subject: [PATCH] initial commit with code pasted from DEBtool_R --- .Rbuildignore | 2 + .gitignore | 4 + DEButilities.Rproj | 20 +++ DESCRIPTION | 10 ++ NAMESPACE | 1 + R/beta0.R | 52 ++++++ R/del.R | 14 ++ R/dget_l_ISO.R | 19 +++ R/dget_l_V1.R | 17 ++ R/dget_lb2.R | 15 ++ R/fnget_lb2.R | 18 ++ R/get_lb.R | 131 ++++++++++++++ R/get_lb2.R | 65 +++++++ R/get_lj.R | 125 ++++++++++++++ R/get_tb.R | 79 +++++++++ R/get_tj.R | 173 +++++++++++++++++++ R/get_tm.R | 151 +++++++++++++++++ R/get_tm_s.R | 147 ++++++++++++++++ R/get_ue0.R | 67 ++++++++ R/initial_scaled_reserve.R | 67 ++++++++ R/nmregr.R | 337 +++++++++++++++++++++++++++++++++++++ R/nmregr_options.R | 108 ++++++++++++ R/repro_rate_metam.R | 89 ++++++++++ R/spline1.R | 60 +++++++ R/tempcorr.R | 45 +++++ 25 files changed, 1816 insertions(+) create mode 100644 .Rbuildignore create mode 100644 .gitignore create mode 100644 DEButilities.Rproj create mode 100644 DESCRIPTION create mode 100644 NAMESPACE create mode 100644 R/beta0.R create mode 100644 R/del.R create mode 100644 R/dget_l_ISO.R create mode 100644 R/dget_l_V1.R create mode 100644 R/dget_lb2.R create mode 100644 R/fnget_lb2.R create mode 100644 R/get_lb.R create mode 100644 R/get_lb2.R create mode 100644 R/get_lj.R create mode 100644 R/get_tb.R create mode 100644 R/get_tj.R create mode 100644 R/get_tm.R create mode 100644 R/get_tm_s.R create mode 100644 R/get_ue0.R create mode 100644 R/initial_scaled_reserve.R create mode 100644 R/nmregr.R create mode 100644 R/nmregr_options.R create mode 100644 R/repro_rate_metam.R create mode 100644 R/spline1.R create mode 100644 R/tempcorr.R diff --git a/.Rbuildignore b/.Rbuildignore new file mode 100644 index 0000000..91114bf --- /dev/null +++ b/.Rbuildignore @@ -0,0 +1,2 @@ +^.*\.Rproj$ +^\.Rproj\.user$ diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..5b6a065 --- /dev/null +++ b/.gitignore @@ -0,0 +1,4 @@ +.Rproj.user +.Rhistory +.RData +.Ruserdata diff --git a/DEButilities.Rproj b/DEButilities.Rproj new file mode 100644 index 0000000..497f8bf --- /dev/null +++ b/DEButilities.Rproj @@ -0,0 +1,20 @@ +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +AutoAppendNewline: Yes +StripTrailingWhitespace: Yes + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source diff --git a/DESCRIPTION b/DESCRIPTION new file mode 100644 index 0000000..3bf0c80 --- /dev/null +++ b/DESCRIPTION @@ -0,0 +1,10 @@ +Package: DEButilities +Type: Package +Title: Utility functions for Dynamic Energy Budget Theory +Version: 0.1.0 +Author: DEBtool maintainers +Maintainer: Philipp Boersch-Supan +Description: Functions copied from the DEBtool_R git repository and repackaged in a CRAN compliant manner to simplify parameter estimation for DEB models. This package is not intended to replace the DEBtool_R translation effort, but rather a quick fix to make functions for initial value calculation accessible in R. +License: What license is it under? +Encoding: UTF-8 +LazyData: true diff --git a/NAMESPACE b/NAMESPACE new file mode 100644 index 0000000..d75f824 --- /dev/null +++ b/NAMESPACE @@ -0,0 +1 @@ +exportPattern("^[[:alpha:]]+") diff --git a/R/beta0.R b/R/beta0.R new file mode 100644 index 0000000..4cc91e4 --- /dev/null +++ b/R/beta0.R @@ -0,0 +1,52 @@ +## beta0 +# particular incomplete beta function + +## +#function f = + +beta0 = function(x0,x1){ + # created 2000/08/16 by Bas Kooijman; modified 2011/04/10 + + ## Syntax + # f = <../beta.m *beta0*> (x0,x1) + + ## Description + # particular incomplete beta function: + # B_x1(4/3,0) - B_x0(4/3,0) = \int_x0^x1 t^(4/3-1) (1-t)^(-1) dt + + # Input + # + # * x0: scalar with lower boundary for integration + # * x1: scalar with uper boundary for integration + # + # Output + # + # * f: scalar with particular incomple beta function + + ## Remarks + # See also <../lib/misc/beta_34_0 *beta_34_0*> + + ## Example of use + # beta0(0.1, 0.2) + + if (x0 < 0 || x0 >= 1 || x1 < 0 || x1 >= 1){ + print('Warning from beta0: argument values outside (0,1) \n') + f = NA + } + + n0 = length(x0) + n1 = length(x1) + if (n0 != n1 && n0 != 1 && n1 != 1) { + print('Warning from beta0: argument sizes do not match \n') + f = NA + } + + x03 = x0 ^ (1/ 3) + x13 = x1 ^ (1/ 3) + a3 = sqrt(3) + f1 = - 3 * x13 + a3 * atan((1 + 2 * x13)/ a3) - log(as.complex(x13 - 1)) + log(as.complex(1 + x13 + x13 ^ 2))/ 2 + f0 = - 3 * x03 + a3 * atan((1 + 2 * x03)/ a3) - log(as.complex(x03 - 1)) + log(as.complex(1 + x03 + x03 ^ 2))/ 2 + f = Re(f1 - f0) + + +} diff --git a/R/del.R b/R/del.R new file mode 100644 index 0000000..640de19 --- /dev/null +++ b/R/del.R @@ -0,0 +1,14 @@ +del <- function(t, el) { + # created 2000/08/21 by Bas Kooijman + # differential equations for change in scaled reserve and length + # during the embryonic period; time is scaled with k_M + # d = [d/dtau e*l^3, d/dtau l]; el = [e*l^3, l] + # the first element is not d/dt e, because e -> infty for t -> 0 + + d = matrix(0,2,1) + + l3 = el[2]^3 + d[2] = (g/3)*(el[1] - el[1]*l3)/(el[1] + g*l3) + d[1] = (3*d[2] - g)*el[1]/el[1] + return(d) +} \ No newline at end of file diff --git a/R/dget_l_ISO.R b/R/dget_l_ISO.R new file mode 100644 index 0000000..fb7d96f --- /dev/null +++ b/R/dget_l_ISO.R @@ -0,0 +1,19 @@ +dget_l_ISO= function(vH, l, pars){ +# vH: scalar with scaled maturity +# l: scalar with scaled length +# dl: scalar with d/dvH l +# called from get_lp, get_lj, get_ls + + k=pars[1] + lT=pars[2] + g=pars[3] + f=pars[4] + sM=pars[5] + + + r = g * (f * sM - lT * sM - l)/ l/ (f + g) # specific growth rate + dl = l * r/ 3 # d/dt l + dvH = f * l^2 * (sM - l * r/ g) - k * vH # d/dt vH + dl = dl/ dvH # dl/ dvH + return(list(dl)) +} diff --git a/R/dget_l_V1.R b/R/dget_l_V1.R new file mode 100644 index 0000000..6887d2d --- /dev/null +++ b/R/dget_l_V1.R @@ -0,0 +1,17 @@ +dget_l_V1= function(vH, l, pars){ + # vH: -, scaled maturity + # l: -, scaled length + # lref: -,scaled length when acceleration starts (can be lb or ls) + # dl: d/dvH l during exponential growth + # called from get_lp, get_lj, get_ls + k=pars[1] + lT=pars[2] + g=pars[3] + f=pars[4] + lref=pars[5] + r = g * (f - lT - lref)/ lref/ (g + f) # specific growth rate + dl = r * l/ 3 # d/dt l + dvH = f * l^3 * (1/ lref - r/ g) - k * vH # d/dt vH + dl = dl/ dvH # dl/ dvH + return(list(dl)) +} \ No newline at end of file diff --git a/R/dget_lb2.R b/R/dget_lb2.R new file mode 100644 index 0000000..8e81acc --- /dev/null +++ b/R/dget_lb2.R @@ -0,0 +1,15 @@ + +dget_lb2=function(x, y, pars){ + lb=Re(pars[1]) + xb=pars[2] + xb3=pars[3] + g=pars[4] + k=pars[5] + # y = x e_H; x = g/(g + e) + # (x,y): (0,0) -> (xb, xb eHb) + + x3 = x^(1/ 3) + l = x3/ (xb3/ lb - beta0(x, xb)/ 3/ g) + dy = l + g - y * (k - x)/ (1 - x) * l/ g/ x + return(list(dy)) +} \ No newline at end of file diff --git a/R/fnget_lb2.R b/R/fnget_lb2.R new file mode 100644 index 0000000..9dda088 --- /dev/null +++ b/R/fnget_lb2.R @@ -0,0 +1,18 @@ + +fnget_lb2=function(lb, xbxb3gvHbk){ + + xb=xbxb3gvHbk[1] + xb3=xbxb3gvHbk[2] + g=xbxb3gvHbk[3] + vHb=xbxb3gvHbk[4] + k=xbxb3gvHbk[5] + + # f = y(x_b) - y_b = 0; x = g/ (e + g); x_b = g/ (e_b + g) + # y(x) = x e_H = x g u_H/ l^3 and y_b = x_b g u_H^b/ l_b^3 + + tspan=seq(1e-10, xb, length=100) + yini=c(y=0) + xy = as.data.frame(ode(yini, tspan, dget_lb2, parms=c(lb, xb, xb3, g, k), method="ode23")) + f = xy$y[length(xy$y)] - xb * g * vHb/ lb^3 + return(f) +} \ No newline at end of file diff --git a/R/get_lb.R b/R/get_lb.R new file mode 100644 index 0000000..b7b590b --- /dev/null +++ b/R/get_lb.R @@ -0,0 +1,131 @@ +## get_lb +# Obtains scaled length at birth, given the scaled reserve density at birth + +## +get_lb = function(p, eb, lb0=NA){ + # created 2007/08/15 by Bas Kooijman; modified 2013/08/19, 2015/01/20 + + ## Syntax + # [lb, info] = <../get_lb.m *get_lb*>(p, eb, lb0) + + ## Description + # Obtains scaled length at birth, given the scaled reserve density at birth. + # + # Input + # + # * p: 3-vector with parameters: g, k, v_H^b (see below) + # * eb: optional scalar with scaled reserve density at birth (default eb = 1) + # * lb0: optional scalar with initial estimate for scaled length at birth (default lb0: lb for k = 1) + # + # Output + # + # * lb: scalar with scaled length at birth + # * info: indicator equals 1 if successful, 0 otherwise + + ## Remarks + # The theory behind get_lb, get_tb and get_ue0 is discussed in + # . + # Solves y(x_b) = y_b for lb with explicit solution for y(x) + # y(x) = x e_H/(1-kap) = x g u_H/ l^3 + # and y_b = x_b g u_H^b/ ((1-kap)l_b^3) + # d/dx y = r(x) - y s(x); + # with solution y(x) = v(x) \int r(x)/ v(x) dx + # and v(x) = exp(- \int s(x) dx). + # A Newton Raphson scheme is used with Euler integration, starting from an optional initial value. + # Replacement of Euler integration by ode23: , + # but that function is much lower. + # Shooting method: . + # Bisection method (via fzero): . + # In case of no convergence, is run automatically as backup. + # Consider the application of for an alternative initial value. + + ## Example of use + # See <../mydata_ue0.m *mydata_ue0*> + + # unpack p + g = p[1] # g = [E_G] * v/ kap * {p_Am}, energy investment ratio + k = p[2] # k = k_J/ k_M, ratio of maturity and somatic maintenance rate coeff + vHb = p[3] # v_H^b = U_H^b g^2 kM^3/ (1 - kap) v^2; U_H^b = M_H^b/ {J_EAm} + + lb0=Re(lb0) + + info = 1 + if (!exists('lb0')){ + lb0 = NA + } + + if (k == 1){ + lb = as.complex(vHb)^(1/ 3) # exact solution for k = 1 + info = 1 + } + + if (is.na(lb0)){ + lb = as.complex(vHb)^(1/3) # exact solution for k = 1 + } else {lb = lb0} + + if (!exists('eb')){ + eb = 1 + } else if (is.na(eb)){ + eb = 1 + } + + n = 1000 + round(1000 * max(0, k - 1)) + xb = g/ (g + eb) + xb3 = xb ^ (1/3) + x = seq(1e-6, xb, length=n) + dx = xb/ n + x3 = x ^ (1/3) + + b = beta0(x, xb)/ (3 * g) + t0 = xb * g * vHb + i = 0 + norm = 1 # make sure that we start Newton Raphson procedure + ni = 100 # max number of iterations + + while (i < ni && norm > 1e-8 && !is.na(norm)){ + l = x3 / (xb3/ lb - b) + s = (k - x) / (1 - x) * l/ g / x + v = exp( - dx * cumsum(s)) + vb = v[n] + r = (g + l) + rv = r / v + t = t0/ lb^3/ vb - dx * sum(rv) + dl = xb3/ lb^2 * l ^ 2 / x3 + dlnl = dl / l + dv = v * exp( - dx * cumsum(s * dlnl)) + dvb = dv[n] + dlnv = dv / v + dlnvb = dlnv[n] + dr = dl + dlnr = dr / r + dt = - t0/ lb^3/ vb * (3/ lb + dlnvb) - dx * sum((dlnr - dlnv) * rv) + # [i lb t dt] # print progress + lb = Re(lb - t/ dt) # Newton Raphson step + norm = Re( t^2 ) + i = i + 1 + } + + if (i == ni || lb < 0 || lb > 1 || is.na(norm)) { # no convergence + # try to recover with a shooting method + if (is.na(lb0)){ + lbinfo = get_lb2(p, eb) + lb=lbinfo[1] + info=lbinfo[2] + } else if (lb0 < 1 && lb0 > 0){ + lbinfo = get_lb2(p, eb, lb0) + lb=lbinfo[1] + info=lbinfo[2] + } else { + lbinfo = get_lb2(p, eb) + lb=lbinfo[1] + info=lbinfo[2] + } + } + + if (info == 0){ + print('warning get_lb: no convergence of l_b') + } + + return(c(lb, info)) + +} \ No newline at end of file diff --git a/R/get_lb2.R b/R/get_lb2.R new file mode 100644 index 0000000..7b8b0a2 --- /dev/null +++ b/R/get_lb2.R @@ -0,0 +1,65 @@ +## get_lb2 +# Obtains scaled length at birth, given the scaled reserve density at birth + +## +get_lb2= function(p, eb, lb0=NA){ +# created 2007/07/26 by Bas Kooijman; modified 2013/08/19, 2015/01/18 + +## Syntax +# [lb info] = <../get_lb2.m *get_lb2*> (p, eb, lb0) + +## Description +# Obtains scaled length at birth, given the scaled reserve density at birth. +# Like get_lb, but using the shooting method, rather than Newton Raphson +# +# Input +# +# * p: 3-vector with parameters: g, k, v_H^b (see below) +# * eb: optional scalar with scaled reserve density at birth (default eb = 1) +# * lb0: optional scalar with initial estimate for scaled length at birth (default lb0: lb for k = 1) +# +# Output +# +# * lb: scalar with scaled length at birth +# * info: indicator equals 1 if successful, 0 otherwise + +## Remarks +# Like , but uses a shooting method in 1 variable + +# unpack p + g = p[1] # g = [E_G] * v/ kap * {p_Am}, energy investment ratio + k = p[2] # k = k_J/ k_M, ratio of maturity and somatic maintenance rate coeff + vHb = p[3] # v_H^b = U_H^b g^2 kM^3/ (1 - kap) v^2; U_H^b = M_H^b/ {J_EAm} + + info = 1 + + + if (is.na(lb0)) { + lb = as.complex(vHb)^(1/ 3) # exact solution for k = 1 + } else { + lb = lb0 + } + if (!exists('eb')){ + eb = 1 + } else if (is.na(eb)){ + eb = 1 + } + + + xb = g/ (eb + g) + xb3 = xb^(1/3) + + xbxb3gvHbk=c(xb, xb3, g, vHb, k) + lb=Re(lb) + + lbflaginfo = fzero(fnget_lb2, lb, maxiter = 100, tol = 10e-8, xbxb3gvHbk) + #lbflaginfo = uniroot(fnget_lb2, c(-10, 10), xbxb3gvHbk) + lb=lbflaginfo$x + info=1 + + if (lb < 0 || lb > 1){ + info = 0 + } + return(c(lb, info)) +} + diff --git a/R/get_lj.R b/R/get_lj.R new file mode 100644 index 0000000..9b36b7d --- /dev/null +++ b/R/get_lj.R @@ -0,0 +1,125 @@ +## get_lj +# Gets scaled length at metamorphosis + +## + #function [lj, lp, lb, info] = + +get_lj=function(p, f, lb0) { + + # created at 2010/02/10 by Bas Kooijman, + # modified 2014/03/03 Starrlight Augustine, 2015/01/18 Bas Kooijman + + ## Syntax + # [lj, lp, lb, info] = <../get_lj.m *get_lj*>(p, f, lb0) + + ## Description + # Type M-acceleration: Isomorph, but V1-morph between vHb and vHj + # This routine obtaines scaled length at metamorphosis lj given scaled muturity at metamorphosis vHj. + # The theory behind get_lj, is discussed in the comments to DEB3. + # If scaled length at birth (third input) is not specified, it is computed (using automatic initial estimate); + # if it is specified. however, is it just copied to the (third) output. + # The code assumes vHb < vHj < vHp (see first input). + # + # Input + # + # * p: 6-vector with parameters: g, k, l_T, v_H^b, v_H^j, v_H^p + # + # if p is a 5-vector, output lp is empty + # + # * f: optional scalar with scaled functional responses (default 1) + # * lb0: optional scalar with scaled length at birth + # + # Output + # + # * lj: scalar with scaled length at metamorphosis + # * lp: scalar with scaled length at puberty + # * lb: scalar with scaled length at birth + # * info: indicator equals 1 if successful + + ## Remarks + # Similar to , which uses root finding, rather than integration + # Scaled length l = L/ L_m with L_m = kap {p_Am}/ [p_M] where {p_Am} is of embryo, + # because the amount of acceleration is food-dependent so the value after metamorphosis is not a parameter. + # {p_Am} and v increase between maturities v_H^b and v_H^j till + # {p_Am} l_j/ l_b and v l_j/ l_b at metamorphosis. + # After metamorphosis l increases from l_j till l_i = f l_j/ l_b - l_T. + # Scaled length l can thus be larger than 1. + # See for foetal development. + + ## Example of use + # get_lj([.5, .1, .1, .01, .2, .3]) + + n_p = length(p) + + # unpack pars + g = p[1] # energy investment ratio + k = p[2] # k_J/ k_M, ratio of maturity and somatic maintenance rate coeff + lT = p[3] # scaled heating length {p_T}/[p_M] + vHb = p[4] # v_H^b = U_H^b g^2 kM^3/ (1 - kap) v^2; U_H^b = M_H^b/ {J_EAm} = E_H^b/ {p_Am} + vHj = p[5] # v_H^j = U_H^j g^2 kM^3/ (1 - kap) v^2; U_H^j = M_H^j/ {J_EAm} = E_H^j/ {p_Am} + if (n_p > 5){ + vHp = p[6] # v_H^p = U_H^j g^2 kM^3/ (1 - kap) v^2; U_B^j = M_H^j/ {J_EAm} + } + + if (!exists('f')) { + f = 1 + } else if (is.na(f)){ + f = 1 + } + + # get lb + if (!exists('lb0')){ + lb0 = NA + } + if (is.na(lb0)){ + lbinfo = get_lb(c(g, k, vHb), f, lb0) + lb=lbinfo[1] + info=lbinfo[2] + } else { + info = 1 + lb = lb0 + } + + if (lT > f - lb) { # no growth is possible at birth + lj = NA + lp = NA + info = 0 + } + + # no acceleration + if (vHb == vHj){ + if (n_p == 5){ # no puberty specified + lbinfo = get_lb(p[c(1, 2, 4)], f, lb0) + lb=lbinfo[1] + info=lbinfo[2] + lj = lb + lp = NA + } # else { # puberty specified +# [lp, lb, info] = get_tp(p([1 2 3 4 6]), f, lb0); +# lj = lb; +# } + } + + + # get lj + library(deSolve) + tspan=seq(vHb, vHj, length=100) + lini=c(l=lb) + out = ode(lini, tspan, dget_l_V1, parms=c(k, lT, g, f, lb), method="ode45") + vH=out[,1] + lj=out[length(out[,2]),2] + sM = lj/ lb + + # get lp + if (n_p > 5){ + tspan=seq(vHj, vHp, length=100) + lini=lj + out = ode(lini, tspan, dget_l_ISO, parms=c(k, lT, g, f, sM), method="ode45") + vH=out[,1] + lp=out[length(out[,2]),2] + } else lp = NA + + +return(c(lj, lp, lb, info)) + +} \ No newline at end of file diff --git a/R/get_tb.R b/R/get_tb.R new file mode 100644 index 0000000..f016391 --- /dev/null +++ b/R/get_tb.R @@ -0,0 +1,79 @@ +## get_tb +# + +## +get_tb= function(p, eb, lb=NA){ +# created at 2007/07/27 by Bas Kooijman; modified 2014/03/17, 2015/01/18 + +## Syntax +# [tb, lb, info] = <../get_tb.m *get_tb*> (p, eb, lb) + +## Description +# Obtains scaled age at birth, given the scaled reserve density at birth. +# Divide the result by the somatic maintenance rate coefficient to arrive at age at birth. +# +# Input +# +# * p: 1 or 3-vector with parameters g, k_J/ k_M, v_H^b +# +# Last 2 values are optional in invoke call to get_lb +# +# * eb: optional scalar with scaled reserve density at birth (default eb = 1) +# * lb: optional scalar with scaled length at birth (default: lb is obtained from get_lb) +# +# Output +# +# * tb: scaled age at birth \tau_b = a_b k_M +# * lb: scalar with scaled length at birth: L_b/ L_m +# * info: indicator equals 1 if successful, 0 otherwise + +## Remarks +# See also for backward integration over maturity and +# for foetal development + +## Example of use +# get_tb([.1;.5;.03]) +# See also <../ mydata_ue0.m *mydata_ue0*> + + if (!exists('eb')) { + eb = 1 # maximum value as juvenile + } + + info = 1 + + if (!exists('lb')) { + if (length(p) < 3){ + print('not enough input parameters, see get_lb \n') + tb = NA + } + lbinfo = get_lb(p, eb) + lb=lbinfo[1] + info=lbinfo[2] + } + if (is.na(lb)){ + lbinfo = get_lb(p, eb) + lb=lbinfo[1] + info=lbinfo[2] + } + + # unpack p + g = p[1] # energy investment ratio + + xb = g/ (eb + g) # f = e_b + ab = 3 * g * xb^(1/ 3)/ lb # \alpha_b + library(pracma) + abxb=c(ab,xb) + tb = 3 * quad(dget_tb, xa= 1e-15, xb=xb, tol = 1.0e-12, trace = FALSE, abxb) + return(c(tb, lb, info)) +} + + # subfunction + +dget_tb= function(x, abxb){ + # called by get_tb + ab=abxb[1] + xb=abxb[2] + f = x ^ (-2/ 3) / (1 - x) / (ab - beta0(x, xb)) + return(f) +} + diff --git a/R/get_tj.R b/R/get_tj.R new file mode 100644 index 0000000..a6e3e7c --- /dev/null +++ b/R/get_tj.R @@ -0,0 +1,173 @@ +## get_tj +# Gets scaled age at metamorphosis + +## + # function [tj, tp, tb, lj, lp, lb, li, rj, rB, info] = + +get_tj = function(p, f, lb0=NA){ + # created at 2011/04/25 by Bas Kooijman, + # modified 2014/03/03 Starrlight Augustine, 2015/01/18 Bas Kooijman + + ## Syntax + # [tj, tp, tb, lj, lp, lb, li, rj, rB, info] = <../get_tj.m *get_tj*> (p, f, lb0) + + ## Description + # Obtains scaled ages at metamorphosis, puberty, birth and the scaled lengths at these ages; + # Food density is assumed to be constant. + # Multiply the result with the somatic maintenance rate coefficient to arrive at unscaled ages. + # Metabolic acceleration occurs between birth and metamorphosis, see also get_ts. + # Notice j-p-b sequence in output, due to the name of the routine + # + # Input + # + # * p: 5 or 6-vector with parameters: g, k, l_T, v_H^b, v_H^j v_H^p + # + # Last value is optional. If ommitted: outputs tp and lp are empty + # + # * f: optional scalar with functional response (default f = 1) + # * lb0: optional scalar with scaled length at birth + # or optional 2-vector with scaled length, l, and scaled maturity, vH + # for a juvenile that is now exposed to f, but previously at another f + # + # Output + # + # * tj: scaled with age at metamorphosis \tau_j = a_j k_M + # + # if length(lb0)==2, tj is the scaled time till metamorphosis + # + # * tp: scaled with age at puberty \tau_p = a_p k_M + # + # if length(lb0)==2, tp is the scaled time till puberty + # + # * tb: scaled with age at birth \tau_b = a_b k_M + # * lj: scaled length at end of V1-stage + # * lp: scaled length at puberty + # * lb: scaled length at birth + # * li: ultimate scaled length + # * rj: scaled exponential growth rate between s and j + # * rB: scaled von Bertalanffy growth rate between b and s and between j and i + # * info: indicator equals 1 if successful, 0 otherwise + + ## Remarks + # See in case of foetal development + + ## Example of use + # get_tj([.5, .1, 0, .01, .05, .2]) + + n_p = length(p) + + # unpack pars + g = p[1] # energy investment ratio + k = p[2] # k_J/ k_M, ratio of maturity and somatic maintenance rate coeff + lT = p[3] # scaled heating length {p_T}/[p_M]Lm + vHb = p[4] # v_H^b = U_H^b g^2 kM^3/ (1 - kap) v^2; U_B^b = M_H^b/ {J_EAm} + vHj = p[5] # v_H^j = U_H^j g^2 kM^3/ (1 - kap) v^2; U_B^j = M_H^j/ {J_EAm} + if (n_p > 5){ + vHp = p[6]} else vHp = 0 # v_H^p = U_H^p g^2 kM^3/ (1 - kap) v^2; U_B^p = M_H^p/ {J_EAm} + + + if (!exists('f')) { + f = 1 + } + if (is.na(f)){ + f = 1 + } + if (!exists('lb0')){ + lb0 = NA + } + + # no acceleration + if (vHb == vHj){ + if (vHp == 0){ # no puberty specified + tblbinfo=get_tb(p[c(1:3)], f, lb0) + tj = tblbinfo[1] + tp = NA + lp = NA + lj = tblbinfo[2] + li = f - lT + rj = 0 + rB = 1/ 3/ (1 + f/ g) + } # else { # puberty specified +# [tp, tb, lp, lb, info] = get_tp(p([1 2 3 4 6]), f, lb0); +# tj = tb; lj = lb; li = f - lT; rj = 0; rB = 1/ 3/ (1 + f/ g); +# } + } + + + # maintenance ratio k = 1: maturity thresholds coincide with length thresholds + if (k == 1 & f * (f - lT)^2 > vHp * k){ # constant maturity density, reprod possible + lb = vHb^(1/3) # scaled length at birth + tblbinfo = get_tb(p[c(1,2,4)], f, lb) # scaled age at birth + tb = tblbinfo[1] + lj = vHj^(1/3) # scaled length at metamorphosis + sM = lj/ lb # acceleration factor + rj = g * (f/ lb - 1 - lT/ lb)/ (f + g) # scaled exponential growth rate between b and j + tj = tb + (log(sM)) * 3/ rj # scaled age at metamorphosis + lp = vHp^(1/3) # scaled length at puberty + li = f * sM - lT # scaled ultimate length + rB = 1/ 3/ (1 + f/ g) # scaled von Bert growth rate between j and i + tp = tj + (log ((li - lj)/ (li - lp)))/ rB # scaled age at puberty + info = 1 + } else if (k == 1 & f * (f - lT)^2 > vHj * k) { # constant maturity density, metam possible + lb = vHb^(1/3) # scaled length at birth + tblbinfo = get_tb(p[c(1,2,4)], f, lb) # scaled age at birth + tb = tblbinfo[1] + lj = vHj^(1/3) # scaled length at metamorphosis + sM = lj/ lb # acceleration factor + rj = g * (f/ lb - 1 - lT/ lb)/ (f + g) # scaled exponential growth rate between b and j + tj = tb + (log(sM)) * 3/ rj # scaled age at metamorphosis + lp = vHp^(1/3) # scaled length at puberty + li = f * sM - lT # scaled ultimate length + rB = 1/ 3/ (1 + f/ g) # scaled von Bert growth rate between j and i + tp = 1e20 # scaled age at puberty + info = 1 + } + + + if (is.na(lb0)){ + tblbinfo = get_tb (p[c(1, 2, 4)], f, lb0) + } else {tblbinfo = get_tb (p[C(1, 2, 4)], f, lb0[1]) } + + tb=tblbinfo[1] + lb=tblbinfo[2] + info_tb=tblbinfo[3] + + ljlplbinfo = get_lj(p, f, lb) + lj=ljlplbinfo[1] + lp=ljlplbinfo[2] + lb=ljlplbinfo[3] + info_tj=ljlplbinfo[3] + + sM = lj/lb # acceleration factor + rj = g * (f/ lb - 1 - lT/ lb)/ (f + g) # scaled exponential growth rate between b and j + tj = tb + log(sM) * 3/ rj # scaled age at metamorphosis + rB = 1/ 3/ (1 + f/ g) # scaled von Bert growth rate between j and i + li = f * sM - lT # scaled ultimate length + + if (is.na(lp)) { # length(p) < 6 + tp = NA + } else if (li <= lp) { # reproduction is not possible + tp = 1e20; # tau_p is never reached + lp = 1; # lp is nerver reached + } else { # reproduction is possible + if (length(lb0) != 2) { # lb0 is absent, empty or a scalar + tp = tj + log((li - lj)/ (li - lp))/ rB + } else{ # lb0 = l and t for a juvenile + tb = NaN + l = lb0[1] + tp = log((li - l)/ (li - lp))/ rB; + } + } + + if (is.na(tp)) info = info_tb + + info = min(info_tb, info_tj) + + if (!is.finite(tp) || !is.finite(tj)){ # tj and tp must be real and positive + info = 0 + } else if (tp < 0 || tj < 0){ + info = 0 + } + + return(c(tj, tp, tb, lj, lp, lb, li, rj, rB, info)) +} diff --git a/R/get_tm.R b/R/get_tm.R new file mode 100644 index 0000000..8d0e4e8 --- /dev/null +++ b/R/get_tm.R @@ -0,0 +1,151 @@ +## get_tm_s +# Obtains scaled mean age at death fro short growth periods + +## + #function [tm, Sb, Sp, info] = + +get_tm_s = function(p, f, lb, lp){ + + # created 2009/02/21 by Bas Kooijman, modified 2014/03/17, 2015/01/18 + + ## Syntax + # [tm, Sb, Sp, info] = <../get_tm_s.m *get_tm_s*>(p, f, lb, lp) + + ## Description + # Obtains scaled mean age at death assuming a short growth period relative to the life span + # Divide the result by the somatic maintenance rate coefficient to arrive at the mean age at death. + # The variant get_tm_foetus does the same in case of foetal development. + # If the input parameter vector has only 4 elements (for [g, lT, ha/ kM2, sG]), + # it skips the calulation of the survival probability at birth and puberty. + # + # Input + # + # * p: 4 or 7-vector with parameters: [g lT ha sG] or [g k lT vHb vHp ha SG] + # * f: optional scalar with scaled reserve density at birth (default eb = 1) + # * lb: optional scalar with scaled length at birth (default: lb is obtained from get_lb) + # * lp: optional scalar with scaled length at puberty + # + # Output + # + # * tm: scalar with scaled mean life span + # * Sb: scalar with survival probability at birth (if length p = 7) + # * Sp: scalar with survival prabability at puberty (if length p = 7) + # * info: indicator equals 1 if successful, 0 otherwise + + ## Remarks + # Theory is given in comments on DEB3 Section 6.1.1. + # See for the general case of long growth period relative to life span + + ## Example of use + # get_tm_s([.5, .1, .1, .01, .2, .1, .01]) + + if (!exists('f')){ + f = 1 + } + else if (is.na(f)){ + f = 1 + } + + if (length(p) >= 7) { + # unpack pars + g = p[1] # energy investment ratio + #k = p[2] # k_J/ k_M, ratio of maturity and somatic maintenance rate coeff + lT = p[3] # scaled heating length {p_T}/[p_M]Lm + #vHb = p[4] # v_H^b = U_H^b g^2 kM^3/ (1 - kap) v^2; U_B^b = M_H^b/ {J_EAm} + #vHp = p[5] # v_H^p = U_H^p g^2 kM^3/ (1 - kap) v^2; U_B^p = M_H^p/ {J_EAm} + ha = p[6] # h_a/ k_M^2, scaled Weibull aging acceleration + sG = p[7] # Gompertz stress coefficient + } + else if (length(p) == 4){ + # unpack pars + g = p[1] # energy investment ratio + lT = p[2] # scaled heating length {p_T}/[p_M]Lm + ha = p[3] # h_a/ k_M^2, scaled Weibull aging acceleration + sG = p[4] # Gompertz stress coefficient + } + + if (abs(sG) < 1e-10){ + sG = 1e-10 + } + + li = f - lT + hW3 = ha * f * g/ 6/ li + hW = hW3^(1/3) # scaled Weibull aging rate + hG = sG * f * g * li^2 + hG3 = hG^3 # scaled Gompertz aging rate + tG = hG/ hW + tG3 = hG3/ hW3 # scaled Gompertz aging rate + + info_lp = 1 + if (length(p) >= 7){ +# if (!exist(lp) & !exist(lb)) { +# [lp lb info_lp] = get_lp(p, f) # !!!!!!!!!!!!!!!!!!!!get_lp? +# } +# elseif exist('lp','var') == 0 || isempty('lp') +# [lp lb info_lp] = get_lp(p, f, lb); +# end +# + # get scaled age at birth, puberty: tb, tp + tblbinfo = get_tb(p, f, lb) + tb=tblbinfo[1] + lb=tblbinfo[2] + info_lp=tblbinfo[3] + irB = 3 * (1 + f/ g) + tp = tb + irB * log((li - lb)/ (li - lp)) + hGtb = hG * tb + Sb = exp((1 - exp(hGtb) + hGtb + hGtb^2/ 2) * 6/ tG3) + hGtp = hG * tp + Sp = exp((1 - exp(hGtp) + hGtp + hGtp^2/ 2) * 6/ tG3) + if (info_lp == 1 & info_tb == 1){ + info = 1 + } + else info = 0 + } + else { # length(p) == 4 + Sb = NaN + Sp = NaN + info = 1 + } + + if (abs(sG) < 1e-10) { + tm = gamma(4/3)/ hW + tm_tail = 0 + } + else if (hG > 0){ + tm = 10/ hG # upper boundary for main integration of S(t) + library(pracma) + tm_tail = expint(exp(tm * hG) * 6/ tG3)/ hG + tm = quad(fnget_tm_s, 0, tm * hW, tol = .Machine$double.eps^0.5, trace = FALSE, tG)/ hW + tm_tail + } + else {# hG < 0 + tm = -1e4/ hG # upper boundary for main integration of S(t) + hw = hW * sqrt( - 3/ tG) # scaled hW + library(pracma) + tm_tail = sqrt(pi)* erfc(tm * hw)/ 2/ hw + tm = quad(fnget_tm_s, 0, tm * hW, tol = .Machine$double.eps^0.5, trace = FALSE, tG)/ hW + tm_tail + } + + +return(c(tm, Sb, Sp, info)) + +} + +# subfunction + +fnget_tm_s = function(t, tG){ + # modified 2010/02/25 + # called by get_tm_s for life span at short growth periods + # integrate ageing surv prob over scaled age + # t: age * hW + # S: ageing survival prob + + hGt = tG * t # age * hG + if (tG > 0) { + S = exp(-t(1 + sum(cumprod(t(hGt) * (1/(4:500)), 2), 2)) * t^3) + } + else { # tG < 0 + S = exp((1 + hGt + hGt^2/ 2 - exp(hGt)) * 6/ tG^3) + } + + return(S) +} \ No newline at end of file diff --git a/R/get_tm_s.R b/R/get_tm_s.R new file mode 100644 index 0000000..e7543a4 --- /dev/null +++ b/R/get_tm_s.R @@ -0,0 +1,147 @@ +## get_tm_s +# Obtains scaled mean age at death fro short growth periods + +## +# function [tm, Sb, Sp, info] = +get_tm_s = function(p, f, lb, lp) { + # created 2009/02/21 by Bas Kooijman, modified 2014/03/17, 2015/01/18 + + ## Syntax + # [tm, Sb, Sp, info] = <../get_tm_s.m *get_tm_s*>(p, f, lb, lp) + + ## Description + # Obtains scaled mean age at death assuming a short growth period relative to the life span + # Divide the result by the somatic maintenance rate coefficient to arrive at the mean age at death. + # The variant get_tm_foetus does the same in case of foetal development. + # If the input parameter vector has only 4 elements (for [g, lT, ha/ kM2, sG]), + # it skips the calulation of the survival probability at birth and puberty. + # + # Input + # + # * p: 4 or 7-vector with parameters: [g lT ha sG] or [g k lT vHb vHp ha SG] + # * f: optional scalar with scaled reserve density at birth (default eb = 1) + # * lb: optional scalar with scaled length at birth (default: lb is obtained from get_lb) + # * lp: optional scalar with scaled length at puberty + # + # Output + # + # * tm: scalar with scaled mean life span + # * Sb: scalar with survival probability at birth (if length p = 7) + # * Sp: scalar with survival prabability at puberty (if length p = 7) + # * info: indicator equals 1 if successful, 0 otherwise + + ## Remarks + # Theory is given in comments on DEB3 Section 6.1.1. + # See for the general case of long growth period relative to life span + + ## Example of use + # get_tm_s([.5, .1, .1, .01, .2, .1, .01]) + + if (!exists('f')) { + f = 1 + } else if (is.na(f)) { + f = 1; + } + + if (length(p) >= 7){ + # unpack pars + g = p[1] # energy investment ratio + #k = p[2] # k_J/ k_M, ratio of maturity and somatic maintenance rate coeff + lT = p[3] # scaled heating length {p_T}/[p_M]Lm + #vHb = p[4] # v_H^b = U_H^b g^2 kM^3/ (1 - kap) v^2; U_B^b = M_H^b/ {J_EAm} + #vHp = p[5] # v_H^p = U_H^p g^2 kM^3/ (1 - kap) v^2; U_B^p = M_H^p/ {J_EAm} + ha = p[6] # h_a/ k_M^2, scaled Weibull aging acceleration + sG = p[7] # Gompertz stress coefficient + } else if (length(p) == 4){ + # unpack pars + g = p[1] # energy investment ratio + lT = p[2] # scaled heating length {p_T}/[p_M]Lm + ha = p[3] # h_a/ k_M^2, scaled Weibull aging acceleration + sG = p[4] # Gompertz stress coefficient + } + + if (abs(sG) < 1e-10) { + sG = 1e-10 + } + + li = f - lT + hW3 = ha * f * g/ 6/ li + hW = hW3^(1/3) # scaled Weibull aging rate + hG = sG * f * g * li^2 + hG3 = hG^3 # scaled Gompertz aging rate + tG = hG/ hW + tG3 = hG3/ hW3 # scaled Gompertz aging rate + + info_lp = 1 + if (length(p) >= 7) { + if (!exists('lp') & !exists('lb')) { + lplbinfo = get_lp(p, f) + lp=lplbinfo[1] + lb=lplbinfo[2] + info_lp=lplbinfo[3] + } else if (!exist('lp') || is.na('lp')){ + lplbinfo = get_lp(p, f, lb) + lp=lplbinfo[1] + lb=lplbinfo[2] + info_lp=lplbinfo[3] + } + + + # get scaled age at birth, puberty: tb, tp + tblbinfo_tb = get_tb(p, f, lb) + tb=tblbinfo_tb[1] + lb=tblbinfo_tb[2] + info_tb=tblbinfo_tb[3] + + irB = 3 * (1 + f/ g) + tp = tb + irB * log((li - lb)/ (li - lp)) + hGtb = hG * tb + Sb = exp((1 - exp(hGtb) + hGtb + hGtb^2/ 2) * 6/ tG3) + hGtp = hG * tp + Sp = exp((1 - exp(hGtp) + hGtp + hGtp^2/ 2) * 6/ tG3) + if (info_lp == 1 && info_tb == 1){ + info = 1 + } else{ + info = 0 + } + } else { # length(p) == 4 + Sb = NaN + Sp = NaN + info = 1 + } + + if (abs(sG) < 1e-10) { + tm = gamma(4/3)/ hW + tm_tail = 0 + } else if (hG > 0) { + tm = 10/ hG # upper boundary for main integration of S(t) + library(pracma) + tm_tail = expint(exp(tm * hG) * 6/ tG3)/ hG + tm = quad(fnget_tm_s, 0, tm * hW, tol = 1.0e-12, trace = FALSE,tG)/ hW + tm_tail + } else {# hG < 0 + tm = -1e4/ hG # upper boundary for main integration of S(t) + hw = hW * sqrt( - 3/ tG) # scaled hW + tm_tail = sqrt(pi)* erfc(tm * hw)/ 2/ hw + tm = quad(fnget_tm_s, 0, tm * hW, tol = 1.0e-12, trace = FALSE, tG)/ hW + tm_tail + } +} + + # subfunction + + #function S = +fnget_tm_s=function(t, tG){ + # modified 2010/02/25 + # called by get_tm_s for life span at short growth periods + # integrate ageing surv prob over scaled age + # t: age * hW + # S: ageing survival prob + + hGt = tG * t # age * hG + if (tG > 0) { + a= matrix(hGt) %*% t(matrix((1/(4:500)))) + S = exp(-(1 + apply((t(apply(a, 1, cumprod))),1, sum)) * t^3) #apply((t(apply(a, 1, cumprod))),2, sum)? + } else {# tG < 0 + S = exp((1 + hGt + hGt^2/ 2 - exp(hGt)) * 6/ tG^3) + } + return(S) +} \ No newline at end of file diff --git a/R/get_ue0.R b/R/get_ue0.R new file mode 100644 index 0000000..bf8e4ee --- /dev/null +++ b/R/get_ue0.R @@ -0,0 +1,67 @@ +## get_ue0 +# gets initial scaled reserve + +## +# function [uE0, lb, info] = +get_ue0=function(p, eb, lb0){ + # created at 2007/07/27 by Bas Kooijman; modified 2010/05/02 + + ## Syntax + # [uE0, lb, info] = <../get_ue0.m *get_ue0*>(p, eb, lb0) + + ## Description + # Obtains the initial scaled reserve given the scaled reserve density at birth. + # Function get_ue0 does so for eggs, get_ue0_foetus for foetuses. + # Specification of length at birth as third input by-passes its computation, + # so if you want to specify an initial value for this quantity, you should use get_lb directly. + # + # Input + # + # * p: 1 or 3 -vector with parameters g, k_J/ k_M, v_H^b, see get_lb + # * eb: optional scalar with scaled reserbe density at birth + # (default: eb = 1) + # * lb0: optional scalar with scaled length at birth + # (default: lb is optained from get_lb) + # + # Output + # + # * uE0: scaled with scaled reserve at t=0: $U_E^0 g^2 k_M^3/ v^2$ + # with $U_E^0 = M_E^0/ \{J_{EAm}\}$ + # * lb: scalar with scaled length at birth + # * info: indicator equals 1 if successful, 0 otherwise + + ## Remarks + # See for foetal development. + # See , for a non-dimensionless scaling. + + ## Example of use + # see <../mydata_ue0.m *mydata_ue0*> + + if (!exists('eb')) { + eb = 1 # maximum value as juvenile + } + + if (!exists('lb0')){ + if (length(p) < 3) { + print('not enough input parameters, see get_lb \n'); + uE0 = NA + lb = NA + info = 0 + } + lbinfo = get_lb(p, eb) + lb=lbinfo[1] + info=lbinfo[2] + } else { + lb = lb0 + info = 1 + } + + + # unpack p + g = p[1] # energy investment ratio + + xb = g/ (eb + g) + uE0 = (3 * g/ (3 * g * xb^(1/ 3)/ lb - beta0(0, xb)))^3 + + return(c(uE0, lb, info)) +} diff --git a/R/initial_scaled_reserve.R b/R/initial_scaled_reserve.R new file mode 100644 index 0000000..a4d7e48 --- /dev/null +++ b/R/initial_scaled_reserve.R @@ -0,0 +1,67 @@ +## initial_scaled_reserve +# Gets initial scaled reserve + +## +# function [U0, Lb, info] + +initial_scaled_reserve=function(f, p, Lb0){ + + # created 2007/08/06 by Bas Kooyman; modified 2009/09/29 + + ## Syntax + # [U0, Lb, info] = <../initial_scaled_reserve.m *initial_scaled_reserve*>(f, p, Lb0) + + ## Description + # Gets initial scaled reserve + # + # Input + # + # * f: n-vector with scaled functional responses + # * p: 5-vector with parameters: VHb, g, kJ, kM, v + # * Lb0: optional n-vector with lengths at birth + # + # Output + # + # * U0: n-vector with initial scaled reserve: M_E^0/ {J_EAm} + # * Lb: n-vector with length at birth + # * info: n-vector with 1's if successful, 0's otherwise + + ## Remarks + # Like , but allows for vector arguments and + # input and output is not downscaled to dimensionless quantities, + + ## Example of use + # p = [.8 .42 1.7 1.7 3.24 .012]; initial_scaled_reserve(1,p) + + # unpack parameters + VHb = p[1] # d mm^2, scaled maturity at birth: M_H^b/[[1-kap]{J_EAm}] + g = p[2] # -, energy investment ratio + kJ = p[3] # 1/d, maturity maintenance rate coefficient + kM = p[4] # 1/d, somatic maintenance rate coefficient + v = p[5] # mm/d, energy conductance + + # if kJ = kM: VHb = g * Lb^3/ v; + + nf = length(f) + U0 = rep(0, length=nf) + Lb = rep(0, length=nf) + info = rep(0, length=nf) + q = c(g, kJ/ kM, VHb * g^2 * kM^3/ v^2) + if (exists('Lb0')){ + lb0 = rep(1, length=nf) * Lb0 * kM * g/ v + } else { + lb0 = ones(nf,1) * get_lb(q,f[1])[1] # initial estimate for scaled length + } + for (i in c(1:nf)){ + lbinfo = get_lb(q, f[i], lb0[i]) + lb=lbinfo[1] + info[i]=lbinfo[2] + ## try get_lb1 or get_lb2 for higher accuracy + Lb[i] = lb * v/ kM/ g; + uE0lbinfo = get_ue0(q, f[i], lb) + uE0=uE0lbinfo[1] + U0[i] = uE0 * v^2/ g^2/ kM^3 + } + + return(c(U0,Lb, info)) +} diff --git a/R/nmregr.R b/R/nmregr.R new file mode 100644 index 0000000..80f9a93 --- /dev/null +++ b/R/nmregr.R @@ -0,0 +1,337 @@ +# nmregr +# Calculates least squares estimates using Nelder Mead's simplex method + +## +#function [q, info] = + +nmregr=function(func, p, varargin){ + + # created 2001/09/07 by Bas Kooijman; modified 2013/10/04 + + ## Syntax + # [q, info] = <../nmregr.m *nmregr*> (func, p, varargin) + + ## Description + # Calculates least squares estimates using Nelder Mead's simplex method + # + # Input + # + # * func: string with name of user-defined function + # + # f = func (p, xyw) with + # p: k-vector with parameters; xyw: (n,c)-matrix; f: n-vector + # [f1, f2, ...] = func (p, xyw1, xyw2, ...) with p: k-vector and + # xywi: (ni,k)-matrix; fi: ni-vector with model predictions + # The dependent variable in the output f; For xyw see below. + # + # * p: (k,2)-matrix with + # + # p(:,1) initial guesses for parameter values + # p(:,2) binaries with yes or no iteration (optional) + # + # * xyzi (read as xyw1, xyw2, .. ): (ni,3) matrix with + # + # xywi(:,1) independent variable i + # xywi(:,2) dependent variable i + # xywi(:,3) weight coefficients i (optional) + # xywi(:,>3) data-pont specific information data (optional) + # The number of data matrices xyw1, xyw2, ... is optional but >0 + # Default for xywi(:,3) is (number of data points in the set i)^-1 + # + # Output + # + # * q: matrix like p, but with least squares estimates + # * info: 1 if convergence has been successful; 0 otherwise + + ## Remarks + # Calls user-defined function 'func' + # Set options with + # See + # for Newton-Raphson method, + # for genetic algorithm, + # for 2 independent variables, and + # for standard deviations proportional to the mean. + # It is usually a good idea to run on the result of nmregr. + + ## Example of use + # See <../mydata_regr.m *mydata_regr*> + + func=match.fun(func) + + i = 1 # initiate data set counter + info = 1 # initiate info setting + ci = as.character(i) # character string with value of i + if (!is.list(varargin)){ + varargin=list(varargin) + } + + nxyw = nrow(summary(varargin)) # number of data sets + + while (i <= nxyw){ # loop across data sets + if (i == 1){ + listxyw = c('xyw1') # initiate list xyw +# listf = c("f1") # initiate list f + } else{ listxyw = paste(listxyw, paste('xyw', ci,sep=""), sep=',') # append list xyw +# listf = c(listf, paste('f', ci,sep="")) +# paste(listf, paste('f', ci,sep=""), sep=',') # append list f + } + i = i + 1 + ci = as.character(i) # character string with value of i + } + +# nl = length(listxyw) +# listxyw = listxyw[1:(nl - 1)] # remove last ',' +# nl = length(listf) +# listf = listf[1:(nl - 1)] # remove last ',' +# +# global_txt = strrep(['global ', listxyw], ',', ' '); +# eval(global_txt); # make data sets global +# + library(pracma) + N = zeros(nxyw, 1) # initiate data counter + Y = vector(length=0) + W = vector(length=0) # initiate observations and weight coefficients + + for (i in c(1:nxyw)) { # loop across data sets + ci = as.character(i) # character string with value of i + # assing unnamed arguments to xyw1, xyw2, etc + eval(parse(text=paste('xyw', ci, ' = varargin[[', ci,']]', sep=""))) + eval(parse(text=paste('N[', ci, '] = nrow(xyw', ci, ')', sep=""))) # number of data points + eval(parse(text=paste('k= ncol(xyw', ci, ')', sep=""))) + eval(parse(text=paste('Y = c(Y,xyw', ci, '[,2])', sep=""))) # append dependent variables + if (k > 2) { + eval(parse(text=paste('W = c(W, xyw', ci, '[,3])', sep=""))) # append weight coefficients + } else {W = c(W, (ones(N[i],1)/ N[i]))} # append weight coefficients + + } + + q = p # copy input parameter matrix to output + info = 1 # convergence has been successful + + np = nrow(p) # k: number of parameters + k = ncol(p) + index = c(1:np) + if (k > 1) { + index = index[p[,2]==1] # indices of iterated parameters + } + n_par = max(length(index)) # number of parameters that must be iterated + + # set options if necessary + if (!exists('max_step_number')){ + nmregr_options(key='max_step_number', 200 * n_par) + } + if (!exists('max_fun_evals')){ + nmregr_options(key='max_fun_evals', 200 * n_par); + } + if (!exists('tol_simplex')){ + nmregr_options(key='tol_simplex', 1e-4); + } + if (!exists('tol_fun')){ + nmregr_options(key='tol_fun', 1e-4) + } + if (!exists('report')){ + nmregr_options(key='report', 1) + } + if (is.na(max_step_number)){ + nmregr_options(key='max_step_number', 200 * n_par) + } + if (is.na(max_fun_evals)){ + nmregr_options(key='max_fun_evals', 200 * n_par) + } + if (is.na(tol_simplex)){ + nmregr_options(key='tol_simplex', 1e-4) + } + if (is.na(tol_fun)){ + nmregr_options(key='tol_fun', 1e-4) + } + if (is.na(report)){ + nmregr_options(key='report', 1) + } + + # Initialize parameters + rho = 1 + chi = 2 + psi = 0.5 + sigma = 0.5 + onesn = ones(1, n_par) + two2np1 = 2:(n_par + 1) + one2n = 1:n_par + np1 = n_par + 1 + + # Set up a simplex near the initial guess. + v = zeros(n_par, np1) + fv = zeros(1,np1) + xin = q[index, 1] # Place input guess in the simplex + v[,1] = xin + eval(parse(text=paste('listf = func(q[,1],', listxyw, ')'))) + if (nxyw == 1) { + fv[,1] = t(W) %*% (listf[[1]] - Y)^2 + } else {fv[,1] = t(W) %*% ((unlist(listf)-Y)^2)} + + # Following improvement suggested by L.Pfeffer at Stanford + usual_delta = 0.05 # 5 percent deltas for non-zero terms + zero_term_delta = 0.00025 # Even smaller delta for zero elements of q + for (j in c(1:n_par)){ + y = xin + if (y[j] != 0){ + y[j] = (1 + usual_delta) * y[j] + } else {y[j] = zero_term_delta} + v[,j+1] = y + q[index,1] = y + eval(parse(text=paste('listf = func(q[,1],', listxyw, ')'))) + if (nxyw == 1) { + fv[1, j + 1] = t(W) %*% ((listf[[1]] - Y) ^ 2) + } else {fv[1, j + 1] = t(W) %*% ((unlist(listf) - Y) ^ 2)} + } + + # sort so v(1,:) has the lowest function value + j=order(fv) + fv=matrix(fv[1,j], nrow=1, byrow=T) + v = v[,j] + + how = 'initial' + itercount = 1 + func_evals = n_par + 1 + if (report == 1){ + print(paste('step ', itercount, ' ssq ', min(fv), '-', + max(fv), ' ', how)) + } + info = 1; + + # Main algorithm --------------------------------------------------------------------------------------------------------------------------------------- + # Iterate until the diameter of the simplex is less than tol_simplex + # AND the function values differ from the min by less than tol_fun, + # or the max function evaluations are exceeded. (Cannot use OR instead of AND.) + + while (func_evals < max_fun_evals && itercount < max_step_number){ + if ( max(max(abs(v[,two2np1] - v[,onesn]))) <= tol_simplex && max(abs(fv[1] - fv[two2np1])) <= tol_fun ){ + break + } + + how = '' + + # Compute the reflection point + + # xbar = average of the n (NOT n+1) best points + xbar = rowSums(v[,one2n])/ n_par + xr = (1 + rho) * xbar - rho * v[,np1] + q[index,1] = xr + eval(parse(text=paste('listf = func(q[,1],', listxyw, ')'))) + if (nxyw == 1){ + fxr = t(W) * ((listf[[1]] - Y) ^ 2) + } else { fxr = t(W) %*% ((unlist(listf) - Y) ^ 2) } + + func_evals = func_evals + 1 + + if (fxr < fv[,1]){ + # Calculate the expansion point + xe = (1 + rho * chi) * xbar - rho * chi * v[, np1] + q[index,1] = xe + eval(parse(text=paste('listf = func(q[,1],', listxyw, ')'))) + if (nxyw == 1){ + fxe = t(W) * ((listf[[1]] - Y) ^ 2) + } else {fxe = t(W) %*% ((unlist(listf) - Y) ^2)} + + func_evals = func_evals + 1 + if (fxe < fxr){ + v[,np1] = xe + fv[,np1] = fxe + how = 'expand' + } else { + v[,np1] = xr + fv[,np1] = fxr + how = 'reflect' + } + } else {# fv[,1] <= fxr + if (fxr < fv[,n_par]){ + v[,np1] = xr + fv[,np1] = fxr + how = 'reflect' + } else{ + # fxr >= fv(:,n_par) + # Perform contraction + if (fxr < fv[,np1]){ + # Perform an outside contraction + xc = (1 + psi * rho) * xbar - psi * rho * v[,np1] + q[index,1] = xc + eval(parse(text=paste('listf = func(q[,1],', listxyw, ')'))) + if (nxyw == 1){ + fxc = t(W) %*% ((listf[[1]] - Y) ^ 2) + } else{fxc = t(W) %*% ((unlist(listf) - Y) ^ 2)} + func_evals = func_evals + 1 + + if (fxc <= fxr){ + v[,np1] = xc + fv[,np1] = fxc + how = 'contract outside' + } else { how = 'shrink' } # perform a shrink + } else { + # Perform an inside contraction + xcc = (1 - psi) * xbar + psi * v[,np1] + q[index,1] = xcc + eval(parse(text=paste('listf = func(q[,1],', listxyw, ')'))) + + if (nxyw == 1){ + fxcc = t(W) %*% ((listf[[1]] - Y) ^ 2) + } else{fxcc = t(W) %*% ((unlist(listf) - Y) ^ 2)} + + func_evals = func_evals + 1 + + if (fxcc < fv[,np1]){ + v[,np1] = xcc + fv[,np1] = fxcc + how = 'contract inside' + } else { + # perform a shrink + how = 'shrink' + } + } + + if ( how == 'shrink' ){ + for (j in two2np1){ + v[,j] = v[,1] + sigma * (v[,j] - v[,1]) + q[index,1] = v[,j] + eval(parse(text=paste('listf = func(q[,1],', listxyw, ')'))) + if (nxyw == 1){ + fv[,j] = t(W) %*% ((listf[[1]] - Y) ^ 2) + } else {fv[,j] = t(W) %*% ((unlist(listf) - Y)^2)} + } + + func_evals = func_evals + n_par + } + } + } + + j=order(fv) + fv=matrix(fv[1,j], nrow=1, byrow=T) + v = v[,j] + itercount = itercount + 1 + if (report == 1){ + print(paste('step ', itercount, ' ssq ', min(fv), '-', max(fv), ' ', how, '\n')) + } + + } # while of main loop + + + q[index,1] = v[,1] + + fval = min(fv) + if (func_evals >= max_fun_evals){ + if (report > 0){ + print(paste('No convergences with ', max_fun_evals, ' function evaluations')) + } + info = 0 + } else if (itercount >= max_step_number ) { + if (report > 0){ + print(paste('No convergences with ', max_step_number, ' steps', sep='')) + } + info = 0 + } else { + if (report > 0){ + print('Successful convergence \n') + } + } + info = 1 + +return(list(q, info)) +} \ No newline at end of file diff --git a/R/nmregr_options.R b/R/nmregr_options.R new file mode 100644 index 0000000..c5c1114 --- /dev/null +++ b/R/nmregr_options.R @@ -0,0 +1,108 @@ +## nmregr_options +# Sets options for function + + ## +nmregr_options = function (key='inexistent', val=NA) { + # created at 2002/02/10 by Bas Kooijman; modified 2015/01/16, 2015/02/27 Goncalo Marques + + ## Syntax + # <../nmregr_options .m *nmregr_options*> (key, val) + + ## Description + # Sets options for function 'nmregr' one by one + # + # Input + # + # * no input: print values to screen + # * one input: + # + # 'default' sets options at default values + # other keys (see below) to print value to screen + # + # * two inputs + # + # 'report': 1 - to report steps to screen; 0 - not to; + # 'max_step_number': maximum number of steps + # 'max_fun_evals': maximum number of function evaluations + # 'tol_simplex': tolerance for how close the simplex points must be + # together to call them the same + # 'tol_tun': tolerance for how close the loss-function values must be + # together to call them the same + # + # Output + # + # * no output, but globals are set to values or values printed to screen + + ## Example of use + # nmregr_options('default'); nmregr_options('report', 0) + + + if (key=='default'){ + report <<- 1 + max_step_number <<- 500 + max_fun_evals <<- 2000 + tol_simplex <<- 1e-4 + tol_fun <<- 1e-4 + + }else if (key=='report'){ + if (!exists('val')){ + if (numel(report) != 0){ + cat('report = ', report) + } else {print('report = unknown \n')} + } else {report <<- val} + + } else if (key=='max_step_number'){ + if (!exists('val')){ + if (numel(max_step_number) != 0){ + cat('max_step_number = ', max_step_number) + } else {print('max_step_number = unknown \n')} + } else{max_step_number <<- val} + + } else if (key== 'max_fun_evals'){ + if (!exists('val')){ + if (numel(max_fun_evals) != 0){ + cat('max_fun_evals = ', max_fun_evals) + } else {print('max_fun_evals = unkown \n')} + } else {max_fun_evals <<- val} + + } else if (key== 'tol_simplex'){ + if (!exists('val')){ + if (numel(tol_simplex) != 0){ + cat('tol_simplex = ', tol_simplex) + } else {print('tol_simplex = unknown \n')} + } else {tol_simplex <<- val} + + } else if (key == 'tol_fun'){ + if (!exists('val')){ + if (numel(tol_fun) != 0){ + cat('tol_fun = ', tol_fun)} + else {print('tol_fun = unknown \n')} + }else {tol_fun <<- val} + + } else { + if (key!='inexistent'){ + cat('key ', key, ' is unkown') + } + if (numel(report) != 0){ + cat('report = ', report) + } else {print('report = unknown \n')} + + if (numel(max_step_number) != 0){ + cat('max_step_number = ', max_step_number) + } else {print('max_step_number = unkown \n')} + + if (numel(max_fun_evals) != 0){ + cat('max_fun_evals = ', max_fun_evals) + } else {print('max_fun_evals = unkown')} + + if (numel(tol_simplex) != 0){ + cat('tol_simplex = ', tol_simplex) + } else {print('tol_simplex = unknown \n')} + + if (numel(tol_fun) != 0){ + cat('tol_fun = ', tol_fun) + } else {print('tol_fun = unknown \n')} + + } + +} \ No newline at end of file diff --git a/R/repro_rate_metam.R b/R/repro_rate_metam.R new file mode 100644 index 0000000..f51bb82 --- /dev/null +++ b/R/repro_rate_metam.R @@ -0,0 +1,89 @@ +## reprod_rate_metam +# gets reproduction rate as function of time for type M acceleration + +## +# function [R, UE0, info] = +reprod_rate_metam = function(L, f, p){ + # created 2009/03/24 by Bas Kooijman, modified 2010/02/23 + + ## Syntax + # [R, UE0, info] = <../reprod_rate_metam.m *reprod_rate_metam*> (L, f, p) + + ## Description + # Calculates the reproduction rate in number of eggs per time + # in the case of acceleration between events b and j + # for an individual of length L and scaled reserve density f + # + # Input + # + # * L: n-vector with length + # * f: scalar with functional response + # * p: 12-vector with parameters: kap kapR g kJ kM LT v UHb UHp Lb Lj Lp + # + # g and v refer to the values for embryo; scaling is always with respect to embryo values + # V1-morphic juvenile between events b and j with E_Hb > E_Hj > E_Hp + # + # * Lb0: optional scalar with length at birth (initial value only) + # + # or optional 2-vector with length, L, and scaled maturity, UH < UHp + # for a juvenile that is now exposed to f, but previously at another f + # + # Output + # + # * R: n-vector with reproduction rates + # * UE0: scalar with scaled initial reserve + # * info: indicator with 1 for success, 0 otherwise + + ## Remarks + # Theory is given in comments to DEB3 section 7.8.2. + # For cumulative reproduction, see , + # *This function will be phased out* and replaced by + + ## Example of use + # p_R = [.8 .95 .42 1.7 1.7 0 0.2 .012 .0366 .01 .2 .5]; + # reprod_rate_j([.05 .1 1 3]', 1, p_R) + + # Explanation of variables: + # R = kapR * pR/ E0 + # pR = (1 - kap) pC - kJ * EHp + # [pC] = [Em] (v/ L + kM (1 + LT/L)) f g/ (f + g); pC = [pC] L^3 + # [Em] = {pAm}/ v + # + # remove energies; now in lengths, time only + # + # U0 = E0/{pAm}; UHp = EHp/{pAm}; SC = pC/{pAm}; SR = pR/{pAm} + # R = kapR SR/ U0 + # SR = (1 - kap) SC - kJ * UHp + # SC = f (g/ L + (1 + LT/L)/ Lm)/ (f + g); Lm = v/ (kM g) + # + # unpack parameters; parameter sequence, cf get_pars_r + kap = p[1] + kapR = p[2] + g = p[3] + kJ = p[4] + kM = p[5] + LT = p[6] + v = p[7] + UHb = p[8] + UHp = p[9] + Lb = p[10] + Lj = p[11] + Lp = p[12] + + Lm = v/ (kM * g) # maximum length + VHb = UHb/ (1 - kap) + + p_UE0 = c(VHb, g, kJ, kM, v) # pars for initial_scaled_reserve + UE0Lbinfo = initial_scaled_reserve(max(1,f), p_UE0, Lb) + UE0=UE0Lbinfo[1] + Lb=UE0Lbinfo[2] + info=UE0Lbinfo[3] + + SC = L ^ 2 * ((g + LT/ Lm) * Lj/ Lb + L/ Lm)/ (1 + g/ f) + SR = (1 - kap) * SC - kJ * UHp + R = (L >= Lp) * kapR * SR/ UE0 # set reprod rate of juveniles to zero + + return(c(R, UE0, info)) +} + + \ No newline at end of file diff --git a/R/spline1.R b/R/spline1.R new file mode 100644 index 0000000..4d08485 --- /dev/null +++ b/R/spline1.R @@ -0,0 +1,60 @@ +spline1 = function(x, knots, Dy1, Dyk) { + # created at 2007/03/29 by Bas Kooijman; modified 2009/09/29 + # + ## Description + # First order splines connect knots by straight lines, and is linear outside the knots. + # Calculates the ordinates and the first derivatives of a first order spline, given abcissa and knot coordinates. + # The spline is interpolating. + # + ## Input + # x: n-vector with abcissa values + # knots: (nk,2)-matrix with coordinates of knots; we must have nk > 3 + # knots(:,1) must be ascending + # Dy1: scalar with first derivative at first knot (optional) + # empty means: zero + # Dyk: scalar with first derivative at last knot (optional) + # empty means: zero + # + ## Ouput + # y: n-vector with spline values (ordinates) + # dy: n-vector with derivatives + # index: n-vector with indices of first knot-abcissa smaller than x + # + ## Remarks + # See ispline1 for integration, rspline1 for roots, espline1 for local extremes. + # + ## Example of use + # x = (1:10)'; y = 3*(x+.1*rand(10,1)).^2; [Y, dY] = spline1([x,y],k); iY = ispline1(x,k); rspline1(k,5) + # See mydata_smooth for further illustration. + + x = as.matrix(x); nx = length(x); nk = nrow(knots); + y = matrix(0,nx,1); dy = y; index = matrix(0,nx,1); # initiate output + + if(exists('Dy1') == FALSE){ # make sure that left clamp is specified + Dy1 = 0; + } + if(exists('Dyk') == FALSE){ # make sure that right clamp is specified + Dyk = 0; + } + + ## derivatives right of knot-abcissa + Dy =c((knots[2:nk,2] - knots[1:nk-1,2]) / + (knots[2:nk,1] - knots[1:nk-1,1]), Dyk); + for(i in 1:nx){ # loop across abcissa values + j = 1; + while(x[i] > knots[min(nk,j),1] & j <= nk){ + j = j + 1; + } + j = j - 1; + if(j == 0){ + y[i] = knots[1,2] - Dy1 * (knots[1,1] - x[i]); + dy[i] = Dy1; + }else{ + y[i] = knots[j,2] - Dy[j] * (knots[j,1] - x[i]); + dy[i] = Dy[j]; + } + index[i] = j; + } + +return(list(y, dy, index)) +} \ No newline at end of file diff --git a/R/tempcorr.R b/R/tempcorr.R new file mode 100644 index 0000000..9366b94 --- /dev/null +++ b/R/tempcorr.R @@ -0,0 +1,45 @@ +tempcorr = function(T,T_1,Tpars){ + # Created at 2002/04/09 by Bas Kooijman; modified 2005/01/24 + # + ## Description + # Calculates the factor with which physiological rates should be multiplied + # to go from a reference temperature to a given temperature. + # + ## Input + # T: vector with new temperatures + # T_1: scalar with reference temperature + # Tpars: 1-, 3- or 5-vector with temperature parameters: + # + ## Output + # TC: vector with temperature correction factor(s) that affect(s) all rates + # + ## Remarks + # shtempcorr shows a graph of this correction factor as function of the temperature. + # + ## Example of use + # tempcorr([330 331 332], 320, [12000 277 318 20000 190000]) and + # shtemp2corr(320, [12000 277 318 20000 190000]). + + ## Code + T_A = Tpars[1]; # Arrhenius temperature + if(length(Tpars) == 1){ + TC = exp(T_A/T_1 - T_A/T); + }else{ + if(length(Tpars) == 3){ + T_L = Tpars[2]; # Lower temp boundary + T_AL = Tpars[3]; # Arrh. temp for lower boundary + TC = exp(T_A/ T_1 - T_A/ T) * (1 + exp(T_AL/ T_1 - T_AL/ T_L)) / + (1 + exp(T_AL/ T - T_AL/ T_L)); + }else{ + T_L = Tpars[2]; # Lower temp boundary + T_H = Tpars[3]; # Upper temp boundary + T_AL = Tpars[4]; # Arrh. temp for lower boundary + T_AH = Tpars[5]; # Arrh. temp for upper boundary + + TC = exp(T_A/ T_1 - T_A/ T) * ... + (1 + exp(T_AL/ T_1 - T_AL/ T_L) + exp(T_AH/ T_H - T_AH/ T_1))/ + (1 + exp(T_AL/ T - T_AL/ T_L) + exp(T_AH/ T_H - T_AH/ T)); + }} + +return(TC) +} \ No newline at end of file