Skip to content

Commit

Permalink
more tidying of namespaces and imported functions
Browse files Browse the repository at this point in the history
  • Loading branch information
pboesu committed Jan 11, 2018
1 parent 5e74f77 commit d4f754a
Show file tree
Hide file tree
Showing 12 changed files with 290 additions and 149 deletions.
5 changes: 3 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,10 @@ 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?
License: MIT + file LICENSE
Encoding: UTF-8
LazyData: true
RoxygenNote: 6.0.1
Suggests: testthat
Imports: pracma
Imports: pracma,
deSolve
2 changes: 2 additions & 0 deletions LICENSE
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
YEAR: 2018
COPYRIGHT HOLDER: Your name goes here
7 changes: 7 additions & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,13 @@
export(beta0)
export(dget_tb)
export(get_lb)
export(get_lj)
export(get_tb)
export(get_tm_s)
export(get_ue0)
importFrom(deSolve,ode)
importFrom(pracma,erfc)
importFrom(pracma,expint)
importFrom(pracma,ones)
importFrom(pracma,quad)
importFrom(pracma,zeros)
28 changes: 14 additions & 14 deletions R/get_lb2.R
Original file line number Diff line number Diff line change
Expand Up @@ -9,18 +9,18 @@ get_lb2= function(p, eb, lb0=NA){
# [lb info] = <../get_lb2.m *get_lb2*> (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.
# 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
# * lb: scalar with scaled length at birth
# * info: indicator equals 1 if successful, 0 otherwise

## Remarks
Expand All @@ -30,12 +30,12 @@ get_lb2= function(p, eb, lb0=NA){
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
lb = as.complex(vHb)^(1/ 3) # exact solution for k = 1
} else {
lb = lb0
}
Expand All @@ -44,22 +44,22 @@ get_lb2= function(p, eb, lb0=NA){
} 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 = pracma::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))
}

71 changes: 45 additions & 26 deletions R/get_lj.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,23 +2,42 @@
# Gets scaled length at metamorphosis

##
#function [lj, lp, lb, info] =
#function [lj, lp, lb, info] =

#' Gets scaled length at metamorphosis
#'
#' 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).
#'
#' @param 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
#' @param f optional scalar with scaled functional responses (default 1)
#' @param lb0 optional scalar with scaled length at birth
#'
#' @return [lj, lp, lb, info] 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
#'
#' @export
#'
#' @importFrom deSolve ode
#'
get_lj=function(p, f, lb0) {

# created at 2010/02/10 by Bas Kooijman,
# 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).
# 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
#
Expand All @@ -35,22 +54,22 @@ get_lj=function(p, f, lb0) {
# * lp: scalar with scaled length at puberty
# * lb: scalar with scaled length at birth
# * info: indicator equals 1 if successful

## Remarks
# Similar to <get_lj1.html get_lj1>, 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,
# 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 <get_lj_foetus.html *get_lj_foetus*> for foetal development.
# See <get_lj_foetus.html *get_lj_foetus*> 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
Expand All @@ -60,13 +79,13 @@ get_lj=function(p, f, lb0) {
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
Expand All @@ -79,37 +98,37 @@ get_lj=function(p, f, lb0) {
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
lp = NA
} # else { # puberty specified
# [lp, lb, info] = get_tp(p([1 2 3 4 6]), f, lb0);
# lj = lb;
# lj = lb;
# }
}


# get lj
library(deSolve)
#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")
out = deSolve::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)
Expand All @@ -122,4 +141,4 @@ get_lj=function(p, f, lb0) {

return(c(lj, lp, lb, info))

}
}
Loading

0 comments on commit d4f754a

Please sign in to comment.