text:
"""
    Selector

Abstract supertype for all selectors.

Selectors are wrappers that indicate that passed values are not the array indices,
but values to be selected from the dimension index, such as `DateTime` objects for
a `Ti` dimension.

Selectors provided in DimensionalData are:

- [`At`](@ref)
- [`Between`](@ref)
- [`Near`](@ref)
- [`Where`](@ref)
- [`Contains`](@ref)

"""
abstract type Selector{T} end

val(sel::Selector) = sel.val
Base.parent(sel::Selector) = sel.val

const SelectorOrInterval = Union{Selector,Interval}

const SelTuple = Tuple{<:SelectorOrInterval,Vararg{<:SelectorOrInterval}}

"""
    At <: Selector

    At(x, atol, rtol)
    At(x; atol=nothing, rtol=nothing)

Selector that exactly matches the value on the passed-in dimensions, or throws an error.
For ranges and arrays, every intermediate value must match an existing value -
not just the end points.

`x` can be any value or `Vector` of values.

`atol` and `rtol` are passed to `isapprox`.
For `Number` `rtol` will be set to `Base.rtoldefault`, otherwise `nothing`,
and wont be used.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(5:7)))
A[X(At(20)), Y(At(6))]

# output

5
```
"""
struct At{T,A,R} <: Selector{T}
    val::T
    atol::A
    rtol::R
end
At(val; atol=nothing, rtol=nothing) = At(val, atol, rtol)

atol(sel::At) = sel.atol
rtol(sel::At) = sel.rtol

struct _True end
struct _False end

selectindices(l::LookupArray, sel::At; kw...) = at(l, sel; kw...)
selectindices(l::LookupArray, sel::At{<:AbstractVector}) = _selectvec(l, sel)

_selectvec(l, sel) = [selectindices(l, rebuild(sel; val=v)) for v in val(sel)]

function at(lookup::NoLookup, sel::At; kw...)
    v = val(sel)
    r = round(Int, v)
    at = atol(sel)
    if isnothing(at)
        v == r || _selvalnotfound(lookup, v)
    else
        at >= 0.5 && error("atol must be small than 0.5 for NoLookup")
        isapprox(v, r; atol=at) || _selvalnotfound(lookup, v)
    end
    r in lookup || throw(BoundsError(lookup, r))
    return r
end
function at(lookup::LookupArray, sel::At; kw...)
    at(order(lookup), lookup, val(sel), atol(sel), rtol(sel); kw...)
end
function at(
    ::Ordered, lookup::LookupArray{<:Union{Number,Dates.TimeType}}, selval, atol, rtol::Nothing;
    err=_True()
)
    x = unwrap(selval)
    i = searchsortedlast(lookup, x)
    # Try the current index
    if i === 0
        i1 = i + 1
        if checkbounds(Bool, lookup, i1) && _is_at(x, lookup[i1], atol)
            return i1
        else
            return _selnotfound_or_nothing(err, lookup, selval)
        end
    elseif _is_at(x, lookup[i], atol)
        return i
    else
        # Try again with the next index
        i1 = i + 1
        if checkbounds(Bool, lookup, i1) && _is_at(x, lookup[i1], atol)
            return i1
        else
            return _selnotfound_or_nothing(err, lookup, selval)
        end
    end
end
# catch-all for an unordered or non-number index
function at(order, lookup::LookupArray, selval, atol, rtol::Nothing; err=_True())
    i = findfirst(x -> _is_at(x, unwrap(selval), atol), parent(lookup))
    if i === nothing
        return _selnotfound_or_nothing(err, lookup, selval)
    else
        return i
    end
end

@inline _is_at(x, y, atol) = x == y
@inline _is_at(x::Real, y::Real, atol::Real) = abs(x - y) <= atol

_selnotfound_or_nothing(err::_True, lookup, selval) = _selvalnotfound(lookup, selval)
_selnotfound_or_nothing(err::_False, lookup, selval) = nothing
@noinline _selvalnotfound(lookup, selval) = throw(ArgumentError("$selval not found in $lookup"))

"""
    Near <: Selector

    Near(x)

Selector that selects the nearest index to `x`.

With [`Points`](@ref) this is simply the index values nearest to the `x`,
however with [`Intervals`](@ref) it is the interval _center_ nearest to `x`.
This will be offset from the index value for `Start` and
[`End`](@ref) loci.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(5:7)))
A[X(Near(23)), Y(Near(5.1))]

# output
4
```
"""
struct Near{T} <: Selector{T}
    val::T
end

selectindices(l::LookupArray, sel::Near) = near(l, sel)
selectindices(l::LookupArray, sel::Near{<:AbstractVector}) = _selectvec(l, sel)

near(lookup::NoLookup, sel::Near{<:Real}) = max(1, min(round(Int, val(sel)), lastindex(lookup)))
function near(lookup::LookupArray, sel::Near)
    span(lookup) isa Union{Irregular,Explicit} && locus(lookup) isa Union{Start,End} &&
        throw(ArgumentError("Near is not implemented for Irregular or Explicit with Start or End loci. Use Contains"))
    near(order(lookup), sampling(lookup), lookup, sel)
end
near(order::Order, ::NoSampling, lookup::LookupArray, sel::Near) = at(lookup, At(val(sel)))
function near(order::Ordered, ::Union{Intervals,Points}, lookup::LookupArray, sel::Near)
    # Unwrap the selector value and adjust it for
    # inderval locus if neccessary
    v = unwrap(val(sel))
    v_adj = _locus_adjust(locus(lookup), v, lookup)
    # searchsortedfirst or searchsortedlast
    searchfunc = _searchfunc(order)
    # Search for the value
    found_i = _inbounds(searchfunc(lookup, v_adj), lookup)

    # Check if this is the lowest possible value allready,
    # and return if so
    if order isa ForwardOrdered
        found_i <= firstindex(lookup) && return found_i
    elseif order isa ReverseOrdered
        found_i >= lastindex(lookup) && return found_i
    end

    # Find which index is nearest,
    # the found index or previous index
    prev_i = found_i - _ordscalar(order)
    dist_to_prev = abs(v_adj - lookup[prev_i])
    dist_to_found = abs(v_adj - lookup[found_i])
    # Compare distance to the found and previous index values
    # We have to use the correct >/>= for Start/End locus
    lessthan = _lt(locus(lookup))
    closest_i = lessthan(dist_to_prev, dist_to_found) ? prev_i : found_i

    return closest_i
end
function near(::Unordered, ::Union{Intervals,Points}, lookup::LookupArray, sel::Near)
    throw(ArgumentError("`Near` has no meaning in an `Unordered` `Sampled` index"))
end

_locus_adjust(locus::Center, v, lookup) = v
_locus_adjust(locus::Start, v, lookup) = v - abs(step(lookup)) / 2
_locus_adjust(locus::End, v, lookup) = v + abs(step(lookup)) / 2
_locus_adjust(locus::Start, v::DateTime, lookup) = v - (v - (v - abs(step(lookup)))) / 2
_locus_adjust(locus::End, v::DateTime, lookup) = v + (v + abs(step(lookup)) - v) / 2

"""
    Contains <: Selector

    Contains(x)

Selector that selects the interval the value is contained by. If the
interval is not present in the index, an error will be thrown.

Can only be used for [`Intervals`](@ref) or [`Categorical`](@ref).

## Example

```jldoctest
using DimensionalData; const DD = DimensionalData
dims_ = X(10:10:20; sampling=DD.Intervals(DD.Center())),
        Y(5:7; sampling=DD.Intervals(DD.Center()))
A = DimArray([1 2 3; 4 5 6], dims_)
A[X(Contains(8)), Y(Contains(6.8))]

# output
3
```
"""
struct Contains{T} <: Selector{T}
    val::T
end

# Filter based on sampling and selector -----------------
selectindices(l::LookupArray, sel::Contains; kw...) = contains(l, sel)
selectindices(l::LookupArray, sel::Contains{<:AbstractVector}) = _selectvec(l, sel)

function contains(l::NoLookup, sel::Contains; kw...)
    i = Int(val(sel))
    i in l || throw(BoundsError(l, i))
    return i
end
contains(l::LookupArray, sel::Contains; kw...) = contains(sampling(l), l, sel; kw...)
# NoSampling (e.g. Categorical) just uses `at`
function contains(::NoSampling, l::LookupArray, sel::Contains; kw...)
    at(l, At(val(sel)); kw...)
end
# Points --------------------------------------
function contains(::Points, l::LookupArray, sel::Contains; err=_True())
    if err isa _True
        throw(ArgumentError("Points LookupArray cannot use `Contains`, use `Near` or `At` for Points."))
    else
        nothing
    end
end
# Intervals -----------------------------------
function contains(sampling::Intervals, l::LookupArray, sel::Contains; err=_True())
    _locus_checkbounds(locus(l), l, sel) || return _boundserror_or_nothing(err)
    contains(order(l), span(l), sampling, locus(l), l, sel; err)
end
# Regular Intervals ---------------------------
function contains(o::Ordered, span::Regular, ::Intervals, locus::Locus, l::LookupArray, sel::Contains;
    err=_True()
)
    v = val(sel)
    i = _searchfunc(locus, o)(l, v)
    return check_regular_contains(span, locus, l, v, i, err)
end
function contains(o::Ordered, span::Regular, ::Intervals, locus::Center, l::LookupArray, sel::Contains;
    err=_True()
)
    v = val(sel) + abs(val(span)) / 2
    i = _searchfunc(locus, o)(l, v)
    return check_regular_contains(span, locus, l, v, i, err)
end

function check_regular_contains(span::Span, locus::Locus, l::LookupArray, v, i, err)
    absstep = abs(val(span))
    if (parent(l) isa AbstractRange) || _lt(locus)(v, l[i] + absstep)
        return i
    else
        return _notcontained_or_nothing(err, v)
    end
end

# Explicit Intervals ---------------------------
function contains(
    o::Ordered, span::Explicit, ::Intervals, locus, l::LookupArray, sel::Contains;
    err=_True()
)
    v = val(sel)
    searchfunc = _searchfunc(_Upper(), o)
    i = searchfunc(view(val(span), 1, :), v; order=ordering(o), lt=_lt(locus))
    if i === 0 || val(span)[2, i] < v
        return _notcontained_or_nothing(err, v)
    else
        return i
    end
end
# Irregular Intervals -------------------------
function contains(
    o::Ordered, span::Irregular, ::Intervals, locus::Locus, l::LookupArray, sel::Contains;
    err=_True()
)
    return _searchfunc(locus, o)(l, val(sel))
end
function contains(
    o::Ordered, span::Irregular, ::Intervals, locus::Center, l::LookupArray, sel::Contains;
    err=_True()
)
    _order_lt(::ForwardOrdered) = (<)
    _order_lt(::ReverseOrdered) = (<=)

    v = val(sel)
    i = searchsortedfirst(l, v)
    i = if i <= firstindex(l)
        firstindex(l)
    elseif i > lastindex(l)
        lastindex(l)
    else
        interval = abs(l[i] - l[i - 1])
        distance = abs(l[i] - v)
        _order_lt(o)(interval / 2, distance) ? i - 1 : i
    end
    return i
end

_boundserror_or_nothing(err::_True) = throw(BoundsError())
_boundserror_or_nothing(err::_False) = nothing

_notcontained_or_nothing(err::_True, selval) = _notcontainederror(selval)
_notcontained_or_nothing(err::_False, selval) = nothing

_notcontainederror(v) = throw(ArgumentError("No interval contains $v"))

_searchfunc(::Locus, ::ForwardOrdered) = searchsortedlast
_searchfunc(::End, ::ForwardOrdered) = searchsortedfirst
_searchfunc(::Locus, ::ReverseOrdered) = searchsortedfirst
_searchfunc(::End, ::ReverseOrdered) = searchsortedlast

"""
    Between <: Selector

    Between(a, b)

Depreciated: use `a..b` instead of `Between(a, b)`. Other `Interval`
objects from IntervalSets.jl, like `OpenInterval(a, b) will also work,
giving the correct open/closed boundaries.

`Between` will e removed in furture to avoid clashes with `DataFrames.Between`.

Selector that retreive all indices located between 2 values,
evaluated with `>=` for the lower value, and `<` for the upper value.
This means the same value will not be counted twice in 2 adjacent
`Between` selections.

For [`Intervals`](@ref) the whole interval must be lie between the
values. For [`Points`](@ref) the points must fall between
the values. Different [`Sampling`](@ref) types may give different
results with the same input - this is the intended behaviour.

`Between` for [`Irregular`](@ref) intervals is a little complicated. The
interval is the distance between a value and the next (for `Start` locus)
or previous (for [`End`](@ref) locus) value.

For [`Center`](@ref), we take the mid point between two index values
as the start and end of each interval. This may or may not make sense for
the values in your indes, so use `Between` with `Irregular` `Intervals(Center())`
with caution.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(5:7)))
A[X(Between(15, 25)), Y(Between(4, 6.5))]

# output

1×2 DimArray{Int64,2} with dimensions:
  X Sampled 20:10:20 ForwardOrdered Regular Points,
  Y Sampled 5:6 ForwardOrdered Regular Points
 4  5
```
"""
struct Between{T<:Union{<:AbstractVector{<:Tuple{Any,Any}},Tuple{Any,Any},Nothing}} <: Selector{T}
    val::T
end
Between(args...) = Between(args)

Base.first(sel::Between) = first(val(sel))
Base.last(sel::Between) = last(val(sel))

abstract type _Side end
struct _Upper <: _Side end
struct _Lower <: _Side end

selectindices(l::LookupArray, sel::Union{Between{<:Tuple},Interval}) = between(l, sel)
function selectindices(lookup::LookupArray, sel::Between{<:AbstractVector})
    inds = Int[]
    for v in val(sel)
        append!(inds, selectindices(lookup, rebuild(sel; val=v)))
    end
end

function between(l::LookupArray, sel::Between)
    a, b = _sorttuple(sel)
    return between(l, a..b)
end
function between(l::NoLookup, sel::Interval)
    x = intersect(sel, first(axes(l, 1))..last(axes(l, 1)))
    return ceil(Int, x.left):floor(Int, x.right)
end
between(l::LookupArray, sel::Interval) = between(sampling(l), l, sel)
# This is the main method called above
function between(sampling::Sampling, l::LookupArray, sel::Interval)
    o = order(l)
    o isa Unordered && throw(ArgumentError("Cannot use an interval or `Between` with Unordered"))
    between(sampling, o, l, sel)
end

_to_closed(T, i::Interval{:closed,:closed}) = i
_to_closed(T, i::Interval{:open,:open}) = (i.left + _minstep(T))..(i.right - _minstep(T))
_to_closed(T, i::Interval{:closed,:open}) = i.left..(i.right - _minstep(T))
_to_closed(T, i::Interval{:open,:closed}) = (i.left + _minstep(T))..i.right

# Maybe not the best way to do this but manipulating <= > etc
# inside `searchsorted` is already too complicated.
_minstep(T) = eps(T)
_minstep(T::Type{<:Integer}) = one(T)
_minstep(T::Type{<:Char}) = 1

function between(sampling::NoSampling, o::Ordered, l::LookupArray, sel::Interval)
    between(Points(), o, l, sel)
end

function between(sampling, o::Ordered, l::LookupArray, sel::Interval)
    lowerbound, upperbound = bounds(l)
    lowsel, highsel = endpoints(sel)
    a = if lowsel > upperbound
        ordered_lastindex(l) + _ordscalar(o)
    elseif lowsel < lowerbound
        ordered_firstindex(l)
    else
        _between_side(_Lower(), o, span(l), sampling, l, sel, lowsel)
    end
    b = if highsel < lowerbound
        ordered_firstindex(l) - _ordscalar(o)
    elseif highsel > upperbound
        ordered_lastindex(l)
    else
        _between_side(_Upper(), o, span(l), sampling, l, sel, highsel)
    end
    # Fix empty range values
    a, b = _maybeflipbounds(o, (a, b))
    # Fix empty range values
    if a > b
        if b < firstindex(l)
            return firstindex(l):(firstindex(l) - 1)
        elseif a > lastindex(l)
            return (lastindex(l) + 1):lastindex(l)
        end
    else
        return a:b
    end
    return a:b
end

# Points -------------------------
function _between_side(side::_Lower, o::Ordered, span, ::Points, l, sel, v)
    i = v <= bounds(l)[1] ? ordered_firstindex(l) : _searchfunc(side, o)(l, v)
    return _close_interval(side, l, sel, l[i], i)
end
function _between_side(side::_Upper, o::Ordered, span, ::Points, l, sel, v)
    i = v >= bounds(l)[2] ? ordered_lastindex(l) : _searchfunc(side, o)(l, v)
    return _close_interval(side, l, sel, l[i], i)
end

# Regular Intervals -------------------------
# Adjust the value for the lookup locus before search
function _between_side(side, o::Ordered, ::Regular, ::Intervals, l, sel, v)
    adj = _locus_adjust(side, l)
    v1 = v + adj
    s = _ordscalar(o)
    i = _searchfunc(side, o)(l, v1)
    fi, li = firstindex(l), lastindex(l)
    v2 = if i > li
        l[li] + adj
    elseif i < fi
        l[fi] + adj
    else
        l[i] - adj
    end
    return _close_interval(side, l, sel, v2, i)
end

_locus_adjust(side, l) = _locus_adjust(side, locus(l), abs(step(span(l))))
_locus_adjust(::_Lower, locus::Start, step) = zero(step)
_locus_adjust(::_Upper, locus::Start, step) = -step
_locus_adjust(::_Lower, locus::Center, step) = step/2
_locus_adjust(::_Upper, locus::Center, step) = -step/2
_locus_adjust(::_Lower, locus::End, step) = step
_locus_adjust(::_Upper, locus::End, step) = -zero(step)

# Explicit Intervals -------------------------
# Rebuild the lookup with the bounds matrix values before searching
function _between_side(side, o::Ordered, span::Explicit, ::Intervals, l, sel, v)
    boundsvec = side isa _Lower ? view(val(span), 1, :) : view(val(span), 2, :)
    l1 = rebuild(l; data=boundsvec)
    i = _searchfunc(side, o)(l1, v)
    return checkbounds(Bool, l1, i) ? (@inbounds _close_interval(side, l1, sel, l1[i], i)) : i
end

# Irregular Intervals -----------------------
#
# This works a little differently to Regular variants,
# as we have to work with unequal step sizes, calculating them
# as we find close values.
#
# Find the inteval the value falls in.
# We need to special-case Center locus for Irregular
_between_side(side, o, span::Irregular, ::Intervals, l, sel, v) = _irreg_side(side, locus(l), o, l, sel, v)

function _irreg_side(side, locus::Union{Start,End}, o, l, sel, v)
    s = _ordscalar(o)
    if v == bounds(l)[1]
        i = ordered_firstindex(l)
        v1 = v
    elseif v == bounds(l)[2]
        i = ordered_lastindex(l)
        v1 = v
    else
        # Search for the value and offset per order/locus/side
        i = _searchfunc(o)(l, v; lt=_lt(side))
        i += - s * (_locscalar(locus) + _sideshift(side))
        # Get the value on the interval edge
        v1 = if i < firstindex(l)
            _maybeflipbounds(l, bounds(l))[1]
        elseif i > lastindex(l)
            _maybeflipbounds(l, bounds(l))[2]
        elseif side isa _Lower && locus isa End
            l[i-s]
        elseif side isa _Upper && locus isa Start
            l[i+s]
        else
            l[i]
        end
    end
    return _close_interval(side, l, sel, v1, i)
end
function _irreg_side(side, locus::Center, o, l, sel, v)
    if v == bounds(l)[1]
        i = ordered_firstindex(l)
        v1 = v
    elseif v == bounds(l)[2]
        i = ordered_lastindex(l)
        v1 = v
    else
        r = _ordscalar(o)
        sh = _sideshift(side)
        i = _searchfunc(o)(l, v; lt=_lt(side))
        (i - r < firstindex(l) ||  i - r > lastindex(l)) && return i
        half_step = abs(l[i] - l[i-r]) / 2
        distance = abs(l[i] - v)
        # Use the right less than </<= to match interval bounds
        i = if _lt(side)(distance, half_step)
            i - sh * r
        else
            i - (1 + sh) * r
        end
        shift = side isa _Lower ? -half_step : half_step
        v1 = l[i] + shift
    end
    return _close_interval(side, l, sel, v1, i)
end

_close_interval(side, l, sel, v, i) = i
function _close_interval(side::_Lower, l, sel::Interval{:open,<:Any}, v, i)
    v in sel ? i : i + _ordscalar(l)
end
function _close_interval(side::_Upper, l, sel::Interval{<:Any,:open}, v, i)
    v in sel ? i : i - _ordscalar(l)
end

_locscalar(::Start) = 1
_locscalar(::End) = 0
_sideshift(::_Lower) = -1
_sideshift(::_Upper) = 1
_ordscalar(l) = _ordscalar(order(l))
_ordscalar(::ForwardOrdered) = 1
_ordscalar(::ReverseOrdered) = -1

_lt(::_Lower) = (<)
_lt(::_Upper) = (<=)

_maybeflipbounds(m::LookupArray, bounds) = _maybeflipbounds(order(m), bounds)
_maybeflipbounds(o::ForwardOrdered, (a, b)) = (a, b)
_maybeflipbounds(o::ReverseOrdered, (a, b)) = (b, a)
_maybeflipbounds(o::Unordered, (a, b)) = (a, b)

"""
    Where <: Selector

    Where(f::Function)

Selector that filters a dimension lookup by any function that
accepts a single value and returns a `Bool`.

## Example

```jldoctest
using DimensionalData

A = DimArray([1 2 3; 4 5 6], (X(10:10:20), Y(19:21)))
A[X(Where(x -> x > 15)), Y(Where(x -> x in (19, 21)))]

# output

1×2 DimArray{Int64,2} with dimensions:
  X Sampled Int64[20] ForwardOrdered Regular Points,
  Y Sampled Int64[19, 21] ForwardOrdered Regular Points
 4  6
```
"""
struct Where{T} <: Selector{T}
    f::T
end

val(sel::Where) = sel.f

# Yes this is everything. `Where` doesn't need lookup specialisation
@inline function selectindices(lookup::LookupArray, sel::Where)
    [i for (i, v) in enumerate(parent(lookup)) if sel.f(v)]
end

"""
    All <: Selector

    All(selectors::Selector...)

Selector that combines the results of other selectors.
The indices used will be the union of all result sorted in ascending order.

## Example

```jldoctest
using DimensionalData, Unitful

dimz = X(10.0:20:200.0), Ti(1u"s":5u"s":100u"s")
A = DimArray((1:10) * (1:20)', dimz)
A[X=All(At(10.0), At(50.0)), Ti=All(1u"s"..10u"s", 90u"s"..100u"s")]

# output

2×4 DimArray{Int64,2} with dimensions:
  X Sampled Float64[10.0, 50.0] ForwardOrdered Regular Points,
  Ti Sampled Quantity{Int64, 𝐓, Unitful.FreeUnits{(s,), 𝐓, nothing}}[1 s, 6 s, 91 s, 96 s] ForwardOrdered Regular Points
 1  2  19  20
 3  6  57  60
```
"""
struct All{S<:Tuple{Vararg{<:SelectorOrInterval}}} <: Selector{S}
    selectors::S
end
All(args::SelectorOrInterval...) = All(args)

@inline function selectindices(lookup::LookupArray, sel::All)
    results = map(s -> selectindices(lookup, s), sel.selectors)
    sort!(union(results...))
end

# selectindices ==========================================================================

"""
    selectindices(lookups, selectors)

Converts [`Selector`](@ref) to regular indices.
"""
function selectindices end
@inline selectindices(lookups::LookupArrayTuple, s1, ss...) = selectindices(lookups, (s1, ss...))
@inline selectindices(lookups::LookupArrayTuple, selectors::Tuple) =
    map((l, s) -> selectindices(l, s), lookups, selectors)
@inline selectindices(lookups::LookupArrayTuple, selectors::Tuple{}) = ()
# @inline selectindices(dim::LookupArray, sel::Val) = selectindices(val(dim), At(sel))
# Standard indices are just returned.
@inline selectindices(::LookupArray, sel::StandardIndices) = sel
@inline function selectindices(l::LookupArray, sel)
    selstr = sprint(show, sel)
    throw(ArgumentError("Invalid index `$selstr`. Did you mean `At($selstr)`? Use stardard indices, `Selector`s, or `Val` for compile-time `At`."))
end
# Vectors are mapped
@inline selectindices(lookup::LookupArray, sel::Selector{<:AbstractVector}) =
    [selectindices(lookup, rebuild(sel; val=v)) for v in val(sel)]

# Unaligned LookupArray ------------------------------------------

# select_unalligned_indices is callled directly from dims2indices

# We use the transformation from the first Transformed dim.
# In practice the others could be empty.
@inline function select_unalligned_indices(lookups::LookupArrayTuple, sel::Tuple{<:Selector,Vararg{<:Selector}})
    coords = [map(val, sel)...]
    transformed = transformfunc(lookups[1])(coords)
    map(_transform2int, sel, transformed)
end

_transform2int(::Near, x) = round(Int, x)
_transform2int(sel::At, x) = _transform2int(sel::At, x, atol(sel))
_transform2int(::At, x, atol::Nothing) = convert(Int, x)
function _transform2int(::At, x, atol)
    i = round(Int, x)
    abs(x - i) <= atol ? i : _transform_notfound(x)
end

@noinline _transform_notfound(x) = throw(ArgumentError("$x not found in Transformed lookups"))

# Shared utils ============================================================================

# Return an inbounds index
_inbounds(is::Tuple, lookup::LookupArray) = map(i -> _inbounds(i, lookup), is)
function _inbounds(i::Int, lookup::LookupArray)
    if i > lastindex(lookup)
        lastindex(lookup)
    elseif i <= firstindex(lookup)
        firstindex(lookup)
    else
        i
    end
end

_sorttuple(sel::Between) = _sorttuple(val(sel))
_sorttuple((a, b)) = a < b ? (a, b) : (b, a)

_lt(::Locus) = (<)
_lt(::End) = (<=)
_gt(::Locus) = (>=)
_gt(::End) = (>)

_locus_checkbounds(loc, lookup::LookupArray, sel::Selector) =  _locus_checkbounds(loc, bounds(lookup), val(sel))
_locus_checkbounds(loc, (l, h)::Tuple, v) = !(_lt(loc)(v, l) || _gt(loc)(v, h))

_searchfunc(::ForwardOrdered) = searchsortedfirst
_searchfunc(::ReverseOrdered) = searchsortedlast

_searchfunc(::_Lower, ::ForwardOrdered) = searchsortedfirst
_searchfunc(::_Lower, ::ReverseOrdered) = searchsortedlast
_searchfunc(::_Upper, ::ForwardOrdered) = searchsortedlast
_searchfunc(::_Upper, ::ReverseOrdered) = searchsortedfirst

hasselection(lookup::LookupArray, sel::At) = at(lookup, sel; err=_False()) === nothing ? false : true
hasselection(lookup::LookupArray, sel::Contains) = contains(lookup, sel; err=_False()) === nothing ? false : true
# Near and Between only fail on Unordered
# Otherwise Near returns the nearest index, and Between an empty range
hasselection(lookup::LookupArray, selnear::Near) = order(lookup) isa Unordered ? false : true
hasselection(lookup::LookupArray, selnear::Union{Interval,Between}) = order(lookup) isa Unordered ? false : true


text:
function MCTS.RandomBandit(mdp::LightDark2D; kwargs...)
    RandomBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        kwargs...)
end
function MCTS.DPWBandit(mdp::LightDark2D; kwargs...)
    DPWBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        exploration_constant=0.75,
        k_action=10.0,
        alpha_action=0.25,
        kwargs...)
end
function MCTS.CBTSBandit(mdp::LightDark2D; kwargs...)
    CBTSBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        exploration_constant=1.0,
        A_max=20,
        n_proposes=100,
        log_length_scale=0.0,
        log_signal_sigma=0.0,
        log_obs_noise=-1.0,
        action_dims=2,
        n_sig=2.0,  #number of standard deviations for GP-UCB
        kwargs...)
end
function MCTS.CBTSDPWBandit(mdp::LightDark2D; kwargs...)
    CBTSDPWBandit(;
        enable_action_pw=true,
        check_repeat_action=true,
        exploration_constant=1.0,
        k_action=5.0,
        alpha_action=0.5,
        n_proposes=100,
        log_length_scale=0.0,
        log_signal_sigma=0.0,
        log_obs_noise=-1.0,
        fit_qs=false,
        action_dims=2,
        n_sig=2.0,  #number of standard deviations for GP-UCB
        kwargs...)
end
function MCTS.ModularSolver(mdp::LightDark2D, b::ModularBandit, seed=0;
                            kwargs...)
    rng = MersenneTwister(seed)
    ModularSolver(;
        n_iterations=100,
        bandit=b,
        depth=10,
        k_state=0.1,
        alpha_state=0.9,
        keep_tree=true,
        check_repeat_state=true,
        rng=rng,
        init_Q=0.0,
        init_N=0,
        next_action=RandomGaussianLD2(rng),
        estimate_value=RolloutEstimator(RandomGaussianLD2(rng)),
        track_best_path=true,
        kwargs...
       )
end
Base.string(::LightDark2D) = "LightDark2D"

ext: jl


text:
using Plots;                            #Use the Package Plots

println("starting ========")

gr();                                              #Set Backend to gr for the package Plots
#plotly();                                         #Set Backend to plotly for the package Plots
dataArr1 =[5,2,3,4,5,5,-6,7,4,2,4,3,5,-2,-4,-5];                    #Create an Array
dataArr2 = [1,2,3,4,5,6,4,2,1,2,4,2,1,1,1,5]

arrLength=length(dataArr1);
xAxisValues=[1:arrLength]

diagramObj = plot(xAxisValues,dataArr1,label="dataArr1", title="Hello!",seriestype = [:bar])
diagramObj = plot!(xAxisValues,dataArr2,label="dataArr2", seriestype = [:line])
display(diagramObj)
println("end ========")

max_stars_repo_name: RainerNeu/juliaExamplesWithJuliaPro

max_issues_repo_name: RainerNeu/juliaExamplesWithJuliaPro

max_forks_repo_name: RainerNeu/juliaExamplesWithJuliaPro


text:
"""
    Field(grid, name, data_type, full_field_on_root)

Creates a distributed array of type `data_type` with dimensions specified by
`grid`. If `full_field_on_root` is true, then the root node allocates memory
the entire array, although in normal operation it only utilizes the portion of
memory that it is assigned to update. Obviously, care must be taken to not use
this option for large simulations.
"""
struct Field{ET, G, N}
    field_name::String
    field_values::Array{ET, N}
    grid::G
    full_field_on_root::Bool
end

function Field(grid::G, name::String, data_type::DataType=FT,
               full_field_on_root::Bool=false
              ) where {G <: UniformGrid{N, FT}} where {N, FT}
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)

    if full_field_on_root && rank == 0
        field_values = Array{data_type, N}(undef, grid.global_num_cells .+
                                           (2*grid.num_guard_cells + 1))
    else
        field_values = Array{data_type, N}(undef, grid.num_cells .+
                                           (2*grid.num_guard_cells + 1))
    end

    return Field{data_type, G, N}(name, field_values, grid, full_field_on_root)
end

Base.eltype(f::Field{ET, G}) where {ET, G} = ET

function Base.show(io::IO, f::Field{ET, G, N}) where {ET, G, N}
    print(io, "Field{$ET, $G, $N}($(f.grid), \"$(f.field_name)\")")
end

function Base.getindex(f::Field{ET, G, N}, indices...) where {ET, G, N}
    @assert length(indices) == N

    return getindex(f.field_values, (indices .+ f.grid.num_guard_cells)...)
end

"""
    communicate_halo!(f::Field, comm)

Communicates the values of each processors guard cells to its neighbors.
"""
function communicate_halo!(f::Field, comm)
    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    size = MPI.Comm_size(comm)

    ngc = f.grid.num_guard_cells
    nc  = f.grid.num_cells

    indices = compute_proc_indices(f.grid.proc_dims, rank)
    onetos = map(Base.OneTo, nc .+ (2*ngc + 1))

    for (i, index) in enumerate(indices)
        lower_indices = tuple(indices[1:i-1]..., indices[i] - 1, indices[i+1:end]...)
        lower_proc_rank = compute_proc_rank(f.grid.proc_dims, lower_indices)
        upper_indices = tuple(indices[1:i-1]..., indices[i] + 1, indices[i+1:end]...)
        upper_proc_rank = compute_proc_rank(f.grid.proc_dims, upper_indices)

        sendrange1 = tuple(onetos[1:i-1]..., ngc + 1:2*ngc + 1,                 onetos[i+1:end]...)
        recvrange1 = tuple(onetos[1:i-1]..., ngc + nc[i] + 1:nc[i] + 2*ngc + 1, onetos[i+1:end]...)

        sview1 = view(f.field_values, sendrange1...)
        rview1 = view(f.field_values, recvrange1...)

        MPI.isend(sview1, lower_proc_rank, 10*i + 1, comm)
        data1, stat1 = MPI.recv(upper_proc_rank, 10*i + 1, comm)
        rview1 .= data1

        sendrange2 = tuple(onetos[1:i-1]..., nc[i] + 1:ngc + nc[i], onetos[i+1:end]...)
        recvrange2 = tuple(onetos[1:i-1]..., 1:ngc,             onetos[i+1:end]...)

        sview2 = view(f.field_values, sendrange2...)
        rview2 = view(f.field_values, recvrange2...)

        MPI.isend(sview2, upper_proc_rank, 10*i + 2, comm)
        data2, stat2 = MPI.recv(lower_proc_rank, 10*i + 2, comm)
        rview2 .= data2
    end
end
"""
    gather_field_on_root!(f::Field)

Gathers the entire field `f` on the root node. Can only be used if `f` was created
with `full_field_on_root = true`.
"""
function gather_field_on_root!(f::Field)
    @assert f.full_field_on_root == true

    comm = MPI.COMM_WORLD
    rank = MPI.Comm_rank(comm)
    size = MPI.Comm_size(comm)

    if rank == 0
        for proc in 1:size-1
            lbs, ubs = compute_proc_bounds(f.grid.global_num_cells, f.grid.proc_dims, proc)
            data, stat = MPI.recv(proc, 1, comm)
            indices = map((lb, ub) -> lb + f.grid.num_guard_cells:ub + f.grid.num_guard_cells, lbs, ubs)
            view(f.field_values, indices...) .= data
        end
    else
        ngc = f.grid.num_guard_cells
        indices = map(nc -> ngc + 1:ngc + nc, f.grid.num_cells)
        MPI.send(view(f.field_values, indices...), 0, 1, comm)
    end
end


text:
using StaticArrays

struct QuadLimbDark{T,VT <: StaticVector{3,T}} <: AbstractLimbDark
    n_max::Int
    u_n::VT
    g_n::VT
    norm::T
end

@doc raw"""
    QuadLimbDark(u::AbstractVector)

A specialized implementation of [`PolynomialLimbDark`](@ref) with a maximum of two terms (quadratic form). This has a completely closed-form solution without any numerical integration. This means there are no intermediate allocations and reduced numerical error.

# Mathematical form
```math
I(\mu) \propto 1 - u_1(1-\mu) - u_2(1-\mu)^2
```

!!! warning "Higher-order terms"
    Higher-order terms will be *ignored*; no error will be thrown

# Examples

```jldoctest quad
ld = QuadLimbDark(Float64[]) # constant term only

b = [0, 1, 2] # impact parameter
r = 0.01 # radius ratio
ld.(b, r)

# output
3-element Vector{Float64}:
 0.9999
 0.9999501061035608
 1.0
```

```jldoctest quad
ld = QuadLimbDark([0.4, 0.26]) # max two terms
ld.(b, r)

# output
3-element Vector{Float64}:
 0.9998785437247428
 0.999974726693709
 1.0
```

# References

See references for [`PolynomialLimbDark`](@ref)
"""
function QuadLimbDark(u::AbstractVector{T}) where T
    n_max = length(u)
    # add constant u_0 term
    if n_max == 0
        u_n = SA[-one(T), zero(T), zero(T)]
    elseif n_max == 1
        u_n = SA[-one(T), u[begin], zero(T)]
    else
        u_n = SA[-one(T), u[begin], u[begin + 1]]
    end
    n_max > 2 && @warn "Higher-order terms will be ignored"
    n_max = min(2, n_max)

    # get Green's basis coefficients
    g_n = compute_gn(u_n)

    # calculate flux normalization factor, which only depends on first two terms
    norm = inv(π * (g_n[begin] + 2/3 * g_n[begin + 1]))

    return QuadLimbDark(n_max, u_n, g_n, norm)
end

function compute(ld::QuadLimbDark, b::S, r) where S
    T = float(S)
    ## check for trivial cases
    if b ≥ 1 + r || iszero(r)
        # completely unobscured
        return one(T)
    elseif r ≥ 1 + b
        # completely obscured
        return zero(T)
    end

    r2 = r^2
    b2 = b^2

    if iszero(b)
        onemr2 = 1 - r2
        sqrt1mr2 = sqrt(onemr2)

        # annular ellipse
        flux = ld.g_n[begin] * onemr2 + 2/3 * ld.g_n[begin + 1] * sqrt1mr2^3
        if ld.n_max > 1
            flux -= ld.g_n[begin + 2] * 2 * r2 * onemr2
        end
        return flux * π * ld.norm
    end

    # take a moment to calculate values used repeatedly in the
    # numerical routines ahead
    onembmr2 = (r + 1 - b) * (1 - r + b)
    onembmr2inv = inv(onembmr2)
    sqonembmr2 = sqrt(onembmr2)
    br = b * r
    fourbr = 4 * br
    fourbrinv = inv(fourbr)
    sqbr = sqrt(br)
    sqbrinv = inv(sqbr)
    onembpr2 = (1 - r - b) * (1 + b + r)
    sqarea = sqarea_triangle(one(T), r, b)
    k2 = max(0, onembpr2 * fourbrinv + 1)
    k = sqrt(k2)
    onemr2mb2 = (1 - r) * (1 + r) - b2
    onemr2pb2 = (1 - r) * (1 + r) + b2
    if k2 > 1
        if k2 > 2
            kc2 = 1 - inv(k2)
        else
            kc2 = onembpr2 * onembmr2inv
        end
    else
        if k2 > 0.5
            kc2 = (r - 1 + b) * (b + r + 1) * fourbrinv
        else
            kc2 = 1 - k2
        end
    end
    kc = sqrt(kc2)

    ## Compute uniform term
    s0, kap0, kite_area2, kck = compute_uniform(b, r; b2, r2, sqarea, fourbrinv)

    flux = ld.g_n[begin] * s0

    if ld.n_max == 0
        return flux * ld.norm
    end

    ## compute linear term
    s1, Eofk, Em1mKdm = compute_linear(b, r; k2, kc, kc2, r2, b2, br, fourbr, sqbr, onembmr2, onembpr2, onembmr2inv)

    flux += ld.g_n[begin + 1] * s1

    if ld.n_max == 1
        return flux * ld.norm
    end

    ## calculate quadratic term
    s2 = compute_quadratic(b, r; s0, r2, b2, kap0, kite_area2, k2)
    flux += ld.g_n[begin + 2] * s2

    return flux * ld.norm
end

function compute_gn(u_n::StaticVector{3,T}) where T
    g_0 = 1 - u_n[begin + 1] - 1.5 * u_n[begin + 2]
    g_1 = u_n[begin + 1] + 2 * u_n[begin + 2]
    g_2 = -0.25 * u_n[begin + 2]
    return SA[g_0, g_1, g_2]
end


text:
#Structures

density,fourier,intensity = createSphericalHarmonicsStructure("$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb", LMAX, KMAX, float(KMAX))
densityCube,fourierCube,intensityCube = createCubicStructure("$(ENV["DETERMINATION_DATA"])/structures/crambin.pdb", 2*KMAX+1, float(KMAX))
# @test_approx_eq_eps similarity(getCube(density), densityCube) 1.0 1e-1
# @test_approx_eq_eps similarity(getCube(intensity), intensityCube) 1.0 1e-2


text:
# game status

struct TransP
    hs::UInt64
    best::Move
    depth::Float64
    flags::Int
    val::Int
end

mutable struct GameStatus
    bookfile::String
    usebook::Bool
    canponder::Bool
    hashsize::UInt32
    board::Board
    lastPV::Array{Move,1}
    lastPVLength::Int
    whiteHeuristics::Array{Int,2}
    blackHeuristics::Array{Int,2}
    inodes::Int32
    nsStart::Int
    moveBufLen::Array{Int,1}
    moveBuf::Array{Move,1}
    triangularLength::Array{Int,1}
    triangularArray::Array{Move,2}
    side::Int
    timedout::Bool
    followpv::Bool
    allownull::Bool
    pvmovesfound::Int
    depth::Float64
    AttackTableNonSlide::Array{BitBoard,2}
    FillRank::Array{BitBoard,1}
    FillFile::Array{BitBoard,1}
    fukyBitsW::BitBoard
    fukyBitsB::BitBoard
    keBitsW::BitBoard
    keBitsB::BitBoard
    MovableKoma::Array{BitBoard,2}
    SenteJin::BitBoard
    GoteJin::BitBoard
    SenteOthers::BitBoard
    GoteOthers::BitBoard
    MoveBeginIndex::Int
    tt::Dict{UInt64,TransP}
    ett::Dict{UInt64,Int64}
    remainTime::Int64
    maxThinkingTime::Int64
    GameStatus() = new()
end

function InitGS()
    gs = GameStatus()
    gs.bookfile = "./Joseki.db"
    gs.usebook  = true
    gs.canponder= true
    gs.hashsize = UInt32(256 * 1024 * 1024)
    gs.board = Board()
    gs.lastPV = [Move(0,0,0,0,0,0) for x = 1:MaxPly]::Array{Move,1}
    gs.lastPVLength = 0
    gs.whiteHeuristics = [0 for x = 1:NumSQ, y = 1:NumSQ]::Array{Int,2}
    gs.blackHeuristics = [0 for x = 1:NumSQ, y = 1:NumSQ]::Array{Int,2}
    gs.inodes = 0
    gs.nsStart = 0
    gs.moveBufLen = [0 for x = 1:MaxPly]::Array{Int,1}
    gs.moveBuf = [Move(0,0,0,0,0,0) for x = 1:MaxMoves]
    gs.triangularLength = [0 for x = 1:MaxPly]::Array{Int,1}
    gs.triangularArray  = [Move(0,0,0,0,0,0) for x = 1:MaxPly, y = 1:MaxPly]::Array{Move,2}
    gs.side = SENTE # 仮に先手として、後で書き換える
    gs.timedout = false
    gs.followpv = true
    gs.allownull = true
    gs.pvmovesfound = 0
    gs.depth = 0.0
    gs.AttackTableNonSlide = [UInt128(0) for piece=MJFU:MJGORY, sq=A9:I1]::Array{BitBoard,2}
    gs.FillRank = [UInt128(0) for r=1:9]::Array{BitBoard,1}
    gs.FillFile = [UInt128(0) for f=1:9]::Array{BitBoard,1}
    gs.fukyBitsW = UInt128(0)::BitBoard
    gs.fukyBitsB = UInt128(0)::BitBoard
    gs.keBitsW   = UInt128(0)::BitBoard
    gs.keBitsB   = UInt128(0)::BitBoard
    gs.MovableKoma = [UInt128(0) UInt128(0) UInt128(0);
                      UInt128(0) UInt128(0) UInt128(0)]::Array{BitBoard,2}
    gs.SenteJin    = UInt128(0)::BitBoard
    gs.GoteJin     = UInt128(0)::BitBoard
    gs.SenteOthers = UInt128(0)::BitBoard
    gs.GoteOthers  = UInt128(0)::BitBoard

    gs.MoveBeginIndex = 0
    gs.tt = Dict{UInt64,TransP}()
    gs.ett = Dict{UInt64,Int64}()

    gs.remainTime = 0 #rewrite by Main.jl
    gs.maxThinkingTime = MAXTHINKINGTIME

    sizehint!(gs.tt,65536*16)
    sizehint!(gs.ett,65536*16)

    #println("$gs")
    gs
end

max_stars_repo_name: kimrin/MechajyoForWindows

max_issues_repo_name: kimrin/MechajyoForWindows

max_forks_repo_name: kimrin/MechajyoForWindows


text:
## utilities for working with Reactive.jl signals
## NOTE: foldp and foldr in Reactive.jl are as in e.g. Haskell
##  https://wiki.haskell.org/Fold

"""
Signal value change function to be used with foldp;
use with (Bool, signal_t=0) tuple as initial fold value
"""
function fsigchange(prev, x)
  prev_val = prev[2]
  # x becomes previous value in the next call
  return (x != prev_val, x)
end

"Bool change signal, true when input signal changes"
function schange{T}(s_inp::Signal{T})
  return map(s -> s[1],
             foldp(fsigchange, (false, s_inp.value), s_inp),
             typ = Bool)
end

"""
Buffer for storing previous signal values to be used with foldp when
indicators are calculated based on signal history.

**IMPORTANT**: Initial value supplied to `foldp` determines buffer window
size, i.e. how many past signal values are retained (rolling window
size). In the case of e.g. SMA that would be moving average window.
Specifying initial value may be tricky: see `test/signals.jl`.

In-place modifies `buffer` argument and returns updated one.
"""
function sighistbuffer!(buffer, valnew)
  push!(buffer, valnew) # add new value to the end of buffer
  shift!(buffer) # remove the earliest element
  return buffer
end

"""
Initialization of `nbuff`-size float-elements buffer
with NaNs and last element `xinit`.
"""
initbuff(nbuff::Int64, xinit::Float64) = [fill(NaN, nbuff-1); xinit]


text:
#plotbuilders.jl: Defines/loads template for building plots
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#
LBLAX_AMPV = "Amplitude [V]"
LBLAX_FREQ = "Frequency [Hz]"
LBLAX_ZVA = "Impedance [V/A]"
LBLAX_MAGZVA = "|Impedance| [V/A]"
LBLAX_POTENTIAL = "Potential [V]"
LBLAX_VOLT = "Voltage [V]"
LBLAX_CURmA = "Current [mA]"
LBLAX_TIME = "Time [s]"

wtrace = 3
line_p1 = cons(:a, line = set(style=:solid, color=:black, width=wtrace))
line_p2 = cons(:a, line = set(style=:dash, color=colorant"red", width=wtrace))
line_Z = cons(:a, line = set(style=:solid, color=:blue, width=wtrace))
line_gnd = cons(:a, line = set(style=:solid, color=colorant"green", width=wtrace))

#Show options
#optleg = EasyPlot.ShowOptions(dim=set(w=350, h=160))
#optnoleg = EasyPlot.ShowOptions(dim=set(w=270, h=160))
optnoleg = EasyPlot.ShowOptions(dim=set(w=480, h=300))

#Plot builder objects:
const plot_builders = Dict{Symbol, EasyPlot.EasyPlotBuilder}()
include("plotbuilders/supplyimpedance.jl")
include("plotbuilders/transientplot.jl")

#==Postproc functions
===============================================================================#
function adjust_inspectdr(gtkplot, tgtmajor)
	DFLTGRID = InspectDR.GridRect(vmajor=true, vminor=true, hmajor=false, hminor=false)
	mplot = gtkplot.src
	mplot.layout[:valloc_title] = 0 #Assume not Multiplot title
	for sp in mplot.subplots
#		sp.layout[:valloc_data] = 300; sp.layout[:halloc_data] = 480
#		sp.layout[:valloc_data] = 150; sp.layout[:halloc_data] = 240
#		sp.layout[:valloc_data] = 100; sp.layout[:halloc_data] = 150
#		sp.layout[:valloc_data] = 150; sp.layout[:halloc_data] = 200
#		sp.layout[:halloc_legend] = WIDTH_LEGEND
		for stripi in sp.strips
			stripi.grid = DFLTGRID
			scale = isa(stripi.yscale, InspectDR.LogScale) ? :log : :lin
			stripi.yscale = InspectDR.AxisScale(scale, tgtmajor=tgtmajor, tgtminor=1)
		end
	end
end

const plotdisplay = Dict{Symbol, EasyPlot.GUIDisplay}(
	:dflt => EasyPlot.GUIDisplay(:InspectDR, postproc=(gtkplot)->adjust_inspectdr(gtkplot, 4)),
	:finey => EasyPlot.GUIDisplay(:InspectDR, postproc=(gtkplot)->adjust_inspectdr(gtkplot, 6)),
)

#Last line


text:
module ParameterJuMP

using SparseArrays

import MutableArithmetics
const MA = MutableArithmetics

using JuMP
export index

export ParameterRef, Param, all_parameters, parametrized_dual_objective_value

# types
# ------------------------------------------------------------------------------
struct ParameterRef <: AbstractVariableRef
    ind::Int64 # local reference
    model::Model
end

# Reference to a constraint in which the parameter has coefficient coef
struct ParametrizedConstraintRef{C}
    cref::ConstraintRef{Model,C}
    coef::Float64
end

const CtrRef{F,S} = ConstraintRef{Model,MOI.ConstraintIndex{F,S},ScalarShape}
const SAF = MOI.ScalarAffineFunction{Float64}

mutable struct _ParameterData
    inds::Vector{Int64}
    next_ind::Int64
    current_values::Vector{Float64}
    future_values::Vector{Float64}
    # Whether the JuMP model is in sync with the value of the parameters
    sync::Bool
    # constraints where it is a RHS
    constraints_map::Dict{Int64,Vector{ParametrizedConstraintRef}}
    parameters_map_saf_in_eq::Dict{
        CtrRef{SAF,MOI.EqualTo{Float64}},
        GenericAffExpr{Float64,ParameterRef},
    }
    parameters_map_saf_in_le::Dict{
        CtrRef{SAF,MOI.LessThan{Float64}},
        GenericAffExpr{Float64,ParameterRef},
    }
    parameters_map_saf_in_ge::Dict{
        CtrRef{SAF,MOI.GreaterThan{Float64}},
        GenericAffExpr{Float64,ParameterRef},
    }
    names::Dict{ParameterRef,String}
    lazy::Bool
    no_duals::Bool
    dual_values::Vector{Float64}

    function _ParameterData()
        return new(
            Int64[],
            1,
            Float64[],
            Float64[],
            true,
            Dict{Int64,Vector{ConstraintRef}}(),
            Dict{CtrRef{SAF,MOI.EqualTo{Float64}},GenericAffExpr{Float64,ParameterRef}}(),
            Dict{CtrRef{SAF,MOI.LessThan{Float64}},GenericAffExpr{Float64,ParameterRef}}(),
            Dict{CtrRef{SAF,MOI.GreaterThan{Float64}},GenericAffExpr{Float64,ParameterRef}}(),
            Dict{ParameterRef,String}(),
            false,
            false,
            Float64[],
        )
    end
end

no_duals(data::_ParameterData) = data.no_duals
no_duals(model::Model) = no_duals(_getparamdata(model))

set_no_duals(model::Model) = set_no_duals(_getparamdata(model))
function set_no_duals(data::_ParameterData)
    if isempty(data.current_values)
        data.no_duals = true
    elseif no_duals(data)
        @warn "No duals mode is already activated"
    else
        error("Parameter JuMP's no duals mode can only be activated in empty models.")
    end
end

"""
    index(p::ParameterRef)::Int64

Return the internal index of the parameter `p`.
"""
JuMP.index(p::ParameterRef) = p.ind

lazy_duals(data::_ParameterData) = data.lazy
lazy_duals(model::Model) = lazy_duals(_getparamdata(model))

set_lazy_duals(model::Model) = set_lazy_duals(_getparamdata(model))
function set_lazy_duals(data::_ParameterData)
    if isempty(data.current_values)
        data.lazy = true
    elseif lazy_duals(data)
        @warn "Lazy mode is already activated"
    else
        error("Parameter JuMP's lazy mode can only be activated in empty models.")
    end
end

set_not_lazy_duals(model::Model) = set_not_lazy_duals(_getparamdata(model))
function set_not_lazy_duals(data::_ParameterData)
    if isempty(data.current_values)
        data.lazy = false
    elseif !lazy_duals(data)
        @warn "Lazy mode is already de-activated"
    else
        error("Parameter JuMP's lazy mode can only be de-activated in empty models.")
    end
end

function _get_param_dict(data::_ParameterData, ::Type{MOI.EqualTo{Float64}})
    return data.parameters_map_saf_in_eq
end

function _get_param_dict(data::_ParameterData, ::Type{MOI.LessThan{Float64}})
    return data.parameters_map_saf_in_le
end

function _get_param_dict(data::_ParameterData, ::Type{MOI.GreaterThan{Float64}})
    return data.parameters_map_saf_in_ge
end

_getmodel(p::ParameterRef) = p.model

function _getparamdata(p::ParameterRef)::_ParameterData
    return _getparamdata(_getmodel(p))::_ParameterData
end

function _getparamdata(model::Model)::_ParameterData
    params = get(model.ext, :ParameterJuMP, nothing)
    if params !== nothing
        return params
    end
    return enable_parameters(model)
end

"""
    is_sync(model::Model)

Test if the JuMP Model is updated to the latest value of all parameters.
"""
is_sync(data::_ParameterData) = data.sync
is_sync(model::Model) = is_sync(_getparamdata(model))

function _add_parameter(model::Model, val::Real)
    params = _getparamdata(model)::_ParameterData
    ind = params.next_ind
    params.next_ind += 1
    push!(params.inds, ind)
    push!(params.current_values, 0.0)
    push!(params.future_values, val)
    if !iszero(val)
        params.sync = false
    end
    push!(params.dual_values, NaN)
    params.constraints_map[ind] = ParametrizedConstraintRef[]
    return ParameterRef(ind, model)
end

_addsizehint!(vec::Vector, len::Integer) = sizehint!(vec, len+length(vec))

"""
    all_parameters(model::AbstractModel)::Vector{ParameterRef}

Returns a list of all parameters currently in the model. The parameters are
ordered by creation time.
"""
function all_parameters(model::AbstractModel)
    data = _getparamdata(model)
    return ParameterRef[ParameterRef(ind, model) for ind in data.inds]
end

# solve
# ------------------------------------------------------------------------------

"""
    enable_parameters(m::Model)

Enables a `Model` to handle `Parameters`.
"""
function enable_parameters(m::Model)
    if haskey(m.ext, :params)
        error("Model already has parameter enabled")
    end
    set_optimize_hook(m, parameter_optimizehook)
    return _initialize_parameter_data(m)
end

function _initialize_parameter_data(model::Model)
    return model.ext[:ParameterJuMP] = _ParameterData()
end

# The constraint has been deleted. We cannot keep `dict` in sync with deletions
# as we return a `ConstraintRef`, not a custom type when the user create a
# parametrized constraint.
function _update_dicts_with_deletion(data::_ParameterData, ::Type{S}) where {S}
    dict = _get_param_dict(data, S)
    to_delete = eltype(keys(dict))[]
    for cref in keys(dict)
        if !is_valid(cref.model, cref)
            push!(to_delete, cref)
        end
    end
    for cref in to_delete
        delete!(dict, cref)
    end
end

function _update_dicts_with_deletion(data::_ParameterData)
    _update_dicts_with_deletion(data, MOI.EqualTo{Float64})
    _update_dicts_with_deletion(data, MOI.LessThan{Float64})
    _update_dicts_with_deletion(data, MOI.GreaterThan{Float64})
    return
end

function parameter_optimizehook(m::Model; kwargs...)
    data = _getparamdata(m)::_ParameterData
    _update_dicts_with_deletion(data)
    sync(data)
    optimize!(m::Model, ignore_optimize_hook = true, kwargs...)
    if !no_duals(data) && has_duals(m)
        _update_duals(data)
    end
    return
end

"""
    sync(model::Model)

Forces the model to update its constraints to the new values of the parameters.
"""
function sync(data::_ParameterData)
    if !is_sync(data)
        _update_constraints(data)
        data.current_values .= data.future_values
        data.sync = true
    end
    return
end
sync(model::Model) = sync(_getparamdata(model))

function parametrized_dual_objective_value(model::AbstractModel)
    params = all_parameters(model)
    linear = [param => dual(param) for param in params]
    obj = dual_objective_value(model)
    constant = obj - sum(value(term.first) * term.second for term in linear)
    return GenericAffExpr(constant, linear)
end

include("constraints.jl")
include("variable_interface.jl")
include("operators.jl")
include("mutable_arithmetics.jl")
include("macros.jl")
include("print.jl")

include("deprecations.jl")

end


text:
#Here we will examine what happens if we set threshold MSE's to continue the process until reached:

include("XOR_MD.jl")
include("Kernals.jl")
include("gaussian_process.jl")

#Set the random seed:
 #Seed for stalzer srand(1234)

"""
We will compare the effect of randomly selecting a learning rate and sigmoid
hyperparameter vs the use of Bayesian Optimization for finding the optimial values wrt the MSE.
Suppose we have limited computing time of 100000 epochs and that we have N tries to
minimise the MSE. Let us say that the learning rate is between a and b
"""

#Initialise Layers and params ==========================================

Layer_1=uniform(0,1,2,2)
Layer_2=uniform(0,1,2,1)

epochs=1000   # was 1000 below, please change back
a=0.001  #Change to 0.001
b=1

c=0.001 #Was the same as above
d=1

N=20  #dont need to change this from the file XOR_Timings as we know it converges within 6 attemps

MSE_Threshold=0.177
 #This is the MSE threshold for random to achieve

ThresholdN=10000000  (#We origionally let the random search have 10,000 attempts and at that stage we stop it)

#Curry the sigmoid functions:
#Curry the sigmoid functions:

function hyper_curry(h)
    return (x->sigmoid(x,h))
end

function hyper_curry_deriv(h)
    return (x->sigmoid_deriv(x,h))
end

Bayesian_Times=zeros(N)

Random_Times=[]

"""
Section 1 -- Run the random section one time to remove compiling timing error
"""

z

# Random Learning Rates First run to remove compiler error!! ========================================

Random_Learning_Rates=uniform(a,b,ThresholdN,1)
Random_Hyperparameters=uniform(c,d,ThresholdN,1)
Random_Mat=cat(2,Random_Learning_Rates,Random_Hyperparameters)

Random_MSE=ones(ThresholdN)

time_so_far=0
#Random_Mat conjoins Random_Learning_Rates and Random_Hyperparameters
# Random_Mat is a Nx2 matrix where Random_Mat[1,:] is the first entry
#with LR_1 and hyperparemeter 1.
for i=1:length(Random_Learning_Rates)

    tic()

    node_function=hyper_curry(Random_Mat[i,2])
    node_deriv=hyper_curry_deriv(Random_Mat[i,2])
    learning_rate=Random_Mat[i,1]
    one_net=Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]
    push!(Random_MSE,one_net)
    if one_net<MSE_Threshold
        println("Threshold value has been reached upon search ",i)

        break
    end
    time_so_far+=toc()
    println("==========================")
    println("Completed pre-run cycle ",i, " of ",length(Random_Learning_Rates))
    println("==========================")
    println("Current Time of pre-run elapsed ",time_so_far)

end

# """
# Section 2-- Run the random selection but this time timing the whole process
# """

srand(1234)

Random_Learning_Rates=uniform(a,b,ThresholdN,1)
Random_Hyperparameters=uniform(c,d,ThresholdN,1)
Random_Mat=cat(2,Random_Learning_Rates,Random_Hyperparameters)

Random_MSE=[]

q=0
non_convergance=0
convergance_val=0
for i=1:length(Random_Learning_Rates)

    tic()

    node_function=hyper_curry(Random_Mat[i,2])
    node_deriv=hyper_curry_deriv(Random_Mat[i,2])
    learning_rate=Random_Mat[i,1]

    one_net=Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]
    push!(Random_MSE,one_net)
    q+=toc()
    push!(Random_Times,q)
    if one_net<MSE_Threshold
        println("Threshold value has been reached upon search ",i)
        convergance_val=i
        println("onenet = ",one_net)
        break
    end

    if i==length(Random_Learning_Rates)
        non_convergance+=1
    end
    println("==========================")
    println("Completed cycle ",i, " of ",length(Random_Learning_Rates))
    println("==========================")
    println("Current time elapsed ",q)
end

if non_convergance==0

    println("Random Learning Rates Training Completed with convergance upon selection ", convergance_val)
else
    println("Random Learning Rates did not converge within the desired threshold value")
end

srand(123)

using PyPlot
# fig = figure("pyplot_subplot_mixed",figsize=(7,7))
ax=axes()

surf(reshape(Random_Learning_Rates[1:1:length(Random_MSE)],size(Random_MSE)),reshape(Random_Hyperparameters[1:1:length(Random_MSE)],size(Random_MSE)),Random_MSE,alpha=0.65,color="#40d5bb")
title("MSE for Randomly Selected")
xlabel("Learning Rate")
ylabel("Hyper-Parameter")
zlabel("Mean Square Error")
grid("on")
show()

# """
# Section 3 - Run Bayesian Opt one time to remove compiler problems
# """
#Initialise Layers and params ==========================================

LR_Test=linspace(a,b,50)
HP_Test=linspace(c,d,50)

#Here is the carteisan product of these written as a vector
Test=gen_points([LR_Test,HP_Test])[1]

#We first have to pick a random point to begin bayesian optimization:

#currently starts with the midpoint, possibly randomise this:
Bayesian_Points=[Test[Int(round(length(Test)/2))]]

#Here we reset the values:
# =====================================================================

#Initialise Layers and params ==========================================

LR_Test=linspace(a,b,50)
HP_Test=linspace(c,d,50)

#Here is the carteisan product of these written as a vector
Test=gen_points([LR_Test,HP_Test])[1]

#This is the break point

#We first have to pick a random point to begin bayesian optimization:

#currently starts with the midpoint, possibly randomise this:
Bayesian_Points=[Test[Int(round(length(Test)/2))]]

#Bayesian_Points is an vector of arrays where in each array first entry is LR second entry is Hyper-Parameters:

#Define hyperparemeter functions:
node_function=hyper_curry(Bayesian_Points[1][2])
node_deriv=hyper_curry_deriv(Bayesian_Points[1][2])

#Define Learning Rate:
learning_rate=Bayesian_Points[1][1]

#Run first train before Bayesian Optimization:
Bayesian_MSE=[Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]]

# =========================================================================

for k=2:N

    D=[(Bayesian_Points[i],Bayesian_MSE[i]) for i=1:length(Bayesian_Points)]
    mu, sigma, D=gaussian_process_chol(std_exp_square_ker,D,1e-6,Test)
    # println("Gaussian Process Complete","\r")
    mu=reshape(mu,length(mu));
    sigma=reshape(sigma,length(sigma))

    new_point=findmin(mu-sigma)[2]

    #Here we will need to change the number 2 to k
    Bayesian_Points=cat(1,Bayesian_Points,[Test[new_point]])

    learning_rate=Bayesian_Points[k][1]

    node_function=hyper_curry(Bayesian_Points[k][2])

    node_deriv=hyper_curry_deriv(Bayesian_Points[k][2])

    value_to_be_appended=Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]

    if value_to_be_appended !=Bayesian_MSE[k-1]
        Bayesian_MSE=cat(1,Bayesian_MSE,[value_to_be_appended])
        println("Epoch Complete")

    else
        println("Found Optimum on the ", k-1, " iteration of ", N, " iterations")
        Bayesian_Points=Bayesian_Points[1:length(Bayesian_Points)-1]
        break
    end

end

# """
# Step 4 -Run the Bayesian Code and Time
# """

srand(1234)

#We first have to pick a random point to begin bayesian optimization:

#currently starts with the midpoint, possibly randomise this:
Bayesian_Points=[Test[Int(round(length(Test)/2))]]

#Bayesian_Points is an vector of arrays where in each array first entry is LR second entry is Hyper-Parameters:

#Define hyperparemeter functions:
node_function=hyper_curry(Bayesian_Points[1][2])
node_deriv=hyper_curry_deriv(Bayesian_Points[1][2])

#Define Learning Rate:
learning_rate=Bayesian_Points[1][1]

tic()
#Run first train before Bayesian Optimization:
Bayesian_MSE=[Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]]
tq=toc()
Bayesian_Times[1]=tq

# =========================================================================

q=0 #preallocate time value at 0
for k=2:N

    tic()
    D=[(Bayesian_Points[i],Bayesian_MSE[i]) for i=1:length(Bayesian_Points)]
    mu, sigma, D=gaussian_process_chol(std_exp_square_ker,D,1e-6,Test)
    # println("Gaussian Process Complete","\r")
    mu=reshape(mu,length(mu));
    sigma=reshape(sigma,length(sigma))

    new_point=findmin(mu-sigma)[2]

    #Here we will need to change the number 2 to k
    Bayesian_Points=cat(1,Bayesian_Points,[Test[new_point]])

    learning_rate=Bayesian_Points[k][1]

    node_function=hyper_curry(Bayesian_Points[k][2])

    node_deriv=hyper_curry_deriv(Bayesian_Points[k][2])

    value_to_be_appended=Train_Neural_Net_Loop(epochs,Layer_1,Layer_2,learning_rate,node_function,node_deriv)[3]

    if value_to_be_appended !=Bayesian_MSE[k-1]
        Bayesian_MSE=cat(1,Bayesian_MSE,[value_to_be_appended])
        println("Epoch Complete")

    else
        println("Found Optimum on the ", k-1, " iteration of ", N, " iterations")
        Bayesian_Points=Bayesian_Points[1:length(Bayesian_Points)-1]
        q+=toc()
        Bayesian_Times[k]=q
        break
    end

    q+=toc()
    Bayesian_Times[k]=q

end

Bayesian_Times2=Bayesian_Times[1:length(Bayesian_MSE)]

# Bayesian Plotting =========================================================

println("Final Times for Random = ",Random_Times[end])
println("Minimum MSE for Random =",minimum(Random_MSE))
println("Final time for Bayes = ",Bayesian_Times2[end])
println("Minimum MSE for Bayes =",minimum(Bayesian_MSE))

LR=[Bayesian_Points[i][1] for i=1:length(Bayesian_Points)]
HP=[Bayesian_Points[i][2] for i=1:length(Bayesian_Points)]

planex=[0,0,1,1]
planey=[0,1,0,1]
planem=[minimum(Random_MSE),minimum(Random_MSE),minimum(Random_MSE),minimum(Random_MSE)]

using PyPlot
# fig = figure("pyplot_subplot_mixed",figsize=(7,7))
# ax=axes()
surf(LR,HP,Bayesian_MSE,alpha=0.65,color="#40d5bb")
surf(planex,planey,planem,alpha=0.3,color="#aa231f")

title("MSE BO - (6 point convergance)")
xlabel("Learning Rate")
ylabel("Hyper-Parameter")
zlabel("Mean Square Error")
grid("off")
show()

if minimum(Random_MSE)<minimum(Bayesian_MSE)
    println(" =========Results=========
        ")
    println("Random selection with ",convergance_val, " attempts achieves", minimum(Random_MSE), "beating Bayesian Optimization which obtained ",minimum(Bayesian_MSE), " in 6 attempts")

else
    if minimum(Random_MSE)<=MSE_Threshold

        println(" =========Results=========
        ")
        println("Bayesian Optimization achieves  ", minimum(Bayesian_MSE)," with 6 attempts, beating Random Search's ",minimum(Random_MSE)," which converged after ", convergance_val ," iterations.")

    else

        println(" =========Results=========
        ")
        println("Bayesian Optimization achieves  ", minimum(Bayesian_MSE)," with 6 attempts, beating Random Search's ",minimum(Random_MSE)," which was did not converge, (cutoff point was",ThresholdN, " searches.)")
    end

end

using PyPlot
# fig = figure("pyplot_subplot_mixed",figsize=(7,7))
# ax=axes()
plot(Bayesian_Times2,Bayesian_MSE,label="Bayesian Optimization",color="#40d5bb")
plot(Random_Times,Random_MSE,label="Random Grid Search",color="#aa231f",alpha=0.45)
title("Threshold MSE = 0.180")
xlabel("Time (s)")
ylabel("MSE")
legend(loc="upper right",fancybox="true")
grid("on")
show()

#0.190 = 2034 searches in time =  121.77159729299996   Minimum MSE for Random =0.18607611627885998
#0.180=  5685   searches in time  =  454.35400524799934
#BO =    seaches in time  =  7562 (s) ( for prerun )


text:
import FFTW: bfft

function interpolate(x, y₁, y₂, RCWAsol)
    sol = RCWAsol
    @assert x == sol.x "grid doesn't match, so likely incorrect"
    O₁ = rotr90(bfft(exp.(-sol.kz₁ * transpose(im * y₁)) .* sol.Õ₁, 1))
    I₂ = rotr90(bfft(exp.(-sol.kz₂ * transpose(im * y₂)) .* sol.Ĩ₂, 1))
    I₁ = rotr90(bfft(exp.( sol.kz₁ * transpose(im * y₁)) .* sol.Ĩ₁, 1))
    O₂ = rotr90(bfft(exp.( sol.kz₂ * transpose(im * y₂)) .* sol.Õ₂, 1))
    (x=x, y₁=y₁, y₂=y₂, I₁=I₁, O₁=O₁, I₂=I₂, O₂=O₂)
end


text:
struct Sequence{Outer<:Reagent,Inner<:Reagent} <: Reagent
    outer::Outer
    inner::Inner
end

then(r::Sequence, actr::Reactable) = then(r.outer, then(r.inner, actr))
then(r::Reagent, actr::Reactable) = Reactor(r, actr)
then(::Identity, actr::Reactable) = actr

Reagents.Reagent(::Commit) = Identity()
Reagents.Reagent(actr::Reactor) = actr.reagent ⨟ Reagent(actr.continuation)

hascas(::Identity) = false

struct Choice{R1<:Reagent,R2<:Reagent} <: Reagent
    r1::R1
    r2::R2
end

hascas(r::Choice) = hascas(r.r1) || hascas(r.r2)

_maysync(r::Reagent) = maysync(then(r, Commit()))
_hascas(r::Reagent) = hascas(then(r, Commit()))

function tryreact!(actr::Reactor{<:Choice}, a, rx::Reaction, offer::Union{Offer,Nothing})
    (; r1, r2) = actr.reagent
    ans1 = tryreact!(then(r1, actr.continuation), a, rx, offer)
    ans1 isa Failure || return ans1
    if offer === nothing
        if _maysync(r1) && _hascas(r2)
            # If the first branch may synchronize, and the second branch has a
            # CAS, we need to simultaneously rescind the offer *and* commit the
            # CASes.
            return Block()
        end
    end
    ans2 = tryreact!(then(r2, actr.continuation), a, rx, offer)
    ans2 isa Failure || return ans2
    if ans1 isa Retry
        return Retry()
    else
        return ans2
    end
end

struct Both{R1<:Reagent,R2<:Reagent} <: Reagent
    r1::R1
    r2::R2
end

hascas(r::Both) = hascas(r.r1) || hascas(r.r2)

function then(r::Both, actr::Reactable)
    (; r1, r2) = r
    function tee1(x)
        function tee2(y1)
            zip12(y2) = Return((y1, y2))
            return Return(x) ⨟ r2 ⨟ Computed(zip12)
        end
        return Return(x) ⨟ r1 ⨟ Computed(tee2)
    end
    return then(Computed(tee1), actr)
end

Base.:∘(inner::Reagent, outer::Sequence) = (inner ∘ outer.inner) ∘ outer.outer
Base.:∘(inner::Reagent, outer::Reagent) = Sequence(outer, inner)

# `|` could be a bit misleading since `Choice` is rather `xor`
Base.:|(r1::Reagent, r2::Reagent) = Choice(r1, r2)
Base.:&(r1::Reagent, r2::Reagent) = Both(r1, r2)

#=
Base.:+(r1::Reagent, r2::Reagent) = Choice(r1, r2)
Base.:>>(r1::Reagent, r2::Reagent) = r2 ∘ r1
Base.:*(r1::Reagent, r2::Reagent) = Both(r1, r2)
=#

struct UntilBegin{R<:Reactor} <: Reagent
    loop::R
end

hascas(r::UntilBegin) = hascas(r.loop)

struct UntilEnd{F,R<:Reactable} <: Reagent
    f::F
    continuation::R  # used only for `hascas`
end

hascas(r::UntilEnd) = hascas(r.continuation)

struct UntilBreak{T,Rx<:Reaction}
    value::T
    rx::Rx
end

then(r::Until, actr::Reactable) =
    then(UntilBegin(then(r.reagent, then(UntilEnd(r.f, actr), Commit()))), actr)

function tryreact!(
    actr::Reactor{<:UntilBegin},
    a,
    rx::Reaction,
    offer::Union{Offer,Nothing},
)
    (; loop) = actr.reagent
    while true
        ans = tryreact!(loop, a, Reaction(), nothing)
        if ans isa UntilBreak
            return tryreact!(actr.continuation, ans.value, combine(rx, ans.rx), offer)
        elseif ans isa Block
            error("`Until(f, reagent)` does not support blocking `reagent`")
        elseif ans isa Failure
            GC.safepoint()
            # TODO: backoff
        else
            ans::Nothing
        end
    end
end

function tryreact!(actr::Reactor{<:UntilEnd}, a, rx::Reaction, offer::Union{Offer,Nothing})
    (; f) = actr.reagent
    b = f(a)
    if b === nothing
        # Commit the reaction
        return tryreact!(actr.continuation, nothing, rx, offer)
    else
        return UntilBreak(b, rx)
    end
end


text:
# module SQLite_Octo.SQL

function execute(adapter::QueryBuilders.AbstractAdapter, commands::String)
    print_with_color(:bold, commands, "\n")
    df = SQLite.query(adapter.db, commands)
    print_with_color(:normal, string(df), "\n")
end


text:
# should be constriants not projection, projection is a function on the contrsrinats
abstract type Constraints end
struct BoxConstraint{T} <: Constraints
    lower::T
    upper::T
end
abstract type Scaling end
struct NoScaling <: Scaling end
scale(S::NoScaling, gx) = gx
struct ConstantPDScaling{T} <: Scaling
    S::T
end
scale(S::ConstantPDScaling, gx) = S.S \ gx

function update_xλ(prob::OptimizationProblem, S::Scaling, λ, x, gx)
    bounds = bounds(prob)
    min.(max.(bounds.lower, x .- λ .* scale(S, gx)), bounds.upper)
end
#update_xλ(p::BoxConstraint, λ, x, gx) = min.(max.(p.lower, x .- λ.*gx), p.upper)

# 5.13 in Kelley IMO
# what's it in Bertsekas?
function sufficient_decrease(fλ, fx, λ, α, x, xλ, gx)
    lhs = fλ
    # We use (22) from Bersekas (1976) to avoid dividing by λ
    rhs = fx - dot(gx, x .- xλ) * α
    lhs <= rhs
end
function sufficient_decrease(fλ, fx, λ, α, gx)
    lhs = fλ
    rhs = fx - norm(gx, 2) * α * λ
    lhs <= rhs
end

function projected_gradient(
    prob::OptimizationProblem,
    f,
    g,
    x0,
    sc::Union{Scaling,NoScaling} = NoScaling();
    itermax,
    ls_itermax,
    β = 0.5,
    α = 1e-6,
    print_trace = true,
)

    f0, g0 = f(x0), g(x0)
    if norm(x0 .- update_xλ(prob, sc, 1, x0, g0), Inf) <= 1e-6
        return x0, f0, g0
    end

    x = x0

    # 5.12 in Kelley IMO
    # what's it in Bertsekas?
    local fx, gx
    for i = 1:itermax
        fx, gx = f(x), g(x)
        if print_trace
            println("i = $i    fx: $(fx)    |gx|: $(norm(gx,Inf))")
            println("x1 = $(x[1])")
            println("x2 = $(x[2])")
        end
        if norm(x .- update_xλ(prob, sc, 1, x, gx)) <= 1e-6
            return x, fx, gx
        end
        for j = 1:ls_itermax
            λ = β^j
            xλ = update_xλ(P, sc, λ, x, gx)
            fλ = f(xλ)
            if sufficient_decrease(fλ, fx, λ, α, x, xλ, gx)
                x = xλ
                break
            end
        end
    end
    x, fx, gx
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
using LeastSquares
using PyPlot

# :) it works for the multivariate case!
test = @afuns (x, y) [1, x + y, y]
test[3](1, 5)

x = [0 0; 1 0; 0 1; 1 1; 0.5 0.5]
y = Float64[0, 1, 2, 0, 1]
plane = lsfit(@afuns((x, y), [1, x, y]), x, y)
plane(5, 4)

# so I need to check the arg list better
badlst = push!([], [1, 2, 3])
tmp = @afuns x badlst
tmp[1](5)

tex = :(hcat([1, 2, 3], [5, 6, 7]))

# so I need to check the arg list better, how do I make sure these are simple functions...
tmp = @afuns x hcat([1, 2, 3], [5, 6, 7])
tmp[1](5)

test = @afuns x [1, 5 + 2*sin(x)]
test[2](6)

test = @afuns x [1, 2*x]
test[1](5)

x = [0.0, 1.0, 3.0, 5.0]
y = [1.0, 0.0, 2.0, 4.0]

feq = lsfit(@afuns(x, [1, x]), x, y)
feq(5)
feq.coef
linreg(x, y)
plot(x, y, "o")
plot(x, map(feq, xi))

feq = lsfit(@afuns(x, [1, x, x^2]), x, y)
feq(5)

x = [-π, -π/2.0, 0, π/2, π]
y = [4.0, 0.0, 1.0, -1.0, -4.0]
feq = lsfit(@afuns(x, [sin(x/2), sin(x), sin(2*x)]), x, y)
feq(5)
plot(x, y, "o")
xi = linspace(-2π, 2π, 100)
plot(xi, map(feq, fx))

v = Float64[5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5]
size(v)
x = 1:length(v)
size(x)
quad = lsfit(@afuns(x, [1, x, x^2]), x, v)
quart = lsfit(@afuns(x, [1, x, x^2, x^3, x^4]), x, v)
plot(x, v, "o")
xvals = linspace(1, length(v), 100)
plot(xvals, map(quad, xvals))
plot(xvals, map(quart, xvals))

max_stars_repo_name: gabrielgellner/LeastSquares

max_issues_repo_name: gabrielgellner/LeastSquares

max_forks_repo_name: gabrielgellner/LeastSquares


text:
export sym_H, orthotropic_H, spd_H, spd_Cholesky, spd_Chol_Orth, spd_zero_to_H

# @doc raw"""
#     sym_H(y::PyObject)
#     sym_H(o::Array)

# Creates a symmetric matrix from 6 parameters
# ```math
# H = \begin{bmatrix}
# y_1 & y_2 & y_3\\
# y_2 & y_4 & y_5 \\
# y_3 & y_5 & y_6
# \end{bmatrix}
# ```
# """
# function sym_H(y::PyObject)
#     y = sym_op(y)
#     z = tf.reshape(y, (-1,3,3))
# end

function sym_H(o::Array)
    [o[1] o[2] o[3];
    o[2] o[4] o[5];
    o[3] o[5] o[6]]
end

@doc raw"""
    orthotropic_H(y::PyObject)
    orthotropic_H(o::Array)

Creates a symmetric matrix from 4 parameters
```math
H = \begin{bmatrix}
y_1 & y_2 & 0\\
y_2 & y_3 & 0 \\
0 & 0 & y_4
\end{bmatrix}
```
"""
function orthotropic_H(y::PyObject)

    y = orthotropic_op(y)

    z = tf.reshape(y, (-1,3,3))

    return z
end

function orthotropic_H(o::Array)
    [o[1] o[2] 0.0;
    o[2] o[3] 0.0;
    0.0 0.0 o[4]]
end

@doc raw"""
    spd_H(o::PyObject, H0::Array{Float64,2})
    spd_H(o::Array, H0::Array)

Creates a SPD matrix from 3 scalars

```math
H = H_0 - \frac{H_0nn'H_0}{1+n'H_0n}
```
"""
function spd_H(o::PyObject, H0::Array{Float64,2})
    if size(o,2)!=3
        error("NNFEM: second dimension of `o` must be 3")
    end
    spd_op = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib", "spd_op")
    ret = spd_op(constant(H0), o)

    # ret.set_shape((-1,3,3))
    return ret
end

function spd_H(o::Array, H0::Array)
    # @show size(o'*H0*o)
    o = o[:]
    H0 - H0*(o*o')*H0/(1+o'*H0*o)
end

@doc raw"""
    spd_Cholesky(o::Array)
    spd_Cholesky(o::PyObject)

Creates a SPD matrix from 6 scalars.

```math
A = LL'
```
where
```math
L = \begin{bmatrix}
o_1 & &  \\
o_2 & o_4 & \\
o_3 & o_5 & o_6
\end{bmatrix}
```
"""
function spd_Cholesky(o::Array)
    # @show size(o'*H0*o)
    [o[1]*o[1] o[1]*o[2] o[1]*o[3];
     o[1]*o[2] o[2]*o[2]+o[4]*o[4] o[2]*o[3]+o[4]*o[5];
     o[1]*o[3] o[2]*o[3]+o[4]*o[5] o[3]*o[3]+o[5]*o[5]+o[6]*o[6]]
end

function spd_Cholesky(o::PyObject)
    if size(o,2)!=6
        error("NNFEM: second dimension of `o` must be 6")
    end

    chol_op = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib", "chol_op")
    ret = chol_op(o)

    tf.reshape(ret, (-1,3,3))
end

@doc raw"""
    spd_Chol_Orth(o::Array)
    spd_Chol_Orth(o::PyObject)

Creates a SPD matrix from 4 scalars.

```math
A = LL'
```
where
```math
L = \begin{bmatrix}
o_1 & & \\
o_2 & o_3 & \\
 &  & o_4
\end{bmatrix}
```
"""
function spd_Chol_Orth(o::Array)
    # @show size(o'*H0*o)
    [o[1]*o[1] o[1]*o[2] 0;
     o[1]*o[2] o[2]*o[2]+o[3]*o[3] 0.0;
     0.0       0.0                 o[4]*o[4]]
end

function spd_Chol_Orth(o::PyObject)
    if size(o,2)!=4
        error("NNFEM: second dimension of `o` must be 4")
    end

    chol_orth_op = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib", "chol_orth_op")
    ret = chol_orth_op(o)

    tf.reshape(ret, (-1,3,3))
end

@doc raw"""
    spd_zero_to_H(o::Array)
    spd_zero_to_H(o::Array, H0inv::Array{Float64,2})

Creates a SPD matrix from 4 scalars.

```math
A = (H_0^{-1} +LL')^{-1}
```
where
```math
L =  \begin{bmatrix}
o_1 & & \\
o_2 & o_3 & \\
 &  & o_4
\end{bmatrix}
```

"""
function spd_zero_to_H(o::PyObject, H0inv::Array{Float64,2})
    if size(o,2)!=4
        error("NNFEM: second dimension of `o` must be 4")
    end

    chol_orth_op = load_op_and_grad("$(@__DIR__)/../../deps/CustomOp/build/libDataLib", "chol_orth_op")
    ret = chol_orth_op(o)

    out = tf.reshape(ret, (-1,3,3))
    inv(H0inv+out)
end

function spd_zero_to_H(o::Array, H0inv::Array{Float64,2})

    # ret = chol_orth_op(o)
    # @show ret
    # out = tf.reshape(ret, (-1,3,3))
    # inv(H0inv+out)

    spd_o = [o[1]*o[1] o[1]*o[2] 0;
         o[1]*o[2] o[2]*o[2]+o[3]*o[3] 0.0;
         0.0       0.0                 o[4]*o[4]]

    inv(H0inv+spd_o)
end


text:
using QuantumGate
using Test
using SparseArrays

@testset "Pauli matrices" begin
    @test σˣ^2 == σ⁰
    @test σʸ^2 == σ⁰
    @test σᶻ^2 == σ⁰
    @test σˣ*σʸ - σʸ*σˣ == 2im*σᶻ
    @test σʸ*σᶻ - σᶻ*σʸ == 2im*σˣ
    @test σᶻ*σˣ - σˣ*σᶻ == 2im*σʸ
    @test σ⁰ == P⁰ + P¹
    @test σᶻ == P⁰ - P¹
end

@testset "Kronecker product" begin
    A = [1 0; 0 1]
    B = [1 0; 0 -1]
    C = [1 0 0 0; 0 -1 0 0; 0 0 1 0; 0 0 0 -1]
    @test A⊗B == C
    @test ⊗(A,A,B,B) == A⊗A⊗B⊗B

    A, B, C, D = rand(2, 2), rand(2, 2), rand(2, 2), rand(2, 2)
    @test A⊗(B+C) ≈ A⊗B + A⊗C
    @test (B+C)⊗A ≈ B⊗A + C⊗A
    @test (A⊗B)⊗C ≈ A⊗(B⊗C)
    @test (A⊗B)*(C⊗D) ≈ (A*C)⊗(B*D)
end

@testset "Circuit identities" begin
    @test Had * σˣ * Had ≈ σᶻ
    @test Had * σʸ * Had ≈ -σʸ
    @test Had * σᶻ * Had ≈ σˣ
end

@testset "Zgate" begin
    @test σᶻ ⊗ σ⁰ == Zgate(1, 2)
    @test σ⁰ ⊗ σᶻ == Zgate(2, 2)
    @test σᶻ ⊗ σ⁰ ⊗ σ⁰ == Zgate(1, 3)
    @test σ⁰ ⊗ σᶻ ⊗ σ⁰ ⊗ σ⁰ == Zgate(2, 4)
end

@testset "CNOTgate" begin
    @test CNOTgate(1, 2, 2) == sparse([
                                    1.0+0.0im 0.0 0.0 0.0;
                                    0.0 1.0 0.0 0.0;
                                    0.0 0.0 0.0 1.0;
                                    0.0 0.0 1.0 0.0
                                ])
    @test CNOTgate(1, 2, 2) == P⁰ ⊗ σ⁰ + P¹ ⊗ σˣ
    @test CNOTgate(1, 2, 3) == P⁰ ⊗ σ⁰ ⊗ σ⁰ + P¹ ⊗ σˣ ⊗ σ⁰
    @test CNOTgate(1, 3, 3) == P⁰ ⊗ σ⁰ ⊗ σ⁰ + P¹ ⊗ σ⁰ ⊗ σˣ
    @test CNOTgate(3, 1, 3) == σ⁰ ⊗ σ⁰ ⊗ P⁰ + σˣ ⊗ σ⁰ ⊗ P¹
end

@testset "SWAP gate" begin
    @test SWAPgate(1, 2, 2) == sparse([
                                    1.0+0.0im 0.0 0.0 0.0;
                                    0.0 0.0 1.0 0.0;
                                    0.0 1.0 0.0 0.0;
                                    0.0 0.0 0.0 1.0
                                ])

    @test SWAPgate(1, 2, 2) == SWAPgate(2, 1, 2)
    @test SWAPgate(1, 3, 3) == SWAPgate(3, 1, 3)
end


text:
using Flux: onecold
using JSON
using CSV
using DataFrames

# default number of bins for continuous features, if unspecified
const default_nbins = 10

# read in features...
atom_data_path = joinpath(dirname(pathof(ChemistryFeaturization)), "..", "data", "pymatgen_atom_data.csv")
const atom_data_df = DataFrame!(CSV.File(atom_data_path))
feature_info_path = joinpath(dirname(pathof(ChemistryFeaturization)), "..", "data", "feature_info.json")
const feature_info = JSON.parsefile(feature_info_path)

const categorical_feature_names = Symbol.(feature_info["categorical"])
const categorical_feature_vals = Dict(fea=>sort(collect(Set(skipmissing(atom_data_df[:, fea])))) for fea in categorical_feature_names)
# but I want blocks to be in my order
categorical_feature_vals[:Block] = ["s", "p", "d", "f"]
const continuous_feature_names = Symbol.(feature_info["continuous"])
const not_features = Symbol.(feature_info["not_features"]) # atomic name, symbol
const avail_feature_names = cat(categorical_feature_names, continuous_feature_names; dims=1)

# compile min and max values of each feature...
const fea_minmax = Dict{Symbol, Tuple{Real, Real}}()
for feature in avail_feature_names
    if !(feature in categorical_feature_names)
        minval = minimum(skipmissing(atom_data_df[:, feature]))
        maxval = maximum(skipmissing(atom_data_df[:, feature]))
        fea_minmax[feature] = (minval, maxval)
    end
end

# Type to store featurization metadata. An array of these specifies a featurization scheme for atoms.
# TODO for Sean, eventually: add PairFeat and BondFeat types (maybe in another file for tidiness)
# TODO, maybe: we could define AbstractFeat that all of these would inherit from, but TBD if that would be useful or not
struct AtomFeat
    name::Symbol # name of feature
    categorical::Bool # whether it's categorical (vs. numerical)
    num_bins::Integer # length of associated subvector
    logspaced::Bool # whether it's logspaced (most relevant for numerical)
    vals::Vector # list of values (length equal to num_bins for categorical, num_bins+1 (specifying bin edges) for numerical)
    # basic standard constructor will check some things...
    function AtomFeat(name::Symbol, categorical::Bool, num_bins::Integer, logspaced::Bool, vals::Vector)
        T = eltype(vals)
        if T==Float64
            T = Float32
        elseif T==Int64
            T = Int32
        end
        if categorical
            if num_bins != length(vals)
                DimensionMismatch("Categorical features should have a number of values equal to the number of bins.")
            end
            if logspaced
                warn("A logarithmically spaced categorical feature doesn't make much sense...are you sure that's what you meant?")
            end
            val_list = vals
        else # numerical feature
            if num_bins + 1 != length(vals)
                DimensionMismatch("Numerical features should have values specifying bin edges. This vector is the wrong length!")
            end
            if !(T<:Real)
                error("Numerical features must have (real) numerical values...") # should figure out how to throw a TypeError propertly
            end
            val_list = sort(vals)
        end
        new(name, categorical, num_bins, logspaced, val_list)
    end
end

# constructor for features where vals gets calculated rather than passed in, works for categorical too if their values are numbers
function AtomFeat(name::Symbol, categorical::Bool, num_bins::Integer, min_val::Real, max_val::Real, logspaced::Bool=false; T=Float32)
    categorical ? len = num_bins : len = num_bins + 1
    if logspaced
        vals = 10 .^ range(log10(min_val), log10(max_val), length=len)
    else
        vals = range(min_val, max_val, length=len)
    end
    AtomFeat(name, categorical, num_bins, logspaced, T.([v for v in vals]))
end

# constructor that will assume categorical features
AtomFeat(name::Symbol, vals::Vector) = AtomFeat(name, true, length(vals), false, vals)

# pretty printing, short form
Base.show(io::IO, f::AtomFeat) = print(io, "$(f.name): AtomFeat with $(f.num_bins) bins")

# pretty printing, long form
Base.show(io::IO, ::MIME"text/plain", f::AtomFeat) = print(io, "AtomFeat\n   name: $(f.name)\n   categorical: $(f.categorical)\n   length: $(f.num_bins)\n   logspaced: $(f.logspaced)\n   bins: $(f.vals)")

"""
Create onehot style vector.

# Arguments
- `feature::AtomFeat`: feature being encoded
- `val`: value of feature to encode

See also: [onecold_bins](@ref), [onehot](@ref Flux.onehot)
"""
function onehot_bins(f::AtomFeat, val)
    if f.categorical
        bin_index = findfirst(isequal(val), f.vals)
    else
        bin_index = searchsorted(f.vals, val).stop
        if bin_index == length(f.vals) # got the max value
            bin_index = bin_index - 1
        elseif isapprox(val, f.vals[1])
            # sometimes it gives 0 otherwise
            bin_index = 1
        end
    end
    onehot_vec = [0.0 for i in 1:f.num_bins]
    onehot_vec[bin_index] = 1.0
    return onehot_vec
end

# TODO, maybe: add signatures for this and onecold version where you give the values rather than the feature object?

"""
    onecold_bins(feature, vec)

Inverse function to onehot_bins, decodes a vector corresponding to a given feature. For a categorical feature, returns a value, for a continuous one, a range.

# Arguments
- `feature::AtomFeat`: feature being encoded
- `vec::Vector`: vector (such as produced by [onehot_bins](@ref))

See also: [onehot_bins](@ref), [onecold](@ref Flux.onecold)

"""
function onecold_bins(f::AtomFeat, vec::Vector)
    if f.categorical
        decoded = onecold(vec, f.vals)
    else
        decoded = decoded = (onecold(vec, f.vals[1:end-1]), onecold(vec, f.vals[2:end]))
    end
    return decoded
end

"Little helper function to check that the logspace vector/boolean is appropriate and convert it to a vector as needed."
function get_logspaced_vec(vec, num_features::Integer)
    if vec==false # default behavior
        logspaced_vec = [false for i in 1:num_features]
    elseif vec==true
        logspaced_vec = [true for i in 1:num_features]
    elseif length(vec) == num_features # specified properly
        logspaced_vec = vec
    elseif length(vec) < num_features
        println("logspaced vector too short. Padding end with falses.")
        logspaced_vec = vcat(vec, [false for i in 1:num_features-size(vec,1)])
    elseif size(vec, 1) > num_features
        println("logspaced vector too long. Cutting off at appropriate length.")
        logspaced_vec = vec[1:num_features]
    end
    return logspaced_vec
end

"""
Function to build a featurization given vectors of metadata.

Note that nbins will be ignored for categorical features.
"""
# TODO: rename this or make_feature_vectors because they sound like they do the same thing...
function build_atom_feats(feature_names::Vector{Symbol}; nbins::Vector{<:Integer}=default_nbins*ones(Int64, size(feature_names,1)), logspaced=false)
    num_features = length(feature_names)

    # figure out spacing for each feature
    logspaced_vec = get_logspaced_vec(logspaced, num_features)

    # build AtomFeat objects
    # TODO: could make this a preallocation since we know the length, probably not critical though
    feature_specs = AtomFeat[]
    #for i in 1:num_features
    for t in zip(feature_names, nbins, logspaced_vec)
        feature_name = t[1]
        nbins = t[2]
        logspaced = t[3]
        if feature_name in categorical_feature_names
            vals = categorical_feature_vals[Symbol(feature_name)]
            feature = AtomFeat(feature_name, vals)
        elseif feature_name in continuous_feature_names
            feature = AtomFeat(feature_name, false, nbins, fea_minmax[feature_name]..., logspaced)
        else
            @warn "$feature_name not found in atom_data_df..."
        end
        push!(feature_specs, feature)
    end
    return feature_specs
end

"""
    make_feature_vectors(featurization)
    make_feature_vectors(feature_names)
    make_feature_vectors(feature_names, nbins)
    make_feature_vectors(feature_names, nbins, logspaced)

Make custom feature vectors, using specified features and numbers of bins. Can be called with an array of AtomFeat objects or by specifying names and other metadata and the objects will be built.

Note that in the latter case, bin numbers will be ignored for categorical features (block, group, and row), but features and nbins vectors should still be the same length. Optionally, feed in vector of booleans with trues at the index of any (continous valued) feature whose bins should be log spaced.

# Arguments
- `featurization::Vector{AtomFeat}`: Featurization scheme in the form of a list of AtomFeat objects. (See below for alternate call signature that generates this list from more primitive inputs)

# Returns
- a dictionary from element symbol => one-hot style feature vector, concatenated in order of feature list.
- A Vector of AtomFeat objects
"""
function make_feature_vectors(featurization::Vector{AtomFeat})
    feature_names = [f.name for f in featurization]
    usable_atom_data =  dropmissing(atom_data_df[:, cat(Symbol.(feature_names), not_features, dims=1)])
    num_dropped = size(atom_data_df)[1] - size(usable_atom_data)[1]
    if num_dropped != 0
        @info "$(num_dropped) elements were dropped so that all features are defined."
    end

    # dict from element symbol to feature vec of that element
    # (if we do any structure-specific features later, e.g. coordination or something,
    # this will have to iterate over every atom in the structure instead...)
    # (but possibly would just want to append those to the end anyway...)
    sym_featurevec = Dict{String, Vector{Float32}}()

    for i in 1:size(usable_atom_data,1)
        el = usable_atom_data.Symbol[i]
        featurevec = []
        # make onehot vector for each feature
        for f in featurization
            val = usable_atom_data[i, Symbol(f.name)]
            subvec = onehot_bins(f, val)
            append!(featurevec, subvec)
        end
        sym_featurevec[el] = featurevec
    end
    return sym_featurevec, featurization
end

"""
alternate call signature...

# Arguments
- `feature_names::Vector{String}`: list of features to be encoded
- `nbins::Vector{Integer}`: number of bins for each feature (in same order)
- `logspaced=false`: (single Bool or vector of them) whether or not to logarithmically space each feature
"""
make_feature_vectors(feature_names::Vector{Symbol}; nbins::Vector{<:Integer}=default_nbins*ones(Int64, size(feature_names,1)), logspaced=false) = make_feature_vectors(build_atom_feats(feature_names; nbins=nbins, logspaced=logspaced))

"""
    chunk_vec(vec, nbins)

Divide up an already-constructed feature vector into "chunks" (presumably one for each feature) of lengths specified by the vector nbins.

Sum of nbins should be equal to the length of vec.

# Examples
```jldoctest
julia> chunk_vec([1,0,0,1,0], [3,2])
2-element Array{Array{Bool,1},1}:
 [1, 0, 0]
 [1, 0]
 ```
"""
function chunk_vec(vec::Vector, nbins::Vector{<:Integer})
    chunks = fill(Bool[], size(nbins, 1))
    @assert length(vec)==sum(nbins) "Total number of bins doesn't match length of feature vector."
    for i in 1:size(nbins,1)
        if i==1
            start_ind = 1
        else
            start_ind = sum(nbins[1:i-1])+1
        end
        chunks[i] = vec[start_ind:start_ind+nbins[i]-1]
    end
    return chunks
end

"""
Check that each subvector in a featurization is a valid one-hot encoding (one true and otherwise all falses).
"""
function vec_valid(vec::Vector, nbins::Vector{<:Integer})
    result = true
    # are there the right number of total bins?
    try
        chunks = chunk_vec(vec, nbins)
    catch e
        result = false
    end
    if result # if we're still good up to now...
        chunks = chunk_vec(vec, nbins) # apparently chunks is local to the previous try block?
        # does each subvector have exactly one true value?
        for i in 1:length(nbins)
            subvec = chunks[i]
            if !(sum(subvec)==1)
                println("Subvector ", i, " is invalid.")
                result = false
            end
        end
    end
    return result
end

"""
Function to invert the binning process. Useful to check that it's working properly, or just to inspect properties once they've been encoded.

Need to feed in a feature vector as well as the lists of features and bin numbers that were used to encode it.
"""
function decode_feature_vector(vec::Vector, featurization::Vector{AtomFeat})
    nbins = [f.num_bins for f in featurization]
    # First, check that the featurization is valid
    if !(vec_valid(vec, nbins))
        error("Vector is invalid!")
    else
        chunks = chunk_vec(vec, nbins)

        # make dict from feature names to corresponding chunks in vector
        fea_chunks = Dict(zip(featurization, chunks))

        return Dict(f.name=>onecold_bins(f, fea_chunks[f]) for f in featurization)
    end
end

# alternate call signature
decode_feature_vector(vec::Vector{<:Real}, feature_names::Vector{String}, nbins::Vector{<:Integer}, logspaced=false) = decode_feature_vector(vec, build_atom_feats(feature_names, nbins, logspaced))


text:
"""
	SimulatorModule
Simulator for lattice.

# modifying _move this when finished_
* simulator_quantum_2d_fractal_initialiteration.jl
  for adding additional iterations for the trotter layer
  * getSingularValues

# Type list
* Simulator,
* ClassicalSimulator,
* Classical2dSimulator,
* Classical3dSimulator,
* Classical2dSquareSimulator,
* QuantumSimulator,
* Quantum2dSimulator,
* Quantum2dSquareSimulator,
* Classical3dSquareSimulator,
* FractalSimulator,
* Classical2dFractalSimulator,
* Quantum2dFractalSimulator,
* Quantum2dFractalInititerSimulator <- working

# Method list
* buildSimulator, #
* getDimM,
* getExpectationValue,
* getFreeEnergy,
* getWholeiteration,
* getNormalizationfactors,
* getData4Energy,
* getCount,
* getNumberOfSites,
* countUp!,
* countDown!,
* isDone,
* setNormalizationfactor!,
* normalizeTensor,
* writeVector,
* isZeroTemperature,
* isClassical,
* getModelname,
* getStates,
* isSymmetricFactorization,
* isPottsModel,
* getHamiltonian,
* getTemperature,
* getExternalfield,
* getEnvParameters,
* getCoareserate,
* getTensorT,
* getTensorW,
* getMeasureOperator,
* setTemperature!,
* setExternalfield!,
* setTensorT!,
* setEnvParameters!,
* initializeCount!,
* getFirstTerm,
* renormalize!,
* getTrotterCount,
* getSpaceCount,
* getTrotterparameter,
* getTrotteriteration,
* getTrotterlayers,
* getBeta, # changed to use trottercount
* setTrotterparameter!,
* setTrotterIteration!,
* renormalize,
* getNewTensorT_2dQ,
* getTensorU,
* getTensorV,
* getTenmatMMd,
* renormalizeX!,
* renormalizeZ!,
* renormalizeY!,
* truncMatrixU,
* matU2tenU,
* simulatorTemperature,
* simulatorQuantum,
* getHausdorffDim,
* getFractalDim,
* getLegextension,
* getTensorP,
* getTensorQ,
* getCoefficient,
* setTensorP!,
* setTensorQ!,
* normalizeTensor!,
* setNorm!,
* constructHalf,
* updateLocalTensors!,
* debug_updateLocalTensors!,
* getNewLegTensor,
* getNewCaretTensor,
* # working_simulator_q2f_re.jl
* renormalizeSpace!,
* renormalizeTrotter!,
* getTensorUy,
* # new version of
* # simulator_quantum_2d_square.jl
* updateCoefficient!,
* setCoefficient!,
* getCoefficient,
* # for debug:
* calculateCoreTensors
* getInititeration(simulator::Quantum2dFractalinititerSimulator)

# Test list
* testClassical2dSquareSimulator,
* testQuantum2dSquareSimulator,
* testSimulateTemp,
* testClassical3dSquareSimulator,
* testClassical2dFractalSimulator,
* testQuantum2dFractalSimulator,
"""
module SimulatorModule
	using HoTRG_lemon.LatticeModule
	using HoTRG_lemon.SpinModule
	using TensorOperations
	using TensorMatrices_lemon

	export
		# simulator.jl
		Simulator,
		buildSimulator, #
		getDimM,
		getExpectationValue,
		getFreeEnergy,
		getWholeiteration,
		getNormalizationfactors,
		getData4Energy,
		getCount,
		getNumberOfSites,
		countUp!,
		countDown!,
		isDone,
		setNormalizationfactor!,
		normalizeTensor,
		writeVector,
			# <inherit from lattice>
		isZeroTemperature,
		isClassical,
		getModelname,
		getStates,
		isSymmetricFactorization,
		isPottsModel,
		getHamiltonian,
		getTemperature,
		getExternalfield,
		getEnvParameters,
		getCoareserate,
		getTensorT,
		getTensorW,
		getMeasureOperator,
		setTemperature!,
		setExternalfield!,
		setTensorT!,
		setEnvParameters!,
		# simulator_classical.jl
		ClassicalSimulator,
		Classical2dSimulator,
		Classical3dSimulator,
		initializeCount!,
		getFirstTerm,
		# simulator_classical_2d_square.jl
		Classical2dSquareSimulator,
		testClassical2dSquareSimulator,
		# simulator_classical_2d_square_renormalize.jl
		renormalize!,
		# simulator_quantum.jl
		QuantumSimulator,
		Quantum2dSimulator,
		getTrotterCount,
		getSpaceCount,
			# <inherit from QuantumLattice>
		getTrotterparameter,
		getTrotteriteration,
		getTrotterlayers,
		getBeta, # changed to use trottercount
		setTrotterparameter!,
		setTrotterIteration!,
		# simulator_quantum_2d_square.jl
		Quantum2dSquareSimulator,
		testQuantum2dSquareSimulator,
		# simulator_quantum_2d_square_renormalize.jl
		renormalize,
		##### to debug
		getNewTensorT_2dQ,
		getTensorU,
		getTensorV,
		getTenmatMMd,
		renormalizeX!,
		renormalizeZ!,
		renormalizeY!,
		truncMatrixU,
		matU2tenU,
		# simulator_temperature.jl
		simulatorTemperature,
		testSimulateTemp,
		# simulator_quantum_data.jl
		simulatorQuantum,
		# simulator_classical_3d_square.jl
		Classical3dSquareSimulator,
		testClassical3dSquareSimulator,
		# simulator_fractal.jl
		FractalSimulator,
			# inherit from FractalLattice
		getHausdorffDim,
		getFractalDim,
		getLegextension,
		getTensorP,
		getTensorQ,
		getCoefficient,
		setTensorP!,
		setTensorQ!,
		normalizeTensor!,
		setNorm!,
		# simulator_classical_2d_fractal.jl
		Classical2dFractalSimulator,
		testClassical2dFractalSimulator,
		# for debugging
		constructHalf,
		updateLocalTensors!,
		debug_updateLocalTensors!,
		getNewLegTensor,
		getNewCaretTensor,
		# working_simulator_q2f.jl
		Quantum2dFractalSimulator,
		testQuantum2dFractalSimulator,
		# working_simulator_q2f_re.jl
		renormalizeSpace!,
		renormalizeTrotter!,
		getTensorUy,
		# new version of
		# simulator_quantum_2d_square.jl
		updateCoefficient!,
		setCoefficient!,
		getCoefficient,
		# for debug:
		calculateCoreTensors,
		# simulator_quantum_2d_fractal_inititer.jl
		Quantum2dFractalInititerSimulator,
		getInititeration,
		printLog,
		#
		getInititeration

	include("simulator.jl")
	include("simulator_classical.jl")
	include("simulator_classical_2d_square.jl")
	include("simulator_classical_2d_square_renormalize.jl")
	include("simulator_quantum.jl")
	include("simulator_quantum_2d_square.jl")
	include("simulator_quantum_2d_square_renormalize.jl")
	include("working_simulator_q2f.jl")
	include("working_simulator_q2f_re.jl")
	#= for the stable old version
	include("simulator_quantum_2d_fractal.jl")
	include("simulator_quantum_2d_fractal_renormalize.jl")
	=#
	include("simulator_classical_temperature.jl")
	include("simulator_quantum_data.jl")
	include("simulator_classical_3d_square.jl")
	include("simulator_classical_3d_square_renormalize.jl")
	include("simulator_classical_2d_fractal_2.jl")
	include("simulator_classical_2d_fractal_renormalize_2.jl")
	 include("simulator_quantum_2d_fractal_inititer.jl")
	 include("renormalize_simulator_quantum_2d_fractal_inititer.jl")
	include("simulator_fractal.jl")
	include("f_normalizeAndSetNorm!.jl")
	include("f_setNorm!.jl")
	include("f_printlog.jl")
	include("f_getInititeration.jl")

end


text:
module ShadowCrushBlock

using ..Ahorn, Maple

@mapdef Entity "brokemiahelper/shadowCrushBlock" ShadowKevin(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, axes::String="both", chillout::Bool=false)

const placements = Ahorn.PlacementDict(
    "Shadow Kevin (Both)" => Ahorn.EntityPlacement(
        ShadowKevin,
        "rectangle"
    ),
    "Shadow Kevin (Vertical)" => Ahorn.EntityPlacement(
        ShadowKevin,
        "rectangle",
        Dict{String, Any}(
            "axes" => "vertical"
        )
    ),
    "Shadow Kevin (Horizontal)" => Ahorn.EntityPlacement(
        ShadowKevin,
        "rectangle",
        Dict{String, Any}(
            "axes" => "horizontal"
        )
    ),
)

const frameImage = Dict{String, String}(
    "none" => "objects/crushblock/block00",
    "horizontal" => "objects/crushblock/block01",
    "vertical" => "objects/crushblock/block02",
    "both" => "objects/crushblock/block03"
)

const smallFace = "objects/crushblock/idle_face"
const giantFace = "objects/crushblock/giant_block00"

const kevinColor = (98, 34, 43) ./ 255

Ahorn.editingOptions(entity::ShadowKevin) = Dict{String, Any}(
    "axes" => Maple.kevin_axes
)

Ahorn.minimumSize(entity::ShadowKevin) = 24, 24
Ahorn.resizable(entity::ShadowKevin) = true, true

Ahorn.selection(entity::ShadowKevin) = Ahorn.getEntityRectangle(entity)

# Todo - Use randomness to decide on Kevin border
function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::ShadowKevin, room::Maple.Room)
    axes = lowercase(get(entity.data, "axes", "both"))
    chillout = get(entity.data, "chillout", false)

    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    giant = height >= 48 && width >= 48 && chillout
    face = giant ? giantFace : smallFace
    frame = frameImage[lowercase(axes)]
    faceSprite = Ahorn.getSprite(face, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    Ahorn.drawRectangle(ctx, 2, 2, width - 4, height - 4, kevinColor)
    Ahorn.drawImage(ctx, faceSprite, div(width - faceSprite.width, 2), div(height - faceSprite.height, 2))

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, 0, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, height - 8, 8, 24, 8, 8)
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, 0, (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, width - 8, (i - 1) * 8, 24, 8, 8, 8)
    end

    Ahorn.drawImage(ctx, frame, 0, 0, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, 0, 24, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, 0, height - 8, 0, 24, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, height - 8, 24, 24, 8, 8)
end

end

max_stars_repo_name: Brokemia/BrokemiaHelper

max_issues_repo_name: Brokemia/BrokemiaHelper

max_forks_repo_name: Brokemia/BrokemiaHelper