text:
# Stella, Themelis, Sopasakis, Patrinos, "A simple and efficient algorithm
# for nonlinear model predictive control", 56th IEEE Conference on Decision
# and Control (2017).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

"""
    PANOCIteration(; <keyword-arguments>)

Instantiate the PANOC algorithm (see [1]) for solving optimization problems
of the form

    minimize f(Ax) + g(x),

where `f` is smooth and `A` is a linear mapping (for example, a matrix).

# Arguments
- `x0`: initial point.
- `f=Zero()`: smooth objective term.
- `A=I`: linear operator (e.g. a matrix).
- `g=Zero()`: proximable objective term.
- `Lf=nothing`: Lipschitz constant of the gradient of x ↦ f(Ax).
- `gamma=nothing`: stepsize to use, defaults to `1/Lf` if not set (but `Lf` is).
- `adaptive=false`: forces the method stepsize to be adaptively adjusted.
- `minimum_gamma=1e-7`: lower bound to `gamma` in case `adaptive == true`.
- `max_backtracks=20`: maximum number of line-search backtracks.
- `H=LBFGS(x0, 5)`: variable metric to use to compute line-search directions.

# References
- [1] Stella, Themelis, Sopasakis, Patrinos, "A simple and efficient algorithm
for nonlinear model predictive control", 56th IEEE Conference on Decision
and Control (2017).
"""

@Base.kwdef struct PANOCIteration{R,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},Tf,TA,Tg,TH}
    f::Tf = Zero()
    A::TA = I
    g::Tg = Zero()
    x0::Tx
    alpha::R = real(eltype(x0))(0.95)
    beta::R = real(eltype(x0))(0.5)
    Lf::Maybe{R} = nothing
    gamma::Maybe{R} = Lf === nothing ? nothing : (alpha / Lf)
    adaptive::Bool = false
    minimum_gamma::R = real(eltype(x0))(1e-7)
    max_backtracks::Int = 20
    H::TH = LBFGS(x0, 5)
end

Base.IteratorSize(::Type{<:PANOCIteration}) = Base.IsInfinite()

Base.@kwdef mutable struct PANOCState{R,Tx,TAx,TH}
    x::Tx             # iterate
    Ax::TAx           # A times x
    f_Ax::R           # value of smooth term
    grad_f_Ax::TAx    # gradient of f at Ax
    At_grad_f_Ax::Tx  # gradient of smooth term
    gamma::R          # stepsize parameter of forward and backward steps
    y::Tx             # forward point
    z::Tx             # forward-backward point
    g_z::R            # value of nonsmooth term (at z)
    res::Tx           # fixed-point residual at iterate (= x - z)
    H::TH             # variable metric
    tau::Maybe{R} = nothing
    x_prev::Tx = zero(x)
    res_prev::Tx = zero(x)
    d::Tx = zero(x)
    Ad::TAx = zero(Ax)
    x_d::Tx = zero(x)
    Ax_d::TAx = zero(Ax)
    f_Ax_d::R = zero(real(eltype(x)))
    grad_f_Ax_d::TAx = zero(Ax)
    At_grad_f_Ax_d::Tx = zero(x)
    z_curr::Tx = zero(x)
end

f_model(iter::PANOCIteration, state::PANOCState) =
    f_model(state.f_Ax, state.At_grad_f_Ax, state.res, iter.alpha / state.gamma)

function Base.iterate(iter::PANOCIteration{R}) where {R}
    x = copy(iter.x0)
    Ax = iter.A * x
    grad_f_Ax, f_Ax = gradient(iter.f, Ax)

    gamma = iter.gamma

    if gamma === nothing
        # compute lower bound to Lipschitz constant of the gradient of x ↦ f(Ax)
        xeps = x .+ R(1)
        grad_f_Axeps, f_Axeps = gradient(iter.f, iter.A * xeps)
        L = norm(iter.A' * (grad_f_Axeps - grad_f_Ax)) / R(sqrt(length(x)))
        gamma = iter.alpha / L
    end

    # compute initial forward-backward step
    At_grad_f_Ax = iter.A' * grad_f_Ax
    y = x - gamma .* At_grad_f_Ax
    z, g_z = prox(iter.g, y, gamma)

    # compute initial fixed-point residual
    res = x - z

    state = PANOCState(
        x=x, Ax=Ax, f_Ax=f_Ax, grad_f_Ax=grad_f_Ax, At_grad_f_Ax=At_grad_f_Ax,
        gamma=gamma, y=y, z=z, g_z=g_z, res=res, H=iter.H,
    )

    return state, state
end

function Base.iterate(
    iter::PANOCIteration{R},
    state::PANOCState{R,Tx,TAx},
) where {R,Tx,TAx}
    Az, f_Az, grad_f_Az, At_grad_f_Az = nothing, nothing, nothing, nothing
    a, b, c = nothing, nothing, nothing

    f_Az_upp = f_model(iter, state)

    # backtrack gamma (warn and halt if gamma gets too small)
    while iter.gamma === nothing || iter.adaptive == true
        if state.gamma < iter.minimum_gamma
            @warn "parameter `gamma` became too small ($(state.gamma)), stopping the iterations"
            return nothing
        end
        Az = iter.A * state.z
        grad_f_Az, f_Az = gradient(iter.f, Az)
        tol = 10 * eps(R) * (1 + abs(f_Az))
        if f_Az <= f_Az_upp + tol
            break
        end
        state.gamma *= 0.5
        state.y .= state.x .- state.gamma .* state.At_grad_f_Ax
        state.g_z = prox!(state.z, iter.g, state.y, state.gamma)
        state.res .= state.x .- state.z
        reset!(state.H)
        f_Az_upp = f_model(iter, state)
    end

    # compute FBE
    FBE_x = f_Az_upp + state.g_z

    # compute direction
    mul!(state.d, state.H, -state.res)

    # store iterate and residual for metric update later on
    state.x_prev .= state.x
    state.res_prev .= state.res

    # backtrack tau 1 → 0
    tau = R(1)
    mul!(state.Ad, iter.A, state.d)

    state.x_d .= state.x .+ state.d
    state.Ax_d .= state.Ax .+ state.Ad
    state.f_Ax_d = gradient!(state.grad_f_Ax_d, iter.f, state.Ax_d)
    mul!(state.At_grad_f_Ax_d, adjoint(iter.A), state.grad_f_Ax_d)

    state.x .= state.x_d
    state.Ax .= state.Ax_d
    state.f_Ax = state.f_Ax_d
    state.grad_f_Ax .= state.grad_f_Ax_d
    state.At_grad_f_Ax .= state.At_grad_f_Ax_d

    copyto!(state.z_curr, state.z)

    sigma = iter.beta * (0.5 / state.gamma) * (1 - iter.alpha)
    tol = 10 * eps(R) * (1 + abs(FBE_x))
    threshold = FBE_x - sigma * norm(state.res)^2 + tol

    for _ in 1:iter.max_backtracks
        state.y .= state.x .- state.gamma .* state.At_grad_f_Ax
        state.g_z = prox!(state.z, iter.g, state.y, state.gamma)
        state.res .= state.x .- state.z
        FBE_x_new = f_model(iter, state) + state.g_z

        if FBE_x_new <= threshold
            # update metric
            update!(state.H, state.x - state.x_prev, state.res - state.res_prev)
            state.tau = tau
            return state, state
        end

        if Az === nothing
            Az = iter.A * state.z_curr
        end

        tau *= 0.5
        state.x .= tau .* state.x_d .+ (1 - tau) .* state.z_curr
        state.Ax .= tau .* state.Ax_d .+ (1 - tau) .* Az

        if ProximalOperators.is_quadratic(iter.f)
            # in case f is quadratic, we can compute its value and gradient
            # along a line using interpolation and linear combinations
            # this allows saving operations
            if grad_f_Az === nothing
                grad_f_Az, f_Az = gradient(iter.f, Az)
            end
            if At_grad_f_Az === nothing
                At_grad_f_Az = iter.A' * grad_f_Az
                c = f_Az
                b = real(dot(state.Ax_d .- Az, grad_f_Az))
                a = state.f_Ax_d - b - c
            end
            state.f_Ax = a * tau^2 + b * tau + c
            state.grad_f_Ax .= tau .* state.grad_f_Ax_d .+ (1 - tau) .* grad_f_Az
            state.At_grad_f_Ax .= tau .* state.At_grad_f_Ax_d .+ (1 - tau) .* At_grad_f_Az
        else
            # otherwise, in the general case where f is only smooth, we compute
            # one gradient and matvec per backtracking step
            state.f_Ax = gradient!(state.grad_f_Ax, iter.f, state.Ax)
            mul!(state.At_grad_f_Ax, adjoint(iter.A), state.grad_f_Ax)
        end
    end

    @warn "stepsize `tau` became too small ($(tau)), stopping the iterations"
    return nothing
end

# Solver

struct PANOC{R, K}
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int
    kwargs::K
end

function (solver::PANOC)(x0; kwargs...)
    stop(state::PANOCState) = norm(state.res, Inf) / state.gamma <= solver.tol
    disp((it, state)) = @printf(
        "%5d | %.3e | %.3e | %.3e\n",
        it,
        state.gamma,
        norm(state.res, Inf) / state.gamma,
        (state.tau === nothing ? 0.0 : state.tau)
    )
    iter = PANOCIteration(; x0=x0, solver.kwargs..., kwargs...)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose
        iter = tee(sample(iter, solver.freq), disp)
    end
    num_iters, state_final = loop(iter)
    return state_final.z, num_iters
end

PANOC(; maxit=1_000, tol=1e-8, verbose=false, freq=10, kwargs...) =
    PANOC(maxit, tol, verbose, freq, kwargs)


text:
using
  PyPlot,
  JLD2,
  Printf,
  FourierFlows

using FFTW: ifft

import PassiveTracerFlows.TracerAdvDiff

# Numerical parameters and time-stepping parameters
nx  = 128         # 2D resolution = nx^2
stepper = "RK4"   # timestepper
dt  = 0.02        # timestep
nsubs  = 200      # number of time-steps for plotting
nsteps = 4nsubs   # total number of time-steps (must be multiple of nsubs)

# Physical parameters
Lx  = 2π      # domain size
kap = 0.002   # diffusivity

gr = TwoDGrid(nx, Lx)
x, y = gridpoints(gr)

# streamfunction (for plotting) and (u,v) flow field
psiampl = 0.2
m, n = 1, 1
psiin = @. psiampl * cos(m*x) * cos(n*y)
uvel(x, y) = +psiampl * n * cos(m*x) * sin(n*y)
vvel(x, y) = -psiampl * m * sin(m*x) * cos(n*y)

prob = TracerAdvDiff.Problem(; steadyflow=true,
    nx=nx, Lx=Lx, kap=kap, u=uvel, v=vvel, dt=dt, stepper=stepper)

sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

# Initial condition c0 = c(x, y, t=0)
amplc0, sigc0 = 0.1, 0.1
c0func(x, y) = amplc0*exp(-(x^2+y^2)/(2sigc0^2))
c0 = @. c0func(x-0.2Lx, y)

TracerAdvDiff.set_c!(prob, c0)

"Plot the concentration field and the (u, v) streamlines."
function plotoutput(prob, fig, axs; drawcolorbar=false)
  cl, v, p, g = prob.clock, prob.vars, prob.params, prob.grid

  t = round(cl.t, digits=2)
  TracerAdvDiff.updatevars!(prob)

  cla()
  pcolormesh(x, y, v.c)

  if drawcolorbar; colorbar(); end

  contour(x, y, psiin, 15, colors="k", linewidths=0.7)
  axis("equal")
  axis("square")
  title("shading: \$c(x, y, t= $t )\$, contours: \$\\psi(x, y)\$")

  pause(0.001)

  nothing
end

fig, axs = subplots(ncols=1, nrows=1, figsize=(8, 8))
plotoutput(prob, fig, axs; drawcolorbar=true)

while cl.step < nsteps
  stepforward!(prob, nsubs)
  plotoutput(prob, fig, axs)
end

ext: jl


text:
function assign_scenarios_to_procs(scen_tree::ScenarioTree
                                   )::Dict{Int,Set{ScenarioID}}
    sp_map = Dict{Int, Set{ScenarioID}}()

    nprocs = nworkers()
    wrks = workers()
    sp_map = Dict(w => Set{ScenarioID}() for w in wrks)

    for (k,s) in enumerate(scenarios(scen_tree))
        w = wrks[(k-1) % nprocs + 1]
        push!(sp_map[w], s)
    end

    return sp_map
end

function _initialize_subproblems(sp_map::Dict{Int,Set{ScenarioID}},
                                 wi::WorkerInf,
                                 scen_tree::ScenarioTree,
                                 constructor::Function,
                                 constructor_args::Tuple,
                                 r::AbstractPenaltyParameter,
                                 warm_start::Bool,
                                 subproblem_callbacks::Vector{SubproblemCallback};
                                 kwargs...
                                 )

    # Send initialization commands
    @sync for (wid, scenarios) in pairs(sp_map)
        @async _send_message(wi,
                             wid,
                             Initialize(constructor,
                                        constructor_args,
                                        (;kwargs...),
                                        typeof(r),
                                        scenarios,
                                        scen_tree,
                                        warm_start,
                                        subproblem_callbacks)
                             )
    end

    # Wait for and process initialization replies
    var_maps = Dict{ScenarioID,Dict{VariableID,VariableInfo}}()
    remaining_maps = copy(scenarios(scen_tree))

    while !isempty(remaining_maps)

        msg = _retrieve_message_type(wi, VariableMap)
        var_maps[msg.scen] = msg.var_info

        delete!(remaining_maps, msg.scen)

    end

    return var_maps
end

function _initialize_lb_subproblems(sp_map::Dict{Int,Set{ScenarioID}},
                                    wi::WorkerInf,
                                    scen_tree::ScenarioTree,
                                    constructor::Function,
                                    constructor_args::Tuple,
                                    warm_start::Bool;
                                    kwargs...
                                    )::Nothing

    @sync for (wid, scenarios) in pairs(sp_map)
        @async _send_message(wi,
                             wid,
                             InitializeLowerBound(constructor,
                                                  constructor_args,
                                                  (;kwargs...),
                                                  scenarios,
                                                  scen_tree,
                                                  warm_start)
                             )
    end

    return
end

function _set_initial_values(phd::PHData,
                             wi::WorkerInf,
                             )::Float64

    # Wait for and process mapping replies
    remaining_init = copy(scenarios(phd.scenario_tree))

    while !isempty(remaining_init)

        msg = _retrieve_message_type(wi, ReportBranch)

        _verify_report(msg)
        _update_values(phd, msg)

        delete!(remaining_init, msg.scen)

    end

    # Start computation of initial PH values -- W values are computed at
    # the end of `_set_penalty_parameter`.
    xhat_res_sq = compute_and_save_xhat(phd)

    return xhat_res_sq
end

function _set_penalty_parameter(phd::PHData,
                                wi::WorkerInf,
                                )::Float64

    if is_subproblem_dependent(typeof(phd.r))

        # Wait for and process penalty parameter messages
        remaining_maps = copy(scenarios(phd.scenario_tree))

        while !isempty(remaining_maps)

            msg = _retrieve_message_type(wi, PenaltyInfo)

            process_penalty_subproblem(phd.r, phd, msg.scen, msg.penalty)

            delete!(remaining_maps, msg.scen)

        end

    end

    if is_initial_value_dependent(typeof(phd.r))

        # Compute penalty parameter values
        process_penalty_initial_value(phd.r, phd)

    end

    if is_variable_dependent(typeof(phd.r))

        # Convert dictionary to subproblem format
        pen_map = Dict{ScenarioID,Dict{VariableID,Float64}}(
            s => Dict{VariableID,Float64}() for s in scenarios(phd)
        )

        for (xhid,penalty) in pairs(penalty_map(phd.r))
            for vid in convert_to_variable_ids(phd, xhid)
                pen_map[scenario(vid)][vid] = penalty
            end
        end

    else

        pen_map = Dict{ScenarioID,Float64}(
            s => get_penalty_value(phd.r) for s in scenarios(phd)
        )

    end

    # Send penalty parameter values to workers
    @sync for (scid, sinfo) in pairs(phd.scenario_map)
        @async begin
            wid = sinfo.pid
            penalty = pen_map[scid]
            _send_message(wi, wid, PenaltyInfo(scid, penalty))
        end
    end

    # Complete computation of initial PH values -- Xhat values are computed at
    # the end of `_set_initial_values`.
    x_res_sq = compute_and_save_w(phd)

    return x_res_sq
end

function initialize(scen_tree::ScenarioTree,
                    model_constructor::Function,
                    r::AbstractPenaltyParameter,
                    user_sp_map::Dict{Int,Set{ScenarioID}},
                    warm_start::Bool,
                    timo::TimerOutputs.TimerOutput,
                    report::Int,
                    lower_bound::Int,
                    subproblem_callbacks::Vector{SubproblemCallback},
                    constructor_args::Tuple,;
                    kwargs...
                    )::Tuple{PHData,WorkerInf}

    # Assign scenarios to processes
    scen_proc_map = isempty(user_sp_map) ? assign_scenarios_to_procs(scen_tree) : user_sp_map
    scen_per_worker = maximum(length.(collect(values(scen_proc_map))))
    n_scenarios = length(scenarios(scen_tree))

    # Start worker loops
    if report > 0
        println("...launching workers...")
    end
    worker_inf = @timeit(timo,
                         "Launch Workers",
                         _launch_workers(2*scen_per_worker, n_scenarios)
                         )

    # Initialize subproblems
    if report > 0
        println("...initializing subproblems...")
        flush(stdout)
    end

    var_map = @timeit(timo,
                      "Initialize Subproblems",
                      _initialize_subproblems(scen_proc_map,
                                              worker_inf,
                                              scen_tree,
                                              model_constructor,
                                              constructor_args,
                                              r,
                                              warm_start,
                                              subproblem_callbacks;
                                              kwargs...)
                      )

    if lower_bound > 0

        if report > 0
            println("...initializing lower-bound subproblems...")
            flush(stdout)
        end

        @timeit(timo,
                "Initialze Lower Bound Subproblems",
                _initialize_lb_subproblems(scen_proc_map,
                                           worker_inf,
                                           scen_tree,
                                           model_constructor,
                                           constructor_args,
                                           warm_start;
                                           kwargs...)
                )

    end

    # Construct master ph object
    phd = @timeit(timo,
                  "Other",
                  PHData(r,
                         scen_tree,
                         scen_proc_map,
                         var_map,
                         timo)
                  )

    # Initial values
    xhat_res_sq = @timeit(timo,
                          "Initial Values",
                          _set_initial_values(phd, worker_inf)
                          )

    # Check for penalty parameter update
    x_res_sq = @timeit(timo,
                       "Penalty Parameter",
                       _set_penalty_parameter(phd,
                                              worker_inf)
                       )

    # Save residual
    _save_residual(phd, -1, xhat_res_sq, x_res_sq, 0.0, 0.0)

    return (phd, worker_inf)
end

function build_extensive_form(tree::ScenarioTree,
                              model_constructor::Function,
                              constructor_args::Tuple,
                              sub_type::Type{S};
                              kwargs...
                              )::JuMP.Model where {S <: AbstractSubproblem}

    model = JuMP.Model()
    JuMP.set_objective_sense(model, MOI.MIN_SENSE)
    JuMP.set_objective_function(model, zero(JuMP.QuadExpr))

    # Below for mapping subproblem variables onto existing extensive form variables
    node_var_map = ef_node_dict_constructor(sub_type)

    for s in scenarios(tree)

        smod = model_constructor(s, constructor_args...; kwargs...)

        snode_var_map = ef_copy_model(model, smod, s, tree, node_var_map)

        @assert(isempty(intersect(keys(node_var_map), keys(snode_var_map))))
        merge!(node_var_map, snode_var_map)
    end

    return model
end


text:
mutable struct Agent
    id::Int64
    opinion::Float64
    inclin_interact::Float64
    perceiv_publ_opinion::Float64
    inactive_ticks::Int16
    feed::Array{Post, 1}
    function Agent(id, opinion)
        # check if opinion value is valid
        if opinion < -1 || opinion > 1
            error("invalid opinion value")
        end

        new(
            id,
            opinion,
            generate_inclin_interact(),
            opinion,
            0,
            Array{Post, 1}(undef, 0)
        )
    end
end

# suppress output of include()
;


text:
function trunc(x)
if typeof(x) == Float64
return 0.0
elseif typeof(x) == Float32
return Float32(0.0)
elseif typeof(x) == Int64
return 0
end
end

max_stars_repo_path: Module 2/Chapter03/truncomp.jl

max_issues_repo_path: Module 2/Chapter03/truncomp.jl

max_forks_repo_path: Module 2/Chapter03/truncomp.jl


text:
function vech(A::AbstractMatrix{T}) where {T}
    m = LinearAlgebra.checksquare(A)
    v = Vector{T}(undef, (m * (m + 1)) >> 1)
    k = 0
    for j in 1:m, i in j:m
        @inbounds v[k += 1] = A[i, j]
    end
    return v
end

function inv_vech(v::AbstractVector{T},
                  n = Int(round(1 / 2 * (sqrt(1 + 8 * length(v) - 1))))) where {T}
    n * (n + 1) / 2 == length(v) || error("length(v) != n(n+1)/2")
    A = LowerTriangular(zeros(T, n, n))
    indices = [0; cumsum(n:-1:1)]
    for j in 1:n, i in j:n
        @inbounds A[i, j] = v[indices[j] + i - (j - 1)]
    end
    return A
end

# inefficient comparision for unit testing.  Not for high performance code
# Does not require the same type
function all_fields_equal(x1::T1, x2::T2, fields) where {T1,T2}
    return all(getfield.(Ref(x1), fields) .== getfield.(Ref(x2), fields))
end

# Helpers for generation

# conditionally call to support functions = nothing
maybe_call_function(f, args...) = f(args...)
maybe_call_function(::Nothing, args...) = nothing

# Rethrow if an exception is non-LAPACK -- we shouldn't fail silently for non-LAPACK errors
# The reason that LAPACK errors are simple rethrown is that they typically mean model is in
# a non-convergence region and new parameters should be drawn
function is_linear_algebra_exception(e)
    if e isa LAPACKException
        return true
    elseif e isa PosDefException
        return true
    else
        return false
    end
end

# See https://github.com/JuliaLang/julia/issues/33409
# Default even lower?
function set_blas_threads(openblas_threads = min(4, Int64(round(Sys.CPU_THREADS / 2))))
    if (BLAS.vendor() == :openblas64)
        println("Setting openblas64 threads = $openblas_threads")
        BLAS.set_num_threads(openblas_threads)
    end
end


text:
"""
This function adds the power limits of generators when there are no CommitmentVariables
"""
function activepower(m::JuMP.AbstractModel, devices::Array{L,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64) where {L <: PowerSystems.ElectricLoad, D <: AbstractControllableLoadForm, S <: PM.AbstractPowerFormulation}

    p_cl = m[:p_cl]
    time_index = m[:p_cl].axes[2]
    name_index = m[:p_cl].axes[1]

    (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent"): true

    pmax_cl = JuMP.JuMPArray(Array{ConstraintRef}(undef, length.(JuMP.axes(p_cl))), name_index, time_index)
    for t in time_index, (ix, name) in enumerate(name_index)
        if name == devices[ix].name
            pmax_cl[name, t] = @constraint(m, p_cl[name, t] <= devices[ix].maxactivepower * values(devices[ix].scalingfactor)[t])
        else
            @error "Bus name in Array and variable do not match"
        end

    end

    JuMP.register_object(m, :loadcontrol_activelimit, pmax_cl)

end

function reactivepower(m::JuMP.AbstractModel, devices::Array{L,1}, device_formulation::Type{D}, system_formulation::Type{S}, time_periods::Int64) where {L <: PowerSystems.ElectricLoad, D <: AbstractControllableLoadForm, S <: AbstractACPowerModel}

    q_cl = m[:q_cl]
    p_cl = m[:p_cl]
    time_index = m[:q_cl].axes[2]
    name_index = m[:q_cl].axes[1]

    (length(time_index) != time_periods) ? @error("Length of time dimension inconsistent") : true

    qmax_cl = JuMP.JuMPArray(Array{ConstraintRef}(undef, length.(JuMP.axes(q_cl))), name_index, time_index)

    for t in time_index, (ix, name) in enumerate(name_index)
        if name == devices[ix].name
            qmax_cl[name, t] = @constraint(m, q_cl[name, t] == p_cl[name, t]*sin(atan((devices[ix].maxactivepower/devices[ix].maxreactivepower))))
        else
            @error "Bus name in Array and variable do not match"
        end
    end

    JuMP.register_object(m, :loadcontrol_reactivelimit, qmax_cl)

end


text:
using LinearAlgebra

include("./kernel.jl")
include("./transform.jl")
include("./incremental.jl")
include("./kernel_system.jl")

@doc raw"""
    Data

A struct containing at most `n` data points stored in column major order. `X` contains data locations, `F` contains data covariates, `Y` contains observations. The data struct may store at most `capacity` points.
"""
mutable struct Data
    capacity::Int
    n::Int
    x::IncrementalColumns{Float64} # Data locations
    fx::IncrementalColumns{Float64} # Covariates
    y::IncrementalVector{Float64} # Observations
end
function Data(capacity, x, fx, y)
    @assert size(x, 2) == size(fx, 2) == size(y, 1)
    n = size(x, 2)
    new_x = incremental_columns(capacity, x)
    new_fx = incremental_columns(capacity, fx)
    new_y = incremental_vector(capacity, y)
    return Data(capacity, n, new_x, new_fx, new_y)
end
function Data(x, fx, y)
    n = size(x, 2)
    return Data(n, x, fx, y)
end
function Data(capacity, d::Int, p::Int)
    new_x = incremental_columns(Float64, capacity, d)
    new_fx = incremental_columns(Float64, capacity, p)
    new_y = incremental_vector(Float64, capacity)
    return Data(capacity, 0, new_x, new_fx, new_y)
end
function add_point!(d::Data, x, fx, y)
    @assert d.n < d.capacity
    add_col!(d.x, x)
    add_col!(d.fx, fx)
    add_element!(d.y, y)
    d.n += 1
    return nothing
end

@doc raw"""
"""
mutable struct ComputeData
    capacity::Int
    n::Int
    d1::Data # unisolvent points
    d2::Data # Reduced dataset - see sjtu notes
    f1::LU{Float64, Matrix{Float64}} # LU factorization of covariates
    w::IncrementalColumns{Float64}
end
function ComputeData(capacity, d1, x, fx, y)
    d2 = Data(capacity, x, fx, y)
    f1 = lu(d1.fx)
    w = incremental_columns(capacity, f1 \ fx)
    return ComputeData(capacity, 0, d1, d2, f1, w)
end

function add_point!(comp::ComputeData, x, fx, y)
    @assert comp.n < comp.capacity
    comp.n += 1
    add_point!(comp.d2, x, fx, y)
    add_col!(comp.w, comp.f1 \ fx) # Add new row to W
    return nothing
end

mutable struct ComputeKernel{K<:AbstractCorrelation,T}
    capacity::Int
    n::Int
    k::K
    θ::T
    k1::Matrix{Float64}
    k12::IncrementalColumns{Float64}
    k2_tilde::IncrementalCholesky{Float64}
end
function ComputeKernel(capacity, comp::ComputeData, k, θ, k0)
    p = size(d1.fx, 1)
    k1 = Array{Float64}(undef, p, p)
    pairwise!(k1, k, d1.x, d1.x, θ...)
    k12 = Array{Float64}(undef, p, capacity)
    k2 = Array{Float64}(undef, capacity, capacity)
    return ComputeKernel(capacity, 0, k, θ, k1, k12, k2)
end

function add_point(c::ComputeKernel, comp::ComputeData, x)
    n = comp.n
    colwise!(c.k12[:, n], c.k, comp.d1.x[:, 1:comp.d1.n], x, c.θ...)
    # Update K2
    colwise!(c.k2[1:n, n], c.k, comp.d2.x[:, 1:comp.d2.n], x, c.θ...)
    # Update Ktilde
    tmp = c.k1 * comp.w[:, n]
    tmp2 = c.k12[:, 1:n]' * comp.w[:, n]
    c.k2[n, n] +=  dot(comp.w[:, n], tmp) - 2 * tmp2[n]
    c.k2[1:n-1, n] .+= (comp.w[:, 1:n-1]' * tmp) .- tmp2[1:n-1]
    # Update Cholesky factor of Ktilde
    # TODO write this more elegantly, this is messy
    v
    add_col!(c.k2_tilde, v)
    return nothing
end

@doc """
"""
struct BTG{G<:AbstractTransform,K<:AbstractCorrelation,C}
    horizon::Int
    n::Int
    g::G
    k::K
    comp::ComputeData
    comp_k::Vector{C}
end
function BTG()
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule NativeFileDialog_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("NativeFileDialog")
JLLWrappers.@generate_main_file("NativeFileDialog", UUID("94d9ae2c-efc7-56f8-9a02-54c47b797961"))
end  # module NativeFileDialog_jll


text:
@testset "argmax" begin
    arg, val = GeometryTypes.argmax(sin, range(0, stop=π, length=101))
    @test arg ≈ π/2
    @test val ≈ 1

    arg, val = GeometryTypes.argmax(cos, range(0, stop=π, length=101))
    @test arg ≈ 0  # 0, not π, because argmax requires a subsequent value to be *strictly better* than the current best
    @test val ≈ 1.0

    @test @allocated(GeometryTypes.argmax(sin, range(0, stop=π, length=101))) == 0

    @test_throws ArgumentError GeometryTypes.argmax(identity, [])
end

@testset "extrema" begin
    xs = [GeometryTypes.Vec2f0(1f0, 1f0), GeometryTypes.Vec2f0(0.5f0, 1.5f0)]
    extr = extrema(xs)

    @test extr[1] ≈ GeometryTypes.Vec2f0(0.5f0, 1f0)
    @test extr[2] ≈ GeometryTypes.Vec2f0(1f0, 1.5f0)

    @testset "minimum" begin
        val = minimum(xs)
        @test val ≈ extr[1]
    end

    @testset "maximum" begin
        val = maximum(xs)
        @test val ≈ extr[2]
    end
end


text:
"""
    color_from_palette(colordescription::String) -> RGB

Pick a similar color to str from within PALETTE. Use for sethue, setbackground
Look for color descriptions in http://juliagraphics.github.io/Colors.jl/stable/namedcolors/#
"""
color_from_palette(colordescription::String) = get(PALETTE, getinverse(PALETTE, parse(RGB, colordescription)))

"The L in HSL is not a good indicator for distinguishing between perceived lightness.
This may provide a better parameter for finding contrasting colors"
function lumin2(colo)
    @warn "lumin2 deprecated"
    c = RGB(colo)
    r, g, b = c.r, c.g, c.b
    γ = 2.2
    0.2126 * r^γ + 0.7152 * g^γ + 0.0722 * b^γ
end

#=
luminance(colo) = begin
    @warn "luminance deprecated"
    HSL(colo).l
end

"""
color_with_luminance(color::T, luminance::Float64) where T
        -> T
"""
function color_with_luminance(color::T, luminance::Float64) where T
    @warn "color_with_luminance deprecated"
    c = HSLA(color)
    T( HSLA(c.h, c.s, luminance, c.alpha)  )
end
=#

get_current_RGB() = RGB(get_current_redvalue(), get_current_greenvalue(), get_current_bluevalue())
get_current_RGBA() = RGBA(get_current_redvalue(), get_current_greenvalue(), get_current_bluevalue(), get_current_alpha())

get_current_lumin2() = begin
    @warn "lumin2 deprecated"
    lumin2(get_current_RGB())
end
get_current_lumin() = lumin(get_current_RGB())

"""
color_with_lumin2(color::T, luminance::Float64) where T
        -> T
The name of this function is somewhat illogical, consider removing. Generally
use the LCHuv colorspace to modify luminance of a color.
"""
function color_with_lumin2(color::T, lumin2::Float64) where T
    c = RGBA(color)
    r, g, b, a = c.r, c.g, c.b, c.alpha
    γ = 2.2
    lr, lg, lb = 0.2126 * r^γ , 0.7152 * g^γ , 0.0722 * b^γ
    l = lr + lg + lb
    lr1, lg1, lb1 = (lumin2 / l ) .* (lr, lg, lb)
    r1, g1, b1 = (0.2126 / lr1)^(-1/γ), (0.7152/lg1)^(-1/γ), (0.0722/lb1)^(-1/γ)
    T(RGBA(r1, g1, b1, a))
end

"""
    color_with_lumin(col, lumin)
Convert the color to have 'luminosity' or similar term, depending on the type of color.
"""
function color_with_lumin(c::C, lumin) where C<:Union{LCHab, LCHuv}
    C(lumin, comp2(c), comp3(c))
end
function color_with_lumin(c::C, lumin) where C<:Union{HSV, HSL}
    C(comp1(c), comp2(c), lumin)
end
function color_with_lumin(c::C, lumin) where C<:Union{LCHabA, LCHuvA}
    C(lumin, comp2(c), comp3(c), comp4(c))
end
function color_with_lumin(c::C, lumin) where C<:Union{HSVA, HSLA}
    C(comp1(c), comp2(c), lumin, comp4(c))
end
function color_with_lumin(c::RGB, lumin)
    convert(RGB, color_with_lumin(convert(LCHuv ,c), lumin))
end
function color_with_lumin(c::RGBA, lumin)
    convert(RGBA, color_with_lumin(convert(LCHuvA ,c), lumin))
end
"""
    lumin(c)

Return luminosity or similar term, depending on color type
"""
function lumin(c::C) where C<:Union{LCHab, LCHuv, LCHabA, LCHuvA}
    comp1(c)
end
function lumin(c::C) where C<:Union{HSV, HSL, HSVA, HSLA}
    comp3(c)
end
function lumin(c::RGB)
    lumin(convert(LCHuv, c))
end
function lumin(c::RGBA)
    lumin(convert(LCHuvA, c))
end

"""
    color_with_hue(col, hue::Angle)
Convert the color to have cylindrical hue.
"""
function color_with_hue(c::C, hue::Angle) where C<:Union{LCHab, LCHuv}
    C(comp1(c), comp2(c), typeof(comp1(c))(hue / °))
end
function color_with_hue(c::C, hue::Angle) where C<:Union{LCHabA, LCHuvA}
    C(comp1(c), comp2(c), typeof(comp1(c))(hue / °), comp4(c))
end

function color_with_hue(c::C, hue::Angle) where C<:Union{HSV, HSL}
    C(typeof(comp1(c))(hue / °), comp2(c), comp3(c))
end
function color_with_hue(c::C, hue::Angle) where C<:Union{HSVA, HSLA}
    C(typeof(comp1(c))(hue / °), comp2(c), comp3(c),  comp4(c))
end

function color_with_hue(c::RGB, hue::Angle)
    convert(RGB, color_with_hue(convert(LCHuv, c), hue ))
end
function color_with_hue(c::RGBA, hue::Angle)
    convert(RGBA, color_with_hue(convert(LCHuvA, c), hue ))
end

"""
    hue(c)

Return angular hue or similar term, depending on color type
"""
function hue(c::C) where C<:Union{LCHab, LCHuv, LCHabA, LCHuvA}
    comp3(c)°
end
function hue(c::C) where C<:Union{HSV, HSL, HSVA, HSLA}
    comp1(c)°
end
function hue(c::RGB)
    hue(convert(LCHuv, c))
end
function hue(c::RGBA)
    hue(convert(LCHuvA, c))
end

"""
    rotate_hue(c, angle::Angle)

In cylindrical colorspace, rotation of hue gives
continuous variation while not chaninging chroma (or similar term)
or perceptual luminance (or similar term).
"""
function rotate_hue(c, ang::Angle)
    modhue = mod(hue(c) + ang , 360°)
    color_with_hue(c, modhue)
end

"""
    color_with_alpha(col, alpha)

Convert col to a color with transparency alpha.
"""
color_with_alpha(col, alpha) = Luxor.coloralpha(col, alpha)

"""
    color_without_alpha(c)

Convert col to a color type without transparency.
"""
color_without_alpha(c::HSVA) = convert(HSV, c)
color_without_alpha(c::HSV) = c

color_without_alpha(c::HSLA) = convert(HSL, c)
color_without_alpha(c::HSL) = c

color_without_alpha(c::LCHabA) = convert(LCHab, c)
color_without_alpha(c::LCHab) = c

color_without_alpha(c::LCHuvA) = convert(LCHuv, c)
color_without_alpha(c::LCHuv) = c

color_without_alpha(c::RGBA) = convert(RGB, c)
color_without_alpha(c::RGB) = c

"""
    chroma(c)

Return chroma or similar term, depending on color type
"""
function chroma(c::C) where C<:Union{LCHab, LCHuv, LCHabA, LCHuvA}
    comp2(c)
end
function chroma(c::C) where C<:Union{HSV, HSL, HSVA, HSLA, RGB, RGBA}
    chroma(convert(LCHuv, c))
end

"""
    saturation(c)

Return chroma or similar term, depending on color type
"""
function saturation(c::C) where C<:Union{LCHab, LCHuv, LCHabA, LCHuvA, RGB, RGBA}
    saturation(convert(HSV, c))
end
function saturation(c::C) where C<:Union{HSV, HSL, HSVA, HSLA}
    comp2(c)
end

"""
    color_with_saturation(col, saturation)
Modify saturation.
"""
function color_with_saturation(c::C, saturation) where C<:Union{LCHab, LCHuv, RGB}
    convert(C, color_with_saturation(convert(HSV, c), saturation ))
end
function color_with_saturation(c::C, saturation) where C<:Union{LCHabA, LCHuvA, RGBA}
    convert(C, color_with_saturation(convert(HSVA, c), saturation ))
end
function color_with_saturation(c::C, saturation) where C<:Union{HSV, HSL}
    C(comp1(c), saturation, comp3(c))
end
function color_with_saturation(c::C, saturation) where C<:Union{HSVA, HSLA}
    C(comp1(c), saturation, comp3(c), comp4(c))
end

# This definition don't really belong anywhere. ColorSchemes is mostly a complete library, but this is useful and missing.
const ColSchemeNoMiddle = ColorSchemes.ColorScheme(SVector([[get(PuOr_8, x, :clamp) for x in range(0, 0.3, length = 10)];
    [get(PuOr_8, 1 - x, :clamp) for x in range(0.3, 0, length = 11)]]...), "MechanicalSketch", "emphasize plus / minus")


text:
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(@__DIR__)
@everywhere using MyMCMCSampler
@everywhere using StatsBase
@everywhere using Combinatorics
@everywhere using Random
@everywhere using PyCall
using Distributions
using DataFrames
using ArgParse
using JLD2, FileIO
using ProgressMeter
using Base.Iterators: enumerate
using UUIDs
using StatsPlots
using BenchmarkTools
@everywhere using Base: @nloops, @ntuple

py"""
import deepdish as dd
def savefile(filepath,data):
    dd.io.save(filepath,data)
"""
rng = MersenneTwister(1234);

@everywhere extensions = quote
    py"""
    import numpy as np
    import networkx as nx
    from networkx.utils import py_random_state
    def compute_count(channel, group):
        divide = channel//group
        remain = channel%group

        out = np.zeros(group, dtype=int)
        out[:remain]=divide+1
        out[remain:]=divide
        return out

    @py_random_state(3)
    def ws_graph(n, k, p, seed=1):
        assert k>=2 and k<=n
        # compute number of edges:
        edge_num = int(round(k*n/2))
        count = compute_count(edge_num, n)
        # print(count)
        G = nx.Graph()
        for i in range(n):
            source = [i]*count[i]
            target = range(i+1,i+count[i]+1)
            target = [node%n for node in target]
            # print(source, target)
            G.add_edges_from(zip(source, target))
        # rewire edges from each node
        nodes = list(G.nodes())
        for i in range(n):
            u = i
            target = range(i+1,i+count[i]+1)
            target = [node%n for node in target]
            for v in target:
                if seed.random() < p:
                    w = seed.choice(nodes)
                    # Enforce no self-loops or multiple edges
                    while w == u or G.has_edge(u, w):
                        w = seed.choice(nodes)
                        if G.degree(u) >= n - 1:
                            break  # skip this rewiring
                    else:
                        G.remove_edge(u, v)
                        G.add_edge(u, w)
        return G

    @py_random_state(4)
    def connected_ws_graph(n, k, p, tries=100, seed=1):
        for i in range(tries):
            # seed is an RNG so should change sequence each call
            G = ws_graph(n, k, p, seed)
            if nx.is_connected(G):
                return G
        raise nx.NetworkXError('Maximum number of tries exceeded')

    def compute_stats(G):
        G_cluster = sorted(list(nx.clustering(G).values()))
        cluster = sum(G_cluster) / len(G_cluster)
        # path
        path = nx.average_shortest_path_length(G)
        return cluster, path

    def norm(x):
        ## to 0-1
        x = x-x.min()
        if x.max()==0:
            x[:]=0
        else:
            x = x/x.max()
        return x

    def get_graph(n,degree,p,repeat=30):
        # row: repeat, column: n,degree,p,seed,cluster,path
        result = np.zeros((repeat, 6))
        for i in range(repeat):
            graph = connected_ws_graph(n=n, k=degree, p=p, seed=i)
            cluster, path = compute_stats(graph)
            result[i] = [n, degree, p, i, cluster, path]
        return result

    def is_connected(adj):
        n = adj.shape[0]
        c = np.arange(n)
        for i in range(n):
            for j in range(i+1,n):
                if adj[i,j]:
                    c[j] = c[i]
        if len(set(c)) > 1:
            return False
        else:
            return True
    """
    end

@everywhere eval(extensions)

### Generate training dataset
s = ArgParseSettings()
@add_arg_table s begin
    "--outdir"
        help = "folder to store dataset file"
        arg_type = String
        default = "data/tmp2"
    "--gamma"
        help = "scaling factor for the potential"
        arg_type = Float64
        default = 1.0
    "--ndata"
        help = "number of data points"
        arg_type = Int
        default = 100
    "--n"
        help = "size of the graph"
        arg_type = Int
        default = 10
end
args = parse_args(ARGS, s)
println(args)
n = args["n"]
ndata = args["ndata"]
γ = args["gamma"]

spin_mode=:spinglass
ks= range(2.0, stop=floor((n-1)*(n-2)/6),length=30)
ps= range(0.0, stop=1.0, length=30)
spin_mode=:spinglass

results = @showprogress pmap(1:ndata) do i
    k = rand(Uniform(2.0, floor((n-1)*(n-2)/6)))
    p = rand(Uniform(0.0, 1.0))
    isconnected = false
    d = nothing
    while !isconnected
        d = random_ws_factor_graph(n,0,k,p,γ=γ,mode=spin_mode,degrees=[3], singleton_dist=:randn)
        A = incident_matrix(d)
        A = (A' * A .>0)
        isconnected = py"is_connected"(A)
    end
    mP = computemarginal(d)[1]
    mP = reshape(mP,:,1)
    return (d,dist2dict(d),mP)
end

# save to file
data2file = map(results) do res
    res[2:end]  # (dict, marginal)
end
outdir = joinpath(args["outdir"], "raw")
mkpath(outdir)
outfilename = joinpath(outdir, "data.dd")
println("Output to: $(outfilename)")
py"savefile"(outfilename,data2file)


text:
using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ImplicitFreeSurface
using Statistics

function geostrophic_adjustment_simulation(free_surface)

    grid = RectilinearGrid(size = (64, 1, 1),
                                  x = (0, 1000kilometers), y = (0, 1), z = (-400meters, 0),
                                  topology = (Bounded, Periodic, Bounded))

    coriolis = FPlane(f=1e-4)

    model = HydrostaticFreeSurfaceModel(grid = grid,
                                        coriolis = coriolis,
                                        free_surface = free_surface)

    gaussian(x, L) = exp(-x^2 / 2L^2)

    U = 0.1 # geostrophic velocity
    L = grid.Lx / 40 # gaussian width
    x₀ = grid.Lx / 4 # gaussian center

    vᵍ(x, y, z) = - U * (x - x₀) / L * gaussian(x - x₀, L)

    g = model.free_surface.gravitational_acceleration

    η₀ = coriolis.f * U * L / g # geostrohpic free surface amplitude

    ηᵍ(x) = η₀ * gaussian(x - x₀, L)

    ηⁱ(x, y) = 2 * ηᵍ(x)

    set!(model, v=vᵍ, η=ηⁱ)
    gravity_wave_speed = sqrt(g * grid.Lz) # hydrostatic (shallow water) gravity wave speed
    wave_propagation_time_scale = model.grid.Δxᶜᵃᵃ / gravity_wave_speed
    simulation = Simulation(model, Δt=2wave_propagation_time_scale, stop_iteration=2)

    return simulation
end

function run_and_analyze(simulation)
    η = simulation.model.free_surface.η
    u, v, w = simulation.model.velocities

    ηx = ComputedField(∂x(η))
    compute!(ηx)

    u₀ = interior(u)[:, 1, 1]
    v₀ = interior(v)[:, 1, 1]
    η₀ = interior(η)[:, 1, 1]
    ηx₀ = interior(ηx)[:, 1, 1]

    run!(simulation)

    compute!(ηx)

    u₁ = interior(u)[:, 1, 1]
    v₁ = interior(v)[:, 1, 1]
    η₁ = interior(η)[:, 1, 1]
    ηx₁ = interior(ηx)[:, 1, 1]

    @show mean(η₀)
    @show mean(η₁)

    Δη = η₁ .- η₀

    return (; η₀, η₁, Δη, ηx₀, ηx₁, u₀, u₁, v₀, v₁)
end

fft_based_free_surface = ImplicitFreeSurface(solver_method=:FFTBased)
pcg_free_surface = ImplicitFreeSurface(solver_method=:PreconditionedConjugateGradient)

free_surfaces = [fft_based_free_surface, pcg_free_surface]
simulations = [geostrophic_adjustment_simulation(free_surface) for free_surface in free_surfaces]
data = [run_and_analyze(sim) for sim in simulations]

fft_data = data[1]
pcg_data = data[2]

#pcg_p = plot([pcg_data.η₀ pcg_data.η₁ pcg_data.Δη], label=["η₀" "ηᵢ" "Δη"], linewidth=2)
#fft_p = plot([fft_data.η₀ fft_data.η₁ fft_data.Δη], label=["η₀" "ηᵢ" "Δη"], linewidth=2)

pcg_p_η = plot([pcg_data.η₀ pcg_data.η₁], label=["η₀" "ηᵢ"], linewidth=2)
fft_p_η = plot([fft_data.η₀ fft_data.η₁], label=["η₀" "ηᵢ"], linewidth=2)

pcg_p_ηx = plot([pcg_data.ηx₀ pcg_data.ηx₁], label=["ηx₀" "ηxᵢ"], linewidth=2)
fft_p_ηx = plot([fft_data.ηx₀ fft_data.ηx₁], label=["ηx₀" "ηxᵢ"], linewidth=2)

pcg_p_u = plot([pcg_data.u₀ pcg_data.u₁], label=["u₀" "uᵢ"], linewidth=2)
fft_p_u = plot([fft_data.u₀ fft_data.u₁], label=["u₀" "uᵢ"], linewidth=2)

p = plot(pcg_p_η, fft_p_η,
         pcg_p_u, fft_p_u,
         pcg_p_ηx, fft_p_ηx,
         layout=(3, 2), titles = ["PCG η" "FFT η" "PCG u" "FFT u" "PCG η_x" "FFT η_x"])

display(p)


text:
"""
Se charge de l'exécution du branch and bound
"""

# retourne la solution du problème, faux si le problème est impossible
function branchAndBound(L::Vector{Elt},A::Matrix{Float64},epsilon;debug=Inf, verbose = false)
    # initialisation de l'arbre
    r = Elt(1,R,2,false,false,0.,"r",0.) # création de l'élément racine
    passedElts = Vector{Elt}(undef,0)
    passedDT = Vector{Float64}(undef,0)
    h = length(L)
    nbElts = length(L)
    dT = 0.
    lowerBound = getLowerBound(L)
    upperBound = Inf
    # Creation de la racine de l'arbre
    root = newNode(passedElts,passedDT,L,lowerBound,r,h,dT,lowerBound)
    # descente récursive dans l'arbre
    finished = false
    lastBestNode = 0
    while !finished && debug > 0
        debug -= 1
        # On calcule le meilleur noeud à explorer de l'arbre. Si le noeud retourné vaux faux c'est qu'il n'y a plus de noeuds à explorer
        tuple = bestNodeToExplore(root,upperBound, verbose=verbose)
        if tuple[1] == false # plus de noeuds à explorer
            finished = true # on s'arrête
        else # on va explorer le noeud et refaire un tour
            aff(root)
            verbose && println("")
            enBas = exploreNode(tuple[1],upperBound,A,nbElts,epsilon, verbose=verbose) # on explore le meilleur noeud à explorer
            # enBas est sensé valoir true sauf si on a atteint le bas de l'arbre
            if enBas != true # on a atteint le bas de l'arbre
                if enBas[2] < upperBound # si la valeur du noeud en bas est meilleure que la borne sup actuelle
                    upperBound = enBas[2]
                    lastBestNode = enBas[1]
                    verbose && println("MEILLEUR NOEUD FINAL -> ",enBas[1].e.name)
                end
            end
        end
    end
    aff(root)
    verbose && println("upperBound = ",upperBound)
    if lastBestNode != 0
    	return lastBestNode.passedElts, lastBestNode.e
    else
    	return false, false
    end
end

# retourne le meilleur noeud à explorer ainsi que sa lowerBound, si le return vaut false,false c'est qu'aucun noeud ne reste à explorer
function bestNodeToExplore(node::Node, upperBound::Float64; str="", verbose = false)
    verbose && println("\n",str,"bestNodeToExplore[",node.e.name,"]")
    if node.h == 0 # si le noeud est une feuille de l'arbre, il est déjà sondé. On ne doit pas le reregarder
        verbose && println("On est en bas de l'arbre ce noeud n'est pas valable")
        return node,Inf
    elseif length(node.children) == 0 # si le noeuf n'a pas d'enfants alors il est candidat du meilleur noeud
        verbose && println(str,"Ce noeud n'a pas d'enfants, on remonte")
        return node,node.lowerBound
    else # sinon c'est un noeuf déjà étudié qui possède des enfants. On va descendre dans ceux ci
        verbose && println(str,"Ce noeud a des enfants et n'est pas tout en bas")
        # on calcule l'enfant qui possède la plus petite lowerBound
        nodeMin = 0
        min = Inf
        str = string(str,"      ")
        for indexChild in 1:length(node.children) # on itère sur les enfants
            child = node.children[indexChild] # récupération de l'enfant
            if child.lowerBound > upperBound # si l'enfant possède une borne inf plus grande que la borne sup connue
                verbose && println(str,"L'enfant est trop nul on le supprime")
                node.children[indexChild] = emptyNode() # on le supprime
            elseif child.empty # l'enfant est déjà supprimé
                verbose && println(str,"Ce noeud a été supprimé")
            else # si les enfants sont encore vivants
                verbose && println(str,"L'enfant est valide, on appelle bestNodeToExplore sur lui")
                tuple = bestNodeToExplore(child, upperBound,str=str, verbose=verbose) # on calcule le meilleur noeuf récursivement
                if tuple[2] < min # si le noeud a une meilleure borne inf que celle connue
                    nodeMin = tuple[1]
                    min = tuple[2]
                end
            end
        end
        if nodeMin == 0 # on n'a trouvé aucun noeud potentiel. On retourne faux faux.
            return false,false
        else
            return nodeMin,min
        end
    end
end

# construit les enfants du noeud passé en paramètres
function exploreNode(node::Node, upperBound::Float64, A::Matrix{Float64}, nbElts::Int, epsilon::Float64; verbose = false)
    verbose && println("exploreNode[",node.e.name," - (",node.dT,",",node.lowerBound,")]")
    node.children = Vector{Node}(undef,node.h) # on alloue le vecteur des enfants de ce noeud
    newUpperBound = Inf
    newOptNode = 0
    for indexChild in 1:length(node.futureElts) # parcours des éléments qui vont devenir les enfants de ce noeud
        e = node.futureElts[indexChild] # élément courant
        # calcul de sa lowerBound
        d = A[node.e.id,e.id] # distance entre le noeud courant et l'enfant
        dT = node.dT + d # distance parcourue depuis la racine
        sumFutureCMin = node.sumFutureCMin - e.cMin
        lowerBound = dT + sumFutureCMin # borne inf de l'enfant
        # vérification des contraintes w et epsilon
        if ((dT > e.limit) && (e.limit != -1.)) || lowerBound > upperBound # si l'enfant est au mieux plus mauvais que la borne sup connue
            node.children[indexChild] = emptyNode()
        else # l'enfant a une chance d'être meilleur que la borne sup connue
            # construction des attributs du nouvel enfant
            nbFutureElts = node.h - 1 # nombre d'elts restants de l'enfant
            nbPassedElts = nbElts - 1 - nbFutureElts # nombre d'elts précédents de l'enfant
            passedElts = newPassedElts(node.passedElts,node.e)
            # on vérifie que l'élément est autorisé (pas de e1 avant s1 typiquement)
            authorized = isOrderRespected(node.futureElts,nbFutureElts+1,e)
            if authorized # le noeud a le droit d'être créé
                timeOK = true
                if e.state == E && e.sourceRemaining
                    distance = 0
                    source = 0
                    if node.e.state == S && node.e.idReq == e.idReq # la source que l'on cherche est juste au dessus de notre noeud
                        source = node.e
                        distance = A[node.e.id,e.id]
                    else # source est contenue dans la liste des elts passés de node
                        sourceFound = false
                        indexElt = 1
                        while !sourceFound && indexElt <= length(node.passedElts) # on parcours la liste des noeuds placés
                            eltBuffer = node.passedElts[indexElt]
                            dTBuffer = node.passedDT[indexElt]
                            if eltBuffer.state == S && eltBuffer.idReq == e.idReq # c'est la source !
                                distance = node.dT + A[node.e.id,e.id] - dTBuffer
                                sourceFound = true
                                source = eltBuffer
                            end
                            indexElt += 1
                        end
                        if !sourceFound
                            verbose && println("ON A PAS TROUVE LA SOURCE WTF CEST PAS NORMAL !!!")
                            exit()
                        end
                    end
                    if distance > (1+epsilon)*A[source.id,e.id] # la contrainte de trajet est respectée
                        timeOK = false
                    end
                end
                if timeOK
                    futureElts = newFutureElts(node.futureElts,nbFutureElts,indexChild) # elements suivants
                    passedDT = newPassedDT(node)
                    node.children[indexChild] = newNode(passedElts,passedDT,futureElts,sumFutureCMin,e,node.h-1,dT,lowerBound) # création de l'enfant
                    if node.h == 1
                        if lowerBound < upperBound
                            newUpperBound = lowerBound
                            newOptNode = node.children[indexChild]
                        end
                    end
                else
                    node.children[indexChild] = emptyNode()
                end
            else # le noeud n'a pas le droit d'être créé
                node.children[indexChild] = emptyNode()
            end
        end
    end
    if node.h == 1 # on est en bas
        return newOptNode,newUpperBound
    else
        return true
    end
end

# construit l'enfant indexChild du noeud root
function newChild(root::Node,indexChild::Int,A::Matrix{Float64})
    e = root.futureElts[indexChild] # élement du noeud enfant
    d = A[root.e.id,e.id] # distance entre le noeud root et le noeud enfant
    dT = root.dT + d # distance du noeud enfant à la racine
    sumFutureCMin = root.sumFutureCMin - e.cMin # somme des distances min des éléments du noeud enfant
    lowerBound = dT + sumFutureCMin # borne inf du noeud enfant
    passedElts = newPassedElts(root.passedElts,root.e)
    passedDT = newPassedDT(root)
    futureElts = newFutureElts(futureElts,length(futureElts)-1,indexChild)
    # on construit le noeud enfant
    return newNode(passedElts,passedDT,futureElts,sumFutureCMin,e,root.h-1,dT,lowerBound)
end

# construit l'attribut passedDT des fils de node
function newPassedDT(node::Node)
    list = Vector{Float64}(undef,length(node.passedDT)+1)
    for indexList in 1:length(node.passedDT)
        list[indexList] = node.passedDT[indexList]
    end
    list[length(node.passedDT)+1] = node.dT
    return list
end

# retourne la somme les cMin de chaque elt
function getLowerBound(v::Vector{Elt})
    lb = 0.
    for e in v
        lb += e.cMin
    end
    return lb
end

# construit la liste des éléments passés des fils de node
function newPassedElts(passedElts,e)
    list = Vector{Elt}(undef,length(passedElts)+1)
    for indexList in 1:length(passedElts)
        list[indexList] = passedElts[indexList]
    end
    list[length(passedElts)+1] = e
    return list
end

# construit la nouvelle liste des éléments restants en supprimant l'élément qui vient d'être utilisé dans le noeud
function newFutureElts(v::Vector{Elt},lengthV::Int,indexChild::Int) # lengthV est la longueur du nouveau tableau attention !
    list = Vector{Elt}(undef,lengthV)
    for indexV in 1:lengthV
        if indexV < indexChild
            list[indexV] = v[indexV]
        else
            list[indexV] = v[indexV+1]
        end
    end
    return list
end

# vérifie que l'ordre des éléments dans le trajet est logique (pas de e1 avant s1)
function isOrderRespected(v::Vector{Elt},lengthV::Int,e::Elt)
    if e.state == S
        return true
    elseif e.state == E
        for indexElt in 1:lengthV # parcours des éléments
            otherE = v[indexElt]
            if otherE.state == S && otherE.idReq == e.idReq
                return false
            end
        end
        return true
    end
end

# fonction d'affichage de l'arbre
function printTree(root::Node;str::String="", verbose = false)
    if root.probed
        verbose && println(str,"PROBED")
    elseif root.empty
        verbose && println(str,"DELETED")
    else
        verbose && println(str,"noeud [",root.e.name," - (",root.dT,",",root.lowerBound,")]")
        str = string(str,"      ")
        for child in root.children
            verbose && println("[",root.e.name,"]")
            verbose && printTree(child,str=str, verbose=verbose)
        end
    end
end

# affiche l'arbre des plannings possibles
function aff(root::Node; verbose = false)
    verbose && println("\nprintTree :")
    verbose && printTree(root, verbose=verbose)
    verbose && println("")
end


text:
struct Pipelines <: AbstractVisualization
    pipelines::AbstractVector
end

function layout(g::SimpleDiGraph)
    xs, ys, _ = solve_positions(Zarate(), g)
    return Point.(ys, -xs)
end

function jsrender(session::Session, pipelines::Pipelines)
    arrow_size = 25
    edge_width = 5
    node_size = 45
    legend_node_size = 20
    colgap = 150

    # set general legend
    palette = vcat(RGB(colorant"black"), Makie.current_default_theme().palette.color[])
    un = reduce(vcat, get_vertex_names.(pipelines.pipelines))
    uc = palette[eachindex(un)]
    scale = AlgebraOfGraphics.CategoricalScale(un, uc, palette, "Step")

    nested_vertices = lift(output(last(pipelines.pipelines))) do _
        return get_vertices.(pipelines.pipelines)
    end
    g = @lift simpledigraph($nested_vertices)

    names = @lift mapreduce(get_vertex_names, append!, $nested_vertices)
    node_color = @lift AlgebraOfGraphics.rescale($names, scale)
    # work around https://github.com/JuliaPlots/GraphMakie.jl/issues/42
    edge_color = @lift ne($g) > 0 ? :black : :transparent
    backgroundcolor = colorant"#F3F4F6"
    width, height = Observable(350), Observable(350)
    fig = Figure(; backgroundcolor)
    ax = Axis(fig[1, 1]; width, height, backgroundcolor)

    on(g) do graph
        points = layout(graph)
        xlims = extrema(first, points)
        ylims = extrema(last, points)
        xlims!(ax, xlims[1] - 0.25,  xlims[2] + 0.25)
        ylims!(ax, ylims[1] - 0.25,  ylims[2] + 0.25)
    end
    notify!(g)
    graphplot!(ax, g;
        arrow_show=true, arrow_size,
        edge_width, edge_color,
        node_size, node_color,
        layout
    )
    hidedecorations!(ax)
    hidespines!(ax)

    Legend(
        fig[1, 2],
        [MarkerElement(; color, marker=:circle, markersize=legend_node_size) for color in scale.plot],
        string.(scale.data), scale.label
    )
    colgap!(fig.layout, colgap)
    on_pixelratio(session, once=true) do pr
        width[] = round(Int, 350pr)
        height[] = round(Int, 350pr)
        AlgebraOfGraphics.resizetocontent!(fig)
    end

    return jsrender(session, scrollable_component(fig))
end


text:
type MOGAFitnessScheme{T} <: ParetoFitnessScheme{T}
  maximise::Vector{Bool}
end

"""
TODO: Describe MOGA.
"""
moga(s::Dict{Any, Any}) =
  MOGAFitnessScheme{eval(Base.parse(Base.get(s, "of", "Float")))}(s["maximise"])

scale!{F}(s::MOGAFitnessScheme, inds::Vector{Tuple{Int, F}}) =
  for (p1_id, p1_f) in inds
    p1_f.rank = 1 + count((p2_id, p2_f) -> dominates(p2_f.scores, p1_f.scores, s.maximise), inds)
  end


text:
module MacroDocs

if VERSION < v"0.4-dev+4319"
    tup(args...) = tuple(args...)
else
    tup(args...) = Tuple{args...}
end

import Docile: @doc, @docstrings, @doc_str, @doc_mstr, meta

using Compat; import Compat.String

@docstrings(manual = ["../../docs/manual.md"])

include(joinpath("MacroDocs", "macro-spec.jl"))

include(joinpath("MacroDocs", "functions.jl"))
include(joinpath("MacroDocs", "globals.jl"))
include(joinpath("MacroDocs", "macros.jl"))
include(joinpath("MacroDocs", "types.jl"))

end


text:
using qccdSimulator.QCCDevControl_Types
using qccdSimulator.QCCDevDes_Types
using Random
using qccdSimulator.QCCDDevControl
using qccdSimulator.QCCDev_Utils

"""
Generates n junctions connected to ZoneInfoDesc (auxZones or gateZones).
repJunc: Repeats a junction ID.
wrongJunctType: Gives a wrong junction type to a ZoneInfoDesc.
isolatedJunc: The first junction is not connected to any ZoneInfoDesc.
repJunc: Repeats a zone ID.
"""
function giveZonesJunctions(nJunctions:: Int64, juncTypes:: Array{String};
            repJunc=false, wrongJuncType=false, isolatedJunc=false, repZone=false)::
            Tuple{Array{ZoneInfoDesc},Array{JunctionInfoDesc}}

    zones = ZoneInfoDesc[]
    junctions = JunctionInfoDesc[]
    sId = 0
    skipAuxZone = wrongJuncType
    for i in 1:nJunctions
        repJunc && i > 1 ? push!(junctions, JunctionInfoDesc(string(i-1), juncTypes[i])) :
        push!(junctions, JunctionInfoDesc(string(i), juncTypes[i]))
        if isolatedJunc
            isolatedJunc = false
            continue
        end
        for j in 1:typesSizes[Symbol(juncTypes[i])]
            if skipAuxZone
                skipAuxZone = false
                continue
            end
            push!(zones, ZoneInfoDesc(string(sId),string(i),string(-1),3))
            if repZone
                repZone = false
                continue
            end
            sId += 1
        end
    end
    return zones, junctions
end

"""
Generates Qubits.
    * qubitID::Int ID of the qubit, if not given it'll be asigned randomly
    * zone::Symbol Desaired zone to create the ion
returns:
    * Qubit struct
"""
function giveQubit(zone::Symbol, qubitID::Int = nothing)
    return qubitID == nothing ? Qubit(rand(), zone) : Qubit(qubitID, zone)
end

"""
Creates some ZoneInfoDesc objects.
"""
function giveZoneInfo(nZones:: Int64;  invZone=false, giveNothing=false)::Array{ZoneInfoDesc}
    zones = ZoneInfoDesc[]
    for i in 1:nZones
        if invZone
            push!(zones,ZoneInfoDesc(string(i),string(i),string(i),2))
            invZone = false
        elseif giveNothing
            push!(zones,ZoneInfoDesc(string(i),"","",2))
            giveNothing = false
        else
            push!(zones,ZoneInfoDesc(string(i),string(i+1),string(i+2),2))
        end
    end
    return zones
end

"""
Creates a struct QCCDevDescription based in the file topology.json
"""
function giveQccDes()::QCCDevDescription
    gateZone:: GateZoneDesc = GateZoneDesc(
        [
            ZoneInfoDesc("1", "", "4", 2),
            ZoneInfoDesc("2", "4", "5", 2),
            ZoneInfoDesc( "3", "7", "8", 2)
        ]
    )
    auxZone:: AuxZoneDesc = AuxZoneDesc(
        [
            ZoneInfoDesc( "4", "1", "2", 2),
            ZoneInfoDesc( "5", "", "9", 2),
            ZoneInfoDesc( "6", "9", "", 2),
            ZoneInfoDesc( "7", "9", "3", 2)
        ]
    )
    junction:: JunctionDesc = JunctionDesc(
        [
            JunctionInfoDesc( "9", "T")
        ]
    )
    loadZone:: LoadZoneDesc = LoadZoneDesc(
        [
            LoadZoneInfoDesc( "8", "3", "")
        ]
    )
    return  QCCDevDescription(gateZone,auxZone,junction,loadZone)
end

function giveQccDes2()::QCCDevDescription
    gateZone:: GateZoneDesc = GateZoneDesc(
        [
            ZoneInfoDesc("1", "", "4", 2),
            ZoneInfoDesc("2", "4", "5", 2),
            ZoneInfoDesc( "3", "7", "8", 8)
        ]
    )
    auxZone:: AuxZoneDesc = AuxZoneDesc(
        [
            ZoneInfoDesc( "4", "1", "2", 2),
            ZoneInfoDesc( "5", "", "9", 2),
            ZoneInfoDesc( "6", "9", "", 2),
            ZoneInfoDesc( "7", "9", "3", 2)
        ]
    )
    junction:: JunctionDesc = JunctionDesc(
        [
            JunctionInfoDesc( "9", "T")
        ]
    )
    loadZone:: LoadZoneDesc = LoadZoneDesc(
        [
            LoadZoneInfoDesc( "8", "3", "")
        ]
    )
    return  QCCDevDescription(gateZone,auxZone,junction,loadZone)
end

"""
Gives adjacency list and corresponding set of shuttles.
faultyEnd0, faultyEnd1: Flags for creating shuttles with wrong connections
"""
function giveShuttlesAdjacency(;faultyEnd0 = false,faultyEnd1 = false)::
    Tuple{Dict{String,Array{Int64}}, Dict{Symbol,Shuttle}}

    adj = Dict{String,Array{Int64}}()
    shuttles = Dict{Symbol,Shuttle}()
    s = 0
    for i in 1:rand(10:30)
        end0 = rand(setdiff(1:100, keys(adj)))
        for j in 1:rand(5:20)
            end1 = rand(setdiff(1:100, [end0]))
            string(end0) in keys(adj) ? push!(adj[string(end0)], end1) :
            adj[string(end0)] = [end1]
            if faultyEnd0
                shuttles[Symbol(s)] = Shuttle(Symbol(s), Symbol(-1), Symbol(end1))
                faultyEnd0 = false
            end
            if faultyEnd0
                shuttles[Symbol(s)] = Shuttle(Symbol(s), Symbol(end0), Symbol(-1))
                faultyEnd1 = false
            end
            shuttles[Symbol(s)] = Shuttle(Symbol(s), Symbol(end0), Symbol(end1))
            s += 1
        end
    end
    return adj, shuttles
end

"""
Creates a struct GateZoneDesc with repeated Ids
"""
function giveGateZoneDescRepeatedId()::GateZoneDesc
    return GateZoneDesc(
        [
            ZoneInfoDesc( "1", "", "4", 2),
            ZoneInfoDesc( "2", "4", "5", 2),
            ZoneInfoDesc( "1", "7", "8", 2)
        ]
    )
end

"""
Creates a struct LoadZoneDesc with repeated Ids
"""
function giveLoadZoneDescRepeatedId()::LoadZoneDesc
    return LoadZoneDesc(
        [
            LoadZoneInfoDesc( "1", "", "4"),
            LoadZoneInfoDesc( "2", "4", "5"),
            LoadZoneInfoDesc( "1", "7", "8")
        ]
    )
end

"""
Creates a struct GateZoneDesc with inexistent connection
"""
function giveGateZoneDescNoConnection()::GateZoneDesc
    return GateZoneDesc(
        [
            ZoneInfoDesc( "1", "", "4", 2),
            ZoneInfoDesc( "2", "9999999999", "5", 2),
            ZoneInfoDesc( "3", "7", "8", 2)
        ]
    )
end

"""
Creates a struct TrapDesc with wrong connected shuttle
"""
function giveGateZoneDescWrongConnectedShuttle()::GateZoneDesc
    return GateZoneDesc(
        [
            ZoneInfoDesc( "1", "", "4", 2),
            ZoneInfoDesc( "2", "7", "5", 2),
            ZoneInfoDesc( "3", "7", "8", 2)
        ]
    )
end

"""
Creates a struct QCCDevControl based in the file giveQccDes()
"""
function giveQccCtrl(;alternateDesc = false)::QCCDevControl
    qccd::QCCDevDescription = alternateDesc ?  giveQccDes2() : giveQccDes()
    gateZones = Dict{Symbol,GateZone}()
    endId = id -> id == "" ? nothing : Symbol(id)
    map(tr -> gateZones[Symbol(tr.id)] = GateZone(Symbol(tr.id), tr.capacity,
                                endId(tr.end0), endId(tr.end1)), qccd.gateZone.gateZones)
    auxZones = Dict{Symbol,AuxZone}()
    map(sh -> auxZones[Symbol(sh.id)] = AuxZone(Symbol(sh.id), sh.capacity,
                                                endId(sh.end0), endId(sh.end1)),
              qccd.auxZone.auxZones)

    loadZones = Dict{Symbol, LoadingZone}()

    endId = id -> id == "" ? nothing : Symbol(id)
    map(aux -> haskey(loadZones, Symbol(aux.id)) ? throw(err(aux.id)) :
               loadZones[Symbol(aux.id)] = LoadingZone(Symbol(aux.id),
               endId(aux.end0), endId(aux.end1)),
               qccd.loadZone.loadZones)

    aux = (zone,id) -> !isnothing(zone) ? filter(x -> x.end0 == id || x.end1 == id, zone) : []
    junctions = Dict{Symbol,Junction}()
    for j ∈ qccd.junction.junctions
        connectedGateZones = aux(qccd.gateZone.gateZones,j.id)
        connectedAuxZones = aux(qccd.auxZone.auxZones,j.id)
        connectedLoadZones = aux(qccd.loadZone.loadZones,j.id)
        ends = Symbol[]
        map(x -> push!(ends,Symbol(x.id)), connectedGateZones)
        map(x -> push!(ends,Symbol(x.id)), connectedAuxZones)
        map(x -> push!(ends,Symbol(x.id)), connectedLoadZones)
        tmpId = j.id == "" ? nothing : Symbol(j.id)
        junctions[tmpId] = Junction(tmpId, Symbol(j.type), ends)
    end

    return QCCDevControl(qccd,:No,false,gateZones,junctions,auxZones,loadZones)
end

max_stars_repo_name: TartuQC/trapped-ion-qccd-simulator

max_issues_repo_name: TartuQC/trapped-ion-qccd-simulator

max_forks_repo_name: TartuQC/trapped-ion-qccd-simulator


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -19.58980225336104 <= q <= 16.332536463331582)

                     add_NL_constraint(m, :(log(1 + exp(0.1097699546178803 + 0.5720816873107069*log(1 + exp(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))) + 0.11808627755560508*log(1 + exp(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))) + -0.12730937872071246*log(1 + exp(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))) + 0.6166642606025916*log(1 + exp(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))) + 0.45996993866072433*log(1 + exp(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))))) + log(1 + exp(-0.8288330741680832 + -0.5580861219922526*log(1 + exp(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))) + -0.8786358511296881*log(1 + exp(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))) + -0.6915444540383735*log(1 + exp(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))) + 0.5554450106282389*log(1 + exp(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))) + -0.20937059223327026*log(1 + exp(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))))) + log(1 + exp(0.13615042538554878 + 0.08502572267265318*log(1 + exp(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))) + -0.9980047609363689*log(1 + exp(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))) + -0.5678523898140919*log(1 + exp(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))) + 0.29618367003382495*log(1 + exp(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))) + -0.7256587322261439*log(1 + exp(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))))) + log(1 + exp(-0.16139551949919406 + -0.5886735107938796*log(1 + exp(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))) + -0.02630736114123211*log(1 + exp(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))) + -0.03441572289730743*log(1 + exp(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))) + 0.3935924617025486*log(1 + exp(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))) + 0.3425056157482871*log(1 + exp(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))))) + log(1 + exp(-0.5329713534130147 + -0.3183971019485696*log(1 + exp(-0.26202143930164956 + 0.3006071626990803*$(x[1]) + 0.5737922637127819*$(x[2]) + -0.4709190185724106*$(x[3]))) + 0.5847814320797045*log(1 + exp(0.1103867136686838 + -0.1701144375741399*$(x[1]) + -0.3114105380777681*$(x[2]) + -0.8433602212441*$(x[3]))) + -0.7330877632484647*log(1 + exp(0.7136961022484924 + 0.7290114146079865*$(x[1]) + -0.34482664008964514*$(x[2]) + 0.3719033643881229*$(x[3]))) + -0.37731199331975374*log(1 + exp(0.7095130860135548 + 0.7761479980944426*$(x[1]) + 0.5578694028199163*$(x[2]) + 0.7897154219898308*$(x[3]))) + -0.08817423127621415*log(1 + exp(-0.2424710608654208 + -0.6645249084724911*$(x[1]) + -0.49648246781073047*$(x[2]) + 0.9115115047829443*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
import LinearAlgebra.BLAS: @blasfunc, BlasInt, BlasReal

export spr!

@static if VERSION < v"1.7-"
    import LinearAlgebra.BLAS: libblas
    thislib = libblas
else
    import LinearAlgebra.BLAS: libblastrampoline
    thislib = libblastrampoline
end

for (fname, elty) in ((:dspr_, :Float64),
                      (:sspr_, :Float32))
    @eval begin
        function spr!(uplo::AbstractChar,
                      n::Integer,
                      α::$elty,
                      x::Union{Ptr{$elty}, AbstractArray{$elty}},
                      incx::Integer,
                      AP::Union{Ptr{$elty}, AbstractArray{$elty}})

            ccall((@blasfunc($fname), thislib), Cvoid,
                  (Ref{UInt8},     # uplo,
                   Ref{BlasInt},   # n,
                   Ref{$elty},     # α,
                   Ptr{$elty},     # x,
                   Ref{BlasInt},   # incx,
                   Ptr{$elty},     # AP,
                   Clong),         # length of uplo
                  uplo,
                  n,
                  α,
                  x,
                  incx,
                  AP,
                  1)
            return AP
        end
    end
end

function spr!(uplo::AbstractChar,
              α::Real, x::Union{DenseArray{T}, AbstractVector{T}},
              AP::Union{DenseArray{T}, AbstractVector{T}}) where {T <: BlasReal}
    require_one_based_indexing(AP, x)
    N = length(x)
    if 2*length(AP) < N*(N + 1)
        throw(DimensionMismatch("Packed symmetric matrix A has size smaller than length(x) = $(N)."))
    end
    return spr!(uplo, N, convert(T, α), x, stride(x, 1), AP)
end

"""
    spr!(uplo, α, x, AP)

Update matrix `A` as `α*A*x*x'`, where `A` is a symmetric matrix provided
in packed format `AP` and `x` is a vector.

With `uplo = 'U'`, the array AP must contain the upper triangular part of the
symmetric matrix packed sequentially, column by column, so that `AP[1]`
contains `A[1, 1]`, `AP[2]` and `AP[3]` contain `A[1, 2]` and `A[2, 2]`
respectively, and so on.

With `uplo = 'L'`, the array AP must contain the lower triangular part of the
symmetric matrix packed sequentially, column by column, so that `AP[1]`
contains `A[1, 1]`, `AP[2]` and `AP[3]` contain `A[2, 1]` and `A[3, 1]`
respectively, and so on.

The scalar input `α` must be real.

The array inputs `x` and `AP` must all be of `Float32` or `Float64` type.
Return the updated `A`.
"""
spr!


text:
# Source file for working with different decays

struct DecayParameter
    # Generic type for configuring sub-module
    val::AbstractFloat
    fixed::Bool
    glob::Bool
end

function param_vector(pmain::AbstractVector{<:AbstractFloat},
    fxparam::AbstractVector{<:Bool},
    gparam::AbstractVector{<:Bool})

    # Converts guesses to DecayParameters

    # Ensure sizes match
    if length(pmain) != length(fxparam) ||
        length(fxparam) != length(gparam)
        error("Parameter vectors must have same length ($(length(pmain))")
    end
    npar = length(pmain)
    param = Vector{DecayParameter}(undef,npar)
    # Fill in vector
    for i = 1:npar
        param[i] = DecayParameter(pmain[i],fxparam[i],gparam[i])
    end
    return param
end

struct DecayComponent
    # Stores amplitude-weighted exponential component
    α::DecayParameter
    τ::DecayParameter
end

struct TcspcFitModel
    # Configuration for fitting data
    τ::AbstractVector{DecayComponent}
    θ::AbstractVector{DecayComponent}
    Q::DecayParameter
    S::DecayParameter
    V::DecayParameter
    Z::DecayParameter
    function TcspcFitModel(pmain::AbstractVector{DecayParameter},
        ntau::Integer,nani::Integer)
        # Builds fit model using parameter vector

        # Ensure that dimensions match specs
        if length(pmain) != ((2 * ntau) + (2 * nani) + 4)
            error("Parameter vector must be $(((2 * ntau) + (2 * nani) + 4)) elements long")
        end
        if ntau < 1
            error("Must provide at least one lifetime")
        end
        if nani < 0
            error("Number of anisotropies must be nonnegative")
        end
        α = pmain[1:2:(2 * ntau)]
        τ = pmain[2:2:(2 * ntau)]
        τvec = Vector{DecayComponent}(undef,ntau)
        for i = 1:ntau
            τvec[i] = DecayComponent(α[i],τ[i])
        end
        if nani > 0
            sa = (2 * ntau) + 1
            ea = sa + (2 * nani)
            β = pmain[sa:2:ea]
            θ = pmain[(sa + 1):2:ea]
            θvec = Vector{DecayComponent}(undef,nani)
            for i = 1:nani
                θvec[i] = DecayComponent(β[i],θ[i])
            end
        else
            θvec = Vector{DecayComponent}(undef,0)
        end
        Q = pmain[(end - 3)]
        S = pmain[(end - 2)]
        V = pmain[(end - 1)]
        Z = pmain[end]
        new(τvec,θvec,Q,S,V,Z)
    end
end

function TcspcFitModel(pmain::AbstractVector{<:AbstractFloat},
    fxparam::AbstractVector{Bool},
    gparam::AbstractVector{Bool},
    ntau::Integer,nani::Integer)
    # Wrapper for TcspcFitModel construction
    param = param_vector(pmain,fxparam,gparam)
    TcspcFitModel(param,ntau,nani)
end

function ncomp(mod::TcspcFitModel)
    # Returns number of variable components as tuple
    τₙ = length(mod.τ)
    θₙ = length(mod.θ)
    return (τₙ,θₙ)
end

function npar(mod::TcspcFitModel)
    # Returns number of parameters in model
    N = sum(ncomp(mod))
    # Don't include color shift in buffer
    return (2 * N) + 3
end

function param_buffer(mod::TcspcFitModel)
    # Groups parameters for easy access in fitting

    pmain = zeros(3,npar(mod))
    # Row 1: Values
    # Row 2: Fixed/free (specified as 1)
    # Row 3: Global/non-global (specified as 1)

    # Start with tau
    τᵢ = 1
    for τ ∈ mod.τ
        pmain[:,2 * τᵢ - 1] = [τ.α.val;τ.α.fixed;τ.α.glob]
        pmain[:,2 * τᵢ] = [τ.τ.val;τ.τ.fixed;τ.τ.glob]
        τᵢ += 1
    end
    if length(mod.θ) > 0
        θᵢ = τᵢ
        for θ ∈ mod.θ
            pmain[:,2 * θᵢ - 1] = [θ.α.val;θ.α.fixed;θ.α.glob]
            pmain[:,2 * θᵢ] = [θ.τ.val;θ.τ.fixed;θ.τ.glob]
            θᵢ += 1
        end
    end
    pmain[:,end - 2] = [mod.S.val;mod.S.fixed;mod.S.glob]
    pmain[:,end - 1] = [mod.V.val;mod.V.fixed;mod.V.glob]
    pmain[:,end] = [mod.Z.val;mod.Z.fixed;mod.Z.glob]
    return transpose(pmain)
end

function isnonnegative(model::TcspcFitModel)
    # Checks whether linear coefficients are nonnegative
    tf = true
    nlcomp = length(model.τ)
    for τ ∈ model.τ
        if τ.α.val < 0
            tf = false
            return tf
        end
    end
    nacomp = length(model.θ)
    for θ ∈ model.θ
        if θ.α.val < 0
            tf = false
            return tf
        end
    end
    tf = tf && model.S.val >= 0 && model.V.val >= 0 && model.Z.val >= 0
end

function taumat(model::TcspcFitModel,tcal::AbstractFloat,
    nbins::Integer)
    # Column for each lifetime in model
    t = -t_vector(tcal,nbins)
    nlcomp = length(model.τ)
    A = zeros(nbins,nlcomp)
    Ai = 1
    for τ ∈ model.τ
        A[:,Ai] = exp.(t / τ.τ.val)
        Ai += 1
    end
    return A
end

function efftau(model::TcspcFitModel)
    # Returns matrix of lifetimes with amplitudes

    nlcomp = length(model.τ)
    x = zeros(nlcomp,2)
    xi = 1
    for τ ∈ model.τ
        x[xi,1] = τ.α.val
        x[xi,2] = τ.τ.val
        xi += 1
    end
    return x
end

function effanitau(model::TcspcFitModel)
    # Returns matrix of effective lifetimes with anisotropy components

    nlcomp = length(model.τ)
    nacomp = length(model.θ)
    ncomp = nlcomp * nacomp
    x = zeros(ncomp,2)
    xi = 1
    for θ ∈ model.θ
        for τ ∈ model.τ
            x[xi,1] = τ.α.val * θ.α.val
            x[xi,2] = ((1 / τ.τ.val) + (1 / θ.τ.val))
            xi += 1
        end
    end
    x[:,2] = x[:,2] .^ (-1)
    return x
end

function animat(model::TcspcFitModel,tcal::AbstractFloat,
    nbins::Integer)

    t = -t_vector(tcal,nbins)
    τeff = effanitau(model)
    ncomp = size(τeff,1)
    A = zeros(nbins,ncomp)
    for i = 1:ncomp
        A[:,i] = exp.(t / τeff[i,2])
    end
    return A
end

function combmat(model::TcspcFitModel,tcal::AbstractFloat,
    nbins::Integer)
    # Returns the combined tau and anisotropy tau matrix
    A1 = taumat(model,tcal,nbins)
    if length(model.θ) > 0
        A2 = animat(model,tcal,nbins)
        return [A1, A2]
    else
        return A1
    end
end

function combmat(model::TcspcFitModel,tax::TcspcTimeAxis)
    # Does the exact same thing, except with a time axis wrapper
    combmat(model,tax.tcal,tax.nbins)
end

function convmat(model::TcspcFitModel,irf::TcspcTransient)
    # Returns the convolved anisotropy matrix from model
    A = combmat(model,irf.t)
    # Convolve each column with the IRF and replace it
    for i = 1:size(A,2)
        temp = conv(irf.data,A[:,i])
        # Try interpolating to see if things change
        A[:,i] = linitpconv(temp,irf.t)
        # A[:,i] = temp[1:irf.t.nbins]
    end
    return A
end

function linmat(model::TcspcFitModel,irf::TcspcTransient,
    bck::TcspcTransient)
    # Returns full unweighted linear matrix
    if irf.t != bck.t
        error("irf and bck must be same!")
    end
    A = convmat(model,irf)
    return [A, irf.data, bck.data, ones(irf.t.nbins)]
end

function linitpconv(temp::AbstractVector{<:AbstractFloat},
    t::TcspcTimeAxis)
    tnew = range(t.tcal,t.tcal * t.nbins,length=length(temp))
    itp = LinearInterpolation(tnew,temp)
    itp(t_vector(t))
end

function wconvmat(model::TcspcFitModel,irf::TcspcTransient)
    # Returns the weighted convolved anisotropy matrix from model

    A = convmat(model,irf)
    # First few: pure lifetime
    Ai = 1
    for τ ∈ model.τ
        A[:,Ai] *= τ.α.val
        Ai += 1
    end
    # Next: combined anisotropy and lifetime
    for θ ∈ model.θ
        for τ ∈ model.τ
            A[:,Ai] *= (θ.α.val * τ.α.val)
            Ai += 1
        end
    end
    return A
end

function wconv_vec(model::TcspcFitModel,irf::TcspcTransient)
    # Returns the weighted vector convolution

    A = wconvmat(model,irf)
    nlcomp = length(model.τ)
    # Start by summing up the pure lifetimes
    I = (1 / 3) * sum(A[:,1:nlcomp],dims=2)
    # Now, do the same except for anisotropies
    I += (η(irf.t) / 3) * sum(A[:,(nlcomp+1):end],dims=2)
end

function model_vec(model::TcspcFitModel,irf::TcspcTransient,
    bck::TcspcTransient)
    # Returns total model vector
    if irf.t != bck.t
        error("irf and bck must have same time axis!")
    end
    I = wconv_vec(model,irf)
    A1 = model.S.val * irf.data
    A2 = model.V.val * bck.data
    A3 = model.Z.val
    I .+ A1 .+ A2 .+ A3
end


text:
## Test ARX model
function arx_mdl(y, u)
    return (0.5*y[1] - 0.5*u[1])
end

@testset "Test on ARX model" begin
    yterms0 = [1]
    uterms0 = [1]
    u = NarmaxLasso.generate_random_input(100, 5)
    y = NarmaxLasso.simulate(u, arx_mdl, yterms0, uterms0;
                             σv=0)
    mdl = NarmaxLasso.generate_all(NarmaxLasso.NarmaxRegressors,
                                   NarmaxLasso.Monomial,
                                   2, 2, 0, 1)
    result = NarmaxLasso.narmax_lasso(y, u, mdl; nλ=5, λ_min_ratio=1e-6)

    @test result.β ≈ [0.0         -0.0398   0.0        0.1998   0.24836;
                      -1.552e-16  -0.4788  -0.49886  -0.4998  -0.499996;
                      0.0          0.0       0.0       -0.1996  -0.2484;
                      0.0          0.43622   0.49973   0.89964   0.99671] rtol=1e-4

    if NarmaxLasso.glmnet_loaded
        result_glmnet = NarmaxLasso.narmax_lasso(y, u, mdl;
            intercept=false, standardize=false, use_glmnet=true, nlambda=5)
        @test result_glmnet.β ≈ [0.0 -0.13385 -0.010741 0.0 0.0349759;
                                 0.0 -0.432755 -0.493341 -0.498908 -0.499592;
                                 0.0 0.0 0.0 0.0 -0.0340471;
                                 0.0 0.289839 0.481891 0.499722 0.569338] rtol=1e-4
    else
        @test_throws(ArgumentError,
         result_glmnet = NarmaxLasso.narmax_lasso(y, u, mdl; use_glmnet=true))
    end
end

## Test ARMAX model
function armax_mdl(y, u, v)
    return (0.5*y[1] - 0.5*u[1] + 0.5*v[1])
end

@testset "Test on ARMAX model" begin
    yterms0 = [1]
    uterms0 = [1]
    vterms0 = [1]
    u = NarmaxLasso.generate_random_input(1000, 5)
    y = NarmaxLasso.simulate(u, armax_mdl, yterms0, uterms0, vterms0;
                             σv=1.0, seed=1)
    mdl = NarmaxLasso.generate_all(NarmaxLasso.NarmaxRegressors,
                                   NarmaxLasso.Monomial,
                                   1, 1, 1, 1)
    result = NarmaxLasso.narmax_lasso(y, u, mdl; nλ=5)
    @test result.β ≈ [0.0  -0.204043  -0.457963  -0.4847    -0.487381;
                      0.0   0.555936   0.478094   0.469134   0.468238;
                      0.0   0.240684   0.514044   0.546466   0.549734] rtol=1e-3
end


text:
using SciPyDiffEq

# Definition of parameters
const a = 1.0
const b = 0.1
const c = 1.5
const d = 0.75

function diffeq(u,p,t)
    du = zeros(2)
    # u1: rabbits
    # u2: foxes
    du[1] = a*u[1] - b*u[1]*u[2]
    du[2] = -c*u[2] + d*b*u[1]*u[2]

    return du

end

# Simulation time
const tspan = (0.0,18.0)

# Initial condition
u0 = zeros(2)
u0[1] = 10
u0[2] = 5

# Solve ODE
prob = ODEProblem(diffeq,u0,tspan)
sol = solve(prob,SciPyDiffEq.BDF(),saveat=(tspan[end]-tspan[1])/999)

#= Visualization of Results
using PyPlot
figure(figsize=(20,8))
plot(sol.t,sol.u)
show()
=#


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE

using AbaqusReader: create_surface_element, create_surface_elements

@testset "create surface element from voluminal element surface" begin
    element = create_surface_element(:Tet4, :S1, [8, 9, 10, 2])
    @test element == (:Tri3, [8, 10, 9])
end

@testset "create surface elements from voluminal element surface" begin
    mesh = Dict(
        "elements" => Dict(16 => [8, 9, 10, 2]),
        "element_types" => Dict(16 => :Tet4),
        "surface_sets" => Dict("LOAD" => [(16, :S1)]))
    elements = create_surface_elements(mesh, "LOAD")
    @test elements[1] == (:Tri3, [8,10,9])
end

@testset "throw error if unknown element or side defined" begin
    @test_throws(Exception, create_surface_element(:Tet5, :S1, [8, 9, 10, 2]))
    @test_throws(Exception, create_surface_element(:Tet4, :S5, [8, 9, 10, 2]))
end


text:
using Plots
using Primes
using Memoize # https://github.com/hymd-research/Memoize.jl.git

@memoize function g(N::T) where {T<:Integer}
    if N == one(T) # meaning N==1
        return one(T)
    else
        fN = sum(e for (_, e) ∈ factor(N).pe)
        ifelse(iseven(fN), 1, -1) + g(N - 1)
    end
end

@time gs = [g(N) for N = 1:1000000]
@assert sum(gs .≥ 0) == 10
scatter(gs)

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
export buffachristiansen

"""
Move the s-th element right after the d-th
"""
function move_after!(p,n,s,d)

    n[d] == s && return

    t1 = n[d]
    t2 = n[s]
    t3 = p[s]

    n[d] = s
    n[s] = t1
    p[s] == 0 || (n[p[s]] = t2)

    p[s] = d
    t1 == 0 || (p[t1] = s)
    t2 == 0 || (p[t2] = t3)
end

function move_before!(p, n, s, d)
    p[d] == s && return

    t1 = p[d]
    t2 = p[s]
    t3 = n[s]

    p[d] = s
    p[s] = t1
    n[s] == 0 || (p[n[s]] = t2)

    n[s] = d
    t1 == 0 || (n[t1] = s)
    t2 == 0 || (n[t2] = t3)
end

function sortneighbors(a, pred)

    n = collect(2:length(a)+1); n[end] = 0
    p = collect(0:length(a)-1); p[1] = 0

    last = 1
    while last != 0
        cand = n[last]
        while cand != 0
            pred(a[last], a[cand]) && break
            cand = n[cand]
        end
        cand == 0 && break
        move_after!(p, n, cand, last)
        last = cand
    end

    first = 1
    while last != 0
        cand = n[last]
        while cand != 0
            pred(a[cand], a[first]) && break
            cand = n[cand]
        end
        cand == 0 && break
        move_before!(p, n, cand, first)
        first = cand
    end

    b = similar(a)
    i, j = last, length(n)
    while true
        b[j] = a[i]
        i = p[i]
        i == 0 && break
        j -= 1
    end

    return b

end

isclosed(a, pred) = length(a)>2 && pred(a[end], a[1])

"""
    buffachristiansen(Γ, γ)

Construct the set of Buffa-Christiansen functions subject to mesh Γ and only
enforcing zero normal components on ∂Γ ∖ γ.
"""
function buffachristiansen(Γ, γ=mesh(coordtype(Γ),1,3))

    @assert CompScienceMeshes.isoriented(Γ)

    T = coordtype(Γ)

    edges = skeleton(Γ, 1)
    fine = barycentric_refinement(Γ)

    in_interior = interior_tpredicate(Γ)
    on_junction = overlap_gpredicate(γ)

    # first pass to determine the number of functions
    numfuncs = 0
    for edge in cells(edges)
        ch = chart(edges, edge)
        !in_interior(edge) && !on_junction(ch) && continue
        numfuncs += 1
    end

    vtof, vton = vertextocellmap(fine)
    jct_pred = overlap_gpredicate(γ)
    bcs, k = Vector{Vector{Shape{T}}}(undef,numfuncs), 1
    for (i,edge) in enumerate(cells(edges))

        ch = chart(edges, edge)
        !in_interior(edge) && !on_junction(ch) && continue

        # index of edge center in fine's vertexbuffer
        p = numvertices(edges) + i

        v = edge[1]
        n = vton[v]
        F = vec(vtof[v,1:n])
        supp = cells(fine)[F]
        bc1 = buildhalfbc(fine, F, v, p, jct_pred)

        v = edge[2]
        n = vton[v]
        F = vec(vtof[v,1:n])
        supp = cells(fine)[F]
        bc2 = buildhalfbc(fine, F, v, p, jct_pred)
        bc2 = Shape{T}[Shape(s.cellid, s.refid, -s.coeff) for s in bc2]

        bcs[k] = [bc1; bc2]
        k += 1
    end

    return RTBasis(fine, bcs)
end

"""
    buildhalfbc(fine, supp::Array{SVector{3,Int},1}, v, p)
"""
function buildhalfbc(fine, S, v, p, onjunction)

    T = coordtype(fine)

    @assert length(S) >= 2
    @assert mod(length(S), 2) == 0
    @assert all(0 .< S .<= numcells(fine))
    @assert all([v in cells(fine)[s] for s in S])

    bf = Shape{T}[]

    n = length(S)
    N = n ÷ 2
    modn(i) = mod1(i,n)

    # sort the fine faces making up the domain
    share_edge(i,j) = length(intersect(cells(fine)[i],cells(fine)[j])) == 2
    S = sortneighbors(S, share_edge)
    @assert all(0 .< S .<= numcells(fine))

    port_faces = findall(i -> p in cells(fine)[i], S)
    port_edges = [something(findfirst(isequal(v), cells(fine)[S[i]]),0) for i in port_faces]

    num_ports = length(port_faces)
    port_fluxes = ones(T,num_ports) / num_ports
    @assert length(port_faces) < 3

    c_on_boundary = !share_edge(S[end], S[1]) || length(S) == 2
    e_on_boundary = length(port_faces) == 1
    @assert c_on_boundary || !e_on_boundary

    # if there are two ports (c must be interior), we want to iterate
    # the support starting with one and ending with the other. This
    # might require reordering the ports we just discovered.
    if !e_on_boundary
        @assert length(port_faces) == 2 "$port_faces, $n"
        if port_faces[2] == modn(port_faces[1]+1)
            reverse!(port_faces)
            reverse!(port_edges)
        end
        @assert port_faces[1] == modn(port_faces[2]+1)
    end

    # Detect the boundary edges
    bnd_faces = Int[]
    bnd_edges = Int[]
    if c_on_boundary
        bnd_faces = [1,n]
        p1 = x -> x in cells(fine)[S[2]] && x != v
        p2 = x -> x in cells(fine)[S[n-1]] && x != v
        e1 = findfirst(p1, cells(fine)[S[1]])
        e2 = findfirst(p2, cells(fine)[S[n]])
        bnd_edges = [e1, e2]
    end

    # Detect the junction edges
    num_junctions = 0
    jct_idxes = Int[]
    if c_on_boundary
        for i in 1:length(bnd_faces)
            f, e = bnd_faces[i], bnd_edges[i]
            a = vertices(fine)[cells(fine)[S[f]][mod1(e+1,3)]]
            b = vertices(fine)[cells(fine)[S[f]][mod1(e+2,3)]]
            seg = simplex(a,b)
            if onjunction(seg)
                push!(jct_idxes, i)
                #push!(jct_edges, e)
            end
        end
        num_junctions = length(jct_idxes)
    end

    # This charge needs to be compensated by interior divergence
    total_charge = (!c_on_boundary || num_junctions == 2) ? 1 : 0
    charges = fill(total_charge/n, n)

    # add the port contributions
    for (f, e, w) in zip(port_faces, port_edges, port_fluxes)
        add!(bf, S[f], e, w)
        charges[f] -= w
    end

    bnd_fluxes = T[]
    if c_on_boundary
        if num_junctions == 0
            bnd_fluxes = T[-(n-port_faces[2])/n, -port_faces[2]/n]
        elseif num_junctions == 1
            bnd_fluxes = -ones(T,2)
            bnd_fluxes[jct_idxes[1]] = 0
        else
            bnd_fluxes = zeros(T,2)
        end
    end

    @assert length(bnd_faces) == length(bnd_edges) == length(bnd_fluxes)

    for (f, e, w) in zip(bnd_faces, bnd_edges, bnd_fluxes)
        add!(bf, S[f], e, w)
        charges[f] -= w
    end

    # The remaining fluxes can be computed based on port and boundary
    # influx, and the knowledge of the charge density by recursion.
    # Iterate beginning with port 0:
    start_face = c_on_boundary ? 1 : port_faces[1]
    for i in 1 : n-1
        # cyclic port_face[1]-based indexing
        j0 = modn(i + start_face - 1)
        j1 = modn(j0+1)
        # get the indices of the faces in fine.faces
        f0, f1 = S[j0], S[j1]
        # what are the local indices of the common edge?
        e0, e1 = getcommonedge(cells(fine)[f0], cells(fine)[f1])
        # add the two half Raviart-Thomas shape functions
        add!(bf, S[j0], abs(e0), +charges[j0])
        add!(bf, S[j1], abs(e1), -charges[j0])
        # update the charge bookkeeping
        charges[j1] += charges[j0]
    end

    return bf
end


text:
using Test
using Quantica

@testset "Quantica.jl" begin
    include("test_lattice.jl")
    include("test_model.jl")
    include("test_hamiltonian.jl")
    include("test_mesh.jl")
    include("test_bandstructure.jl")
end


text:
using ConstantQ
using Base.Test

import ConstantQ: rawdata, _speckernel, _tempkernel, Frequency
import DSP: hamming, hanning

immutable DummyFrequency <: Frequency
end

let
    freq = DummyFrequency()
    @test_throws Exception nfreqs(freq)
    @test_throws Exception freqs(freq)
end

let
    GeometricFrequency(60, 61)

    # keyword arguments
    GeometricFrequency(min=60, max=5000)
    GeometricFrequency(min=60, max=5000, bins=24)
end

let
    GeometricFrequency(60, 61)
    GeometricFrequency(60, 60)
    @test_throws ArgumentError GeometricFrequency(60, 59)

    GeometricFrequency(60, 5000, 24)
    @test_throws ArgumentError GeometricFrequency(60, 5000, -1)
end

let
    for bins in [24, 48]
        freq = GeometricFrequency(min=55, max=2000, bins=bins)
        @test nbins_per_octave(freq) == bins
    end
end

let
    freq = GeometricFrequency(min=60, max=5000, bins=24)
    @test isapprox(q(freq), 34.12708770892056)
end

let
    @test nfreqs(GeometricFrequency(60, 120, 20)) == 20
    @test nfreqs(GeometricFrequency(60, 240, 20)) == 40
    @test nfreqs(GeometricFrequency(60, 120, 30)) == 30
    @test nfreqs(GeometricFrequency(60, 240, 30)) == 60
end

# constant diff in log-frequency domain
let
    f = freqs(GeometricFrequency(min=60, max=5000))
    d = diff(log.(f))
    for e in d
        @test isapprox(d[1], e)
    end
end

function invalid_win()
    println("invalid")
end

let
    fs = 16000
    freq = GeometricFrequency(min=60, max=5000)
    KernelProperty(fs, freq, hamming)
    @test_throws Exception KernelProperty(fs, freq, invalid_win)
end

# SpectralKernelMatrix
let
    fs = 16000
    freq = GeometricFrequency(min=60, max=5000)
    kp = KernelProperty(fs, freq, hamming)

    S = spzeros(Complex{Float64}, 5, 5)
    K = SpectralKernelMatrix(S, kp)
    @test property(K) == kp
    @test rawdata(K) == S
    @test issparse(rawdata(K))
    @test issparse(K)

    S = zeros(Complex{Float64}, 5, 5)
    K = SpectralKernelMatrix(S, kp)
    @test rawdata(K) == S
    @test !issparse(rawdata(K))
    @test !issparse(K)
end

# TemporalKernelMatrix
let
    fs = 16000
    freq = GeometricFrequency(min=60, max=5000)
    kp = KernelProperty(fs, freq, hamming)

    S = spzeros(Complex{Float64}, 5, 5)
    K = TemporalKernelMatrix(S, kp)
    @test property(K) == kp
    @test rawdata(K) == S
    @test issparse(rawdata(K))
    @test issparse(K)

    S = zeros(Complex{Float64}, 5, 5)
    K = TemporalKernelMatrix(S, kp)
    @test rawdata(K) == S
    @test !issparse(rawdata(K))
    @test !issparse(K)
end

# _speckernel
let
    fs = 16000
    freq = GeometricFrequency(174.5, fs/2)
    K = _speckernel(Float64, fs, freq, hamming, 0.005)
    @test issparse(K)
    @test eltype(K) == Complex{Float64}

    K = _speckernel(Float32, fs, freq, hamming, 0.005)
    @test issparse(K)
    @test eltype(K) == Complex{Float32}
end

# speckernel
let
    fs = 16000
    freq = GeometricFrequency(174.5, fs/2)
    K = speckernel(Float64, fs, freq, hamming, 0.005)
    @test isa(K, SpectralKernelMatrix)
    @test issparse(K)

    K = speckernel(Float64, fs)
    @test isa(K, SpectralKernelMatrix)
    @test issparse(K)

    rawK = ConstantQ.rawdata(K)
    @test size(K) == size(rawK)
    @test length(K) == length(rawK)

    # getindex
    for i=1:10
       @test rawK[i] == K[i]
    end

    @test isa(full(K), DenseMatrix)
end

# _tempkernel
let
    fs = 16000
    freq = GeometricFrequency(174.5, fs/2)
    K = _tempkernel(Float64, fs, freq, hamming)
    @test !issparse(K)
    @test eltype(K) == Complex{Float64}

    K = _tempkernel(Float32, fs, freq, hamming)
    @test !issparse(K)
    @test eltype(K) == Complex{Float32}
end

# tempkernel
let
    fs = 16000
    freq = GeometricFrequency(min=220, max=440)

    # compute speckernel
    K = speckernel(Float64, fs, freq, hamming, 0.0)
    K = rawdata(K)
    fftlen = size(K, 1)

    # back to tempkernel
    k = conj(ifft(full(conj(K .* fftlen)), 1))

    # check correctness
    expected = rawdata(tempkernel(Float64, fs, freq, hamming))
    @test isapprox(k, expected)
end

# cqt user-friendly interface
let
    srand(98765)
    x = rand(Float64, 60700)
    fs = 16000

    freq = GeometricFrequency(min=60, max=5000)
    K = speckernel(Float64, fs, freq)
    X, timeaxis, freqaxis = cqt(x, fs, K)
    @test isa(X, Matrix{Complex{Float64}})

    @test first(timeaxis) == 0.0
    @test last(timeaxis) <= length(x)/fs
    @test first(freqaxis) == 60
    @test last(freqaxis) <= 5000
end

# check correctness
let
    srand(98765)
    x = rand(Float64, 60700)
    fs = 16000
    hopsize = convert(Int, round(fs * 0.01))
    freq = GeometricFrequency(min=220, max=fs/2)

    # CQT in frequency-domain
    # make sure threshold to be zero
    K = speckernel(Float64, fs, freq, hamming, 0.0)
    @test isa(K, SpectralKernelMatrix)
    X1, t1, f1 = cqt(x, fs, K, hopsize=hopsize)

    # CQT in time-domain
    k = tempkernel(Float64, fs, freq, hamming)
    @test isa(k, TemporalKernelMatrix)
    X2, t2, f2 = cqt(x, fs, k, hopsize=hopsize)

    @test isapprox(norm(X1), norm(X2), atol=1.0e-7)
    @test isapprox(X1, X2, atol=1.0e-7)
    @test isapprox(t1, t2)
    @test isapprox(f1, f2)
end