CompartmentalSystems package¶

Subpackages¶

CompartmentalSystems.bins package

Submodules¶

CompartmentalSystems.BlockIvp module¶

class CompartmentalSystems.BlockIvp.BlockIvp(time_str: str, start_blocks: List[Tuple[str, numpy.ndarray]], functionss: List[List[Tuple[Callable, List[str]]]], disc_times: Tuple[float] = ())[source]¶

Bases: object

Helper class to build initial value systems from functions that operate on blocks of the state_variables.

block_solve(t_span, first_step=None, **kwargs)[source]¶

block_solve_functions(t_span, first_step=None, **kwargs)[source]¶

classmethod build_rhss(time_str: str, start_blocks: List[Tuple[str, numpy.ndarray]], functionss: List[List[Tuple[Callable, List[str]]]]) → Callable[[numpy.float64, numpy.ndarray], numpy.ndarray][source]¶

The function returns a function dot_X=f(t,X) suitable as the right-hand side for the ode solver scipy.solve_ivp from a collection of array valued functions that compute blocks of dot_X from time and blocks of X rather than from single equations.

A special application is the creation of block triangular systems, to integrate variables whose time derivative depends on the solution of an original system instantaneously along with it.

Assume that X_1(t) is the solution of the initial value problem (ivp)

ivp_1: dot_X_1=f_1(t,X), X_1(t_0)

and X_2(t) the solution of another ivp

ivp_2: dot_X_2=f_2(t,X_1,X_2), X_2(t_0) whose right-hand side depends on x_1.

Then we can obtain the solution of both ivps simultaneously by combining them into one.

(dot_X_1, dox_X_2)^t = (f_1(t,X_1),f_2(t,X_1,X_2))^t

For n instead of 2 variables one has: (dot_X_1, dot_X_2,…, dot_X_n)^t

= (f_1(t,X_1), f_2(t,X_1,X_2),…, f_n(t,X_1,…X_n))^t

For a full lower triangular system the block derivative dot_X_i depends on t and ALL the blocks X_1,…,X_i but often it will only depend on SOME of the previous blocks so that f_i has a considerably smaller argument list.

This function therefore allows to specify WHICH blocks the f_i depend on. Consider the following 5+2*2 = 9 -dimensional block diagonal example:

b_s=block_rhs(: time_str=’t’

,start_blocks=[(‘X1’,np.ones((5,1)),(‘X2’,np.ones((2,2)))] ,functions=[

((lambda x : x*2 ), [‘X1’] )

,((lambda t,x : t*x ), [‘t’ ,’X2’])

])

The first argument ‘time_str’ denotes the alias for the t argument to be used later in the signature of the block functions. The second argument ‘start_blocks’ describes the decomposition of X into blocks by a list of tuples of the form (‘Name’,array). The third argument ‘functions’ is a list of tuples of the function itself and the list of the names of its block arguments as specified in the ‘start_blocks’ argument. Order is important for the ‘start_blocks’ and the ‘functions’. It is assumed that the i-th function computes the derivative of the i-th block. The names of the blocks itself are arbitrary and have no meaning apart from their correspondence in the start_blocks and functions argument.

check_block_exists(block_name)[source]¶

CompartmentalSystems.BlockOde module¶

class CompartmentalSystems.BlockOde.BlockOde(time_str: str, block_names_and_shapes: List[Tuple[str, Tuple[int]]], functionss: List[List[Tuple[Callable, List[str]]]], disc_times: Tuple[int] = ())[source]¶

Bases: object

Helper class to build a system from functions that operate on blocks of the state_variables.

blockIvp(start_blocks: List[Tuple[str, numpy.ndarray]]) → CompartmentalSystems.BlockIvp.BlockIvp [source]¶: Extends the system to an initial value problem by adding startvalues for the blocks. It checks that the names of the names of the blocks coincide.

check_block_exists(block_name)[source]¶

CompartmentalSystems.Cache module¶

class CompartmentalSystems.Cache.Cache(keys, values, smr_hash)[source]¶

Bases: object

end_time_from_phi_ind(ind)[source]¶

classmethod from_file(filename)[source]¶

phi_ind(tau)[source]¶

save(filename)[source]¶

start_time_from_phi_ind(ind)[source]¶

CompartmentalSystems.cs_plotter module¶

class CompartmentalSystems.cs_plotter.CSPlotter(state_vector, inputs, outputs, internal_fluxes, pipe_colors={'linear': 'blue', 'no state dependence': 'red', 'nonlinear': 'green'}, visible_pool_names=True, pool_size_scale_factor=1, pool_color='blue', pool_alpha=0.3, pipe_alpha=0.5, connectionstyle='arc3, rad=0.1', arrowstyle='simple', mutation_scale=50, fontsize=24)[source]¶

Bases: object

legend(ax)[source]¶

plot_pools_and_fluxes(ax)[source]¶

class CompartmentalSystems.cs_plotter.Pool(x, y, size, pool_color, pool_alpha, pipe_alpha, connectionstyle, arrowstyle, mutation_scale)[source]¶

Bases: object

plot(ax)[source]¶

plot_input_flux(ax, color)[source]¶

plot_internal_flux_to(ax, pool_to, color)[source]¶

plot_name(ax, name, fontsize)[source]¶

plot_output_flux(ax, color)[source]¶

CompartmentalSystems.discrete_model_run module¶

exception CompartmentalSystems.discrete_model_run.DMRError[source]¶

Bases: Exception

Generic error occurring in this module.

class CompartmentalSystems.discrete_model_run.DiscreteModelRun(times, Bs, xs)[source]¶

Bases: object

acc_external_output_vector()[source]¶

acc_internal_flux_matrix()[source]¶

acc_net_external_input_vector()[source]¶

acc_net_external_output_vector()[source]¶

acc_net_internal_flux_matrix()[source]¶

age_densities_1_single_value_func(start_age_densities)[source]¶

age_densities_2_single_value_func()[source]¶

age_densities_single_value_func(start_age_densities)[source]¶

age_moment_vector(order, start_age_moments)[source]¶

age_quantiles_at_time(q, t, pools, start_age_densities)[source]¶

backward_transit_time_moment(order, start_age_moments)[source]¶

compute_start_age_moments(max_order, time_index=0)[source]¶

compute_start_m_factorial_moment(order, time_index=0)[source]¶

property dts¶: The lengths of the time intervals.

classmethod from_Bs_and_net_Us(start_values, times, Bs, net_Us)[source]¶: Bs State transition operators for one time step

classmethod from_SmoothModelRun(smr, data_times=None)[source]¶

classmethod from_fluxes(start_values, times, net_Us, net_Fs, net_Rs)[source]¶

classmethod load_from_file(filename)[source]¶

property net_Us¶

property nr_pools¶

pool_age_densities_func(start_age_densities, coarsity=1)[source]¶

classmethod reconstruct_B(x, F, R)[source]¶

classmethod reconstruct_Bs(xs, Fs, Rs)[source]¶

classmethod reconstruct_Bs_without_xs(start_values, Us, Fs, Rs)[source]¶

classmethod reconstruct_from_fluxes_and_solution(data_times, xs, Fs, Rs)[source]¶

save_to_file(filename)[source]¶

solve()[source]¶

start_age_densities_func()[source]¶

property start_values¶

system_age_moment(order, start_age_moments)[source]¶

CompartmentalSystems.discrete_model_run_14C module¶

class CompartmentalSystems.discrete_model_run_14C.DiscreteModelRun_14C(dmr, start_values_14C, net_Us_14C, decay_rate=0.0001209681)[source]¶

Bases: CompartmentalSystems.discrete_model_run.DiscreteModelRun

Construct and return a DiscreteModelRun_14C instance that: models the 14C component of the original model run.

Parameters

dmr (DiscreteModelRun) – original model run
start_values_14C (numpy.nd_array, nr_pools) – 14C start values.
Fa_func (func(t)) – returns atmospheric fraction to be multiplied with the input vector
rate (decay) – The decay rate to be used, defaults to 0.0001209681 (daily).

acc_net_external_input_vector_Delta_14C(alpha=None)[source]¶

acc_net_external_output_vector()[source]¶

acc_net_external_output_vector_Delta_14C(alpha=None)[source]¶

acc_net_internal_flux_matrix_Delta_14C(alpha=None)[source]¶

property external_output_vector¶

solve_Delta_14C(alpha=None)[source]¶

exception CompartmentalSystems.discrete_model_run_14C.Error[source]¶

Bases: Exception

Generic error occurring in this module.

CompartmentalSystems.discrete_model_run_with_gross_fluxes module¶

exception CompartmentalSystems.discrete_model_run_with_gross_fluxes.DMRError[source]¶

Bases: Exception

Generic error occurring in this module.

class CompartmentalSystems.discrete_model_run_with_gross_fluxes.DiscreteModelRunWithGrossFluxes(times, Bs, xs, gross_Us, gross_Fs, gross_Rs)[source]¶

Bases: CompartmentalSystems.discrete_model_run.DiscreteModelRun, CompartmentalSystems.model_run.ModelRun

acc_external_input_vector()[source]¶

acc_gross_external_input_vector()[source]¶

acc_gross_external_output_vector()[source]¶

acc_gross_internal_flux_matrix()[source]¶

classmethod from_SmoothModelRun(smr, data_times=None)[source]¶

classmethod reconstruct_from_fluxes_and_solution(data_times, xs, net_Fs, net_Rs, gross_Us, gross_Fs, gross_Rs)[source]¶

CompartmentalSystems.example_smooth_model_runs module¶

CompartmentalSystems.example_smooth_model_runs.critics()[source]¶

CompartmentalSystems.example_smooth_model_runs.emanuel_1()[source]¶

CompartmentalSystems.example_smooth_model_runs.nonlinear_two_pool()[source]¶

CompartmentalSystems.example_smooth_reservoir_models module¶

CompartmentalSystems.example_smooth_reservoir_models.critics(symbs)[source]¶

CompartmentalSystems.example_smooth_reservoir_models.emanuel(symbs)[source]¶

CompartmentalSystems.example_smooth_reservoir_models.minimal(symbs)[source]¶

CompartmentalSystems.example_smooth_reservoir_models.nonlinear_two_pool(symbs)[source]¶

CompartmentalSystems.helpers_reservoir module¶

CompartmentalSystems.helpers_reservoir.F_Delta_14C(C12, C14, alpha=None)[source]¶

CompartmentalSystems.helpers_reservoir.MH_sampling(N, PDF, start=1.0)[source]¶

CompartmentalSystems.helpers_reservoir.arrange_subplots(n)[source]¶

CompartmentalSystems.helpers_reservoir.array_integration_by_ode(f: Callable[[float], numpy.ndarray], a: float, b: float, *args, **kwargs) → numpy.ndarray[source]¶

CompartmentalSystems.helpers_reservoir.array_integration_by_values(f: Callable[[float], numpy.ndarray], taus: numpy.ndarray, *args, **kwargs) → numpy.ndarray[source]¶

CompartmentalSystems.helpers_reservoir.array_quad_result(f: Callable[[float], numpy.ndarray], a: float, b: float, epsrel=0.001, *args, **kwargs) → numpy.ndarray[source]¶

CompartmentalSystems.helpers_reservoir.check_parameter_dict_complete(model, parameter_dict, func_set)[source]¶

Check if the parameter set the function set are complete: to enable a model run.

Parameters

model (SmoothReservoirModel) – The reservoir model on which the model run bases.
parameter_dict (dict) – {x: y} with x being a SymPy symbol and y being a numerical value.
func_set (dict) – {f: func} with f being a SymPy symbol and func being a Python function. Defaults to dict().

Returns

set of free symbols, parameter_dict is complete if: free_symbols is the empty set

Return type

free_symbols (set)

CompartmentalSystems.helpers_reservoir.const_of_t_maker(const)[source]¶

CompartmentalSystems.helpers_reservoir.custom_lru_cache_wrapper(maxsize=None, typed=False, stats=False)[source]¶

CompartmentalSystems.helpers_reservoir.deprecation_warning(txt)[source]¶

CompartmentalSystems.helpers_reservoir.draw_rv(CDF)[source]¶

CompartmentalSystems.helpers_reservoir.f_of_t_maker(sol_funcs, ol)[source]¶

CompartmentalSystems.helpers_reservoir.factor_out_from_matrix(M)[source]¶

CompartmentalSystems.helpers_reservoir.flux_dict_string(d, indent=0)[source]¶

CompartmentalSystems.helpers_reservoir.func_subs(t, Func_expr, func, t0)[source]¶: returns the function part_func where part_func(_,_,…) =func(_,t=t0,_..) (func partially applied to t0) The position of argument t in the argument list is found by examining the Func_expression argument. :param t: the symbol to be replaced by t0 :type t: sympy.symbol :param t0: the value to which the function will be applied :type t0: value :param func: a python function :type func: function :param Func_exprs: An expression for an undefined Function :type Func_exprs: sympy.Function

CompartmentalSystems.helpers_reservoir.generalized_inverse_CDF(CDF, u, start_dist=0.0001, tol=1e-08)[source]¶

CompartmentalSystems.helpers_reservoir.has_pw(expr)[source]¶

CompartmentalSystems.helpers_reservoir.integrate_array_func_for_nested_boundaries(f: Callable[[float], numpy.ndarray], integrator: Callable[[Callable[[float], numpy.ndarray], float, float], numpy.ndarray], tuples: Sequence[Tuple[float, float]]) → Sequence[float][source]¶

CompartmentalSystems.helpers_reservoir.is_DiracDelta(expr)[source]¶: Check if expr is a Dirac delta function.

CompartmentalSystems.helpers_reservoir.is_compartmental(M)[source]¶

CompartmentalSystems.helpers_reservoir.jacobian(vec, state_vec)[source]¶

CompartmentalSystems.helpers_reservoir.load_csv(filename)[source]¶

CompartmentalSystems.helpers_reservoir.make_cut_func_set(func_set)[source]¶

CompartmentalSystems.helpers_reservoir.melt(ndarr, identifiers=None)[source]¶

CompartmentalSystems.helpers_reservoir.net_Fs_from_discrete_Bs_and_xs(Bs, xs)[source]¶

CompartmentalSystems.helpers_reservoir.net_Rs_from_discrete_Bs_and_xs(Bs, xs, decay_corr=None)[source]¶

CompartmentalSystems.helpers_reservoir.net_Us_from_discrete_Bs_and_xs(Bs, xs)[source]¶

CompartmentalSystems.helpers_reservoir.numerical_function_from_expression(expr, tup, parameter_dict, func_set)[source]¶

CompartmentalSystems.helpers_reservoir.numerical_rhs(state_vector, time_symbol, rhs, parameter_dict, func_dict)[source]¶

CompartmentalSystems.helpers_reservoir.numerical_rhs_old(state_vector, time_symbol, rhs, parameter_dict, func_set, times)[source]¶

CompartmentalSystems.helpers_reservoir.numsol_symbolic_system_old(state_vector, time_symbol, rhs, parameter_dict, func_set, start_values, times)[source]¶

CompartmentalSystems.helpers_reservoir.numsol_symbolical_system(state_vector, time_symbol, rhs, parameter_dicts, func_dicts, start_values, times, disc_times=())[source]¶

CompartmentalSystems.helpers_reservoir.parse_input_function(u_i, time_symbol)[source]¶

Return an ordered list of jumps in the input function u.

Parameters: u (SymPy expression) – input function in \(\dot{x} = B\,x + u\)
Returns: ascending list of jumps in u

CompartmentalSystems.helpers_reservoir.phi_tmax(s, t_max, block_ode, x_s, x_block_name, phi_block_name)[source]¶

CompartmentalSystems.helpers_reservoir.print_quantile_error_statisctics(qs_ode, qs_pi)[source]¶

CompartmentalSystems.helpers_reservoir.save_csv(filename, melted, header)[source]¶

CompartmentalSystems.helpers_reservoir.stochastic_collocation_transform(M, CDF)[source]¶

CompartmentalSystems.helpers_reservoir.stride(data, strides)[source]¶

CompartmentalSystems.helpers_reservoir.tup2str(tup)[source]¶

CompartmentalSystems.helpers_reservoir.warning(txt)[source]¶

CompartmentalSystems.helpers_reservoir.x_phi_ode(srm, parameter_dicts, func_dicts, x_block_name='x', phi_block_name='phi', disc_times=())[source]¶

CompartmentalSystems.helpers_reservoir.x_phi_tmax(s, t_max, block_ode, x_s, x_block_name, phi_block_name)[source]¶

CompartmentalSystems.helpers_reservoir.x_tmax(s, t_max, block_ode, x_s, x_block_name, phi_block_name)[source]¶

CompartmentalSystems.model_run module¶

class CompartmentalSystems.model_run.ModelRun[source]¶

Bases: object

abstract acc_gross_external_input_vector()[source]¶

abstract acc_gross_external_output_vector()[source]¶

abstract acc_gross_internal_flux_matrix()[source]¶

abstract acc_net_external_input_vector()[source]¶

abstract acc_net_external_output_vector()[source]¶

abstract acc_net_internal_flux_matrix()[source]¶

abstract property dts¶

abstract property nr_pools¶

abstract solve(alternative_start_values=None)[source]¶

CompartmentalSystems.model_run.plot_attributes(mrs, file_name)[source]¶

CompartmentalSystems.model_run.plot_stocks_and_fluxes(mrs, file_name, labels=None)[source]¶

CompartmentalSystems.myOdeResult module¶

CompartmentalSystems.myOdeResult.get_sub_t_spans(t_span, disc_times)[source]¶

class CompartmentalSystems.myOdeResult.myOdeResult(y, t, sol)[source]¶: Bases: scipy.integrate._ivp.ivp.OdeResult

CompartmentalSystems.myOdeResult.solve_ivp_pwc(rhss, t_span, y0, disc_times=(), **kwargs)[source]¶

CompartmentalSystems.picklegzip module¶

CompartmentalSystems.picklegzip.dump(object, filename, protocol=- 1)[source]¶: Save an object to a compressed disk file. Works well with huge objects.

CompartmentalSystems.picklegzip.load(filename)[source]¶: Loads a compressed object from disk.

CompartmentalSystems.pwc_model_run module¶

exception CompartmentalSystems.pwc_model_run.Error[source]¶

Bases: Exception

Generic error occurring in this module.

class CompartmentalSystems.pwc_model_run.PWCModelRun(model, parameter_dicts, start_values, times, disc_times, func_dicts=None)[source]¶

Bases: CompartmentalSystems.model_run.ModelRun

B_func(vec_sol_func=None)[source]¶

Phi(T, S)[source]¶

acc_gross_external_input_vector(data_times=None)[source]¶

acc_gross_external_output_vector(data_times=None)[source]¶

acc_gross_internal_flux_matrix(data_times=None)[source]¶

acc_net_external_input_vector(data_times=None)[source]¶

acc_net_external_output_vector(data_times=None)[source]¶

acc_net_internal_flux_matrix(data_times=None)[source]¶

property boundaries¶

property dts¶

external_input_flux_funcss()[source]¶

external_input_vector_func(cut_off=True)[source]¶

external_output_flux_funcss()[source]¶

fake_discretized_Bs(data_times=None)[source]¶

initialize_state_transition_operator_cache(lru_maxsize, lru_stats=False, size=1)[source]¶

internal_flux_funcss()[source]¶

join_flux_funcss_rc(flux_funcss, key)[source]¶

join_functions_rc(funcs)[source]¶

myhash()[source]¶: Compute a hash considering SOME but NOT ALL properties of a model run. The function’s main use is to detect saved state transition operator cashes that are no longer compatible with the model run object that wants to use them. This check is useful but NOT COMPREHENSIVE.

property nr_intervals¶

property nr_pools¶

sol_funcs()[source]¶

solve(alternative_start_values=None)[source]¶

solve_func(alternative_start_values=None)[source]¶

CompartmentalSystems.pwc_model_run_14C module¶

exception CompartmentalSystems.pwc_model_run_14C.Error[source]¶

Bases: Exception

Generic error occurring in this module.

class CompartmentalSystems.pwc_model_run_14C.PWCModelRun_14C(pwc_mr, start_values_14C, Fa_func, decay_rate=0.0001209681)[source]¶

Bases: CompartmentalSystems.pwc_model_run.PWCModelRun

acc_gross_external_input_vector_Delta_14C(data_times=None, alpha=None)[source]¶

acc_gross_external_output_vector(data_times=None)[source]¶

acc_gross_external_output_vector_Delta_14C(data_times=None, alpha=None)[source]¶

acc_gross_internal_flux_matrix_Delta_14C(data_times=None, alpha=None)[source]¶

acc_net_external_input_vector_Delta_14C(data_times=None, alpha=None)[source]¶

acc_net_external_output_vector(data_times=None)[source]¶

acc_net_external_output_vector_Delta_14C(data_times=None, alpha=None)[source]¶

acc_net_internal_flux_matrix_Delta_14C(data_times=None, alpha=None)[source]¶

external_output_flux_funcss()[source]¶

property external_output_vector¶

solve_Delta_14C(alpha=None)[source]¶

CompartmentalSystems.pwc_model_run_14C.pfile_C14Atm_NH()[source]¶

CompartmentalSystems.pwc_model_run_fd module¶

class CompartmentalSystems.pwc_model_run_fd.PWCModelRunFD(Bs, us, pwc_mr)[source]¶

Bases: CompartmentalSystems.model_run.ModelRun

B_func(vec_sol_func=None)[source]¶

classmethod B_pwc(times, Bs)[source]¶

acc_gross_external_input_vector()[source]¶

acc_gross_external_output_vector()[source]¶

acc_gross_internal_flux_matrix()[source]¶

acc_net_external_input_vector()[source]¶

acc_net_external_output_vector()[source]¶

acc_net_internal_flux_matrix()[source]¶

classmethod create_B_generic(Bs)[source]¶

classmethod create_srm_generic(time_symbol, Bs, us)[source]¶

classmethod create_u_generic(us)[source]¶

property dts¶

external_input_vector_func(cut_off=True)[source]¶

fake_discretized_Bs(data_times=None)[source]¶

classmethod from_Bs_and_us(time_symbol, data_times, start_values, Bs, us)[source]¶

classmethod from_gross_fluxes(time_symbol, data_times, start_values, gross_Us, gross_Fs, gross_Rs)[source]¶

property model¶

property nr_pools¶

classmethod reconstruct_B_surrogate(dt, x0, gross_U, gross_F, gross_R, B0)[source]¶

classmethod reconstruct_Bs(times, start_values, gross_Us, gross_Fs, gross_Rs)[source]¶

solve(alternative_start_values=None)[source]¶

property start_values¶

property times¶

classmethod u_pwc(times, us)[source]¶

exception CompartmentalSystems.pwc_model_run_fd.PWCModelRunFDError[source]¶

Bases: Exception

Generic error occurring in this module.

CompartmentalSystems.smooth_model_run module¶

Module for numerical treatment of piece-wise continuous reservoir models.

An abstract SmoothReservoirModel is filled with life by giving initial values, a parameter set, a time grid, and potentially additional involved functions to it.

The model can then be run and as long as the model is linear, based on the state transition operator age and transit time distributions can be computed.

Nonlinear models can be linearized along a solution trajectory.

Counting of compartment/pool/reservoir numbers start at zero and the total number of pools is \(d\).

exception CompartmentalSystems.smooth_model_run.Error[source]¶

Bases: Exception

Generic error occurring in this module.

class CompartmentalSystems.smooth_model_run.SmoothModelRun(model, parameter_dict, start_values, times, func_set=None)[source]¶

Bases: CompartmentalSystems.model_run.ModelRun

Class for a model run based on a SmoothReservoirModel.

model¶

The reservoir model on which the model run bases.

Type: SmoothReservoirModel

parameter_dict¶

{x: y} with x being a SymPy symbol and y being a numerical value.

Type: dict

start_values¶

The vector of start values.

Type: numpy.array

times¶

The time grid used for the simulation. Typically created by numpy.linspace.

Type: numpy.array

func_set¶

{f: func} with f being a SymPy symbol and func being a Python function. Defaults to dict().

Type: dict

Pool counting starts with 0. In combined structures for pools and system, the system is at the position of a (d+1) st pool.

B_func(vec_sol_func=None)[source]¶

Phi(T, S)[source]¶

Phi_func()[source]¶

acc_gross_external_input_vector(data_times=None)[source]¶

Return the grid of accumulated external input vectors.

Returns: len(times) x nr_pools
Return type: numpy.ndarray

acc_gross_external_output_vector(data_times=None)[source]¶

Return the vectors of accumulated external outputs.

Returns: len(times)-1 x nr_pools
Return type: numpy.ndarray

acc_gross_internal_flux_matrix(data_times=None)[source]¶

Return the grid of flux matrices.

Returns: len(times) x nr_pools x nr_pools
Return type: numpy.ndarray

acc_net_external_input_vector(data_times=None)[source]¶

acc_net_external_output_vector(data_times=None)[source]¶

acc_net_internal_flux_matrix(data_times=None)[source]¶

add_constant_age_distribution_surface_plotly(fig, opacity=0.7, index=0)[source]¶

Add a grey and transparent density surface to an existing Plotly density plot.

If index is not specified it is assumed to be 0 and the values correspond to the first time in the times porperty of the model run (the age distribution at the beginning) and repeated for all times. The plotted surface represents an age distribution that is constant in time. It is intended to increase the visibility of changes in the age distribution with time. Note that this constant age distribution does NOT necessarily correspond to a possible (constant) development of the system. This would only be true if the system was in equilibrium and the age distribution was the equilibrium age distribution. While this special case is a very interesting application this function does not assertain that such an equlibrium situation is even possible.

Parameters

fig (Plotly figure) – The existing density plot to which the surface is added.
opacity (between 0 and 1, optional) – The opacity of the new surface. Defaults to 0.9. Unfortunately, the opacity option does not seem to work properly.
index (int, optional) – The time index from which the age distribution data is taken. Defaults to 0 such that the constant distribution is computed at time \(t_0\).

Returns

None. Instead fig is changed in place.

add_line_to_density_plot_plotly(fig, data, color, name, time_stride=1, width=5, on_surface=True, bottom=True, legend_on_surface=False, legend_bottom=False)[source]¶

Add a line to an existing Plotly density plot.

Parameters

fig (Plotly figure) – Contains already a density plot to which the new line is added.
data (numpy.array len(times)) – The age data of the new line.
color (#RRGGBB) – The color of the new line.
name (str) – The name of the new line for the legend.
time_stride (int, optional) – Coarsity of the plot in the time direction to save memory. Defaults to 1 meaning that all times are plotted and no memory is saved.
width (int, optional) – Width of the new line. Defaults to 5.
on_surface (bool, optional) – If True, a new line with the given age data is plotted on top of the given density. Defaults to True.
bottom (bool optional) – If True, a new line with the given age data is plotted in the xy-plane. Defaults to True.
legend_on_surface (bool, optional) – If True, the line on the surface is mentioned in the legend. Has no effect if on_surface is False. Defaults to False.
legend_bottom (bool, optional) – If True, the line in the xy-plane is mentioned in the legend. Has no effect if bottom is False. Defaults to False.

Returns

None. Instead fig is changed in place.

age_densities(pool_age_densities, system_age_density)[source]¶

Combine pool and system age densities to one numpy.array.

Parameters

pool_age_densites (numpy.ndarray len(ages) x len(times) x nr_pools) – The pool age density values.
system_age_density (numpy.ndarray len(ages) x len(times)) – The system age density values.

Returns

(len(ages) x len(times) x (nr_pools+1)). The system age density values are appended to the end of the pool density values (system = pool d+1 with d = nr_pools).

Return type

numpy.ndarray

age_moment_vector(order, start_age_moments=None)[source]¶

Compute the order th pool age moment vector over the time grid by an ODE system.

This function solves an ODE system to obtain the pool age moments very fast. If the system has empty pools at the beginning, the semi-explicit formula is used until all pools are non-empty. Then the ODE system starts.

Parameters

order (int) – The order of the pool age moments to be computed.
start_age_moments (numpy.ndarray order x nr_pools, optional) – Given initial age moments up to the order of interest. Can possibly be computed by moments_from_densities(). Defaults to None assuming zero initial ages.

Returns

len(times) x nr_pools. The order th pool age moments over the time grid.

Return type

numpy.ndarray

age_moment_vector_from_densities(order, start_age_densities)[source]¶

Compute the order th moment of the pool ages by integration.

This function is extremely slow, since for each pool the integral over the density is computed based on the singe-valued functions. It is implemented only for the sake of completeness and to test the results obtained by faster methods.

Parameters

order (int) – The order of the moment to be computed.
start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\).

Returns

len(times) x nr_pools. Contains the order th moment of the pool ages over the time grid.

Return type

numpy.ndarray

age_moment_vector_semi_explicit(order, start_age_moments=None, times=None)[source]¶

Compute the order th moment of the pool ages by a semi-explicit formula.

This function bases on a semi-explicit formula such that no improper integrals need to be computed.

Parameters

order (int) – The order of the age moment to be computed.
start_age_moments (numpy.ndarray order x nr_pools, optional) – Given initial age moments up to the order of interest. Can possibly be computed by moments_from_densities(). Defaults to None assuming zero initial ages.
times (numpy.array, optional) – Time grid. Defaults to None and the original time grid is used.

Returns

len(times) x nr_pools. The order th pool age moments over the time grid.

Return type

numpy.ndarray

apply_to_forward_transit_time_simulation(f_dict={'mean': <function mean>}, N=10000, M=2, k=5, MH=False)[source]¶

This is just a tentative approach.

To be honest, the problem of an unkown infinite future cannot be solved this way, since the densities used to simulate forward transit time random variables are biased by the cut-off at the end of the time grid.

apply_to_forward_transit_time_simulation_its(f_dict, times, N=1000, k=5)[source]¶

This is just a tentative approach.

To be honest, the problem of an unkown infinite future cannot be solved this way, since the densities used to simulate forward transit time random variables are biased by the cut-off at the end of the time grid.

backward_transit_time_density(pool_age_densities)[source]¶

Compute the backward transit time based on given pool age densities.

The result is obtained by computing a weighted sum of the pool age densities according to output rates.

Parameters: pool_age_densites (numpy.ndarray len(ages) x len(times) x nr_pools) – The pool age density values.
Returns: len(ages) x len(times). Mass leaving the system with the respective age at the respective time.
Return type: numpy.ndarray

backward_transit_time_density_single_value_func(start_age_densities)[source]¶

Return a function that returns a single value for the backward transit time density.

Parameters: start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\).
Returns: Python function p_sv: p_sv(a, t) returns the mass that leaves the system at time t with age a.

backward_transit_time_moment(order, start_age_moments=None)[source]¶

Compute the order th backward transit time moment based on the age_moment_vector().

Parameters

order (int) – The order of the backward transit time moment that is to be computed.
start_age_moments (numpy.ndarray order x nr_pools, optional) – Given initial age moments up to the order of interest. Can possibly be computed by moments_from_densities(). Defaults to None assuming zero initial ages.

Returns

The order th backward transit time moment over the time grid.

Return type

numpy.array

backward_transit_time_moment_from_density(order, start_age_densities)[source]¶

Compute the order th backward transit time moment based on an improper integral over the density.

This function is extremely slow and implemented only for the sake of completeness and for testing results from faster approaches.

Parameters

order (int) – The order of the backward transit time moment that is to be computed.
start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\).

Returns

The order th backward transit time moment over the time grid.

Return type

numpy.array

combine_pools_and_system_values(pools_values, system_values)[source]¶

Append the system values to the pool values as if they belonged to another pool.

Parameters

pools_values (numpy.ndarray len(times) x nr_pools) – The values to be saved over the time-pool grid.
system_values (numpy.array) – The system values to be saved over the time grid.

Returns

len(times) x (nr_pools+1) The pool and system values over the time-pool grid with the system added at the end as another pool.

Return type

numpy.ndarray

cumulative_backward_transit_time_distribution_single_value(start_age_densities=None, F0=None)[source]¶

Return a function for the cumulative backward transit time distribution.

Parameters

start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\). Defaults to None.
F0 (Python function) – A function of age that returns a numpy.array containing the masses with age less than or equal to the age at time \(t_0\). Defaults to None.

Raises

Error – If both start_age_densities and F0 are None. One must be given. It is fastest to provide F0, otherwise F0 will be computed by numerical integration of start_age_densities.

Returns

Python function F_sv: F_sv(a, t) is the mass leaving the system at time t with age less than or equal to a.

cumulative_forward_transit_time_distribution_single_value_func(cut_off=True)[source]¶

Return a function for the cumulative forward transit time distribution.

Parameters: cut_off (bool, optional) – If True, no density values are going to be computed after the end of the time grid, instead numpy.nan will be returned. Defaults to True. False might lead to unexpected behavior.
Returns: Python function F_sv: F_sv(a, t) is the mass leaving the system at time t+a with age less than or equal to a.

cumulative_pool_age_distributions_single_value(start_age_densities=None, F0=None)[source]¶

Return a function for the cumulative pool age distributions.

Parameters

start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\). Defaults to None.
F0 (Python function) – A function of age that returns a numpy.array containing the masses with age less than or equal to the age at time \(t_0\). Defaults to None.

Raises

Error – If both start_age_densities and F0 are None. One must be given. It is fastest to provide F0, otherwise F0 will be computed by numerical integration of start_age_densities.

Returns

Python function F_sv: F_sv(a,t) is the vector of pool masses (numpy.array) with age less than or equal to a at time t.

cumulative_system_age_distribution_single_value(start_age_densities=None, F0=None)[source]¶

Return a function for the cumulative system age distribution.

Parameters

start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\). Defaults to None.
F0 (Python function) – A function of age that returns a numpy.array containing the masses with age less than or equal to the age at time \(t_0\). Defaults to None.

Raises

Error – If both start_age_densities and F0 are None. One must be given. It is fastest to provide F0, otherwise F0 will be computed by numerical integration of start_age_densities.

Returns

Python function F_sv: F_sv(a, t) is the mass in the system with age less than or equal to a at time t.

density_values(density_sv, values)[source]¶

Compute the density value over the time grid at ages given in values.

This function can be used to obtain the density values at mean or median values to draw a curve on top of the density surface. But actually this is now implemented in a much faster way based on the surface itself.

Parameters

density_sv (Python function) – A function that takes a, t as arguments and returns density value with age a at time t. Potentially coming from system_age_density_single_value().
values (numpy.array) – The ages over the time grid at which the density values are to be computed.

Returns

The density values over the time grid based on the given age values.

Return type

numpy.array

density_values_for_pools(pool_densities_sv, pool_age_values)[source]¶

Compute the pool age densities over the time grid at ages given in pool_age_values.

This function can be used to obtain the density values at mean or median values to draw a curve on top of the density surface. But actually this is now implemented in a much faster way based on the surface itself.

Parameters

pool_age_densites_sv (Python function) – A function that takes a, t as arguments and returns a vector of pool contents with mass a at time t. Potentially coming from pool_age_densities_single_value().
pool_age_values (numpy.ndarray len(times) x nr_pools) – The ages over the time-pool grid at which the density values are to be computed.

Returns

(len(times) x nr_pools) The pool density values over the time-pool grid based on the given age values.

Return type

numpy.ndarray

static distribution_quantile(quantile, F, norm_const=None, start_value=None, method='brentq', tol=1e-08)[source]¶

Return distribution quantile (one single value) of a given distribution.

The compuation is done by computing the generalized inverse of the respective cumulative distribution using a nonlinear root search algorithm.

Parameters

quantile (between 0 and 1) – The relative share of mass that is considered to be left of the computed value. A value of 0.5 leads to the computation of the median of the distribution.
F (Python function) – A function of age a that returns the mass that is of age less than or equal to a.
norm_const (numpy.array, optional) – The amount of total mass of the distribution. Defaults to one assuming a given probability distribution.
start_value (float, optional) – A start value for the nonlinear search. A good values is slighty greater than the solution value. Defaults to zero.
method (str) – The method that is used for finding the roots of a nonlinear function. Either ‘brentq’ or ‘newton’. Defaults to ‘brentq’.
tol (float) – The tolerance used in the numerical root search algorithm. A low tolerance decreases the computation speed tremendously, so a value of 1e-01 might already be fine. Defaults to 1e-08.

Returns

The computed quantile value of the distribution.

Return type

float

distribution_quantiles(quantile, F_sv, norm_consts=None, start_values=None, times=None, method='brentq', tol=1e-08)[source]¶

Return distribution quantiles over the time grid of a given distribution.

The compuation is done by computing the generalized inverse of the respective cumulative distribution using a nonlinear root search algorithm. Depending on how slowly the cumulative distribution can be computed, this can take quite some time.

Parameters

quantile (between 0 and 1) – The relative share of mass that is considered to be left of the computed value. A value of 0.5 leads to the computation of the median of the distribution.
F_sv (Python function) – A function of age a and time t that returns the mass that is of age less than or equal to a at time t.
norm_consts (numpy.array, optional) – An array over the time grid of total masses over all ages. Defaults to an array of ones assuming given probability distributions.
start_values (numpy.array, optional) – An array over the time grid of start values for the nonlinear search. Good values are slighty greater than the solution values. Must have the same length as times. Defaults to an array of zeros.
times (numpy.array, optional) – Time grid on which to compute the quantiles. Defaults to None in which case the orignal time grid is used.
method (str) – The method that is used for finding the roots of a nonlinear function. Either ‘brentq’ or ‘newton’. Defaults to ‘brentq’.
tol (float) – The tolerance used in the numerical root search algorithm. A low tolerance decreases the computation speed tremendously, so a value of 1e-01 might already be fine. Defaults to 1e-08.

Returns

The computed quantile values over the time grid.

Return type

numpy.array

property dts¶: The lengths of the time intervals.

external_input_flux_funcs()[source]¶

Return a dictionary of the external input fluxes. The resulting functions base on sol_funcs and are linear interpolations.

Returns: {key: func} with key representing the pool which receives the input and func a function of time that returns a float.
Return type: dict

property external_input_vector¶

Return the grid of external input vectors.

Returns: len(times) x nr_pools
Return type: numpy.ndarray

external_input_vector_func(cut_off=True)[source]¶

Return a vector valued function for the external inputs.

Returns: Python function u: u(t) is a numpy.array containing the external inputs at time t. Note: If the required (future) values for the input exceed the maximum of times they are assumed to be zero if cut_off is True. If cut_off is False then the input function is assumed to be valid everywhere which might be dangerous if they are extrapolated from data.

external_output_flux_funcs()[source]¶

Return a dictionary of the external output fluxes.

Returns: {key: func} with key representing the pool from which the output comes and func a function of time that returns a float.
Return type: dict

property external_output_vector¶

Return the grid of external output vectors.

Returns: len(times) x nr_pools
Return type: numpy.ndarray

fake_discretized_Bs(data_times=None)[source]¶

fake_gross_discretized_output(data_times)[source]¶

fake_net_discretized_output(data_times)[source]¶

forward_transit_time_density_func(cut_off=True, times=None)[source]¶

Return a function based on an age array for the forward transit time density.

Parameters

cut_off (bool, optional) – If True, no density values are going to be computed after the end of the time grid, instead numpy.nan will be returned. Defaults to True and False might lead to unexpected behavior.
times (numpy.array, optional) – Time grid. Defaults to None and the original time grid is used.

Returns

Python function p: p(ages) is a numpy.ndarray len(ages) x len(times) that gives the mass that will leave the system with the respective age when it came in at time t, where ages is a numpy.array.

forward_transit_time_density_single_value_func(cut_off=True, my_B_func=None)[source]¶

Return a function that returns a single value for the forward transit time density.

Parameters: cut_off (bool, optional) – If True, no density values are going to be computed after the end of the time grid, instead numpy.nan will be returned. Defaults to True and False might lead to unexpected behavior.
Returns: Python function p_sv: p_sv(a, t) returns how much mass will leave the system with age a when it came in at time t.

initialize_state_transition_operator_cache(lru_maxsize, lru_stats=False, size=1)[source]¶

internal_flux_funcs()[source]¶

Return a dictionary of the internal fluxes.

Returns: {key: func} with key=(pool_from, pool_to) representing the pools involved and func a function of time that returns a float.
Return type: dict

linearize()[source]¶

Return a linearized SmoothModelRun instance.

Linearization happens along the solution trajectory. Only for linear systems all functionality is guaranteed, this is why nonlinear systems should be linearized first.

Returns: A linearized version of the original SmoothModelRun, with the solutions now being part of func_set.
Return type: SmoothModelRun

linearize_old()[source]¶

Return a linearized SmoothModelRun instance.

Linearization happens along the solution trajectory. Only for linear systems all functionality is guaranteed, this is why nonlinear systems should be linearized first.

Returns: A linearized version of the original SmoothModelRun, with the solutions now being part of func_set.
Return type: SmoothModelRun

load_density_csv(filename, ages, times=None)[source]¶

Load density values from a csv file.

Attention: It is assumed that the data were saved before with the very same ages, times, and pools.

Parameters

filename (str) – The csv file from which the data are to be read.
ages (numpy.array) – The ages corresponding to the age indices. What is needed here is in fact only the length of the age grid.
times (numpy.array, optional) – The ages corresponding to the time indices. What is needed here is in fact only the length of the time grid. Defaults to None and the original times are being used.

Returns

len(ages) x len(times) The density values over the ages-times grid.

Return type

numpy.ndarray

load_pools_and_system_densities_csv(filename, ages)[source]¶

Load pool and system age densities from a csv file.

Attention: It is assumed that the data were saved before with the very same ages, times, and pools. Furthermore, it is assumed that the system value always follows the pool values.

Parameters

filename (str) – The csv file from which the data are to be read.
ages (numpy.array) – The ages corresponding to the age indices. What is needed here is in fact only the length of the age grid.

Returns

len(ages) x len(times) x (nr_pools+1) The density values for the pools and the system over the ages-times-(pools+system) grid.

Return type

numpy.ndarray

load_pools_and_system_value_csv(filename)[source]¶

Load pool and system values from a csv file.

Values could be the mean/median age, for example. One dimension less than a density.

Attention: It is assumed that the data were saved before with the very same ages, times, and pools. Furthermore, it is assumed that the system value always follows the pool values.

Parameters: filename (str) – The csv file from which the data are to be read.
Returns: len(times) x (nr_pools+1) The values for the pools and the system over the times-(pools+system) grid.
Return type: numpy.ndarray

load_state_transition_operator_cache(filename)[source]¶

load_value_csv(filename)[source]¶

static moments_from_densities(max_order, densities)[source]¶

Compute the moments up to max_order of the given densities.

Parameters

max_order (int) – The highest order up to which moments are to be computed.
densities (numpy.array) – Each entry is a Python function of one variable (age) that represents a probability density function.

Returns

moments x pools, containing the moments of the given densities.

Return type

numpy.ndarray

myhash()[source]¶: Compute a hash considering SOME but NOT ALL properties of a model run. The function’s main use is to detect saved state transition operator cashes that are no longer compatible with the model run object that wants to use them. This check is useful but NOT COMPREHENSIVE.

property no_input_model¶

property nr_pools¶

Return the number of pools involved in the model.

Type: int

property output_rate_vector¶

Return the grid of output rate vectors.

Returns: len(times) x nr_pools, solution/output_vector
Return type: numpy.ndarray

output_rate_vector_at_t(t)[source]¶

Return a vector of output rates at time t.

Parameters: t (float) – The time at which the output rates are computed.
Returns: The ith entry contains the output rate of pool i at time t.
Return type: numpy.array

output_vector_func(t)[source]¶

Return a vector of the external output fluxes at time t.

Returns: The i th entry is the output from pool i at time t.
Return type: numpy.array

plot_3d_density_plotly(title, density_data, ages, age_stride=1, time_stride=1)[source]¶

Create a 3-dimendional density plot with Plotly.

The colors are created such that they are constant along the age-time diagonal. Thus, equal colors mark equal entry time into the system.

Parameters

title (str) – The title of the new figure.
density_data (numpy.ndarray len(ages) x len(times)) – The density data to be plotted.
age_stride (int, optional) – Coarsity of the plot in the age direction to save memory. Defaults to 1 meaning that all times are plotted and no memory is saved.
time_stride (int, optional) – Coarsity of the plot in the time direction to save memory. Defaults to 1 meaning that all times are plotted and no memory is saved.

Returns

Plotly figure.

plot_external_input_fluxes(fig, fontsize=10)[source]¶

Plot all external inpput fluxes.

For each external input flux a new subplot is added.

Parameters

fig (Matplotlib figure) – The fig to which the subplots are added.
fontsize (float, optional) – Defaults to 10.

Returns

None. Instead fig is changed in place.

plot_external_output_fluxes(fig, fontsize=10)[source]¶

Plot all external output fluxes.

For each external output flux a new subplot is added.

Parameters

fig (Matplotlib figure) – The fig to which the subplots are added.
fontsize (float, optional) – Defaults to 10.

Returns

None. Instead fig is changed in place.

plot_internal_fluxes(fig, fontsize=10)[source]¶

Plot all internal fluxes.

For each internal flux a new subplot is added.

Parameters

fig (Matplotlib figure) – The fig to which the subplots are added.
fontsize (float, optional) – Defaults to 10.

Returns

None. Instead fig is changed in place.

plot_mean_ages(fig, start_mean_ages)[source]¶

Plot the time evolution of the mean ages for all pools and the system.

For each pool and the system a separate subplot is created.

Parameters

fig (Matplotlib figure) – The fig to which the subplots are added.
start_mean_ages (numpy.array) – Contains the start mean ages of the pools at time \(t_0\).

Returns

None. Instead fig is changed in place.

plot_mean_backward_transit_time(ax, start_mean_ages)[source]¶

Plot the time evolution of the mean backward transit time.

For each pool and the system a separate subplot is created.

Parameters

ax (Matplotlib axis) – The ax onto which the plot is done.
start_mean_ages (numpy.array) – Contains the start mean ages of the pools at time \(t_0\).

Returns

None. Instead ax is changed in place.

plot_phase_plane(ax, i, j, fontsize=10)[source]¶

Plot one single phase plane.

Parameters

ax (Matplotlib axis) – The axis onto which the phase plane is plotted.
i (int) – The numbers of the pools for which the phase plane is plotted.
j (int) – The numbers of the pools for which the phase plane is plotted.
fontsize (float, optional) – Defaults to 10.

Returns

None. Instead ax is changed in place.

plot_phase_planes(fig, fontsize=10)[source]¶

Plot all phase planes.

For each (i,j)-phase plane a new subplot is added.

Parameters

fig (Matplotlib figure) – The fig to which the subplots are added.
fontsize (float, optional) – Defaults to 10.

Returns

None. Instead fig is changed in place.

plot_solutions(fig, fontsize=10)[source]¶

Plot the solution trajectories.

For each trajectory (nr_pools+1) a new subplot is created.

Parameters

fig (Matplotlib figure) – The fig to which the subplots are added.
fontsize (float, optional) – Defaults to 10.

Returns

None. Instead fig is changed in place.

pool_age_densities_func(start_age_densities=None)[source]¶

Return a function that takes an array of ages and returns the pool age densities.

Parameters: start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\) for every pool. Defaults to None, meaning that all initial mass is considered to have zero age.
Returns: Python function p: p(ages) returns a numpy.ndarray len(ages) x len(times) x nr_pools containing the pool contents with the respective ages at the respective times, where ages is a numpy.array.

pool_age_densities_single_value(start_age_densities=None)[source]¶

Return a function for the pool age densities.

Parameters: start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\). Defaults to None, meaning that all initial mass is considered to have zero age.
Returns: Python function p_sv: p_sv(a, t) returns a numpy.array containing the pool contents with age a at time t.

pool_age_distribution_quantiles_pool_by_ode(quantile, pool, start_age_densities, F0=None, check_time_indices=None, **kwargs)[source]¶

Return pool age distribution quantile over the time grid for one single pool.

The compuation is done by solving an ODE as soon as the pool is nonempty. The initial value is obtained by computing the generalized inverse of the pool age distribution by a numerical root search algorithm.

Parameters

quantile (between 0 and 1) – The relative share of mass that is considered to be left of the computed value. A value of 0.5 leads to the computation of the median of the distribution.
pool (int) – The number of the pool for which the age quantile is to be computed.
start_age_densities (Python function) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\).
F0 (Python function, optional) – A function of age that returns a numpy.array containing the masses with age less than or equal to the age at time \(t_0\). It is fastest to provide F0, otherwise F0 will be computed by numerical integration of start_age_densities. Defaults to None.
check_time_indices (numpy.array, optional) –

Indices of the tiime
grid on which the ODE result are checked against an explicit solution computed by the pseudo-inverse of the cumulative

distribution function. Defaults to None in which case no check is performed.
kwargs – Passed to the solve_ivp, e.g., method or max_step.

Raises

Error – If start_age_densities is None.

Returns

(len(times)) The computed quantile values over the time grid.

Return type

numpy.ndarray

pool_age_distributions_quantiles(quantile, start_values=None, start_age_densities=None, F0=None, method='brentq', tol=1e-08)[source]¶

Return pool age distribution quantiles over the time grid.

The compuation is done by computing the generalized inverse of the respective cumulative distribution using a nonlinear root search algorithm. Depending on how slowly the cumulative distribution can be computed, this can take quite some time.

Parameters

quantile (between 0 and 1) – The relative share of mass that is considered to be left of the computed value. A value of 0.5 leads to the computation of the median of the distribution.
start_values (numpy.ndarray, len(times) x nr_pools, optional) – For each pool an array over the time grid of start values for the nonlinear search. Good values are slighty greater than the solution values. Defaults to an array of zeros for each pool
start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\). Defaults to None.
F0 (Python function) – A function of age that returns a numpy.array containing the masses with age less than or equal to the age at time \(t_0\). Defaults to None.
method (str) – The method that is used for finding the roots of a nonlinear function. Either ‘brentq’ or ‘newton’. Defaults to ‘brentq’.
tol (float) – The tolerance used in the numerical root search algorithm. A low tolerance decreases the computation speed tremendously, so a value of 1e-01 might already be fine. Defaults to 1e-08.

Raises

Error – If both start_age_densities and F0 are None. One must be given. It is fastest to provide F0, otherwise F0 will be computed by numerical integration of start_age_densities.

Returns

(len(times) x nr_pools) The computed quantile values over the time-pool grid.

Return type

numpy.ndarray

pool_age_distributions_quantiles_by_ode(quantile, start_age_densities, F0=None, check_time_indices=None, **kwargs)[source]¶

Return pool age distribution quantiles over the time grid.

The compuation is done by solving an ODE for each pool as soon as the pool is nonempty. The initial value is obtained by computing the generalized inverse of the pool age distribution by a numerical root search algorithm.

Parameters

quantile (between 0 and 1) – The relative share of mass that is considered to be left of the computed value. A value of 0.5 leads to the computation of the median of the distribution.
start_age_densities (Python function) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\).
F0 (Python function, optional) – A function of age that returns a numpy.array containing the masses with age less than or equal to the age at time \(t_0\). Defaults to None. It is fastest to provide F0, otherwise F0 will be computed by numerical integration of start_age_densities.
check_time_indices (numpy.array, optional) –

Indices of the tiime
grid on which the ODE result are checked against an explicit solution computed by the pseudo-inverse of the cumulative

distribution function. Defaults to None in which case no check is performed.
kwargs – Passed to the solve_ivp, e.g., method or max_step.

Returns

(len(times) x nr_pools) The computed quantile values over the time-pool grid.

Return type

numpy.ndarray

save_density_csv(filename, density, ages, times=None)[source]¶

Save density values over age-time grid to csv file.

Parameters

filename (str) – The name of the csv file to be written.
density (numpy.ndarray len(ages) x len(times)) – The density values to be saved over the age-time grid.
ages (numpy.array) – The ages corresponding to the indices in the zeroth dimension of the density ndarray.
times (numpy.array, optional) – An alternative time grid to be used. Defaults to None which means that the original time grid is going to be used.

Returns:

save_pools_and_system_density_csv(filename, pool_age_densities, system_age_density, ages)[source]¶

Save the pool and system age densities to a csv file.

The system value will be coded into pool number -1.

Parameters

filename (str) – The name of the csv file to be written.
pool_age_densites (numpy.ndarray len(ages) x len(times) x nr_pools) – The pool age density values.
system_age_density (numpy.ndarray len(ages) x len(times)) – The system age density values.
ages (numpy.array) – The ages that correspond to the indices in the zeroth dimension of the density arrays.

Returns

None

save_pools_and_system_value_csv(filename, pools_ndarr, system_arr)[source]¶

Save pool and system values to a csv file.

Values could be the mean age, for example. One dimension less than a density. The system value will be coded into pool number -1.

Parameters

filename (str) – The name of the csv file to be written.
pools_ndarr (numpy.ndarray len(times) x nr_pools) – The values to be saved over the time-pool grid.
system_arr (numpy.array) – The values over the time grid corresponding to the system.

Returns

Return type

None

save_state_transition_operator_cache(filename)[source]¶

save_value_csv(filename, arr, times=None)[source]¶

Save values over the time grid to a csv file.

Parameters

filename (str) – The name of the csv file to be written.
arr (numpy.array) – The values to be saved over the time grid.
times (np.array, optional) – The time grid to be used. Defaults to None in which case the orginal time grid is used.

Returns:

sol_funcs()[source]¶

Returns list of linearly interpolated solution functions per pool. :returns:

List of Python functions [f[i]], where f[i](t) returns
pool i’s content at time t.

sol_funcs_dict_by_name()[source]¶: Return linearly interpolated solution functions. as a dictionary indexed by the name (string) of the state variables

sol_funcs_dict_by_symbol()[source]¶: Return linearly interpolated solution functions as a dictionary indexed by the symbols of the state variables

sol_funcs_old()[source]¶

Return linearly interpolated solution functions.

Returns: Python function f: f(t) returns a numpy.array containing the pool contents at time t.

solve(alternative_start_values=None)[source]¶

Solve the model and return a solution grid. If the solution has been computed previously (even by other methods) the cached result will be returned.

Parameters: alternative_start_values (numpy.array) – If not given, the original start_values are used.
Returns: len(times) x nr_pools, contains the pool contents at the times given in the time grid.
Return type: numpy.ndarray

solve_func(alternative_start_values=None)[source]¶

Solve the model and return a function of time.

Parameters: alternative_start_values (numpy.array, optional) – If not given, the original start_values are used.
Returns: Python function f: f(t) is a numpy.array that containts the pool contents at time t.

solve_old(alternative_times=None, alternative_start_values=None)[source]¶

Solve the model and return a solution grid.

Parameters

alternative_times (numpy.array) – If not given, the original time grid is used.
alternative_start_values (numpy.array) – If not given, the original start_values are used.

Returns

len(times) x nr_pools, contains the pool contents at the times given in the time grid.

Return type

numpy.ndarray

solve_single_value_old(alternative_start_values=None)[source]¶

Solve the model and return a function of time.

Parameters: alternative_start_values (numpy.array, optional) – If not given, the original start_values are used.
Returns: Python function f: f(t) is a numpy.array that containts the pool contents at time t.

system_age_density(pool_age_densities)[source]¶

Return the system age density based on the given pool age densities.

Parameters: pool_age_densites (numpy.ndarray len(ages) x len(times) x nr_pools) – The pool age density values.
Returns: (len(ages) x len(times)) The sum of the pool age contents over all pools.
Return type: numpy.ndarray

system_age_density_single_value(start_age_densities=None)[source]¶

Return a function for the system age density.

Parameters: start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\). Defaults to None, meaning that all initial mass is considered to have zero age.
Returns: Python function sys_p_sv: sys_p_sv(a, t) returns the system content with age a at time t.

system_age_distribution_quantiles(quantile, start_values=None, start_age_densities=None, F0=None, method='brentq', tol=1e-08)[source]¶

Return system age distribution quantiles over the time grid.

The compuation is done by computing the generalized inverse of the respective cumulative distribution using a nonlinear root search algorithm. Depending on how slowly the cumulative distribution can be computed, this can take quite some time.

Parameters

quantile (between 0 and 1) – The relative share of mass that is considered to be left of the computed value. A value of 0.5 leads to the computation of the median of the distribution.
start_values (numpy.array, optional) – An array over the time grid of start values for the nonlinear search. Good values are slighty greater than the solution values. Must have the same length as times. Defaults to an array of zeros.
start_age_densities (Python function, optional) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\). Defaults to None.
F0 (Python function) – A function of age that returns a numpy.array containing the masses with age less than or equal to the age at time \(t_0\). Defaults to None.
method (str) – The method that is used for finding the roots of a nonlinear function. Either ‘brentq’ or ‘newton’. Defaults to ‘brentq’.
tol (float) – The tolerance used in the numerical root search algorithm. A low tolerance decreases the computation speed tremendously, so a value of 1e-01 might already be fide. Defaults to 1e-08.

Raises

Error – If both start_age_densities and F0 are None. One must be given. It is fastest to provide F0, otherwise F0 will be computed by numerical integration of start_age_densities.

Returns

The computed quantile values over the time grid.

Return type

numpy.array

system_age_distribution_quantiles_by_ode(quantile, start_age_densities, F0=None, check_time_indices=None, **kwargs)[source]¶

Return system age distribution quantile over the time grid.

The compuation is done by solving an ODE as soon as the system is nonempty. The initial value is obtained by computing the generalized inverse of the system age distribution by a numerical root search algorithm.

Parameters

quantile (between 0 and 1) – The relative share of mass that is considered to be left of the computed value. A value of 0.5 leads to the computation of the median of the distribution.
pool (int) – The number of the pool for which the age quantile is to be computed.
start_age_densities (Python function) – A function of age that returns a numpy.array containing the masses with the given age at time \(t_0\).
F0 (Python function, optional) – A function of age that returns a numpy.array containing the masses with age less than or equal to the age at time \(t_0\). It is fastest to provide F0, otherwise F0 will be computed by numerical integration of start_age_densities. Defaults to None.
check_time_indices (numpy.array, optional) –

Indices of the tiime
grid on which the ODE result are checked against an explicit solution computed by the pseudo-inverse of the cumulative

distribution function. Defaults to None in which case no check is performed.
kwargs – Passed to the solve_ivp, e.g., method or max_step.

Raises

Error – If start_age_densities is None.

Returns

The computed quantile values over the time grid.

Return type

numpy.ndarray

system_age_moment(order, start_age_moments=None)[source]¶

Compute the order th system age moment vector over the time grid by an ODE system.

The pool age moments are computed by age_moment_vector() and then weighted corresponding to the pool contents.

Parameters

order (int) – The order of the pool age moments to be computed.
start_age_moments (numpy.ndarray order x nr_pools, optional) – Given initial age moments up to the order of interest. Can possibly be computed by moments_from_densities(). Defaults to None assuming zero initial ages.

Returns

The order th system age moment over the time grid.

Return type

numpy.array

to_14C_explicit(start_values_14C, Fa_func, decay_rate=0.0001209681)[source]¶

Construct and return a SmoothModelRun instance that: models the 14C component additional to the original model run.

Parameters

start_values_14C (numpy.nd_array, nr_pools) – 14C start values.
Fa_func (func(t)) – returns atmospheric fraction to be multiplied with the input vector
rate (decay) – The decay rate to be used, defaults to 0.0001209681 (daily).

Returns

SmoothModelRun

x_solve_func_skew()[source]¶

CompartmentalSystems.smooth_model_run_14C module¶

exception CompartmentalSystems.smooth_model_run_14C.Error[source]¶

Bases: Exception

Generic error occurring in this module.

class CompartmentalSystems.smooth_model_run_14C.SmoothModelRun_14C(smr, start_values_14C, Fa_func, decay_rate=0.0001209681)[source]¶

Bases: CompartmentalSystems.smooth_model_run.SmoothModelRun

Construct and return a SmoothModelRun_14C instance that: models the 14C component of the original model run.

Parameters

smr (SmoothModelRun) – original model run
start_values_14C (numpy.nd_array, nr_pools) – 14C start values.
Fa_func (func(t)) – returns atmospheric fraction to be multiplied with the input vector
rate (decay) – The decay rate to be used, defaults to 0.0001209681 (daily).

acc_gross_external_input_vector_Delta_14C(data_times=None, alpha=None)[source]¶

acc_gross_external_output_vector(data_times=None)[source]¶

Return the vectors of accumulated external outputs.

Returns: len(times)-1 x nr_pools
Return type: numpy.ndarray

acc_gross_external_output_vector_Delta_14C(data_times=None, alpha=None)[source]¶

acc_gross_internal_flux_matrix_Delta_14C(data_times=None, alpha=None)[source]¶

acc_net_external_input_vector_Delta_14C(data_times=None, alpha=None)[source]¶

acc_net_external_output_vector(data_times=None)[source]¶

acc_net_external_output_vector_Delta_14C(data_times=None, alpha=None)[source]¶

acc_net_internal_flux_matrix_Delta_14C(data_times=None, alpha=None)[source]¶

external_output_flux_funcs()[source]¶

Return a dictionary of the external output fluxes.

Returns: {key: func} with key representing the pool from which the output comes and func a function of time that returns a float.
Return type: dict

property external_output_vector¶

Return the grid of external output vectors.

Returns: len(times) x nr_pools
Return type: numpy.ndarray

solve_Delta_14C(alpha=None)[source]¶

CompartmentalSystems.smooth_model_run_14C.pfile_C14Atm_NH()[source]¶

CompartmentalSystems.smooth_reservoir_model module¶

Module for symbolical treatment of smooth reservoir models.

This module handles the symbolic treatment of compartmental/reservoir/pool models. It does not deal with numerical computations and model simulations, but rather defines the underlying structure of the respective model.

All fluxes or matrix entries are supposed to be SymPy expressions. Smooth means that no Piecewise or DiracDelta functions should be involved in the model description.

Counting of compartment/pool/reservoir numbers starts at zero and the total number of pools is \(d\).

exception CompartmentalSystems.smooth_reservoir_model.Error[source]¶

Bases: Exception

Generic error occurring in this module.

class CompartmentalSystems.smooth_reservoir_model.SmoothReservoirModel(state_vector, time_symbol, input_fluxes={}, output_fluxes={}, internal_fluxes={})[source]¶

Bases: object

General class of smooth reservoir models.

state_vector¶

The model’s state vector \(x\). Its entries are SymPy symbols.

Type: SymPy dx1-matrix

state_variables¶

Names of the variables in the state vector. Its entries are of type str.

Type: list of str

time_symbol¶

The model’s time symbol.

Type: SymPy symbol

input_fluxes¶

The model’s external input fluxes. {key1: flux1, key2: flux2} with key the pool number and flux a SymPy expression for the influx.

Type: dict

output_fluxes¶

The model’s external output fluxes. {key1: flux1, key2: flux2} with key the pool number and flux a SymPy expression for the outflux.

Type: dict

internal_fluxes¶

The model’s internal_fluxes. {key1: flux1, key2: flux2} with key = (pool_from, pool_to) and flux a SymPy expression for the flux.

Type: dict

property F¶

The right hand side of the differential equation \(\dot{x}=B\,x+u\).

Type: SymPy dx1-matrix

age_moment_system(max_order)[source]¶

Return the age moment system of the model.

Parameters

max_order (int) – The maximum order up to which the age moment system is created (1 for the mean).

Returns

extended_state (SymPy d*(max_order+1)x1-matrix): The extended
state vector of the age moment system.
extended_rhs (SymPy d*(max_order+1)x1-matrix): The extended right
hand side of the age moment ODE.

Return type

tuple

all_fluxes()[source]¶

property compartmental_matrix¶

\(B\) from \(\dot{x} = B\,x+u\).

Returns: \(B = \xi\,T\,N\)
Return type: SymPy dxd-matrix
Type: SymPy Matrix

property external_inputs¶

Return the vector of external inputs.

Type: SymPy dx1-matrix

property external_outputs¶

Return the vector of external outputs.

Type: SymPy dx1-matrix

figure(figure_size=7, 7, logo=False, thumbnail=False)[source]¶

Return a figure representing the reservoir model.

Parameters

figure_size (2-tuple, optional) – Width and height of the figure. Defaults to (7,7).
logo (bool, optional) – If True, figure_size set to (3,3), no legend, smaller font size. Defaults to False.
thumbnail (bool, optional) – If True, produce a very small version, no legend. Defaults to False.

Returns

Figure representing the reservoir model.

Return type

Matplotlib figure

property free_symbols¶: Returns the superset of the free symbols of the flux expressions including the state variables.

classmethod from_B_u(state_vector, time_symbol, B, u) → CompartmentalSystems.smooth_reservoir_model.SmoothReservoirModel [source]¶

Construct and return a SmoothReservoirModel instance from: \(\dot{x}=B\,x+u\)

Parameters

state_vector (SymPy dx1-matrix) – The model’s state vector \(x\). Its entries are SymPy symbols.
time_symbol (SymPy symbol) – The model’s time symbol.
B (SymPy dxd-matrix) – The model’s compartmental matrix.
u (SymPy dx1-matrix) – The model’s external input vector.

Returns

SmoothReservoirModel

classmethod from_state_variable_indexed_fluxes(state_vector, time_symbol, input_fluxes, output_fluxes, internal_fluxes) → CompartmentalSystems.smooth_reservoir_model.SmoothReservoirModel [source]¶

Return an instance of SmoothReservoirModel.

Parameters

state_vector (SymPy dx1-matrix) – The model’s state vector \(x\). Its entries are SymPy symbols.
time_symbol (SymPy symbol) – The model’s time symbol.
input_fluxes (dict) – The model’s external input fluxes. {key1: flux1, key2: flux2} with key the symbol of the target pool. (as used in the state vector) and flux a SymPy expression for the influx.
output_fluxes (dict) – The model’s external output fluxes. {key1: flux1, key2: flux2} with key the symbols for the source pool (as used in the state vector) and flux a SymPy expression for the outflux.
internal_fluxes (dict) – The model’s internal_fluxes. {key1: flux1, key2: flux2} with key = (source pool symbol, target pool symbol) and flux a SymPy expression for the flux.

Returns

SmoothReservoirModel

property function_expressions¶: Returns the superset of the free symbols of the flux expressions.

property internal_inputs¶

Return the vector of internal inputs.

Type: SymPy dx1-matrix

property internal_outputs¶

Return the vector of internal outputs.

Type: SymPy dx1-matrix

property is_compartmental¶

Returns checks that all fluxes are nonnegative at the time of implementation this functionality sympy did not support relations in predicates yet. So while the following works:

with assuming(Q.positive(x) & Q.positive(y)):: print(ask(Q.positive(2*x+y)

it is not possible yet to get a meaningful answer to:

with assuming(Q.is_true(x>0) & Q.is_true(y>0)):: print(ask(Q.positive(2*x+y)

We therefore cannot implement more elaborate assumptions like k_1-(a_12+a_32)>=0 but still can assume all the state_variables and the time_symbol to be nonnegative Therefore we can check the compartmental_property best after all paramater value have been substituted. At the moment the function throws an exception if this is not the case.

property is_linear¶

Returns True if we can make SURE that the model is linear by checking that the jacobian is not state dependent. Note that external numerical functions of state variables are represented as sympy.Function f(x_1,x_2,…,t) Sympy will consider the derivative of math:df/dx_i with respect to state variable math:x_i as function math:g(x_1, x_2,…) too, since it can not exclude this possibility if we know f only numerically. In consequence this method will return False even if the numerical implementation of f IS linear in math:x_1,x_2,… . To avoid this situation you can just reformulate linear external functions math:f(x_1,x_2,…,t) as linear combinations of state independent external functions math:f(x_1,x_2,…,t)=g_1(t)x_1+g_2(t)x_2+… so that sympy can detect the linearity.

Returns: ‘True’, ‘False’
Return type: bool

is_state_dependent(expr)[source]¶

property jacobian¶

property no_input_model¶

property nr_pools¶

Return the number of pools involved in the model.

Type: int

plot_pools_and_fluxes(ax, mutation_scale=50, fontsize=24, thumbnail=False, legend=True)[source]¶

port_controlled_Hamiltonian_representation()[source]¶

tuple: \(J, R, N, x, u\) from \(\dot{x} = [J(x)-R(x)] \frac{\partial}{\partial x}H+u\). with \(H=\sum_i x_i \implies \frac{\partial}{\partial x}H =(1,1,...,1)\)

Returns

J (skew symmetric SymPy dxd-matrix) of internal fluxbalances: \(J_{i,j}=r_{j,i}-r_{i,j}\)
Q (SymPy dxd-matrix): Diagonal matrix describing the dissipation rates (outfluxes).
x (SymPy dx1-matrix): The model’s state vector.
u (SymPy dx1-matrix): The model’s external input vector.

Return type

tuple

property state_variable_set¶

steady_states(par_set=None)[source]¶

subs(parameter_dict)[source]¶: Returns a new instance of class: SmoothReservoirModel with all parameters in the parameter_dict replaced by their values by calling subs on all the flux expressions. :param parameter_dict: A dictionary with the structure {parameter_symbol:parameter_value,….}

xi_T_N_u_representation(factor_out_xi=True)[source]¶

tuple: \(\xi, T, N, x, u\) from \(\dot{x} = \xi\,T\,N\,x+u\).

Parameters

factor_out_xi (bool) – If true, xi is extracted from the matrix, otherwise \(xi=1\) will be returned. (Defaults to True.)

Returns

xi (SymPy number): Environmental coefficient.
T (SymPy dxd-matrix): Internal fluxes. Main diagonal contains
-1 entries.
N (SymPy dxd-matrix): Diagonal matrix containing the decomposition
rates.
x (SymPy dx1-matrix): The model’s state vector.
u (SymPy dx1-matrix): The model’s external input vector.

Return type

tuple

CompartmentalSystems.smooth_reservoir_model_14C module¶

Module for symbolical treatment of smooth 14C reservoir models.

This module handles the symbolic treatment of compartmental/reservoir/pool models. It does not deal with numerical computations and model simulations, but rather defines the underlying structure of the respective model.

All fluxes or matrix entries are supposed to be SymPy expressions. Smooth means that no Piecewise or DiracDelta functions should be involved in the model description.

Counting of compartment/pool/reservoir numbers starts at zero and the total number of pools is \(d\).

exception CompartmentalSystems.smooth_reservoir_model_14C.Error[source]¶

Bases: Exception

Generic error occurring in this module.

class CompartmentalSystems.smooth_reservoir_model_14C.SmoothReservoirModel_14C(srm, decay_symbol, Fa)[source]¶

Bases: CompartmentalSystems.smooth_reservoir_model.SmoothReservoirModel

General class of smooth 14C reservoir models.

decay_symbol¶

The model’s decay symbol.

Type: SymPy symbol

property output_fluxes_corrected_for_decay¶

CompartmentalSystems.start_distributions module¶

Module for computing age distributions or moments thereof to be used as start distributions in subsequent simulations.

The age distribution at the start \(t_0\) is NOT defined by the reservoir model and the initial values. In fact EVERY age distribution can be chosen. The implemented algorithms will correctly project it to any time \(t\). This module provides several ways to generate such a distribution.

The functions containing the word ‘distributions’ usually return a vector valued function of age representing pool wise age distributions that are NOT normalized (Integrating of a vector component over all ages yields the mass of the corresponding pool.) , and in some cases a start vector that should be used in the subsequent simulation for which the start age distribution is computed.

The functions containing the word ‘moments’ usually return an array: moments x pools, containing the moments of the pool ages . representing the initial values needed for the moment systems e.g. the mean of a start distribution to be used as initial value for the mean age system. The length of the list is determined by the maximum order of the moment system to be solved. In some cases a consistent start vector is also provided.

Zero start age distributions¶

The distributions eaisiest to imagine are those that start with zero age:

The one with all pools empty provided by: :py:meth:`start_age_distributions_from_zero_initial_content and the respective moments by: :py:meth:`start_age_moments_from_zero_initial_content
The one where all initial mass has age zero, provided by: :py:meth:`start_age_distributions_from_zero_age_initial_content

Established distributions¶

However for many applications one is interested in the CHANGE of an age distribution that has been established over a (possibly infinitely) long period of time.

Spinup

If we start the computation with all pools empty at time \(0\) and run it till time \(t = a_{max}\), the resulting distribution will be non zero only in the interval \([0,a_{max}]\). Such a distribution is provided by: start_age_distributions_from_empty_spinup() and the moments by: :py:meth:`start_age_moments_from_empty_spinup

Note that the finiteness of the spin up time has to be considered in the choice of questions that can be asked. For instance do not mistake the fact that the percentage of material older than \(a_{max}\) will increase over time for a property of the system, where it is actually a property of the start distribution resulting from the finiteness of the spin up time.
Distributions induced by steady states of the autonomous system, if those exist.

If the term ‘established’ is taken to the limit of infinity one can look for a related system that has persisted unchanged for all times \(t<t_0\) and start with the age distribution created by such a system. Such a distribution can only occur if the system has been in a steady state. For a general non-autonomous system this is very unlikely that such a steady state exist at all. However we can consider a related autonomous system resulting from ‘freezing’ the general non-autonomous system at a time \(t_0\). Even for such an autonomous system it is uncertain if and where equilibria exist. This has to be checked before an equilibrium age distribution can be computed. Actually the following steps have to be performed:
1. Transform the general nonlinear non-autonomous system into a nonlinear autonomous system by freezing it at time \(t=t_0\):
2. Compute \(u_0(x)=u(t_0,x_0)\) and \(B_0(x)=B(t_0,x_0)\)
3. Look for an equilibrium \(x_{fix}\) of the frozen system such that \(0=B_0(x_{fix})+u_0(x_{fix})\). If the frozen system is linear the we can compute the fixed point explicitly : \(x_{fix}=B_0^{-1}u_0\). In general the frozen system will be nonlinear and we will have to look for the fixed point numerically.
4. Compute the age distribution of the system at equilibrium \(x_{fix}\). This can be done using the formulas for linear autonomous pool models. (At the fixed point the nonlinear system is identical to a linear one. This is a special case of the general idea of linearizing a nonlinear model with respect to a trajectory, which in case of a fixed point is constant.)
All these steps are performed by start_age_distributions_from_steady_state() and start_age_moments_from_steady_state().

Note that \(x_{fix}\) is the compatible start vector that has to be used along with this start distributions for the following computation.

CompartmentalSystems.start_distributions.compute_fixedpoint_numerically(srm, t0, x0, parameter_dict, func_set)[source]¶

CompartmentalSystems.start_distributions.lapm_for_steady_state(srm, t0, parameter_dict, func_set, x0=None)[source]¶

If a fixedpoint of the frozen system can be found, create a linear autonomous model as an equivalent for the frozen (generally nonlinear) system there.

The function performs the following steps:

Substitute symbols and symbolic functions with the parameters and numeric functions.
Transform the general nonlinear non-autonomous system into a nonlinear autonomous system by freezing it at time \(t=t_0\):
Compute \(u_0(x)=u(t_0,x_0)\) and \(B_0(x)=B(t_0,x_0)\)
Look for an equilibrium \(x_{fix}\) of the frozen system such that \(0=B_0(x_{fix})+u_0(x_{fix})\). If the frozen system is linear the we can compute the fixed point explicitly : \(x_{fix}=B_0^{-1}u_0\). In general the frozen system will be nonlinear and we will have to look for the fixed point numerically.
Create a linear autonomous pool model. that can be investigated with the package LAPM

This is a special case of the general linearization of a nonlinear model along a trajectory, which in case of a fixed point is constant. At the fixed point the age distribution and solution of the nonlinear system are identical to those of a linear one.

Parameters

srm (SmoothReservoirModel) – The (symbolic) model
par_set (dict) – The parameter set that transforms the symbolic model into a numeric one. The keys are the sympy symbols, the values are the values used for the simulation.
func_set (dict) – The keys are the symbolic sympy expressions for external functions the values are the numeric functions to be used in the simulation.
t0 (float) – The time where the non-autonomous system is frozen.
x0 (numpy.ndarray) – An initial guess to start the fixed point iteration. If the frozen model is linear it will be ignored and can therefore be omitted.

Returns

lapm is an instance of (LAPM.linear_autonomous_pool_model.LinearAutonomousPoolModel) representing the linearization with respect to the (constant) fixed point trajectory. It yields the same age distribution as the frozen possibly nonlinear system at the fixed point. \(x_{fix}\) is a one dimensional vector representing the equilibrium. This is returned since it is very likely needed as start vector in the simulation for which the start distributions has been computed. (The computed distribution assumes the system to be in this state.)

Return type

(lapm, x_fix) (tuple)

CompartmentalSystems.start_distributions.start_age_distributions_from_empty_spinup(srm, t_max, parameter_dict, func_set)[source]¶

Finite age spin up from empty pools

Creates a SmoothModelRun object with empty pools at \(t=0\), runs it until \(t=t_{max}\) an returns the age distribution at \(t_{max}\)

Parameters

srm (SmoothReservoirModel) – The (symbolic) model
t_max (float) – The end of the spinup (which starts at t=0 and runs until t=t_max)
parameter_dict (dict) – The parameter set that transforms the symbolic model into a numeric one. The keys are the sympy symbols, the values are the values used for the simulation.
func_set (dict) – The keys are the symbolic sympy expressions for external functions the values are the numeric functions to be used in the simulation

Returns

a_dist_at_end_of_spinup is a vector valued function of age. a_dist_at_end_of_spinup(a)[i] reflects the mass of age \(a\) in pool i.

sol_t_max is a one dimensional vector representing the pool contents at the end of the spinup. This is returned since it is very likely needed as start vector in the simulation for which the start distributions has been computed.

Return type

(a_dist_at_end_of_spinup, sol_t_max) (tuple)

CompartmentalSystems.start_distributions.start_age_distributions_from_steady_state(srm, t0, parameter_dict, func_set, x0=None)[source]¶

Compute the age distribution of the system at equilibrium \(x_{fix}\) , by means of a linear autonomous pool model with identical age distributions. The fixed point and the linear model are provided by: lapm_for_steady_state()

Parameters

srm (SmoothReservoirModel) – The (symbolic) model
parameter_dict (dict) – The parameter set that transforms the symbolic model into a numeric one. The keys are the sympy symbols, the values are the values used for the simulation.
func_set (dict) – The keys are the symbolic sympy expressions for external functions, the values are the numeric functions to be used in the simulation
t0 (float) – The time where the non-autonomous system is frozen.
x0 (numpy.ndarray) – An initial guess to start the fixed point iteration. If the frozen model is linear it will be ignored and can therefore be omitted.

Returns

a_dist_function is a vector valued function of age. a_dist_function(a)[i] reflects the mass of age \(a\) in pool i. \(x_{fix}\) is a one dimensional vector representing the

equilibrium.

This is returned since it is very likely needed as start vector in the simulation for which the start distributions has been computed. (The computed distribution assumes the system to be in this state.)

Return type

(a_dist_function, x_fix) (tuple)

CompartmentalSystems.start_distributions.start_age_distributions_from_zero_age_initial_content(srm, x0)[source]¶

Returns the age distribution (function) for a system into which all initial mass is injected instantaneous at \(t=0\).

The function returns a vector valued function dist(a) of the age a that returns the start vector for a=0 and a zero vector of the same size everywhere else. This represents the age distribution of mass in a compartmental system where the age of all pool contents is zero. This means that the distribution is NOT normalized to 1.

Parameters

srm (SmoothReservoirModel) – The (symbolic) model
x0 (numpy.ndarray) – The contents of the pools at \(t=0\)

Returns

a vector valued function of age. dist(a)[i] reflects the mass of age \(a\) in pool i.

Return type

dist (callable)

CompartmentalSystems.start_distributions.start_age_distributions_from_zero_initial_content(srm)[source]¶

Returns the age distribution (function) for an empty system.

The function returns a vector valued function dist(a) of the age a that is zero everywhere dist(a)=0 for all a. This represents the age distribution of a compartmental system with all pools empty.

Parameters: srm (SmoothReservoirModel) – The (symbolic) model
Returns: dist, a vector valued function of age. dist(a)[i] reflects the mass of age \(a\) in pool i.
Return type: callable

CompartmentalSystems.start_distributions.start_age_moments_from_empty_spinup(srm, t_max, parameter_dict, func_set, max_order)[source]¶

CompartmentalSystems.start_distributions.start_age_moments_from_steady_state(srm, t0, parameter_dict, func_set, max_order)[source]¶

Compute the age moments of the system at equilibrium \(x_{fix}\) , by means of a linear autonomous pool model with identical age distributions. The fixed point and the linear model are provided by: lapm_for_steady_state()

Parameters

srm (SmoothReservoirModel) – The (symbolic) model
parameter_dict (dict) – The parameter set that transforms the symbolic model into a numeric one. The keys are the sympy symbols, the values are the values used for the simulation.
func_set (dict) – The keys are the symbolic sympy expressions for external functions the values are the numeric functions to be used in the simulation
t0 (float) – The time where the non-autonomous system is frozen.
x0 (numpy.ndarray) – An initial guess to start the fixed point iteration. If the frozen model is linear it will be ignored and can therefore be omitted.
max_order (int) – The highest order up to which moments are to be computed.

Returns

moments x pools, containing the moments of the pool ages in equilibrium.

Return type

numpy.ndarray

CompartmentalSystems.start_distributions.start_age_moments_from_zero_initial_content(srm, max_order)[source]¶

The function returns an array of shape (max_order, srm.nr_pools)

The values are set to numpy.nan to be consistent with other parts of the code. For instance the mean age (first moment) of an empty pool is undefined ) :param srm: The (symbolic) model :type srm: SmoothReservoirModel :param max_order: The highest order up to which moments are to be

computed.

CompartmentalSystems package¶

Subpackages¶

Submodules¶

CompartmentalSystems.BlockIvp module¶

CompartmentalSystems.BlockOde module¶

CompartmentalSystems.Cache module¶

CompartmentalSystems.cs_plotter module¶

CompartmentalSystems.discrete_model_run module¶

CompartmentalSystems.discrete_model_run_14C module¶

CompartmentalSystems.discrete_model_run_with_gross_fluxes module¶

CompartmentalSystems.example_smooth_model_runs module¶

CompartmentalSystems.example_smooth_reservoir_models module¶

CompartmentalSystems.helpers_reservoir module¶

CompartmentalSystems.model_run module¶

CompartmentalSystems.myOdeResult module¶

CompartmentalSystems.picklegzip module¶

CompartmentalSystems.pwc_model_run module¶

CompartmentalSystems.pwc_model_run_14C module¶

CompartmentalSystems.pwc_model_run_fd module¶

CompartmentalSystems.smooth_model_run module¶

CompartmentalSystems.smooth_model_run_14C module¶

CompartmentalSystems.smooth_reservoir_model module¶

CompartmentalSystems.smooth_reservoir_model_14C module¶

CompartmentalSystems.start_distributions module¶

Zero start age distributions¶

Established distributions¶

Module contents¶

CompartmentalSystems

Navigation

Related Topics