text:
# do some profiling
using Plots, LaTeXStrings
gr(colorbar=false,size=(600,150),legend=false,grid=false,xticks=false,yticks=false,axis=false)
msg =:mediumseagreen

using Revise, FourierGPE

# ==== Initialize simulation
L = (60.0,)
N = (512,)
sim = Sim(L,N)
@unpack_Sim sim

# ==== set simulation parameters
μ = 25.0
g = -0.01
γ = 0.0
Ns = 200
ξs = 2/abs(g)/Ns
us = 20
tf = π |> Float64
Nt = 150
ti = 0.0
t = LinRange(ti,tf,Nt)

# ==== soliton wavefunction
x = X[1]
ψs(x) = sqrt(Ns/2ξs)*sech(x/ξs)*exp(im*us*x)
ψi = ψs.(x)
ϕi = kspace(ψi,sim)

# ==== Set all fields
@pack_Sim! sim

# ====== Evolve in k space ==========
sol = runsim(sim)
# ===================================

import FourierGPE:Lgp!, nlin!, V
V(x,t) = 0.0

function nlin!(dϕ,ϕ,sim::Sim{1},t)
    @unpack g,X,V0 = sim; x = X[1]
    dϕ .= ϕ
    xspace!(dϕ,sim)
    @. dϕ *= V0 + V(x,[t]) + g*abs2(dϕ)
    kspace!(dϕ,sim)
    return nothing
end

function Lgp!(dϕ,ϕ,sim,t)
    @unpack γ,μ,espec = sim
    nlin!(dϕ,ϕ,sim,t)
    @. dϕ = -im*(1.0 - im*γ)*(dϕ + (espec - μ)*ϕ)
    return nothing
end

using BenchmarkTools

dϕi = copy(ϕi)
@btime Lgp!(dϕi,ϕi,sim,t)
# factor of 10 speed up in MKL in 1D
# if there isn't a speed up, there must be another bottleneck

@btime kspace!(dϕi,sim)
@btime xspace!(dϕi,sim)
@btime @. dϕi *= g*abs2(dϕ)

# ===============================
# 2d profiling: vortex precession example

using Revise, FourierGPE

# ==== Initialize simulation
L = (18.0,18.0)
N = (128,128)
sim = Sim(L,N)
@unpack_Sim sim

# ==== set simulation parameters
μ = 12.0

# Time dependent potential function (here trivial t dep)
import FourierGPE.V
V(x,y,t) = 0.5*(x^2 + y^2)
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)

# ==== make initial state
x,y = X
ψi = ψ0.(x,y',μ,g)
ϕi = kspace(ψi,sim)
@pack_Sim! sim

# ==== evolve
@time sol = runsim(sim)

# ==== ground state

ϕg = sol[end]
ψg = xspace(ϕg,sim)
showpsi(x,y,ψg)

# ==== set simulation parameters
γ = 0.0
t = LinRange(ti,tf,Nt)
ϕi = kspace(ψg,sim)
# reltol = 1e-7
# alg = DP5()ß

#---
using VortexDistributions
R(w) = sqrt(2*μ/w^2)
R(1)
rv = 3.
healinglength(x,y,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,μ,g)))
ξ0 = healinglength.(0.,0.,μ,g)
ξ = healinglength(rv,0.,μ,g)
#---

pv = PointVortex(rv,0.,1)
vi = ScalarVortex(ξ,pv)
psi = Torus(copy(ψg),x,y)
vortex!(psi,vi)
showpsi(x,y,psi.ψ)
ψi .= psi.ψ
ϕi = kspace(ψi,sim)

# ===== compare with Fetter JLTP 2010
ξ = 1/sqrt(μ)
Rtf = R(1)
Ωm = 3*log(Rtf/ξ/sqrt(2))/2/Rtf^2
Ωv = Ωm/(1-rv^2/Rtf^2)

# or a precession period
Tv = 2*pi/Ωv

ti = 0.; tf = 0.1*Tv
t = LinRange(ti,tf,Nt)

@pack_Sim! sim
#---
# ==== evolve
@time solv = runsim(sim)

# FFTW 1.8s
# MKL 1.1s

import FourierGPE:Lgp!, nlin!, V
V(x,t) = 0.0

function nlin!(dϕ,ϕ,sim::Sim{1},t)
    @unpack g,X,V0 = sim; x = X[1]
    dϕ .= ϕ
    xspace!(dϕ,sim)
    @. dϕ *= V0 + V(x,[t]) + g*abs2(dϕ)
    kspace!(dϕ,sim)
    return nothing
end

function Lgp!(dϕ,ϕ,sim,t)
    @unpack γ,μ,espec = sim
    nlin!(dϕ,ϕ,sim,t)
    @. dϕ = -im*(1.0 - im*γ)*(dϕ + (espec - μ)*ϕ)
    return nothing
end

using BenchmarkTools

dϕi = copy(ϕi)
@btime Lgp!(dϕi,ϕi,sim,t)
# factor of ?? speed up in MKL in 2D
# if there isn't a speed up, there must be another bottleneck

@btime kspace!(dϕi,sim)
@btime xspace!(dϕi,sim)
@btime @. dϕi *= g*abs2(dϕ)


text:
import ArgParse, Gurobi, DataFrames
using ArgCheck

const arg_parse_settings = ArgParse.ArgParseSettings(prog="var", fromfile_prefix_chars='@',
                                                     description="Calculates the steady state as function of ξ.")

@ArgParse.add_arg_table arg_parse_settings begin
    "--sto"
        help = "stoichiometric matrix"
        arg_type = String
        required = true
    "--mets"
        help = "metabolites"
        arg_type = String
        required = true
    "--rxns"
        help = "reactions"
        arg_type = String
        required = true
    "--A"
        help = "biomass synthesis cost"
        arg_type = Float64
        default = 0.
        range_tester = A -> 0 ≤ A < Inf
    "--c"
        help = "media"
        arg_type = String
        required = true
    "--xi"
        help = "File with a list of ξ values to use. Execution stops at the first infeasible value of ξ found."
        arg_type = String
        required = true
    "--out"
        help = "output file"
        arg_type = String
        required = true
    "--prefix"
        help = "all paths arguments (--sto, --mets, --rxns, --c, --xi, --out) will be prefixed with this string."
        arg_type = String
        default = ""
end

const args = ArgParse.parse_args(arg_parse_settings)

set_multiobj_n!(gurobi_model::Gurobi.Model, n::Int) = Gurobi.set_intattr!(gurobi_model, "NumObj", n)

function set_multiobj!(gurobi_model::Gurobi.Model, i::Int, c, priority, weight)
    @assert length(c) == Gurobi.num_vars(gurobi_model)
    oldi = Gurobi.get_int_param(gurobi_model, "ObjNumber")
    Gurobi.set_int_param!(gurobi_model, "ObjNumber", i)
    Gurobi.set_dblattrarray!(gurobi_model, "ObjN", 1, Gurobi.num_vars(gurobi_model), c)
    Gurobi.set_intattr!(gurobi_model, "ObjNPriority", priority)
    Gurobi.set_dblattr!(gurobi_model, "ObjNWeight", weight)
    Gurobi.set_int_param!(gurobi_model, "ObjNumber", oldi)
end

function set_ub!(gurobi_model::Gurobi.Model, u::Vector{Float64})
    @assert length(u) == Gurobi.num_vars(gurobi_model)
    Gurobi.set_dblattrarray!(gurobi_model, "UB", 1, Gurobi.num_vars(gurobi_model), u)
end

function set_rhs!(gurobi_model::Gurobi.Model, rhs::Vector{Float64})
    @assert Gurobi.num_constrs(gurobi_model) == length(rhs)
    Gurobi.set_dblattrarray!(gurobi_model, "RHS", 1, Gurobi.num_constrs(gurobi_model), rhs)
end

function GurobiEnv()
    env = Gurobi.Env()
    # suppress Gurobi logging output
    Gurobi.setparam!(env, "OutputFlag", 0)
    Gurobi.setparam!(env, "Quad", 1)
    Gurobi.setparam!(env, "Method", 4)
    #Gurobi.setparam!(env, "MultiObjMethod", 1)
    Gurobi.setparam!(env, "FeasibilityTol", 1e-9)
    Gurobi.setparam!(env, "OptimalityTol", 1e-9)
    Gurobi.setparam!(env, "MarkowitzTol", 0.9)
    #Gurobi.setparam!(env, "LogToConsole", 0)
    #Gurobi.setparam!(env, "LogFile", path_prefix * "gurobi.log"))
    return env
end

function solve!(gurobi_model::Gurobi.Model)::Bool
    Gurobi.update_model!(gurobi_model)
    Gurobi.optimize(gurobi_model)
    return Gurobi.get_status(gurobi_model) == :optimal
end

function U(V, c, ξ)::Float64
    @argcheck V ≥ 0 && c ≥ 0 && 0 ≤ ξ < Inf
    return c == 0 ? 0. : min(V, c / ξ)
end

function main()
    path_prefix = strip(args["prefix"])

    ξlist = readdlm(path_prefix * strip(args["xi"]))
    @assert all(0 .≤ ξlist .< Inf)

    sto = DataFrames.readtable(path_prefix * strip(args["sto"]); allowcomments=true, separator=' ')
    S = sparse(sto[:i], sto[:k], sto[:s])
    m, n = Base.size(S)
    @assert m < n

    mets = DataFrames.readtable(path_prefix * strip(args["mets"]); allowcomments=true, separator=' ')
    rxns = DataFrames.readtable(path_prefix * strip(args["rxns"]); allowcomments=true, separator=' ')
    @assert Base.size(mets, 1) == m && Base.size(rxns, 1) == n

    @assert all(rxns[:lb] .≤ 0 .≤ rxns[:ub])
    @assert all(rxns[:an] .≥ 0) && all(rxns[:ap] .≥ 0) && all(mets[:bp] .≥ 0) && all(mets[:bn] .≥ 0)
    @assert all(mets[:L] .≤ 0 .≤ mets[:V])
    @assert all(-Inf .< mets[:e] .< Inf) && all(-Inf .< mets[:y] .< Inf)

    c = DataFrames.readtable(path_prefix * strip(args["c"]); allowcomments=true, separator=' ')
    @assert all(c[:mM] .≥ 0)
    @assert all(c[:id] == mets[:id])

    x = zeros(2n + 2m + 2)   # [r+, r-, u+, u-, μ, ϕ]
    upper::Vector{Float64} = [rxns[:ub]; -rxns[:lb]; mets[:V]; -mets[:L]; Inf; 1.]     # upper bounds
    @assert length(upper) == length(x) == 2m + 2n + 2

    # equality constrains
    W = sparse([ S            -S              speye(m)          -speye(m)     -mets[:y]       spzeros(m)
                 rxns[:ap]'    rxns[:an]'     mets[:bp]'        mets[:bn]'    args["A"]       -1 ])

    @assert Base.size(W) == (m + 1, length(x))

    # the problem has the form:
    # max μ s.t. W⋅x = [e,0] and 0 ≤ x ≤ upper.
    # with the secondary objective: min ϕv

    LP = Gurobi.gurobi_model(GurobiEnv(), sense = :maximize,
                             Aeq = W, beq = [collect(mets[:e]); 0],
                             lb = zeros(x), ub = upper, f = zeros(x))

    set_multiobj_n!(LP, 2)
    set_multiobj!(LP, 0, [zeros(2n + 2m); 1;  0], 10, 1.0)  # max μ
    set_multiobj!(LP, 1, [zeros(2n + 2m); 0; -1], 5,  1.0)  # min ϕ

    r = zeros(n)
    u = zeros(m)
    s = zeros(m)

    xout = open(path_prefix * strip(args["out"]) * ".x", "w")
    uout = open(path_prefix * strip(args["out"]) * ".u", "w")
    rout = open(path_prefix * strip(args["out"]) * ".r", "w")
    sout = open(path_prefix * strip(args["out"]) * ".s", "w")

    write(xout, "xi\tmu\tphi\n")
    write(uout, "xi"); for met in mets[:id]; write(uout, "\t" * met); end; write(uout, "\n")
    write(rout, "xi"); for rxn in rxns[:id]; write(rout, "\t" * rxn); end; write(rout, "\n")
    write(sout, "xi"); for met in mets[:id]; write(sout, "\t" * met); end; write(sout, "\n")

    for ξ in ξlist
        # update uptake bounds
        view(upper, 2n + (1:m)) .= U.(mets[:V], c[:mM], ξ)
        set_ub!(LP, upper)

        solve!(LP) || return
        x .= Gurobi.get_solution(LP)
        r .= x[1:n] .- x[n + (1:n)]
        u .= x[2n + (1:m)] .- x[2n + m + (1:m)]
        s .= c[:mM] .- u * ξ
        μ = x[2n + 2m + 1]
        ϕ = x[2n + 2m + 2]

        write(xout, "$ξ\t$μ\t$ϕ\n")
        @printf "ξ=%-7.13f μ=%-7.13f ϕ=%-7.13f\n" ξ μ ϕ

        write(uout, "$ξ"); for i = 1:m; write(uout, "\t$(u[i])"); end; write(uout, "\n")
        write(rout, "$ξ"); for k = 1:n; write(rout, "\t$(r[k])"); end; write(rout, "\n")
        write(sout, "$ξ"); for i = 1:m; write(sout, "\t$(s[i])"); end; write(sout, "\n")
    end

    close(xout)
    close(uout)
    close(rout)
    close(sout)
end

main()

ext: jl

max_stars_repo_name: cossio/DFBAeq

max_issues_repo_name: cossio/DFBAeq

max_forks_repo_name: cossio/DFBAeq


text:
"""
System block
"""

export systemequation, odeproblem, scope
export SystemBlock, blk!, connect!, assignvar!
export SimulationResult

mutable struct SystemBlock <: AbstractBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}
    env::Dict{Symbol,Any}
    states::Dict{Symbol,AbstractBlock}

    function SystemBlock()
        new(Vector{InPort}(), Vector{OutPort}(), Dict{Symbol,Any}(), Dict{Symbol,AbstractBlock}())
    end
end

"""
IO
"""

function Base.show(io::IO, blk::SystemBlock)
    println(io, "SystemBlock()")
end

"""
functions
"""

function blk!(system::SystemBlock, label::Symbol, blk::AbstractBlock)::AbstractBlock
    system.env[label] = blk
end

function blk!(system::SystemBlock, label::Symbol, blk::AbstractIntegratorBlock)::AbstractBlock
    system.states[blk.state] = blk
    system.env[label] = blk
end

function blk!(system::SystemBlock, label::Symbol, blk::SystemBlock)::AbstractBlock
    for (k,v) in blk.states
        system.states[k] = v
    end
    for (k,v) in blk.env
        k2 = Symbol(label, "_", k)
        system.env[k2] = blk.env[k]
    end
    system.env[label] = blk
end

function blk!(system::SystemBlock, label::Symbol, blk::InBlock)::AbstractBlock
    push!(system.inport, blk.inport[1])
    system.env[label] = blk
end

function blk!(system::SystemBlock, label::Symbol, blk::OutBlock)::AbstractBlock
    push!(system.outport, blk.outport[1])
    system.env[label] = blk
end

"""
connect
"""

function connect!(src::AbstractSignalLine{Tv}, dest::InPort, ::Type{Tv} = Float64) where Tv
    dest.line = src
    nothing
end

function connect!(src::OutPort, dest::InPort, ::Type{Tv} = Float64) where Tv
    connect!(signal(src, Tv), dest, Tv)
end

function connect!(src::Vector{OutPort}, dest::InPort, ::Type{Tv} = Float64) where Tv
    connect!(src[1], dest, Tv)
end

function connect!(src::AbstractBlock, dest::InPort, ::Type{Tv} = Float64) where Tv
    connect!(src.outport[1], dest, Tv)
end

function connect!(src::Any, dest::Vector{InPort}, ::Type{Tv} = Float64) where Tv
    connect!(src, dest[1], Tv)
end

function connect!(src::Any, dest::AbstractBlock, ::Type{Tv} = Float64) where Tv
    connect!(src, dest.inport[1], Tv)
end

"""
assign var
"""

function assignvar!(system::SystemBlock, label::Symbol, x::Any)
    system.env[label] = x
end

"""
system equation
"""

function systemequation(system::SystemBlock, label::Symbol = :model)
    u = gensym()
    du = Symbol("d", u)
    states = sort(collect(keys(system.states)))
    n = length(states)
    dict = Dict(zip(states, [:($u[$i]) for i = 1:n]))
    expr = [Expr(:(=), :($(du)[$i]), _replace(dict, _tosystemexpr(system.states[states[i]]))) for i = 1:n]
    initvec = [_toexpr(system.states[x].initialcondition) for x = states]
    init = Symbol(label, "init")
    func = Symbol(label, "!")
    quote
        $init = $initvec
        function $func($(du), $(u), p, t)
            $(Expr(:block, expr...))
        end
    end
end

function _systemequation_func(system::SystemBlock)
    u = gensym()
    du = Symbol("d", u)
    states = sort(collect(keys(system.states)))
    n = length(states)
    dict = Dict(zip(states, [:($u[$i]) for i = 1:n]))
    expr = [Expr(:(=), :($(du)[$i]), _replace(dict, _tosystemexpr(system.states[states[i]]))) for i = 1:n]
    func = quote
        ($(du), $(u), p, t) -> begin
            $(Expr(:block, expr...))
        end
    end
    eval(func)
end

function _systemequation_init(system::SystemBlock)
    states = sort(collect(keys(system.states)))
    [eval(_toexpr(system.states[x].initialcondition)) for x = states]
end

struct NoSolution
    t
end

function odeproblem(system::SystemBlock; tspan)
    if length(system.states) > 0
        func = _systemequation_func(system)
        initvec = _systemequation_init(system)
        DifferentialEquations.ODEProblem(func, initvec, tspan)
    else
        NoSolution([x for x = tspan])
    end
end

function DifferentialEquations.solve(x::NoSolution, args...; kargs...)
    x
end

function scope(system::SystemBlock, var::Symbol, solution)
    u = :solvedvalue
    states = sort(collect(keys(system.states)))
    n = length(states)
    dict = Dict(zip(states, [:($u(t)[$i]) for i = 1:n]))
    varx = signal(system.env[var])
    expr = _replace(dict, _toexpr(varx))
    func = eval(quote
        (t, solvedvalue) -> begin
            $(expr)
        end
            end)
    t -> begin
        func(t, solution)
        end
end

function _replace(dict, expr::Expr)
    args = [_replace(dict, x) for x = expr.args]
    Expr(expr.head, args...)
end

function _replace(dict, expr::Symbol)
    if expr in keys(dict)
        :($(dict[expr]))
    else
        expr
    end
end

function _replace(dict, expr::Any)
    expr
end


text:
# A line segment in 3D space defined by its two endpoints.
# For ray tracing purposes, the line starts at points[1] and goes to points[2]
struct LineSegment_3D{T <: AbstractFloat}
    points::NTuple{2, Point_3D{T}}
end

# Constructors
# -------------------------------------------------------------------------------------------------
LineSegment_3D(p₁::Point_3D, p₂::Point_3D) = LineSegment_3D((p₁, p₂))

# Base
# -------------------------------------------------------------------------------------------------
Base.broadcastable(l::LineSegment_3D) = Ref(l)

# Methods
# -------------------------------------------------------------------------------------------------
function (l::LineSegment_3D{T})(r::R) where {T <: AbstractFloat, R <: Real}
    return l.points[1] + T(r) * (l.points[2] - l.points[1])
end

#arc_length(l::LineSegment_3D) = distance(l.points[1], l.points[2])
#
#function intersect(l₁::LineSegment_3D{T}, l₂::LineSegment_3D{T}) where {T <: AbstractFloat}
#    # NOTE: Doesn't work for colinear lines. (v⃗ × u⃗ = 0⃗)
#    #
#    # Using the equation of a line in parametric form
#    # For l₁ = x⃗₁ + rv⃗ and l₂ = x⃗₂ + su⃗
#    # x⃗₁ + rv⃗ = x⃗₂ + su⃗                             subtracting x⃗₁ from both sides
#    # rv⃗ = (x⃗₂-x⃗₁) + su⃗                             w⃗ = x⃗₂-x⃗₁
#    # rv⃗ = w⃗ + su⃗                                   cross product with u⃗ (distributive)
#    # r(v⃗ × u⃗) = w⃗ × u⃗ + s(u⃗ × u⃗)                   u⃗ × u⃗ = 0
#    # r(v⃗ × u⃗) = w⃗ × u⃗                              dot product v⃗ × u⃗ to each side
#    # r(v⃗ × u⃗) ⋅(v⃗ × u⃗) = (w⃗ × u⃗) ⋅(v⃗ × u⃗)          divide by (v⃗ × u⃗) ⋅(v⃗ × u⃗)
#    # r = [(w⃗ × u⃗) ⋅(v⃗ × u⃗)]/[(v⃗ × u⃗) ⋅(v⃗ × u⃗)]
#    # Note that if the lines are parallel or collinear, v⃗ × u⃗ = 0
#    # We need to ensure r, s ∈ [0, 1]. Verifying this condition for r is simple, but we need to
#    # solve for s as well.
#    # x⃗₂ + su⃗ = x⃗₁ + rv⃗                              subtracting x⃗₂ from both sides
#    # su⃗ = -w⃗ + rv⃗                                   we see that each element must satisfy
#    # s(u⃗ ⋅ u⃗) = (-w⃗ + rv⃗) ⋅ u⃗                       hence
#    # s = (rv⃗ - w⃗) ⋅ u⃗/(u⃗ ⋅ u⃗)
#    # If the lines are skew, s and r represent the parameters of the points of closest
#    # approach - Intersection of two lines in three-space, Ronald Goldman, in Graphics
#    # Gems by Andrew S. Glassner.
#    v⃗ = l₁.points[2] - l₁.points[1]
#    u⃗ = l₂.points[2] - l₂.points[1]
#    w⃗ = l₂.points[1] - l₁.points[1]
#    if norm(v⃗ × u⃗) > 5.0e-5
#        r = ((w⃗ × u⃗) ⋅ (v⃗ × u⃗))/((v⃗ × u⃗) ⋅ (v⃗ × u⃗))
#        p = l₁(r)
#        s = (r*v⃗ - w⃗) ⋅ u⃗/(u⃗ ⋅ u⃗)
#        return (0 ≤ s ≤ 1) && (0 ≤ r ≤ 1) ? (true, p) : (false, p)
#    else
#        return (false, Point_3D(T, 0))
#    end
#end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
obcolor(::Obstacle) = (0,0.6,0)
obcolor(::Union{RandomWall, RandomDisk}) = (149/255, 88/255, 178/255)
obcolor(::Union{FiniteSplitterWall, SplitterWall, Antidot, Ellipse}) = (0.8,0.0,0)
obcolor(::PeriodicWall) = (0.8,0.8,0)
obalpha(::Obstacle) = 0.5
obalpha(::Union{Antidot, Ellipse}) = 0.1
obls(::Obstacle) = "solid"
obls(::Union{FiniteSplitterWall, SplitterWall, Antidot, Ellipse}) = "dashed"
obls(::PeriodicWall) = "dotted"

"""
    plot(obst::Obstacle; kwargs...)
Plot given obstacle on the current `PyPlot` axes.

The default arguments for each type of obstacle have been chosen for maximum
clarity and consistency.

The `kwargs...` given by the user are keywords passed directly into PyPlot's
constructors. For `Wall` obstacles, kwargs are passed into `PyPlot.plot()`. For
`Circular` obstacles, kwargs are passed into `matplotlib.patches.Circle` or `Arc`.
"""
function plot(d::Obstacle) end

function plot(d::Circular; kwargs...)
    edgecolor = obcolor(d)
    facecolor = (edgecolor..., obalpha(d))
    circle1 = PyPlot.plt.Circle(d.c, d.r;
        edgecolor = edgecolor, facecolor = facecolor,
        linestyle = obls(d), lw = 2.0, kwargs...)
    PyPlot.gca().add_artist(circle1)
end

function plot(d::Semicircle; kwargs...)
    theta1 = atan(d.facedir[2], d.facedir[1])*180/π + 90
    theta2 = theta1 + 180
    edgecolor = obcolor(d)
    s1 = PyPlot.matplotlib.patches.Arc(d.c, 2d.r, 2d.r, theta1 = theta1, theta2 = theta2, edgecolor = edgecolor,
    lw = 2.0, kwargs...)
    PyPlot.gca().add_artist(s1)
end

function plot(w::Wall; kwargs...)
    if typeof(w) <: FiniteWall &&  w.isdoor
       PyPlot.plot([w.sp[1],w.ep[1]],[w.sp[2],w.ep[2]];
       color="black", linestyle = "-", lw = 2.0, kwargs...)
       PyPlot.plot([w.sp[1],w.ep[1]],[w.sp[2],w.ep[2]];
       color=(0, 0.9, 0.9), linestyle = "--", lw = 2.0, kwargs...)
    else
        PyPlot.plot([w.sp[1],w.ep[1]],[w.sp[2],w.ep[2]];
        color=obcolor(w),
        linestyle = obls(w), lw = 2.0, kwargs...)
    end
end

function plot(e::Ellipse; kwargs...)
    edgecolor = obcolor(e)
    facecolor = (edgecolor..., obalpha(e))
    ellipse = PyPlot.matplotlib.patches.Ellipse(e.c, 2e.a, 2e.b;
        edgecolor = edgecolor, facecolor = facecolor,
        linestyle = obls(e), lw = 2.0, kwargs...)
    PyPlot.gca().add_artist(ellipse)
end


text:
import Base.vec, Base.repmat, Base.reshape
export vec, repmat, reshape

function vec(x::Number)
  return x
end

function repmat(x::Number, m::Int64, n::Int64)
  return repmat([x], m, n)
end

function reshape(x::Number, m::Int64, n::Int64)
  return x
end


text:
using RegularizedRegression: SlopeRegression, learn, coefficients

@testset "Slope regression" begin
    for _ in 1:NumTestRepetitions
        test_learner(SlopeRegression, 15.0) # Seems this needs to be higher for SLOPE than for SCAD and MCP
    end
end


text:
using KernelAbstractions
using AMDGPU
using ROCKernels
using CUDA
using CUDAKernels

# One day for Intel GPUs
# using oneAPI
# using OneKernels

using LinearAlgebra
using Printf
using PowerModels
using LazyArtifacts
using ExaTronKernels
const KA = KernelAbstractions
include("admm/opfdata.jl")
include("admm/environment.jl")
include("admm/generator_kernel.jl")
include("admm/eval_kernel.jl")
include("admm/polar_kernel.jl")
include("admm/bus_kernel.jl")
include("admm/tron_kernel.jl")
include("admm/acopf_admm_gpu.jl")

# `datafile`: the name of the test file of type `String`
# here: MATPOWER case2868rte.m in ExaData project Artifact
datafile = joinpath(artifact"ExaData", "ExaData", "matpower", "case9.m")
# datafile = joinpath(artifact"ExaData", "ExaData", "matpower", "case9241pegase.m")
# datafile = joinpath(artifact"ExaData", "ExaData", "matpower", "case_ACTIVSg70k.m")
# `rho_pq`: ADMM parameter for power flow of type `Float64`
rho_pq = parse(Float64, ARGS[1])
# `rho_va`: ADMM parameter for voltage and angle of type `Float64`
rho_va = parse(Float64, ARGS[2])
# `max_iter`: maximum number of iterations of type `Int`
max_iter = parse(Int, ARGS[3])
# Indicate which GPU device to use
device = CPU()
# device = CUDADevice()
# device = ROCDevice()
# device = OneDevice()
# Use polar formulation for branch problems
use_polar = true
# verbose = 0: No output
# verbose = 1: Final result metrics
# verbose = 2: Iteration output
verbose = 1

println("Running case file: $datafile with")
println("rho_pq = $rho_pq,")
println("rho_va = $rho_va")
println("max_iter = $max_iter")

env = admm_rect_gpu(
    datafile;
    iterlim=max_iter,
    rho_pq=rho_pq,
    rho_va=rho_va,
    scale=1e-4,
    use_polar=use_polar,
    device=device,
    verbose=verbose
)


text:
@testset "context data" begin

    sentences = [["hi", "tests", "world"], ["tests", "are", "sure", "fun"]]
    cs = [Context(3, ["hi", "tests", "world"]), Context(3, ["tests", "are", "sure", "fun"])]
    true_skipgram = word2vecdata(:skipgram)
    true_cbow = word2vecdata(:cbow)
    for context in cs
        push!(true_skipgram, context)
        push!(true_cbow, context)
    end

    skipgram_data = (process_input(sentences, 1, 1, :skipgram)).data

    @test length(skipgram_data) == length(true_skipgram)
    for i in 1:length(skipgram_data)
        @test iterate(skipgram_data, i) == iterate(true_skipgram, i)
    end

    cbow_data = (process_input(sentences, 1, 1, :cbow)).data

    @test length(cbow_data) == length(true_cbow)
    for i in 1:length(cbow_data)
        @test iterate(cbow_data, i) == iterate(true_cbow, i)
    end

    # cbow = CBOWData(:cbow, contexts)
    # skipgram = CBOWData(:skipgram, contexts)
    #
    # @test iterate(skipgram) == (("tests", ["hi", "world"]), 2)
    # @test iterate(cbow) == ((["hi", "world"], "tests"), 2)

end


text:
export plotGrid

"""
function plotGrid

	plots nodal and cell-centered grids

	Requirement: PyPlot must be installed.

	Usage:

	plotGrid(M::AbstractTensorMesh) - plots nodal grid M
	plotGrid(x,M)                   - plots grid defined by x (nodal or cc)

	plotGrid is based on PyPlot.plot and PyPlot.plot3D. All keywords from those functions
	are available and forwarded to those methods.

"""
function plotGrid(y,M;spacing=[1,1,1],color="blue",kwargs...)

    if length(y)==M.dim*prod(M.n)
        nn = M.n
    elseif length(y)==M.dim*prod(M.n .+ 1)
        nn = M.n .+ 1
    else
        error("plotGrid - unknown grid type, length(y)=$(length(y)), n=$(M.n), dim=$(M.dim)")
    end
    Y = reshape(y,(prod(nn),M.dim))
    yi(d) = reshape(Y[:,d],tuple(nn...))
    if M.dim==2
        J1 = 1:spacing[1]:nn[1]; y1 = reshape(yi(1),tuple(nn...));
        J2 = 1:spacing[2]:nn[2]; y2 = reshape(yi(2),tuple(nn...));
        p1 = plot!(y1[:,J2],y2[:,J2]; color=color,legend=nothing,  kwargs...);
        p2 = plot!(y1[J1,:]',y2[J1,:]';color=color,legend=nothing,  kwargs...);
    elseif M.dim==3
        pt1 = p3!(yi(1),yi(2),yi(3),nn,[1 2 3],spacing; color=color,legend=nothing, kwargs...)
		pt2 = p3!(yi(1),yi(2),yi(3),nn,[2,1,3],spacing; color=color,legend=nothing, kwargs...);
        pt3 = p3!(yi(1),yi(2),yi(3),nn,[3,1,2],spacing; color=color,legend=nothing, kwargs...);
	else
		error("plotGrid - cannot handle $(M.dim)-dimensional grids")
    end
end

plotGrid(M::AbstractMesh;kwargs...) = plotGrid(getNodalGrid(M),M;kwargs...)

"""
function p3

	helper function for plotGrid

"""
function p3!(y1,y2,y3,m,order,s;kwargs...)
    m  = m[order];
    s  = s[order];
    y(x)  =  reshape(permutedims(x,order),(m[1],m[2]*m[3]));
    y1    = y(y1); y2 = y(y2); y3 = y(y3);
    dx = Int.(round.(m./s))
    J2 = 1:s[2]:m[2];
    J3 = 1:s[3]:m[3];
    K2 = collect(J2)*ones(Int64,1,length(J3));
    K3 = ones(Int64,length(J2),1)*(m[2]*(collect(J3).-1))'
    K  = vec(K2)+vec(K3);
    for k=1:length(K)
        plot3D!(y1[:,K[k]],y2[:,K[k]], y3[:,K[k]];kwargs...)
    end
end


text:
#= This file is a part of URIValidator. License is MIT:
https://github.com/rafpyprog/URIValidator.jl/blob/master/LICENSE.md =#

__precompile__()

include("urlsplit.jl")

type ValidationError <: Exception
    var::String
end

function scheme_is_valid(scheme)
    #= Validates the Scheme component (RFC3986 item 3.1)

    Scheme names consist of a sequence of characters beginning with a letter
    and  followed  by any combination of letters, digits, plus ("+"), period
    ("."), or hyphen ("-")
    =#
    scheme = lowercase(scheme) # the canonical form is lowercase
    starts_with_letter = ismatch(r"^[a-z]", scheme)
    combination_letters_digits_plus_period_hyphen = ismatch(r"^[a-z0-9+\.-]+$",
                                                            scheme)

    if starts_with_letter & combination_letters_digits_plus_period_hyphen
        is_valid = true
    else
        is_valid = false
    end
   return is_valid
end

function authority_is_valid(authority)
    #= Validates the Scheme component (RFC3986 item 3.2)

    The authority component is preceded by a double slash ("//") and is
   terminated by the next slash ("/"), question mark ("?"), or number
   sign ("#") character, or by the end of the URI.

   authority   = [ userinfo "@" ] host [ ":" port ]
    =#
    x = 1
end

function path_is_valid(path)
    x = 2
    #= If a URI contains an authority component, then the path component
   must either be empty or begin with a slash ("/") character.  =#
end

type Validator
    URI
    scheme
    authority

    function Validator(URI)
        this = new()
        this.URI = URI

        URI_components = URISplit(this.URI)

        # scheme validation
        this.scheme = URI_components.scheme

        if !scheme_is_valid(this.scheme)
            throw(ValidationError("Scheme '" * this.scheme * "' is invalid."))
        end

        #authority validation'
        this.authority = URI_components.authority

        return this
    end
end


text:
# ndgrid-unitr.jl

using LazyGrids: ndgrid, ndgrid_array, GridUR
using Test: @test, @testset, @test_throws, @inferred

@testset "unitr" begin
    (L, M, N) = (5, 6, 7)
    (xa, ya, za) = @inferred ndgrid_array(1:L, 1:M, 1:N)
    (xl, yl, zl) = @inferred ndgrid(1:L, 1:M, 1:N)

    @test eltype(xl) === Int
    @test size(xl) === (L, M, N)
    @test xl isa GridUR
    @test yl isa GridUR
    @test zl isa GridUR
    @test xl == xa
    @test yl == ya
    @test zl == za
    @test xl[2,3,4] == xa[2,3,4] == 2
    @test yl[2,3,4] == ya[2,3,4] == 3
    @test zl[2,3,4] == za[2,3,4] == 4

    x, y, z = (-1:2, -3:4, 5:9)
    (xa, ya, za) = @inferred ndgrid_array(x, y, z)
    (xl, yl, zl) = @inferred ndgrid(x, y, z)

    @test size(xl) === (length(x), length(y), length(z))
    @test xl isa GridUR
    @test yl isa GridUR
    @test zl isa GridUR
    @test xl == xa
    @test yl == ya
    @test zl == za
    @test xl[2,3,4] == xa[2,3,4] == x[2]
    @test yl[2,3,4] == ya[2,3,4] == y[3]
    @test zl[2,3,4] == za[2,3,4] == z[4]
end


text:
const LLVM_TYPES = IdDict{Type{<:NativeTypes},String}(
    Float32 => "float",
    Float64 => "double",
    Int8 => "i8",
    Int16 => "i16",
    Int32 => "i32",
    Int64 => "i64",
    UInt8 => "i8",
    UInt16 => "i16",
    UInt32 => "i32",
    UInt64 => "i64",
    Bool => "i8",
    Bit => "i1"
)
const LLVM_TYPES_SYM = IdDict{Symbol,String}(
    :Float32 => "float",
    :Float64 => "double",
    :Int8 => "i8",
    :Int16 => "i16",
    :Int32 => "i32",
    :Int64 => "i64",
    :UInt8 => "i8",
    :UInt16 => "i16",
    :UInt32 => "i32",
    :UInt64 => "i64",
    :Bool => "i8",
    :Bit => "i1",
    :Nothing => "void"
)

const JULIAPOINTERTYPE = 'i' * string(8sizeof(Int))

vtype(W, typ) = isone(abs(W)) ? typ : "<$W x $typ>"
vtype(W, T::DataType) = vtype(W, LLVM_TYPES[T])
julia_type(W, T) = isone(abs(W)) ? T : _Vec{W,T}

suffix(W::Int, s::String) = W == -1 ? s : 'v' * string(W) * s
suffix(W::Int, T) = suffix(W, suffix(T))
suffix(::Type{Ptr{T}}) where {T} = "p0" * suffix(T)
function suffix(@nospecialize(T))
    if T === Float32 || T === Float64
        t = 'f'
    elseif T === Bool
        return "i1"
    else
        t = 'i'
    end
    string(t, 8sizeof(T))
end

# Type-dependent LLVM constants
function llvmconst(T, val)
    iszero(val) && return "zeroinitializer"
    "$(LLVM_TYPES[T]) $val"
end
function llvmconst(::Type{Bool}, val)
    Bool(val) ? "i1 1" : "zeroinitializer"
end
function llvmconst(W::Integer, T, val)
    isa(val, Number) && iszero(val) && return "zeroinitializer"
    typ = LLVM_TYPES[T]
    '<' * join(("$typ $(val)" for _ in Base.OneTo(W)), ", ") * '>'
end
function llvmconst(W::Integer, ::Type{Bool}, val)
    Bool(val) || return "zeroinitializer"
    typ = "i1"
    '<' * join(("$typ $(Int(val))" for _ in Base.OneTo(W)), ", ") * '>'
end
function llvmconst(W::Int, v::String)
    '<' * join((v for _ in Base.OneTo(W)), ", ") * '>'
end
# function llvmtypedconst(T, val)
#     typ = LLVM_TYPES[T]
#     iszero(val) && return "$typ zeroinitializer"
#     "$typ $val"
# end
# function llvmtypedconst(::Type{Bool}, val)
#     Bool(val) ? "i1 1" : "i1 zeroinitializer"
# end

function llvmname(op, WR, WA, T, TA)
    lret = LLVM_TYPES[T]
    ln = "llvm.$op.$(suffix(WR,T))"
    (isone(abs(WR)) || T !== first(TA)) ? ln * '.' * suffix(maximum(WA),first(TA)) : ln
end

function llvmcall_expr(op, WR, R, WA, TA, ::Nothing = nothing)
    ff = llvmname(op, WR, WA, R, TA)
    argt = Expr(:tuple)
    foreach(WT -> push!(argt.args, julia_type(WT[1], WT[2])), zip(WA,TA))
    call = Expr(:call, :ccall, ff, :llvmcall, julia_type(WR, R), argt)
    foreach(n -> push!(call.args, Expr(:call, :data, Symbol(:v, n))), 1:length(TA))
    Expr(:block, Expr(:meta, :inline), isone(abs(WR)) ? call : Expr(:call, :Vec, call))
end
function llvmcall_expr(op, WR, R, WA, TA, flags::String)
    lret = LLVM_TYPES[R]
    lvret = vtype(WR, lret)
    lop = llvmname(op, WR, WA, R, TA)
    instr = "$lvret $flags @$lop"
    larg_types = vtype.(WA, TA)
    decl = "declare $lvret @$(lop)(" * join(larg_types, ", ") * ')'
    args_for_call = ("$T %$(n-1)" for (n,T) ∈ enumerate(larg_types))
    instrs = """%res = call $flags $lvret @$(lop)($(join(args_for_call, ", ")))
        ret $lvret %res"""
    args = Expr(:curly, :Tuple); foreach(WT -> push!(args.args, julia_type(WT[1], WT[2])), zip(WA,TA))
    arg_syms = [Expr(:call, :data, Symbol(:v,n)) for n ∈ 1:length(TA)]
    # println(instrs)
    llvmcall_expr(decl, instrs, julia_type(WR, R), args, lvret, larg_types, arg_syms)
end

@static if VERSION ≥ v"1.6.0-DEV.674"
    function llvmcall_expr(decl::String, instr::String, @nospecialize(ret), @nospecialize(args), lret::String, largs, arg_syms, callonly::Bool = false)
        mod = """
            $decl

            define $lret @entry($(join(largs, ", "))) alwaysinline {
            top:
                $instr
            }
        """
        # attributes #0 = { alwaysinline }
        call = Expr(:call, :llvmcall, (mod, "entry"), ret, args)
        foreach(arg -> push!(call.args, arg), arg_syms)
        if first(lret) === '<'
            call = Expr(:call, :Vec, call)
        end
        callonly && return call
        Expr(:block, Expr(:meta, :inline), call)
    end
else
    function llvmcall_expr(decl::String, instr::String, @nospecialize(ret), @nospecialize(args), lret::String, largs, arg_syms, callonly::Bool = false)
        call = Expr(:call, :llvmcall, (decl, instr), ret, args)
        foreach(arg -> push!(call.args, arg), arg_syms)
        if first(lret) === '<'
            call = Expr(:call, :Vec, call)
        end
        callonly && return call
        Expr(:block, Expr(:meta, :inline), call)
    end
end


text:
module LinkedListTest

using Test
using Tower.Nodes

@testset "Create single Nodes with the correct values" begin

    node1 = Node{Int64}(1)
    @test value(node1) == 1
    @test next(node1) == nothing

    node2 = Node{Int128}(Int128(10))
    @test value(node2) == Int128(10)
    @test next(node2) == nothing

    node3 = Node{Vector{Int64}}([1, 2, 3])
    @test value(node3) == [1, 2, 3]
    @test next(node3) == nothing

    node4 = Node{Tuple{Int64,Int64}}((11, 12))
    @test value(node4) == (11, 12)
    @test next(node4) == nothing
end

@testset "Create nodes with the next node and the correct value" begin
    node3 = Node{Int64}(30)
    @test value(node3) == 30
    @test next(node3) == nothing

    node2 = Node{Int64}(20, node3)
    @test value(node2) == 20
    @test next(node2) == node3

    node1 = Node{Int64}(10, node2)
    @test value(node1) == 10
    @test next(node1) == node2

    println(node1)
end

@testset "function value returns the correct value" begin
    node1 = Node{Float64}(12.3)
    @test value(node1) == 12.3

    node2 = Node{Int128}(Int128(123))
    @test value(node2) == Int128(123)
end

end


text:
module EFT
export twoSum, twoDiff, twoProd, twoDiv

# * Documentation of functions

@doc raw"""
    twoSum(a, b)

Error-Free Transformation for the sum of two numbers:

(x,y) = twoSum(a,b)

  ``\Leftrightarrow``

``x = fl(a+b)``  and  ``x+y = a+b``
"""
function twoSum end

@doc raw"""
    twoProd(a, b)

Error-Free Transformation for the product of two numbers:

(x,y) = twoProd(a, b)

  ``\Leftrightarrow``

``x = fl(a\,b)``  and  ``x+y = a\,b``
"""
function twoProd end

@doc raw"""
    twoDiv(a, b)

Approximate transformation for the division. The transformation is not exact
(the error is not representable) but at least the sign of y should be correct:

(x, y) = twoDiv(a, b)

  ``\Leftrightarrow``

``x = fl(a/b)``  and  ``y ≈ a/b - x``
"""
function twoDiv end

# * Generic implementation

# Generic implementation of twoSum. Uses the algorithm by D. Knuth.
function twoSum(a, b)
    x = a + b
    z = x - a
    y = (a - (x-z)) + (b-z)
    return (x,y)
end

twoDiff(a, b) = twoSum(a,-b)

@doc raw"""
    split(a)

Error-free splitting of a FP number:

(x,y) = split(a)

  ``\Leftrightarrow``

``x+y = a``, and the representations of x and y do not overlap

This implementation uses the algorithm by Dekker.
"""
function split(a)
    # New methods should be added to compute the correct
    # constant for each FP representation
    #
    # factor = 1+2^s, where s is the mantissa length
    factor(x::Float64) = 134217729
    factor(x::Float32) = 4097

    c = factor(a) * a
    x = c - (c - a)
    y = a - x
    return (x, y)
end

# Generic implementation of TwoProd. Uses the algorithm by Veltkamp.
function twoProd(a, b)
    x = a * b
    (a1, a2) = split(a)
    (b1, b2) = split(b)
    y = a2*b2 - (((x - a1*b1) - a2*b1) - a1*b2)
    return (x,y)
end

function twoDiv(a, b)
    x = a / b
    (p1, p2) = twoProd(x, b)
    y = (a-p1-p2)/b
    return (x, y)
end

# * Specific algorithms for Float32

# Since double precision (Float64) is available to run higher-precision
# computations, it should be used to compute more efficiently transformations of
# Float32 operations.

macro float32_transformation(operator, name)
    quote
        function $(esc(name))(a::Float32, b::Float32)
            x = $operator(a, b)
            r = $operator(Float64(a), Float64(b)) - Float64(x)
            y = Float32(r)
            return (x,y)
        end
    end
end

@float32_transformation(+, twoSum)
@float32_transformation(*, twoProd)
@float32_transformation(/, twoDiv)

end


text:
using Test, LinearMaps, LinearAlgebra, SparseArrays, InteractiveUtils
using LinearMaps: FiveArg

@testset "block maps" begin
    @testset "hcat" begin
        m = 3
        n = 4
        for elty in (Float32, ComplexF64), n2 in (0, 2)
            A11 = rand(elty, m, n)
            A12 = rand(elty, m, n2)
            a = rand(elty, m)
            L = @inferred hcat(LinearMap(A11), LinearMap(A12))
            @test L.maps isa Tuple
            Lv = @inferred LinearMaps.BlockMap{elty}([LinearMap(A11), LinearMap(A12)], (2,))
            @test Lv.maps isa Vector
            @test L == Lv == LinearMaps.BlockMap([LinearMap(A11), LinearMap(A12)], (2,))
            @test occursin("$m×$(n+n2) LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L))
            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
            @test L isa LinearMaps.BlockMap{elty}
            if elty <: Complex
                @test_throws ErrorException LinearMaps.BlockMap{Float64}((LinearMap(A11), LinearMap(A12)), (2,))
            end
            A = [A11 A12]
            x = rand(n+n2)
            @test size(L) == size(A) == size(Lv)
            @test Matrix(L) == A == Matrix(Lv) == mul!(copy(A), L, 1, true, false)
            @test L * x ≈ A * x ≈ Lv * x
            L = @inferred hcat(LinearMap(A11), LinearMap(A12), LinearMap(A11))
            A = [A11 A12 A11]
            @test Matrix(L) == A == mul!(zero(A), L, 1, true, false)
            A = [I I I A11 A11 A11 a]
            @test (@which [A11 A11 A11]).module != LinearMaps
            @test (@which [I I I A11 A11 A11]).module != LinearMaps
            @test (@which hcat(I, I, I)).module != LinearMaps
            @test (@which hcat(I, I, I, LinearMap(A11), A11, A11)).module == LinearMaps
            maps = @inferred LinearMaps.promote_to_lmaps(ntuple(i->m, 7), 1, 1, I, I, I, LinearMap(A11), A11, A11, a)
            @inferred LinearMaps.rowcolranges(maps, (7,))
            L = @inferred hcat(I, I, I, LinearMap(A11), A11, A11, a)
            @test L == [I I I LinearMap(A11) LinearMap(A11) LinearMap(A11) LinearMap(a)]
            x = ones(elty, size(L, 2))
            @test L isa LinearMaps.BlockMap{elty}
            @test L * x ≈ A * x
            L = @inferred hcat(I, I, I, LinearMap(A11), A11, A11, a, a, a, a)
            @test occursin("$m×$(3m+3n+4) LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L))
            L = @inferred hcat(I, I, I, LinearMap(A11), A11, A11, a, a, a, a, a, a, a)
            @test occursin("$m×$(3m+3n+7) LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L))
        end
        A11 = zeros(m+1, n)
        A12 = zeros(m, n)
        @test_throws DimensionMismatch hcat(LinearMap(A11), LinearMap(A12))
    end

    @testset "vcat" begin
        m = 2
        n = 3
        for elty in (Float32, ComplexF64)
            A11 = rand(elty, m, n)
            v = rand(elty, n)
            L = @inferred vcat(LinearMap(A11))
            @test L == [LinearMap(A11);]
            @test Matrix(L) ≈ A11
            A21 = rand(elty, 2m, n)
            L = @inferred vcat(LinearMap(A11), LinearMap(A21))
            @test L.maps isa Tuple
            @test L isa LinearMaps.BlockMap{elty}
            @test occursin("$(3m)×$n LinearMaps.BlockMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), L))
            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
            Lv = LinearMaps.BlockMap{elty}([LinearMap(A11), LinearMap(A21)], (1,1))
            @test Lv.maps isa Vector
            @test L == Lv
            @test (@which [A11; A21]).module != LinearMaps
            A = [A11; A21]
            x = rand(elty, n)
            @test size(L) == size(A)
            @test Matrix(L) == A == mul!(copy(A), L, 1, true, false)
            @test Matrix(Lv) == A == mul!(copy(A), Lv, 1, true, false)
            @test L * x ≈ Lv * x ≈ A * x
            A = [I; I; I; A11; A11; A11; reduce(hcat, fill(v, n))]
            @test (@which [I; I; I; A11; A11; A11; v v v v v v v v v v]).module != LinearMaps
            L = @inferred vcat(I, I, I, LinearMap(A11), LinearMap(A11), LinearMap(A11), reduce(hcat, fill(v, n)))
            @test L == [I; I; I; LinearMap(A11); LinearMap(A11); LinearMap(A11); reduce(hcat, fill(v, n))]
            @test L isa LinearMaps.BlockMap{elty}
            @test L * x ≈ A * x
        end
        A11 = zeros(m, n+1)
        A21 = zeros(2m, n)
        @test_throws DimensionMismatch vcat(LinearMap(A11), LinearMap(A21))
    end

    @testset "hvcat" begin
        m1 = 2
        m2 = 3
        n = 3
        for elty in (Float32, ComplexF64)
            A11 = rand(elty, m1, m1)
            A12 = ones(elty, m1, m2)
            A21 = rand(elty, m2, m1)
            A22 = ones(elty, m2, m2)
            A = [A11 A12; A21 A22]
            @test (@which [A11 A12; A21 A22]).module != LinearMaps
            @inferred hvcat((2,2), LinearMap(A11), LinearMap(A12), LinearMap(A21), LinearMap(A22))
            L = [LinearMap(A11) LinearMap(A12); LinearMap(A21) LinearMap(A22)]
            @test L.maps isa Tuple
            Lv = @inferred LinearMaps.BlockMap{elty}([LinearMap(A11), LinearMap(A12), LinearMap(A21), LinearMap(A22)], (2,2))
            @test Lv.maps isa Vector
            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
            @test @inferred !issymmetric(L)
            @test @inferred !ishermitian(L)
            x = rand(m1+m2)
            @test L isa LinearMaps.BlockMap{elty}
            @test size(L) == size(A)
            @test L * x ≈ Lv * x ≈ A * x
            @test Matrix(L) == Matrix(Lv) == A == mul!(zero(A), L, 1)
            for α in (false, 1, rand()), β in (false, 1, rand())
                @test mul!(copy(A), L, 2, α, β) ≈ A*(2α + β)
            end
            @test convert(AbstractMatrix, L) == A
            A = [I A12; A21 I]
            @test (@which [I A12; A21 I]).module != LinearMaps
            @inferred hvcat((2,2), I, LinearMap(A12), LinearMap(A21), I)
            L = @inferred hvcat((2,2), I, LinearMap(A12), LinearMap(A21), I)
            @test L isa LinearMaps.BlockMap{elty}
            @test size(L) == (m1+m2, m1+m2)
            @test Matrix(L) ≈ A
            @test L * x ≈ A * x
            y = randn(elty, size(L, 1))
            for α in (0, 1, rand(elty)), β in (0, 1, rand(elty))
                @test mul!(copy(y), L, x, α, β) ≈ y*β .+ A*x*α
            end
            X = rand(elty, m1+m2, n)
            Y = randn(elty, size(L, 1), n)
            for α in (0, 1, rand(elty)), β in (0, 1, rand(elty))
                @test mul!(copy(Y), L, X, α, β) ≈ Y*β .+ A*X*α
            end
            A = ones(elty, m1, m1); LA = LinearMap(A)
            B = zeros(elty, m2, 3m1); LB = LinearMap(B)
            @test [LA LA LA; LB] isa LinearMaps.BlockMap{elty}
            @test Matrix([LA LA LA; LB]) ≈ [A A A; B]
            @test [LB; LA LA LA] isa LinearMaps.BlockMap{elty}
            @test Matrix([LB; LA LA LA]) ≈ [B; A A A]
            @test [I; LA LA LA] isa LinearMaps.BlockMap{elty}
            @test Matrix([I; LA LA LA]) ≈ [I; A A A]
            A12 = LinearMap(zeros(elty, m1, m2+1))
            A21 = LinearMap(zeros(elty, m2, m1))
            @test_throws DimensionMismatch [I A12; A21 I]
            @test_throws DimensionMismatch [I A21; A12 I]
            @test_throws DimensionMismatch [A12 A12; A21 A21]
            @test_throws DimensionMismatch [A12 A21; A12 A21]

            # basic test of "misaligned" blocks
            M = ones(elty, 3, 2) # non-square
            A = LinearMap(M)
            B = [I A; A I]
            C = [I M; M I]
            @test B isa LinearMaps.BlockMap{elty}
            @test Matrix(B) == C
            @test Matrix(transpose(B)) == transpose(C)
            @test Matrix(adjoint(B)) == C'
        end
    end

    @testset "adjoint/transpose" begin
        for elty in (Float32, ComplexF64), transform in (transpose, adjoint)
            A12 = rand(elty, 10, 10)
            A = [I A12; transform(A12) I]
            L = [I LinearMap(A12); transform(LinearMap(A12)) I]
            @test @inferred(LinearMaps.MulStyle(L)) === FiveArg()
            if elty <: Complex
                if transform == transpose
                    @test @inferred issymmetric(L)
                else
                    @test @inferred ishermitian(L)
                end
            end
            if elty <: Real
                @test @inferred ishermitian(L)
                @test @inferred issymmetric(L)
            end
            x = rand(elty, 20)
            @test L isa LinearMaps.LinearMap{elty}
            @test size(L) == size(A)
            @test L * x ≈ A * x
            @test Matrix(L) == A
            @test convert(AbstractMatrix, L) == A
            @test sparse(L) == sparse(A)
            Lt = @inferred transform(L)
            @test Lt isa LinearMaps.LinearMap{elty}
            @test Lt * x ≈ transform(A) * x
            @test convert(AbstractMatrix, Lt) == transform(A)
            @test sparse(transform(L)) == transform(A)
            for α in (false, 1, rand()), β in (false, 1, rand())
                @test mul!(copy(transform(A)), Lt, 2, α, β) ≈ transform(A)*(2α + β)
            end
            Lt = @inferred transform(LinearMap(L))
            @test Lt * x ≈ transform(A) * x
            @test Matrix(Lt) == Matrix(transform(A))
            A21 = rand(elty, 10, 10)
            A = [I A12; A21 I]
            L = [I LinearMap(A12); LinearMap(A21) I]
            Lt = @inferred transform(L)
            @test Lt isa LinearMaps.LinearMap{elty}
            @test Lt * x ≈ transform(A) * x
            @test Matrix(Lt) ≈ Matrix(transform(LinearMap(L))) ≈ Matrix(transform(A))
            @test Matrix(transform(LinearMap(L))+transform(LinearMap(L))) ≈ 2Matrix(transform(A))
            X = rand(elty, size(L, 1), 10)
            Y = randn(elty, size(L, 2), 10)
            for α in (0, 1, rand(elty)), β in (0, 1, rand(elty))
                @test mul!(copy(Y), Lt, X, α, β) ≈ Y*β .+ transform(A)*X*α
            end
        end
    end

    @testset "block diagonal maps" begin
        for elty in (Float32, ComplexF64)
            m = 2; n = 3
            M1 = 10*(1:m) .+ (1:(n+1))'; L1 = LinearMap(M1)
            M2 = randn(elty, m, n+2); L2 = LinearMap(M2)
            M3 = randn(elty, m, n+3); L3 = LinearMap(M3)

            if elty <: Complex
                @test_throws ErrorException LinearMaps.BlockDiagonalMap{Float64}((L1, L2, L3, L2, L1))
            end
            # Md = diag(M1, M2, M3, M2, M1) # unsupported so use sparse:
            Md = Matrix(blockdiag(sparse.((M1, M2, M3, M2, M1))...))
            @test (@which blockdiag(sparse.((M1, M2, M3, M2, M1))...)).module != LinearMaps
            @test (@which cat(M1, M2, M3, M2, M1; dims=(1,2))).module != LinearMaps
            x = randn(elty, size(Md, 2))
            Bd = @inferred blockdiag(L1, L2, L3, L2, L1)
            @test Bd.maps isa Tuple
            Bdv = @inferred LinearMaps.BlockDiagonalMap{elty}([L1, L2, L3, L2, L1])
            @test Bdv.maps isa Vector
            @test @inferred(LinearMaps.MulStyle(Bd)) === FiveArg()
            @test occursin("$(5m)×$(5n+9) LinearMaps.BlockDiagonalMap{$elty}", sprint((t, s) -> show(t, "text/plain", s), Bd))
            @test Matrix(Bd) == Md == mul!(zero(Md), Bd, 1)
            for α in (false, rand()), β in (false, rand())
                @test mul!(copy(Md), Bd, 2, α, β) ≈ Md*(2α + β)
            end
            @test convert(AbstractMatrix, Bd) isa SparseMatrixCSC
            @test sparse(Bd) == Md
            @test Matrix(@inferred blockdiag(L1)) == M1
            @test Matrix(@inferred blockdiag(L1, L2)) == blockdiag(sparse.((M1, M2))...)
            Bd2 = @inferred cat(L1, L2, L3, L2, L1; dims=(1,2))
            @test_throws ArgumentError cat(L1, L2, L3, L2, L1; dims=(2,2))
            @test Bd == Bdv == Bd2
            @test Bd == blockdiag(L1, M2, M3, M2, M1)
            @test size(Bd) == (5m, 5n+9)
            @test !issymmetric(Bd)
            @test !ishermitian(Bd)
            @test (@inferred Bd * x) ≈ Bdv * x ≈ Md * x
            for transform in (identity, adjoint, transpose)
                @test Matrix(@inferred transform(Bd)) == transform(Md)
                @test Matrix(@inferred transform(LinearMap(Bd))) == transform(Md)
            end
            y = randn(elty, size(Md, 1))
            for α in (0, 1, rand(elty)), β in (0, 1, rand(elty))
                @test mul!(copy(y), Bd, x, α, β) ≈ y*β .+ Md*x*α
                @test mul!(copy(y), Bdv, x, α, β) ≈ y*β .+ Md*x*α
            end
            X = randn(elty, size(Md, 2), 3)
            Y = randn(elty, size(Md, 1), 3)
            for α in (0, 1, rand(elty)), β in (0, 1, rand(elty))
                @test mul!(copy(Y), Bd, X, α, β) ≈ Y*β .+ Md*X*α
                @test mul!(copy(Y), Bdv, X, α, β) ≈ Y*β .+ Md*X*α
            end
        end
    end

    @testset "function block map" begin
        N = 5
        T = ComplexF64
        CS! = LinearMap{T}(cumsum!,
            (y, x) -> (copyto!(y, x); reverse!(cumsum!(y, reverse!(y)))), N;
            ismutating=true)
        A = rand(T, N, N)
        B = rand(T, N, N)
        LT = LowerTriangular(ones(T, N, N))
        L1 = [CS! CS! CS!; CS! CS! CS!; CS! CS! CS!]
        M1 = [LT LT LT; LT LT LT; LT LT LT]
        L2 = [CS! LinearMap(A) CS!; LinearMap(B) CS! CS!; CS! CS! CS!]
        M2 = [LT A LT
              B LT LT
              LT LT LT]
        u = rand(T, 3N)
        v = rand(T, 3N)
        for α in (false, true, rand(T)), β in (false, true, rand(T))
            for transform in (identity, adjoint), (L, M) in ((L1, M1), (L2, M2))
                # @show α, β, transform
                @test mul!(copy(v), transform(L), u, α, β) ≈ transform(M)*u*α + v*β
                @test mul!(copy(v), transform(LinearMap(L)), u, α, β) ≈ transform(M)*u*α + v*β
                @test mul!(copy(v), LinearMap(transform(L)), u, α, β) ≈ transform(M)*u*α + v*β
                if transform != adjoint
                    transL = transform(L)
                    alloc = @allocated similar(v)
                    if L == L2 && α != false
                        @test_broken (@allocated mul!(v, transL, u, α, β)) <= alloc
                    else
                        @test (@allocated mul!(v, transL, u, α, β)) <= alloc
                    end
                end
            end
        end
    end
end


text:
# utilities for overapprox_nd_relational.jl and overest_nd.jl
EVAL_CONSTS = true # some solvers can handle expressions like "log(2)" or "1/6" but others can only handle 0.16666*x
# When EVAL_CONSTS is true, it will evaluate an expression like "1/6" and turn it into 0.16666
# TODO: add this functionality to be turned on or off for expressions like log(2)

# list of operators and fucntions that are not counted as variables
special_oper = [:+, :-, :/, :*, :^]
special_func = [:exp, :log, :log10,
                :sin, :cos, :tan,
                :sinh, :cosh, :tanh,
                :asin, :acos, :atan,
                :asinh, :acosh, :atanh]

increasing_special_func = [:exp, :log, :log10,
                         :tan, :sinh, :tanh,
                         :asin, :atan,
                         :asinh, :atanh, :acosh]
special_consts = [π]

N_VARS = 0 # number of variables; has to be defined globally.
@debug("N_VARS := 0")

add_var(bound) = add_var()

function add_var()
    # bound.nvars += 1
    # @ is the symbol preceding A in ascii
    #return Symbol('v'*('@'+bound.nvars))

    # the counter has to be a global variable.
    global N_VARS
    N_VARS += 1
    return Symbol("v_$N_VARS")
end

function reset_NVARS()
    global N_VARS
    N_VARS = 0
end

function to_pairs(B)
    """
    This function converts the output of overest, a tuple of (x points, y points) in the form ready
        for closed form generator: [(x1,y1), (x2, y2), ...]
    B is an array of points.
    """
    xp,yp = B
    x = vcat(xp...)
    y = vcat(yp...)
    pairs = collect(zip(x,y))
    return pairs
end

# todo: does this do the same thing as SymEngine, free_symbols or Tomer's autoline.get_symbols ?
function find_variables(expr::Expr)
    """
    given an expression expr, this function finds all the variables.
    it is useful to identify 1d functions that can be directly implemented
    in the overest algorithm, without any composition hacking.
    Example: find_variables(:(x+1))     = [:x]
             find_variables(:(-x+2y-z)) = [:x, :y, :z]
             find_variables(:(log(x)))  = [:x]
             find_variables(:(x+ x*z))   = [:x, :z]
    """
    all_vars = []
    for arg in expr.args
        if arg isa Expr
            all_vars = vcat(all_vars, find_variables(arg))
        elseif arg isa Symbol
            if !(arg in special_oper) && !(arg in special_func) && !(arg in special_consts)
                all_vars = vcat(all_vars, arg)
            end
        end
    end
    # # use symengine function instead
    # all_vars = Symbol.(free_symbols(Basic(expr))) # this bungles the case of 0*x because Basic makes it just 0 and then the expression looks like it has no free vars
    return unique(all_vars)
end

function is_affine(expr)
    """
    given an expression expr, this function determines if the expression
    is an affine function.

    Example: is_affine(:(x+1))       = true
             is_affine(:(-x+(2y-z))) = true
             is_affine(:(log(x)))    = false
             is_affine(:(x + x*z))   = false
             is_affine(:(x/6))       = true
             is_affine(:(5*x))       = true
             is_affine(:(log(2)*x))  = true
             is_affine(:(-x))        = true
     """
    # it is number
    if is_number(expr)
        return true
    elseif expr isa Symbol # symbol
        return true
    elseif expr isa Expr
        check_expr_args_length(expr)
        func = expr.args[1]::Symbol
        if func ∉ [:+, :-, :*, :/] # only these operations are allowed
            return false
        else  # func ∈ [:+, :-, :*, :/]
            if func == :* # one of args has to be a number
                option1 =  is_number(expr.args[2]) && is_affine(expr.args[3])
                option2 =  is_number(expr.args[3]) && is_affine(expr.args[2])
                return (option1 || option2)
            elseif func == :/ # second arg has to be a number
                return is_affine(expr.args[2]) && is_number(expr.args[3])
            else # func is + or -
                return all(is_affine.(expr.args[2:end]))
            end
        end
    else
        @error "This case should not be reached. Expression of the wrong type was passed to is_affine()."
        return false # if not a number, symbol, or Expr, return false
    end
end

is_outer_affine(s::Symbol) = true
is_outer_affine(r::Real) = true
function is_outer_affine(expr::Expr)
    """
    5*sin(x) - 3*cos(y) is _outer_ affine. It can be re-written:
    5*z - 3*w    where z = sin(x)   w = cos(y)
    """
    if is_number(expr)
        return true
    else
        check_expr_args_length(expr)
        func = expr.args[1]
        if func ∈ [:+, :-]
            # check args
            return all(is_outer_affine.(expr.args[2:end]))
        elseif func == :* # one of args has to be a number
            option1 =  is_number(expr.args[2]) && is_outer_affine(expr.args[3])
            option2 =  is_number(expr.args[3]) && is_outer_affine(expr.args[2])
            return (option1 || option2)
        elseif func == :/ # second arg has to be a number
            return is_number(expr.args[3])
        else
            return false
        end
    end
end

function add_ϵ(points, ϵ::T where T <: Real)
    `Add ϵ to the y values of all points in a container`
    @debug "Added $ϵ to all points"
    new_points = Tuple{Float64, Float64}[]
    for p in points
        push!(new_points, (p[1], p[2] + ϵ))
    end
    return new_points
end

function rewrite_division_by_const(e)
    return e
end
function rewrite_division_by_const(expr::Expr)
    if expr.args[1] == :/ && !is_number(expr.args[2]) && is_number(expr.args[3])
        if EVAL_CONSTS
            c = eval(1/expr.args[3])
            return :($c * $(expr.args[2]))
        else
            return :( (1/$(expr.args[3])) * $(expr.args[2]) )
        end
    else
        return expr
    end
end

function get_sincos_regions(a,b; offset=0)
    """
    Return inflection points for sin, cos
    """

    n̂_a = ceil((a - offset) / π)
    n̂_b = floor((b - offset) / π)
    n_array = [i for i in n̂_a:n̂_b]

    if length(n_array) == 0 # no inflection points within interval
        @assert sin(a + offset) != 0.0
        return n_array, sin(a + offset) < 0 # true if convex, false implies concave
    else # greater than 0 length
        return [offset + n*π for n in n_array], nothing # nothing denotes mixed convexity
    end
end

function get_regions_unary(func::Symbol, a, b)
    """
    Return the zeros of the second derivative of function func and/or whether it is convex or not

    d2f_zeros, convex = get_regions(func, a, b)
    """
    if func == :cos
        d2f_zeros, convex = get_sincos_regions(a, b, offset=π/2)
    elseif func == :sin
        d2f_zeros, convex = get_sincos_regions(a, b, offset=0)
    elseif func == :exp
        d2f_zeros, convex = [], true
    elseif func == :log
        @assert a > 0
        d2f_zeros, convex = [], false
    elseif func == :tanh
        if b ≤ 0
            d2f_zeros, convex = [], true
        elseif a ≥ 0
            d2f_zeros, convex = [], false
        else # spans 0
            d2f_zeros, convex = [0], nothing
        end
    else
        d2f_zeros, convex = nothing, nothing
    end

    return d2f_zeros, convex
end

function division_d2f_regions(e, arg, a, b)
    # function defined over the interval [a,b]
    c = e.args[2]
    @assert is_number(c)
    @assert (a <= b)
    if eval(c) > 0
        if a > 0
            d2f_zeros, convex = [], true
        elseif b < 0
            d2f_zeros, convex = [], false
        else
            error("ERROR: interval ["*string(a)*","*string(b)*"] straddles zero for function c/x.The function is discontinuous and cannot be bounded.")
        end
    elseif eval(c) < 0
        if a > 0
            d2f_zeros, convex = [], false
        elseif b < 0
            d2f_zeros, convex = [], true
        else
            error("ERROR: interval ["*string(a)*","*string(b)*"] straddles zero for function c/x.The function is discontinuous and cannot be bounded.")
        end
    end
    return d2f_zeros, convex
end

function exponent_d2f_regions(e, arg, a, b)
    if is_number(e.args[2]) # c^x
        @assert eval(e.args[2]) > 0 # only real valued over real arguments for c > 0
        d2f_zeros, convex = [], true
    elseif is_number(e.args[3]) # x^c (polynomials)
        x = e.args[2]
        c = e.args[3]
        # a few cases
        # 1) c is fractional. only valid for x >= 0. check. and convex (increasing) no inflection points
        if (c % 1) != 0
            @assert a >= 0 # ensures whole interval is > 0
            if (c > 0) && (c < 1)
                d2f_zeros, convex = [], false
            else
                d2f_zeros, convex = [], true
            end
        # 2) c is odd (3, 5, 7): inflection point at zero may be applicable. convex for x>0, concave for x<0
        elseif (c % 2) == 1
            if a > 0
                d2f_zeros, convex = [], true
            elseif b < 0
                d2f_zeros, convex = [], false
            else # interval overlaps 0
                d2f_zeros, convex = [0], nothing
            end
        # 3) c is even (2, 4, 6): convex, no inflection points
        elseif (c % 2) == 0
            if c < 0
                # discontinuous at zero for negative even numbers
                @assert ((a > 0) || (b < 0)) && (a ≤ b)
            end
            d2f_zeros, convex = [], true
        end
    else
        # TODO: handle x^y
        error("Expression type not handled.")
    end
    return d2f_zeros, convex
end

function get_regions_1arg(e::Expr, arg::Symbol, a, b)
    # TODO:
    # check if c/x or c^x or x^c
    # multiplication between two variables is expanded using log and exp
    # TODO: add case to catch x^1 which is affine
    func = e.args[1]
    if func == :/
        d2f_zeros, convex = division_d2f_regions(e, arg, a, b)
    elseif func == :^
        d2f_zeros, convex = exponent_d2f_regions(e, arg, a, b)
    else
        @debug "Expression not supported."
        d2f_zeros, convex = nothing, nothing
    end
    @debug "Expression $e is convex: $convex and has inflection points: $d2f_zeros"
    return d2f_zeros, convex
end

function find_UB(func, a, b, N; lb=false, plot=false, ϵ=0.0, d2f_zeros=nothing, convex=nothing)

    """
    This function finds the piecewise linear upperbound (lowerbound) of
        a given function func in internval [a,b] with N
        sampling points over each concave/convex region.
        see the overest_new.jl for more details.

    Return values are points (UB_points) and the min-max closed form (UB_sym)
    """
    UB = bound(func, a, b, N; lowerbound=lb, d2f_zeros=d2f_zeros, convex=convex, plot=plot)
    UB_points = unique(sort(to_pairs(UB), by = x -> x[1]))
    # if !lb
    #     println("Is upper bound.")
    # else
    #     println("Is lower bound.")
    # end
    #println("Max of points: ", maximum(UB_points))
    if abs(ϵ) > 0
        # treat epsilon as relative percentage e.g. 0.01 is 1%
        # xpts, ypts = UB
        ypts = vcat(UB[2]...)
        pt_range = maximum(ypts) - minimum(ypts)
        rel_ϵ = pt_range*ϵ
        @debug("Relative ϵ percent is $(ϵ*100)% which for function in range [$(minimum(ypts)), $(maximum(ypts))] is $rel_ϵ")
        UB_points = add_ϵ(UB_points, rel_ϵ) # return new points shifted by relative epsilon up or down
    end
    #println("Max of points after ϵ adjustment: ", maximum(UB_points))
    UB_sym = closed_form_piecewise_linear(UB_points)
    #UB_eval = SymEngine.lambdify(:(x -> $UB_sym), [:x])
    return UB_points::Vector{Tuple{Float64, Float64}}, UB_sym #, UB_eval
end

function find_1d_range(B)
    """
    given a set of points B representing a piecewise linear bound function
        this function finds the range of the bound.
    """

    y_pnts = [point[2] for point in B]
    min_B = minimum(y_pnts)
    max_B = maximum(y_pnts)
    return min_B, max_B
end

function check_expr_args_length(expr)
    if length(expr.args) > 3
        throw(ArgumentError("""
        Operation $(expr.args[1]) has $(length(expr.args)-1) arguments.
        Use parantheses to make this into multiple operations, each with two arguments,
        For example, change (x+y+z)^2 to ((x+y)+z)^2.
                          """))
    end
end

# Think we could use stuff from https://github.com/JuliaIntervals/IntervalArithmetic.jl
# but whose to say if it's better tested?
function find_affine_range(expr, range_dict)
    """
    given a an affine expression expr, this function finds the
        lower and upper bounds of the range.
    """
    @assert is_affine(expr)
    @debug "$expr is affine"

    @debug "testing if $expr can be eval'ed to a number"
    try (return eval(expr), eval(expr)) catch; nothing end # works if expr is a number
    @debug "testing if $expr is a symbol"
    expr isa Symbol ? (return range_dict[expr]) : nothing

    check_expr_args_length(expr)
    all_vars = find_variables(expr)

    if length(all_vars) == 1
        a, b = range_dict[all_vars[1]]
        expr_copy = copy(expr)
        c = eval(substitute!(expr_copy, all_vars[1], a))
        expr_copy = copy(expr)
        d = eval(substitute!(expr_copy, all_vars[1], b))
        return min(c,d), max(c,d)
    end

    func = expr.args[1]
    a, b = find_affine_range(expr.args[2], range_dict)
    c, d = find_affine_range(expr.args[3], range_dict)
    if func == :+
        return (a+c), (b+d)
    elseif func == :-
        return (a-d), (b-c)
    elseif func == :*
        # Does this assume anything about the order of the multiplication being const*variable ?
        try
            a = eval(expr.args[2])
            a > 0 ? (return a*c, a*d) : return a*d, a*c
        catch
            c = eval(expr.args[3])
            c > 0 ? (return c*a, c*b) : (return c*b, c*a)
        end
    elseif func == :/
        c = eval(expr.args[3])
        c > 0 ? (return a/c, b/c) : return( b/c, a/c)
    end
end

mutable struct MyException <: Exception
    var::String
end

function find_range(expr, range_dict)
    """
    Find range of PWL function.
    """
    if is_affine(expr)
        l,u = find_affine_range(expr, range_dict)
        return [l,u]
    elseif is_relu(expr)
        @debug expr " is relu"
        if expr.args[1] == :relu
            inner_expr = expr.args[2]
        elseif expr.args[1] == :max # max with 0
            if expr.args[2] == 0.0
                inner_expr = expr.args[3]
            elseif expr.args[3] == 0.0
                inner_expr = expr.args[2]
            end
        end
        l,u = find_range(inner_expr, range_dict)
        @debug "l = " l ", u = " u
        return [0, max(0, u)]
    elseif expr.args[1] == :min
        x = expr.args[2]
        y = expr.args[3]
        xl, xu = find_range(x, range_dict)
        yl, yu = find_range(y, range_dict)
        return [min(xl,yl), min(xu,yu)]
    elseif expr.args[1] == :max
        x = expr.args[2]
        y = expr.args[3]
        xl, xu = find_range(x, range_dict)
        yl, yu = find_range(y, range_dict)
        return [max(xl,yl), max(xu,yu)]
    else
        throw(MyException("not implemented yet"))
    end
end

function is_relu(expr)
    if length(expr.args) < 2
        return false
    end
    c1 = (expr.args[1] == :max && expr.args[2] == 0.)
    c2 = length(expr.args) > 2 ? (expr.args[1] == :max && expr.args[3] == 0.) : false
    c3 = expr.args[1] == :relu
    return c1 || c2 || c3
end

function is_min(expr)
    return expr.args[1] == :min
end

# todo: pretty sure we can use SymEngine.subs.
# maybe better tested? but subs is also overly complicated...
function substitute!(expr, old, new)

    """
    This function substitutes the old value `old' with the
         new value `new' in the expression expr

    Example: substitute!(:(x^2+1), :x, :(y+1)) = :((y+1)^2+1))
    """

    for (i,arg) in enumerate(expr.args)
       if arg==old
           expr.args[i] = new
       elseif arg isa Expr
           substitute!(arg, old, new)
       end
    end
    return expr
end

function substitute!(expr::Expr, old_list::Vector{Any}, new_list::Array{Any})
    for (k, v) in zip(old_list, new_list)
        substitute!(expr, k, v)
    end
    return expr
end

∉(e, set) = !(e ∈ set)

function reduce_args_to_2(f::Symbol, arguments::Array)
    """
    reduce_args_to_2(x+y) = x+y
    reduce_args_to_2(x+y+z) = x+(y+z)
    reduce_args_to_2(0.5+y+z) = (0.5+y)+z
    """
    if (length(arguments) <= 2) | (f ∉ [:+, :*])
        e = Expr(:call)
        e.args = [f, map(reduce_args_to_2, arguments)...]
        return e
    elseif length(arguments) ==3 #length(args) = 3 and f ∈ [:+, :*]
        if is_number(arguments[1])
            a = reduce_args_to_2(arguments[2])
            fbc = reduce_args_to_2(f, arguments[1:3 .!= 2])
        else
            a = reduce_args_to_2(arguments[1])
            fbc = reduce_args_to_2(f, arguments[2:3])
        end
        return :($f( $a, $fbc))
    else
        first_half = reduce_args_to_2(f, arguments[1:2])
        second_half = reduce_args_to_2(f, arguments[3:end])
        return reduce_args_to_2(f, [first_half, second_half])
        error("""
            operations with more than 3 arguments are not supported.
            use () to break the arguments.
            Example= x+y+z+w -> (x+y) + (y+z)
            """)
    end
end

reduce_args_to_2(x::Number) = x
reduce_args_to_2(x::Symbol) = x

function reduce_args_to_2(expr::Expr)
    #println(expr)
    f = expr.args[1]
    arguments = expr.args[2:end]
    return reduce_args_to_2(f::Symbol, arguments::Array)
end

function is_number(expr)
    try
        eval(expr)
        return true
    catch
        return false
    end
end

function is_unary(expr::Expr)
    # one for function, one for single argument to function
    return length(expr.args) == 2
end

function is_1d(expr::Expr)
    return length(find_variables(expr)) == 1
end

is_binary(expr::Expr) = length(expr.args) == 3

function multiply_interval(range, constant)
    S = [range[1]*constant, range[2]*constant]
    return [min(S...), max(S...)]
end


text:
__precompile__()

module ModifiedStokesSolver

using Compat
using Compat.Libdl
using Compat.LinearAlgebra
using Compat.MathConstants
using Compat.@info

using PyCall

using CurveDiscretization
using BarycentricLagrange

export DiscreteCurve
export discretize, multiply_connected
export doublelayer, doublelayer_fast, doublelayer_self, doublelayer_matrix
export system_matrix, system_matvec
export correction_block_matrix
export stokeslet
export pressure
export stresslet, stresslet_split
export dblkernel, dblkernel_log, dblkernel_block!
export fmm_stokeslet_direct, fmm_stokeslet_targ
export fmm_stresslet_direct, fmm_stresslet_prep
export fmm_stresslet_targ, fmm_stresslet_self
export ndgrid, interior_points
export fs_modstokes, per_modstokes, per_gradient

include("../../external/mbh2dfmm/julia/ModStokesFMM2D.jl")
include("types.jl")
include("helsing.jl")
include("layerpotential.jl")
include("integralequation.jl")
include("kernels.jl")
include("kernel_gradient.jl")
include("legendre.jl")
include("sparse.jl")
include("freespace.jl")
include("periodic.jl")
include("matlab_compat.jl")
include("largealpha.jl")

function __init__()
    # Load libgomp for FMM
    Libdl.dlopen("libgomp", Libdl.RTLD_GLOBAL)
end

end


text:
"""
ResNet Models

Pre trained ResNet{18,34,50,101,152} weights are available.
See below to see how to use it!

```julia
  julia> include(KnetLayers.dir("examples","resnet.jl"))

  julia> using .ResNetLib

  julia> m = ResNet{50}(trained=true)
  ResNet{50}()

  julia> topK(m.labels, m("./gray-wolf_sam-parks.png");K=5)
  5-element Array{String,1}:
  "timber wolf, grey wolf, gray wolf, Canis lupus"
  "white wolf, Arctic wolf, Canis lupus tundrarum"
  "red wolf, maned wolf, Canis rufus, Canis niger"
  "dingo, warrigal, warragal, Canis dingo"
  "coyote, prairie wolf, brush wolf, Canis latrans"
```
"""
module ResNetLib

using Statistics, Knet, KnetLayers, BSON, ImageMagick, Images

  struct ResNet{N}
    layers::Chain
    labels::Array{String}
  end

  struct BasicV2
    layers::Chain
    downsample::Union{Conv,Nothing}
  end

  function BasicV2(channels, stride; downsample=false, in_channels=0, kwargs...)
    layers = Chain(
                  BatchNorm(in_channels),
                  ReLU(),
                  Conv(height=3, width=3, inout=in_channels=>channels,
                       padding = 1, stride = stride, binit=nothing, mode=1),
                  BatchNorm(channels÷4),
                  ReLU(),
                  Conv(height=3, width=3, inout=channels=>channels,
                       padding = 1, stride = 1, binit=nothing, mode=1)
                  )
    if downsample
      return BasicV2(layers, Conv(height=1, width=1, inout=in_channels=>channels,
                                  padding = 0, stride = stride, mode=1, binit=nothing))
    else
      return BasicV2(layers,nothing)
    end
  end

  function (m::BasicV2)(x)
        residual = x
        x = m.layers[1:2](x)
        if m.downsample !== nothing
            residual = m.downsample(x)
        end
        x = m.layers[3:length(m.layers.layers)](x)
        return x + residual
  end

  struct BottleneckV2
    layers::Chain
    downsample::Union{Conv,Nothing}
  end

  function BottleneckV2(channels, stride; downsample=false, in_channels=0, kwargs...)
    layers = Chain(
                  BatchNorm(in_channels),
                  ReLU(),
                  Conv(height=1, width=1, inout=in_channels=>channels÷4,
                       padding = 0, stride = 1, binit=nothing, mode=1),
                  BatchNorm(channels÷4),
                  ReLU(),
                  Conv(height=3, width=3, inout=channels÷4=>channels÷4,
                       padding = 1, stride = stride, binit=nothing, mode=1),
                  BatchNorm(channels÷4),
                  ReLU(),
                  Conv(height=1, width=1, inout=channels÷4=>channels,
                       padding = 0, stride = 1, binit=nothing, mode=1)
                  )
    if downsample
      return BottleneckV2(layers,
                          Conv(height=1, width=1, inout=in_channels=>channels,
                          padding = 0, stride = stride, mode=1, binit=nothing))
    else
      return BottleneckV2(layers,nothing)
    end
  end

  function (m::BottleneckV2)(x)
        residual = x
        x = m.layers[1:2](x)
        if m.downsample !== nothing
            residual = m.downsample(x)
        end
        x = m.layers[3:length(m.layers.layers)](x)
        return x + residual
  end

  function _make_layer(block, layers, channels, stride, stage_index; in_channels=0)
    layer = [block(channels, stride, downsample=(channels != in_channels), in_channels=in_channels)]
    for _ in 1:layers-1
           push!(layer, block(channels, 1, downsample=false, in_channels=channels))
    end
    return Chain(layer...)
  end

  function _init(block, layers, channels; classes=1000, N=50, stage=0)
    @assert length(layers) == length(channels) - 1 "error"
    top = Chain(BatchNorm(3),
                  Conv(height=7, width=7, inout=3=>channels[1],
                       padding = 1, stride = 2, binit=nothing, mode=1),
                  BatchNorm(channels[1]), ReLU(),
                  Pool(window=3, padding = 1, stride = 2))

      stages = Chain[]
      in_channels = channels[1]
      for (i, num_layer) in enumerate(layers)
            stride = i == 1 ? 1 : 2
            push!(stages,_make_layer(block, num_layer, channels[i+1], stride, i+1, in_channels=in_channels))
            in_channels = channels[i+1]
            stage==i && return Chain(top,Chain(stages...))
      end
      bottom =  Chain(BatchNorm(channels[end]), ReLU(), Pool(window=(7,7), mode=1))
      stage==5 && return Chain(top,Chain(stages...), bottom)
      classifier = Chain(mat,Linear(input=in_channels,output=classes))
      return Chain(top,Chain(stages...),bottom,classifier)
  end

  configs = Dict(18  => (BasicV2, [2, 2, 2, 2], [64, 64, 128, 256, 512]),
                 34  => (BasicV2, [3, 4, 6, 3], [64, 64, 128, 256, 512]),
                 50  => (BottleneckV2, [3, 4, 6, 3], [64, 256, 512, 1024, 2048]),
                 101 => (BottleneckV2, [3, 4, 23, 3],[64, 256, 512, 1024, 2048]),
                 152 => (BottleneckV2, [3, 8, 36, 3],[64, 256, 512, 1024, 2048]))

  @inline (m::ResNet)(x) = m.layers(x)

  (m::ResNet)(x::Union{AbstractMatrix,AbstractString}) where {T,N} = m(preprocess(x))
  topK(labels::Vector{String},y;K=5) = labels[sortperm(vec(y);rev=true)[1:K]]

  function ResNet{N}(;trained=false, stage=0, mfile=KnetLayers.dir("examples","resnet$(N)v2.bson")) where N
     resnet = ResNet{N}(_init(configs[N]...; stage=stage),getLabels())
     if trained
       if !isfile(mfile)
          download(mfile,mfile) #FIXME
       end
       weights = BSON.load(mfile)
       loadResNet!(resnet,weights)
     end
     return resnet
  end

  Base.show(io::IO, ::ResNet{N}) where N = print(io, "ResNet{$N}()")

  ###
  #### Utils
  ###
  const atype = KnetLayers.arrtype
  transfer!(p::Param, x)  = transfer!(p.value,x)
  transfer!(p::KnetArray, x::AbstractArray) = p .= KnetArray(x)
  transfer!(p,x) =  p .= x
  to4D(x) = reshape(convert(atype,x),1,1,length(x),1)
  toArrType(x) = convert(atype,x)

  function getLabels(labels=KnetLayers.dir("examples","imagenet_labels.txt"))
    if !isfile(labels)
        download("https://github.com/ekinakyurek/KnetLayers.jl/releases/download/v0.2.0/imagenet_labels.txt",labels) #FIXME
    end
    return readlines(labels)
  end

  import Base: /
  /(a::RGB, b::RGB) = RGB(a.r/b.r, a.g/b.g, a.b/b.b)

  function preprocess(img::AbstractMatrix)
    # Resize such that smallest edge is 256 pixels long
      img = resize_smallest_dimension(img, 256)
      im  = (center_crop(img, 224) .- RGB(0.485, 0.456, 0.406)) ./ RGB(0.229, 0.224, 0.225)
      z   =  channelview(im)
      z1  = Float32.(permutedims(z, (3, 2, 1))[:,:,:,:]);
  end

  preprocess(img::AbstractString)  = preprocess(RGB.(load(img)))

  # Resize an image such that its smallest dimension is the given length
  function resize_smallest_dimension(im, len)
    reduction_factor = len/minimum(size(im)[1:2])
    new_size = size(im)
    new_size = (
        round(Int, size(im,1)*reduction_factor),
        round(Int, size(im,2)*reduction_factor),
    )
    if reduction_factor < 1.0
      # Images.jl's imresize() needs to first lowpass the image, it won't do it for us
      im = imfilter(im, KernelFactors.gaussian(0.75/reduction_factor), Inner())
    end
    return imresize(im, new_size)
  end

  # Take the len-by-len square of pixels at the center of image `im`
  function center_crop(im, len)
    l2 = div(len,2)
    adjust = len % 2 == 0 ? 1 : 0
    return im[div(end,2)-l2:div(end,2)+l2-adjust,div(end,2)-l2:div(end,2)+l2-adjust]
  end

  # Load Utils
  function loadBNLayer!(m::BatchNorm, weights, prefix)
    m.moments.var  = weights[Symbol(prefix*"running_var")]  |> to4D
    m.moments.mean = weights[Symbol(prefix*"running_mean")] |> to4D
    m.params       = vcat(weights[Symbol(prefix*"gamma")],weights[Symbol(prefix*"beta")]) |> toArrType
  end

  function loadConvLayer!(m::Conv, weights, prefix; binit=false)
    transfer!(m.weight, weights[Symbol(prefix*"weight")])
    if binit
      transfer!(m.bias.b,  weights[Symbol(prefix*"bias")])
    end
  end

  function loadDenseLayer!(m::Linear, weights, prefix; binit=true)
    transfer!(m.mult.weight, permutedims(weights[Symbol(prefix*"weight")],(2,1)))
    if binit
      transfer!(m.bias.b,  weights[Symbol(prefix*"bias")])
    end
  end

  function loadTop!(top::Chain, weights, prefix)
    loadBNLayer!(top[1], weights, string(prefix,"batchnorm0_"))
    loadConvLayer!(top[2], weights, string(prefix,"conv0_"))
    loadBNLayer!(top[3], weights, string(prefix,"batchnorm1_"))
  end

  function loadBottom!(bottom::Chain, weights, prefix)
    loadBNLayer!(bottom[1], weights, string(prefix,"batchnorm2_"))
  end

  function loadFinal!(bottom::Chain, weights, prefix)
    loadDenseLayer!(bottom[2], weights, string(prefix,"dense0_"))
  end

  function loadBlock!(block::BottleneckV2, weights, prefix, bn, conv)
    loadBNLayer!(block.layers[1], weights, string(prefix,"batchnorm$(bn)_"))
    loadConvLayer!(block.layers[3], weights, string(prefix,"conv$(conv)_"))
    bn+=1; conv+=1
    loadBNLayer!(block.layers[4], weights, string(prefix,"batchnorm$(bn)_"))
    loadConvLayer!(block.layers[6], weights, string(prefix,"conv$(conv)_"))
    bn+=1; conv+=1
    loadBNLayer!(block.layers[7], weights, string(prefix,"batchnorm$(bn)_"))
    loadConvLayer!(block.layers[9], weights, string(prefix,"conv$(conv)_"))
    bn+=1; conv+=1
    if block.downsample !== nothing
      loadConvLayer!(block.downsample, weights, string(prefix,"conv$(conv)_"))
      conv+=1
    end
    return bn, conv
  end

  function loadBlock!(block::BasicV2, weights, prefix, bn, conv)
    loadBNLayer!(block.layers[1], weights, string(prefix,"batchnorm$(bn)_"))
    loadConvLayer!(block.layers[3], weights, string(prefix,"conv$(conv)_"))
    bn+=1; conv+=1
    loadBNLayer!(block.layers[4], weights, string(prefix,"batchnorm$(bn)_"))
    loadConvLayer!(block.layers[6], weights, string(prefix,"conv$(conv)_"))
    bn+=1; conv+=1
    if block.downsample !== nothing
      loadConvLayer!(block.downsample, weights, string(prefix,"conv$(conv)_"))
      conv+=1
    end
    return bn, conv
  end

  function loadStage!(stage::Chain, weights, prefix)
    bn   = 0
    conv = 0
    for (i,block) in enumerate(stage.layers)
        bn, conv = loadBlock!(block,weights,prefix,bn,conv)
    end
  end

  const idmap = Dict{Int,Int}(18=>2,34=>3,50=>4,101=>5,152=>7)
  function loadResNet!(resnet::ResNet{N}, weights; prefix="resnetv2") where N
    prefix=prefix*string(idmap[N],"_");
    loadTop!(resnet.layers[1],weights,prefix)
    for (i,stage) in enumerate(resnet.layers[2])
        loadStage!(resnet.layers[2][i], weights, "$(prefix)stage$(i)_")
    end
    if length(resnet.layers) > 2
      loadBottom!(resnet.layers[3],weights,prefix)
    end
    if length(resnet.layers) > 3
      loadFinal!(resnet.layers[4],weights,prefix)
    end
    return resnet
  end

  export ResNet, topK
end


text:
@testset "Mathematics" begin
    @testset "Combinatorics" begin include("Combinatorics.jl") end
    @testset "VectorSpaces" begin include("VectorSpaces.jl") end
    @testset "AlgebraOverFields" begin include("AlgebraOverFields.jl") end
    @testset "QuantumNumbers" begin include("QuantumNumbers.jl") end
end


text:
# 2018 CODATA value
# https://en.wikipedia.org/wiki/Fine-structure_constant#Measurement
const c = 137.035999084
const α = inv(c)


text:
using YaoLocations
using MLStyle
using YaoHIR
using Test

module TestIntrinsic
using YaoHIR: @intrinsic

@intrinsic X
@intrinsic R(theta::T) where {T <: Real}
@intrinsic SWAP

end

@test YaoHIR.routine_name(TestIntrinsic.X) == :X
@test YaoHIR.routine_name(TestIntrinsic.R(1.0)) == :R
@test TestIntrinsic.R(1.0).theta == 1.0

display(YaoHIR.X)
display(YaoHIR.Rx(1.0))
display(YaoHIR.Operation(YaoHIR.X, 2.0))

circ = Chain(
    Gate(YaoHIR.X, Locations(1)),
    Core.SSAValue(1),
    Ctrl(Gate(Core.SSAValue(1), Locations(3)), CtrlLocations(2))
)

print(circ)

@test YaoHIR.leaves(circ) == [Gate(YaoHIR.X, Locations(1)),
    Core.SSAValue(1),
    Ctrl(Gate(Core.SSAValue(1), Locations(3)), CtrlLocations(2))
]


text:
# # Example
#
#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/example.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/example.ipynb)
#md #
# This example shows part of the basic functionality of the module.
#
#nb # If you prefer JupyterLab's interface over Jupyter Notebook's, replace `notebooks`
#nb # with `lab/tree` in the URL.
#nb #-
# First off, let's import the package:
using Scats
#md #+
# Let's set up an instance of the interface of the package to obtain access to its
# capabilities.
s = Scats.API();
#md #+
# We don't have a time series yet, but we can use the built-in generator to get one. To do
# this, we should specify time series generator parameters using a text file, but we do not
# have one either. Fortunately, the package provides a template for each built-in structure,
# so we'll use one of them.
## Create a temporary file for generator parameters
file, _ = mktemp()

## Fill this file with the generator settings template
s.Gen.example(file)

## Finally, read these parameters into the appropriate structure
s.Gen(file)
#md #+
# Let's generate a time series:
## Call the generation
s.gen!()
#md #+
# Done. The time series is now contained in the instance of [Input](@ref Input)' structure.
#
# Let's render it for validation. But first, we need a package for that.
using Pkg
Pkg.add("PyPlot")
#md #+
# And now let's import it:
using PyPlot
#md #+
# It's optional, but it's never harmful to make graphs prettier:
## Get access to rcParams
rcP = PyPlot.PyDict(PyPlot.matplotlib."rcParams")

## Set DPI
rcP["savefig.dpi"] = 300
rcP["figure.dpi"] = 300

## Activate TeX support
rcP["text.usetex"] = true

## Activate languages support
rcP["text.latex.preamble"] = raw"\usepackage[main=english]{babel}"

## Switch to Computer Modern font set
rcP["mathtext.fontset"] = "cm"

## Set font sizes
rcP["font.size"] = 18
rcP["legend.fontsize"] = 12;
#md #+
# Now let's get a plot:
## Create a plot
plot(s.Input.t, s.Input.x, color = "#425378")

## Add a title
title(raw"\textrm{Time series}")

## Add labels
xlabel(raw"\textrm{Time}")
ylabel(raw"\textrm{Value}")

## Create a directory for figures
!isdir("figures") && mkdir("figures")

## Save the figure
savefig("figures/input.pdf", bbox_inches = "tight")
gcf()

#-
## Close the figure
close()


text:
include("PDI_exact_resolution.jl")

include("InstanceLoader.jl")
include("Helper.jl")
include("GraphHelper.jl")

INSTANCE_PATH = "../PRP_instances/A_014_#ABS1_15_1.prp"
#INSTANCE_PATH = "../PRP_instances/A_050_ABS14_50_1.prp"

function BranchAndCutPDI(params, nodes, demands, costs)

	model = createPDI_Bard_Nananukul_compacte(params, nodes, demands, costs)

	function lazySep(c_data)

		n = params["n"]
		l = params["l"]
		Q = params["Q"]

		for t=1:l

			fullHisto = Set{Int64}()
			sousTours = Set{Int64}[]

			for i in 1:n

				#Si noeud pas déjà dans un circuit trouvé
				if !( i in fullHisto ) && ( round( callback_value(c_data, variable_by_name(model, "z[$i,$t]")) ) == 1 )

					foundZero = false

					current = i

					##Initiliser historique avec i
					histo = Set{Int64}(current)
					#push!(histo, current)

					##Pour faire les deux sens
					##INUTILE ENFAIT
					#for reverse in range 1:2

					while true

						#On cherche le noeud précédent
						for j in 0:n
							if current != j

								#Si on trouve le noeud précédent
								if round( callback_value(c_data, variable_by_name(model, "x[$j,$current,$t]")) ) == 1
									current = j
									break
								end

							end
						end

						#Fin du circuit avec passage en 0
						if current == 0
							foundZero = true
							break
						end

						#Fin du circuit mais pas de passage en 0
						if current in histo
							break
						end

						push!(histo, current)

					end

					#Si pas de zero trouvé, sous-tour
					if !foundZero
						push!(sousTours, histo)
					end

					#Ajouter le cycle trouver a l'historique pour ne pas refaire l'algorithme dessus
					fullHisto = union(fullHisto, histo)

				end
			end

			#Set complet
			fullSet = Set(0:n)

			for sousTour in sousTours

				#Set sans le sous tour
				setNotInSousTour = setdiff(fullSet, sousTour)

				#FCCs
				if length(sousTour)>=1
					con = @build_constraint(sum( variable_by_name(model, "x[$i,$j,$t]") for i in setNotInSousTour, j in sousTour) >= sum( variable_by_name(model, "q[$i,$t]") for i in sousTour)/Q )
					MOI.submit(model, MOI.LazyConstraint(c_data), con)
				end

				#FCCs Reversed Q
				#con = @build_constraint(Q * sum( variable_by_name(model, "x[$i,$j,$t]") for i in setNotInSousTour, j in sousTour) >= sum( variable_by_name(model, "q[$i,$t]") for i in sousTour) )
				#MOI.submit(model, MOI.LazyConstraint(c_data), con)

				#GFSECs
				#if length(sousTour)>=2
					#con = @build_constraint(sum( i != j ? variable_by_name(model, "x[$i,$j,$t]") : 0 for i in sousTour, j in sousTour) <= length(sousTour) - sum( variable_by_name(model, "q[$i,$t]") for i in sousTour)/Q )
					#MOI.submit(model, MOI.LazyConstraint(c_data), con)
				#end
				#GFSECs Reversed Q by Adulyasak
				#if length(sousTour)>=2
				#con = @build_constraint( Q * sum( i != j ? variable_by_name(model, "x[$i,$j,$t]") : 0 for i in sousTour, j in sousTour) <= sum(Q*variable_by_name(model, "z[$i,$t]") - variable_by_name(model, "q[$i,$t]") for i in sousTour) )
				#MOI.submit(model, MOI.LazyConstraint(c_data), con)
				#end

				#Pour user cut
				#MOI.submit(LP, MOI.UserCut(cb_data), con)
			end

		end
	end

	#Mauvais algorithme
	#=
	function userSep(c_data)

		n = params["n"]
		l = params["l"]
		Q = params["Q"]

		for t=1:l

			fullHisto = Set{Int64}()
			sousTours = Set{Int64}[]

			for i in 1:n

				#Si noeud pas déjà dans un circuit trouvé
				if !( i in fullHisto ) && ( callback_value(c_data, variable_by_name(model, "z[$i,$t]")) > 0 )

					# println("z[$i,$t]", callback_value(c_data, variable_by_name(model, "z[$i,$t]")))
					###DEBUG

					# for j in 0:n
					# 	if i != j
					# 		println("x[$j,$i,$t]", callback_value(c_data, variable_by_name(model, "x[$j,$i,$t]")))
					# 	end
					# end
					# continue

					foundZero = false

					current = i

					##Initiliser historique avec i
					histo = Set{Int64}(current)
					#push!(histo, current)

					##Pour faire les deux sens
					##INUTILE ENFAIT
					#for reverse in range 1:2

					while true

						#On cherche le noeud précédent
						for j in 0:n
							if current != j

								#Si on trouve le noeud précédent
								if callback_value(c_data, variable_by_name(model, "x[$j,$current,$t]")) > 0
									current = j
									break
								end

							end
						end

						#Fin du circuit avec passage en 0
						if current == 0
							#println("ZERO")
							foundZero = true
							break
						end

						#Fin du circuit mais pas de passage en 0
						if current in histo
							#println("DANS HISTO")
							break
						end

						push!(histo, current)

					end

					#Si pas de zero trouvé, sous-tour
					if !foundZero
						push!(sousTours, histo)
					end

					#Ajouter le cycle trouver a l'historique pour ne pas refaire l'algorithme dessus
					fullHisto = union(fullHisto, histo)

				end
			end

			#Set complet
			fullSet = Set(0:n)

			for sousTour in sousTours

				#Set sans le sous tour
				setNotInSousTour = setdiff(fullSet, sousTour)

				#con = @build_constraint(sum( variable_by_name(model, "x[$i,$j,$t]") for i in setNotInSousTour, j in sousTour) >= sum( variable_by_name(model, "q[$i,$t]") for i in sousTour)/Q )
				#MOI.submit(model, MOI.UserCut(c_data), con)

				con = @build_constraint( Q * sum( i != j ? variable_by_name(model, "x[$i,$j,$t]") : 0 for i in sousTour, j in sousTour) <= sum(Q*variable_by_name(model, "z[$i,$t]") - variable_by_name(model, "q[$i,$t]") for i in sousTour) )
				MOI.submit(model, MOI.UserCut(c_data), con)

			end

		end

	end
	=#

	# our userSep_ViolatedAcyclicCst function sets a LazyConstraintCallback of CPLEX
	MOI.set(model, MOI.LazyConstraintCallback(), lazySep)

	#MOI.set(model, MOI.UserCutCallback(), userSep)
	#set_time_limit_sec(model,10800) #time limit de 3h
	optimize!(model)

	return model
end

function PDItoCircuits(model, params, nodes, demands, costs)

	allCircuits = []

	for t in 1:params["l"]
		push!(allCircuits, vrpToCircuit(model, params, true, t))
	end

	return allCircuits

end


text:
module Chebyshev

using FFTW
using Future: copy!

import Base: *, \
struct Transform{T,I}
    dct!::FFTW.r2rFFTWPlan{T,(3,),true,1}
end

function plan_transform!(x::Vector{T}) where T
    Transform{T,true}(FFTW.plan_r2r!(x, FFTW.REDFT00))
end

function plan_transform(x::Vector{T}) where T
    Transform{T,false}(FFTW.plan_r2r!(x, FFTW.REDFT00))
end

(C::Transform{T,true}  * x::Vector{T}) where T = transform!(x, C.dct!)
(C::Transform{T,false} * x::Vector{T}) where T = transform(x, C.dct!)
(C::Transform{T,true}  \ A::Vector{T}) where T = inv_transform!(A, C.dct!)
(C::Transform{T,false} \ A::Vector{T}) where T = inv_transform(A, C.dct!)

"""
    Chebyshev.nodes(N, T = Float64)

Gives `N` extrema Chebyshev points of type `T` in [-1 ,1].
"""
function nodes(N, T = Float64)
    T[cos(θ) for θ in range(π, 0, length=N)]
end

"""
    Chebyshev.firstkind(C, s[, offset = 0])

Evaluate a Chebyshev series of the first kind.

# Arguments

- `C`: a vector with the coefficients of the series
- `s`: a scalar or an array of evaluation points ∈ [-1, 1]
- `offset`: an optional starting index of the series summation
"""
firstkind(A, ss, offset = 0) = [firstkind(A, s, offset) for s in ss]
function firstkind(A::AbstractArray{C}, s::S, offset = 0)::C where {C,S <: Real}
    -1 ≤ s ≤ 1 || throw(DomainError("s ∉ [-1,1]"))
    T₋ = s
    T  = one(S)

    for i in 1:offset
        T₋, T = T, 2s*T - T₋
    end

    N = length(A)
    f = zero(C)
    for i in 1+offset:N
        f += A[i]*T
        T₋, T = T, 2s*T - T₋
    end
    f
end

"""
    Chebyshev.secondkind(C, s[, offset = 0])

Evaluate a Chebyshev series of the second kind.

# Arguments

- `C`: a vector with the coefficients of the series
- `s`: a scalar or an array of evaluation points ∈ [-1, 1]
- `offset`: an optional starting index of the series summation
"""
secondkind(C, ss, offset = 0) = [secondkind(C, s, offset) for s in ss]
function secondkind(A::AbstractArray{C}, s::S, offset = 0)::C where {C,S <: Real}
    -1 ≤ s ≤ 1 || throw(DomainError("s ∉ [-1,1]"))
    U₋ = zero(s)
    U  = one(S)

    for i in 1:offset
        U₋, U = U, 2s*U - U₋
    end

    N = length(A)
    f = zero(C)
    for i in 1+offset:N
        f += A[i]*U
        U₋, U = U, 2s*U - U₋
    end
    f
end

"""
    Chebyshev.transform(x[, plan!])

Compute the discrete Chebyshev transform of `x`.

# Arguments

- `x`: samples of a function distributed along [extrema Chebyshev nodes](@ref Chebyshev.nodes)
- `plan!`: an optional pre-planned in-place DCT-I used to compute the transform
"""
function transform(x, plan! = FFTW.plan_r2r!(x, FFTW.REDFT00))
    C = copy(x)
    transform!(C, plan!)
end

"""
    Chebyshev.transform!(A,[ x, plan!])

Compute the discrete Chebyshev transform of `x` in-place.

# Arguments

- `A`: the output vector (also the input vector if `x` is omitted)
- `x`: an optional input vector with samples of a function distributed
  along [extrema Chebyshev nodes](@ref Chebyshev.nodes)
- `plan!`: an optional pre-planned in-place DCT-I used to compute the transform
"""
function transform!(x, plan! = FFTW.plan_r2r!(x, FFTW.REDFT00))
    N = length(x)
    if N != length(plan!)
        throw(BoundsError("`x` must have the same size as the preplanned DCT"))
    end

    plan!*x

    s = 1/(N-1)

    @inbounds begin
        for n in 2:2:N
            x[n-1] *= s
            x[n]   *= -s
        end
        if isodd(N)
            x[N] *= s
        end

        x[1] /= 2
        x[N] /= 2
    end

    x
end

function transform!(A::T, x::T, plan! = FFTW.plan_r2r!(x, FFTW.REDFT00)) where T
    copy!(A, x)
    transform!(A, plan!)
end

"""
    Chebyshev.inv_transform(A[, plan!])

Perform the inverse Chebyshev transform.

This is the inverse of [`Chebyshev.transform`](@ref).

# Fields

- `A`: the coefficients of a Chebyshev series
- `plan!`: an optional pre-planned in-place DCT-I used to compute the transform
"""
function inv_transform(A, plan! = FFTW.plan_r2r!(A, FFTW.REDFT00))
    x = copy(A)
    inv_transform!(x, plan!)
end

"""
    Chebyshev.inv_transform!(x, [C, plan!)

Perform the inverse Chebyshev transform in place.

This is the inverse of [`Chebyshev.transform!`](@ref)

# Arguments

- `x`: the output vector to store the reconstructed function with
  points distributed along extrema Chebyshev nodes (also the input
  vector if `A` is omitted)
- `A`: an optional input vector with coefficients of a Chebyshev series
- `plan!`: an optional pre-planned in-place DCT-I used to compute the transform
"""
function inv_transform!(A, plan! = FFTW.plan_r2r!(A, FFTW.REDFT00))
    N = length(A)
    if N != length(plan!)
        throw(BoundsError("`A` must have the same size as the preplanned DCT"))
    end

    @inbounds begin
        for n in 2:2:N-2
            A[n]   *= -0.5
            A[n+1] *=  0.5
        end
        if iseven(N)
            A[N] *= -1
        else
            A[N-1] *= -0.5
        end
    end

    plan!*A
end

function inv_transform!(x::T, A::T, plan! = FFTW.plan_r2r!(A, FFTW.REDFT00)) where T
    copy!(x, A)
    Chebyshev.inv_transform!(x, plan!)
end

"""
    clenshaw_curtis_weights(N)

Returns the `N` weights used in the Clenshaw Curtis quadrature
"""
function clenshaw_curtis_weights(N)
    n = N - 1
    w = zeros(N)

    # Unlike the documentation, we move the 1/(N-1)
    # factor inside the r vector.
    w[1] = 1.0/n
    for k in 2:2:n
        w[k+1] = (1.0 + (-1)^k)/(1 - k^2)/n
    end

    # Scale by S
    w[1] *= 2; w[N] *= 2
    # Multiply by D
    FFTW.r2r!(w, FFTW.REDFT00)
    # Scale by S⁻¹
    w[1] /= 2; w[N] /= 2

    return w
end

end


text:
using Plots, LaTeXStrings
gr(grid=false,legend=false,titlefontsize=12,size=(500,300),transpose=true,colorbar=false)

using FourierGPE

## set simulation parameters ====
L = (40.0,)
N = (512,)
sim = Sim(L,N)
@unpack_Sim sim

μ = 25.0
## potential
import FourierGPE.V
V(x,t) = 0.5*x^2

## TF state
ψ0(x,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,0.0)/μ,0.0)+im*0.0)

## initial state
x = X[1]
ψi = ψ0.(x,μ,g)
ϕi = kspace(ψi,sim)
@pack_Sim! sim

## imaginary time to find ground state
sol = runsim(sim)

## pull out and check the ground state
ϕg = sol[end]
ψg = xspace(ϕg,sim)
plot(x,one.(x)*μ,ls=:solid,c=:gray,w=1)
plot!(x,V.(x,0.0),c=:red,w=4,alpha=0.4)
plot!(x,g*abs2.(ψg),c=c3,w=2,fill=(0,0.4,c3),size=(600,200))
xlims!(-10,10); ylims!(0,1.3*μ)
title!(L"\textrm{local}\; \mu(x)")
xlabel!(L"x/a_x"); ylabel!(L"\mu(x)/\hbar\omega_x")

## imprint dark soliton
ψf = xspace(sol[end],sim)
c = sqrt(μ)
ξ = 1/c
v = 0.5*c
xs = 0.
f = sqrt(1-(v/c)^2)

ψs = ψf.*(f*tanh.(f*(x .-xs)/ξ).+im*v/c);
showpsi(x,ψs)
xlims!(-10,10)

## set simulation parameters
γ = 0.0
tf = 8*pi/sqrt(2); t = LinRange(ti,tf,Nt)

# define new simulation using previous parameters
simSoliton = Sim(sim;γ=γ,tf=tf,t=t)
ϕi = kspace(ψs,simSoliton)
@pack! simSoliton = ϕi

## evolve in k space
sols = runsim(simSoliton)

## plot

ϕf = sols[end-3]
ψf = xspace(ϕf,simSoliton)
showpsi(x,ψf)

anim = @animate for i in 1:length(t)-3
    ψ = xspace(sols[i],simSoliton)
    y = g*abs2.(ψ)
    plot(x,y,c=c3,w=2,fill=(0,0.4,c3),size=(600,300),grid=false)
    xlims!(-10,10); ylims!(0,1.3*μ)
    title!(L"\textrm{local}\; \mu(x)")
    xlabel!(L"x/a_x"); ylabel!(L"\mu(x)/\hbar\omega_x")
end

gif(anim,"./examples/soliton.gif",fps=25)