text:
# global cache
# ============

const ANALYZED_LINFOS  = IdSet{MethodInstance}() # keeps `MethodInstance`s analyzed by JET

const JET_GLOBAL_CACHE = IdDict{MethodInstance,Vector{InferenceErrorReportCache}}()

function CC.code_cache(interp::JETInterpreter)
    cache  = JETGlobalCache(interp, code_cache(interp.native))
    worlds = WorldRange(get_world_counter(interp))
    return WorldView(cache, worlds)
end

struct JETGlobalCache
    interp::JETInterpreter
    native::WorldView{InternalCodeCache}
end

function CC.haskey(wvc::WorldView{JETGlobalCache}, mi::MethodInstance)
    ret = CC.haskey(WorldView(wvc.cache.native, wvc.worlds), mi)
    if ret && !(mi in ANALYZED_LINFOS)
        add_jet_callback!(mi) # XXX: forcibly register a callback for caches in system image
    end
    return ret
end

function CC.get(wvc::WorldView{JETGlobalCache}, mi::MethodInstance, default)
    # ignore code cache for a `MethodInstance` that is not yet analyzed by JET;
    # this happens when the `MethodInstance` is cached within the system image
    mi in ANALYZED_LINFOS || return default # force inference

    ret = CC.get(WorldView(wvc.cache.native, wvc.worlds), mi, default)
    if isa(ret, CodeInstance)
        # cache hit, now we need to append cached reports associated with this `MethodInstance`
        global_cache = get(JET_GLOBAL_CACHE, mi, nothing)
        if isa(global_cache, Vector{InferenceErrorReportCache})
            interp = wvc.cache.interp
            for cached in global_cache
                restored = restore_cached_report!(cached, interp)
                push!(interp.to_be_updated, restored) # should be updated in `abstract_call` (after exiting `typeinf_edge`)
                # # TODO make this hold
                # @assert first(cached.st).linfo === mi "invalid global restoring"
            end
        end
    end
    return ret
end

function CC.getindex(wvc::WorldView{JETGlobalCache}, mi::MethodInstance)
    r = CC.get(wvc, mi, nothing)
    r === nothing && throw(KeyError(mi))
    return r::CodeInstance
end

function CC.setindex!(wvc::WorldView{JETGlobalCache}, ci::CodeInstance, mi::MethodInstance)
    CC.setindex!(WorldView(wvc.cache.native, wvc.worlds), ci, mi)
    add_jet_callback!(mi)
    return nothing
end

function add_jet_callback!(linfo)
    @static BACKEDGE_CALLBACK_ENABLED || return nothing

    if !isdefined(linfo, :callbacks)
        linfo.callbacks = Any[invalidate_jet_cache!]
    else
        if !any(@nospecialize(cb)->cb===invalidate_jet_cache!, linfo.callbacks)
            push!(linfo.callbacks, invalidate_jet_cache!)
        end
    end
    return nothing
end

function invalidate_jet_cache!(replaced, max_world, depth = 0)
    delete!(ANALYZED_LINFOS, replaced)
    delete!(JET_GLOBAL_CACHE, replaced)

    if isdefined(replaced, :backedges)
        for mi in replaced.backedges
            mi = mi::MethodInstance
            (mi in ANALYZED_LINFOS || haskey(JET_GLOBAL_CACHE, mi)) || continue # otherwise infinite loop
            invalidate_jet_cache!(mi, max_world, depth+1)
        end
    end
    return nothing
end

# local cache
# ===========

struct JETLocalCache
    interp::JETInterpreter
    cache::Vector{InferenceResult}
end

CC.get_inference_cache(interp::JETInterpreter) = JETLocalCache(interp, get_inference_cache(interp.native))

function CC.cache_lookup(linfo::MethodInstance, given_argtypes::Vector{Any}, cache::JETLocalCache)
    inf_result = cache_lookup(linfo, given_argtypes, cache.cache)

    isa(inf_result, InferenceResult) || return nothing

    # cache hit, try to restore local report caches
    interp = cache.interp
    sv = interp.current_frame::InferenceState

    # constant prop' hits a cycle (recur into same non-constant analysis), we just bail out
    isa(inf_result.result, InferenceState) && return inf_result

    analysis_result = jet_cache_lookup(linfo, given_argtypes, interp.cache)

    # corresponds to report throw away logic in `_typeinf(interp::JETInterpreter, frame::InferenceState)`
    filter!(!is_from_same_frame(sv.linfo, linfo), interp.reports)

    if isa(analysis_result, AnalysisResult)
        for cached in analysis_result.cache
            restored = restore_cached_report!(cached, interp)
            push!(interp.to_be_updated, restored) # should be updated in `abstract_call_method_with_const_args`
            # # TODO make this hold
            # @assert first(cached.st).linfo === linfo "invalid local restoring"
        end
    end

    return inf_result
end

# should sync with the implementation of `cache_lookup(linfo::MethodInstance, given_argtypes::Vector{Any}, cache::Vector{InferenceResult})`
function jet_cache_lookup(linfo::MethodInstance, given_argtypes::Vector{Any}, cache::Vector{AnalysisResult})
    method = linfo.def::Method
    nargs::Int = method.nargs
    method.isva && (nargs -= 1)
    length(given_argtypes) >= nargs || return nothing
    for cached_result in cache
        cached_result.linfo === linfo || continue
        cache_match = true
        cache_argtypes = cached_result.argtypes
        cache_overridden_by_const = cached_result.overridden_by_const
        for i in 1:nargs
            if !is_argtype_match(given_argtypes[i],
                                 cache_argtypes[i],
                                 CC.getindex(cache_overridden_by_const, i))
                cache_match = false
                break
            end
        end
        if method.isva && cache_match
            cache_match = is_argtype_match(tuple_tfunc(given_argtypes[(nargs + 1):end]),
                                           cache_argtypes[end],
                                           CC.getindex(cache_overridden_by_const, CC.length(cache_overridden_by_const)))
        end
        cache_match || continue
        return cached_result
    end
    return nothing
end

CC.push!(cache::JETLocalCache, inf_result::InferenceResult) = CC.push!(cache.cache, inf_result)


text:
function random_walk(steps::Int64)::Array{Float64, 1}
    SQRT_STEPS = sqrt(steps)
    walk = ones(steps)
    for idx in 2:steps
        if rand() < .5
            y = -1
        else
            y = 1
        end
        walk[idx] = walk[idx-1] + (y / SQRT_STEPS)
    end
    walk
end

random_walk(100)

multiple = 1
for i in 1:8
    global multiple
    t = @elapsed random_walk(multiple)
    println("jl_random_walk\t$t")
    multiple *= 10
end

ext: jl

max_stars_repo_name: monkeyusage/montecarlo

max_issues_repo_name: monkeyusage/montecarlo

max_forks_repo_name: monkeyusage/montecarlo


text:
# some commentary

struct WeaveFeaturization <: AbstractFeaturization
    atom_features::Vector{<:AbstractAtomFeatureDescriptor}
    pair_features::Vector{<:AbstractPairFeatureDescriptor}
end

function encodable_elements(fzn::WeaveFeaturization)
    # TODO: implement me!
end


text:
function rot_n_stride!(traj::ContactTraj, cache::ContactTraj, stride::Vector{T}, window::Vector{Int}) where T
    rotate!(traj, cache)
    mpc_stride!(traj, stride, window)
    return nothing
end

function rotate!(traj::ContactTraj, cache::ContactTraj)
    H = traj.H

	cache.q[1] .= traj.q[1]
    cache.u[1] .= traj.u[1]
    cache.w[1] .= traj.w[1]
    cache.γ[1] .= traj.γ[1]
    cache.b[1] .= traj.b[1]
    cache.z[1] .= traj.z[1]
    cache.θ[1] .= traj.θ[1]

    for t = 1:H+1
        traj.q[t] .= traj.q[t+1]
    end

    for t = 1:H-1
        traj.u[t] .= traj.u[t+1]
        traj.w[t] .= traj.w[t+1]
        traj.γ[t] .= traj.γ[t+1]
        traj.b[t] .= traj.b[t+1]
        traj.z[t] .= traj.z[t+1]
        traj.θ[t] .= traj.θ[t+1]
    end

    traj.q[end] .= cache.q[1]
    traj.u[end] .= cache.u[1]
    traj.w[end] .= cache.w[1]
    traj.γ[end] .= cache.γ[1]
    traj.b[end] .= cache.b[1]
    traj.z[end] .= cache.z[1]
    traj.θ[end] .= cache.θ[1]

    return nothing
end

function set_window!(traj::ContactTraj, ref::ContactTraj, window)

	for (i, t) in enumerate(window)
		traj.q[i] .= ref.q[t]
	end

	for (i, t) in enumerate(window[1:end-2])
		traj.u[i] .= ref.u[t]
		traj.w[i] .= ref.w[t]
		traj.γ[i] .= ref.γ[t]
		traj.b[i] .= ref.b[t]
		traj.z[i] .= ref.z[t]
		traj.θ[i] .= ref.θ[t]
	end

    return nothing
end

function set_trajectory!(traj::ContactTraj, ref_traj::ContactTraj)
	H = ref_traj.H
	for t = 1:H
		traj.q[t] .= ref_traj.q[t]
		traj.u[t] .= ref_traj.u[t]
		traj.w[t] .= ref_traj.w[t]
		traj.γ[t] .= ref_traj.γ[t]
		traj.b[t] .= ref_traj.b[t]
		traj.z[t] .= ref_traj.z[t]
		traj.θ[t] .= ref_traj.θ[t]
	end
	traj.q[H+1] .= ref_traj.q[H+1]
	traj.q[H+2] .= ref_traj.q[H+2]
	return nothing
end

"""
    Update the last two cofiguaraton to be equal to the fisrt two up to an constant offset.
"""
function mpc_stride!(traj::ContactTraj, stride::Vector{T}, window::Vector{Int}) where T
    H = traj.H

    for t = H+1:H+2
		traj.q[t] .= traj.q[t-H]
		traj.q[t] .+= stride

		# update_z!(traj, t - 2)
		τ = t - 2
		traj.z[τ][traj.iq2] = traj.q[τ+2]

		# update_θ!(traj, t - 2)
		traj.θ[τ][traj.iq0] = traj.q[τ]
		traj.θ[τ][traj.iq1] = traj.q[τ+1]

		# update_θ!(traj, t - 3)
		τ = t - 3
		traj.θ[τ][traj.iq0] = traj.q[τ]
		traj.θ[τ][traj.iq1] = traj.q[τ+1]
    end

    return nothing
end

function get_stride(model::Model, traj::ContactTraj) #TODO: dispatch over environment / model
    stride = zeros(model.nq)
    stride[1] = traj.q[end-1][1] - traj.q[1][1]
    return stride
end

function update_altitude!(alt::Vector{T}, ϕ::Vector{T}, s::Simulation{T}, traj::Trajectory{T}, t::Int, nc::Int, N_sample::Int;
	threshold = 1.0, verbose = false) where T

	idx1 = max(0, t - 1 - N_sample) + 1

	for i = 1:nc
		γ_max = 0.0
		idx_max = 0

		for j = idx1:t-1
			if traj.γ[j][i] > γ_max
				γ_max = traj.γ[j][i]
				idx_max = j
			end
		end

		if γ_max > threshold
			s.ϕ(ϕ, traj.q[idx_max+2])
			alt[i] = ϕ[i]
			verbose && println(" ")
			verbose && println("point $i in contact")
			verbose && println("sim_step : $idx_max")
			verbose && println("alt      : $ϕ")
			verbose && println("force    : $(traj.γ[idx_max][i])")
		end
	end
end

function update_altitude!(alt::Vector{T}, ϕ::Vector{T}, s::Simulation{T}, x::Vector{T}, nq::Int, nc::Int;
	threshold = 1.0,
	verbose = false) where T

	s.ϕ(ϕ, x[1:2nq])
	γ = x[2nq .+ (1:nc)]

	for i = 1:nc
		if γ[i] > threshold
			alt[i] = ϕ[i]
			verbose && println(" ")
			verbose && println("point $i in contact")
			verbose && println("sim_step : $idx_max")
			verbose && println("alt      : $(ϕ[i])")
			verbose && println("force    : $(γ[i])")
		end
	end
end

function live_plotting(model::Model{T}, ref_traj::ContactTraj,
		sim_traj::Trajectory, newton::Newton, q0::AbstractVector{T},
		q1::AbstractVector{T}, t::Int) where T
	nq = model.nq
	nu = model.nu

	ql = 1
	qu = nq
	ul = 1
	uu = nu
	plt = plot(layout=grid(2,1,heights=[0.7, 0.3], figsize=[(1000, 1000),(400,400)]), legend=:bottomright, xlims=(0,20))
	plot!(plt[1,1], hcat([Vector(x[ql:qu]) for x in newton.traj.q]...)', color=:blue, linewidth=1.0, label=nothing)
	plot!(plt[1,1], hcat([Vector(x[ql:qu]) for x in ref_traj.q]...)', linestyle=:dot, color=:red, linewidth=3.0, label=nothing)

	scatter!((2-1/N_sample)ones(nq), sim_traj.q[t+0], markersize=8.0, color=:lightgreen, label="simul. q0 q1")
	scatter!(2*ones(nq), sim_traj.q[t+1], markersize=8.0, color=:lightgreen, label=nothing)

	scatter!(plt[1,1], 1*ones(nq), q0, markersize=6.0, color=:blue, label="newton q0 q1")
	scatter!(plt[1,1], 2*ones(nq), q1, markersize=6.0, color=:blue, label=nothing)

	scatter!(plt[1,1], 1*ones(nq), ref_traj.q[1], markersize=4.0, color=:red, label="refer. q0 q1 q2")
	scatter!(plt[1,1], 2*ones(nq), ref_traj.q[2], markersize=4.0, color=:red, label=nothing)
	scatter!(plt[1,1], 3*ones(nq), ref_traj.q[3], markersize=4.0, color=:red, label=nothing)

	plot!(plt[2,1], hcat([Vector(x[ul:uu]) for x in newton.traj.u]...)', color=:blue, linewidth=1.0, label=nothing)#"u newton")
	plot!(plt[2,1], hcat([Vector(x[ul:uu]) for x in ref_traj.u[1:end]]...)', linestyle=:dot, color=:red, linewidth=3.0, label=nothing)#"u refer.")
	display(plt)
end


text:
module LibSpatialIndex

    include("capi.jl")

    version() = unsafe_string(C.SIDX_Version())

    function _checkresult(result::C.RTError, message::String)
        if result == C.RT_Failure || result == C.RT_Fatal
            error(message)
        end
    end

    """
    The RTree index [guttman84] is a balanced tree structure that consists of
    index nodes, leaf nodes and data.

    Every node (leaf and index) has a fixed capacity of entries, (the node
    capacity) chosen at index creation An `RTree` abstracts the data with their
    Minimum Bounding Region (`MBR`) and clusters these MBRs according to various
    heuristics in the leaf nodes. Queries are evaluated from the root of the
    tree down the leaves. Since the index is balanced nodes can be under full.
    They cannot be empty though. A `fillfactor` specifies the minimum number of
    entries allowed in any node. The fill factor is usually close to `70%`.

    # Options

    * `ndim`: Dimensionality of the data that will be inserted.
    * `indextype`: one of `RT_RTree` (default), `RT_MVRTree`, or `RT_TPRTree`.
    * `variant`: one of `RT_Linear`, `RT_Quadratic`, or `RT_Star` (default).
    * `storage`: one of `RT_Memory` (default), `RT_Disk`, or `RT_Custom`.
    * `indexcapacity`: The index node capacity. Default is `100`.
    * `leafcapacity`: The leaf node capacity. Default is `100`.
    * `leafpoolcapacity`: Default is `100`.
    * `indexpoolcapacity`: Default is `100`.
    * `regionpoolcapacity`: Default is `1000`.
    * `pointpoolcapacity`: Default is `500`.
    * `tightMBR`: Default is `true`.
    * `nearminimumoverlapfactor`: Default is `32`.
    * `fillfactor`: The fill factor. Default is `0.7`.
    * `splitdistributionfactor`: Default is `0.4`.
    * `reinsertfactor`: Default is `0.3`.

    # Performance

    Dataset size, data density, etc. have nothing to do with capacity and page
    size. What you are trying to achieve is fast reads from the disk and take
    advantage of disk buffering and prefetching.

    Normally, you select the page size to be equal to the disk page size
    (depends on how you format the drive). Then you choose the node capacity to
    be enough to fill the whole page (including data entries if you have any).

    There is a tradeoff though in making node capacity too large. The larger the
    capacity, the longer the “for loops” for inserting, deleting, locating node
    entries take (CPU time). On the other hand, the larger the capacity the
    shorter the tree becomes, reducing the number of random IOs to reach the
    leaves. Hence, you might want to fit multiple nodes (of smaller capacity)
    inside a single page to balance I/O and CPU cost, in case your disk page
    size is too large.

    Finally, if you have enough buffer space to fit most of the index nodes in
    main memory, then large capacities do not make too much sense, because the
    height of the tree does not matter any more. Of course, the smaller the
    capacity, the larger the number of leaf nodes you will have to retrieve from
    disk for range queries (point queries and nearest neighbor queries will not
    suffer that much). So very small capacities hurt as well.

    There is another issue when trying to fit multiple nodes per page: Leftover
    space. You might have leftover space due to data entries that do not have a
    fixed size. Unfortunately, in that case, leftover space per page is lost,
    negatively impacting space usage.

    Selecting the right page size is easy; make it equal to the disk page size.
    Selecting the right node capacity is an art.

    # References
    [guttman84] “R-Trees: A Dynamic Index Structure for Spatial Searching”
        Antonin Guttman, Proc. 1984 ACM-SIGMOD Conference on Management of
        Data (1985), 47-57.
    """
    mutable struct RTree
        index::C.IndexH
        indextype::C.RTIndexType
        ndim::UInt32
        variant::C.RTIndexVariant
        storage::C.RTStorageType
        indexcapacity::UInt32
        leafcapacity::UInt32
        leafpoolcapacity::UInt32
        indexpoolcapacity::UInt32
        regionpoolcapacity::UInt32
        pointpoolcapacity::UInt32
        tightMBR::Bool
        nearminimumoverlapfactor::UInt32
        fillfactor::Float64
        splitdistributionfactor::Float64
        reinsertfactor::Float64

        function RTree(ndim::Integer;
                indextype::C.RTIndexType = C.RT_RTree,
                variant::C.RTIndexVariant = C.RT_Star,
                storage::C.RTStorageType = C.RT_Memory,
                indexcapacity::Integer = 100,
                leafcapacity::Integer = 100,
                leafpoolcapacity::Integer = 100,
                indexpoolcapacity::Integer = 100,
                regionpoolcapacity::Integer = 1000,
                pointpoolcapacity::Integer = 500,
                tightMBR::Bool = true,
                nearminimumoverlapfactor::Integer = 32,
                fillfactor::Real = 0.7,
                splitdistributionfactor::Real = 0.4,
                reinsertfactor::Real = 0.3
            )
            p = C.IndexProperty_Create()
            C.IndexProperty_SetIndexType(p, indextype)
            C.IndexProperty_SetDimension(p, UInt32(ndim))
            C.IndexProperty_SetIndexVariant(p, variant)
            C.IndexProperty_SetIndexStorage(p, storage)
            C.IndexProperty_SetIndexCapacity(p, UInt32(indexcapacity))
            C.IndexProperty_SetLeafCapacity(p, UInt32(leafcapacity))
            C.IndexProperty_SetLeafPoolCapacity(p, UInt32(leafpoolcapacity))
            C.IndexProperty_SetIndexPoolCapacity(p, UInt32(indexpoolcapacity))
            C.IndexProperty_SetRegionPoolCapacity(p, UInt32(regionpoolcapacity))
            C.IndexProperty_SetPointPoolCapacity(p, UInt32(pointpoolcapacity))
            # IndexProperty_SetBufferingCapacity(iprop::IndexPropertyH, value::UInt32)
            C.IndexProperty_SetEnsureTightMBRs(p, UInt32(tightMBR))
            # IndexProperty_SetOverwrite(iprop::IndexPropertyH, value::UInt32)
            C.IndexProperty_SetNearMinimumOverlapFactor(p, UInt32(nearminimumoverlapfactor))
            # IndexProperty_SetWriteThrough(iprop::IndexPropertyH, value::UInt32)
            C.IndexProperty_SetFillFactor(p, fillfactor)
            C.IndexProperty_SetSplitDistributionFactor(p, splitdistributionfactor)
            # IndexProperty_SetTPRHorizon(iprop::IndexPropertyH, value::Cdouble)
            C.IndexProperty_SetReinsertFactor(p, reinsertfactor)
            # IndexProperty_SetFileName(iprop::IndexPropertyH, value)
            # IndexProperty_SetFileNameExtensionDat(iprop::IndexPropertyH, value)
            # IndexProperty_SetFileNameExtensionIdx(iprop::IndexPropertyH, value)
            # IndexProperty_SetCustomStorageCallbacksSize(iprop::IndexPropertyH, value::UInt32)
            # IndexProperty_SetCustomStorageCallbacks(iprop::IndexPropertyH, value)
            # IndexProperty_SetIndexID(iprop::IndexPropertyH, value::Int64)
            index = C.Index_Create(p)
            Bool(C.Index_IsValid(index)) || error("Invalid index")
            C.IndexProperty_Destroy(p)
            rtree = new(index, indextype, ndim, variant, storage, indexcapacity,
                leafcapacity, leafpoolcapacity, indexpoolcapacity,
                regionpoolcapacity, pointpoolcapacity, tightMBR,
                nearminimumoverlapfactor, fillfactor, splitdistributionfactor,
                reinsertfactor)
            finalizer( x -> (C.Index_Destroy(x.index); x.index = C_NULL), rtree)
            rtree
        end
    end

    """
    Inserts an item into the `rtree` with given `id` and boundingbox specified
    by `minvalues` and `maxvalues`, where the item lies within the interval
    `[minvalues[i],maxvalues[i]]` for each axis `i` in 1, ..., `ndim`.
    """
    function insert!(
            rtree::RTree,
            id::Integer,
            minvalues::Vector{Float64},
            maxvalues::Vector{Float64}
        )
        C.Index_InsertData(rtree.index, Int64(id), pointer(minvalues),
            pointer(maxvalues), UInt32(length(minvalues)), Ptr{UInt8}(0), Cint(0)
        )
    end

    """
    Returns a vector of `id`s corresponding to items in `rtree` that intersects
    the box specified by `minvalues` and `maxvalues`.

    Each item intersects the interval `[minvalues[i],maxvalues[i]]` for each
    axis `i` in 1, ..., `ndim`.
    """
    function intersects(
            rtree::RTree,
            minvalues::Vector{Float64},
            maxvalues::Vector{Float64}
        )
        items = Ref{Ptr{Int64}}()
        nresults = Ref{UInt64}()
        result = C.Index_Intersects_id(rtree.index, pointer(minvalues),
            pointer(maxvalues), UInt32(length(minvalues)), items, nresults
        )
        _checkresult(result, "Index_Intersects_id: Failed to evaluate")
        unsafe_wrap(Array, items[], nresults[])
    end

    """
    Returns a vector of `id`s corresponding to items in `rtree` that intersects
    the coordinates specified by `point`.
    """
    intersects(rtree::RTree, point::Vector{Float64}) = intersects(rtree, point, point)

    """
    Returns a vector of `id`s corresponding to the `k` items in `rtree`
    that are nearest to the box specified by `minvalues` and `maxvalues`.

    Each item intersects the interval `[minvalues[i],maxvalues[i]]` for each
    axis `i` in 1, ..., `ndim`. If there are fewer than `k` items in `rtree`,
    it will return less than `k` items. On the other hand, if there are ties
    between some of the items, it might return more than `k` items.
    """
    function knn(
            rtree::RTree,
            minvalues::Vector{Float64},
            maxvalues::Vector{Float64},
            k::Integer
        )
        items = Ref{Ptr{Int64}}()
        nresults = Ref{UInt64}(k)
        result = C.Index_NearestNeighbors_id(rtree.index, pointer(minvalues),
            pointer(maxvalues), UInt32(length(minvalues)), items, nresults)
        _checkresult(result, "Index_NearestNeighbors_id: Failed to evaluate")
        unsafe_wrap(Array, items[], nresults[])
    end

    """
    Returns a vector of `id`s corresponding to the `k` items in `rtree`
    that are nearest to the box specified by `minvalues` and `maxvalues`.

    If there are fewer than `k` items in `rtree`, it will return less than `k`
    items. On the other hand, if there are ties between some of the items,
    it might return more than `k` items.
    """
    knn(rtree::RTree, point::Vector{Float64}, k::Integer) = knn(rtree, point, point, k)

end # module


text:
# Copyright (c) 2019 Arpit Bhatia and contributors                               #src
#                                                                                #src
# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
# of this software and associated documentation files (the "Software"), to deal  #src
# in the Software without restriction, including without limitation the rights   #src
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
# copies of the Software, and to permit persons to whom the Software is          #src
# furnished to do so, subject to the following conditions:                       #src
#                                                                                #src
# The above copyright notice and this permission notice shall be included in all #src
# copies or substantial portions of the Software.                                #src
#                                                                                #src
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
# SOFTWARE.                                                                      #src

# # Solvers and Solutions

# **Originally Contributed by**: Arpit Bhatia

# The purpose of this part of the tutorial is to introduce you to solvers and
# how to use them with JuMP. We'll also learn what to do with a problem after
# the solver has finished optimizing it.

# ## What is a Solver?

# A solver is a software package that incorporates algorithms for finding
# solutions to one or more classes of problem. For example, GLPK, which we used
# in the previous tutorials is a solver for linear programming (LP) and mixed
# integer programming (MIP) problems. It incorporates algorithms such as the
# simplex method, interior-point method etc. JuMP currently supports a number of
# open-source and commercial solvers which can be viewed in the [Supported-solvers](@ref)
# table.

# ## What is MathOptInterface?

# Each mathematical optimization solver API has its own concepts and data
# structures for representing optimization models and obtaining results.
# However, it is often desirable to represent an instance of an optimization
# problem at a higher level so that it is easy to try using different solvers.

# MathOptInterface (MOI) is an abstraction layer designed to provide an
# interface to mathematical optimization solvers so that users do not need to
# understand multiple solver-specific APIs. MOI can be used directly, or through
# a higher-level modeling interface like JuMP.

# Note that JuMP re-exports MathOptInterface and you can use the shortcut `MOI`
# to refer to MathOptInterface in your code.

# ## Constructing a model

# JuMP models can be created in three different modes: `AUTOMATIC`, `MANUAL` and
# `DIRECT`. We'll use the following LP to illustrate them.

# ```math
# \begin{aligned}
# & \max_{x,y} & x + 2y \\
# & \;\;\text{s.t.} & x + y &\leq 1 \\
# & & 0\leq x, y &\leq 1 \\
# \end{aligned}
# ```

using JuMP
using GLPK

# ### `AUTOMATIC` Mode

# #### With Optimizer

# This is the easiest method to use a solver in JuMP. In order to do so, we
# simply set the solver inside the Model constructor.

model_auto = Model(GLPK.Optimizer)
@variable(model_auto, 0 <= x <= 1)
@variable(model_auto, 0 <= y <= 1)
@constraint(model_auto, x + y <= 1)
@objective(model_auto, Max, x + 2y)
optimize!(model_auto)
objective_value(model_auto)

# #### No Optimizer (at first)

# It is also possible to create a JuMP model with no optimizer attached. After
# the model object is initialized empty and all its variables, constraints and
# objective are set, then we can attach the solver at `optimize!` time.

model_auto_no = Model()
@variable(model_auto_no, 0 <= x <= 1)
@variable(model_auto_no, 0 <= y <= 1)
@constraint(model_auto_no, x + y <= 1)
@objective(model_auto_no, Max, x + 2y)
set_optimizer(model_auto_no, GLPK.Optimizer)
optimize!(model_auto_no)
objective_value(model_auto_no)

# Note that we can also enforce the automatic mode by passing
# `caching_mode = MOIU.AUTOMATIC` in the Model function call.

# ### `MANUAL` Mode

# This mode is similar to the `AUTOMATIC` mode, but there are less protections
# from the user getting errors from the solver API. On the other side, nothing
# happens silently, which might give the user more control. It requires
# attaching the solver before the solve step using the `MOIU.attach_optimizer()`
# function.

model_manual = Model(GLPK.Optimizer, caching_mode = MOIU.MANUAL)
@variable(model_manual, 0 <= x <= 1)
@variable(model_manual, 0 <= y <= 1)
@constraint(model_manual, x + y <= 1)
@objective(model_manual, Max, x + 2y)
MOIU.attach_optimizer(model_manual)
optimize!(model_manual)
objective_value(model_manual)

# ### `DIRECT` Mode

# Some solvers are able to handle the problem data directly. This is common for
# LP/MIP solver but not very common for open-source conic solvers. In this case
# we do not set a optimizer, we set a backend which is more generic and is able
# to hold data and not only solving a model.

model_direct = direct_model(GLPK.Optimizer())
@variable(model_direct, 0 <= x <= 1)
@variable(model_direct, 0 <= y <= 1)
@constraint(model_direct, x + y <= 1)
@objective(model_direct, Max, x + 2y)
optimize!(model_direct)
objective_value(model_direct)

# ### Solver Options

# Many of the solvers also allow options to be passed in. However, these options
# are solver-specific. To find out the various options available, please check
# out the individual solver packages. Some examples for the GLPK solver are
# given below.

using GLPK

# To turn off printing (i.e. silence the solver),

model = Model(optimizer_with_attributes(GLPK.Optimizer, "msg_lev" => 0));

# To increase the maximum number of simplex iterations:

model = Model(optimizer_with_attributes(GLPK.Optimizer, "it_lim" => 10_000));

# To set the solution timeout limit (in milliseconds):

model = Model(optimizer_with_attributes(GLPK.Optimizer, "tm_lim" => 5_000));

# ## How to querying the solution

# So far we have seen all the elements and constructs related to writing a JuMP
# optimization model. In this section we reach the point of what to do with a
# solved problem. JuMP follows closely the concepts defined in MathOptInterface
# to answer user questions about a finished call to `optimize!(model)`. The
# three main steps in querying a solution are given below. We'll use the model
# we created in `AUTOMATIC` mode with an optimizer attached in this section.

# ### The termination status

# Termination statuses are meant to explain the reason why the optimizer stopped
# executing in the most recent call to `optimize!`.

termination_status(model_auto)

# You can view the different termination status codes by referring to the docs
# or though checking the possible types using the below command.

display(typeof(MOI.OPTIMAL))

# ### The primal and dual status

# These statuses indicate what kind of result is available to be queried from
# the model. It's possible that no result is available to be queried. We shall
# discuss more on the dual status and solutions in the Duality tutorial.

primal_status(model_auto)

#-

dual_status(model_auto)

# As we saw before, the result (solution) status codes can be viewed directly
# from Julia.

display(typeof(MOI.FEASIBLE_POINT))

# ### Getting the primal solution

# Provided the primal status is not `MOI.NO_SOLUTION`, we can inspect the
# solution values and optimal cost.

value(x)

#-

value(y)

#-

objective_value(model_auto)

# Since it is possible that no solution is available to be queried from the
# model, calls to [`value`](@ref) may throw errors. Hence, it is recommended to
# check for the presence of solutions.

model_no_solution = Model(GLPK.Optimizer)
@variable(model_no_solution, 0 <= x <= 1)
@variable(model_no_solution, 0 <= y <= 1)
@constraint(model_no_solution, x + y >= 3)
@objective(model_no_solution, Max, x + 2y)

optimize!(model_no_solution)

try #hide
if termination_status(model_no_solution) == MOI.OPTIMAL
    optimal_solution = value(x)
    optimal_objective = objective_value(model_no_solution)
elseif termination_status(model_no_solution) == MOI.TIME_LIMIT && has_values(model_no_solution)
    suboptimal_solution = value(x)
    suboptimal_objective = objective_value(model_no_solution)
else
    error("The model was not solved correctly.")
end
catch err; showerror(stderr, err); end  #hide

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ..Helper: init

function simulate!(x, m, y)
    dim, steps = size(x)
    for t = 2:steps
        A_mul_B!(y, m, view(x, :, t-1))

        # x[:, t] .= tanh.(y)
        for i = 1:dim
            x[i, t] = tanh(y[i])
        end
    end
end

function prepare(name; m=nothing)
    res = init(name)
    if m != nothing
        res[:m] = m
    end
    function run()
        simulate!(res[:x], res[:m], res[:y])
        res
    end
end

max_stars_repo_name: tkf/comparatist

max_issues_repo_name: tkf/comparatist

max_forks_repo_name: tkf/comparatist


text:
# ---
# title: 706. Design HashMap
# id: problem706
# author: Qling
# date: 2021-03-14
# difficulty: Easy
# categories: Hash Table, Design
# link: <https://leetcode.com/problems/design-hashmap/description/>
# hidden: true
# ---
#
# Design a HashMap without using any built-in hash table libraries.
#
# To be specific, your design should include these functions:
#
#   * `put(key, value)` : Insert a (key, value) pair into the HashMap. If the value already exists in the HashMap, update the value.
#   * `get(key)`: Returns the value to which the specified key is mapped, or -1 if this map contains no mapping for the key.
#   * `remove(key)` : Remove the mapping for the value key if this map contains the mapping for the key.
#
#
# **Example:**
#
#
#
#     MyHashMap hashMap = new MyHashMap();
#     hashMap.put(1, 1);
#     hashMap.put(2, 2);
#     hashMap.get(1);            // returns 1
#     hashMap.get(3);            // returns -1 (not found)
#     hashMap.put(2, 1);          // update the existing value
#     hashMap.get(2);            // returns 1
#     hashMap.remove(2);          // remove the mapping for 2
#     hashMap.get(2);            // returns -1 (not found)
#
#
#
# **Note:**
#
#   * All keys and values will be in the range of `[0, 1000000]`.
#   * The number of operations will be in the range of `[1, 10000]`.
#   * Please do not use the built-in HashMap library.
#
#
## @lc code=start
using LeetCode

struct MyHashMap
    buckets::Int
    table::Vector{Vector{T}} where {T}
    MyHashMap() = new(1009, [Tuple{Int,Int}[] for _ = 1:1009])
end

hash_706(hash_map::MyHashMap, key::Int) = mod1(key, hash_map.buckets)

function put_706!(hash_map::MyHashMap, key::Int, value::Int)
    hash_key = hash_706(hash_map, key)

    if !isempty(hash_map.table[hash_key])
        for (i, item) in enumerate(hash_map.table[hash_key])
            if item[1] == key
                hash_map.table[hash_key][i] = (key, value)
                return
            end
        end
    end

    push!(hash_map.table[hash_key], (key, value))
    return
end

function remove_706!(hash_map::MyHashMap, key)
    hash_key = hash_706(hash_map, key)

    if !isempty(hash_map.table[hash_key])
        for (i, item) in enumerate(hash_map.table[hash_key])
            if item[1] == key
                deleteat!(hash_map.table[hash_key], i)
                return
            end
        end
    end
end

function get_706(hash_map::MyHashMap, key)
    hash_key = hash_706(hash_map, key)

    for item in hash_map.table[hash_key]
        !isempty(item) && item[1] == key && return item[2]
    end

    return -1
end
## @lc code=end


text:
# Instead of using cudaStreamSynchronize, this implements
# a version that waits on a condition. This allows you to
# use this with Julia Task interface. A practical application
# is that one Julia process can feed multiple devices, with the jobs
# running asynchronously.

immutable StreamHandle
    handle::Ptr{Void}
end

convert(::Type{Ptr{Void}}, shnd::StreamHandle) = shnd.handle

abstract AbstractStream

type Stream <: AbstractStream
    handle::StreamHandle
    c::Condition
end
function Stream()
    p = Array(Ptr{Void}, 1)
    rt.cudaStreamCreate(p)
    hnd = StreamHandle(p[1])
    Stream(hnd, Condition())
end
NullStream() = Stream(StreamHandle(C_NULL), Condition())

convert(::Type{Ptr{Void}}, s::Stream) = convert(Ptr{Void}, s.handle)

const null_stream = NullStream()

destroy(s::Stream) = rt.cudaStreamDestroy(s.ptr)

synchronize(s::Stream) = rt.cudaStreamSynchronize(s)

function wait(s::Stream)
    runnotify = data -> notify(s.c)
    notifyasync = Base.SingleAsyncWork(runnotify)
    rt.cudaStreamAddCallback(s, c_async_send_cudastream, notifyasync.handle, 0)
    wait(s.c)
    rt.checkerror(rt.cudaGetLastError())
end

async_send_cudastream(hnd, status, data) = (ccall(:uv_async_send,Cint,(Ptr{Void},),data); nothing)
const c_async_send_cudastream = cfunction(async_send_cudastream, Void, (rt.cudaStream_t, rt.cudaError_t, Ptr{Void}))


text:
using LinearAlgebra ## "norm"
using StatsBase ## "sample"
using Statistics ## "mean"

## Some parts have been vectorized for better performance.
"""
**BFO** is a function that computes the **classical (non-adaptive) BFO (Bacterial Foraging Optimization)**

**Inputs**:
1. **J**: a function with domain R^n
2. **Range**: exploration range: Range^n

3. **S** = 10, number of bacteria
4. **Sr** = 4,  number of bacteria removed in reproductive step
5. **Nc** = 20, number of chemotactic steps
6. **Ns** = 5, number of swim steps
7. **Nre** = 5, number of reproductive steps
8. **Ned** = 2, elimination and dispersal steps
9. **Ped** = 0.3, probability of elimination
10. **Ci** = (mean(Range[:,2].-Range[:,1])/S), the run-length unit
**Output**: a dictionary that stores
1. the minimum value of J
2. the point achieving this minimu value
3. the path of each bacterium (for plotting illustration)
"""
function BFO(J, Range::Array{Float64,2}; S = 10::Int, Sr = 4::Int, Nc = 20::Int, Ns = 5::Int,
        Nre = 5::Int, Ned = 2::Int, Ped = 0.3::Float64, Ci = (mean(Range[:,2].-Range[:,1])/S)::Float64)

    n = size(Range,1) ## dimension of the input of J
    ## randomly generate S bacteria in Range
    B_loc = RandPts(Range,S)
    ## a dictionary recording the path of bacterium i
    Path_Dict = Dict{Int,Array{Float64,2}}(i=>[zeros(n,0) B_loc[:,i]] for i=1:S)
    for l = 1:Ned ## index of elimination-dispersal steps
        for k = 1:Nre ## index of reproductive steps
            for j = 1:Nc ## index of chemotactic steps
                ## Chemotactic Step; variable used: Range, S, B_loc, Ci, J
                println("Chemotactic step $j")
                @time Chemo = ChemotacticStep(J, Range::Array{Float64,2}, S::Int, Ns::Int, Ci::Float64, B_loc::Array{Float64,2})
                B_loc = Chemo[1]
                for i=1:S
                    Path_Dict[i] = [Path_Dict[i] Chemo[2][i]]
                end
            end ## end of chemotactic steps

            if k<Nre
                ## Reproductive Step; variable used: B_loc, Sr, S, J
                println("Reproduction step $k")
                @time Repro = ReproductionStep(J, B_loc::Array{Float64,2}, Sr::Int, S::Int, n::Int)
                B_loc = Repro[1]
                for i=1:S
                    Path_Dict[i] = [Path_Dict[i] Repro[2][i]]
                end
            end
        end
        ## Elimination-Dispersal Step; variable used: Ped, Range, B_loc
        println("Elimination-Dispersal step $l")
        if l<Ned
            @time EliDis = ElimDispStep(B_loc, Ped, Range, S, n)
            B_loc = EliDis[1]
            for i=1:S
                Path_Dict[i] = [Path_Dict[i] EliDis[2][i]]
            end
        end
    end
    B_best = Int(sortslices([[J(B_loc[:,i]) for i=1:S] collect(1:S)], dims = 1)[1,2]) ## best bacterium
    X_best = B_loc[:,B_best] ## best location for minimizing J
    J_best = J(X_best) ## best (minimum) J value
    return Dict("Minimum"=>J_best, "Minimum Point"=>X_best, "Path_Dict"=>Path_Dict)
end

"""
**ChemotacticStep** is a function performing the chemotactic step in BFO.
"""
function ChemotacticStep(J, Range::Array{Float64,2}, S::Int, Ns::Int, Ci::Float64, B_loc::Array{Float64,2})
    Path = Dict{Int, Array{Float64,2}}(i=>zeros(size(Range,1),2) for i=1:S)
    ToSwim = collect(1:S) ## monitor the swimming bacteria
    Tumble = RandUnit(Range,S)
    m = 0 ## index of swimming
    while (length(ToSwim)!=0)&(m<Ns) ## tumble/swim
        m = m + 1
        J_old = [J(B_loc[:,i]) for i in ToSwim]

        B_loc_new = copy(B_loc)
        B_loc_new[:,ToSwim] = B_loc[:,ToSwim].+Ci*Tumble[:,ToSwim]

        ## Mirror back the out-of-range bacteria
        ## Range_tail = Range[:,2]
        ## Range_head = Range[:,1]
        ## 1. mirror back bacteria tumbling out from below
        B_loc_new[:,ToSwim] = (abs.(B_loc_new[:,ToSwim].-Range[:,1])).+Range[:,1]
        ## 2. mirror back bacteria tumbling out from above
        B_loc_new[:,ToSwim] = Range[:,2].-(abs.(Range[:,2].-B_loc_new[:,ToSwim]))

        ## Evaluate J at new locations
        J_new = [J(B_loc_new[:,i]) for i in ToSwim]

        ## Find out improved bacteria and update the swimming ones
        ToSwim = ToSwim[findall(J_new.<J_old)]

        for i in ToSwim
            Path[i] = [Path[i] B_loc_new[:,i]]
        end
        B_loc[:,ToSwim] = B_loc_new[:,ToSwim]
    end ## end of tumble/swim
    return (B_loc,Path)
end

"""
**ReproductionStep** is a function performing the reproduction step in BFO.
"""
function ReproductionStep(J, B_loc::Array{Float64,2}, Sr::Int, S::Int, n::Int)
    Path = Dict{Int,Array{Float64,2}}(i=>zeros(n,2) for i=1:S)
    ## (!!!) I define the health of a bacterium as the J value of its current location
    Health = [J(B_loc[:,i]) for i=1:S]
    Health_sort = sortslices([Health collect(1:S)], dims = 1) ## sort the bacteria according to health
    B_survive = Array{Int,1}(Health_sort[1:S-Sr,2]) ## ## pick out the most healthy (S-Sr) bacteria
    B_die = setdiff(collect(1:S), B_survive)
    B_loc[:,B_die] = B_loc[:,B_survive[1:Sr]] ## reproduce the most healthy Sr bacteria
    for i in B_die
        Path[i] = [Path[i] B_loc[:,i]]
    end
    return (B_loc, Path)
end

"""
**ElimDispStep** is a function performing the elimination-dispersal step in BFO.
"""
function ElimDispStep(B_loc::Array{Float64,2}, Ped::Float64, Range::Array{Float64,2}, S::Int, n::Int)
    Path = Dict{Int,Array{Float64,2}}(i=>zeros(n,0) for i=1:S)
    Alive = findall([sample([false, true], aweights([Ped,1-Ped])) for i=1:S]) # true = alive, false = kill
    Killed = setdiff(collect(1:S),Alive) ## kiilled bacteria
    ## Rebuild the killed bacteria in random locations in Range
    B_loc[:,Killed] = RandPts(Range,length(Killed))
    ## update the path dictionary
    for i in Killed
        Path[i] = [Path[i] B_loc[:,i]]
    end
    return (B_loc, Path)
end

"""
**RandPts** is a function that randomly uniformly generates m points in dimension size(Range,2) inside Range.

Inputs:

m: number of random points to generate

Range: dim x 2 matrix; the d-th row is the range of the d-th coordinate

Output:

an dim x m matrix; each column is a desired random vector.

"""
function RandPts(Range::Array{Float64,2}, m::Int)
    ## dim = size(Range,1)
    ## randPts = rand(dim,m)
    ## Coeff = (Range[:,2].-Range[:,1])
    ## Intercept = Range[:,1]
    ## randPts = ((Range[:,2].-Range[:,1]).*rand(size(Range,1),m)).+Range[:,1]
    return ((Range[:,2].-Range[:,1]).*rand(size(Range,1),m)).+Range[:,1]
end

"""
**RandUnit** is a function generating m random unit vectors in dimension d. The vector generated is uniform with respect to a Range where random walk is to be applied.

**Inputs**:

1. Range: a d x m matrix; the kth row is the range for the kth variable.
2. m: number of vectors to generate

**Output**:

an d x m matrix whose columns are the desired vectors.

"""
function RandUnit(Range::Array{Float64,2}, m::Int)
    ## d = size(Range,1)
    ## Scales = Range[:,2].-Range[:,1]
    Rand = (rand(size(Range,1),m).-1/2).*(Range[:,2].-Range[:,1])
    for k = 1:m
        Rand[:,k]./=norm(Rand[:,k])
    end
    return Rand
end


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

using DynamicalSystems
using DelayEmbeddings
using DelimitedFiles
using PyPlot

include("../../src/pecuzal_method.jl")
include("../../src/data_analysis_functions.jl")

## This script was for testing reasons only and had no impact of any results
# shown in the paper

## We analyze the reconstruction from standard time delay embedding, the
# MDOP embedding algorithm and the Garcia & Almeida method and compare it to
# our propsed method: The pecuzal method.
# For evaluation of the different reconstructions we consider Uzal's L-statistic,
# the fraction of recurrence rates from the JRP of the reconstruction and the
# reference and only the reference. As a third criterion we look at the
# recurrence time entropy, since it is related to the Kolmogorov-entropy as a
# dynamical invariant. We also look at two other RQA-measures and, moreover, the
# generalized mutual false nearest neighbors.

# Here we look at a chaotic system: The Roessler attractor in the funnel regime

## Integrate system and determine theiler window
using Random
using DifferentialEquations
Random.seed!(234)
# set time interval for integration
N = 15000 # number of samples
transients = 2000
dt_tra = 0.01
tspan = (dt_tra*N) + (dt_tra*transients)

function roessler!(du,u,p,t)
    x,y,z = u
    a,b,c = p
    du[1] = -(y+z) # dx
    du[2] = x + (a*y) # dy
    du[3] = b+z*(x-c) # dz
end

function prob_func_roe(prob,i,repeat)
  @. prob.u0 = [2*rand(); 2*rand(); 2*rand()]
  prob
end

# parameter values
a = .2925
b = .1
c = 8.5
p = [a, b, c]

u0 = [2*rand(); 2*rand(); 2*rand()]
roe = ODEProblem(roessler!,u0,tspan,p)
tr = solve(roe; saveat = dt_tra)
tr = tr.u

roe_ensemble = EnsembleProblem(roe; prob_func = prob_func_roe)
sim = solve(roe_ensemble; saveat = dt_tra, trajectories = 1000)

pygui(true)

figure()
plot3D(tr[:,1], tr[:,2], tr[:,3])

a = .36
b = .01
c = 14
u0 = [rand(); 5+rand(); 50+rand()]
roe = Systems.roessler(u0 ; a=a, b=b, c=c)

# set time interval for integration
N = 10000 # number of samples
dt_tra = 0.01
t = 0:dt_tra:(dt_tra*N)

tr = trajectory(roe, (N*dt_tra); dt = dt_tra, Ttr = (2000*dt_tra))
tr = regularize(tr)

w1 = estimate_delay(tr[:,1], "mi_min")
w2 = estimate_delay(tr[:,2], "mi_min")

w = maximum(hcat(w1,w2))
taus = 0:100

## Perform reconstructions and compute according L-value

# L-value of reference
L_ref = uzal_cost(tr, Tw = (4*w), w = w, samplesize=1)

#Standard TDE
@time Y_tde, _ = standard_embedding_cao(tr[:,2])
L_tde = uzal_cost(Y_tde, Tw = (4*w), w = w, samplesize=1)

#MDOP embedding
println("Computation time MDOP:")
@time begin
    tw1 = mdop_maximum_delay(tr[:,1])
    tw2 = mdop_maximum_delay(tr[:,2])
    lm1 = DelayEmbeddings.findlocalminima(tw1[2])
    if lm1[1] ≤ 2
        lmm1 = try lm1[2]
        catch
            lm1[1]
        end
    else
        lmm1 = lm1[1]
    end
    lm2 = DelayEmbeddings.findlocalminima(tw2[2])
    if lm2[1] ≤ 2
        lmm2 = try lm2[2]
        catch
            lm2[1]
        end
    else
        lmm2 = lm2[1]
    end
    tau_max = maximum(hcat(lmm1, lmm2))
    taus_mdop = 0:tau_max

    Y_mdop, τ_vals_mdop, ts_vals_mdop, FNNs_mdop , βs_mdop = mdop_embedding(tr[:,1:2];
                                                        τs = taus_mdop , w = w)
    L_mdop = uzal_cost(Y_mdop, Tw = (4*w), w = w, samplesize=1)
end

#Garcia & Almeida
println("Computation time Garcia & Almeida method:")
@time Y_GA, τ_vals_GA, ts_vals_GA, FNNs_GA , ns_GA = garcia_almeida_embedding(tr[:,1:2];
                                                    τs = taus , w = w, T = w)
L_GA = uzal_cost(Y_GA, Tw = (4*w), w = w, samplesize=1)

#Pecuzal
println("Computation time pecuzal method:")
@time Y_pec, τ_vals_pec, ts_vals_pec, Ls_pec , εs_pec = pecuzal_embedding_update(tr[:,1:2];
                                                            τs = taus , w = w)
L_pec = minimum(Ls_pec)

@time Y_pec, τ_vals_pec, ts_vals_pec, Ls_pec , εs_pec = pecuzal_embedding_update(tr[:,1:2];
                                                            τs = taus , w = w)

## compute other evaluation statistics

@time mfnn1, mfnn2, mfnn3, mfnn4, f1, f2, f3, f4, RQA_ref, RQA1, RQA2,
                                            RQA3, RQA4, R_ref, R1, R2, R3, R4 =
        perform_recurrence_analysis(tr, Dataset(Y_tde), Dataset(Y_mdop),
                    Dataset(Y_GA), Dataset(Y_pec); ε = 0.08, w = w, kNN = 10)

NN = 1000
figure()
subplot(231)
RR_ref = grayscale(R_ref[1:NN,1:NN])
imshow(RR_ref, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("reference")
subplot(232)
RR1 = grayscale(R1[1:NN,1:NN])
imshow(RR1, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("TDE")
subplot(233)
RR2 = grayscale(R2[1:NN,1:NN])
imshow(RR2, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("MDOP")
subplot(234)
RR3 = grayscale(R3[1:NN,1:NN])
imshow(RR3, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("GA")
subplot(235)
RR4 = grayscale(R4[1:NN,1:NN])
imshow(RR4, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("PECUZAL")

figure()
subplot(131)
RR_ref = grayscale(R_ref[1:NN,1:NN])
imshow(RR_ref, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("Reference")

subplot(132)
RR4 = grayscale(R4[1:NN,1:NN])
imshow(RR4, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("Reconstruction")

subplot(133)
JRP = grayscale(elementwise_product(R4[1:NN,1:NN],R_ref[1:NN,1:NN]))
imshow(JRP, cmap = "binary_r", extent = (1, size(RR_ref)[1], 1, size(RR_ref)[2]))
title("Joint Recurrence Plot (JRP)")

figure()
plot3D(tr[:,1], tr[:,2], tr[:,3])
title("reference")
figure()
plot3D(Y_tde[:,1], Y_tde[:,2], Y_tde[:,3])
title("TDE")
figure()
plot3D(Y_GA[:,1], Y_GA[:,2], Y_GA[:,3])
title("GA")
figure()
plot3D(Y_mdop[:,1], Y_mdop[:,2], Y_mdop[:,3])
title("MDOP")
figure()
plot3D(Y_pec[:,1], Y_pec[:,2], Y_pec[:,3])
title("PECUZAL")

display("mfnn_tde: $mfnn1")
display("mfnn_mdop: $mfnn2")
display("mfnn_GA: $mfnn3")
display("mfnn_pec: $mfnn4")

display("*****")

display("f_tde: $f1")
display("f_mdop: $f2")
display("f_GA: $f3")
display("f_pec: $f4")

display("*****")

display("DET_ref: $(abs(RQA_ref.DET - RQA_ref.DET)/RQA_ref.DET)")
display("DET_tde: $(abs(RQA_ref.DET - RQA1.DET)/RQA_ref.DET)")
display("DET_mdop: $(abs(RQA_ref.DET - RQA2.DET)/RQA_ref.DET)")
display("DET_GA: $(abs(RQA_ref.DET - RQA3.DET)/RQA_ref.DET)")
display("DET_pec: $(abs(RQA_ref.DET - RQA4.DET)/RQA_ref.DET)")

display("*****")

display("RTE_ref: $(abs(RQA_ref.RTE - RQA_ref.RTE)/RQA_ref.RTE)")
display("RTE_tde: $(abs(RQA_ref.RTE - RQA1.RTE)/RQA_ref.RTE)")
display("RTE_mdop: $(abs(RQA_ref.RTE - RQA2.RTE)/RQA_ref.RTE)")
display("RTE_GA: $(abs(RQA_ref.RTE - RQA3.RTE)/RQA_ref.RTE)")
display("RTE_pec: $(abs(RQA_ref.RTE - RQA4.RTE)/RQA_ref.RTE)")

display("*****")

display("Lmax_ref: $(abs(RQA_ref.Lmax - RQA_ref.Lmax)/RQA_ref.Lmax)")
display("Lmax_tde: $(abs(RQA_ref.Lmax - RQA1.Lmax)/RQA_ref.Lmax)")
display("Lmax_mdop: $(abs(RQA_ref.Lmax - RQA2.Lmax)/RQA_ref.Lmax)")
display("Lmax_GA: $(abs(RQA_ref.Lmax - RQA3.Lmax)/RQA_ref.Lmax)")
display("Lmax_pec: $(abs(RQA_ref.Lmax - RQA4.Lmax)/RQA_ref.Lmax)")

display("*****")

display("ENTR_ref: $(abs(RQA_ref.ENTR - RQA_ref.ENTR)/RQA_ref.ENTR)")
display("ENTR_tde: $(abs(RQA_ref.ENTR - RQA1.ENTR)/RQA_ref.ENTR))")
display("ENTR_mdop: $(abs(RQA_ref.ENTR - RQA2.ENTR)/RQA_ref.ENTR))")
display("ENTR_GA: $(abs(RQA_ref.ENTR - RQA3.ENTR)/RQA_ref.ENTR))")
display("ENTR_pec: $(abs(RQA_ref.ENTR - RQA4.ENTR)/RQA_ref.ENTR))")

display("*****")

display("L_ref: $(abs(RQA_ref.L - RQA_ref.L)/RQA_ref.L)")
display("L_tde: $(abs(RQA_ref.L - RQA1.L)/RQA_ref.L)")
display("L_mdop: $(abs(RQA_ref.L - RQA2.L)/RQA_ref.L)")
display("L_GA: $(abs(RQA_ref.L - RQA3.L)/RQA_ref.L)")
display("L_pec: $(abs(RQA_ref.L - RQA4.L)/RQA_ref.L)")

display("*****")

display("MRT_ref: $(abs(RQA_ref.MRT - RQA_ref.MRT)/RQA_ref.MRT)")
display("MRT_tde: $(abs(RQA_ref.MRT - RQA1.MRT)RQA_ref.MRT)")
display("MRT_mdop: $(abs(RQA_ref.MRT - RQA2.MRT)RQA_ref.MRT)")
display("MRT_GA: $(abs(RQA_ref.MRT - RQA3.MRT)RQA_ref.MRT)")
display("MRT_pec: $(abs(RQA_ref.MRT - RQA4.MRT)RQA_ref.MRT)")

display("*****")

display("LAM_ref: $(abs(RQA_ref.LAM - RQA_ref.LAM)/RQA_ref.LAM)")
display("LAM_tde: $(abs(RQA_ref.LAM - RQA1.LAM)RQA_ref.LAM)")
display("LAM_mdop: $(abs(RQA_ref.LAM - RQA2.LAM)RQA_ref.LAM)")
display("LAM_GA: $(abs(RQA_ref.LAM - RQA3.LAM)RQA_ref.LAM)")
display("LAM_pec: $(abs(RQA_ref.LAM - RQA4.LAM)RQA_ref.LAM)")

display("*****")

display("RTE-MRT_ref: $(abs(RQA_ref.RTE/RQA_ref.MRT - RQA_ref.RTE/RQA_ref.MRT)/(RQA_ref.RTE/RQA_ref.MRT))")
display("RTE-MRT_tde: $(abs(RQA_ref.RTE/RQA_ref.MRT - RQA1.RTE/RQA1.MRT)/(RQA_ref.RTE/RQA_ref.MRT))")
display("RTE-MRT_mdop: $(abs(RQA_ref.RTE/RQA_ref.MRT - RQA2.RTE/RQA2.MRT)/(RQA_ref.RTE/RQA_ref.MRT))")
display("RTE-MRT_GA: $(abs(RQA_ref.RTE/RQA_ref.MRT - RQA3.RTE/RQA3.MRT)/(RQA_ref.RTE/RQA_ref.MRT))")
display("RTE-MRT_pec: $(abs(RQA_ref.RTE/RQA_ref.MRT - RQA4.RTE/RQA4.MRT)/(RQA_ref.RTE/RQA_ref.MRT))")

display("*****")


text:
using Test, Julianda

foo = Julianda.Config.loadConfig("../config")

@test length(Julianda.Account.listAccounts(foo)) != 0

bar = Julianda.Account.listAccounts(foo)
foo2 = Julianda.Config.changeAccount(foo, bar[1].id)
@test foo2 == foo
@test Julianda.Config.saveConfig("../config_test", foo)

@test typeof(Julianda.Account.getAccount(foo)) == Julianda.Account.account
@test typeof(Julianda.Account.getAccountSummary(foo)) == Julianda.Account.account
@test length(Julianda.Account.getAccountInstruments(foo)) != 0
@test Julianda.Account.setAccountConfig(foo, "Testing", "0.05")

max_stars_repo_name: CalebDepatie/Julianda

max_forks_repo_name: CalebDepatie/Julianda


text:
# ================================================================
# Reach set interfaces
# ================================================================

"""
    AbstractReachSet{N}

Abstract type for all reach-sets types.

### Notes

A reach-set is a set representation `X` associated to a given time span `Δt`.

In its simplest form, we represent reach-sets with a struct that wraps the tuple
`(X, Δt)`, where `X` corresponds to a geometric set, eg. a polytope, and `Δt` is
the interval with the time span associated to this reach-set.

This type is parametric in `N`. The parameter `N` represents for the type of numerical
coefficient used in the representation (typically, `N = Float64`).

Although concrete subtypes of `AbstractReachSet` may represent the set `X` in
different ways, or carry additional information as an extra type field, they
should all implement the `AbstractReachSet` interface to enable shared functionality
for reach-set types. In particular, each concrete subtype should implement the
following methods:

- `set`    -- return the geometric set
- `setrep` -- return the type of the set representation
- `tspan`  -- return the time interval span
- `tstart` -- return the initial time
- `tend`   -- return the final time
- `dim`    -- return ambient dimension of the reach-set
"""
abstract type AbstractReachSet{N} end

# convenience union for dispatch on structs that behave like a set
const SetOrReachSet = Union{LazySet, UnionSet, UnionSetArray, IA.Interval, IA.IntervalBox, AbstractReachSet}

"""
    basetype(T::Type{<:AbstractReachSet})

Return the base type of the given reach-set type (i.e., without type parameters).

### Input

- `T` -- reach-set type, used for dispatch

### Output

The base type of `T`.
"""
basetype(T::Type{<:AbstractReachSet}) = Base.typename(T).wrapper

"""
    set(R::AbstractReachSet)

Return the geometric set represented by this reach-set.

### Input

- `R` -- reach-set

### Output

The set wrapped by the given reach-set.
"""
function set(::AbstractReachSet) end

# no-op
set(X::LazySet) = X

"""
    setrep(R::AbstractReachSet)

Return the type of the set representation of this reach-set.

### Input

- `R` -- reach-set

### Output

Type of the set representation of the given reach-set.
"""
setrep(::AbstractReachSet)

"""
    tstart(R::AbstractReachSet)

Return the initial time of this reach-set.

### Input

- `R` -- reach-set

### Output

A float representing the initial time of the given reach-set.
"""
function tstart(::AbstractReachSet) end

"""
    tend(R::AbstractReachSet)

Return the final time of this reach-set.

### Input

- `R` -- reach-set

### Output

A float representing the final time of the given reach-set.
"""
function tend(::AbstractReachSet) end

"""
    tspan(R::AbstractReachSet)

Return time span of this reach-set.

### Input

- `R` -- reach-set

### Output

The interval representing the time span of the given reach-set.
"""
function tspan(::AbstractReachSet) end

# if `contiguous` is true, we assume that the time span of each reach-set in this slice
# is sorted, e.g. F[1:4] so it suffices to look at the first and last reach-sets
function tstart(F::VRT; contiguous=false) where {RT<:AbstractReachSet, VRT<:AbstractVector{RT}}
    if contiguous
        t0 = tstart(first(F))
    else
        t0 = minimum(tstart, F)
    end
    return t0
end

# if `contiguous` is true, we assume that the time span of each reach-set in this slice
# is sorted, e.g. F[1:4] so it suffices to look at the first and last reach-sets
function tend(F::VRT; contiguous=false) where {RT<:AbstractReachSet, VRT<:AbstractVector{RT}}
    if contiguous
        tf = tend(last(F))
    else
        tf = maximum(tend, F)
    end
    return tf
end

# if `contiguous` is true, we assume that the time span of each reach-set in this slice
# is sorted, e.g. F[1:4] so it suffices to look at the first and last reach-sets
function tspan(F::VRT; contiguous=false) where {RT<:AbstractReachSet, VRT<:AbstractVector{RT}}
    if contiguous
        t0 = tstart(first(F))
        tf = tend(last(F))
    else
        t0 = minimum(tstart, F)
        tf = maximum(tend, F)
    end
    return TimeInterval(t0, tf)
end

"""
    dim(R::AbstractReachSet)

Return the ambient dimension of the reach-set.

### Input

- `R` -- reach-set

### Output

An integer corresponding to the ambient dimension of the given reach-set.
"""
function LazySets.dim(::AbstractReachSet) end

# for internal use (dimensional checks in solve)
_dim(X::AbstractReachSet) = dim(X)

"""
    copy(R::AbstractReachSet)

Return a copy of the given reach-set.

### Input

- `R`  -- reach-set

### Output

A new reach-set of the sam type and the same field values as `R`.
"""
Base.copy(R::AbstractReachSet) = deepcopy(R)

"""
    shift(R::AbstractReachSet, t0::Number)

Perform a time-shift of the given reach-set.

### Input

- `R`  -- reach-set
- `t0` -- number that corresponds to the time-shift

### Output

A new reach-set of the same type of `R` such that its time-span has been shifted
by `t0`.
"""
function shift(R::AbstractReachSet, t0::Number) end

# ------------------------------
# Methods to check disjointness
# ------------------------------

# fallback uses  internal function _is_intersection_empty, which admit a pre-processing
# step for the reach-set / algorithm choice
@commutative isdisjoint(R::AbstractReachSet, Y::SetOrReachSet, method::AbstractDisjointnessMethod=FallbackDisjointness()) =_is_intersection_empty(R, Y, method)

# disambiguations
isdisjoint(R1::AbstractReachSet, R2::AbstractReachSet, method::AbstractDisjointnessMethod=FallbackDisjointness()) =_is_intersection_empty(R1, R2, method)

# vector of reach-sets
@commutative isdisjoint(R::AbstractVector{<:AbstractReachSet}, Y::SetOrReachSet, method::AbstractDisjointnessMethod=FallbackDisjointness()) = all(Ri -> _is_intersection_empty(Ri, Y, method), R)

# internal implementation
_is_intersection_empty(R::AbstractReachSet, Y::AbstractReachSet, method::FallbackDisjointness) = isdisjoint(set(R), set(Y))

@commutative _is_intersection_empty(R::AbstractReachSet, Y::LazySet, method::FallbackDisjointness) = isdisjoint(set(R), Y)

@commutative function _is_intersection_empty(R::AbstractReachSet, Y::LazySet, method::ZonotopeEnclosure)
    Z = overapproximate(R, Zonotope)
    return isdisjoint(set(Z), Y)
end

@commutative function _is_intersection_empty(R::AbstractReachSet, Y::LazySet, method::BoxEnclosure)
    H = overapproximate(R, Hyperrectangle)
    return isdisjoint(set(H), Y)
end

# used for disjointness check in continuos post-operators
function _is_intersection_empty(P::HPolyhedron,  Z::Zonotope, ::BoxEnclosure)
    return isdisjoint(P, box_approximation(Z))
end

# in this method we assume that the intersection is non-empty
@commutative function _is_intersection_empty(R::AbstractReachSet, Y::LazySet, method::Dummy)
    return false
end

# ------------------------------
# Concrete intersection methods
# ------------------------------

intersection(R::AbstractReachSet, S::AbstractReachSet) = _intersection(set(R), set(S), FallbackIntersection())
_intersection(R::AbstractReachSet, S::AbstractReachSet, method::AbstractIntersectionMethod) = _intersection(set(R), set(S), method)

# fallback methods for reach-sets
@commutative function _intersection(R::AbstractReachSet, X::LazySet, method::AbstractIntersectionMethod)
    return _intersection(set(R), X, method)
end

# ------------------------------
# Methods to check inclusion
# ------------------------------


text:
function _pwf2pm_info!(pm_data::Dict, pwf_data::Dict, option::String, status::Char, section::String)
    key = lowercase(option)
    value = status == 'L' ? true : status == 'D' ? false : error("Execution option $key status not defined")
    pm_data["info"]["actions"][key] = value
end

function _pwf2pm_info!(pm_data::Dict, pwf_data::Dict, option::String, value::Real, section::String)
    key = lowercase(option)
    pm_data["info"]["parameters"][key] = value
end

function _pwf2pm_info!(pm_data::Dict, pwf_data::Dict)
    pm_data["info"] = Dict("actions" => Dict{String,Any}(), "parameters" => Dict{String,Any}())
    info_sections = ["DOPC", "DOPC IMPR", "DCTE"]

    for section in info_sections
        if haskey(pwf_data, section)

            for (option,value) in pwf_data[section]
                _pwf2pm_info!(pm_data, pwf_data, option, value, section)
            end
        end
    end
end


text:
using ModelTES

"""Return a BiasedTES object containing parameters for a HOLMES notional TES"""
function pholmes(L=50e-9, RbiasFraction=0.2)
    n = 3.25; Tc = 0.1; Tbath = 0.07
    k = 2.33e-8; C = 0.5e-12
    Rl = 0.3e-3; Rn = 10e-3; Rpara=0.0; R0=Rn*RbiasFraction
    alpha = 200.0; beta = 2.0
    model = ModelTES.ShankRIT(alpha, beta, n, Tc, Tbath, k, R0, Rn);
    tes_param = TESParams(n,Tc,Tbath,k,C,L,Rl,Rpara,Rn,model)
    bt = BiasedTES(tes_param, R0)
end

"""Return a BiasedTES object containing parameters for Kelsey's low-E TES model"""
function lowEpix()
    Tc=0.094; Tbath = 0.075
    R0 = 1.64e-3; Rl = 0.35e-3; Rn = 8.2e-3; Rpara=0.0
    n=3.15; k = 1.48e-9; C = 0.1e-12
    alpha = 175.0; beta = 2.28; L=50e-9
    model = ModelTES.ShankRIT(alpha, beta, n, Tc, Tbath, k, R0, Rn);
    tes_param = TESParams(n,Tc,Tbath,k,C,L,Rl,Rpara,Rn,model)
    bt = BiasedTES(tes_param, R0)
end

"""Return a BiasedTES object containing parameters for Kelsey's high-E TES model"""
function highEpix()
    Tc=0.107; Tbath = 0.065
    R0 = 1.55e-3; Rl = 0.35e-3; Rn = 10.3e-3; Rpara=0.0
    n=3.3; k = 3.50e-9; C = 0.72e-12
    alpha = 92.0; beta = 1.4; L=300e-9
    model = ModelTES.ShankRIT(alpha, beta, n, Tc, Tbath, k, R0, Rn);
    tes_param = TESParams(n,Tc,Tbath,k,C,L,Rl,Rpara,Rn,model)
    bt = BiasedTES(tes_param, R0)
end

"""Return a BiasedTES object containin parameters for the Oct 2016 LCLS-II model"""
function LCLSII_2016(L=52e-9)
    Tc=0.075; Tbath = 0.055
    Rl = 0.34e-3; Rn = 13.5e-3; Rpara=0.0; R0=0.15*Rn
    n=3.5; G=130e-12; k = G/(n*Tc^(n-1)); C = 0.11e-12
    alpha = 120.0; beta = 1.69
    model = ModelTES.ShankRIT(alpha, beta, n, Tc, Tbath, k, R0, Rn);
    tes_param = TESParams(n,Tc,Tbath,k,C,L,Rl,Rpara,Rn,model)
    bt = BiasedTES(tes_param, R0)
end

"""Return a BiasedTES object containin parameters for the April 2018 LCLS-II model"""
function LCLSII(L=67e-9)
    Tc=0.050; Tbath = 0.030
    Rl = 0.35e-3; Rn = 13.3e-3; Rpara=0.0; R0=0.20*Rn
    n=3.25; G=170e-12; k = G/(n*Tc^(n-1)); C = 0.132e-12
    alpha = 197.0; beta = 1.6
    model = ModelTES.ShankRIT(alpha, beta, n, Tc, Tbath, k, R0, Rn);
    tes_param = TESParams(n,Tc,Tbath,k,C,L,Rl,Rpara,Rn,model)
    bt = BiasedTES(tes_param, R0)
end


text:
struct Cost{T}
    #TODO: types for methods
    evaluate
    gradient_state
    gradient_action
    hessian_state_state
    hessian_action_action
    hessian_action_state
    evaluate_cache::Vector{T}
    gradient_state_cache::Vector{T}
    gradient_action_cache::Vector{T}
    hessian_state_state_cache::Matrix{T}
    hessian_action_action_cache::Matrix{T}
    hessian_action_state_cache::Matrix{T}
end

function Cost(f::Function, num_state::Int, num_action::Int; num_parameter::Int=0)
    #TODO: option to load/save methods
    @variables x[1:num_state], u[1:num_action], w[1:num_parameter]

    evaluate = f(x, u, w)
    gradient_state = Symbolics.gradient(evaluate, x)
    gradient_action = Symbolics.gradient(evaluate, u)
    hessian_state_state = Symbolics.jacobian(gradient_state, x)
    hessian_action_action = Symbolics.jacobian(gradient_action, u)
    hessian_action_state = Symbolics.jacobian(gradient_action, x)

    evaluate_func = eval(Symbolics.build_function([evaluate], x, u, w)[2])
    gradient_state_func = eval(Symbolics.build_function(gradient_state, x, u, w)[2])
    gradient_action_func = eval(Symbolics.build_function(gradient_action, x, u, w)[2])
    hessian_state_state_func = eval(Symbolics.build_function(hessian_state_state, x, u, w)[2])
    hessian_action_action_func = eval(Symbolics.build_function(hessian_action_action, x, u, w)[2])
    hessian_action_state_func = eval(Symbolics.build_function(hessian_action_state, x, u, w)[2])

    return Cost(evaluate_func,
        gradient_state_func, gradient_action_func,
        hessian_state_state_func, hessian_action_action_func, hessian_action_state_func,
        zeros(1),
        zeros(num_state), zeros(num_action),
        zeros(num_state, num_state), zeros(num_action, num_action), zeros(num_action, num_state))
end

Objective{T} = Vector{Cost{T}} where T

function cost(costs::Vector{Cost{T}}, states, actions, parameters) where T
    J = 0.0
    for (t, cost) in enumerate(costs)
        cost.evaluate(cost.evaluate_cache, states[t], actions[t], parameters[t])
        J += cost.evaluate_cache[1]
    end
    return J
end

function cost_gradient!(gradient_states, gradient_actions, costs::Vector{Cost{T}}, states, actions, parameters) where T
    H = length(costs)
    for (t, cost) in enumerate(costs)
        cost.gradient_state(cost.gradient_state_cache, states[t], actions[t], parameters[t])
        @views gradient_states[t] .= cost.gradient_state_cache
        fill!(cost.gradient_state_cache, 0.0) # TODO: confirm this is necessary
        t == H && continue
        cost.gradient_action(cost.gradient_action_cache, states[t], actions[t], parameters[t])
        @views gradient_actions[t] .= cost.gradient_action_cache
        fill!(cost.gradient_action_cache, 0.0) # TODO: confirm this is necessary
    end
end

function cost_hessian!(hessian_state_state, hessian_action_action, hessian_action_state, costs::Vector{Cost{T}}, states, actions, parameters) where T
    H = length(costs)
    for (t, cost) in enumerate(costs)
        cost.hessian_state_state(cost.hessian_state_state_cache, states[t], actions[t], parameters[t])
        @views hessian_state_state[t] .+= cost.hessian_state_state_cache
        fill!(cost.hessian_state_state_cache, 0.0) # TODO: confirm this is necessary
        t == H && continue
        cost.hessian_action_action(cost.hessian_action_action_cache, states[t], actions[t], parameters[t])
        cost.hessian_action_state(cost.hessian_action_state_cache, states[t], actions[t], parameters[t])
        @views hessian_action_action[t] .+= cost.hessian_action_action_cache
        @views hessian_action_state[t]  .+= cost.hessian_action_state_cache
        fill!(cost.hessian_action_action_cache, 0.0) # TODO: confirm this is necessary
        fill!(cost.hessian_action_state_cache, 0.0) # TODO: confirm this is necessary
    end
end


text:
using StochasticDiffEq, Random
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
using DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_2Dlinear
Random.seed!(100)

#Let the solver determine the initial stepsize for you!
sol = solve(prob_sde_2Dlinear, SRI())

# using Plots
# plot(sol,plot_analytic=true)
# gui()

#Make sure it does a good job
sol.t[2] > 1e-7


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
config = MOIT.TestConfig()

@testset "NormInfinity" begin
    bridged_mock = MOIB.Constraint.NormInfinity{Float64}(mock)

    MOIT.basic_constraint_tests(bridged_mock, config,
        include = [(F, MOI.NormInfinityCone) for F in [
            MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}, MOI.VectorQuadraticFunction{Float64}
        ]])

    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 0.5, 1.0],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[0.0, 1.0, 0.0, 0.0]],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1], [-1]])

    MOIT.norminf1vtest(bridged_mock, config)
    MOIT.norminf1ftest(bridged_mock, config)

    @testset "Test mock model" begin
        var_names = ["x", "y", "z"]
        MOI.set(mock, MOI.VariableName(), MOI.get(mock, MOI.ListOfVariableIndices()), var_names)
        nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())
        @test length(nonneg) == 1
        MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")
        zeros = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Zeros}())
        @test length(zeros) == 2
        MOI.set(mock, MOI.ConstraintName(), zeros[1], "x_eq")
        MOI.set(mock, MOI.ConstraintName(), zeros[2], "y_eq")

        s = """
        variables: x, y, z
        nonneg: [x + -1.0y, x + -1.0z, x + y, x + z] in MathOptInterface.Nonnegatives(4)
        x_eq: [-1.0 + x] in MathOptInterface.Zeros(1)
        y_eq: [-0.5 + y] in MathOptInterface.Zeros(1)
        maxobjective: y + z
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(mock, model, var_names, ["nonneg", "x_eq", "y_eq"])
    end

    @testset "Test bridged model" begin
        var_names = ["x", "y", "z"]
        MOI.set(bridged_mock, MOI.VariableName(), MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)
        norminf = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormInfinityCone}())
        @test length(norminf) == 1
        MOI.set(bridged_mock, MOI.ConstraintName(), norminf[1], "norminf")
        zeros = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Zeros}())
        @test length(zeros) == 2
        MOI.set(bridged_mock, MOI.ConstraintName(), zeros[1], "x_eq")
        MOI.set(bridged_mock, MOI.ConstraintName(), zeros[2], "y_eq")

        s = """
        variables: x, y, z
        norminf: [1.0x, y, z] in MathOptInterface.NormInfinityCone(3)
        x_eq: [-1.0 + x] in MathOptInterface.Zeros(1)
        y_eq: [-0.5 + y] in MathOptInterface.Zeros(1)
        maxobjective: y + z
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(bridged_mock, model, var_names, ["norminf", "x_eq", "y_eq"])
    end

    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormInfinityCone}()))

    @testset "$attr" for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = [4.0, 1.0, -2.0]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
    end

    test_delete_bridge(bridged_mock, ci, 3, ((MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0),))
end

@testset "NormOne" begin
    bridged_mock = MOIB.Constraint.NormOne{Float64}(mock)

    MOIT.basic_constraint_tests(bridged_mock, config,
        include = [(F, MOI.NormOneCone) for F in [
            MOI.VectorOfVariables, MOI.VectorAffineFunction{Float64}, MOI.VectorQuadraticFunction{Float64}
        ]])

    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 0.5, 0.5, 0.5, 0.5],
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) => [[1.0, 1.0, 0.0, 0.0]],
        (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [[1.0]],
        (MOI.VectorAffineFunction{Float64}, MOI.Zeros) => [[-1], [0]])

    MOIT.normone1vtest(bridged_mock, config)
    MOIT.normone1ftest(bridged_mock, config)

    var_names = ["x", "y", "z"]
    MOI.set(bridged_mock, MOI.VariableName(), MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)

    nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives}())
    greater = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}}())
    u, v = MOI.get(mock, MOI.ListOfVariableIndices())[4:5]
    @testset "Test mock model" begin
        MOI.set(mock, MOI.VariableName(), u, "u")
        MOI.set(mock, MOI.VariableName(), v, "v")
        @test length(nonneg) == 1
        MOI.set(mock, MOI.ConstraintName(), nonneg[1], "nonneg")
        @test length(greater) == 1
        MOI.set(mock, MOI.ConstraintName(), greater[1], "greater")
        zeros = MOI.get(mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Zeros}())
        @test length(zeros) == 2
        MOI.set(mock, MOI.ConstraintName(), zeros[1], "x_eq")
        MOI.set(mock, MOI.ConstraintName(), zeros[2], "y_eq")

        s = """
        variables: x, y, z, u, v
        nonneg: [u + -1.0y, v + -1.0z, u + y, v + z] in MathOptInterface.Nonnegatives(4)
        greater: x + -1.0u + -1.0v >= 0.0
        x_eq: [-1.0 + x] in MathOptInterface.Zeros(1)
        y_eq: [-0.5 + y] in MathOptInterface.Zeros(1)
        maxobjective: y + z
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(mock, model, [var_names; "u"; "v"], ["nonneg", "greater", "x_eq", "y_eq"])
    end

    @testset "Test bridged model" begin
        normone = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormOneCone}())
        @test length(normone) == 1
        MOI.set(bridged_mock, MOI.ConstraintName(), normone[1], "normone")
        zeros = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.Zeros}())
        @test length(zeros) == 2
        MOI.set(bridged_mock, MOI.ConstraintName(), zeros[1], "x_eq")
        MOI.set(bridged_mock, MOI.ConstraintName(), zeros[2], "y_eq")

        s = """
        variables: x, y, z
        normone: [1.0x, y, z] in MathOptInterface.NormOneCone(3)
        x_eq: [-1.0 + x] in MathOptInterface.Zeros(1)
        y_eq: [-0.5 + y] in MathOptInterface.Zeros(1)
        maxobjective: y + z
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(bridged_mock, model, var_names, ["normone", "x_eq", "y_eq"])
    end

    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.NormOneCone}()))

    @testset "$attr" for attr in [MOI.ConstraintPrimalStart(), MOI.ConstraintDualStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = [4.0, 1.0, -2.0]
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
        if attr isa MOI.ConstraintPrimalStart
            @test MOI.get(mock, MOI.VariablePrimalStart(), u) == 1
            @test MOI.get(mock, MOI.VariablePrimalStart(), v) == 2
            @test MOI.get(mock, attr, nonneg[1]) == [0.0, 4.0, 2.0, 0.0]
            @test MOI.get(mock, attr, greater[1]) == 1
        else
            @test MOI.get(mock, attr, nonneg[1]) == [0.0, 2.0, 1.0, 0.0]
            @test MOI.get(mock, attr, greater[1]) == 4
        end
    end

    test_delete_bridge(bridged_mock, ci, 3, (
        (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
        (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, 0)))
end


text:
metadata = fill("metadata", 4)
an1 = ArrayNode(rand(Float32, 3, 4))
b = BagNode(an1, [1:4, 0:-1], metadata)
an2 = ArrayNode(NGramMatrix(["test", "skunk", "mill", "julia"], 3, 10, 17))
wb = WeightedBagNode(an2, [1:2,3:4], Float32[1,2,3,4], metadata)
n1 = ProductNode((b=b,wb=wb))
an3 = ArrayNode(SparseMatrixCSC(rand(Float32, 10, 2)))
n2 = ProductNode((n1, an3))

n2m = reflectinmodel(n2)
n1m, an3m = n2m.ms
bm, wbm = n1m.ms
an1m = bm.im
an2m = wbm.im

ss = ["GGGCGGCGA", "CCTCGCGGG", "TTTTCGCTATTTATGAAAATT", "TTCCGGTTTAAGGCGTTTCCG"]
base_arr = collect("ACGT")
possibilities = product(base_arr, base_arr, base_arr) |> collect |> x->reshape(x, 1, :) .|> (x->reduce(*, x)) |> x->x[:]

struct NumberNode
    n::Int
    chs::Vector{NumberNode}
end
import HierarchicalUtils: NodeType, LeafNode, InnerNode, noderepr, children
NodeType(::Type{NumberNode}) = InnerNode()
noderepr(n::NumberNode) = string(n.n)
children(n::NumberNode) = n.chs

function Mill.unpack2mill(ds::LazyNode{:Codons})
    s = ds.data
    ss = map(x -> reduce.(*, partition(x, 3)),s)
    x = reduce(hcat, map(x->Flux.onehotbatch(x, possibilities), ss))
    BagNode(ArrayNode(x), Mill.length2bags(length.(ss)))
end

# specification of printing
NodeType(::Type{<:LazyNode{:Codons}}) = InnerNode()
children(n::LazyNode{:Codons}) = (n.data,)
NodeType(n::Vector{<:AbstractString}) = InnerNode()
noderepr(n::Vector{<:AbstractString}) = "Array $(length(n)) items"
NodeType(n::AbstractString) = LeafNode()
children(n::Vector{<:AbstractString}) = (n...,)

t2 = treemap(n2, n2m) do (k1, k2), chs
    NumberNode(rand(1:10), collect(chs))
end

function buf_printtree(data; kwargs...)
    buf = IOBuffer()
    printtree(buf, data; kwargs...)
    String(take!(buf))
end

@testset "list traversal" begin
    for (n1, n2) in NodeIterator(n2, n2m)
        @test list_traversal(n1) == list_traversal(n2)
    end
end

@testset "getindex on strings" begin
    @test n2[""] === n2
    @test n2["E"] === n1
    @test n2["I"] === b
    @test n2["K"] === an1
    @test n2["M"] === wb
    @test n2["O"] === an2
    @test n2["U"] === an3

    @test n2m[""] === n2m
    @test n2m["E"] === n1m
    @test n2m["I"] === bm
    @test n2m["K"] === an1m
    @test n2m["M"] === wbm
    @test n2m["O"] === an2m
    @test n2m["U"] === an3m
end

@testset "children" begin
    @test Set(children(n2)) == Set([n1, an3])
    @test Set(children(n1)) == Set([b, wb])
    @test Set(children(b)) == Set([an1])
    @test Set(children(wb)) == Set([an2])

    @test Set(children(n2m)) == Set([n1m, an3m])
    @test Set(children(n1m)) == Set([bm, wbm])
    @test Set(children(bm)) == Set([an1m])
    @test Set(children(wbm)) == Set([an2m])
end

@testset "nchildren" begin
    @test nchildren(an1) == nchildren(an2) == nchildren(an3) == 0
    @test nchildren(b) == nchildren(wb) == 1
    @test nchildren(n1) == nchildren(n2) == 2

    @test nchildren(an1m) == nchildren(an2m) == nchildren(an3m) == 0
    @test nchildren(bm) == nchildren(wbm) == 1
    @test nchildren(n1m) == nchildren(n2m) == 2
end

@testset "nnodes" begin
    @test nnodes(an1) == nnodes(an2) == nnodes(an3) == 1
    @test nnodes(wb) == nnodes(b) == 2
    @test nnodes(n1) == nnodes(b) + nnodes(wb) + 1
    @test nnodes(n2) == nnodes(n1) + nnodes(an3) + 1

    @test nnodes(an1m) == nnodes(an2m) == nnodes(an3m) == 1
    @test nnodes(wbm) == nnodes(bm) == 2
    @test nnodes(n1m) == nnodes(bm) + nnodes(wbm) + 1
    @test nnodes(n2m) == nnodes(n1m) + nnodes(an3m) + 1
end

@testset "nleafs" begin
    @test nleafs(an1) == nleafs(an2) == nleafs(an3) == 1
    @test nleafs(wb) == nleafs(b) == 1
    @test nleafs(n1) == nleafs(b) + nleafs(wb)
    @test nleafs(n2) == nleafs(n1) + nleafs(an3)

    @test nleafs(an1m) == nleafs(an2m) == nleafs(an3m) == 1
    @test nleafs(wbm) == nleafs(bm) == 1
    @test nleafs(n1m) == nleafs(bm) + nleafs(wbm)
    @test nleafs(n2m) == nleafs(n1m) + nleafs(an3m)
end

@testset "NodeIterator" begin
    @test Set(NodeIterator(n2)) == Set([an1, an2, an3, b, wb, n1, n2])

    @test Set(NodeIterator(n2m)) == Set([an1m, an2m, an3m, bm, wbm, n1m, n2m])
end

@testset "LeafIterator" begin
    @test Set(LeafIterator(n2)) == Set([an1, an2, an3])

    @test Set(LeafIterator(n2m)) == Set([an1m, an2m, an3m])
end

@testset "TypeIterator" begin
    @test Set(TypeIterator(AbstractMillNode, n2)) == Set(NodeIterator(n2))
    @test Set(TypeIterator(AbstractBagNode, n2)) == Set([b, wb])
    @test Set(TypeIterator(AbstractProductNode, n2)) == Set([n1, n2])

    @test Set(TypeIterator(AbstractMillModel, n2m)) == Set(NodeIterator(n2m))
    @test Set(TypeIterator(BagModel, n2m)) == Set([bm, wbm])
    @test Set(TypeIterator(ProductModel, n2m)) == Set([n1m, n2m])
end

@testset "Iteration over multiple trees" begin
    @test NodeIterator(n2, n2m) |> collect == collect(zip(NodeIterator(n2) |> collect, NodeIterator(n2m) |> collect))
end

@testset "printtree" begin
    @test buf_printtree(n2, trav=true) ==
        """
        ProductNode with 2 obs [""]
          ├── ProductNode with 2 obs ["E"]
          │     ├─── b: BagNode with 2 obs ["I"]
          │     │         └── ArrayNode(3×4 Array with Float32 elements) ["K"]
          │     └── wb: WeightedBagNode with 2 obs ["M"]
          │               └── ArrayNode(17×4 NGramMatrix with Int64 elements) ["O"]
          └── ArrayNode(10×2 SparseMatrixCSC with Float32 elements) ["U"]
        """

    @test buf_printtree(n2m, trav=true) ==
        """
        ProductModel … ↦ ArrayModel(Dense(20, 10)) [""]
          ├── ProductModel … ↦ ArrayModel(Dense(20, 10)) ["E"]
          │     ├─── b: BagModel … ↦ BagCount([SegmentedMean(10); SegmentedMax(10)]) ↦ ArrayModel(Dense(21, 10)) ["I"]
          │     │         └── ArrayModel(Dense(3, 10)) ["K"]
          │     └── wb: BagModel … ↦ BagCount([SegmentedMean(10); SegmentedMax(10)]) ↦ ArrayModel(Dense(21, 10)) ["M"]
          │               └── ArrayModel(Dense(17, 10)) ["O"]
          └── ArrayModel(Dense(10, 10)) ["U"]
        """
end

@testset "LazyNode" begin
	ds = LazyNode{:Codons}(ss)
	m = Mill.reflectinmodel(ds)
	@test nchildren(ds) == 1
	@test nleafs(ds) == 4

    @test buf_printtree(ds, trav=true) ==
        """
        LazyNode{Codons} with 4 obs [""]
          └── Array 4 items ["U"]
                ├── "GGGCGGCGA" ["Y"]
                ├── "CCTCGCGGG" ["c"]
                ├── "TTTTCGCTATTTATGAAAATT" ["g"]
                └── "TTCCGGTTTAAGGCGTTTCCG" ["k"]
        """

    @test buf_printtree(m, trav=true) ==
        """
        LazyModel{Codons} [""]
          └── BagModel … ↦ BagCount([SegmentedMean(10); SegmentedMax(10)]) ↦ ArrayModel(Dense(21, 10)) ["U"]
                └── ArrayModel(Dense(64, 10)) ["k"]
        """
end


text:
using BenchmarkTools
using CategoricalArrays
using Printf

const SUITE = BenchmarkGroup()

SUITE["isequal"] = BenchmarkGroup()

function sumequals(A::AbstractArray, v::Any)
    n = 0
    @inbounds for x in A
        n += isequal(x, v)
    end
    n
end

ca = CategoricalArray(repeat(string.('A':'J'), outer=1000))
nca = CategoricalArray(repeat([missing; string.('A':'J')], outer=1000))

SUITE["isequal"]["CategoricalArray{String}", "String"] =
    @benchmarkable  sumequals(ca, "D")
SUITE["isequal"]["CategoricalArray{Union{String, Missing}}", "String"] =
    @benchmarkable sumequals(nca, "D")

function sumequals(A::AbstractArray, v::CategoricalValue)
    n = 0
    @inbounds for x in A
        n += isequal(x, v)
    end
    n
end

ca = CategoricalArray(repeat(string.('A':'J'), outer=1000))
nca = CategoricalArray(repeat([missing; string.('A':'J')], outer=1000))

SUITE["isequal"]["CategoricalArray{String}", "CategoricalValue"] =
    @benchmarkable sumequals(ca, ca[1])
SUITE["isequal"]["CategoricalArray{Union{String, Missing}}", "CategoricalValue"] =
    @benchmarkable sumequals(nca, nca[1])

SUITE["isequal"]["CategoricalArray{String}", "CategoricalValue different pool"] =
    @benchmarkable sumequals(ca, nca[1])
SUITE["isequal"]["CategoricalArray{Union{String, Missing}}", "CategoricalValue  different pool"] =
    @benchmarkable sumequals(nca, ca[1])

# With many levels, checking whether some levels have been dropped can be very slow (#93)
SUITE["many levels"] = BenchmarkGroup()

a = rand([@sprintf("id%010d", k) for k in 1:1000], 10000)

SUITE["many levels"]["CategoricalArray(::Vector{String})"] =
    @benchmarkable CategoricalArray(a)

a = rand([@sprintf("id%010d", k) for k in 1:1000], 10000)
ca = CategoricalArray(a)

levs = levels(ca)
SUITE["many levels"]["levels! with original levels"] =
    @benchmarkable levels!(ca, levs)

levs = reverse(levels(ca))
SUITE["many levels"]["levels! with resorted levels"] =
    @benchmarkable levels!(ca, levs)

levs = [levels(ca); [@sprintf("id2%010d", k) for k in 1:1000]]
SUITE["many levels"]["levels! with many additional levels"] =
    @benchmarkable levels!(ca, levs)

SUITE["repeated assignment"] = BenchmarkGroup()

function mycopy!(dest, src)
    @inbounds for i in eachindex(dest, src)
        dest[i] = src[i]
    end
end

a = categorical(rand([@sprintf("id%010d", k) for k in 1:1000], 10000))
b = CategoricalArray{String}(undef, 10000)
c = categorical(rand([@sprintf("id%010d", k) for k in 1:1000], 10000))
d = categorical(rand([@sprintf("id%010d", k) for k in 1001:2000], 10000))

SUITE["repeated assignment"]["empty dest"] =
    @benchmarkable mycopy!(b2, a) setup = b2=copy(b)
SUITE["repeated assignment"]["same levels dest"] =
    @benchmarkable mycopy!(c2, a) setup = c2=copy(c)
SUITE["repeated assignment"]["many levels dest"] =
    @benchmarkable mycopy!(d2, a) setup = d2=copy(d)

orig_vec = (x -> repeat(x, 32)).(string.([x % 1000 for x in 1:1000000]))
cat2merge_vec = (x -> repeat(x, 32)).(string.([x % 1000 for x in 1:100000]))
SUITE["recode"] = BenchmarkGroup()
SUITE["recode"]["vectors"] = @benchmarkable recode(orig_vec, cat2merge_vec => "None")
SUITE["recode"]["categorical_vectors"] = @benchmarkable recode(categorical(orig_vec), cat2merge_vec => "None")
SUITE["recode"]["matrices"] = @benchmarkable recode(reshape(orig_vec, :, 1), reshape(cat2merge_vec, :, 1) => "None")
SUITE["recode"]["categorical_matrices"] = @benchmarkable recode(categorical(reshape(orig_vec, :, 1)), reshape(cat2merge_vec, :, 1) => "None")


text:
#=using Pkg
Pkg.activate(joinpath(dirname(@__FILE__), "../../../"))
using Test, OrdinaryDiffEq, HDF5, ModelConstructors, Plots, NLsolve, Statistics, Roots
include(joinpath(dirname(@__FILE__), "../../../src/includeall.jl"))
=#
start_time = time_ns()
# This only uses q for the input
default = false
nlsolve_newton = false
nlsolve_trustreg = false
forwarddiff = false
use_find_zero = true # 3 times as slow
verbose = false
if default
    nlsolve_trustreg = false
    use_find_zero = false
    nlsolve_newton = false
end

## Parameters
ρₑ = 0.03
ρₕ = 0.03
τ  = .005
aₑ = .15
aₕ = .7 * aₑ
σ  = .03
γ  = 2.
χ₁ = 1.
χ₂ = 5.
δ  = 0.03
Δ  = .8
pnorm = Inf
T = 1000
tol = 1e-6
update_hyp = .9 # 1 is all weight on the ratio of a rolling window of error averages, 0 is none
reduce_hyp = .95 # proportion by which to reduce Δ if the error grows in size
discounting = 8

## Grid
# uneven grid on [0, 1] with 1000 points, with more points near 0 and 1

N   = 300 # 1000
zz  = collect(range(0.001, stop = 0.999, length = N))
# η   = 3*zz.^2  - 2*zz.^3
η   = zz
dη  = diff(η)

default_err = fill(NaN, N, 3)
nlsolve_newton_err = fill(NaN, N)
nlsolve_trustreg_err = fill(NaN, N)
nlsolve_findzero_err = fill(NaN, N)

# default_vₑ = copy(vₑ)
# default_vₕ = copy(vₕ)
## Terminal conditions for value functions V and Vₕ
vₑ = aₑ^(-γ) .* η .^ (1. - γ)
vₕ = aₑ^(-γ) .* (1. .- η) .^ (1. - γ)
vₑ₀ = similar(vₑ)
vₕ₀ = similar(vₕ)
vₑ₁ = similar(vₑ)
vₕ₁ = similar(vₕ)
Q   =  ones(N)
Qp  = zeros(N)
SSQ = zeros(N)
Psi = zeros(N)
Iota   = zeros(N)
Φ   = zeros(N)
A   = zeros(N)
S   = zeros(N)
MU  = zeros(N)
Mvₑ  = zeros(N)
Mvₕ = zeros(N)
Mvₑfd  = zeros(N)
Mvₕfd = zeros(N)
G   = zeros(N)
Svₑ  = zeros(N)
Svₕ = zeros(N)
CN  = zeros(N)
C_N = zeros(N)

F = zeros(1) # For static step

function staticstep!(F, x, η, dη, q_old, GS, VVlp, region)
    if region == 1
        q   = x[1]
        ψ   = (q^(1 / γ) * GS - (aₕ - (χ₁ * q - 1) / χ₂)) / (aₑ - aₕ)

        ssq = sqrt((aₑ - aₕ) / q / (ψ - η) / VVlp)
        F[1] = ssq * (q - (q - q_old) * (ψ - η) / dη) - σ * q
        if verbose
            @show q, ψ, η, ssq, F[1]
        end
        # F[1] = log(q) / γ + log(GS) - log(aₑ * ψ + aₕ * (1. - ψ) - (χ₁ * q - 1) / χ₂)
        # F[2] = ssq * (q - (q - q_old) * (ψ - η) / dη) - σ * q
        # F[3] = aₑ - aₕ - q * (ψ - η) * ssq^2 * VVlp
    elseif region == 2
        q    = x[1]
        F[1] = log(q) / γ + log(GS) - log(aₑ - (χ₁ * q - 1.) / χ₂)
    end
    nothing
end

function staticstep(x, η, dη, q_old, GS, VVlp, region)
    F = zeros(1)
    staticstep!(F, x, η, dη, q_old, GS, VVlp, region)
    return F[1]
end

t_err_vec = fill(NaN, T)
Δ_vec = fill(NaN, T)
dobreak = Bool[false]
Δ_vec[1:2] .= Δ
for t in 1:T

    ## Initialization
    GS1 = (η        ./ vₑ  ) .^ (1. / γ)
    GS2 = ((1 .- η) ./ vₕ ) .^ (1. / γ)
    GS  =  GS1 + GS2

    vₑp = diff(vₑ)  ./ dη
    vₕp = diff(vₕ) ./ dη

    vₑpl  = vₑp  ./ vₑ[2:N]
    vₕpl = vₕp ./ vₕ[2:N]

    VVlp = vₕpl - vₑpl + 1. ./ (η[2:N] .* (1. .- η[2:N]))

    ## Initial condition at η[1]
    Psi[1] = 0.
    SSQ[1] = σ
    qL     = [0.]
    qR     = [aₑ * χ₂ + 1.]
    q₀     = [0.]
    for k = 1:30
        q₀[1] = (qL[1] + qR[1])/2
        ι     = (χ₁ * q₀[1] - 1.) / χ₂
        A₀   = aₕ - ι

        if log(q₀[1]) / γ + log(GS[1]) > log(A₀)
            qR[1] = q₀[1]
        else
            qL[1] = q₀[1]
        end
    end
    Q[1]  = q₀[1]

    ## Newton method
    #  find q, psi and ssq = σ + σ^q from value functions
    break_pt = Int[0]
    EN = zeros(3)
    minq = zeros(1)
    maxq = zeros(1)
    for n = 2:N
        q     = Q[n - 1]
        q_old = Q[n - 1]
        ψ     = Psi[n - 1]
        ssq   = SSQ[n - 1]

        # errors given guess
        if nlsolve_newton || nlsolve_trustreg || use_find_zero
            minq[1] = find_zero(x -> η[n] - (x^(1 / γ) * GS[n] - (aₕ - (χ₁ * x - 1) / χ₂)) / (aₑ - aₕ), minq[1])
            maxq[1] = find_zero(x -> 1. - (x^(1 / γ) * GS[n] - (aₕ - (χ₁ * x - 1) / χ₂)) / (aₑ - aₕ), maxq[1])
            if verbose
                @show minq, maxq
            end
            # xₙ₋₁  = [q]
            # xₙ₋₁ = q
            # xₙ₋₁  = [(minq + maxq) / 2.]
            xₙ₋₁  = (minq[1] + maxq[1]) / 2.
            qₙ = if forwarddiff
                nlsolve((F, x) -> staticstep!(F, x, η[n], dη[n - 1], q_old, GS[n], VVlp[n - 1], 1), xₙ₋₁,
                         method = nlsolve_newton ? :newton : :trust_region, autodiff = :forward).zero
            elseif use_find_zero
                find_zero(x -> staticstep(x, η[n], dη[n - 1], q_old, GS[n], VVlp[n - 1], 1), xₙ₋₁) # (minq + maxq) / 2.)
            else
                nlsolve((F, x) -> staticstep!(F, x, η[n], dη[n - 1], q_old, GS[n], VVlp[n - 1], 1), xₙ₋₁,
                         method = nlsolve_newton ? :newton : :trust_region).zero
            end

            staticstep!(F, qₙ, η[n], dη[n - 1], q_old, GS[n], VVlp[n - 1], 1)
            if nlsolve_newton
                nlsolve_newton_err[n] = F[1]
            elseif nlsolve_trustreg
                nlsolve_trustreg_err[n] = F[1]
            elseif use_find_zero
                nlsolve_findzero_err[n] = F[1]
            end

            EN[1] = qₙ[1]
            EN[2] = (qₙ[1]^(1 / γ) * GS[n] - (aₕ - (χ₁ * qₙ[1] - 1.) / χ₂)) / (aₑ - aₕ)
            EN[3] = σ * qₙ[1] / (qₙ[1] - (qₙ[1] - q_old) * (EN[2] - η[n]) / dη[n - 1])
        elseif default
            xₙ₋₁  = [q, ψ, ssq]

            ER = [log(q) / γ + log(GS[n]) - log(aₑ * ψ + aₕ * (1. - ψ) - (χ₁ * q - 1) / χ₂),
                  ssq * (q - (q - q_old) * (ψ - η[n]) / dη[n - 1]) - σ * q,
                  aₑ - aₕ - q * (ψ - η[n]) * ssq^2 * VVlp[n - 1]]

            # matrix of derivatives of errors
            # could shorten it since q_old = q

            QN = zeros(3,3)

            QN[1, :] = [1 / (q * γ) + 1/((aₑ - aₕ) * ψ + aₕ - (q - 1.) / χ₂) / χ₂,
                        -(aₑ - aₕ) / ((aₑ - aₕ) * ψ + aₕ - (q - 1.) / χ₂), 0.]

            QN[2, :] = [ssq * (1. - (ψ - η[n]) / dη[n - 1]) - σ,
                        -ssq * (q - q_old) / dη[n-1],  q - (q - q_old) * (ψ - η[n]) / dη[n-1]]

            QN[3, :] = [- (ψ - η[n]) * ssq^2 * VVlp[n-1],
                        -q * ssq^2 * VVlp[n - 1], -2 * q * (ψ - η[n]) * ssq * VVlp[n - 1]]

            EN = [q, ψ, ssq] - QN \ ER
            staticstep!(F, EN, η[n], dη[n - 1], q_old, GS[n], VVlp[n - 1], 1)
            default_err[n, :] .= F
        end

#=            use_nlsolve_newton = false
            use_nlsolve_trustreg = false
            use_nlsolve_newton_qonly = false
            use_nlsolve_trustreg_qonly = false=#

        # if the boundary of the crisis regime has been reached, we have
        # ψ = 1 from now on
        if EN[2] > 1.
            break_pt[1] = n
            break
        end

        # save results
        Q[n]   = EN[1]
        Psi[n] = EN[2]
        SSQ[n] = EN[3]
        Qp[n]  = (Q[n] - q_old) / dη[n-1]
    end

    ## Newton method for ψ = 1, for remaining eta
    for n in break_pt[1]:N
        q = Q[n - 1]

        ER = log(q) / γ + log(GS[n]) - log(aₑ - (χ₁ * q - 1.) / χ₂)
        QN = 1. / (q * γ)  + 1. / (aₑ - (q - 1.) / χ₂) / χ₂
        EN1 = q - ER / QN

        Qp[n]  = (EN1 - q) / dη[n-1]
        Q[n]   =  EN1
        Psi[n] = 1.
        SSQ[n] = 1. / (1. - (1. - η[n]) * Qp[n] / Q[n]) * σ
    end

    ## Computing the PDE
    # Volatility of η
    S     = (Psi - η) .* SSQ # This is σ_η * η
    S[1]  = 0.
    S[N]  = 0.
    Iota .= (χ₁ .* Q .- 1.) ./ χ₂
    Φ    .= log.(χ₁ .* Q) / χ₂
    A    .= aₑ .* Psi + aₕ .* (1. .- Psi) - Iota

    # for the law of motion of η, C/N = (η Q)^(1/γ - 1)/V^(1/γ)
    # recall that GS1 = (η./V).^(1/γ)
    CN  = GS1 .* Q .^ (1. / γ - 1.) ./ η # C/N*η
    C_N = GS2 .* Q.^(1. / γ - 1.) ./ (1. .- η)

    Svₑ[2:N]  = vₑpl .* S[2:N]
    Svₕ[2:N] = vₕpl .* S[2:N]

    VarSig  = -Svₑ + S ./ η        + SSQ .+ (γ - 1) * σ
    VarSig_ = -Svₕ + S ./ (1. .- η) + SSQ .+ (γ - 1) * σ

    ind = 2:(N - 1)
    MU[ind] .= ((aₑ .- Iota[ind]) ./ Q[ind] - CN[ind] .- τ) .* η[ind] + S[ind] .* (VarSig[ind] - SSQ[ind]) # This is μ_η * η

    Mvₑ = ρₑ .- CN  - (1 - γ) .* (Φ .- δ + Svₑ .* σ .- γ * σ^2 / 2.)
    Mvₕ = ρₕ .- C_N - (1 - γ) .* (Φ .- δ + Svₕ .* σ .- γ * σ^2 / 2.)

    for ind in 2:(N - 1)
        Mvₑfd[ind] = ((vₑ[ind + 1] - vₑ[ind]) / dη[ind] * max(0, MU[ind]) +
            (vₑ[ind] - vₑ[ind - 1]) / dη[ind - 1] * min(0, MU[ind])) / vₑ[ind] +
            (vₑ[ind + 1] - 2 * vₑ[ind] + vₑ[ind - 1]) / (dη[ind] * dη[ind - 1]) / 2 / vₑ[ind] * S[ind]^2;
        Mvₕfd[ind] = ((vₕ[ind + 1] - vₕ[ind]) / dη[ind] * max(0, MU[ind]) +
            (vₕ[ind] - vₕ[ind - 1]) / dη[ind - 1] * min(0, MU[ind])) / vₕ[ind] +
            (vₕ[ind + 1] - 2 * vₕ[ind] + vₕ[ind - 1]) / (dη[ind] * dη[ind - 1]) / 2 / vₕ[ind] * S[ind]^2;
    end

    t_err_vec[t] = norm(vcat(Mvₑfd[2:end - 1] - Mvₑ[2:end - 1], Mvₕfd[2:end - 1] - Mvₕ[2:end - 1]), pnorm)

    if dobreak[1]
        break
    end
    ## Updating vₑ and vₕ
    # last parameter is dt*ρₑ if dt is small, can be at most 1 (1 corresponds to policy iteration)
    # it is more agressive to set the last parameter closer to 1, but code may not converge

    vₑ₁ .= upwind_parabolic_pde(η, Mvₑ, MU, S.^2, G, vₑ, Δ_vec[t])
    vₕ₁ .= upwind_parabolic_pde(η, Mvₕ, MU, S.^2, G, vₕ, Δ_vec[t])

    if t_err_vec[t] < tol
        dobreak[1] = true
    end
    if t > 1 && t < T
        if t_err_vec[t] / t_err_vec[t - 1] < 1.
            Δ_weight = (1 - update_hyp) + update_hyp * mean(t_err_vec[max(2, t - discounting):t] ./ t_err_vec[max(1, t - discounting - 1):t - 1])
            Δ_vec[t + 1] = (1 - Δ_weight) + Δ_weight * Δ_vec[t]

            vₑ₀ .= vₑ
            vₕ₀ .= vₕ
            vₑ .= vₑ₁
            vₕ .= vₕ₁
        elseif t_err_vec[t] / t_err_vec[t - 1] > 1
            Δ_vec[t + 1] = reduce_hyp * Δ_vec[t];
            vₑ .= vₑ₀
            vₕ .= vₕ₀
        end
    else
        vₑ .= vₑ₁
        vₕ .= vₕ₁
    end
end

end_time = time_ns()

elasped_time = (end_time - start_time) / 1e9
@show elasped_time

breakT = findfirst(isnan.(t_err_vec))
if isnothing(breakT)
    breakT = T
else
    breakT -= 1
    t_err_vec = t_err_vec[1:breakT]
    Δ_vec = Δ_vec[1:breakT]
end
@show breakT

nothing