text:
using Test
using JuMP
using MathOptInterface
using MIPVerify
using MIPVerify:
    relu,
    get_target_indexes,
    set_max_indexes,
    get_max_index,
    matmul,
    tight_upperbound,
    tight_lowerbound,
    abs_ge,
    masked_relu,
    is_constant,
    get_tightening_algorithm,
    mip,
    lp,
    interval_arithmetic,
    DEFAULT_TIGHTENING_ALGORITHM,
    TighteningAlgorithm,
    MIPVerifyExt,
    upper_bound,
    lower_bound
@isdefined(TestHelpers) || include("../TestHelpers.jl")

function count_binary_variables(m::Model)
    JuMP.num_constraints(m, JuMP.VariableRef, MathOptInterface.ZeroOne)
end

TestHelpers.@timed_testset "core_ops.jl" begin
    @testset "is_constant" begin
        @testset "JuMP.AffExpr" begin
            m = TestHelpers.get_new_model()
            @test is_constant(zero(JuMP.VariableRef))
            @test is_constant(one(JuMP.VariableRef))
            x = @variable(m)
            y = @variable(m)
            z = 2 * x + 3 * y
            @test !is_constant(z)
        end
        @testset "JuMP.VariableRef" begin
            m = TestHelpers.get_new_model()
            x = @variable(m)
            @test !is_constant(x)
        end
    end

    @testset "owner_model" begin
        m = TestHelpers.get_new_model()
        y1 = @variable(m)
        x1 = one(JuMP.VariableRef) * 1
        x2 = one(JuMP.VariableRef) * 2
        @test MIPVerify.owner_model([x1, y1]) == m
        @test_throws MethodError MIPVerify.owner_model([x1, x2])
    end

    @testset "get_tightening_algorithm" begin
        m = TestHelpers.get_new_model()

        tightening_algorithms = [interval_arithmetic, mip, lp]

        @testset "if variable known to be constant, always use interval_arithmetic" begin
            x = one(JuMP.VariableRef) # is_constant(x)==true
            for alg in tightening_algorithms
                @test get_tightening_algorithm(x, alg) == interval_arithmetic
            end
            @test get_tightening_algorithm(x, nothing) == interval_arithmetic
        end

        @testset "if variable not known to be constant" begin
            @testset "use tightening algorithm if specified" begin
                m = TestHelpers.get_new_model()
                y = @variable(m)
                for alg in tightening_algorithms
                    @test get_tightening_algorithm(y, alg) == alg
                end
            end
            @testset "fall back to model-level tightening algorithm if specified" begin
                m = TestHelpers.get_new_model()
                y = @variable(m)
                for alg in tightening_algorithms
                    m.ext[:MIPVerify] = MIPVerifyExt(alg)
                    @test get_tightening_algorithm(y, nothing) == alg
                end
            end
            @testset "fall back to package default tightening algorithm as last resort" begin
                m = TestHelpers.get_new_model()
                y = @variable(m)
                @test get_tightening_algorithm(y, nothing) == DEFAULT_TIGHTENING_ALGORITHM
            end
        end
    end

    @testset "maximum(xs)" begin
        @testset "single variable to maximize over" begin
            m = TestHelpers.get_new_model()
            x1 = @variable(m, lower_bound = 0, upper_bound = 3)
            xmax = MIPVerify.maximum([x1])

            # no binary variables need to be introduced
            @test count_binary_variables(m) == 0

            @objective(m, Max, x1)
            optimize!(m)
            solve_output = JuMP.value(xmax)
            @test solve_output ≈ 3
        end
        @testset "multiple variables to maximize over, some constant" begin
            m = TestHelpers.get_new_model()
            x0 = one(JuMP.VariableRef) # add constant variable at start
            x1 = @variable(m, lower_bound = 0, upper_bound = 3)
            x2 = @variable(m, lower_bound = 4, upper_bound = 5)
            x3 = @variable(m, lower_bound = 2, upper_bound = 7)
            x4 = @variable(m, lower_bound = -1, upper_bound = 1)
            x5 = @variable(m, lower_bound = -3, upper_bound = 1)
            xmax = MIPVerify.maximum([x0, x1, x2, x3, x4, x5])

            # an efficient implementation does not add binary variables for x1, x4 and x5
            @test count_binary_variables(m) <= 2

            # elements of the input array are made to take their maximum value
            @objective(m, Max, x1 + x2 + x3 + x4 + x5)
            optimize!(m)

            solve_output = JuMP.value(xmax)
            @test solve_output ≈ 7
        end
        @testset "single variable to maximize over, constant" begin
            m = TestHelpers.get_new_model()
            x1 = one(JuMP.VariableRef) * 3
            xmax = MIPVerify.maximum([x1])

            # no binary variables need to be introduced
            @test count_binary_variables(m) == 0

            solve_output = JuMP.value(xmax)
            @test solve_output ≈ 3
        end
        @testset "multiple variables to maximize over, all constant" begin
            m = TestHelpers.get_new_model()
            x1 = one(JuMP.VariableRef)
            x2 = one(JuMP.VariableRef) * 2
            xmax = MIPVerify.maximum([x1, x2])

            # no binary variables need to be introduced
            @test count_binary_variables(m) == 0

            solve_output = JuMP.value(xmax)
            @test solve_output ≈ 2
        end
        @testset "regression test to deal with indexing issue in v0.8.0" begin
            m = TestHelpers.get_new_model()
            x1 = @variable(m, lower_bound = -2, upper_bound = 2) # upper_bound of this variable is low enough that it gets filtered away
            x2 = @variable(m, lower_bound = 2.5, upper_bound = 100)
            x3 = @variable(m, lower_bound = 3, upper_bound = 3.3)
            xmax = MIPVerify.maximum([x1, x2, x3])

            # an efficient implementation does not add binary variables for x1
            @test count_binary_variables(m) <= 2

            # elements of the input array are made to take their maximum value
            @objective(m, Max, xmax)
            optimize!(m)

            solve_output = JuMP.value(xmax)
            @test solve_output ≈ 100
        end
        @testset "lower_bound on one matches upper_bound on another; output expected to be constant" begin
            m = TestHelpers.get_new_model()
            x1 = @variable(m, lower_bound = -6, upper_bound = 2)
            x2 = @variable(m, lower_bound = 2, upper_bound = 2)

            xmax = MIPVerify.maximum([x1, x2])

            # no binary variables need to be introduced
            @test count_binary_variables(m) == 0

            optimize!(m)

            solve_output = JuMP.value(xmax)
            @test solve_output ≈ 2
        end
    end

    @testset "maximum(xs, ls, us)" begin
        @testset "single variable to maximize over" begin
            m = TestHelpers.get_new_model()
            x1 = @variable(m, lower_bound = 0, upper_bound = 3)
            xmax = MIPVerify.maximum([x1], [0], [1])

            # no binary variables need to be introduced
            @test count_binary_variables(m) == 0

            @objective(m, Max, x1)
            optimize!(m)
            solve_output = JuMP.value(xmax)
            @test solve_output ≈ 3
        end
        @testset "multiple variables to maximize over, all constant" begin
            m = TestHelpers.get_new_model()
            x1 = one(JuMP.VariableRef)
            x2 = one(JuMP.VariableRef) * 2
            xmax = MIPVerify.maximum([x1, x2], [1, 2], [1, 2])

            # no binary variables need to be introduced
            @test count_binary_variables(m) == 0

            solve_output = JuMP.value(xmax)
            @test solve_output ≈ 2
        end
    end

    @testset "maximum_ge" begin
        @testset "multiple variables to maximize over, all constant" begin
            m = TestHelpers.get_new_model()
            x1 = one(JuMP.VariableRef)
            x2 = one(JuMP.VariableRef) * 2
            xmax = MIPVerify.maximum_ge([x1, x2])

            solve_output = JuMP.value(xmax)
            @test solve_output ≈ 2
        end
    end

    @testset "abs_ge" begin
        @testset "positive constant input" begin
            m = TestHelpers.get_new_model()
            x = one(JuMP.VariableRef)
            x_abs = MIPVerify.abs_ge(x)

            solve_output = JuMP.value(x_abs)
            @test solve_output ≈ 1
        end

        @testset "negative constant input" begin
            m = TestHelpers.get_new_model()
            x = one(JuMP.VariableRef) * -1
            x_abs = MIPVerify.abs_ge(x)

            solve_output = JuMP.value(x_abs)
            @test solve_output ≈ 1
        end
    end

    @testset "relu(x)" begin
        @testset "Numerical Input" begin
            @test relu(5) == 5
            @test relu(0) == 0
            @test relu(-1) == 0
        end

        @testset "Variable Input" begin
            @testset "constant" begin
                x = one(JuMP.VariableRef)
                x_r = relu(x)
                @test MIPVerify.is_constant(x_r)
                @test x_r.constant == x.constant
            end
            @testset "strictly non-negative" begin
                m = TestHelpers.get_new_model()
                x = @variable(m, lower_bound = 0, upper_bound = 1)
                x_r = relu(x)

                # no binary variables should be introduced
                @test count_binary_variables(m) == 0

                @objective(m, Max, 2 * x_r - x)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 1
            end
            @testset "strictly non-positive" begin
                m = TestHelpers.get_new_model()
                x = @variable(m, lower_bound = -1, upper_bound = 0)
                x_r = relu(x)

                # no binary variables should be introduced
                @test count_binary_variables(m) == 0

                @objective(m, Max, 2 * x_r - x)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 1
            end
            @testset "regular" begin
                m = TestHelpers.get_new_model()
                x = @variable(m, lower_bound = -1, upper_bound = 2)
                x_r = relu(x)

                # at most one binary variable to be introduced
                @test count_binary_variables(m) <= 1

                @objective(m, Max, 2 * x_r - x)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 2
            end
        end
    end

    @testset "relu(x, l, u)" begin
        @testset "Variable Input" begin
            @testset "strictly non-negative" begin
                m = TestHelpers.get_new_model()
                x = @variable(m)
                x_r = relu(x, 1, 2)
                @test count_binary_variables(m) == 0

                @objective(m, Max, x_r - x)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 0
            end
            @testset "strictly non-positive" begin
                m = TestHelpers.get_new_model()
                x = @variable(m)
                x_r = relu(x, -2, -1)
                @test upper_bound(x_r) == 0
                @test lower_bound(x_r) == 0
            end
            @testset "constant" begin
                m = TestHelpers.get_new_model()
                x = @variable(m)
                x_r = relu(x, 2, 2)
                @test upper_bound(x_r) == 2
                @test lower_bound(x_r) == 2
            end
        end
    end

    @testset "masked_relu" begin
        @testset "Numerical Input" begin
            @test masked_relu(5, -1) == 0
            @test masked_relu(0, -1) == 0
            @test masked_relu(-5, -1) == 0
            @test masked_relu(5, 0) == 5
            @test masked_relu(0, 0) == 0
            @test masked_relu(-5, 0) == 0
            @test masked_relu(5, 1) == 5
            @test masked_relu(0, 1) == 0
            @test masked_relu(-5, 1) == -5
        end

        @testset "Variable Input, single" begin
            @testset "mask is negative" begin
                m = TestHelpers.get_new_model()
                x = @variable(m, lower_bound = -1, upper_bound = 2)
                x_r = masked_relu(x, -1)

                # no binary variables to be introduced
                @test count_binary_variables(m) == 0

                @objective(m, Max, 2 * x_r - x)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 1
                @test JuMP.value(x) ≈ -1
                @test JuMP.value(x_r) ≈ 0

                @objective(m, Min, 2 * x_r - x)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ -2
                @test JuMP.value(x) ≈ 2
                @test JuMP.value(x_r) ≈ 0
            end
            @testset "mask is 0" begin
                m = TestHelpers.get_new_model()
                x = @variable(m, lower_bound = -1, upper_bound = 2)
                x_r = masked_relu(x, 0)

                # at most one binary variable to be introduced
                @test count_binary_variables(m) <= 1

                @objective(m, Max, 2 * x_r - x)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 2
                @test JuMP.value(x) ≈ 2
                @test JuMP.value(x_r) ≈ 2

                @objective(m, Min, 2 * x_r - x)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 0
                @test JuMP.value(x) ≈ 0
                @test JuMP.value(x_r) ≈ 0
            end
            @testset "mask is positive" begin
                m = TestHelpers.get_new_model()
                x = @variable(m, lower_bound = -1, upper_bound = 2)
                x_r = masked_relu(x, 1)

                # no binary variables to be introduced
                @test count_binary_variables(m) == 0

                @objective(m, Max, 2 * x_r - x)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 2
                @test JuMP.value(x) ≈ 2
                @test JuMP.value(x_r) ≈ 2

                @objective(m, Min, 2 * x_r - x)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ -1
                @test JuMP.value(x) ≈ -1
                @test JuMP.value(x_r) ≈ -1
            end
        end

        @testset "Variable Input, array" begin
            @testset "invalid mask" begin
                m = TestHelpers.get_new_model()
                @variable(m, x[1:4], lower_bound = -1, upper_bound = 2)

                @test_throws AssertionError masked_relu(x, [-1, 0, 1])
            end
            @testset "valid mask" begin
                m = TestHelpers.get_new_model()
                @variable(m, x[1:3], lower_bound = -1, upper_bound = 2)

                x_r = masked_relu(x, [-1, 0, 1])

                @objective(m, Max, sum(2 * x_r - x))
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 5
                @test JuMP.value.(x) ≈ [-1, 2, 2]
                @test JuMP.value.(x_r) ≈ [0, 2, 2]

                @objective(m, Min, sum(2 * x_r - x))
                optimize!(m)
                @test JuMP.objective_value(m) ≈ -3
                @test JuMP.value.(x) ≈ [2, 0, -1]
                @test JuMP.value.(x_r) ≈ [0, 0, -1]
            end
        end
    end

    @testset "abs_ge" begin
        @testset "strictly non-negative" begin
            m = TestHelpers.get_new_model()
            x = @variable(m, lower_bound = 0, upper_bound = 1)
            x_a = abs_ge(x)

            # no binary variables should be introduced
            @test count_binary_variables(m) == 0

            @objective(m, Max, 2 * x_a - x)
            optimize!(m)
            @test JuMP.objective_value(m) ≈ 1
        end
        @testset "strictly non-positive" begin
            m = TestHelpers.get_new_model()
            x = @variable(m, lower_bound = -1, upper_bound = 0)
            x_a = abs_ge(x)

            # no binary variables should be introduced
            @test count_binary_variables(m) == 0

            @objective(m, Max, 2 * x_a - x)
            optimize!(m)
            @test JuMP.objective_value(m) ≈ 3
        end
        @testset "regular" begin
            m = TestHelpers.get_new_model()
            x = @variable(m, lower_bound = -2, upper_bound = 2)
            x_a = abs_ge(x)

            # no binary variables should be introduced
            @test count_binary_variables(m) == 0

            @objective(m, Max, 2 * x_a - x)
            optimize!(m)
            @test JuMP.objective_value(m) ≈ 6
        end
        @testset "abs_ge is not strict" begin
            # in particular we only need to satisfy the property |x_a| > x
            m = TestHelpers.get_new_model()
            x = @variable(m, lower_bound = -4, upper_bound = 2)
            x_a = abs_ge(x)

            # no binary variables should be introduced
            @test count_binary_variables(m) == 0

            @objective(m, Min, x_a - x)
            optimize!(m)
            @test JuMP.objective_value(m) ≈ 0
        end
    end

    @testset "get_target_indexes" begin
        @test get_target_indexes(1, 5) == [1]
        @test get_target_indexes(1, 5, invert_target_selection = true) == [2, 3, 4, 5]
        @test get_target_indexes([2, 4], 5) == [2, 4]
        @test get_target_indexes([2, 4], 5, invert_target_selection = true) == [1, 3, 5]
        @test_throws AssertionError get_target_indexes(6, 5)
        @test_throws AssertionError get_target_indexes([1, 6], 5)
    end

    @testset "set_max_indexes" begin
        @testset "single target index" begin
            @testset "vanilla" begin
                m = TestHelpers.get_new_model()
                x = @variable(m, [i = 1:3])
                @constraint(m, x[2] == 5)
                @constraint(m, x[3] == 1)
                set_max_indexes(m, x, [1])
                @objective(m, Min, x[1])
                optimize!(m)
                @test JuMP.value(x[1]) ≈ 5
            end
            @testset "with margin" begin
                margin = 3
                m = TestHelpers.get_new_model()
                x = @variable(m, [i = 1:3])
                @constraint(m, x[2] == 5)
                @constraint(m, x[3] == 1)
                set_max_indexes(m, x, [1], margin = margin)
                @objective(m, Min, x[1])
                optimize!(m)
                @test JuMP.value(x[1]) ≈ 5 + margin
            end
        end
        @testset "multiple target indexes" begin
            @testset "vanilla" begin
                m = TestHelpers.get_new_model()
                x = @variable(m, [i = 1:3])
                @constraint(m, x[1] == 5)
                set_lower_bound(x[2], 0)
                set_upper_bound(x[2], 10)
                set_lower_bound(x[3], -1)
                set_upper_bound(x[3], 10)
                set_max_indexes(m, x, [2, 3])
                @objective(m, Min, x[2] + x[3])
                optimize!(m)
                @test JuMP.value(x[2]) ≈ 5
                @test JuMP.value(x[3]) ≈ -1
            end
            @testset "with margin" begin
                margin = 3
                m = TestHelpers.get_new_model()
                x = @variable(m, [i = 1:3])
                @constraint(m, x[1] == 5)
                set_lower_bound(x[2], 0)
                set_upper_bound(x[2], 10)
                set_lower_bound(x[3], -1)
                set_upper_bound(x[3], 10)
                set_max_indexes(m, x, [2, 3], margin = margin)
                @objective(m, Min, x[2] + x[3])
                optimize!(m)
                @test JuMP.value(x[2]) ≈ 5 + margin
                @test JuMP.value(x[3]) ≈ -1
            end
            @testset "first JuMPLinearType is constant" begin
                @testset "selected variable has non-constant value, and can take the maximum value" begin
                    m = TestHelpers.get_new_model()
                    x1 = one(JuMP.VariableRef)
                    x2 = @variable(m, lower_bound = 4, upper_bound = 5)
                    set_max_indexes(m, [x1, x2], [2])
                    @objective(m, Min, x2)
                    optimize!(m)
                    @test JuMP.value(x2) ≈ 4
                end
                @testset "selected variable has non-constant value, and cannot take the maximum value" begin
                    m = TestHelpers.get_new_model()
                    x1 = one(JuMP.VariableRef)
                    x2 = @variable(m, lower_bound = -5, upper_bound = -4)
                    set_max_indexes(m, [x1, x2], [2])
                    @objective(m, Min, x2)
                    optimize!(m)
                    solve_status = JuMP.termination_status(m)
                    @test (
                        solve_status in
                        [MathOptInterface.INFEASIBLE, MathOptInterface.INFEASIBLE_OR_UNBOUNDED]
                    )
                end
                @testset "selected variable has constant value, and can take the maximum value" begin
                    m = TestHelpers.get_new_model()
                    x1 = one(JuMP.VariableRef)
                    x2 = @variable(m, lower_bound = -5, upper_bound = -4)
                    set_max_indexes(m, [x1, x2], [1])
                    @objective(m, Min, x2)
                    optimize!(m)
                    @test JuMP.value(x2) ≈ -5
                end
                @testset "selected variable has constant value, and cannot take the maximum value" begin
                    m = TestHelpers.get_new_model()
                    x1 = one(JuMP.VariableRef)
                    x2 = @variable(m, lower_bound = 4, upper_bound = 5)
                    set_max_indexes(m, [x1, x2], [1])
                    @objective(m, Min, x2)
                    optimize!(m)
                    solve_status = JuMP.termination_status(m)
                    @test (
                        solve_status in
                        [MathOptInterface.INFEASIBLE, MathOptInterface.INFEASIBLE_OR_UNBOUNDED]
                    )
                end
            end
        end
    end

    @testset "Bounds" begin
        @testset "Bounds on variables" begin
            m = TestHelpers.get_new_model()
            x = @variable(m, [i = 1:2], lower_bound = -1, upper_bound = 1)

            A1 = [1 -0.5; -0.5 1]
            b1 = [0, 0]
            p1 = Linear(A1, b1)

            A2 = [1 -1; 1 -1]
            b2 = [0, 0]
            p2 = Linear(A2, b2)

            test_cases = [(interval_arithmetic, -3.0, 3.0), (lp, -2.0, 2.0), (mip, -1.5, 1.5)]

            for (algorithm, l, u) in test_cases
                @testset "tightening with $(algorithm)" begin
                    m.ext[:MIPVerify] = MIPVerify.MIPVerifyExt(algorithm)
                    output = (x |> p1 |> ReLU() |> p2)

                    @test tight_upperbound(output[1], nta = algorithm) ≈ u
                    @test tight_lowerbound(output[2], nta = algorithm) ≈ l
                end
            end
        end

        @testset "Bounds on constants" begin
            x1 = one(JuMP.VariableRef)
            @test tight_upperbound(x1) == 1
            @test tight_lowerbound(x1) == 1
        end

    end

    @testset "get_relu_type" begin
        @test MIPVerify.get_relu_type(-3, -1) == MIPVerify.zero_output
        @test MIPVerify.get_relu_type(0, 0) == MIPVerify.zero_output # special case; this can go either way
        @test MIPVerify.get_relu_type(1, 1) == MIPVerify.constant_output
        @test MIPVerify.get_relu_type(-3, 5) == MIPVerify.split
        @test MIPVerify.get_relu_type(3, 5) == MIPVerify.linear_in_input
    end
end

max_stars_repo_name: vtjeng/MIPVerify

max_issues_repo_name: vtjeng/MIPVerify

max_forks_repo_name: vtjeng/MIPVerify


text:
import Base: split,
             ∈
import Random: rand
import IntervalArithmetic: inf, sup

"""
    AbstractIntervalMatrix{IT} <: AbstractMatrix{IT}

Abstract supertype for interval matrix types.
"""
abstract type AbstractIntervalMatrix{IT} <: AbstractMatrix{IT} end

"""
    IntervalMatrix{T, IT<:Interval{T}, MT<:AbstractMatrix{IT}} <: AbstractIntervalMatrix{IT}

An interval matrix i.e. a matrix whose coefficients are intervals. This type is
parametrized in the number field, the interval type, and the matrix type.

### Fields

- `mat` -- matrix whose entries are intervals

### Examples

```jldoctest
julia> A = IntervalMatrix([-0.9±0.1 0±0; 0±0 -0.9±0.1])
2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}:
 [-1, -0.799999]   [0, 0]
  [0, 0]          [-1, -0.799999]
```
"""
struct IntervalMatrix{T, IT<:Interval{T}, MT<:AbstractMatrix{IT}} <: AbstractIntervalMatrix{IT}
    mat::MT
end

import Base:size, IndexStyle, getindex, setindex!, copy

IndexStyle(::Type{<:IntervalMatrix}) = IndexLinear()
size(M::IntervalMatrix) = size(M.mat)
getindex(M::IntervalMatrix, i::Int) = getindex(M.mat, i)
setindex!(M::IntervalMatrix, X, inds...) = setindex!(M.mat, X, inds...)
copy(M::IntervalMatrix) = IntervalMatrix(copy(M.mat))

"""
    opnorm(A::IntervalMatrix, p::Real=Inf)

The matrix norm of an interval matrix.

### Input

- `A` -- interval matrix
- `p` -- (optional, default: `Inf`) the class of `p`-norm

### Notes

The matrix ``p``-norm of an interval matrix ``A`` is defined as

```math
    ‖A‖_p := ‖\\max(|\\text{inf}(A)|, |\\text{sup}(A)|)‖_p
```

where ``\\max`` and ``|·|`` are taken elementwise.
"""
function LinearAlgebra.opnorm(A::IntervalMatrix, p::Real=Inf)
    if p == Inf || p == 1
        return LinearAlgebra.opnorm(max.(abs.(inf(A)), abs.(sup(A))), p)
    else
        error("the interval matrix norm for this value of p=$p is not implemented")
    end
end

"""
    inf(A::IntervalMatrix{T}) where {T}

Return the infimum of this interval matrix, which corresponds to taking the
element-wise infimum of `A`.

### Input

- `A` -- interval matrix

### Output

A scalar matrix whose coefficients are the infima of each element in `A`.
"""
function inf(A::IntervalMatrix{T}) where {T}
    return map(inf, A)
end

"""
    sup(A::IntervalMatrix{T}) where {T}

Return the supremum of this interval matrix, which corresponds to taking the
element-wise supremum of `A`.

### Input

- `A` -- interval matrix

### Output

A scalar matrix whose coefficients are the suprema of each element in `A`.
"""
function sup(A::IntervalMatrix{T}) where {T}
    return map(sup, A)
end

"""
    split(A::IntervalMatrix{T}) where {T}

Split an interval matrix ``A`` into two scalar matrices ``C`` and ``S``
such that ``A = C + [-S, S]``.

### Input

- `A` -- interval matrix

### Output

A pair `(C, S)` such that the entries of `C` are the central points and the
entries of `S` are the (nonnegative) radii of the intervals in `A`.
"""
function split(A::IntervalMatrix{T}) where {T}
    m, n = size(A)
    C = Matrix{T}(undef, m, n)
    S = Matrix{T}(undef, m, n)

    @inbounds for j in 1:n
        for i in 1:m
            itv = A[i, j]
            radius = (sup(itv) - inf(itv)) / T(2)
            C[i, j] = inf(itv) + radius
            S[i, j] = radius
        end
    end

    return C, S
end

"""
    ∈(M::AbstractMatrix, A::AbstractIntervalMatrix)

Check whether a concrete matrix is an instance of an interval matrix.

### Input

- `M` -- concrete matrix
- `A` -- interval matrix

### Output

`true` iff `M` is an instance of `A`

### Algorithm

We check for each entry in `M` whether it belongs to the corresponding interval
in `A`.
"""
function ∈(M::AbstractMatrix, A::AbstractIntervalMatrix)
    @assert size(M) == size(A) "incompatible matrix sizes (M: $(size(M)), A: " *
                               "$(size(A)))"

    m, n = size(A)
    @inbounds for j in 1:n
        for i in 1:m
            if M[i, j] ∉ A[i, j]
                return false
            end
        end
    end
    return true
end

# random interval
function rand(::Type{Interval}; N::Type{<:Real}=Float64,
              rng::AbstractRNG=GLOBAL_RNG)::Interval{N}
    x, y = randn(rng, N), randn(rng, N)
    return x < y ? Interval(x, y) : Interval(y, x)
end

# create a random interval for the given numeric type and random-number generator
@inline function _rand_interval(; N=Float64, rng::AbstractRNG=GLOBAL_RNG)
    x, y = randn(rng, N), randn(rng, N)
    return x < y ? Interval(x, y) : Interval(y, x)
end

"""
    rand(::Type{IntervalMatrix}, m::Int=2, n::Int=2;
         N=Float64, rng::AbstractRNG=GLOBAL_RNG)

Return a random interval matrix of the given size and numeric type.

### Input

- `IntervalMatrix` -- type, used for dispatch
- `m`              -- (optional, default: `2`) number of rows
- `n`              -- (optional, default: `2`) number of columns
- `rng`            -- (optional, default: `GLOBAL_RNG`) random-number generator

### Output

An interval matrix of size ``m × n`` whose coefficients are normally-distributed
intervals of type `N` with mean `0` and standard deviation `1`.
"""
function rand(::Type{IntervalMatrix}, m::Int=2, n::Int=2;
              N=Float64, rng::AbstractRNG=GLOBAL_RNG)
    B = Matrix{Interval{N}}(undef, m, n)
    for j in 1:n
        for i in 1:m
            @inbounds B[i, j] = _rand_interval(N=N, rng=rng)
        end
    end
    return IntervalMatrix(B)
end

"""
    sample(A::IntervalMatrix{T}; rng::AbstractRNG=GLOBAL_RNG) where {T}

Return a sample of the given random interval matrix.

### Input

- `A`   -- interval matrix
- `m`   -- (optional, default: `2`) number of rows
- `n`   -- (optional, default: `2`) number of columns
- `rng` -- (optional, default: `GLOBAL_RNG`) random-number generator

### Output

An interval matrix of size ``m × n`` whose coefficients are normally-distributed
intervals of type `N` with mean `0` and standard deviation `1`.
"""
function sample(A::IntervalMatrix{T}; rng::AbstractRNG=GLOBAL_RNG) where {T}
    m, n = size(A)
    B = Matrix{T}(undef, m, n)
    @inbounds for j in 1:n
        for i in 1:m
            itv = A[i, j]
            B[i, j] = (sup(itv) - inf(itv)) * rand(rng) + inf(itv)
        end
    end
    return B
end

ext: jl


text:
using Statistics, SpecialFunctions

"""
    tnmom2c(c, a, b)

Computes ⟨(x - c)^2⟩ for the standard normal distribution truncated to the
interval [a, b].
"""
function tnmom2c(c::Real, a::Real, b::Real)
    if !(a ≤ b) || isnan(c)
        return oftype(middle(a, b), NaN)
    elseif a == b
        return oftype(middle(a, b), (a - c) * (b - c))
    elseif isinf(c)
        return c^2
    elseif abs(a) > abs(b)
        return tnmom2c(-c, -b, -a)
    elseif isinf(a) && isinf(b)
        return oftype(middle(a, b), 1 + c^2)
    elseif c ≤ 0
        m1 = tnmean(a, b)
        @assert m1 ≥ 0
        #= at this point m1 ≥ 0, so the subtraction is fine =#
        return tnmom2(a, b) + c^2 - 2c * m1
    elseif isfinite(a) && isinf(b)
        return 1 + c^2 + √(2/π) * (a - 2c) / erfcx(a / √2)
    end

    @assert a < b
    @assert 0 ≤ a ≤ b || a ≤ 0 ≤ b

    exp_a = exp(-a^2 / 2)
    exp_b = exp(-b^2 / 2)
    erfcx_abs_a = erfcx(abs(a) / √2)
    erfcx_b = erfcx(b / √2)
    erfc_abs_a = exp_a * erfcx_abs_a
    erfc_b = exp_b * erfcx_b

    ra = √(2/π) * (a - 2c) / (1 + c^2)    # a / (1 + c^2) - 2 / (1/c + c)
    rb = √(2/π) * (b - 2c) / (1 + c^2)

    Δ = (b - a) * middle(a, b)

    if a ≤ 0 ≤ b
        erfc_a = 2 - erfc_abs_a
        erfc_ma = erfc_abs_a
        erfc_mb = 2 - erfc_b
        erf_a = erfc_abs_a - 1
        erf_b = 1 - erfc_b
        @assert erf_b ≥ erf_a

        γa = ra * exp_a - erf_a
        γb = rb * exp_b - erf_b
        @assert γa ≥ γb

        if γb ≤ 0 ≤ γa
            m = γa - γb
            @assert m ≥ 0
        elseif γa ≤ 0 && γb ≤ 0
            m = (ra * exp_a + erfc_a) - (rb * exp_b + erfc_b)
            @assert m ≥ 0
        elseif γa ≥ 0 && γb ≥ 0
            m = (ra * exp_a - erfc_ma) - (rb * exp_b - erfc_mb)
            @assert m ≥ 0
        end
        m /= (erf_b - erf_a)
    elseif 0 ≤ a ≤ b
        erfc_a = erfc_abs_a
        erfc_ma = 2 - erfc_abs_a
        erfc_mb = 2 - erfc_b
        erf_a = 1 - erfc_a
        erf_b = 1 - erfc_b
        erfcx_a = erfcx_abs_a

        γa = ra * exp_a - erf_a
        γb = rb * exp_b - erf_b
        @assert γa ≥ γb

        if γa ≥ 0 ≥ γb
            m = (γa - γb) / (erfc_a - erfc_b)
            @assert m ≥ 0
        elseif γa ≤ 0 && γb ≤ 0 # ζ ≤ a ≤ b
            exp_Δ = exp(-Δ)
            m = (ra + erfcx_a) - (rb + erfcx_b) * exp_Δ
            @assert !isnan(m)
            @assert m ≥ 0
            m /= erfcx_a - erfcx_b * exp_Δ
            @assert m ≥ 0
        elseif γa ≥ 0 && γb ≥ 0
            m = (ra * exp_a - erfc_ma) - (rb * exp_b - erfc_mb)
            m /= erfc_a - erfc_b
            @assert m ≥ 0
        else
            error("trouble")
        end
    else
        error("trouble")
    end
    @assert m ≥ 0
    return (1 + c^2) * m

    # γa = √(2/π) * (a - 2c) / (1 + c^2) * exp(-a^2 / 2) - ea
    # γb = √(2/π) * (b - 2c) / (1 + c^2) * exp(-b^2 / 2) - eb
    #
    # ea = erf(a / √2)
    # eb = erf(b / √2)
    #
    # if γb ≤ 0 ≤ γa
    #     return (1 + c^2) * (γa - γb) / (eb - ea)
    # else
    #     ea = erfcx(a / √2)
    #     eb = erfcx(b / √2)
    #     γa = √(2/pi) * (a - 2c) / (1 + c^2) + ea
    #     γb = √(2/pi) * (b - 2c) / (1 + c^2) + eb
    #     @assert 0 ≤ Δ ≤ 1
    #     #@assert γb ≤ γa
    #     @assert eb ≤ ea
    #     return (1 + c^2) * (γa - γb * Δ) / (ea - eb * Δ)
    # end
end

"""
    tnmom2c(c, a, b, μ, σ)

Computes ⟨(x - c)^2⟩ for the normal distribution with mean μ and standard
deviation σ, truncated to the interval [a, b].
"""
function tnmom2c(c, a, b, μ, σ)
    α = (a - μ) / σ
    β = (b - μ) / σ
    γ = (c - μ) / σ
    return σ^2 * tnmom2c(γ, α, β)
end


text:
# Extend later to multiple dimensions by using multple dispatch, probably we'll now use a Vector of Tuples of Ints
# or by using generated functions
function init_derivatives!(m::AbstractNLCTModel, instructions::AbstractDict{Symbol, Vector{Int}},
                           statevars::Vector{Symbol} = collect(keys(get_state_variables(m))))

    derivs = get_derivatives(m)
    state = statevars[1]

    n = length(derivs) # Get existing length of derivatives that have already been added
    for (k, v) in instructions
        for i in v
            n += 1
            if i == 1
                derivs[Symbol(:∂, k, :_∂, state)] = n
            elseif i == 2
                derivs[Symbol(:∂², k, :_∂, state, "²")] = n
            end
        end
    end
end

function init_derivatives!(m::AbstractNLCTModel, instructions::AbstractDict{Symbol, Vector{Tuple{Int, Int}}},
                           statevars::Vector{Symbol} = collect(keys(get_state_variables(m))))

    derivs = get_derivatives(m)
    state1 = statevars[1]
    state2 = statevars[2]

    n = length(derivs) # Get existing length of derivatives that have already been added
    for (k, v) in instructions
        for i in v
            n += 1
            if i == (1, 0)
                derivs[Symbol(:∂, k, :_∂, state1)] = n
            elseif i == (2, 0)
                derivs[Symbol(:∂², k, :_∂, state1, "²")] = n
            elseif i == (0, 1)
                derivs[Symbol(:∂, k, :_∂, state2)] = n
            elseif i == (0, 2)
                derivs[Symbol(:∂², k, :_∂, state2, "²")] = n
            elseif i == (1, 1)
                derivs[Symbol(:∂², k, :_∂, state1, :∂, state2)] = n
            end
        end
    end
end

"""
```
function standard_derivs(dims::Int)
```

returns instructions for which derivatives
to calculate for standard continuous time models,
depending on the dimension `dims` of the state space.

For a one-dimensional model, we request the first and second derivatives.

For a multi-dimensional model, we request the first, second,
and (first) mixed partial derivatives, e.g.
for a 2D model, we want ∂f_∂x, ∂f²_∂x2, ∂f²_∂x∂y, ∂f_∂y, ∂f²_∂y2.

The instructions are returned as a Vector of Ints or Vector of Tuples of Ints.
"""
function standard_derivs(dims::Int)
    if dims == 1
        v = Vector{Int}(undef, 2)
        v[1] = 1
        v[2] = 2
        return v
    elseif dims == 2
        v = Vector{Tuple{Int, Int}}(undef, 5)
        v[1] = (1, 0)
        v[2] = (0, 1)
        v[3] = (2, 0)
        v[4] = (1, 1)
        v[5] = (0, 2)
        return v
    end
end


text:
export LSFManager, addprocs_lsf

struct LSFManager <: ClusterManager
    np::Integer
    bsub_flags::Cmd
end

function launch(manager::LSFManager, params::Dict, launched::Array, c::Condition)
    try
        dir = params[:dir]
        exename = params[:exename]
        exeflags = params[:exeflags]

        np = manager.np

        jobname = `julia-$(getpid())`

        cmd = `$exename $exeflags $(worker_arg())`
        bsub_cmd = `bsub -I $(manager.bsub_flags) -cwd $dir -J $jobname "$cmd"`

        stream_proc = [open(bsub_cmd) for i in 1:np]

        for i in 1:np
            config = WorkerConfig()
            config.io = stream_proc[i]
            push!(launched, config)
            notify(c)
        end

    catch e
        println("Error launching workers")
        println(e)
    end
end

manage(manager::LSFManager, id::Int64, config::WorkerConfig, op::Symbol) = nothing

kill(manager::LSFManager, id::Int64, config::WorkerConfig) = kill(config.io)

addprocs_lsf(np::Integer, bsub_flags::Cmd=``; params...) =
        addprocs(LSFManager(np, bsub_flags); params...)


text:
# Example that computes the full inverse wavelet transform using all the
# scales from Jmin=1, to Jmax=log2(n)-1 the coarsest scale

#push path
push!(LOAD_PATH,pwd()*"/src") # for julia 0.4 you need to put your path
#import libraries
using FWT
using PyPlot
using Images

n = 1024;
name = "data_samples/ernst_reuter_haus.bmp";
f = load_image(name, n);
f = rescale(sum(f,3));
f = f[:,:,1];

# filter call
h=filt_gen(WT.db2)
# high pass filter
g = cat(1, 0, h[length(h):-1:2]) .* ( (-1).^(1:length(h)) )
g = cat(1, 0, h[length(h):-1:2]) .* ( (-1).^(1:length(h)) )

# First compute wavelet transform
# Perform of the wavelet transform
Jmax = round(Int64,log2(n))-1;
Jmin = 1;
fW = copy(f);
clf;
for j=Jmax:-1:Jmin
    A = fW[1:2^(j+1),1:2^(j+1)];
    for d=1:2
        Coarse = subsampling(cconvol(A,h,d),d);
        Detail = subsampling(cconvol(A,g,d),d);
        A = cat(d, Coarse, Detail );
    end
    fW[1:2^(j+1),1:2^(j+1)] = A;
    j1 = Jmax-j;
    if j1<4
        imageplot(A[1:2^j,2^j+1:2^(j+1)], "Horizontal, j=$j", 3,4, j1 + 1);
        imageplot(A[2^j+1:2^(j+1),1:2^j], "Vertical, j=$j", 3,4, j1 + 5);
        imageplot(A[2^j+1:2^(j+1),2^j+1:2^(j+1)], "Diagonal, j=$j", 3,4, j1 + 9);
    end
end

# Now the inverse

f1 = copy(fW);
clf;
for j=Jmin:Jmax
    A = f1[1:2^(j+1),1:2^(j+1)];
    for d=1:2
        if d==1
            Coarse = A[1:2^j,:];
            Detail = A[2^j+1:2^(j+1),:];
        else
            Coarse = A[:,1:2^j];
            Detail = A[:,2^j+1:2^(j+1)];
        end
        Coarse = cconvol(upsampling(Coarse,d),FWT.reverse(h),d);
        Detail = cconvol(upsampling(Detail,d),FWT.reverse(g),d);
        A = Coarse + Detail;
        j1 = Jmax-j;
        if j1>0 && j1<5
            imageplot(A, "Partial reconstruction, j=$j", 2,2,j1);
        end
    end
    f1[1:2^(j+1),1:2^(j+1)] = A;
end


text:
# TODO
# Only for 1dof joints
mutable struct TrackingLQR{T,N,NK} <: Controller
    K::Vector{Vector{SMatrix{1,NK,T,NK}}} # for each time step and each eqc

    xd::Vector{Vector{SVector{3,Float64}}} # for each time step and each eqc
    vd::Vector{Vector{SVector{3,Float64}}}# for each time step and each eqc
    qd::Vector{Vector{UnitQuaternion{T}}} # for each time step and each eqc
    ωd::Vector{Vector{SVector{3,Float64}}} # for each time step and each eqc

    eqcids::Vector{Integer}
    Fτd::Vector{Vector{SVector{1,Float64}}} # for each time step and each eqc

    control!::Function

    function TrackingLQR(mechanism::Mechanism{T,Nn,Nb}, storage::Storage{T,N}, Fτ::Vector{<:Vector{<:AbstractVector{T}}}, eqcids::AbstractVector{<:Integer},
            Q::Vector{<:AbstractMatrix{T}}, R::Vector{<:AbstractMatrix{T}};
            controlfunction::Function = control_trackinglqr!
        ) where {T, Nn, Nb, N}

        Q = cat(Q...,dims=(1,2))*mechanism.Δt
        R = cat(R...,dims=(1,2))*mechanism.Δt

        xd = [[SA{T}[0; 0; 0] for i=1:Nb] for j=1:N]
        vd = [[SA{T}[0; 0; 0] for i=1:Nb] for j=1:N]
        qd = [[one(UnitQuaternion{T}) for i=1:Nb] for j=1:N]
        ωd = [[SA{T}[0; 0; 0] for i=1:Nb] for j=1:N]

        for k = 1:N
            for i = 1:Nb
                xd[k][i] = storage.x[i][k]
                vd[k][i] = storage.v[i][k]
                qd[k][i] = storage.q[i][k]
                ωd[k][i] = storage.ω[i][k]
            end
        end

        # calculate K
        Ku = dlqr(mechanism, xd, vd, qd, ωd, Fτ, eqcids,Q,R,N)

        new{T, N, size(Ku[1][1])[2]}(Ku, xd, vd, qd, ωd, eqcids, Fτ, controlfunction)
    end
end

function control_trackinglqr!(mechanism::Mechanism{T,Nn,Nb}, lqr::TrackingLQR{T,N}, k) where {T,Nn,Nb,N}
    Δz = zeros(T,Nb*12)
    qvm = QuatVecMap()
    for (id,body) in pairs(mechanism.bodies)
        colx = (id-1)*12+1:(id-1)*12+3
        colv = (id-1)*12+4:(id-1)*12+6
        colq = (id-1)*12+7:(id-1)*12+9
        colω = (id-1)*12+10:(id-1)*12+12

        state = body.state
        Δz[colx] = state.xsol[2]-lqr.xd[k][id]
        Δz[colv] = state.vsol[2]-lqr.vd[k][id]
        Δz[colq] = rotation_error(state.qsol[2],lqr.qd[k][id],qvm)
        Δz[colω] = state.ωsol[2]-lqr.ωd[k][id]
    end

    if k<N
        for (i,id) in enumerate(lqr.eqcids)
            u = lqr.Fτd[k][i] - lqr.K[k][i]*Δz
            setForce!(mechanism, geteqconstraint(mechanism, id), u)
        end
    end

    return
end

function dlqr(mechanism::Mechanism{T,Nn,Nb,Ne}, xd, vd, qd, ωd, Fτd, eqcids,Q,R,N) where {T,Nn,Nb,Ne}
    bodyids = getid.(mechanism.bodies)
    mx = Nb*6*2
    mu = size(eqcids)[1]
    mλ = 0
    for eqc in mechanism.eqconstraints
        ConstrainedDynamics.isinactive(eqc) && continue
        mλ += ConstrainedDynamics.length(eqc)
    end
    Ku = [[zeros(1,size(Q)[1]) for j=1:mu] for i=1:(N-1)]
    Kλ = [[zeros(1,size(Q)[1]) for j=1:mλ] for i=1:(N-1)]
    Pk = Q

    k = 0
    for outer k=N-1:-1:1
        A, Bu, Bλ, G = linearsystem(mechanism, xd[k], vd[k], qd[k], ωd[k], Fτd[k], bodyids, eqcids)
        D = Bu - Bλ/(G*Bλ)*G*Bu
        M11 = R + D'*Pk*Bu
        M12 = D'*Pk*Bλ
        M21 = G*Bu
        M22 = G*Bλ

        M = [M11 M12;M21 M22]
        b = [D'*Pk;G]*A

        Kk = M\b

        for i=1:mu
            Ku[k][i] = Kk[i:i,:]
        end

        Kuk = Kk[1:mu,:]
        Kλk = Kk[mu+1:mu+mλ,:]

        Abar = A-Bu*Kuk-Bλ*Kλk
        Pkp1 = Q + Kuk'*R*Kuk + Abar'*Pk*Abar

        if norm(Pk-Pkp1) < 1e-5
            break
        end

        Pk = Pkp1
    end

    for k2=k-1:-1:1
        Ku[k2] = Ku[k2+1]
    end

    return Ku
end


text:
"""
Vanilla RNN, i.e. ⤦\n
    hᵗ = f(w*xᵗ + u*hᵗ⁻¹ .+ b)
"""
mutable struct RNN <: Block
    w::VarOrNil # input to hidden weights
    b::VarOrNil # bias of hidden units
    u::VarOrNil # recurrent weights
    f::Function # activation function
    h::Any      # hidden variable
    function RNN(inputSize::Int, hiddenSize::Int, fn::Function=relu; type::Type=Array{Float32})
        T = eltype(type)
        w = randn(T, hiddenSize, inputSize) .* sqrt( T(2/inputSize) )
        b = zeros(T, hiddenSize, 1)
        u = randdiagonal(T, hiddenSize; from=-0.1990, to=0.1997)
        new(Variable{type}(w,true,true,true),
            Variable{type}(b,true,true,true),
            Variable{type}(u,true,true,true), fn, nothing)
    end
    function RNN(fn::Function)
        new(nothing, nothing, nothing, fn, nothing)
    end
end

function clone(this::RNN; type::Type=Array{Float32})
    cloned = RNN(this.f)
    cloned.w = clone(this.w, type=type)
    cloned.b = clone(this.b, type=type)
    cloned.u = clone(this.u, type=type)
    return cloned
end

mutable struct RNNs <: Block
    layers::Vector{RNN}
    function RNNs(topology::Vector{Int}, fn::Array{F}; type::Type=Array{Float32}) where F
        n = length(topology) - 1
        layers = Vector{RNN}(undef, n)
        for i = 1:n
            layers[i] = RNN(topology[i], topology[i+1], fn[i]; type=type)
        end
        new(layers)
    end
end

Base.getindex(m::RNNs,     k...) =  m.layers[k...]
Base.setindex!(m::RNNs, v, k...) = (m.layers[k...] = v)
Base.length(m::RNNs)       = length(m.layers)
Base.lastindex(m::RNNs)    = length(m.layers)
Base.firstindex(m::RNNs)   = 1
Base.iterate(m::RNNs, i=firstindex(m)) = i>length(m) ? nothing : (m[i], i+1)

function Base.show(io::IO, m::RNN)
    SIZE = size(m.w)
    TYPE = typeof(m.w.value)
    print(io, "RNN($(SIZE[2]), $(SIZE[1]), $(m.f); type=$TYPE)")
end

function Base.show(io::IO, m::RNNs)
    print(io, "RNNs\n      (\n          ")
    join(io, m.layers, ",\n          ")
    print(io, "\n      )")
end

function resethidden(m::RNN)
    m.h = nothing
end

function resethidden(model::RNNs)
    for m in model
        resethidden(m)
    end
end

function forward(m::RNN, x::Variable{T}) where T
    f = m.f  # activition function
    w = m.w  # input's weights
    b = m.b  # input's bias
    u = m.u  # memory's weights
    h = m.h ≠ nothing ? m.h : Variable{T}(Zeros(T, size(w,1), size(x,2)))
    x = f(matAddVec(w*x + u*h, b))
    m.h = x
    return x
end

function forward(model::RNNs, x::Variable)
    for m in model
        x = forward(m, x)
    end
    return x
end

function predict(m::RNN, x::T) where T
    f = m.f        # activition function
    w = m.w.value  # input's weights
    b = m.b.value  # input's bias
    u = m.u.value  # memory's weights
    h = m.h ≠ nothing ? m.h : Zeros(T, size(w,1), size(x,2))
    x = f(w*x + u*h .+ b)
    m.h = x
    return x
end

function predict(model::RNNs, x)
    for m in model
        x = predict(m, x)
    end
    return x
end

"""
    unbiasedof(m::RNN)

unbiased weights of RNN block
"""
function unbiasedof(m::RNN)
    weights = Vector(undef, 1)
    weights[1] = m.w.value
    return weights
end

function weightsof(m::RNN)
    weights = Vector(undef,3)
    weights[1] = m.w.value
    weights[2] = m.b.value
    weights[3] = m.u.value
    return weights
end

"""
    unbiasedof(model::RNNs)

unbiased weights of RNNs block
"""
function unbiasedof(model::RNNs)
    weights = Vector(undef, 0)
    for m in model
        append!(weights, unbiasedof(m))
    end
    return weights
end

function weightsof(m::RNNs)
    weights = Vector(undef,0)
    for i = 1:length(m)
        append!(weights, weightsof(m[i]))
    end
    return weights
end

function gradsof(m::RNN)
    grads = Vector(undef,3)
    grads[1] = m.w.delta
    grads[2] = m.b.delta
    grads[3] = m.u.delta
    return grads
end

function gradsof(m::RNNs)
    grads = Vector(undef,0)
    for i = 1:length(m)
        append!(grads, gradsof(m[i]))
    end
    return grads
end

function zerograds!(m::RNN)
    for v in gradsof(m)
        v .= zero(v)
    end
end

function zerograds!(m::RNNs)
    for v in gradsof(m)
        v .= zero(v)
    end
end

function paramsof(m::RNN)
    params = Vector{Variable}(undef,3)
    params[1] = m.w
    params[2] = m.b
    params[3] = m.u
    return params
end

function xparamsof(m::RNN)
    xparams = Vector{XVariable}(undef,3)
    xparams[1] = ('w', m.w)
    xparams[2] = ('b', m.b)
    xparams[3] = ('u', m.u)
    return xparams
end

function paramsof(m::RNNs)
    params = Vector{Variable}(undef,0)
    for i = 1:length(m)
        append!(params, paramsof(m[i]))
    end
    return params
end

function xparamsof(m::RNNs)
    xparams = Vector{XVariable}(undef,0)
    for i = 1:length(m)
        append!(xparams, xparamsof(m[i]))
    end
    return xparams
end

function nparamsof(m::RNN)
    lw = length(m.w)
    lb = length(m.b)
    lu = length(m.u)
    return (lw + lb + lu)
end

function bytesof(model::RNN, unit::String="MB")
    n = nparamsof(model) * elsizeof(model.w)
    return blocksize(n, uppercase(unit))
end

function nparamsof(m::RNNs)
    num = 0
    for i = 1:length(m)
        num += nparamsof(m[i])
    end
    return num
end

function bytesof(model::RNNs, unit::String="MB")
    n = nparamsof(model) * elsizeof(model[1].w)
    return blocksize(n, uppercase(unit))
end

function to(type::Type, m::RNN)
    m.w = to(type, m.w)
    m.b = to(type, m.b)
    m.u = to(type, m.u)
    return m
end

function to!(type::Type, m::RNN)
    m = to(type, m)
    return nothing
end

function to(type::Type, m::RNNs)
    for layer in m
        layer = to(type, layer)
    end
    return m
end

function to!(type::Type, m::RNNs)
    for layer in m
        to!(type, layer)
    end
end

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
"""
    STD_Linear_RR!(X::Vector{Interval{Float64}},opt,UBD::Float64)

Performs standard range reduction using linear underestimator forms via
McCormick relaxations and contracts `X::Vector{Interval{Float64}}` in place.
The `opt` object has the fields `.numVar`, `.numConstr`, `.f`, `.g`, `.gL`, `.gU`,
`.gL_Loc`, `.gU_Loc` and `solver.LP_solver`. `f(x)` is the objective, `g(x)` is
the constraint function, `numVar` is the number of decision variables, `numConstr`
is the number of constraints, `gL` is the lower bound, `gU` is the upper bound,
`gL_Loc` is an index to indicate the lower bound is finite, and `gU_Loc` is an index
that indicates the upper bound is finite. The upper bound is `UBD::Float64`.
"""
function STD_Linear_RR!(X::Vector{Interval{Float64}},opt,UBD::Float64)

    # gets prior McCormick relaxation setting
    mu_temp::Int64 = copy(EAGOSmoothMcCormickGrad.MC_param.mu)
    set_diff_relax(0)

    # sets variable bounds
    l::Vector{Float64} = [X[i].lo for i=1:opt[1].numVar]
    u::Vector{Float64} = [X[i].hi for i=1:opt[1].numVar]
    x0::Vector{Float64} = mid.(X)

    # computes relaxation of g and f
    x_mc::Vector{SMCg{opt[1].numVar,Float64}} = [SMCg{opt[1].numVar,Float64}(x0[i],x0[i],seed_g(Float64,i,opt[1].numVar),seed_g(Float64,i,opt[1].numVar),X[i],false,X,x0) for i=1:opt[1].numVar]
    f_mc::SMCg{opt[1].numVar,Float64} = opt[1].f(x_mc)
    if opt[1].numConstr>0
        c::Vector{SMCg{opt[1].numVar,Float64}} = opt[1].g(x_mc)
        dcdx::SparseMatrixCSC{Float64,Int64} = spzeros(length(opt[1].gL_loc)+length(opt[1].gU_loc)+1,opt[1].numVar)
    else
        dcdx = spzeros(1,opt[1].numVar)
    end
    # forms coefficients of linear relaxations
    dcdx[1,:] = f_mc.cv_grad
    if opt[1].numConstr>0
        cx_ind1::Int64 = 2
        for i in opt[1].gL_loc
            for j=1:opt[1].numVar
                if (c[i].cv_grad[j] != 0.0)
                    dcdx[cx_ind1,j] = c[i].cv_grad[j]
                end
            end
            cx_ind1 += 1
        end
        for i in opt[1].gU_loc
            for j=1:opt[1].numVar
                if (c[i].cc_grad[j] != 0.0)
                    dcdx[cx_ind1,j] = -c[i].cc_grad[j]
                end
            end
            cx_ind1 += 1
        end
    end
    # forms rhs of linear relaxations
    if opt[1].numConstr>0
        rhs::Vector{Float64} = zeros(Float64,length(opt[1].gL_loc)+length(opt[1].gU_loc)+1)
    else
        rhs = zeros(Float64,1)
    end
    rhs[1] = UBD + f_mc.cv - sum(x0.*f_mc.cv_grad)

    if opt[1].numConstr>0
        cx_ind2::Int64 = 2
        for i in opt[1].gU_loc
            rhs[cx_ind2] = sum(x0[:].*c[i].cv_grad[:])+opt[1].gU[i]-c[i].cv
            cx_ind2 += 1
        end
        for i in opt[1].gL_loc
            rhs[cx_ind2] = sum(-x0[:].*c[i].cc_grad[:])-opt[1].gL[i]+c[i].cc
            cx_ind2 += 1
        end
    end

    # solve nx lower & upper bounding problems
    for i=1:opt[1].numVar

        # solve lower bounding problems
        c_obj::Vector{Float64} = [j == i ? 1.0 : 0.0 for j=1:opt[1].numVar]
        model = buildlp(c_obj, dcdx, '<', rhs, l, u, opt[1].solver.LP_solver)
        MathProgBase.optimize!(model)
        stat_flag::Symbol = MathProgBase.status(model)
        if (stat_flag == :Optimal)
            l[i] = MathProgBase.getobjval(model)
        elseif (stat_flag == :Infeasible)
            feas = false
            return feas
        end

        # solve upper bounding problems
        c_obj[i] = -1.0
        model = buildlp(c_obj, dcdx, '<', rhs, l, u, opt[1].solver.LP_solver)
        MathProgBase.optimize!(model)
        stat_flag = MathProgBase.status(model)
        if (stat_flag == :Optimal)
            u[i] = -MathProgBase.getobjval(model)
        elseif (stat_flag == :Infeasible)
            feas = false
            return feas
        end
    end

    # reset McCormick relaxation to original value
    set_diff_relax(mu_temp)

    # stores outputs
    X[:] = [Interval(l[i],u[i]) for i=1:opt[1].numVar]
    return true
end

"""
    Imp_Linear_RR!(X::Vector{Interval{Float64}},opt,UBD::Float64)

Performs standard range reduction using linear underestimator forms via
McCormick relaxations and contracts `X::Vector{Interval{Float64}}` in place. However,
it only contracts the p variables. The `opt` object has the fields `.numVar`, `.numConstr`, `.f`, `.g`, `.gL`, `.gU`,
`.gL_Loc`, `.gU_Loc`, `solver.LP_solver`, and `solver.Implicit_Opts.nx`. `f(x)` is the objective, `g(x)` is
the constraint function, `numVar` is the number of decision variables, `numConstr`
is the number of constraints, `gL` is the lower bound, `gU` is the upper bound,
`gL_Loc` is an index to indicate the lower bound is finite, and `gU_Loc` is an index
that indicates the upper bound is finite. The upper bound is `UBD::Float64`.
"""
function Imp_Linear_RR!(X::Vector{Interval{Float64}},opt, UBD::Float64)
    # gets prior McCormick relaxation setting
    mu_temp::Int64 = copy(EAGOSmoothMcCormickGrad.MC_param.mu)
    set_diff_relax(0)

    # sets variable bounds
    nx::Int64 = opt[1].solver.Implicit_Options.nx
    np::Int64 = opt[1].numVar - nx

    l::Vector{Float64} = [X[i].lo for i=(nx+1):(nx+np)]
    u::Vector{Float64} = [X[i].hi for i=(nx+1):(nx+np)]
    x0::Vector{Float64} = (l+u)/2.0

    # computes relaxation of g and f
    p_mc::Vector{SMCg{np,Float64}} = [SMCg{np,Float64}(x0[i],x0[i],seed_g(Float64,opt[1].numVar,i),seed_g(Float64,opt[1].numVar,i),X[i],false,X,x0) for i=1:opt[1].numVar]
    if (opt[1].solver.Implicit_Options.Inplace)
        out,z,x_mc = InGenExpansionParams(opt[1].solver.Implicit_Options.h,
                                          opt[1].solver.Implicit_Options.hj,
                                          X[1:nx],
                                          X[(nx+1):(nx+np)],
                                          p0,
                                          opt[1].solver.Implicit_Options.Param)
    else
        out,z,x_mc = GenExpansionParams(opt[1].solver.Implicit_Options.h,
                                        opt[1].solver.Implicit_Options.hj,
                                        X[1:nx],
                                        X[(nx+1):(nx+np)],
                                        p0,
                                        opt[1].solver.Implicit_Options.Param)
    end
    f_mc::SMCg{np,Float64} = opt[1].solver.ImplicitOpts.f(x_mc[end],p_SMC)
    c::Vector{SMCg{np,Float64}} = opt[1].solver.ImplicitOpts.g(x_mc[end],p_SMC)

    # UPDATE gL_loc & gU_loc
    gL_temp = filter(x->(x>nx),opt[1].gL_loc)
    gU_temp = filter(x->(x>nx),opt[1].gU_loc)

    c_cv::Vector{Float64} = [c[i].cv for i=1:opt[1].numConstr]
    c_cc::Vector{Float64} = [c[i].cc for i=1:opt[1].numConstr]
    dcdx_cv::Array{Float64,2} = Float64[c[i].cv_grad[j] for i=1:opt[1].numConstr,j=1:np]
    dcdx_cc::Array{Float64,2} = Float64[-c[i].cc_grad[j] for i=1:opt[1].numConstr,j=1:np]
    cval = vcat(c_cv[gU_temp-nx]-opt[1].gU[gU_temp-nx],-c_cc[gL_temp-nx]+opt[1].gL[gL_temp-nx])
    dcdx::Array{Float64,2} = vcat([f_mc.cv_grad[i] for i=1:np]',
                                  dcdx_cv[opt[1].gU_loc,:],
                                  dcdx_cc[opt[1].gL_loc,:])

    rhs1::Vector{Float64} = [sum([x0[j]*dcdx_cv[i,j] for j=1:np]) for i=1:opt[1].numConstr]
    rhs2::Vector{Float64} = [sum([x0[j]*dcdx_cc[i,j] for j=1:np]) for i=1:opt[1].numConstr]
    rhs = vcat([UBD + f_mc.cv - sum([x0[i]*f_mc.cv_grad[i] for i=1:opt[1].numVar])],
               rhs1[gU_temp-nx]+opt[1].gU[gU_temp]-c_cv[gU_temp-nx],
               rhs2[gL_temp-nx]-opt[1].gL[gL_temp]+c_cc[gL_temp-nx])

    # solve nx lower & upper bounding problems
    for i=1:opt[1].numVar
        # solve lower bounding problems
        c_obj::Vector{Float64} = [j == i ? 1.0 : 0.0 for j=1:np]
        model = buildlp(c_obj, dcdx, '<', rhs, l, u, opt[1].solver.LP_solver)
        MathProgBase.optimize!(model)
        stat_flag::Symbol = MathProgBase.status(model)
        if (stat_flag == :Optimal)
            l[i] = MathProgBase.getobjval(model)
        elseif (stat_flag == :Infeasible)
            feas = false
            return feas
        end
        # solve upper bounding problems
        c_obj[i] = -1.0
        model = buildlp(c_obj, dcdx, '<', rhs, l, u, opt[1].solver.LP_solver)
        MathProgBase.optimize!(model)
        stat_flag = MathProgBase.status(model)
        if (stat_flag == :Optimal)
            u[i] = -MathProgBase.getobjval(model)
        elseif (stat_flag == :Infeasible)
            feas = false
            return feas
        end
    end

    # reset McCormick relaxation to original value
    set_diff_relax(mu_temp)

    X[(nx+1):(nx+np)] = [Interval(l[i],u[i]) for i=1:np]
    return true
end


text:
# Spinning sideband intensity calculations based on
# J. Herzfeld and A. E. Berger, J. Chem. Phys. 1980, 73, 6021–6030.

module SolidStateNMR

using Cubature

# τ functions
τ₋(α,β,θ) = 1/24*cos(2α)*(3 + cos(2β))*sin(2θ) -
            1/6*sin(2α)*cos(β)*cos(2θ) +
            √2/6*cos(2α)sin(2β)sin(θ) -
            √2/3*sin(2α)*sin(β)*cos(θ)
τ₊(β,θ) = 1/24*(cos(2β)-1)*sin(2θ) +
          √2/6*sin(2β)*sin(θ)

# Δ functions, σ values in ppm
Δ₋(Ω,ωᵣ,σxx,σyy) = -Ω/ωᵣ*(σxx - σyy)*1.0e-6
Δ₊(Ω,ωᵣ,σiso,σzz) = -3Ω/ωᵣ*(σiso - σzz)*1.0e-6

# Integrand to F function
f(θ,N,α,β,Δ₊,Δ₋) = 1/(2π)*exp(im*(-N*θ + Δ₋*τ₋(α,β,θ) + Δ₊*τ₊(β,θ)))

function F(N,α,β,Δ₊,Δ₋;maxevals=150)
    (val1,err1) = pquadrature(θ->real(f(θ,N,α,β,Δ₊,Δ₋)),0.0,2π;maxevals=maxevals)
    (val2,err2) = pquadrature(θ->imag(f(θ,N,α,β,Δ₊,Δ₋)),0.0,2π;maxevals=maxevals)
    return val1+im*val2
end

# Relative intensity of the Nth sideband
intensity(N,Δ₊,Δ₋;maxevals=1000) = hcubature(x->sin(x[2])*abs2(F(N,x[1],x[2],Δ₊,Δ₋)),[0,0],[2π,π];maxevals=maxevals)[1]/4π

# List of sideband intensities:
# -maxN, -maxN + 1, ..., -1, 0, 1, ..., maxN - 1, maxN
# σ: eigenvalues of the shielding tensor in ppm
# ωᵣ: rotor speed in Hz
# Ω: B₀ in Hz
# maxN: number of sidebands calculated on either side of the main peak
function sidebands(Ω,ωᵣ,σ; maxN=3, maxevals=500)
    σiso = (σ[1] + σ[2] + σ[3])/3
    D₊ = Δ₊(Ω,ωᵣ,σiso,σ[3])
    D₋ = Δ₋(Ω,ωᵣ,σ[1],σ[2])
    Is = [intensity(N, D₊, D₋;maxevals=maxevals) for N=0:maxN]
    [reverse(Is[2:end]);Is]  # I[-n] = I[n]
end

# Lorentzian function
# a: peak maximum
# Γ: full width at half height
# x₀: center
# Maximum: 2a/(π*Γ)
L(a,x₀,Γ) = x -> Γ*a/(2π*((x-x₀)^2+(Γ/2)^2))

# Composition of Lorentzian functions
# a: amplitudes
# f₀: center frequencies
# w: peak widths at half height
function compose(a,f₀,w)
    funcs = [L(x[1],x[2],x[3]) for x=zip(a,f₀,w)]
    f -> sum(l(f) for l=funcs)
end

# Solid state NMR spectrum for a single atom
# σ: Eigenvalues of the shielding tensor in ppm
# ref: chemical shift reference (used for spectrum alignment)
# Ω: B₀ in Hz
# ωᵣ: rotor speed in Hz
# w: peak width at half height
# maxN: number of sidebands calculated on either side of the main peak
function spectrum{T<:Number}(σ::Array{T,1},ref,Ω,ωᵣ,w; maxN=3)
    ref_freq = ref*1.0e-6*Ω # reference in Hz
    center = ref - mean(σ)*1.0e-6*Ω # main peak location in Hz
    a = sidebands(Ω,ωᵣ,σ; maxN=maxN)
    f₀ = [center+(ωᵣ*n) for n=-maxN:maxN]
    Γ = w * ones(2*maxN+1)
    compose(a,f₀,Γ)
end

# Solid state NMR spectrum for a collection of atoms
# (This is the function you most likely want to use)
# atoms: [("atom_label1", σ₁), ("atom_label1", σ₂) , ...]
# where σᵢ is the 3-element array of the shielding tensor
# eigenvalues for atom i
# ref: chemical shift reference (used for spectrum alignment)
# Ω: B₀ in Hz
# ωᵣ: rotor speed in Hz
# w: peak width at half height
# maxN: number of sidebands calculated on either side of the main peak
function spectrum{T<:Number,S<:AbstractString}(atoms::Array{Tuple{S,Array{T,1}},1},ref_ppm,Ω,ωᵣ,w; maxN=3)
    spects = pmap(atom->spectrum(atom[2],ref_ppm,Ω,ωᵣ,w; maxN=maxN),atoms)
    x->sum(s(x) for s=spects)
end
end


text:
include("general/pogm_restart.jl")
export pogm_restart

include("general/ncg.jl")
export ncg

include("general/ogm_ls.jl")
export ogm_ls

include("general/poweriter.jl")
export poweriter


text:
function Well1024a()
    return Well1024a(())
end

function Well1024a(arg0::Vector{jint})
    return Well1024a((Vector{jint},), arg0)
end

function Well1024a(arg0::jint)
    return Well1024a((jint,), arg0)
end

function Well1024a(arg0::jlong)
    return Well1024a((jlong,), arg0)
end

function next_int(obj::Well1024a)
    return jcall(obj, "nextInt", jint, ())
end


text:
abstract type fitness_functions end

"""
    mutable struct exponential_fitness <: fitness_functions

Fitness landscape that has the form of an exponential function.
"""
mutable struct exponential_fitness <: fitness_functions
    l::Int64
    beta::Float64
    epsilon::Float64
    f0::Float64
    fl::Float64
    l0::Int64
    E_Star::Function
end

# Setting some defaults
exponential_fitness(;l=10, beta=1, epsilon=2, f0=1, fl=0, l0=5, E_Star=Est) = exponential_fitness(l, beta, epsilon, f0, fl, l0, E_Star)

"""
    function fitness(E::Float64, p::fermi_fitness)

Evaluate fermi fitness function.
"""
function _fitness(E::Real, p::exponential_fitness)::Float64
    fitness = p.f0 * (1 - exp(p.beta * (E - p.E_Star(p.l, p.epsilon, p.l0)))) - p.fl * p.l
    return fitness
end

"""
    function fitness(E::Float64, l::Int64, p::fermi_fitness)

Evaluate fermi fitness function.
"""
function _fitness(E::Real, l::Int64, p::exponential_fitness)::Float64
    fitness = p.f0 * (1 - exp(p.beta * (E - p.E_Star(l, p.epsilon, p.l0)))) - p.fl * l
    return fitness
end

"""
    mutable struct fermi_fitness <: fitness_functions

Fitness landscape that has the form of a fermi function.

# Defaults

```
fermi_fitness(;l=l, beta=beta, f0=f0) = fermi_fitness(l, beta, f0, 0, Est)
fermi_fitness(;beta=beta, f0=f0) = fermi_fitness(10, beta, f0, 0, Est)
fermi_fitness(;l=l, beta=beta, f0=f0, fl=fl) = fermi_fitness(l, beta, f0, fl, Est)
```
"""
mutable struct fermi_fitness <: fitness_functions
    l::Int64
    beta::Float64
    epsilon::Float64
    f0::Float64
    fl::Float64
    l0::Int64
    E_Star::Function
end

# Setting some defaults
fermi_fitness(;l=10, beta=1, epsilon=2, f0=1, fl=0, l0=10, E_Star=Est) = fermi_fitness(l, beta, epsilon, f0, fl, l0, E_Star)

"""
    function fitness(E::Float64, p::fermi_fitness)

Evaluate fermi fitness function.
"""
function _fitness(E::Real, p::fermi_fitness)::Float64
    fitness = (p.f0 * (1/(1 + exp(p.beta * (E - p.E_Star(p.l, p.epsilon, p.l0))))) - p.fl * p.l)
    return fitness
end

"""
    function fitness(E::Float64, p::fermi_fitness)

Evaluate fermi fitness function.
"""
function _dE_fitness(E::Real, p::fermi_fitness)
    dfitness = -p.f0 * 1/(1 + exp(p.beta * (E - p.E_Star(p.l, p.epsilon, p.l0))))^2 * p.beta * exp(p.beta * (E - p.E_Star(p.l, p.epsilon, p.l0)))
    return dfitness
end

"""
    function fitness(E::Float64, l::Int64, p::fermi_fitness)

Evaluate fermi fitness function.
"""
function _fitness(E::Real, l::Int64, p::fermi_fitness)::Float64
    fitness = (p.f0 * (1/(1 + exp(p.beta * (E - p.E_Star(l, p.epsilon, p.l0))))) - p.fl * l)
    return fitness
end

"""
    function fitness(E::Float64, p::fermi_fitness)

Evaluate fermi fitness function.
"""
function _dE_fitness(E::Real, l::Int64, p::fermi_fitness)
    dfitness = -p.f0 * 1/(1 + exp(p.beta * (E - p.E_Star(l, p.epsilon, p.l0))))^2 * p.beta * exp(p.beta * (E - p.E_Star(l, p.epsilon, p.l0)))
    return dfitness
end

"""
    Est(l)

Threshold of the fermi landscape.
"""
function Est(l::Int, epsilon::Real, l0::Int)::Real
    return (3l/4 - l0) * epsilon
end

mutable struct quadratic_fitness <: fitness_functions
    l::Int64
    c::Real
    l0::Int64
    epsilon::Real
    E_Star::Function
end

# Setting some defaults
quadratic_fitness(;l=10, c=1, l0=10, epsilon=1, E_Star=Est) = quadratic_fitness(l, c, l0, epsilon, E_Star)

"""
    function fitness(E::Float64, p::quadratic_fitness)

Evaluate quadratic fitness function.
"""
function _fitness(E::Float64, p::quadratic_fitness)
    fitness = -p.c * (E - p.E_Star(p.l, p.epsilon, p.l0))^2
    return fitness
end

"""
    function fitness(E::Float64, l::Int64, p::quadratic_fitness)

Evaluate fermi fitness function.
"""
function _fitness(E::Float64, l::Int64, p::quadratic_fitness)
    fitness = -p.c * (E - p.E_Star(l, p.epsilon, p.l0))^2
    return fitness
end

"""
    mutable struct linear_fitness <: fitness_functions

Fitness landscape that is linear.
"""
mutable struct linear_fitness <: Jevo.fitness_functions
    l::Int64
    s::Real
    fl::Real
    n::Int64
    epsilon::Real
    l0::Int64
end

# Setting some defaults
linear_fitness(;l=10, s=1, fl=0, n=4, epsilon=1, l0=10) = linear_fitness(l, s, fl, n, epsilon, l0)

function _fitness(E::Real, p::linear_fitness)::Float64
    fitness = -p.fl * p.l + p.s/(p.l * p.epsilon) * ((p.n-1)/p.n * (p.l * p.epsilon)  - p.l0  - E )
    return fitness
end

function _fitness(E::Real, l::Int64, p::linear_fitness)::Float64
    fitness = -p.fl * l + p.s/(p.l * p.epsilon) * ((p.n-1)/p.n * (p.l * p.epsilon)  - p.l0 - E )
    return fitness
end

#=
"""
    mutable struct linear_fitness <: fitness_functions

Fitness landscape that has the form of a fermi function.
"""
mutable struct semi_linear_fitness <: Jevo.fitness_functions
    l::Int64
    f0::Real
    fl::Float64
    n::Int64
    epsilon::Float64
end

# Setting some defaults
semi_linear_fitness(;l=10, f0=1, fl=0, n=4, epsilon=1) = semi_linear_fitness(l, f0, fl, n, epsilon)

function _fitness(E::Real, p::semi_linear_fitness)::Float64
    if E < p.epsilon * ((p.n-1)/p.n *p.l  - 5)
        fitness = -p.fl * p.l + p.f0
    else
        fitness = -p.fl * p.l - p.f0 / (5p.epsilon) * (E - p.l * p.epsilon * (p.n - 1)/p.n)
    end
    return fitness
end

function _fitness(E::Real, l::Int64, p::semi_linear_fitness)::Float64
    if E < p.epsilon * ((p.n-1)/p.n * l  - 5)
        fitness = -p.fl * l + p.f0
    else
        fitness = -p.fl * l - p.f0 / (5p.epsilon) * (E - l * p.epsilon * (p.n - 1)/p.n)
    end
    return fitness
end
=#

# Add some beautiful broadcasting
fitness(E, l, p) = _fitness(E, l, p)
Broadcast.broadcasted(::typeof(fitness), E, l, p) = broadcast(_fitness, E, l, Ref(p))

fitness(E, p) = _fitness(E, p)
Broadcast.broadcasted(::typeof(fitness), E, p) = broadcast(_fitness, E, Ref(p))

dE_fitness(E, l, p) = _dE_fitness(E, l, p)
Broadcast.broadcasted(::typeof(dE_fitness), E, l, p) = broadcast(_dE_fitness, E, l, Ref(p))

dE_fitness(E, p) = _dE_fitness(E, p)
Broadcast.broadcasted(::typeof(dE_fitness), E, p) = broadcast(_dE_fitness, E, Ref(p))

mutable struct num_fermi <: Jevo.fitness_functions
    n::Real
    l_0::Real
    gap::Real
    f0::Real
    fl::Real
end

num_fermi(;n=4, l_0=10, gap=7, f0=400, fl=10) = num_fermi(n, l_0, gap, f0, fl)

# Neutral Expectation
γ_0(n) = (n-1)/n

# Binding Threshold
γ_1(l, n, l_0) = γ_0(n) - l_0/l

# Binding probability
pb(γ, l, n, l_0, gap) = 1 / (1 + exp(gap * (l / l_0) * (γ - γ_1(l, n, l_0))))

# Fitness component for functional binding
F_b(γ, l, n, l_0, gap, f0) = f0 * pb(γ, l, n, l_0, gap)

# Fitness component of genomic constraint
F_c(l, l_0, fl) = - fl * l / l_0

# Total Fitness
F(γ, l, n, l_0, gap, f0, fl) = F_b(γ, l, n, l_0, gap, f0) + F_c(l, l_0, fl)

Jevo._fitness(γ::Real, l, p::num_fermi)::Float64 = F(γ/(l * p.gap/p.l_0), l, p.n, p.l_0, p.gap, p.f0, p.fl)

mutable struct const_fitness <: Jevo.fitness_functions
    s::Real
end

const_fitness(;s=1) = const_fitness(s)

function _fitness(p::const_fitness)::Real
    return p.s
end

function _fitness(E::Real, p::const_fitness)::Real
    return p.s
end

function _fitness(E::Real, l::Int64, p::const_fitness)::Real
    return p.s
end

fitness(p::const_fitness) = _fitness(p::const_fitness)


text:
##########################################################################################
# Basic imports and dimension definitions
##########################################################################################
using DimensionalData
using DimensionalData.Dimensions, DimensionalData.LookupArrays
using DimensionalData: basetypeof, broadcast_dims
using DimensionalData.Dimensions: setdims
using Dates

Time = DimensionalData.Ti

AbDimArray = DimensionalData.AbstractDimArray

"""
    dimindex(A::ClimArray, d) → i
Get the index of dimension `d`.
"""
dimindex(A, i::Int) = i
dimindex(A, dim) = DimensionalData.dimnum(A, dim)

"""
    gnv(object) → x
Short for "get numeric value", this function will return the pure numeric value(s)
of the given object. Convenience function for quickly getting the numeric data of
dimensional arrays or dimensions.
"""
gnv(x) = x
gnv(x::Union{AbDimArray, LookupArray}) = parent(x)
gnv(x::Dimension) = parent(parent(x))

export At, (..), Between, Near # Selectors from DimensionalArrays.jl
export hasdim, dims, dimindex
export Time, Lon, Lat, dims, Coord, Hei, Pre, Ti
export CoordinateSpace, OrthogonalSpace, spacestructure
export DimensionalData # for accessing its functions
export setdims
export gnv

@dim Lon IndependentDim "Longitude"
@dim Lat IndependentDim "Latitude"
@dim Hei IndependentDim "Height"
@dim Pre IndependentDim "Pressure"

STANDARD_DIMS = (Lon, Lat, Time, Hei, Pre, Coord)

"""
    COMMONNAMES
A dictionary of common names of dimensions (as strings) to actual dimension types.
"""
const COMMONNAMES = Dict(
    "lat" => Lat,
    "latitude" => Lat,
    "y" => Lat,
    "yc" => Lat,
    "lon" => Lon,
    "long" => Lon,
    "x" => Lon,
    "xc" => Lon,
    "longitude" => Lon,
    "time" => Time,
    "t" => Time,
    "age" => Time,
    "height" => Hei,
    "altitude" => Hei,
    "pressure" => Pre,
    "level" => Pre,
)

##########################################################################################
# Space types
##########################################################################################

# the following traits for the the way space is configured. currently the options are
# CoordinateSpace, which is for points/coordinates
# while the OrthogonalSpace is for standard Longitude x Latitude dimensions.
# Dispatch on `CoordinateSpace` or other types is type-stable
# because it comes from the underlying dimension

abstract type SpaceType end

"""
Space information is represented by two orthogonal dimensions `Lon, Lat`,
one being longitude and the other being latitude.
"""
struct OrthogonalSpace <: SpaceType end

"""
Space information is represented by a single dimension `Coord`, whose
elements are coordinates, i.e. 2-element `SVector(longitude, latitude)`.
Each coordinate represents the center of an arbitrary polygon in space.
The actual limits of each polygon are not included in the dimension for performance reasons.

In statistical functions such as [`spaceagg`](@ref), it is assumed that entry of the
coordinates covers **an equal amount of area**. If this is not the case, you can simply
provide an additional weights vector which would correspond to the area covered.

This dimension also allows indexing by latitude ranges, e.g. you can do
```julia
A # some `ClimArray` with a `Coord` dimension
A[Coord(Lat(-30..30)))]
```

Most functions of ClimateBase.jl implicitly assume that the coordinates are
sorted by latidude. You can achieve this with the following code:
```julia
A # some `ClimArray` with a `Coord` dimension
coords = gnv(dims(A, Coord))
si = sortperm(coords; by = reverse)
A = A[Coord(si)]
```
**This is done automatically by [`ncread`](@ref).**
"""
struct CoordinateSpace <: SpaceType end

spacestructure(a::AbDimArray) = spacestructure(dims(a))
function spacestructure(dims)
    if hasdim(dims, Coord)
        CoordinateSpace()
    elseif hasdim(dims, Lon) || hasdim(dims, Lat)
        OrthogonalSpace()
    else
        error("Array does not have any spatial dimensions: `Lon`, `Lat`, or `Coord`.")
    end
end

##########################################################################################
# ClimArray definition and DimensionalData.jl extensions
##########################################################################################
# TODO: remove `refdims` field
export ClimArray
struct ClimArray{T,N,D<:Tuple,R<:Tuple,A<:AbstractArray{T,N},Me} <: AbDimArray{T,N,D,A}
    data::A
    dims::D
    refdims::R
    name::Symbol
    attrib::Me
end
ClimArray(A::AbDimArray) = ClimArray(A.data, A.dims, A.refdims, A.name, A.metadata)
ClimArray(A::ClimArray, dims::Tuple = A.dims; name = A.name, attrib = A.attrib) =
ClimArray(A.data, format(dims, A.data), A.refdims, Symbol(name), attrib)

"""
    ClimArray(A::Array, dims::Tuple; name = "", attrib = nothing)

`ClimArray` is a structure that contains numerical array data bundled with dimensional
information, a name and an `attrib` field (typically a dictionary) that holds general
attributes.
You can think of `ClimArray` as a in-memory representation of a CFVariable.

At the moment, a `ClimArray` is using `DimArray` from DimensionalData.jl, and
all basic handling of `ClimArray` is offered by `DimensionalData` (see below).

`ClimArray` is created by passing in standard array data `A` and a
tuple of dimensions `dims`. See [`ncread`](@ref) to automatically create a `ClimArray`
from a .nc file. For obtaining the raw numeric values of a `ClimArray` or any of its
dimensions, use the function [`gnv`](@ref).

## Example
```julia
using ClimateBase, Dates
Time = ClimateBase.Ti # more intuitive name for time dimension
lats = -90:5:90
lons = 0:10:359
t = Date(2000, 3, 15):Month(1):Date(2020, 3, 15)
# dimensional information:
dimensions = (Lon(lons), Lat(lats), Time(t))
data = rand(36, 37, 241) # numeric data
A = ClimArray(data, dimensions)
```
"""
ClimArray(A::AbstractArray, dims::Tuple; refdims=(), name="", attrib=nothing) =
    ClimArray(A, format(dims, A), refdims, Symbol(name), attrib)
ClimArray(A::AbstractArray, dims::Tuple, name; refdims=(), attrib=nothing) =
    ClimArray(A, format(dims, A), refdims, Symbol(name), attrib)

Base.parent(A::ClimArray) = A.data
Base.@propagate_inbounds Base.setindex!(A::ClimArray, x, I::Vararg{DimensionalData.StandardIndices}) =
    setindex!(A.data, x, I...)

DimensionalData.metadata(A::ClimArray) = A.attrib
DimensionalData.basetypeof(::ClimArray) = ClimArray

DimensionalData.rebuild(
    ::ClimArray, data, dims, refdims, name, metadata
) = ClimArray(data, format(dims, data), refdims, name, metadata)
DimensionalData.rebuild(
    A::ClimArray;
    data=gnv(A), dims=dims(A), refdims=refdims(A), name=name(A), metadata=metadata(A)
) = ClimArray(data, format(dims, data), refdims, name, metadata)

# The following basic methods allow indexing with tuples, (Time(5), Lon(3))
Base.getindex(A::ClimArray, i::Tuple) = A[i...]
Base.setindex!(A::ClimArray, x, i::Tuple) = setindex!(A, x, i...)
Base.view(A::ClimArray, i::Tuple) = view(A, i...)

# Convenience
Base.ones(A::AbDimArray) = basetypeof(A)(ones(eltype(A), size(A)), dims(A))
Base.zeros(A::AbDimArray) = basetypeof(A)(zeros(eltype(A), size(A)), dims(A))


text:
using Photon
import Knet

# This file contains an example how to train a network on the MNIST dataset and
# plot some metrics afterwards (if you run this example in Juno or a
# Jupyter notebook).
#
# The MNIST dataset contains small (28x28) gray scale images that depict a
# handwritten digit. So the goal is to predict what the digit is.

# Load the MNIST data set
include(Knet.dir("data", "mnist.jl"))
trndata, tstdata = mnistdata()

# Create a simple Convolutional network.
model = Sequential(
      Conv2D(16, 3, relu),
      Conv2D(16, 3, relu),
      MaxPool2D(),
      Dense(64, relu),
      Dense(10)
)

# Create a workout containing the model, a loss function and
# the optimizer (default SGD) and register the accuracy.
workout = Workout(model, nll, acc=OneHotBinaryAccuracy())

# Run the training for 10 epochs and we don't need a convertor since
# mnistdata function already does that work.
fit!(workout, trndata, tstdata; epochs=10, convertor=identity)

println("\nTrained the model in $(workout.epochs) epochs.")

# Now let's plot some results. If you haven't installed Plots yet, you'll
# need to run:  using Pkg; Pkg.add("Plots")
import Plots
plotmetrics(Plots, workout, [:loss, :val_loss, :val_acc])


text:
using FractionalSystems, Plots

function RF(du, u, p, t)
    du[1] = u[2]*(u[3]-1+u[1]*u[1])+0.1*u[1]
    du[2] = u[1]*(3*u[3]+1-u[1]*u[1])+0.1*u[2]
    du[3] = -2*u[3]*(0.98+u[1]*u[2])
end

LE, tspan=FOLyapunov(RF, 0.999, 0, 0.02, 300, [0.1; 0.1; 0.1], 0.005, 1000)

plot(tspan, LE[1, :])
plot!(tspan, LE[2, :])
plot!(tspan, LE[3, :])


text:
"""
    conv_bn(kernelsize, inplanes, outplanes, activation = relu;
            rev = false, preact = true,
            stride = 1, pad = 0, dilation = 1, groups = 1, [bias, weight, init],
            initβ = Flux.zeros32, initγ = Flux.ones32, ϵ = 1f-5, momentum = 1f-1)

Create a convolution + batch normalization pair with activation.

# Arguments

  - `kernelsize`: size of the convolution kernel (tuple)
  - `inplanes`: number of input feature maps
  - `outplanes`: number of output feature maps
  - `activation`: the activation function for the final layer
  - `rev`: set to `true` to place the batch norm before the convolution
  - `preact`: set to `true` to place the activation function before the batch norm
    (only compatible with `rev = false`)
  - `stride`: stride of the convolution kernel
  - `pad`: padding of the convolution kernel
  - `dilation`: dilation of the convolution kernel
  - `groups`: groups for the convolution kernel
  - `bias`, `weight`, `init`: initialization for the convolution kernel (see [`Flux.Conv`](#))
  - `initβ`, `initγ`: initialization for the batch norm (see [`Flux.BatchNorm`](#))
  - `ϵ`, `momentum`: batch norm parameters (see [`Flux.BatchNorm`](#))
"""
function conv_bn(kernelsize, inplanes, outplanes, activation = relu;
                 rev = false, preact = false,
                 initβ = Flux.zeros32, initγ = Flux.ones32, ϵ = 1.0f-5, momentum = 1.0f-1,
                 kwargs...)
    layers = []
    if rev
        activations = (conv = activation, bn = identity)
        bnplanes = inplanes
    else
        activations = (conv = identity, bn = activation)
        bnplanes = outplanes
    end
    if preact
        rev ? throw(ArgumentError("preact and rev cannot be set at the same time")) :
        activations = (conv = activation, bn = identity)
    end
    push!(layers,
          Conv(kernelsize, Int(inplanes) => Int(outplanes), activations.conv; kwargs...))
    push!(layers,
          BatchNorm(Int(bnplanes), activations.bn;
                    initβ = initβ, initγ = initγ, ϵ = ϵ, momentum = momentum))
    return rev ? reverse(layers) : layers
end

"""
    depthwise_sep_conv_bn(kernelsize, inplanes, outplanes, activation = relu;
                          rev = false,
                          stride = 1, pad = 0, dilation = 1, [bias, weight, init],
                          initβ = Flux.zeros32, initγ = Flux.ones32,
                          ϵ = 1f-5, momentum = 1f-1)

Create a depthwise separable convolution chain as used in MobileNet v1.
This is sequence of layers:

  - a `kernelsize` depthwise convolution from `inplanes => inplanes`
  - a batch norm layer + `activation`
  - a `kernelsize` convolution from `inplanes => outplanes`
  - a batch norm layer + `activation`

See Fig. 3 in [reference](https://arxiv.org/abs/1704.04861v1).

# Arguments

  - `kernelsize`: size of the convolution kernel (tuple)
  - `inplanes`: number of input feature maps
  - `outplanes`: number of output feature maps
  - `activation`: the activation function for the final layer
  - `rev`: set to `true` to place the batch norm before the convolution
  - `stride`: stride of the first convolution kernel
  - `pad`: padding of the first convolution kernel
  - `dilation`: dilation of the first convolution kernel
  - `bias`, `weight`, `init`: initialization for the convolution kernel (see [`Flux.Conv`](#))
  - `initβ`, `initγ`: initialization for the batch norm (see [`Flux.BatchNorm`](#))
  - `ϵ`, `momentum`: batch norm parameters (see [`Flux.BatchNorm`](#))
"""
function depthwise_sep_conv_bn(kernelsize, inplanes, outplanes, activation = relu;
                               rev = false,
                               initβ = Flux.zeros32, initγ = Flux.ones32,
                               ϵ = 1.0f-5, momentum = 1.0f-1,
                               stride = 1, kwargs...)
    return vcat(conv_bn(kernelsize, inplanes, inplanes, activation;
                        rev = rev, initβ = initβ, initγ = initγ,
                        ϵ = ϵ, momentum = momentum,
                        stride = stride, groups = Int(inplanes), kwargs...),
                conv_bn((1, 1), inplanes, outplanes, activation;
                        rev = rev, initβ = initβ, initγ = initγ,
                        ϵ = ϵ, momentum = momentum))
end

"""
    skip_projection(inplanes, outplanes, downsample = false)

Create a skip projection
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:

  - `inplanes`: the number of input feature maps
  - `outplanes`: the number of output feature maps
  - `downsample`: set to `true` to downsample the input
"""
function skip_projection(inplanes, outplanes, downsample = false)
    return downsample ?
           Chain(conv_bn((1, 1), inplanes, outplanes, identity; stride = 2, bias = false)) :
           Chain(conv_bn((1, 1), inplanes, outplanes, identity; stride = 1, bias = false))
end

# array -> PaddedView(0, array, outplanes) for zero padding arrays
"""
    skip_identity(inplanes, outplanes[, downsample])

Create a identity projection
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:

  - `inplanes`: the number of input feature maps
  - `outplanes`: the number of output feature maps
  - `downsample`: this argument is ignored but it is needed for compatibility with [`resnet`](#).
"""
function skip_identity(inplanes, outplanes)
    if outplanes > inplanes
        return Chain(MaxPool((1, 1); stride = 2),
                     y -> cat(y,
                              zeros(eltype(y),
                                    size(y, 1),
                                    size(y, 2),
                                    outplanes - inplanes, size(y, 4)); dims = 3))
    else
        return identity
    end
end
skip_identity(inplanes, outplanes, downsample) = skip_identity(inplanes, outplanes)

"""
    squeeze_excite(channels, reduction = 4)

Squeeze and excitation layer used by MobileNet variants
([reference](https://arxiv.org/abs/1905.02244)).

# Arguments

  - `channels`: the number of input/output feature maps
  - `reduction = 4`: the reduction factor for the number of hidden feature maps
    (must be >= 1)
"""
function squeeze_excite(channels, reduction = 4)
    @assert (reduction>=1) "`reduction` must be >= 1"
    return SkipConnection(Chain(AdaptiveMeanPool((1, 1)),
                                conv_bn((1, 1), channels, channels ÷ reduction, relu;
                                        bias = false)...,
                                conv_bn((1, 1), channels ÷ reduction, channels, hardσ)...),
                          .*)
end

"""
    invertedresidual(kernel_size, inplanes, hidden_planes, outplanes, activation = relu;
                     stride, reduction = nothing)

Create a basic inverted residual block for MobileNet variants
([reference](https://arxiv.org/abs/1905.02244)).

# Arguments

  - `kernel_size`: The kernel size of the convolutional layers
  - `inplanes`: The number of input feature maps
  - `hidden_planes`: The number of feature maps in the hidden layer
  - `outplanes`: The number of output feature maps
  - `activation`: The activation function for the first two convolution layer
  - `stride`: The stride of the convolutional kernel, has to be either 1 or 2
  - `reduction`: The reduction factor for the number of hidden feature maps
    in a squeeze and excite layer (see [`squeeze_excite`](#)).
    Must be >= 1 or `nothing` for no squeeze and excite layer.
"""
function invertedresidual(kernel_size, inplanes, hidden_planes, outplanes,
                          activation = relu; stride, reduction = nothing)
    @assert stride in [1, 2] "`stride` has to be 1 or 2"
    pad = @. (kernel_size - 1) ÷ 2
    conv1 = (inplanes == hidden_planes) ? identity :
            Chain(conv_bn((1, 1), inplanes, hidden_planes, activation; bias = false))
    selayer = isnothing(reduction) ? identity : squeeze_excite(hidden_planes, reduction)
    invres = Chain(conv1,
                   conv_bn(kernel_size, hidden_planes, hidden_planes, activation;
                           bias = false, stride, pad = pad, groups = hidden_planes)...,
                   selayer,
                   conv_bn((1, 1), hidden_planes, outplanes, identity; bias = false)...)
    return (stride == 1 && inplanes == outplanes) ? SkipConnection(invres, +) : invres
end

function invertedresidual(kernel_size::Integer, args...; kwargs...)
    return invertedresidual((kernel_size, kernel_size), args...; kwargs...)
end


text:
using BeastUtils.Simulation, BeastUtils.TreeUtils, BeastUtils.SlowLikelihoods,
      BeastUtils.MatrixUtils
using LinearAlgebra, Test, PhyloNetworks

import Random
seed =666
Random.seed!(seed)

const var_diff_tol = 1e-10

function check_simulation(tsm::Simulation.TraitSimulationModel,
                            μ_true::Vector{Float64},
                            V_true::Matrix{Float64};
                            reps::Int = 100_000,
                            μ_tol::Float64 = 1e-1,
                            v_tol::Float64 = 1e-1)

    V = zeros(n * p, n * p)
    μ = zeros(n * p)

    for i = 1:reps
        y = vec(Simulation.simulate(tsm))
        μ .+= y
        V .+= y * y'
    end

    μ ./= reps

    V ./= reps
    V = V - μ * μ'

    V_diff = abs.(V - V_true)
    μ_diff = abs.(μ - μ_true)

    v_max = maximum(V_diff)
    μ_max = maximum(μ_diff)

    @show v_max
    @show μ_max

    return v_max < v_tol && μ_max < μ_tol
end

n = 5
p = 3
taxa = ["taxon$i" for i = 1:n]

## Tree only

# Known tree
newick = "((A:1.0,(B:0.5,C:0.8):0.1):2.0,(D:0.5,E:2.0):1.5);"
tree = readTopology(newick)
Σ = randpd(p, rescale=true)
Ψ1 = [3.0 2.0 2.0 0.0 0.0;
      2.0 2.6 2.1 0.0 0.0;
      2.0 2.1 2.9 0.0 0.0;
      0.0 0.0 0.0 2.0 1.5;
      0.0 0.0 0.0 1.5 3.5]

taxa1 = ["A", "B", "C", "D", "E"]
μ_true = zeros(n * p)
V_true = kron(Σ, Ψ1)
tdm = Simulation.TreeDiffusionModel(tree, Σ)
tsm = Simulation.TraitSimulationModel(taxa1, tdm)

model_params = SlowLikelihoods.sim_to_params(tsm, pss=Inf)
V_true2 = var(model_params, tree, taxa1)

@test maximum(abs.(V_true - V_true2)) < var_diff_tol

@test check_simulation(tsm, μ_true, V_true)

perm = [3, 2, 5, 4, 1] # scramble order
taxa2 = taxa1[perm]
Ψ2 = Ψ1[perm, perm]

μ_true = zeros(n * p)
V_true = kron(Σ, Ψ2)
tdm = Simulation.TreeDiffusionModel(tree, Σ)
tsm = Simulation.TraitSimulationModel(taxa2, tdm)

model_params = SlowLikelihoods.sim_to_params(tsm, pss=Inf)
V_true2 = var(model_params, tree, taxa2)
@test maximum(abs.(V_true - V_true2)) < var_diff_tol

@test check_simulation(tsm, μ_true, V_true)

# Random Tree
tree = TreeUtils.rtree(n, labels = taxa, keep_order = false, ultrametric = true)
standardize_height!(tree)
Σ = randpd(p, rescale=true)
tdm = Simulation.TreeDiffusionModel(tree, Σ)
tsm = Simulation.TraitSimulationModel(taxa, tdm)

μ_true = zeros(n * p)
Ψ = TreeUtils.vcv(tree, taxa)
V_true = kron(Σ, Ψ)

model_params = SlowLikelihoods.sim_to_params(tsm, pss=Inf)
V_true2 = var(model_params, tree, taxa)
@test maximum(abs.(V_true - V_true2)) < var_diff_tol

@test check_simulation(tsm, μ_true, V_true)

## Residual variance

tree = TreeUtils.rtree(n, labels = taxa, keep_order = false, ultrametric = false)

Σ = randn(p, p)
Γ = randn(p, p)

Σ = Σ * Σ' # needs to be positive definite
Γ = Γ * Γ' # needs to be positive definite

tdm = Simulation.TreeDiffusionModel(tree, Σ)
rvm = Simulation.ResidualVarianceModel(Γ)
tsm = Simulation.TraitSimulationModel(taxa, tdm, rvm)

μ_true = zeros(n * p)
Ψ = TreeUtils.vcv(tdm.tree, taxa)
V_true = kron(Σ, Ψ) + kron(Γ, Diagonal(ones(n)))

model_params = SlowLikelihoods.sim_to_params(tsm, pss=Inf)
V_true2 = var(model_params, tree, taxa)
@test maximum(abs.(V_true - V_true2)) < var_diff_tol

@test check_simulation(tsm, μ_true, V_true)

## Latent factor model

k = 2

tree = TreeUtils.rtree(n, labels = taxa, keep_order = false, ultrametric = false)

Σ = Diagonal(ones(k))
Λ = Diagonal(rand(p))
L = randn(k, p)

tdm = Simulation.TreeDiffusionModel(tree, Σ)
lfm = Simulation.LatentFactorModel(L, Λ)
tsm = Simulation.TraitSimulationModel(taxa, tdm, lfm)

μ_true = zeros(n * p)
Ψ = TreeUtils.vcv(tdm.tree, taxa)
V_true = kron(L' * L, Ψ) + kron(Λ, Diagonal(ones(n)))

@test check_simulation(tsm, μ_true, V_true)


text:
# This file was generated, do not modify it. # hide
XGBR = @load XGBoostRegressor

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
# ------------------------------------------------------------------------------
#                              VML SERVICE FUNCTION
# ------------------------------------------------------------------------------

if debug println("MKL : VML SERVICE") end
for fname in (  :vml_get_mode, :vml_set_mode,
                :vml_get_accuracy, :vml_set_accuracy,
                :vml_check_error)

                if debug println("export : ", string(fname)) end
                @eval export $fname
end

# ------------------------------------------------------------------------------

# credit : Simon Kornblith (VML.jl)
immutable VMLAccuracy
    mode::UInt
end

const VML_LA = VMLAccuracy(0x00000001)
const VML_HA = VMLAccuracy(0x00000002)
const VML_EP = VMLAccuracy(0x00000003)

Base.show(io::IO, m::VMLAccuracy) = print(io, m == VML_LA ? "VML_LA" :
                                              m == VML_HA ? "VML_HA" : "VML_EP")

vml_get_mode() = ccall((:_vmlGetMode, libvml), Cuint, ())
vml_set_mode(mode::Integer) = (ccall((:_vmlSetMode, libvml), Cuint, (UInt,), mode); nothing)

vml_set_accuracy(m::VMLAccuracy) = vml_set_mode((vml_get_mode() & ~0x03) | m.mode)
vml_get_accuracy() = VMLAccuracy(vml_get_mode() & 0x3)

vml_set_mode((vml_get_mode() & ~0x0000FF00))

function vml_check_error()
    vml_error = ccall((:_vmlClearErrStatus, libvml), Cint, ())
    if vml_error != 0
        if vml_error == 1
            error(DomainError())
        elseif vml_error == 2 || vml_error == 3 || vml_error == 4
            # Singularity, overflow, or underflow
            # I don't think Base throws on these
        elseif vml_error == 1000
            warn("VML does not support $(vml_get_accuracy); lower accuracy used instead")
        else
            error("an unexpected error occurred in VML ($vml_error)")
        end
    end
end

#=
/* file: mkl_vml_defines.h */
/*
//  Copyright(C) 2006-2011 Intel Corporation. All Rights Reserved.
//
//  The source code, information  and  material ("Material") contained herein is
//  owned  by Intel Corporation or its suppliers or licensors, and title to such
//  Material remains  with Intel Corporation  or its suppliers or licensors. The
//  Material  contains proprietary information  of  Intel or  its  suppliers and
//  licensors. The  Material is protected by worldwide copyright laws and treaty
//  provisions. No  part  of  the  Material  may  be  used,  copied, reproduced,
//  modified, published, uploaded, posted, transmitted, distributed or disclosed
//  in any way  without Intel's  prior  express written  permission. No  license
//  under  any patent, copyright  or  other intellectual property rights  in the
//  Material  is  granted  to  or  conferred  upon  you,  either  expressly,  by
//  implication, inducement,  estoppel or  otherwise.  Any  license  under  such
//  intellectual  property  rights must  be express  and  approved  by  Intel in
//  writing.
//
//  *Third Party trademarks are the property of their respective owners.
//
//  Unless otherwise  agreed  by Intel  in writing, you may not remove  or alter
//  this  notice or  any other notice embedded  in Materials by Intel or Intel's
//  suppliers or licensors in any way.
*/
/*
//++
//  Macro definitions visible on user level.
//--
*/

#ifndef __MKL_VML_DEFINES_H__
#define __MKL_VML_DEFINES_H__

#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */

/*
//++
//  MACRO DEFINITIONS
//  Macro definitions for VML mode and VML error status.
//
//  VML mode controls VML function accuracy, floating-point settings (rounding
//  mode and precision) and VML error handling options. Default VML mode is
//  VML_HA | VML_ERRMODE_DEFAULT, i.e. VML high accuracy functions are
//  called, and current floating-point precision and the rounding mode is used.
//
//  Error status macros are used for error classification.
//--
*/

/*
//  VML FUNCTION ACCURACY CONTROL
//  VML_HA - when VML_HA is set, high accuracy VML functions are called
//  VML_LA - when VML_LA is set, low accuracy VML functions are called
//  VML_EP - when VML_EP is set, enhanced performance VML functions are called
//
//  NOTE: VML_HA, VML_LA and VML_EP must not be used in combination
*/
#define VML_LA 0x00000001
#define VML_HA 0x00000002
#define VML_EP 0x00000003

/*
//  SETTING OPTIMAL FLOATING-POINT PRECISION AND ROUNDING MODE
//  Definitions below are to set optimal floating-point control word
//  (precision and rounding mode).
//
//  For their correct work, VML functions change floating-point precision and
//  rounding mode (if necessary). Since control word changing is typically
//  expensive operation, it is recommended to set precision and rounding mode
//  to optimal values before VML function calls.
//
//  VML_FLOAT_CONSISTENT  - use this value if the calls are typically to single
//                          precision VML functions
//  VML_DOUBLE_CONSISTENT - use this value if the calls are typically to double
//                          precision VML functions
//  VML_RESTORE           - restore original floating-point precision and
//                          rounding mode
//  VML_DEFAULT_PRECISION - use default (current) floating-point precision and
//                          rounding mode
//  NOTE: VML_FLOAT_CONSISTENT, VML_DOUBLE_CONSISTENT, VML_RESTORE and
//        VML_DEFAULT_PRECISION must not be used in combination
*/
#define VML_DEFAULT_PRECISION 0x00000000
#define VML_FLOAT_CONSISTENT  0x00000010
#define VML_DOUBLE_CONSISTENT 0x00000020
#define VML_RESTORE           0x00000030

/*
//  VML ERROR HANDLING CONTROL
//  Macros below are used to control VML error handler.
//
//  VML_ERRMODE_IGNORE   - ignore errors
//  VML_ERRMODE_ERRNO    - errno variable is set on error
//  VML_ERRMODE_STDERR   - error description text is written to stderr on error
//  VML_ERRMODE_EXCEPT   - exception is raised on error
//  VML_ERRMODE_CALLBACK - user's error handler function is called on error
//  VML_ERRMODE_DEFAULT  - errno variable is set, exceptions are raised and
//                         user's error handler is called on error
//  NOTE: VML_ERRMODE_IGNORE must not be used in combination with
//        VML_ERRMODE_ERRNO, VML_ERRMODE_STDERR, VML_ERRMODE_EXCEPT,
//        VML_ERRMODE_CALLBACK and VML_ERRMODE_DEFAULT.
*/
#define VML_ERRMODE_IGNORE   0x00000100
#define VML_ERRMODE_ERRNO    0x00000200
#define VML_ERRMODE_STDERR   0x00000400
#define VML_ERRMODE_EXCEPT   0x00000800
#define VML_ERRMODE_CALLBACK 0x00001000
#define VML_ERRMODE_DEFAULT  \
VML_ERRMODE_ERRNO | VML_ERRMODE_CALLBACK | VML_ERRMODE_EXCEPT

/*
//  OpenMP(R) number of threads mode macros
//  VML_NUM_THREADS_OMP_AUTO   - Maximum number of threads is determined by
//                               environmental variable OPM_NUM_THREADS or
//                               omp_set_num_threads() function
//  VML_NUM_THREADS_OMP_FIXED  - Number of threads is determined by
//                               environmental variable OPM_NUM_THREADS
//                               omp_set_num_threads() functions
*/
#define VML_NUM_THREADS_OMP_AUTO   0x00000000
#define VML_NUM_THREADS_OMP_FIXED  0x00010000

/*
//  FTZ & DAZ mode macros
//  VML_FTZDAZ_ON   - FTZ & DAZ MXCSR mode enabled
//                    for faster (sub)denormal values processing
//  VML_FTZDAZ_OFF  - FTZ & DAZ MXCSR mode disabled
//                    for accurate (sub)denormal values processing
*/
#define VML_FTZDAZ_ON   0x00280000
#define VML_FTZDAZ_OFF  0x00140000

/*
//  ACCURACY, FLOATING-POINT CONTROL, FTZDAZ AND ERROR HANDLING MASKS
//  Accuracy, floating-point and error handling control are packed in
//  the VML mode variable. Macros below are useful to extract accuracy and/or
//  floating-point control and/or error handling control settings.
//
//  VML_ACCURACY_MASK           - extract accuracy bits
//  VML_FPUMODE_MASK            - extract floating-point control bits
//  VML_ERRMODE_MASK            - extract error handling control bits
//                                (including error callback bits)
//  VML_ERRMODE_STDHANDLER_MASK - extract error handling control bits
//                                (not including error callback bits)
//  VML_ERRMODE_CALLBACK_MASK   - extract error callback bits
//  VML_NUM_THREADS_OMP_MASK    - extract OpenMP(R) number of threads mode bits
//  VML_FTZDAZ_MASK             - extract FTZ & DAZ bits
*/
#define VML_ACCURACY_MASK           0x0000000F
#define VML_FPUMODE_MASK            0x000000F0
#define VML_ERRMODE_MASK            0x0000FF00
#define VML_ERRMODE_STDHANDLER_MASK 0x00000F00
#define VML_ERRMODE_CALLBACK_MASK   0x0000F000
#define VML_NUM_THREADS_OMP_MASK    0x00030000
#define VML_FTZDAZ_MASK             0x003C0000

/*
//  ERROR STATUS MACROS
//  VML_STATUS_OK        - no errors
//  VML_STATUS_BADSIZE   - array dimension is not positive
//  VML_STATUS_BADMEM    - invalid pointer passed
//  VML_STATUS_ERRDOM    - at least one of arguments is out of function domain
//  VML_STATUS_SING      - at least one of arguments caused singularity
//  VML_STATUS_OVERFLOW  - at least one of arguments caused overflow
//  VML_STATUS_UNDERFLOW - at least one of arguments caused underflow
//  VML_STATUS_ACCURACYWARNING - function doesn't support set accuracy mode,
//                               lower accuracy mode was used instead
*/
#define VML_STATUS_OK                    0
#define VML_STATUS_BADSIZE              -1
#define VML_STATUS_BADMEM               -2
#define VML_STATUS_ERRDOM                1
#define VML_STATUS_SING                  2
#define VML_STATUS_OVERFLOW              3
#define VML_STATUS_UNDERFLOW             4
#define VML_STATUS_ACCURACYWARNING       1000

#ifdef __cplusplus
}
#endif /* __cplusplus */

#endif /* __MKL_VML_DEFINES_H__ */
=#


text:
using Distributions, MCMC, ZVSimulator

function C(n::Int, c::Float64)
  X = eye(n)
  [(j <= n-i) ? X[i+j, i] = X[i, i+j] = c^j : nothing for i = 1:(n-1), j = 1:(n-1)]
  X
end

df, npars, a = 5., 4, 0.25
c = ((df-2)/df)*C(npars, a)

mcmodel = model(MvTDist(df, zeros(npars), c), init=zeros(npars))
samplers = [HMCDA() for i in 1:5]
tasks = mcmodel * samplers * SerialMC(steps=10000, burnin=1000)

results = psim_serialmc_vrf(tasks)

# writedlm_vrf(results, "./", ' ')
# writedlm_reformat("_means.txt", "./_means_pprint.txt", fill("%.3e", 3))