From d4f754a5d46b13760d0ecd3d0da41dee9d4a890a Mon Sep 17 00:00:00 2001 From: Philipp Boersch-Supan Date: Thu, 11 Jan 2018 17:07:22 -0500 Subject: [PATCH] more tidying of namespaces and imported functions --- DESCRIPTION | 5 +- LICENSE | 2 + NAMESPACE | 7 +++ R/get_lb2.R | 28 +++++----- R/get_lj.R | 71 +++++++++++++++--------- R/get_tm_s.R | 98 ++++++++++++++++++++------------- R/initial_scaled_reserve.R | 20 +++---- R/nmregr.R | 110 ++++++++++++++++++++----------------- R/tempcorr.R | 20 +++---- man/get_lj.Rd | 26 +++++++++ man/get_tm_s.Rd | 31 +++++++++++ man/nmregr.Rd | 21 +++++++ 12 files changed, 290 insertions(+), 149 deletions(-) create mode 100644 LICENSE create mode 100644 man/get_lj.Rd create mode 100644 man/get_tm_s.Rd create mode 100644 man/nmregr.Rd diff --git a/DESCRIPTION b/DESCRIPTION index b56da7e..c58a83c 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -5,9 +5,10 @@ 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? +License: MIT + file LICENSE Encoding: UTF-8 LazyData: true RoxygenNote: 6.0.1 Suggests: testthat -Imports: pracma +Imports: pracma, + deSolve diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..d456648 --- /dev/null +++ b/LICENSE @@ -0,0 +1,2 @@ +YEAR: 2018 +COPYRIGHT HOLDER: Your name goes here diff --git a/NAMESPACE b/NAMESPACE index b5caa03..7ee99c6 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -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) diff --git a/R/get_lb2.R b/R/get_lb2.R index 7b8b0a2..34dcd41 100644 --- a/R/get_lb2.R +++ b/R/get_lb2.R @@ -9,7 +9,7 @@ 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 @@ -17,10 +17,10 @@ get_lb2= function(p, eb, lb0=NA){ # * 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 @@ -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 } @@ -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)) } - + diff --git a/R/get_lj.R b/R/get_lj.R index 9b36b7d..c45f234 100644 --- a/R/get_lj.R +++ b/R/get_lj.R @@ -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 # @@ -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 , 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 for foetal development. - + # 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 @@ -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 @@ -79,13 +98,13 @@ 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 @@ -93,23 +112,23 @@ get_lj=function(p, 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) @@ -122,4 +141,4 @@ get_lj=function(p, f, lb0) { return(c(lj, lp, lb, info)) -} \ No newline at end of file +} diff --git a/R/get_tm_s.R b/R/get_tm_s.R index e7543a4..abea014 100644 --- a/R/get_tm_s.R +++ b/R/get_tm_s.R @@ -2,47 +2,71 @@ # Obtains scaled mean age at death fro short growth periods ## -# function [tm, Sb, Sp, info] = +# function [tm, Sb, Sp, info] = + +#' Obtains scaled mean age at death fro short growth periods +#' +#' 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. +#' +#' @param p 4 or 7-vector with parameters: [g lT ha sG] or [g k lT vHb vHp ha SG] +#' @param f optional scalar with scaled reserve density at birth (default eb = 1) +#' @param lb optional scalar with scaled length at birth (default: lb is obtained from get_lb) +#' @param lp optional scalar with scaled length at puberty +#' +#' @return [tm, Sb, Sp, info] +#' * 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 +#' +#' @export +#' +#' @importFrom pracma quad expint erfc +#' 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 - # + # * p: + # * f: + # * lb: + # * lp: + # # 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 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 @@ -59,11 +83,11 @@ 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 @@ -71,7 +95,7 @@ get_tm_s = function(p, f, lb, lp) { 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')) { @@ -79,26 +103,26 @@ get_tm_s = function(p, f, lb, lp) { lp=lplbinfo[1] lb=lplbinfo[2] info_lp=lplbinfo[3] - } else if (!exist('lp') || is.na('lp')){ + } else if (!exists('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) + Sp = exp((1 - exp(hGtp) + hGtp + hGtp^2/ 2) * 6/ tG3) if (info_lp == 1 && info_tb == 1){ info = 1 } else{ @@ -109,33 +133,33 @@ 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) - 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 + #library(pracma) + tm_tail = pracma::expint(exp(tm * hG) * 6/ tG3)/ hG + tm = pracma::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 + tm_tail = sqrt(pi)* pracma::erfc(tm * hw)/ 2/ hw + tm = pracma::quad(fnget_tm_s, 0, tm * hW, tol = 1.0e-12, trace = FALSE, tG)/ hW + tm_tail } } - + # subfunction - - #function S = + + #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 + # t: age * hW # S: ageing survival prob - + hGt = tG * t # age * hG if (tG > 0) { a= matrix(hGt) %*% t(matrix((1/(4:500)))) @@ -143,5 +167,5 @@ fnget_tm_s=function(t, tG){ } else {# tG < 0 S = exp((1 + hGt + hGt^2/ 2 - exp(hGt)) * 6/ tG^3) } - return(S) -} \ No newline at end of file + return(S) +} diff --git a/R/initial_scaled_reserve.R b/R/initial_scaled_reserve.R index a4d7e48..4c6db35 100644 --- a/R/initial_scaled_reserve.R +++ b/R/initial_scaled_reserve.R @@ -7,10 +7,10 @@ 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 # @@ -25,23 +25,23 @@ initial_scaled_reserve=function(f, p, Lb0){ # * 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 + + ## 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) @@ -50,7 +50,7 @@ initial_scaled_reserve=function(f, p, Lb0){ 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 + lb0 = pracma::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]) @@ -62,6 +62,6 @@ initial_scaled_reserve=function(f, p, Lb0){ 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 index 80f9a93..7775cf0 100644 --- a/R/nmregr.R +++ b/R/nmregr.R @@ -1,16 +1,26 @@ # nmregr -# Calculates least squares estimates using Nelder Mead's simplex method +# ## -#function [q, info] = +#function [q, info] = +#' Calculates least squares estimates using Nelder Mead's simplex method +#' +#' @param func string with name of user-defined function +#' @param p (k,2)-matrix with p(:,1) initial guesses for parameter values; p(:,2) binaries with yes or no iteration (optional) +#' @param varargin MATLAB idiom... +#' +#' @return [q, info] +#' +#' @importFrom pracma zeros ones +#' 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 # @@ -42,31 +52,31 @@ nmregr=function(func, p, varargin){ # # * 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 Newton-Raphson method, # for genetic algorithm, - # for 2 independent variables, and + # 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 @@ -78,17 +88,17 @@ nmregr=function(func, p, varargin){ 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 +# + #library(pracma) + N = pracma::zeros(nxyw, 1) # initiate data counter Y = vector(length=0) W = vector(length=0) # initiate observations and weight coefficients @@ -104,10 +114,10 @@ nmregr=function(func, p, varargin){ } 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) @@ -115,7 +125,7 @@ nmregr=function(func, p, varargin){ 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) @@ -147,20 +157,20 @@ nmregr=function(func, p, varargin){ 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) + onesn = pracma::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) + v = pracma::zeros(n_par, np1) + fv = pracma::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, ')'))) @@ -182,13 +192,13 @@ nmregr=function(func, p, varargin){ 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 + } + + # 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 @@ -197,7 +207,7 @@ nmregr=function(func, p, varargin){ 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, @@ -207,11 +217,11 @@ nmregr=function(func, p, varargin){ 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] @@ -222,7 +232,7 @@ nmregr=function(func, p, varargin){ } 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] @@ -231,7 +241,7 @@ nmregr=function(func, p, varargin){ 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 @@ -242,13 +252,13 @@ nmregr=function(func, p, varargin){ fv[,np1] = fxr how = 'reflect' } - } else {# fv[,1] <= fxr + } else {# fv[,1] <= fxr if (fxr < fv[,n_par]){ v[,np1] = xr fv[,np1] = fxr how = 'reflect' } else{ - # fxr >= fv(:,n_par) + # fxr >= fv(:,n_par) # Perform contraction if (fxr < fv[,np1]){ # Perform an outside contraction @@ -259,7 +269,7 @@ nmregr=function(func, p, varargin){ 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 @@ -270,13 +280,13 @@ nmregr=function(func, p, varargin){ 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 @@ -286,7 +296,7 @@ nmregr=function(func, p, varargin){ how = 'shrink' } } - + if ( how == 'shrink' ){ for (j in two2np1){ v[,j] = v[,1] + sigma * (v[,j] - v[,1]) @@ -296,12 +306,12 @@ nmregr=function(func, p, varargin){ 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] @@ -311,10 +321,10 @@ nmregr=function(func, p, varargin){ } } # while of main loop - - + + q[index,1] = v[,1] - + fval = min(fv) if (func_evals >= max_fun_evals){ if (report > 0){ @@ -328,10 +338,10 @@ nmregr=function(func, p, varargin){ info = 0 } else { if (report > 0){ - print('Successful convergence \n') + print('Successful convergence \n') } } info = 1 return(list(q, info)) -} \ No newline at end of file +} diff --git a/R/tempcorr.R b/R/tempcorr.R index 9366b94..030ea0f 100644 --- a/R/tempcorr.R +++ b/R/tempcorr.R @@ -2,8 +2,8 @@ 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. + # 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 @@ -12,13 +12,13 @@ tempcorr = function(T,T_1,Tpars){ # ## 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. + # 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]). + # 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 @@ -28,7 +28,7 @@ tempcorr = function(T,T_1,Tpars){ 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)) / + 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 @@ -36,10 +36,10 @@ tempcorr = function(T,T_1,Tpars){ 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))/ + 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 +} diff --git a/man/get_lj.Rd b/man/get_lj.Rd new file mode 100644 index 0000000..38ef607 --- /dev/null +++ b/man/get_lj.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_lj.R +\name{get_lj} +\alias{get_lj} +\title{Gets scaled length at metamorphosis} +\usage{ +get_lj(p, f, lb0) +} +\arguments{ +\item{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} + +\item{f}{optional scalar with scaled functional responses (default 1)} + +\item{lb0}{optional scalar with scaled length at birth} +} +\value{ +[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 +} +\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). +} diff --git a/man/get_tm_s.Rd b/man/get_tm_s.Rd new file mode 100644 index 0000000..b650dca --- /dev/null +++ b/man/get_tm_s.Rd @@ -0,0 +1,31 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_tm_s.R +\name{get_tm_s} +\alias{get_tm_s} +\title{Obtains scaled mean age at death fro short growth periods} +\usage{ +get_tm_s(p, f, lb, lp) +} +\arguments{ +\item{p}{4 or 7-vector with parameters: [g lT ha sG] or [g k lT vHb vHp ha SG]} + +\item{f}{optional scalar with scaled reserve density at birth (default eb = 1)} + +\item{lb}{optional scalar with scaled length at birth (default: lb is obtained from get_lb)} + +\item{lp}{optional scalar with scaled length at puberty} +} +\value{ +[tm, Sb, Sp, info] +* 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 +} +\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. +} diff --git a/man/nmregr.Rd b/man/nmregr.Rd new file mode 100644 index 0000000..69c615d --- /dev/null +++ b/man/nmregr.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/nmregr.R +\name{nmregr} +\alias{nmregr} +\title{Calculates least squares estimates using Nelder Mead's simplex method} +\usage{ +nmregr(func, p, varargin) +} +\arguments{ +\item{func}{string with name of user-defined function} + +\item{p}{(k,2)-matrix with p(:,1) initial guesses for parameter values; p(:,2) binaries with yes or no iteration (optional)} + +\item{varargin}{MATLAB idiom...} +} +\value{ +[q, info] +} +\description{ +Calculates least squares estimates using Nelder Mead's simplex method +}