API reference

Importing Model

nemf.load_model(model_path, ref_data_path=None)[source]

Loads a model by reading its model configuration from from file

Parameters:
Returns:

model – class objects that contains the model
Return type:

nemf_model class

Model Functions

class nemf.models.model_class(model_path, ref_data_path=None)[source]
load_reference_data(ref_data_path=None, **kwargs)[source]

Loads reference data used in model optimization from file

Either, the path to the reference data is provided in the yaml configuration file or passed to this function. Latter overwrites the path in the configuration file.

Parameters:ref_data_path (string (optional)) – path to the file containing the reference data
nemf.models.import_interaction_functions(func)[source]

Adds the functions from ‘func’ to the globals in models

func : list
List containing the functions that will be added to globals

Visualizing Model

nemf.interaction_graph(model)[source]

shows a graph/network of all compartments and their interactions

Forcasting Model

nemf.forward_model(model, method='Radau', verbose=False, t_eval=None, **kwargs)[source]

Runs the time integration for a provided model configuration.

Parameters:
  • model (model_class object) – class object containing the model configuration and its related methods. See load_configuration
  • method (string, optional) –

    Type of solver used for the initial-value problem aka forecasting. Should be on of:

    • ’RK45’ (default): Explicit Runge-Kutta method of order 5(4) [1]. The error is controlled assuming accuracy of the fourth-order method, but steps are taken using the fifth-order accurate formula (local extrapolation is done). A quartic interpolation polynomial is used for the dense output [2]. Can be applied in the complex domain.
    • ’RK23’: Explicit Runge-Kutta method of order 3(2) [3]. The error is controlled assuming accuracy of the second-order method, but steps are taken using the third-order accurate formula (local extrapolation is done). A cubic Hermite polynomial is used for the dense output. Can be applied in the complex domain.
    • ’DOP853’: Explicit Runge-Kutta method of order 8 [13]. Python implementation of the “DOP853” algorithm originally written in Fortran [14]. A 7-th order interpolation polynomial accurate to 7-th order is used for the dense output. Can be applied in the complex domain.
    • ’Radau’: Implicit Runge-Kutta method of the Radau IIA family of order 5 [4]. The error is controlled with a third-order accurate embedded formula. A cubic polynomial which satisfies the collocation conditions is used for the dense output.
    • ’BDF’: Implicit multi-step variable-order (1 to 5) method based on a backward differentiation formula for the derivative approximation [5]. The implementation follows the one described in [6]. A quasi-constant step scheme is used and accuracy is enhanced using the NDF modification. Can be applied in the complex domain.
    • ’LSODA’: Adams/BDF method with automatic stiffness detection and switching [7], [8]. This is a wrapper of the Fortran solver from ODEPACK.

    Explicit Runge-Kutta methods (‘RK23’, ‘RK45’, ‘DOP853’) should be used for non-stiff problems and implicit methods (‘Radau’, ‘BDF’) for stiff problems [9]. Among Runge-Kutta methods, ‘DOP853’ is recommended for solving with high precision (low values of rtol and atol). If not sure, first try to run ‘RK45’. If it makes unusually many iterations, diverges, or fails, your problem is likely to be stiff and you should use ‘Radau’ or ‘BDF’. ‘LSODA’ can also be a good universal choice, but it might be somewhat less convenient to work with as it wraps old Fortran code.

  • verbose (bool, optional) – Flag for extra verbosity during runtime
  • t_eval (1d-array, optional) – contains time stamps in posix time for which a solution shall be found and returned.
Returns:

model – class object containing the model configuration, model run results, and its related methods
Return type:

model_class object

References

[1](1, 2, 3) J. R. Dormand, P. J. Prince, “A family of embedded Runge-Kutta formulae”, Journal of Computational and Applied Mathematics, Vol. 6, No. 1, pp. 19-26, 1980.
[2](1, 2) L. W. Shampine, “Some Practical Runge-Kutta Formulas”, Mathematics of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
[3](1, 2, 3) P. Bogacki, L.F. Shampine, “A 3(2) Pair of Runge-Kutta Formulas”, Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
[4](1, 2) E. Hairer, G. Wanner, “Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems”, Sec. IV.8.
[5](1, 2) Backward Differentiation Formula on Wikipedia.
[6](1, 2) L. F. Shampine, M. W. Reichelt, “THE MATLAB ODE SUITE”, SIAM J. SCI. COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
[7](1, 2) A. C. Hindmarsh, “ODEPACK, A Systematized Collection of ODE Solvers,” IMACS Transactions on Scientific Computation, Vol 1., pp. 55-64, 1983.
[8](1, 2) L. Petzold, “Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations”, SIAM Journal on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148, 1983.
[9](1, 2, 3) Stiff equation on Wikipedia.
[13](1, 2, 3) E. Hairer, S. P. Norsett G. Wanner, “Solving Ordinary Differential Equations I: Nonstiff Problems”, Sec. II.
[14](1, 2) Page with original Fortran code of DOP853.

Fitting Model

nemf.inverse_model(model, nlp_method='SLSQP', ivp_method='Radau', sample_sets=3, maxiter=1000, seed=137, verbose=False, debug=False)[source]

Fits the model to data.

Optimizes a set of randomly generated free parameters and returns their optimized values and the corresponding fit-model and cost- function output

Parameters:
  • model (model_class object) – class object containing the model configuration and its related methods. See load_configuration()
  • nlp_method (string, optional) –

    Type of solver for the non-linear-programming problem aka fitting. Should be one of:

    • ‘trust-constr’
    • ‘SLSQP’
    • ’L-BFGS-B’
    • ’TNC
    • ’Powell’

    For problems with constraints use one of:

    • ‘trust-constr’
    • ‘SLSQP’
  • ivp_method (string, optional) –

    Type of solver used for the initial-value problem aka forecasting. Should be on of:

    • ’Radau’ (default):
      Implicit Runge-Kutta method of the Radau IIA family
    • ’RK45’: Explicit Runge-Kutta method of order 5(4) [1].
    • ’RK23’: Explicit Runge-Kutta method of order 3(2) [3].
    • ’DOP853’: Explicit Runge-Kutta method of order 8 [13].
    • ’BDF’: Implicit multi-step variable-order (1 to 5) method
    • ’LSODA’: Adams/BDF method with automatic stiffness detection

    Explicit Runge-Kutta methods (‘RK23’, ‘RK45’, ‘DOP853’) should be used for non-stiff problems and implicit methods (‘Radau’, ‘BDF’) for stiff problems [9]. Among Runge-Kutta methods, ‘DOP853’ is recommended for solving with high precision (low values of rtol and atol). If not sure, first try to run ‘RK45’. If it makes unusually many iterations, diverges, or fails, your problem is likely to be stiff and you should use ‘Radau’ or ‘BDF’. ‘LSODA’ can also be a good universal choice, but it might be somewhat less convenient to work with as it wraps old Fortran code.

  • sample_sets (positive integer, optional) – Amount of randomly generated sample sets used as initial free parameters
  • maxiter (positive integer, optional) – Maximal amount of iterations allowed in the gradient descent algorithm.
  • seed (positive integer, optional) – Initializes the random number generator. Used to recreate the same set of pseudo-random numbers. Helpfull when debugging.
  • verbose (boo, optional) – Flag for extra verbosity during runtime
Returns:

model – class object containing the model configuration, model run results (parameters, model, prediction, cost), and its related methods
Return type:

model_class object

Notes

Non-linear-programming solvers aka minimizers

Bound-Constrained minimization

Method L-BFGS-B uses the L-BFGS-B algorithm [B6], [B7] for bound constrained minimization.

Method Powell is a modification of Powell’s method [B3], [B4] which is a conjugate direction method. It performs sequential one-dimensional minimizations along each vector of the directions set (direc field in options and info), which is updated at each iteration of the main minimization loop. The function need not be differentiable, and no derivatives are taken. If bounds are not provided, then an unbounded line search will be used. If bounds are provided and the initial guess is within the bounds, then every function evaluation throughout the minimization procedure will be within the bounds. If bounds are provided, the initial guess is outside the bounds, and direc is full rank (default has full rank), then some function evaluations during the first iteration may be outside the bounds, but every function evaluation after the first iteration will be within the bounds. If direc is not full rank, then some parameters may not be optimized and the solution is not guaranteed to be within the bounds.

Method TNC uses a truncated Newton algorithm [B5], [B8] to minimize a function with variables subject to bounds. This algorithm uses gradient information; it is also called Newton Conjugate-Gradient. It differs from the Newton-CG method described above as it wraps a C implementation and allows each variable to be given upper and lower bounds.

Constrained Minimization

Method SLSQP uses Sequential Least SQuares Programming to minimize a function of several variables with any combination of bounds, equality and inequality constraints. The method wraps the SLSQP Optimization subroutine originally implemented by Dieter Kraft [B12]. Note that the wrapper handles infinite values in bounds by converting them into large floating values. Method trust-constr is a trust-region algorithm for constrained optimization. It swiches between two implementations depending on the problem definition. It is the most versatile constrained minimization algorithm implemented in SciPy and the most appropriate for large-scale problems. For equality constrained problems it is an implementation of Byrd-Omojokun Trust-Region SQP method described in [B17] and in [B5], p. 549. When inequality constraints are imposed as well, it swiches to the trust-region interior point method described in [B16]. This interior point algorithm, in turn, solves inequality constraints by introducing slack variables and solving a sequence of equality-constrained barrier problems for progressively smaller values of the barrier parameter. The previously described equality constrained SQP method is used to solve the subproblems with increasing levels of accuracy as the iterate gets closer to a solution.

Initial-Value-Problem solvers aka forecasting

The available options are:

  • ‘RK45’ (default): Explicit Runge-Kutta method of order 5(4) [1].
    The error is controlled assuming accuracy of the fourth-order method, but steps are taken using the fifth-order accurate formula (local extrapolation is done). A quartic interpolation polynomial is used for the dense output [2]. Can be applied in the complex domain.
  • ‘RK23’: Explicit Runge-Kutta method of order 3(2) [3]. The error
    is controlled assuming accuracy of the second-order method, but steps are taken using the third-order accurate formula (local extrapolation is done). A cubic Hermite polynomial is used for the dense output. Can be applied in the complex domain.
  • ‘DOP853’: Explicit Runge-Kutta method of order 8 [13].
    Python implementation of the “DOP853” algorithm originally written in Fortran [14]. A 7-th order interpolation polynomial accurate to 7-th order is used for the dense output. Can be applied in the complex domain.
  • ‘Radau’: Implicit Runge-Kutta method of the Radau IIA family of
    order 5 [4]. The error is controlled with a third-order accurate embedded formula. A cubic polynomial which satisfies the collocation conditions is used for the dense output.
  • ‘BDF’: Implicit multi-step variable-order (1 to 5) method based
    on a backward differentiation formula for the derivative approximation [5]. The implementation follows the one described in [6]. A quasi-constant step scheme is used and accuracy is enhanced using the NDF modification. Can be applied in the complex domain.
  • ‘LSODA’: Adams/BDF method with automatic stiffness detection and
    switching [7], [8]. This is a wrapper of the Fortran solver from ODEPACK.

Explicit Runge-Kutta methods (‘RK23’, ‘RK45’, ‘DOP853’) should be used for non-stiff problems and implicit methods (‘Radau’, ‘BDF’) for stiff problems [9]. Among Runge-Kutta methods, ‘DOP853’ is recommended for solving with high precision (low values of rtol and atol). If not sure, first try to run ‘RK45’. If it makes unusually many iterations, diverges, or fails, your problem is likely to be stiff and you should use ‘Radau’ or ‘BDF’. ‘LSODA’ can also be a good universal choice, but it might be somewhat less convenient to work with as it wraps old Fortran code.

References

[B3]Powell, M J D. 1964. An efficient method for finding the minimum of a function of several variables without calculating derivatives. The Computer Journal 7: 155-162.
[B4]Press W, S A Teukolsky, W T Vetterling and B P Flannery. Numerical Recipes (any edition), Cambridge University Press.
[B5](1, 2) Nocedal, J, and S J Wright. 2006. Numerical Optimization. Springer New York.
[B6]Byrd, R H and P Lu and J. Nocedal. 1995. A Limited Memory Algorithm for Bound Constrained Optimization. SIAM Journal on Scientific and Statistical Computing 16 (5): 1190-1208.
[B7]Zhu, C and R H Byrd and J Nocedal. 1997. L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization. ACM Transactions on Mathematical Software 23 (4): 550-560.
[B8]Nash, S G. Newton-Type Minimization Via the Lanczos Method. 1984. SIAM Journal of Numerical Analysis 21: 770-778.
[B12]Kraft, D. A software package for sequential quadratic programming. 1988. Tech. Rep. DFVLR-FB 88-28, DLR German Aerospace Center – Institute for Flight Mechanics, Koln, Germany.
[B16]Byrd, Richard H., Mary E. Hribar, and Jorge Nocedal. 1999. An interior point algorithm for large-scale nonlinear programming. SIAM Journal on Optimization 9.4: 877-900.
[B17]Lalee, Marucha, Jorge Nocedal, and Todd Plantega. 1998. On the implementation of an algorithm for large-scale equality constrained optimization. SIAM Journal on Optimization 8.3: 682-706.

Visualizing Results

nemf.output_summary(model)[source]

reads the data saved in the model class and depending on this data chooses a visualization method with the help of draw_output_summary to present the results

Interaction Functions

Here a list of all currently implemented interaction functions is presented. These can be used in the YAML model configuration to describe interactions between to compartments.

Note

The compartments are referenced via indices in the implementation of the interaction functions. This is simply an implementation detail. The user can (and should) use the compartment names when referencing them in the YAML configuration file. The framework handles the mapping of the names to the corresponding indeces internally.

nemf.interaction_functions.excretion(X, idx_A, idx_B, coefficient)

linear response with respect to origin/prey compartment

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • coefficient (float) – governs the slope of the linear response
Returns:

df – change in the origin and destitnation compartment. Calculated by coefficient*origin_compartment
Return type:

float
nemf.interaction_functions.exudation(X, idx_A, idx_B, coefficient)

linear response with respect to origin/prey compartment

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • coefficient (float) – governs the slope of the linear response
Returns:

df – change in the origin and destitnation compartment. Calculated by coefficient*origin_compartment
Return type:

float
nemf.interaction_functions.holling_type_0(X, idx_A, coefficient)[source]

linear response with respect to destination/predator compartment

For examples see: Examples

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • coefficient (float) – governs the slope of the linear response
Returns:

df – change in the origin and destitnation compartment. Calculated by coefficient*destination_compartment
Return type:

float
nemf.interaction_functions.holling_type_I(X, idx_A, idx_B, coefficient)[source]
linear response with respect to both origin/pray and
destination/predator compartment.

For examples see: Examples

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • coefficient (float) – governs the slope of the linear response
Returns:

df – change in the origin and destitnation compartment. Calculated by coefficient * origin_compartment * destination_compartment
Return type:

float
nemf.interaction_functions.holling_type_II(X, idx_A, idx_B, food_processing_time, hunting_rate)[source]
non-linear response with respect to origin/pray with linear response
with respect to destination/predator compartment

The response with respect to the origin compartment ‘B’ is approximately linear for small ‘B’ and converges towards an upper limit governed by the ‘food_processing_time’ for large ‘B’. For examples see: Examples

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • coefficient (float) – governs the slope of the linear response
Returns:

df – change in the origin and destitnation compartment. Calculated by consumption_rate = ((hunting_rate * origin_compartment) / (1 + hunting_rate * food_processing_time * origin_compartment)) * destination_compartment
Return type:

float
nemf.interaction_functions.holling_type_III(X, idx_A, idx_B, saturation_rate, consumption_rate_limit)[source]
non-linear response with respect to origin/pray with linear response
with respect to destination/predator compartment

The response with respect to the origin compartment ‘B’ is approximately quadratic for small ‘B’ and converges towards an upper limit governed by the ‘food_processing_time’ for large ‘B’. For examples see: Examples

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • saturation_rate (float) – first parameter of the interaction. governs the slope of the non lineare response for small pray populations.
  • consumption_rate_limit (float) – second parameter of the interaction. governs the upper limit of the response.
Returns:

df – change in the origin and destitnation compartment. Calculated by consumption_rate = ((consumption_rate_limit * saturation_rate * B**2)/

(consumption_rate_limit + (saturation_rate*B**2)))*A
Return type:

float
nemf.interaction_functions.inverse_type_0(X, idx_A, idx_B, coefficient)[source]

linear response with respect to origin/prey compartment

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • coefficient (float) – governs the slope of the linear response
Returns:

df – change in the origin and destitnation compartment. Calculated by coefficient*origin_compartment
Return type:

float
nemf.interaction_functions.linear_mortality(X, idx_A, idx_B, coefficient)

linear response with respect to origin/prey compartment

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • coefficient (float) – governs the slope of the linear response
Returns:

df – change in the origin and destitnation compartment. Calculated by coefficient*origin_compartment
Return type:

float
nemf.interaction_functions.non_negative_dual_nutrient_limited_growth(X, idx_A, idx_B, growth_rate, half_saturation, second_nutrient, second_half_saturation)[source]
nemf.interaction_functions.non_negative_excretion(X, idx_A, idx_B, coefficient)
nemf.interaction_functions.non_negative_exudation(X, idx_A, idx_B, coefficient)
nemf.interaction_functions.non_negative_holling_type_0(X, idx_A, coefficient)[source]
nemf.interaction_functions.non_negative_holling_type_I(X, idx_A, idx_B, coefficient)[source]
nemf.interaction_functions.non_negative_holling_type_II(X, idx_A, idx_B, food_processing_time, hunting_rate)[source]
nemf.interaction_functions.non_negative_inverse_type_0(X, idx_A, idx_B, coefficient)[source]
nemf.interaction_functions.non_negative_linear_mortality(X, idx_A, idx_B, coefficient)
nemf.interaction_functions.non_negative_nutrient_limited_growth(X, idx_A, idx_B, growth_rate, half_saturation)[source]
nemf.interaction_functions.non_negative_remineralisation(X, idx_A, idx_B, coefficient)
nemf.interaction_functions.non_negative_stress_dependant_exudation(X, idx_A, idx_B, coefficient)
nemf.interaction_functions.nutrient_limited_growth(X, idx_A, idx_B, growth_rate, half_saturation)[source]

non-linear response with respect to destination/predator compartment

Similar to holling_type_II and is a reparameterization of holling II. The response with respect to the origin compartment ‘B’ is approximately linear for small ‘B’ and converges towards an upper limit governed by the ‘growth_rate’ for large ‘B’. For examples see: Examples

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • growth_rate (float) – first parameter of the interaction. governs the upper limit of the response.
  • half_saturation (float) – second parameter of the interaction. governs the slope of the response.
Returns:

df – change in the origin and destitnation compartment. Calculated by consumption_rate = ((hunting_rate * origin_compartment) / (1 + hunting_rate * food_processing_time * origin_compartment)) * destination_compartment
Return type:

float
nemf.interaction_functions.remineralisation(X, idx_A, idx_B, coefficient)
linear response with respect to both origin/pray and
destination/predator compartment.

For examples see: Examples

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • coefficient (float) – governs the slope of the linear response
Returns:

df – change in the origin and destitnation compartment. Calculated by coefficient * origin_compartment * destination_compartment
Return type:

float
nemf.interaction_functions.sloppy_feeding(holling_type, coeff, *args)[source]

calls holling_type functions with an extra “efficiency” coefficient. the inverse of the efficiency is then supposed to flow into a different compartment

nemf.interaction_functions.stress_dependant_exudation(X, idx_A, idx_B, coefficient)
linear response with respect to both origin/pray and
destination/predator compartment.

For examples see: Examples

Parameters:
  • X (np.array) – containing the current state of the contained quantity of each compartment
  • idx_A (integer) – index of the element representing the destination/predator compartment
  • idx_B (integer) – index of the element representing the origin/pray compartment
  • coefficient (float) – governs the slope of the linear response
Returns:

df – change in the origin and destitnation compartment. Calculated by coefficient * origin_compartment * destination_compartment
Return type:

float