text:
#### Exceptions ####

"""
Exception indicating that the specified method is not implemented for an interface.
"""
struct UnimplementedError <: Exception
    msg::String
end

#### Internal Types and Methods ####

struct Indexer
    next_index::Dict{NodeID, Index}
    indices::Dict{NodeID, Dict{String, Index}}
end

function Indexer()::Indexer
    idxr = Indexer(Dict{NodeID, Index}(),
                   Dict{NodeID, Dict{String,Index}}())
    return idxr
end

function _retrieve_and_advance_index(idxr::Indexer, nid::NodeID)::Index
    if !haskey(idxr.next_index, nid)
        idxr.next_index[nid] = Index(zero(INDEX))
    end
    idx = idxr.next_index[nid]
    idxr.next_index[nid] = _increment(idx)
    return idx
end

function index(idxr::Indexer, nid::NodeID, name::String)::Index

    if !haskey(idxr.indices, nid)
        idxr.indices[nid] = Dict{String, Index}()
    end
    node_vars = idxr.indices[nid]

    if haskey(node_vars, name)

        idx = node_vars[name]

    else

        idx = _retrieve_and_advance_index(idxr, nid)
        node_vars[name] = idx

    end

    return idx
end

struct VariableData
    name::String
    xhat_id::XhatID
end

mutable struct ProblemData
    obj::Float64
    sts::MOI.TerminationStatusCode
    time::Float64
    lb_obj::Float64
    lb_sts::MOI.TerminationStatusCode
    lb_time::Float64
end

ProblemData() = ProblemData(0.0, MOI.OPTIMIZE_NOT_CALLED, 0.0,
                            0.0, MOI.OPTIMIZE_NOT_CALLED, 0.0)

struct ScenarioInfo
    pid::Int
    prob::Float64
    branch_vars::Dict{VariableID, Float64}
    leaf_vars::Dict{VariableID, Float64}
    w_vars::Dict{VariableID, Float64}
    xhat_vars::Dict{VariableID, Float64}
    problem_data::ProblemData
end

function ScenarioInfo(pid::Int,
                      prob::Float64,
                      branch_ids::Set{VariableID},
                      leaf_ids::Set{VariableID}
                      )::ScenarioInfo

    branch_map = Dict{VariableID, Float64}()
    w_dict = Dict{VariableID, Float64}()
    x_dict = Dict{VariableID, Float64}()
    for vid in branch_ids
        branch_map[vid] = 0.0
        w_dict[vid] = 0.0
        x_dict[vid] = 0.0
    end

    leaf_map = Dict{VariableID, Float64}(vid => 0.0 for vid in leaf_ids)

    return ScenarioInfo(pid,
                        prob,
                        branch_map,
                        leaf_map,
                        w_dict,
                        x_dict,
                        ProblemData())
end

function create_ph_dicts(sinfo::ScenarioInfo)
    return (sinfo.w_vars, sinfo.xhat_vars)
end

function objective_value(sinfo::ScenarioInfo)::Float64
    return sinfo.problem_data.obj
end

function retrieve_variable_value(sinfo::ScenarioInfo, vid::VariableID)::Float64
    if haskey(sinfo.branch_vars, vid)
        vi = sinfo.branch_vars[vid]
    else
        vi = sinfo.leaf_vars[vid]
    end
    return vi
end

struct PHIterate
    xhat::Dict{XhatID,Float64}
    x::Dict{VariableID,Float64}
    w::Dict{VariableID,Float64}
end

struct PHLowerBound
    lower_bound::Float64
    gap::Float64
    rel_gap::Float64
end

struct PHResidual
    abs_res::Float64
    rel_res::Float64
    xhat_sq::Float64
    x_sq::Float64
end

struct PHHistory
    iterates::Dict{Int, PHIterate}
    residuals::Dict{Int, PHResidual}
    lower_bounds::Dict{Int, PHLowerBound}
end

function PHHistory()
    return PHHistory(Dict{Int, PHIterate}(),
                     Dict{Int, PHResidual}(),
                     Dict{Int, PHLowerBound}(),
                     )
end

function _save_iterate(phh::PHHistory,
                       iter::Int,
                       phi::PHIterate
                       )::Nothing
    phh.iterates[iter] = phi
    return
end

function _save_lower_bound(phh::PHHistory,
                           iter::Int,
                           phlb::PHLowerBound
                           )::Nothing
    phh.lower_bounds[iter] = phlb
    return
end

function _save_residual(phh::PHHistory, iter::Int, res::PHResidual)::Nothing
    phh.residuals[iter] = res
    return
end

#### User Facing Types ####

"""
Struct for user callbacks.

**Fields**

* `name::String` : User's name for the callback. Defaults to `string(h)`.
* `h::Function` : Callback function. See notes below for calling signature.
* `initialize::Function` : Function to initialize `ext` *after* subproblem creation has occurred.
* `ext::Dict{Symbol,Any}` : Dictionary to store data between callback calls or needed parameters.

The callback function `h` must have the signature
    `h(ext::Dict{Symbol,Any}, phd::PHData, winf::WorkerInf, niter::Int)::Bool`
where `ext` is the same dictionary given to the `Callback` constructor, `phd` is the standard PH data structure (see `PHData`), `winf` is used for communicating with subproblems (see `apply_to_subproblem`)  and `niter` is the current iteration. The callback may return `false` to stop PH.

The `initialize` function must have the signature
    `initialize(ext::Dict{Symbol,Any}, phd::PHData)`
where `ext` is the same dictionary given to the `Callback` constructor and `phd` is the standard PH data structure (see `PHData`).
"""
struct Callback
    name::String
    h::Function
    initialize::Function
    ext::Dict{Symbol,Any}
end

function Base.show(io::IO, cb::Callback)
    print(io, cb.name)
    return
end

"""
Type representing a consensus variable.

The following functions are available to the user to interact with consensus variables
* [`is_integer`](@ref)
* [`value`](@ref)
* [`variables`](@ref)
"""
mutable struct HatVariable
    value::Float64 # Current value of variable
    vars::Set{VariableID} # All nonhat variable ids that contribute to this variable
    is_integer::Bool # Flag indicating that this variable is an integer (includes binary)
end

HatVariable(is_int::Bool)::HatVariable = HatVariable(0.0, Set{VariableID}(),is_int)

function HatVariable(val::Float64,vid::VariableID, is_int::Bool)
    return HatVariable(val, Set{VariableID}([vid]), is_int)
end

## Primary PH Data Structure ##

"""
Data structure used to store information and results for a stochastic programming problem.

See the following functions make use of this object:
* [`apply_to_subproblem`](@ref)
* [`branch_value`](@ref)
* [`consensus_variables`](@ref)
* [`convert_to_variable_ids`](@ref)
* [`convert_to_xhat_id`](@ref)
* [`get_callback`](@ref)
* [`get_callback_ext`](@ref)
* [`is_leaf`](@ref)
* [`name`](@ref)
* [`probability`](@ref)
* [`scenario_bundle`](@ref)
* [`scenarios`](@ref)
* [`stage_id`](@ref)
* [`value`](@ref)
* [`w_value`](@ref)
* [`xhat_value`](@ref)

The following post solution functions are also available:
* [`leaf_value`](@ref)
* [`retrieve_soln`](@ref)
* [`retrieve_aug_obj_value`](@ref)
* [`retrieve_obj_value`](@ref)
* [`retrieve_no_hats`](@ref)
* [`retrieve_w`](@ref)

If the corresponding save options are enabled, the saved terms may be accessed with one of the following:
* [`lower_bounds`](@ref)
* [`residuals`](@ref)
* [`retrieve_xhat_history`](@ref)
* [`retrieve_no_hat_history`](@ref)
* [`retrieve_w_history`](@ref)
"""
struct PHData
    r::AbstractPenaltyParameter
    scenario_tree::ScenarioTree
    scenario_map::Dict{ScenarioID, ScenarioInfo}
    callbacks::Vector{Callback}
    xhat::Dict{XhatID, HatVariable}
    variable_data::Dict{VariableID, VariableData}
    history::PHHistory
    time_info::TimerOutputs.TimerOutput
end

function PHData(r::AbstractPenaltyParameter,
                tree::ScenarioTree,
                scen_proc_map::Dict{Int, Set{ScenarioID}},
                var_map::Dict{ScenarioID, Dict{VariableID, VariableInfo}},
                time_out::TimerOutputs.TimerOutput
                )::PHData

    var_data = Dict{VariableID,VariableData}()
    xhat_dict = Dict{XhatID, HatVariable}()
    idxr = Indexer()

    scenario_map = Dict{ScenarioID, ScenarioInfo}()
    for (pid, scenarios) in pairs(scen_proc_map)
        for scen in scenarios

            branch_ids = Set{VariableID}()
            leaf_ids = Set{VariableID}()

            for (vid, vinfo) in pairs(var_map[scen])

                vnode = node(tree, vid.scenario, vid.stage)

                if isnothing(vnode)
                    error("Unable to locate scenario tree node for variable '$(vinfo.name)' occuring in scenario $(vid.scenario) and stage $(vid.stage).")
                end

                idx = index(idxr, vnode.id, vinfo.name)
                xhid = XhatID(vnode.id, idx)
                vdata = VariableData(vinfo.name, xhid)
                var_data[vid] = vdata

                if is_leaf(vnode)

                    push!(leaf_ids, vid)

                else

                    push!(branch_ids, vid)

                    if haskey(xhat_dict, xhid)
                        if is_integer(xhat_dict[xhid]) != vinfo.is_integer
                            error("Variable '$(vinfo.name)' must be integer or non-integer in all scenarios in which it is used.")
                        end
                    else
                        xhat_dict[xhid] = HatVariable(vinfo.is_integer)
                    end
                    add_variable(xhat_dict[xhid], vid)

                end
            end

            scenario_map[scen] = ScenarioInfo(pid,
                                              tree.prob_map[scen],
                                              branch_ids,
                                              leaf_ids,
                                              )
        end
    end

    return PHData(r,
                  tree,
                  scenario_map,
                  Vector{Callback}(),
                  xhat_dict,
                  var_data,
                  PHHistory(),
                  time_out,
                  )
end

function Base.show(io::IO, phd::PHData)
    nscen = length(scenarios(phd))
    print(io, "PH structure for a stochastic program with $(nscen) scenarios.")
    return
end

#### User Facing Functions ####

## Callback Functions ##

"""
    Callback(f::Function)

Creates a `Callback` structure for function `f`.
"""
function Callback(f::Function)
    return Callback(string(f), f, (::Dict{Symbol,Any},::PHData)->(), Dict{Symbol,Any}())
end

"""
    Callback(f::Function, ext::Dict{Symbol,Any})

Creates a `Callback` structure for function `f` with the external data dictionary `ext`.
"""
function Callback(f::Function, ext::Dict{Symbol,Any})
    return Callback(string(f), f, (::Dict{Symbol,Any},::PHData)->(), ext)
end

"""
    Callback(f::Function, initialize::Function)

Creates a `Callback` structure for function `f` with initializer `initialize`.
"""
function Callback(f::Function, initialize::Function)
    return Callback(string(f), f, initialize, Dict{Symbol,Any}())
end

"""
    Callback(name::String, f::Function, ext::Dict{Symbol,Any})

Creates a `Callback` structure for function `f` with the name `name` and the external data dictionary `ext`.
"""
function Callback(name::String, f::Function, ext::Dict{Symbol,Any})
    return Callback(name, f, (::Dict{Symbol,Any},::PHData)->(), ext)
end

"""
    Callback(f::Function, initialize::Function, ext::Dict{Symbol,Any})

Creates a `Callback` structure for function `f` with the external data dictionary `ext` which will be initialized with `initialize`.
"""
function Callback(f::Function, initialize::Function, ext::Dict{Symbol,Any})
    return Callback(string(f), f, initialize, ext)
end

"""
    cb(f::Function)
    cb(f::Function, ext::Dict{Symbol,Any})
    cb(f::Function, initialize::Function)
    cb(f::Function, initialize::Function, ext::Dict{Symbol,Any})
    cb(name::String, f::Function, ext::Dict{Symbol,Any})

Shorthand for [`Callback`](@ref) functions with the same signature.
"""
cb(f::Function) = Callback(f)
cb(f::Function, ext::Dict{Symbol,Any}) = Callback(f, ext)
cb(f::Function, initialize::Function) = Callback(f, initialize)
cb(f::Function, initialize::Function, ext::Dict{Symbol,Any}) = Callback(f, initialize, ext)
cb(name::String, f::Function, ext::Dict{Symbol,Any}) = Callback(name, f, ext)
cb(name::String, f::Function, initialize::Function, ext::Dict{Symbol,Any}) = Callback(name, f, initialize, ext)

## Consensus Variable Functions ##

function add_variable(a::HatVariable, vid::VariableID)
    push!(a.vars, vid)
    return
end

"""
    is_integer(a::HatVariable)::Bool

Returns true if the consensus variable is an integer variable. The consensus variable is an integer if the contributing subproblem variables are all integer variables.
"""
function is_integer(a::HatVariable)::Bool
    return a.is_integer
end

"""
    set_value(a::HatVariable, v::Float64)

Sets the current value of `a` to `v`.
"""
function set_value(a::HatVariable, v::Float64)::Nothing
    a.value = v
    return
end

"""
   value(a::HatVariable)::Float64

Returns the current value of `a`.
"""
function value(a::HatVariable)::Float64
    return a.value
end

"""
    variables(a::HatVariable)::Set{VariableID}

Returns the variable ids for all subproblem variables contributing to this variable.
"""
function variables(a::HatVariable)::Set{VariableID}
    return a.vars
end

## PHData Interaction Functions ##
# NOTE: Additional, more complicated functions are in utils.jl

"""
    add_callback(phd::PHData,
                 cb::Callback
                 )::Nothing

Adds the callback `cb` to the given PH problem.
"""
function add_callback(phd::PHData,
                      cb::Callback
                      )::Nothing

    push!(phd.callbacks, cb)
    cb.initialize(cb.ext, phd)

    return
end

"""
    apply_to_subproblem(to_apply::Function,
                        phd::PHData,
                        winf::WorkerInf,
                        scid::ScenarioID,
                        args::Tuple=(),
                        kwargs::NamedTuple=NamedTuple(),
                        )

Applies the function `to_apply` to the subproblem with scenario id `scid`.
"""
function apply_to_subproblem(to_apply::Function,
                             phd::PHData,
                             winf::WorkerInf,
                             scid::ScenarioID,
                             args::Tuple=(),
                             kwargs::NamedTuple=NamedTuple(),
                             )

    _send_message(winf,
                  phd.scenario_map[scid].pid,
                  SubproblemAction(scid,
                                   to_apply,
                                   args,
                                   kwargs)
                  )

    return
end

"""
    branch_value(phd::PHData, vid::VariableID)::Float64

Returns the value of the variable associated with `vid`. Must be a branch variable.

See also: [`leaf_value`](@ref), [`value`](@ref)
"""
function branch_value(phd::PHData, vid::VariableID)::Float64
    return phd.scenario_map[scenario(vid)].branch_vars[vid]
end

"""
    branch_value(phd::PHData, scen::ScenarioID, stage::StageID, idx::Index)::Float64

Returns the value of the variable associated with with scenario `scen`, stage `stage` and index `idx`. Must be a branch variable.

See also: [`leaf_value`](@ref), [`value`](@ref)
"""
function branch_value(phd::PHData, scen::ScenarioID, stage::StageID, idx::Index)::Float64
    return branch_value(phd, VariableID(scen, stage, idx))
end

"""
    consensus_variables(phd::PHData)::Dict{XhatID,HatVariable}

Returns the collection of consensus variables for the problem.
"""
function consensus_variables(phd::PHData)::Dict{XhatID,HatVariable}
    return phd.xhat
end

"""
    convert_to_variable_ids(phd::PHData, xid::XhatID)::Set{VariableID}

Convert the given consensus variable id to the contributing individual subproblem variable ids.

**Arguments**

* `phd::PHData` : PH data structure for the corresponding problem
* `xid::XhatID` : consensus variable id to convert
"""
function convert_to_variable_ids(phd::PHData, xid::XhatID)::Set{VariableID}
    return variables(phd.xhat[xid])
end

"""
    convert_to_xhat_id(phd::PHData, vid::VariableID)::XhatID

Convert the given `VariableID` to the consensus variable id (`XhatID`).

**Arguments**

* `phd::PHData` : PH data structure for the corresponding problem
* `vid::VariableID` : variable id to convert
"""
function convert_to_xhat_id(phd::PHData, vid::VariableID)::XhatID
    return phd.variable_data[vid].xhat_id
end

"""
    get_callback(phd::PHData, name::String)::Callback

Retrieve the callback with name `name`.
"""
function get_callback(phd::PHData, name::String)::Callback
    return_cb = nothing
    for cb in phd.callbacks
        if cb.name == name
            return_cb = cb
        end
    end

    if isnothing(return_cb)
        error("Unable to find callback $name.")
    end

    return return_cb
end

"""
    get_callback_ext(phd::PHData, name::String)::Dict{Symbol,Any}

Retrieve the external dictionary for callback `name`.
"""
function get_callback_ext(phd::PHData, name::String)::Dict{Symbol,Any}
    cb = get_callback(phd, name)
    return cb.ext
end

"""
    is_leaf(phd::PHData, xhid::XhatID)::Bool

Returns true if the given consensus variable id belongs to a leaf vertex in the scenario tree.
"""
function is_leaf(phd::PHData, xhid::XhatID)::Bool
    return is_leaf(phd.scenario_tree, xhid.node)
end

"""
    leaf_value(phd::PHData, vid::VariableID)::Float64
    leaf_value(phd::PHData, scen::ScenarioID, stage::StageID, idx::Index)::Float64

Returns the value of the variable associated with `vid` or with scenario `scen`, stage `stage` and index `idx`. Must be a leaf variable.

**WARNING:** For computational efficiency, leaf values are collected only at the end of a PH run. Therefore, using this function in a callback will result in an error.

See also: [`branch_value`](@ref), [`value`](@ref)
"""
function leaf_value(phd::PHData, vid::VariableID)::Float64
    return phd.scenario_map[scenario(vid)].leaf_vars[vid]
end

function leaf_value(phd::PHData, scen::ScenarioID, stage::StageID, idx::Index)::Float64
    return leaf_value(phd, VariableID(scen, stage, idx))
end

"""
    name(phd::PHData, vid::VariableID)::String

Returns the name of the consensus variable for the given `VariableID`.
"""
function name(phd::PHData, vid::VariableID)::String
    if !haskey(phd.variable_data, vid)
        error("No name available for variable id $vid")
    end
    return phd.variable_data[vid].name
end

"""
    name(phd::PHData, xid::XhatID)::String

Returns the name of the consensus variable for the given [`XhatID`](@ref). The name is the same given to the individual scenario variables.
"""
function name(phd::PHData, xid::XhatID)::String
    return name(phd, first(convert_to_variable_ids(phd, xid)))
end

"""
    probability(phd::PHData, scenario::ScenarioID)::Float64

Returns the probability of the given scenario.
"""
function probability(phd::PHData, scenario::ScenarioID)::Float64
    return phd.scenario_map[scenario].prob
end

"""
    scenario_bundle(phd::PHData, xid::XhatID)::Set{ScenarioID}

Returns the scenarios contributing to the consensus variable associated with `xid`.
"""
function scenario_bundle(phd::PHData, xid::XhatID)::Set{ScenarioID}
    return scenario_bundle(phd.scenario_tree, xid.node)
end

"""
    scenario_tree(phd::PHData)::ScenarioTree

Returns the scenario tree for created stochastic programming problem.
"""
function scenario_tree(phd::PHData)::ScenarioTree
    return phd.scenario_tree
end

"""
    scenarios(phd::PHData)::Set{ScenarioID}

Returns the set of all scenarios for the stochastic problem.
"""
function scenarios(phd::PHData)::Set{ScenarioID}
    return scenarios(phd.scenario_tree)
end

"""
    stage_id(phd::PHData, xid::XhatID)::StageID

Returns the [`StageID`](@ref) in which the given consensus variable is.
"""
function stage_id(phd::PHData, xid::XhatID)::StageID
    return phd.scenario_tree.tree_map[xid.node].stage
end

"""
    print_timing([io::IO], phd::PHData)

Prints timining information from the Progressive Hedging solve.
"""
function print_timing(phd::PHData)
    println(phd.time_info)
    return
end

function print_timing(io::IO, phd::PHData)
    println(io, phd.time_info)
    return
end

"""
   value(phd::PHData, vid::VariableID)

Returns the value of the variable associated with `vid`.

See also: [`branch_value`](@ref), [`leaf_value`](@ref)
"""
function value(phd::PHData, vid::VariableID)::Float64
    return retrieve_variable_value(phd.scenario_map[scenario(vid)], vid)
end

"""
   value(phd::PHData, scen::ScenarioID, stage::StageID, idx::Index)::Float64

Returns the value of the variable associated with scenario `scen`, stage `stage` and index `idx`.

See also: [`branch_value`](@ref), [`leaf_value`](@ref)
"""
function value(phd::PHData, scen::ScenarioID, stage::StageID, idx::Index)::Float64
    vid = VariableID(scen, stage, idx)
    return value(phd, vid)
end

"""
   w_value(phd::PHData, vid::VariableID)

Returns the value of the variable associated with `vid`. Only available for branch variables.
"""
function w_value(phd::PHData, vid::VariableID)::Float64
    return phd.scenario_map[scenario(vid)].w_vars[vid]
end

"""
    w_value(phd::PHData, scen::ScenarioID, stage::StageID, idx::Index)::Float64

Returns the value of the variable associated with scenario `scen`, stage `stage` and index `idx`. Only available for branch variables.
"""
function w_value(phd::PHData, scen::ScenarioID, stage::StageID, idx::Index)::Float64
    return w_value(phd, VariableID(scen, stage, idx))
end

"""
   xhat_value(phd::PHData, xhid::VariableID)::Float64

Returns the value of the consensus variable associated with `xhid`. Only available for leaf variables after calling `solve`.  Available for branch variables at any time.
"""
function xhat_value(phd::PHData, xhat_id::XhatID)::Float64
    return value(phd.xhat[xhat_id])
end

"""
   xhat_value(phd::PHData, vid::VariableID)::Float64

Returns the value of the consensus variable associated with `vid`. Only available for leaf variables after calling `solve`.  Available for branch variables at any time.
"""
function xhat_value(phd::PHData, vid::VariableID)::Float64
    return xhat_value(phd, convert_to_xhat_id(phd, vid))
end

"""
   xhat_value(phd::PHData, scen::ScenarioID, stage::StageID, idx::Index)::Float64

Returns the value of the consensus variable associated with scenario `scen`, stage `stage` and index `idx`. Only available for leaf variables after calling `solve`.  Available for branch variables at any time.
"""
function xhat_value(phd::PHData, scen::ScenarioID, stage::StageID, idx::Index)::Float64
    return xhat_value(phd, VariableID(scen, stage, idx))
end


text:
# Model
include_model("rocket")
model_sl = free_time_model(additive_noise_model(model_slosh))

# Horizon
T = 41

# Bounds
ul, uu = control_bounds(model_sl, T, [-5.0; 0.0; 0.01], [5.0; 100.0; 1.0])

# Initial and final states
x1_slosh = [15.0;
	  model_sl.l1 + 10.0;
	  -45.0 * pi / 180.0;
	  0.0;
	  -10.0;
	  -10.0;
	  -1.0 * pi / 180.0;
	  0.0]

xT_slosh = [0.0;
	  model_sl.l1;
	  0.0;
	  0.0;
	  0.0;
	  0.0;
	  0.0;
	  0.0]

# Radius of landing pad
r_pad = 0.25

# xl <= x <= xl
_xl_slosh = -Inf * ones(model_sl.n)
_xl_slosh[2] = model_sl.l1
_xu_slosh = Inf * ones(model_sl.n)

xl_slosh, xu_slosh = state_bounds(model_sl, T, _xl_slosh, _xu_slosh,
 	x1 = x1_slosh, xT = xT_slosh)

xl_slosh[T][1] = -1.0 * r_pad
xu_slosh[T][1] = r_pad
xl_slosh[T][2] = xT_slosh[2] - 0.001
xu_slosh[T][2] = xT_slosh[2] + 0.001
xl_slosh[T][3] = -1.0 * pi / 180.0
xu_slosh[T][3] = 1.0 * pi / 180.0
xl_slosh[T][4] = -0.001
xu_slosh[T][4] = 0.001
xl_slosh[T][5] = -0.001
xu_slosh[T][5] = 0.001
xl_slosh[T][6] = -0.01 * pi / 180.0
xu_slosh[T][6] = 0.01 * pi / 180.0

# Objective
obj_slosh = quadratic_time_tracking_objective(
	[(t != T ? Diagonal([1.0; 10.0; 1.0; 0.1; 1.0; 10.0; 1.0; 0.1])
		: Diagonal([10.0; 100.0; 10.0; 0.1; 10.0; 100.0; 10.0; 0.1])) for t = 1:T],
	[Diagonal([1.0e-1 * ones(2); 0.0]) for t = 1:T-1],
    [xT_slosh for t = 1:T],
	[zeros(model_sl.m) for t = 1:T],
	1.0)

# Time step constraints
con_free_time = free_time_constraints(T)

# Problem
prob_slosh = trajectory_optimization_problem(model_sl,
					obj_slosh,
					T,
                    xl = xl_slosh,
                    xu = xu_slosh,
                    ul = ul,
                    uu = uu,
                    con = con_free_time)

# Trajectory initialization
x0_slosh = linear_interpolation(x1_slosh, xT_slosh, T) # linear interpolation on state
u0 = [ones(model_sl.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0_slosh = pack(x0_slosh, u0, prob_slosh)

# Solve
if false # set to true to reoptimize
	include_snopt()
    @time z̄_slosh , info = solve(prob_slosh, copy(z0_slosh),
		nlp = :SNOPT7,
		time_limit = 60 * 10)
    @save joinpath(@__DIR__, "sol_slosh.jld2") z̄_slosh
else
    println("Loading solution...")
    @load joinpath(@__DIR__, "sol_slosh.jld2") z̄_slosh
end

x̄_slosh, ū_slosh = unpack(z̄_slosh, prob_slosh)
@show sum([ū_slosh[t][end] for t = 1:T-1]) # ~2.74

ext: jl


text:
export nonspatial, nonspatial!

function birth!(tumor::Vector{Cell}, mutations::Vector{Mutation}, parent, cur_id, cur_mutation, μ, t)

    pos = Float64[]

	new = Cell(cur_id, pos, copy(parent.mutations), parent.index, parent.b, t, SVector{0,Float64}())
	push!(tumor, new)

	m = 0
	if rand() < μ/2
		m += 1
		push!(mutations, Mutation(parent.index, isempty(parent.mutations) ? 0 : last(parent.mutations), t, length(tumor), SVector{0,Float64}()))
		push!(parent.mutations, cur_mutation+m)
	end
	if rand() < μ/2
		m += 1
		push!(mutations, Mutation(new.index, isempty(new.mutations) ? 0 : last(new.mutations), t, length(tumor), SVector{0,Float64}()))
		push!(new.mutations, cur_mutation+m)
	end
	return m
end

function nonspatial!( tumor::Vector{Cell}, mutations::Vector{Mutation}, until; b, d, μ, cur_id = 1, cur_mutation = 0, t = 0.0, seed=nothing, showprogress=true)

    isnothing(seed) || Random.seed!(seed)

    N = length(tumor)

    prog = showprogress ? ProgressUnknown("Progress: Tumor size ") : nothing
    while loop_condition(N,t, until)

        N==0 && error("Tumor died")
        row = rand(1:N)
        parent = tumor[row]

        p = rand()*(b+d) - b
        t += randexp()/((d+b)*N)

        if p < 0.
            cur_id += 1
            N += 1
            cur_mutation += birth!(tumor, mutations, parent, cur_id, cur_mutation, μ, t)
        elseif p < d
            N -= 1
            deleteat!(tumor, row)
        end
        showprogress && ProgressMeter.update!(prog, N )
    end
    return Dict{Symbol, Any}(:index => cur_id, :mutation => cur_mutation, :time => t)
end
nonspatial!( until; tumor::Vector{Cell}, mutations::Vector{Mutation}, simparams...) =  nonspatial!( tumor, mutations, until; simparams...)

function nonspatial( until; cur_mutation = 0, b, simparams...)
    p₀ = Float64[]
    tumor = [Cell(index=1, position=p₀, parent=0, mutations=collect(1:cur_mutation), b=b, t_birth=0.0, p_birth=SVector{0,Float64}(p₀))]

	mutations = [Mutation(; origin=1, ancestor=0, t_birth=0., N_birth=1, p_birth=SVector{0,Float64}(p₀)) for m=1:cur_mutation]

    output = nonspatial!(tumor, mutations, until; cur_mutation=cur_mutation, b=b, simparams...)
	output[:tumor] = tumor
	output[:mutations] = mutations

	return output
end

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
## connect to Optim.jl as backend
function sem_wrap_optim(par, F, G, H, sem::AbstractSem)
    sem(par, !isnothing(F), !isnothing(G), !isnothing(H))
    if !isnothing(G) G .= gradient(sem) end
    if !isnothing(H) H .= hessian(sem) end
    if !isnothing(F) return objective(sem)[1] end
end

function SemFit(optimization_result::Optim.MultivariateOptimizationResults, model::AbstractSem, start_val)
    return SemFit(
        optimization_result.minimum,
        optimization_result.minimizer,
        start_val,
        model,
        optimization_result
    )
end

optimizer(res::Optim.MultivariateOptimizationResults) = Optim.summary(res)
n_iterations(res::Optim.MultivariateOptimizationResults) = Optim.iterations(res)
convergence(res::Optim.MultivariateOptimizationResults) = Optim.converged(res)

function sem_fit(model::Sem{O, I, L, D}; start_val = start_val, kwargs...) where {O, I, L, D <: SemDiffOptim}

    if !isa(start_val, Vector)
        start_val = start_val(model; kwargs...)
    end

    result = Optim.optimize(
                Optim.only_fgh!((F, G, H, par) -> sem_wrap_optim(par, F, G, H, model)),
                start_val,
                model.diff.algorithm,
                model.diff.options)
    return SemFit(result, model, start_val)

end

function sem_fit(model::SemFiniteDiff{O, I, L, D}; start_val = start_val, kwargs...) where {O, I, L, D <: SemDiffOptim}

    if !isa(start_val, Vector)
        start_val = start_val(model; kwargs...)
    end

    result = Optim.optimize(
                Optim.only_fgh!((F, G, H, par) -> sem_wrap_optim(par, F, G, H, model)),
                start_val,
                model.diff.algorithm,
                model.diff.options)
    return SemFit(result, model, start_val)

end

function sem_fit(model::SemEnsemble{N, T , V, D, S}; start_val = start_val, kwargs...) where {N, T, V, D <: SemDiffOptim, S}

    if !isa(start_val, Vector)
        start_val = start_val(model; kwargs...)
    end

    result = Optim.optimize(
                Optim.only_fgh!((F, G, H, par) -> sem_wrap_optim(par, F, G, H, model)),
                start_val,
                model.diff.algorithm,
                model.diff.options)
    return SemFit(result, model, start_val)

end


text:
import ACEbase: acquire!, release!
using Base.Threads: threadid, nthreads
using DataStructures: Stack

struct ArrayCache{T}
   vecs::Vector{Stack{Vector{T}}}
   mats::Vector{Stack{Matrix{T}}}
end

struct CachedArray{N, T} <: AbstractArray{T, N}
   A::Array{T, N}
   pool::ArrayCache{T}
end

release!(A::Any) = nothing
release!(pA::CachedArray) = release!(pA.pool, pA)

using Base: @propagate_inbounds

@propagate_inbounds function Base.getindex(pA::CachedArray, I...)
   @boundscheck checkbounds(pA.A, I...)
   @inbounds pA.A[I...]
end

@propagate_inbounds function Base.setindex!(pA::CachedArray, val, I...)
   @boundscheck checkbounds(pA.A, I...)
   @inbounds pA.A[I...] = val
end

# Base.getindex(pA::CachedArray, args...) = getindex(pA.A, args...)

# Base.setindex!(pA::CachedArray, args...) = setindex!(pA.A, args...)

Base.length(pA::CachedArray) = length(pA.A)

Base.eltype(pA::CachedArray) = eltype(pA.A)

Base.size(pA::CachedArray, args...) = size(pA.A, args...)

Base.parent(pA::CachedArray) = pA.A

function ArrayCache{T}() where {T}
   nt = nthreads()
   vecs = [ Stack{Vector{T}}() for _=1:nt ]
   mats = [ Stack{Matrix{T}}() for _=1:nt ]
   return ArrayCache(vecs, mats)
end

acquire!(c::ArrayCache{T}, ::Type{T}, len::Integer) where {T} =
         acquire!(c, len)

acquire!(c::ArrayCache{T}, len::Integer, ::Type{T}) where {T} =
         acquire!(c, len)

acquire!(c::ArrayCache{T}, ::Type{S}, len::Integer) where {T, S} =
         Vector{S}(undef, len)

acquire!(c::ArrayCache{T}, len::Integer, ::Type{S}) where {T, S} =
         Vector{S}(undef, len)

function acquire!(c::ArrayCache{T}, len::Integer) where {T}
   stack = c.vecs[threadid()]
   if isempty(stack)
      A = Vector{T}(undef, len)
   else
      A = pop!(stack)
      resize!(A, len)
   end
   return CachedArray(A, c)
end

release!(c::ArrayCache, cA::CachedArray{1}) =
      push!(c.vecs[threadid()], cA.A)

## Alternative Array Cache Implementation

struct GenArrayCache
   vecs::Vector{Stack{Vector{UInt8}}}
end

# struct GenCachedArray{T, N, AT <: AbstractArray{T, N}} <: AbstractArray{T, N}
#    A::AT
#    pool::GenArrayCache
# end

using UnsafeArrays

struct GenCachedArray{T, N} <: AbstractArray{T, N}
   A::UnsafeArray{T, N}
   _A::Vector{UInt8}
   pool::GenArrayCache
end

# release!(A::Any) = nothing
release!(pA::GenCachedArray) = release!(pA.pool, pA)

using Base: @propagate_inbounds

@propagate_inbounds function Base.getindex(pA::GenCachedArray, I...)
   @boundscheck checkbounds(pA.A, I...)
   @inbounds pA.A[I...]
end

@propagate_inbounds function Base.setindex!(pA::GenCachedArray, val, I...)
   @boundscheck checkbounds(pA.A, I...)
   @inbounds pA.A[I...] = val
end

# Base.getindex(pA::CachedArray, args...) = getindex(pA.A, args...)

# Base.setindex!(pA::CachedArray, args...) = setindex!(pA.A, args...)

Base.length(pA::GenCachedArray) = length(pA.A)

Base.eltype(pA::GenCachedArray) = eltype(pA.A)

Base.size(pA::GenCachedArray, args...) = size(pA.A, args...)

Base.parent(pA::GenCachedArray) = pA.A

function GenArrayCache() where {T}
   nt = nthreads()
   vecs = [ Stack{Vector{UInt8}}() for _=1:nt ]
   return GenArrayCache(vecs)
end

function acquire!(c::GenArrayCache, ::Type{T}, len::Integer) where {T}
   szofT = sizeof(T)
   szofA = len * szofT
   stack = c.vecs[threadid()]
   if isempty(stack)
      _A = Vector{UInt8}(undef, szofA)
   else
      _A = pop!(stack)
      resize!(_A, szofA)
   end
   # A = reinterpret(T, _A)
   # return GenCachedArray(A, c)

   ptr = Base.unsafe_convert(Ptr{T}, _A)
   # A = Base.unsafe_wrap(Array, ptr, len)
   A = UnsafeArray(ptr, (len,))
   return GenCachedArray(A, _A, c)
end

# release!(c::GenArrayCache, cA::GenCachedArray) =
#       push!(c.vecs[threadid()], parent(cA.A))

release!(c::GenArrayCache, cA::GenCachedArray) =
      push!(c.vecs[threadid()], cA._A)


text:
__precompile__()

"""
    MMD Kernels
"""
module Kernels

export expect, RBFKernel, loss

abstract type AbstractKernel end

struct RBFKernel <: AbstractKernel
    sigmas::Vector{Float64}
    matrix::Matrix{Float64}
end

function RBFKernel(nqubits::Int, sigmas::Vector{Float64}, isbinary::Bool)
    basis = collect(0:(1<<nqubits - 1))
    return RBFKernel(sigmas, rbf_kernel_matrix(basis, basis, sigmas, isbinary))
end

expect(kernel::RBFKernel, px::Vector{Float64}, py::Vector{Float64}) = px' * kernel.matrix * py

# RBF Kernel

function rbf_kernel_matrix(x::Vector, y::Vector, sigmas::Vector{Float64}, isbinary::Bool)
    if length(sigmas) == 0
        throw(DimensionMismatch("At least 1 sigma prameters is required for RBM kernel!"))
    end

    if isbinary
        dx2 = map(count_ones, xor.(x, y'))
    else
        dx2 = (x .- y').^2
    end

    K = 0
    for sigma in sigmas
        gamma = 1.0 / (2 * sigma)
        K = K + exp.(-gamma * dx2)
    end
    return K
end

function loss(px::AbstractVecOrMat{Float64}, kernel::RBFKernel, py::AbstractVecOrMat{Float64})
    pxy = px - py
    return expect(kernel, pxy, pxy)
end

end

max_stars_repo_name: Roger-luo/QCBM

max_issues_repo_name: Roger-luo/QCBM

max_forks_repo_name: Roger-luo/QCBM


text:
using JSONSchema, JSON
import BinaryProvider

using Test

### load the "json-schema-org/JSON-Schema-Test-Suite" project from github
tsurl = "https://github.com/json-schema-org/JSON-Schema-Test-Suite/archive/master.tar.gz"

destdir = mktempdir()
dwnldfn = joinpath(destdir, "test-suite.tar.gz")
BinaryProvider.download(tsurl, dwnldfn, verbose=true)

unzipdir = joinpath(destdir, "test-suite")
BinaryProvider.unpack(dwnldfn, unzipdir)

jsonTestFilesDirectory = joinpath(unzipdir, "JSON-Schema-Test-Suite-master", "tests")
localReferenceTestsDirectory = mktempdir(jsonTestFilesDirectory)

"""
Write test files for locally referenced schema files.
These files have the same format as JSON Schema org test files. They are written
to a sibling directory to JSON-Schema-Test-Suite-master/tests/draft* directories
so they can be consumed the same way as the draft*/*.json test files.
"""
function writeLocalReferenceTestFiles()
    # sibling directory for testing a relative path containing "../"
    referenceComponentsDirectory = mktempdir(jsonTestFilesDirectory)

    write(joinpath(referenceComponentsDirectory, "localReferenceSchemaOne.json"), """
    {
        "type": "object",
        "properties": {
            "localRefOneResult": { "type": "string" }
        }
    }
    """)

    write(joinpath(referenceComponentsDirectory, "localReferenceSchemaTwo.json"), """
    {
        "type": "object",
        "properties": {
            "localRefTwoResult": { "type": "number" }
        }
    }
    """)

    write(joinpath(referenceComponentsDirectory, "nestedLocalReference.json"), """
    {
        "type": "object",
        "properties": {
            "result": { "\$ref": "file:localReferenceSchemaOne.json#/properties/localRefOneResult" }
        }
    }
    """)

    write(joinpath(localReferenceTestsDirectory, "localReferenceTest.json"), """[
    {
        "description": "test locally referenced schemas",
        "schema": {
            "type": "object",
            "properties": {
                "result1": { "\$ref": "file:../$(basename(abspath(referenceComponentsDirectory)))/localReferenceSchemaOne.json#/properties/localRefOneResult" },
                "result2": { "\$ref": "../$(basename(abspath(referenceComponentsDirectory)))/localReferenceSchemaTwo.json#/properties/localRefTwoResult" }
            },
            "oneOf": [
                {
                    "required": ["result1"]
                },
                {
                    "required": ["result2"]
                }
            ]
        },
        "tests": [
            {
                "description": "reference only local schema 1",
                "data": {"result1": "some text" },
                "valid": true
            },
            {
                "description": "reference only local schema 2",
                    "data": {"result2": 1234 },
                    "valid": true
            },
            {
                "description": "incorrect reference to local schema 1",
                    "data": { "result1": true },
                    "valid": false
            },
            {
                "description": "reference neither local schemas",
                    "data": { "result": true },
                    "valid": false
            },
            {
                "description": "reference both local schemas",
                "data": {"result1": "some text", "result2": 500 },
                "valid": false
            }
        ]
    }
    ]""")

    write(joinpath(localReferenceTestsDirectory, "nestedLocalReferenceTest.json"), """[
    {
        "description": "test locally referenced schemas",
        "schema": {
            "type": "object",
            "properties": {
                "result": { "\$ref": "file:../$(basename(abspath(referenceComponentsDirectory)))/nestedLocalReference.json#/properties/result" }
            }
        },
        "tests": [
            {
                "description": "nested reference, correct type",
                "data": {"result": "some text" },
                "valid": true
            },
            {
                "description": "nested reference, incorrect type",
                "data": {"result": 1234 },
                "valid": false
            }
        ]
    }
    ]""")
    # return directory name (not path) of tests for use in testing below
    return
end

writeLocalReferenceTestFiles()

@testset begin

    ################################################################################
    ### Applying test suites for draft 4/6 specifications, and local ref tests   ###
    ################################################################################

    # add custom directory containing tests for locally referenced schema files
    localRefTestDirectoryName = basename(abspath(localReferenceTestsDirectory))
    @testset "Test suite for $draftfn" for draftfn in ["draft4", "draft6", localRefTestDirectoryName]
        tsdir = joinpath(jsonTestFilesDirectory, draftfn)

        # the test suites use the 'remotes' folder to simulate remote refs with the
        #  'http://localhost:1234' url.  To have tests cope with this, the id dictionary
        # is preloaded with the files in ''../remotes'
        idmap0 = Dict{String, Any}()
        remfn = joinpath(tsdir, "../../remotes")
        for rn in ["integer.json", "name.json", "subSchemas.json", "folder/folderInteger.json"]
            idmap0["http://localhost:1234/" * rn] = Schema(JSON.parsefile(joinpath(remfn, rn))).data
        end

        @testset "$tfn" for tfn in filter(n -> occursin(r"\.json$",n), readdir(tsdir))
            fn = joinpath(tsdir, tfn)
            schema = JSON.parsefile(fn)
            @testset "- $(subschema["description"])" for subschema in (schema)
                spec = subschema["schema"] isa Bool ?
                    Schema(subschema["schema"]; idmap0=idmap0) :
                    Schema(subschema["schema"]; idmap0=idmap0, parentFileDirectory = dirname(fn))

                @testset "* $(subtest["description"])" for subtest in subschema["tests"]
                    @test isvalid(subtest["data"], spec) == subtest["valid"]
                end
            end
        end
    end

end


text:
using Metal
using Test
using Random
using LinearAlgebra

# All M-series chips have a unified memory model (as of this commit).
# This means the CPU, GPU, and other hardware accelerators share the
# exact same memory pool. This is unlike most other GPU systems
# where the CPU and GPU have different physical memory that must
# be synced (usually across a PCI-E) connection.
# A unified memory system is beneficial for many workloads because
# these transfers are unnecessary and allow for a more seamless
# use of the available heterogeneous compute architecture.

# This document is meant to showcase potential usecases allowed by unified memory.

function simple_kernel(arr)
    idx = thread_position_in_grid_1d()
    arr[idx] = cos(arr[idx])
    return
end

### Unsafe wrapping a GPU array by a CPU array and altering eachother's (same) data
# Note that since the GPU array allocation requires more assumptions about the
# underlying data, in a unified memory architecture, the GPU array should first
# be allocated, then wrapped by a CPU array...not the other way around.

dims = tuple(16,16)
# Create a Metal array with a default storage mode of shared (both CPU and GPU get access)
arr_mtl = MtlArray{Float32}(undef, dims)
# Unsafe wrap the contents of the Metal array with a CPU array
arr_cpu = unsafe_wrap(Array{Float32}, arr_mtl, dims)

# Alter the CPU array data
arr_cpu .= pi

# These changes are reflected in the Metal array
Metal.@allowscalar @test arr_mtl[1] == Float32(pi)

# Now launch a kernel altering the Metal array
Metal.@sync @metal threads=1024 grid=1024 simple_kernel(arr_mtl)

# These changes are reflected in the wrapped CPU array
@test all(arr_cpu .== -1.0)

### Using CPU functions that are unavailable to the GPU in the middle of a series of kernels
# This would otherwise require a set of copies back and forth.

round.(rand!(arr_cpu)*100)
# Example 1: Calculate the deteminant using CPU implementation
det(arr_cpu)
# Example 2: Singular Value Decomposition
svd(arr_cpu)

# Care must still be taken to make sure that any GPU work is finished before altering data on the CPU
# GPU kernel launches are still ASYNCHRONOUS!
# Relevant doc page: https://developer.apple.com/documentation/metal/mtlstoragemode/mtlstoragemodeshared?language=objc
# TODO: Come up with simultaneous launch of GPU and CPU work that results in undesired behavior

function long_kernel(arr, dummy)
    idx = thread_position_in_grid_1d()
    for i in 1:100000
        dummy[1] += Float32(0.3)
    end
    arr[idx] = Float32(pi)
    arr[idx] = cos(arr[idx])
    return
end

# Make larger arrays to make the kernel take non-trivial time
dims = 1024*1024
# Create a Metal array with a default storage mode of shared (both CPU and GPU get access)
arr_mtl = MtlArray{Float32}(undef, dims)
# Unsafe wrap the contents of the Metal array with a CPU array
arr_cpu = unsafe_wrap(Array{Float32}, arr_mtl, dims)
dummy_mtl = MtlArray{Float32}(undef, 1)

rand!(arr_cpu)
# Now launch a kernel altering the Metal array
@metal threads=1024 grid=1024 long_kernel(arr_mtl, dummy_mtl)

# we need to synchronize the device as the kernel may not have finished yet
synchronize()
@test all(arr_cpu .== -1.0)


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

module MicroPerf
import Perftests: @perf, meta

## recursive fib ##
fib(n) = n < 2 ? n : fib(n-1) + fib(n-2)
@perf fib(20) meta("fib", "Recursive fibonacci")

## parse integer ##
function parseintperf(t)
    local n, m
    for i=1:t
        n = rand(UInt32)
        s = hex(n)
        m = UInt32(parse(Int64,s,16))
    end
    return n
end

@perf parseintperf(1000) meta("parse_int", "Integer parsing")

## mandelbrot set: complex arithmetic and comprehensions ##

function mandel(z)
    c = z
    maxiter = 80
    for n = 1:maxiter
        if abs(z) > 2
            return n-1
        end
        z = z^2 + c
    end
    return maxiter
end

mandelperf() = [mandel(complex(r,i)) for i=-1.:.1:1., r=-2.0:.1:0.5]
@perf mandelperf() meta("mandel", "Calculation of mandelbrot set")

## numeric vector sort ##
function qsort!(a,lo,hi)
    i, j = lo, hi
    while i < hi
        pivot = a[(lo+hi)>>>1]
        while i <= j
            while a[i] < pivot; i += 1; end
            while a[j] > pivot; j -= 1; end
            if i <= j
                a[i], a[j] = a[j], a[i]
                i, j = i+1, j-1
            end
        end
        if lo < j; qsort!(a,lo,j); end
        lo, j = i, hi
    end
    return a
end

sortperf(n) = qsort!(rand(n), 1, n)
@perf sortperf(5000) meta("quicksort", "Sorting of random numbers using quicksort")

## slow pi series ##
function pisum()
    sum = 0.0
    for j = 1:500
        sum = 0.0
        for k = 1:10000
            sum += 1.0/(k*k)
        end
    end
    sum
end
@perf pisum() meta("pi_sum", "devectorized", "Devectorized summation of a power series")

## slow pi series, vectorized ##
function pisumvec()
    s = 0.0
    a = [1:10000...]
    for j = 1:500
        s = sum(1./(a.^2))
    end
    s
end

@perf pisumvec() meta("pi_sum", "vectorized", "Vectorized summation of a power series")

## random matrix statistics ##
function randmatstat(t)
    n = 5
    v = zeros(t)
    w = zeros(t)
    for i=1:t
        a = randn(n,n)
        b = randn(n,n)
        c = randn(n,n)
        d = randn(n,n)
        P = [a b c d]
        Q = [a b; c d]
        v[i] = trace((P.'*P)^4)
        w[i] = trace((Q.'*Q)^4)
    end
    return (std(v)/mean(v), std(w)/mean(w))
end

@perf randmatstat(1000) meta("rand_mat_stat", "Statistics on a random matrix")

## largish random number gen & matmul ##

@perf rand(1000,1000)*rand(1000,1000) meta("rand_mat_mul", "Multiplication of random matrices")

## printfd ##
@unix_only begin
    function printfd(n)
        open("/dev/null","w") do io
            for i = 1:n
                @printf(io,"%d %d\n",i,i+1)
            end
        end
    end

    printfd(1)
    @perf printfd(100000) meta("printfd", "Printing to a file descriptor")
end

end # module


text:
using EchelleCCFs
using Test

    @testset "CCF plan" begin
        Δv = 7e3
        λ1 = 5025.0
        d1 = 0.5
        d2 = 0.25
        λ2 = 5075.0
        @test_nowarn BasicLineList([λ1,λ2], [d1,d2])
        ll = BasicLineList([λ1,λ2], [d1,d2])
        @test_nowarn BasicCCFPlan(line_list=ll, mask_shape=TopHatCCFMask(Δv))
        ccfpl = BasicCCFPlan(line_list=ll, mask_shape=TopHatCCFMask(Δv))
        @test_nowarn calc_ccf_v_grid(ccfpl)
        v_grid = calc_ccf_v_grid(ccfpl)
        @test length(v_grid) == calc_length_ccf_v_grid(ccfpl)
    end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

include("strings/substring.jl")
include("strings/search.jl")
include("strings/unicode.jl")

import .Unicode: textwidth, islowercase, isuppercase, isalpha, isdigit, isnumeric, iscntrl, ispunct,
    isspace, isprint, isxdigit, lowercase, uppercase, titlecase, lowercasefirst, uppercasefirst

include("strings/util.jl")
include("strings/io.jl")

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
using PyPlot, MRIReco

##### simple example ####

N = 256
numCoils = 8
I = shepp_logan(N)
I = circularShutterFreq!(I,1)

coilsens = birdcageSensitivity(N, 8, 1.5)

# simulation parameters
params = Dict{Symbol, Any}()
params[:simulation] = "fast"
params[:trajName] = "Spiral"
params[:numProfiles] = 6
params[:numSamplingPerProfile] = div(N*N,16)
params[:windings] = div(N,16)
params[:AQ] = 2.0e-2
params[:senseMaps] = coilsens

# do simulation
acqData = simulation(I, params)

# reco parameters
params = Dict{Symbol, Any}()
params[:reco] = "multiCoil"
params[:reconSize] = (N,N)
params[:regularization] = "L2"
params[:λ] = 1.e-3
params[:iterations] = 40
params[:solver] = "cgnr"
params[:senseMaps] = coilsens
Ireco = reconstruction(acqData, params)

# use PyPlot for interactive display
figure(1)
clf()
subplot(1,2,1)
title("Phantom")
imshow(abs.(I))
subplot(1,2,2)
title("Reconstruction")
imshow(abs.(Ireco[:,:,1,1,1]))

# export images
filename = joinpath(dirname(pathof(MRIReco)),"../docs/src/assets/senseReco.png")
exportImage(filename, abs.(Ireco[:,:,1,1,1]) )


text:
### Core functionality for method deletion
using Core: MethodInstance
using Base: MethodList

const themethod = Ref{Any}(nothing)

function delete_method(m::Method)
    # While methods(getfield(m.module, m.name)) works in many cases, it fails when
    # a nonex-ported function is extended by a different module.
    # Invalidate the cache of the callers
    try
    invalidate_callers(m.specializations)
    # Remove from method list
    mt = get_methodtable(m.sig)
    ml = MethodList(mt)
    drop_from_tme!(mt, :defs, m)
    # Delete compiled instances of this method from the cache
    drop_from_tme!(mt, :cache, m)
    # Delete this signature
    deleteat!(ml.ms, findfirst(x->x==m, ml.ms))
    catch err
        themethod[] = m
        showerror(STDERR, err)
        rethrow(err)
    end
nothing
end

delete_method(::Void) = nothing

get_methodtable(u::UnionAll) = get_methodtable(u.body)
get_methodtable(sig) = _get_methodtable(sig.parameters[1])
_get_methodtable(u::UnionAll) = (@show u; _get_methodtable(u.body))
_get_methodtable(f) = f.name.mt

invalidate_callers(::Void) = nothing

function invalidate_callers(tml::TypeMapLevel)
    global themethod
    println("here")
    themethod[] = tml
    nothing
end

function invalidate_callers(tme::TypeMapEntry)
    while isa(tme, TypeMapEntry)
        invalidate_callers(tme.func)
        tme = tme.next
    end
    nothing
end

function invalidate_callers(mi::MethodInstance, iddict=ObjectIdDict())
    iddict[mi] = true
    if isdefined(mi, :backedges)
        for c in mi.backedges
            haskey(iddict, c) && continue
            invalidate_callers(c, iddict)
        end
    end
    mtc = get_methodtable(mi.def.sig)
    drop_from_tme!(mtc, :cache, mi)
    drop_from_tme!(mi.def, :specializations, mi)
end

function drop_from_tme!(mt::Union{MethodTable,Method}, fn::Symbol, m::Union{Method,MethodInstance})
    nodeprev, node = nothing, getfield(mt, fn)
    while isa(node, TypeMapEntry)
        if node.func === m
            if nodeprev == nothing
                setfield!(mt, fn, node.next)
            else
                nodeprev.next = node.next
            end
        end
        nodeprev, node = node, node.next
    end
    mt
end

### Parsing expressions to determine which method to delete
const ExLike = Union{Expr,Revise.RelocatableExpr}
# Much is taken from ExpressionUtils.jl but generalized to work with ExLike

function get_signature(ex::E) where E <: ExLike
    while ex.head == :macrocall && isa(ex.args[end], E) || is_trivial_block_wrapper(ex)
        ex = ex.args[end]::E
    end
    if ex.head == :function
        return ex.args[1]
    elseif ex.head == :(=) && isa(ex.args[1], E)
        ex = ex.args[1]::E
        if ex.head == :where || ex.head == :call
            return ex
        end
    end
    nothing
end

function get_method(mod::Module, sig::ExLike)
    t = split_sig(mod, convert(Expr, sig))
    mths = Base._methods_by_ftype(t, -1, typemax(UInt))
    if !isempty(mths)
        # There might be many methods, but the one that should match should be the
        # last one, since methods are ordered by specificity
        return mths[end][3]
    end
    warn("Revise failed to find any methods for signature ", t, "\n  Most likely it was already deleted.")
    nothing
end

function split_sig(mod::Module, ex::Expr)
    t = split_sig_expr(mod, ex)
    return eval(mod, t) # fex), eval(mod, argex)
end

function split_sig_expr(mod::Module, ex::Expr, wheres...)
    if ex.head == :where
        return split_sig_expr(mod, ex.args[1], ex.args[2:end], wheres...)
    end
    fex = ex.args[1]
    sigex = Expr(:curly, :Tuple, :(typeof($fex)), argtypeexpr.(ex.args[2:end])...)
    for w in wheres
        sigex = Expr(:where, sigex, w...)
    end
    sigex
end

function is_trivial_block_wrapper(ex::ExLike)
    if ex.head == :block
        return length(ex.args) == 1 ||
            (length(ex.args) == 2 && (is_linenumber(ex.args[1]) || ex.args[1]===nothing))
    end
    false
end
is_trivial_block_wrapper(Compat.@nospecialize arg) = false

function is_linenumber(Compat.@nospecialize stmt)
    (isa(stmt, ExLike) && (stmt).head == :line) || isa(stmt, LineNumberNode)
end

argtypeexpr(s::Symbol) = :Any
function argtypeexpr(ex::ExLike)
    if ex.head == :...
        return :(Vararg{$(argtypeexpr(ex.args[1]))})
    end
    if ex.head != :(::)
        @show ex.head ex.args
    end
    ex.head == :(::) || throw(ArgumentError("expected :(::) expression, got ex.head = $(ex.head)"))
    1 <= length(ex.args) <= 2 || throw(ArgumentError("expected 1 or 2 args, got $(ex.args)"))
    ex.args[end]
end


text:
using DocOpt

doc="""WORLD-based mel-cepstrum extraction for audio signals.
Usage:
    mcep.jl [options] <src_dir> <dst_dir>
    mcep.jl --version
    mcep.jl -h | --help

Options:
    -h --help        show this message
    --period=PERIOD  frame period in msec [default: 5.0]
    --order=ORDER    order of mel cepsrum [default: 40]
    --alpha=ALPHA    all-pass constant [default: 0.0]
    --max=MAX        Maximum number that will be processed [default: 200]
"""

using VoiceConversion
using WAV
using WORLD
using HDF5, JLD

import MelGeneralizedCepstrums: mcepalpha

using Logging
@Logging.configure(level=DEBUG, output=STDOUT)

searchdir(path, key) = filter(x -> contains(x, key), readdir(path))

function mkdir_if_not_exist(dir)
    if !isdir(dir)
        println("Create $(dir)")
        run(`mkdir -p $dir`)
    end
end

# process one file
function _mcep(path, period, order, α, savepath)
    x, fs = wavread(path)
    @assert size(x, 2) == 1 "The input data must be monoral."
    x = map(Float64, vec(x))
    fs = Int(fs)

    if α == 0.0
        α = mcepalpha(fs)
    end

    f0, timeaxis = dio(x, fs, DioOption(period=period))
    f0 = stonemask(x, fs, timeaxis, f0)
    spectrogram = cheaptrick(x, fs, timeaxis, f0)
    mc = sp2mc(spectrogram, order, α)

    save(savepath,
         "description", "WORLD-based Mel-cepstrum",
         "type", "MelCepstrum",
         "fs", fs,
         "period", period,
         "order", order,
         "fftlen", get_fftsize_for_cheaptrick(fs),
         "alpha", α,
         "feature_matrix", mc
         )
end

let
    args = docopt(doc, version=v"0.0.2")

    srcdir = args["<src_dir>"]
    dstdir = args["<dst_dir>"]
    mkdir_if_not_exist(dstdir)

    period = parse(Float64, args["--period"])
    order = parse(Int, args["--order"])
    nmax = parse(Int, args["--max"])
    α = parse(Float64, args["--alpha"])

    files = searchdir(srcdir, ".wav")
    @info("$(length(files)) data found.")

    # perform feature extraction for each file
    for (count, filename) in enumerate(files)
        path = joinpath(srcdir, filename)
        savepath = joinpath(dstdir, string(splitext(basename(path))[1], "_wmcep.jld"))

        @info("Start processing $(path)")
        elapsed = @elapsed _mcep(path, period, order, α, savepath)
        @info("Elapsed time in feature extraction is $(elapsed) sec.")
        @info("Dumped to $(savepath)")

        count += 1
        count > nmax && break
    end

    println("Finished")
end


text:
sa2y(dims::Tuple,is::Int,ia::Int)=sub2ind(dims,is,ia)
y2sa(dims::Tuple,iy::Int)=ind2sub(dims,iy)

function sa2y(dims::Tuple,sa::Array{Int,2})
	T=size(sa,1)
	y=zeros(Int,T)
	for t=1:T
		y[t]=sa2y(dims,sa[t,1],sa[t,2])
	end
	y
end

sa2y(dims::Tuple,sa::Array{Array{Int,2},1})=map(sai->sa2y(dims,sai),sa)

function y2sa(dims::Tuple,y::Array{Int,1})
	T=length(y)
	sa=zeros(Int,T,2)
	for t=1:T
		sai=y2sa(dims,y[t])
		sa[t,1]=sai[1]
		sa[t,2]=sai[2]
	end
	sa
end

y2sa(dims::Tuple,y::Array{Array{Int,1},1})=map(yi->y2sa(dims,yi),y)

k2xs=y2sa
xs2k=sa2y


text:
#!/usr/bin/env julia

using Luxor, Colors
w, h = 600, 600

fname = "/tmp/stars.png"
Drawing(w, h, fname)
origin()
cols = [RGB(rand(3)...) for i in 1:50]
background("grey20")
x = -w/2
for y in 100 * randn(h, 1)
    setcolor(cols[rand(1:end)])
    # star(x::Real, y::Real, radius::Real, npoints::Int64=5, ratio::Real=0.5, orientation=0, action=:nothing; vertices = false, reversepath=false)
    star(x, y, 10, rand(4:7), rand(3:7)/10, 0, :fill)
    x += 2
end
finish()
println("finished test: output in $(fname)")


text:
using Random, Distributions, StatsPlots; pyplot()
Random.seed!(2)

alpha  = 0.05
L(obs) = obs - (1-sqrt(alpha))
U(obs) = obs + (1-sqrt(alpha))

mu = 5.57
triDist = TriangularDist(mu-1,mu+1,mu)

N = 100
hitBounds, missBounds = zeros(N, 2), zeros(N,2)
for i in 1:N
    observation = rand(triDist)
    LL, UU = L(observation), U(observation)
    if LL <= mu && mu <= UU
        hitBounds[i,:] = [LL  UU-LL]
    else
        missBounds[i,:] = [LL  UU-LL]
    end
end

groupedbar(hitBounds, bar_position=:stack,
    c=:blue, la=0, fa=[0 1], label="", ylims=(3,8))
groupedbar!(missBounds, bar_position=:stack,
    c=:red, la=0, fa=[0 1], label="", ylims=(3,8))
plot!([0,N+1],[mu,mu],
    c=:black, xlims=(0,N+1),
    ylims=(3,8), label="Parameter value", ylabel="Value Estimate")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using LightQuery
import Documenter
using Documenter: deploydocs, makedocs
makedocs(
    sitename = "LightQuery.jl",
    modules = [LightQuery],
    doctest = false,
    pages = [
        "Usage and performance notes" => "index.md",
        "Beginner tutorial" => "beginner_tutorial.md",
        "Reshaping tutorial" => "reshaping_tutorial.md",
        "Interface" => "interface.md",
    ],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
)
deploydocs(repo = "github.com/bramtayl/LightQuery.jl.git")


text:
context("faces") do

context("constructors") do
    f1 = Face(1,2,3)
    @fact f1 --> Face{3,Int,0}(1,2,3)
    @fact Face{3,Int,0}(f1) --> f1
    @fact Face{3,UInt8,0}(f1) --> Face{3,UInt8,0}(1,2,3)
    @fact Face{3,Int,-1}(f1) --> Face{3,Int,-1}(0,1,2)
end

context("getindex") do
    a = [1,2,3,4]
    @fact a[Face{3,Int,0}(1,2,3)] --> (1,2,3)
    @fact a[Face{3,Int,-1}(0,1,2)] --> (1,2,3)
    @fact a[Face{3,Int,1}(2,3,4)] --> (1,2,3)
end

context("setindex") do
    a = [1,2,3,4]
    a[Face(1,2,3)] = [7,6,5]
    @fact a --> [7,6,5,4]
end

end


text:
using NCDatasets

"""
    Find indices of particle with lowest mse.
"""
function particle_diags_min_index(ds::NCDataset; metric::String = "mse_full")
    mse_mat = ds.group["particle_diags"][metric][:] # (n_particles, n_iterations)
    replace!(mse_mat, NaN => Inf)
    part_min_i, iter_min_i = Tuple(argmin(mse_mat))
    return part_min_i, iter_min_i
end

"""
    optimal_parameters(ds_path::String; method::String = "best_nn_particle_mean", metric::String = "mse_full")
Given path to Diagnostics.nc file, return optimal parameter set across all iterations and corresponding parameter names.

Inputs:
 - `ds_path`: path to Diagnostics.nc file
 - `method`: method for computing optimal parameters. Returns parameters of:
        "best_particle" - particle with lowest mse in training (`metric` = "mse_full") or validation (`metric` = "val_mse_full") set.
        "best_nn_particle_mean" - particle nearest to ensemble mean for the iteration with lowest mse.
 - `metric` : mse metric to find the minimum of {"mse_full", "val_mse_full"}.
Returns:
    - `u_names`: vector of parameter names
    - `u`: vector of optimal parameter values
"""

function optimal_parameters(ds_path::String; method::String = "best_nn_particle_mean", metric::String = "mse_full")

    valid_metrics = ("mse_full", "val_mse_full")
    @assert metric in valid_metrics "Invalid metric: $metric"

    valid_methods = ("best_nn_particle_mean", "best_particle")
    @assert method in valid_methods "Invalid method: $method"

    if method == "best_nn_particle_mean"
        optimal = NCDataset(ds_path) do ds
            mse_arr = ds.group["metrics"][string(metric, "_nn_mean")][:]
            replace!(mse_arr, NaN => Inf)
            mse_min_i = argmin(mse_arr)
            nn_mean_index = ds.group["metrics"]["nn_mean_index"][mse_min_i]
            u = ds.group["particle_diags"]["phi"][nn_mean_index, :, mse_min_i] # (n_particles, n_params, n_iterations)
            u_names = ds.group["ensemble_diags"]["param"][:]
            @info "Optimal parameters found at: iteration = $mse_min_i ; nn particle = $nn_mean_index"
            (; u_names, u)
        end
    elseif method == "best_particle"
        optimal = NCDataset(ds_path) do ds
            phi = ds.group["particle_diags"]["phi"][:] # (n_particles, n_params, n_iterations)
            part_min_i, iter_min_i = particle_diags_min_index(ds; metric = metric)
            u = phi[part_min_i, :, iter_min_i]
            u_names = ds.group["particle_diags"]["param"][:]
            @info "Optimal parameters found at: iteration = $iter_min_i ; particle = $part_min_i"
            (; u_names, u)
        end
    end

    return optimal.u_names, optimal.u
end

"""
    optimal_mse(ds_path::String; method::String = "best_nn_particle_mean", metric::String = "mse_full")
Given path to Diagnostics.nc file, return minimum mse across iterations.

Inputs:
 - `ds_path`: path to Diagnostics.nc file
 - `method`: method for computing minimum mse. Returns training and validation mse of:
        "best_particle" - particle with lowest mse in training (`metric` = "mse_full") or validation (`metric` = "mse_full_val") set.
        "best_nn_particle_mean" - particle nearest to ensemble mean for the iteration with lowest mse.
 - `metric` : mse metric to find the minimum of {"mse_full", "val_mse_full"}.

Returns:
    - `mse_train`: mse for training set
    - `mse_val`: mse for validation set
"""

function optimal_mse(ds_path::String; method::String = "best_nn_particle_mean", metric::String = "mse_full")

    valid_metrics = ("mse_full", "val_mse_full")
    @assert metric in valid_metrics "Invalid metric: $metric"

    valid_methods = ("best_nn_particle_mean", "best_particle")
    @assert method in valid_methods "Invalid method: $method"

    if method == "best_nn_particle_mean"
        optimal = NCDataset(ds_path) do ds
            mse_arr = ds.group["metrics"][string(metric, "_nn_mean")][:]
            replace!(mse_arr, NaN => Inf)
            mse_min_i = argmin(mse_arr)
            mse_train = ds.group["metrics"]["mse_full_nn_mean"][mse_min_i]
            # if val dataset provided, return val mse
            if ("val_mse_full_nn_mean" in keys(ds.group["metrics"]))
                mse_val = ds.group["metrics"]["val_mse_full_nn_mean"][:][mse_min_i]
            else
                mse_val = nothing
            end
            @info "Optimal mse found at: iteration = $mse_min_i"
            (; mse_train, mse_val)
        end
    elseif method == "best_particle"
        optimal = NCDataset(ds_path) do ds
            part_min_i, iter_min_i = particle_diags_min_index(ds; metric = metric)
            mse_mat = ds.group["particle_diags"]["mse_full"][:] # (n_particles, n_iterations)
            mse_train = mse_mat[part_min_i, iter_min_i]
            # if val dataset provided, return val mse
            if ("val_mse_full" in keys(ds.group["particle_diags"]))
                mse_mat_val = ds.group["particle_diags"]["val_mse_full"][:]
                mse_val = mse_mat_val[part_min_i, iter_min_i]
            else
                mse_val = nothing
            end
            @info "Optimal mse found at: iteration = $iter_min_i ; particle = $part_min_i"
            (; mse_train, mse_val)
        end
    end

    return optimal.mse_train, optimal.mse_val
end

"""
    Given path to Diagnostics.nc file, return mean(mse across particles) in the final iteration for training and validation set.
"""
function final_mse(ds_path::String; metric::String = "mse_full_mean")
    final = NCDataset(ds_path) do ds
        mse_mean = ds.group["metrics"][metric][:]
        mse_train = mse_mean[length(mse_mean) - 1]
        val_metric = string("val_", metric)
        # if val dataset provided, return val mse
        if (val_metric in keys(ds.group["metrics"]))
            ds_mse_val_mean = ds.group["metrics"][val_metric][:]
            mse_val = ds_mse_val_mean[length(ds_mse_val_mean) - 1]
        else
            mse_val = nothing
        end
        (; mse_train, mse_val)
    end

    return final.mse_train, final.mse_val
end

"""
    Given path to Diagnostics.nc file, return parameter names and mean values in final iteration.
"""
function final_parameters(ds_path::String)
    final = NCDataset(ds_path) do ds
        ds_phi_mean = ds.group["ensemble_diags"]["phi_mean"]
        u = ds_phi_mean[:, size(ds_phi_mean, 2) - 1]
        u_names = ds.group["ensemble_diags"]["param"][:]
        (; u_names, u)
    end
    return final.u_names, final.u
end


text:
# dense Nauty graphs

WORDSIZE == 64 || error("WORDSIZE is not 64. Probably all hell will break loose.")

export DenseNautyGraph, DenseNautyDiGraph

"""
A dense Nauty graph is a bit matrix: if the graph has n vertices,
then it is coded as an (WORDSIZE*m)×n bit matrix, with m=(n+WORDSIZE-1)÷WORDSIZE.
Each column of that matrix is thus a chunk of memory of size m, and the ith column represents the adjacency list of the the ith vertex. The bits are reversed in each word, so to know if the ith vertex is connected to vertex j, check data[setpos(j),i].
"""
struct DenseNautyGraph <: AbstractGraph{Int}
    data::BitMatrix
end
LightGraphs.is_directed(::Type{DenseNautyGraph}) = false
LightGraphs.is_directed(::DenseNautyGraph) = false

struct DenseNautyDiGraph <: AbstractGraph{Int}
    data::BitMatrix
end
LightGraphs.is_directed(::Type{DenseNautyDiGraph}) = true
LightGraphs.is_directed(::DenseNautyDiGraph) = true

DenseNautyXGraph = Union{DenseNautyGraph,DenseNautyDiGraph}

num_setwords(n::Int) = (n+WORDSIZE-1)÷WORDSIZE
setpos(n::Int) = 1+((n-1)⊻0x3f)

DenseNautyGraph(n::Int) = (data = BitMatrix(undef,WORDSIZE*num_setwords(n),n); fill!(data,false); DenseNautyGraph(data))

function DenseNautyGraph(mat::Union{Matrix,BitMatrix})
    n = size(mat,1)
    mat == mat' || error("mat should be a square, symmetric matrix, $mat, $(mat')")

    data = BitMatrix(undef,WORDSIZE*num_setwords(n),n)
    fill!(data,false)
    for s=1:n, d=1:n
        data[setpos(d),s] = mat[s,d]
    end
    DenseNautyGraph(data)
end

function DenseNautyDiGraph(mat::Union{Matrix{Int},Matrix{Bool},BitMatrix})
    n = size(mat,1)
    n == size(mat,2) || error("mat should be a square matrix")

    data = BitMatrix(undef,WORDSIZE*num_setwords(n),n)
    fill!(data,false)
    for s=1:n, d=1:n
        data[setpos(d),s] = mat[s,d]
    end
    DenseNautyDiGraph(data)
end

function DenseNautyGraph(g::AbstractGraph{T}) where T <: Integer
    ng = DenseNautyGraph(nv(g))
    for e=edges(g)
        add_edge!(ng,e)
    end
    ng
end

function DenseNautyDiGraph(g::AbstractGraph{T}) where T <: Integer
    ng = DenseNautyDiGraph(nv(g))
    for e=edges(g)
        add_edge!(ng,e)
    end
    ng
end

convert(::Type{DenseNautyGraph}, g) = DenseNautyGraph(g)
convert(::Type{DenseNautyDiGraph}, g) = DenseNautyDiGraph(g)

Base.:(==)(g::DenseNautyXGraph, h::DenseNautyXGraph) = g.data == h.data
Base.hash(g::DenseNautyXGraph) = hash(g.data)
Base.copy(g::DenseNautyXGraph) = typeof(g)(copy(g.data))

Base.eltype(g::DenseNautyXGraph) = Int
LightGraphs.edgetype(g::DenseNautyXGraph) = SimpleEdge{eltype(g)}

LightGraphs.nv(g::DenseNautyXGraph) = size(g.data,2)

function LightGraphs.ne(g::DenseNautyGraph)
    n = 0
    for i=1:nv(g), j=1:i
        if g.data[setpos(i),j] || g.data[setpos(j),i]
            n += 1
        end
    end
    n
end
LightGraphs.ne(g::DenseNautyDiGraph) = count(g.data)

Base.zero(g::DenseNautyGraph) = DenseNautyGraph(BitMatrix(undef,0,0))
Base.zero(g::DenseNautyDiGraph) = DenseNautyDiGraph(BitMatrix(undef,0,0))

Base.reverse(g::DenseNautyGraph) = g
function Base.reverse(g::DenseNautyDiGraph)
    n = nv(g)
    data = BitMatrix(undef,WORDSIZE*num_setwords(n),n)
    fill!(data,false)
    for s=1:n, d=1:n
        data[setpos(d),s] = g.data[setpos(s),d]
    end
    DenseNautyDiGraph(data)
end

LightGraphs.edges(g::DenseNautyXGraph) = DenseNautyEdgeIter(g)
struct DenseNautyEdgeIter <: AbstractEdgeIter
    g::DenseNautyXGraph
end
function Base.iterate(eit::DenseNautyEdgeIter, state=(1,0))
    g = eit.g
    n = nv(g)
    s, d = state
    while true
        d += 1
        if d > n
            s += 1
            d = 1
        end
        if s > n
            return nothing
        end
        if g.data[setpos(d),s]
            return Edge(s,d), (s,d)
        end
    end
end

LightGraphs.has_edge(g::DenseNautyXGraph, s::Int, d::Int) = g.data[setpos(d),s]

function LightGraphs.inneighbors(g::DenseNautyXGraph, v::Int)
    n = Int[]
    for i=1:nv(g)
        g.data[setpos(v),i] && push!(n,i)
    end
    n
end

function LightGraphs.outneighbors(g::DenseNautyXGraph, v::Int)
    n = Int[]
    for i=1:nv(g)
        g.data[setpos(i),v] && push!(n,i)
    end
    n
end

LightGraphs.add_edge!(g::DenseNautyGraph, s::Int, d::Int) = (g.data[setpos(d),s] = g.data[setpos(s),d] = true)
LightGraphs.add_edge!(g::DenseNautyDiGraph, s::Int, d::Int) = (g.data[setpos(d),s] = true)
LightGraphs.add_edge!(g::DenseNautyXGraph, e::Edge) = add_edge!(g, e.src, e.dst)

LightGraphs.rem_edge!(g::DenseNautyGraph, s::Int, d::Int) = (g.data[setpos(d),s] = g.data[setpos(s),d] = false; true)
LightGraphs.rem_edge!(g::DenseNautyDiGraph, s::Int, d::Int) = (g.data[setpos(d),s] = false; true)
LightGraphs.rem_edge!(g::DenseNautyXGraph, e::Edge) = rem_edge!(g, e.src, e.dst)


text:
@testset "ExternalSamplingMCCFR" begin
    env = OpenSpielEnv("kuhn_poker";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = ExternalSamplingMCCFRPolicy(rng=StableRNG(123))
    run(p, env, StopAfterStep(1000))
    @test_broken RLZoo.nash_conv(p, env) < 0.05

    env = OpenSpielEnv("leduc_poker";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = ExternalSamplingMCCFRPolicy(rng=StableRNG(123))
    run(p, env, StopAfterStep(1000))
    @test RLZoo.nash_conv(p, env) < 2.5

    env = OpenSpielEnv("liars_dice";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = ExternalSamplingMCCFRPolicy(rng=StableRNG(123))
    run(p, env, StopAfterStep(100))
    @test RLZoo.nash_conv(p, env) < 1.6

    env = OpenSpielEnv("kuhn_poker(players=3)";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = ExternalSamplingMCCFRPolicy(rng=StableRNG(123))
    run(p, env, StopAfterStep(100))
    @info "nash_conv for kuhn_poker(players=3)" RLZoo.nash_conv(p, env)
end

@testset "OutcomeSamplingMCCFR" begin
    env = OpenSpielEnv("kuhn_poker";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = OutcomeSamplingMCCFRPolicy(rng=StableRNG(123))
    run(p, env, StopAfterStep(10000))
    @test_broken RLZoo.nash_conv(p, env) < 0.04

    env = OpenSpielEnv("leduc_poker";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = OutcomeSamplingMCCFRPolicy(rng=StableRNG(123))
    run(p, env, StopAfterStep(10000))
    @test RLZoo.nash_conv(p, env) < 3

    env = OpenSpielEnv("liars_dice";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = OutcomeSamplingMCCFRPolicy(rng=StableRNG(123))
    run(p, env, StopAfterStep(1000))
    @test RLZoo.nash_conv(p, env) < 1.7

    env = OpenSpielEnv("kuhn_poker(players=3)";default_state_style=RLBase.Information{String}(), is_chance_agent_required=true)
    p = OutcomeSamplingMCCFRPolicy(rng=StableRNG(123))
    run(p, env, StopAfterStep(1000))
    @info "nash_conv for kuhn_poker(players=3)" RLZoo.nash_conv(p, env)
end


text:
using JLD2
using FileIO
using QXTns

@testset "Test tensor cache" begin
    tc = QXTools.TensorCache()
    a = rand(Float64, 3, 4, 5)
    sym = push!(tc, a)
    @test sym == push!(tc, a)
    @test sym != push!(tc, rand(Float64, 3, 4, 5))
    # test approximate matching
    @test sym == push!(tc, a .+ 0.1*eps(Float64))
    @test sym != push!(tc, a .+ 1.1*eps(Float64))

    @test tc[sym] == a
    @test length(tc) == 3

    # test saving tensor cache
    mktempdir() do path
        save_cache(tc, joinpath(path, "tmp.jld2"))
        loaded_data = load(joinpath(path, "tmp.jld2"))
        @test all([(loaded_data[x] == tc[Symbol(x)]) for x in keys(loaded_data)])
    end
end

@testset "Test build contraction tree" begin
    circ = create_test_circuit()
    tnc = convert_to_tnc(circ)
    plan = min_fill_contraction_plan(tnc)
    cg = build_compute_graph(tnc, plan)
    # test that number of nodes in compute graph is
    # tensors in network + length of plan + 1 for save node
    @test length(cg.root) == length(tnc) + length(plan) + 1
end


text:
using Base.Test

randop() = rand(["-->", "→",
                 "||",
                 "&&",
                 "<", "==", "<:", ">:",
                 "<|", "|>",
                 ":",
                 "+", "-",
                 ">>", "<<",
                 "*", "/",
                 "//",
                 "^", "↑",
                 "::",
                 ".", "->"])

function test_expr(str)
    x, ps = CSTParser.parse(ParseState(str))

    x0 = Expr(x)
    x1 = remlineinfo!(Base.parse(str))
    !ps.errored && x0 == x1# && isempty(span(x))
end

@testset "All tests" begin
@testset "Operators" begin
    @testset "Binary Operators" begin
        for iter = 1:250
            str = join([["x$(randop())" for i = 1:19];"x"])

            @test test_expr(str)
        end
    end
    @testset "Conditional Operator" begin
        strs = ["a ? b : c"
                "a ? b:c : d"
                "a ? b:c : d:e"]
        for str in strs
            @test test_expr(str)
        end
    end

    @testset "Dot Operator" begin
        @test "a.b"  |> test_expr
        @test "a.b.c"  |> test_expr
        @test "(a(b)).c"  |> test_expr
        @test "(a).(b).(c)"  |> test_expr
        @test "(a).b.(c)"  |> test_expr
        @test "(a).b.(c+d)"  |> test_expr
    end

    @testset "Unary Operator" begin
        @test "+" |> test_expr
        @test "-" |> test_expr
        @test "!" |> test_expr
        @test "~" |> test_expr
        @test "&" |> test_expr
        @test "::" |> test_expr
        @test "<:" |> test_expr
        @test ">:" |> test_expr
        @test "¬" |> test_expr
        @test "√" |> test_expr
        @test "∛" |> test_expr
        @test "∜" |> test_expr
    end

    @testset "Unary Operator" begin
        @test "a=b..." |> test_expr
        @test "a-->b..." |> test_expr
        @test "a&&b..." |> test_expr
        @test "a||b..." |> test_expr
        @test "a<b..." |> test_expr
        @test "a:b..." |> test_expr
        @test "a+b..." |> test_expr
        @test "a<<b..." |> test_expr
        @test "a*b..." |> test_expr
        @test "a//b..." |> test_expr
        @test "a^b..." |> test_expr
        @test "a::b..." |> test_expr
        @test "a where b..." |> test_expr
        @test "a.b..." |> test_expr
    end

    @testset "unary op calls" begin
        @test "+(a,b)" |> test_expr
        @test "-(a,b)" |> test_expr
        @test "!(a,b)" |> test_expr
        @test "¬(a,b)" |> test_expr
        @test "~(a,b)" |> test_expr
        @test_broken "<:(a,b)" |> test_expr
        @test "√(a,b)" |> test_expr
        @test "\$(a,b)" |> test_expr
        @test ":(a,b)" |> test_expr
        @test "&a" |> test_expr
        @test "&(a,b)" |> test_expr
        @test "::a" |> test_expr
        @test "::(a,b)" |> test_expr
    end

    @testset "where precedence" begin
        @test "a = b where c = d" |> test_expr
        @test "a = b where c" |> test_expr
        @test "b where c = d" |> test_expr

        @test "a ? b where c : d" |> test_expr

        @test "a --> b where c --> d" |> test_expr
        @test "a --> b where c" |> test_expr
        @test "b where c --> d" |> test_expr

        @test "a || b where c || d" |> test_expr
        @test "a || b where c" |> test_expr
        @test "b where c || d" |> test_expr

        @test "a && b where c && d" |> test_expr
        @test "a && b where c" |> test_expr
        @test "b where c && d" |> test_expr

        @test "a <: b where c <: d" |> test_expr
        @test "a <: b where c" |> test_expr
        @test "b where c <: d" |> test_expr

        @test "a <| b where c <| d" |> test_expr
        @test "a <| b where c" |> test_expr
        @test "b where c <| d" |> test_expr

        @test "a : b where c : d" |> test_expr
        @test "a : b where c" |> test_expr
        @test "b where c : d" |> test_expr

        @test "a + b where c + d" |> test_expr
        @test "a + b where c" |> test_expr
        @test "b where c + d" |> test_expr

        @test "a << b where c << d" |> test_expr
        @test "a << b where c" |> test_expr
        @test "b where c << d" |> test_expr

        @test "a * b where c * d" |> test_expr
        @test "a * b where c" |> test_expr
        @test "b where c * d" |> test_expr

        @test "a // b where c // d" |> test_expr
        @test "a // b where c" |> test_expr
        @test "b where c // d" |> test_expr

        @test "a ^ b where c ^ d" |> test_expr
        @test "a ^ b where c" |> test_expr
        @test "b where c ^ d" |> test_expr

        @test "a :: b where c :: d" |> test_expr
        @test "a :: b where c" |> test_expr
        @test "b where c :: d" |> test_expr

        @test "a.b where c.d" |> test_expr
        @test "a.b where c" |> test_expr
        @test "b where c.d" |> test_expr

        @test "a where b where c" |> test_expr
    end
end

@testset "Type Annotations" begin
    @testset "Curly" begin
        @test "x{T}" |> test_expr
        @test "x{T,S}" |> test_expr
        @test "a.b{T}" |> test_expr
        @test "a(b){T}" |> test_expr
        @test "(a(b)){T}" |> test_expr
        @test "a{b}{T}" |> test_expr
        @test "a{b}(c){T}" |> test_expr
        @test "a{b}.c{T}" |> test_expr
        @test """x{T,
        S}""" |> test_expr
    end
end

@testset "Tuples" begin
    @test "1," |> test_expr
    @test "1,2" |> test_expr
    @test "1,2,3" |> test_expr
    @test "()" |> test_expr
    @test "(==)" |> test_expr
    @test "(1)" |> test_expr
    @test "(1,)" |> test_expr
    @test "(1,2)" |> test_expr
    @test "(a,b,c)" |> test_expr
    @test "(a...)" |> test_expr
    @test "((a,b)...)" |> test_expr
    @test "a,b = c,d" |> test_expr
    @test "(a,b) = (c,d)" |> test_expr
end

@testset "Function Calls" begin
    @testset "Simple Calls" begin
        @test "f(x)" |> test_expr
        @test "f(x,y)" |> test_expr
        @test "f(g(x))" |> test_expr
        @test "f((x,y))" |> test_expr
        @test "f((x,y), z)" |> test_expr
        @test "f(z, (x,y), z)" |> test_expr
        @test "f{a}(x)" |> test_expr
        @test "f{a<:T}(x::T)" |> test_expr
    end

    @testset "Keyword Arguments" begin
        @test "f(x=1)" |> test_expr
        @test "f(x=1,y::Int = 1)" |> test_expr
    end

    @testset "Compact Declaration" begin
        @test "f(x) = x" |> test_expr
        @test "f(x) = g(x)" |> test_expr
        @test "f(x) = (x)" |> test_expr
        @test "f(x) = (x;y)" |> test_expr
        @test "f(g(x)) = x" |> test_expr
        @test "f(g(x)) = h(x)" |> test_expr
    end

    @testset "Standard Declaration" begin
        @test "function f end" |> test_expr
        @test "function f(x) x end" |> test_expr
        @test "function f(x); x; end" |> test_expr
        @test "function f(x) x; end" |> test_expr
        @test "function f(x); x end" |> test_expr
        @test "function f(x) x;y end" |> test_expr
        @test """function f(x)
            x
        end""" |> test_expr
        @test """function f(x,y =1)
            x
        end""" |> test_expr
        @test """function f(x,y =1;z =2)
            x
        end""" |> test_expr
    end
    @testset "Anonymous" begin
        @test "x->y" |> test_expr
        @test "(x,y)->x*y" |> test_expr
        @test """function ()
            return
        end""" |> test_expr
    end
end

@testset "Types" begin
    @testset "Abstract" begin
        @test "abstract t" |> test_expr
        @test "abstract t{T}" |> test_expr
        @test "abstract t <: S" |> test_expr
        @test "abstract t{T} <: S" |> test_expr
    end

    @testset "Bitstype" begin
        @test "bitstype 64 Int" |> test_expr
        @test "bitstype 4*16 Int" |> test_expr
    end

    @testset "Typealias" begin
        @test "typealias name fsd" |> test_expr
    end

    @testset "Structs" begin
        @test "type a end" |> test_expr
        @test "type a; end" |> test_expr
        @test "type a; b;end" |> test_expr
        @test """type a
            arg1
        end""" |> test_expr
        @test """type a <: T
            arg1::Int
            arg2::Int
        end""" |> test_expr
        @test """type a
            arg1::T
        end""" |> test_expr
        @test """type a{T}
            arg1::T
            a(args) = new(args)
        end""" |> test_expr
        @test """type a <: Int
            arg1::Vector{Int}
        end""" |> test_expr
        @test """immutable a <: Int
            arg1::Vector{Int}
        end""" |> test_expr
    end
end

@testset "Modules" begin
    @testset "Imports " begin
        @test "import ModA" |> test_expr
        @test "import .ModA" |> test_expr
        @test "import ..ModA.a" |> test_expr
        @test "import ModA.subModA" |> test_expr
        @test "import ModA.subModA: a" |> test_expr
        @test "import ModA.subModA: a, b" |> test_expr
        @test "import ModA.subModA: a, b.c" |> test_expr
        @test "import .ModA.subModA: a, b.c" |> test_expr
        @test "import ..ModA.subModA: a, b.c" |> test_expr
    end
    @testset "Export " begin
        @test "export ModA" |> test_expr
        @test "export a, b, c" |> test_expr
    end
end

@testset "Generators" begin
    @test "(y for y in X)" |> test_expr
    @test "((x,y) for x in X, y in Y)" |> test_expr
    @test "(y.x for y in X)" |> test_expr
    @test "((y) for y in X)" |> test_expr
    @test "(y,x for y in X)" |> test_expr
    @test "((y,x) for y in X)" |> test_expr
    @test "[y for y in X]" |> test_expr
    @test "[(y) for y in X]" |> test_expr
    @test "[(y,x) for y in X]" |> test_expr
    @test "Int[y for y in X]" |> test_expr
    @test "Int[(y) for y in X]" |> test_expr
    @test "Int[(y,x) for y in X]" |> test_expr
    @test """
    [a
    for a = 1:2]
    """ |> test_expr
    @test "[ V[j][i]::T for i=1:length(V[1]), j=1:length(V) ]" |> test_expr
    @test "all(d ≥ 0 for d in B.dims)" |> test_expr
end

@testset "Macros " begin
    @test "@mac" |> test_expr
    @test "@mac a b c" |> test_expr
    @test "@mac f(5)" |> test_expr
    @test "(@mac x)" |> test_expr
    @test "Mod.@mac a b c" |> test_expr
    # @test "[@mac a b]" |> test_expr
    @test "@inline get_chunks_id(i::Integer) = _div64(Int(i)-1)+1, _mod64(Int(i) -1)" |> test_expr
    @test "@inline f() = (), ()" |> test_expr
end

@testset "Square " begin
    @testset "vect" begin
        @test "[x]" |> test_expr
        @test "[(1,2)]" |> test_expr
        @test "[x...]" |> test_expr
        @test "[1,2,3,4,5]" |> test_expr
    end

    @testset "ref" begin
        @test "t[i]" |> test_expr
        @test "t[i, j]" |> test_expr
    end

    @testset "vcat" begin
        @test "[x;]" |> test_expr
        @test "[x;y;z]" |> test_expr
        @test """[x
                  y
                  z]""" |> test_expr
        @test """[x
                  y;z]""" |> test_expr
        @test """[x;y
                  z]""" |> test_expr
        @test "[x,y;z]" |> test_expr
    end

    @testset "typed_vcat" begin
        @test "t[x;]" |> test_expr
        @test "t[x;y]" |> test_expr
        @test """t[x
                   y]""" |> test_expr
        @test "t[x;y]" |> test_expr
        @test "t[x y; z]" |> test_expr
        @test "t[x, y; z]" |> test_expr
    end

    @testset "hcat" begin
        @test "[x y]" |> test_expr
    end

    @testset "typed_hcat" begin
        @test "t[x y]" |> test_expr
    end

    @testset "Comprehension" begin
        @test "[i for i = 1:10]" |> test_expr
        @test "Int[i for i = 1:10]" |> test_expr
    end
end

@testset "Keyword Blocks" begin
    @testset "If" begin
        @test "if cond end" |> test_expr
        @test "if cond; a; end" |> test_expr
        @test "if cond a; end" |> test_expr
        @test "if cond; a end" |> test_expr
        @test """if cond
            1
            1
        end""" |> test_expr
        @test """if cond
        else
            2
            2
        end""" |> test_expr
        @test """if cond
            1
            1
        else
            2
            2
        end""" |> test_expr
        @test "if 1<2 end" |> test_expr
        @test """if 1<2
            f(1)
            f(2)
        end""" |> test_expr
        @test """if 1<2
            f(1)
        elseif 1<2
            f(2)
        end""" |> test_expr
        @test """if 1<2
            f(1)
        elseif 1<2
            f(2)
        else
            f(3)
        end""" |> test_expr
        @test "if cond a end" |> test_expr
    end

    @testset "Try" begin
        @test "try f(1) end" |> test_expr
        @test "try; f(1) end" |> test_expr
        @test "try; f(1); end" |> test_expr
        @test "try; f(1); catch e; e; end" |> test_expr
        @test "try; f(1); catch e; e end" |> test_expr
        @test "try; f(1); catch e e; end" |> test_expr
        @test """try
            f(1)
        catch
        end""" |> test_expr
        @test """try
            f(1)
        catch
            error(err)
        end""" |> test_expr
        @test """try
            f(1)
        catch err
            error(err)
        end""" |> test_expr
        @test """try
            f(1)
        catch
            error(err)
        finally
            stop(f)
        end""" |> test_expr
        @test """try
            f(1)
        catch err
            error(err)
        finally
            stop(f)
        end""" |> test_expr
        @test """try
            f(1)
        finally
            stop(f)
        end""" |> test_expr
    end
    @testset "For" begin
        @test """for i = 1:10
            f(i)
        end""" |> test_expr
        @test """for i = 1:10, j = 1:20
            f(i)
        end""" |> test_expr
    end

    @testset "Let" begin
        @test """let x = 1
            f(x)
        end""" |> test_expr
        @test """let x = 1, y = 2
            f(x)
        end""" |> test_expr
    end

    @testset "Do" begin
        @test """f(X) do x
            return x
        end""" |> test_expr
        @test """f(X,Y) do x,y
            return x,y
        end""" |> test_expr
    end
end

@testset "No longer broken things" begin
    @test "[ V[j][i]::T for i=1:length(V[1]), j=1:length(V) ]" |> test_expr
    @test "all(d ≥ 0 for d in B.dims)" |> test_expr
    @test ":(=)" |> test_expr
    @test ":(1)" |> test_expr
    @test ":(a)" |> test_expr
    @test "(@_inline_meta(); f(x))" |> test_expr
    @test "isa(a,b) != c" |> test_expr
    @test "isa(a,a) != isa(a,a)" |> test_expr
    @test "@mac return x" |> test_expr
    @test "a,b," |> test_expr
    @test "m!=m" |> test_expr
    @test "+(x...)" |> test_expr
    @test "+(promote(x,y)...)" |> test_expr
    @test "\$(x...)" |> test_expr
    @test "ccall(:gethostname, stdcall, Int32, ())" |> test_expr
    @test "@inbounds @ncall a b c" |> test_expr
    @test "(a+b)``" |> test_expr
    @test "(-, ~)" |> test_expr
    @test """function +(x::Bool, y::T)::promote_type(Bool,T) where T<:AbstractFloat
                return ifelse(x, oneunit(y) + y, y)
            end""" |> test_expr
    @test """finalizer(x,x::GClosure->begin
                    ccall((:g_closure_unref,Gtk.GLib.libgobject),Void,(Ptr{GClosure},),x.handle)
                end)""" |> test_expr
    @test "function \$A end" |> test_expr
    @test "&ctx->exe_ctx_ref" |> test_expr
    @test ":(\$(docstr).\$(TEMP_SYM)[\$(key)])" |> test_expr
    @test "SpecialFunctions.\$(fsym)(n::Dual)" |> test_expr
    @test "(Base.@_pure_meta;)" |> test_expr
    @test "@M a b->(@N c = @O d e f->g)" |> test_expr
    @test "! = f" |> test_expr
    @test "[a=>1, b=>2]" |> test_expr
    @test "a.\$(b)" |> test_expr
    @test "a.\$f()" |> test_expr
    @test "4x/y" |> test_expr
    @test """
    ccall(:jl_finalize_th, Void, (Ptr{Void}, Any,),
                Core.getptls(), o)
    """ |> test_expr
    @test """
    A[if n == d
        i
    else
        (indices(A,n) for n = 1:nd)
    end...]
    """ |> test_expr
    @test """Base.@__doc__(bitstype \$(sizeof(basetype) * 8) \$(esc(typename)) <: Enum{\$(basetype)})""" |> test_expr
    @test """
    @spawnat(p,
        let m = a
            isa(m, Exception) ? m : nothing
        end)
    """ |> test_expr
    @test "[@spawn f(R, first(c), last(c)) for c in splitrange(length(R), nworkers())]" |> test_expr
    @test "M.:(a)" |> test_expr
    @test """
            begin
                for i in I for j in J
                    if cond
                        a
                    end
                end end
            end""" |> test_expr
    @test "-f.(a.b + c)" |> test_expr
    @test ":(import Base: @doc)" |> test_expr
    @test "[a for a in A for b in B]" |> test_expr
    @test "+(a,b,c...)" |> test_expr
    @test """@testset a for t in T
        t
    end""" |> test_expr
    @test "import Base.==" |> test_expr
    @test "a`text`" |> test_expr
    @test "a``" |> test_expr
    @test "a`text`b" |> test_expr
    @test "[a; a 0]" |> test_expr
    @test "[a, b; c]" |> test_expr
    @test "t{a; b} " |> test_expr
    @test "a ~ b + c -d" |> test_expr
    @test "y[j=1:10,k=3:2:9; isodd(j+k) && k <= 8]" |> test_expr
    @test "(8=>32.0, 12=>33.1, 6=>18.2)" |> test_expr
    @test "(a,b = c,d)" |> test_expr
    @test "[ -1 -2;]" |> test_expr
    @test "-2y" |> test_expr # precedence
    @test "'''" |> test_expr #tokenize
    @test """
        if j+k <= deg +1
        end
        """ |> test_expr
    @test "function f() ::T end" |> test_expr # ws closer
    @test "import Base: +, -, .+, .-" |> test_expr
    @test "[a +   + l]" |> test_expr #ws closer
    @test "@inbounds C[i,j] = - α[i] * αjc" |> test_expr
    @test "@inbounds C[i,j] = - n * p[i] * pj" |> test_expr
    @test """
        if ! a
            b
        end
        """ |> test_expr # ws closer
    @test "[:-\n:+]" |> test_expr
    @test "::a::b" |> test_expr
    @test "-[1:nc]" |> test_expr
    @test "f() where {a} = b" |> test_expr
    @test "@assert .!(isna(res[2]))" |> test_expr # v0.6
    @test "-((attr.rise / PANGO_SCALE)pt).value" |> test_expr
    @test "!(a = b)" |> test_expr
    @test "-(1)a" |> test_expr
    @test "(Base.@_pure_meta;)" |> test_expr
    @test "!(a)::T" |> test_expr
    @test "a::b where T<:S" |> test_expr
    @test "+(x::Bool, y::T)::promote_type(Bool,T) where T<:AbstractFloat" |> test_expr
    @test "T where V<:(T where T)" |> test_expr
    @test "function ^(z::Complex{T}, p::Complex{T})::Complex{T} where T<:AbstractFloat end" |> test_expr
    @test "function +(a) where T where S end" |> test_expr
    @test "function -(x::Rational{T}) where T<:Signed end" |> test_expr
    @test "\$(a)(b)" |> test_expr
    @test "if !(a) break end" |> test_expr
    @test "module a() end" |> test_expr
    @test "M.r\"str\" " |> test_expr
    @test "f(a for a in A if cond)" |> test_expr
    @test "\"dimension \$d is not 1 ≤ \$d ≤ \$nd\" " |> test_expr
    @test "-(-x)^1" |> test_expr
end

@testset "Broken things" begin
    @test_broken "\$(a) * -\$(b)" |> test_expr
    @test_broken "function(f, args...; kw...) end" |> test_expr
    @test_broken """
        "\\\\\$ch"
        """ |> test_expr
    @test_broken "µs" |> test_expr # normalize unicode
    @test_broken """let f = ((; a = 1, b = 2) -> ()),
                m = first(methods(f))
                @test DSE.keywords(f, m) == [:a, :b]
            end""" |> test_expr
    @test_broken "-1^a" |> test_expr
end

test_fsig_decl(str) = (x->x.id).(CSTParser._get_fsig(CSTParser.parse(str)).defs)
@testset "func-sig variable declarations" begin
    @test test_fsig_decl("f(x) = x") == [:x]
    @test test_fsig_decl("""function f(x)
        x
    end""") == [:x]

    @test test_fsig_decl("f{T}(x::T) = x") == [:T, :x]
    @test test_fsig_decl("""function f{T}(x::T)
        x
    end""") == [:T, :x]

    @test test_fsig_decl("f(x::T) where T = x") == [:T, :x]
    @test test_fsig_decl("""function f(x::T) where T
        x
    end""") == [:T, :x]

    @test test_fsig_decl("f(x::T{S}) where T where S = x") == [:T, :S, :x]
    @test test_fsig_decl("""function f(x::T{S}) where T where S
        x
    end""") == [:T, :S, :x]
end

end


text:
#=
Motivation:
From GoF book, a class TCPConnection responds different to requests
depending on its current state of Established, Listening, or Closed.
=#
module StateExample

export LISTENING, ESTABLISHED, CLOSED

abstract type AbstractState end

struct ListeningState <: AbstractState end
struct EstablishedState <: AbstractState end
struct ClosedState <: AbstractState end

const LISTENING = ListeningState()
const ESTABLISHED = EstablishedState()
const CLOSED = ClosedState()

"""
A Connection object contains the current `state` and a reference
to a connection `conn`
"""
struct Connection{T <: AbstractState,S}
    state::T
    conn::S
end

# Use multiple dispatch
send(c::Connection, msg) = send(c.state, c.conn, msg)

# Implement `send` method for each state
send(::ListeningState, conn, msg) = error("No connection yet")
send(::EstablishedState, conn, msg) = write(conn, msg * "\n")
send(::ClosedState, conn, msg) = error("Connection already closed")

# test functions
function test(state, msg)
    c = Connection(state, stdout)
    try
        send(c, msg)
    catch ex
        println("$(ex) for message '$msg'")
    end
    return nothing
end

function test()
    test(LISTENING, "hello world 1")
    test(CLOSED, "hello world 2")
    test(ESTABLISHED, "hello world 3")
end

end #module

using .StateExample
StateExample.test()