tidied up pracma dependency, docs

DESCRIPTION

@@ -10,3 +10,4 @@ Encoding: UTF-8
 LazyData: true
 RoxygenNote: 6.0.1
 Suggests: testthat
+Imports: pracma

NAMESPACE

@@ -1,5 +1,7 @@
 # Generated by roxygen2: do not edit by hand
 
 export(dget_tb)
+export(get_lb)
 export(get_tb)
 export(get_ue0)
+importFrom(pracma,quad)

R/get_lb.R

@@ -1,66 +1,77 @@
-## get_lb
-# Obtains scaled length at birth, given the scaled reserve density at birth
-
-##
-get_lb = function(p, eb, lb0=NA){
+#' get_lb
+#'
+#' Obtains scaled length at birth, given the scaled reserve density at birth
+#'
+#'
+#'
+#' @param p 3-vector with parameters: g, k, v_H^b (see below)
+#' @param eb optional scalar with scaled reserve density at birth (default eb = 1)
+#' @param lb0 optional scalar with initial estimate for scaled length at birth (default lb0 = NA, will use lb for k = 1)
+#'
+#' @author Bas Kooijman
+#'
+#' @return [lb, info] = lb: scalar with scaled length at birth; info: indicator equals 1 if successful, 0 otherwise
+#' @export
+#'
+get_lb = function(p, eb = 1, 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. 
+  # 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 
+  # * 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 
+  # 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. 
+  # 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 (!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     
+    lb = as.complex(vHb)^(1/3) # exact solution for k = 1
   } else {lb = lb0}
 
   if (!exists('eb')){
@@ -75,13 +86,13 @@ get_lb = function(p, eb, lb0=NA){
   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
@@ -104,7 +115,7 @@ get_lb = function(p, eb, lb0=NA){
     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)){
@@ -121,11 +132,11 @@ get_lb = function(p, eb, lb0=NA){
       info=lbinfo[2]
     }
   }
-  
+
   if (info == 0){
-    print('warning get_lb: no convergence of l_b')
+    warning('get_lb: no convergence of l_b')
   }
-  
+
   return(c(lb, info))
 
-}
+}

R/get_tb.R

@@ -11,6 +11,7 @@
 #' @author Bas Kooijman
 #'
 #' @return 3-element vector containing 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
+#' @importFrom pracma quad
 #' @export
 #'
 #' @examples
@@ -75,7 +76,7 @@ get_tb <- function(p, eb = 1, lb=NA){
 
   xb = g/ (eb + g)  # f = e_b
   ab = 3 * g * xb^(1/ 3)/ lb  # \alpha_b
-  library(pracma)
+  #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))

R/get_tm.R

@@ -2,50 +2,50 @@
 # Obtains scaled mean age at death fro short growth periods
 
 ##
-  #function [tm, Sb, Sp, info] = 
+  #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. 
+  # 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.   
+  # 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. 
+  # Theory is given in comments on DEB3 Section 6.1.1.
   # See <get_tm.html *get_tm*> 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
@@ -63,19 +63,19 @@ get_tm_s = function(p, f, lb, lp){
     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)) {
@@ -84,7 +84,7 @@ get_tm_s = function(p, f, lb, 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]
@@ -95,7 +95,7 @@ get_tm_s = function(p, f, lb, 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) 
+    Sp = exp((1 - exp(hGtp) + hGtp + hGtp^2/ 2) * 6/ tG3)
     if (info_lp == 1 & info_tb == 1){
       info = 1
     }
@@ -106,39 +106,39 @@ get_tm_s = function(p, f, lb, lp){
   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)
+    #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
+    tm = pracma::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)
+    #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
+    tm = pracma::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 
+  # 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)
@@ -148,4 +148,4 @@ fnget_tm_s = function(t, tG){
   }
 
   return(S)
-}
+}

R/get_ue0.R

@@ -74,7 +74,7 @@ get_ue0=function(p, eb = 1, lb0 = NA){
   xb = g/ (eb + g)
   uE0 = (3 * g/ (3 * g * xb^(1/ 3)/ lb - beta0(0, xb)))^3
 
-  print(paste("uE0", uE0))
+  #print(paste("uE0", uE0))
 
   return(c(uE0, lb, info))
 }