initial commit with code pasted from DEBtool_R

.Rbuildignore

@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

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

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 <[email protected]>
+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

NAMESPACE

@@ -0,0 +1 @@
+exportPattern("^[[:alpha:]]+")

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)
+
+
+}

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)
+}

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))
+}

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))
+}

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))
+}

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)
+}

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 
+  #    <http://www.bio.vu.nl/thb/research/bib/Kooy2009b.html Kooy2009b>.
+  # 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: <get_lb1.html *get_lb1*>,
+  #  but that function is much lower.
+  # Shooting method: <get_lb2.html *get_lb2*>.
+  # Bisection method (via fzero): <get_lb3.html *get_lb3*>.
+  # In case of no convergence, <get_lb2.html *get_lb2*> is run automatically as backup.
+  # Consider the application of <get_lb_foetus.html *get_lb_foetus*> 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))
+
+}

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 <get_lb.html *get_lb*>, 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))
+}
+