text:
module VariableTemplates

debug() = false
using StaticArrays

export Vars, Grad

include("vars.jl")
include("varsindex.jl")
include("concretize.jl")
include("varsize.jl")
include("flattened_names.jl")
include("flattened_tup_chain.jl")
include("flattened_nt.jl")

abstract type AbstractVars{S, A, TC, TH} end
struct Vars{
    S,   # Variable index map
    A,   # Array
    TC,  # Tuple chain (tuple of Union{Symbol,Int})
    TH   # Tuple hash (NamedTuple): keys=hash of all tup-chains, values=(gid,len)
    } <: AbstractVars{S, A, TC, TH}
    array::A
end

struct Grad{
    S,   # Variable index map
    A,   # Array
    TC,  # Tuple chain (tuple of Union{Symbol,Int})
    TH   # Tuple hash (NamedTuple): keys=hash of all tup-chains, values=(gid,len)
    } <: AbstractVars{S, A, TC, TH}
    array::A
end

Base.parent(v::AbstractVars) = getfield(v, :array)
Base.eltype(v::AbstractVars) = eltype(parent(v))
Base.similar(v::AbstractVars{S, A, TC, TH}) where {S, A, TC, TH} =
    typeof(v)(similar(parent(v)))
Base.propertynames(v::AbstractVars{S}) where {S} = propertynames(S)

function Vars{S}(array) where {S}
    if S isa NamedTuple
        S_used = S
    else
        S_used = concretize(S)
    end
    TH = get_hash(S_used)
    A = typeof(array)
    TC = ()
    return Vars{S_used, A, TC, TH}(array)
end
function Grad{S}(array) where {S}
    if S isa NamedTuple
        S_used = S
    else
        S_used = concretize(S)
    end
    TH = get_hash(S_used)
    A = typeof(array)
    TC = ()
    return Grad{S_used, A, TC, TH}(array)
end

function get_hash(S)
    return flattened(
        S,
        :_,
        (pn,p) -> pn==:gid || pn == :type,
        (pn,p)->p,
        (nt)-> :gid in propertynames(nt) || :type in propertynames(nt))
end

#####
##### Internal methods
#####

get_tup_chain(v::AbstractVars{S,A,TC}) where {S,A,TC} = TC

get_val(::Val{T}) where {T} = T
get_val(::Type{Val{T}}) where {T} = T

# get length of result
get_len(::Type{Tuple{}}) = 0
get_len(::Type{Tuple{T}}) where {T<:NamedTuple} = get_len(parameters(T)[2])
get_len(::Type{Tuple{T}}) where {T} = get_len(T)
get_len(::Type{T}) where {T<:Real} = 1
get_len(::Type{T}) where {N,FT<:Real,T<:SVector{N,FT}} = N
get_len(::Type{T}) where {N,FT,T<:SArray{N,FT}} = N
get_len(::Type{T}) where {N,FT<:Real,M,T<:SHermitianCompact{N, FT, M}} = M

get_result(a, gid, ::Type{T}, ::Vars) where {T<:Real} = T(a[gid])
get_result(a, gid, ::Type{T}, ::Vars) where {S,FT<:Real,N,L,T<:SArray{S,FT,N,L},Vars} = T(a[gid:gid+L-1])
function get_result(a, gid, ::Type{T}, ::Vars) where {N,FT<:Real,M,T<:SHermitianCompact{N, FT, M}}
    return T(a[gid:gid+M-1])
end

get_result(a, gid, ::Type{T}, g::Grad) where {T<:Real} = T(a[1:grad_dim(g), gid])
function get_result(a, gid, ::Type{T}, g::Grad) where {S,FT<:Real,N,L,T<:SArray{S,FT,N,L}}
    M = grad_dim(g)
    MMatrix{M,L,FT}(a[1:M, gid:gid+L-1])
end
function get_result(a, gid, ::Type{T}, g::Grad) where {N,FT<:Real,P,T<:SHermitianCompact{N, FT, P}}
    M = grad_dim(g)
    return T(a[1:M, gid:gid+P-1])
end

function set_result!(a, gid, ::Type{T}, val, ::Vars) where {T<:Real}
    a[gid] = T(val)
end
function set_result!(a, gid, ::Type{T}, val, ::Vars) where {S,FT<:Real,N,L,T<:SArray{S,FT,N,L}}
    a[gid:gid+L-1] .= T(val)
end
function set_result!(a, gid, ::Type{T}, val, ::Vars) where {N,FT<:Real,P,T<:SHermitianCompact{N, FT, P}}
    a[gid:gid+P-1] .= T(val).lowertriangle
end

function set_result!(a, gid, ::Type{T}, val, g::Grad) where {T<:Real}
    M = grad_dim(g)
    a[1:M, gid] = T(val)
end
function set_result!(a, gid, ::Type{T}, val, g::Grad) where {S,FT<:Real,N,L,T<:SArray{S,FT,N,L}}
    M = grad_dim(g)
    a[1:M, gid:gid+L-1] .= MMatrix{M,L,FT}(val)
end
function set_result!(a, gid, ::Type{T}, val, g::Grad) where {N,FT<:Real,P,T<:SHermitianCompact{N, FT, P}}
    M = grad_dim(g)
    a[1:M, gid:gid+P-1] .= T(val).lowertriangle
end

Base.@propagate_inbounds function getproperty_generic(
    v::AbstractVars{S, A, TC, TH},
    tup_chain,
    S_new
) where {S, A, TC, TH}
    joined_tup_chain = join_tup_chain(tup_chain, :_)
    a = parent(v)
    if hasproperty(TH, joined_tup_chain)
        nt = getproperty(TH, joined_tup_chain)
        gid = get_val(nt.gid)
        type = get_val(nt.type)
        return get_result(a, gid, type, v)
    else
        if v isa Vars
            return Vars{S_new,A,tup_chain,TH}(a)
        else
            return Grad{S_new,A,tup_chain,TH}(a)
        end
    end
end

Base.@propagate_inbounds Base.getproperty(t::Tuple, tup_chain::Tuple) =
    Base.getproperty(Base.getindex(t, tup_chain[1]), tup_chain[2:end])

Base.@propagate_inbounds Base.getproperty(nt::NamedTuple, tup_chain::Tuple) =
    Base.getproperty(Base.getproperty(nt, tup_chain[1]), tup_chain[2:end])

Base.@propagate_inbounds Base.getproperty(nt::NamedTuple, syms::Tuple{S}) where {S} =
    Base.getproperty(nt, syms[1])

update_S(S, sym::Symbol) = getproperty(S, sym)
update_S(S, i::Int) = getindex(S, i)
update_S(S, tup_chain) = getproperty(S, tup_chain)

Base.@propagate_inbounds Base.getproperty(v::AbstractVars{S}, tup_chain::Tuple) where {S} =
    getproperty_generic(v, tup_chain, update_S(S, tup_chain))
Base.@propagate_inbounds Base.getproperty(v::AbstractVars{S, A, TC}, sym::Symbol) where {S, A, TC} =
    getproperty_generic(v, (TC..., sym), update_S(S, sym))

grad_dim(v::Grad{S,A}) where {S,A<:SubArray} = size(fieldtype(A, 1), 1)
grad_dim(v::Grad{S,A}) where {S,A} = size(A, 1)

Base.@propagate_inbounds function Base.setproperty!(
    v::AbstractVars{S, A, TC, TH},
    sym::Symbol,
    val
) where {S, A, TC, TH}
    sep = :_
    new_tup_chain = (TC..., sym)
    joined_tup_chain = join_tup_chain(new_tup_chain, sep)
    a = parent(v)
    if hasproperty(TH, joined_tup_chain)
        nt = getproperty(TH, joined_tup_chain)
        gid = get_val(nt.gid)
        type = get_val(nt.type)
        return set_result!(a, gid, type, val, v)
    else
        error("Cannot set property of non-concrete variable")
    end
end

Base.@propagate_inbounds @inline function Base.getindex(
    v::Vars{S, A, TC, TH},
    i::Int,
) where {S, A, TC, TH}
    sep = :_
    new_tup_chain = (TC..., i)
    joined_tup_chain = join_tup_chain(new_tup_chain, sep)
    a = parent(v)
    if hasproperty(TH, joined_tup_chain)
        nt = getproperty(TH, joined_tup_chain)
        gid = get_val(nt.gid)
        type = get_val(nt.type)
        return get_result(a, gid, type, v)
    else
        S_new = update_S(S, i)
        return Vars{S_new,A,new_tup_chain,TH}(a)
    end
end

Base.@propagate_inbounds @inline function Base.setindex!(
    v::AbstractVars{S, A, TC, TH},
    i::Int,
    val,
) where {S, A, TC, TH}
    sep = :_
    new_tup_chain = (TC..., i)
    joined_tup_chain = join_tup_chain(new_tup_chain, sep)
    a = parent(v)
    if hasproperty(TH, joined_tup_chain)
        nt = getproperty(TH, joined_tup_chain)
        gid = get_val(nt.gid)
        type = get_val(nt.type)
        return set_result!(a, gid, type, val, v)
    else
        error("Cannot set index of non-concrete variable")
    end
end

end # module


text:
# all the packages used in Julia
Pkg.update()
Pkg.add("MathProgBase")
Pkg.add("Mosek")
Pkg.add("ECOS")
Pkg.add("SCS")
Pkg.add("Ipopt")
Pkg.add("KNITRO")
Pkg.add("NLopt")
Pkg.add("Cuba")
Pkg.add("CPUTime")
Pkg.add("JSON")
Pkg.update()
#EOF

ext: jl

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
using Random
using RecipesBase

import MinimalRLCore
using Distributions: Uniform

module OneDTmazeConst

const UP = 1
const RIGHT = 2
const DOWN = 3
const LEFT = 4

const ACTIONS = [UP, RIGHT, DOWN, LEFT]

const EPSILON = 0.05
const ACTION_STEP = 0.08
const MOVE_NOISE = 0.01

end

"""
    1d TMaze

"""
Base.@kwdef struct OneDTMaze <: MinimalRLCore.AbstractEnvironment

    pos::Vector{Float64} = [0.0, 0.0]

    cumulant_schedule::CumulantSchedule = TMazeCumulantSchedules.Constant(1.0)
    starts::String = "beg"
    env_reward::Float64 = -0.01
    move_noise::Float64 = OneDTmazeConst.MOVE_NOISE

end

OneDTMaze(starts::String) = OneDTMaze(starts=starts)
OneDTMaze(cs::CumulantSchedule, starts::String) = OneDTMaze(cumulant_schedule=cs, starts=starts)
OneDTMaze(cs::CumulantSchedule, starts::String, r::Float64) = OneDTMaze(cumulant_schedule=cs, starts=starts, env_reward = r)
OneDTMaze(cs::CumulantSchedule, starts::String, r::Float64, move_noise::Float64) = OneDTMaze(cumulant_schedule=cs, starts=starts, env_reward = r, move_noise = move_noise)

volume(::OneDTMaze) = 0.8 + 1.0 + 2*0.4

MinimalRLCore.get_actions(env::OneDTMaze) = OneDTmazeConst.ACTIONS
valid_action(env::OneDTMaze, action) = action in MinimalRLCore.get_actions(env)
MinimalRLCore.get_reward(env::OneDTMaze) = env.env_reward # -> determines if the agent_state is terminal
MinimalRLCore.is_terminal(env::OneDTMaze, pos=env.pos) = any(check_goal.([env], 1:4, [pos]))

check_goal(env::OneDTMaze, goal, pos=env.pos) = check_goal(OneDTMaze, goal, pos)

function check_goal(env::Type{OneDTMaze}, goal, pos, epsilon=OneDTmazeConst.EPSILON)
    cur_x = pos[1]
    cur_y = pos[2]
    if goal == (goal isa String ? "G1" : 1)
        (cur_x == 0.0 && range_check(cur_y, 1.0-epsilon, 2.0)) # G1
    elseif goal == (goal isa String ? "G2" : 2)
        (cur_x == 0.0 && range_check(cur_y, 0.0, 0.6+epsilon)) # G2
    elseif goal == (goal isa String ? "G3" : 3)
        (cur_x == 1.0 && range_check(cur_y, 1.0-epsilon, 2.0)) # G3
    elseif goal == (goal isa String ? "G4" : 4)
        (cur_x == 1.0 && range_check(cur_y, 0.0, 0.6+epsilon)) # G4
    end
end

function get_cumulants(env::OneDTMaze, cs::CumulantSchedule)
    cumulants = zeros(4)
    if check_goal(env, 1)
        cumulants[1] = get_cumulant(cs, "G1")
    elseif check_goal(env, 2)
        cumulants[2] = get_cumulant(cs, "G2")
    elseif check_goal(env, 3)
        cumulants[3] = get_cumulant(cs, "G3")
    elseif check_goal(env, 4)
        cumulants[4] = get_cumulant(cs, "G4")
    end
    cumulants
end

function MinimalRLCore.get_state(env::OneDTMaze)
    cumulants = get_cumulants(env, env.cumulant_schedule)
    return vcat(env.pos, cumulants)
end

function get_random_state(env::OneDTMaze)

    if env.starts == "beg"
        [0.5, rand(Uniform(0.0, 0.1))]
    elseif env.starts == "const"
        [0.5, 0.0]
    elseif env.starts == "whole"
        h = rand()
        if h < 0.8/volume(env)
            [0.5, rand(Uniform(0.0, 0.8))]
        elseif h < 1.8/volume(env)
            [rand(Uniform(0.0, 1.0)), 0.8]
        else
            if rand() > 0.5
                [0.0, rand(Uniform(0.65, 0.95))]
            else
                [1.0, rand(Uniform(0.65, 0.95))]
            end
        end
    end
end

function MinimalRLCore.reset!(env::OneDTMaze)
    ns = get_random_state(env)
    env.pos .= ns
end

function MinimalRLCore.reset!(env::OneDTMaze,
                              start_state::Vector{Float64})
   env.pos .= start_state
end

function MinimalRLCore.environment_step!(env::OneDTMaze, action, rng::AbstractRNG=Random.GLOBAL_RNG)
    ODTMC = OneDTmazeConst
    @boundscheck valid_action(env, action)
    if env.move_noise == 0.0
        rand_mov = 0
    else
        rand_mov = rand(rng, Uniform(-env.move_noise, env.move_noise))
    end

    x_mov, y_mov = if action == ODTMC.UP
        (0.0, ODTMC.ACTION_STEP + rand_mov)
    elseif action == ODTMC.DOWN
        (0.0, -ODTMC.ACTION_STEP + rand_mov)
    elseif action == ODTMC.RIGHT
        (ODTMC.ACTION_STEP + rand_mov, 0.0)
    elseif action == ODTMC.LEFT
        (-ODTMC.ACTION_STEP + rand_mov, 0.0)
    end

    cur_x = env.pos[1]
    cur_y = env.pos[2]

    # checking for hallways
    if (cur_x == 0.0 && cur_y == 0.8)
        # In left junction
        # Can't move left
        if (x_mov < 0)
            x_mov = 0
        end
    elseif (cur_x == 0.5 && cur_y == 0.8)
        # In middle junction
        # Can't move up
        if (y_mov > 0)
            y_mov = 0
        end
    elseif (cur_x == 1.0 && cur_y == 0.8)
        # in right junction
        # can't move right
        if (x_mov > 0)
            x_mov = 0
        end
    elseif (cur_x == 0.5 || cur_x == 0.0 || cur_x == 1.0)
        # in one of the vertical hallways
        # can't move horizontally
        x_mov = 0
    elseif (cur_y == 0.8)
        # in horizontal hallway
        # can't move vertically
        y_mov = 0
    else
        println("SOMETHING IS BROKEN WITH 1D TMaze")
    end

    # Moving
    cur_x += x_mov
    cur_y += y_mov

    # clamping between known bounds
    cur_x = clamp(cur_x, 0.0, 1.0)
    cur_y = clamp(cur_y, 0.0, 1.0)

    # snapping to junctions
    if (range_check(cur_x, 0.0, 0.0 + ODTMC.EPSILON)
        && range_check(cur_y, 0.8 - ODTMC.EPSILON, 0.8 + ODTMC.EPSILON))
        # Snapping to left junction
        cur_x, cur_y = 0.0, 0.8
    elseif (range_check(cur_x, 0.5 - ODTMC.EPSILON, 0.5 + ODTMC.EPSILON)
        && range_check(cur_y, 0.8 - ODTMC.EPSILON , 0.8 + ODTMC.ACTION_STEP + env.move_noise))
        # second condition checks for whether the agent somehow makes too big of a step upwards.
        # In middle junction
        # Snapping to middle junction
        cur_x, cur_y = 0.5, 0.8
    elseif (range_check(cur_x, 1.0 - ODTMC.EPSILON, 1.0)
        && range_check(cur_y, 0.8 - ODTMC.EPSILON, 0.8 + ODTMC.EPSILON))
        # in right junction
        # Snapping to right junction
        cur_x, cur_y = 1.0, 0.8
    end

    env.pos[1] = cur_x
    env.pos[2] = cur_y

    update!(env.cumulant_schedule, env.pos)

end

module OneDTMazePlotParams

using Colors
import ColorSchemes

const SIZE = 10
const BG = Colors.RGB(1.0, 1.0, 1.0)
const WALL = Colors.RGB(0.3, 0.3, 0.3)
const AC = Colors.RGB(0.69921875, 0.10546875, 0.10546875)
const GOAL = Colors.RGB(0.796875, 0.984375, 0.76953125)
const GOAL_PALETTE = ColorSchemes.tol_muted
const AGENT = fill(AC, SIZE, SIZE)

end

@recipe function f(env::OneDTMaze)
    ticks := nothing
    foreground_color_border := nothing
    grid := false
    legend := false
    aspect_ratio := 1
    xaxis := false
    yaxis := false

    PP = OneDTMazePlotParams

    screen = fill(PP.WALL, PP.SIZE*12, PP.SIZE*12)

    # Paths
    screen[Int(PP.SIZE-PP.SIZE//2+1):Int(9*PP.SIZE + PP.SIZE//2), (Int(11//2*PP.SIZE) + 1):(Int((13//2)*PP.SIZE))] .= PP.BG
    screen[(Int(17//2*PP.SIZE) + 1):Int(19//2*PP.SIZE), Int(PP.SIZE-PP.SIZE//2+1):Int(11*PP.SIZE + PP.SIZE//2)] .= PP.BG

    screen[Int(14//2*PP.SIZE + 1):Int(22//2*PP.SIZE), Int(1//2*PP.SIZE+1):Int(3//2*PP.SIZE)] .= PP.BG
    screen[Int(14//2*PP.SIZE + 1):Int(22//2*PP.SIZE), Int(21//2*PP.SIZE+1):Int(23//2*PP.SIZE)] .= PP.BG

    # GOALS
    ϵ = OneDTmazeConst.EPSILON
    if OneDTmazeConst.EPSILON == 0.0
        screen[Int(13//2*PP.SIZE+1):Int(14//2*PP.SIZE), Int(1//2*PP.SIZE+1):Int(3//2*PP.SIZE)] .= PP.GOAL_PALETTE[2]
        screen[Int(22//2*PP.SIZE+1):Int(23//2*PP.SIZE), Int(1//2*PP.SIZE+1):Int(3//2*PP.SIZE)] .= PP.GOAL_PALETTE[1]
        screen[Int(13//2*PP.SIZE+1):Int(14//2*PP.SIZE), Int(21//2*PP.SIZE+1):Int(23//2*PP.SIZE)] .= PP.GOAL_PALETTE[4]
        screen[Int(22//2*PP.SIZE+1):Int(23//2*PP.SIZE), Int(21//2*PP.SIZE+1):Int(23//2*PP.SIZE)] .= PP.GOAL_PALETTE[3]
    else
        ϵ_min_size = maximum([Int(floor(PP.SIZE*ϵ)),1])
        screen[Int(13//2*PP.SIZE+1):Int(14//2*PP.SIZE+ϵ_min_size), Int(1//2*PP.SIZE+1):Int(3//2*PP.SIZE)] .= PP.GOAL_PALETTE[2]
        screen[Int(22//2*PP.SIZE-ϵ_min_size+1):Int(23//2*PP.SIZE), Int(1//2*PP.SIZE+1):Int(3//2*PP.SIZE)] .= PP.GOAL_PALETTE[1]
        screen[Int(13//2*PP.SIZE+1):Int(14//2*PP.SIZE+ϵ_min_size), Int(21//2*PP.SIZE+1):Int(23//2*PP.SIZE)] .= PP.GOAL_PALETTE[4]
        screen[Int(22//2*PP.SIZE-ϵ_min_size+1):Int(23//2*PP.SIZE), Int(21//2*PP.SIZE+1):Int(23//2*PP.SIZE)] .= PP.GOAL_PALETTE[3]
    end

    state = env.pos
    x, y = Int(floor((env.pos[1]*10 + 1)*PP.SIZE)), Int(floor((env.pos[2]*10 + 1)*PP.SIZE))

    screen[Int(y-PP.SIZE//5 + 1):Int(y+PP.SIZE//5), Int(x-PP.SIZE//2 + 1):Int(x+PP.SIZE//2)] .= PP.AC
    screen[Int(y-PP.SIZE//2 + 1):Int(y+PP.SIZE//2), Int(x-PP.SIZE//5 + 1):Int(x+PP.SIZE//5)] .= PP.AC

    screen[end:-1:1, :]
end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
include(joinpath(@__DIR__, "local_optimizer.jl"))
if !isdefined(Main, :FCFP)
    include(joinpath(@__DIR__, "FCFP.jl"))
end

module DSSP

    using Printf
    using JuMP
    using Clp
    using Random
    import Main.FCFP.FCFPInstance
    import Main.LocalOptimizer.optimize_locally

    # Debug Data
    #= network_objective = 450
    V = [1,2,3,4,5,6,7]
    E = [(1,2,0), (1,3,0), (1,4,0), (2,3,0), (4,3,0),(2,6,0),(2,7,0),(3,7,0),(4,5,0)]
    production_nodes = [1]
    building_nodes = [5,6,7]
    junctions_nodes = [2,3,4]
    costs = Dict(e => 1 for e=E)
    heat_demand = Dict(v => 200 for v=building_nodes)
    production = Dict(v => 500 for v=production_nodes)
    capacities = Dict{Tuple, Float32}() =#

    function edge_cost_func(e::Tuple, inst::FCFPInstance)
        if e in keys(inst.fixed_costs)
            return inst.fixed_costs[e]
        else
            return inst.fixed_costs[e[2], e[1], e[3]]
        end
    end

    function write_optimization_LP_model(inst::FCFPInstance)
        ### Create model ###
        m = Model(with_optimizer(Clp.Optimizer, LogLevel=0))
        ### Variables
        @variable(m, 0 <= x[inst.all_edges] <= inst.network_objective)
        # Capacities of old pipes
        for e=keys(inst.capacities)
            setupperbound(x[e], inst.capacities[e])
        end
        ### Constraints
        # Source objective flow:
        @constraint(m, sum(x[e] for e=inst.all_edges if e[1] in inst.production_nodes) >= inst.network_objective)
        # Well objective flow:
        @constraint(m, sum(x[e] for e=inst.all_edges if e[2] in inst.building_nodes) >= inst.network_objective)
        # Each node follows the nodal rule, i.e. sum(flow_in) = sum(flow_out)
        for n=inst.building_nodes
            @constraint(m, sum(x[e] for e=inst.all_edges if e[2]==n) <= inst.heat_demand[n])
        end
        for n=inst.old_buildings
            @constraint(m, sum(x[e] for e=inst.all_edges if e[2]==n) >= inst.heat_demand[n]-1e-5)
        end
        for n=inst.production_nodes
            @constraint(m, sum(x[e] for e=inst.all_edges if e[1]==n) <= inst.production[n])
        end
        for n=inst.junctions_nodes
            @constraint(m, sum(x[e] for e=inst.all_edges if e[2]==n) == sum(x[e] for e=inst.all_edges if e[1]==n))
        end
        ### Objective
        @objective(m, Min, sum(x[e] * edge_cost_func(e, inst) for e=inst.all_edges))

        return m, x
    end

    function update_costs(inst::FCFPInstance, current_solution, last_solution, edge_costs, history_costs, tol, schema = "II" )
            nb_modified = 0
            for e=inst.all_edges
                    flow = current_solution[e]
                    if abs(flow - last_solution[e]) > tol
                            nb_modified += 1
                    end
                    if flow > 0
                            new_cost = edge_cost_func(e, inst) / flow
                            edge_costs[e] = new_cost
                            if schema == "I"
                                    history_costs[e] = max(new_cost, history_costs[e])
                            else
                                    history_costs[e] = new_cost
                            end
                    else
                            edge_costs[e] = history_costs[e]
                    end
            end
            return edge_costs, history_costs, nb_modified
    end

    function disrupt_costs(inst::FCFPInstance, edge_costs)
            for e=inst.all_edges
                    edge_costs[e] = max(0, edge_costs[e] + randn() * edge_costs[e]/2)
            end
            return edge_costs
    end

    function optimize_with_DSSP(network_objective, V, costs, heat_demand,
                                production, capacities, old_buildings, print_function = print,
                                postprocess=true, schema="II", maxiter=1000,
                                tol=1e-2)
        simulated_annealing = false
        simulated_annealing_max_iterations = 30
        simulated_annealing_it = 0
        Random.seed!(0)
        max_iter_without_improvement = 20
        inst = FCFPInstance(network_objective, V, costs, heat_demand, production, capacities, old_buildings)
        m, x = write_optimization_LP_model(inst)
        edge_costs =  Dict{Tuple , Float32}()
        for e=inst.all_edges
                if e in keys(capacities)
                        edge_costs[e] =  0 + edge_cost_func(e, inst) / capacities[e]
                else
                        edge_costs[e] =  0
                end
        end
        current_solution = Dict(e => 0 for e=inst.all_edges)
        history_costs = copy(edge_costs)
        best_solution = copy(current_solution)
        best_cost = Inf
        print_function("  Iteration  |  Current Obj.  |  Real Obj.  |  Best Obj.  |  Nb modifed Cost  | Solve Status\n")
        print_function(string(repeat("_", 92), "\n"))
        counter = 0
        nb_iter_without_improvement = 0
        while true
                last_solution = copy(current_solution)
                # Set new objective
                @objective(m, Min, sum(x[e] * edge_costs[e] for e=inst.all_edges))
                # Solve current problem
                optimize!(m)
                status = termination_status(m)
                current_solution = Dict(e => value(x[e]) for e=inst.all_edges)
                current_cost = objective_value(m)
                # Update costs
                edge_costs, history_costs, nb_modified = update_costs(inst, current_solution, last_solution, edge_costs, history_costs, tol, schema)
                nb_iter_without_improvement += 1
                # Save best solution
                real_current_cost = sum(edge_cost_func(e, inst) for e=inst.all_edges if current_solution[e] > 0)
                if real_current_cost < best_cost
                        best_solution = copy(current_solution)
                        best_cost = real_current_cost
                        nb_iter_without_improvement = 0
                end
                # Print information
                print_function("$(@sprintf("%13d|%16f|%13f|%13f|%19d|%14s", counter, current_cost, real_current_cost, best_cost, nb_modified, status))")
                counter += 1
                # Stopping conditions
                converged = all(abs(current_solution[e] - last_solution[e]) < tol for e=inst.all_edges)
                converged = converged || (nb_iter_without_improvement > max_iter_without_improvement)
                if converged
                        if simulated_annealing && (simulated_annealing_it < simulated_annealing_max_iterations)
                                print_function("=== Disrupting step ===")
                                edge_costs = disrupt_costs(inst, edge_costs)
                                simulated_annealing_it += 1
                                nb_iter_without_improvement = 0
                        else
                                break
                        end
                end
                if counter > maxiter
                        break
                end
        end
        # Post process : local optimization
        if postprocess
                best_solution = optimize_locally(best_solution, inst, print_function)
        end
        # Return final solution
        return best_solution, best_cost
    end

    export optimize_with_DSSP

end

max_stars_repo_name: Planheat/Planheat-Tool

max_issues_repo_name: Planheat/Planheat-Tool

max_forks_repo_name: Planheat/Planheat-Tool


text:
## THE CONSTANT REGRESSOR

const MMI = MLJModelInterface
export ConstantClassifier, ConstantRegressor,
       DeterministicConstantClassifier,
       ProbabilisticConstantClassifer

import Distributions

"""
ConstantRegressor(; distribution_type=Distributions.Normal)

A regressor that, for any new input pattern, predicts the univariate
probability distribution best fitting the training target data. Use
`predict_mean` to predict the mean value instead.
"""
struct ConstantRegressor{D} <: MMI.Probabilistic end

function ConstantRegressor(; distribution_type=Distributions.Normal)
    model   = ConstantRegressor{distribution_type}()
    message = clean!(model)
    isempty(message) || @warn message
    return model
end

function MMI.clean!(model::ConstantRegressor{D}) where D
    message = ""
    D <: Distributions.Sampleable ||
        error("$model.distribution_type is not a valid distribution_type.")
    return message
end

MMI.reformat(::ConstantRegressor, X) = (MMI.matrix(X),)
MMI.reformat(::ConstantRegressor, X, y) = (MMI.matrix(X), y)
MMI.selectrows(::ConstantRegressor, I, A) = (view(A, I, :),)
MMI.selectrows(::ConstantRegressor, I, A, y) = (view(A, I, :), y[I])

function MMI.fit(::ConstantRegressor{D}, verbosity::Int, A, y) where D
    fitresult = Distributions.fit(D, y)
    cache     = nothing
    report    = NamedTuple()
    return fitresult, cache, report
end

MMI.fitted_params(::ConstantRegressor, fitresult) =
    (target_distribution=fitresult,)

MMI.predict(::ConstantRegressor, fitresult, Xnew) =
    fill(fitresult, nrows(Xnew))

##
## THE CONSTANT DETERMINISTIC REGRESSOR (FOR TESTING)
##

struct DeterministicConstantRegressor <: MMI.Deterministic end

function MMI.fit(::DeterministicConstantRegressor, verbosity::Int, X, y)
    fitresult = mean(y)
    cache     = nothing
    report    = NamedTuple()
    return fitresult, cache, report
end

MMI.reformat(::DeterministicConstantRegressor, X) = (MMI.matrix(X),)
MMI.reformat(::DeterministicConstantRegressor, X, y) = (MMI.matrix(X), y)
MMI.selectrows(::DeterministicConstantRegressor, I, A) = (view(A, I, :),)
MMI.selectrows(::DeterministicConstantRegressor, I, A, y) =
    (view(A, I, :), y[I])

MMI.predict(::DeterministicConstantRegressor, fitresult, Xnew) =
    fill(fitresult, nrows(Xnew))

##
## THE CONSTANT CLASSIFIER
##

"""
ConstantClassifier()

A classifier that, for any new input pattern, `predict`s the
`UnivariateFinite` probability distribution `d` best fitting the
training target data. So, `pdf(d, level)` is the proportion of levels
in the training data coinciding with `level`. Use `predict_mode` to
obtain the training target mode instead.
"""
mutable struct ConstantClassifier <: MMI.Probabilistic
    testing::Bool
    bogus::Int
end

ConstantClassifier(; testing=false, bogus=0) =
    ConstantClassifier(testing, bogus)

function MMI.reformat(model::ConstantClassifier, X)
    model.testing && @info "reformatting X"
    return (MMI.matrix(X),)
end

function MMI.reformat(model::ConstantClassifier, X, y)
    model.testing && @info "reformatting X, y"
    return (MMI.matrix(X), y)
end

function MMI.reformat(model::ConstantClassifier, X, y, w)
    model.testing && @info "reformatting X, y, w"
    return (MMI.matrix(X), y, w)
end

function MMI.selectrows(model::ConstantClassifier, I, A)
    model.testing && @info "resampling X"
    return (view(A, I, :),)
end

function MMI.selectrows(model::ConstantClassifier, I, A, y)
    model.testing && @info "resampling X, y"
    return (view(A, I, :), y[I])
end

function MMI.selectrows(model::ConstantClassifier, I, A, y, ::Nothing)
    model.testing && @info "resampling X, y, nothing"
    return (view(A, I, :), y[I], nothing)
end

function MMI.selectrows(model::ConstantClassifier, I, A, y, w)
    model.testing && @info "resampling X, y, nothing"
    return (view(A, I, :), y[I], w[I])
end

# here `args` is `y` or `y, w`:
function MMI.fit(::ConstantClassifier, verbosity::Int, A, y, w=nothing)
    fitresult = Distributions.fit(MLJBase.UnivariateFinite, y, w)
    cache     = nothing
    report    = NamedTuple
    return fitresult, cache, report
end

MMI.fitted_params(::ConstantClassifier, fitresult) =
    (target_distribution=fitresult,)

MMI.predict(::ConstantClassifier, fitresult, Xnew) =
    fill(fitresult, nrows(Xnew))

##
## DETERMINISTIC CONSTANT CLASSIFIER (FOR TESTING)
##

struct DeterministicConstantClassifier <: MMI.Deterministic end

function MMI.fit(::DeterministicConstantClassifier, verbosity::Int, X, y)
    # dump missing target values and make into a regular array:
    fitresult = mode(skipmissing(y) |> collect) # a CategoricalValue
    cache     = nothing
    report    = NamedTuple()
    return fitresult, cache, report
end

MMI.reformat(::DeterministicConstantClassifier, X) = (MMI.matrix(X),)
MMI.reformat(::DeterministicConstantClassifier, X, y) = (MMI.matrix(X), y)
MMI.selectrows(::DeterministicConstantClassifier, I, A) = (view(A, I, :),)
MMI.selectrows(::DeterministicConstantClassifier, I, A, y) =
    (view(A, I, :), y[I])

MMI.predict(::DeterministicConstantClassifier, fitresult, Xnew) =
    fill(fitresult, nrows(Xnew))

#
# METADATA
#

metadata_pkg.((ConstantRegressor, ConstantClassifier,
               DeterministicConstantRegressor, DeterministicConstantClassifier),
              name="MLJModels",
              uuid="d491faf4-2d78-11e9-2867-c94bc002c0b7",
              url="https://github.com/alan-turing-institute/MLJModels.jl",
              julia=true,
              license="MIT",
              is_wrapper=false)

metadata_model(ConstantRegressor,
               input=MMI.Table,
               target=AbstractVector{MMI.Continuous},
               weights=false,
               descr="Constant regressor (Probabilistic).",
               path="MLJModels.ConstantRegressor")

metadata_model(DeterministicConstantRegressor,
               input=MMI.Table,
               target=AbstractVector{MMI.Continuous},
               weights=false,
               descr="Constant regressor (Deterministic).",
               path="MLJModels.DeterministicConstantRegressor")

metadata_model(ConstantClassifier,
               input=MMI.Table,
               target=AbstractVector{<:MMI.Finite},
               weights=true,
               descr="Constant classifier (Probabilistic).",
               path="MLJModels.ConstantClassifier")

metadata_model(DeterministicConstantClassifier,
               input=MMI.Table,
               target=AbstractVector{<:MMI.Finite},
               weights=false,
               descr="Constant classifier (Deterministic).",
               path="MLJModels.DeterministicConstantClassifier")


text:
#=

=#
macro RExpr(body)

end

function generate_expr(arg::Int64)
	return :(kRPC.Remote.KRPC.ConstantInt($(convert(Int32, arg)))
end

function generate_expr(arg::Float64)
	return :(kRPC.Remote.KRPC.ConstantDouble($arg))
end

function generate_expr(arg::Float32)
	return :(kRPC.Remote.KRPC.ConstantFloat($arg))
end

function generate_expr(arg::Bool)
	return :(kRPC.Remote.KRPC.ConstantBool($arg))
end

function generate_expr(arg::String)
	return :(kRPC.Remote.KRPC.ConstantString($arg))
end

function generate_expr()


text:
type LMMDriftCalculator
  numberOfRates::Int
  numberOfFactors::Int
  isFullFactor::Bool
  numeraire::Int
  alive::Int
  displacements::Vector{Float64}
  oneOverTaus::Vector{Float64}
  C::Matrix{Float64}
  pseudo::Matrix{Float64}
  tmp::Vector{Float64}
  e::Matrix{Float64}
  downs::Vector{Int}
  ups::Vector{Int}
end

function LMMDriftCalculator(pseudo::Matrix{Float64}, displacements::Vector{Float64}, taus::Vector{Float64}, numeraire::Int, alive::Int)
  numberOfRates = length(taus)
  prows, pcols = size(pseudo)
  numberOfFactors = pcols
  isFullFactor = numberOfFactors == numberOfRates
  tmp = zeros(numberOfRates)
  e = zeros(pcols, prows)
  downs = Vector{Int}(numberOfRates)
  ups = Vector{Int}(numberOfRates)

  # checking reqs
  numberOfRates > 0 || error("Dim out of range")
  length(displacements) == numberOfRates || error("displacements out of range")
  prows == numberOfRates || error("pseudo rows not consistent with dim")
  (pcols > 0 && pcols <= numberOfRates) || error("pseudo rows not consistent with pseudo cols")
  alive <= numberOfRates || error("Alive out of bounds")
  numeraire >= alive || error("Numeraire $(numeraire) smaller than alive $(alive)")

  # precompute 1 / taus
  oneOverTaus = 1.0 ./ taus

  # Compute covariance matrix from pseudo root
  pT = pseudo' # transpose
  C = pseudo * pT

  # Compute lower and upper extrema for (non reduced) drift calculation
  for i = alive:numberOfRates
    downs[i] = min(i+1, numeraire)
    ups[i] = max(i+1, numeraire)
  end

  return LMMDriftCalculator(numberOfRates, numberOfFactors, isFullFactor, numeraire, alive, displacements, oneOverTaus, C, pseudo,
                            tmp, e, downs, ups)
end

function compute_reduced!(lmm::LMMDriftCalculator, forwards::Vector{Float64}, drifts::Vector{Float64})
  # Compute drifts with factor reduction using pseudo square root of covariance matrix

  # precompute forwards factor
  for i = lmm.alive:lmm.numberOfRates
    lmm.tmp[i] = (forwards[i] + lmm.displacements[i]) / (lmm.oneOverTaus[i] + forwards[i])
  end

  # Enforce initialization
  for r = 1:lmm.numberOfFactors
    lmm.e[r, max(1, lmm.numeraire-1)] = 0.0
  end
  # Now compute drifts: take the numeraire P_N (numeraire = N) as the reference point
  # divide the summation into 3 steps

  # 1st Step: the drift corresponding to the numeraire P_N is zero
  # (if N=1 no drift is null, if N = numberOfRates the last drift is null)
  if lmm.numeraire > 1
    drifts[lmm.numeraire - 1] = 0.0
  end

  # 2nd Step: then, move backwards from N-2 (included) back to alive(included)
  # (if N=1, jump to 3rd Step, if N = numberOfRates, the e[r, N-1] are correctly initialized)
  for i = lmm.numeraire - 2:-1:lmm.alive
    drifts[i] = 0.0
    for r = 1:lmm.numberOfFactors
      lmm.e[r, i] = lmm.e[r, i+1] + lmm.tmp[i+1] * lmm.pseudo[i+1, r]
      drifts[i] -= lmm.e[r, i] * lmm.pseudo[i, r]
    end
  end

  # 3rd step: now move forward from N (included) up to n (excluded)
  # (if N = 0 this is the only relevant computation)
  for i = lmm.numeraire:lmm.numberOfRates
    drifts[i] = 0.0
    for r = 1:lmm.numberOfFactors
      if i == 1
        lmm.e[r, i] = lmm.tmp[i] * lmm.pseudo[i, r]
      else
        lmm.e[r, i] = lmm.e[r, i-1] + lmm.tmp[i] * lmm.pseudo[i, r]
      end
      drifts[i] += lmm.e[r, i] * lmm.pseudo[i, r]
    end
  end
  return lmm, drifts
end

function compute!(lmm::LMMDriftCalculator, forwards::Vector{Float64}, drifts::Vector{Float64})
  # TODO check for full factor
  compute_reduced!(lmm, forwards, drifts)

  return lmm, drifts
end


text:
#include("../ReadWrite/Header.jl")

"""
**SeisWindow**
*Window a seis file using header words.*
**IN**
* in
* out
* key=[]
* minval=[]
* maxval=[]
* note that windowing along the time axis is achieved by using the key "t".
**OUT**
*Credits: AS, 2015*
"""

function SeisWindow(in::String,out::String;key=[],minval=[],maxval=[])
    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
	extent = ReadTextHeader(in)
	tmin = extent.o1
	tmax = extent.o1 + extent.d1*(extent.n1-1)
	for ikey=1:length(key)
		if key[ikey] == "t"
			tmin = minval[ikey]
			tmax = maxval[ikey]
		end
	end
	itmin = convert(Int32,round((tmin - extent.o1)/extent.d1) + 1)
	if (itmin < 1)
		itmin = 1
	end
	itmax = convert(Int32,round((tmax - extent.o1)/extent.d1) + 1)
	if (itmax) > extent.n1
		itmax = extent.n1
	end
	SeisWindowHeaders(in,out;key=key,minval=minval,maxval=maxval,tmin=tmin,tmax=tmax)
	FetchTraces(in,out;itmin=itmin,itmax=itmax)
	tmp = join(["tmp_SeisWindow_",string(round(Int,rand()*100000))])
	@compat SeisProcessHeaders(out,tmp,[UpdateHeader],[Dict(:itmin=>itmin,:itmax=>itmax)])
	filename_h_tmp = join([DATAPATH tmp "@headers@"])
	filename_h_out = join([DATAPATH out "@headers@"])
	cp(filename_h_tmp,filename_h_out,remove_destination=true);
	rm(filename_h_tmp);
	#rm(tmp);
end

function FetchTraces(in::String,out::String;ntrace=500,itmin=round(Int,1),itmax=round(Int,9999999999))

	NX = GetNumTraces(out)
	itrace = round(Int,1)
	filename_data_in = ParseDataName(in)
    DATAPATH = get(ENV,"DATAPATH",join([pwd(),"/"]))
	filename_data_out = join([DATAPATH out "@data@"])
	stream_in  = open(filename_data_in)
	stream_out  = open(filename_data_out,"w")
	extent = ReadTextHeader(in)
	nt = extent.n1
	if itmin < 1
		itmin = round(Int,1)
	end
	if itmax > nt
		itmax = round(Int,nt)
	end
	while itrace <= NX
		if (itrace > 1)
			stream_out = open(filename_data_out,"a")
		end
		nx = NX - itrace + 1
		ntrace = nx > ntrace ? ntrace : nx
		h = SeisReadHeaders(out,group="some",itrace=itrace,ntrace=ntrace)
		d = zeros(Float32,itmax-itmin+1,ntrace)
		SeekTraces!(d,stream_in,h,itmin,itmax,nt,round(Int,ntrace))
		write(stream_out,d)
		close(stream_out)
		itrace += ntrace
	end
	close(stream_in)
end

function SeekTraces!{T}(d::AbstractArray{T,2},stream_in::IOStream,h::Array{Header,1},itmin,itmax,nt,ntrace)

	d1 = zeros(Float32,nt)
	for ix = 1 : ntrace
			position = 4*nt*(h[ix].tracenum-1)
			seek(stream_in,position)
			d1 = read(stream_in,Float32,nt)
			for it = itmin : itmax
				d[it - itmin + 1,ix] = d1[it]
			end
	end
	nothing
end

function UpdateHeader(h;itrace=1,itmin=1,itmax=9e9)
	ot = (itmin-1)*h[1].d1 + h[1].o1
	nt = itmax - itmin + 1
	for j = 1 : length(h)
		h[j].o1 = ot
		h[j].n1 = nt
		h[j].tracenum = j + itrace - 1
	end

	return h
end


text:
myParse(x) = parse.(Int, x)
heightmap = Matrix(hcat(myParse.(collect.(readlines("input.txt")))...)')

nRows, nCols = size(heightmap)

function surrounding(row, col)
    arr = [(row - 1, col), (row, col - 1), (row + 1, col), (row, col + 1)]
    return [heightmap[x, y] for (x, y) in arr if x in 1:nRows && y in 1:nCols]
end

function coordfromIndex(index)
    row = ((index - 1) % nRows) + 1
    col = ((index - 1) ÷ nRows) + 1
    return row, col
end

function indexfromCoord(row, col)
    return nRows * (col - 1) + row
end

function surrounding(index)
    row, col = coordfromIndex(index)
    return surrounding(row, col)
end

function surroundingIndices(index)
    row, col = coordfromIndex(index)
    arr = [(row - 1, col), (row, col - 1), (row + 1, col), (row, col + 1)]
    return [indexfromCoord(x, y) for (x, y) in arr if x in 1:nRows && y in 1:nCols]
end

#! Part 1

function checkLower(height, index)
    val = min(surrounding(index)...)
    return height < val
end

function sumP1()
    sum = 0

    for (i, height) in enumerate(heightmap)
        if checkLower(height, i)
            sum += height + 1
        end
    end

    return sum
end

sumP1()

#! Part 2

function getProductPt2()
    visited = [heightmap[i] == 9 ? true : false for i in 1:length(heightmap)]

    basinSizes = []

    function basinSearch(index, basin = [])
        if !visited[index]
            visited[index] = true
            push!(basin, index)
            for i in surroundingIndices(index)
                basinSearch(i, basin)
            end
        end
        return length(basin)
    end

    for (i, _) in enumerate(heightmap)
        size = basinSearch(i)
        if size != 0
            push!(basinSizes, size)
        end
    end

    sort!(basinSizes, rev=true)

    return prod(basinSizes[1:3])

end

getProductPt2()

max_stars_repo_path: Day 9/main.jl

max_issues_repo_path: Day 9/main.jl

max_forks_repo_path: Day 9/main.jl


text:
module ALASMPCCmod

using OutputALASmod
using ActifMPCCmod
using UnconstrainedMPCCActif
using MPCCmod
using Relaxation

using NLPModels

"""
Type ALASMPCC :

liste des constructeurs :
ALASMPCC(mod::MPCCmod.MPCC,r::Float64,s::Float64,t::Float64, prec::Float64)
ALASMPCC(mod::MPCCmod.MPCC,r::Float64,s::Float64,t::Float64, prec::Float64, rho::Vector)

liste des méthodes :

liste des accesseurs :

liste des fonctions :
solvePAS(alas::ALASMPCC)

EndingTest(alas::ALASMPCC,Armijosuccess::Bool,small_step::Bool,feas::Float64,dual_feas::Float64,k::Int64)
SlackComplementarityProjection(alas::ALASMPCC)
LagrangeInit(alas::ALASMPCC,rho::Vector,xj::Vector)
LagrangeCompInit(alas::ALASMPCC,rho::Vector,xj::Vector)
LagrangeUpdate(alas::ALASMPCC,rho::Vector,xjk::Vector,uxl::Vector,uxu::Vector,ucl::Vector,ucu::Vector,usg::Vector,ush::Vector)
Penaltygen(alas::ALASMPCCmod.ALASMPCC,x::Vector,yg::Vector,yh::Vector,rho::Vector,usg::Vector,ush::Vector,uxl::Vector,uxu::Vector,ucl::Vector,ucu::Vector)
RhoUpdate(rho::Vector,err::Vector,alas::ALASMPCC)
RhoDetail(rho::Vector,alas::ALASMPCC)
"""

type ALASMPCC
 mod::MPCCmod.MPCC
 #paramètres du pb
 r::Float64
 s::Float64
 t::Float64
 tb::Float64

 #paramètres algorithmiques
 prec::Float64 #precision à 0
 rho_init::Vector #nombre >= 0
end

function ALASMPCC(mod::MPCCmod.MPCC,r::Float64,s::Float64,t::Float64, prec::Float64, rho::Vector)

 rho_init=rho
 tb=mod.paramset.tb(r,s,t)

 return ALASMPCC(mod,r,s,t,tb,prec,rho_init)
end

"""
solvePAS() : méthode de penalisation/activation de contraintes DHKM 16'

Méthode avec gestion à l'intérieur du paramètre de pénalité
"""
function solvePAS(alas::ALASMPCC; verbose::Bool=true)

#Initialisation paramètres :
 rho=alas.rho_init
 k_max=alas.mod.paramset.ite_max_alas
 l_max=alas.mod.paramset.ite_max_viol
 rho_update=alas.mod.paramset.rho_update
 obj_viol=alas.mod.paramset.goal_viol

 n=length(alas.mod.mp.meta.x0)

 gradpen=Vector(n)
 gradpen_prec=Vector(n)
 xj=Vector(n+2*alas.mod.nb_comp)
 xjk=Vector(n+2*alas.mod.nb_comp)
 xjkl=Vector(n+2*alas.mod.nb_comp)

# S0 : initialisation du problème avec slack (projection sur _|_ )
 xj=SlackComplementarityProjection(alas)

 #variables globales en sortie du LineSearch
 step=1.0
 Armijosuccess=true
 wnew=zeros(Bool,0,0)

 #Initialisation multiplicateurs de la contrainte de complémentarité
 lg,lh,lphi=LagrangeCompInit(alas,rho,xj)

 # Initialisation multiplicateurs : uxl,uxu,ucl,ucu :
 uxl,uxu,ucl,ucu,usg,ush=LagrangeInit(alas,rho,xj)

# S1 : Initialisation du MPCC_Actif
 #objectif du problème avec pénalité lagrangienne sans contrainte :
 pen_mpcc=CreatePenaltyNLP(alas,xj,rho,usg,ush,uxl,uxu,ucl,ucu)

 #initialise notre problème pénalisé avec des contraintes actives w
 ma=ActifMPCCmod.MPCC_actif(pen_mpcc,alas.r,alas.s,alas.t,
                            alas.mod.nb_comp,alas.mod.paramset,
                            alas.mod.algoset.direction,
                            alas.mod.algoset.linesearch)

# S2 : Major Loop Activation de contrainte
 #initialisation
 k=0
 xjk=xj
 small_step=false
 ht,gradpen=NLPModels.objgrad(ma.nlp,xjk)

 #direction initial
 dj=zeros(n+2*alas.mod.nb_comp)

 l_negative=findfirst(x->x<0,[lg;lh;lphi])!=0

 feas=MPCCmod.viol_contrainte_norm(alas.mod,xjk)
 feasible=feas<=alas.prec

 dual_feas=norm(ActifMPCCmod.grad(ma,xjk,gradpen),Inf)
 multi_norm=alas.mod.algoset.scaling_dual(usg,ush,uxl,uxu,ucl,ucu,
                                          lg,lh,lphi,alas.prec,
                                          alas.mod.paramset.precmpcc,
                                          norm(rho,Inf),dual_feas)
 dual_feasible=dual_feas/multi_norm<=alas.prec

 oa=OutputALASmod.OutputALAS(xjk,dj,feas,dual_feas,rho,ht)

 while k<k_max && (l_negative || !((feasible || minimum(rho)==alas.mod.paramset.rho_max*ma.crho ) && dual_feasible)) && Armijosuccess && !small_step

  l=0
  xjkl=xjk
  feask=feas
  gradpen_prec=gradpen

  #boucle pour augmenter le paramètre de pénalisation tant qu'on a pas baissé suffisament la violation
  while l<l_max && !dual_feasible && Armijosuccess && !small_step
   #Unconstrained Solver modifié pour résoudre le sous-problème avec contraintes de bornes
 xjkl,ma.w,dj,step,wnew,outputArmijo,small_step,ols,gradpen,ht=UnconstrainedMPCCActif.LineSearchSolve(ma,xjkl,dj,step,gradpen,ht)
   feas=MPCCmod.viol_contrainte_norm(alas.mod,xjkl);
   feasible=feas<=alas.prec

   dual_feas=norm(ActifMPCCmod.grad(ma,xjkl,gradpen),Inf)
   dual_feasible=dual_feas<=alas.mod.algoset.unconstrained_stopping(alas.prec,rho)

   OutputALASmod.Update(oa,xjkl,rho,feas,dual_feas,dj,step,ols,ht,n)

   Armijosuccess=(outputArmijo==0)
   if outputArmijo==2
    @show "Unbounded subproblem"
    return xj,EndingTest(alas,Armijosuccess,small_step,feas,dual_feas,k),rho,oa
   end
   l+=1
  end

  #On met à jour rho si ça n'a pas été:
  if (l==l_max || small_step || !Armijosuccess || dual_feasible) && (feas>obj_viol*feask && !feasible)
   verbose && print_with_color(:red, "Max ité Unc. Min. l=$l |x|=$(norm(xjkl,Inf)) |c(x)|=$feas |L'|=$dual_feas Arm=$Armijosuccess small_step=$small_step rho=$(norm(rho,Inf))  \n")
   rho=RhoUpdate(alas,rho,ma.crho,abs(MPCCmod.viol_contrainte(alas.mod,xjkl)))
   ActifMPCCmod.setcrho(ma,alas.mod.algoset.crho_update(feas,rho)) #on change le problème donc on réinitialise beta

   #met à jour la fonction objectif après la mise à jour de rho
   ma.nlp,ht,gradpen=UpdatePenaltyNLP(alas,rho,xjkl,usg,ush,uxl,uxu,ucl,ucu,ma.nlp,objpen=ht,gradpen=gradpen,crho=ma.crho)
   ActifMPCCmod.setbeta(ma,0.0) #on change le problème donc on réinitialise beta
   #ActifMPCCmod.sethess(ma,H) #on change le problème donc on réinitialise Hess
   dual_feas=norm(ActifMPCCmod.grad(ma,xjkl,gradpen),Inf)
   dual_feasible=dual_feas<=alas.prec
  end
  xjk=xjkl

  #Mise à jour des paramètres de la pénalité Lagrangienne (uxl,uxu,ucl,ucu)
  uxl,uxu,ucl,ucu,usg,ush=LagrangeUpdate(alas,rho,xjk,uxl,uxu,ucl,ucu,usg,ush)

  #met à jour la fonction objectif après la mise à jour des multiplicateurs
  temp,ht,gradpen=UpdatePenaltyNLP(alas,rho,xjk,usg,ush,uxl,uxu,ucl,ucu,
                                   ma.nlp,objpen=ht,gradpen=gradpen)

  #Mise à jour des multiplicateurs de la complémentarité
  if alas.mod.nb_comp>0
   lg,lh,lphi=ActifMPCCmod.LSQComputationMultiplier(ma,gradpen,xjk)
   l_negative=findfirst(x->x<0,[lg;lh;lphi])!=0 #vrai si un multiplicateur est negatif
  end
  multi_norm=alas.mod.algoset.scaling_dual(usg,ush,uxl,uxu,ucl,ucu,lg,lh,lphi,alas.prec,alas.mod.paramset.precmpcc,norm(rho,Inf),dual_feas)

  feasible=feas<=alas.prec
  dual_feasible=dual_feas/multi_norm<=alas.prec

  #si on bloque mais qu'un multiplicateur est <0 on continue
  if l_negative && (!Armijosuccess || small_step)
   Armijosuccess=true
   small_step=false
  end

  #Relaxation rule si on a fait un pas de Armijo-Wolfe et si un multiplicateur est du mauvais signe
  if k!=0 && (dot(gradpen,dj)>=alas.mod.paramset.tau_wolfe*dot(gradpen_prec,dj) || step==0.0) && findfirst(x->x<0,[lg;lh;lphi])!=0

   ActifMPCCmod.RelaxationRule(ma,xjk,lg,lh,lphi,wnew)
   if findfirst(x->x<0,[lg;lh;lphi])!=0
    dual_feasible=false #si on a relaché des contraintes, on force à faire une itération
   end
   verbose && print_with_color(:yellow, "Active set: $(ma.wcc) \n")
  end

  k+=1
  verbose && k>=k_max && print_with_color(:red, "Max ité. Lagrangien \n")
 end

#Traitement finale :
 xj=xjk

 stat=EndingTest(alas,Armijosuccess,small_step,feas,dual_feas,k)

 return xj,stat,rho,oa
end

function EndingTest(alas::ALASMPCC,Armijosuccess::Bool,small_step::Bool,feas::Float64,dual_feas::Float64,k::Int64)

 stat=0
 if !Armijosuccess
  print_with_color(:red, "Failure : Armijo Failure\n")
  stat=1
 end
 if small_step
  print_with_color(:red, "Failure : too small step\n")
  #stat=1
 end
 if feas>alas.prec
  print_with_color(:red, "Failure : Infeasible Solution. norm: $feas\n")
  #stat=1
 end
 if dual_feas>alas.prec
  if k>=alas.mod.paramset.ite_max_alas
   print_with_color(:red, "Failure : Non-optimal Sol. norm: $dual_feas\n")
   stat=1
  else
   print_with_color(:red, "Inexact : Fritz-John Sol. norm: $dual_feas\n")
   #stat=0
  end
 end

 return stat
end

function Feasible(alas,x::Vector)
 return MPCCmod.viol_contrainte_norm(alas.mod,x)<=alas.prec
end
function DualFeasible(alas,gradpen::Vector)
 return norm(gradpen[1:n],Inf)<=alas.prec
end

"""
SlackComplementarityProjection :
projete le point x0 sur la contrainte de complémentarité avec les slack
pré-requis : x0 doit être de taille n+2q
"""
function SlackComplementarityProjection(alas::ALASMPCC)

 nb_comp=alas.mod.nb_comp

 if nb_comp==0
  return alas.mod.xj
 end

 #initialisation :
 #x=[alas.mod.xj;alas.mod.G(alas.mod.xj);alas.mod.H(alas.mod.xj)]
 x=[alas.mod.xj;NLPModels.cons(alas.mod.G,alas.mod.xj);NLPModels.cons(alas.mod.H,alas.mod.xj)]
 n=length(x)-2*nb_comp

 #projection sur les contraintes de positivité relaxés : yG>=tb et yH>=tb
 x[n+1:n+nb_comp]=max.(x[n+1:n+nb_comp],ones(nb_comp)*alas.tb)
 x[n+nb_comp+1:n+2*nb_comp]=max.(x[n+nb_comp+1:n+2*nb_comp],ones(nb_comp)*alas.tb)

 #projection sur les contraintes papillons : yG<=psi(yH,r,s,t) et yH<=psi(yG,r,s,t)
 for i=1:nb_comp
  psiyg=Relaxation.psi(x[n+i],alas.r,alas.s,alas.t)
  psiyh=Relaxation.psi(x[n+nb_comp+i],alas.r,alas.s,alas.t)

  if x[n+i]-psiyh>0 && x[n+nb_comp+i]-psiyg>0
   x[n+i]>=x[n+nb_comp+i] ? x[n+nb_comp+i]=psiyg : x[n+i]=psiyh
  end
 end

 return x
end

"""
Initialisation des multiplicateurs de Lagrange
"""
function LagrangeInit(alas::ALASMPCC,rho::Vector,xj::Vector)

  n=length(alas.mod.mp.meta.x0)
  rho_eqg,rho_eqh,rho_ineq_lvar,rho_ineq_uvar,rho_ineq_lcons,rho_ineq_ucons=RhoDetail(alas,rho)

  uxl=max.(rho_ineq_lvar.*(alas.mod.mp.meta.lvar-xj[1:n]),zeros(n))
  uxu=max.(rho_ineq_uvar.*(xj[1:n]-alas.mod.mp.meta.uvar),zeros(n))
  if alas.mod.mp.meta.ncon!=0
   c(z)=NLPModels.cons(alas.mod.mp,z)
   nc=length(alas.mod.mp.meta.y0) #nombre de contraintes

   ucl=max.(rho_ineq_lcons.*(alas.mod.mp.meta.lcon-c(xj[1:n])),zeros(nc))
   ucu=max.(rho_ineq_ucons.*(c(xj[1:n])-alas.mod.mp.meta.ucon),zeros(nc))
  else
   ucl=[];ucu=[];
  end

  usg=zeros(alas.mod.nb_comp)
  ush=zeros(alas.mod.nb_comp)

 return uxl,uxu,ucl,ucu,usg,ush
end

function LagrangeCompInit(alas::ALASMPCC,rho::Vector,xj::Vector)

 n=length(alas.mod.mp.meta.x0)
 nb_comp=alas.mod.nb_comp
 rho_eqg,rho_eqh,rho_ineq_lvar,rho_ineq_uvar,rho_ineq_lcons,rho_ineq_ucons=RhoDetail(alas,rho)

 psiyg=Relaxation.psi(xj[n+1:n+nb_comp],alas.r,alas.s,alas.t)
 psiyh=Relaxation.psi(xj[n+1+nb_comp:n+2*nb_comp],alas.r,alas.s,alas.t)
 phi=(xj[n+1:n+nb_comp]-psiyh).*(xj[n+1+nb_comp:n+2*nb_comp]-psiyg)

 lphi=max.(phi,zeros(nb_comp))
 lg=max.(-rho_eqg.*(xj[n+1:n+nb_comp]-alas.tb),zeros(nb_comp))
 lh=max.(-rho_eqh.*(xj[n+1+nb_comp:n+2*nb_comp]-alas.tb),zeros(nb_comp))

 return lg,lh,lphi
end

"""
Mise à jour des multiplicateurs de Lagrange
"""
function LagrangeUpdate(alas::ALASMPCC,rho::Vector,xjk::Vector,
                        uxl::Vector,uxu::Vector,
                        ucl::Vector,ucu::Vector,
                        usg::Vector,ush::Vector)

  n=length(alas.mod.mp.meta.x0)
   rho_eqg,rho_eqh,rho_ineq_lvar,rho_ineq_uvar,rho_ineq_lcons,rho_ineq_ucons=RhoDetail(alas,rho)

  uxl=uxl+max.(rho_ineq_lvar.*(xjk[1:n]-alas.mod.mp.meta.uvar),-uxl)
  uxu=uxu+max.(rho_ineq_uvar.*(alas.mod.mp.meta.lvar-xjk[1:n]),-uxu)
  if alas.mod.mp.meta.ncon!=0
   c=NLPModels.cons(alas.mod.mp,xjk[1:n])
   ucl=ucl+max(rho_ineq_lcons.*(c-alas.mod.mp.meta.ucon),-ucl)
   ucu=ucu+max(rho_ineq_ucons.*(alas.mod.mp.meta.lcon-c),-ucu)
  end

  if alas.mod.nb_comp!=0
   G(x)=NLPModels.cons(alas.mod.G,x)
   H(x)=NLPModels.cons(alas.mod.H,x)
   usg=usg+rho_eqg.*(G(xjk[1:n])-xjk[n+1:n+alas.mod.nb_comp])
   ush=ush+rho_eqh.*(H(xjk[1:n])-xjk[n+alas.mod.nb_comp+1:n+2*alas.mod.nb_comp])
  end

 return uxl,uxu,ucl,ucu,usg,ush
end

"""
Fonction de pénalité générique :
"""
function Penaltygen(alas::ALASMPCCmod.ALASMPCC,
                    x::Vector,yg::Vector,yh::Vector,
                    rho::Vector,usg::Vector,ush::Vector,
                    uxl::Vector,uxu::Vector,
                    ucl::Vector,ucu::Vector)
 nb_comp=alas.mod.nb_comp

 err=MPCCmod.viol_contrainte(alas.mod,x,yg,yh)

 return alas.mod.algoset.penalty(alas.mod.mp,err,alas.mod.nb_comp,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
end

function GradPenaltygen(alas::ALASMPCCmod.ALASMPCC,
                        x::Vector,yg::Vector,yh::Vector,
                        rho::Vector,usg::Vector,ush::Vector,
                        uxl::Vector,uxu::Vector,
                        ucl::Vector,ucu::Vector)
 g=Vector(x)
 err=MPCCmod.viol_contrainte(alas.mod,x,yg,yh)

 return alas.mod.algoset.penalty(alas.mod.mp,alas.mod.G,alas.mod.H,
                                 err,alas.mod.nb_comp,
                                 x,yg,yh,rho,
                                 usg,ush,uxl,uxu,ucl,ucu,g)
end

function HessPenaltygen(alas::ALASMPCCmod.ALASMPCC,
                        x::Vector,yg::Vector,yh::Vector,
                        rho::Vector,usg::Vector,ush::Vector,
                        uxl::Vector,uxu::Vector,
                        ucl::Vector,ucu::Vector)
 n=length(alas.mod.mp.meta.x0)
 nb_comp=alas.mod.nb_comp

 Hess=zeros(0,0)
 Hf=tril([NLPModels.hess(alas.mod.mp,x) zeros(n,2*nb_comp);zeros(2*nb_comp,n) eye(2*nb_comp)])

 err=MPCCmod.viol_contrainte(alas.mod,x,yg,yh)

 return Hf+alas.mod.algoset.penalty(alas.mod.mp,alas.mod.G,alas.mod.H,
                                 err,nb_comp,
                                 x,yg,yh,rho,
                                 usg,ush,uxl,uxu,ucl,ucu,Hess)
end

function ObjGradPenaltygen(alas::ALASMPCCmod.ALASMPCC,
                            x::Vector,yg::Vector,yh::Vector,
                            rho::Vector,usg::Vector,ush::Vector,
                            uxl::Vector,uxu::Vector,
                            ucl::Vector,ucu::Vector)
 n=length(x)
 nb_comp=alas.mod.nb_comp

 err=MPCCmod.viol_contrainte(alas.mod,x,yg,yh)

 f=NLPModels.obj(alas.mod.mp,x)
 f+=alas.mod.algoset.penalty(alas.mod.mp,err,alas.mod.nb_comp,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
 g=vec([NLPModels.grad(alas.mod.mp,x)' zeros(2*alas.mod.nb_comp)'])
 g+=alas.mod.algoset.penalty(alas.mod.mp,alas.mod.G,alas.mod.H,
                                 err,alas.mod.nb_comp,
                                 x,yg,yh,rho,
                                 usg,ush,uxl,uxu,ucl,ucu,Array{Float64}(0))

 return f,g
end

function ObjGradHessPenaltygen(alas::ALASMPCCmod.ALASMPCC,
                            x::Vector,yg::Vector,yh::Vector,
                            rho::Vector,usg::Vector,ush::Vector,
                            uxl::Vector,uxu::Vector,
                            ucl::Vector,ucu::Vector)
 n=length(x)
 nb_comp=alas.mod.nb_comp

 err=MPCCmod.viol_contrainte(alas.mod,x,yg,yh)

 f=NLPModels.obj(alas.mod.mp,x)
 f+=alas.mod.algoset.penalty(alas.mod.mp,err,alas.mod.nb_comp,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
 g=vec([NLPModels.grad(alas.mod.mp,x)' zeros(2*alas.mod.nb_comp)'])
 g+=alas.mod.algoset.penalty(alas.mod.mp,alas.mod.G,alas.mod.H,
                                 err,alas.mod.nb_comp,
                                 x,yg,yh,rho,
                                 usg,ush,uxl,uxu,ucl,ucu,g)
 Hess=tril([NLPModels.hess(alas.mod.mp,x) zeros(n,n);zeros(2*nb_comp,n) eye(2*nb_comp)])
 Hess+=alas.mod.algoset.penalty(alas.mod.mp,alas.mod.G,alas.mod.H,
                                 err,nb_comp,
                                 x,yg,yh,rho,
                                 usg,ush,uxl,uxu,ucl,ucu,Array{Float64}(0,0))

 return f,g,Hess
end

function CreatePenaltyNLP(alas::ALASMPCCmod.ALASMPCC,
                          xj::Vector,rho::Vector,
                          usg::Vector,ush::Vector,
                          uxl::Vector,uxu::Vector,
                          ucl::Vector,ucu::Vector)
 n=length(alas.mod.mp.meta.x0)
 nb_comp=alas.mod.nb_comp

 penf(x,yg,yh)=NLPModels.obj(alas.mod.mp,x)+Penaltygen(alas,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
 penf(x)=penf(x[1:n],x[n+1:n+alas.mod.nb_comp],x[n+alas.mod.nb_comp+1:n+2*alas.mod.nb_comp])

 gpenf(x,yg,yh)=vec([NLPModels.grad(alas.mod.mp,x)' zeros(2*alas.mod.nb_comp)'])+GradPenaltygen(alas,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
 gpenf(x)=gpenf(x[1:n],x[n+1:n+alas.mod.nb_comp],x[n+alas.mod.nb_comp+1:n+2*alas.mod.nb_comp])

 gfpenf(x,yg,yh)=ObjGradPenaltygen(alas,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
 gfpenf(x)=gfpenf(x[1:n],x[n+1:n+alas.mod.nb_comp],x[n+alas.mod.nb_comp+1:n+2*alas.mod.nb_comp])

 Hpenf(x,yg,yh)=HessPenaltygen(alas,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
 Hpenf(x)=Hpenf(x[1:n],x[n+1:n+alas.mod.nb_comp],x[n+alas.mod.nb_comp+1:n+2*alas.mod.nb_comp])

 return NLPModels.SimpleNLPModel(x->penf(x), xj,
                                lvar=[alas.mod.mp.meta.lvar;alas.tb*ones(2*alas.mod.nb_comp)],
                                uvar=[alas.mod.mp.meta.uvar;Inf*ones(2*alas.mod.nb_comp)],
				g=x->gpenf(x),H=x->Hpenf(x),fg=x->gfpenf(x))
end

function UpdatePenaltyNLP(alas::ALASMPCCmod.ALASMPCC,
                          rho::Vector,xj::Vector,
                          usg::Vector,ush::Vector,
                          uxl::Vector,uxu::Vector,
                          ucl::Vector,ucu::Vector,
                          pen_nlp::NLPModels.AbstractNLPModel;gradpen::Vector=[],objpen::Float64=zeros(0),crho::Float64=1.0)
 n=length(alas.mod.mp.meta.x0)

 penf(x,yg,yh)=NLPModels.obj(alas.mod.mp,x)+Penaltygen(alas,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
 penf(x)=penf(x[1:n],x[n+1:n+alas.mod.nb_comp],x[n+alas.mod.nb_comp+1:n+2*alas.mod.nb_comp])

 gpenf(x,yg,yh)=vec([NLPModels.grad(alas.mod.mp,x)' zeros(2*alas.mod.nb_comp)'])+GradPenaltygen(alas,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
 gpenf(x)=gpenf(x[1:n],x[n+1:n+alas.mod.nb_comp],x[n+alas.mod.nb_comp+1:n+2*alas.mod.nb_comp])

 gfpenf(x,yg,yh)=ObjGradPenaltygen(alas,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
 gfpenf(x)=gfpenf(x[1:n],x[n+1:n+alas.mod.nb_comp],x[n+alas.mod.nb_comp+1:n+2*alas.mod.nb_comp])

 Hpenf(x,yg,yh)=HessPenaltygen(alas,x,yg,yh,rho,usg,ush,uxl,uxu,ucl,ucu)
 Hpenf(x)=Hpenf(x[1:n],x[n+1:n+alas.mod.nb_comp],x[n+alas.mod.nb_comp+1:n+2*alas.mod.nb_comp])

 pen_nlp.f=x->penf(x)
 pen_nlp.g=x->gpenf(x)
 pen_nlp.fg=x->gfpenf(x)
 pen_nlp.H=x->Hpenf(x)

 if isempty(objpen)
  return pen_nlp
 else
   if minimum(rho)==maximum(rho)
    #if minimum(rho)==alas.mod.paramset.rhomax
    #fobj=NLPModels.obj(alas.mod.mp,xj)
    #f=fobj+alas.mod.paramset.rho_update*(objpen-fobj)
    #g=alas.mod.paramset.rho_max*crho*(gradpen-vec([NLPModels.grad(alas.mod.mp,xj)' zeros(2*alas.mod.nb_comp)']))
    f,g=gfpenf(xj)
   else
    f,g=gfpenf(xj)
   end

  return pen_nlp,f,g
 end
end

"""
Mise à jour de rho:
structure de rho :
2*mod.nb_comp
length(mod.mp.meta.lvar)
length(mod.mp.meta.uvar)
length(mod.mp.meta.lcon)
length(mod.mp.meta.ucon)
"""
function RhoUpdate(alas::ALASMPCC,rho::Vector,crho::Float64,err::Vector;l::Int64=0)
 #rho[find(x->x>0,err)]*=alas.rho_update
 #rho*=alas.mod.paramset.rho_update
 rho=min.(rho*alas.mod.paramset.rho_update,ones(length(rho))*alas.mod.paramset.rho_max*crho)

 return rho
end

"""
Renvoie : rho_eq, rho_ineq_var, rho_ineq_cons
"""
function RhoDetail(alas::ALASMPCC,rho::Vector)
 nb_ineq_lvar=length(alas.mod.mp.meta.lvar)
 nb_ineq_uvar=length(alas.mod.mp.meta.uvar)
 nb_ineq_lcons=length(alas.mod.mp.meta.lcon)
 nb_ineq_ucons=length(alas.mod.mp.meta.ucon)
 nb_comp=alas.mod.nb_comp

 return rho[1:nb_comp],rho[nb_comp+1:2*nb_comp],rho[2*nb_comp+1:2*nb_comp+nb_ineq_lvar],
        rho[2*nb_comp+nb_ineq_lvar+1:2*nb_comp+nb_ineq_lvar+nb_ineq_uvar],
        rho[2*nb_comp+nb_ineq_lvar+nb_ineq_uvar+1:2*nb_comp+nb_ineq_lvar+nb_ineq_uvar+nb_ineq_lcons],
        rho[2*nb_comp+nb_ineq_lvar+nb_ineq_uvar+nb_ineq_lcons+1:2*nb_comp+nb_ineq_lvar+nb_ineq_uvar+nb_ineq_lcons+nb_ineq_ucons]
end

#end of module
end


text:
mutable struct API
    source
    structs::OrderedDict
    funcs::OrderedDict
    consts::OrderedDict
    enums::OrderedDict
    eval
end

API(source, structs, funcs, consts, enums) = API(source, structs, funcs, consts, enums, nothing)

extract_fields(apis, sym) = OrderedDict(vcat(collect.(getproperty.(apis, sym))...))
merge(apis::API...) = API(getproperty.(apis, :source), extract_fields.(Ref(apis), [:structs, :funcs, :consts, :enums])..., nothing)

Base.show(io::IO, pf::API) = print(io, "API with $(length(pf.structs)) structs, $(length(pf.funcs)) functions, $(length(pf.consts)) consts and $(length(pf.enums)) enums from $(pf.source)")

function definition_begins(line, ::Type{FDefinition})
    !isnothing(match(r"\w+\(.*\)\s+=(?!=)", line)) || !isnothing(match(r"^function ", line))
end
definition_ends(line, ::Type{FDefinition}) = strip(line) == "end" || (definition_begins(line, FDefinition) && isnothing(match(r"function \w+", line)))

definition_begins(line, ::Type{SDefinition}) = split(line, " ")[1] ∈ ["mutable", "struct"]
definition_ends(line, ::Type{SDefinition}) = (definition_begins(line, SDefinition) && split(line, " ")[end] == "end") || strip(line) == "end"

definition_begins(line, ::Type{CDefinition}) = startswith(line, "const ")
definition_ends(line, ::Type{CDefinition}) = !endswith(strip(line), ",")

definition_begins(line, ::Type{EDefinition}) = strip(first(split(line, " "))) ∈ ["@enum", "@cenum"]
definition_ends(line, ::Type{EDefinition}) = (definition_begins(line, EDefinition) && strip(split(line, " ")[end]) != "begin") || "end" ∈ strip.(split(line, " "))

function detect_definition(io, T)
    lines = String[]
    line = readline(io)
    if eof(io)
        definition_begins(line, T) && definition_ends(line, T) && return [line]
        return
    end
    while !definition_begins(line, T)
        line = readline(io)
        eof(io) && definition_begins(line, T) && definition_ends(line, T) && return [line]
        eof(io) && return
    end
    push!(lines, line)
    while !definition_ends(line, T)
        line = readline(io)
        eof(io) && !definition_ends(line, T) && error("EOF while scanning for end of definition for type $T")
        push!(lines, line)
    end
    lines
end

function parse_for_definition(file, ::Type{T}) where {T <: Declaration}
    defs = T[]
    open(file) do io
        while !eof(io)
            lines = detect_definition(io, T)
            if !isnothing(lines)
                decl = T(join(lines, "\n"))
                !(decl isa FDefinition && last(decl.signature.args).name == "fun_ptr") && !isalias(decl.name) ? push!(defs, decl) : nothing
            end
        end
    end
    defs
end

function convert_constptr_to_struct!(api)
    for cdef ∈ values(api.consts)
        if is_opaque(api.eval(cdef.value)) && !startswith(cdef.name, "PFN") && !isalias(cdef.name) # ignore function pointers
            sdef = parse_ptr(cdef.name)
            api.structs[sdef.name] = sdef
        end
    end
end

function parse_ptr(sym)
    base = remove_prefix(sym)
    SDefinition(sym, false, OrderedDict("handle" => "Ptr{Nothing}"))
end

function API(file, eval_symbol)
    defs = []
    for decl ∈ [SDefinition, FDefinition, CDefinition, EDefinition]
        push!(defs, parse_for_definition(file, decl))
    end
    api = API(file, OrderedDict.(defs)...)
    api.eval = eval_symbol
    convert_constptr_to_struct!(api)
    api
end

API(files::AbstractArray{<: AbstractString}, eval_symbol) = merge(API.(files, eval_symbol)...)


text:
## Abstract interface for planners and solutions ##
export Planner

"Abstract planner type, which defines the interface for planners."
abstract type Planner end

(planner::Planner)(domain::Domain, state::State, spec) =
    solve(planner, domain, state, spec)

solve(planner::Planner, domain::Domain, state::State, spec::Specification) =
    error("Not implemented.")
solve(planner::Planner, domain::Domain, state::State, goals::Vector{<:Term}) =
    solve(planner, domain, state, Specification(goals))
solve(planner::Planner, domain::Domain, state::State, goal::Term) =
    solve(planner, domain, state, Specification(goal))

# Path-search planning algorithms
include("planners/path_search.jl")
include("planners/bfs.jl")
include("planners/forward.jl")
include("planners/backward.jl")

# Policy-based planning algorithms
include("planners/rtdp.jl")
include("planners/mcts.jl")

# External planners
include("planners/external.jl")


text:
@doc replace(ParallelKernel.HIDE_COMMUNICATION_DOC, "@init_parallel_kernel" => "@init_parallel_stencil") macro hide_communication(args...) check_initialized(); esc(:(ParallelStencil.ParallelKernel.@hide_communication($(args...)))); end
@doc replace(ParallelKernel.ZEROS_DOC,              "@init_parallel_kernel" => "@init_parallel_stencil") macro zeros(args...)              check_initialized(); esc(:(ParallelStencil.ParallelKernel.@zeros($(args...)))); end
@doc replace(ParallelKernel.ONES_DOC,               "@init_parallel_kernel" => "@init_parallel_stencil") macro ones(args...)               check_initialized(); esc(:(ParallelStencil.ParallelKernel.@ones($(args...)))); end
@doc replace(ParallelKernel.RAND_DOC,               "@init_parallel_kernel" => "@init_parallel_stencil") macro rand(args...)               check_initialized(); esc(:(ParallelStencil.ParallelKernel.@rand($(args...)))); end
@doc replace(ParallelKernel.PARALLEL_INDICES_DOC,   "@init_parallel_kernel" => "@init_parallel_stencil") macro parallel_indices(args...)   check_initialized(); esc(:(ParallelStencil.ParallelKernel.@parallel_indices($(args...)))); end
@doc replace(ParallelKernel.PARALLEL_ASYNC_DOC,     "@init_parallel_kernel" => "@init_parallel_stencil") macro parallel_async(args...)     check_initialized(); esc(:(ParallelStencil.ParallelKernel.@parallel_async($(args...)))); end
@doc replace(ParallelKernel.SYNCHRONIZE_DOC,        "@init_parallel_kernel" => "@init_parallel_stencil") macro synchronize(args...)        check_initialized(); esc(:(ParallelStencil.ParallelKernel.@synchronize($(args...)))); end
# NOTE: @parallel does not appear here as it is extended and therefore defined in parallel.jl
@doc replace(ParallelKernel.GRIDDIM_DOC,            "@init_parallel_kernel" => "@init_parallel_stencil") macro gridDim(args...)            check_initialized(); esc(:(ParallelStencil.ParallelKernel.@gridDim($(args...)))); end
@doc replace(ParallelKernel.BLOCKIDX_DOC,           "@init_parallel_kernel" => "@init_parallel_stencil") macro blockIdx(args...)           check_initialized(); esc(:(ParallelStencil.ParallelKernel.@blockIdx($(args...)))); end
@doc replace(ParallelKernel.BLOCKDIM_DOC,           "@init_parallel_kernel" => "@init_parallel_stencil") macro blockDim(args...)           check_initialized(); esc(:(ParallelStencil.ParallelKernel.@blockDim($(args...)))); end
@doc replace(ParallelKernel.THREADIDX_DOC,          "@init_parallel_kernel" => "@init_parallel_stencil") macro threadIdx(args...)          check_initialized(); esc(:(ParallelStencil.ParallelKernel.@threadIdx($(args...)))); end
@doc replace(ParallelKernel.SYNCTHREADS_DOC,        "@init_parallel_kernel" => "@init_parallel_stencil") macro sync_threads(args...)       check_initialized(); esc(:(ParallelStencil.ParallelKernel.@sync_threads($(args...)))); end
@doc replace(ParallelKernel.SHAREDMEM_DOC,          "@init_parallel_kernel" => "@init_parallel_stencil") macro sharedMem(args...)          check_initialized(); esc(:(ParallelStencil.ParallelKernel.@sharedMem($(args...)))); end
@doc replace(replace(ParallelKernel.PKSHOW_DOC,     "@init_parallel_kernel" => "@init_parallel_stencil"), "pk_show"    => "ps_show")    macro ps_show(args...)     check_initialized(); esc(:(ParallelStencil.ParallelKernel.@pk_show($(args...)))); end
@doc replace(replace(ParallelKernel.PKPRINTLN_DOC,  "@init_parallel_kernel" => "@init_parallel_stencil"), "pk_println" => "ps_println") macro ps_println(args...)  check_initialized(); esc(:(ParallelStencil.ParallelKernel.@pk_println($(args...)))); end

"""
    @init_parallel_stencil(package, numbertype, ndims)

Initialize the package ParallelStencil, giving access to its main functionality. Creates a module `Data` in the module where `@init_parallel_stencil` is called from. The module `Data` contains the types `Data.Number`, `Data.Array` and `Data.DeviceArray` (type `?Data` *after* calling `@init_parallel_stencil` to see the full description of the module).

# Arguments
- `package::Module`: the package used for parallelization (CUDA or Threads).
- `numbertype::DataType`: the type of numbers used by @zeros, @ones and @rand and in all array types of module `Data` (e.g. Float32 or Float64). It is contained in `Data.Number` after @init_parallel_stencil.
- `ndims::Integer`: the number of dimensions used for the stencil computations in the kernels (1, 2 or 3).

See also: [`Data`](@ref)
"""
macro init_parallel_stencil(args...)
    posargs, kwargs_expr = split_args(args)
    if (length(args) > 3)            @ArgumentError("too many arguments.")
    elseif (0 < length(posargs) < 3) @ArgumentError("there must be either three or zero positional arguments.")
    end
    kwargs = split_kwargs(kwargs_expr)
    if (length(posargs) == 3) package, numbertype_val, ndims_val = extract_posargs_init(__module__, posargs...)
    else                      package, numbertype_val, ndims_val = extract_kwargs_init(__module__, kwargs)
    end
    if (package == PKG_NONE) @ArgumentError("the package argument cannot be ommited.") end #TODO: this error message will disappear, once the package can be defined at runtime.
    if (ndims == NDIMS_NONE) @ArgumentError("the ndims argument cannot be ommited.") end #TODO: this error message will disappear, once the ndims can be defined at runtime.
    check_already_initialized(package, numbertype_val, ndims_val)
    esc(init_parallel_stencil(__module__, package, numbertype_val, ndims_val))
end

function init_parallel_stencil(caller::Module, package::Symbol, numbertype::DataType, ndims::Integer)
    datadoc_call = :(@doc replace(ParallelStencil.ParallelKernel.DATA_DOC, "@init_parallel_kernel" => "@init_parallel_stencil") Data)
    ParallelKernel.init_parallel_kernel(caller, package, numbertype; datadoc_call=datadoc_call)
    set_package(package)
    set_numbertype(numbertype)
    set_ndims(ndims)
    set_initialized(true)
    return nothing
end

macro is_initialized() is_initialized() end
macro get_package() get_package() end
macro get_numbertype() get_numbertype() end
macro get_ndims() get_ndims() end
let
    global is_initialized, set_initialized, set_package, get_package, set_numbertype, get_numbertype, set_ndims, get_ndims, check_initialized, check_already_initialized
    _is_initialized::Bool       = false
    package::Symbol             = PKG_NONE
    numbertype::DataType        = NUMBERTYPE_NONE
    ndims::Integer              = NDIMS_NONE
    set_initialized(flag::Bool) = (_is_initialized = flag)
    is_initialized()            = _is_initialized
    set_package(pkg::Symbol)    = (package = pkg)
    get_package()               = package
    set_numbertype(T::DataType) = (numbertype = T)
    get_numbertype()            = numbertype
    set_ndims(n::Integer)       = (ndims = n)
    get_ndims()                 = ndims
    check_initialized()         = if !is_initialized() @NotInitializedError("no macro or function of the module can be called before @init_parallel_stencil.") end

    function check_already_initialized(package::Symbol, numbertype::DataType, ndims::Integer)
        if is_initialized()
            if package==get_package() && numbertype==get_numbertype() && ndims==get_ndims()
                @warn "ParallelStencil has already been initialized, with the same arguments. If you are using ParallelStencil interactively in the REPL, then you can ignore this message. If you are using ParallelStencil non-interactively, then you are likely using ParallelStencil in an inconsistent way: @init_parallel_stencil should only be called once, right after 'using ParallelStencil'."
            else
                @IncoherentCallError("ParallelStencil has already been initialized, with different arguments. If you are using ParallelStencil interactively in the REPL and want to avoid restarting Julia, then you can call ParallelStencil.@reset_parallel_stencil() and rerun all parts of your code that use ParallelStencil features (including kernel definitions and array allocations). If you are using ParallelStencil non-interactively, then you are using ParallelStencil in an invalid way: @init_parallel_stencil should only be called once, right after 'using ParallelStencil'.")
            end
        end
    end
end

function extract_posargs_init(caller::Module, package, numbertype, ndims) # NOTE: this function takes not only symbols: numbertype can be anything that evaluates to a type in the caller and for package will be checked wether it is a symbol in check_package and a proper error message given if not.
    package, numbertype_val = extract_posargs_init(caller, package, numbertype)
    ndims_val = eval_arg(caller, ndims)
    check_ndims(ndims_val)
    return package, numbertype_val, ndims_val
end

function extract_kwargs_init(caller::Module, kwargs::Dict)
    package, numbertype_val = ParallelKernel.extract_kwargs_init(caller, kwargs)
    if (:ndims in keys(kwargs)) ndims_val = eval_arg(caller, kwargs[:ndims]); check_ndims(ndims_val)
    else                        ndims_val = NUMBERTYPE_NONE
    end
    return package, numbertype_val, ndims_val
end


text:
module EAGOParametricInterval

using EAGOIntervalArithmetic
using IntervalArithmetic

type Param_Bisect_Opts
  DAGflag::Bool
  LPflag::Bool
  kmax_main
  kmax_cntr
  style
  display
  ptol
  etol
  rtol
  DAGpass
  p_rel_bisect
  DAGh
  DAGg
  DAGsym
end
Param_Bisect_Opts() = Param_Bisect_Opts(false, #DAGflag::Bool
                                        false, #LPflag::Bool
                                        1E2, #kmax_main
                                        1E2, #kmax_cntr
                                        "KrawczykCW", #style
                                        "Full", #display
                                        1E-3, #ptol
                                        1E-3, #etol
                                        1E-6, #rtol
                                        5, #DAGpass
                                        false, # p_rel_bisect
                                        [],
                                        [],
                                        [])

include("lib/Sparse_Conditioner.jl")
include("lib/Parametric_Utility.jl")
include("lib/Parametric_Contractor.jl")
#include("lib/Parametric_Test.jl")
#include("lib/Parametric_Bisection.jl") # Used for generalized bisection
#include("lib/Parametric_Main.jl")

#=
export Param_Bisect_Opts,setprec, Miranda, MirandaExc, partialIncTop,
       partialIncBot, Strict_XinY, isEqual, extDivide, extProcess,
=#

export  PI_NewtonGS, PI_KrawczykCW,
        PIn_NewtonGS, PIn_KrawczykCW,
        PId_NewtonGS, PId_KrawczykCW,
        Sparse_Precondition!, SparseInSto

end


text:
using Test
using Rectangle

include("bstutils.jl")

@testset "Binary Search Trees" begin
    t = BinarySearchTree{Int, Int}()
    for i in 1:100
        insert!(t, i, 10i)
    end
    @test depth(t) == 100

    a = [4, 5, 3, 1, 10, 2, 7, 6, 8, 9]
    t = BinarySearchTree{Int, Int}()

    @test get!(t, 1, 2) == get(t, 1, 3)
    @test endswith(string(t), "Tree with 1 nodes.\nRoot at: 1.\n")
    @test get!(t, 2, 3) == get(t, 2, 4)
    @test get!(t, 2, 4) == 3
    @test get(t, 2, 4) == 3
    @test delete!(t, 1) == (1 => 2)
    @test delete!(t, 2) == (2 => 3)
    @test isempty(t)

    for i in a
        insert!(t, i, i*10)
    end
    @test collect(Iterator(t)) == [1=>10, 2=>20, 3=>30, 4=>40, 5=>50,
                                   6=>60, 7=>70, 8=>80, 9=>90, 10=>100]

    b = [i for i in 1:100]
    for i = 1:50
        n1 = rand(b)
        n2 = rand(b)
        b[n1], b[n2] = b[n2], b[n1]
    end
    s = BinarySearchTree{Int, Int}()
    for i in b
       insert!(s, i, 10i)
    end
    bb = [i for i in 1:10]
    @test collect(Iterator(s, -3, 10)) == Pair.(bb, 10bb)
    bb = [i for i in 11:20]
    @test collect(Iterator(s, bb[1], bb[end])) == Pair.(bb, 10bb)
    bb = [i for i in 101:103]
    @test collect(Iterator(s, bb[1], bb[end])) == []
    bb = [i for i in 99:103]
    @test collect(Iterator(s, bb[1], bb[end])) == [99 => 990, 100 => 1000]
    bb = [i for i in -3:2]
    @test collect(Iterator(s, bb[1], bb[end])) == [1 => 10, 2 => 20]
    bb = [i for i in -3:0]
    @test collect(Iterator(s, bb[1], bb[end])) == []
    for i = 1:10
        insert!(s, i, 10i)
    end
    @test collect(Iterator(s, 1, 2)) == [1 => 10, 1=>10, 2=>20, 2 => 20]
    @test length(t) == 10
    @test !isempty(t)
    @test maximum(t) == (10 =>100)
    @test minimum(t) == (1 => 10)
    data = tree_get_data(t)
    @test data[1] == sort(a)
    @test delete!(t, 5) == (5 => 50)
    @test delete!(t, 25) === nothing
    data = tree_get_data(t)
    @test data[2] == [10, 20, 30, 40, 60, 70, 80, 90, 100]
    empty!(t)
    @test isempty(t)
    @testset "Non-unique" begin
        t.unique = false
        n = 1000
        a = []
        for i = 1:n
            j = rand(1:10000)
            push!(a, j, i)
        end
        for i = 1:n
            insert!(t, a[i], i)
        end
        @test parentvalidity(t)
        @test bstvalidity(t)
        for i = 1:div(n, 2)
            delete!(t, a[i])
        end
        @test parentvalidity(t)
        @test bstvalidity(t)
        for i = 1:div(n, 2)
            insert!(t, a[i], i)
        end
        @test parentvalidity(t)
        @test bstvalidity(t)
    end
    @testset "Unique" begin
        t = BinarySearchTree{Int, Int}()
        t.unique = true
        n = 1000
        b = [i for i in 1:n]
        for i = 1:div(n, 2)
            x = rand(1:n)
            y = rand(1:n)
            b[x], b[y] = b[y], b[x]
        end
        a = Tuple{Int, Int}[]
        for i = 1:n
            insert!(t, b[i], i)
            push!(a, (b[i], i))
        end
        ka, kv = tree_get_data(t)
        aa = sort(a, lt= (x, y) -> x[1] < y[1])
        @test aa == collect(zip(ka, kv))
        @test bstvalidity(t)
        @test parentvalidity(t)
        @test_throws ErrorException insert!(t, 10, 10)

        # Randomly delete 50% of the data. Data in a are not sorted

        n_2 = div(n, 2)
        for i = 1:n_2
            Base.delete!(t, b[i])
        end
        bb = b[n_2+1:end]
        a = a[n_2+1:end]
        ka, kv = tree_get_data(t)
        aa = sort(a, lt= (x, y) -> x[1] < y[1])
        @test aa == collect(zip(ka, kv))
        @test bstvalidity(t)
        @test parentvalidity(t)
    end
end

@testset "Red and Black Trees" begin

    a = [4, 5, 3, 1, 10, 2, 7, 6, 8, 9]
    t = RBTree{Int, Int}()

    @test isempty(t)

    for i in a
        insert!(t, i, i*10)
    end
    @test collect(Iterator(t)) == [1=>10, 2=>20, 3=>30, 4=>40, 5=>50,
                                   6=>60, 7=>70, 8=>80, 9=>90, 10=>100]
    @test collect(Iterator(t, 3, 8)) == [3=>30, 4=>40, 5=>50,
                                         6=>60, 7=>70, 8=>80]
    @test parentvalidity(t)
    @test bstvalidity(t)
    @test rbvalidity(t)
    @test length(t) == 10
    @test !isempty(t)
    @test maximum(t) == (10 =>100)
    @test minimum(t) == (1 => 10)
    data = tree_get_data(t)
    @test data[1] == sort(a)
    @test delete!(t, 5) == (5 => 50)
    @test parentvalidity(t)
    @test bstvalidity(t)
    @test rbvalidity(t)
    data = tree_get_data(t)
    @test data[2] == [10, 20, 30, 40, 60, 70, 80, 90, 100]
    empty!(t)
    @test isempty(t)

    @testset "Non-unique" begin
        t.unique = false
        n = 10000
        a = Int[]
        for i = 1:n
            j = rand(1:1000)
            push!(a, j)
        end
        for i = 1:n
            insert!(t, a[i], i)
        end
        @test depth(t) <= 2*log2(n)
        @test parentvalidity(t)
        @test bstvalidity(t)
        @test rbvalidity(t)
        delete!(t, t.root.k)
        for i = 1:div(n, 2)
            delete!(t, a[i])
        end
        @test rbvalidity(t)
        @test parentvalidity(t)
        @test bstvalidity(t)

        for i = 1:div(n, 2)
            insert!(t, a[i], i)
        end
        @test rbvalidity(t)
        @test parentvalidity(t)
        @test bstvalidity(t)
        @test depth(t) <= 2*log2(n)
    end
    empty!(t)
    @testset "Unique" begin
        t.unique = true
        n = 10000
        b = [i for i in 1:n]

        for i = div(n, 2):-1:1
            x = rand(1:n)
            y = rand(1:n)
            b[x], b[y] = b[y], b[x]
        end
        a = Tuple{Int, Int}[]
        pa = Pair{Int, Int}[]
        for i = 1:n
            insert!(t, b[i], i)
            push!(a, (b[i], i))
            push!(pa, (b[i] => i))
        end
        @test depth(t) <= 2*log2(n)
        ka, kv = tree_get_data(t)
        aa = sort(a, lt= (x, y) -> x[1] < y[1])
        pa = sort(pa, lt= (x, y) -> x[1] < y[1])
        @test aa == collect(zip(ka, kv))
        @test pa == collect(Iterator(t))
        @test rbvalidity(t)
        @test bstvalidity(t)
        @test parentvalidity(t)
        #@test_throws ErrorException insert!(t, 10, 10)

        # Randomly delete 50% of the data. Data in a are not sorted
        n_2 = div(n, 2)
        for i = 1:n_2
            Base.delete!(t, b[i])
            ka, kv = tree_get_data(t)
        end
        bb = b[n_2+1:end]
        a  = a[n_2+1:end]
        ka, kv = tree_get_data(t)
        aa = sort(a, lt= (x, y) -> x[1] < y[1])
        @test aa == collect(zip(ka, kv))
        @test rbvalidity(t)
        @test bstvalidity(t)
        @test parentvalidity(t)
    end
end


text:
"""
    ThreadedScans.dac!(op, xs::AbstractVector) -> xs

Compute inclusive scan. Use a divide-and-conquer strategy for spawning the
tasks.  Intermediate reductions are propagated sequentially.  It is slightly
better than `ThreadedScans.linear!`.

# Keyword arguments
* `ntasks = Threads.nthreads()`: number of tasks
"""
function ThreadedScans.dac!(
    op,
    xs;
    ntasks::Integer = Threads.nthreads(),
    spin::Union{Nothing,Integer} = nothing,
)
    if length(xs) <= 2
        if length(xs) == 2
            xs[begin+1] = op(xs[begin], xs[begin+1])
        end
        return xs
    end
    if ntasks < 2
        scan!(op, xs)
        return xs
    end
    nchunks = min(length(xs), ntasks + 1)
    ntasks = nchunks - 1
    @assert ntasks > 1

    chunks = partitionto(xs, nchunks)
    promises = [OneWayPromise{eltype(xs)}() for _ in 1:nchunks]
    dac_spawn_foreach(1:nchunks) do i
        if i == 1
            put!(promises[1], scan!(op, chunks[1]))
        else
            local accr = foldl(op, chunks[i])
            local accl = take!(promises[i-1], spin)
            put!(promises[i], op(accl, accr))
            scan!(op, accl, chunks[i])
        end
    end

    return xs
end


text:
struct UnscentedKF
	xk::Vector
	xp::Vector
	P::AbstractMatrix
	K::AbstractMatrix
	w_mean::Vector
	w_cov::Vector
	X::AbstractMatrix
	w_mean_::Vector
	w_cov_::Vector
	X_::AbstractMatrix
	Pyy::AbstractMatrix
	Pxy::AbstractMatrix
	filter::SquareRootFilter
end

"""
	sol = unscented_kalman(
        y, x₀, f, g, P, R, Q;
        α=1e-3, κ=1.0, β=2.0, sqrt_filter=CholeskySqrt(),
        )

Performs the non-linear adaptive unscented Kalman filter algorithm for the estimation of
a state y. This is the additive (zero mean) noise case.

y:  Vector with one measurement (observation) of m = length(x) variables
x₀: A priori estimation of the variables
f:	State Model
g:	Measurement Model
P:  Initial predicted process covariance matrix
R:  Measurement (observable) covariance matrix
Q:  Process noise covariance matrix

α:  Distance between the sample point and mean point: 1 ≤ α ≤ 1e-3 (default)
κ:  Scaling factor (default 1.0)
β:  Prior knowledge of the distribution (default 2.0, Gaussian process)

sqrt_filter: Square root filtering technique for the covariance matrix.
			 Accepted values are: CholeskyModSqrt() (default), SVDSqrt(), CholeskySqrt()

sol: UnscentedKF type
	xk:  	 Kalman filter state estimation
	xp:  	 Predicted state estimate
    P:   	 Process covariance matrix a posteriori
    K:   	 Kalman gain
	w_mean:  Weights of the mean
	w_cov: 	 Weights of the covariance
	S:		 Sigma matrix, size = (num_variables, 2*num_variables+1)
	w_mean_: Weights of the mean expanded (updated process)
	w_cov_:  Weights of the covariance expanded (updated process)
	S_:      Sigma matrix expanded (updated process), size = (#variables, 2*(2*num_variables+1))
	Pyy:     Covariance matrix measurement upate
	Pxy:	 Cross-covariance matrix measurement upate
	filter:  Filter used for the covariance matrix

"""
function unscented_kalman(
	x::Vector,
    x₀::Vector,
    f::Function,
    g::Function,
    P::AbstractMatrix,
	R::AbstractMatrix,
	Q::AbstractMatrix;
	α::Real=1e-3,
	κ::Real=3-length(x),
	β::Real=2.0,
	sqrt_filter::SquareRootFilter=CholeskyModSqrt(),
	)

	m = length(x₀)

	## Predict process
	# Sigma points and weights
	_X = _sigma_points(x₀, P, sqrt_filter, α, κ, β, m)
	w_m, w_c = _mean_covariance_weights(m; α=α, β=β, κ=κ)
	# Propagate sigma points through the non-linear function
	X_ = _propagate_sigma_nlf(f, _X)
	# Predicted state estimation, a priori,
	x̂ = _estimation(X_, w_m)
	# Predicted process covariance matrix, a priori
	P_ = _covariance_matrix_ukf(P, X_, x̂, w_c, Q)

	## Update process
	# Sigma points a posteriori
	_Y = _sigma_points(X_, R, sqrt_filter, α, κ, β, size(X_, 1))
	w_m_, w_c_ = _mean_covariance_weights(2*m; α=α, β=β, κ=κ)
	# Measurement estimation
	Y_ = _propagate_sigma_nlf(g, _Y)
	ŷ = _estimation(Y_, w_m_)
	# Covariance matrix
	Pyy = _covariance_matrix_ukf(P, Y_, ŷ, w_c_, R)
	# Cross-covariance matrix
	Pxy = _covariance_matrix_ukf(P, _Y, x̂, Y_, ŷ, w_c_)
	# Kalman gain
	K = _kalman_gain(P, Pxy, Pyy)
	# State estimation measurement update (a posteriori)
	xk = muladd(K, (x .- ŷ), x̂)
	# Covariance matrix measurement update (a posteriori)
	P_ .= P_ .- K*Pyy*K'

	return UnscentedKF(
		vec(xk), vec(x̂),
		P_, K,
		vec(w_m), vec(w_c), X_,
		vec(w_m_), vec(w_c_), Y_,
		Pyy, Pxy,
		sqrt_filter,
		)
end

function _sigma_points(
	x::Vector, P::AbstractMatrix, filter, α::Real, β::Real, κ::Real, m::Int64,
	)
	λ = α^2 * (m + κ) - m
	Pi = _square_root_filter(filter, sqrt(m + λ) .* P)
	return hcat(x, x .+ Pi, x .- Pi)
end

function _sigma_points(
	X::Matrix, P::AbstractMatrix, filter, α::Real, β::Real, κ::Real, m::Int64,
	)
	λ = α^2 * (m + κ) - m
	Pi = _square_root_filter(filter, sqrt(m + λ).*P)
	return hcat(X, X[:,1] .+ Pi, X[:,1] .- Pi)
end

function _mean_covariance_weights(m::Int64; α::Real=1e-3, β::Real=2.0, κ::Real=1.0)
	λ = α^2 * (m + κ) - m
	w_m = ones(2*m + 1) .* 0.5 / (m + λ) # mean weights
	w_c = copy(w_m) # covariance weights
	w_m[1] = λ / (m + λ)
	w_c[1] += one(β) + β - α^2
	return w_m, w_c
end

_propagate_sigma_nlf(f::Function, χ::Array) = (x -> f(x))(χ)

_estimation_product(X::Matrix, w::Vector) = ((A, x) -> x' .* A)(X, w)
function _estimation(X::Matrix, w::Vector)
	return vec(sum(_estimation_product(X::Matrix, w::Vector), dims=2))
end

function _covariance_matrix_ukf(P::AbstractMatrix, X::Matrix, x::Vector, w::Vector, N)
	A = X .- x
	return _return_type(P, w' .* A*A' .+ N)
end

function _covariance_matrix_ukf(
	P::AbstractMatrix, X::Matrix, x::Vector, Y::Matrix, y::Vector, w::Vector,
	)
	return _return_type(P, (w' .* (X .- x)) * (Y .- y)')
end

function _kalman_gain(P::AbstractMatrix, Pxy::AbstractMatrix, Pyy::AbstractMatrix)
	F = qr(Pyy)
	return _return_type(P, F.R \ (F.Q' * Pxy))
end


text:
"""
Tuning parameter for metropolis step of a continuous parameter
with support on the real line.
"""
mutable struct TuningParam{T}
  value::T
  acceptanceCount::Int
  currentIter::Int
  batchsize::Int
  delta::Function
  targetAcc::Float64

  TuningParam(v::V) where V = new{V}(v, 0, 0, 50,
                                     n::Int -> max(n^(-0.5), 0.01),
                                     0.44)
end

"""
Performs update of the tuning parameter supplied, given whether
the proposed value was accepted or not.
"""
function update_tuning_param_default(param::TuningParam, accept::Bool)
  if accept
    param.acceptanceCount += 1
  end

  param.currentIter += 1

  if param.currentIter % param.batchsize == 0
    n = Int(floor(param.currentIter / param.batchsize))
    factor = exp(param.delta(n))
    if acceptanceRate(param) > param.targetAcc
      param.value *= factor
    else
      param.value /= factor
    end

    param.acceptanceCount = 0
  end

  return
end

function acceptanceRate(param::TuningParam)
  return param.acceptanceCount / param.batchsize
end


text:
module ContextualDispatch

using MacroTools: @capture, postwalk, unblock, rmlines
using Mixtape
import Mixtape: allow,
                transform,
                optimize!,
                CompilationContext,
                @load_abi,
                jit
using CodeInfoTools

# Controls when the transform on lowered code applies.
# I don't want to apply the recursive transform in type inference.
# So when stacklevel > 1, don't apply.
# I'll just dispatch to `Mixtape.call` at runtime.
mutable struct Mix{T} <: CompilationContext
    stacklevel::Int
    ctx::T
end
Mix(ctx::T) where T = Mix(1, ctx)
Mix{T}() where T = Mix(1, T())

# Basically == a Cassette context.
abstract type Context end

# Our version of overdub.
overdub(::Context, f, args...) = f(args...)
allow(::Mix, m::Module, args...) = true

# The transform inserts state, then wraps calls in (overdub).
# Then, anytime there's a return value --
# create a tuple of (ret, state)
# and returns that.
# Consider monadic lifting f: R -> R => trans => R -> (R, state).
swap(r, e) = e
function swap(r, e::Expr)
    if e.head == :(=)
        return Expr(:(=), e.args[1], swap(r, e.args[2]))
    elseif e.head == :call
        return Expr(:call, overdub, r, e.args[1:end]...)
    end
    return e
end

prehook!(ctx::Context, b, sig) = CodeInfoTools.identity(b)
posthook!(ctx::Context, b, sig) = CodeInfoTools.identity(b)

# Potentially can be sped up. Profile.
function transform(mix::Mix{T}, src, sig) where T
    mix.stacklevel == 1 || return src
    prebuilder = CodeInfoTools.Builder(src)
    prehook!(mix.ctx, prebuilder, sig)
    new = CodeInfoTools.finish(prebuilder)
    b = CodeInfoTools.Builder(new)
    q = push!(b, Expr(:call, T))
    rets = Any[]
    for (v, st) in b
        b[v] = swap(q, st)
        st isa Core.ReturnNode && push!(rets, v => st)
    end
    for (n, ret) in rets
        v = insert!(b, n, Expr(:call, Base.tuple, ret.val, q))
        b[n] = Core.ReturnNode(v)
    end
    new = CodeInfoTools.finish(b)
    postbuilder = CodeInfoTools.Builder(new)
    posthook!(mix.ctx, postbuilder, sig)
    new = CodeInfoTools.finish(postbuilder)
    mix.stacklevel += 1
    return new
end

macro load_call()
    expr = quote
        using Mixtape
        ContextualDispatch.@load_abi()
        call(ctx::T, fn, args...) where T <: Context = call(fn, args...; ctx = Mix(ctx))
    end
    esc(expr)
end

export overdub, Context, @load_call,
       Mix, jit, prehook!, posthook!

end # module