Evaluation Functions

eval_ode_run

Evaluation functions for the EnsembleProblem. Given a set of measures the solution sol is evaluated seperatly per dimension. An additional set of global measures take in the complete solution and return a single number or a matrix. Handing over this function to DEMCBBProblem (and thus also to EnsembleProblem) the expected signature is (sol, i::Int) -> (results, repeat::Bool). Here, there are several more general versions that can be adjusted to the experiment.

eval_ode_run(sol, i, state_filter::Array{Int64,1}, eval_funcs::Array{<:Function,1}, matrix_eval_funcs::Union{AbstractArray, Nothing}=nothing, global_eval_funcs::Union{AbstractArray, Nothing}=nothing; failure_handling::Symbol=:None, cyclic_setback::Bool=false, replace_inf=nothing, flag_past_measures=false, state_filter_only_per_dim=false)

• sol: solution of one of the EnsembleProblem runs, should have only timesteps with constant time intervals between them
• i: Int, number of iteration/run
• state_filter: Array with indicies of all dimensions (of the solutions) that should be evaluated
• eval_funcs: Array of functions that should be applied to every dimension of the solution. Signature: (sol::AbstractArray) -> measure or (sol::AbstractArray, previous_results::AbstractArray) -> measure depending on the value of flag_past_measures.
• matrix_eval_funcs: Array of functions that should be applied to the complete Ndimensional solution and return a matrix (or vector), like e.g covariance or correlation, signature (N-Dim Array w/ Samples ::AbstractArray) -> Measure::AbstractArray (technical detail: length(measure)!=Ndim (system dimension))
• global_eval_funcs: Array of functions that should be applied to the complete N-dimensional solution, signature (N-Dim Array w/ Samples ::AbstractArray) -> Measure::Number
• failure_handling: How failure of integration is handled. Should be :None (do no checks), :Inf (If retcode==:DtLessThanMin: return Inf) or :Repeat (If no succes, repeat the trial (only works with random initial conditions))
• cyclic_setback: Bool, if true $N*2\pi$ is substracted from the solution so that the first element of the solution that is analysed is within $[-\pi, \pi]$. Usefull e.g. for phase oscillators.
• replace_inf: Number or Nothing, if a number replaces all Infs in the solution with this number. Can be usefull if one still wants to distinguish between different solutions containing Infs, +Inf is replaced by the Number, -Inf by (-1)*Number.
• flag_past_measures::Bool: If true als function within eval_funcs also receive the previous results (of the other measures for the same dimension) as an extra arguments. Thus all functions need to have a signature (sol::AbstractArray, previous_results::AbstractArray) -> measure. If false the functions only receive the solution vector, thus the function should have the signature (sol::AbstractArray) -> measure
• state_filter_only_per_dim: Only per-dimension measures are affected by the state_filter.

Example function

In order to derive a valid evaluation function from this for the MCBBProblem one can define a function similar to this:

function my_eval_ode_run(sol, i)
    N_dim = length(sol.prob.u0)
    state_filter = collect(1:N_dim)
    eval_funcs = [mean, std]
    eval_ode_run(sol, i, state_filter, eval_funcs)
end

Utilizing previous results

If one wants to utilze the previous results (and don't compute measures twice), one has to use the flag_past_measures=true option. This is only possible for the per dimension measures. An example could read:

function my_eval_ode_run(sol, i)
    N_dim = length(sol.prob.u0)
    state_filter = collect(1:N_dim)
    meanval(u::AbstractArray, past_measures::AbstractArray) = StatsBase.mean(u)
    standarddev(u::AbstractArray, past_measures::AbstractArray) = StatsBase.std(u; mean=past_measures[1], corrected=false)
    eval_funcs = [meanval, standarddev, empirical_1D_KL_divergence_hist]
    eval_ode_run(sol, i, state_filter, eval_funcs; flag_past_measures=true)
end

Latter function is also already available as a default eval_ode_run in this library. The order of the functions is important. In this example meanval will always get an empty array as the second argument, standarddev will get an array with the result from meanval as the only value and empirical_1D_KL_divergence_hist will get an additional array with the results from meanval and standarddev.

eval_ode_run(sol, i)

Default eval_ode_run, identical to the code above.

Continue Integration / Response Analysis

 eval_ode_run(sol, i, state_filter::Array{Int64,1}, eval_funcs::Array{<:Function,1},  matrix_eval_funcs::Union{AbstractArray, Nothing}, global_eval_funcs::Union{AbstractArray, Nothing}, par_var::OneDimParameterVar, eps::Float64, par_bounds::AbstractArray, distance_matrix_func; failure_handling::Symbol=:None, cyclic_setback::Bool=false, flag_past_measures::Bool=false, N_t::Int=200, alg=nothing, debug::Bool=false, return_pm::Bool, new_tspan::Union{Nothing, AbstractArray}, kwargs...)

Evaluation function that continues each integration and computes the same measures for par+eps and par-eps. Returns the results of the usual eval_ode_run (all measures) and additionally the response of the distance function to the paramater increase/decrease.

• par_var: ParameterVar struct, same as handed over to the problem type.

• eps: Number, response at par+/-eps

• distance_matrix_func: Same distance matrix functions that will also be used for the later analysis/clustering, expected signature: (sol::MCBBSol, prob::MCBBProblem) -> D::AbstractArray, . Attension: if the weight vector is provided this version of the distance it needs to have one less element as the function later used before clustering because the result of the response analysis is an additional measure.

• N_t: Time steps for the continued integration

• alg: Algorithm for solve()

• debug: If true, also returns the DifferentialEquations problem solved for the continuation.

• return_pm: If true, returns D(p+dp) AND D(p-dp). If false returns the mean of these.

• new_tspan: timespan for continued integeration, default: 15% of the original timespan

• all further keyword arguments will be handed over to solve(prob::DEMCBBProblem, ...)

empirical_1D_KL_divergence_hist(u::AbstractArray, mu::Number, sig::Number, hist_bins::Int=31, n_stds::Number=3, sig_tol=1e-4::Number)

One measure that can be used with eval_ode_run. Computes the empirical Kullback-Leibler Divergence of the input to a normal distribution with the same mean and std as the input, thus it is a measure of normality. This version does this with histograms.

• u: Input Array
• mu: Mean of u
• sig: Std of u
• hist_bins: number of bins of the histogram to estimate the empirical pdf of the data
• n_stds: Interval that the histogram covers in numbers of stds (it covers mean +/- n_stds*std)
• sig_tol: At times, all KL div methods run into difficulties when sig gets really small, for sig<sig_tol 0 is returned as a result because in the limit of sig -> 0 the reference distribution is a delta distribution and the data is constant thus also a delta distribution. hence the distributions are identical and the KL div should be zero.

empirical_1D_KL_divergence_pc(u::AbstractArray, mu::Number, sig::Number)

One measure that can be used with eval_ode_run. Computes the empirical Kullback-Leibler Divergence of the input to a normal distribution with the same mean and std as the input, thus it is a measure of normality. This version does this based on a linearly interpolated emperical CDF, see Perez-Cruz (IEEE, 2008). This version can run into numerical difficulties for discrete systems with alternating inputs like [a,-a,a,a] and large 2a. For reasonable continous input it is a better and parameter free approximation to the KL divergence than the histogram estimator.

• u: Input Array
• mu: Mean of u
• sig: Std of u

wasserstein_ecdf(u::AbstractArray, mu::Number, sig::Number)

One measure that can be used with eval_ode_run. Computes the 1-wasserstein distance based on ECDFs.

• u: Input Array
• mu: Mean of u
• sig: Std of u

correlation(sol::AbstractArray)

Example function for matrix_eval_funcs. This routine calculates the absolute value of the Pearson correlation between the time series of all system dimensions and returns it as a matrix.

correlation_hist(sol::AbstractArray, nbins::Int=30)

Example function for matrix_eval_funcs. This routine calculates the absolute value of the Pearson correlation between the time series of all system dimensions and returns the weights of histogram fitted to all of these values. It uses the same binning for all calculations with the edges calculated by 0:1. /nbins:1.

correlation_ecdf(sol::AbstractArray, nbins::Int=30)

Example function for matrix_eval_funcs. This routine calculates the absolute value of the Pearson correlation between the time series of all system dimensions and the ECDF of a histogram fitted to all of these values. It uses the same binning for all calculations with the edges calculated by 0:1. /nbins:1.

check_inf_nan(sol::myMCSol)

Checks if any of the results is inf or nan and returns the indices in a dictionary with keys Inf and NaN