text:
######### CmdStan program example  ###########

using CmdStan, StanMCMCChain, MCMCChain, Test, JLD, Statistics

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  bernoullimodel = "
  data {
    int<lower=1> N;
    int<lower=0,upper=1> y[N];
  }
  parameters {
    real<lower=0,upper=1> theta;
  }
  model {
    theta ~ beta(1,1);
    y ~ bernoulli(theta);
  }
  "

  observeddata = [
    Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 1, 0, 0, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 0, 0, 0, 1, 0, 1, 1]),
    Dict("N" => 10, "y" => [0, 0, 0, 1, 0, 0, 0, 1, 0, 1])
  ]

  global stanmodel, rc, chains, cnames, d

  stanmodel = Stanmodel(num_samples=1200, thin=1, name="bernoulli",
    model=bernoullimodel, output_format=:mcmcchain);

 @time rc, chains, cnames = stan(stanmodel, observeddata, ProjDir, diagnostics=false,
    CmdStanDir=CMDSTAN_HOME);

  describe(chains)
  println()

  @time save("tmp/chains.jld",
    "range", chains.range,
    "a3d", chains.value,
    "names", chains.names,
    "chains", chains.chains)
  println()

  @time d = load(joinpath(ProjDir, "tmp", "chains.jld"))
  println()

  chn2 = MCMCChain.Chains(d["a3d"], names=d["names"])

  describe(chn2)

  @test 0.1 <  mean(chn2.value[:, 8, :] ) < 0.6

end # cd


text:
function setup_TriangleMesh(::Type{T}, ::Type{U}) where {T, U}
    name = "tri"
    vertices = [Point{2,T}(   0, 0),
                Point{2,T}(   1, 0),
                Point{2,T}(1//2, 1),
                Point{2,T}(3//2, 1)]
    polytopes = [Triangle{U}(1, 2, 3), Triangle{U}(2, 4, 3)]
    groups = Dict{String, BitSet}()
    groups["A"] = BitSet([1])
    groups["B"] = BitSet([2])
    return PolytopeVertexMesh(name, vertices, polytopes, groups)
end

function setup_QuadrilateralMesh(::Type{T}, ::Type{U}) where {T, U}
    name = "quad"
    vertices = [Point{2,T}(0, 0),
                Point{2,T}(1, 0),
                Point{2,T}(1, 1),
                Point{2,T}(0, 1),
                Point{2,T}(2, 0),
                Point{2,T}(2, 1)]
    polytopes = [Quadrilateral{U}(1, 2, 3, 4), Quadrilateral{U}(2, 5, 6, 3)]
    groups = Dict{String, BitSet}()
    groups["A"] = BitSet([1])
    groups["B"] = BitSet([2])
    return PolytopeVertexMesh(name, vertices, polytopes, groups)
end

function setup_MixedPolygonMesh(::Type{T}, ::Type{U}) where {T, U}
    name = "tri_quad"
    vertices = [Point{2,T}(0, 0),
                Point{2,T}(1, 0),
                Point{2,T}(1, 1),
                Point{2,T}(0, 1),
                Point{2,T}(2, 0)]
    polytopes = [Quadrilateral{U}(1, 2, 3, 4), Triangle{U}(2, 5, 3)]
    groups = Dict{String, BitSet}()
    groups["A"] = BitSet([1])
    groups["B"] = BitSet([2])
    return PolytopeVertexMesh(name, vertices, polytopes, groups)
end

function setup_QuadraticTriangleMesh(::Type{T}, ::Type{U}) where {T, U}
    name = "tri6"
    vertices = [Point{2,T}(   0,    0),
                Point{2,T}(   1,    0),
                Point{2,T}(   0,    1),
                Point{2,T}(1//2,    0),
                Point{2,T}(2//5, 2//5),
                Point{2,T}(   0, 1//2),
                Point{2,T}(   1,    1),
                Point{2,T}(   1, 1//2),
                Point{2,T}(1//2,    1)]
    polytopes = [QuadraticTriangle{U}(1, 2, 3, 4, 5, 6),
                 QuadraticTriangle{U}(2, 7, 3, 8, 9, 5)]
    groups = Dict{String, BitSet}()
    groups["A"] = BitSet([1])
    groups["B"] = BitSet([2])
    return PolytopeVertexMesh(name, vertices, polytopes, groups)
end

function setup_QuadraticQuadrilateralMesh(::Type{T}, ::Type{U}) where {T, U}
    name = "quad8"
    vertices = [Point{2,T}(    0,    0),
                Point{2,T}(    1,    0),
                Point{2,T}(    1,    1),
                Point{2,T}(    0,    1),
                Point{2,T}(    2,    0),
                Point{2,T}(    2,    1),
                Point{2,T}( 1//2,    0),
                Point{2,T}(7//10, 1//2),
                Point{2,T}( 1//2,    1),
                Point{2,T}(    0, 1//2),
                Point{2,T}( 3//2,    0),
                Point{2,T}(    2, 1//2),
                Point{2,T}( 3//2,    1),
               ]
    polytopes = [QuadraticQuadrilateral{U}(1, 2, 3, 4, 7, 8, 9, 10),
                 QuadraticQuadrilateral{U}(2, 5, 6, 3, 11, 12, 13, 8)]
    groups = Dict{String, BitSet}()
    groups["A"] = BitSet([1])
    groups["B"] = BitSet([2])
    return PolytopeVertexMesh(name, vertices, polytopes, groups)
end

function setup_MixedQuadraticPolygonMesh(::Type{T}, ::Type{U}) where {T, U}
    name = "tri6_quad8"
    vertices = [Point{2,T}(    0,    0),
                Point{2,T}(    1,    0),
                Point{2,T}(    1,    1),
                Point{2,T}(    0,    1),
                Point{2,T}(    2,    0),
                Point{2,T}( 1//2,    0),
                Point{2,T}(7//10, 1//2),
                Point{2,T}( 1//2,    1),
                Point{2,T}(    0, 1//2),
                Point{2,T}( 3//2,    0),
                Point{2,T}( 3//2, 1//2),
               ]
    polytopes = [QuadraticQuadrilateral{U}(1, 2, 3, 4, 6, 7, 8, 9),
                 QuadraticTriangle{U}(2, 5, 3, 10, 11, 7)]
    groups = Dict{String, BitSet}()
    groups["A"] = BitSet([1])
    groups["B"] = BitSet([2])
    return PolytopeVertexMesh(name, vertices, polytopes, groups)
end

ext: jl

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
"""
    ConstantAnisotropicDiffusivity{FT, KH, KV}

Parameters for constant anisotropic diffusivity models.
"""
struct ConstantAnisotropicDiffusivity{FT, KH, KV} <: TensorDiffusivity{FT}
    νh :: FT
    νv :: FT
    κh :: KH
    κv :: KV
    function ConstantAnisotropicDiffusivity{FT}(νh, νv, κh, κv) where FT
        return new{FT, typeof(κh), typeof(κv)}(νh, νv, convert_diffusivity(FT, κh), convert_diffusivity(FT, κv))
    end
end

"""
    ConstantAnisotropicDiffusivity(; νh, νv, κh, κv)

Returns parameters for a constant anisotropic diffusivity closure with constant horizontal
and vertical viscosities `νh`, `νv` and constant horizontal and vertical tracer
diffusivities `κh`, `κv`. `κh` and `κv` may be `NamedTuple`s with fields corresponding
to each tracer, or a single number to be a applied to all tracers.

By default, a viscosity of `ν = 1.05×10⁻⁶` m² s⁻¹ is used for both the horizontal
and vertical viscosity, and a diffusivity of `κ = 1.46×10⁻⁷` m² s⁻¹ is used
for the horizontal and vertical diffusivities applied to every tracer.
These values are the approximate viscosity and thermal diffusivity for seawater at 20°C
and 35 psu, according to Sharqawy et al., "Thermophysical properties of seawater: A review
of existing correlations and data" (2010).
"""
ConstantAnisotropicDiffusivity(FT=Float64; νh=ν₀, νv=ν₀, κh=κ₀, κv=κ₀) =
    ConstantAnisotropicDiffusivity{FT}(νh, νv, κh, κv)

function with_tracers(tracers, closure::ConstantAnisotropicDiffusivity{FT}) where FT
    κh = tracer_diffusivities(tracers, closure.κh)
    κv = tracer_diffusivities(tracers, closure.κv)
    return ConstantAnisotropicDiffusivity{FT}(closure.νh, closure.νv, κh, κv)
end

calculate_diffusivities!(K, arch, grid, closure::ConstantAnisotropicDiffusivity, args...) = nothing

@inline ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure::ConstantAnisotropicDiffusivity, U, args...) =
    ∂ⱼνᵢⱼ∂ᵢu(i, j, k, grid, closure.νh, closure.νh, closure.νv, U.u)

@inline ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, closure::ConstantAnisotropicDiffusivity, U, args...) =
    ∂ⱼνᵢⱼ∂ᵢv(i, j, k, grid, closure.νh, closure.νh, closure.νv, U.v)

@inline ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, closure::ConstantAnisotropicDiffusivity, U, args...) =
    ∂ⱼνᵢⱼ∂ᵢw(i, j, k, grid, closure.νh, closure.νh, closure.νv, U.w)

@inline function ∇_κ_∇c(i, j, k, grid, closure::ConstantAnisotropicDiffusivity,
                        c, ::Val{tracer_index}, args...) where tracer_index

    @inbounds κh = closure.κh[tracer_index]
    @inbounds κv = closure.κv[tracer_index]

    return ∂ⱼκᵢⱼ∂ᵢc(i, j, k, grid, κh, κh, κv, c)
end


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

## thread/task locking abstraction

"""
    AbstractLock

Abstract supertype describing types that
implement the synchronization primitives:
[`lock`](@ref), [`trylock`](@ref), [`unlock`](@ref), and [`islocked`](@ref).
"""
abstract type AbstractLock end
function lock end
function unlock end
function trylock end
function islocked end
unlockall(l::AbstractLock) = unlock(l) # internal function for implementing `wait`
relockall(l::AbstractLock, token::Nothing) = lock(l) # internal function for implementing `wait`
assert_havelock(l::AbstractLock) = assert_havelock(l, Threads.threadid())
assert_havelock(l::AbstractLock, tid::Integer) =
    (islocked(l) && tid == Threads.threadid()) ? nothing : error("concurrency violation detected")
assert_havelock(l::AbstractLock, tid::Task) =
    (islocked(l) && tid === current_task()) ? nothing : error("concurrency violation detected")
assert_havelock(l::AbstractLock, tid::Nothing) = error("concurrency violation detected")

"""
    AlwaysLockedST

This struct does not implement a real lock, but instead
pretends to be always locked on the original thread it was allocated on,
and simply ignores all other interactions.
It also does not synchronize tasks; for that use a real lock such as [`RecursiveLock`](@ref).
This can be used in the place of a real lock to, instead, simply and cheaply assert
that the operation is only occurring on a single cooperatively-scheduled thread.
It is thus functionally equivalent to allocating a real, recursive, task-unaware lock
immediately calling `lock` on it, and then never calling a matching `unlock`,
except that calling `lock` from another thread will throw a concurrency violation exception.
"""
struct AlwaysLockedST <: AbstractLock
    ownertid::Int16
    AlwaysLockedST() = new(Threads.threadid())
end
assert_havelock(l::AlwaysLockedST) = assert_havelock(l, l.ownertid)
lock(l::AlwaysLockedST) = assert_havelock(l)
unlock(l::AlwaysLockedST) = assert_havelock(l)
trylock(l::AlwaysLockedST) = l.ownertid == Threads.threadid()
islocked(::AlwaysLockedST) = true

## condition variables

"""
    GenericCondition

Abstract implementation of a condition object
for synchonizing tasks objects with a given lock.
"""
struct GenericCondition{L<:AbstractLock}
    waitq::Vector{Any}
    lock::L

    GenericCondition{L}() where {L<:AbstractLock} = new{L}([], L())
    GenericCondition{L}(l::L) where {L<:AbstractLock} = new{L}([], l)
    GenericCondition(l::AbstractLock) = new{typeof(l)}([], l)
end

assert_havelock(c::GenericCondition) = assert_havelock(c.lock)
lock(c::GenericCondition) = lock(c.lock)
unlock(c::GenericCondition) = unlock(c.lock)
trylock(c::GenericCondition) = trylock(c.lock)
islocked(c::GenericCondition) = islocked(c.lock)

"""
    wait([x])

Block the current task until some event occurs, depending on the type of the argument:

* [`Channel`](@ref): Wait for a value to be appended to the channel.
* [`Condition`](@ref): Wait for [`notify`](@ref) on a condition.
* `Process`: Wait for a process or process chain to exit. The `exitcode` field of a process
  can be used to determine success or failure.
* [`Task`](@ref): Wait for a `Task` to finish. If the task fails with an exception, the
  exception is propagated (re-thrown in the task that called `wait`).
* [`RawFD`](@ref): Wait for changes on a file descriptor (see the `FileWatching` package).

If no argument is passed, the task blocks for an undefined period. A task can only be
restarted by an explicit call to [`schedule`](@ref) or [`yieldto`](@ref).

Often `wait` is called within a `while` loop to ensure a waited-for condition is met before
proceeding.
"""
function wait(c::GenericCondition)
    ct = current_task()
    assert_havelock(c)
    push!(c.waitq, ct)
    token = unlockall(c.lock)

    try
        return wait()
    catch
        filter!(x->x!==ct, c.waitq)
        rethrow()
    finally
        relockall(c.lock, token)
    end
end

"""
    notify(condition, val=nothing; all=true, error=false)

Wake up tasks waiting for a condition, passing them `val`. If `all` is `true` (the default),
all waiting tasks are woken, otherwise only one is. If `error` is `true`, the passed value
is raised as an exception in the woken tasks.

Return the count of tasks woken up. Return 0 if no tasks are waiting on `condition`.
"""
notify(c::GenericCondition, @nospecialize(arg = nothing); all=true, error=false) = notify(c, arg, all, error)
function notify(c::GenericCondition, @nospecialize(arg), all, error)
    assert_havelock(c)
    cnt = 0
    if all
        cnt = length(c.waitq)
        for t in c.waitq
            schedule(t, arg, error=error)
        end
        empty!(c.waitq)
    elseif !isempty(c.waitq)
        cnt = 1
        t = popfirst!(c.waitq)
        schedule(t, arg, error=error)
    end
    return cnt
end

notify_error(c::GenericCondition, err) = notify(c, err, true, true)

n_waiters(c::GenericCondition) = length(c.waitq)

"""
    isempty(condition)

Return `true` if no tasks are waiting on the condition, `false` otherwise.
"""
isempty(c::GenericCondition) = isempty(c.waitq)

# default (Julia v1.0) is currently single-threaded
# (although it uses MT-safe versions, when possible)
"""
    Condition()

Create an edge-triggered event source that tasks can wait for. Tasks that call [`wait`](@ref) on a
`Condition` are suspended and queued. Tasks are woken up when [`notify`](@ref) is later called on
the `Condition`. Edge triggering means that only tasks waiting at the time [`notify`](@ref) is
called can be woken up. For level-triggered notifications, you must keep extra state to keep
track of whether a notification has happened. The [`Channel`](@ref) and [`Event`](@ref) types do
this, and can be used for level-triggered events.

This object is NOT thread-safe. See [`Threads.Condition`](@ref) for a thread-safe version.
"""
const Condition = GenericCondition{AlwaysLockedST}

## scheduler and work queue

global const Workqueue = Task[]

function enq_work(t::Task)
    t.state == :runnable || error("schedule: Task not runnable")
    ccall(:uv_stop, Cvoid, (Ptr{Cvoid},), eventloop())
    push!(Workqueue, t)
    t.state = :queued
    return t
end

schedule(t::Task) = enq_work(t)

"""
    schedule(t::Task, [val]; error=false)

Add a [`Task`](@ref) to the scheduler's queue. This causes the task to run constantly when the system
is otherwise idle, unless the task performs a blocking operation such as [`wait`](@ref).

If a second argument `val` is provided, it will be passed to the task (via the return value of
[`yieldto`](@ref)) when it runs again. If `error` is `true`, the value is raised as an exception in
the woken task.

# Examples
```jldoctest
julia> a5() = sum(i for i in 1:1000);

julia> b = Task(a5);

julia> istaskstarted(b)
false

julia> schedule(b);

julia> yield();

julia> istaskstarted(b)
true

julia> istaskdone(b)
true
```
"""
function schedule(t::Task, arg; error=false)
    # schedule a task to be (re)started with the given value or exception
    if error
        t.exception = arg
    else
        t.result = arg
    end
    return enq_work(t)
end

# fast version of `schedule(t, arg); wait()`
function schedule_and_wait(t::Task, arg=nothing)
    t.state == :runnable || error("schedule: Task not runnable")
    if isempty(Workqueue)
        return yieldto(t, arg)
    else
        t.result = arg
        push!(Workqueue, t)
        t.state = :queued
    end
    return wait()
end

"""
    yield()

Switch to the scheduler to allow another scheduled task to run. A task that calls this
function is still runnable, and will be restarted immediately if there are no other runnable
tasks.
"""
yield() = (enq_work(current_task()); wait())

"""
    yield(t::Task, arg = nothing)

A fast, unfair-scheduling version of `schedule(t, arg); yield()` which
immediately yields to `t` before calling the scheduler.
"""
function yield(t::Task, @nospecialize x = nothing)
    t.state == :runnable || error("schedule: Task not runnable")
    t.result = x
    enq_work(current_task())
    return try_yieldto(ensure_rescheduled, Ref(t))
end

"""
    yieldto(t::Task, arg = nothing)

Switch to the given task. The first time a task is switched to, the task's function is
called with no arguments. On subsequent switches, `arg` is returned from the task's last
call to `yieldto`. This is a low-level call that only switches tasks, not considering states
or scheduling in any way. Its use is discouraged.
"""
function yieldto(t::Task, @nospecialize x = nothing)
    t.result = x
    return try_yieldto(identity, Ref(t))
end

function try_yieldto(undo, reftask::Ref{Task})
    try
        ccall(:jl_switchto, Cvoid, (Any,), reftask)
    catch
        undo(reftask[])
        rethrow()
    end
    ct = current_task()
    exc = ct.exception
    if exc !== nothing
        ct.exception = nothing
        throw(exc)
    end
    result = ct.result
    ct.result = nothing
    return result
end

# yield to a task, throwing an exception in it
function throwto(t::Task, @nospecialize exc)
    t.exception = exc
    return yieldto(t)
end

function ensure_rescheduled(othertask::Task)
    ct = current_task()
    if ct !== othertask && othertask.state == :runnable
        # we failed to yield to othertask
        # return it to the head of the queue to be scheduled later
        pushfirst!(Workqueue, othertask)
        othertask.state = :queued
    end
    if ct.state == :queued
        # if the current task was queued,
        # also need to return it to the runnable state
        # before throwing an error
        i = findfirst(t->t===ct, Workqueue)
        i === nothing || deleteat!(Workqueue, i)
        ct.state = :runnable
    end
    nothing
end

@noinline function poptask()
    t = popfirst!(Workqueue)
    if t.state != :queued
        # assume this somehow got queued twice,
        # probably broken now, but try discarding this switch and keep going
        # can't throw here, because it's probably not the fault of the caller to wait
        # and don't want to use print() here, because that may try to incur a task switch
        ccall(:jl_safe_printf, Cvoid, (Ptr{UInt8}, Int32...),
            "\nWARNING: Workqueue inconsistency detected: popfirst!(Workqueue).state != :queued\n")
        return
    end
    t.state = :runnable
    return Ref(t)
end

function wait()
    while true
        if isempty(Workqueue)
            c = process_events(true)
            if c == 0 && eventloop() != C_NULL && isempty(Workqueue)
                # if there are no active handles and no runnable tasks, just
                # wait for signals.
                pause()
            end
        else
            reftask = poptask()
            if reftask !== nothing
                result = try_yieldto(ensure_rescheduled, reftask)
                process_events(false)
                # return when we come out of the queue
                return result
            end
        end
    end
    # unreachable
end

if Sys.iswindows()
    pause() = ccall(:Sleep, stdcall, Cvoid, (UInt32,), 0xffffffff)
else
    pause() = ccall(:pause, Cvoid, ())
end

## async event notifications

"""
    AsyncCondition()

Create a async condition that wakes up tasks waiting for it
(by calling [`wait`](@ref) on the object)
when notified from C by a call to `uv_async_send`.
Waiting tasks are woken with an error when the object is closed (by [`close`](@ref).
Use [`isopen`](@ref) to check whether it is still active.
"""
mutable struct AsyncCondition
    handle::Ptr{Cvoid}
    cond::Condition
    isopen::Bool

    function AsyncCondition()
        this = new(Libc.malloc(_sizeof_uv_async), Condition(), true)
        associate_julia_struct(this.handle, this)
        finalizer(uvfinalize, this)
        err = ccall(:uv_async_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
            eventloop(), this, uv_jl_asynccb::Ptr{Cvoid})
        if err != 0
            #TODO: this codepath is currently not tested
            Libc.free(this.handle)
            this.handle = C_NULL
            throw(_UVError("uv_async_init", err))
        end
        return this
    end
end

"""
    AsyncCondition(callback::Function)

Create a async condition that calls the given `callback` function. The `callback` is passed one argument,
the async condition object itself.
"""
function AsyncCondition(cb::Function)
    async = AsyncCondition()
    waiter = Task(function()
        while isopen(async)
            success = try
                wait(async)
                true
            catch exc # ignore possible exception on close()
                isa(exc, EOFError) || rethrow()
            end
            success && cb(async)
        end
    end)
    # must start the task right away so that it can wait for the AsyncCondition before
    # we re-enter the event loop. this avoids a race condition. see issue #12719
    yield(waiter)
    return async
end

## timer-based notifications

"""
    Timer(delay; interval = 0)

Create a timer that wakes up tasks waiting for it (by calling [`wait`](@ref) on the timer object).

Waiting tasks are woken after an initial delay of `delay` seconds, and then repeating with the given
`interval` in seconds. If `interval` is equal to `0`, the timer is only triggered once. When
the timer is closed (by [`close`](@ref) waiting tasks are woken with an error. Use [`isopen`](@ref)
to check whether a timer is still active.
"""
mutable struct Timer
    handle::Ptr{Cvoid}
    cond::Condition
    isopen::Bool

    function Timer(timeout::Real; interval::Real = 0.0)
        timeout ≥ 0 || throw(ArgumentError("timer cannot have negative timeout of $timeout seconds"))
        interval ≥ 0 || throw(ArgumentError("timer cannot have negative repeat interval of $interval seconds"))

        this = new(Libc.malloc(_sizeof_uv_timer), Condition(), true)
        err = ccall(:uv_timer_init, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), eventloop(), this)
        if err != 0
            #TODO: this codepath is currently not tested
            Libc.free(this.handle)
            this.handle = C_NULL
            throw(_UVError("uv_timer_init", err))
        end

        associate_julia_struct(this.handle, this)
        finalizer(uvfinalize, this)

        ccall(:uv_update_time, Cvoid, (Ptr{Cvoid},), eventloop())
        ccall(:uv_timer_start,  Cint,  (Ptr{Cvoid}, Ptr{Cvoid}, UInt64, UInt64),
              this, uv_jl_timercb::Ptr{Cvoid},
              UInt64(round(timeout * 1000)) + 1, UInt64(round(interval * 1000)))
        return this
    end
end

unsafe_convert(::Type{Ptr{Cvoid}}, t::Timer) = t.handle
unsafe_convert(::Type{Ptr{Cvoid}}, async::AsyncCondition) = async.handle

function wait(t::Union{Timer, AsyncCondition})
    isopen(t) || throw(EOFError())
    stream_wait(t, t.cond)
end

isopen(t::Union{Timer, AsyncCondition}) = t.isopen

function close(t::Union{Timer, AsyncCondition})
    if t.handle != C_NULL && isopen(t)
        t.isopen = false
        isa(t, Timer) && ccall(:uv_timer_stop, Cint, (Ptr{Cvoid},), t)
        ccall(:jl_close_uv, Cvoid, (Ptr{Cvoid},), t)
    end
    nothing
end

function uvfinalize(t::Union{Timer, AsyncCondition})
    if t.handle != C_NULL
        disassociate_julia_struct(t.handle) # not going to call the usual close hooks
        close(t)
        t.handle = C_NULL
    end
    t.isopen = false
    nothing
end

function _uv_hook_close(t::Union{Timer, AsyncCondition})
    uvfinalize(t)
    notify_error(t.cond, EOFError())
    nothing
end

function uv_asynccb(handle::Ptr{Cvoid})
    async = @handle_as handle AsyncCondition
    notify(async.cond)
    nothing
end

function uv_timercb(handle::Ptr{Cvoid})
    t = @handle_as handle Timer
    if ccall(:uv_timer_get_repeat, UInt64, (Ptr{Cvoid},), t) == 0
        # timer is stopped now
        close(t)
    end
    notify(t.cond)
    nothing
end

"""
    sleep(seconds)

Block the current task for a specified number of seconds. The minimum sleep time is 1
millisecond or input of `0.001`.
"""
function sleep(sec::Real)
    sec ≥ 0 || throw(ArgumentError("cannot sleep for $sec seconds"))
    wait(Timer(sec))
    nothing
end

# timer with repeated callback
"""
    Timer(callback::Function, delay; interval = 0)

Create a timer that wakes up tasks waiting for it (by calling [`wait`](@ref) on the timer object) and
calls the function `callback`.

Waiting tasks are woken and the function `callback` is called after an initial delay of `delay` seconds,
and then repeating with the given `interval` in seconds. If `interval` is equal to `0`, the timer
is only triggered once. The function `callback` is called with a single argument, the timer itself.
When the timer is closed (by [`close`](@ref) waiting tasks are woken with an error. Use [`isopen`](@ref)
to check whether a timer is still active.

# Examples

Here the first number is printed after a delay of two seconds, then the following numbers are printed quickly.

```julia-repl
julia> begin
           i = 0
           cb(timer) = (global i += 1; println(i))
           t = Timer(cb, 2, interval = 0.2)
           wait(t)
           sleep(0.5)
           close(t)
       end
1
2
3
```
"""
function Timer(cb::Function, timeout::Real; interval::Real = 0.0)
    t = Timer(timeout, interval = interval)
    waiter = Task(function()
        while isopen(t)
            success = try
                wait(t)
                true
            catch exc # ignore possible exception on close()
                isa(exc, EOFError) || rethrow()
                false
            end
            success && cb(t)
        end
    end)
    # must start the task right away so that it can wait for the Timer before
    # we re-enter the event loop. this avoids a race condition. see issue #12719
    yield(waiter)
    return t
end

max_stars_repo_name: EricForgy/julia

max_issues_repo_name: EricForgy/julia

max_forks_repo_name: EricForgy/julia


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:5] <= 1)

                     @variable(m, -36.453305806680525 <= q <= 37.480198521751326)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.5060529352779994 + 0.34605798298901425*1/(1 + exp(-(-0.49627455608242954 + -0.46375151869975495*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + 0.4918746175115887*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + 0.38157591143072933*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.5439511060803288*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + 0.8504066842806357*1/(1 + exp(-(-0.2321501407718678 + 0.5600091754803458*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.0034071169835434567*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + -0.6765465160367254*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + 0.8220752699067027*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + 0.462211279142271*1/(1 + exp(-(-0.9807289687619161 + -0.5968595287514287*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.3897253830589995*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + 0.5397807246416386*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.24558591797668194*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + -0.0428411186394575*1/(1 + exp(-(0.7899732054652109 + -0.5363355140855166*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.6727460735415476*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + -0.4129142282119118*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.4150508411017726*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))))))))) + 1/(1 + exp(-(0.02992678336456356 + -0.8781954986928167*1/(1 + exp(-(-0.49627455608242954 + -0.46375151869975495*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + 0.4918746175115887*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + 0.38157591143072933*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.5439511060803288*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + -0.84297522225665*1/(1 + exp(-(-0.2321501407718678 + 0.5600091754803458*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.0034071169835434567*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + -0.6765465160367254*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + 0.8220752699067027*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + 0.44937020244474857*1/(1 + exp(-(-0.9807289687619161 + -0.5968595287514287*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.3897253830589995*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + 0.5397807246416386*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.24558591797668194*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + -0.012040065223148932*1/(1 + exp(-(0.7899732054652109 + -0.5363355140855166*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.6727460735415476*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + -0.4129142282119118*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.4150508411017726*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))))))))) + 1/(1 + exp(-(-0.25761034404046734 + -0.19911102915318102*1/(1 + exp(-(-0.49627455608242954 + -0.46375151869975495*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + 0.4918746175115887*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + 0.38157591143072933*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.5439511060803288*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + 0.7457041629864394*1/(1 + exp(-(-0.2321501407718678 + 0.5600091754803458*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.0034071169835434567*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + -0.6765465160367254*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + 0.8220752699067027*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + 0.7638178262546895*1/(1 + exp(-(-0.9807289687619161 + -0.5968595287514287*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.3897253830589995*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + 0.5397807246416386*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.24558591797668194*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + 0.07889233901952375*1/(1 + exp(-(0.7899732054652109 + -0.5363355140855166*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.6727460735415476*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + -0.4129142282119118*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.4150508411017726*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))))))))) + 1/(1 + exp(-(0.1153458350428731 + 0.9297971624965538*1/(1 + exp(-(-0.49627455608242954 + -0.46375151869975495*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + 0.4918746175115887*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + 0.38157591143072933*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.5439511060803288*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + -0.7017891371439648*1/(1 + exp(-(-0.2321501407718678 + 0.5600091754803458*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.0034071169835434567*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + -0.6765465160367254*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + 0.8220752699067027*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + -0.8093911002605614*1/(1 + exp(-(-0.9807289687619161 + -0.5968595287514287*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.3897253830589995*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + 0.5397807246416386*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.24558591797668194*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5]))))))) + 0.43763294749909276*1/(1 + exp(-(0.7899732054652109 + -0.5363355140855166*1/(1 + exp(-(0.8994638843665332 + -0.21460729883264662*$(x[1]) + -0.6173390778114793*$(x[2]) + 0.5769355957045894*$(x[3]) + -0.29054193627721014*$(x[4]) + 0.9844022485972483*$(x[5])))) + -0.6727460735415476*1/(1 + exp(-(0.2763051493368929 + 0.8533339132011375*$(x[1]) + -0.36082426348518304*$(x[2]) + -0.694925986364535*$(x[3]) + -0.6948010174882064*$(x[4]) + 0.21743382793841226*$(x[5])))) + -0.4129142282119118*1/(1 + exp(-(0.044944779340077456 + 0.03687124531051067*$(x[1]) + -0.5175326375722826*$(x[2]) + -0.9986359847691255*$(x[3]) + -0.33739875692716526*$(x[4]) + 0.025638388186982386*$(x[5])))) + -0.4150508411017726*1/(1 + exp(-(0.8609963866663906 + -0.4203074235446964*$(x[1]) + 0.91553501702727*$(x[2]) + 0.8830101091308449*$(x[3]) + 0.07593641537869811*$(x[4]) + -0.2012556317025167*$(x[5])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function f()
    print('>')
    readline()
    f()
end
f()

max_stars_repo_name: fuzzynoise/ono

max_issues_repo_name: fuzzynoise/ono

max_forks_repo_name: fuzzynoise/ono


text:
# wraps variables found in UNITS.INC for use in Julia

struct UNITS_s
    UNITL::Array{Float64,1}
    UNITM::Array{Float64,1}
    UNITT::Array{Float64,1}
    UNITF::Array{Float64,1}
    UNITP::Array{Float64,1}
    UNITH::Array{Float64,1}
    UNITE::Array{Float64,1}
    UNITV::Array{Float64,1}
    UNITA::Array{Float64,1}
    UNITS::Array{Float64,1}
    UNITI::Array{Float64,1}
    UNITB::Array{Float64,1}
    UNITD::Array{Float64,1}
    UNITR::Array{Float64,1}
    UNITQ::Array{Float64,1}
    UNITG::Array{Float64,1}
    UNITC::Array{Float64,1}
    UNCHL::Array{Cchar,1}
    UNCHM::Array{Cchar,1}
    UNCHT::Array{Cchar,1}
    UNCHF::Array{Cchar,1}
    UNCHP::Array{Cchar,1}
    UNCHH::Array{Cchar,1}
    UNCHE::Array{Cchar,1}
    UNCHV::Array{Cchar,1}
    UNCHA::Array{Cchar,1}
    UNCHS::Array{Cchar,1}
    UNCHI::Array{Cchar,1}
    UNCHB::Array{Cchar,1}
    UNCHD::Array{Cchar,1}
    UNCHR::Array{Cchar,1}
    UNCHQ::Array{Cchar,1}
    UNCHG::Array{Cchar,1}
    UNCHC::Array{Cchar,1}
    NUL::Array{Int32,1}
    NUM::Array{Int32,1}
    NUT::Array{Int32,1}
    NUF::Array{Int32,1}
    NUP::Array{Int32,1}
    NUH::Array{Int32,1}
    NUE::Array{Int32,1}
    NUV::Array{Int32,1}
    NUA::Array{Int32,1}
    NUS::Array{Int32,1}
    NUI::Array{Int32,1}
    NUB::Array{Int32,1}
    NUD::Array{Int32,1}
    NUR::Array{Int32,1}
    NUQ::Array{Int32,1}
    NUG::Array{Int32,1}
    NUC::Array{Int32,1}
    function UNITS_s()
        as_unr = cglobal((:as_unr_, libaswing), Float64); as_unr_idx = as_unr
        UNITL  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITM  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITT  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITF  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITP  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITH  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITE  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITV  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITA  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITS  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITI  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITB  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITD  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITR  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITQ  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITG  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)
        UNITC  = unsafe_wrap(Array, as_unr_idx, 1); as_unr_idx += sizeof(Float64)

        as_unc = cglobal((:as_unc_, libaswing), Cchar); as_unc_idx = as_unc
        UNCHL  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHM  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHT  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHF  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHP  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHH  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHE  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHV  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHA  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHS  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHI  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHB  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHD  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHR  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHQ  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHG  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32
        UNCHC  = unsafe_wrap(Array, as_unc_idx, 32); as_unc_idx += sizeof(Cchar)*32

        as_uni = cglobal((:as_uni_, libaswing), Int32); as_uni_idx = as_uni
        NUL  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUM  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUT  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUF  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUP  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUH  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUE  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUV  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUA  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUS  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUI  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUB  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUD  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUR  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUQ  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUG  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)
        NUC  = unsafe_wrap(Array, as_uni_idx, 1); as_uni_idx += sizeof(Int32)

        return new(UNITL, UNITM, UNITT, UNITF, UNITP, UNITH, UNITE,
            UNITV, UNITA, UNITS, UNITI, UNITB, UNITD, UNITR,
            UNITQ, UNITG, UNITC,
            UNCHL, UNCHM, UNCHT, UNCHF, UNCHP, UNCHH, UNCHE,
            UNCHV, UNCHA, UNCHS, UNCHI, UNCHB, UNCHD, UNCHR,
            UNCHQ, UNCHG, UNCHC,
            NUL,   NUM,   NUT,   NUF,   NUP,   NUH,   NUE,
            NUV,   NUA,   NUS,   NUI,   NUB,   NUD,   NUR,
            NUQ,   NUG,   NUC)
    end
end


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 960703e4-98b3-11eb-327d-e3640b19b7e8
begin
	import Pkg; Pkg.activate(@__DIR__)
	using Plots, PlutoUI
end

# ╔═╡ 688ba03f-a052-421d-9c0c-91b33a27d02c
begin
	#@gif for p ∈ plots
#		plot(p)
#	end
end;

# ╔═╡ 664f2bc3-bd52-4240-9443-b2fcd7500f04
begin
	T₁ = 200
	T₂ = 80
	ω₀ = 0.25
	tmax = 1000
	N = 10000
	ts = range(0,tmax,length=N)
	M0x = 0.5
	M0y = 0.5
	M0z = 0
	M = 1
	fMx(t) = (M0x * cos(ω₀ * t) - M0y * sin(ω₀ * t)) * exp(-t/T₂)
	fMy(t) = (M0y * cos(ω₀ * t) + M0y * sin(ω₀ * t)) * exp(-t/T₂)
	fMz(t) = M0z * exp(-t/T₁) + M * (1 - exp(-t/T₁))
	Mxs = fMx.(ts)
	Mys = fMy.(ts)
	Mzs = fMz.(ts)
end;

# ╔═╡ bfd1b44e-b63e-4644-bb91-f2e7d3a6d2eb
plots = map(1:100:N) do i
	@views plot(Mxs[1:i], Mys[1:i], Mzs[1:i], xlims=(-1,1), ylims=(-1,1),zlims=(0,1), leg=:none, aspect_ratio=:equal)
end;

# ╔═╡ 8daf0654-5b8e-4e53-8655-6c52173cc0d5
@bind i Slider(1:length(plots), show_value=true)

# ╔═╡ 6cc56a8a-f937-4735-9f5f-4162aec6c97c
begin
	plots[i]
end

# ╔═╡ Cell order:
# ╠═688ba03f-a052-421d-9c0c-91b33a27d02c
# ╟─8daf0654-5b8e-4e53-8655-6c52173cc0d5
# ╠═6cc56a8a-f937-4735-9f5f-4162aec6c97c
# ╠═bfd1b44e-b63e-4644-bb91-f2e7d3a6d2eb
# ╠═664f2bc3-bd52-4240-9443-b2fcd7500f04
# ╠═960703e4-98b3-11eb-327d-e3640b19b7e8


text:
"A type for the ensemble average scattering coefficients.
Here they are discretised in terms of the depth x of the halfspace"
mutable struct AverageWave{T<:AbstractFloat}
    hankel_order::Int # largest hankel order
    x::Vector{T} # spatial mesh
    amplitudes::Array{Complex{T}} # a matrix of the scattering amplitudes, size(A_mat) = (length(x), 2hankel_order +1)
    # Enforce that the dimensions are correct
    function AverageWave{T}(hankel_order::Int, x::Vector{T}, amplitudes::Array{Complex{T}}) where T <: AbstractFloat
        if (length(x), 2*hankel_order+1) != size(amplitudes)[1:2]
            error("The amplitudes of AverageWave does not satisfy size(amplitudes)[1:2] == (length(X), 2*hankel_order+1)")
        end
        new(hankel_order,x,amplitudes)
    end
end

AverageWave(M::Int, x::AbstractVector{T}, as::AbstractArray{Complex{T}}) where T<:AbstractFloat = AverageWave{T}(M,collect(x),collect(as))

function AverageWave(x::AbstractVector{T}, A_mat::Array{Complex{T}}) where T<:Number
    AverageWave(Int((size(A_mat,2)-1)/2), collect(x), A_mat)
end

"Approximates the error in AverageWave.amplitudes due to the mesh AverageWave.x."
function average_error(avg_w::AverageWave)
    ddf = circshift(avg_w.amplitudes, (2,0,0)) - 2*circshift(avg_w.amplitudes, (1,0,0)) + avg_w.amplitudes
    h = (avg_w.x[2] - avg_w.x[1])
    ddf = ddf[3:end,:,:]/(h^2)
    max_ddf = maximum(abs.(ddf))

    # trapezoidal error
    return max_ddf*h^3/12
end

"Calculates an AverageWave from one EffectiveWave"
function AverageWave(xs::AbstractVector{T}, wave_eff::EffectiveWave{T}) where T<:Number

    amps = wave_eff.amplitudes
    ho = wave_eff.hankel_order
    θ_eff = wave_eff.θ_eff

    S = size(amps,2)

    average_amps = [
        im^T(m)*exp(-im*m*θ_eff)*amps[m+ho+1,s]*exp(im*wave_eff.k_eff*cos(θ_eff)*x)
    for x in xs, m=-ho:ho, s=1:S]

    return AverageWave(ho,xs,average_amps)
end

"Calculates an AverageWave from a vector of EffectiveWave"
function AverageWave(xs::AbstractVector{T}, wave_effs::Vector{EffectiveWave{T}}) where T<:Number

    ho = wave_effs[1].hankel_order
    avg_wave_effs = [AverageWave(xs, wave) for wave in wave_effs]
    amps = sum(avg_wave_effs[i].amplitudes[:,:,:] for i in eachindex(avg_wave_effs))

    return AverageWave(ho, xs, amps)
end

"Numerically solved the integral equation governing the average wave. Optionally can use wave_eff to approximate the wave away from the boundary."
function AverageWave(ω::T, medium::Medium{T}, specie::Specie{T};
        radius_multiplier::T = 1.005,
        x::AbstractVector{T} = [zero(T)],
        tol::T = T(1e-4),
        wave_effs::Vector{EffectiveWave{T}} = [zero(EffectiveWave{T})],
        max_size::Int = 1000,
        kws...
    ) where T<:Number

    k = real(ω/medium.c)

    if x == [zero(T)]
        if maximum(abs(w.k_eff) for w in wave_effs) == zero(T)
            wave_effs = effective_waves(real(ω/medium.c), medium, [specie];
                radius_multiplier=radius_multiplier, tol=tol, mesh_points=2, kws...)
        end
        # estimate a large coarse non-dimensional mesh based on the lowest attenuating effective wave
        a12 = T(2)*radius_multiplier*specie.r
        x = x_mesh(wave_effs[1]; tol = tol,  a12 = a12, max_size=max_size)
    end

    X = x.*k
    (MM_quad,b_mat) = average_wave_system(ω, X, medium, specie; tol = tol, kws...);

    M = Int( (size(b_mat,2) - 1)/2 )
    J = length(collect(X)) - 1

    len = (J + 1) * (2M + 1)
    MM_mat = reshape(MM_quad, (len, len));
    b = reshape(b_mat, (len));

    As = MM_mat\b
    As_mat = reshape(As, (J+1, 2M+1, 1))

    return AverageWave(M, collect(X)./k, As_mat)
end

"note that this uses the non-dimensional X = k*depth"
function average_wave_system(ω::T, X::AbstractVector{T}, medium::Medium{T}, specie::Specie{T};
        θin::Float64 = 0.0, tol::T = 1e-6,
        radius_multiplier::T = 1.005,
        scheme::Symbol = :trapezoidal,
        hankel_order::Int = maximum_hankel_order(ω, medium, [specie]; tol = tol),
        t_vecs::Vector{Vector{Complex{T}}} = t_vectors(ω, medium, [specie]; hankel_order = hankel_order),
        kws...
    ) where T<:AbstractFloat

    k = real(ω/medium.c)
    a12k = radius_multiplier*T(2)*real(k*specie.r);
    M = hankel_order;

    J = length(X) - 1
    h = X[2] - X[1]

    PQ_quad = intergrand_kernel(X, a12k; M = M, θin = θin, scheme=scheme);

    MM_quad = [
        (specie.num_density/(k^2))*t_vecs[1][m+M+1]*PQ_quad[l,m+M+1,j,n+M+1] - ( (m==n && j==l) ? 1.0+0.0im : 0.0+0.0im)
    for  l=1:(J+1), m=-M:M, j=1:(J+1), n=-M:M];

    b_mat = [
        -t_vecs[1][m+M+1]*exp(im*X[l]*cos(θin))*exp(im*m*(pi/2.0 - θin))
    for l = 1:(J+1), m = -M:M]

    return (MM_quad,b_mat)
end

"Returns x the mesh used to discretise the integral equations."
function x_mesh(wave_eff_long::EffectiveWave{T}, wave_eff_short::EffectiveWave{T} = wave_eff_long;
        tol::T = T(1e-5),  a12::T = T(Inf),
        max_size::Int = 1000,
        min_size::Int = 5,
        max_x::T = (-log(tol))/abs(cos(wave_eff_long.θ_eff)*imag(wave_eff_long.k_eff))
) where T<:AbstractFloat

    #= The default min_X result in:
        abs(exp(im*min_X*cos(θ_effs[end])*k_effs[end]/k)) < tol
    =#
    # estimate a reasonable derivative based on more rapidly varying wave_eff_short.
    df = abs(wave_eff_short.k_eff * cos(wave_eff_short.θ_eff))

    # Based on Simpson's rule
        # dX  = (tol*90 / (df^4))^(1/5)
    # Based on trapezoidal integration
        dx  = (tol * T(30) / (df^2))^(1/3)

    # prioritise a mesh fine enough to give accurate discrete integrals
    if dx > a12 dx = a12 end

    if max_x/dx + 1 > max_size
        if dx == a12
            max_x = (max_size - 1)*dx
            @warn("The mesh max_size = $max_size which was too small for tol = $tol. Will shrink meshed region.")
        else # otherwise priortise shrink max_x and increase dx proportionally
            a = sqrt(max_x/((max_size - 1)*dx))
            dx = min(dx*a,a12)
            max_x = (max_size - 1)*dx
            @warn("The mesh max_size = $max_size which was too small for tol = $tol. Will make a smaller and coarser mesh.")
        end
    elseif max_x/dx + 1 < min_size
        max_x = (min_size - 1)*dx
    end

    # if whole correction length a12k was given, then make dX/a12k = integer
    if a12 != T(Inf) && dx < a12
        n = ceil(a12 / dx)
        dx = a12/n
    end

    return zero(T):dx:max_x
end


text:
# This code tests the autodiffed Drift/Kepler step

include("../src/kepler_drift_gamma.jl")

@testset "kepler_drift_gamma" begin

function test_kepler_drift(dlnq::BigFloat,drift_first::Bool)
# Call as: save,jac_num,jacobian=test_kepler_drift(1e-4,true)
# This routine runs a test of the kep_elliptic_jacobian function in kepler_solver.jl
# Define the central force constant in terms of AU and days:
k = (2.0*pi/365.25)^2
# Initial position at 1 AU:
x0 = [0.01,1.0,0.01]
# Circular velocity:
r0 = sqrt(x0[1]*x0[1]+x0[2]*x0[2]+x0[3]*x0[3])
vcirc = sqrt(k/r0)
# Define initial velocity at apastron:
v0 = [.9*vcirc,0.01*vcirc,0.01*vcirc]  # The eccentricity is about ~2(1-v0/vcirc).
#h = 18.0 # 18-day timesteps
h = 0.000005 # 18-day timesteps
#h = 0.05 # 18-day timesteps
hbig = big(h)

s0::Float64 = 0.0
x = zeros(Float64,3)
v = zeros(Float64,3)
t = 0.0
nsteps = 1
xsave=zeros(Float64,12,nsteps)  # First is time (1,:); next three are position (2-4,:); next three velocity (5-7,:); then r (8,:), drdt (9,:); then beta (10,:); finally s & ds (11:12,:)
# Create a variable to store the state at the end of a step:
state=zeros(Float64,12)
state_diffbig=zeros(BigFloat,12)

xsave[1,1]=0.0
for jj=1:3
  xsave[jj+1,1]=x0[jj]
  xsave[jj+4,1]=v0[jj]
end
# Save beta:
xsave[8,1]=r0
beta0 = 2.0*k/r0-dot(v0,v0)
xsave[10,1]=beta0
# First, compute it with the old version:
jacobian_old =zeros(Float64,7,7)
d_old = Derivatives(Float64)
iter = kep_drift_ell_hyp!(x0,v0,k,h,s0,state,jacobian_old,drift_first,d_old)
#println("s old: ",state[11])

# Check that old and new are giving the same answer:
delxv = jac_delxv_gamma!(x0,v0,k,h,drift_first;debug=true)
println("state: ",state)
println("delxv: ",delxv)
println("Old-new: ",state[2:7]-delxv[1:6])

# Now compute algebraic Jacobian:
delxv,jac_alg = jac_delxv_gamma!(x0,v0,k,h,drift_first;grad=true,auto=false,debug=true)

# Next, compute autodiff Jacobian:
delxv,jacobian = jac_delxv_gamma!(x0,v0,k,h,drift_first;grad=true,auto=true,debug=true)

# Now, do finite differences at higher precision:
kbig = big(k)
x0big = big.(x0); v0big = big.(v0)
# Now compute big-float precision autodiff Jacobian:
delxv_big,jacobian_big = jac_delxv_gamma!(x0big,v0big,kbig,hbig,drift_first;grad=true,auto=true,debug=true)
#println("Auto diff Jacobian: ",jacobian)
jac_frac = jac_alg[1:6,:]./convert(Array{Float64,2},jacobian_big[1:6,:]).-1.0
println("Fractional Jacobian difference: ",maxabs(jac_frac[.~isnan.(jac_frac)]))

# Now compute finite-difference Jacobian:
delxv,jac_num = jac_delxv_gamma!(x0big,v0big,kbig,hbig,drift_first;grad=true,auto=false,dlnq=dlnq,debug=true)
name = ["x01","x02","x03","v01","v02","v03","gamma","r","fm1","fdot","gmhgdot","gdotm1"]
for i=1:12
  println(i," ",name[i])
  println("Fini diff Jacobian: ",convert(Array{Float64,1},jac_num[i,:]))
  println("Algebraic Jacobian: ",convert(Array{Float64,1},jac_alg[i,:]))
  println("Auto diff Jacobian: ",convert(Array{Float64,1},jacobian_big[i,:]))
  println("Frac diff Jacobian: ",jac_alg[i,:]./convert(Array{Float64,1},jacobian_big[i,:]).-1.0)
end
println("gamma alg  gradient: ",jac_alg[7,:])
println("gamma auto gradient: ",convert(Array{Float64,1},jacobian_big[7,:]))
println("gamma diff gradient: ",convert(Array{Float64,1},jac_num[7,:]))
println("r alg  gradient: ",jac_alg[8,:])
println("r auto gradient: ",convert(Array{Float64,1},jacobian_big[8,:]))
println("r diff gradient: ",convert(Array{Float64,1},jac_num[8,:]))
println("fm1 gradient  alg: ",jac_alg[9,:])
println("fm1 gradient auto: ",convert(Array{Float64,1},jacobian_big[9,:]))
println("fm1 gradient diff: ",convert(Array{Float64,1},jac_num[9,:]))
println("fdot gradient  alg: ",jac_alg[10,:])
println("fdot gradient auto: ",convert(Array{Float64,1},jacobian_big[10,:]))
println("fdot gradient diff: ",convert(Array{Float64,1},jac_num[10,:]))
println("g-h gdot gradient   alg: ",jac_alg[11,:])
println("g-h gdot gradient  auto: ",convert(Array{Float64,1},jacobian_big[11,:]))
println("g-h gdot gradient  diff: ",convert(Array{Float64,1},jac_num[11,:]))
println("gdot-1 gradient  alg: ",jac_alg[12,:])
println("gdot-1 gradient auto: ",convert(Array{Float64,1},jacobian_big[12,:]))
println("gdot-1 gradient diff: ",convert(Array{Float64,1},jac_num[12,:]))

println("Dim of jacobian_big: ",size(jacobian_big)," dim of jac_num: ",size(jac_num))
println("Maximum jac_big-jac_num: ",maxabs(convert(Array{Float64,2},jacobian_big-jac_num)))
println("Maximum jac_big-jac_alg: ",maxabs(convert(Array{Float64,2},jacobian_big-jac_alg)))
return xsave,jac_num,jac_alg
end

# First try:
xsave,jac_num1,jacobian=test_kepler_drift(big(1e-20),true)
emax = 0.0; imax = 0; jmax = 0
#for i=1:6, j=1:8
for i=1:12, j=1:8
  if jacobian[i,j] != 0.0
    diff = abs(convert(Float64,jac_num1[i,j])/jacobian[i,j]-1.0)
    if  diff > emax
      emax = diff; imax = i; jmax = j
#      println("i ",i," j ",j," jac_num1: ",convert(Float64,jac_num1[i,j])," jacobian: ",jacobian[i,j])
    end
  end
end
println("Maximum fractional error: ",emax," ",imax," ",jmax)
println("Maximum jacobian-jac_num: ",maxabs(jacobian-convert(Array{Float64,2},jac_num1)))

# Second try:
xsave,jac_num1,jacobian=test_kepler_drift(big(1e-20),false)

emax = 0.0; imax = 0; jmax = 0
for i=1:6, j=1:8
  if jacobian[i,j] != 0.0
    diff = abs(convert(Float64,jac_num1[i,j])/jacobian[i,j]-1.0)
    if  diff > emax
      emax = diff; imax = i; jmax = j
#      println("i ",i," j ",j," jac_num1: ",convert(Float64,jac_num1[i,j])," jacobian: ",jacobian[i,j])
    end
  end
end
println("Maximum fractional error: ",emax," ",imax," ",jmax)
println("Maximum jacobian-jac_num: ",maxabs(jacobian-convert(Array{Float64,2},jac_num1)))

@test isapprox(jacobian,jac_num1;norm=maxabs)
end


text:
struct CellFinder{Tv,Ti}
    xgrid::ExtendableGrid{Tv,Ti}
    xCellFaces::Adjacency{Ti}
    xFaceCells::Adjacency{Ti}
    xCellGeometries::GridEGTypes
    facetogo::Array{Array{Ti,1},1}
    previous_cells::Array{Ti,1}
    EG::Array{Type{<:AbstractElementGeometry},1}
    L2G4EG::Vector{L2GTransformer{Tv,Ti,EG,CS} where {EG<:AbstractElementGeometry,CS<:AbstractCoordinateSystem}} # each Geometry has its own L2GTransformer
    invA::Matrix{Tv}
    xreftest::Array{Tv,1}
    cx::Vector{Tv}
end

## postprocessing function that calcualte conditions to decide over which face to go next
function postprocess_xreftest!(xreftest::Array{Tv}, ::Type{<:Edge1D}) where {Tv}
    xreftest[2] = 1 - xreftest[1]
    return nothing
end
function postprocess_xreftest!(xreftest::Array{Tv}, ::Type{<:Triangle2D}) where {Tv}
    xreftest[3] = 1 - xreftest[1] - xreftest[2]
    return nothing
end
function postprocess_xreftest!(xreftest::Array{Tv}, ::Type{<:Tetrahedron3D}) where {Tv}
    xreftest[4] = 1 - xreftest[1] - xreftest[2] - xreftest[3]
    return nothing
end
function postprocess_xreftest!(xreftest::Array{Tv}, ::Type{<:Parallelogram2D}) where {Tv}
    xreftest[3] = 1 - xreftest[1]
    xreftest[4] = 1 - xreftest[2]
    return nothing
end
function postprocess_xreftest!(xreftest::Array{Tv}, ::Type{<:Parallelepiped3D}) where {Tv}
    xreftest[4] = 1 - xreftest[1]
    xreftest[5] = 1 - xreftest[2]
    xreftest[6] = 1 - xreftest[3]
    return nothing
end

function CellFinder(xgrid::ExtendableGrid{Tv,Ti}) where {Tv,Ti}
    CS = xgrid[CoordinateSystem]
    EG = xgrid[UniqueCellGeometries]
    L2G4EG = Vector{L2GTransformer{Tv,Ti,EG,CS} where {EG<:AbstractElementGeometry,CS<:AbstractCoordinateSystem}}(undef,length(EG))
    facetogo = Array{Array{Ti,1},1}(undef,length(EG))
    xreftestlength::Int = 0
    for j = 1 : length(EG)
        L2G4EG[j] = L2GTransformer(EG[j], xgrid, ON_CELLS)
        if EG[j] <: AbstractElementGeometry1D
            facetogo[j] = [1, 2]
            xreftestlength = max(xreftestlength,2)
        elseif EG[j] <: Triangle2D
            facetogo[j] = [3, 1, 2]
            xreftestlength = max(xreftestlength,3)
        elseif EG[j] <: Tetrahedron3D
            facetogo[j] = [4, 2, 1, 3]
            xreftestlength = max(xreftestlength,4)
        elseif EG[j] <: Parallelogram2D
            facetogo[j] = [4, 1, 2, 3]
            xreftestlength = max(xreftestlength,4)
        elseif EG[j] <: Parallelepiped3D
            facetogo[j] = [5, 2, 1, 3, 4, 6]
            xreftestlength = max(xreftestlength,6)
        else
            @error "ElementGeometry not supported by CellFinder"
        end
    end

    edim = dim_element(EG[1])
    A = zeros(Tv,edim,edim)
    xreftest = zeros(Tv,xreftestlength)
    return CellFinder{Tv,Ti}(xgrid, xgrid[CellFaces], xgrid[FaceCells], xgrid[CellGeometries], facetogo, zeros(Ti,3), EG, L2G4EG, A, xreftest, zeros(Tv,edim))
end

# modification of pdelib function in gfind.cxx
function gFindLocal!(xref, CF::CellFinder{Tv,Ti}, x; icellstart::Int = 1, eps = 1e-14) where{Tv,Ti}

    # works for convex domainsand simplices only !
    xCellFaces::Adjacency{Ti} = CF.xCellFaces
    xFaceCells::Adjacency{Ti} = CF.xFaceCells
    xCellGeometries::GridEGTypes = CF.xCellGeometries
    EG::GridEGTypes = CF.EG
    cx::Vector{Tv} = CF.cx
    cEG::Int = 0
    facetogo::Array{Array{Ti,1},1} = CF.facetogo
    icell::Int = icellstart
    previous_cells::Array{Ti,1} = CF.previous_cells
    fill!(previous_cells,0)
    xreftest::Array{Tv,1} = CF.xreftest
    L2G::L2GTransformer{Tv,Ti} = CF.L2G4EG[1]
    L2Gb::Vector{Tv} = L2G.b

    invA::Matrix{Tv} = CF.invA
    imin::Int = 0

    while (true)
        # find current cell geometry index
        cEG = 1
        while xCellGeometries[icell] != EG[cEG]
            cEG += 1
        end

        # update local 2 global map
        L2G = CF.L2G4EG[cEG]
        update_trafo!(L2G, icell)
        L2Gb = L2G.b

        # compute barycentric coordinates of node
        for j = 1 : length(x)
            cx[j] = x[j] - L2Gb[j]
        end
        mapderiv!(invA,L2G,xref)
        fill!(xreftest,0)
        for j = 1 : length(x), k = 1 : length(x)
            xreftest[k] += invA[j,k] * cx[j]
        end
        postprocess_xreftest!(xreftest,xCellGeometries[icell])

        # find minimal barycentric coordinate with
        imin = 1
        for i = 2 : length(facetogo[cEG])
            if xreftest[imin] >= xreftest[i]
                imin = i
            end
        end

        # if all barycentric coordinates are within [0,1] the including cell is found
        if xreftest[imin] >= -eps
            xref .= view(xreftest,1:length(xref))
            return icell
        end

        # otherwise: go into direction of minimal barycentric coordinates
        for j = 1 : length(previous_cells)-1
            previous_cells[j] = previous_cells[j+1]
        end
        previous_cells[end] = icell
        icell = xFaceCells[1,xCellFaces[facetogo[cEG][imin],icell]]
        if icell == previous_cells[end]
            icell = xFaceCells[2,xCellFaces[facetogo[cEG][imin],icell]]
            if icell == 0
                @debug  "could not find point in any cell and ended up at boundary of domain (maybe x lies outside of the domain ?)"
                return 0
            end
        end

        if icell == previous_cells[end-1]
            @debug  "could not find point in any cell and ended up going in circles (better try brute force search)"
            return 0
        end
    end

    return 0
end

function gFindBruteForce!(xref, CF::CellFinder{Tv,Ti}, x; eps = 1e-14) where {Tv,Ti}

    cx::Vector{Tv} = CF.cx
    cEG::Int = 0
    EG::GridEGTypes = CF.EG
    xreftest::Array{Tv,1} = CF.xreftest
    xCellGeometries::GridEGTypes = CF.xCellGeometries
    facetogo::Array{Array{Ti,1},1} = CF.facetogo
    L2Gb::Vector{Tv} = CF.L2G4EG[1].b
    invA::Matrix{Tv} = CF.invA
    imin::Int = 0

    for icell = 1 : num_sources(CF.xgrid[CellNodes])

        # find current cell geometry index
        cEG = 1
        while xCellGeometries[icell] != EG[cEG]
            cEG += 1
        end

        # update local 2 global map
        L2G = CF.L2G4EG[cEG]
        update_trafo!(L2G, icell)
        L2Gb = L2G.b

        # compute barycentric coordinates of node
        for j = 1 : length(x)
            cx[j] = x[j] - L2Gb[j]
        end
        mapderiv!(invA,L2G,xref)
        fill!(xreftest,0)
        for j = 1 : length(x), k = 1 : length(x)
            xreftest[k] += invA[j,k] * cx[j]
        end
        postprocess_xreftest!(xreftest,xCellGeometries[icell])

        # find minimal barycentric coordinate with
        imin = 1
        for i = 2 : length(facetogo[cEG])
            if xreftest[imin] >= xreftest[i]
                imin = i
            end
        end

        # if all barycentric coordinates are within [0,1] the including cell is found
        if xreftest[imin] >= -eps
            xref .= view(xreftest,1:length(xref))
            return icell
        end
    end

    @debug "gFindBruteForce did not find any cell that contains x = $x (make sure that x is inside the domain, or try reducing $eps)"

    return 0
end


text:
using Turing, Random, Test
using DynamicPPL: @varname, getsym
using Distributions: Normal
using StatsFuns

dir = splitdir(splitdir(pathof(Turing))[1])[1]
include(dir*"/test/test_utils/AllUtils.jl")

@turing_testset "util.jl" begin
    i = 1
    vn = @varname s
    @test getsym(vn) === :s
    @test vn.indexing == ""

    vn = @varname x[1,2][1+5][45][3][i]
    @test getsym(vn) === :x
    @test vn.indexing == "[1,2][6][45][3][1]"
    @test StatsFuns.logistic(1.1) == 1.0 / (exp(-1.1) + 1.0)
    @test StatsFuns.logit(0.3) ≈ -0.8472978603872036 atol=1e-9
    @test !isnan(StatsFuns.logit(1.0))
    @test isinf(StatsFuns.logit(1.0))
    @test !isnan(StatsFuns.logit(0.0))
    @test isinf(StatsFuns.logit(0.0))
end


text:
import HTTP

export method, GET, route, page, probability, query

# Request type

method(m::AbstractString, app...) = branch(req -> req[:method] == m, app...)
method(ms, app...) = branch(req -> req[:method] in ms, app...)

GET(app...) = method("GET", app...)

# Path routing

splitpath(p::AbstractString) = split(p, "/", keepempty=false)
splitpath(p) = p

function matchpath(target, path)
  length(target) > length(path) && return
  params = d()
  for i = 1:length(target)
    if startswith(target[i], ":")
      params[Symbol(target[i][2:end])] = path[i]
    else
      target[i] == path[i] || return
    end
  end
  return params
end

function matchpath!(target, req)
  ps = matchpath(target, req[:path])
  ps == nothing && return false
  merge!(params!(req), ps)
  splice!(req[:path], 1:length(target))
  return true
end

route(p::Vector, app...) = branch(req -> matchpath!(p, req), app...)
route(p::AbstractString, app...) = route(splitpath(p), app...)
route(app...) = route([], app...)
route(app::Base.Callable, p) = route(p, app)
route(app1::Base.Callable, app2::Base.Callable) = route([], app1, app2)

page(p::Vector, app...) = branch(req -> length(p) == length(req[:path]) && matchpath!(p, req), app...)
page(p::AbstractString, app...) = page(splitpath(p), app...)
page(app...) = page([], app...)
page(app::Base.Callable, p) = page(p, app)
page(app1::Base.Callable, app2::Base.Callable) = page([], app1, app2)

# Query routing

function matchquery(q, req)
  qdict = HTTP.URIs.queryparams(req[:query])
  length(q) != length(qdict) && return false
  for (key, value) in q
    if haskey(qdict, key) && (value == "" || value == qdict[key])
      continue
    else
      return false
    end
  end
  return true
end

query(q::Dict{<:AbstractString, <:AbstractString}, app...) =
  branch(req -> matchquery(q, req), app...)

# Misc

probability(x, app...) = branch(_->rand()<x, app...)
# Old typo
@deprecate probabilty(x, app...) probability(x, app...)


text:
abstract type CondMixture end

struct Tsederberg <: CondMixture end

function k_fluid(Mud::DrillFluid, model::Tsederberg)

end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 3f54a71c-06aa-11ec-3712-0dbdbf56b590
# edit the code below to set your name and Cornell ID (i.e. email without @cornell.edu)

student = (name = "Kyra Frederick", cornell_id = "kcf29")

# ╔═╡ ecc6dc76-32bc-46cd-9d81-67fe85afc975
md""" make sure to check if this has been updated"""

# ╔═╡ daf6f050-6fe5-474b-b83c-bc14870b6efb
md"""
# Assignment 5

A river is receiving waste discharges from two sources that are 15km apart, as shown in the Figure below. Characteristics of the river inflow and waste discharges are given in the table, and the values of ``k_a``, ``k_c``, and ``k_n`` are 0.55, 0.35, and 0.25 per day, respectively.

"""

# ╔═╡ 409babef-def0-49a6-a82f-552e3dfb0a3a
md"""
**Table 1**
"""

# ╔═╡ 4651f993-b5da-4c04-9d0e-2c021b3d2732
md"""
| **Parameter** | **River Inflow** | **Waste Stream 1** | **Waste Stream 2 ** |
|:-------------:|:------------:|:-------------:|:---------:|
| ``Q`` 		 		| 100,000 ``m^3d`` | 10,000 ``m^3d`` | 15,000 ``m^3d`` |
| ``DO`` 		 		| 350 mg/L 		   | 5 mg/L 		 | 5 mg/L  		   |
| ``CBOD``		 		| 150 mg/L 		   | 50 mg/L 		 | 45 mg/L  	   |
| ``NBOD``		 		| 150 mg/L 		   | 35 mg/L 		 | 35 mg/L  	   |
| ``DO_{Saturated}``		| 150 mg/L 		   |     -  		 |    -   		   |
"""

# ╔═╡ 43b15cf7-0b54-4873-8ef4-828e63de6a95
md"""
## Part 1

Determine how long (in distance) it takes for the DO of the river to recover to the regulated level of 6mg/L if the waste is untreated.
WS1 = Waste stream 1
WS2= waste stream 2

``[DO]=[DO_{0}]e^{-k_{a}t_{W1}}``

``DO_{@x=0}= (Q_{river}DO_{river}+Q_{WS1}DO_{WS1})/(Q_{river}+Q_{WS1})``

``DO_{@x=14.99999}= [DO_{@x=0}]e^{-k_{a}t}``

* ``t_{W1}=A_{c}*L/Q=A_{c}*15,000m/(100,000+10,000m^3/d)=``

``DO_{@x=15}= ((Q_{river}+Q_{WS1}) DO_{@x=14.99999}+Q_{WS2}DO_{WS2})/((Q_{river}+Q_{WS1})+Q_{WS2})``

``DO_{@x>15}= [DO_{@x=15}]e^{-k_{a}t_{W2}}``

Ans: ``t_{W2}=t_{W1} * ln(6mg/L/DO_{@x=15})/-k_{a}``

"""

# ╔═╡ 9f1fcbbb-c856-42e9-a53d-fbfd7a0846de
md"""
## Part 2
What is the minimum level of treatment (% removal) for waste stream 2 that will ensure that the ``DO`` level never drops below 4 mg/L, with waste stream 1 untreated?

"""

# ╔═╡ b4bda1ca-25d4-4384-a054-f7a18cd2af29
md"""
## Part 3:
What is the minimum level of equal treatment (% removal) of the two sources that will ensure that the ``DO`` level never drops below 4 mg/L?

"""

# ╔═╡ c682b9cf-1c5f-43d3-aa5d-abe34469fc5d
md"""
## Part 4:
If the ``CBOD`` values entering the river upstream are known to vary uniformly between ``\pm 20\%`` of the value given in the table, what is the (approximate) likelihood that the ``DO`` will drop below 4 mg/L under the treatment strategy you proposed for Part 3?
"""

# ╔═╡ 0cbd3384-1ff8-434e-b081-8a0ec326e1a3
md"""
Use this cell to discuss your results. Feel free to add new cells to query any quantities of interest from the model.
"""

# ╔═╡ b979c498-3335-4766-9025-cc402ff91a37
md"""
## Citations

If you used any external resource in the process of solving this problem, provide a reference below. These could be other people, websites, datasets, or code examples.
"""

# ╔═╡ 459f424f-3700-4832-a665-c9c5659de2b1
md"""
List external resources here.
"""

# ╔═╡ Cell order:
# ╠═3f54a71c-06aa-11ec-3712-0dbdbf56b590
# ╠═ecc6dc76-32bc-46cd-9d81-67fe85afc975
# ╟─daf6f050-6fe5-474b-b83c-bc14870b6efb
# ╟─409babef-def0-49a6-a82f-552e3dfb0a3a
# ╟─4651f993-b5da-4c04-9d0e-2c021b3d2732
# ╠═43b15cf7-0b54-4873-8ef4-828e63de6a95
# ╟─9f1fcbbb-c856-42e9-a53d-fbfd7a0846de
# ╟─b4bda1ca-25d4-4384-a054-f7a18cd2af29
# ╟─c682b9cf-1c5f-43d3-aa5d-abe34469fc5d
# ╠═0cbd3384-1ff8-434e-b081-8a0ec326e1a3
# ╠═b979c498-3335-4766-9025-cc402ff91a37
# ╠═459f424f-3700-4832-a665-c9c5659de2b1

max_stars_repo_name: kyfrederick/hmfaysal-omega-theme

max_issues_repo_name: kyfrederick/hmfaysal-omega-theme

max_forks_repo_name: kyfrederick/hmfaysal-omega-theme


text:
val(x::CompoundPeriod) = isempty(x) ? 0 : canonicalize(x)
val(x) = isempty(x) ? 0 : value(x)

function utc_delta()
    utctime, localtime = now(UTC), now()
    delta_minutes = floor(localtime - utctime, Minute)
    delta_hour = floor(delta_minutes, Hour)
    delta_hour_value = value(delta_hour)
    delta_minute = canonicalize(delta_minutes - delta_hour)
    delta_minute_value = val(delta_minute)
    delta_sign = '+'
    if delta_hour_value < 0
        delta_sign = '-'
    end
    delta_hour_str = string(abs(delta_hour_value))
    if abs(delta_hour_value) < 10
       delta_hour_str = "0" * delta_hour_str
    end
    delta_minute_str = string(delta_minute_value)
    if length(delta_minute_str) < 2
        delta_minute_str = "0" * delta_minute_str
    end
    delta_time_str = delta_sign * delta_hour_str * ":" * delta_minute_str
    delta_time_str
end

const LOCAL_TZ_DELTA = utc_delta()

function timestamp(nd::NanoDate; sep="", utctime=true, localtime=false, postfix=true)
    if localtime == true
        utctime = false
        postfix = true
    elseif postfix == true
        utctime = true
    else
        utctime = false
    end

    if utctime
        postfix = 'Z'
    elseif localtime
        postfix = LOCAL_TZ_DELTA
    else
        postfix = ""
    end

    tstamp = string(nd; sep) * postfix
    tstamp
end


text:
"""
    StateVectorType

A trait defined on an `AbstractFunction`, describing how state vectors are composed
and how the error state, or the "difference" between two vectors, is computed.
These "error states" should be Euclidean vectors, and can have a different dimension than
the original state vectors.

The following types are provided by RobotDynamics:

* [`EuclideanState`](@ref)
* [`RotationState`](@ref)

## Defining a new state vector type
To define a custom state vector type (e.g. `NewStateType`), you need to implement
the function that calculates the difference between state vectors, as well the
"error state Jacobian" and it's derivative.
First, you must specify the dimension of the error state by defining

    errstate_dim(::NewStateType, fun::AbstractFunction)

Then you define method that calculates the error state between two state vectors:

    state_diff!(::NewStateType, model, dx, x, x0)

where `dx` should have dimension `e = output_dim(model)` and `x` and `x0` should
have dimension `n = state_dim(model)`.

The other piece of information we need is the "error state Jacobian." This is the
first-order approximation of the mapping from local error state values to the true
state vector. Let `x ⊗ x0` be the composition of two state vectors, and let `φ` be
the function that maps a Euclidean error state `δx` to a state vector. The error
state Jacobian is the Jacobian of `x ⊗ φ(δx)` with respect to `δx`, taking the limit
as `δx → 0`. The resulting function is solely a function of the state `x`, and is
implemented with the following function:

    errstate_jacobian!(::NewStateType, model, G, x)

where `G` has dimensions `(n,e)`.

We'll also need the derivative of this function. For computational efficiency,
we compute the Jacobian of the error state Jacobian transposed with some vector `x̄`,
where `x̄` is state vector of size `(n,)`. More precisely, it is the Hessian of

```math
\\bar{x}^T (x \\otimes \\varphi(\\delta x))
````
This function is calculated using

    ∇errstatef_jacobian!(::NewStateType, model, ∇G, x, xbar)

where `∇G` is of size `(e,e)`.
"""
abstract type StateVectorType end

"""
    EuclideanState

The space of standard Euclidean vectors, where composition is element-wise addition,
and the error state is just element-wise subtraction.
"""
struct EuclideanState <: StateVectorType end

"""
    RotationState

The space of vectors composes of Euclidean states intermixed with 3D rotations,
represented using any 3 or 4-parameter representation (usually unit quaternions).
The error state Jacobians and error state calculation are all calculating using
[Rotations.jl](https://github.com/JuliaGeometry/Rotations.jl).
"""
struct RotationState <: StateVectorType end

"""
    statevectortype

Queries the StateVectorType trait on a function.
"""
statevectortype(::Type{<:AbstractFunction}) = EuclideanState()
statevectortype(fun::AbstractFunction) = statevectortype(typeof(fun))

# Dispatch on `statevectortype` trait
"""
    errstate_dim

The length of the Euclidean error state for the given state vector type.
"""
errstate_dim(fun::AbstractFunction) = errstate_dim(statevectortype(fun), fun)

"""
    state_diff!(fun, dx, x, x0)

Calculate the error state between `x` and `x0`. The Euclidean error state `dx` should
have dimension equal to `output_dim(fun)` while `x` and `x0` should have dimension
equal to `state_dim(fun)`. For Euclidean states, this is just

```math
\\delta x = x - x_0
```
"""
state_diff!(fun::AbstractFunction, dx, x, x0) = state_diff!(statevectortype(fun), fun, dx, x, x0)

state_diff(fun::AbstractFunction, x, x0) = state_diff(statevectortype(fun), fun, x, x0)

"""
    errstate_jacobian!

The Jacobian of

```math
x \\otimes \\varphi(\\delta x)
```
with respect to ``\\delta x``.
Here ``\\otimes`` is the composition function for two elements of the state space,
and ``\\varphi(\\delta x)`` is the function that elements of the Euclidean error state to
the state vector space.
"""
errstate_jacobian!(fun::AbstractFunction, J, z::AbstractKnotPoint) =
    errstate_jacobian!(statevectortype(fun), fun, J, state(z))
errstate_jacobian!(fun::AbstractFunction, J, x) =
    errstate_jacobian!(statevectortype(fun), fun, J, x)

"""
    ∇errstate_jacobian!

The Hessian of

```math
\\bar{x}^T (x \\otimes \\varphi(\\delta x))
```
with respect to ``\\delta x``.
Here ``\\otimes`` is the composition function for two elements of the state space,
and ``\\varphi(\\delta x)`` is the function that elements of the Euclidean error state to
the state vector space. The vector ``\\bar{x}`` is some element of the state space.
"""
∇errstate_jacobian!(fun::AbstractFunction, J, z::AbstractKnotPoint, dx) =
    ∇errstate_jacobian!(statevectortype(fun), fun, J, state(z), dx)

∇errstate_jacobian!(fun::AbstractFunction, J, x, dx) =
    ∇errstate_jacobian!(statevectortype(fun), fun, J, x, dx)

# Euclidean state vectors
errstate_dim(::EuclideanState, fun::AbstractFunction) = state_dim(fun)
state_diff!(::EuclideanState, fun::AbstractFunction, dx, x, x0) = dx .= x .- x0
state_diff(::EuclideanState, fun::AbstractFunction, x, x0) = x - x0
∇errstate_jacobian!(::EuclideanState, fun::AbstractFunction, ∇G, x, dx) = ∇G .= 0

function errstate_jacobian!(::EuclideanState, fun::AbstractFunction, J, x)
    J .= 0
    for i = 1:size(J,1)
        J[i,i] = one(eltype(J))
    end
    return nothing
end


text:
################################################
# Sections of the code below were proposed by Simon Danisch:
#     https://gist.github.com/SimonDanisch/ae046f3a0c78b26242e78fa9b139aa11#file-benchmark-jl
################################################

using BenchmarkTools

"""
    regular_time_step(u::Matrix)

Find the Jacobi numerical approximation using un-optimized loops.
"""
function regular_time_step(u::Matrix)
    n, m = size(u)
    err = 0.0
    for i in 2:n-1
        for j in 2:m-1
            tmp = u[i,j]
            u[i,j] = ((u[i-1, j] + u[i+1, j] +
                       u[i, j-1] + u[i, j+1])*4.0 +
                       u[i-1,j-1] + u[i-1,j+1] +
                       u[i+1,j-1] + u[i+1,j+1])/20.0

            diff = u[i,j] - tmp
            err = err + diff*diff
        end
    end
    return sqrt(err)
end

"""
    optimized_time_step{T}(u::Matrix{T})

Find the Jacobi numerical approximation using optimized loops.
"""
function optimized_time_step{T}(u::Matrix{T})
    n, m = size(u)
    err = T(0.0) # not an optimization, but makes it work for all element types
    @inbounds for i = 2:n-1
        for j = 2:m-1
            tmp = u[j, i];
            u[j, i] = ( T(4) *(u[j-1, i] + u[j+1, i] +
                        u[j, i-1] + u[j, i+1]) +
                        u[j-1, i-1] + u[j+1, i+1] +
                        u[j+1, i-1] + u[j-1, i+1]) / T(20)
            diff = u[j, i] - tmp
            err += diff*diff;
        end
    end
    return sqrt(err)
end

"""
    optimized_time_step{T}(u::Matrix{T})

Find the Jacobi numerical approximation using optimized loops and SIMD.
"""
function optimized_time_step_simd{T}(u::Matrix{T})
    n, m = size(u)
    err = T(0.0) # not an optimization, but makes it work for all element types
    @inbounds for i = 2:n-1
        @simd for j = 2:m-1
            tmp = u[j, i];
            u[j, i] = ( T(4) *(u[j-1, i] + u[j+1, i] +
                        u[j, i-1] + u[j, i+1]) +
                        u[j-1, i-1] + u[j+1, i+1] +
                        u[j+1, i-1] + u[j-1, i+1]) / T(20)
            diff = u[j, i] - tmp
            err += diff*diff;
        end
    end
    return sqrt(err)
end

"""
    optimized_vectorized_time_step(u::Matrix)

Find the Jacobi numerical approximation using optimized vectorization.
"""
function optimized_vectorized_time_step(u::Matrix)
    (n, m) = size(u)
    u_old = copy(u)
    #@views u[2:n-1, 2:m-1] = ((u[1:n-2, 2:m-1] + u[3:n,   2:m-1] +
    u[2:n-1, 2:m-1] = ((u[1:n-2, 2:m-1] + u[3:n,   2:m-1] +
                        u[2:n-1, 1:m-2] + u[2:n-1, 3:m])*4.0 +
                        u[1:n-2, 1:m-2] + u[1:n-2, 3:m] +
                        u[3:n,   1:m-2] + u[3:n,   3:m])/20.0
    v = vec(u - u_old)
    return sqrt(dot(v, v))
end

# C version
const clib_name5 = tempname()
source = """
#include <math.h>
double c_time_step(int N, double *u[])
{
    double tmp, diff;
    int i, j;
    double err;
    err = 0.0;
    for (i = 1; i < N-1; i++) {
        for (j = 1; j < N-1; j++) {
            tmp = u[j][i];
            u[j][i] = ( 4.0 *(u[j-1][i] + u[j+1][i] +
                        u[j][i-1] + u[j][i+1]) +
                        u[j-1][i-1] + u[j+1][i+1] +
                        u[j+1][i-1] + u[j-1][i+1])/20.;
            diff = u[j][i] - tmp;
            err += diff*diff;
        }
    }
    return sqrt(err);
}
"""
run(`gcc --version`)
# gcc (Ubuntu 5.4.1-2ubuntu1~16.04) 5.4.1 20160904
open(`gcc -fPIC -Ofast -O3 -xc -shared -o $(clib_name5).so -`, "w") do pipe
    print(pipe, source)
end
function c_time_step(X::Vector{Vector{Float64}})
    ccall(("c_time_step", clib_name5), Float64, (Cint, Ptr{Ptr{Float64}}), length(X), X)
end

function solver(u, time_step)
    iter = 0
    err = 2.0
    while iter < 100000 && err > 1e-6
        err = time_step(u)
        iter += 1
    end
    u, err, iter
end

#########################################
# Get the number of points from the command line
#########################################
n, = size(ARGS)
if (n < 1)
   println("Usage: jacobiDriver4-benchmark.jl N")
   println("       --->  Please specify the number of points.")
   exit()
end

num_points = parse(Int, ARGS[1])

println("-------------------------------")
println(@sprintf "Jacobi Iterations: %d" num_points)
println("-------------------------------")

#function init_array(numPoints = 100)
function init_array(numPoints = num_points)
    pi_c = 4.0 * atan(1.0)
    v = zeros(Float64, numPoints, numPoints)
    x = [(i-1)*pi_c/(numPoints-1) for i in 1:numPoints]
    v[1,:] = [sin(a) for a in x]
    v[numPoints,:] = [sin(a)*exp(-pi_c) for a in x]
    v
end
init_c() = (v = init_array(); v_2 = [v[i, :] for i=1:size(v, 1)])

funcs = (regular_time_step, optimized_time_step, optimized_vectorized_time_step, optimized_time_step_simd) #
benchmarks = map(funcs) do time_step
    v = init_array()
    solver(v, time_step) # compile
    v = init_array()
    # should be enough iterations already to be okay with just a simple elapsed
    result = @elapsed begin
        u, err, it = solver(v, time_step)
    end
    Symbol(time_step) => (result, it)
end

v = init_c()
result = @elapsed begin
    u, err, it = solver(v, c_time_step)
end
benchmarks = (benchmarks..., :c_time_step => (result, it))
for (name, (t, sweeps)) in benchmarks
    println(name, ":")
    println("          Number of sweeps: ", sweeps)
    println("              Elapsed Time: $(t) Seconds")
end


text:
struct SOSPolynomialBridge{T, F <: MOI.AbstractVectorFunction,
                           DT <: AbstractSemialgebraicSet,
                           VBS <: AbstractVariableBridge, MCT,
                           BT <: PolyJuMP.AbstractPolynomialBasis,
                           MT <: MP.AbstractMonomial,
                           MVT <: AbstractVector{MT}} <: MOIB.AbstractBridge
    variable_bridge::VBS
    certificate_monomials::MVT
    zero_constraint::MOI.ConstraintIndex{F, PolyJuMP.ZeroPolynomialSet{DT, BT, MT, MVT}}
end

function SOSPolynomialBridge{T, F, DT, VBS, MCT, BT, MT, MVT}(
    model::MOI.ModelLike, f::MOI.AbstractVectorFunction,
    s::SOSPolynomialSet{<:AbstractAlgebraicSet}) where {
        # Need to specify types to avoid ambiguity with the default constructor
        T, F <: MOI.AbstractVectorFunction, DT <: AbstractSemialgebraicSet,
        VBS <: AbstractVariableBridge, MCT,
        BT <: PolyJuMP.AbstractPolynomialBasis, MT <: MP.AbstractMonomial,
        MVT <: AbstractVector{MT}
    }
    @assert MOI.output_dimension(f) == length(s.monomials)
    p = MP.polynomial(collect(MOIU.eachscalar(f)), s.monomials)
    # FIXME convert needed because the coefficient type of `r` is `Any` otherwise if `domain` is `AlgebraicSet`
    r = convert(typeof(p), rem(p, ideal(s.domain)))
    X = monomials_half_newton_polytope(MP.monomials(r), s.newton_polytope)
    Q, variable_bridge = add_matrix_variable_bridge(model, MCT, length(X), T)
    g = build_gram_matrix(Q, X)
    q = r - g
    set = PolyJuMP.ZeroPolynomialSet(s.domain, s.basis, MP.monomials(q))
    coefs = MOIU.vectorize(MP.coefficients(q))
    zero_constraint = MOI.add_constraint(model, coefs, set)
    return SOSPolynomialBridge{T, F, DT, VBS, MCT, BT, MT, MVT}(
        variable_bridge, X, zero_constraint)
end

function MOI.supports_constraint(::Type{SOSPolynomialBridge{T}},
                                 ::Type{<:MOI.AbstractVectorFunction},
                                 ::Type{<:SOSPolynomialSet{<:AbstractAlgebraicSet}}) where T
    return true
end
function MOIB.added_constraint_types(::Type{SOSPolynomialBridge{T, F, DT, VBS, MCT, BT, MT, MVT}}) where {T, F, DT, VBS, MCT, BT, MT, MVT}
    added = [(F, PolyJuMP.ZeroPolynomialSet{DT, BT, MT, MVT})]
    return append_added_constraint_types(added, MCT, T)
end
function MOIB.concrete_bridge_type(::Type{<:SOSPolynomialBridge{T}},
                                   F::Type{<:MOI.AbstractVectorFunction},
                                   ::Type{<:SOSPolynomialSet{DT, CT, <:PolyJuMP.MonomialBasis, MT, MVT}}) where {T, DT<:AbstractAlgebraicSet, CT, MT, MVT}
    # promotes VectorOfVariables into VectorAffineFunction, it should be enough
    # for most use cases
    G = MOIU.promote_operation(-, T, F, MOI.VectorOfVariables)
    MCT = matrix_cone_type(CT)
    VBS = union_vector_bridge_types(MCT, T)
    return SOSPolynomialBridge{T, G, DT, VBS, MCT, PolyJuMP.MonomialBasis, MT, MVT}
end

# Attributes, Bridge acting as an model
function MOI.get(bridge::SOSPolynomialBridge, attr::MOI.NumberOfVariables)
    return MOI.get(bridge.variable_bridge, attr)
end
function MOI.get(::SOSPolynomialBridge{T, F, DT, VBS, MCT, BT, MT, MVT},
                 ::MOI.NumberOfConstraints{F, PolyJuMP.ZeroPolynomialSet{DT, BT, MT, MVT}}) where {
        # Need to specify types to avoid ambiguity with the method redirecting
        # to `variable_bridge`
        T, F <: MOI.AbstractVectorFunction, DT <: AbstractSemialgebraicSet,
        VBS <: AbstractVariableBridge, MCT <: MOI.AbstractVectorSet,
        BT <: PolyJuMP.AbstractPolynomialBasis, MT <: MP.AbstractMonomial,
        MVT <: AbstractVector{MT}
    }
    return 1
end
function MOI.get(b::SOSPolynomialBridge{T, F, DT, VBS, MCT, BT, MT, MVT},
                 ::MOI.ListOfConstraintIndices{F, PolyJuMP.ZeroPolynomialSet{DT, BT, MT, MVT}}) where {
        # Need to specify types to avoid ambiguity with the method redirecting
        # to `variable_bridge`
        T, F <: MOI.AbstractVectorFunction, DT <: AbstractSemialgebraicSet,
        VBS <: AbstractVariableBridge, MCT <: MOI.AbstractVectorSet,
        BT <: PolyJuMP.AbstractPolynomialBasis, MT <: MP.AbstractMonomial,
        MVT <: AbstractVector{MT}
    }
    return [b.zero_constraint]
end
function MOI.get(bridge::SOSPolynomialBridge,
                 attr::MOI.NumberOfConstraints)
    return MOI.get(bridge.variable_bridge, attr)
end
function MOI.get(bridge::SOSPolynomialBridge,
                 attr::MOI.ListOfConstraintIndices)
    return MOI.get(bridge.variable_bridge, attr)
end

# Indices
function MOI.delete(model::MOI.ModelLike, bridge::SOSPolynomialBridge)
    # First delete the constraints in which the Gram matrix appears
    MOI.delete(model, bridge.zero_constraint)
    # Now we delete the Gram matrix
    MOI.delete(model, bridge.variable_bridge)
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike,
                 attr::MOI.ConstraintDual,
                 bridge::SOSPolynomialBridge)
    return MOI.get(model, attr, bridge.zero_constraint)
end
function MOI.get(::MOI.ModelLike, ::CertificateMonomials,
                 bridge::SOSPolynomialBridge)
    return bridge.certificate_monomials
end
function MOI.get(model::MOI.ModelLike,
                 attr::GramMatrixAttribute,
                 bridge::SOSPolynomialBridge)
    return build_gram_matrix(MOI.get(model, attr, bridge.variable_bridge),
                             bridge.certificate_monomials)
end
function MOI.get(model::MOI.ModelLike,
                 attr::MomentMatrixAttribute,
                 bridge::SOSPolynomialBridge)
    return build_moment_matrix(MOI.get(model, attr, bridge.variable_bridge),
                               bridge.certificate_monomials)
end


text:
@testset "SME" begin
    res = BRS.sme(collect(1:10:100), collect(1:10:100))
    @test res.value == 1
    @test res.nUp == 1
    @test res.nDown == 0
end

@testset "RMSSD ratio" begin
    res = BRS.rmssdr(collect(1:10:100), collect(1:10:100))
    @test res == 1
end

@testset "xBRS" begin
    res = BRS.xbrs(collect(1000:1000:10000), collect(1:10:100))
    @test res.value ≈ 100
end

@testset "tfBRS" begin
    res = BRS.tfbrs(collect(100:100:10000), collect(1:1:100), n = 40)
    @test res.value ≈ 100
end

@testset "prsaBRS" begin
    res = BRS.prsabrs(collect(100:100:1000000), collect(1:1:10000))
    @test isapprox(res.value, 100, atol = 0.01)
end

@testset "arBRS" begin
    res = BRS.arbrs(collect(1.0:1000.0), collect(1.0:1000.0))
    @test isapprox(res.αHF, 1)
    @test isnan(res.αLF)
end

@testset "getBRS" begin

    res = getBRS(collect(1.0:1000.0), collect(1.0:1000.0))
    @test length(res) == 7
end


text:
using Revise
using StatsBase, KernelDensity, Distributions
dir = "C:/Users/stravaglino3/Downloads/julia_code/src/"
includet(dir*"models.jl")
includet(dir*"models_dissociation.jl")
includet(dir*"helper_functions.jl")
includet(dir*"binning_helper_fs.jl")
includet(dir*"mle_fitting.jl")
includet(dir*"mle_plotting.jl")
using DataFrames, Plots, ColorSchemes
data_dir = "C:/Users/stravaglino3/Downloads/julia_code/CD3gd/data/"
CD3_data = CSV.read(data_dir*"dissociations_cleaned.csv",DataFrame)

γϵ_binned = CSV.read(data_dir*"data/rebinned/cd3ge-tcr-rebin binned bond lifetime.csv",DataFrame)
δϵ_binned = CSV.read(data_dir*"cd3de-tcr-unbinned.txtbinned bond lifetime.csv",DataFrame)
mix_binned = CSV.read(data_dir*"data/rebinned/cd3mix-tcr-rebin.txtbinned bond lifetime.csv",DataFrame)

γϵ = delete_missing(CD3_data[:,2:3])
δϵ = delete_missing(CD3_data[:,8:9])
mix = delete_missing(CD3_data[:,14:15])

pyplot()
gr()
s1=Plots.scatter(γϵ[1],γϵ[2], label="γϵ");
s2=Plots.scatter(δϵ[1],δϵ[2], label="δϵ");
s3=Plots.scatter(mix[1],mix[2], label="mix");
plot(s1,s2,s3)

γϵ  = F_t_Data(γϵ[1],γϵ[2],   InterpKDE(kde(γϵ[1])))
δϵ  = F_t_Data(δϵ[1],δϵ[2],   InterpKDE(kde(δϵ[1])))
mix = F_t_Data(mix[1],mix[2], InterpKDE(kde(mix[1])))

h1=histogram2d(γϵ.F, γϵ.t, nbins=(40,50),title="γϵ")#, c=:jet);#,aspect_ratio=2)
h2=histogram2d(δϵ.F, δϵ.t, nbins=(40,60),title="δϵ")#, c=:turbo);#,aspect_ratio=2)
h3=histogram2d(mix.F, mix.t, nbins=(40,50),title="mix")#, c=:turbo);#,aspect_ratio=2)
plot(h1,h2,h3)

#________________________________________________________________________________________________
## try fitting with 2 species model
# rates from Dunn's Science paper:
r_ini = [14.0,0.1,
        0.14,0.4,
        5.0,0.2,
        20.0,4.0]

# γϵ.t = γϵ.t .- minimum(γϵ.t)

model = bimolecular_diss_bell!
# u₀ = [1.0,0.0]
u₀ = [] # don't provide initial states
u₀_type = "equilibrium" # calculate initial states based on equilibrium
# u₀ = "detailed_balance"
# u₀ = "force"
# u₀ = []
tspan = (0.0,maximum(γϵ.t))
γϵ_ft_M = M(model,u₀,u₀_type,tspan,[],[],0.0)

pinit = [14.1,0.01,
        0.1,0.6,
        3.,0.2,
        150.0,4.0]

opt =  [(14.05,14.15),(0.,10.),
        (0.,1.),(0.,10.),
        (0.,10.),(0.,10.),
        (0.,500.),(0.,10.)]

# pinit = rand(8)*10
@time mle_loss(pinit,γϵ,γϵ_ft_M)
#=
res = so.minimize(x->mle_loss(10 .^ x,γϵ,γϵ_ft_M), log10.(pinit),method="COBYLA",
                        bounds=opt,
                        options=Dict("maxiter"=>10000,
                                 "disp" => true) )
res["x"]
=#   # the random seed

@time rates_γϵ, l_γϵ, hes_γϵ = optimization(γϵ_ft_M,γϵ,opt)
using LinearAlgebra
errors = sqrt.(abs.( diag(inv(hes_γϵ)) ))
# ks_γϵ = 10 .^ rates_γϵ
ks_γϵ = rates_γϵ
ks_γϵ = pinit

# pyplot()
plotly()
dt = 0.1
ts = 0:dt:maximum(γϵ.t)
Fs = 0:1:maximum(γϵ.F)
# Fs = 0:5:150
ks_γϵ = [14.1,0.0,
        0.1,0.6,
        3.,0.2,
        150.0,4.0]

sols_γϵ = plot_ft_fit(γϵ_ft_M,ks_γϵ,γϵ)

plot(Fs,mean_lifetime.(sols_γϵ),ylabel="<t>",xlabel="F",title="γϵ")
plot_F_vs_t(clean_dissociation_data(γϵ_binned))

us = map(x -> initial_u(ks_γϵ,γϵ_ft_M,x), 0:30)
scatter(0:30,reduce(hcat,us)')

plot_bell_equations([11.0,0.0,
                     0.14,0.4,
                     3.0,0.2,
                     20.0,4.0])
plt = plot_bell_equations(vcat(k₋ᵒf_γϵ, ks_γϵ))
# display(plt)
pyplot()
u₀_γϵ = initial_state(ks_γϵ[1],ks_γϵ[3],ks_γϵ[5],ks_γϵ[7])
initial_state(ks_γϵ[5],ks_γϵ[7])
a
________________________________________________________________________________________________
## δϵ fitting with 2 species model
model = bimolecular_diss_bell!
# u₀ = [1.0,0.0]
u₀ = "equilibrium"
tspan = (0.0,maximum(δϵ.t))
δϵ_ft_M = M(model,u₀,tspan,[],[],0.0)
pinit = [11.1,0.1,
        0.1,0.5,
        0.1,0.5,
        10.0,5.0]
opt =  [(11.05,11.15),(0.,10.),
        (0.,1.),(0.,10.),
        (0.,10.),(0.,10.),
        (0.,100.),(0.,10.)]
@time rates_δϵ, l_δϵ, hes_δϵ = optimization(δϵ_ft_M,δϵ,opt)

ks_δϵ = rates_δϵ

ks_δϵ = [11.087136651210209, 0.0026318767020595732,
        0.02256859725182868, 0.5066206884832165,
        0.21328859839615166, 0.058733617801501145,
        1.5349514553341914, 0.6116955164634462]

plotly()
ts = 0:dt:maximum(δϵ.t)
Fs = 0:1:maximum(δϵ.F)
sols_δϵ = plot_ft_fit(δϵ_ft_M,ks_δϵ,δϵ)
pyplot()
plotly()
plot(Fs,mean_lifetime.(sols_δϵ),ylabel="<t>",xlabel="F",title="δϵ")
plot_F_vs_t(clean_dissociation_data(δϵ_binned))

pyplot()
plot_bell_equations(ks_δϵ)
u₀_δϵ = initial_state(ks_δϵ[1],ks_δϵ[3],ks_δϵ[5],ks_δϵ[7])

________________________________________________________________________________________________
## mix fitting with 2 species model

model = cd3_tri_dissociation_bell!
u₀ = vcat(u₀_γϵ./3, u₀_δϵ./3, 1/3,0.0)
u₀ = vcat(u₀_γϵ./2, u₀_δϵ./2, 0.0,0.0)
tspan = (0.0,maximum(mix.t))
mix_ft_M = M(model,u₀,tspan,[],[])
r1 = ks_γϵ
r2 = ks_δϵ
pinit = [0.6,0.1,
        0.1,0.5,
        0.1,0.5,
        10.0,5.0]
opt =  [(0.59,0.61),(0.,10.),
        (0.,1.),(0.,1.),
        (0.,10.),(0.,1.),
        (0.,100.),(0.,10.)]
@time rates_mix, l_mix = optimization(mix_ft_M,mix,opt)

ks_mix = rates_mix
ks_mix = pinit
ts = 0:dt:maximum(mix.t)
Fs = 0:1:35
h2 = plot_ft_fit(mix_ft_M,ks_mix,mix)

histogram(γϵ.F,bins=10)

h = fit(Histogram, rand(100))
plot(h)
U = kde(γϵ.F;kernel=Distributions.Normal)
ik = InterpKDE(U)

pdf(U, [1,2,3])
plot(U.x,U.density)
scatter!(γϵ.F, zeros(size(γϵ.F)), markersize=7,markeralpha=0.3)

size(γϵ.F)

X = randn(100)
Y = randn(100)

k = kde(X)

pdf(k, k.x) ≈ k.density

ik = InterpKDE(k)

pdf(ik,2)