text:
""" Finds the linear min/max bounds of JuMP.VariableRefs."""
function find_linear_bounds!(gm::GlobalModel; bbls::Array{BlackBoxLearner} = gm.bbls, all_bounds::Bool = false)
    unbounds = get_unbounds(bbls)
    if all_bounds
        unbounds = get_bounds(bbls)
    end
    if isnothing(unbounds)
        return
    end
    clear_tree_constraints!(gm)
    m = gm.model
    orig_objective = JuMP.objective_function(m)
    new_bounds = copy(unbounds)
    @showprogress 0.5 "Finding bounds..." for var in collect(keys(unbounds))
        if isinf(unbounds[var][1]) || all_bounds
            @objective(m, Min, var);
            JuMP.optimize!(m);
            if termination_status(m) == MOI.OPTIMAL
                new_bounds[var][1] = getvalue(var);
            end
        end
        if isinf(unbounds[var][2]) || all_bounds
            @objective(m, Max, var);
            JuMP.optimize!(m);
            if termination_status(m) == MOI.OPTIMAL
                new_bounds[var][2] = getvalue(var);
            end
        end
    end
    # Revert objective, bounds and return new_bounds
    @objective(m, Min, orig_objective)
    bound!(gm, new_bounds)
    return new_bounds
end

"""
    find_bounds!(gm::GlobalModel; bbls::Array{BlackBoxLearner} = [], all_bounds::Bool = false)

Finds the outer variable bounds of GlobalModel by solving only over the linear constraints
and listed bbls.
TODO: improve! Only find bounds of non-binary variables.
"""
function find_bounds!(gm::GlobalModel; bbls::Array{BlackBoxLearner} = gm.bbls, all_bounds::Bool = false)
    linear_bounds = find_linear_bounds!(gm, bbls = bbls, all_bounds = all_bounds)
    return linear_bounds
end

"""
    ridge_regress(X::DataFrame, Y::Array, solver::MOI.OptimizerWithAttributes; rho::Float64 = 0., weights = ones(length(Y)))

Performs ridge regression on data.
"""
function ridge_regress(X::DataFrame, Y::Array, solver::MOI.OptimizerWithAttributes; rho::Float64 = 0., weights = ones(length(Y)))
    m = JuMP.Model(with_optimizer(solver))
    normalized_X, lbs, ubs = normalized_data(X);
    @variable(m, x[1:size(X,2)])
    @variable(m, offset)
    @objective(m, Min, sum((Y - normalized_X*x .- offset).^2) + rho*sum(x.^2))
    status = optimize!(m)
    return getvalue(offset), getvalue.(x)./(ubs-lbs)
end

"""
    u_regress(X::DataFrame, Y::Array, solver::MOI.OptimizerWithAttributes)

Finds upper regressors of data that are conservative.
"""
function u_regress(X::DataFrame, Y::Array, solver::MOI.OptimizerWithAttributes)
    if size(X, 1) < 2*size(X, 2)
        @warn("Upper regression doesn't have enough data, thus returning constant bounds. ")
        return maximum(Y), zeros(size(X,2))
    end
    m = JuMP.Model(with_optimizer(solver))
    @variable(m, α[1:size(X, 2)])
    @variable(m, α0)
    @variable(m, err[1:length(Y)] >= 0)
    @constraint(m, [i=1:length(Y)], sum(Array{Float64}(X[i, :]).* α) + α0 >= Y[i])
    @constraint(m, err .>= (Matrix(X) * α .+ α0 - Y))
    @constraint(m, err .>= -(Matrix(X) * α .+ α0 - Y))
    @objective(m, Min, sum(err.^2))
    try
        optimize!(m)
        return getvalue(α0), getvalue.(α)
    catch
        @warn("Infeasible u_regress, returning constant bounds.")
        return maximum(Y), zeros(size(X,2))
    end
end

"""
    l_regress(X::DataFrame, Y::Array, solver::MOI.OptimizerWithAttributes)

Finds lower regressors of data that are conservative.
"""
function l_regress(X::DataFrame, Y::Array, solver::MOI.OptimizerWithAttributes)
    if size(X, 1) < 2*size(X, 2)
        @warn("Lower regression doesn't have enough data, thus returning constant bounds. ")
        return minimum(Y), zeros(size(X,2))
    end
    m = JuMP.Model(with_optimizer(solver))
    @variable(m, β[1:size(X, 2)])
    @variable(m, β0)
    @variable(m, err[1:length(Y)] >= 0)
    @constraint(m, [i=1:length(Y)], sum(Array{Float64}(X[i, :]).* β) + β0 <= Y[i])
    @constraint(m, err .>= (Matrix(X) * β .+ β0 - Y))
    @constraint(m, err .>= -(Matrix(X) * β .+ β0 - Y))
    @objective(m, Min, sum(err.^2))
    try
        optimize!(m)
        return getvalue(β0), getvalue.(β)
    catch
        @warn("Infeasible l_regress, returning constant bounds.")
        return minimum(Y), zeros(size(X,2))
    end
end

"""
    svm(X::DataFrame, Y::Array, solver::MOI.OptimizerWithAttributes, threshold = 0)

Finds the unregularized SVM split, where threshold is the allowable error.
"""
function svm(X::Matrix, Y::Array, solver::MOI.OptimizerWithAttributes, threshold = 0)
    m = JuMP.Model(with_optimizer(solver))
    @variable(m, error[1:length(Y)] >= 0)
    @variable(m, β[1:size(X, 2)])
    @variable(m, β0)
    for i=1:length(Y)
        @constraint(m, threshold + error[i] >= Y[i] - β0 - sum(X[i,:] .* β))
        @constraint(m, threshold + error[i] >= -Y[i] + β0 + sum(X[i,:] .* β))
    end
    @objective(m, Min, sum(error))
    optimize!(m)
    return getvalue(β0), getvalue.(β)
end

"""
    reweight(X::Matrix, mag::Float64 = 10)

Gaussian reweighting of existing data by proximity to previous solution.
Note: mag -> Inf results in uniform reweighting.
Returns:
- weights: weights of X rows, by Euclidian distance
"""
function reweight(bbl::BlackBoxLearner, sol::DataFrame, mag::Float64 = 10.)
    n_samples, n_features = size(bbl.X);
    bounds = get_bounds(bbl.vars)
    bound_dist = [abs(bounds[var][1] - bounds[var][2]) for var in bbl.vars]
    vecsol = [sol[1,strvar] for strvar in string.(bbl.vars)]
    distances = []
    for i = eachrow(bbl.X)
        push!(distances, sum((values(i) .- vecsol).^2 ./ bound_dist))
    end
    weights = exp.(-1/mag*distances);
    return weights
end

"""
    add_infeasibility_cuts(gm::GlobalModel)

Adds cuts reducing infeasibility of BBC inequalities.
"""
function add_infeasibility_cuts!(gm::GlobalModel)
    #TODO: CHECK REGRESSION EQUALITIES.
    var_vals = solution(gm)
    cut_count = 0
    feas_tol = get_param(gm, :tighttol)
    for i=1:length(gm.bbls)
        # Inequality Classifiers
        if get_param(gm.bbls[i], :gradients) && gm.bbls[i] isa BlackBoxClassifier && !gm.bbls[i].equality
            bbc = gm.bbls[i]
            if bbc.feas_gap[end] <= -feas_tol
                leaf = bbc.active_leaves[1]
                @assert length(bbc.active_leaves) == 1
                rel_vals = var_vals[:, string.(bbc.vars)]
                eval!(bbc, rel_vals)
                Y = bbc.Y[end]
                update_gradients(bbc, [size(bbc.X, 1)])
                cut_grad = bbc.gradients[end, :]
                push!(bbc.mi_constraints[leaf],
                    @constraint(gm.model, sum(Array(cut_grad) .* (bbc.leaf_variables[leaf][2] .-
                                        (Array(rel_vals)' .* bbc.leaf_variables[leaf][1]))) +
                                        Y * bbc.leaf_variables[leaf][1] + bbc.relax_var >= 0))
                cut_count += 1
            end
            for ll in bbc.lls
                if ll.feas_gap[end] <= -feas_tol
                    leaf = ll.active_leaves[1]
                    @assert length(ll.active_leaves) == 1
                    rel_vals = DataFrame(Array(var_vals[:, string.(ll.vars)]), string.(bbc.vars))
                    eval!(bbc, rel_vals)
                    Y = bbc.Y[end]
                    update_gradients(bbc, [size(bbc.X, 1)])
                    cut_grad = bbc.gradients[end, :]
                    push!(ll.mi_constraints[leaf],
                    @constraint(gm.model, sum(Array(cut_grad) .* (ll.leaf_variables[leaf][2] .-
                                        (Array(rel_vals)' .* ll.leaf_variables[leaf][1]))) +
                                        Y * ll.leaf_variables[leaf][1] + ll.relax_var >= 0))
                    cut_count += 1
                end
            end
        # Convex Regressors
        elseif get_param(gm.bbls[i], :gradients) && gm.bbls[i] isa BlackBoxRegressor && gm.bbls[i].convex
            bbr = gm.bbls[i]
            if bbr.feas_gap[end] <= -feas_tol
                rel_vals = var_vals[:, string.(bbr.vars)]
                eval!(bbr, rel_vals)
                Y = bbr.Y[end]
                update_gradients(bbr, [size(bbr.X, 1)])
                cut_grad = bbr.gradients[end, :]
                push!(bbr.mi_constraints[1],
                    @constraint(gm.model, bbr.dependent_var + bbr.relax_var >=
                        sum(Array(cut_grad) .* (bbr.vars .- Array(rel_vals)')) + Y))
                    cut_count += 1
            end
            for ll in bbr.lls
                if ll.feas_gap[end] <= -feas_tol
                    rel_vals = DataFrame(Array(var_vals[:, string.(ll.vars)]), string.(bbr.vars))
                    eval!(bbr, rel_vals)
                    Y = bbr.Y[end]
                    update_gradients(bbr, [size(bbr.X, 1)])
                    cut_grad = bbr.gradients[end, :]
                    push!(ll.mi_constraints[1],
                        @constraint(gm.model, ll.dependent_var + ll.relax_var >=
                            sum(Array(cut_grad) .* (ll.vars .- Array(rel_vals)')) + Y))
                        cut_count += 1
                end
            end
        # elseif get_param(gm.bbls[i], :gradients) && gm.bbls[i] isa BlackBoxRegressor
        #     bbr = gm.bbls[i]
        #     if bbr.feas_gap[end] <= -feas_tol
        #         rel_vals = var_vals[:, string.(bbr.vars)]
        #         eval!(bbr, rel_vals)
        #         Y = bbr.Y[end]
        #         update_gradients(bbr, [size(bbr.X, 1)])
        #         cut_grad = bbr.gradients[end, :]
        #         push!(bbr.mi_constraints[1],
        #             @constraint(gm.model, bbr.dependent_var + bbr.relax_var >=
        #                 sum(Array(cut_grad) .* (bbr.vars .- Array(rel_vals)')) + Y))
        #             cut_count += 1
        #     end
        #     for ll in bbr.lls
        #         if ll.feas_gap[end] <= -feas_tol
        #             rel_vals = DataFrame(Array(var_vals[:, string.(ll.vars)]), string.(bbr.vars))
        #             eval!(bbr, rel_vals)
        #             Y = bbr.Y[end]
        #             update_gradients(bbr, [size(bbr.X, 1)])
        #             cut_grad = bbr.gradients[end, :]
        #             push!(ll.mi_constraints[1],
        #                 @constraint(gm.model, ll.dependent_var + ll.relax_var >=
        #                     sum(Array(cut_grad) .* (ll.vars .- Array(rel_vals)')) + Y))
        #                 cut_count += 1
        #         end
        #     end
        end
        # TODO: add infeasibility cuts for Regressors that are locally convex.
        # TODO: add infeasibility cuts for equalities as well.
    end
    @info "$(cut_count) infeasibility cuts added."
    return cut_count
end


text:
using Test

using CompScienceMeshes
using BEAST
using StaticArrays

ω = 1.0

l1 = meshsegment(1.0,1/2)
l2 = meshsegment(1.0,1/4)

idcs = BEAST.interior_and_junction_vertices(l1, boundary(l1))
@test length(idcs) == 3

lag1 = lagrangec0d1(l1, boundary(l1))
lag2 = lagrangecxd0(l2)

id = Identity()
G = assemble(id, lag1, lag2)
H = assemble(id, lag2, lag1)

Gt = [
    3 1 0 0
    1 3 3 1
    0 0 1 3] // 16

@test norm(G-H') ≈ 0.0  atol = sqrt(eps()) #G == H'
@test norm(G-Gt) ≈ 0.0  atol = sqrt(eps())

# Test whether localoperator and localoperator 2 give the same result
# for test and basis functions defined on the same mesh.
id = Identity()
nc = NCross()

#fn = Pkg.dir("BEAST","test","sphere2.in")
fn = joinpath(dirname(@__FILE__),"assets","sphere316.in")
m = readmesh(fn)
rt = raviartthomas(m)

G1 = zeros(ComplexF64, numfunctions(rt), numfunctions(rt))
BEAST.assemble_local_matched!(id, rt, rt, (v,m,n)->(G1[m,n]+=v))

G2 = zeros(ComplexF64, numfunctions(rt), numfunctions(rt))
BEAST.assemble_local_mixed!(id, rt, rt, (v,m,n)->(G2[m,n]+=v))

# Test wether the gram matrix assembled works for
# different but overlapping meshes
m1 = meshrectangle(1.0,1.0,1.0)
m2 = translate(m1, point(0.5,0.5,0.0))
I = Identity()
x1 = lagrangecxd0(m1)
x2 = lagrangecxd0(m2)
G = assemble(I,x1,x2)
@test sum(G) ≈ 1/4

ext: jl


text:
using Random
function parameter_selection()
        # maximum parameter

        max_Interceptioncapacity_Grass = 2.0
        max_Interceptioncapacity_Rip = 3.0
        max_Ks = 0.07
        max_Kf = 0.6
        max_Soilstoaragecapacity_Grass = 200.0
        max_Soilstoaragecapacity_Rip = 200.0
        max_Soilstoaragecapacity_Bare = 50.0

        precission = 0.001
        precission_soilcap = 0.01

        beta_Bare = rand(0.1:precission: 2.0)
        beta_Forest = rand(0.1: precission: 2.0)
        beta_Grass = rand(0.1: precission: 2.0)
        beta_Rip = rand(0.1: precission: 2.0)
        Ce = rand(0.4: precission :0.8)
        Drainagecapacity = 0.0
        Interceptioncapacity_Bare = 0.0
        Interceptioncapacity_Forest = rand(1.0:precission:3.0)
        # parameter constraint Interception interceptioncapacity grass and rip lower than Interceptioncapacity_Forest
        if Interceptioncapacity_Forest < max_Interceptioncapacity_Grass
                Interceptioncapacity_Grass = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Grass = rand(0.0:precission: max_Interceptioncapacity_Grass)
        end

        if Interceptioncapacity_Forest < max_Interceptioncapacity_Rip
                Interceptioncapacity_Rip = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Rip = rand(0.0:precission: max_Interceptioncapacity_Rip)
        end
        # parameter constraints on fast reservoir coefficients
        Kf_Rip = rand(0.2:precission:3.0)
        if Kf_Rip < max_Kf
                Kf = rand(0.1:precission:Kf_Rip - precission)
        else
                Kf = rand(0.1:precission: max_Kf)
        end

        if Kf < max_Ks
                Ks = rand(0.001:precission * 0.1: Kf - precission)
        else
                Ks = rand(0.001:precission * 0.1: max_Ks)
        end
        Meltfactor = rand(1.0:precission:5.0)
        Mm = rand(0.001:precission * 0.1:1.0)
        Precipitation_Gradient = 0.0
        #Precipitation_Gradient = round(random_parameter(0, 0.0045), precission= 5)
        Ratio_Pref = rand(0.1:precission:0.9)
        # Parameter Constrain SOilstoragecapacity Forest >= Grass >= Rip/Bare
        Soilstoaragecapacity_Forest = rand(50.0:precission_soilcap:400.0)
        if Soilstoaragecapacity_Forest < max_Soilstoaragecapacity_Grass
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap:Soilstoaragecapacity_Forest - precission_soilcap)
        else
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap: max_Soilstoaragecapacity_Grass)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Rip
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap: max_Soilstoaragecapacity_Rip)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Bare
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap: max_Soilstoaragecapacity_Bare)
        end

        Temp_Thresh = rand(-2.0:precission:2.0)
        Ratio_Riparian = rand(0.05:precission:0.5)
         # based on calculation of recession curve

        bare_parameters = Parameters(beta_Bare, Ce, 0, Interceptioncapacity_Bare, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Bare, Temp_Thresh)
        forest_parameters = Parameters(beta_Forest, Ce, 0, Interceptioncapacity_Forest, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Forest, Temp_Thresh)
        grass_parameters = Parameters(beta_Grass, Ce, 0, Interceptioncapacity_Grass, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Grass, Temp_Thresh)
        rip_parameters = Parameters(beta_Rip, Ce, Drainagecapacity, Interceptioncapacity_Rip, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Rip, Temp_Thresh)
        slow_parameters = Slow_Paramters(Ks, Ratio_Riparian)

        parameters_array = [beta_Bare, beta_Forest, beta_Grass, beta_Rip, Ce, Interceptioncapacity_Forest, Interceptioncapacity_Grass, Interceptioncapacity_Rip, Kf_Rip, Kf, Ks, Meltfactor, Mm, Ratio_Pref, Ratio_Riparian, Soilstoaragecapacity_Bare, Soilstoaragecapacity_Forest, Soilstoaragecapacity_Grass, Soilstoaragecapacity_Rip, Temp_Thresh]

        return [bare_parameters, forest_parameters, grass_parameters, rip_parameters], slow_parameters, parameters_array
end

function parameter_selection_palten()
        # maximum parameter

        max_Interceptioncapacity_Grass = 2.0
        max_Interceptioncapacity_Rip = 3.0
        max_Ks = 0.07
        max_Kf = 0.5
        max_Soilstoaragecapacity_Grass = 250.0
        max_Soilstoaragecapacity_Rip = 200.0
        max_Soilstoaragecapacity_Bare = 50.0

        precission = 0.001
        precission_soilcap = 0.01

        beta_Bare = rand(0.1:precission: 2.0)
        beta_Forest = rand(0.1: precission: 2.0)
        beta_Grass = rand(0.1: precission: 2.0)
        beta_Rip = rand(0.1: precission: 2.0)
        Ce = rand(0.4: precission :0.8)
        Drainagecapacity = 0.0
        Interceptioncapacity_Bare = 0.0
        Interceptioncapacity_Forest = rand(1.0:precission:3.0)
        # parameter constraint Interception interceptioncapacity grass and rip lower than Interceptioncapacity_Forest
        if Interceptioncapacity_Forest < max_Interceptioncapacity_Grass
                Interceptioncapacity_Grass = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Grass = rand(0.0:precission: max_Interceptioncapacity_Grass)
        end

        if Interceptioncapacity_Forest < max_Interceptioncapacity_Rip
                Interceptioncapacity_Rip = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Rip = rand(0.0:precission: max_Interceptioncapacity_Rip)
        end
        # parameter constraints on fast reservoir coefficients
        Kf_Rip = rand(0.2:precission:3.0)
        if Kf_Rip < max_Kf
                Kf = rand(0.05:precission:Kf_Rip - precission)
        else
                Kf = rand(0.05:precission: max_Kf)
        end

        if Kf < max_Ks
                Ks = rand(0.001:precission * 0.1: Kf - precission)
        else
                Ks = rand(0.001:precission * 0.1: max_Ks)
        end
        Meltfactor = rand(1.0:precission:5.0)
        Mm = rand(0.001:precission * 0.1:1.0)
        Precipitation_Gradient = 0.0
        #Precipitation_Gradient = round(random_parameter(0, 0.0045), precission= 5)
        Ratio_Pref = rand(0.2:precission:0.8)
        # Parameter Constrain SOilstoragecapacity Forest >= Grass >= Rip/Bare
        Soilstoaragecapacity_Forest = rand(50.0:precission_soilcap:500.0)
        if Soilstoaragecapacity_Forest < max_Soilstoaragecapacity_Grass
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap:Soilstoaragecapacity_Forest - precission_soilcap)
        else
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap: max_Soilstoaragecapacity_Grass)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Rip
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap: max_Soilstoaragecapacity_Rip)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Bare
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap: max_Soilstoaragecapacity_Bare)
        end

        Temp_Thresh = rand(-2.0:precission:2.0)
        Ratio_Riparian = rand(0.05:precission:0.5)
         # based on calculation of recession curve

        bare_parameters = Parameters(beta_Bare, Ce, 0, Interceptioncapacity_Bare, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Bare, Temp_Thresh)
        forest_parameters = Parameters(beta_Forest, Ce, 0, Interceptioncapacity_Forest, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Forest, Temp_Thresh)
        grass_parameters = Parameters(beta_Grass, Ce, 0, Interceptioncapacity_Grass, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Grass, Temp_Thresh)
        rip_parameters = Parameters(beta_Rip, Ce, Drainagecapacity, Interceptioncapacity_Rip, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Rip, Temp_Thresh)
        slow_parameters = Slow_Paramters(Ks, Ratio_Riparian)

        parameters_array = [beta_Bare, beta_Forest, beta_Grass, beta_Rip, Ce, Interceptioncapacity_Forest, Interceptioncapacity_Grass, Interceptioncapacity_Rip, Kf_Rip, Kf, Ks, Meltfactor, Mm, Ratio_Pref, Ratio_Riparian, Soilstoaragecapacity_Bare, Soilstoaragecapacity_Forest, Soilstoaragecapacity_Grass, Soilstoaragecapacity_Rip, Temp_Thresh]

        return [bare_parameters, forest_parameters, grass_parameters, rip_parameters], slow_parameters, parameters_array
end

function parameter_selection_feistritz()
        # maximum parameter

        max_Interceptioncapacity_Grass = 2.0
        max_Interceptioncapacity_Rip = 3.0
        max_Ks = 0.05
        max_Kf = 0.6
        max_Soilstoaragecapacity_Grass = 200.0
        max_Soilstoaragecapacity_Rip = 150.0
        max_Soilstoaragecapacity_Bare = 50.0

        precission = 0.001
        precission_soilcap = 0.01

        beta_Bare = rand(0.1:precission: 2.0)
        beta_Forest = rand(0.1: precission: 2.0)
        beta_Grass = rand(0.1: precission: 2.0)
        beta_Rip = rand(0.1: precission: 2.0)
        Ce = rand(0.4: precission :0.8)
        Drainagecapacity = 0.0
        Interceptioncapacity_Bare = 0.0
        Interceptioncapacity_Forest = rand(1.0:precission:3.0)
        # parameter constraint Interception interceptioncapacity grass and rip lower than Interceptioncapacity_Forest
        if Interceptioncapacity_Forest < max_Interceptioncapacity_Grass
                Interceptioncapacity_Grass = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Grass = rand(0.0:precission: max_Interceptioncapacity_Grass)
        end

        if Interceptioncapacity_Forest < max_Interceptioncapacity_Rip
                Interceptioncapacity_Rip = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Rip = rand(0.0:precission: max_Interceptioncapacity_Rip)
        end
        # parameter constraints on fast reservoir coefficients
        Kf_Rip = rand(0.2:precission:3.0)
        if Kf_Rip < max_Kf
                Kf = rand(0.05:precission:Kf_Rip - precission)
        else
                Kf = rand(0.05:precission: max_Kf)
        end

        if Kf < max_Ks
                Ks = rand(0.001:precission * 0.1: Kf - precission)
        else
                Ks = rand(0.001:precission * 0.1: max_Ks)
        end
        Meltfactor = rand(1.0:precission:5.0)
        Mm = rand(0.001:precission * 0.1:1.0)
        Precipitation_Gradient = 0.0
        #Precipitation_Gradient = round(random_parameter(0, 0.0045), precission= 5)
        Ratio_Pref = rand(0.1:precission:0.7)
        # Parameter Constrain SOilstoragecapacity Forest >= Grass >= Rip/Bare
        Soilstoaragecapacity_Forest = rand(50.0:precission_soilcap:250.0)
        if Soilstoaragecapacity_Forest < max_Soilstoaragecapacity_Grass
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap:Soilstoaragecapacity_Forest - precission_soilcap)
        else
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap: max_Soilstoaragecapacity_Grass)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Rip
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap: max_Soilstoaragecapacity_Rip)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Bare
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap: max_Soilstoaragecapacity_Bare)
        end

        Temp_Thresh = rand(-2.0:precission:2.0)
        Ratio_Riparian = rand(0.05:precission:0.5)
         # based on calculation of recession curve

        bare_parameters = Parameters(beta_Bare, Ce, 0, Interceptioncapacity_Bare, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Bare, Temp_Thresh)
        forest_parameters = Parameters(beta_Forest, Ce, 0, Interceptioncapacity_Forest, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Forest, Temp_Thresh)
        grass_parameters = Parameters(beta_Grass, Ce, 0, Interceptioncapacity_Grass, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Grass, Temp_Thresh)
        rip_parameters = Parameters(beta_Rip, Ce, Drainagecapacity, Interceptioncapacity_Rip, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Rip, Temp_Thresh)
        slow_parameters = Slow_Paramters(Ks, Ratio_Riparian)

        parameters_array = [beta_Bare, beta_Forest, beta_Grass, beta_Rip, Ce, Interceptioncapacity_Forest, Interceptioncapacity_Grass, Interceptioncapacity_Rip, Kf_Rip, Kf, Ks, Meltfactor, Mm, Ratio_Pref, Ratio_Riparian, Soilstoaragecapacity_Bare, Soilstoaragecapacity_Forest, Soilstoaragecapacity_Grass, Soilstoaragecapacity_Rip, Temp_Thresh]

        return [bare_parameters, forest_parameters, grass_parameters, rip_parameters], slow_parameters, parameters_array
end

function parameter_selection_pitztal()
        # maximum parameter

        max_Interceptioncapacity_Grass = 2.0
        max_Interceptioncapacity_Rip = 3.0
        max_Ks = 0.07
        max_Kf = 0.5
        max_Soilstoaragecapacity_Grass = 250.0
        max_Soilstoaragecapacity_Rip = 200.0
        max_Soilstoaragecapacity_Bare = 50.0

        precission = 0.001
        precission_soilcap = 0.01

        beta_Bare = rand(0.1:precission: 2.0)
        beta_Forest = rand(0.1: precission: 2.0)
        beta_Grass = rand(0.1: precission: 2.0)
        beta_Rip = rand(0.1: precission: 2.0)
        Ce = rand(0.4: precission :0.8)
        Drainagecapacity = 0.0
        Interceptioncapacity_Bare = 0.0
        Interceptioncapacity_Forest = rand(1.0:precission:3.0)
        # parameter constraint Interception interceptioncapacity grass and rip lower than Interceptioncapacity_Forest
        if Interceptioncapacity_Forest < max_Interceptioncapacity_Grass
                Interceptioncapacity_Grass = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Grass = rand(0.0:precission: max_Interceptioncapacity_Grass)
        end

        if Interceptioncapacity_Forest < max_Interceptioncapacity_Rip
                Interceptioncapacity_Rip = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Rip = rand(0.0:precission: max_Interceptioncapacity_Rip)
        end
        # parameter constraints on fast reservoir coefficients
        Kf_Rip = rand(0.2:precission:3.0)
        if Kf_Rip < max_Kf
                Kf = rand(0.05:precission:Kf_Rip - precission)
        else
                Kf = rand(0.05:precission: max_Kf)
        end

        if Kf < max_Ks
                Ks = rand(0.001:precission * 0.1: Kf - precission)
        else
                Ks = rand(0.001:precission * 0.1: max_Ks)
        end
        Meltfactor = rand(2.0:precission:6.0)
        Mm = rand(0.001:precission * 0.1:1.0)
        Precipitation_Gradient = 0.0
        #Precipitation_Gradient = round(random_parameter(0, 0.0045), precission= 5)
        Ratio_Pref = rand(0.2:precission:0.9)
        # Parameter Constrain SOilstoragecapacity Forest >= Grass >= Rip/Bare
        Soilstoaragecapacity_Forest = rand(50.0:precission_soilcap:500.0)
        if Soilstoaragecapacity_Forest < max_Soilstoaragecapacity_Grass
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap:Soilstoaragecapacity_Forest - precission_soilcap)
        else
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap: max_Soilstoaragecapacity_Grass)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Rip
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap: max_Soilstoaragecapacity_Rip)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Bare
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap: max_Soilstoaragecapacity_Bare)
        end

        Temp_Thresh = rand(-2.5:precission:2.5)
        Ratio_Riparian = rand(0.05:precission:0.5)
         # based on calculation of recession curve

        bare_parameters = Parameters(beta_Bare, Ce, 0, Interceptioncapacity_Bare, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Bare, Temp_Thresh)
        forest_parameters = Parameters(beta_Forest, Ce, 0, Interceptioncapacity_Forest, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Forest, Temp_Thresh)
        grass_parameters = Parameters(beta_Grass, Ce, 0, Interceptioncapacity_Grass, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Grass, Temp_Thresh)
        rip_parameters = Parameters(beta_Rip, Ce, Drainagecapacity, Interceptioncapacity_Rip, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Rip, Temp_Thresh)
        slow_parameters = Slow_Paramters(Ks, Ratio_Riparian)

        parameters_array = [beta_Bare, beta_Forest, beta_Grass, beta_Rip, Ce, Interceptioncapacity_Forest, Interceptioncapacity_Grass, Interceptioncapacity_Rip, Kf_Rip, Kf, Ks, Meltfactor, Mm, Ratio_Pref, Ratio_Riparian, Soilstoaragecapacity_Bare, Soilstoaragecapacity_Forest, Soilstoaragecapacity_Grass, Soilstoaragecapacity_Rip, Temp_Thresh]

        return [bare_parameters, forest_parameters, grass_parameters, rip_parameters], slow_parameters, parameters_array
end

function parameter_selection_silbertal()
        # maximum parameter

        max_Interceptioncapacity_Grass = 2.0
        max_Interceptioncapacity_Rip = 3.0
        max_Ks = 0.1
        max_Kf = 1.0
        max_Soilstoaragecapacity_Grass = 250.0
        max_Soilstoaragecapacity_Rip = 250.0
        max_Soilstoaragecapacity_Bare = 50.0

        precission = 0.001
        precission_soilcap = 0.01

        beta_Bare = rand(0.1:precission: 2.0)
        beta_Forest = rand(0.1: precission: 2.0)
        beta_Grass = rand(0.1: precission: 2.0)
        beta_Rip = rand(0.1: precission: 2.0)
        Ce = rand(0.4: precission :0.8)
        Drainagecapacity = 0.0
        Interceptioncapacity_Bare = 0.0
        Interceptioncapacity_Forest = rand(1.0:precission:3.0)
        # parameter constraint Interception interceptioncapacity grass and rip lower than Interceptioncapacity_Forest
        if Interceptioncapacity_Forest < max_Interceptioncapacity_Grass
                Interceptioncapacity_Grass = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Grass = rand(0.0:precission: max_Interceptioncapacity_Grass)
        end

        if Interceptioncapacity_Forest < max_Interceptioncapacity_Rip
                Interceptioncapacity_Rip = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Rip = rand(0.0:precission: max_Interceptioncapacity_Rip)
        end
        # parameter constraints on fast reservoir coefficients
        Kf_Rip = rand(0.2:precission:3.0)
        if Kf_Rip < max_Kf
                Kf = rand(0.05:precission:Kf_Rip - precission)
        else
                Kf = rand(0.05:precission: max_Kf)
        end

        if Kf < max_Ks
                Ks = rand(0.001:precission * 0.1: Kf - precission)
        else
                Ks = rand(0.001:precission * 0.1: max_Ks)
        end
        Meltfactor = rand(1.0:precission:6.0)
        Mm = rand(0.001:precission * 0.1:1.0)
        Precipitation_Gradient = 0.0
        #Precipitation_Gradient = round(random_parameter(0, 0.0045), precission= 5)
        Ratio_Pref = rand(0.1:precission:0.9)
        # Parameter Constrain SOilstoragecapacity Forest >= Grass >= Rip/Bare
        Soilstoaragecapacity_Forest = rand(50.0:precission_soilcap:400.0)
        if Soilstoaragecapacity_Forest < max_Soilstoaragecapacity_Grass
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap:Soilstoaragecapacity_Forest - precission_soilcap)
        else
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap: max_Soilstoaragecapacity_Grass)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Rip
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap: max_Soilstoaragecapacity_Rip)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Bare
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap: max_Soilstoaragecapacity_Bare)
        end

        Temp_Thresh = rand(-2.:precission:2.)
        Ratio_Riparian = rand(0.05:precission:0.5)
         # based on calculation of recession curve

        bare_parameters = Parameters(beta_Bare, Ce, 0, Interceptioncapacity_Bare, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Bare, Temp_Thresh)
        forest_parameters = Parameters(beta_Forest, Ce, 0, Interceptioncapacity_Forest, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Forest, Temp_Thresh)
        grass_parameters = Parameters(beta_Grass, Ce, 0, Interceptioncapacity_Grass, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Grass, Temp_Thresh)
        rip_parameters = Parameters(beta_Rip, Ce, Drainagecapacity, Interceptioncapacity_Rip, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Rip, Temp_Thresh)
        slow_parameters = Slow_Paramters(Ks, Ratio_Riparian)

        parameters_array = [beta_Bare, beta_Forest, beta_Grass, beta_Rip, Ce, Interceptioncapacity_Forest, Interceptioncapacity_Grass, Interceptioncapacity_Rip, Kf_Rip, Kf, Ks, Meltfactor, Mm, Ratio_Pref, Ratio_Riparian, Soilstoaragecapacity_Bare, Soilstoaragecapacity_Forest, Soilstoaragecapacity_Grass, Soilstoaragecapacity_Rip, Temp_Thresh]

        return [bare_parameters, forest_parameters, grass_parameters, rip_parameters], slow_parameters, parameters_array
end

function parameter_selection_defreggental()
        # maximum parameter

        max_Interceptioncapacity_Grass = 2.0
        max_Interceptioncapacity_Rip = 3.0
        max_Ks = 0.07
        max_Kf = 0.5
        max_Soilstoaragecapacity_Grass = 250.0
        max_Soilstoaragecapacity_Rip = 250.0
        max_Soilstoaragecapacity_Bare = 50.0

        precission = 0.001
        precission_soilcap = 0.01

        beta_Bare = rand(0.1:precission: 2.0)
        beta_Forest = rand(0.1: precission: 2.0)
        beta_Grass = rand(0.1: precission: 2.0)
        beta_Rip = rand(0.1: precission: 2.0)
        Ce = rand(0.4: precission :0.8)
        Drainagecapacity = 0.0
        Interceptioncapacity_Bare = 0.0
        Interceptioncapacity_Forest = rand(1.0:precission:3.0)
        # parameter constraint Interception interceptioncapacity grass and rip lower than Interceptioncapacity_Forest
        if Interceptioncapacity_Forest < max_Interceptioncapacity_Grass
                Interceptioncapacity_Grass = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Grass = rand(0.0:precission: max_Interceptioncapacity_Grass)
        end

        if Interceptioncapacity_Forest < max_Interceptioncapacity_Rip
                Interceptioncapacity_Rip = rand(0.0:precission:Interceptioncapacity_Forest - precission)
        else
                Interceptioncapacity_Rip = rand(0.0:precission: max_Interceptioncapacity_Rip)
        end
        # parameter constraints on fast reservoir coefficients
        Kf_Rip = rand(0.2:precission:3.0)
        if Kf_Rip < max_Kf
                Kf = rand(0.05:precission:Kf_Rip - precission)
        else
                Kf = rand(0.05:precission: max_Kf)
        end

        if Kf < max_Ks
                Ks = rand(0.001:precission * 0.1: Kf - precission)
        else
                Ks = rand(0.001:precission * 0.1: max_Ks)
        end
        Meltfactor = rand(1.0:precission:6.0)
        Mm = rand(0.001:precission * 0.1:1.0)
        Precipitation_Gradient = 0.0
        #Precipitation_Gradient = round(random_parameter(0, 0.0045), precission= 5)
        Ratio_Pref = rand(0.2:precission:0.9)
        # Parameter Constrain SOilstoragecapacity Forest >= Grass >= Rip/Bare
        Soilstoaragecapacity_Forest = rand(50.0:precission_soilcap:400.0)
        if Soilstoaragecapacity_Forest < max_Soilstoaragecapacity_Grass
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap:Soilstoaragecapacity_Forest - precission_soilcap)
        else
                Soilstoaragecapacity_Grass = rand(5.0:precission_soilcap: max_Soilstoaragecapacity_Grass)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Rip
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Rip = rand(4.9:precission_soilcap: max_Soilstoaragecapacity_Rip)
        end

        if Soilstoaragecapacity_Grass < max_Soilstoaragecapacity_Bare
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap:Soilstoaragecapacity_Grass - precission_soilcap)
        else
                Soilstoaragecapacity_Bare = rand(1.0:precission_soilcap: max_Soilstoaragecapacity_Bare)
        end

        Temp_Thresh = rand(-0.:precission:2.5)
        Ratio_Riparian = rand(0.05:precission:0.5)
         # based on calculation of recession curve

        bare_parameters = Parameters(beta_Bare, Ce, 0, Interceptioncapacity_Bare, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Bare, Temp_Thresh)
        forest_parameters = Parameters(beta_Forest, Ce, 0, Interceptioncapacity_Forest, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Forest, Temp_Thresh)
        grass_parameters = Parameters(beta_Grass, Ce, 0, Interceptioncapacity_Grass, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Grass, Temp_Thresh)
        rip_parameters = Parameters(beta_Rip, Ce, Drainagecapacity, Interceptioncapacity_Rip, Kf, Meltfactor, Mm, Ratio_Pref, Soilstoaragecapacity_Rip, Temp_Thresh)
        slow_parameters = Slow_Paramters(Ks, Ratio_Riparian)

        parameters_array = [beta_Bare, beta_Forest, beta_Grass, beta_Rip, Ce, Interceptioncapacity_Forest, Interceptioncapacity_Grass, Interceptioncapacity_Rip, Kf_Rip, Kf, Ks, Meltfactor, Mm, Ratio_Pref, Ratio_Riparian, Soilstoaragecapacity_Bare, Soilstoaragecapacity_Forest, Soilstoaragecapacity_Grass, Soilstoaragecapacity_Rip, Temp_Thresh]

        return [bare_parameters, forest_parameters, grass_parameters, rip_parameters], slow_parameters, parameters_array
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
#
# This file is part of the SlowActors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#
"""
    SlowActors

A minimal Julia actors library for using the `Actors` interface.

It is not registered since it has demonstrative purpose.
You can install it by:

```julia
pkg> add "https://github.com/JuliaActors/SlowActors.jl"
```
"""
module SlowActors

"Gives the package version."
const version = v"0.1.0"

using Reexport, DataStructures, Distributed

@reexport using Actors

include("types.jl")
include("actor.jl")
include("com.jl")

export newLink, spawn

end


text:
# Generate documentation with this command:
# (cd docs && julia --color=yes make.jl)

push!(LOAD_PATH, "..")

using Documenter
using FastSphericalHarmonics

makedocs(; sitename="FastSphericalHarmonics", format=Documenter.HTML(),
         modules=[FastSphericalHarmonics])

deploydocs(; repo="github.com/eschnett/FastSphericalHarmonics.jl.git",
           devbranch="main", push_preview=true)


text:
"""
    Problem{P, Q}(network::Network, input::P, output::Q)

Problem definition for neural verification.

The verification problem consists of: for all  points in the input set,
the corresponding output of the network must belong to the output set.
"""
struct Problem{P, Q}
    network::Network
    input::P
    output::Q
end

"""
    Result
Supertype of all result types.

See also: [`BasicResult`](@ref), [`CounterExampleResult`](@ref), [`AdversarialResult`](@ref), [`ReachabilityResult`](@ref)
"""
abstract type Result end

status(result::Result) = result.status

function validate_status(st::Symbol)
    @assert st ∈ (:holds, :violated, :unknown) "unexpected status code: `:$st`.\nOnly (:holds, :violated, :unknown) are accepted"
    return st
end

"""
    BasicResult(status::Symbol)

Result type that captures whether the input-output constraint is satisfied.
Possible status values:\n
    :holds (io constraint is satisfied always)\n
    :violated (io constraint is violated)\n
    :unknown (could not be determined)
"""
struct BasicResult <: Result
    status::Symbol
end

"""
    CounterExampleResult(status, counter_example)

Like `BasicResult`, but also returns a `counter_example` if one is found (if status = :violated).
The `counter_example` is a point in the input set that, after the NN, lies outside the output set.
"""
struct CounterExampleResult <: Result
    status::Symbol
    counter_example::Vector{Float64}
    CounterExampleResult(s, ce) = new(validate_status(s), ce)
end

"""
    AdversarialResult(status, max_disturbance)

Like `BasicResult`, but also returns the maximum allowable disturbance in the input (if status = :violated).
"""
struct AdversarialResult <: Result
	status::Symbol
	max_disturbance::Float64
    AdversarialResult(s, md) = new(validate_status(s), md)
end

"""
    ReachabilityResult(status, reachable)

Like `BasicResult`, but also returns the output reachable set given the input constraint (if status = :violated).
"""
struct ReachabilityResult <: Result
	status::Symbol
	reachable::Vector{<:AbstractPolytope}
    ReachabilityResult(s, r) = new(validate_status(s), r)
end

# Additional constructors:
CounterExampleResult(s) = CounterExampleResult(s, Float64[])
AdversarialResult(s)    = AdversarialResult(s, -1.0)
ReachabilityResult(s)   = ReachabilityResult(s, AbstractPolytope[])


text:
import Zygote, ForwardDiff
using Random, Test
using YaoBlocks, YaoArrayRegister
using ChainRulesCore: Tangent

@testset "recursive_create_tangent" begin
    c = chain(put(5, 2 => chain(Rx(1.4), Rx(0.5))), cnot(5, 3, 1), put(5, 3 => Rx(-0.5)))
    tc = YaoBlocks.AD.recursive_create_tangent(c)
    @test tc isa Tangent
end

@testset "construtors" begin
    @test Zygote.gradient(x -> x.list[1].blocks[1].theta, sum([chain(1, Rz(0.3))]))[1] == (
        list = NamedTuple{
            (:blocks,),
            Tuple{Vector{NamedTuple{(:block, :theta),Tuple{Nothing,Float64}}}},
        }[(blocks = [(block = nothing, theta = 1.0)],)],
    )
    @test_broken Zygote.gradient(x -> getfield(getfield(x, :content), :theta), Daggered(Rx(0.5)))[1] ==
                 (content = (block = nothing, theta = 1.0),)
end

@testset "rules" begin
    h = put(5, 3 => Z) + put(5, 2 => X)
    c = chain(put(5, 2 => chain(Rx(1.4), Rx(0.5))), cnot(5, 3, 1), put(5, 3 => Rx(-0.5)))
    r = rand_state(5)
    g0 = reinterpret(
        ComplexF64,
        ForwardDiff.gradient(
            x -> real(expect(h, ArrayReg([Complex(x[2i-1], x[2i]) for i in 1:length(x)÷2]))),
            reinterpret(Float64, r.state),
        ),
    )
    @test Zygote.gradient(x -> real(expect(h, ArrayReg(x))), r.state)[1] ≈ g0
    @test Zygote.gradient(x -> real(expect(h, ArrayReg{1}(reshape(statevec(x), :, 1)))), r)[1].state ≈
          g0
    @test Zygote.gradient(x -> real(expect(h, ArrayReg(reshape(state(x), :, 1)))), r)[1].state ≈ g0
    @test Zygote.gradient(x -> real(expect(h, copy(x))), r)[1].state ≈ g0
    @test Zygote.gradient(x -> real(expect(h, parent(x'))), r)[1].state ≈ g0

    g1 = reinterpret(
        ComplexF64,
        ForwardDiff.gradient(
            x -> real(sum(abs2, [Complex(x[2i-1], x[2i]) for i in 1:length(x)÷2])),
            reinterpret(Float64, r.state),
        ),
    )
    @test Zygote.gradient(x -> real(sum(abs2, state(x'))), r)[1].state ≈ g1
    @test Zygote.gradient(x -> real(sum(abs2, statevec(x'))), r)[1].state ≈ g1
    # zygote does not work if `sin` is not here,
    # because it gives an adjoint of different type as the output matrix type.
    @test AD.extract_circuit_gradients!(
        Zygote.gradient(x -> real(sum(sin, Matrix(x))), c)[1].blocks,
        Float64[],
    ) ≈ ForwardDiff.gradient(x -> real(sum(sin, Matrix(dispatch(c, x)))), parameters(c))
end

@testset "adwith zygote" begin
    c = chain(put(5, 2 => chain(Rx(0.4), Rx(0.5))), cnot(5, 3, 1), put(5, 3 => Rx(-0.5)))
    dispatch!(c, :random)

    function loss(reg::AbstractRegister, circuit::AbstractBlock{N}) where {N}
        reg = apply(copy(reg), circuit)
        st = state(reg)
        sum(real(st .* st))
    end

    # apply
    reg0 = zero_state(5)
    params = rand!(parameters(c))
    paramsδ = Zygote.gradient(params -> loss(reg0, dispatch(c, params)), params)[1]
    regδ = Zygote.gradient(reg -> loss(reg, c), reg0)[1]
    fparamsδ = ForwardDiff.gradient(
        params ->
            loss(ArrayReg(Matrix{Complex{eltype(params)}}(reg0.state)), dispatch(c, params)),
        params,
    )
    fregδ = ForwardDiff.gradient(
        x -> loss(
            ArrayReg([Complex(x[2i-1], x[2i]) for i in 1:length(x)÷2]),
            dispatch(c, Vector{real(eltype(x))}(parameters(c))),
        ),
        reinterpret(Float64, reg0.state),
    )
    @test fregδ ≈ reinterpret(Float64, regδ.state)
    @test fparamsδ ≈ paramsδ

    # expect and fidelity
    c = chain(
        put(5, 5 => Rx(1.5)),
        put(5, 1 => Rx(0.4)),
        put(5, 4 => Rx(0.2)),
        put(5, 2 => chain(Rx(0.4), Rx(0.5))),
        cnot(5, 3, 1),
        put(5, 3 => Rx(-0.5)),
    )
    h = chain(repeat(5, X, 1:5))
    reg = rand_state(5)
    function loss2(reg::AbstractRegister, circuit::AbstractBlock{N}) where {N}
        return 5 * real(expect(h, copy(reg) => circuit) + fidelity(reg, apply(reg, circuit)))
    end
    params = rand!(parameters(c))
    fδc = ForwardDiff.gradient(
        params ->
            loss2(ArrayReg(Matrix{Complex{eltype(params)}}(reg.state)), dispatch(c, params)),
        params,
    )
    δr, δc = Zygote.gradient((reg, params) -> loss2(reg, dispatch(c, params)), reg, params)
    @test δc ≈ fδc

    fregδ = ForwardDiff.gradient(
        x -> loss2(
            ArrayReg([Complex(x[2i-1], x[2i]) for i in 1:length(x)÷2]),
            dispatch(c, Vector{real(eltype(x))}(params)),
        ),
        reinterpret(Float64, reg.state),
    )
    @test fregδ ≈ reinterpret(Float64, δr.state)

    # operator fidelity
    c = chain(
        put(5, 5 => Rx(1.5)),
        put(5, 1 => Rx(0.4)),
        put(5, 4 => Rx(0.2)),
        put(5, 2 => chain(Rx(0.4), Rx(0.5))),
        cnot(5, 3, 1),
        put(5, 3 => Rx(-0.5)),
    )
    h = chain(repeat(5, X, 1:5))
    function loss3(circuit::AbstractBlock{N}, h) where {N}
        return operator_fidelity(circuit, h)
    end
    params = rand!(parameters(c))
    fδc = ForwardDiff.gradient(params -> loss3(dispatch(c, params), h), params)
    δc, = Zygote.gradient(p -> loss3(dispatch(c, p), h), params)
    @test δc ≈ fδc

    # NOTE: operator back propagation in expect is not implemented!
    # to differentiate operators, we need to use the expensive `mat_back` function.
end

@testset "add block" begin
    H = sum([chain(5, put(k => Z)) for k in 1:5])
    c = chain(put(5, 2 => chain(Rx(0.4), Rx(0.5))), cnot(5, 3, 1), put(5, 3 => Rx(-0.5)))
    dispatch!(c, :random)
    function loss(reg::AbstractRegister, circuit::AbstractBlock{N}) where {N}
        reg = apply(copy(reg), circuit)
        st = state(reg)
        reg2 = apply(copy(reg), H)
        st2 = state(reg2)
        sum(real(st .* st2))
    end

    reg0 = zero_state(5)
    params = rand!(parameters(c))
    paramsδ = Zygote.gradient(params -> loss(reg0, dispatch(c, params)), params)[1]
    fparamsδ = ForwardDiff.gradient(
        params ->
            loss(ArrayReg(Matrix{Complex{eltype(params)}}(reg0.state)), dispatch(c, params)),
        params,
    )
    @test fparamsδ ≈ paramsδ

    regδ = Zygote.gradient(reg -> loss(reg, c), reg0)[1]
    fregδ = ForwardDiff.gradient(
        x -> loss(
            ArrayReg([Complex(x[2i-1], x[2i]) for i in 1:length(x)÷2]),
            dispatch(c, Vector{real(eltype(x))}(parameters(c))),
        ),
        reinterpret(Float64, reg0.state),
    )
    @test fregδ ≈ reinterpret(Float64, regδ.state)
end

max_stars_repo_name: QuantumBFS/YaoStdCircuit

max_issues_repo_name: QuantumBFS/YaoStdCircuit


text:
module SuccessiveConvexifications
"""
module SuccessiveConvexifications a.k.a. SCvx

[Ref] Mao, Y., Szmuk, M., Xu, X., & Acikmese, B. (2018).
Successive Convexification:
A Superlinearly Convergent Algorithm for Non-convex Optimal Control
Problems.

[Note]
N > 0: discrete-time horizon length
x ∈ R^(n_x, 1): state
u ∈ R^(n_u, 1): input
X ∈ R^(N, n_x): concatenated state, X = [x1 x2 ... xN]
U ∈ R^((N-1), n_u): concatenated input, U = [u1 u2 ... u(N-1)]
D ∈ R^(N, n_x): optimisation variables of convex subproblem
w.r.t state, D = [d1 d2 ... dN]
W ∈ R^((N-1), n_u): optimisation variables of convex subproblem
w.r.t input, W = [w1 w2 ... w(N-1)]
"""

using Parameters
using Convex
const Cvx = Convex  # to clarify the code
using Mosek
using MosekTools
const SOLVER = Mosek
# using Debugger

using Reexport
include("Linearisers.jl")
@reexport using .Linearisers

export SCvx
export initial_guess!, solve!, flush!, append!

##################### Extension ####################
function Convex.pos(x::Array)
    y = deepcopy(x)
    for i in 1:length(y)
        if y[i] < 0
            y[i] = 0.0
        end
    end
    return y
end

##################### SuccessiveConvexifications ####################
"""
# Note
`objs` and `consts` is an array containing obj. funcs. and constraint funcs.,
resp.
1) The objective function is the sum of obj. funcs. in `objs` arrays.
2) The constraint function is combination of all const. funcs. in `consts` arrays.
    1) Inequality constraint: f ≤ 0
    2) Equality constraint: f = 0
# Arguments
- `objs_(option)::Array`: contains objective functions.
    - option: `path` or `terminal`
- `consts_(option1)_(option2)::Array`: contains constraint functions.
    - option1: `path`, `initial` or `terminal
    - option2: `eq` or `ineq`
"""
@with_kw mutable struct SCvx
    # hyperparameters
    N::Int64;  @assert N > 2
    i::Int64 = 0
    i_max::Int64 = 200
    n_x::Int64
    n_u::Int64
    r_k::Float64 = 1e0
    λ::Float64 = 1e5
    ρ0::Float64 = 0.0
    ρ1::Float64 = 0.25
    ρ2::Float64 = 0.7
    rl::Float64 = 1.0e-1
    α::Float64 = 2.0e0
    β::Float64 = 2.0e0
    ϵ::Float64 = 1e-3
    # problem
    objs_path::Array{Function} = []
    objs_terminal::Array{Function} = []
    consts_path_ineq::Array{Function} = []
    consts_path_eq::Array{Function} = []
    consts_initial_ineq::Array{Function} = []
    consts_initial_eq::Array{Function} = []
    consts_terminal_ineq::Array{Function} = []
    consts_terminal_eq::Array{Function} = []
    # default initial guess
    X_k::Array = zeros(N, n_x)
    U_k::Array = zeros(N-1, n_u)
    X_k_sol = [X_k]  # contain kth solutions including zeroth sol
    U_k_sol = [U_k]  # contain kth solutions including zeroth sol
    obj_k::Array{Number} = []
    flag::Array{String} = []
end

function initial_guess!(scvx::SCvx, X_k::Array, U_k::Array)
    if scvx.N != size(X_k)[1] || scvx.N != (size(U_k)[1]+1)
        error("Invalid length")
    elseif scvx.n_x != size(X_k)[2] || scvx.n_u != size(U_k)[2]
        error("Dimension disagreement")
    end
    scvx.X_k, scvx.U_k = X_k, U_k
    return scvx
end

"""
    flush!(scvx)

Flush the initial values of solution.
You would need this when applying this algorithm to fixed final time problem.
"""
function flush!(scvx::SCvx)
    scvx.N = scvx.N - 1
    scvx.X_k = scvx.X_k[2:end, :]
    scvx.U_k = scvx.U_k[2:end, :]
end

function Base.append!(scvx::SCvx, symbol, func::Array{<:Function})
    symbol = Symbol(symbol)
    append!(getproperty(scvx, symbol), func)
end

function Cvx.solve!(scvx::SCvx; verbose::Bool=false)
    is_not_skipped = true
    is_not_stopped = true
    scvx.i = 0
    while is_not_stopped
        # step 1: solve convex subproblem
        D, W, jacob_dict = solve_cvx_subprob(scvx, verbose=verbose)
        # step 2: compute the actual and predicted change
        J_k, obj_k = get_J(scvx, scvx.X_k, scvx.U_k, with_obj=true)
        diff_J = J_k - get_J(scvx, scvx.X_k+D.value, scvx.U_k+W.value)
        diff_L = J_k - get_L(scvx, D.value, W.value, jacob_dict)
        # if diff_J == 0.0  # TODO: should be diff_L == 0.0 ?
        if diff_L == 0.0
            flag = "diff_L = 0"
            break
        else
            ρ_k = diff_J / diff_L
        end
        # step 3: compare the ratio
        if ρ_k <= scvx.ρ0
            scvx.r_k = scvx.r_k /scvx.α  # reject and contracted
            is_not_skipped = false  # go to step 1
            flag = "rejected"
        else
            is_not_skipped = true
        end
        if is_not_skipped
            scvx.X_k, scvx.U_k = (
                                  scvx.X_k + D.value, scvx.U_k + W.value
                                 )  # accept
            if ρ_k < scvx.ρ1
                r_k = scvx.r_k / scvx.α
                flag = "accept but contracted"
            elseif ρ_k < scvx.ρ2
                r_k = scvx.r_k
                flag = "accept and maintain"
            else
                r_k = scvx.r_k * scvx.β
                flag = "accept and enlarged"
            end
            scvx.r_k = max(r_k, scvx.rl)
            is_not_stopped = diff_J > scvx.ϵ && scvx.i < scvx.i_max
            if !is_not_stopped
                if scvx.i >= scvx.i_max
                    flag = "maximum iteration"
                elseif diff_J <= scvx.ϵ
                    flag = "success"
                end
            end
        end
        update!(scvx, obj_k, flag)
    end
    return scvx
end

function update!(scvx, obj_k, flag; verbose=false)
    scvx.i += 1
    push!(scvx.obj_k, obj_k)
    push!(scvx.flag, flag)
    push!(scvx.X_k_sol, scvx.X_k)
    push!(scvx.U_k_sol, scvx.U_k)
    if verbose && scvx.i % 10 == 0
        println("iter = $(scvx.i)")
    end
end

function solve_cvx_subprob(scvx::SCvx; verbose::Bool=false)
    N = scvx.N
    n_x, n_u = scvx.n_x, scvx.n_u
    D, W = Cvx.Variable(N, n_x), Cvx.Variable(N-1, n_u)
    ## problem formulation
    # constraints
    constraints = [norm([reshape(D, prod(size(D)), 1);
                         reshape(W, prod(size(W)), 1)], 2) <= scvx.r_k]
    # linearised and penalised objective
    obj_linearised, jacob_dict = get_obj_linearised(scvx,
                                                    D, W,
                                                    jacob_dict=Dict())
    (
     obj_const_penalty_linearised,
     jacob_dict
    ) = get_obj_const_penalty_linearised(scvx,
                                         D, W, jacob_dict=jacob_dict)
    # TODO: add convex objective and constraints for path and terminal
    obj_linearised_penalised = (obj_linearised
                                + obj_const_penalty_linearised)
    prob = minimize(obj_linearised_penalised, constraints)
    if SOLVER == Mosek
        Cvx.solve!(prob, SOLVER.Optimizer(LOG=verbose))  # for Mosek
    else
        Cvx.solve!(prob, SOLVER.Optimizer(verbose=verbose))  # for others
    end
    return D, W, jacob_dict
end

################## get_J ##################
function get_J(scvx::SCvx, X, U; with_obj=false)
    obj = get_obj(scvx, X, U)
    const_penalty = get_obj_const_penalty(scvx, X, U)
    J = obj + const_penalty
    if with_obj
        return J, obj
    else
        return J
    end
end

## get_obj
function get_obj(scvx::SCvx, X, U)
    obj = 0.0
    obj += calculate_objs_path(scvx, X, U)
    obj += calculate_objs_terminal(scvx, X)
    return obj
end

function calculate_objs_path(scvx::SCvx, X, U)::Number
    # path obj
    N = scvx.N
    obj = zeros(1)
    if scvx.objs_path == []
    else
        for n in 1:N-1
            for obj_path in scvx.objs_path
                X_n, U_n = X[n, :], U[n, :]
                obj += obj_path(X_n, U_n)
            end
        end
    end
    return obj[1]
end

function calculate_objs_terminal(scvx::SCvx, X)::Number
    # terminal obj
    obj = zeros(1)
    if scvx.objs_terminal == []
    else
        for obj_terminal in scvx.objs_terminal
            X_N = X[end, :]
            obj += obj_terminal(X_N)
        end
    end
    return obj[1]
end

## get_obj_const_penalty
function get_obj_const_penalty(scvx::SCvx, X, U)
    obj = 0.0
    obj += calculate_const_path_penalty(scvx, X, U)
    obj += calculate_const_initial_penalty(scvx, X, U)
    obj += calculate_const_terminal_penalty(scvx, X)
    return obj
end

function calculate_const_path_penalty(scvx::SCvx, X, U)
    N = scvx.N
    obj = 0.0
    for n in 1:N-1
        X_n, U_n, X_n_next = X[n, :], U[n, :], X[n+1, :]
        if scvx.consts_path_ineq == []
        else
            for const_ineq in scvx.consts_path_ineq
                obj += calculate_ineq_penalty(scvx, const_ineq, X_n, U_n)
            end
        end
        if scvx.consts_path_eq == []
        else
            for const_eq in scvx.consts_path_eq
                obj += calculate_eq_penalty(scvx,
                                            const_eq, X_n, U_n, X_n_next)
            end
        end
    end
    return obj
end

function calculate_const_initial_penalty(scvx::SCvx, X, U)
    X_1, U_1 = X[1, :], U[1, :]
    obj = 0.0
    if scvx.consts_initial_ineq == []
    else
        for const_ineq in scvx.consts_initial_ineq
            obj += calculate_ineq_penalty(scvx, const_ineq, X_1, U_1)
        end
    end
    if scvx.consts_initial_eq == []
    else
        for const_eq in scvx.consts_initial_eq
            obj += calculate_eq_penalty(scvx, const_eq, X_1, U_1)
        end
    end
    return obj
end

function calculate_const_terminal_penalty(scvx::SCvx, X)
    X_N = X[end, :]
    obj = 0.0
    if scvx.consts_terminal_ineq == []
    else
        for const_ineq in scvx.consts_terminal_ineq
            obj += calculate_ineq_penalty(scvx, const_ineq, X_N)
        end
    end
    if scvx.consts_terminal_eq == []
    else
        for const_eq in scvx.consts_terminal_eq
            obj += calculate_eq_penalty(scvx, const_eq, X_N)
        end
    end
    return obj
end

################## get_L ##################
function get_L(scvx::SCvx, D::Array, W::Array, jacob_dict)::Number
    obj_linearised, _ = get_obj_linearised(scvx,
                                           D, W,
                                           jacob_dict=jacob_dict)
    (
     obj_const_penalty_linearised, _
    ) = get_obj_const_penalty_linearised(scvx,
                                         D, W, jacob_dict=jacob_dict)
    L = obj_linearised + obj_const_penalty_linearised
    return L
end

## get_obj_linearised
function get_obj_linearised(scvx::SCvx, D, W; jacob_dict=Dict())
    obj = 0.0
    obj_path, jacob_dict = calculate_objs_path_linearised(scvx,
                                                          D, W,
                                                          jacob_dict)
    obj += obj_path  # path
    (
     obj_terminal, jacob_dict
    ) = calculate_objs_terminal_linearised(scvx, D, W, jacob_dict)
    obj += obj_terminal  # terminal
    return obj, jacob_dict
end

function calculate_objs_path_linearised(scvx::SCvx, D, W, jacob_dict)
    name = "objs_path"
    if !haskey(jacob_dict, name)
        jacob_dict[name] = []
    end
    N = scvx.N
    X_k, U_k = scvx.X_k, scvx.U_k  # linearisation point
    obj = 0.0
    DW = [D[1:end-1, :] W]
    if scvx.objs_path == []
    else
        for n in 1:N-1
            X_k_n, U_k_n = X_k[n, :], U_k[n, :]
            for obj_path in scvx.objs_path
                DW_n = DW[n, :]
                (
                 jacob, jacob_dict
                ) = use_or_calculate_jacobian(name, size(jacob_dict[name])[1] < N-1, jacob_dict, obj_path, X_k_n, U_k_n)
                obj_path_tmp = (obj_path(X_k_n, U_k_n)
                                + jacob * reshape(DW_n,
                                                  prod(size(DW_n)), 1)
                               )[1]  # reshape for Cvx.Variable
                obj += obj_path_tmp
            end
        end
    end
    return obj, jacob_dict
end

function calculate_objs_terminal_linearised(scvx::SCvx, D, W, jacob_dict)
    name = "objs_terminal"
    if !haskey(jacob_dict, name)
        jacob_dict[name] = []
    end
    obj = 0.0
    X_k = scvx.X_k
    if scvx.objs_terminal == []
    else
        for obj_terminal in scvx.objs_terminal
            X_k_N = X_k[end, :]
            D_N = D[end, :]
            name = "objs_terminal"
            (
             jacob, jacob_dict
            ) = use_or_calculate_jacobian(name, size(jacob_dict[name])[1] == 0, jacob_dict, obj_terminal, X_k_N)
            obj_tmp = (obj_terminal(X_k_N)
                       + jacob * reshape(D_N, prod(size(D_N)), 1)
                      )[1]  # reshape for Cvx.Variable
            obj += obj_tmp
        end
    end
    return obj, jacob_dict
end

## get_obj_const_penalty_linearised
function get_obj_const_penalty_linearised(scvx::SCvx, D, W;
                                          jacob_dict=Dict())
    obj = 0.0
    (
     obj_path, jacob_dict
    ) = calculate_const_path_penalty_linearised(scvx, D, W, jacob_dict)
    obj += obj_path  # path
    (
     obj_initial, jacob_dict
    ) = calculate_const_initial_penalty_linearised(scvx, D, W, jacob_dict)
    obj += obj_initial  # initial
    (
     obj_terminal, jacob_dict
    ) = calculate_const_terminal_penalty_linearised(scvx,
                                                    D, W, jacob_dict)
    obj += obj_terminal  # terminal
    return obj, jacob_dict
end

function calculate_const_path_penalty_linearised(scvx::SCvx,
                                                 D, W, jacob_dict)
    name_ineq = "consts_path_ineq"
    if !haskey(jacob_dict, name_ineq)
        jacob_dict[name_ineq] = []
    end
    name_eq = "consts_path_eq"
    if !haskey(jacob_dict, name_eq)
        jacob_dict[name_eq] = []
    end
    obj = 0.0
    N = scvx.N
    X_k, U_k = scvx.X_k, scvx.U_k  # linearisation point
    DW = [D[1:end-1, :] W]
    DW_extended = [D[1:end-1, :] W D[2:end, :]]

    for n in 1:N-1
        X_k_n, U_k_n, X_k_n_next = X_k[n, :], U_k[n, :], X_k[n+1, :]
        DW_n = DW[n, :]
        D_n, W_n, D_n_next = D[n, :], W[n, :], D[n+1, :]
        DW_n_extended = DW_extended[n, :]
        if scvx.consts_path_ineq == []
        else
            for const_ineq in scvx.consts_path_ineq
                (
                 jacob, jacob_dict
                ) = use_or_calculate_jacobian(name_ineq, size(jacob_dict[name_ineq])[1] < N-1, jacob_dict, const_ineq, X_k_n, U_k_n, n=n)
                const_ineq_tmp(
                               X_k_n, U_k_n, DW_n, jacob
                              ) = (
                                   const_ineq(X_k_n, U_k_n)
                                   + jacob * reshape(DW_n,
                                                     prod(size(DW_n)), 1)
                                  )  # reshape for Cvx.Variable
                obj += calculate_ineq_penalty(scvx, const_ineq_tmp,
                                              X_k_n, U_k_n, DW_n, jacob)
            end
        end
        if scvx.consts_path_eq == []
        else
            for const_eq in scvx.consts_path_eq
                (
                 jacob, jacob_dict
                ) = use_or_calculate_jacobian(name_eq, size(jacob_dict[name_ineq])[1] < N-1, jacob_dict, const_eq, X_k_n, U_k_n, X_k_n_next, n=n)
                const_eq_tmp(
                             X_k_n, U_k_n, X_k_n_next,
                             DW_n_extended, jacob
                            ) = (
                                 const_eq(X_k_n, U_k_n, X_k_n_next)
                                 + jacob * reshape(DW_n_extended, prod(size(DW_n_extended)), 1)
                                )  # reshape for Cvx.Variable
                obj += calculate_eq_penalty(scvx, const_eq_tmp,
                                            X_k_n, U_k_n, X_k_n_next,
                                            DW_n_extended, jacob)
            end
        end
    end
    return obj, jacob_dict
end

function calculate_const_initial_penalty_linearised(scvx::SCvx,
                                                    D, W, jacob_dict)
    name_ineq = "consts_initial_ineq"
    if !haskey(jacob_dict, name_ineq)
        jacob_dict[name_ineq] = []
    end
    name_eq = "consts_initial_eq"
    if !haskey(jacob_dict, name_eq)
        jacob_dict[name_eq] = []
    end
    obj = 0.0
    X_k, U_k = scvx.X_k, scvx.U_k  # linearisation point
    X_k_1, U_k_1 = X_k[1, :], U_k[1, :]
    DW = [D[1:end-1, :] W]
    DW_1 = DW[1, :]
    D_1, W_1 = D[1, :], W[1, :]
    if scvx.consts_initial_ineq == []
    else
        for const_ineq in scvx.consts_initial_ineq
            (
             jacob, jacob_dict
            ) = use_or_calculate_jacobian(name_ineq, size(jacob_dict[name_ineq])[1] == 0, jacob_dict, const_ineq, X_k_1, U_k_1)
            const_ineq_tmp(
                           X_k_1, U_k_1, DW_1, jacob
                          ) = const_ineq(X_k_1, U_k_1) + jacob * DW_1
            obj += calculate_ineq_penalty(scvx, const_ineq_tmp,
                                          X_k_1, U_k_1, DW_1, jacob)
        end
    end
    if scvx.consts_initial_eq == []
    else
        for const_eq in scvx.consts_initial_eq
            (
             jacob, jacob_dict
            ) = use_or_calculate_jacobian(name_eq, size(jacob_dict[name_eq])[1] == 0, jacob_dict, const_eq, X_k_1, U_k_1)
            const_eq_tmp(
                         X_k_1, U_k_1, DW_1, jacob
                        ) = (
                             const_eq(X_k_1, U_k_1)
                             + jacob * reshape(DW_1, prod(size(DW_1)), 1)
                            )  # reshape for Cvx.Variable
            obj += calculate_eq_penalty(scvx, const_eq_tmp,
                                        X_k_1, U_k_1, DW_1, jacob)
        end
    end
    return obj, jacob_dict
end

function calculate_const_terminal_penalty_linearised(scvx::SCvx,
                                                     D, W, jacob_dict)
    name_ineq = "consts_terminal_ineq"
    if !haskey(jacob_dict, name_ineq)
        jacob_dict[name_ineq] = []
    end
    name_eq = "consts_terminal_eq"
    if !haskey(jacob_dict, name_eq)
        jacob_dict[name_eq] = []
    end
    obj = 0.0
    X_k, U_k = scvx.X_k, scvx.U_k  # linearisation point
    X_k_N = X_k[end, :]
    D_N = D[end, :]
    if scvx.consts_terminal_ineq == []
    else
        for const_ineq in scvx.consts_terminal_ineq
            (
             jacob, jacob_dict
            ) = use_or_calculate_jacobian(name_ineq, size(jacob_dict[name_ineq])[1] == 0, jacob_dict, const_ineq, X_k_N)
            const_ineq_tmp(
                           X_k_N, D_N, jacob
                          ) = (
                               const_ineq(X_k_N)
                               + jacob * reshape(D_N, prod(size(D_N)), 1)
                              )  # reshape for Cvx.Variable
            obj += calculate_ineq_penalty(scvx, const_ineq_tmp,
                                          X_k_N, D_N, jacob)
        end
    end
    if scvx.consts_terminal_eq == []
    else
        for const_eq in scvx.consts_terminal_eq
            (
             jacob, jacob_dict
            ) = use_or_calculate_jacobian(name_eq, size(jacob_dict[name_eq])[1] == 0, jacob_dict, const_eq, X_k_N)
            const_eq_tmp(
                         X_k_N, D_N, jacob
                        ) = (
                             const_eq(X_k_N)
                             + jacob * reshape(D_N, prod(size(D_N)), 1)
                            )  # reshape for Cvx.Variable
            obj += calculate_eq_penalty(scvx, const_eq_tmp,
                                        X_k_N, D_N, jacob)
        end
    end
    return obj, jacob_dict
end

################## etc ##################
function calculate_ineq_penalty(scvx::SCvx, func, args...)
    return scvx.λ * sum(pos(func(args...)))
end

function calculate_eq_penalty(scvx::SCvx, func, args...)
    return scvx.λ * norm(func(args...), 1)
end

function use_or_calculate_jacobian(name, condition, jacob_dict, func,
                                   args...; n=1)
    # TODO: add convex check
    if condition
        jacob = get_jacobian(func, args...)
        jacob_dict[name] = [jacob]
    else
        jacob = jacob_dict[name][n]
    end
    return jacob, jacob_dict
end

end


text:
test_gage_principal_axes(
    Array[
     [1., 0., 0.],
     [0., 1., 0.],
     [1., 1., 0.]],
    [sqrt(2)/2. -sqrt(2)/2. 0.;
     sqrt(2)/2.  sqrt(2)/2. 0.;
     0.          0.         1.], "test using array of 3 vectors as input")


text:
@testset "Project Rosalind" begin

    @testset "Project Rosalind: Count Nucleotides" begin
        @test count_nucleotides("AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC") == "20 21 17 12"
    end

    @testset "Project Rosalind: DNA to RNA" begin
        @test dna2rna("GATGGAACTTGACTACGTAAATT") == "GAUGGAACUUGACUACGUAAAUU"
    end

    @testset "Project Rosalind: Reverse Complement" begin
        @test reverse_complement("AAAACCCGGT") == "ACCGGGTTTT"
    end

end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia


text: eval(load_forwardbackward_discretization())


text:
#############################################################################
#############################################################################
#
# A script that runs all unit tests in the project.
#
#############################################################################
#############################################################################

include("../poly_factorization_project.jl")

println("Starting tests")
x = x_poly()
####
# Execute unit tests for integers
###
include("integers_test.jl")
test_euclid_ints()
test_ext_euclid_ints()

####
# Execute unit tests for polynomials
####
include("polynomials_test.jl")
prod_test_poly()
prod_test_poly_mod_p()
prod_derivative_test_poly()
prod_derivative_test_poly_mod_p()
ext_euclid_test_poly_mod_p()
division_test_poly_mod_p()
power_test_poly()
power_test_poly_modp()

####
# Execute unit tests for polynomial factorization
####
include("factorization_test.jl")
@time factor_test_poly_mod_p()

###################################
# Benchmark tests for polynomials #
###################################
include("polynomials_test.jl")
prod_benchmark()
prod_modP_benchmark()
crt_benchmark()


text:
@everywhere module SingleParticle
using StatsBase
using KernelDensity
using Distributions

export Particle, init_particle, init_particle2, bounds_check, bounds_check2, judge, MCstep, MCstep, MCstep_parallel, MCstep_parallel2

mutable struct Particle
    x_present::Array{Float64, 1}
    obj_func::Function # return relative weight
    obj_func_val::Float64
    weight::Float64
    cov::Array{Float64, 2}
end

function init_particle(lbs, ubs, obj_func, weight)
    """
    obj_func: must return -1 as a flag if has an error.
    """
    @assert length(lbs) == length(ubs)
    num_param = length(lbs)
    obj_func_val = -1
    kernel = Uniform(0, 10)
    x_present = zeros(length(lbs))   # generate parameter
    while obj_func_val < 0
        for i in 1:num_param
            x_present[i] = rand(kernel)
        end
        obj_func_val = obj_func(x_present)
    end
    return Particle(x_present, obj_func, obj_func_val, weight, eye(num_param))
end

function bounds_check(x_next)
    flag = true
    if x_next < 0 || x_next > 10
        flag = false
    end
    return flag
end

function judge(particle, proposed, F_obj_new, prior, ϵ)
    flag = false
    prob_present = prod([pdf(p, particle.x_present[i]) for (i, p) in enumerate(prior)])
    prob_proposed = prod([pdf(p, proposed[i]) for (i, p) in enumerate(prior)])
    MH_ratio = (F_obj_new < ϵ) * prob_proposed / ((particle.obj_func_val < ϵ) * prob_present)
    if rand() < MH_ratio
        flag = true
    else
        nothing
    end
    return flag
end

function MCstep(particle, prior, ϵ)
    """
    particle transition based on MH algorithm.
    """
    proposed = zeros(particle.x_present)
    kernel = MultivariateNormal(particle.x_present, particle.cov)
    F_obj_new = -1
    flag_MH = 0

    while F_obj_new < 0
        flag_bound = false
        while flag_bound == false
            proposed = rand(kernel)
            # proposed = np.random[:multivariate_normal](particle.x_present, np.array(particle.cov))
            for (i, x) in enumerate(proposed)
                flag_bound = bounds_check(x)
                if flag_bound == false
                    break
                end
            end
        end
        F_obj_new = particle.obj_func(proposed)
    end

    if judge(particle, proposed, F_obj_new, prior, ϵ) == true
        particle.x_present = proposed
        particle.obj_func_val = F_obj_new
        flag_MH = 1
    end

    return particle, flag_MH
end

function MCstep_parallel(particle, prior, ϵ)
    """
    particle transition based on MH algorithm.
    """
    paramset = particle.x_present
    obj_func_val = particle.obj_func_val
    proposed = zeros(particle.x_present)
    kernel = MultivariateNormal(particle.x_present, particle.cov)
    F_obj_new = -1
    flag_MH = 0

    while F_obj_new < 0
        flag_bound = false
        while flag_bound == false
            proposed = rand(kernel)
            # proposed = np.random[:multivariate_normal](particle.x_present, np.array(particle.cov))
            for (i, x) in enumerate(proposed)
                flag_bound = bounds_check(x)
                if flag_bound == false
                    break
                end
            end
        end
        F_obj_new = particle.obj_func(proposed)
    end

    if judge(particle, proposed, F_obj_new, prior, ϵ) == true
        paramset = proposed
        obj_func_val = F_obj_new
        flag_MH = 1
    end

    return paramset, obj_func_val, flag_MH
end

end


text:
#=------------------------------------------------------------------------------
   Helper Functions used throughout the package. This is a good place for code
   that shouldn't be exposed to the user.

  Operations Glossary
  -------------------
    * alterfilename - ()
    * reduce_dictionaries - ()
    * multiplicity_factor - ()
    * add_perm! - ()
    * string_as_varname_function - ()
    * matrix_to_dictionary - ()
    * reduce_edges - ()
    * parse - ()
    * redo_indexing - ()
    * remap_indices - ()
    *

------------------------------------------------------------------------------=#
"""-----------------------------------------------------------------------------
-----------------------------------------------------------------------------"""
function alterFilename(s::String,insert::String;
                       delimiter::String=".", keep_postfix::Bool=true)
    suffix,postfix = split(s,delimiter)
    if keep_postfix
        return suffix * insert * delimiter * keep_postfix
    else
        return suffix * insert
    end
end

"""-----------------------------------------------------------------------------
    reduce_dictionaries!(D1,D2)

This function takes in two dictionaries assumed to have the same value type and
reduces the contents of the second dicionary into the first. When two keys exist
in both dictionaries, their values are summed together.
-----------------------------------------------------------------------------"""
function reduce_dictionaries!(D1::Dict{Array{Int,1},N},
                              D2::Dict{Array{Int,1},N}) where N <: AbstractFloat
    for (key,val) in D2
        if haskey(D1,key)
            D1[key] += val
        else
            D1[key] = val
        end
    end
end

"""-----------------------------------------------------------------------------
    multiplicity_factor(indices)

  This function takes in a list of indices and returns the multinomial
coefficient computed by the frequency of the values in the indices. Works as a
helper function for computing the number of non-zeros the edge represents in the
vector contraction routines.

Input:
------
* indices -(Array{Int,1}):

  The indices associated with the hyper edge.

Output:
-------
* multinomial_factor - (Int)

   The number of non-zeros this edge represents in the original tensor.

Note
----

TODO: Change to include tuples too.
-----------------------------------------------------------------------------"""
function multiplicity_factor(indices::Array{Int,1},multiplicities::Dict{Int,Int},
                             final_counts::Array{Int,1})

    #reset multiplicities dictionary
    for key in keys(multiplicities)
        delete!(multiplicities,key)
    end

    #reset final_counts
    for i in 1:length(indices)
        final_counts[i] = 0
    end

    for index in indices

        if haskey(multiplicities,index)
            multiplicities[index] += 1
        else
            multiplicities[index] = 1
        end
    end

    #copy into format that can be passed to multinomial
    i = 1
    for (_,val) in multiplicities
        final_counts[i] = val
        i += 1
    end

    return multinomial(final_counts...)
end

function multiplicity_factor(indices::Array{Int,1})

    multiplicities = Dict{Int,Int}()

    for index in indices

        if haskey(multiplicities,index)
            multiplicities[index] += 1
        else
            multiplicities[index] = 1
        end
    end

    #copy into format that can be passed to multinomial
    final_counts = zeros(Int,length(indices))
    i = 1
    for (_,val) in multiplicities
        final_counts[i] = Int(val)
        i += 1
    end

    return multinomial(final_counts...)
end

"""-----------------------------------------------------------------------------
    nnz_count(indices)

  Computes the number of non-zeros found in the symmetric tensor.

Input:
------
* indices -(Array{Int,2}):

  An array with all the indices of the hyper edges found in the tensor.

Output:
-------
* nnz_count - (Int)

   The number of non-zeros in the original tensor.

-----------------------------------------------------------------------------"""
function nnz_count(indices::Array{Int,2})

    nnz_coutn = 0
    edge_count, _  = size(indices)
    for i =1:edge_count
        nnz_count += multiplicity_factor(indices[i,:])
    end

    nnz_count
end

"""
helper function for testing and adding in the permutation as an edge into the
dictionary used to build the symmetric tensor from the

TODO: add in types
"""
function add_perm!(edges,A,order,indices)
    val = A[indices...]
    sum = 0.0

    tol = 1e-12
    for p in unique(permutations(indices))
        sum += A[p...]
        if abs(A[p...] - val) > tol
            error(string("tensor is not symmetric: index ",p," has different values"))
        end
    end
    #add in the average of the values of the permutations of the tensor
    if sum > 0.0
        edges[collect(indices)] = sum/multiplicity_factor(collect(indices))
    end
end

#creates a new variable, bounds the last loop, used in meta code generators
function string_as_varname_function(s::AbstractString, v::Any)
   s = Symbol(s)
   @eval (($s) = ($v))
end

"""-----------------------------------------------------------------------------
    matrix_to_dictionary()

  This function takes in a matrix containing the indices of hyperedges as
columns and converts it into a dictionary linking the sorted indices to values.
If no values are passed in, then the weights of each hyperedge are assumed to be
1.0.

Inputs:
-------
* A - (k x n Array{Int,2}):

    Each column corresponds to the indices of a hyper edges.

Output:
-------
* D - (Dictionary{Array{Int,1},Float64}):

    The dictionary used to create the super symmetric tensor.

Note:
-----

Add in weights to link to hyper edges.
-----------------------------------------------------------------------------"""
function matrix_to_dictionary(A::Array{Int,2})
    D = Dict{Array{Int,1},Float64}()
    _,n = size(A)

    for j in 1:n
        sorted_indices = sort(A[:,j])
        if !haskey(D,sorted_indices)
            D[sorted_indices] = 1.0
        end
    end

    return D
end

"""-----------------------------------------------------------------------------
    reduce_edges(edges)

  This function takes in the edges to passed to create a new SSSTensor and reduces
the list into a dictionary where the edges that have the same edges are added
together.

Input:
------
* edges -(Array{{Array{Int,1},Number},1}):

    2 Tuples which contain the sorted indices and an edge weight associated with
    it.

Output:
-------
* edge_dict - (Dict{Array{Int,1},Number}):

    The resulting dictionary which has the edges aggregated together.
-----------------------------------------------------------------------------"""
function reduce_edges(edges::Array{Tuple{Array{Int,1},N},1}) where N <: AbstractFloat
    edge_dict = Dict{Array{Int,1},N}()

    for (indices, weight) in edges
        if haskey(edge_dict,indices)
            edge_dict[indices] += weight
        else
            edge_dict[indices] = weight
        end
    end

    return edge_dict
end

function reduce_edges(edges::Dict{Array{Int,1},N}) where N <: AbstractFloat
    return edges
end

function reduce_edges!(edge_dict::Dict{Array{Int,1},N},
                       edges::Array{Tuple{Array{Int,1},N},1}) where N <: AbstractFloat
    for (indices, weight) in edges
        if haskey(edge_dict,indices)
            edge_dict[indices] += weight
        else
            edge_dict[indices] = weight
        end
    end
    return edge_dict
end

"""-----------------------------------------------------------------------------
    redo_indexing!(hyperedges)

  rewrites the indices in the hyper edges to index from 1 instead of 0.
-----------------------------------------------------------------------------"""
function redo_indexing!(hyperedges::Array{Tuple{Array{Int,1},N},1}) where N <: AbstractFloat
    order = length(hyperedges[1][1])

    for i = 1:length(hyperedges)
        for j = 1:order
            hyperedges[i][1][j] += 1
        end
    end
end

"""-----------------------------------------------------------------------------
    remap_indices!(hyperedges)

  rewrites the vertices indices in the
-----------------------------------------------------------------------------"""
function remap_indices!(hyperedges::Array{Tuple{Array{Int,1},N},1}) where N <: AbstractFloat
    order = length(hyperedges[1][1])

    new_id = 1
    remapping_dict = Dict{Int,Int}()

    for i =1:length(hyperedges)

        for j =1:order
            if !haskey(remapping_dict,hyperedges[i][1][j])
                remapping_dict[hyperedges[i][1][j]] = new_id
                new_id += 1
            end
            hyperedges[i][1][j] = remapping_dict[hyperedges[i][1][j]]
        end
        sort!(hyperedges[i][1])
    end
end


text:
module PlantHydraulics

using CLIMAParameters
using CLIMAParameters.Planet: T_freeze, grav, ρ_cloud_liq
using ConstrainedRootSolvers: NewtonBisectionMethod, ReduceStepMethodND,
            ResidualTolerance, SolutionTolerance, SolutionToleranceND,
            find_peak, find_zero
using DocStringExtensions: TYPEDEF, TYPEDFIELDS
using Parameters
using Statistics
using WaterPhysics

# define global constants
struct EarthParameterSet <: AbstractEarthParameterSet end
const EARTH = EarthParameterSet()
GRAVITY(FT) = FT( grav(EARTH) );
K_25(FT)    = FT( T_freeze(EARTH) ) + 25;
R_GAS(FT)   = FT( gas_constant() );
ρ_H₂O(FT)   = FT( ρ_cloud_liq(EARTH) );
ρg_MPa(FT)  = ρ_H₂O(FT) * GRAVITY(FT) * FT(1e-6);

# export public types --- soil vulnerability
export AbstractSoilVC,
       BrooksCorey,
       VanGenuchten

# export public types --- xylem vulnerability
export AbstractXylemVC,
       WeibullDual,
       WeibullSingle

# export public types --- pressure volume curve
export AbstractCapacity,
       PVCurveLinear,
       PVCurveSegmented

# export public types --- flow mode
export AbstractFlowMode,
       NonSteadyStateMode,
       SteadyStateMode

# export public types --- hydraulic tissue
export AbstractHydraulicOrgan,
       LeafHydraulics,
       RootHydraulics,
       StemHydraulics

# export public types --- hydraulic system
export AbstractPlantOrganism,
       GrassLikeOrganism,
       PalmLikeOrganism,
       TreeLikeOrganism,
       TreeSimple

# export public functions --- initialize plant
export create_grass,
       create_palm,
       create_soil_VC,
       create_tree,
       fit_soil_VC!

# export public functions --- curves related
export p_from_volume,
       soil_erwc,
       soil_k_ratio_erwc,
       soil_k_ratio_p25,
       soil_k_ratio_rwc,
       soil_k_ratio_swc,
       soil_p_25_erwc,
       soil_p_25_rwc,
       soil_p_25_swc,
       soil_rwc,
       soil_swc,
       xylem_k_ratio,
       xylem_p_crit

# export public functions
export pressure_profile!,
       inititialize_legacy!,
       critical_flow,
       xylem_risk,
       plant_conductances!,
       roots_flow!,
       xylem_flow,
       update_PVF!,
       temperature_effects!,
       end_pressure,
       fit_xylem_VC

include("types/curves.jl")
include("types/flow.jl"  )
include("types/organ.jl" )
include("types/plant.jl" )

include("initialize/legacy.jl")
include("initialize/soil.jl"  )
include("initialize/plant.jl" )

include("curves/capacity.jl")
include("curves/soil.jl"    )
include("curves/xylem.jl"   )

include("hydraulics/conductance.jl")
include("hydraulics/flow.jl"       )
include("hydraulics/pressure.jl"   )
include("hydraulics/temperature.jl")

include("hydraulics/capacitance.jl")

end # module

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
"""
Generate constant signal.

# Parameters:
- `k`: Constant output value
"""
function Constant(;name, k=1)
    @named output = RealOutput()
    pars = @parameters k=k
    eqs = [
        output.u ~ k
    ]
    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate sine signal.

# Parameters:
- `frequency`: [Hz] Frequency of sine wave
- `amplitude`: Amplitude of sine wave
- `phase`: [rad] Phase of sine wave
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
"""
function Sine(;name,
    frequency,
    amplitude=1,
    phase=0,
    offset=0,
    start_time=0)

    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
    eqs = [
        output.u ~ offset + ifelse(t < start_time, 0, amplitude* sin(2*pi*frequency*(t - start_time) + phase))
    ]
    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate cosine signal.

# Parameters:
- `frequency`: [Hz] Frequency of sine wave
- `amplitude`: Amplitude of sine wave
- `phase`: [rad] Phase of sine wave
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
"""
function Cosine(;name,
    frequency,
    amplitude=1,
    phase=0,
    offset=0,
    start_time=0)

    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase
    eqs = [
        output.u ~ offset + ifelse(t < start_time, 0, amplitude* cos(2*pi*frequency*(t - start_time) + phase))
    ]
    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate current time signal.

# Parameters:
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
"""
function ContinuousClock(;name, offset=0, start_time=0)
    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time
    eqs = [
        output.u ~ offset + ifelse(t < start_time, 0, t - start_time)
    ]
    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate ramp signal.

# Parameters:
- `height`: Height of ramp
- `duration`: [s] Duration of ramp (= 0.0 gives a Step)
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
"""
function Ramp(;name,
    offset=0,
    height=1,
    duration=1,
    start_time=0)

    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time height=height duration=duration
    eqs = [
        output.u ~ offset + ifelse(t < start_time, 0,
            ifelse(t < (start_time + duration), (t - start_time) * height / duration, height))
    ]
    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate step signal.

# Parameters:
- `height`: Height of step
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
"""
function Step(;name, offset=0, height=1, start_time=0)
    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time height=height
    eqs = [
        output.u ~ offset + ifelse(t < start_time, 0, height)
    ]
    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

"""
Generate exponentially damped sine signal.

# Parameters:
- `frequency`: [Hz] Frequency of sine wave
- `amplitude`: Amplitude of sine wave
- `damping`: [1/s] Damping coefficient of sine wave
- `phase`: [rad] Phase of sine wave
- `offset`: Offset of output signal
- `start_time`: [s] Output `y = offset` for `t < start_time`
"""
function ExpSine(;name,
    frequency,
    amplitude=1,
    damping=0.1,
    phase=0,
    offset=0,
    start_time=0)

    @named output = RealOutput()
    pars = @parameters offset=offset start_time=start_time amplitude=amplitude frequency=frequency phase=phase damping=damping
    eqs = [
        output.u ~ offset + ifelse(t < start_time, 0, amplitude * exp(-damping * (t - start_time)) * sin(2*pi*frequency*(t - start_time) + phase))
    ]
    compose(ODESystem(eqs, t, [], pars; name=name), [output])
end

# TODO:
# - Exponentials    Generate a rising and falling exponential signal
# - Pulse   Generate pulse signal of type Real
# - SawTooth    Generate saw tooth signal
# - Trapezoid   Generate trapezoidal signal of type Real


text:
# Tests for empty unions of intervals
# I = ∅
I = IntervalUnion()
# J = [0,1] ∪ ]3,5]
J = IntervalUnion([Interval(0,1),Interval(3,true,5)])
@test empty(I)
@test number_of_components(I) == 0
@test string(I) == "∅"
@test cardinal(I) == 0
@test left(I) == []
@test right(I) == []
@test limits(I) == []
@test compact(I) == IntervalUnion()
@test compactness(I) == 0
@test 0 ∉ I
@test 1.2 ∉ I
@test I ∩ J == IntervalUnion()
@test J ∩ I == IntervalUnion()
@test I ∪ J == J
@test J ∪ I == J
@test I ⊆ J
@test J ⊈ I
@test I ⊆ I
@test complement(I) == IntervalUnion([Interval(-Inf,true,Inf,true)])
@test setdiff(I,J) == IntervalUnion()
@test setdiff(J,I) == J
@test symdiff(I,J) == J
@test symdiff(J,I) == J
@test jaccard(I,J) == 0
@test jaccard(J,I) == 0
@test overlap_coefficient(I,J) == 0
@test overlap_coefficient(J,I) == 0
@test dice_coefficient(I,J) == 0
@test dice_coefficient(J,I) == 0

# I = [0,1]
I = IntervalUnion([Interval(0,1)])
II = IntervalUnion(0,1)
@test I == II
@test string(I) == "[0,1]"
@test string(II) == "[0,1]"

# J = ]1,2.5]
J = IntervalUnion([Interval(1,true,2.5)])
JJ = IntervalUnion(1,true,2.5)
@test J == JJ
@test string(J) == "]1,2.5]"
@test string(JJ) == "]1,2.5]"

# K = ]0,1[ ∪ ]1,2.6] ∪ [2.7,4]
K = IntervalUnion([Interval(0,true,1,true),Interval(1,true,2.6),Interval(2.7,4)])
@test string(K) == "]0,1[ ∪ ]1,2.6] ∪ [2.7,4]"

# L = [-1,-0.998[
L = IntervalUnion([Interval(-1,-0.998,true)])
LL = IntervalUnion(-1,-0.998,true)
@test L == LL
@test string(L) == "[-1,-0.998["
@test string(LL) == "[-1,-0.998["

# L = ]-2.3,1.67[
M = IntervalUnion([Interval(-2.3,true,1.67,true)])
MM = IntervalUnion(-2.3,true,1.67,true)
@test M == MM
@test string(M) == "]-2.3,1.67["
@test string(MM) == "]-2.3,1.67["

# Tests for equality
@test I == I
@test I != J
@test K != J

# Tests for number of components
@test number_of_components(I) == 1
@test number_of_components(J) == 1
@test number_of_components(K) == 3

# Tests for limits
@test left(I) == [0]
@test right(I) == [1]
@test left(J) == right(I)
@test left(K) == [0,1,2.7]
@test right(K) == [1,2.6,4]
@test limits(I) == [0,1]
@test limits(J) == [1,2.5]
@test limits(K) == [0,1,2.6,2.7,4]

# Tests for ∈
@test 0.55 ∈ I
@test 0 ∈ I
@test 1 ∉ J
@test 1.000000001 ∈ J
@test 1 ∉ K
@test 2.77 ∈ K

# Tests for cardinal
@test cardinal(I) == 1.0
@test cardinal(J) == 1.5
@test cardinal(K) == 3.9

# Tests for ⊆
@test I ⊈ K
@test J ⊆ K
@test I ⊆ I
@test K ⊆ K

# Tests for ∪
# [0,1] ∪ ]1,2.5] = [0,2.5]
@test I ∪ J == IntervalUnion([Interval(0,2.5)])
# [0,1] ∪ ]0,1[ ∪ ]1,2.6] ∪ [2.7,4] = [0,2.6] ∪ [2.7,4]
@test I ∪ K == IntervalUnion([Interval(0,2.6),Interval(2.7,4)])
# ]1,2.5] ∪ ]0,1[ ∪ ]1,2.6] ∪ [2.7,4] = ]0,1[ ∪ ]1,2.6] ∪ [2.7,4]
@test J ∪ K == K
@test I ∪ I == I
@test J ∪ J == J
@test K ∪ K == K

# Tests for ∩
# [0,1] ∩ ]1,2.5] = ]0,0[
@test I ∩ J == IntervalUnion()
# [0,1] ∩ (]0,1[ ∪ ]1,2.6] ∪ [2.7,4]) = ]0,1[
@test I ∩ K == IntervalUnion([Interval(0,true,1,true)])
# ]1,2.5] ∩ (]0,1[ ∪ ]1,2.6] ∪ [2.7,4]) = ]1,2.5]
@test J ∩ K == J
@test I ∩ I == I
@test J ∩ J == J
@test K ∩ K == K

# Tests for complement
# comp([0,1]) = ]-inf,0[ ∪ ]1,inf[
@test complement(I) == IntervalUnion([Interval(-Inf,true,0,true),Interval(1,true,Inf,true)])
# comp([1,2.5]) = ]-inf,1] ∪ ]2.5,inf[
@test complement(J) == IntervalUnion([Interval(-Inf,true,1),Interval(2.5,true,Inf,true)])
# comp(]0,1[ ∪ ]1,2.6] ∪ [2.7,4]) = ]-inf,0] ∪ [1,1] ∪ ]2.6,2.7[ ∪ ]4,inf[
@test complement(K) == IntervalUnion([Interval(-Inf,true,0),Interval(1,1),Interval(2.6,true,2.7,true),Interval(4,true,Inf,true)])
# comp([1,1]) = ]-inf,1[ ∪ ]1,inf[
@test complement(IntervalUnion([Interval(1,1)])) == IntervalUnion([Interval(-Inf,true,1,true),Interval(1,true,Inf,true)])
@test complement(complement(I)) == I
@test complement(complement(J)) == J
@test complement(complement(K)) == K
# complement(]-Inf,Inf[) = ∅
@test complement(IntervalUnion([Interval(-Inf,true,Inf,true)])) == IntervalUnion()
# complement(]-Inf,2]) = ]2,Inf[
@test complement(IntervalUnion([Interval(-Inf,true,2)])) == IntervalUnion([Interval(2,true,Inf,true)])
# complement(]-Inf,-1.23[) = [-1.23,Inf[
@test complement(IntervalUnion([Interval(-Inf,true,-1.23,true)])) == IntervalUnion([Interval(-1.23,Inf,true)])
# complement(]1,Inf[) = ]-Inf,1]
@test complement(IntervalUnion([Interval(1,true,Inf,true)])) == IntervalUnion([Interval(-Inf,true,1)])
# complement([-2.76,Inf[) = ]-Inf,-2.76[
@test complement(IntervalUnion([Interval(-2.76,Inf,true)])) == IntervalUnion([Interval(-Inf,true,-2.76,true)])

# Tests for setdiff
# setdiff([0,2] ∪ [3,4], [1,3] ∪ [3.5,5]) = [0,1[ ∪ ]3,3.5[
@test setdiff(IntervalUnion([Interval(0,2),Interval(3,4)]),IntervalUnion([Interval(1,3),Interval(3.5,5)]))==IntervalUnion([Interval(0,1,true),Interval(3,true,3.5,true)])
# setdiff([0,1], ]0,1[ ∪ ]1,2.6] ∪ [2.7,4]) = [0,0] ∪ [1,1]
@test setdiff(I,K) == IntervalUnion([Interval(0,0),Interval(1,1)])

# Tests for symdiff
# symdiff([0,1], ]0,1[ ∪ ]1,2.6] ∪ [2.7,4]) =[0,0] ∪ [1,2.6] ∪ [2.7,4]
@test symdiff(I,K) == IntervalUnion([Interval(0,0),Interval(1,2.6),Interval(2.7,4)])

# Tests for compact
# compact(]0,1] ∪ [2,3])=[0,3]
@test compact(IntervalUnion([Interval(0,true,1),Interval(2,3)]))==IntervalUnion([Interval(0,3)])

# Tests for compactness
# compactness(]0,1] ∪ [2,3])=2/3
@test compactness(IntervalUnion([Interval(0,true,1),Interval(2,3)]))==2/3

# Tests for superset
# superset(]0,1] ∪ [2,3])=]0,3]
@test superset(IntervalUnion([Interval(0,true,1),Interval(2,3)]))==IntervalUnion([Interval(0,true,3)])

# Tests for super_complement
# super_complement(]0,1] ∪ [2,3]) = ]-Inf,0] ∪ ]3,Inf[
@test super_complement(IntervalUnion([Interval(0,true,1),Interval(2,3)]))==IntervalUnion([Interval(-Inf,true,0),Interval(3,true,Inf,true)])

# Tests for restricted_complement
# restricted_complement([0,1] ∪ [4,5]) = ]1,4[
@test restricted_complement(IntervalUnion([Interval(0,1),Interval(4,5)]))==IntervalUnion([Interval(1,true,4,true)])
# restricted_complement(]0,1[ ∪ ]2.5,3] ∪ [3.5,6]) = [1,2.5] ∪ ]3,3.5[
@test restricted_complement(IntervalUnion([Interval(0,true,1,true),Interval(2.5,true,3),Interval(3.5,6)]))==IntervalUnion([Interval(1,2.5),Interval(3,true,3.5,true)])

# Tests for complement_cardinal
# complement_cardinal([0,1] ∪ [4,5]) = 3
@test complement_cardinal(IntervalUnion([Interval(0,1),Interval(4,5)]))==3
# complement_cardinal(]0,1[ ∪ ]2.5,3] ∪ [3.5,6]) = 2
@test complement_cardinal(IntervalUnion([Interval(0,true,1,true),Interval(2.5,true,3),Interval(3.5,6)]))==2

# Tests for sampling unions of intervals
# I = ]-2.3,0[ ∪ ]0,1]
I = IntervalUnion([Interval(-2.3,true,0,true),Interval(0,true,1)])
s = sample(I)
@test (-2.3 < s < 0) || (0 < s <= 1)
samples = sample(I,100)
@test length(samples) == 100
for s in samples
	@test (-2.3 < s < 0) || (0 < s <= 1)
end
# I = [1,1] ∪ [3.4,3.4] ∪ [6,6]
I = IntervalUnion([Interval(1,1),Interval(3.4,3.4),Interval(6,6)])
s = sample(I)
@test s ∈ [1,3.4,6]