text:
# create a dictionary with all the possible models
function get_car_models(env::UrbanEnv, get_model::Function)
    d = Dict{SVector{2, LaneTag}, DriverModel}()

    #
    r1 = SVector(LaneTag(1,1), LaneTag(2,1))
    d[r1] = get_model(env, r1)

    r2 = SVector(LaneTag(1,1), LaneTag(5,1))
    d[r2] = get_model(env, r2)

    r3 = SVector(LaneTag(3,1), LaneTag(4,1))
    d[r3] = get_model(env, r3)

    r4 = SVector(LaneTag(3,1), LaneTag(5,1))
    d[r4] = get_model(env, r4)

    return d
end

function get_stop_model(env::UrbanEnv, r::SVector{2, LaneTag})
    route = [env.roadway[tag] for tag in find_route(env, r)]
    intersection_entrances = get_start_lanes(env.roadway)
    if !(route[1] ∈ intersection_entrances)
        intersection = Lane{Float64}[]
        intersection_exits = Lane{Float64}[]
    else
        intersection_exits = get_exit_lanes(env.roadway)
        intersection=Lane{Float64}[route[1], route[2]]
    end
    navigator = RouteFollowingIDM(route=route, a_max=3., σ=1.)
    intersection_driver = StopIntersectionDriver(navigator= navigator,
                                                intersection=intersection,
                                                intersection_entrances = intersection_entrances,
                                                intersection_exits = intersection_exits,
                                                stop_delta=maximum(env.params.crosswalk_width),
                                                accel_tol=0.,
                                                priorities = env.priorities)
    crosswalk_drivers = Vector{CrosswalkDriver}(length(env.crosswalks))
    for i=1:length(env.crosswalks)
        cw_conflict_lanes = get_conflict_lanes(env.crosswalks[i], env.roadway)
        crosswalk_drivers[i] = CrosswalkDriver(navigator = navigator,
                                crosswalk = env.crosswalks[i],
                                conflict_lanes = cw_conflict_lanes,
                                intersection_entrances = intersection_entrances,
                                yield=!isempty(intersect(cw_conflict_lanes, route))
                                )
        # println(" yield to cw ", i, " ", crosswalk_drivers[i].yield)
    end
    return UrbanDriver(navigator=navigator,
                    intersection_driver=intersection_driver,
                    crosswalk_drivers=crosswalk_drivers
                    )
end

function get_ttc_model(env::UrbanEnv, r::SVector{2, LaneTag})
    route = [env.roadway[tag] for tag in find_route(env, r)]
    intersection_entrances = get_start_lanes(env.roadway)
    if !(route[1] ∈ intersection_entrances)
        intersection = Lane{Float64}[]
        intersection_exits = Lane{Float64}[]
    else
        intersection_exits = get_exit_lanes(env.roadway)
        intersection=Lane{Float64}[route[1], route[2]]
    end
    navigator = RouteFollowingIDM(route=route, a_max=2., σ=1.)
    intersection_driver = TTCIntersectionDriver(navigator = navigator,
                                            intersection = intersection,
                                            intersection_pos = VecSE2(mdp.env.params.inter_x,
                                                                        mdp.env.params.inter_y),
                                            stop_delta = maximum(mdp.env.params.crosswalk_width),
                                            accel_tol = 0.,
                                            priorities = mdp.env.priorities,
                                            ttc_threshold = (mdp.env.params.x_max - mdp.env.params.inter_x)/mdp.env.params.speed_limit
                                            )
    crosswalk_drivers = Vector{CrosswalkDriver}(length(env.crosswalks))
    for i=1:length(env.crosswalks)
        cw_conflict_lanes = get_conflict_lanes(env.crosswalks[i], env.roadway)
        crosswalk_drivers[i] = CrosswalkDriver(navigator = navigator,
                                crosswalk = env.crosswalks[i],
                                conflict_lanes = cw_conflict_lanes,
                                intersection_entrances = intersection_entrances,
                                yield=!isempty(intersect(cw_conflict_lanes, route))
                                )
        # println(" yield to cw ", i, " ", crosswalk_drivers[i].yield)
    end
    return UrbanDriver(navigator=navigator,
                    intersection_driver=intersection_driver,
                    crosswalk_drivers=crosswalk_drivers
                    )
end

function AutomotiveDrivingModels.observe!(model::UrbanDriver,
                                          ego::VehicleState,
                                          car::VehicleState,
                                          ped::VehicleState,
                                          roadway::Roadway)
    AutomotiveDrivingModels.observe!(model.navigator, ego, car, roadway)

    for driver in model.crosswalk_drivers
        AutomotiveDrivingModels.observe!(driver, ego, car, ped, roadway)
    end
    a_lon_crosswalks = minimum([driver.a.a_lon for driver in model.crosswalk_drivers])
    a_lon = min(a_lon_crosswalks, model.navigator.a)

    model.a = LonAccelDirection(a_lon, model.navigator.dir)
    return model.a
end

function AutomotiveDrivingModels.observe!(model::TTCIntersectionDriver,
                                          ego::VehicleState,
                                          car::VehicleState,
                                          roadway::Roadway)

    model.priority = ttc_check(model, car, ego, roadway)

end

function AutomotiveDrivingModels.observe!(model::CrosswalkDriver,
                                          ego::VehicleState,
                                          car::VehicleState,
                                          ped::VehicleState,
                                          roadway::Roadway)
    a_lon =0.
    a_lon_idm = model.navigator.a_max
    if !model.yield
        a_lon = a_lon_idm
    else
        model.priority = !is_crossing(ped, model.crosswalk, model.conflict_lanes, roadway) ||
        has_passed(model, car, roadway)
        if !model.priority
            dist_to_cw = get_distance_to_crosswalk(model, car, roadway, -model.stop_delta)
            a_lon = min(a_lon_idm, stop_at_dist(model, car, dist_to_cw))
            # println("dist to cw $dist_to_cw, crosswalk $(model.crosswalk.tag), $(get_lane(roadway, car) ∈ model.intersection_entrances)")
            # println("priority $(model.priority), crossing $(is_crossing(ped, model.crosswalk, model.conflict_lanes, roadway)), a_lon $a_lon, v $(car.v)")
        else
            a_lon = a_lon_idm # just idm
        end
    end
    model.a = LonAccelDirection(a_lon, model.navigator.dir)
end

function AutomotiveDrivingModels.observe!(model::RouteFollowingIDM,
                                          ego::VehicleState,
                                          car::VehicleState,
                                          roadway::Roadway)
    Δv = model.v_des - car.v
    acc = Δv*model.k_spd
    model.a = clamp(acc, -model.d_max, model.a_max)
    # fore = is_neighbor_fore_along_lane(car, ego, roadway)
    # v_car = car.v
    # if fore.ind != nothing
    #     headway, v_oth = fore.Δs, car.v
    # else
    #     headway, v_oth = NaN, NaN
    # end
    # track_longitudinal!(model, v_car, v_oth, headway)
    cur_lane = get_lane(roadway, car)
    set_direction!(model, cur_lane, roadway)
    return model
end

function is_crossing(ped::VehicleState, crosswalk::Lane{T}, conflict_lanes::Vector{Lane{T}}, roadway::Roadway) where T<:Real
    ped_lane = get_lane(roadway, ped)
    if ped_lane.tag != crosswalk.tag
        return false
    end
    if 4.0 <= ped.posF.s <= 10.
        return true
    end
    # at this point, the pedestrian is not on the road
    # check if the pedestrian is going to cross or not
    if AutomotivePOMDPs.direction_from_center(ped, crosswalk) > 0. && get_lane(roadway, ped).tag == crosswalk.tag
        return true
    end
    return false
end

function ttc_check(model::TTCIntersectionDriver, ego::VehicleState, car::VehicleState, roadway::Roadway)
    posF = car.state.posF
    int_x, int_y, int_θ = model.intersection_pos
    lane = get_lane(roadway, car)
    int_proj = Frenet(model.intersection_pos, lane, roadway)
    if normsquared(VecE2(model.intersection_pos - car.state.posG)) < inter_width^2 # vehicle is in the middle
        ttc = 0.
    else
        ttc = (int_proj.s - posF.s)/car.state.v
    end

    if 0 <= ttc < model.ttc_threshold
        return false
    else
        return true
    end

end

function is_neighbor_fore_along_lane(ego::VehicleState, car::VehicleState, roadway::Roadway, ego_def::VehicleDef = VehicleDef())
    tag_start = ego.posF.roadind.tag
    s_base = ego.posF.s + ego_def.length/2*cos(ego.posF.ϕ)
    tag_start = ego.posF.roadind.tag
    targetpoint_primary = VehicleTargetPointRear()
    targetpoint_valid = VehicleTargetPointFront()

    is_neighbor_fore_along_lane(ego, car, roadway, tag_start, s_base)
end

function is_neighbor_fore_along_lane(
    ego::VehicleState,
    veh::VehicleState,
    roadway::Roadway,
    tag_start::LaneTag,
    s_base::Float64;
    max_distance_fore::Float64 = 250.0, # max distance to search forward [m]
    )

    best_ind = nothing
    best_dist = max_distance_fore
    tag_target = tag_start
    lane = roadway[tag_target]

    dist_searched = 0.0
    while dist_searched < max_distance_fore

        s_adjust = NaN
        if veh.posF.roadind.tag == tag_target
            s_adjust = 0.0

        elseif is_between_segments_hi(veh.posF.roadind.ind, lane.curve) &&
                is_in_entrances(roadway[tag_target], veh.posF.roadind.tag)

            distance_between_lanes = norm(VecE2(roadway[tag_target].curve[1].pos - roadway[veh.posF.roadind.tag].curve[end].pos))
            s_adjust = -(roadway[veh.posF.roadind.tag].curve[end].s + distance_between_lanes)

        elseif is_between_segments_lo(veh.posF.roadind.ind) &&
                is_in_exits(roadway[tag_target], veh.posF.roadind.tag)

            distance_between_lanes = norm(VecE2(roadway[tag_target].curve[end].pos - roadway[veh.posF.roadind.tag].curve[1].pos))
            s_adjust = roadway[tag_target].curve[end].s + distance_between_lanes
        end

        if !isnan(s_adjust)
            s_valid = veh.posF.s + VehicleDef().length/2*cos(veh.posF.ϕ) + s_adjust
            dist_valid = s_valid - s_base + dist_searched
            if dist_valid ≥ 0.0
                s_primary = veh.posF.s - VehicleDef().length/2*cos(veh.posF.ϕ) + s_adjust
                dist = s_primary - s_base + dist_searched
                if dist < best_dist
                    best_dist = dist
                    best_ind = 1
                end
            end
        end
        if best_ind != nothing
            break
        end

        if !has_next(lane) ||
           (tag_target == tag_start && dist_searched != 0.0) # exit after visiting this lane a 2nd time
            break
        end

        dist_searched += (lane.curve[end].s - s_base)
        s_base = -norm(VecE2(lane.curve[end].pos - next_lane_point(lane, roadway).pos)) # negative distance between lanes
        tag_target = next_lane(lane, roadway).tag
    end
    NeighborLongitudinalResult(best_ind, best_dist)
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
mutable struct DefaultConstraint <: Constraint
    rhs::Float64;
    type::Symbol;
    variablesPositiveCoefficients::Dict{Int,Float64};
    variablesNegativeCoefficients::Dict{Int,Float64};
end

"""
    constructConstraint(rhs, t, vars, coefs)

Constructs a `DefaultConstraint` with right-hand side `rhs`, type `t` and left-hand side given by the variables in `vars` and coefficients in `coefs`.

**Note:** Both, the rhs and coefs must be of type Float64.

# Example
```jldoctest
julia> OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 12, 33], [2.0, 5.0, 3.0, 1.0]);
```
"""
function constructConstraint(rhs::Float64,
                             constraintType::Symbol,
                             variables::Array{Int,1},
                             coefficients::Array{Float64,1})
    dictPositive = Dict{Int,Float64}();
    dictNegative = Dict{Int,Float64}();
    @inbounds for i in 1:length(variables)
        if coefficients[i] > 0.0
            dictPositive[variables[i]] = coefficients[i];
        else
            dictNegative[variables[i]] = coefficients[i];
        end
    end
    return DefaultConstraint(rhs, constraintType, dictPositive, dictNegative);
end

"""
    get_rhs(c)

Returns the right-hand side value of constraint `c`

# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 12, 33], [2.0, 5.0, 3.0, 1.0]);
julia> OREnvironment.get_rhs(c)
12.0
```
"""
get_rhs(c::Constraint)::Float64 = c.rhs;

"""
    set_rhs!(c, val)

Sets the right-hand side value of constraint `c` to `val`

# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 12, 33], [2.0, 5.0, 3.0, 1.0]);
julia> OREnvironment.set_rhs!(c, 25.0);
julia> OREnvironment.get_rhs(c)
25.0
```
"""
set_rhs!(c::Constraint, value::Float64) = c.rhs = value;

"""
    get_type(c)

Returns the constraint type of `c`. Admisible values are :lessOrEq, :equal and :greaterOrEq.

# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 12, 33], [2.0, 5.0, 3.0, 1.0]);
julia> OREnvironment.get_type(c)
:lessOrEq
```
"""
get_type(c::Constraint)::Symbol = c.type;

"""
    set_type!(c, val)

Sets the constraint type of `c` to `val`. Admisible values are :lessOrEq, :equal and :greaterOrEq.

# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 12, 33], [2.0, 5.0, 3.0, 1.0]);
julia> OREnvironment.set_type!(c, :greaterOrEq);
julia> OREnvironment.get_type(c)
:greaterOrEq
```
"""
set_type!(c::Constraint, value::Symbol) = c.type = value;

"""
    is_variable(c, var)

Returns true if variable `var` is in constraint `c`.

# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 12, 33], [2.0, 5.0, 3.0, 1.0]);
julia> OREnvironment.is_variable(c, 33)
true
julia> OREnvironment.is_variable(c, 13)
false
```
"""
function is_variable(c::Constraint, variable::Int)::Bool
    local answer::Bool = haskey(c.variablesPositiveCoefficients, variable) ||
                         haskey(c.variablesNegativeCoefficients, variable);
    return answer;
end

"""
    get_coefficient(c, var)

Returns the coefficient of variable `var` in constraint `c`. If the variable is not in the constraint it returns 0.0;

# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 12, 33], [2.0, 5.0, 3.0, 1.0]);
julia> OREnvironment.get_coefficient(c, 33)
1.0
julia> OREnvironment.get_coefficient(c, 3)
5.0
```
"""
function get_coefficient(c::Constraint, variable::Int)::Float64
    local value::Float64 = get(c.variablesPositiveCoefficients, variable, 0.0);
    if value == 0.0
        value = get(c.variablesNegativeCoefficients, variable, 0.0);
    end
    return value;
end

"""
    set_coefficient!(c, var, coef)

Sets in constraint `c` coefficient `coef` for variable `var`.

**Note:** If the variable was not in the constraint, it is added to it. Note that the vector that relates variables and constraints in the `Problem` struct is not updated when calling this adding a variable using this method!

# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 12, 33], [2.0, 5.0, 3.0, 1.0]);
julia> OREnvironment.set_coefficient!(c, 33, 123.0);
julia> OREnvironment.get_coefficient(c, 33)
123.0
```
"""
function set_coefficient!(c::Constraint, variable::Int, value::Float64)
    if value > 0.0
        c.variablesPositiveCoefficients[variable] = value;
    else
        c.variablesNegativeCoefficients[variable] = value;
    end
    return nothing;
end

"""
    get_relationship_variables_constraints(vconstraints, nVariables)

Given an Array of constraints `vconstraints` and the number of variables in the problem, returns an Array where the i-th element is an Array with the indexes of the constraints where the i-th variable appears.

This function is employed in the `Problem` struct so it is possible to quickly identify what constraints in the problem are affected when a variable changes its value. This permits faster updates of the left-hand side value associated to a given solution.

# Example
```jldoctest
julia> c1 = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 4], [2.0, 5.0, 3.0]);
julia> c2 = OREnvironment.constructConstraint(5.0, :equal, [2, 3, 5], [9.0, 1.0, 7.0]);
julia> constraints = [c1, c2];
julia> nVariables = 5;
julia> OREnvironment.get_relationship_variables_constraints(constraints, nVariables)
5-element Array{Array{Int64,1},1}:
 [1]
 [2]
 [1, 2]
 [1]
 [2]
 # So variable 1 appears in constraint 1, variable 2 in constraint 2,
 # variable 3 in constraints 1 and 2...
```
"""
function get_relationship_variables_constraints(constraints::Array{<:Constraint,1},
                                                nVariables::Int)::Array{Array{Int,1}, 1}
    if length(constraints) == 0 return Array{Array{Int,1}, 1}() end
    variablesConstraints = [Int[] for i in 1:nVariables];
    @inbounds for i in 1:length(constraints)
      add_constraint_index_to_variables!(constraints[i], i, variablesConstraints);
    end
    return variablesConstraints;
end

function add_constraint_index_to_variables!(c::Constraint,
                                            idx::Int,
                                            variablesConstraints::Array{Array{Int,1}, 1})
    @inbounds for variable in keys(c.variablesPositiveCoefficients)
        push!(variablesConstraints[variable], idx);
    end
    @inbounds for variable in keys(c.variablesNegativeCoefficients)
        push!(variablesConstraints[variable], idx);
    end
    return nothing;
end

"""
    read_constraints(file)

Returns an Array of constraints. The constraints are read from `file`.

See example in the file with the tests: *./test/Constraints.jl*.
"""
function read_constraints(file::String)
    fileStream = open(file, "r");
    lines = readlines(fileStream);
    constraints = [constructConstraint(0.0, :equal, Int[], Array{Float64,1}())
                   for i in 1:length(lines)];
    @inbounds Threads.@threads for i in 1:length(lines)
      parse_constraint!(constraints[i], lines[i], i);
    end
    return constraints;
end

function parse_constraint!(c::Constraint, line::String, idx::Int)
    local sign = parse_and_set_constraint_sign!(c, line, idx);
    elements = split(line, sign);

    set_rhs!(c, parse(Float64, elements[2]));

    lhs = split(elements[1], "+");
    @inbounds for i in 1:length(lhs)
      terms = split(lhs[i], "x");
      coef = parse(Float64, terms[1]);
      var  = parse(Int, terms[2]);
      set_coefficient!(c, var, coef);
    end
end

function parse_and_set_constraint_sign!(c::Constraint, line::String, idx::Int)
    local sign = "";
    if occursin("<=", line)
      sign = "<=";
      set_type!(c, :lessOrEq);
    elseif occursin(">=", line)
      sign = ">=";
      set_type!(c, :greaterOrEq);
    elseif occursin("=", line) && !occursin("<", line) && !occursin(">", line) && !occursin("==", line)
      sign = "=";
      set_type!(c, :equal);
    else
      error("Incorrect sign in constraint ", idx, ". Maybe you wrote => instead of >=, or == instead of =");
    end
    return sign;
end

"""
    is_feasible(c, lhs)

Given a constraint `c`, returns if the left-hand side consumption `lhs` is feasible or not.

# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 4], [2.0, 5.0, 3.0]);
julia> OREnvironment.is_feasible(c, 6.0)
true
julia> OREnvironment.is_feasible(c, 16.0)
false
```
"""
is_feasible(c::Constraint, lhs::Float64)::Bool = is_feasible(get_rhs(c), lhs, Val(get_type(c)));
is_feasible(rhs::Float64, lhs::Float64, ::Val{:lessOrEq})::Bool    = lhs <= rhs;
is_feasible(rhs::Float64, lhs::Float64, ::Val{:equal})::Bool       = lhs ≈ rhs;
is_feasible(rhs::Float64, lhs::Float64, ::Val{:greaterOrEq})::Bool = lhs >= rhs;

"""
    is_active(c, lhs)

Given a constraint `c`, returns if the left-hand side consumption `lhs` makes the constraint active or not.

# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 4], [2.0, 5.0, 3.0]);
julia> OREnvironment.is_active(c, 6.0)
false
julia> OREnvironment.is_active(c, 12.0)
true
```
"""
is_active(c::Constraint, lhs::Float64)::Bool = lhs ≈ get_rhs(c);

"""
    compute_lhs_after_increment(var, Δ, currentLHS, c)

For a change of `Δ` in the value of variable `var`, returns the new left-hand side consumption in constraint `c` given the current consumption `currentLHS`.

# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 4], [2.0, 5.0, 3.0]);
julia> Δ = 3; currentLHS = 9.0;
julia> OREnvironment.compute_lhs_after_increment(1, Δ, currentLHS, c)
15.0 # = 2*Δ + currentLHS
```
"""
function compute_lhs_after_increment(variable::Int,
                                     Δvariable::Real,
                                     currentLHS::Float64,
                                     c::Constraint)::Float64
    local Δconsumption::Float64 = Δvariable*get_coefficient(c, variable);
    local lhs::Float64 = currentLHS + Δconsumption;
    return lhs;
end

"""
    is_increment_feasible(var, Δ, currentLHS, c)

For a change of `Δ` in the value of variable `var`, returns if new left-hand side consumption in constraint `c` is feasible given the current consumption `currentLHS`.
# Example
```jldoctest
julia> c = OREnvironment.constructConstraint(12.0, :lessOrEq, [1, 3, 4], [2.0, 5.0, 3.0]);
julia> Δ = 3; currentLHS = 9.0;
julia> OREnvironment.is_increment_feasible(1, Δ, currentLHS, c)
false # new lhs would be 15 which larger than 12
```
"""
function is_increment_feasible(variable::Int,
                               Δvariable::Real,
                               currentLHS::Float64,
                               c::Constraint)::Bool
    local lhs::Float64 = compute_lhs_after_increment(variable, Δvariable, currentLHS, c);
    local answer::Bool = is_feasible(c, lhs);
    return answer;
end

############################
# General methods for Constraints when dealing with Solutions
############################
"""
    is_increment_feasible(s, vconstraints, var, Δ, idxConstraints)

Given solution `s`, the constraints of the problem `vconstraints`, a change `Δ` (respect to the current solution `s`) in variable `var` and an Array with the constraint indexes where variable `var` appears, returns if the change is feasible or not.

**Note:** That the change is feasible DOES NOT mean that the solution is feasible. Given a feasible solution, a feasible change implies that the new solution is still feasible. But given an infeasible solution, a feasible change means that the constraints where `var` appears are not violated by the current change. However, as the original solution was infeasible, the new one is still infeasible. If the return of the function is `false`, then the solution after the change is infeasible regardless of the original feasibility situation.

See examples in the file with the tests: *./test/Constraints.jl*.
"""
function is_increment_feasible(s::Solution,
                               constraints::Array{<:Constraint,1},
                               variable::Int,
                               Δvariable::Real,
                               idxConstraints::Array{Int,1})::Bool
    @inbounds for i in idxConstraints
        local currentLHS::Float64 = get_constraint_consumption(s, i);
        local feasible::Bool = is_increment_feasible(variable, Δvariable, currentLHS, constraints[i]);
        if !feasible
            return false;
        end
    end
    return true;
end

"""
    is_current_consumption_feasible(s, vconstraints)

Returns true if the current consumption associated with solution `s` is feasible respect to constraints `vconstraints`.

**Note:** This function DOES NOT compute the current consumption of solution `s`. It just checks the last values saved in memory and compare them with the right-hand side values of the constraints. This is for efficiency reasons. Usually, given a solution `s`, when you locally modified it, it is possible to efficiently recompute the current consumption of the solution by just checking the constraints affected by the change. Therefore, to check feasibility it would be a waste of time to recompute the left-hand sides again instead of checking the values in memory.

If the lhs must be recomputed to check feasibility, use function [`is_feasible(s::Solution, constraints::Array{<:Constraint,1})`](@ref) instead.

See examples in the file with the tests: *./test/Constraints.jl*.
"""
function is_current_consumption_feasible(s::Solution,
                                         constraints::Array{<:Constraint,1})::Bool
    local N::Int = length(constraints);
    @inbounds for i in 1:N
        local lhs::Float64 = get_constraint_consumption(s, i);
        local feasible::Bool = is_feasible(constraints[i], lhs);
        if !feasible
            return false;
        end
    end
    return true;
end

"""
    compute_lhs(c, s)

Computes the consumption (left-hand side) of constraint `c` by solution `s`.

**Note:** It returns the value of the lhs, but it DOES NOT modified the consumption status of solution `s`.

See examples in the file with the tests: *./test/Constraints.jl*.
"""
function compute_lhs(c::Constraint, s::Solution)::Float64
    local lhs::Float64 = 0.0;
    for (variable, coef) in c.variablesPositiveCoefficients
         lhs += coef*get_solution(s, variable);
    end
    for (variable, coef) in c.variablesNegativeCoefficients
         lhs += coef*get_solution(s, variable);
    end
    return lhs;
end

"""
    is_feasible(s, vconstraints)

Returns true if the consumption incurred by solution `s` is feasible respect to constraints `vconstraints`.

**Note:** This function computes the left-hand sides of constraints using `s`, but it DOES NOT update the current consumption of solution `s`.

**Tip:** If the lhs does not have to be recomputed, then check function [`is_current_consumption_feasible`](@ref) for a better performance. Actually, that function should be preferred over `is_feasible(s, vconstraints)`.

See examples in the file with the tests: *./test/Constraints.jl*.
"""
function is_feasible(s::Solution,
                     constraints::Array{<:Constraint,1})::Bool
    for c in constraints
        local lhs::Float64 = compute_lhs(c, s);
        local feasible::Bool = is_feasible(c, lhs);
        if !feasible
            return false;
        end
    end
    return true;
end

"""
    is_active(c, s)

Returns true if constraint `c` is active under the consumption incurred by solution `s`.

**Note:** This function computes the left-hand side of the constraint using `s`, but it DOES NOT update the current consumption of solution `s`.

**Tip:** If the lhs does not have to be recomputed, then check function [`is_active_under_current_consumption`](@ref) for a better performance. Actually, that function should be preferred over this one.

See examples in the file with the tests: *./test/Constraints.jl*.
"""
is_active(c::Constraint, s::Solution)::Bool = is_active(c, compute_lhs(c, s));

"""
    is_active_under_current_consumption(c, idx, s)

Returns true if the `idx`-th constraint (`c`) is active under the consumption incurred by solution `s`.

**Note:** This function DOES NOT compute the current consumption of solution `s`. It just checks the last value saved in memory and compare it with the right-hand side values of the constraint. That is why `idx` must be provided.

If the lhs must be recomputed, use function [`is_active(c::Constraint, s::Solution)`](@ref) instead.

See examples in the file with the tests: *./test/Constraints.jl*.
"""
function is_active_under_current_consumption(c::Constraint, idxConstraint::Int, s::Solution)::Bool
    local lhs::Float64 = get_constraint_consumption(s, idxConstraint);
    return is_active(c, lhs);
end

ext: jl


text:
const NAUX = 10
const MAXTEXTURE = 1000

@enum mjtGridPos begin
   GRID_TOPLEFT
   GRID_TOPRIGHT
   GRID_BOTTOMLEFT
   GRID_BOTTOMRIGHT
end

@enum mjtFramebuffer begin
   FB_WINDOW
   FB_OFFSCREEN
end

@enum mjtFontScale begin
   FONTSCALE_50  = 50
   FONTSCALE_100 = 100
   FONTSCALE_150 = 150
   FONTSCALE_200 = 200
   FONTSCALE_250 = 250
   FONTSCALE_300 = 300
end

@enum mjtFont begin
   FONT_NORMAL
   FONT_SHADOW
   FONT_BIG
end

struct mjrRect
   left::Cint
   bottom::Cint
   width::Cint
   height::Cint
end

struct mjrContext
   lineWidth::Cfloat
   shadowClip::Cfloat
   shadowScale::Cfloat
   fogStart::Cfloat
   fogEnd::Cfloat
   fogRGBA::SVector{4, Cfloat}
   shadowSize::Cint
   offWidth::Cint
   offHeight::Cint
   offSamples::Cint
   fontScale::Cint
   auxWidth::SVector{NAUX, Cint}
   auxHeight::SVector{NAUX, Cint}
   auxSamples::SVector{NAUX, Cint}

   offFBO::UInt32
   offFBO_r::UInt32
   offColor::UInt32
   offColor_r::UInt32
   offDepthStencil::UInt32
   offDepthStencil_r::UInt32
   shadowFBO::UInt32
   shadowTex::UInt32

   auxFBO::SVector{NAUX, UInt32}
   auxFBO_r::SVector{NAUX, UInt32}
   auxColor::SVector{NAUX, UInt32}
   auxColor_r::SVector{NAUX, UInt32}

   ntexture::Cint
   textureType::SVector{100, Cint}
   texture::SVector{100, UInt32}
   basePlane::UInt32
   baseMesh::UInt32
   baseHField::UInt32
   baseBuiltin::UInt32
   baseFontNormal::UInt32
   baseFontShadow::UInt32
   baseFontBig::UInt32
   rangePlane::Cint
   rangeMesh::Cint
   rangeHField::Cint
   rangeBuiltin::Cint
   rangeFont::Cint

   nskin::Cint
   skinvertVBO::Ptr{UInt32}
   skinnormalVBO::Ptr{UInt32}
   skintexcoordVBO::Ptr{UInt32}
   skinfaceVBO::Ptr{UInt32}

   charWidth::SVector{127, Cint}
   charWidthBig::SVector{127, Cint}
   charHeight::Cint
   charHeightBig::Cint

   glewInitialized::Cint
   windowAvailable::Cint
   windowSamples::Cint
   windowStereo::Cint
   windowDoublebuffer::Cint
   currentBuffer::Cint
end


text:
#=
  Langford's number problem L(2,N) in Julia ConstraintSolver.

  Langford's number problem (CSP lib problem 24)
  http://www.csplib.org/prob/prob024/
  """
  Arrange 2 sets of positive integers 1..k to a sequence,
  such that, following the first occurence of an integer i,
  each subsequent occurrence of i, appears i+1 indices later
  than the last.
  For example, for k=4, a solution would be 41312432
  """

  * John E. Miller: Langford's Problem
    http://www.lclark.edu/~miller/langford.html

  * Encyclopedia of Integer Sequences for the number of solutions for each k
    http://www.research.att.com/cgi-bin/access.cgi/as/njas/sequences/eisA.cgi?Anum=014552

  Note: For k=4 there are two different solutions:
     solution:[4,1,3,1,2,4,3,2]
     position:[2,5,3,1,4,8,7,6]
  and
     solution:[2,3,4,2,1,3,1,4]
     position:[5,1,2,3,7,4,6,8]

  Which this symmetry breaking

     Solution[1] #< Solution[K2],

  then just the second solution is shown.

  Note: There are only solutions when K mod 4 == 0 or K mod 4 == 3.

  - For a generalized version L(k,n) where k >= 2, see http://hakank.org/langford_generalized.pi
  - For Nickerson's variant, see http://hakank.org/langford_nickerson.pi

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# Note: my_element is defined in constraints_utils.jl
#
function langford(k=5,symmetry_breaking=true,all_solutions=true,print_solutions=true)

    if !((k % 4 == 0) || (k % 4 == 3))
        println("There is no solutions for k ($k) unless K mod 4 == 0 or K mod 4 == 3")
        return
    end

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                           # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                           "branch_strategy" => :IMPS, # default
                                                           # "branch_strategy" => :ABS, # Activity Based Search
                                                           # "activity.decay" => 0.999, # default 0.999
                                                           # "activity.max_probes" => 10, # default, 10
                                                           # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>16,

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                        ))
    k2 = 2*k

    @variable(model, 1 <= position[1:k2] <= k2, Int)
    @variable(model, 1 <= solution[1:k2] <= k, Int)

    @constraint(model, position in CS.AllDifferentSet())

    # symmetry breaking
    if symmetry_breaking
        @constraint(model, solution[1] <= solution[k2])
    end

    for i in 1:k
        @constraint(model, position[k+i] == position[i] + i+1)
        # @constraint(model, solution[position[i]] == i )
        my_element(model, position[i], solution, i)
        my_element(model, position[k+i], solution, i)
    end

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                positionx = convert.(Integer,JuMP.value.(position; result=sol))
                solutionx = convert.(Integer,JuMP.value.(solution; result=sol))
                # println("position:$positionx\n")
                println("solution:$solutionx position:$positionx\n")
            end
        end
    else
        println("status:$status")
    end

    return status
end

all_solutions = true
print_solutions = true
symmetry_breaking = true
@time langford(7,symmetry_breaking,all_solutions,print_solutions)

for n in 1:20
    if n % 4 == 3  || n % 4 == 0
        println("n:$n")
        @time status = langford(n,true,all_solutions,false)
        if status != MOI.OPTIMAL
            break
        end
    end
end

println("\nTime to first solution:")
all_solutions = false
print_solutions = true
symmetry_breaking = true
for n in 1:20
    if n % 4 == 3  || n % 4 == 0
        println("n:$n")
        @time status = langford(n,symmetry_breaking,all_solutions,print_solutions)
        if status != MOI.OPTIMAL
            break
        end
    end
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
using DataFrames, Dates, Plots, StatsPlots, Distributions, LinearAlgebra, Printf

function importPQ(csvprefix::String,nfiles::Int)
    P = fill(DataFrame(),nfiles)
    Q = fill(DataFrame(),nfiles)
    for i in 1:nfiles
        P[i] = CSV.read(csvprefix*"P$(i).csv", DataFrame, comment="#", types=[Time, Float64], normalizenames=true)
        Q[i] = CSV.read(csvprefix*"Q$(i).csv", DataFrame, comment="#", types=[Time, Float64], normalizenames=true)
        println("start (P,Q): ($(P[i].timestamp[1]),$(Q[i].timestamp[1]))\tlength: ($(length(P[i].vb_va)),$(length(Q[i].vb_va)))")
        Millisecond(-5) < (P[i].timestamp[1]-Q[i].timestamp[1]) < Millisecond(5) || @warn "t0 differs by >5ms between P$(i) and Q$(i)"
        length(P[i].vb_va) == length(Q[i].vb_va) || @warn "P Q lengths not equal. Samples may be missing."
    end
    return P,Q
end

function checksampling(labels::Array{String,1}, datasets::Array{DataFrame,1}...; boxplt=false)
    length(labels) == length(datasets) || throw(ArgumentError("Length of labels does not match number of datasets."))
    rateplts = Dict();
    for (data,label) in zip(datasets,labels)
        δts = [data[i].timestamp[2:end].-data[i].timestamp[1:end-1] for i in 1:length(data)]
        hs = [unique(δt) for δt in δts]
        println(label*" set of sampling rates")
        display(hs)
        inrange = Millisecond(4) .<= vcat(hs...) .<= Millisecond(6)
        .&(inrange...) || @warn "sampling rate outside range 4-6 milliseconds"
        if boxplt
            rateplts[label] = boxplot([[δt|>Dates.value|>x->x/1000 for δt in δts_] for δts_ in δts], ylabel = "δt [μs]")
            xticks!(:none)
        else
            rateplts[label] = histogram([[δt|>Dates.value|>x->x/1000 for δt in δts_] for δts_ in δts], xlabel = "δt [μs]", ylabel = "n", bins=20, alpha=0.5)
        end
        title!(label)
    end
    plot(values(rateplts)..., size = (length(rateplts)*300,400))
end

function plotPQ(P::Array{DataFrame,1}, Q::Array{DataFrame,1}, title::String; scatterplt=true, size=(1200,400))
    pltP = plot();
    pltQ = plot();
    length(P) == length(Q) || @warn "unequal number of dataframes in P and Q."
    for i in 1:length(P)
        nrow(P[i]) == nrow(Q[i]) || @warn "unequal number of rows in P[$(i)] and Q[$(i)]."
        scatter!(pltP, P[i].timestamp, P[i].vb_va, label = "P$(i)", size = size)
        xlabel!("Time [s]")
        ylabel!("P [mmHg]")
        scatter!(pltQ, Q[i].timestamp, Q[i].vb_va, label = "Q$(i)", size = size)
        xlabel!("Time [s]")
        ylabel!("Q [L/min]")
    end
    plot(pltP,pltQ, layout=(2,1), title=[title ""], xlabel = ["" "Time [s]"])
end

"""
Takes time aligned vectors pressure `p` and flow `q` in L/min with samling rate h and returns volume vector `v` and a PV plot.
"""
function pvloop(p,q,h)
    #FIXME: peak detect and plot beat by beat.
    v = 1000/60*h.*cumsum(q) # L/min ⇒ milliliters
    pltPV = plot(v, p, ylabel="Pressure [mmhg]", xlabel="Volume [mL]", legend=:none);
    return v, pltPV
end

"""
Takes dataframes p and q (with columns "timestamp" and "vb_va") of pressure in mmHg and flow in L/min and returns a plot including p vs t, q vs t, and p vs volume.
"""
function plotPQV(p::DataFrame, q::DataFrame; h=0.005, pheight= 0.7, size=(800,400))
    nrow(p) == nrow(q) || @warn "unequal number of rows in P and Q."
    pltP = plot(p.timestamp, p.vb_va)
    ylabel!("Pressure [mmHg]")
    pltQ = plot(q.timestamp, q.vb_va)
    ylabel!("Flow [L/min]")
    pltPQ = plot(pltP,pltQ, layout = grid(2, 1, heights=[pheight,1-pheight]), xlabel = ["" "Time [s]"], legend=:none)
    v = 1000/60*h.*cumsum(q.vb_va) # L/min ⇒ milliliters
    pltPV = plot(v, p.vb_va, xlabel="Volume [mL]", legend=:none, ylabel="Pressure [mmHg]");
    plt = plot(pltPQ,pltPV,layout=(1,2),size=size,grid=true)
    return plt, v
end

function plot_result(θhat, x0hat, y, u, h)
    yhat, xhat = simulate(θhat,x0hat,u,h,xout=true)
    t = range(0,step=h,length=length(u))
    plty = plot([t,t],[y,yhat], xticks=:none, label = ["Experiment" "Model"])
    title!(@sprintf "Model Fit: Rp=%0.2f, C=%0.2f, Rc=%0.2f L=%0.2f (P̄/Q̄=%0.2f)" vcat(θhat, mean(y)/mean(u))...)
    ylabel!("AoP [mmHg]")
    pltu = plot(t,u,label=:none)
    ylabel!("AoQ [L/min]")
    xlabel!("Time [s]")
    yhat, xhat, t, plot(plty, pltu, grid=true, layout = grid(2, 1, heights=[0.67,0.33]))
end

function id_multiple_init(u, y, x00, θ0, n=100; h=0.005, θid=[1,1,1,1], x0id=true,x0solve=false, G_=G, paramdist=[true,true,true,true], σy=0, σu=0,)
    # fixme add check on paradist length
    id = IdConditions(θ0,x00,θid,x0id,x0solve,G_)

    θdists = fill(Uniform(0,1),length(θ0))
    for (i,dist) in enumerate(paramdist)
        if (dist)
            θdists[i] = Uniform(0.001*θ0[i],10*θ0[i])
        end
    end

    θ0s = []
    θhats = []
    x0hats = []
    res = []

    for _ in 1:n
        y_in = σy == 0 ? y : y + randn(length(u),1).*σy
        u_in = σu == 0 ? u : u + randn(length(u),1).*σu
        θ0_ = [dist ? rand(θdist) : θ0_ for (dist, θdist, θ0_) in zip(paramdist, θdists, θ0)]
        x00_ = solvex0(G_(θ0_,h), u_in)
        @time phat, res_ = identify(θ0_,x00_,u_in,y_in,h,θid=θid,x0id=x0id,x0solve=x0solve, G_=G_);
        θhat_, x0hat_ = deal(phat,id)
        if id.x0solve x0hat_ = solvex0(G_(θhat_,h),u_in) end

        push!(θ0s, θ0_)
        push!(θhats, θhat_)
        push!(x0hats, x0hat_)
        push!(res, res_)
    end
    θ0s, θhats, x0hats, res
end

"""
Summarize results from multiple initialization.
"""
function multiinitresults(θhats, x0hats, res)
    # fixme uses global u
    θlabels = ["Rp", "C", "Rc", "L"]
    # generalize for any θ:
    # θlabels = ["θ$(n)" for n in 1:length(θhats[1])]
    x0labels = ["x0$(n)" for n in 1:length(x0hats[1])]

    results = DataFrame([Symbol(label)=>Float64[] for label in vcat(θlabels[1:length(θhats[1])],x0labels)]...)

    for (θhat, x0hat) in zip(θhats, x0hats)
        push!(results, vcat(θhat, x0hat))
    end

    last_iters = [Optim.trace(res_)[end] for res_ in res]; # use other than end?
    results.MSE = [last_iter.value for last_iter in last_iters]./length(u);
    results.g = [last_iter.metadata["g(x)"] for last_iter in last_iters];
    results.hess = [last_iter.metadata["h(x)"] for last_iter in last_iters];
    results.poshess = [isposdef(last_iter.metadata["h(x)"]) for last_iter in last_iters];
    results.hesscond = [cond(last_iter.metadata["h(x)"]) for last_iter in last_iters];
    results.hessSVD = [svd(last_iter.metadata["h(x)"]) for last_iter in last_iters];
    results.idx = [1:length(θhats)...];

    @show sort!(results, [:MSE,:hesscond], rev=(false, false));

    nanratio = sum(isnan.(results[!,:MSE]))/length(results[!, :MSE])

    delete!(results, isnan.(results[!,:MSE]))

    return results, nanratio
end

"""
plot parameter histograms before (θ0) and after (θhat) multiple initialization fits.
"""
function fithistograms(θ0s,θhats,θlabels)
    hists0 = []
    hists = []
    for i in 1:length(θlabels)
        push!(hists0, histogram([θ0[i] for θ0 in θ0s], bins=100, ylabel=θlabels[i]))
        hist = histogram([θhat[i] for θhat in θhats], bins=100)
        #vline!(hist, [θ0[i]], c=:red, lw=2, ls=:dash, label=θlabels[i]*" true")
        push!(hists, hist)
    end

    title!(hists0[1],"initial θ");
    hist0_plt = plot(hists0..., layout=(length(hists0),1), legend=:none);
    title!(hists[1],"fit θ");
    hist_plt = plot(hists..., layout=(length(hists0),1), legend=:none);

    plot(hist0_plt, hist_plt, layout=(1,2))
end

max_stars_repo_name: hpigot/windkessel-id

max_issues_repo_name: hpigot/windkessel-id

max_forks_repo_name: hpigot/windkessel-id


text:
using Random, Statistics, DataFrames, Distributed

# Run multiple instances of a learner
# (type and parameter values defined in 'settings')
# on Gaussian environments (Section 2.3.1)
# (defined in 'specs')
function comparison(settings, specs::SpecsGaussian)

    df = @distributed vcat for setting in settings
        # @show setting
        # @show (setting[1][1].type, setting[1][2], setting[1][3])
        parameterset = getparamset(setting)
        # @show parameterset
        callback = run!(specs.envtype, setting[1][1].type, parameterset,
                        sigma = setting[1][2],
                        changeprobability = setting[1][3],
                        nsteps = specs.nsteps,
                        iseedenvpolicy = specs.seedenvpolicy[setting[2]],
                        iseedlearner = specs.seedlearner[setting[2]]
                        )
        merror = computesquareddiff(callback.true_history[2:end],
                                    callback.estimated_history[2:end])
        DataFrame(learner = setting[1][1].name,
                type = setting[1][1].type,
                param = parameterset,
                sigma = setting[1][2],
                changeprobability = setting[1][3],
                seed = specs.seedenvpolicy[setting[2]],
                seedlearner = specs.seedlearner[setting[2]],
                merrors = [merror],
                states = [callback.state_history],
                switches = [callback.switchflag],
                #nposteriors = [callback.nposteriors[2:end]]
                # nresamplings = [callback.nresamplings_history[end]]
                )
        end
    df
end


text:
abstract type AbstractNode end

abstract type AbstractScene <: AbstractNode end

abstract type AbstractIOEvent end

abstract type AbstractMotion end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
##############################################################################
#
# Text output
#
##############################################################################

function escapedprint(io::IO, x::Any, escapes::String)
    print(io, x)
end

function escapedprint(io::IO, x::String, escapes::String)
    print_escaped(io, x, escapes)
end

function printtable(io::IO,
                    df::AbstractDataFrame;
                    header::Bool = true,
                    separator::Char = ',',
                    quotemark::Char = '"',
                    nastring::String = "NA")
    n, p = size(df)
    etypes = eltypes(df)
    if header
        cnames = _names(df)
        for j in 1:p
            print(io, quotemark)
            print(io, cnames[j])
            print(io, quotemark)
            if j < p
                print(io, separator)
            else
                print(io, '\n')
            end
        end
    end
    for i in 1:n
        for j in 1:p
            if ! (isna(df[j],i))
                if ! (etypes[j] <: Real)
		    print(io, quotemark)
		    escapedprint(io, df[i, j], "\"'")
		    print(io, quotemark)
                else
		    print(io, df[i, j])
                end
            else
		print(io, nastring)
            end
            if j < p
		print(io, separator)
            else
                print(io, '\n')
            end
        end
    end
    return
end

function printtable(df::AbstractDataFrame;
                    header::Bool = true,
                    separator::Char = ',',
                    quotemark::Char = '"',
                    nastring::String = "NA")
    printtable(STDOUT,
               df,
               header = header,
               separator = separator,
               quotemark = quotemark,
               nastring = nastring)
    return
end

# Infer configuration settings from filename
function writetable(filename::String,
                    df::AbstractDataFrame;
                    header::Bool = true,
                    separator::Char = getseparator(filename),
                    quotemark::Char = '"',
                    nastring::String = "NA",
                    append::Bool = false)

    if endswith(filename, ".bz") || endswith(filename, ".bz2")
        throw(ArgumentError("BZip2 compression not yet implemented"))
    end

    if append && isfile(filename) && filesize(filename) > 0
        file_df = readtable(filename, header = false, nrows = 1)

        # Check if number of columns matches
        if size(file_df, 2) != size(df, 2)
            throw(DimensionMismatch("Number of columns differ between file and DataFrame"))
        end

        # When 'append'-ing to a nonempty file,
        # 'header' triggers a check for matching colnames
        if header
            if any(i -> symbol(file_df[1, i]) != index(df)[i], 1:size(df, 2))
                throw(KeyError("Column names don't match names in file"))
            end

            header = false
        end
    end

    openfunc = endswith(filename, ".gz") ? gzopen : open

    openfunc(filename, append ? "a" : "w") do io
        printtable(io,
                   df,
                   header = header,
                   separator = separator,
                   quotemark = quotemark,
                   nastring = nastring)
    end

    return
end

##############################################################################
#
# HTML output
#
##############################################################################

function html_escape(cell::String)
    cell = replace(cell, "&", "&")
    cell = replace(cell, "<", "<")
    cell = replace(cell, ">", ">")
    return cell
end

function Base.writemime(io::IO,
                        ::MIME"text/html",
                        df::AbstractDataFrame)
    n = size(df, 1)
    cnames = _names(df)
    write(io, "
")
    write(io, "
")
    write(io, "	")
    for column_name in cnames
        write(io, "$column_name	")
    end
    write(io, "
")
    tty_rows, tty_cols = Base.tty_size()
    mxrow = min(n,tty_rows)
    for row in 1:mxrow
        write(io, "
")
        write(io, "$row	")
        for column_name in cnames
            cell = string(df[row, column_name])
            write(io, "$(html_escape(cell))	")
        end
        write(io, "
")
    end
    if n > mxrow
        write(io, "
")
        write(io, "⋮	")
        for column_name in cnames
            write(io, "⋮	")
        end
        write(io, "
")
    end
    write(io, "
")
end

##############################################################################
#
# MIME
#
##############################################################################

function Base.writemime(io::IO,
                        ::MIME"text/csv",
                        df::AbstractDataFrame)
    printtable(io, df, true, ',')
end

function Base.writemime(io::IO,
                        ::MIME"text/tab-separated-values",
                        df::AbstractDataFrame)
    printtable(io, df, true, '\t')
end


text:
@defcomp welfare begin
    CEMUTOTPER      = Variable(index=[time])    #Period utility
    CUMCEMUTOTPER   = Variable(index=[time])    #Cumulative period utility
    PERIODU         = Variable(index=[time])    #One period utility function
    UTILITY         = Variable()                #Welfare Function

    CPC             = Parameter(index=[time])   #Per capita consumption (thousands 2010 USD per year)
    l               = Parameter(index=[time])   #Level of population and labor
    rr              = Parameter(index=[time])   #Average utility social discount rate
    elasmu          = Parameter()               #Elasticity of marginal utility of consumption
    scale1          = Parameter()               #Multiplicative scaling coefficient
    scale2          = Parameter()               #Additive scaling coefficient

    function run_timestep(p, v, d, t)
        # Define function for PERIODU
        v.PERIODU[t] = (p.CPC[t] ^ (1 - p.elasmu) - 1) / (1 - p.elasmu) - 1

        # Define function for CEMUTOTPER
        v.CEMUTOTPER[t] = v.PERIODU[t] * p.l[t] * p.rr[t]

        # Define function for CUMCEMUTOTPER
        v.CUMCEMUTOTPER[t] = v.CEMUTOTPER[t] + (!is_first(t) ? v.CUMCEMUTOTPER[t-1] : 0)

        # Define function for UTILITY
        if t == TimestepIndex(100)
            v.UTILITY = 5 * p.scale1 * v.CUMCEMUTOTPER[t] + p.scale2

            utility = 5 * p.scale1 * v.CUMCEMUTOTPER[t] + p.scale2
        end
    end
end


text:
using OperatorFlux
using Test
using Flux
using LinearAlgebra

@testset "SpectralCovariance - FourierTransform" begin
    # auxiliary data
    M = 32

    # test constructor
    trafo = FourierTransform(modes = (3, 3))
    ch = 3 => 1
    rank = 2
    sc = SpectralCovariance(trafo, rank, ch)
    @test sc.transform == trafo
    @test sc.in_channels == 3
    @test sc.out_channels == 1
    @test ndims(sc.weights_μ) == 2
    @test ndims(sc.weights_d) == 2
    @test ndims(sc.weights_w) == 3
    @test size(sc.weights_μ) == (3, 1)
    @test size(sc.weights_d) == (3, 1)
    @test size(sc.weights_w) == (3, 2, 1)

    # test covariance
    trafo = FourierTransform(modes = (3, 3))
    ch = 3 => 1
    rank = 2
    sc = SpectralCovariance(trafo, rank, ch)
    M = 32
    z = rand(3, M, M, 12)
    μ, D, V = sc(z)
    @test size(μ) == (1, M, M, 12)
    @test size(D) == (1, M, M, 12)
    @test size(V) == (1, M, M, rank, 12)
    @test all(D .> 0)
end

@testset "SpectralCovariance - ChebyshevTransform" begin
    # auxiliary data
    M = 32

    # test constructor
    trafo = ChebyshevTransform(modes = (3, 3))
    ch = 3 => 1
    rank = 2
    sc = SpectralCovariance(trafo, rank, ch)
    @test sc.transform == trafo
    @test sc.in_channels == 3
    @test sc.out_channels == 1
    @test ndims(sc.weights_μ) == 2
    @test ndims(sc.weights_d) == 2
    @test ndims(sc.weights_w) == 3
    @test size(sc.weights_μ) == (3, 1)
    @test size(sc.weights_d) == (3, 1)
    @test size(sc.weights_w) == (3, 2, 1)

    # test covariance
    trafo = ChebyshevTransform(modes = (3, 3))
    ch = 3 => 1
    rank = 2
    sc = SpectralCovariance(trafo, rank, ch)
    M = 32
    z = rand(3, M, M, 12)
    μ, D, V = sc(z)
    @test size(μ) == (1, M, M, 12)
    @test size(D) == (1, M, M, 12)
    @test size(V) == (1, M, M, rank, 12)
    @test all(D .> 0)
end

@testset "SpectralCovariance - LegendreTransform" begin
    # auxiliary data
    M = 32

    # test constructor
    trafo = LegendreTransform(modes = (3, 3))
    ch = 3 => 1
    rank = 2
    sc = SpectralCovariance(trafo, rank, ch)
    @test sc.transform == trafo
    @test sc.in_channels == 3
    @test sc.out_channels == 1
    @test ndims(sc.weights_μ) == 2
    @test ndims(sc.weights_d) == 2
    @test ndims(sc.weights_w) == 3
    @test size(sc.weights_μ) == (3, 1)
    @test size(sc.weights_d) == (3, 1)
    @test size(sc.weights_w) == (3, 2, 1)

    # test covariance
    trafo = LegendreTransform(modes = (3, 3))
    ch = 3 => 1
    rank = 2
    sc = SpectralCovariance(trafo, rank, ch)
    M = 32
    z = rand(3, M, M, 12)
    μ, D, V = sc(z)
    @test size(μ) == (1, M, M, 12)
    @test size(D) == (1, M, M, 12)
    @test size(V) == (1, M, M, rank, 12)
    @test all(D .> 0)
end


text:
# Primal infeasible test problems
# x >= 0, elements in A >=0, elements in b <0 and s in R+

# module PrimalInfeasibleTest
# using COSMO, Test, Random, SparseArrays, LinearAlgebra
# include("./../COSMOTestUtils.jl")

# This test is precision agnostic and will be run with T precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

@testset "Primal infeasible QP problems - Testset 1 (T = $(T))" for T in UnitTestFloats
    # for iii = 1:nn

    if precision(T) >= precision(Float64)
        nn = 1
        rng = Random.MersenneTwister(74747)

        # choose size of problem
        n = rand(rng,  5:50)
        m = 2 * n
        A = sprand(rng, T, m, n, 0.8)
        b = -rand(rng, T, m)
        A = [A; -spdiagm(0 => ones(T, n))]
        b = [b; zeros(T, n)]

        # create dual feasible problem
        P = generate_pos_def_matrix(rng, n; MT = T)
        ytrue = rand(rng, T, m + n, 1)
        xtrue = rand(rng, T, n, 1)
        q = (-P * xtrue -  A' * ytrue)[:]

        constraint = COSMO.Constraint(-A, b, COSMO.Nonnegatives)

        settings = COSMO.Settings{T}(max_iter=10000, eps_abs = 1e-5, eps_rel=1e-5)
        model = COSMO.Model{T}()
        assemble!(model, P, q, constraint, settings = settings)
        res1 = COSMO.optimize!(model);

       @test res1.status == :Primal_infeasible
   end
 end


text:
export get_leaf_nodes, get_leaf_nodes!

"""
    get_leaf_nodes(jlt::AbstractJLBoostTreeModel)

jlt - The JLBoostTree

"""
function get_leaf_nodes(jlt::AbstractJLBoostTree)
    get_leaf_nodes!(AbstractJLBoostTree[], jlt)
end

function get_leaf_nodes!(leaf_nodes::Vector{AbstractJLBoostTree}, jlt::AbstractJLBoostTree)
    if has_children(jlt)
        for child in children(jlt)
            get_leaf_nodes!(leaf_nodes, child)
        end
    else
        push!(leaf_nodes, jlt)
    end
    leaf_nodes
end


text:
struct DoubleFloat{T} <: MultipartFloat
   hi::T
   lo::T

   # these forms ensure hi, lo are canonically valued
   function DoubleFloat{T}(hi::T, lo::T) where {T<:IEEEFloat}
       hi, lo = add_2(hi, lo)
       return new{T}(hi, lo)
   end

   function DoubleFloat{DoubleFloat{T}}(hi::DoubleFloat{T}, lo::DoubleFloat{T}) where {T<:IEEEFloat}
     hi1, lo1 = HILO(hi)
     hi2, lo2 = HILO(lo)
     hihi, hilo, lohi, lolo = add_4(hi1, hi2, lo1, lo2)
     zhi = DoubleFloat{T}(hihi, hilo)
     zlo = DoubleFloat{T}(lohi, lolo)
     return new{DoubleFloat{T}}(zhi, zlo)
   end
   # this form does not alter hi, lo
   function DoubleFloat(hi::T, lo::T) where {T<:AbstractFloat}
       return new{T}(hi, lo)
   end
end

const Double64 = DoubleFloat{Float64}
const Double32 = DoubleFloat{Float32}
const Double16 = DoubleFloat{Float16}

const QuadrupleF64 = DoubleFloat{DoubleFloat{Float64}}
const QuadrupleF32 = DoubleFloat{DoubleFloat{Float32}}
const QuadrupleF16 = DoubleFloat{DoubleFloat{Float16}}

@inline HI(x::T) where {T<:IEEEFloat} = x
@inline LO(x::T) where {T<:IEEEFloat} = zero(F)
@inline HILO(x::T) where {T<:IEEEFloat} = x, zero(F)

@inline HI(x::T) where {F<:IEEEFloat, T<:DoubleFloat{F}} = x.hi
@inline LO(x::T) where {F<:IEEEFloat, T<:DoubleFloat{F}} = x.lo
@inline HILO(x::T) where {F<:IEEEFloat, T<:DoubleFloat{F}} = x.hi, x.lo

@inline HI(x::T) where {F<:IEEEFloat, G<:DoubleFloat{F}, T<:DoubleFloat{G}} = x.hi
@inline LO(x::T) where {F<:IEEEFloat, G<:DoubleFloat{F}, T<:DoubleFloat{G}} = x.lo
@inline HILO(x::T) where {F<:IEEEFloat, G<:DoubleFloat{F}, T<:DoubleFloat{G}} = x.hi, x.lo

@inline HI(x::Tuple{T,T}) where {T<:IEEEFloat} = x[1]
@inline LO(x::Tuple{T,T}) where {T<:IEEEFloat} = x[2]
@inline HILO(x::Tuple{T,T}) where {T<:IEEEFloat} = x

@inline HI(x::Tuple{T,T}) where {F<:IEEEFloat, T<:DoubleFloat{F}} = x[1]
@inline LO(x::Tuple{T,T}) where {F<:IEEEFloat, T<:DoubleFloat{F}} = x[2]
@inline HILO(x::Tuple{T,T}) where {F<:IEEEFloat, T<:DoubleFloat{F}} = x

@inline HI(x::Tuple{T,T}) where {F<:IEEEFloat, T<:DoubleFloat{DoubleFloat{F}}} = x[1]
@inline LO(x::Tuple{T,T}) where {F<:IEEEFloat, T<:DoubleFloat{DoubleFloat{F}}} = x[2]
@inline HILO(x::Tuple{T,T}) where {F<:IEEEFloat, T<:DoubleFloat{DoubleFloat{F}}} = x

@inline Double64(x::Tuple{Float64,Float64}) = DoubleFloat(x[1], x[2])
@inline Double32(x::Tuple{Float32,Float32}) = DoubleFloat(x[1], x[2])
@inline Double16(x::Tuple{Float16,Float16}) = DoubleFloat(x[1], x[2])
@inline Double64(x::Tuple{T,T}) where {T<:IEEEFloat} = DoubleFloat(Float64(x[1]), Float64(x[2]))
@inline Double32(x::Tuple{T,T}) where {T<:IEEEFloat} = DoubleFloat(Float32(x[1]), Float32(x[2]))
@inline Double16(x::Tuple{T,T}) where {T<:IEEEFloat} = DoubleFloat(Float16(x[1]), Float16(x[2]))
@inline DoubleFloat(x::Tuple{T,T}) where {T<:AbstractFloat} = DoubleFloat{T}(x[1], x[2])

@inline Double64(hi::T) where {T<:IEEEFloat} = Double64(Float64(hi), zero(Float64))
@inline Double32(hi::T) where {T<:IEEEFloat} = Double32(Float32(hi), zero(Float32))
@inline Double16(hi::T) where {T<:IEEEFloat} = Double16(Float16(hi), zero(Float16))
@inline Double64(hi::T) where {T<:Integer} = Double64(Float64(hi))
@inline Double32(hi::T) where {T<:Integer} = Double32(Float32(hi))
@inline Double16(hi::T) where {T<:Integer} = Double16(Float16(hi))

for (F,D,U) in ((:Float64, :Double64, :(Union{Float32, Float16})),
                (:Float32, :Double32, :(Union{Float64, Float16})),
                (:Float16, :Double16, :(Union{Float64, Float32})))
  @eval begin
    @inline $D(hi::T, lo::T) where {T<:$U} = $D($F(hi), $F(lo))
    @inline $D(hi::T, lo::I) where {T<:$U, I<:Integer} = $D($F(hi), $F(lo))
    @inline $D(hi::I, lo::T) where {T<:$U, I<:Integer} = $D($F(hi), $F(lo))
    @inline $D(hi::I, lo::I) where {I<:Integer} = $D($F(hi), $F(lo))
  end
end

DoubleFloat{T}(x::DoubleFloat{T}, y::T) where {T<:IEEEFloat} = DoubleFloat(x, DoubleFloat(y))
DoubleFloat{T}(x::T, y::DoubleFloat{T}) where {T<:IEEEFloat} = DoubleFloat(DoubleFloat(x), y)
DoubleFloat{T}(x::DoubleFloat{T}, y::DoubleFloat{T}) where {T<:IEEEFloat} = DoubleFloat(x,y)

DoubleFloat{T}(x::DoubleFloat{DoubleFloat{T}}, y::T) where {T<:IEEEFloat} = DoubleFloat(x, DoubleFloat(DoubleFloat(y),DoubleFloat(zero(T))))
DoubleFloat{T}(x::T, y::DoubleFloat{DoubleFloat{T}}) where {T<:IEEEFloat} = DoubleFloat(DoubleFloat(DoubleFloat(x), DoubleFloat(zero(T))), y)
DoubleFloat{T}(x::DoubleFloat{DoubleFloat{T}}, y::DoubleFloat{T}) where {T<:IEEEFloat} = DoubleFloat(x, DoubleFloat(y, DoubleFloat(zero(T))))
DoubleFloat{T}(x::DoubleFloat{T}, y::DoubleFloat{DoubleFloat{T}}) where {T<:IEEEFloat} = DoubleFloat(DoubleFloat(x, DoubleFloat(zero(T))), y)
DoubleFloat{T}(x::DoubleFloat{DoubleFloat{T}}, y::DoubleFloat{DoubleFloat{T}}) where {T<:IEEEFloat} = DoubleFloat(x,y)

Double64(x::Double32) = Double64(add_2(Float64(HI(x)), Float64(LO(x))))
Double64(x::Double16) = Double64(add_2(Float64(HI(x)), Float64(LO(x))))
Double32(x::Double16) = Double32(add_2(Float32(HI(x)), Float32(LO(x))))
Double32(x::Double64) = Double32(BigFloat(x))
Double16(x::Double64) = Double16(BigFloat(x))
Double16(x::Double32) = Double16(BigFloat(x))

# a type specific hash function helps the type to 'just work'
const hash_double_lo = (UInt === UInt64) ? 0x9bad5ebab034fe78 : 0x72da40cb
const hash_0_double_lo = hash(zero(UInt), hash_double_lo)

@inline Base.hash(z::DoubleFloat{T}, h::UInt) where {T<:IEEEFloat} =
    hash(z.hi, h) ⊻ hash(z.lo)

@inline Base.hash(z::DoubleFloat{T}, h::UInt) where {F<:IEEEFloat, T<:DoubleFloat{F}} =
    hash(z.hi) ⊻ hash(z.lo)

@inline Base.hash(z::DoubleFloat{T}, h::UInt) where {F<:IEEEFloat, G<:DoubleFloat{F}, T<:DoubleFloat{G}} =
    hash(z.hi) ⊻ hash(z.lo)

Base.precision(::DoubleFloat{T}) where {T<:IEEEFloat} = 2 * precision(T)


text:
module GeometryTests

using Test

@testset "GridGraphs" begin include("GridGraphsTests.jl") end

@testset "Grids" begin include("GridsTests.jl") end

@testset "FaceLabels" begin include("FaceLabelsTests.jl") end

@testset "DiscreteModels" begin include("DiscreteModelsTests.jl") end

@testset "UnstructuredGeometry" begin include("UnstructuredGeometryTests.jl") end

@testset "GeometryWrappers" begin include("GeometryWrappersTests.jl") end

@testset "CartesianGeometry" begin include("CartesianGeometryTests.jl") end

@testset "GridPortions" begin include("GridPortionsTests.jl") end

@testset "BoundaryGrids" begin include("BoundaryGridsTests.jl") end

@testset "SkeletonGrids" begin include("SkeletonGridsTests.jl") end

@testset "SimplexifyTests" begin include("SimplexifyTests.jl") end

end # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -14.06559781669544 <= q <= 13.37751741699121)

                     add_NL_constraint(m, :(softplus(-0.7556338259975872 + -0.2850076752045205*softplus(0.7432962759699406 + 0.616365769922441*softplus(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*softplus(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*softplus(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.5819748667402531*softplus(0.6367931081099356 + 0.08862555157308405*softplus(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*softplus(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*softplus(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + 0.4871037750371179*softplus(0.12597161846431115 + 0.10770543162761159*softplus(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*softplus(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*softplus(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) + softplus(0.9186311390475046 + -0.8559712370251158*softplus(0.7432962759699406 + 0.616365769922441*softplus(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*softplus(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*softplus(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.2512416329037399*softplus(0.6367931081099356 + 0.08862555157308405*softplus(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*softplus(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*softplus(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + 0.7592758779705528*softplus(0.12597161846431115 + 0.10770543162761159*softplus(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*softplus(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*softplus(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) + softplus(-0.7151630404444922 + -0.15417874070690862*softplus(0.7432962759699406 + 0.616365769922441*softplus(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + -0.9563519707270385*softplus(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.21036883565748976*softplus(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.016895941959414795*softplus(0.6367931081099356 + 0.08862555157308405*softplus(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.5030828808953012*softplus(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + -0.4316150506981855*softplus(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5]))) + -0.17047547651198514*softplus(0.12597161846431115 + 0.10770543162761159*softplus(-0.4636854780986108 + 0.31515700957585446*$(x[1]) + 0.3880835610738931*$(x[2]) + 0.5564947842471106*$(x[3]) + -0.6707203467350809*$(x[4]) + 0.9006732364765586*$(x[5])) + 0.46799585735795457*softplus(0.08793534010659565 + 0.5972587053731475*$(x[1]) + 0.5980877634648243*$(x[2]) + 0.5697895882823776*$(x[3]) + -0.6631393507177727*$(x[4]) + -0.46226121533035425*$(x[5])) + 0.45029698767672954*softplus(-0.49150226791798435 + -0.92821223742826*$(x[1]) + -0.3666783267897977*$(x[2]) + -0.7302712314714839*$(x[3]) + 0.924313656245149*$(x[4]) + 0.3686570879245057*$(x[5])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#
# Function that wraps the coordinates to obtain minimum images
# around a defined center
#
# It modifies the x,y,z input vectors
#

function wrap!(sides,x,y,z;center=[0.,0.,0.],sel=[-1])

  # If no selection was set, wrap all
  if sel[1] == -1
    n = length(x)
    for i in 1:n
      wrapone!(i,sides,x,y,z,center)
    end

  # Wrap only the selection
  else
    n = length(sel)
    for i in 1:n
      wrapone!(sel[i],sides,x,y,z,center)
    end
  end

end

function wrapone!(i,sides,x,y,z,center)

  x[i] = x[i] - center[1]
  y[i] = y[i] - center[2]
  z[i] = z[i] - center[3]

  x[i] = x[i]%sides[1]
  y[i] = y[i]%sides[2]
  z[i] = z[i]%sides[3]

  if x[i] > sides[1]/2 ; x[i] = x[i] - sides[1] ; end
  if y[i] > sides[2]/2 ; y[i] = y[i] - sides[2] ; end
  if z[i] > sides[3]/2 ; z[i] = z[i] - sides[3] ; end

  if x[i] < -sides[1]/2 ; x[i] = x[i] + sides[1] ; end
  if y[i] < -sides[2]/2 ; y[i] = y[i] + sides[2] ; end
  if z[i] < -sides[3]/2 ; z[i] = z[i] + sides[3] ; end

  x[i] = x[i] + center[1]
  y[i] = y[i] + center[2]
  z[i] = z[i] + center[3]

end


text:
module day9

import ..aoc19: run_program!
using DelimitedFiles

readInput() = readdlm(joinpath(@__DIR__, "input.txt"), ',', Int) |> vec

function part1(data, programInput = 1)
    output = Channel{Int}(Inf)
    input = Channel{Int}(Inf)
    put!(input, programInput)
    run_program!(data, input, output)
    close(output)
    [o for o in output]
end

part2(data) = part1(data, 2)

end


text:
using IonSim: _alaguerre
export two_ion_ms, rabi_flop

"""
    two_ion_ms(tspan, Ω::Real, ν::Real, δ::Real, η::Real, n̄::Real)
[ref](https://doi.org/10.1103/PhysRevA.62.022311)

Assumes vibrational mode starts in a thermal state with: ``\\langle a^\\dagger a\\rangle = n̄``
and ions start in doubly ground state. Returns `(ρgg, ρee)`, the population in the doubly
ground and doubly excited state, respectively. ``[Ω], [ν], [δ] = Hz``
"""
function two_ion_ms(tspan, Ω::Real, ν::Real, δ::Real, η::Real, n̄::Real)
    ρgg = Float64[]; ρee = Float64[]
    n̄ *= 1.0; Ω *= 2π; ν *= 2π; δ *= 2π
    for t in tspan
        F = -√2 * η * Ω * sin((ν - δ) * t) / (ν - δ)
        G =  -√2 * η * Ω * (1 - cos((ν - δ) * t)) / (ν - δ)
        A = -((η * Ω)^2 / (ν - δ)) * (t - (1/(2 * (ν - δ))) * sin(2 * (ν - δ) * t))
        term1 = (1/2) * exp(-(F^2 + G^2)/4)
        term2 = (1/8) * exp(-(F^2 + G^2))
        ρee_t = 0.0
        ρgg_t  = 0.0
        for n in 0:200
            pn = _Pn(n̄, n)
            isnan(pn) && continue
            t1_1 = _laguerre((F^2 + G^2) / 2, n) * cos(A + (1/2) * F * G) * term1
            t2_2 = _laguerre(2 * (F^2 + G^2), n) * term2
            ρgg_t += pn * (t1_1 + t2_2)
            ρee_t += pn * (-t1_1 + t2_2)
        end
        push!(ρgg, ρgg_t)
        push!(ρee, ρee_t)
    end
    3/8 .+ ρgg, 3/8 .+ ρee
end

"""
    rabi_flop(tspan, Ω::Real, η::Real, n̄::Real; s::Int=0)

Single ion rabi flop. Returns:
``\\sum_{n=0}^∞ p_n sin^2(\\Omega_n t)``
 with
``\\Omega_n = Ωe^{-η^2/2}η^s\\sqrt{\\frac{n!}{(n+s)!}}L_{n}^{s}(η^2)``

where ``s`` is the order of the (blue) sideband that we are driving and ``L_{n}^{s}`` is the
associated Laguerre polynomial. [ref](https://doi.org/10.1103/RevModPhys.75.281)

"""
function rabi_flop(tspan, Ω::Real, η::Real, n̄::Real; s::Int=0)
    n̄ *= 1.0; Ω *= 2π
    p = Vector{Float64}(undef, 0)
    for t in tspan
        pi = 0.0
        for n in 0:300
            pi += _Pn(n̄, n) * sin(Ω/2 * exp(-η^2/2) * η^s * sqrt(1/prod(n+1:n+s)) * _alaguerre(η^2, n, s) * t)^2
        end
        push!(p, pi)
    end
    p
end

function _laguerre(x, n)
    L = 1.0, -x + 1
    if n < 2
        return L[n+1]
    end
    for i in 2:n
        L = L[2], ((2i - 1 - x) * L[2] - (i - 1) * L[1]) / i
    end
    L[2]
end

_Pn(n̄::Real, n::Int) = (n̄ / (n̄ + 1))^n / (n̄ + 1)


text:
using DataDrivenControl
using LinearAlgebra

@testset "irl" begin
    @testset "linear_irl" begin
        # See [1, Example 2. Continuous-Time Optimal Adaptive Control Using IRL]
        # Refs
        # [1] “Reinforcement Learning and Feedback Control: Using Natural Decision Methods to Design Optimal Adaptive Controllers,” IEEE Control Syst., vol. 32, no. 6, pp. 76–105, Dec. 2012, doi: 10.1109/MCS.2012.2214134.
        n, m = 2, 1
        Q = Matrix(I, n, n)
        R = Matrix(I, m, m)
        B = [0 2]'
        irl = DataDrivenControl.LinearIRL(Q, R, B)
        x = rand(2)
        P = [0.0500 0.0039;
             0.0039 0.2085]  # true solution
        w = [P[1, 1], P[2, 1]+P[1, 2], P[2, 2]]
        @test DataDrivenControl.optimal_input(irl, x, w) ≈ - inv(R) * B' * P * x
        @test DataDrivenControl.value(irl, w, x) ≈ x' * P * x
        # data buffer
        ŵ = rand(Int(n*(n+1)/2))
        û = DataDrivenControl.optimal_input(irl, x, ŵ)
        buffer = DataDrivenControl.DataBuffer()
        cost = DataDrivenControl.QuadraticCost(Q, R)
        ts = 0:irl.T:1
        for t in ts
            push!(buffer, irl, cost; t=t, x=x, u=û, w=ŵ)
        end
        @test length(buffer.data_array) == length(ts)
        ŵ_new = deepcopy(ŵ)
        i = deepcopy(irl.i)
        DataDrivenControl.policy_iteration!(irl, buffer, ŵ_new)
        DataDrivenControl.value_iteration!(irl, buffer, ŵ_new)
        @test ŵ_new != ŵ  # updated?
        @test i + 1 == irl.i
        DataDrivenControl.reset!(irl)
        @test irl.i == irl.i_init
    end
end


text:
module Sum
export Σ

function Σ(list) # List summation
  sum = 0
  for i in list
    sum = sum + list[i]
  end
  return sum
end

end


text:
export delay

import Base: close
import Base: show

"""
    delay(delay::Int)

Creates a delay operator, which delays the emission of items from the source Observable
by a given timeout.

# Arguments:
- `delay::Int`: the delay duration in milliseconds (a number) until which the emission of the source items is delayed.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref)
"""
delay(delay::Int) = DelayOperator(delay)

struct DelayOperator <: InferableOperator
    delay :: Int
end

function on_call!(::Type{L}, ::Type{L}, operator::DelayOperator, source) where L
    return proxy(L, source, DelayProxy(operator.delay))
end

operator_right(operator::DelayOperator, ::Type{L}) where L = L

struct DelayProxy <: ActorSourceProxy
    delay :: Int
end

actor_proxy!(::Type{L}, proxy::DelayProxy, actor::A)   where { L, A } = DelayActor{L, A}(proxy.delay, actor)
source_proxy!(::Type{L}, proxy::DelayProxy, source::S) where { L, S } = DelayObservable{L, S}(source)

struct DelayDataMessage{L}
    data :: L
end

struct DelayErrorMessage
    err
end

struct DelayCompleteMessage end

const DelayMessage{L} = Union{DelayDataMessage{L}, DelayErrorMessage, DelayCompleteMessage}

struct DelayQueueItem{L}
    message    :: DelayMessage{L}
    emmited_at :: Float64
end

struct DelayCompletionException <: Exception end

mutable struct DelayActor{L, A} <: Actor{L}
    is_cancelled :: Bool
    delay        :: Int
    actor        :: A
    channel      :: Channel{DelayQueueItem{L}}

    DelayActor{L, A}(delay::Int, actor::A) where { L, A } = begin
        channel = Channel{DelayQueueItem{L}}(Inf)
        self    = new(false, delay, actor, channel)

        task = @async begin
            try
                while !self.is_cancelled
                    item = take!(channel)::DelayQueueItem{L}
                    sleepfor = (item.emmited_at + convert(Float64, self.delay / MILLISECONDS_IN_SECOND)) - time()
                    if sleepfor > 0.0
                        sleep(sleepfor)
                    else
                        yield()
                    end
                    if !self.is_cancelled
                        __process_delayed_message(self, item.message)
                    end
                end
            catch err
                if !(err isa DelayCompletionException)
                    __process_delayed_message(self, DelayErrorMessage(err))
                end
            end
        end

        bind(channel, task)

        return self
    end
end

__process_delayed_message(actor::DelayActor{L}, message::DelayDataMessage{L}) where L = next!(actor.actor, message.data)
__process_delayed_message(actor::DelayActor,    message::DelayErrorMessage)           = begin error!(actor.actor, message.err); close(actor); end
__process_delayed_message(actor::DelayActor,    message::DelayCompleteMessage)        = begin complete!(actor.actor); close(actor); end

on_next!(actor::DelayActor{L}, data::L) where L = put!(actor.channel, DelayQueueItem{L}(DelayDataMessage{L}(data), time()))
on_error!(actor::DelayActor{L}, err)    where L = put!(actor.channel, DelayQueueItem{L}(DelayErrorMessage(err), time()))
on_complete!(actor::DelayActor{L})      where L = put!(actor.channel, DelayQueueItem{L}(DelayCompleteMessage(), time()))

Base.close(actor::DelayActor) = close(actor.channel, DelayCompletionException())

struct DelayObservable{L, S} <: Subscribable{L}
    source :: S
end

function on_subscribe!(observable::DelayObservable, actor::DelayActor)
    return DelaySubscription(actor, subscribe!(observable.source, actor))
end

struct DelaySubscription <: Teardown
    actor
    subscription
end

as_teardown(::Type{<:DelaySubscription}) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::DelaySubscription)
    if !subscription.actor.is_cancelled
        close(subscription.actor)
    end
    subscription.actor.is_cancelled = true
    unsubscribe!(subscription.subscription)
    return nothing
end

Base.show(io::IO, ::DelayOperator)              = print(io, "DelayOperator()")
Base.show(io::IO, ::DelayProxy)                 = print(io, "DelayProxy()")
Base.show(io::IO, ::DelayActor{L})      where L = print(io, "DelayActor($L)")
Base.show(io::IO, ::DelayObservable{L}) where L = print(io, "DelayObservable($L)")
Base.show(io::IO, ::DelaySubscription)          = print(io, "DelaySubscription()")


text:
#=

function T_scale(model::GERG2008,x=SA[1.])
    return mixing_rule_asymetric(
        (a,b)->sqrt(a * b),
        _gerg_asymetric_mix_rule,
        x,
        model.Tc,
        model.gamma_T,
        model.beta_T,
    )
end

# function T_scale(model::GERG2008)
#     return only(model.Tc)
# end

function T_scales(model::GERG2008,x=SA[1.])
    return model.Tc
end

# function p_scale(model::GERG2008)
#     return only(model.Pc)
# end

function p_scale(model::GERG2008,x=SA[1.])
    return dot(x,model.Pc)
end

# function p_scales(model::GERG2008,x=SA[1.])
#     return only(model.Pc)
# end

function _v_scale(model::GERG2008,x=SA[1.])
    return mixing_rule_asymetric(
        (a,b) -> ((cbrt(a) + cbrt(b))*0.5)^3,
        _gerg_asymetric_mix_rule,
        x,
        model.vc,
        model.gamma_v,
        model.beta_v,
    )
end

# function _v_scale(model::GERG2008)
#     return only(model.vc)
# end

# function lb_volume(model::GERG2008)
#     return only(model.lb_v)
# end

function lb_volume(model::GERG2008,x=SA[1.];phase=:unknown)
    return dot(x,model.lb_v)
end

function _delta(model::GERG2008, rho, T, x=SA[1.])
    vcmix = _v_scale(model,x)
    return rho * vcmix
end

function _tau(model::GERG2008, rho, T, x=SA[1.])
    Tcmix  = T_scale(model,x)
    return Tcmix / T
end

# function _delta(model::GERG2008, rho, T)
#     vcmix = _v_scale(model)
#     return rho * vcmix
# end

# function _tau(model::GERG2008, rho, T)
#     Tcmix  = T_scale(model)
#     return Tcmix / T
# end

function _f0(model::GERG2008, ρ, T, x)
    RR = 8.314472 / 8.314510
    common_type = promote_type(typeof(ρ), typeof(T), eltype(x))
    _0 = zero(common_type)
    res = _0
    x0 = _0  #for comparison
    ao2 = _0
    ao_zero = _0

    ρc = model.rhoc
    Tc = model.Tc
    nr = model.nr
    zeta = model.zeta
    for i in model.ideal_iters[1]
        x[i] != x0 && begin
            ao2 = ao_zero
            ao3 = ao_zero
            δ = ρ / ρc[i]
            τ = Tc[i] / T
            ao1 = nr[1, i] + nr[2, i] * τ + nr[3, i] * log(τ)
            ao2 =
            nr[4, i] * log(abs(sinh(zeta[1, i] * τ))) -
            nr[5, i] * log(cosh(zeta[2, i] * τ)) +
            nr[6, i] * log(abs(sinh(zeta[3, i] * τ))) -
            nr[7, i] * log(cosh(zeta[4, i] * τ))

            ao3 = log(δ)
            a0 = RR * (ao1 + ao2) + ao3
            res += x[i] * (a0 + log(x[i]))
        end
    end

    for i in model.ideal_iters[2]
        x[i] != x0 && begin
            ao2 = ao_zero
            ao3 = ao_zero
            δ = ρ / ρc[i]
            τ = Tc[i] / T
            ao1 = nr[1, i] + nr[2, i] * τ + nr[3, i] * log(τ)
            ao2 =
                nr[4, i] * log(abs(sinh(zeta[1, i] * τ))) -
                nr[5, i] * log(cosh(zeta[2, i] * τ)) +
                nr[6, i] * log(abs(sinh(zeta[3, i] * τ)))
            ao3 = log(δ)
            a0 = RR * (ao1 + ao2) + ao3
            res += x[i] * (a0 + log(x[i]))
        end
    end

    for i in model.ideal_iters[3]
        x[i] != x0 && begin
            ao2 = ao_zero
            ao3 = ao_zero
            δ = ρ / ρc[i]
            τ = Tc[i] / T
            ao1 = nr[1, i] + nr[2, i] * τ + nr[3, i] * log(τ)
            ao2 =
                nr[4, i] * log(abs(sinh(zeta[1, i] * τ))) -
                nr[5, i] * log(cosh(zeta[2, i] * τ))
            ao3 = log(δ)
            a0 = RR * (ao1 + ao2) + ao3
            res += x[i] * (a0 + log(x[i]))
        end
    end

    for i in model.ideal_iters[4]
        x[i] != x0 && begin
            ao2 = ao_zero
            ao3 = ao_zero
            δ = ρ / ρc[i]
            τ = Tc[i] / T
            ao1 = nr[1, i] + nr[2, i] * τ + nr[3, i] * log(τ)
            ao3 = log(δ)
            a0 = RR * (ao1) + ao3
            res += x[i] * (a0 + log(x[i]))
        end
    end

    return res
end

function _f0(model::GERG2008, rho, T)

    RR = 8.314472 / 8.314510
    common_type = promote_type(typeof(rho), typeof(T))

    res = zero(common_type)
    ao1 = zero(common_type)
    ao2 = zero(common_type)
    ao_zero = zero(common_type)

    delta = rho / model.rhoc[1]
    tau = model.Tc[1] / T

    if !iszero(length(model.ideal_iters[1]))
        ao1 = model.nr[1, 1] + model.nr[2, 1] * tau + model.nr[3, 1] * log(tau)
        ao2 =
            model.nr[4, 1] * log(abs(sinh(model.zeta[1, 1] * tau))) -
            model.nr[5, 1] * log(cosh(model.zeta[2, 1] * tau)) +
            model.nr[6, 1] * log(abs(sinh(model.zeta[3, 1] * tau))) -
            model.nr[7, 1] * log(cosh(model.zeta[4, 1] * tau))

        ao3 = log(delta)
        a0 = RR * (ao1 + ao2) + ao3
        res +=  a0
    end

    if !iszero(length(model.ideal_iters[2]))
        ao1 = model.nr[1, 1] + model.nr[2, 1] * tau + model.nr[3, 1] * log(tau)
        ao2 =
            model.nr[4, 1] * log(abs(sinh(model.zeta[1, 1] * tau))) -
            model.nr[5, 1] * log(cosh(model.zeta[2, 1] * tau)) +
            model.nr[6, 1] * log(abs(sinh(model.zeta[3, 1] * tau)))
        ao3 = log(delta)
        a0 = RR * (ao1 + ao2) + ao3
        res +=  a0
    end

    if !iszero(length(model.ideal_iters[3]))
        ao1 = model.nr[1, 1] + model.nr[2, 1] * tau + model.nr[3, 1] * log(tau)
        ao2 =
            model.nr[4, 1] * log(abs(sinh(model.zeta[1, 1] * tau))) -
            model.nr[5, 1] * log(cosh(model.zeta[2, 1] * tau))
        ao3 = log(delta)
        a0 = RR * (ao1 + ao2) + ao3
        res +=  a0
    end

    if !iszero(length(model.ideal_iters[4]))
        ao1 = model.nr[1, 1] + model.nr[2, 1] * tau + model.nr[3, 1] * log(tau)
        ao3 = log(delta)
        a0 = RR * (ao1) + ao3
        res +=  a0
    end

    return res
end

function _fr1(model::GERG2008, delta, tau)
    common_type = promote_type(typeof(delta), typeof(tau))
    res = zero(common_type)
    res0 = zero(common_type)
    res1 = zero(common_type)

    res1 = res0
    for k = 1:model.k_pol_ik[1]
        res1 += model.n0ik[1][k] * (delta^model.d0ik[1][k]) * (tau^model.t0ik[1][k])
    end

    for k = (model.k_pol_ik[1]+1):(model.k_exp_ik[1]+model.k_pol_ik[1])
        res1 +=
            model.n0ik[1][k] *
            (delta^model.d0ik[1][k]) *
            (tau^model.t0ik[1][k]) *
            exp(-delta^model.c0ik[1][k-model.k_pol_ik[1]])
    end
    res += res1
    return res
end

function _fr1(model::GERG2008, delta, tau,x)
    common_type = promote_type(typeof(delta), typeof(tau), eltype(x))
    res = zero(common_type)
    res0 = zero(common_type)
    res1 = zero(common_type)
    x0 = zero(common_type)

    for i = 1:length(model.components)
        x[i] != x0 && begin
            res1 = res0
            for k = 1:model.k_pol_ik[i]
                res1 += model.n0ik[i][k] * (delta^model.d0ik[i][k]) * (tau^model.t0ik[i][k])
            end

            for k = (model.k_pol_ik[i]+1):(model.k_exp_ik[i]+model.k_pol_ik[i])
                res1 +=
                    model.n0ik[i][k] *
                    (delta^model.d0ik[i][k]) *
                    (tau^model.t0ik[i][k]) *
                    exp(-delta^model.c0ik[i][k-model.k_pol_ik[i]])
            end
            res += x[i] * res1

        end
    end

    return res
end

function _fr2(model::GERG2008, delta, tau, x)
    common_type = promote_type(typeof(delta), typeof(tau), eltype(x))
    res = zero(common_type)
    res0 = zero(common_type)
    x0 = zero(common_type)

    for kk = 1:length(model.Aij_indices)
        i1, i2, i0 = model.Aij_indices[kk]
        x[i1] != x0 && x[i2] != x0 && begin
            res1 = res0
            for j = 1:model.k_pol_ijk[i0]
                res1 +=
                    model.nijk[i0][j] * (delta^model.dijk[i0][j]) * (tau^model.tijk[i0][j])
            end

            for j = (model.k_pol_ijk[i0]+1):(model.k_pol_ijk[i0]+model.k_exp_ijk[i0])

                idx = j - model.k_pol_ijk[i0]

                res1 +=
                    model.nijk[i0][j] *
                    (delta^model.dijk[i0][j]) *
                    (tau^model.tijk[i0][j]) *
                    exp(
                        -model.etaijk[i0][idx] * (delta - model.epsijk[i0][idx])^2 -
                        model.betaijk[i0][idx] * (delta - model.gammaijk[i0][idx]),
                    )
            end
            res += res1 * x[i1] * x[i2] * model.fij[i0]

        end
    end
    return res
end

function a_ideal(model::GERG2008, v, T, z=SA[1.0])
    N = sum(z)
    len = length(z)
    x = z/N
    rho = 1.0e-3 / (v/N)
    if len == 1
        return _f0(model, rho, T)
    else
        return _f0(model, rho, T, x)
    end
end

function a_res(model::GERG2008, v, T, z=SA[1.0])
    N = sum(z)
    len = length(z)
    rho = 1.0e-3 / (v/N)
    x = z/N
    delta = _delta(model, rho, T, x)
    tau = _tau(model, rho, T, x)
    if len == 1
        return _fr1(model, delta, tau)
    else
        return (_fr1(model, delta, tau, x) + _fr2(model, delta, tau, x))
    end
end

function eos(model::GERG2008, v, T, z=SA[1.0])
    negative_vt(v,T) && return nan_num(v,T,z)
    N = sum(z)
    len = length(z)
    rho = 1.0e-3 / (v/N)
    x = z/N
    R = R̄
    delta = _delta(model, rho, T, x)
    tau = _tau(model, rho, T, x)
    if len == 1
        return N * R * T * (_f0(model, rho, T) + _fr1(model, delta, tau))
    else
        return N * R * T * (_f0(model, rho, T, x)+ _fr1(model, delta, tau, x) + _fr2(model, delta, tau, x))
    end
end
=#


text:
to_observable(o::Observable) = o
to_observable(o::AbstractObservable) = Observables.observe(o)
to_observable(o) = Observable(o)

const JSONContext = JSON.Writer.StructuralContext
const JSONSerialization = JSON.Serializations.CommonSerialization

struct MithrilSerialization <: JSONSerialization end

JSON.show_json(io::JSONContext, ::MithrilSerialization, x::JSString) = print(io, x)

printjs(io::IO, val) = JSON.show_json(io, MithrilSerialization(), val)

print_component(io::IO, x) = printjs(io, x)
print_component(io::IO, f::Function) = f(io)
print_component(io::IO, s::JSString) = print(io, s)

function print_component(io::IO, n::Node)
    print(io, "m(")
    printjs(io, instanceof(n).tag)
    print(io, ", {")
    properties = props(n)
    attributes = get(properties, :attributes, Dict())
    events = get(properties, :events, Dict())
    pairs = Iterators.flatten([attributes, properties, ("on$key" => val for (key, val) in events)])
    for (key, val) in pairs
        key in (:attributes, :events) && continue
        printjs(io, key)
        print(io, ": ")
        printjs(io, val)
        print(io, ",")
    end
    print(io, "}")
    print(io, ", [")
    for child in children(n)
        print_component(io, child)
        print(io, ", ")
    end
    print(io, "])")
end

struct MithrilComponent{T}
    template::T
    data::NamedTuple
    function MithrilComponent(template::T, data::NamedTuple) where {T}
        new{T}(template, map(to_observable, data))
    end
end

template(m::MithrilComponent) = getfield(m, :template)
data(m::MithrilComponent) = getfield(m, :data)
is_vnode(m::MithrilComponent) = getfield(m, :is_vnode)

Base.propertynames(m::MithrilComponent) = propertynames(data(m))
Base.getproperty(m::MithrilComponent, s::Symbol) = getproperty(data(m), s)[]
function Base.setproperty!(m::MithrilComponent, s::Symbol, val)
    setindex!(getproperty(data(m), s), val)
end

WebIO.render(m::MithrilComponent) = WebIO.render(Scope(m))

function print_component(io::IO, m::MithrilComponent{<:Node})
    print(io, "{view: () => ")
    print_component(io, template(m))
    print(io, "}")
end

print_component(io::IO, m::MithrilComponent) = print_component(io, template(m))

function WebIO.Scope(m::MithrilComponent)
    s = Scope(imports =
        [
         "https://unpkg.com/mithril@next/mithril.js",
         "https://www.gitcdn.xyz/repo/piever/InteractResources/v0.4.0/bulma/main_confined.min.css"
        ]
    )
    datanames = String[]
    for (key, val) in pairs(data(m))
        skey = string(key)
        setobservable!(s, skey, to_observable(val))
        push!(datanames, skey)
        onjs(s[skey], js"function (value) {this.m.redraw()}")
    end
    io = IOBuffer()
    print(io, js"""
    function (m) {
        this.m = m;
        function addProperty(obj, name) {
            Object.defineProperty(obj, name, {
                get: function() {return _webIOScope.getObservableValue(name);},
                set: function(val) {_webIOScope.setObservableValue(name, val);}
            });
        }
        function Data(names) {
            var key;
            for (key of names) {addProperty(this, key);}
        }
        var data = new Data($datanames);
    """)
    print(io, "var template = ")
    print_component(io, m)
    print(io, ";\n")
    print(io, "m.mount(_webIOScope.dom, {view: function () {return m('div.interact-widget', m(template));}});\n}")
    onimport(s, JSString(String(take!(io))))
    return s
end

struct MithrilWidget{T, S}<:AbstractWidget{T, S}
    m::MithrilComponent
    observe::AbstractObservable{S}
    function MithrilWidget{T}(m::MithrilComponent, observe::AbstractObservable{S}) where {T, S}
        new{T, S}(m, observe)
    end
end

component(mw::MithrilWidget) = getfield(mw, 1)
Observables.observe(mw::MithrilWidget) = getfield(mw, 2)

Base.propertynames(mw::MithrilWidget) = propertynames(component(mw))
Base.getproperty(mw::MithrilWidget, s::Symbol) = getproperty(component(mw), s)
Base.setproperty!(mw::MithrilWidget, s::Symbol, val) = setproperty!(component(mw), s, val)
Widgets.render(mw::MithrilWidget) = WebIO.render(component(mw))


text:
using NURBS
using Base.Test

@testset "B-Spline curve generator" begin

    @testset "bspline" begin
        b = [1. 2 4 3; 1 3 3 1; 0 0 0 0]
        @test typeof(bspline(4,4,5,b)[1]) == Array{Float64,2}
        @test typeof(bspline(4,4,5,b)[2]) == Array{Array{Int64,1},1}
        @test isapprox(bspline(4,4,5,b)[1],[1.0 1.875 2.75 3.25 3.0; 1.0 2.125 2.5 2.125 1.0; 0.0 0.0 0.0 0.0 0.0],atol=1e-5)
    end

    @testset "bsplineu" begin
        b = [0. 3 6 9; 0 10 3 6; 0 0 0 0]
        @test typeof(bsplineu(4,3,5,b)[1]) == Array{Float64,2}
        @test typeof(bsplineu(4,3,5,b)[2]) == Array{Array{Int64,1},1}
        @test isapprox(bsplineu(4,3,5,b)[1],[1.5 3.0 4.5 6.0 7.5; 5.0 7.875 6.5 4.25 4.5; 0.0 0.0 0.0 0.0 0.0],atol=1e-5)
    end

    @testset "dbspline" begin
        b = [1. 2 4 3; 1 3 3 1; 0 0 0 0]
        @test typeof(dbspline(4,4,5,b)[1]) == Array{Float64,2}
        @test typeof(dbspline(4,4,5,b)[2]) == Array{Float64,2}
        @test typeof(dbspline(4,4,5,b)[3]) == Array{Float64,2}
        @test typeof(dbspline(4,4,5,b)[4]) == Array{Array{Int64,1},1}
        @test isapprox(dbspline(4,4,5,b)[1],[1.0 1.875 2.75 3.25 3.0; 1.0 2.125 2.5 2.125 1.0; 0.0 0.0 0.0 0.0 0.0],atol=1e-5)
        @test isapprox(dbspline(4,4,5,b)[2],[3.0 3.75 3.0 0.75 -3.0; 6.0 3.0 0.0 -3.0 -6.0; 0.0 0.0 0.0 0.0 0.0],atol=1e-5)
        @test isapprox(dbspline(4,4,5,b)[3],[6.0 0.0 -6.0 -12.0 -18.0; -12.0 -12.0 -12.0 -12.0 -12.0; 0.0 0.0 0.0 0.0 0.0],atol=1e-5)
    end

    @testset "dbsplineu" begin
        b = [0. 3 6 9; 0 10 3 6; 0 0 0 0]
        @test typeof(dbsplineu(4,3,5,b)[1]) == Array{Float64,2}
        @test typeof(dbsplineu(4,3,5,b)[2]) == Array{Float64,2}
        @test typeof(dbsplineu(4,3,5,b)[3]) == Array{Float64,2}
        @test typeof(dbsplineu(4,3,5,b)[4]) == Array{Array{Int64,1},1}
        @test isapprox(dbsplineu(4,3,5,b)[1],[1.5 3.0 4.5 6.0 7.5; 5.0 7.875 6.5 4.25 4.5; 0.0 0.0 0.0 0.0 0.0],atol=1e-5)
        @test isapprox(dbsplineu(4,3,5,b)[2],[3.0 3.0 3.0 3.0 3.0; 10.0 1.5 -7.0 -2.0 3.0; 0.0 0.0 0.0 0.0 0.0],atol=1e-5)
        @test isapprox(dbsplineu(4,3,5,b)[3],[0.0 0.0 0.0 0.0 0.0; -17.0 -17.0 10.0 10.0 10.0; 0.0 0.0 0.0 0.0 0.0],atol=1e-5)
    end

    @testset "nmatrix" begin
        @test typeof(nmatrix(3)[1]) == Array{Float64,2}
        @test typeof(nmatrix(3)[2]) == Float64
        @test nmatrix(3)[1] == [1.0 -2.0 1.0; -2.0 2.0 0.0; 1.0 1.0 0.0]
        @test nmatrix(3)[2] == 0.5
        @test nmatrix(4)[1] == [-1.0 3.0 -3.0 1.0; 3.0 -6.0 3.0 0.0; -3.0 0.0 3.0 0.0; 1.0 4.0 1.0 0.0]
        @test isapprox(nmatrix(4)[2],0.1666666666,atol=1e-10)
    end

    @testset "matpbspl" begin
        b = [0. 3 6 9; 0 10 3 6; 0 0 0 0]
        @test typeof(matpbspl(4,3,5,b)[1]) == Array{Float64,2}
        @test typeof(matpbspl(4,3,5,b)[2]) == Array{Array{Int64,1},1}
        @test isapprox(matpbspl(4,3,4,b)[1],[1.5 2.5 3.5 4.5 5.5 6.5 7.5 7.83333 6.83333 4.5 2.16667 1.16667 1.5 2.5 3.5; 5.0 7.38889 7.88889 6.5 4.72222 4.05556 4.5 5.0 4.5 3.0 1.88889 2.55556 5.0 7.38889 7.88889; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0],atol=1e-5)
    end
end


text:
import UncertainData: RandomSequences

x = rand(100)
y = rand(100)

test_ccm = ConvergentCrossMappingTest(timeseries_lengths = [45, 50], n_reps = 20)
test_cm = CrossMappingTest(n_reps = 10)
test_vf = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = 1:5)
test_tog = TransferOperatorGridTest(binning = RectangularBinning(5), ηs = 1:5)
test_jdd = JointDistanceDistributionTest()
test_jddt = JointDistanceDistributionTTest()
predtest = VisitationFrequencyTest(binning = RectangularBinning(5), ηs = -5:5)
test_pa = PredictiveAsymmetryTest(predictive_test = predtest)
test_npa = NormalisedPredictiveAsymmetryTest(predtest, f = 1.0)
test_NNmi = NearestNeighbourMITest(ηs = -5:5)

Nseq = 80
nchunks = 3
rtest_cm = RandomSequencesTest(test_cm, RandomSequences(nchunks, Nseq))
rtest_ccm = RandomSequencesTest(test_ccm, RandomSequences(nchunks, Nseq))
rtest_vf = RandomSequencesTest(test_vf, RandomSequences(nchunks, Nseq))
rtest_tog = RandomSequencesTest(test_tog, RandomSequences(nchunks, Nseq))
rtest_jdd = RandomSequencesTest(test_jdd, RandomSequences(nchunks, Nseq))
rtest_jddt = RandomSequencesTest(test_jddt, RandomSequences(nchunks, Nseq))
rtest_pa = RandomSequencesTest(test_pa, RandomSequences(nchunks, Nseq))
rtest_npa = RandomSequencesTest(test_npa, RandomSequences(nchunks, Nseq))
rtest_NNmi = RandomSequencesTest(test_NNmi, RandomSequences(nchunks, Nseq))

@test rtest_cm isa RandomSequencesTest
@test rtest_ccm isa RandomSequencesTest
@test rtest_vf isa RandomSequencesTest
@test rtest_tog isa RandomSequencesTest
@test rtest_jdd isa RandomSequencesTest
@test rtest_jddt isa RandomSequencesTest
@test rtest_pa isa RandomSequencesTest
@test rtest_npa isa RandomSequencesTest
@test rtest_NNmi isa RandomSequencesTest

import CausalityTools.CausalityTests.get_return_type
import HypothesisTests: OneSampleTTest

T = typeof(1.0)

@test causality(x, y, rtest_cm) isa Vector{Vector{T}}
@test causality(x, y, rtest_ccm) isa Vector{Vector{Vector{T}}}
@test causality(x, y, rtest_vf) isa Vector{Vector{T}}
@test causality(x, y, rtest_tog) isa Vector{Vector{T}}
@test causality(x, y, rtest_jdd) isa Vector{Vector{T}}
@test causality(x, y, rtest_jddt) isa Vector{OneSampleTTest}
@test causality(x, y, rtest_pa) isa Vector{Vector{T}}
@test causality(x, y, rtest_npa) isa Vector{Vector{T}}
@test causality(x, y, rtest_NNmi) isa Vector{Vector{T}}

@test causality(x, y, rtest_cm) |> length == 3
@test causality(x, y, rtest_ccm) |> length == 3
@test causality(x, y, rtest_vf) |> length == 3
@test causality(x, y, rtest_tog) |> length == 3
@test causality(x, y, rtest_jdd) |> length == 3
@test causality(x, y, rtest_jddt) |> length == 3
@test causality(x, y, rtest_pa) |> length == 3
@test causality(x, y, rtest_npa) |> length == 3
@test causality(x, y, rtest_NNmi) |> length == 3