text:
# some basic definitions for generic matches

execute(k::Config, m, s, i) = error("$m did not expect to be called with state $s")
success(k::Config, m, s, t, i, r) = error("$m did not expect to receive state $s, result $r")
failure(k::Config, m, s) = error("$m did not expect to receive state $s, failure")

execute(k::Config, m::Matcher, s::Dirty, i) = FAILURE

# many matchers delegate to a child, making only slight modifications.
# we can describe the default behaviour just once, here.
# child matchers then need to implement (1) state creation (typically on
# response) and (2) anything unusual (ie what the matcher actually does)

# assume this has a matcher field
abstract type Delegate<:Matcher end

# assume this has a state field
abstract type DelegateState<:State end

execute(k::Config, m::Delegate, s::Clean, i) = Execute(m, s, m.matcher, CLEAN, i)

execute(k::Config, m::Delegate, s::DelegateState, i) = Execute(m, s, m.matcher, s.state, i)

failure(k::Config, m::Delegate, s) = FAILURE

# various weird things for completeness

@auto_hash_equals mutable struct Epsilon<:Matcher
    name::Symbol
    Epsilon() = new(:Epsilon)
end

execute(k::Config, m::Epsilon, s::Clean, i) = Success(DIRTY, i, EMPTY)

@auto_hash_equals mutable struct Insert<:Matcher
    name::Symbol
    text
    Insert(text) = new(:Insert, text)
end

execute(k::Config, m::Insert, s::Clean, i) = Success(DIRTY, i, Any[m.text])

@auto_hash_equals mutable struct Dot<:Matcher
    name::Symbol
    Dot() = new(:Dot)
end

function execute(k::Config, m::Dot, s::Clean, i)
    if iterate(k.source, i) == nothing
        FAILURE
    else
        c, i = iterate(k.source, i)
        Success(DIRTY, i, Any[c])
    end
end

@auto_hash_equals mutable struct Fail<:Matcher
    name::Symbol
    Fail() = new(:Fail)
end

execute(k::Config, m::Fail, s::Clean, i) = FAILURE

# evaluate the sub-matcher, but replace the result with EMPTY

@auto_hash_equals mutable struct Drop<:Delegate
    name::Symbol
    matcher::Matcher
    Drop(matcher) = new(:Drop, matcher)
end

@auto_hash_equals struct DropState<:DelegateState
    state::State
end

success(k::Config, m::Drop, s, t, i, r::Value) = Success(DropState(t), i, EMPTY)

# exact match

@auto_hash_equals mutable struct Equal<:Matcher
    name::Symbol
    string
    Equal(string) = new(:Equal, string)
end

always_print(::Equal) = true
function print_field(m::Equal, ::Type{Val{:string}})
    if isa(m.string, AbstractString)
        "\"$(m.string)\""
    else
        :string
    end
end

function execute(k::Config, m::Equal, s::Clean, i)
    for x in m.string
        if iterate(k.source, i) == nothing
            return FAILURE
        end
        y, i = iterate(k.source, i)
        if x != y
            return FAILURE
        end
    end
    Success(DIRTY, i, Any[m.string])
end

# repetition

# in the first two cases (Depth and Breadth) we perform a tree search
# of all possible states (limited by the maximum number of matches),
# yielding when we have a result within the lo/hi range.

abstract type Repeat_<:Matcher end   # _ to avoid conflict with function in 0.3

ALL = typemax(Int)

abstract type RepeatState<:State end

function Repeat(m::Matcher, lo, hi; flatten=true, greedy=true, backtrack=true)
    if greedy
        backtrack ? Depth(m, lo, hi; flatten=flatten) : Depth!(m, lo, hi; flatten=flatten)
    else
        backtrack ? Breadth(m, lo, hi; flatten=flatten) : Breadth!(m, lo, hi; flatten=flatten)
    end
end
Repeat(m::Matcher, lo; flatten=true, greedy=true, backtrack=true) = Repeat(m, lo, lo; flatten=flatten, greedy=greedy, backtrack=backtrack)
Repeat(m::Matcher; flatten=true, greedy=true, backtrack=true) = Repeat(m, 0, ALL; flatten=flatten, greedy=greedy, backtrack=backtrack)

print_field(m::Repeat_, ::Type{Val{:lo}}) = "lo=$(m.lo)"
print_field(m::Repeat_, ::Type{Val{:hi}}) = "hi=$(m.hi)"
print_field(m::Repeat_, ::Type{Val{:flatten}}) = "flatten=$(m.flatten)"

function repeat_success(m::Repeat_, results::Vector{Value})
    if m.flatten
        flatten(results)
    else
        Any[results...]
    end
end

# depth-first (greedy) state and logic

@auto_hash_equals mutable struct Depth<:Repeat_
    name::Symbol
    matcher::Matcher
    lo::Integer
    hi::Integer
    flatten::Bool
    Depth(m, lo, hi; flatten=true) = new(:Depth, m, lo, hi, flatten)
end

# greedy matching is effectively depth first traversal of a tree where:
# * performing an additional match is moving down to a new level
# * performaing an alternate match (backtrack+match) moves across
# the traversal requires a stack.  the DepthState instances below all
# store that stack - actually three of them.  the results stack is
# unusual / neat in that it is also what we need to return.

# unfortunately, things are a little more complex, because it's not just
# DFS, but also post-order.  which means there's some extra messing around
# so that the node ordering is correct.

abstract type DepthState<:RepeatState end

# an arbitrary iter to pass to a state where it's not needed (typically from
# failure)
arbitrary(s::DepthState) = s.iters[1]

@auto_hash_equals struct DepthSlurp{I}<:DepthState
    # there's a mismatch in lengths here because the empty results is
    # associated with an iter and state
    results::Vector{Value} # accumulated.  starts []
    iters::Vector{I}       # at the end of the result.  starts [i].
    states::Vector{State}  # at the end of the result.  starts [DIRTY],
                           # since [] at i is returned last.
end

@auto_hash_equals struct DepthYield{I}<:DepthState
    results::Vector{Value}
    iters::Vector{I}
    states::Vector{State}
end

@auto_hash_equals struct DepthBacktrack{I}<:DepthState
    results::Vector{Value}
    iters::Vector{I}
    states::Vector{State}
end

# when first called, create base state and make internal transition

execute(k::Config{S,I}, m::Depth, s::Clean, i::I) where {S,I} = execute(k, m, DepthSlurp{I}(Vector{Value}(), I[i], State[DIRTY]), i)

# repeat matching until at bottom of this branch (or maximum depth)

max_depth(m::Depth, results) = m.hi == length(results)

function execute(k::Config, m::Depth, s::DepthSlurp, i)
    if max_depth(m, s.results)
        execute(k, m, DepthYield(s.results, s.iters, s.states), i)
    else
        Execute(m, s, m.matcher, CLEAN, i)
    end
end

function success(k::Config, m::Depth, s::DepthSlurp, t, i, r::Value)
    results = Value[s.results..., r]
    iters = vcat(s.iters, i)
    states = vcat(s.states, t)
    if max_depth(m, results)
        execute(k, m, DepthYield(results, iters, states), i)
    else
        Execute(m, DepthSlurp(results, iters, states), m.matcher, CLEAN, i)
    end
end

function failure(k::Config, m::Depth, s::DepthSlurp)
    # the iter passed on is arbitrary and unused
    execute(k, m, DepthYield(s.results, s.iters, s.states), arbitrary(s))
end

# yield a result and set state to backtrack

function execute(k::Config, m::Depth, s::DepthYield, unused)
    n = length(s.results)
    if n >= m.lo
        Success(DepthBacktrack(s.results, s.iters, s.states), s.iters[end], repeat_success(m, s.results))
    else
        # we need to continue searching in case there's some other weird
        # case that gets us back into valid matches
        execute(k, m, DepthBacktrack(s.results, s.iters, s.states), arbitrary(s))
    end
end

# backtrack once and then move down again if possible.  we cannot repeat a
# path because we always advance child state.

function execute(k::Config, m::Depth, s::DepthBacktrack, unused)
    if length(s.iters) == 1  # we've exhausted the search space
        FAILURE
    else
        # we need the iter from *before* the result
        Execute(m, DepthBacktrack(s.results[1:end-1], s.iters[1:end-1], s.states[1:end-1]), m.matcher, s.states[end], s.iters[end-1])
    end
end

# backtrack succeeded so move down
function success(k::Config, m::Depth, s::DepthBacktrack, t, i, r::Value)
    execute(k, m, DepthSlurp(Value[s.results..., r], vcat(s.iters, i), vcat(s.states, t)), i)
end

# we couldn't move down, so yield this point
function failure(k::Config, m::Depth, s::DepthBacktrack)
    execute(k, m, DepthYield(s.results, s.iters, s.states), arbitrary(s))
end

# breadth-first specific state and logic

@auto_hash_equals mutable struct Breadth<:Repeat_
    name::Symbol
    matcher::Matcher
    lo::Integer
    hi::Integer
    flatten::Bool
    Breadth(m, lo, hi; flatten=true) = new(:Breadth, m, lo, hi, flatten)
end

# minimal matching is effectively breadth first traversal of a tree where:
# * performing an additional match is moving down to a new level
# * performaing an alternate match (backtrack+match) moves across
# the traversal requires a queue.  unfortunately, unlike with greedy,
# that means we need to store the entire result for each node.

# on the other hand, because the results are pre-order, the logic is simpler
# than for the greedy match (wikipedia calls this "level order" so my
# terminology may be wrong).

@auto_hash_equals struct Entry{I}
    iter::I
    state::State
    results::Vector{Value}
end

abstract type BreadthState<:RepeatState end

arbitrary(s::BreadthState) = s.start

@auto_hash_equals struct BreadthGrow{I}<:BreadthState
    start::I  # initial iter
    queue::Vector{Entry{I}}  # this has to be immutable for caching
end

@auto_hash_equals struct BreadthYield{I}<:BreadthState
    start::I  # initial iter
    queue::Vector{Entry{I}}  # this has to be immutable for caching
end

# when first called, create base state and make internal transition

execute(k::Config{S,I}, m::Breadth, s::Clean, i::I) where {S,I} = execute(k, m, BreadthYield{I}(i, Entry{I}[Entry{I}(i, CLEAN, Any[])]), i)

# yield the top state

function execute(k::Config, m::Breadth, s::BreadthYield, unused)
    q = s.queue[1]
    n = length(q.results)
    if n >= m.lo
        Success(BreadthGrow(s.start, s.queue), q.iter, repeat_success(m, q.results))
    else
        execute(k, m, BreadthGrow(s.start, s.queue), unused)
    end
end

# add the next row

function execute(k::Config, m::Breadth, s::BreadthGrow, unused)
    if length(s.queue[1].results) > m.hi
        FAILURE
    else
        Execute(m, s, m.matcher, CLEAN, s.queue[1].iter)
    end
end

success(k::Config, m::Breadth, s::BreadthGrow, t, i, r::Value) = Execute(m, BreadthGrow(s.start, vcat(s.queue, Entry(i, t, Value[s.queue[1].results..., r]))), m.matcher, t, i)

# discard what we have yielded and grown

function failure(k::Config, m::Breadth, s::BreadthGrow)
    if (length(s.queue) > 1)
        execute(k, m, BreadthYield(s.start, s.queue[2:end]), arbitrary(s))
    else
        FAILURE
    end
end

# largest first (greedy) repetition without backtracking of child
# matchers

@auto_hash_equals mutable struct Depth!<:Repeat_
    name::Symbol
    matcher::Matcher
    lo::Integer
    hi::Integer
    flatten::Bool
    Depth!(m, lo, hi; flatten=true) = new(:Depth!, m, lo, hi, flatten)
end

mutable struct DepthSlurp!{I}<:RepeatState
    result::Vector{Value}
    iters::Vector{I}
end
hash(::DepthSlurp!, h::UInt) = throw(CacheException())

arbitrary(s::DepthSlurp!) = s.iters[1]

@auto_hash_equals struct DepthYield!{I}<:RepeatState
    result::Vector{Value}
    iters::Vector{I}
end

execute(k::Config{S,I}, m::Depth!, s::Clean, i::I) where {S,I} = execute(k, m, DepthSlurp!{I}(Value[], I[i]), i)

function execute(k::Config, m::Depth!, s::DepthSlurp!, i)
    if length(s.result) < m.hi
        Execute(m, s, m.matcher, CLEAN, i)
    else
        execute(k, m, DepthYield!(s.result, s.iters), i)
    end
end

function success(k::Config, m::Depth!, s::DepthSlurp!, t::State, i, r::Value)
    push!(s.result, r)
    push!(s.iters, i)
    execute(k, m, s, i)
end

failure(k::Config, m::Depth!, s::DepthSlurp!) = execute(k, m, DepthYield!(s.result, s.iters), arbitrary(s))

function execute(k::Config, m::Depth!, s::DepthYield!, unused)
    if length(s.result) == m.lo
        Success(DIRTY, s.iters[end], repeat_success(m, s.result))
    else
        Success(DepthYield!(s.result[1:end-1], s.iters[1:end-1]), s.iters[end], repeat_success(m, s.result))
    end
end

# smallest first (non-greedy) repetition without backtracking of child
# matchers

@auto_hash_equals mutable struct Breadth!<:Repeat_
    name::Symbol
    matcher::Matcher
    lo::Integer
    hi::Integer
    flatten::Bool
    Breadth!(m, lo, hi; flatten=true) = new(:Breadth!, m, lo, hi, flatten)
end

@auto_hash_equals struct BreadthState!{I}<:RepeatState
    result::Vector{Value}
    iter::I
end

function execute(k::Config, m::Breadth!, s::Clean, i)
    if m.lo == 0
        Success(BreadthState!(Value[], i), i, EMPTY)
    else
        execute(k, m, BreadthState!(Value[], i), i)
    end
end

function execute(k::Config, m::Breadth!, s::BreadthState!, unused)
    if length(s.result) == m.hi
        FAILURE
    else
        Execute(m, s, m.matcher, CLEAN, s.iter)
    end
end

function success(k::Config, m::Breadth!, s::BreadthState!, t::State, i, r::Value)
    result = Value[s.result..., r]
    if length(result) >= m.lo
        Success(BreadthState!(result, i), i, repeat_success(m, result))
    else
        execute(k, m, BreadthState!(result, i), i)
    end
end

failure(k::Config, m::Breadth!, s::BreadthState!) = FAILURE

# match all in a sequence with backtracking

# there are two nearly identical matchers here - the only difference is
# whether results are merged (Seq/+) or not (And/&).

# we need two different types so that we can define + and & appropriately.
# to make the user API more conssistent we also define Series (similar to
# Repeat) that takes a flatten argument.  finally, both are so similar
# that they can share the same state.

function Series(m::Matcher...; flatten=true, backtrack=true)
    if flatten
        backtrack ? Seq(m...) : Seq!(m...)
    else
        backtrack ? And(m...) : And!(m...)
    end
end

abstract type Series_<:Matcher end

# first, the backtracking version

@auto_hash_equals mutable struct Seq<:Series_
    name::Symbol
    matchers::Vector{Matcher}
    Seq(m::Matcher...) = new(:Seq, [m...])
    Seq(m::Vector{Matcher}) = new(:Seq, m)
end

serial_success(m::Seq, results::Vector{Value}) = flatten(results)

@auto_hash_equals mutable struct And<:Series_
    name::Symbol
    matchers::Vector{Matcher}
    And(m::Matcher...) = new(:And, Matcher[m...])
    And(m::Vector{Matcher}) = new(:And, m)
end

# copy to get type right (Array{Value,1} -> Array{Any,1})
serial_success(m::And, results::Vector{Value}) = Any[results...]

@auto_hash_equals struct SeriesState{I}<:State
    results::Vector{Value}
    iters::Vector{I}
    states::Vector{State}
end

# when first called, call first matcher

function execute(k::Config, m::Series_, s::Clean, i)
    if length(m.matchers) == 0
        Success(DIRTY, i, EMPTY)
    else
        Execute(m, SeriesState(Value[], [i], State[]), m.matchers[1], CLEAN, i)
    end
end

# if the final matcher matched then return what we have.  otherwise, evaluate
# the next.

function success(k::Config, m::Series_, s::SeriesState, t, i, r::Value)
    n = length(s.iters)
    results = Value[s.results..., r]
    iters = vcat(s.iters, i)
    states = vcat(s.states, t)
    if n == length(m.matchers)
        Success(SeriesState(results, iters, states), i, serial_success(m, results))
    else
        Execute(m, SeriesState(results, iters, states), m.matchers[n+1], CLEAN, i)
    end
end

# if the first matcher failed, fail.  otherwise backtrack

function failure(k::Config, m::Series_, s::SeriesState)
    n = length(s.iters)
    if n == 1
        FAILURE
    else
        Execute(m, SeriesState(s.results[1:end-1], s.iters[1:end-1], s.states[1:end-1]), m.matchers[n-1], s.states[end], s.iters[end-1])
    end
end

# try to advance the current match

function execute(k::Config, m::Series_, s::SeriesState, i)
    @assert length(s.states) == length(m.matchers)
    Execute(m, SeriesState(s.results[1:end-1], s.iters[1:end-1], s.states[1:end-1]), m.matchers[end], s.states[end], s.iters[end-1])
end

# next, the non-backtracking version

abstract type Series!<:Matcher end

@auto_hash_equals mutable struct Seq!<:Series!
    name::Symbol
    matchers::Vector{Matcher}
    Seq!(m::Matcher...) = new(:Seq!, Matcher[m...])
    Seq!(m::Vector{Matcher}) = new(:Seq!, m)
end

serial_success(m::Seq!, results::Vector{Value}) = flatten(results)

@auto_hash_equals mutable struct And!<:Series!
    name::Symbol
    matchers::Vector{Matcher}
    And!(m::Matcher...) = new(:And!, Matcher[m...])
    ANd!(m::Vector{Matcher}) = new(:And!, m)
end

serial_success(m::And!, results::Vector{Value}) = Any[results...]

@auto_hash_equals struct SeriesState!<:State
    results::Vector{Value}
    i  # index into current alternative
end

execute(k::Config, m::Series!, s::Clean, i) = execute(k, m, SeriesState!(Value[], 0), i)

function execute(k::Config, m::Series!, s::SeriesState!, i)
    if s.i == length(m.matchers)
        return Success(DIRTY, i, serial_success(m, s.results))
    else
        return Execute(m, SeriesState!(s.results, s.i+1), m.matchers[s.i+1], CLEAN, i)
    end
end

success(k::Config, m::Series!, s::SeriesState!, t, i, r::Value) = execute(k, m, SeriesState!(Value[s.results..., r], s.i), i)

failure(k::Config, m::Series!, s::SeriesState!) = FAILURE

# alternatives (these need to be separate matchers so that | works ok)

function Alternatives(m::Matcher...; backtrack=true)
    backtrack ? Alt(m...) : Alt!(m...)
end

abstract type Alternatives_<:Matcher end

# first, the backtracking version

@auto_hash_equals mutable struct Alt<:Alternatives_
    name::Symbol
    matchers::Vector{Matcher}
    Alt(matchers::Matcher...) = new(:Alt, Matcher[matchers...])
    Alt(matchers::Vector{Matcher}) = new(:Alt, matchers)
end

@auto_hash_equals struct AltState{I}<:State
    state::State
    iter::I
    i::Int  # index into current alternative
end

function execute(k::Config, m::Alt, s::Clean, i)
    if length(m.matchers) == 0
        FAILURE
    else
        execute(k, m, AltState(CLEAN, i, 1), i)
    end
end

function execute(k::Config, m::Alt, s::AltState, i)
    Execute(m, s, m.matchers[s.i], s.state, s.iter)
end

function success(k::Config, m::Alt, s::AltState, t, i, r::Value)
    Success(AltState(t, s.iter, s.i), i, r)
end

function failure(k::Config, m::Alt, s::AltState)
    if s.i == length(m.matchers)
        FAILURE
    else
        execute(k, m, AltState(CLEAN, s.iter, s.i+1), s.iter)
    end
end

# next, the non-backtracking version (for some sense of backtracking -
# it does not re-match children again, but does try other
# alternatives)

@auto_hash_equals mutable struct Alt!<:Alternatives_
    name::Symbol
    matchers::Vector{Matcher}
    Alt!(matchers::Matcher...) = new(:Alt!, Matcher[matchers...])
    Alt!(matchers::Vector{Matcher}) = new(:Alt!, matchers)
end

@auto_hash_equals struct AltState!{I}<:State
    iter::I
    i::Int  # index into current alternative
end

execute(k::Config, m::Alt!, s::Clean, i) = execute(k, m, AltState!(i, 0), i)

function execute(k::Config, m::Alt!, s::AltState!, i)
    if s.i == length(m.matchers)
        FAILURE
    else
        Execute(m, AltState!(i, s.i+1), m.matchers[s.i+1], CLEAN, i)
    end
end

success(k::Config, m::Alt!, s::AltState!, t, i, r::Value) = Success(s, i, r)

failure(k::Config, m::Alt!, s::AltState!) = execute(k, m, AltState!(s.iter, s.i), s.iter)

# evaluate the child, but discard values and do not advance the iter

@auto_hash_equals mutable struct Lookahead<:Delegate
    name::Symbol
    matcher::Matcher
    Lookahead(matcher) = new(:Lookahead, matcher)
end

always_print(::Delegate) = true

@auto_hash_equals struct LookaheadState<:DelegateState
    state::State
    iter
end

execute(k::Config, m::Lookahead, s::Clean, i) = Execute(m, LookaheadState(s, i), m.matcher, CLEAN, i)

success(k::Config, m::Lookahead, s, t, i, r::Value) = Success(LookaheadState(t, s.iter), s.iter, EMPTY)

# if the child matches, fail; if the child fails return EMPTY
# repeated calls fail (the internal matcher is not backtracked).

@auto_hash_equals mutable struct Not<:Matcher
    name::Symbol
    matcher::Matcher
    Not(matcher) = new(:Not, matcher)
end

@auto_hash_equals struct NotState<:State
    iter
end

execute(k::Config, m::Not, s::Clean, i) = Execute(m, NotState(i), m.matcher, CLEAN, i)

execute(k::Config, m::Not, s::NotState, i) = FAILURE

success(k::Config, m::Not, s::NotState, t, i, r::Value) = FAILURE

failure(k::Config, m::Not, s::NotState) = Success(s, s.iter, EMPTY)

# match a regular expression.

# because Regex match against strings, this matcher works only against
# string sources.

# for efficiency, we need to know the offset where the match finishes.
# we do this by adding r"(.??)" to the end of the expression and using
# the offset from that.

# we also prepend ^ to anchor the match

@auto_hash_equals mutable struct Pattern<:Matcher
    name::Symbol
    text::AbstractString
    regex::Regex
    groups::Tuple
    Pattern(r::Regex, group::Int...) = new(:Pattern, r.pattern, Regex("^(?:" * r.pattern * ")(.??)"), group)
    Pattern(s::AbstractString, group::Int...) = new(:Pattern, s, Regex("^(?:" * s * ")(.??)"), group)
    Pattern(s::AbstractString, flags::AbstractString, group::Int...) = new(:Pattern. s, Regex("^(?:" * s * ")(.??)", flags), group)
end

print_field(m::Pattern, ::Type{Val{:text}}) = "text=\"$(m.text)\""
print_field(m::Pattern, ::Type{Val{:regex}}) = "regex=r\"$(m.regex.pattern)\""

function execute(k::Config, m::Pattern, s::Clean, i)
    x = match(m.regex, forwards(k.source, i))
    if x == nothing
        FAILURE
    else
        i = discard(k.source, i, x.offsets[end]-1)
        if length(m.groups) > 0
            Success(DIRTY, i, Any[x.captures[i] for i in m.groups])
        else
            Success(DIRTY, i, Any[x.match])
        end
    end
end

# support loops

mutable struct Delayed<:Matcher
    name::Symbol
    matcher::Nullable{Matcher}
    Delayed() = new(:Delayed, Nullable{Matcher}())
end

function print_matcher(m::Delayed, known::Set{Matcher})
    function producer(c::Channel)
        tag = "$(m.name)"
        if (isnull(m.matcher))
            put!(c, "$(tag) OPEN")
        elseif m in known
            put!(c, "$(tag)...")
        else
            put!(c, "$(tag)")
            push!(known, m)
            for (i, line) in enumerate(print_matcher(get(m.matcher), known))
                put!(c, i == 1 ? "`-$(line)" : "  $(line)")
            end
        end
    end
    Channel(c -> producer(c))
end

function execute(k::Config, m::Delayed, s::Dirty, i)
    Response(DIRTY, i, FAILURE)
end

function execute(k::Config, m::Delayed, s::State, i)
    if isnull(m.matcher)
        error("set Delayed matcher")
    else
        execute(k, get(m.matcher), s, i)
    end
end

# end of stream / string

@auto_hash_equals mutable struct Eos<:Matcher
    name::Symbol
    Eos() = new(:Eos)
end

function execute(k::Config, m::Eos, s::Clean, i)
    if iterate(k.source, i) == nothing
        Success(DIRTY, i, EMPTY)
    else
        FAILURE
    end
end

# this is general, but usually not much use with backtracking

mutable struct ParserError{I}<:Exception
    msg::AbstractString
    iter::I
end

@auto_hash_equals mutable struct Error<:Matcher
    name::Symbol
    msg::AbstractString
    Error(msg::AbstractString) = new(:Error, msg)
end

function execute(k::Config, m::Error, s::Clean, i::I) where {I}
    msg = m.msg
    try
        msg = diagnostic(k.source, i, m.msg)
    catch x
        println("cannot generate diagnostic for $(typeof(k.source))")
        # continue with original message
    end
    throw(ParserError{I}(msg, i))
end


text:
struct WeightedLeastSquares
    estimate::Array{Float64,1}
    rmse::Float64
    mae::Float64
    wrss::Float64
end

struct ErrorMetric
    rmse::Float64
    mae::Float64
    wrss::Float64
end

struct ErrorMetricWiden
    rmse::Float64
    mae::Float64
    wrss::Float64
    rmseGBPWLS::Float64
    maeGBPWLS::Float64
end

########## Compute the GBP error metrics ##########
function errorMetric(gbp::Union{ContinuousModel, ContinuousTreeModel})
    observationGBP = gbp.system.jacobian * gbp.inference.mean
    NactiveRows = activeRows(gbp)

    rmse = 0.0; mae = 0.0; wrss = 0.0
    @inbounds for i = 1:length(gbp.system.observation)
        rmse += (gbp.system.observation[i] - observationGBP[i])^2
        mae += abs(gbp.system.observation[i] - observationGBP[i])
        wrss += (gbp.system.observation[i] - observationGBP[i])^2 / gbp.system.variance[i]
    end
    rmse = (rmse / NactiveRows)^(1/2)
    mae = mae / NactiveRows

    return ErrorMetric(rmse, mae, wrss)
end

########## Compute the GBP error metrics and error according to the WLS ##########
@inline function errorMetric(gbp::Union{ContinuousModel, ContinuousTreeModel}, wls::WeightedLeastSquares)
    observationGBP = gbp.system.jacobian * gbp.inference.mean
    NactiveRows = activeRows(gbp)

    rmse = 0.0; mae = 0.0; wrss = 0.0
    @inbounds for i = 1:length(gbp.system.observation)
        rmse += (gbp.system.observation[i] - observationGBP[i])^2
        mae += abs(gbp.system.observation[i] - observationGBP[i])
        wrss += (gbp.system.observation[i] - observationGBP[i])^2 / gbp.system.variance[i]
    end
    rmse = (rmse / NactiveRows)^(1/2)
    mae = mae / NactiveRows

    rmseGBPWLS = 0.0; maeGBPWLS = 0.0
    @inbounds for i = 1:gbp.graph.Nvariable
        rmseGBPWLS += (wls.estimate[i] - gbp.inference.mean[i])^2
        maeGBPWLS += abs(wls.estimate[i] - gbp.inference.mean[i])
    end
    rmseGBPWLS = (rmseGBPWLS / gbp.graph.Nvariable)^(1/2)
    maeGBPWLS = maeGBPWLS / gbp.graph.Nvariable

    return ErrorMetricWiden(rmse, mae, wrss, rmseGBPWLS, maeGBPWLS)
end

########## Compute WLS solution and error metrics ##########
function wls(gbp::Union{ContinuousModel, ContinuousTreeModel})
    W = spdiagm(0 =>  @. 1.0 / sqrt(gbp.system.variance))
    A = W * gbp.system.jacobian
    G = A' * A
    b = A' * W * gbp.system.observation
    x = G \ b

    observationWLS = gbp.system.jacobian * x
    NactiveRows = activeRows(gbp)

    rmse = 0.0; mae = 0.0; wrss = 0.0
    @inbounds for i = 1:length(gbp.system.observation)
        rmse += (gbp.system.observation[i] - observationWLS[i])^2
        mae += abs(gbp.system.observation[i] - observationWLS[i])
        wrss += (gbp.system.observation[i] - observationWLS[i])^2 / gbp.system.variance[i]
    end
    rmse = (rmse / NactiveRows)^(1/2)
    mae = mae / NactiveRows

    return WeightedLeastSquares(x, rmse, mae, wrss)
end

########## Find number of nonzeros rows ##########
@inline function activeRows(gbp)
    NactiveRows = length(gbp.system.observation)
    @inbounds for (k, i) in enumerate(gbp.system.observation)
        if i == 0.0
            if all(gbp.system.jacobianTranspose[:, k] .== 0.0)
                NactiveRows -= 1
            end
        end
    end

    return NactiveRows
end

########## Check keyword arguments ##########
function checkKeywords(prob, alpha, variance)
    #### Check convergence parameters
    if prob <= 0.0 || prob >= 1.0
        error("Invalid prob value.")
    end
    if alpha <= 0.0 || alpha >= 1.0
        error("Invalid alpha value.")
    end

    #### Check initial variance
    if variance < 0.0
        error("Invalid variance value.")
    end
end

ext: jl


text:
include("scripts.jl")
Random.seed!(7891)

function scr_ml(ARGS)

    # PARSING ARGUMENTS
    # id_scr = 0; n_distr = 10; n_bio = 3; n_sample = 20; mxcb = 100
    # inp = "../data/data_f_cs.csv"; out = "../res"
    id_scr = parse(Int32, ARGS[1])
    # number of biotracers
    n_bio = parse(Int64, ARGS[2])
    # number of indiv used for building distribution
    n_distr = parse(Int64, ARGS[3])
    # number of sample used to test
    n_sample = parse(Int64, ARGS[4])
    # number of replicates
    n_rep = parse(Int64, ARGS[5])
    # maximum column combintaions
    mxcb = parse(Int64, ARGS[6])
    # path to data (pca-transformed or not)
    pca = (parse(Int32, ARGS[7]) == 1)
    println(pca)
    inp = pca ? "../data/data_f_pca.csv" : "../data/data_f_cs.csv"
    if pca
        println("Only one set of columns will be used")
    end

    # FIXED PARAMETERS
    # path to res folder
    # number of regions considered
    n_reg = 3
    n_tot_indiv = 30
    n_tot_bio = 17
    # data augmentation parameters
    aug = Int64(1e3)
    sig = .01

    df = get_data(false, inp)

    out = "res"*ARGS[2]*"_"*ARGS[7]*"/"
    mkdir(out)

    println("Data loaded ✓")

    if id_scr == 1
        # NB n_distr used as n_distr_max
        scr_ndistr(df, n_bio, n_distr, n_sample, mxcb, n_rep, n_tot_bio, n_tot_indiv,
            n_reg, aug, sig, out, pca)
    elseif id_scr == 2
        scr_nsample(df, n_bio, n_distr, n_sample, mxcb, n_rep, n_tot_bio,
        n_tot_indiv, n_reg, aug, sig, out, pca)
    elseif id_scr == 3
        scr_noise(df, n_bio, n_distr, n_sample, mxcb, n_rep, n_tot_bio,
        n_tot_indiv, n_reg, aug, sig, out, pca)
    else
        # NB n_bio used as n_bio_max
        scr_nbio_v2(df, n_bio, n_distr, n_sample, mxcb, n_rep, n_tot_bio,
        n_tot_indiv, n_reg, aug, sig, out, pca)
    end

    return 0
end

println("All scripts loaded ✓")
scr_ml(ARGS)

# df = get_data(false, inp)
# inp = "../data/data_f_cs.csv"
# scr_train(df)


text:
using Test

@timed_testset "sequential/" begin
    include("sequential/generated_weights.jl")
    if Base.find_package("Gurobi") != nothing
        # Skip these tests if Gurobi is not installed.
        # The corresponding networks are too large for CBC to deal with.
        include("sequential/trained_weights.jl")
    end
end


text:
export Weights, create_weights, compute_weights, component,
        extractcoeff, rearrange_weights,
        rearrange, ncoeff

import Base: +, *, size

struct Weights
        Nx::Int64
        p::Union{Int64, Array{Int64,1}}
        Ne::Int64
        woff::Array{Float64,2}
        wdiag::Array{Float64,2}
        w∂k::Array{Float64,2}
end

function Base.show(io::IO, W::Weights)
  println(io,"Weights for a transport map with Nx = $(W.Nx), p = $(W.p) and Ne = $(W.Ne)")
end

function create_weights(T::RadialMap)
        Nx = T.Nx
        p = T.p
        if p==0
                woff  = zeros(Nx-1)
                wdiag = zeros(2*Nx)
                w∂k   = zeros(Nx)
        else
                woff  = zeros((Nx-1)*(p+1))
                wdiag = zeros(Nx*(p+3))
                w∂k   = zeros(Nx*(p+2))
        end
        return woff, wdiag, w∂k
end

function create_weights(T::RadialMap, X::AbstractMatrix{Float64})
        NxX, Ne = size(X)
        Nx = T.Nx
        @assert NxX == Nx "Wrong dimension of the ensemble matrix `X`"
        p = T.p
        if p==0
                woff  = zeros(Ne, Nx-1)
                wdiag = zeros(Ne, 2*Nx)
                w∂k   = zeros(Ne, Nx)
        else
                woff  = zeros(Ne, (Nx-1)*(p+1))
                wdiag = zeros(Ne, Nx*(p+3))
                w∂k   = zeros(Ne, Nx*(p+2))
        end
        return Weights(Nx, p, Ne, woff, wdiag, w∂k)
end

function compute_weights(u::ui, z::Float64, woff)
        p = u.p
        # ChecNx that woff has the right size
        if p==-1
                #No weights
        elseif p==0
                @assert size(woff, 1) == 1
                woff[1] = z
        else
                @assert size(woff, 1) == p+1

                woff[1] = z
                for j=1:p
                woff[j+1] = rbf(u.ξi[j], u.σi[j])(z)
                end
        end
        # woff
end

function compute_weights(u::uk, z::Float64, wdiag, w∂k; withconstant::Bool = true)
        p = u.p
        # Check that wdiag and w∂k have the right size
        if withconstant == true
        if p==-1
                # No weights
        elseif p==0
                @assert size(wdiag, 1) == 2
                @assert size(w∂k, 1) == 1

                wdiag[1:2] .= [1.0; z]
                w∂k[1] = 1.0
        else
                @assert size(wdiag, 1) == p+3
                @assert size(w∂k, 1) == p+2

                wdiag[1] = 1.0
                # ψ₀ and ψ₀′
                wdiag[2] = ψ₀(u.ξk[1], u.σk[1])(z)
                w∂k[1] = ψ₀′(u.ξk[1], u.σk[1])(z)

                # ψₚ₊₁ and ψₚ₊₁′
                wdiag[p+3] = ψpp1(u.ξk[p+2], u.σk[p+2])(z)
                w∂k[p+2] = ψpp1′(u.ξk[p+2], u.σk[p+2])(z)

                # ψⱼ and ψⱼ′
                @inbounds for j=2:p+1
                wdiag[j+1] = ψj(u.ξk[j], u.σk[j])(z)
                w∂k[j] = rbf(u.ξk[j], u.σk[j])(z)
                end
        end
        else
        if p==-1
                # No weights
        elseif p==0
                @assert size(wdiag, 1) == 1
                @assert size(w∂k, 1) == 1

                wdiag[1] = z
                w∂k[1] = 1.0
        else
                @assert size(wdiag, 1) == p+2
                @assert size(w∂k, 1) == p+2

                # ψ₀ and ψ₀′
                wdiag[1] = ψ₀(u.ξk[1], u.σk[1])(z)
                w∂k[1] = ψ₀′(u.ξk[1], u.σk[1])(z)

                # ψₚ₊₁ and ψₚ₊₁′
                wdiag[p+2] = ψpp1(u.ξk[p+2], u.σk[p+2])(z)
                w∂k[p+2] = ψpp1′(u.ξk[p+2], u.σk[p+2])(z)

                # ψⱼ and ψⱼ′
                @inbounds for j=2:p+1
                        wdiag[j] = ψj(u.ξk[j], u.σk[j])(z)
                        w∂k[j] = rbf(u.ξk[j], u.σk[j])(z)
                end
        end
        end

        # wdiag, w∂k
end

function compute_weights(u::uk, z::Float64, wdiag; withconstant::Bool = true)
        p = u.p
        # Check that wdiag and w∂k have the right size
        if withconstant ==true
        if p==-1
                # No weights
        elseif p==0
                @assert size(wdiag, 1) == 2
                wdiag[1:2] .= [1.0; z]
        else
                @assert size(wdiag, 1) == p+3
                wdiag[1] = 1.0
                # ψ₀ and ψ₀′
                wdiag[2] = ψ₀(u.ξk[1], u.σk[1])(z)

                # ψₚ₊₁ and ψₚ₊₁′
                wdiag[p+3] = ψpp1(u.ξk[p+2], u.σk[p+2])(z)

                # ψⱼ and ψⱼ′
                @inbounds for j=2:p+1
                wdiag[j+1] = ψj(u.ξk[j], u.σk[j])(z)
                end
        end
        else
                if p==-1
                        # No weights
                elseif p==0
                        @assert size(wdiag, 1) == 1
                        wdiag[1] = z
                else
                        @assert size(wdiag, 1) == p+2
                        # ψ₀ and ψ₀′
                        wdiag[1] = ψ₀(u.ξk[1], u.σk[1])(z)

                        # ψₚ₊₁ and ψₚ₊₁′
                        wdiag[p+2] = ψpp1(u.ξk[p+2], u.σk[p+2])(z)

                        # ψⱼ and ψⱼ′
                        @inbounds for j=2:p+1
                        wdiag[j] = ψj(u.ξk[j], u.σk[j])(z)
                        end
                end
        end
end

function compute_weights(T::RadialMap, z::AbstractVector{Float64}, woff::AbstractVector{Float64}, wdiag::AbstractVector{Float64}, w∂k::AbstractVector{Float64})
        @get T (Nx, p)
        if p==0
                if Nx==1
                wd = view(wdiag,:)
                w∂ = view(w∂k,Nx:Nx)
                weights(component(T.C[1],1), z[1], wd, w∂)
                else
                # Fill wdiag and w∂k
                @inbounds for i=1:Nx
                        wd = view(wdiag,(i-1)*2+1:i*2)
                        w∂ = view(w∂k,i:i)
                       compute_weights(component(T.C[i],i), z[i], wd, w∂)
                end
                # Fill woff
                @inbounds for i=1:Nx-1
                        wo = view(woff,i:i)
                       compute_weights(component(T.C[end],i), z[i], wo)
                end
                end
        else
                if Nx==1
                wd = view(wdiag,:)
                w∂ = view(w∂k,:)
               compute_weights(component(T.C[1],1), z[1], wd, w∂)
                else
                        # Fill wdiag and w∂k
                        @inbounds for i=1:Nx
                                wd = view(wdiag, (i-1)*(p+3)+1:i*(p+3))
                                w∂ = view(w∂k, (i-1)*(p+2)+1:i*(p+2))
                               compute_weights(component(T.C[i],i), z[i], wd, w∂)
                        end
                        # Fill woff
                        @inbounds for i=1:Nx-1
                                wo = view(woff, (i-1)*(p+1)+1:i*(p+1))
                               compute_weights(component(T.C[end],i), z[i], wo)
                        end
                end
        end
end

function compute_weights(T::RadialMap, X::AbstractMatrix{Float64}, woff::Array{Float64,2}, wdiag::Array{Float64,2}, w∂k::Array{Float64,2})
        NxX, Ne = size(X)
        @get T (Nx, p)

        @assert Nx == NxX "Error in dimension of the ensemble and size of the KR map"

        if p==0
                if Nx==1
                utmp = component(T.C[1],1)
                @inbounds for l=1:Ne
                        wd = view(wdiag,l,:)
                        w∂ = view(w∂k,l,1:1)
                       compute_weights(utmp, X[1,l], wd, w∂)
                end
                else
                # Fill wdiag and w∂k
                @inbounds for i=1:Nx
                utmp = component(T.C[i],i)
                        for l=1:Ne
                        wd = view(wdiag, l, (i-1)*2+1:i*2)
                        w∂ = view(w∂k,l,i:i)
                       compute_weights(utmp, X[i,l], wd, w∂)
                        end
                end
                # Fill woff
                @inbounds  for i=1:Nx-1
                        utmp = component(T.C[end],i)
                        for l=1:Ne
                        wo = view(woff,l, i:i)
                       compute_weights(utmp, X[i,l], wo)
                        end
                end
                end
        else
                if Nx==1
                utmp = component(T.C[1],1)
                @inbounds for l=1:Ne
                        wd = view(wdiag,l,:)
                        w∂ = view(w∂k,l,:)
                       compute_weights(utmp, X[1,l], wd, w∂)
                end
                else
                        # Fill wdiag and w∂k
                        @inbounds for i=1:Nx
                                utmp = component(T.C[i],i)
                                for l=1:Ne
                                wd = view(wdiag, l, (i-1)*(p+3)+1:i*(p+3))
                                w∂ = view(w∂k, l, (i-1)*(p+2)+1:i*(p+2))
                               compute_weights(utmp, X[i,l], wd, w∂)
                                end
                        end
                        # Fill woff
                        @inbounds for i=1:Nx-1
                                utmp = component(T.C[end],i)
                                for l=1:Ne
                                wo = view(woff, l, (i-1)*(p+1)+1:i*(p+1))
                               compute_weights(utmp, X[i,l], wo)
                                end
                        end
                end
        end
end

function compute_weights(T::RadialMap, X, W::Weights)
        Nx, Ne = size(X)
        @assert T.p == W.p "Error value of p, can't return weights"
        @assert T.Nx == W.Nx "Error value of Nx, can't return weights"
        @assert T.Nx == Nx "Error value of Nx, can't return weights"

        return compute_weights(T, X, W.woff, W.wdiag, W.w∂k)
end

# Extract the weights on
function rearrange_weights(W::Weights, l::Int64)
        @get W (Nx, p, Ne, woff, wdiag)
        if p==0
                if l==1
                        wquad = wdiag[:,1:2]
                else
                        wquad = vcat(woff[:,1:l-1], wdiag[:,2*(l-1)+1:2*l])
                end
        else
                if l==1
                        wquad = wdiag[:,1:p+3]
                else
                        wquad = vcat(woff[:,1:(l-1)*(p+1)], wdiag[:,(l-1)*(p+3)+1:l*(p+3)])
                end
        end
        return  wquad#, w∂
end

# This function extracts the right weights for the cost function, wquad will be the part of the weights
# associated with the quadratic function, while w∂ will be associated with the log barrier function
# rearrange_weights does not use the line of 1 in wdiag
function rearrange(W::Weights, l::Int64)
        @get W (Nx, p, Ne, woff, wdiag, w∂k)
        if p==0
                if l==1
                        wo = zeros(0,0)
                        wd = wdiag[:,2:2]
                        w∂ = w∂k[:,1:1]
                else
                        wo = woff[:,1:l-1]
                        wd = wdiag[:,2*l:2*l]
                        w∂ = w∂k[:,l:l]
                end
        else
                if l==1
                        wo = zeros(0,0)
                        wd = wdiag[:,2:p+3]
                        w∂ = w∂k[:,1:p+2]
                else
                        wo = woff[:,1:(l-1)*(p+1)]
                        wd = wdiag[:,(l-1)*(p+3)+2:l*(p+3)]
                        w∂ = w∂k[:,(l-1)*(p+2)+1:l*(p+2)]
                end
        end
        return  wo, wd, w∂
end

## Procedure for a sparse map, weights cannot be reused from one line to the other,
#  need to treat each k-th component individually

function ncoeff(k::Int64, p::Array{Int64})
        noff = 0
        ndiag = 0
        n∂k = 0
        for i=1:k-1
                if p[i]==-1
                elseif p[i]==0
                        noff += 1
                else
                        noff += deepcopy(p[i])+1
                end
        end

        if p[k]==-1
        elseif p[k]==0
                ndiag = 1
                n∂k = 1
        else
                ndiag = p[k]+2
                n∂k = p[k]+2
        end

        return noff, ndiag, n∂k
end

# This function can be called directly into optimize to get the right weights

function compute_weights(C::SparseRadialMapComponent, z::AbstractVector{Float64})
        @get C (Nx, p)
        # Determine the number of coefficients for each ( ncoeff)
        noff, ndiag, n∂k = ncoeff(Nx, p)

        woff  = zeros(noff)
        # The constant of the diagonal term is omitted
        wdiag = zeros(ndiag)
        w∂k   = zeros(n∂k)
        count = 0
        # Fill off-diagonal components
        @inbounds for i=1:Nx-1
                pidx = p[i]
                if pidx == -1
                        # No weights to compute
                elseif pidx == 0
                        # Linear function
                        wo = view(woff,count+1)
                        compute_weights(component(C,i), z[i], wo)
                        count +=1
                else
                        # Linear term + pidx rbf functions
                        wo = view(woff,count+1:count+pidx+1)
                        compute_weights(component(C,i), z[i], wo)
                        count += deepcopy(pidx) + 1
                end
        end
        @assert count == noff "The size is correct for the off-diagonal weights"

        # Fill diagonal components
        if p[Nx] ==-1
                # No weights to compute
        else
                # Affine function for p[k]=0
                # Constant plus (p[k]+2) ψ functions for p[k]>0
                wd = view(wdiag,:)
                w∂ = view(w∂k,:)
               compute_weights(component(C,Nx), z[Nx], wd, w∂,  withconstant = false)
        end
        return woff, wdiag, w∂k
end

function compute_weights(C::SparseRadialMapComponent, X::AbstractMatrix{Float64})
        NxX, Ne = size(X)
        @get C (Nx, p)
        @assert NxX == Nx
        # Determine the number of coefficients for each ( ncoeff)
        noff, ndiag, n∂k = ncoeff(Nx, p)
        if noff ==0
        woff  = zeros(0,0)
        else
        woff  = zeros(Ne, noff)
        end
        # The constant of the diagonal term is omitted
        if ndiag ==0
        wdiag = zeros(0,0)
        else
        wdiag = zeros(Ne, ndiag)
        end
        if n∂k==0
        w∂k = zeros(0,0)
        else
        w∂k   = zeros(Ne, n∂k)
        end

        count = 0
        # Fill off-diagonal components
        @inbounds for i=1:Nx-1
                pidx = p[i]
                if pidx == -1
                        # No weights to compute
                elseif pidx == 0
                        # Linear function
                        for l=1:Ne
                                wo = view(woff,l:l,count+1)
                                compute_weights(component(C,i), X[i,l], wo)
                        end
                        count +=1
                else
                        for l=1:Ne
                                # Linear term + pidx rbf functions
                                wo = view(woff,l,count+1:count+pidx+1)
                                compute_weights(component(C,i), X[i,l], wo)
                        end
                        count += deepcopy(pidx) + 1
                end
        end
        @assert count == noff "The size is correct for the off-diagonal weights"

        # Fill diagonal components
        if p[Nx] ==-1
                # No weights to compute
        else
                for l=1:Ne
                # Affine function for p[Nx]=0
                # Constant plus (p[Nx]+2) ψ functions for p[Nx]>0
                        wd = view(wdiag,l,:)
                        w∂ = view(w∂k,l,:)
                        compute_weights(component(C,Nx), X[Nx,l], wd, w∂, withconstant = false)
                end
        end
        return woff, wdiag, w∂k
end


text:
module BIOES

function decode(ids::Vector{Int}, tagdict::Dict{String,Int})
    id2tag = Array{String}(length(tagdict))
    for (k,v) in tagdict
        id2tag[v] = k
    end

    spans = Tuple{Int,Int,String}[]
    bpos = 0
    for i = 1:length(ids)
        tag = id2tag[ids[i]]
        tag == "O" && continue
        startswith(tag,"B") && (bpos = i)
        startswith(tag,"S") && (bpos = i)
        nexttag = i == length(ids) ? "O" : id2tag[ids[i+1]]
        if (startswith(tag,"S") || startswith(tag,"E")) && bpos > 0
            tag = id2tag[ids[bpos]]
            basetag = length(tag) > 2 ? tag[3:end] : ""
            push!(spans, (bpos,i,basetag))
            bpos = 0
        end
    end
    spans
end

end


text:
"""
The functions provided by YAXArrays are supposed to work on different types of cubes. This module defines the interface for all
Data types that
"""
module Cubes
using DiskArrays: DiskArrays, eachchunk
using Distributed: myid
using Dates: TimeType
using IntervalSets: Interval, (..)
using Base.Iterators: take, drop
using ..YAXArrays: workdir, YAXDefaults
using YAXArrayBase: YAXArrayBase, iscompressed, dimnames, iscontdimval
import YAXArrayBase: getattributes, iscontdim, dimnames, dimvals, getdata

export concatenatecubes, caxes, subsetcube, readcubedata, renameaxis!, YAXArray

"""
This function calculates a subset of a cube's data
"""
function subsetcube end

"Returns the axes of a Cube"
function caxes end

include("Axes.jl")
using .Axes:
    CubeAxis,
    RangeAxis,
    CategoricalAxis,
    findAxis,
    getAxis,
    axVal2Index,
    axname,
    axsym,
    axVal2Index_lb,
    axVal2Index_ub,
    renameaxis,
    axcopy

mutable struct CleanMe
    path::String
    persist::Bool
    function CleanMe(path::String, persist::Bool)
        c = new(path, persist)
        finalizer(clean, c)
        c
    end
end
function clean(c::CleanMe)
    if !c.persist && myid() == 1
        if !isdir(c.path) && !isfile(c.path)
            @warn "Cube directory $(c.path) does not exist. Can not clean"
        else
            rm(c.path, recursive = true)
        end
    end
end

"""
    YAXArray{T,N}

An array labelled with named axes that have values associated with them.
It can wrap normal arrays or, more typically DiskArrays.

### Fields

* `axes` a `Vector{CubeAxis}` containing the Axes of the Cube
* `data` N-D array containing the data
"""
struct YAXArray{T,N,A<:AbstractArray{T,N},AT}
    axes::AT
    data::A
    properties::Dict{String}
    cleaner::Vector{CleanMe}
    function YAXArray(axes, data, properties, cleaner)
        if ndims(data) != length(axes)
            throw(
                ArgumentError(
                    "Can not construct YAXArray, supplied data dimension is $(ndims(data)) while the number of axes is $(length(axes))",
                ),
            )
        elseif ntuple(i -> length(axes[i]), ndims(data)) != size(data)
            throw(
                ArgumentError(
                    "Can not construct YAXArray, supplied data size is $(size(data)) while axis lenghts are $(ntuple(i->length(axes[i]),ndims(data)))",
                ),
            )
        else
            return new{eltype(data),ndims(data),typeof(data),typeof(axes)}(
                axes,
                data,
                properties,
                cleaner,
            )
        end
    end
end

YAXArray(axes, data, properties = Dict{String,Any}(); cleaner = CleanMe[]) =
    YAXArray(axes, data, properties, cleaner)
function YAXArray(x::AbstractArray)
    ax = caxes(x)
    props = getattributes(x)
    YAXArray(ax, x, props)
end
Base.size(a::YAXArray) = size(getdata(a))
Base.size(a::YAXArray, i::Int) = size(getdata(a), i)
function Base.getproperty(a::YAXArray, s::Symbol)
    ax = axsym.(caxes(a))
    i = findfirst(isequal(s), ax)
    if i === nothing
        return getfield(a, s)
    else
        return caxes(a)[i]
    end
end
function Base.propertynames(a::YAXArray, private::Bool = false)
    if private
        (axsym.(caxes(a))..., :axes, :data, :properties)
    else
        (axsym.(caxes(a))..., :axes, :data)
    end
end

Base.ndims(a::YAXArray{<:Any,N}) where {N} = N
Base.eltype(a::YAXArray{T}) where {T} = T
Base.permutedims(c::YAXArray, p) =
    YAXArray(caxes(c)[collect(p)], permutedims(getdata(c), p), c.properties, c.cleaner)
caxes(c::YAXArray) = getfield(c, :axes)
function caxes(x)
    map(enumerate(dimnames(x))) do a
        i, s = a
        v = YAXArrayBase.dimvals(x, i)
        iscontdim(x, i) ? RangeAxis(s, v) : CategoricalAxis(s, v)
    end
end

"""
    readcubedata(cube)

Given any array implementing the YAXArray interface it returns an in-memory [`YAXArray`](@ref) from it.
"""
function readcubedata(x)
    YAXArray(collect(CubeAxis, caxes(x)), getindex_all(x), getattributes(x))
end

cubechunks(c::YAXArray) = common_size(eachchunk(getdata(c)))
cubechunks(x) = common_size(eachchunk(x))
common_size(a::DiskArrays.GridChunks) = a.chunksize
function common_size(a)
    ntuple(ndims(first(a))) do idim
        otherdims = setdiff(1:ndims(first(a)), idim)
        allengths = map(i -> length(i.indices[idim]), a)
        for od in otherdims
            allengths = unique(allengths, dims = od)
        end
        length(allengths) < 3 ? allengths[1] : allengths[2]
    end
end

getindex_all(a) = getindex(a, ntuple(_ -> Colon(), ndims(a))...)
Base.getindex(x::YAXArray, i...) = getdata(x)[i...]
chunkoffset(c::YAXArray) = common_offset(eachchunk(getdata(c)))
chunkoffset(x) = common_offset(eachchunk(x))
common_offset(a::DiskArrays.GridChunks) = a.offset
function common_offset(a)
    ntuple(ndims(a)) do idim
        otherdims = setdiff(1:ndims(a), idim)
        allengths = map(i -> length(i[idim]), a)
        for od in otherdims
            allengths = unique(allengths, dims = od)
        end
        @assert length(allengths) == size(a, idim)
        length(allengths) < 3 ? 0 : allengths[2] - allengths[1]
    end
end

# Implementation for YAXArrayBase interface
YAXArrayBase.dimvals(x::YAXArray, i) = caxes(x)[i].values

function YAXArrayBase.dimname(x::YAXArray, i)
    axsym(caxes(x)[i])
end

YAXArrayBase.getattributes(x::YAXArray) = x.properties

YAXArrayBase.iscontdim(x::YAXArray, i) = isa(caxes(x)[i], RangeAxis)

YAXArrayBase.getdata(x::YAXArray) = getfield(x, :data)

function YAXArrayBase.yaxcreate(::Type{YAXArray}, data, dimnames, dimvals, atts)
    axlist = map(dimnames, dimvals) do dn, dv
        iscontdimval(dv) ? RangeAxis(dn, dv) : CategoricalAxis(dn, dv)
    end
    YAXArray(axlist, data, atts)
end
YAXArrayBase.iscompressed(c::YAXArray) = _iscompressed(getdata(c))
_iscompressed(c::DiskArrays.PermutedDiskArray) = iscompressed(c.a.parent)
_iscompressed(c::DiskArrays.SubDiskArray) = iscompressed(c.v.parent)
_iscompressed(c) = YAXArrayBase.iscompressed(c)

function renameaxis!(c::YAXArray, p::Pair)
    axlist = caxes(c)
    i = findAxis(p[1], axlist)
    axlist[i] = renameaxis(axlist[i], p[2])
    c
end
function renameaxis!(c::YAXArray, p::Pair{<:Any,<:CubeAxis})
    i = findAxis(p[1], caxes(c))
    i === nothing && throw(ArgumentError("$(p[1]) Axis not found"))
    length(caxes(c)[i].values) == length(p[2].values) ||
        throw(ArgumentError("Length of replacement axis must equal length of old axis"))
    caxes(c)[i] = p[2]
    c
end

function _subsetcube end

function subsetcube(z::YAXArray{T}; kwargs...) where {T}
    newaxes, substuple = _subsetcube(z, collect(Any, map(Base.OneTo, size(z))); kwargs...)
    newdata = view(getdata(z), substuple...)
    YAXArray(newaxes, newdata, z.properties, cleaner = z.cleaner)
end

sorted(x, y) = x < y ? (x, y) : (y, x)

#TODO move everything that is subset-related to its own file or to axes.jl
interpretsubset(subexpr::Union{CartesianIndices{1},LinearIndices{1}}, ax) =
    subexpr.indices[1]
interpretsubset(subexpr::CartesianIndex{1}, ax) = subexpr.I[1]
interpretsubset(subexpr, ax) = axVal2Index(ax, subexpr, fuzzy = true)
function interpretsubset(subexpr::NTuple{2,Any}, ax)
    x, y = sorted(subexpr...)
    Colon()(sorted(axVal2Index_lb(ax, x), axVal2Index_ub(ax, y))...)
end
interpretsubset(subexpr::NTuple{2,Int}, ax::RangeAxis{T}) where {T<:TimeType} =
    interpretsubset(map(T, subexpr), ax)
interpretsubset(subexpr::UnitRange{<:Integer}, ax::RangeAxis{T}) where {T<:TimeType} =
    interpretsubset(T(first(subexpr))..T(last(subexpr) + 1), ax)
interpretsubset(subexpr::Interval, ax) = interpretsubset((subexpr.left, subexpr.right), ax)
interpretsubset(subexpr::AbstractVector, ax::CategoricalAxis) =
    axVal2Index.(Ref(ax), subexpr, fuzzy = true)

function _subsetcube(z, subs; kwargs...)
    kwargs = Dict{Any,Any}(kwargs)
    for f in YAXDefaults.subsetextensions
        f(kwargs)
    end
    newaxes = deepcopy(collect(caxes(z)))
    foreach(kwargs) do kw
        axdes, subexpr = kw
        axdes = string(axdes)
        iax = findAxis(axdes, caxes(z))
        if isa(iax, Nothing)
            throw(ArgumentError("Axis $axdes not found in cube"))
        else
            oldax = newaxes[iax]
            subinds = interpretsubset(subexpr, oldax)
            subs2 = subs[iax][subinds]
            subs[iax] = subs2
            if !isa(subinds, AbstractVector) && !isa(subinds, AbstractRange)
                newaxes[iax] = axcopy(oldax, oldax.values[subinds:subinds])
            else
                newaxes[iax] = axcopy(oldax, oldax.values[subinds])
            end
        end
    end
    substuple = ntuple(i -> subs[i], length(subs))
    inewaxes = findall(i -> isa(i, AbstractVector), substuple)
    newaxes = newaxes[inewaxes]
    @assert length.(newaxes) ==
            map(length, filter(i -> isa(i, AbstractVector), collect(substuple)))
    newaxes, substuple
end

Base.getindex(a::YAXArray; kwargs...) = subsetcube(a; kwargs...)

Base.read(d::YAXArray) = getindex_all(d)

function formatbytes(x)
    exts = ["bytes", "KB", "MB", "GB", "TB"]
    i = 1
    while x >= 1024
        i = i + 1
        x = x / 1024
    end
    return string(round(x, digits = 2), " ", exts[i])
end
cubesize(c::YAXArray{T}) where {T} = (sizeof(T)) * prod(map(length, caxes(c)))
cubesize(::YAXArray{T,0}) where {T} = sizeof(T)

getCubeDes(::CubeAxis) = "Cube axis"
getCubeDes(::YAXArray) = "YAXArray"
getCubeDes(::Type{T}) where {T} = string(T)
Base.show(io::IO, c::YAXArray) = show_yax(io, c)

function show_yax(io::IO, c)
    println(io, getCubeDes(c), " with the following dimensions")
    for a in caxes(c)
        println(io, a)
    end
    foreach(getattributes(c)) do p
        if p[1] in ("labels", "name", "units")
            println(io, p[1], ": ", p[2])
        end
    end
    println(io, "Total size: ", formatbytes(cubesize(c)))
end

include("TransformedCubes.jl")
include("Slices.jl")
end


text:
# Code qui determine un jeu de poids soit local a 2 points, soit global a tous les points

using PyPlot
using LinearAlgebra, Printf
function testidee()
    c1 = rand(0:100,20)
    c2 = rand(0:100,20)
    x  = rand(0.0:1.0,20)
    for i=1:5
        x[rand(1:20)]=rand()
    end
    Lfrac=[]
    for i=1:20
        if x[i]!=0.0 && x[i]!=1.0
            push!(Lfrac,i)
        end
    end
    @show x
    @show Lfrac

    # --------------------------------------------------------------------------
    P=[]
    push!(P,[rand(1:200) rand(800:1000) ])
    while P[end][1]<900 && P[end][2]>100
        push!(P,[rand(P[end][1]+1:min(1000,P[end][1]+500)) rand(max(1,P[end][2]-500):P[end][2]-1)])
    end
    # ---
    fig = figure("test idees",figsize=(8,8))
    xlim(0,1500)
    ylim(0,1500)
    for i=1:length(P)
        scatter(P[i][1],P[i][2], color="black")
    end
    i=rand(1:length(P)-1)
    a=P[i]
    x1,y1 = a
    b=P[i+1]
    x2,y2 = b
    Δx = abs(x1-x2)
    Δy = abs(y1-y2)
    λ1 =  1 - Δx / (Δx+Δy)
    λ2 =  1 - Δy / (Δx+Δy)
    @show P
    @show a
    @show b
    @show Δx
    @show Δy
    println("directions deduites localement")
    @show λ1
    @show λ2
    xlim(0,1000)
    ylim(0,1000)
    scatter(x1, y1, color="red", marker="x")
    scatter(x2, y2, color="red", marker="x")
    zλ = λ1 * dot(c1,x) + λ2 * dot(c2,x)
    #----
    n1=P[end][1]
    n2=P[1][2]
    xm=(x1+x2)/2.0
    ym=(y1+y2)/2.0
    Δx = abs(n1-xm)
    Δy = abs(n2-ym)
    λ1 =  1 - Δx / (Δx+Δy)
    λ2 =  1 - Δy / (Δx+Δy)
    @show Δx
    @show Δy
    println("directions deduites globalement")
    @show λ1
    @show λ2
    plot(n1, n2, xm, ym, linestyle="-", color="blue", marker="+")
    annotate("",
	xy=[xm;ym],# Arrow tip
	xytext=[n1;n2], # Text offset from tip
	arrowprops=Dict("arrowstyle"=>"->"))

    # --------------------------------------------------------------------------
    for f in Lfrac  @printf("%5d ",f) end ; @printf("\n")
    for f in Lfrac  @printf("%5d ",c1[f]) end; @printf("\n")
    for f in Lfrac  @printf("%5d ",c2[f]) end; @printf("\n")
    for f in Lfrac  @printf("%5.2f ",λ1*c1[f]+λ2*c2[f]) end; @printf("\n")
    for f in Lfrac  @printf("%5.2f ",x[f]) end; @printf("\n")

end

testidee()

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
# # Pandemic Control
# A stochastic optimal control problem were we seek to control the spread of a
# contagion with uncertain parameters.

# ## Problem Statement and Model
# We wish to determine optimal social distancing policy to control spread of a
# contagion and minimize the economic impact of doing so. We'll model the spread
# of the virus through a given population using the SEIR Model which considers 4
# population subsets that follow:
# ```math
# \text{Susceptible} \rightarrow \text{Exposed} \rightarrow \text{Infectious} \rightarrow \text{Recovered}
# ```
# This model is formalized as:
# ```math
# \begin{gathered}
# \frac{ds(t)}{dt} = (u(t) - 1)\beta si(t) \\
# \frac{de(t)}{dt} = (1 - u(t))\beta si(t) - \xi e(t) \\
# \frac{di(t)}{dt} = \xi e(t) - \gamma i(t)\\
# \frac{dr(t)}{dt} = \gamma i(t) \\
# si(t) = s(t) i(t)
# \end{gathered}
# ```
# where ``s(t)`` is the susceptible population, ``e(t)`` is the exposed population,
# ``i(t)`` is the infectious population, ``r(t)`` is the recovered population,
# and ``u(t) \in [0, 1]`` is the enforced population isolation (social distancing).
# The other values denote model parameters that will be specific to the contagion
# in question.

# For our case study, we'll consider the incubation constant ``\xi`` to be an
# uncertain parameter ``\xi \sim \mathcal{U}(\underline{\xi}, \overline{\xi})``.
# This introduces to infinite dependencies into our model: time ``t`` and
# uncertain incubation ``\xi``. Moreover we'll seek to minimize the isolation
# measures ``u(t)`` that are implemented while limiting the amount of infected
# individuals ``i(t)`` to be below a threshold ``i_{max}``. Thus, the optimization
# problem becomes:
# ```math
# \begin{aligned}
# &&\min_{} &&& \int_{t \in \mathcal{D}_{t}} u(t) dt \\
# && \text{s.t.} &&& \frac{\partial s(t, \xi)}{\partial t} = (u(t) - 1)\beta si(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
# &&&&& \frac{\partial e(t, \xi)}{\partial t} = (1 - u(t))\beta si(t, \xi) - \xi e(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
# &&&&& \frac{\partial i(t, \xi)}{\partial t} = \xi e(t, \xi) - \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
# &&&&& \frac{\partial r(t, \xi)}{\partial t} = \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
# &&&&& si(t, \xi) = s(t, \xi) i(t, \xi), && \forall \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
# &&&&& s(0, \xi) = s_0, e(0, \xi) = e_0, i(0, \xi) = i_0, r(0, \xi) = r_0, && \forall \xi \in \mathcal{D}_{\xi} \\
# &&&&& i(t, \xi) \leq i_{max}, && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
# &&&&& u(t) \in [0, 0.8] \\
# &&&&& \xi \sim \mathcal{U}(\underline{\xi}, \overline{\xi})
# \end{aligned}
# ```
# Notice that the SEIR model variables now all depend on both ``t`` and ``\xi``,
# except ``u(t)`` which means we need to decide our policy before knowing the
# true value of ``\xi``.

# ## Modeling in InfiniteOpt
# Traditional modeling frameworks like JuMP cannot be used to solve this model
# directly because it is infinite, contains partial difference equations, and
# contains a time expectation. We would first have to transform it into a
# transcripted (discretized) variant by applying all the necessary techniques
# (e.g., orthogonal collocation over finite elements, trapezoid rule, etc.)
# which is very combersome and nontrivial in this case.

# However, we can directly model the above form in InfiniteOpt and it will take
# care of the rest! Let's get started by importing the needed packages and
# specifying the model parameters that we'll need.

using InfiniteOpt, Ipopt, Distributions, Plots

## Set the SEIR parameters
γ = 0.303
β = 0.727
N = 1e5
ξ_min = 0.1
ξ_max = 0.6

## Set the domain information
i_max = 0.02
ϵ = 0.005
t0 = 0
tf = 200
extra_ts = [0.001, 0.002, 0.004, 0.008, 0.02, 0.04, 0.08, 0.2, 0.4, 0.8]
num_samples = 5

## Set the intial condition values
e0 = 1 / N
i0 = 0
r0 = 0
s0 = 1 - 1 / N;

# ## Model Initialization
# Now let's setup our infinite model and select Ipopt as our optimizer that will
# be used to solve it. This is accomplished making an [`InfiniteModel`](@ref):
model = InfiniteModel(Ipopt.Optimizer);

# ## Infinite Parameter Definition
# We now need to define the infinite parameters ``t \in [t_0, t_f]`` and
# ``\xi \sim \mathcal{U}(\underline{\xi}, \overline{\xi})``. This is accomplished
# with [`@infinite_parameter`](@ref). We'll also include the following
# specifications:
# - use 51 equidistant time points
# - specify that orgothonal collocation using 2 nodes should be used for time derivatives
# - specify that the number of random scenarios should equal `num_samples`
# - add `extra_ts` as extra time points
@infinite_parameter(model, t ∈ [t0, tf], num_supports = 51,
                    derivative_method = OrthogonalCollocation(2))
@infinite_parameter(model, ξ ~ Uniform(ξ_min, ξ_max), num_supports = num_samples)
add_supports(t, extra_ts)

# ## Infinite Variable Definition
# With our infinite parameters defined, we can now define our infinite variables:
# - ``s(t, \xi) \geq 0``
# - ``e(t, \xi) \geq 0``
# - ``i(t, \xi) \geq 0``
# - ``r(t, \xi) \geq 0``
# - ``si(t, \xi)``
# - ``0 \leq u(t) \leq 0.8``
@variable(model, s ≥ 0, Infinite(t, ξ))
@variable(model, e ≥ 0, Infinite(t, ξ))
@variable(model, i ≥ 0, Infinite(t, ξ))
@variable(model, r ≥ 0, Infinite(t, ξ))
@variable(model, si, Infinite(t, ξ))
@variable(model, 0 ≤ u ≤ 0.8, Infinite(t), start = 0.2)

# ## Objective Definition
# Now its time to add the objective
# ``\min \ \int_{t \in \mathcal{D}_{t}} u(t) dt`` using `@objective`:
@objective(model, Min, ∫(u, t))

# ## Constraint Definition
# The last step now to defining our model is to define the constraints using
# `@constraint`. This will involve defining the initial conditions:
# ```math
# s(0, \xi) = s_0, e(0, \xi) = e_0, i(0, \xi) = i_0, r(0, \xi) = r_0, \ \forall \xi \in \mathcal{D}_{\xi}
# ```
# the model equations:
# ```math
# \begin{aligned}
# &&& \frac{\partial s(t, \xi)}{\partial t} = (u(t) - 1)\beta si(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
# &&& \frac{\partial e(t, \xi)}{\partial t} = (1 - u(t))\beta si(t, \xi) - \xi e(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
# &&& \frac{\partial i(t, \xi)}{\partial t} = \xi e(t, \xi) - \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
# &&& \frac{\partial r(t, \xi)}{\partial t} = \gamma i(t, \xi), && \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi} \\
# &&& si(t, \xi) = s(t, \xi) i(t, \xi), && \forall \forall t \in \mathcal{D}_{t}, \xi \in \mathcal{D}_{\xi}, \\
# \end{aligned}
# ```
# and the infection limit constraint:
# ```math
# i(t, \xi) \leq i_{max}, \ \forall t \in \mathcal{D}_t, \xi \in \mathcal{D}_{\xi}.
# ```
## Define the initial conditions
@constraint(model, s == s0, DomainRestrictions(t => 0))
@constraint(model, e == e0, DomainRestrictions(t => 0))
@constraint(model, i == i0, DomainRestrictions(t => 0))
@constraint(model, r == r0, DomainRestrictions(t => 0))

## Define the SEIR equations
@constraint(model, s_constr, ∂(s, t) == -(1 - u) * β * si)
@constraint(model, e_constr, ∂(e, t) == (1 - u) * β * si - ξ * e)
@constraint(model, i_constr, ∂(i, t) == ξ * e - γ * i)
@constraint(model, r_constr, ∂(r, t) == γ * i)
@constraint(model, si == s * i)

## Define the infection rate limit
@constraint(model, imax_constr, i ≤ i_max)

# ## Display the Infinite Model
# Let's display `model` now that it is fully defined:
print(model)

# ## Optimize the Model
# Let's solve our model and find the optimal policy. All we have to do is
# invoke `optimize!` and the model will automatically be transformed solved behind
# the scenes:
optimize!(model)

# ## Retrieve and Plot the Results
# Now we can retrieve the optimal results and plot them to visualize the optimal
# policy. Note that the values of infinite variables will be returned as arrays
# corresponding to how the supports were used to discretize our model. We can
# retrieve our values using `value`.

# Get the results:
r_opt = value(r, ndarray = true) * 100 # make the population fractions into percentages
s_opt = value(s, ndarray = true) * 100
i_opt = value(i, ndarray = true) * 100
e_opt = value(e, ndarray = true) * 100
u_opt = value(u)
obj_opt = objective_value(model)
ts = value(t)
ξs = value(ξ);

# Plot the values of ``r(t, \xi)`` and ``s(t, \xi)`` over time with
# confidence bounds:
r_mean = mean(r_opt, dims = 2)
r_std = std(r_opt, dims = 2)
plot(ts, r_mean, label = "r(t, ξ)", linecolor = :red, background_color = :transparent)
plot!(ts, r_mean + r_std, linecolor = :red, linestyle = :dash, linealpha = 0.4, label = "")
plot!(ts, r_mean - r_std, linecolor = :red, linestyle = :dash, linealpha = 0.4, label = "")

s_mean = mean(s_opt, dims = 2)
s_std = std(s_opt, dims = 2)
plot!(ts, s_mean, label = "s(t, ξ)", linecolor = :blue)
plot!(ts, s_mean + s_std, linecolor = :blue, linestyle = :dash, linealpha = 0.4, label = "")
plot!(ts, s_mean - s_std, linecolor = :blue, linestyle = :dash, linealpha = 0.4, label = "")
ylabel!("Pop. (%)")
xlabel!("Time (Days)")

# Plot the values of ``i(t, \xi)`` and ``e(t, \xi)`` over time with
# confidence bounds:
i_mean = mean(i_opt, dims = 2)
i_std = std(i_opt, dims = 2)
plot(ts, i_mean, label = "i(t, ξ)", linecolor = :green, background_color = :transparent)
plot!(ts, i_mean + i_std, linecolor = :green, linestyle = :dash, linealpha = 0.4, label = "")
plot!(ts, i_mean - i_std, linecolor = :green, linestyle = :dash, linealpha = 0.4, label = "")

e_mean = mean(e_opt, dims = 2)
e_std = std(e_opt, dims = 2)
plot!(ts, e_mean, label = "e(t, ξ)", linecolor = :purple)
plot!(ts, e_mean + e_std, linecolor = :purple, linestyle = :dash, linealpha = 0.4, label = "")
plot!(ts, e_mean - e_std, linecolor = :purple, linestyle = :dash, linealpha = 0.4, label = "")
ylabel!("Pop. (%)")
xlabel!("Time (Days)")

# Plot the values of ``u(t)`` over time:
plot(ts, u_opt, linecolor = :orange, label = "u(t)", ylims = (-0.02, 1.02),
     background_color = :transparent)
xlabel!("Time (Days)")
ylabel!("Distancing Ratio")

# ### Maintenance Tests
# These are here to ensure this example stays up to date.
using Test
@test termination_status(model) == MOI.LOCALLY_SOLVED
@test has_values(model)
@test u_opt isa Vector{<:Real}


text:
using FixedPointNumbers, Base.Test

for f in ["normed.jl", "fixed.jl"]
    println("Testing $f")
    include(f)
end

if v"0.5.0" <= VERSION < v"0.6.0-dev"
    @test isempty(detect_ambiguities(FixedPointNumbers, Base, Core))
end


text:
using YaoBase, SparseArrays, BitBasis, YaoArrayRegister
using ..SymEngine
export @ket_str, @bra_str
export SymReg, AdjointSymReg, SymRegOrAdjointSymReg, expand
export szero_state

YaoArrayRegister._warn_type(raw::AbstractArray{Basic}) = nothing

const SymReg{B,MT} = ArrayReg{B,Basic,MT} where {MT<:AbstractMatrix{Basic}}
const AdjointSymReg{B,MT} = AdjointArrayReg{B,Basic,MT}
const SymRegOrAdjointSymReg{B,MT} = Union{SymReg{B,MT},AdjointSymReg{B,MT}}

function SymReg{B,MT}(r::ArrayReg{B,<:Number}) where {B,MT<:AbstractMatrix{Basic}}
    return ArrayReg{B,Basic,MT}(MT(Basic.(r.state)))
end

function SymReg(r::ArrayReg{B,<:Number}) where {B}
    smat = SparseMatrixCSC(Basic.(_pretty_basic.(r.state)))
    return ArrayReg{B,Basic,SparseMatrixCSC{Basic,Int}}(SparseMatrixCSC(smat))
end

_pretty_basic(x) = x
_pretty_basic(x::Real) = isinteger(x) ? Int(x) : x
_pretty_basic(x::Complex) = isreal(x) ? _pretty_basic(real(x)) : x

function ket_m(s)
    v, N = parse_str(s)
    st = spzeros(Basic, 1 << N, 1)
    st[v+1] = 1
    return ArrayReg{1}(st)
end

const MAX_SYM_QUBITS = 10

function Base.show(io::IO, r::SymRegOrAdjointSymReg{1})
    if nqubits(r) < MAX_SYM_QUBITS
        print_sym_state(io, r)
    else
        summary(io, r)
        print(io, "\n    active qubits: ", nactive(r), "/", nqubits(r))
    end
end

Base.:(*)(x::SymReg{B,MT}, y::SymReg{B,MT}) where {B,MT} = SymReg{B,MT}(kron(state(x), state(y)))
Base.:(^)(x::SymReg{B,MT}, n::Int) where {B,MT} = SymReg{B,MT}(kron(state(x) for _ in 1:n))

Base.:(*)(x::AdjointSymReg{B,MT}, y::AdjointSymReg{B,MT}) where {B,MT} =
    adjoint(parent(x) * parent(y))
Base.:(^)(x::AdjointSymReg{B,MT}, n::Int) where {B,MT} = adjoint(parent(x)^n)

SymEngine.expand(x::SymReg{B}) where {B} = ArrayReg{B}(expand.(state(x)))

"""
    szero_state(n; nbatch=1)

Create a symbolic zero state, same as `ket"000"`, but allows you use an integer.
"""
szero_state(args...; kwargs...) = zero_state(Basic, args...; kwargs...)

function YaoBase.partial_tr(r::SymReg{B}, locs) where {B}
    orders = setdiff(1:nqubits(r), locs)
    focus!(r, orders)
    state = sum(rank3(r); dims = 2)
    relax!(r, orders)
    return ArrayReg(state)
end


text:
#Contain the Zone_target structure to load hold and manipulate the targets in a zone
"""loads the targets in the zone and keep track of whether a newzone in needed to be loaded
"""
struct Zone_Target
    #Loads the targets relvant for a particular zone
    #Any thing stored in two element array will always asume body first and head second
    max_target::Integer #maximum number of target
    max_usp::Integer #total maximum number of unique priority for the full tile

    zones::Array{Integer,1}
    #first entry for actual number of target and second entry for number of unique priority
    ntarget::Array{Integer,1}

    index::Array{Integer,1}
    X_FP::Array{Float64,1} #To store the x-position on focal plane of target
    Y_FP::Array{Float64,1} #To store y-position on focal plane of target
    PRIORITY::Array{Integer,1}

    #boolean array
    pot_ass::Array{Bool,1} #To set for the potential assignment
    pot_coll::Array{Bool,1} #To set for potential collision

    #To store theta,phi for a particular positioner
    theta_pos::Array{Float64,1}
    phi_pos::Array{Float64,1}

    #unique sorted priority for this tile
    usp::Array{Integer,1}
    usp_bool_sel::Array{Integer,2} #To store the boolean selection of each priority

    #To store the possible assignments
    possible_ass::Array{Int32,1}

    function  Zone_Target(maximum_target::Integer,nunique_priority::Integer)
        new(maximum_target,nunique_priority,zeros(3),zeros(2),zeros(maximum_target),zeros(maximum_target),zeros(maximum_target),
        zeros(maximum_target),zeros(Bool,maximum_target),zeros(Bool,maximum_target),
        zeros(maximum_target),zeros(maximum_target),
        zeros(nunique_priority),zeros(Bool,maximum_target,nunique_priority),
        zeros(Int32,maximum_target))
    end
end

function load_target_zone(zone_target::Zone_Target,zone_x,zone_y,fp_bound,targets_dic)
    #This zone is already loaded so nothing to be done
    if((zone_target.zones[1]==zone_x) & (zone_target.zones[2]==zone_y))
        return
    else
        zone_target.zones[1]=zone_x
        zone_target.zones[2]=zone_y
        zone_target.zones[3]=get_zonexy(zone_x,zone_y,fp_bound)
    end

    #indsel=zeros(Int32,1000)

    count_targ=1
    for zx in maximum([zone_x-1,1]):minimum([zone_x+1,fp_bound.nx])
        for zy in maximum([zone_y-1,1]):minimum([zone_y+1,fp_bound.ny])
            zone_xy=get_zonexy(zx,zy,fp_bound)

            if(!(zone_xy in targets_dic["zone_map_keys"]))
                continue
            end
            imin=targets_dic["zone_map"][zone_xy][1]
            imax=targets_dic["zone_map"][zone_xy][2]
            #@show zx,zy,zone_xy,targets_dic["zone_map"][zone_xy]
            #@show targets_dic["zone_x"][imin:imax]
            count_this=imax-imin
            zone_target.index[count_targ:count_targ+count_this] .= collect(imin:imax)
            count_targ=count_targ+count_this+1
        end
    end

    #set everything else to 0
    if(count_targ<=zone_target.max_target)
        zone_target.index[count_targ:end] .= 0
    end
    #Due to Julia indexing convetion 1 need to be subtracte
    count_targ=count_targ-1

    tmp_ind=view(zone_target.index,1:count_targ)
    #sort the indices
    sort!(tmp_ind)

    zone_target.ntarget[1]=count_targ

    zone_target.X_FP[1:count_targ]=targets_dic["X_FP"][tmp_ind]
    zone_target.Y_FP[1:count_targ]=targets_dic["Y_FP"][tmp_ind]
    zone_target.PRIORITY[1:count_targ]=targets_dic["PRIORITY"][tmp_ind]

    us_priority=sort(unique(zone_target.PRIORITY[1:zone_target.ntarget[1]]),rev=true)
    nup=size(us_priority,1)
    zone_target.ntarget[2]=nup
    zone_target.usp[1:nup] .= us_priority

    imin=1;imax=zone_target.ntarget[1]
    for uu in 1:zone_target.ntarget[2]
        zone_target.usp_bool_sel[imin:imax,uu] .= (zone_target.PRIORITY[imin:imax].== us_priority[uu])
    end

end

max_stars_repo_name: shadaba/MockFBA

max_issues_repo_name: shadaba/MockFBA

max_forks_repo_name: shadaba/MockFBA


text:
"Non deterministic variable"
abstract type NonDetVar end


text:
import .Tree

struct WeightedGraph <: Graph
    idx::IncidenceIndex
    lam::Vector{Float64}
end

@inline num_nodes(g::WeightedGraph) =
    num_nodes(g.idx)

@inline num_edges(g::WeightedGraph) =
    num_edges(g.idx)

@inline function enumerate_edges(f::F, g::WeightedGraph) where {F<:Function}
    for u in 1:num_nodes(g)
        for (v, ei) in g.idx[u]
            if u < v
                f(ei, u, v, g.lam[ei])
            end
        end
    end
end

@inline function enumerate_edges(f::F, t::Tree.WeightedTree{W}) where {F<:Function, W}
    for (u, v) in enumerate(t.tree.parent)
        if u == t.tree.root
            continue
        end
        f(u, u, v, t.lam[u])
    end
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
function mm_spec_kj3_dimspec()::MPCCDimSpec
    dimspec = MPCCDimSpec(2, 1, 2, 0, 0, 1, 0)
    return dimspec
end

function mm_spec_kj3_defn(x, pr, ps)

    f = (x[1] - pr[1])^2 + (x[2] - pr[1])^2

    ce = Vector{Num}()

    ci = Vector{Num}()

    F = Matrix{Num}(undef, 2, 1)
    F[1,1] = x[1]
    F[2,1] = x[2]

    return MPCCDefinition(f, ce, ci, F)
end

function mm_spec_kj3_nzmask()::MPCCModelNZMask
    nz_mask = MPCCModelNZMask(
        # ce
        [],
        # ci
        [   Set{Int64}([1, 2]),
            Set{Int64}([1]) ],
        # F
        [
            # F1
            [   Set{Int64}([1]) ],
            # F2
            [   Set{Int64}([2]) ]
        ] )
    return nz_mask
end

function mm_spec_kj3_testvectors()
    return Vector{MPCCModelTestVector}(undef, 0)
end

function mm_spec_kj3_knownsols()

    function local_kj3_knownsol_1(pr::Vector{R}, ps::Vector{Int64}) where {R <: Real}
        @assert ( 1 == length(pr) && 0 == length(ps) )
        t = pr[1]

        if ( t <= 0 )
            x_sol = [0, 0]
            f_sol = 2(-t)^2
        else
            x_sol = [t, 0]
            f_sol = (x_sol[1] - t)^2 + (x_sol[2] - t)^2
        end

        return (x_sol, f_sol)
    end

    function local_kj3_knownsol_2(pr::Vector{R}, ps::Vector{Int64}) where {R <: Real}
        @assert ( 1 == length(pr) && 0 == length(ps) )
        t = pr[1]

        if ( t <= 0 )
            x_sol = [0, 0]
            f_sol = 2(-t)^2
        else
            x_sol = [0, t]
            f_sol = (x_sol[1] - t)^2 + (x_sol[2] - t)^2
        end

        return (x_sol, f_sol)
    end

    return [ local_kj3_knownsol_1, local_kj3_knownsol_2 ]
end

function mm_spec_kj3_parameterisations(t)
    defn1 = MPCCParameterisationDefn(
            Vector{Num}([ t ]),
            (-5.0, 5.0),
            "Standard"
        )

    return  Vector{MPCCParameterisationDefn}( [ defn1 ] )
end


text:
"""
For converting the concentrations into molecules

Arguments:  cell_radius = the radius of the modelled cell (m)
            baseline,
            indicator,
            endogeneous,
            immobile,
            peak,
            f_i,
            f_e,
            f_im
Returns:    m_baseline, m_indicator, m_endogeneous, m_immobile, m_peak, m_f_i,
            m_f_e, m_f_im
"""
function convertUnits(cell_radius::Float64, baseline::Float64,  indicator::Float64, endogeneous::Float64,
                      immobile::Float64,    peak::Float64,      f_i::Float64,       f_e::Float64,
                      f_im::Float64)
  cell_volume = sphereVolume(cell_radius)
  m_baseline = molarsToMolecules(baseline, cell_volume)
  m_indicator = molarsToMolecules(indicator, cell_volume)
  m_endogeneous = molarsToMolecules(endogeneous, cell_volume)
  m_immobile = molarsToMolecules(immobile, cell_volume)
  m_peak = molarsToMolecules(peak, cell_volume)
  m_f_i = perMolarToPerMolecule(f_i, cell_volume)
  m_f_e = perMolarToPerMolecule(f_e, cell_volume)
  m_f_im = perMolarToPerMolecule(f_im, cell_volume)
  return cell_volume, m_baseline, m_indicator, m_endogeneous, m_immobile, m_peak, m_f_i, m_f_e, m_f_im
end


text:
function finite_diff(f::Function, args::Tuple, i::Int, dx::Float64)
    pos_args = Any[args...]
    pos_args[i] += dx
    neg_args = Any[args...]
    neg_args[i] -= dx
    return (f(pos_args...) - f(neg_args...)) / (2. * dx)
end

const dx = 1e-6

@testset "bernoulli" begin
    f = (x::Bool, prob::Float64) -> logpdf(Gen.Bernoulli(), x, prob)
    args = (false, 0.3,)
    actual = logpdf_grad(Gen.Bernoulli(), args...)
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    args = (true, 0.3,)
    actual = logpdf_grad(Gen.Bernoulli(), args...)
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
end

@testset "beta" begin
    f = (x, alpha, beta) -> logpdf(Gen.Beta(), x, alpha, beta)
    args = (0.4, 0.2, 0.3)
    actual = logpdf_grad(Gen.Beta(), args...)
    @test isapprox(actual[1], finite_diff(f, args, 1, dx))
    @test isapprox(actual[2], finite_diff(f, args, 2, dx))
    @test isapprox(actual[3], finite_diff(f, args, 3, dx))
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
using Images: Gray, channelview, rawview
using Colors
using FileIO
import Printf

using Plots

function load_image_as_grayscale(imagename)
    oriimg = load(imagename)
    grayimg = Gray.(oriimg)
    rawview(channelview(grayimg)) / 255
end

function demo_scalespace(akaze)
    filenames = map(1:length(akaze.evolution_)) do n "akazetest-$(lpad(n,2,'0')).png" end
    for (fn, x) in zip(filenames, akaze.evolution_)
        save(fn, Gray.(x.Lt))
    end
    (j,k) = size(akaze.evolution_)
    run(`montage -filter point -geometry 512x512 -tile $(j)x$k $filenames akazetest.png`)
    imshow(load("akazetest.png"))
end

function plot_features(kpts, cmap="Blues", style=:dash)
    mycirc = hcat([[cos(t), sin(t)] for t in (2 * π * (0:53) / 53)]...)

    data = map(kpts) do kp
        xy = [kp.pt.x; kp.pt.y]
        xx = xy .+ kp.size * [0 cos(kp.angle); 0 sin(kp.angle)]
        circ = xy .+ kp.size * mycirc
        [circ[1,:], circ[2,:], xx[1,:], xx[2,:], xx[1,1:1], xx[2,1:1]]
    end
    cx = [d[1] for d in data]
    cy = [d[2] for d in data]
    rx = [d[3] for d in data]
    ry = [d[4] for d in data]
    px = [d[5] for d in data]
    py = [d[6] for d in data]

    plot!(cx, cy, l=2, color = Colors.colormap(cmap)[90], label = "", style=style, alpha=0.85)
    plot!(rx, ry, l=2, color = Colors.colormap(cmap)[75], label = "", alpha=0.85)
    scatter!(px,py, color = Colors.colormap(cmap)[50], shape = :o, m=5, label = "", msw=0)
end

function original_akaze_features(imagename; diffusivity=0.001, nsublevels=4, omax=4, dthreshold=0.001)
    dthreshold_str = Printf.@sprintf("%.15f", dthreshold)

    cmd = `/home/user/src/akaze/build/bin/akaze_features $imagename --diffusivity $(Int(diffusivity)) --show_results 0 --dthreshold $dthreshold_str --descriptor 5 --nsublevels $nsublevels --omax $omax`
    display(cmd)
    run(cmd)

    data = open("keypoints.txt") do x
        el = eachline(x)
        iterate(el)
        iterate(el)
        ncols = 4
        reshape([parse(Float64, st) for r in el for st in split(r)[1:ncols]], ncols, :)
    end
    keypoints = map(eachcol(data)) do (x,y,size,angle)
        KeyPoint(pt=Point(x,y), size=size, angle=angle)
    end
    desc = open("keypoints.txt") do x
        el = eachline(x)
        rows = parse(Int, iterate(el)[1])
        cols = parse(Int, iterate(el)[1])
        reshape([parse(UInt8, st) for r in el for st in split(r)[end-rows+1:end]], rows, cols)
    end

    keypoints, desc
end

function dump_keypoints_text(kpts, desc)
    thedata = map(kpts, eachcol(desc)) do kp, dd
        @Printf.sprintf("%f %f %f %f %s", kp.pt.x, kp.pt.y, kp.size, kp.angle, join(dd, " "))
    end
    [[@Printf.sprintf("%d",x) for x in size(desc)]...,thedata...]
end

function print_keypoints_for_benchmark(kpts, desc)
    @Printf.printf("%d\n%d", size(desc)...)
    foreach(kpts, eachcol(desc)) do kp, dd
        sc = 1.0 / (kp.size)^2
        @Printf.printf("%f %f %f 0 %f %s\n", kp.pt.x, kp.pt.y, sc, sc, join(dd, " "))
    end
end


text:
__precompile__()

module GNSSTools

    include("greeting.jl")
    export logo, greeting
    # printstyled(logo, color=:default, bold=false)
    # println(greeting)

    # using Reexport

    using Random
    using FFTW
    using Base.Threads
    # Set the number of threads for FFTW to use. It is set to the number of
    # threads that Julia is set to use. User should launch Julia with the number
    # of threads they wish to use.
    FFTW.set_num_threads(nthreads())
    # FFTW.set_num_threads(1)
    using ProgressMeter
    using Statistics
    using LinearAlgebra
    using CPUTime
    using SatelliteToolbox
    using Interpolations
    using QuadGK
    using HDF5
    using StatsBase
    using SpecialFunctions
    using Optim
    using PyPlot

    pygui(true)

    include("constants.jl")
    include("signal_types.jl")
    include("general_tools.jl")
    include("generate_phase_noise.jl")
    include("calc_code_val.jl")
    include("definesignal.jl")
    include("generatesignal.jl")
    include("codes/l1ca_code_generator.jl")
    include("codes/l5_code_generator.jl")
    include("gnss_data_tools.jl")
    include("course_acquisition_tools.jl")
    include("fine_acquisition_tools.jl")
    include("rinex_tools.jl")
    include("tracking_tools.jl")
    include("dkalman.jl")
    include("calcdoppler.jl")
    include("orbit_tools.jl")
    include("tle_tools.jl")
    include("data_process.jl")
    include("doppler_tools.jl")
    include("processing_tools.jl")
    include("demos.jl")
    include("benchmarks.jl")

    export Rₑ
    export L1_freq
    export h_parms_tcxo, h_parms_ocxo, h_parms_cesium
    export l1ca_code_length
    export l1ca_db_chipping_rate
    export l1ca_chipping_rate
    export l1ca_codes, define_l1ca_code_type
    export L5_freq
    export L5_code_length
    export nh10_code_length
    export nh20_code_length
    export L5_chipping_rate
    export nh_chipping_rate
    export L5_db_chipping_rate
    export l5q_codes
    export nh10
    export nh20
    export l5i_codes
    export definesignal, definesignal!, definereplica, definereplica!
    export generatesignal!, generatereplica!
    export fft_correlate
    export GNSSData
    export loaddata, reloaddata!, signal_to_sc_data_file
    export gencorrresult
    export courseacquisition, courseacquisition!
    export v_t2p_fa, p_fa2v_t, false_alarm_pdf, detection_pdf, v_t2p_d, p_d2v_t
    export fineacquisition
    export gnsstypes
    export calcinitcodephase
    export calccodeidx
    export calctvector, calc_t_at_i
    export calc_snr
    export loadrinex
    export trackprn
    export calcA, calcC, calcQ
    export plotresults
    export dkalman
    export demo
    export courseacquisitiontest
    export meshgrid
    export data_info_from_name
    export initorbitinfo
    export calcdoppler, calcdopplerrate
    export getGPSSatnums
    export getTLEs
    export GPSData
    export dataprocess
    export doppler2chips, get_chips_and_ϕ
    export define_constellation
    export calcelevation
    export course_acq_est
    export plot_satellite_orbit
    export doppler_distribution
    export definecodetype, definesignaltype
    export get_eop
    export generate_phase_noise, generate_phase_noise!, plot_spectrum
    export calc_code_val
    export process
    export CN0_monte_carlo
    export define_l5_code_type, define_l1ca_code_type
    export get_distribution_bounds
    export binary2OneOrNegativeOne
    export make3dplot
    export v_t2p_fa, p_fa2v_t, v_t2p_d, p_d2v_t
    export noncoherent_loss, processing_gain
    export phase_noise_variance

end # module


text:
using Base.Test

using DataReaders

@testset "google_daily" begin
    dr = DataReader("google");

    dt_start = DateTime("2015-04-01");
    dt_end = DateTime("2015-04-15");

    @testset "one symbol" begin
        symb = DataSymbol("MSFT");

        response = get(dr, symb, dt_start, dt_end);

        ## DataFrame
        df = DataFrame(response)
        println(df)

        ## TimeArray
        ta_price, ta_volume = TimeArray(response)
        println(ta_price)
        println(ta_volume)
    end

    @testset "several symbols" begin
        symbols = DataSymbol.(["IBM", "MSFT"])

        response = get(dr, symbols, dt_start, dt_end);

        data = DataFrame(response);
        println(data)
    end
end


text:
sbmlfile = joinpath("data", "ecoli_core.xml")
download_data_file(
    "http://systemsbiology.ucsd.edu/sites/systemsbiology.ucsd.edu/files/Attachments/Images/downloads/Ecoli_core/ecoli_core_model.xml",
    sbmlfile,
    "78692f8509fb36534f4f9b6ade23b23552044f3ecd8b48d84d484636922ae907",
)

@testset "Conversion from and to SBML model" begin
    sbmlm = load_sbml_model(sbmlfile)
    sm = convert(StandardModel, sbmlm)
    sbmlm2 = convert(SBMLModel, sm)

    @test Set(reactions(sbmlm)) == Set(reactions(sbmlm2))
    @test Set(reactions(sbmlm)) == Set(reactions(sm))
    @test Set(metabolites(sbmlm)) == Set(metabolites(sbmlm2))
    @test all([
        sbmlm.sbml.reactions[i].stoichiometry == sbmlm2.sbml.reactions[i].stoichiometry for
        i in reactions(sbmlm2)
    ])
end


text:
"""
cgs_efield!(incident_electric_field, fft_tmp, G, vector_potential_3d, dielectric, vector_potential_1d, AToE, Ig, pfft, pifft, x, p, r, rt, u, v, q, uq; tol = 1e-15, maxit = 18)\n
Inplace version for solving the VIE. the result will be in x.\n
"""
function cgs_efield!(
    incident_electric_field,
    fft_tmp,
    G,
    vector_potential_3d,
    dielectric,
    vector_potential_1d,
    AToE,
    Ig,
    pfft,
    pifft,
    x,
    p,
    r,
    rt,
    u,
    v,
    q,
    uq;
    tol = 1e-15,
    maxit = 18
)
    T = ComplexF32
    n2b = norm(incident_electric_field)
    n2b = n2b*tol
    #copyto!(dst,src)
    copyto!(x, incident_electric_field)
    _total_minus_scattered_electric_field!(p, x, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, pfft, pifft, -one(T))
    # p = b - Ax
    axpy!(one(T), incident_electric_field, p)
    # r= b - Ax
    copyto!(r, p)
    # u = b - Ax
    copyto!(u, p)
    # v = A(b - Ax)
    _total_minus_scattered_electric_field!(v, p, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, pfft, pifft, one(T))
    copyto!(rt, r)
    ρ = dot(rt, r)
    for _ in 1:maxit
        σ = dot(rt, v)
        α = ρ/σ
        # q = u-α*v
        caxpy!(q, -α, v, u)
        # uq = u + q
        caxpy!(uq, one(T), u, q)
        # x = x + α*uq
        axpy!(α, uq, x)

        _total_minus_scattered_electric_field_for_residual!(r, uq, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, pfft, pifft, -α)

        if norm(r) <= n2b
          break #return xx,xy,xz
        end

        ρold = ρ
        ρ = dot(rt, r)
        β = ρ/ρold
        # u = r + β*q
        caxpy!(u, β, q, r)
        #p = u + β*(q+β*p)
        axpby!(one(T), q, β, p)
        axpby!(one(T), u, β, p)
        _total_minus_scattered_electric_field!(v, p, vector_potential_1d, vector_potential_3d, G, dielectric, Ig, AToE, fft_tmp, pfft, pifft, one(T))
    end

    return nothing
end


text:
using DamagedShearBand
using Documenter

makedocs(;
    modules=[DamagedShearBand],
    authors="Léo Petit <leo.petit@live.fr> and contributors",
    repo="https://github.com/Leooop/DamagedShearBand.jl/blob/{commit}{path}#L{line}",
    sitename="DamagedShearBand.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Leooop.github.io/DamagedShearBand.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Leooop/DamagedShearBand.jl",
)


text:
using Documenter, CmdStan

makedocs(
    modules = [CmdStan],
    format = :html,
    sitename = "StanJulia/CmdStan.jl",
    pages = Any[
        "Home" => "INTRO.md",
        "Installation" => "INSTALLATION.md",
        "Walkthrough" => "WALKTHROUGH.md",
        "Versions" => "VERSIONS.md",
        "Index" => "index.md"
    ]
)

deploydocs(
    repo = "github.com/StanJulia/CmdStan.jl.git",
 )


text:
struct RozemaAnisotropicMinimumDissipation{FT, K} <: AbstractAnisotropicMinimumDissipation{FT}
     C :: FT
    Cb :: FT
     ν :: FT
     κ :: K
    function RozemaAnisotropicMinimumDissipation{FT}(C, Cb, ν, κ) where FT
        return new{FT, typeof(κ)}(C, Cb, ν, convert_diffusivity(FT, κ))
    end
end

"""
    RozemaAnisotropicMinimumDissipation(FT=Float64; C=0.33, ν=1.05e-6, κ=1.46e-7)

Returns a `RozemaAnisotropicMinimumDissipation` closure object of type `FT` with

    * `C` : Poincaré constant
    * `ν` : 'molecular' background viscosity
    * `κ` : 'molecular' background diffusivity for each tracer

See Rozema et al., " (2015)
"""
RozemaAnisotropicMinimumDissipation(FT=Float64; C=0.33, Cb=0.0, ν=ν₀, κ=κ₀) =
    RozemaAnisotropicMinimumDissipation{FT}(C, Cb, ν, κ)

function with_tracers(tracers, closure::RozemaAnisotropicMinimumDissipation{FT}) where FT
    κ = tracer_diffusivities(tracers, closure.κ)
    return RozemaAnisotropicMinimumDissipation{FT}(closure.C, closure.Cb, closure.ν, κ)
end

# Bindings
const RAMD = RozemaAnisotropicMinimumDissipation

"Return the filter width for Anisotropic Minimum Dissipation on a Regular Cartesian grid."
@inline Δx(i, j, k, grid::RegularCartesianGrid, ::RAMD) = grid.Δx
@inline Δy(i, j, k, grid::RegularCartesianGrid, ::RAMD) = grid.Δy
@inline Δz(i, j, k, grid::RegularCartesianGrid, ::RAMD) = grid.Δz

# We only have regular grids for now. When we have non-regular grids this will need to be changed.
const Δx_ccc = Δx
const Δy_ccc = Δy
const Δz_ccc = Δz
const Δx_ccf = Δx
const Δy_ccf = Δy
const Δz_ccf = Δz

function TurbulentDiffusivities(arch::AbstractArchitecture, grid::AbstractGrid, tracers, ::RAMD)
    νₑ = CellField(arch, grid)
    κₑ = TracerFields(arch, grid, tracers)
    return (νₑ=νₑ, κₑ=κₑ)
end

@inline function ν_ccc(i, j, k, grid::AbstractGrid{FT}, closure::RAMD, buoyancy, U, C) where FT
    q = tr_∇u_ccc(i, j, k, grid, U.u, U.v, U.w)

    if q == 0
        νˢᵍˢ = zero(FT)
    else
        r = Δ²ₐ_uᵢₐ_uⱼₐ_Σᵢⱼ_ccc(i, j, k, grid, closure, U.u, U.v, U.w)
        ζ = Δ²ᵢ_wᵢ_bᵢ_ccc(i, j, k, grid, closure, buoyancy, U.w, C)
        νˢᵍˢ = -closure.C * (r - closure.Cb * ζ) / q
    end

    return max(zero(FT), νˢᵍˢ) + closure.ν
end

@inline function κ_ccc(i, j, k, grid::AbstractGrid{FT}, closure::RAMD, c, ::Val{tracer_index},
                       U) where {FT, tracer_index}

    @inbounds κ = closure.κ[tracer_index]

    σ = θᵢ²_ccc(i, j, k, grid, c)

    if σ == 0
        κˢᵍˢ = zero(FT)
    else
        ϑ =  Δ²ⱼ_uᵢⱼ_cⱼ_cᵢ_ccc(i, j, k, grid, closure, U.u, U.v, U.w, c)
        κˢᵍˢ = - closure.C * ϑ / σ
    end

    return max(zero(FT), κˢᵍˢ) + κ
end

#####
##### The *** 30 terms *** of AMD
#####

@inline function Δ²ₐ_uᵢₐ_uⱼₐ_Σᵢⱼ_ccc(i, j, k, grid, closure, u, v, w)
    Δx = Δx_ccc(i, j, k, grid, closure)
    Δy = Δy_ccc(i, j, k, grid, closure)
    Δz = Δz_ccc(i, j, k, grid, closure)

    ijk = (i, j, k, grid)
    uvw = (u, v, w)
    ijkuvw = (i, j, k, grid, u, v, w)

    Δx²_uᵢ₁_uⱼ₁_Σ₁ⱼ = Δx^2 * (
         Σ₁₁(ijkuvw...) * ∂x_u(ijk..., u)^2
      +  Σ₂₂(ijkuvw...) * ▶xy_cca(ijk..., ∂x_v², uvw...)
      +  Σ₃₃(ijkuvw...) * ▶xz_cac(ijk..., ∂x_w², uvw...)

      +  2 * ∂x_u(ijkuvw...) * ▶xy_cca(ijk..., ∂x_v_Σ₁₂, uvw...)
      +  2 * ∂x_u(ijkuvw...) * ▶xz_cac(ijk..., ∂x_w_Σ₁₃, uvw...)
      +  2 * ▶xy_cca(ijk..., ∂x_v, uvw...) * ▶xz_cac(ijk..., ∂x_w, uvw...) * ▶yz_acc(ijk..., Σ₂₃, uvw...)
    )

    Δy²_uᵢ₂_uⱼ₂_Σ₂ⱼ = Δy^2 * (
      + Σ₁₁(ijkuvw...) * ▶xy_cca(ijk..., ∂y_u², uvw...)
      + Σ₂₂(ijkuvw...) * ∂y_v(ijk..., v)^2
      + Σ₃₃(ijkuvw...) * ▶yz_acc(ijk..., ∂y_w², uvw...)

      +  2 * ∂y_v(ijkuvw...) * ▶xy_cca(ijk..., ∂y_u_Σ₁₂, uvw...)
      +  2 * ▶xy_cca(ijk..., ∂y_u, uvw...) * ▶yz_acc(ijk..., ∂y_w, uvw...) * ▶xz_cac(ijk..., Σ₁₃, uvw...)
      +  2 * ∂y_v(ijkuvw...) * ▶yz_acc(ijk..., ∂y_w_Σ₂₃, uvw...)
    )

    Δz²_uᵢ₃_uⱼ₃_Σ₃ⱼ = Δz^2 * (
      + Σ₁₁(ijkuvw...) * ▶xz_cac(ijk..., ∂z_u², uvw...)
      + Σ₂₂(ijkuvw...) * ▶yz_acc(ijk..., ∂z_v², uvw...)
      + Σ₃₃(ijkuvw...) * ∂z_w(ijk..., w)^2

      +  2 * ▶xz_cac(ijk..., ∂z_u, uvw...) * ▶yz_acc(ijk..., ∂z_v, uvw...) * ▶xy_cca(ijk..., Σ₁₂, uvw...)
      +  2 * ∂z_w(ijkuvw...) * ▶xz_cac(ijk..., ∂z_u_Σ₁₃, uvw...)
      +  2 * ∂z_w(ijkuvw...) * ▶yz_acc(ijk..., ∂z_v_Σ₂₃, uvw...)
    )

    return Δx²_uᵢ₁_uⱼ₁_Σ₁ⱼ + Δy²_uᵢ₂_uⱼ₂_Σ₂ⱼ + Δz²_uᵢ₃_uⱼ₃_Σ₃ⱼ
end

#####
##### trace(∇u) = uᵢⱼ uᵢⱼ
#####

@inline function tr_∇u_ccc(i, j, k, grid, uvw...)
    ijk = (i, j, k, grid)

    return (
        # ccc
        ∂x_u²(ijk..., uvw...)
      + ∂y_v²(ijk..., uvw...)
      + ∂z_w²(ijk..., uvw...)

        # ffc
      + ▶xy_cca(ijk..., ∂x_v², uvw...)
      + ▶xy_cca(ijk..., ∂y_u², uvw...)

        # fcf
      + ▶xz_cac(ijk..., ∂x_w², uvw...)
      + ▶xz_cac(ijk..., ∂z_u², uvw...)

        # cff
      + ▶yz_acc(ijk..., ∂y_w², uvw...)
      + ▶yz_acc(ijk..., ∂z_v², uvw...)
    )
end

@inline function Δ²ᵢ_wᵢ_bᵢ_ccc(i, j, k, grid, closure, buoyancy, w, C)
    ijk = (i, j, k, grid)

    Δx = Δx_ccc(ijk..., closure)
    Δy = Δy_ccc(ijk..., closure)
    Δz = Δz_ccc(ijk..., closure)

    Δx²_wx_bx = Δx^2 * (▶xz_cac(ijk..., ∂x_faa, w)
                          * ▶x_caa(ijk..., ∂x_faa, buoyancy_perturbation, buoyancy, C))

    Δy²_wy_by = Δy^2 * (▶yz_acc(ijk..., ∂y_afa, w)
                          * ▶y_aca(ijk..., ∂y_afa, buoyancy_perturbation, buoyancy, C))

    Δz²_wz_bz = Δz^2 * (∂z_aac(ijk..., w)
                          * ▶z_aac(ijk..., ∂z_aaf, buoyancy_perturbation, buoyancy, C))

    return Δx²_wx_bx + Δy²_wy_by + Δz²_wz_bz
end

@inline function Δ²ⱼ_uᵢⱼ_cⱼ_cᵢ_ccc(i, j, k, grid, closure, u, v, w, c)
    ijk = (i, j, k, grid)

    Δx = Δx_ccc(ijk..., closure)
    Δy = Δy_ccc(ijk..., closure)
    Δz = Δz_ccc(ijk..., closure)

    Δx²_cx_ux = Δx^2 * (
                 ∂x_caa(ijk..., u) * ▶x_caa(ijk..., ∂x_c², c)
        + ▶xy_cca(ijk..., ∂x_v, v) * ▶x_caa(ijk..., ∂x_faa, c) * ▶y_aca(ijk..., ∂y_afa, c)
        + ▶xz_cac(ijk..., ∂x_w, w) * ▶x_caa(ijk..., ∂x_faa, c) * ▶z_aac(ijk..., ∂z_aaf, c)
    )

    Δy²_cy_uy = Δy^2 * (
          ▶xy_cca(ijk..., ∂y_u, u) * ▶y_aca(ijk..., ∂y_afa, c) * ▶x_caa(ijk..., ∂x_faa, c)
        +        ∂y_aca(ijk..., v) * ▶y_aca(ijk..., ∂y_c², c)
        + ▶xz_cac(ijk..., ∂y_w, w) * ▶y_aca(ijk..., ∂y_afa, c) * ▶z_aac(ijk..., ∂z_aaf, c)
    )

    Δz²_cz_uz = Δz^2 * (
          ▶xz_cac(ijk..., ∂z_u, u) * ▶z_aac(ijk..., ∂z_aaf, c) * ▶x_caa(ijk..., ∂x_faa, c)
        + ▶yz_acc(ijk..., ∂z_v, v) * ▶z_aac(ijk..., ∂z_aaf, c) * ▶y_aca(ijk..., ∂y_afa, c)
        +        ∂z_aac(ijk..., w) * ▶z_aac(ijk..., ∂z_c², c)
    )

    return Δx²_cx_ux + Δy²_cy_uy + Δz²_cz_uz
end

@inline θᵢ²_ccc(i, j, k, grid, c) = (
      ▶x_caa(i, j, k, grid, ∂x_c², c)
    + ▶y_aca(i, j, k, grid, ∂y_c², c)
    + ▶z_aac(i, j, k, grid, ∂z_c², c)
)