text:
import Distributions.entropy

immutable TSBPMM
    α::Float64
    qv::Vector{Beta}
    π::Vector{Float64}
    z::Matrix{Float64}

    cluster_update::Function
    cluster_loglikelihood::Function
    object_loglikelihood::Function
    cluster_entropy::Function

    function TSBPMM(N::Int64, T::Int64, α::Float64,
            logpi::Vector{Float64}, z::Matrix{Float64},
            cluster_update::Function,
            cluster_loglikelihood::Function,
            object_loglikelihood::Function,
            cluster_entropy::Function)
        @assert length(logpi) == T "logpi must have size T"
        @assert size(z) == (N, T) "z has incostistent size"

        return new(α, [Beta(1., α) for i=1:T-1],
                logpi,
                z,
                cluster_update,
                cluster_loglikelihood,
                object_loglikelihood,
                cluster_entropy)
    end
end

function TSBPMM(N::Int64, T::Int64, α::Float64,
            cluster_update::Function,
            cluster_loglikelihood::Function,
            object_loglikelihood::Function,
            cluster_entropy::Function; random_init=false)
    logpi = zeros(T)
    z = zeros(N, T)

    if random_init
        rand!(logpi)
        ps = sum(logpi)
        @devec logpi ./= ps
        log!(logpi)

        rand!(z)
        for i=1:N
            zs = sum(z[i, :])
            z[i, :] ./= zs
        end
    else
        logpi[:] = -log(T)
        z[:, :] = 1. / T
    end

    return TSBPMM(N, T, α,
            logpi,
            z,
            cluster_update,
            cluster_loglikelihood,
            object_loglikelihood,
            cluster_entropy)
end

N(mix::TSBPMM) = size(mix.z, 1)
truncation_level(mix::TSBPMM) = size(mix.z, 2)

function infer(mix::TSBPMM, niter::Int64, ltol::Float64;
        iter_callback::Function = (oksa...) -> begin end,
        disable_dp = false)
    prev_lb = variational_lower_bound(mix)
    for iter=1:niter
        variational_update(mix; disable_dp=false)

        lb = variational_lower_bound(mix)

        iter_callback(mix, iter, lb)

        if abs(lb - prev_lb) < ltol
            return iter
        end

        prev_lb = lb
    end

    return niter
end

function variational_update(mix::TSBPMM; disable_dp=false)
    if !disable_dp
        ts = 0.
        for k=truncation_level(mix):-1:1
            zk = view(mix.z, :, k)
            mix.cluster_update(k, zk)
            zs = sum(zk)
            if k < truncation_level(mix)
                mix.qv[k] = Beta(1. + zs, mix.α + ts)
            end
            ts += zs
        end

        logpi!(mix.π, mix)
    else
        mix.π[:] = 0.
    end

    z = zeros(truncation_level(mix))
    for i=1:N(mix)
        for k=1:truncation_level(mix)
            z[k] = mix.π[k] + mix.object_loglikelihood(k, i)
        end

        @devec z[:] -= maximum(z)
        exp!(z)
        @devec z[:] ./= sum(z)

        assert(abs(sum(z) - 1.) < 1e-7)

        mix.z[i, :] = z
    end
end

function variational_lower_bound(mix::TSBPMM)
    return loglikelihood(mix) + entropy(mix)
end

meanlog(beta::Beta) = digamma(beta.α) - digamma(beta.α + beta.β)
meanlogmirror(beta::Beta) = digamma(beta.β) - digamma(beta.α + beta.β)
meanmirror(beta::Beta) = beta.β/ (beta.α + beta.β)
logmeanmirror(beta::Beta) = log(beta.β) - log(beta.α + beta.β)

function logpi!(π::Vector{Float64}, mix::TSBPMM)
    r = 0.
    for k=1:truncation_level(mix)-1
        π[k] = meanlog(mix.qv[k]) + r
        r += meanlogmirror(mix.qv[k])
    end
    π[truncation_level(mix)] = r
end

function loglikelihood(mix::TSBPMM)
    ll = 0.

    ts = 0.
    for k=truncation_level(mix):-1:1
        zk = view(mix.z, :, k)
#        zk = mix.z[:, k]
        ll += mix.cluster_loglikelihood(k, zk)
        assert(!isnan(ll))

        zs = sum(zk)
        if k <= truncation_level(mix) - 1
            qv = mix.qv[k]
            ll += zs * meanlog(qv) + (mix.α+ts-1) * meanlogmirror(qv) - lbeta(1., mix.α)
            assert(!isnan(ll))
        end

        ts += zs
    end

    return ll
end

function entropy(mix::TSBPMM)
    ee = 0.
    ee += entropy(mix.z)

    for k=1:truncation_level(mix)
        if k < truncation_level(mix)
            ee += entropy(mix.qv[k])
        end
        ee += mix.cluster_entropy(k)
    end

    return ee
end

function map_assignments(mix::TSBPMM)
    z = zeros(Int64, N(mix))
    for i=1:N(mix)
        z[i] = indmax(sub(mix.z, i, :))
    end
    return z
end

function pi!(π::Vector{Float64}, mix::TSBPMM)
    r = 0.
    for k=1:truncation_level(mix)-1
        qv = mix.qv[k]
        π[k] = exp(log(mean(qv)) + r)
        r += logmeanmirror(qv)
    end
    π[truncation_level(mix)] = exp(r)
    assert(abs(sum(π) - 1.) < 1e-7)
end

export TSBPMM, infer, variational_lower_bound, map_assignments, pi!, T


text:
macro _change_bounds_fn(model_type, idx_type, args...)
    body = last(args)
    typeof(body) == Expr || throw(DomainError(body, "missing function body"))
    plural = :plural in args
    plural_s = plural ? "s" : ""
    inplace = :inplace in args
    fname = Symbol(:change_bound, plural_s, inplace ? "!" : "")
    idx_var = Symbol(
        :rxn,
        idx_type == :Int ? "_idx" :
        idx_type == :String ? "_id" :
        throw(DomainError(idx_type, "unsupported index type for change_bound macro")),
        plural_s,
    )
    example_idx =
        plural ? (idx_type == :Int ? [123, 234] : ["ReactionA", "ReactionC"]) :
        (idx_type == :Int ? 123 : "\"ReactionB\"") #= unquoting is hard =#
    example_val = plural ? [4.2, 100.1] : 42.3
    missing_default = plural ? :((nothing for _ in $idx_var)) : nothing

    bound_type = Float64
    if plural
        idx_type = AbstractVector{eval(idx_type)}
        bound_type = AbstractVector{bound_type}
    end

    docstring = """
        $fname(
            model::$model_type,
            $idx_var::$idx_type;
            lower = $missing_default,
            upper = $missing_default,
        )

    Change the specified reaction flux bound$(plural_s) in the model
    $(inplace ? "in-place" : "and return the modified model").

    # Example
    ```
    $(inplace ? "new_model = " : "")$fname(model, $example_idx, lower=$(-0.5 .* example_val), upper=$example_val)
    ```
    """

    return :(@doc $docstring $fname(
        model::$model_type,
        $idx_var::$idx_type;
        lower = $missing_default,
        upper = $missing_default,
    ) = $body)
end

ext: jl


text:
"""Solve the initializationproblem: construct matrices of indirect utility and
initial choice depending on exogenous variables."""
initialendogstatevars(p::DDP, meshValFun::Array{Float64}) =
    _initialendogstatevars(
        p.options.initialize.problem, meshValFun,
        p.tStateVectors, getchoicevarszero(p), getnonchoicevarszero(p))

function _initialendogstatevars(initializationproblem, meshValFun::Array{Float64},
    tStateVectors,
    tChoiceVectors::Tuple{Vector{Float64}}, tExogStateVectors)

    exogtNodes = length.(tExogStateVectors)
    mV0 = zeros(exogtNodes)
    mChoice0 = zeros(exogtNodes)

    vActions::Vector{Real} = tChoiceVectors[1]

    for iexog in CartesianIndices(exogtNodes)

        interimvalue = -Inf
        interimChoice = 0

        for iChoice = 1:length(vActions)

                value = initializationproblem(meshValFun[iChoice,iexog.I...], vActions[iChoice])

                if value > interimvalue
                    interimvalue = value
                    interimChoice = vActions[iChoice]
                end

        end

        mV0[iexog.I...] = interimvalue
        mChoice0[iexog.I...] = interimChoice

    end

    return mV0, (mChoice0,)

end

function _initialendogstatevars(initializationproblem, meshValFun::Array{Float64},
    tStateVectors, tChoiceVectors::NTuple{2,Vector{Float64}}, tExogStateVectors)

    exogtNodes = length.(tExogStateVectors)
    mV0 = zeros(exogtNodes)
    mChoiceOne0 = zeros(exogtNodes)
    mChoiceTwo0 = zeros(exogtNodes)

    vActionsOne, vActionsTwo = tChoiceVectors

    for iexog in CartesianIndices(exogtNodes)

        interimvalue = -Inf
        interimChoiceOne = 0
        interimChoiceTwo = 0

        for iChoiceOne = 1:length(vActionsOne), iChoiceTwo = 1:length(vActionsTwo)

                value = initializationproblem(meshValFun[iChoiceOne, iChoiceTwo, iexog.I...],
                            (vActionsOne[iChoiceOne], vActionsTwo[iChoiceTwo]))

                if value > interimvalue
                    interimvalue = value
                    interimChoiceOne = vActionsOne[iChoiceOne]
                    interimChoiceTwo = vActionsTwo[iChoiceTwo]
                end

        end

        mV0[iexog.I...] = interimvalue
        mChoiceOne0[iexog.I...] = interimChoiceOne
        mChoiceTwo0[iexog.I...] = interimChoiceTwo

    end

    return mV0, (mChoiceOne0, mChoiceTwo0)

end


text:
using LinearAlgebra
@testset "degopt" begin
    x1a = [3; 4.0]
    x1b = [0.1; 0.3]
    x2a = [0.1; 0.1; -0.1]
    x2b = [0.3; 0; -0.1]
    x = [(x1a, x1b), (x2a, x2b)]

    z = [1.0; -0.1; 0.1; 0.01]

    degopt1 = Degopt(deepcopy(x), deepcopy(z))
    # testing graph_degopt
    (g1, _) = graph_degopt(degopt1)

    # Testing grow
    degopt2 = grow!(deepcopy(degopt1))
    (g2, _) = graph_degopt(degopt1)
    @test eval_graph(g1, 0.3) == eval_graph(g2, 0.3)

    # Test get_degopt_crefs
    (x_crefs, y_crefs) = get_degopt_crefs(2)
    set_coeffs!(g1, 5.0, [x_crefs[1][2][2]])

    x[1][2][2] = 5.0
    degopt3 = Degopt(deepcopy(x), deepcopy(z))
    (g3, _) = graph_degopt(degopt3)

    @test eval_graph(g1, -3.0) == eval_graph(g3, -3.0)

    # Test scale!

    degopt4 = deepcopy(degopt1)
    scale!(degopt4, 2)
    (g4, _) = graph_degopt(degopt4)
    (g1, _) = graph_degopt(degopt1)
    @test eval_graph(g4, 0.3) == eval_graph(g1, 2 * 0.3)

    # Test square!
    degopt5 = deepcopy(degopt1)
    square!(degopt5)
    (g5, _) = graph_degopt(degopt5)
    @test eval_graph(g5, 2 + 2im) == eval_graph(g1, 2 + 2im)^2

    # Test initialize as degopt
    degopt6 = Degopt(g1)
    (g6, _) = graph_degopt(degopt6)
    @test eval_graph(g6, 10.0) == eval_graph(g1, 10.0)

    # Normalization

    degopt8 = normalize!(deepcopy(degopt1))
    (g8, _) = graph_degopt(degopt8)
    (g1, _) = graph_degopt(degopt1)
    @test eval_graph(g8, pi) ≈ eval_graph(g1, pi)

    # get_degopt_coeff
    (HA, HB, y) = get_degopt_coeffs(g8)
    (g9, _) = graph_degopt(Degopt(HA, HB, y))
    @test eval_graph(g8, pi) == eval_graph(g9, pi)
end


text:
struct MyModel <: AbstractMCMC.AbstractModel end

struct MySample{A,B}
    a::A
    b::B
end

struct MySampler <: AbstractMCMC.AbstractSampler end
struct AnotherSampler <: AbstractMCMC.AbstractSampler end

struct MyChain{A,B,S} <: AbstractMCMC.AbstractChains
    as::Vector{A}
    bs::Vector{B}
    stats::S
end

MyChain(a, b) = MyChain(a, b, NamedTuple())

function AbstractMCMC.step(
    rng::AbstractRNG,
    model::MyModel,
    sampler::MySampler,
    state::Union{Nothing,Integer} = nothing;
    sleepy = false,
    loggers = false,
    kwargs...
)
    # sample `a` is missing in the first step
    a = state === nothing ? missing : rand(rng)
    b = randn(rng)

    loggers && push!(LOGGERS, Logging.current_logger())
    sleepy && sleep(0.001)

    _state = state === nothing ? 1 : state + 1

    return MySample(a, b), _state
end

function AbstractMCMC.bundle_samples(
    samples::Vector{<:MySample},
    model::MyModel,
    sampler::MySampler,
    ::Any,
    ::Type{MyChain};
    stats = nothing,
    kwargs...
)
    as = [t.a for t in samples]
    bs = [t.b for t in samples]

    return MyChain(as, bs, stats)
end

function isdone(
    rng::AbstractRNG,
    model::MyModel,
    s::MySampler,
    samples,
    state,
    iteration::Int;
    kwargs...
)
    # Calculate the mean of x.b.
    bmean = mean(x.b for x in samples)
    return abs(bmean) <= 0.001 || iteration >= 10_000 || state >= 10_000
end

# Set a default convergence function.
function AbstractMCMC.sample(model, sampler::MySampler; kwargs...)
    return sample(Random.GLOBAL_RNG, model, sampler, isdone; kwargs...)
end

function AbstractMCMC.chainscat(
    chain::Union{MyChain,Vector{<:MyChain}},
    chains::Union{MyChain,Vector{<:MyChain}}...
)
    return vcat(chain, chains...)
end

# Conversion to NamedTuple
Base.convert(::Type{NamedTuple}, x::MySample) = (a = x.a, b = x.b)


text:
module Covar

export @covar_system, ⟂, val, var, err, numvals, CorrelatedVar, CovariantVar, DerivedVar,
       CovariantSystem

using StaticArrays

"""
Represents a system of covariant variables

    vals:   Values of the covariant variables
    covar:  The covariance matrix
"""
struct CovariantSystem{T<:AbstractFloat,AV<:AbstractVector{T},AM<:AbstractMatrix{T}}
    vals::AV
    covar::AM

    function CovariantSystem(vals::AV, covar::AM) where
        {T<:AbstractFloat,AV<:AbstractVector{T},AM<:AbstractMatrix{T}}
        @assert(size(covar) == length.((vals, vals)),
                "CovariantSystem has to satisfy: vals ∈ ℝ^n, covar ∈ ℝ^(n×n) where n ∈ ℕ")
        return new{T,AV,AM}(vals, covar)
    end
end

numvals(x::CovariantSystem) = length(x.vals)

function Base.convert(::Type{CovariantSystem{_,AV,AM}}, cs) where {_,AV,AM}
    return CovariantSystem(convert(AV, cs.vals), convert(AM, cs.covar))
end

function Base.convert(::Type{CovariantSystem{T,AV,AM}},
                      cs::CovariantSystem{T,AV,AM}) where {T,AV,AM}
    return cs
end

function Base.promote_rule(::Type{CovariantSystem{T1,AV1,AM1}}, ::Type{CovariantSystem{T2,AV2,AM2}}) where
    {T1,AV1,AM1,T2,AV2,AM2}
    return CovariantSystem{promote_type(T1, T2),promote_type(AV1, AV2),promote_type(AM1, AM2)}
end

abstract type CorrelatedVar <: Number end

"""
Represents a variable of a CovariantSystem

    index:  Index of the variable in the CovariantSystem
    system: CovariantSystem of the variable
"""
struct CovariantVar{T<:AbstractFloat,AV<:AbstractVector{T},AM<:AbstractMatrix{T}} <: CorrelatedVar
    index::Int
    system::CovariantSystem{T,AV,AM}

    function CovariantVar(index::Int, system::CovariantSystem{T,AV,AM}) where
        {T<:AbstractFloat,AV<:AbstractVector{T},AM<:AbstractMatrix{T}}
        @assert index ≤ length(system.vals) "Index of CovariantVar not in range"
        return new{T,AV,AM}(index, system)
    end
end

system(x::CovariantVar) = x.system
val(x::CovariantVar) = system(x).vals[x.index]
var(x::CovariantVar) = system(x).covar[x.index, x.index]

⟂(x::CovariantVar, y::CovariantVar) = x.system !== y.system

"""
Creates a CovariantSystem and returns an Array containing the CovariantVars

See also: `CovariantSystem`, `CovariantVar`
"""
macro covar_system(vals, covar)
    return quote
        begin
            cs = CovariantSystem($(esc(vals)), $(esc(covar)))
            CovariantVar.(keys(cs.vals), Ref(cs))
        end
    end
end

function Base.show(io::IO, x::CovariantVar)
    println(io, "x_$(x.index) = $(val(x)) ± $(err(x))")
    print(io, "  Vₘ = $(x.system.covar)")
end

include("gradients.jl")

"""
Represents a variable derived from one or multiple CovariantSystems

    val:     Value of the Variable
    grads:   Gradients with respect to the corresponding CovariantSystems
    systems: The underlying CovariantSystems
"""
struct DerivedVar{T<:AbstractFloat,G<:Gradients} <: CorrelatedVar
    val::T
    grads::G
end

val(x::DerivedVar) = x.val

function var(x::DerivedVar{T}) where T
    var = zero(T)
    for (system,grad) in x.grads
        var += grad' * (system.covar * grad)
    end
    return var
end

err(x::CorrelatedVar) = √var(x)

function Base.show(io::IO, x::DerivedVar)
    print(io, "f({xᵢ}) = $(val(x)) ± $(err(x))")
end

function onehot(::Type{AV}, val, i::Int, len::Int) where
    {T<:AbstractFloat,AV<:AbstractVector{T}}
    a = zeros(T, len)
    @inbounds a[i] = val
    return a
end

function onehot(SV::Type{SVector{N,T}}, val, i::Int, len::Int) where {N,T<:AbstractFloat}
    a = zeros(MVector{N,T})
    @inbounds a[i] = val
    return SVector(a)
end

function onehot(SV::Type{SVector{1,T}}, val, i::Int, len::Int) where {T<:AbstractFloat}
    return SV(val)
end

function DerivedVar(x::CovariantVar{T,AV}) where {T,AV}
    return DerivedVar(val(x), Gradients(system(x), onehot(AV, one(T), x.index, numvals(system(x)))))
end

include("rules.jl")

end


text:
#!/usr/bin/env bash
# -*- mode: julia -*-
#=
curdir="$(dirname "$(readlink -f "$0")")"
parentdir=$(dirname "$curdir")
echo "$parentdir"
exec julia-1.1.1 --project="$parentdir" "${BASH_SOURCE[0]}"
=#
module ROSWrapper
#TODO strongly and consistently type (esp w/ gpu to avoid Float64/32 conv)
#TODO get rid of branching in cb's
using RoboLib.Util: binarize, takeN, rangebearing2point
using CoordinateTransformations, Rotations
using RoboLib.Geom: Pose2D, project2D, Scale2D
using StaticArrays
using LinearAlgebra
using Statistics
using Random
using mushr_pf.ParticleFilter: setup_pf, update_act!, update_obs!, reset!, resample!, OccMap

# ROS Imports
using RobotOS
@rosimport sensor_msgs.msg: LaserScan
@rosimport nav_msgs.msg: Odometry
@rosimport nav_msgs.srv: GetMap
@rosimport geometry_msgs.msg: PoseStamped, PoseWithCovarianceStamped, PoseArray, Quaternion
@rosimport std_msgs.msg: Header, Float64
@rosimport ackermann_msgs.msg: AckermannDriveStamped
@rosimport vesc_msgs.msg: VescStateStamped

rostypegen(@__MODULE__)
using .sensor_msgs.msg
using .nav_msgs.msg
using .nav_msgs.srv
using .geometry_msgs.msg
using .std_msgs.msg: Header, Float64Msg
using .ackermann_msgs.msg
using .vesc_msgs.msg: VescStateStamped
# -----

include("settings.jl")
include("util.jl")

function laser_cb(msg::LaserScan, rospf::ROSPF)
    @debugtask begin
        h = msg.header
        trecv = time()
        if rospf.last_laser == -1
            rospf.last_laser = trecv

            # a wise man once said "buffers are the way to happiness"
            angles = LinRange(msg.angle_min, msg.angle_max, length(msg.ranges))
            resize!(rospf.angles, length(angles))
            resize!(rospf.laseridxs, length(angles))
            rospf.angles .= angles
        else
            dt = trecv - rospf.last_laser
            rospf.last_laser = trecv

            #println((1/dt, size())

            idxs = validindices(msg.ranges, rospf)
            angles = takeN(view(rospf.angles, idxs), rospf.nrays)
            ranges = takeN(view(msg.ranges, idxs), rospf.nrays)

            @assert length(angles) == length(ranges) == rospf.nrays

            update_obs!(rospf.pf, (angles, ranges), dt)
            resample!(rospf.pf)

            meanpose = mean(rospf.pf.belief.particles)
            publish(rospf.odom_pub, rosodometrystamped(meanpose, "map"))
            publish(rospf.pose_pub, rosposestamped(meanpose, "map"))

            # Publish a subset of the particles
            particles = view(rospf.pf.belief.particles, rand!(rospf.rand_part_idxs, eachindex(rospf.rand_part_idxs)))
            publish(rospf.particle_pub, rosposearraystamped(particles, "map"))

            # Publish the range measurements used in the reweighting step
            laserposes!(rospf, angles, ranges, meanpose)
            laserposes = rosposearraystamped(rospf.laserposes, "map")
            publish(rospf.laserpose_pub, laserposes)
        end
    end
end

function clicked_pose_cb(msg::PoseWithCovarianceStamped, rospf::ROSPF)
    @debugtask begin
        p = msg.pose.pose.position
        q = msg.pose.pose.orientation

        r2d = project2D(Quat(q.w, q.x, q.y, q.z))
        t2d = SVector{2}(p.x, p.y)
        reset!(rospf.pf, Pose2D(r2d, t2d))

        rospf.initialized = true
    end
end

function servo_cb(msg, rospf::ROSPF)
    rospf.last_servo = msg.data
end

function vesc_cb(msg, rospf::ROSPF)
    @debugtask begin
        trecv = to_sec(msg.header.stamp)
        if rospf.last_servo == -1
            return
        elseif rospf.last_vesc == -1
            rospf.last_vesc = trecv
        else
            dt = trecv - rospf.last_vesc
            rospf.last_vesc = trecv

            v = (msg.state.speed - rospf.speed2erpm_offset) / rospf.speed2erpm_gain
            delta = (rospf.last_servo - rospf.steering2servo_offset) / rospf.steering2servo_gain

            update_act!(rospf.pf, SVector(v, delta), dt)
        end
    end
end

function ackermann_cb(msg, rospf::ROSPF)
    @debugtask begin
        h = msg.header
        d = msg.drive
        trecv = time()
        if rospf.last_acker == -1
            rospf.last_acker = trecv
        else
            dt = trecv - rospf.last_acker
            rospf.last_acker = trecv

            ctrl = SVector(d.speed, d.steering_angle)
            update_act!(rospf.pf, ctrl, dt)
        end
    end
end

#TODO(cxs): nested params
#TODO(cxs): handle array/dict params
#TODO(cxs): handle unknown
function get_params(d::Dict{Symbol, Any})
    for (param_name, value) in d
        if isa(value, Dict)
            get_params(value)
        elseif isa(value, Union{Number, String})
            if has_param('~' * string(param_name))
                new_value = oftype(value, get_param('~' * string(param_name), value))
                @assert typeof(value) == typeof(new_value)
                loginfo("Found param \"$param_name\". Default value: \"$value\", new value: \"$new_value\"")
                d[param_name] = new_value
            else
                loginfo("Param \"$param_name\" not set. Using default value: \"$value\"")
            end
        end
    end
end

function process_ros_map(map_msg)
    info = map_msg.map.info

    s = info.resolution
    p = info.origin.position
    q = info.origin.orientation

    # something something row major vs column major
    occmap = reshape(map_msg.map.data', (info.width, info.height))'
    occmap = OccMap((el)->el!=0, occmap, SVector(p.x, p.y), project2D(Quat(q.w, q.x, q.y, q.z)), s)

    return occmap
end

function main()
    init_node("particle_filter")

    get_params(conf)

    loginfo("Waiting for map service...")
    wait_for_service(rosconf[:static_map_service])
    map_msg = ServiceProxy(rosconf[:static_map_service], GetMap)(GetMapRequest())
    occmap = process_ros_map(map_msg)

    pose_pub      = Publisher(rosconf[:pose_topic], PoseStamped, queue_size = 1)
    particle_pub  = Publisher(rosconf[:particle_topic], PoseArray, queue_size = 1)
    laserpose_pub = Publisher(rosconf[:laserpose_topic], PoseArray, queue_size = 1)
    odom_pub = Publisher(rosconf[:odom_topic], Odometry, queue_size = 1)

    # setup PF
    pfconf[:rng] = Random.MersenneTwister(pfconf[:rng_seed])
    pf = setup_pf(occmap, pfconf)
    rospf = ROSPF(conf, pf, pose_pub, particle_pub, odom_pub, occmap, laserpose_pub)

    pose_sub  = Subscriber(rosconf[:initialpose_topic], PoseWithCovarianceStamped, clicked_pose_cb, (rospf,), queue_size=1)

    loginfo("Waiting for initial pose...")
    while !rospf.initialized
        rossleep(Duration(0.001))
    end

    laser_sub = Subscriber(rosconf[:scan_topic], LaserScan, laser_cb, (rospf,), queue_size=1)
    servo_sub = Subscriber(rosconf[:servo_pos_topic], Float64Msg, servo_cb, (rospf,), queue_size=1)
    vesc_sub = Subscriber(rosconf[:vesc_topic], VescStateStamped, vesc_cb, (rospf,), queue_size=1)
    #ackermann_sub = Subscriber(rosconf[:ackermann_topic], AckermannDriveStamped, ackermann_cb, (rospf,), queue_size=1)

    loginfo("Initialization complete")
    d = Duration(0.001)
    while !is_shutdown() rossleep(d) end
end

end # module

using .ROSWrapper: main
if !isinteractive()
    main()
end


text:
function hypergraph_density(H, number_of_species)
    temp = [count(h -> length(h) == i, H)/binomial(number_of_species, i) for i=1:number_of_species]

    return temp
end

function hypergraph_dimension(H)
    return maximum(length.(H))
end

function get_descendents(G, h)
    return [union(h,v) for v in neighbors(G, last(h))]
end

function missing_edges(hypergraph)
    return setdiff(minimal_simplicial_complex(hypergraph), hypergraph)
end

"""
read_hypergraph(filename)
"""
function read_hypergraph(filename)
    Hp = readlines(filename)
    temp = [ vcat([parse(Int, ss) for ss in split.(h, ",")]...) for h in Hp]
    return temp
end

"""
hypergraph_subdivide(H; expansion = false)
"""
function hypergraph_subdivide(H; expansion = false)

    # Build graph of inclusion between hyperedges.
    # convention is [i,j] for (i -> j)
    G = []
    for size = reverse(2:hypergraph_dimension(H))
            nodes_s = findall(h -> length(h) == size , H) # v ∈ V such that length(H[v]) == size
            nodes_sm1 = findall(h -> length(h) < size, H) # v ∈ V such that length(H[v]) < size

            for p in nodes_s, c in nodes_sm1
                    if H[c] ⊆ H[p]
                            push!(G, [p,c])
                    end
            end
    end

    if expansion .== false # return just G
            return G
    else # expand G to the Vietoris-Rips comples

            # Expand G into K_H
            roots = setdiff(1:length(H), last.(G)) # roots of G, all nodes without incoming edges (characterized by last item of pairs [i,j])
            tops = setdiff(1:length(H), first.(G)) # tops of G, all nodes without outgoing edges

            list = [[i] for i=1:length(H)] #initialize list of hyperedges
            let temp = [[h] for h in roots]
                    @inbounds @showprogress "Subdividing: " for i =1:hypergraph_dimension(H)-1
                            #temp = [union(h, v) for h in temp for v in neighbors(G, last(h)) ]
                            temp = vcat(map(h -> get_descendents(G, h), temp)...)
                            union!(list, temp)
                            # remove all nodes that are already at the top of G
                            temp = filter(x -> .!issubset(last(x), tops), temp)
                    end
            end
            return list
    end
end

"""
betti_hypergraph_ripscomplex(H; max_dim = 3)

Computes Betti numbers of hypergraph.
Based on computing the hypergraph subdivision without expansion, which returns the inclusion graph
"""
function betti_hypergraph_ripscomplex(H; max_dim = 3)

    if length(H) == 0
        return Int.(zeros(max_dim + 1))
    elseif length(H) == 1
        return Int.(vcat(1, zeros(max_dim)))
    end

    _max_dim = min(max_dim, length(H) - 2);

    G = hypergraph_subdivide(H; expansion = false)
    # build distance matrix to give as input to Ripserer.
    M = 2.0*ones(length(H), length(H))
    [ M[g[1], g[2]] = 1 for g in G ]
    [ M[g[2], g[1]] = 1 for g in G ]
    M[diagind(M)] .= 0;

    #call Ripser
    flt = Rips(M; threshold=1.5, sparse=true)
    barcodes = ripserer(flt; dim_max = _max_dim)
    # interpret the outout. Barcodes is an array of tuples of dim = max_dim + 1.
    # All tuples where the second entry is Inf are the persistent ones.
    # Betti numbers just count the number of those persistent ones.
    betti = [ count(x -> x[2] .== Inf, barcodes[i]) for i in 1:length(barcodes)]

    padding = (max_dim+1) - length(betti);
    return vcat(betti, Int.(zeros(padding)))

    ## old function using Ripser library: https://github.com/mtsch/Ripser.jl#master
    # if length(H) == 0
    #     return zeros(max_dim + 1)
    # else
    #
    #     G = hypergraph_subdivide(H; expansion = false)
    #     # build distance matrix to give as input to Ripser.
    #     M = 2.0*ones(length(H), length(H))
    #     [ M[g[1], g[2]] = 1 for g in G ]
    #
    #     #call Ripser
    #     barcodes = ripser(M; dim_max = max_dim, threshold = 1.5)
    #     # interpret the outout. Barcodes is an array of tuples of dim = max_dim + 1.
    #     # All tuples where the second entry is Inf are the persistent ones.
    #     # Betti numbers just count the number of those persistent ones.
    #     betti = [ count(x -> x[2] .== Inf, barcodes[i]) for i in 1:length(barcodes)]
    #
    #     return betti
    # end
end

"""
minimal_simplicial_complex(hypergraph)

Computes the minimal simplicial complex containig the hypergraph H
"""
function minimal_simplicial_complex(hypergraph)
    K = Array{Int64, 1}[]
    #K = []
    for h in hypergraph
        union!(K, collect(combinations(h)))
    end
    return unique(Set, K) # delete duplicated simplices
end

"""
disassembly_hypergraph(hypergraph)

"""
function disassembly_hypergraph(hypergraph)
    M = missing_edges(hypergraph)
    temp = filter(h -> length(h) == 1, hypergraph) # hyperedges of dimension = 1 (i.e., vertices)

    return union(temp, M)
end

max_issues_repo_name: mtangulo/coexistenceholes

max_forks_repo_name: mtangulo/coexistenceholes


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

JLLWrappers.@generate_main_file_header("Rclone")
JLLWrappers.@generate_main_file("Rclone", UUID("d3707c4a-c138-5086-90bb-c65de4e50315"))
end  # module Rclone_jll


text:
# Time selection & restriction functions
abstract SnapMode
type SnapExact    <: SnapMode; end
type SnapPrevious <: SnapMode; end
type SnapNext     <: SnapMode; end
type SnapNearest  <: SnapMode; end

time2idx(sig::Signal, t::SecondT) = time2idx(sig, t, SnapExact())
function time2idx(sig::Signal, t::SecondT, ::SnapExact)
    ts = sig.time
    ts[1] <= t <= ts[end] || throw(BoundsError())
    r = searchsorted(ts, t)
    isempty(r) && throw(InexactError())
    return r[1]
end
function time2idx(sig::Signal, t::SecondT, ::SnapNearest)
    ts = sig.time
    t <= ts[1] && return 1 # No bounds errors; these are easy cases
    t >= ts[end] && return length(ts)
    i = searchsortedlast(ts, t)
    return ifelse(t - ts[i] > ts[i+1] - t, i+1, i)
end
function time2idx(sig::Signal, t::SecondT, ::SnapPrevious)
    ts = sig.time
    t < ts[1] && throw(BoundsError())
    searchsortedlast(ts, t)
end
function time2idx(sig::Signal, t::SecondT, ::SnapNext)
    ts = sig.time
    t > ts[end] && throw(BoundsError())
    searchsortedfirst(ts, t)
end

# With regular signals, we can compute and allow out-of-range indices
# TODO: "lift" these computations as FloatRange does?
_time2idxf(sig::RegularSignal, t::SecondT) = (t - sig.time[1])*samplingfreq(sig)
time2idx(sig::RegularSignal, t::SecondT, ::SnapExact)    = 1+Int(_time2idxf(sig,t))
# TODO: fix inference here
time2idx(sig::RegularSignal, t::SecondT, ::SnapNearest)  = 1+Int(round(Int, _time2idxf(sig,t)))
time2idx(sig::RegularSignal, t::SecondT, ::SnapPrevious) = 1+Int(floor(Int, _time2idxf(sig,t)))
time2idx(sig::RegularSignal, t::SecondT, ::SnapNext)     = 1+Int(ceil(Int, _time2idxf(sig,t)))

# Time restriction
before(sig::Signal, t::SecondT) = before(sig, time2idx(sig, t, SnapPrevious()))
before(sig::Signal, i::Real)    = within(sig, 1, i)

after(sig::Signal, t::SecondT) = after(sig, time2idx(sig, t, SnapNext()))
after(sig::Signal, i::Real)    = within(sig, i, length(sig.time))

within(sig::Signal, i::(Real, Real))          = within(sig, t[1], t[2])
within(sig::Signal, t::(SecondT, SecondT))    = within(sig, time2idx(sig, t[1], SnapNext()), time2idx(sig, t[2], SnapPrevious()))
within(sig::Signal, t1::SecondT, t2::SecondT) = within(sig, time2idx(sig, t1,   SnapNext()), time2idx(sig, t2,   SnapPrevious()))
function within(sig::Signal, i1::Real, i2::Real) # The real (and only) workhorse
    data = sub(sig.data, i1:i2, :)
    data = reshape(data, tuple(length(i1:i2), size(sig)...))
    Signal(sig.time[i1:i2], data)
end

# Windowing a regular signal returns a signal of a signals: one for each channel
# each with a timebase of the window size and length(at) repetitions
# Window defaults to windowing all channels
# convert seconds to relative indices ahead of time for regular signals

# Note that, given locations in times that are not exactly on a sample point,
# this *snaps* the time vector of the resulting Signal to the nearest sample.
# Were this not done, there are two things that become problematic:
#  1. Fencepost problems: If the window is exactly `m` multiples of the sampling
#     rate (or is specified in indices), the result of windowing may have m or
#     m+1 samples depending on the event location (m+1 when exactly on a sample
#     point, m otherwise). This would be very difficult to handle downstream.
#  2. We could no longer combine multiple repetitions into one Signal object as
#     the times of each repetition would vary by as much as ±0.5dt. The only
#     sensible way to combine such signals is to upsample them all to a common
#     timebase. This is something that the user could do ahead of time if they
#     are concerned about jitter on the scale of ±0.5dt.
function window{S<:SecondT}(sig::RegularSignal, at::AbstractVector{S}, within::(Real, Real))
    window(sig, [time2idx(sig, a, SnapNearest()) for a in at], within)
end
# Can't use within::(Union(SecondT,Real), Union(SecondT,Real)) due to ambiguity
function window{S<:SecondT}(sig::RegularSignal, at::AbstractVector{S}, within::(SecondT, SecondT))
    window(sig, [time2idx(sig, a, SnapNearest()) for a in at], within)
end
function window{R<:Real}(sig::RegularSignal, at::AbstractVector{R}, within::(SecondT, SecondT))
    window(sig, at, (time2idx(sig, within[1], SnapNext()), time2idx(sig, within[2], SnapPrevious())))
end
function window{R<:Real,R2<:Real}(sig::RegularSignal, at::AbstractVector{R}, within::(R2, R2))
    r = within[1]:within[2]
    dt = float(samplingrate(sig)) # TODO: Math with SIUnits is very annoying.
    t = inseconds((within[1]*dt):dt:(within[2]*dt))

    vi = similar(sig.data, tuple(length(r), length(at), size(sig)...))
    for chan=1:length(sig)
        for rep in 1:length(at)
            vi[:, rep, chan] = sig.data[r+at[rep], chan]
        end
    end
    Signal(sig.time[at], Signal(t, vi))
end

# Windowing an irregular signal cannot aggregate the Signals together into a
# common time-base across all the repetitions. We could aggregate *channels*
# together within each repetition, but that is the effective transpose of what
# we need. Furthermore, specifying indices vs. time can behave very differently
# because we don't snap to the nearest sample like RegularSignals. So there are
# four different cases, each constructing a matrix of one-channel Signals:
function window{S<:SecondT}(sig::Signal, at::AbstractVector{S}, within::(Real, Real))
    # Set number of indexes within each window, about time
    r1, r2 = within
    r = r1:r2
    cs = [(i = time2idx(sig, t, SnapNearest()); signal(time(sig, i+r) - t, c[i+r])) for t in at, c in sig]
    signal(at, cs)
end
function window{R<:Real}(sig::Signal, at::AbstractVector{R}, within::(Real, Real))
    # Set number of indexes within each window, about indices
    r1, r2 = within
    r = r1:r2
    cs = [signal(time(sig, i+r) - time(sig, i), c[i+r]) for i in at, c in sig]
    signal(time(sig, at), cs)
end
function window{S<:SecondT}(sig::Signal, at::AbstractVector{S}, within::(SecondT, SecondT))
    # Time windows, about times
    t1, t2 = within
    cs = [signal(timewithin(sig, t+t1, t+t2)-t, chanwithin(sig, c, t+t1, t+t2)) for t in at, c in 1:length(sig)]
    signal(at, cs)
end
function window{R<:Real}(sig::Signal, at::AbstractVector{R}, within::(SecondT, SecondT))
    # time windows, about indices
    t1, t2 = within
    cs = [(t = time(sig, i); signal(timewithin(sig, t+t1, t+t2)-t, chanwithin(sig, c, t+t1, t+t2))) for i in at, c in 1:length(sig)]
    signal(time(sig, at), cs)
end


text:
"""
    _outlength(A, I, J)
    _outlength(u, I)
Determine the size of the output for an operation like extract or range-based indexing.
"""
function _outlength(A, I, J)
    if I == ALL
        Ilen = size(A, 1)
    else
        Ilen = length(I)
    end
    if J == ALL
        Jlen = size(A, 2)
    else
        Jlen = length(J)
    end
    return Ilen, Jlen
end

function _outlength(u, I)
    if I == ALL
        wlen = size(u)
    else
        wlen = length(I)
    end
    return wlen
end

"""
    extract!(C::GBMatrix, A::GBMatOrTranspose, I, J; kwargs...)::GBMatrix
    extract!(C::GBVector, A::GBVector, I; kwargs...)::GBVector

Extract a submatrix or subvector from `A` into `C`.

# Arguments
- `C::Union{GBVector, GBMatrix}`: the submatrix or subvector extracted from `A`.
- `A::GBArray`
- `I` and `J`: A colon, scalar, vector, or range indexing A.

# Keywords
- `mask::Union{Nothing, GBArray} = nothing`: mask where
    `size(M) == (max(I), max(J))`.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc::Union{Nothing, Descriptor} = nothing`

# Returns
- `GBMatrix` or `GBVector`: the modified array `C`, now containing the matrix `A[I, J]` or
    `A[I]` for a vector.

# Throws
- `GrB_DIMENSION_MISMATCH`: If `size(C) != (max(I), max(J))` or `size(C) != size(mask)`.
"""
extract!

function extract!(
    C::GBMatrix, A::GBMatOrTranspose, I, J;
    mask = nothing, accum = nothing, desc = nothing
)
    I, ni = idx(I)
    J, nj = idx(J)
    I isa Number && (I = UInt64[I])
    J isa Number && (J = UInt64[J])
    mask === nothing && (mask = C_NULL)
    desc = _handledescriptor(desc; in1 = A)
    libgb.GrB_Matrix_extract(C, mask, getaccum(accum, eltype(C)), parent(A), I, ni, J, nj, desc)
    return C
end

function extract!(
    C::GBMatrix, A::GBMatOrTranspose, ::Colon, J;
    mask = nothing, accum = nothing, desc = nothing
)
    return extract!(C, A, ALL, J; mask, accum, desc)
end

function extract!(
    C::GBMatrix, A::GBMatOrTranspose, I, ::Colon;
    mask = nothing, accum = nothing, desc = nothing
)
    return extract!(C, A, I, ALL; mask, accum, desc)
end

function extract!(
    C::GBMatrix, A::GBMatOrTranspose, ::Colon, ::Colon;
    mask = nothing, accum = nothing, desc = nothing
)
    return extract!(C, A, ALL, ALL; mask, accum, desc)
end

"""
    extract(A::GBMatOrTranspose, I, J; kwargs...)::GBMatrix
    extract(A::GBVector, I; kwargs...)::GBVector

Extract a submatrix or subvector from `A`

# Arguments
- `A::GBArray`: the array being indexed.
- `I` and `J`: A colon, scalar, vector, or range indexing A.

# Keywords
- `mask::Union{Nothing, GBArray} = nothing`: mask where
    `size(M) == (max(I), max(J))`.
- `accum::Union{Nothing, AbstractBinaryOp} = nothing`: binary accumulator operation
    where `C[i,j] = accum(C[i,j], T[i,j])` where T is the result of this function before accum is applied.
- `desc::Descriptor = nothing`

# Returns
- `GBMatrix`: the submatrix `A[I, J]`.

# Throws
- `GrB_DIMENSION_MISMATCH`: If `(max(I), max(J)) != size(mask)`.
"""
extract

function extract(
    A::GBMatOrTranspose, I, J;
    mask = nothing, accum = nothing, desc = nothing
)
    Ilen, Jlen = _outlength(A, I, J)
    C = similar(A, Ilen, Jlen)
    return extract!(C, A, I, J; mask, accum, desc)
end

function extract(
    A::GBMatOrTranspose, ::Colon, J;
    mask = nothing, accum = nothing, desc = nothing
)
    return extract(A, ALL, J; mask, accum, desc)
end

function extract(
    A::GBMatOrTranspose, I, ::Colon;
    mask = nothing, accum = nothing, desc = nothing
)
    return extract(A, I, ALL; mask, accum, desc)
end

function extract(
    A::GBMatOrTranspose, ::Colon, ::Colon;
    mask = nothing, accum = nothing, desc = nothing
)
    return extract(A, ALL, ALL; mask, accum, desc)
end

function Base.getindex(
    A::GBMatOrTranspose, ::Colon, j;
    mask = nothing, accum = nothing, desc = nothing
)
    return extract(A, ALL, j; mask, accum, desc)
end

function extract!(
    w::GBVector, u::GBVector, I;
    mask = nothing, accum = nothing, desc = nothing
)
    I, ni = idx(I)
    desc = _handledescriptor(desc)
    mask === nothing && (mask = C_NULL)
    libgb.GrB_Matrix_extract(w, mask, getaccum(accum, eltype(w)), u, I, ni, UInt64[1], 1, desc)
    return w
end

function extract!(
    w::GBVector, u::GBVector, ::Colon;
    mask = nothing, accum = nothing, desc = nothing
)
    return extract!(w, u, ALL; mask, accum, desc)
end

function extract(
    u::GBVector, I;
    mask = nothing, accum = nothing, desc = nothing
)
    wlen = _outlength(u, I)
    w = similar(u, wlen)
    return extract!(w, u, I; mask, accum, desc)
end

function extract(u::GBVector, ::Colon; mask = nothing, accum = nothing, desc = nothing)
    extract(u, ALL; mask, accum, desc)
end


text:
function param_values()::Vector{Float64}
    p::Vector{Float64} = zeros(C.NUM)

    p[C.InhCDK] = 2.0
    p[C.CycBCdk1T] = 8.1808
    p[C.CycACdk2T] = 1.0000
    p[C.PP1T] = 1.0000
    p[C.kapp1] = 0.0115
    p[C.kapp1a] = 0.7054
    p[C.kipp1] = 0.0018
    p[C.kipp1C] = 0.7549
    p[C.kPP1Gw] = 18.4724
    p[C.ENSAtot] = 1.0000
    p[C.B55tot] = 0.2500
    p[C.SubT] = 1.0000
    p[C.kass] = 617.2807
    p[C.kdis] = 0.0088
    p[C.kcatB55] = 1.0338
    p[C.kGwENSA] = 20.8811
    p[C.kppxGwl] = 0.1560
    p[C.kcBc1Sub] = 0.0080
    p[C.kcBc1G] = 0.2393
    p[C.Gwtot] = 1.0000
    p[C.kB55G] = 496.5636
    p[C.kB55Sub] = 0.0593
    p[C.kcAc2G] = 0.1916
    p[C.k25S] = 0.0050
    p[C.k25F] = 0.9411
    p[C.kweeS] = 0.0050
    p[C.kweeF] = 47.2937
    p[C.kcBc1W1] = 1.3132
    p[C.kcBc125] = 1.3132
    p[C.kppxY15] = 0.0050
    p[C.kcAc2W1] = 0.1096
    p[C.kcAc225] = 0.1096
    p[C.kB55W1] = 0.5511
    p[C.kB5525] = 0.5511
    p[C.Kd] = 0.025

    return p
end

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
############################################################################
### test sorting
############################################################################

sort!(partable)

model_ml_sorted = Sem(
    specification = partable,
    data = dat
)

@testset "graph sorting" begin
    @test model_ml_sorted.imply.I_A isa LowerTriangular
end

@testset "ml_solution_sorted" begin
    solution_ml_sorted = sem_fit(model_ml_sorted)
    update_estimate!(partable, solution_ml_sorted)
    @test SEM.compare_estimates(par_ml, partable, 0.01)
end


text:
include("./imports.jl")
@testset ExtendedTestSet "calciumscoring" begin
    ii1 = [1 2 0 1
           1 5 0 1]
    ii2 = [1 2 0 1
           1 5 0 1]
    img = cat(ii1, ii2, dims=3)

    spc = [0.625, 0.625, 0.5]

    mm1 = [1 0 0 1
           1 1 0 1]
    mm2 = [1 0 0 1
           1 1 0 1]
    msk = cat(mm1, mm2, dims=3)
    @testset ExtendedTestSet "calciumscoring" begin
        answer = (0.6510416666666667, 1.953125)
        test = compute_calcium_scores(img, spc, msk)
        @test (test[1] ≈ answer[1]) && (test[2] ≈ answer[2])
    end
    @testset ExtendedTestSet "calciumscoring" begin
        answer = (0.2604166666666667, 0.78125)
        test = compute_calcium_scores(img, spc, msk, 0, 1)
        @test (test[1] ≈ answer[1]) && (test[2] ≈ answer[2])
    end
end


text:
println("Dj) Test of the MultiPhotonDeExcitation module with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-MultiPhotonDeExcitation.sum")
setDefaults("method: continuum, Galerkin")           ## setDefaults("method: continuum, Galerkin")  "method: continuum, asymptotic Coulomb"
                                                     ## setDefaults("method: normalization, Ong-Russek")
setDefaults("method: normalization, pure sine")      ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")
setDefaults("unit: rate", "a.u.")

## grid = Radial.Grid(Radial.Grid(false), rnt = 4.0e-6, h = 5.0e-2, hp = 1.0e-2, rbox = 10.0)
grid = Radial.Grid(Radial.Grid(true), rnt = 4.0e-6, h = 5.0e-2, rbox = 10.0)
name = "two-photon absorption or emission of lithium-like xenon"

if  false
    # Green function for two-photon absorption or emission of lithium-like xenon
    refConfigs       = [Configuration("[He] 2p")]
    levelSymmetries  = [LevelSymmetry(1//2, Basics.minus), LevelSymmetry(3//2, Basics.minus),
                        LevelSymmetry(1//2, Basics.plus),  LevelSymmetry(3//2, Basics.plus)]
    greenSettings    = GreenSettings(3, [0, 1], 0.01, true, LevelSelection())

    greenRep         = Representation(name, Nuclear.Model(54.), grid, refConfigs,
                                      GreenExpansion( AtomicState.DampedSpaceCI(), Basics.DeExciteSingleElectron(), levelSymmetries, 3, greenSettings) )
                                   ## GreenExpansion( AtomicState.CoreSpaceCI(), Basics.DeExciteSingleElectron(),
                                   ## GreenExpansion( AtomicState.DampedSpaceCI(), Basics.DeExciteSingleElectron()
    greenOut         = generate(greenRep, output=true)
    mpGreen          = greenOut["Green channels"]

elseif  false
    #
    println("Calculate two-photon absorption cross sections by monochromatic and equally-polarized photons from the same beam")
    mpSettings       = MultiPhotonDeExcitation.Settings(MultiPhotonDeExcitation.TwoPhotonAbsorptionMonochromatic(),
                                                        [E1,M1], [UseCoulomb], mpGreen, 0, true, LineSelection())
    wa = Atomic.Computation(Atomic.Computation(), name=name, grid=grid, nuclearModel=Nuclear.Model(54.),
                            initialConfigs= [Configuration("1s^2 2s")],
                            finalConfigs  = [Configuration("1s^2 3s")],
                            process = MultiPhotonDE(), processSettings= mpSettings  )
    wb = perform(wa)

elseif  true
    #
    println("Calculate two-photon emission cross sections: energy-differential + total cross sections")
    mpSettings       = MultiPhotonDeExcitation.Settings(MultiPhotonDeExcitation.TwoPhotonEmission(),
                                                        [E1,M1], [UseCoulomb], mpGreen, 4, true, LineSelection())

    wc = Atomic.Computation(Atomic.Computation(), name=name, grid=grid, nuclearModel=Nuclear.Model(54.),
                            initialConfigs= [Configuration("1s^2 3s")],
                            finalConfigs  = [Configuration("1s^2 2s")],
                            process = MultiPhotonDE(), processSettings= mpSettings    )
    wd = perform(wc)

end

setDefaults("print summary: close", "")


text:
using JuMP
using GLPKMathProgInterface
using Ipopt

n = 2;
lowRes = -4;
hiRes  = 4;
m = Model()
@defVar(m, lowRes <= α[i=1:n]  <=hiRes)
@defVar(m, lowRes <= β1[i=1:n] <=hiRes)
@defVar(m, lowRes <= β2[i=1:n] <=hiRes)
@defVar(m, lowRes <= β3[i=1:n] <=hiRes)
@defVar(m, lowRes <= β4[i=1:n] <=hiRes)

@defVar(m, lowRes <= A011 <= hiRes)
@defVar(m, lowRes <= A012 <= hiRes)
@defVar(m, lowRes <= A021 <= hiRes)
@defVar(m, lowRes <= A022 <= hiRes)

@defVar(m, lowRes <= A111 <= hiRes)
@defVar(m, lowRes <= A112 <= hiRes)
@defVar(m, lowRes <= A121 <= hiRes)
@defVar(m, lowRes <= A122 <= hiRes)

@defVar(m, lowRes <= B011 <= hiRes)
@defVar(m, lowRes <= B012 <= hiRes)
@defVar(m, lowRes <= B021 <= hiRes)
@defVar(m, lowRes <= B022 <= hiRes)

@defVar(m, lowRes <= B111 <= hiRes)
@defVar(m, lowRes <= B112 <= hiRes)
@defVar(m, lowRes <= B121 <= hiRes)
@defVar(m, lowRes <= B122 <= hiRes)

A0 = [A011 A012
     A021 A022]
A1 = [A111 A112
     A121 A122]
B0 = [B011 B012
      B021 B022]
B1 = [B111 B112
      B121 B122]
e = [1 1]'

@setObjective(m, Max, 0)
@addConstraints(m,begin
                  sum{α[i], i=1:n}==1
                  sum{β1[i], i=1:n}==1
                  sum{β1[i], i=1:n}==1
                  sum{β2[i], i=1:n}==0
                  sum{β3[i], i=1:n}==0
                  sum{β4[i], i=1:n}==0
                  sum{β4[i], i=1:n}==0
                end)
@addNLConstraints(m,begin
                  B111*β1[1]+B112*β1[1]+B121*β1[2]+B122*β1[2]==0
                  B111*β2[1]+B112*β2[1]+B121*β2[2]+B122*β2[2]==1
                  B111*β3[1]+B112*β3[1]+B121*β3[2]+B122*β3[2]==0
                  B111*β4[1]+B112*β4[1]+B121*β4[2]+B122*β4[2]==0
                  A011*α[1] +A012*α[1] +A021*α[2] +A022*α[2] ==1//2
                  B011*α[1] +B012*α[1] +B021*α[2] +B022*α[2] ==1
                  A111*β1[1]+A112*β1[1]+A121*β1[2]+A122*β1[2]==1
                  A111*β2[1]+A112*β2[1]+A121*β2[2]+A122*β2[2]==0
                  A111*β3[1]+A112*β3[1]+A121*β3[2]+A122*β3[2]==-1
                  A111*β4[1]+A112*β4[1]+A121*β4[2]+A122*β4[2]==0
                end)

@addNLConstraint(m,α[1] * (B011+B012)*(B011+B012)  + α[2] * (B021+B022)*(B021+B022)  == 3//2)
@addNLConstraint(m,β1[1] * (B111+B112)*(B111+B112)  + β1[2] * (B121+B122)*(B121+B122)  == 1)
@addNLConstraint(m,β2[1] * (B111+B112)*(B111+B112)  + β2[2] * (B121+B122)*(B121+B122)  == 0)
@addNLConstraint(m,β3[1] * (B111+B112)*(B111+B112)  + β3[2] * (B121+B122)*(B121+B122)  == -1)
@addNLConstraint(m,β4[1] * (B111+B112)*(B111+B112)  + β4[2] * (B121+B122)*(B121+B122)  == 2)
@addNLConstraint(m,β1[1] * (B111*(B111+B112) + B112*(B121+B122))  + β1[2] * (B121*(B111+B112) + B122*(B121+B122))  == 0)
@addNLConstraint(m,β2[1] * (B111*(B111+B112) + B112*(B121+B122))  + β2[2] * (B121*(B111+B112) + B122*(B121+B122))  == 0)
@addNLConstraint(m,β3[1] * (B111*(B111+B112) + B112*(B121+B122))  + β3[2] * (B121*(B111+B112) + B122*(B121+B122))  == 0)
@addNLConstraint(m,β4[1] * (B111*(B111+B112) + B112*(B121+B122))  + β4[2] * (B121*(B111+B112) + B122*(B121+B122))  == 1)

@addNLConstraint(m,(1/2)*((A111*(B011+B012) + A112*(B021+B022))*β1[1] +
      (A121*(B011+B012) + A122*(B021+B022))β1[2]) +
      (1/3)*((A111*(B011+B012) + A112*(B021+B022))β3[1] +
      (A121*(B011+B012) + A122*(B021+B022))β3[2]) == 0)

print(m)

status = solve(m)

println("Objective value: ", getObjectiveValue(m))
println("β1 = ", getValue(β1))
println("β2 = ", getValue(β2))
println("β3 = ", getValue(β3))
println("β4 = ", getValue(β4))


text:
for script in readdir("./script/quick_adjacency_decomposition")
    # run only julia files in ./script/quick_adjacency_decomposition/ directory
    if occursin(r"^.*\.jl$", script) && !occursin(r"^run_all\.jl$", script)
        println("./quick_adjacency_decomposition/$script")
        @time include("./$script")
    end
end