text:
function apply_branching!(node::Node, bu::BoundUpdate{LT})
    node.depth += 1
    node.bound_update = bu
    push!(node.lt_bounds, bu)
    return nothing
end

function apply_branching!(node::Node, bu::BoundUpdate{GT})
    node.depth += 1
    node.bound_update = bu
    push!(node.gt_bounds, bu)
    return nothing
end

function apply_branching!(node::Node, ac::AffineConstraint{LT})
    node.depth += 1
    # NOTE: Can potentially modify bd in-place; would need to note this in the
    # contract, though.
    f = CSAF(ac.f.coeffs, ac.f.indices, 0.0)
    s = LT(ac.s.upper - ac.f.constant)
    push!(node.lt_general_constrs, AffineConstraint(f, s))
    return nothing
end

function apply_branching!(node::Node, ac::AffineConstraint{GT})
    node.depth += 1
    # NOTE: Can potentially modify bd in-place; would need to note this in the
    # contract, though.
    f = CSAF(ac.f.coeffs, ac.f.indices, 0.0)
    s = GT(ac.s.lower - ac.f.constant)
    push!(node.gt_general_constrs, AffineConstraint(f, s))
    return nothing
end

function down_branch(node::Node, branch_cvi::CVI, val::Float64)
    return _branch(node, branch_cvi, LT(floor(Int, val)))
end

function up_branch(node::Node, branch_cvi::CVI, val::Float64)
    return _branch(node, branch_cvi, GT(ceil(Int, val)))
end

function _branch(node::Node, branch_cvi::CVI, set::S) where {S<:Union{LT,GT}}
    # TODO: Can likely reuse this memory instead of copying
    new_node = copy(node)
    apply_branching!(new_node, BoundUpdate(branch_cvi, set))
    return new_node
end

"""
Abstract type representing a possible way the algorithm could branch. Allows us
to understand the value in taking this particular branching choice without
instantiating all of the data structures necessary to actually represent it.
"""
abstract type AbstractBranchingCandidate end

"""
Represent a variable branching candidate in terms of: 1) the variable being
branched on, and 2) its fractional value at the LP solution.
"""
struct VariableBranchingCandidate <: AbstractBranchingCandidate
    index::CVI
    value::Float64
end

"""
    branching_candidates(form::DMIPFormulation, nr::NodeResult, config::AlgorithmConfig{B}) where {B <: AbstractBranchingRule}

Given the formulation, the result from the node LP solve, and the branching
rule (of type `B` as encoded in `config`), returns a vector of potential
branching candidates.

For example, if `B <: AbstractVariableBranchingRule`, there is a method which
will inspect the current LP solution and return branching candidates for each
integer variable taking a fractional value.
"""
function branching_candidates end

function branching_candidates(
    form::DMIPFormulation,
    parent_result::NodeResult,
    config::AlgorithmConfig{B},
)::Vector{VariableBranchingCandidate} where {B<:AbstractVariableBranchingRule}
    candidates = VariableBranchingCandidate[]
    for cvi in all_variables(form)
        var_set = get_variable_kind(form, cvi)
        if typeof(var_set) <: Union{ZO,GI}
            xi = parent_result.x[index(cvi)]
            xi_f = _approx_floor(xi, config.int_tol)
            xi_c = _approx_ceil(xi, config.int_tol)
            frac_val = min(xi - xi_f, xi_c - xi)
            if frac_val > config.int_tol
                # Otherwise, we're integral up to tolerance.
                push!(candidates, VariableBranchingCandidate(cvi, xi))
            end
        end
    end
    return candidates
end

"""
Abstract type representing a numeric scoring for a branching candidate.
"""
abstract type AbstractBranchingScore end

"""
Generic representation for scoring a variable branching candidate in terms of
a score for the down branch, a score for the up branch, and an aggregate
for the two.
"""
struct VariableBranchingScore <: AbstractBranchingScore
    down_branch_score::Float64
    up_branch_score::Float64
    aggregate_score::Float64
end

"""
    aggregate_score(::AbstractBranchingScore)::Float64

Represent a scoring for a branching candidate using a scalar numeric value.
Conventionally, this will be a value between 0 and 1.
"""
function aggregate_score end

aggregate_score(vbs::VariableBranchingScore) = vbs.aggregate_score

"""
    branch_on(parent_node::Node, candidate::AbstractBranchingCandidate, score::AbstractBranchingScore)::NTuple{Node}

Given a node, a selected branching candidate, and a score for that candidate,
apply that branching. The return value will be a tuple of `Node`s. These nodes
will be proved in increasing order of preference, i.e. the _last_ entry in the
return tuple will be the favorite.

Note:
* The memory in `parent_node` should remain untouched unaliased.
"""
function branch_on end

function branch_on(
    parent_node::Node,
    candidate::VariableBranchingCandidate,
    score::VariableBranchingScore,
)::Tuple{Node,Node}
    db = down_branch(parent_node, candidate.index, candidate.value)
    ub = up_branch(parent_node, candidate.index, candidate.value)
    if score.down_branch_score > score.up_branch_score
        return ub, db
    else
        return db, ub
    end
end

"""
    branching_score(state::CurrentState, bc::AbstractBranchingCandidate, nr::NodeResult, config::AlgorithmConfig)::AbstractBranchingScore

Scores a given branching candidate that could potentially be taken at the
current node.

Note:
* The two methods are added to hide the
"""
function branching_score end

function branching_score(
    ::CurrentState,
    bc::VariableBranchingCandidate,
    ::Node,
    parent_result::NodeResult,
    config::AlgorithmConfig{MostInfeasible},
)
    xi = parent_result.x[index(bc.index)]
    xi_f = _approx_floor(xi, config.int_tol)
    xi_c = _approx_ceil(xi, config.int_tol)
    f⁺ = xi_c - xi
    f⁻ = xi - xi_f
    return VariableBranchingScore(f⁻, f⁺, min(f⁻, f⁺))
end

function branching_score(
    ::CurrentState,
    bc::VariableBranchingCandidate,
    parent_node::Node,
    parent_result::NodeResult,
    config::AlgorithmConfig{PseudocostBranching},
)
    up_hist = config.branching_rule.upward_pseudocost_hist
    down_hist = config.branching_rule.downward_pseudocost_hist

    ψ⁺ = bc.index in keys(up_hist) ? up_hist[bc.index].ψ : config.branching_rule.ψ⁺_average
    ψ⁻ = bc.index in keys(down_hist) ? down_hist[bc.index].ψ : config.branching_rule.ψ⁻_average

    xi = parent_result.x[index(bc.index)]
    xi_f = _approx_floor(xi, config.int_tol)
    xi_c = _approx_ceil(xi, config.int_tol)
    f⁺ = xi_c - xi
    f⁻ = xi - xi_f

    μ = config.branching_rule.μ
    ψ_aggregate = (1 - μ) * min(f⁺*ψ⁺, f⁻*ψ⁻) + μ * max(f⁺*ψ⁺, f⁻*ψ⁻)
    return VariableBranchingScore(f⁻*ψ⁻, f⁺*ψ⁺, ψ_aggregate)
end

# TODO: Make this function more efficient by caching sb_model somewhere (likely
# in config object, rather than as a function argument).
function branching_score(
    state::CurrentState,
    bc::VariableBranchingCandidate,
    ::Node,
    parent_result::NodeResult,
    config::AlgorithmConfig{StrongBranching},
)
    sb_model = config.lp_solver_factory(state, config)
    MOI.copy_to(sb_model, state.gurobi_model)

    c = parent_result.cost
    cvi = bc.index
    xi = parent_result.x[index(cvi)]

    low_bound = _approx_floor(xi, config.int_tol)
    up_bound = _approx_ceil(xi, config.int_tol)
    ci = state.constraint_state.base_state.var_constrs[index(cvi)]

    c_up = _branch_cost(sb_model, GT(up_bound), ci)
    c_down = _branch_cost(sb_model, LT(low_bound), ci)

    μ = config.branching_rule.μ
    f⁻ = c_down - c
    f⁺ = c_up - c
    f = (1 - μ) * min(f⁺, f⁻) + μ * max(f⁺, f⁻)
    return VariableBranchingScore(f⁻, f⁺, f)
end

"""
    _branch_cost(model::Gurobi.Optimizer, constraint::Union{LT, GT}, ci::CI)

Computes the optimal objective cost of a branch.

The model is added with the new constraint, solved to optimal, and then
set back to the original.
"""
function _branch_cost(model::Gurobi.Optimizer, constraint::Union{LT,GT}, ci::CI)
    interval = MOI.get(model, MOI.ConstraintSet(), ci)
    temp_interval = IN(
        if constraint isa LT
            interval.lower
        else
            max(constraint.lower, interval.lower)
        end,
        if constraint isa GT
            interval.upper
        else
            min(constraint.upper, interval.upper)
        end,
    )

    MOI.set(model, MOI.ConstraintSet(), ci, temp_interval)
    MOI.optimize!(model)
    term_status = MOI.get(model, MOI.TerminationStatus())

    # NOTE: If the parent node LP was unbounded, we should never enter this
    # function. Therefore, we should never see an unbounded LP result.
    cost = (
        if term_status == MOI.OPTIMAL
            MOI.get(model, MOI.ObjectiveValue())
        elseif term_status == MOI.INFEASIBLE ||
               term_status == MOI.INFEASIBLE_OR_UNBOUNDED
            Inf
        else
            error("Unexpected termination status $term_status at node LP
                   when performing strong branching.")
        end
    )
    MOI.set(model, MOI.ConstraintSet(), ci, interval)
    return cost
end

"""
Given a current node and node LP result, branch. The return value will be a
tuple of `Node`s corresponding to the childen created by the branching.

This is a generic implementation which should support multiway branching and
general constraint branching. However, the "easy" case with binary variable
branching should be type stable and efficient.

Note:
* This method will error if branching is not valid. This is encoded via
`branching_candidates`: if no candidates are returned, then this function will
error.
"""
function branch(
    state::CurrentState,
    form::DMIPFormulation,
    parent_node::Node,
    parent_result::NodeResult,
    config::AlgorithmConfig,
)::NTuple
    # NOTE: The following calls dispatch on `config`. In the "simple" case
    # (i.e. variable branching), this should be type stable. In the more
    # complex case (general branching), it may not be.
    candidates = branching_candidates(form, parent_result, config)
    if isempty(candidates)
        error("No branching candidates--are you sure you want to be branching?")
    end
    scores = Dict(
        candidate =>
            branching_score(state, candidate, parent_node, parent_result, config) for
        candidate in candidates
    )
    # TODO: Can make all of this more efficient, if bottleneck.
    agg_scores = Dict(
        candidate => aggregate_score(scores[candidate]) for
        candidate in candidates
    )
    best_agg_score, best_candidate = findmax(agg_scores)
    return branch_on(parent_node, best_candidate, scores[best_candidate])
end


text:
# # [Select](@id man-select)

using Julio: Julio, Events
using Test

# `Julio.select` can be used for executing exactly one ready *event* from a set
# of events.

# ## Selecting a queue

function test_simple_select()
    #=
    To demonstrate how `Julio.select` works, suppose that we have multiple
    queues and waiting for an element from them.
    =#
    qin1, qout1 = Julio.queue()  # input/output endpoints for the first queue
    qin2, qout2 = Julio.queue()  # input/output endpoints for the second queue
    #=
    Suppose there is an element in the second queue:
    =#
    put!(qin2, 222)
    #=
    Then waiting on `take!` event of an empty queue and a nonempty queue
    does not block. It executes the `take!` event on the nonempty queue:
    =#
    selected = nothing
    Julio.select(
        (take!, qout1) => item -> begin
            selected = item
        end,
        (take!, qout2) => item -> begin
            selected = item
        end,
    )
    @test selected == 222
    #=
    ---
    =#
end

# ## Selecting an arbitrary event

# `Julio.select` can also be used with various events.

function test_mixed_select()
    #=
    It supports, for example, unbuffered channel:
    =#
    send_endpoint, _ = Julio.channel()
    #=
    ...and read/write on `IO` objects such as a pipe:
    =#
    Julio.open(`echo "hello"`) do output
        #=
        ...and acquiring locks:
        =#
        lck = Julio.Lock()
        # t = @task lock(lck)   #src
        # yield(t)              #src

        #=
        When the `Base` API `$f` has keyword arguments, you can use
        `Julio.Events.$f` to create an event.  Note that `($f, args...)` is
        equivalent to `Julio.Events.$f(args...)`.
        =#
        selected = nothing
        Julio.select(
            Events.readline(output; keep = true) => item -> begin
                selected = item
            end,
            Events.lock(lck) => _ -> begin
                unlock(lck)
                selected = :lock
            end,
            (put!, send_endpoint, 1) => _ -> begin
                selected = :put_1
            end,
        )
        #=
        In the above example, since the `output` pipe and the `lock` are
        both ready, `Julio.select` can select any one of them. However,
        since there is no other task `take!`ing the element from the channel
        `cho`, the `put!` event can not be selected.
        =#
        @test selected in ("hello\n", :lock)
        #=
        If `Events.lock(lck)` was selected, the `output` is not consumed:
        =#
        if selected === :lock
            @test readline(output) == "hello"
        end
        #=
        ---
        =#
    end
end

# ## [Example: bounding search results](@id ex-bounding-search)
#
# Suppose we need to move at least `minitems` items from one channel to another
# while filtering them using a predicate function `f`. Furthermore, we don't
# want to loose any items. That is to say, once an item is taken from the input
# channel, it must be put into the output channel (unless `f` evaluates to
# false).  Note that the we cannot use the [cancel scope](@ref man-cancel-scope)
# (naively) due to the last requirement; i.e., it is not correct to cancel the
# task when it's blocked while putting the item to the output channel. While we
# can still use the cancel scope by surrounding post-`take!` code in a
# `Julio.shield` block, the folowing code demonstrates more straightforward
# approach based on `Julio.select`.
#
# To setup cancellation specific to one event (`take!`), we can use explicit
# "cancellation token" and combine it with the original event.

using Julio: maybetake!, tryput!

function channel_filter!(f, output, input, minitems; ntasks = 4 * Threads.nthreads())
    nitems = Threads.Atomic{Int}(0)
    done = Julio.Promise{Nothing}()  # cancellation token
    Julio.withtaskgroup() do tg
        for _ in 1:ntasks
            Julio.spawn!(tg) do
                while true
                    m = Julio.select(
                        (fetch, done) => Returns(nothing),  # return nothing when done
                        (maybetake!, input),  # return Some(x) if we took x
                    )
                    x = @something(m, break)  # break if done
                    if f(x)
                        put!(output, x)
                        # If enough items have been sent, signal other tasks to finish.
                        if Threads.atomic_add!(nitems, 1) + 1 >= minitems
                            tryput!(done, nothing)
                            break
                        end
                    end
                end
            end
        end
    end
end

# In the above example, we use `(Julio.maybetake!, input)` event instead of
# `(take!, input)` event.  This is for handling the case `input` is closed.
# That is to say, `m === nothing` if `input` is closed or `tryput!(done,
# nothing)` has been executed.

function test_channel_filter()
    Julio.withtaskgroup() do tg
        send_endpoint1, receive_endpoint1 = Julio.channel()
        send_endpoint2, receive_endpoint2 = Julio.channel()
        Julio.spawn!(tg) do
            try
                for i in 1:15
                    put!(send_endpoint1, i)
                end
            finally
                close(send_endpoint1)
            end
        end
        Julio.spawn!(tg) do
            try
                channel_filter!(isodd, send_endpoint2, receive_endpoint1, 3; ntasks = 2)
            finally
                close(send_endpoint2)
            end
        end
        try
            out2 = collect(receive_endpoint2)
            out1 = collect(receive_endpoint1)
            sort!(out2)
            sort!(out1)
            @test length(out2) >= 3  # `channel_filter!` produced at least 3 elements
            @test out2 == (1:length(out2)) .* 2 .- 1
            @test out1 == out1[1]:out1[end]
        finally
            close(receive_endpoint1)
            close(receive_endpoint2)
        end
    end
end

ext: jl


text:
using Enzyme
mutable struct Model
    F::Vector{Float64}
    F_H::Vector{Float64}
    t::Float64
    h::Float64
end

function func_U(t)
    e = exp(1)
    return 2.0*((e^(3.0*t))-(e^3))/((e^(3.0*t/2.0))*(2.0+(e^3)))
end

function func(F, X,t)
    F[2] = X[1]*X[1]+0.5*(func_U(t)*func_U(t))
    F[1] = 0.5*X[1]+ func_U(t)
    return nothing
end

function advance(model)
    F_H = model.F_H
    F = model.F
    t = model.t
    h = model.h
    func(F, F_H,t)
    F[1] = F_H[1] + h/2.0*F[1]
    F[2] = F_H[2] + h/2.0*F[2]
    func(F,F,t+h/2.0)
    model.F[1] = F_H[1] + h*F[1]
    model.F[2] = F_H[2] + h*F[2]
    return nothing
end

function opt_sol(Y,t)
    e = exp(1)
    Y[1] = (2.0*e^(3.0*t)+e^3)/(e^(3.0*t/2.0)*(2.0+e^3))
    Y[2] = (2.0*e^(3.0*t)-e^(6.0-3.0*t)-2.0+e^6)/((2.0+e^3)^2)
    return
end

function opt_lambda(L,t)
    e = exp(1)
    L[1] = (2.0*e^(3-t)-2.0*e^(2.0*t))/(e^(t/2.0)*(2+e^3))
    L[2] = 1.0
    return
end


text:
"""
`loc_dev(g) -> deviance from local indices`

Evaluate the sum of the squared deviance residuals on the local indices of `y`
"""
function loc_dev!{T<:AbstractFloat}(g::PGLM{T})
    dev = zero(T)
    @inbounds for j in localindexes(g.y)
        sm = zero(T)
        @simd for k in 1:length(g.βs)
            sm += g.Xt[k,j]*g.βs[k]
        end
        g.η[j] = sm
        g.μ[j] = invlink(g.l,sm)
        dev += g.wt[j] * devresid2(g.d,g.y[j],g.μ[j])
    end
    g.dev[g.y.pidx] = dev
end

Base.size(g::PGLM) = ((p,n) = size(g.Xt); (n,p,length(g.dev)))

"""
Evaluate the deviance for `PGLM` model `g` using step factor `s`

This method also updates the `η` and `μ` members of `g`
"""
function StatsBase.deviance{T<:AbstractFloat}(g::PGLM{T},s::T)
    @simd for k in 1:length(g.β)
        @inbounds g.βs[k] = g.β[k] + s * g.δβ[k]
    end
    fill!(g.dev,zero(T))
    pmap(loc_dev!,fill(g,nworkers()))
    sum(g.dev)
end
StatsBase.deviance{T<:AbstractFloat}(g::PGLM{T}) = deviance(g,zero(T))

"""
`loc_updateXtW(g) -> nothing`

local function to update the `XtW` array
"""
function loc_updateXtW!{T<:AbstractFloat}(g::PGLM{T})
    n,p,npr = size(g)
    k = procs(g.y)[g.y.pidx]
    if g.l == canonical(g.d)
        @inbounds for ii in localindexes(g.y)
            W = g.wt[ii] * varfunc(g.d,g.μ[ii])
            for j in 1:p
                Wj = g.Xt[j,ii]
                g.XtWr[j,k] += Wj * (g.y[ii] - g.μ[ii])
                Wj *= W
                @simd for i in j:p
                    g.XtWX[i,j,k] += g.Xt[i,ii] * Wj
                end
            end
        end
    else
        for ii in localindexes(g.y)
            mueta = μη(g.l,g.η[ii])
            W = g.wt[ii] * abs2(mueta)/max(eps(T),varfunc(g.d,g.μ[ii]))
            for j in 1:p
                @simd for i in j:p
                    g.XtWX[i,j,k] += g.Xt[i,ii] * W * g.Xt[j,ii]
                end
                g.XtWr[j,k] += g.Xt[j,ii]*W*(g.y[ii] - g.μ[ii])/max(eps(T),mueta)
            end
        end
    end
end

"""
`updateXtW(g) -> δβ`

update the `XtW`, `XtWr` and `δβ` arrays in `g`
"""
function updateXtW!{T<:AbstractFloat}(g::PGLM{T})
    n,p,npr = size(g)
    fill!(g.XtWX,zero(T))
    fill!(g.XtWr,zero(T))
    pmap(loc_updateXtW!, fill(g,nworkers()))
    for k in procs(g.y)
        for j in 1:p
            for i in j:p
                g.XtWX[i,j,1] += g.XtWX[i,j,k]
            end
            g.δβ[j] += g.XtWr[j,k]
        end
    end
    A_ldiv_B!(cholfact!(sub(sdata(g.XtWX),:,:,1),:L),g.δβ)
end

function StatsBase.fit(g::GLM;verbose::Bool=false, maxIter::Integer=30,
                       minStepFac::Real=0.001, convTol::Real=1.e-6)
    g.fit && return g
    maxIter > zero(maxIter) || error("maxIter must be positive")
    zero(minStepFac) < minStepFac < one(minStepFac) || error("minStepFac must be in (0,1)")
    T = eltype(g.y)

    cvg = false
    devold = deviance(g,zero(T))
    for i in 1:maxIter
        updateXtW!(g)
        s = one(T)                      # step factor
        dev = deviance(g,s)
        while dev > devold
            s *= convert(T,0.5)         # halve the step factor
            s > minStepFac || error("step-halving failed at β₀ = $(g.β)")
            dev = deviance(g,s)
        end
        copy!(g.β,g.βs)
        crit = (devold - dev)/dev
        verbose && println("$i: $dev, $crit")
        if crit < convTol
            cvg = true
            break
        end
        devold = dev
    end
    cvg || error("failure to converge in $maxIter iterations")
    g.fit = true
    g
end

function loc_initμη!(g::PGLM)
    for i in localindexes(g.y)
        g.μ[i] = mustart(g.d,g.y[i],g.wt[i])
        g.η[i] = link(g.l,g.μ[i])
    end
end

"""
`initμη(g) -> g`

Initialize the `μ` and `η` members of `g`
"""
initμη!(g::PGLM) = (pmap(loc_initμη!,fill(g,nworkers())), g)

function initμη!(g::SGLM)
    for i in 1:length(g.y)
        g.μ[i] = mustart(g.d,g.y[i],g.wt[i])
        g.η[i] = link(g.l,g.μ[i])
    end
    g
end

function updateXtW!{T<:AbstractFloat}(g::SGLM{T})
    p,n = size(g.Xt)
    fill!(g.XtWX,zero(T))
    fill!(g.δβ,zero(T))
    if g.l == canonical(g.d)
        @inbounds for ii in 1:length(g.y)
            W = sqrt(g.wt[ii] * varfunc(g.d,g.μ[ii]))
            g.wtres[ii] = g.y[ii] - g.μ[ii]
            @simd for j in 1:p
                g.XtW[j,ii] = W * g.Xt[j,ii]
            end
        end
        return A_ldiv_B!(cholfact(BLAS.syrk!('L','N',one(T),g.XtW,zero(T),g.XtWX),:L),
                         BLAS.gemv!('N',one(T),g.Xt,g.wtres,zero(T),g.δβ))
    end
    for ii in 1:length(g.y)
        mueta = μη(g.l,g.η[ii])
        W = g.wt[ii] * abs2(mueta)/max(eps(T),varfunc(g.d,g.μ[ii]))
        for j in 1:p
            @simd for i in j:p
                g.XtWX[i,j] += g.Xt[i,ii] * W * g.Xt[j,ii]
            end
            g.δβ[j] += g.Xt[j,ii]*W*(g.y[ii] - g.μ[ii])/max(eps(T),mueta)
        end
    end
    A_ldiv_B!(cholfact!(g.XtWX,:L),g.δβ)
end

"""
Evaluate the deviance for an `SGLM` model `g` using step factor `s`

This method also updates the `η` and `μ` members of `g`
"""
function StatsBase.deviance{T<:AbstractFloat}(g::SGLM{T},s::T)
    ## g.βs := g.β + s * g.δβ; g.η := g.Xt'*g.βs
    BLAS.gemv!('T',one(T),g.Xt,BLAS.axpy!(s,g.δβ,copy!(g.βs,g.β)),zero(T),g.η)
    dev = zero(T)
    @inbounds for j in 1:length(g.y)
        g.μ[j] = invlink(g.l,g.η[j])
        dev += g.wt[j] * devresid2(g.d,g.y[j],g.μ[j])
    end
    dev
end


text:
using Flux, Adapt

to_device(use_cuda::Val{:cpu}, x) = to_host(x)
to_device(use_cuda::Val{:gpu}, x) = to_gpu(x)

to_host(m) = fmap(x -> adapt(Array, x), m)
to_host(x::Array) = x

to_gpu(x) = fmap(CuArrays.cu, x)
# gpu_if_avail(x, use_cuda::Val{true}) = fmap(CuArrays.cu, x)


text:
# This file was generated, do not modify it. # hide
linear = @load LinearRegressor pkg=MLJLinearModels
ridge = @load RidgeRegressor pkg=MultivariateStats; ridge.lambda = 0.01
knn = @load KNNRegressor; knn.K = 4
tree = @load DecisionTreeRegressor; min_samples_leaf=1
forest = @load RandomForestRegressor pkg=DecisionTree; forest.n_trees=500
svm = @load SVMRegressor;

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
module SeisPGVPGA

# module to be used over modules
using  SeisIO, FileIO, Dates, DSP, Statistics, DataFrames, CSV, Trapz, Distributed

export

    map_seisPGVPGA,
    map_readPGVPGA,
    get_requeststations,
    get_starttimelist

# NOTE: Do not change the order of include.
include("func_SeisPGVPGA.jl")
include("map_SeisPGVPGA.jl")

end


text:
using Displaz

function example(N)
    t = range(0, stop=1, length=N)'
    P = (10 .+ 0.2*rand(Float64, size(t))) .* vcat(exp.(5*t).*cos.(100*pi*t),
                                                   exp.(5*t).*sin.(100*pi*t),
                                                   100*t);
    C = vcat(t, 1 .- t, zero(t));
    sz  = map(Float32, t)
    shape = mod.(rand(UInt8, N), 6)
    plot3d!(P, color = C, label="Top series", markershape=shape)
    P[3,:] = -P[3,:]
    plot3d!(P, color = C, markersize = sz, markershape="-", label="Bottom series")
end

example(1_000_000)


text:
abstract type AbstractImageManip end
function call(im::AbstractImageManip, img) end

struct NoManip <: AbstractImageManip end
(im::NoManip)(img) = img

mutable struct ImageBuffer{CB<:CircularBuffer, IM<:AbstractImageManip}
    buffer::CB
    image_manip::IM
    curr_step::Int64
end

ImageBuffer(size::Int64, im, T::Type) =
    ImageBuffer(
        CircularBuffer(size, (Int64, T), (:step, :image)),
        im, 0)

# ImageBuffer(size) = ImageBuffer(size, NoManip())

function add!(imb::IMB, img) where {IMB<:ImageBuffer}
    imb.curr_step += 1
    add!(imb.buffer, (imb.curr_step, imb.image_manip(img)))
end

@forward ImageBuffer.buffer getindex


text:
using PyPlot
using OffsetArrays

export showImage, plotRadius, plotEllipse

function showImage(image::OffsetArray{Float64,2,Array{Float64,2}})
        clf()
        imshow(image, aspect = "equal", extent = (minimum(axes(image)[1]), maximum(axes(image)[1]), minimum(axes(image)[2]), maximum(axes(image)[2])))
        colorbar()
        gcf()
end

function showImage(image::Array{Float64,2})
        clf()
        imshow(image, aspect = "equal", extent = (minimum(axes(image)[1]), maximum(axes(image)[1]), minimum(axes(image)[2]), maximum(axes(image)[2])))
        colorbar()
        gcf()
end

function showImage(x::Array{Float64,2}, y::Array{Float64,2}, z::Array{Float64,2})
        clf()
        pcolormesh(x, y, z)
        xlim([minimum(x), maximum(x)]);
        ylim([minimum(y), maximum(y)]);
        gca().set_aspect("equal")
        colorbar()
        gcf()
end

function plotRadius(radius::Prec64)
        plt.gca().add_artist(plt.Circle((0.0, 0.0), radius = radius, color = "r", fill=false))
        gcf()
end

function plotEllipse(center::Tuple{Float64,Float64}, width::Float64, height::Float64, angle::Float64)
        ax = gca()
        ell = matplotlib.patches.Ellipse(center, width, height, angle, color = "r", fill = false)
        ax.add_patch(ell)
        gcf()
end

function plotEllipse(ell::ParametricFormEllipse)
        ax = gca()
        center = tuple(ell.center...)
        width = ell.axis_lengths[1]
        height = ell.axis_lengths[2]
        angle = ell.ccw_angle
        ell = matplotlib.patches.Ellipse(center, width, height, angle, color = "r", fill = false)
        ax.add_patch(ell)
        gcf()
end

function plotEllipse(ell::ConicFormEllipse)
        plotEllipse(conic2parametric(ell))
end


text:
# Optimizes the code for computing the derivatives
# of I_v and J_v with respect to k.
include("../src/IJv_derivative.jl")

function optimize_IJv_derivative(k2::Float64)
v_max = 20
Iv_lower = zeros(Float64,v_max+1); Jv_lower = zeros(Float64,v_max+1)
Iv_raise = zeros(Float64,v_max+1); Jv_raise = zeros(Float64,v_max+1)
dIvdk_lower = zeros(Float64,v_max+1); dJvdk_lower = zeros(Float64,v_max+1)
dIvdk_raise = zeros(Float64,v_max+1); dJvdk_raise = zeros(Float64,v_max+1)
Iv_big_lower = zeros(BigFloat,v_max+1); Jv_big_lower = zeros(BigFloat,v_max+1)
Iv_big_raise = zeros(BigFloat,v_max+1); Jv_big_raise = zeros(BigFloat,v_max+1)
dIvdk_big_lower = zeros(BigFloat,v_max+1); dJvdk_big_lower = zeros(BigFloat,v_max+1)
dIvdk_big_raise = zeros(BigFloat,v_max+1); dJvdk_big_raise = zeros(BigFloat,v_max+1)
# This computes I_v for the largest v, and then works down to smaller values:
if k2 <= 1
  kc = convert(Float64,sqrt(1.-big(k2)))
  k = convert(Float64,sqrt(big(k2)))
  kck = convert(Float64,sqrt((1.-big(k2))*big(k2)))
  kap = convert(Float64,2*asin(sqrt(big(k2))))
  k2_big = big(k2)
  k_big = sqrt(k2_big)
  kc_big = sqrt(1-big(k2))
  kck_big = k_big*kc_big
  kap_big = 2*asin(k_big)
else
  kc = convert(Float64,sqrt(1.-inv(big(k2))))
  k = convert(Float64,sqrt(big(k2)))
  kap = convert(typeof(k2),pi); kck=zero(k2)
  k2_big = big(k2)
  k_big = sqrt(k2_big)
  kc_big = sqrt(big(1.0)-inv(big(k2)))
  kck_big = big(0.0)
  kap_big = big(pi)
end
if k2 > 0
  dIJv_lower_dk!(v_max,k2,kck,kc,kap,Iv_lower,Jv_lower,dIvdk_lower,dJvdk_lower)
  dIJv_lower_dk!(v_max,k2_big,kck_big,kc_big,kap_big,Iv_big_lower,Jv_big_lower,dIvdk_big_lower,dJvdk_big_lower)
  dIJv_raise_dk!(v_max,k2,kck,kc,kap,Iv_raise,Jv_raise,dIvdk_raise,dJvdk_raise)
  dIJv_raise_dk!(v_max,k2_big,kck_big,kc_big,kap_big,Iv_big_raise,Jv_big_raise,dIvdk_big_raise,dJvdk_big_raise)
#  for v=0:v_max
#    println("v: ",v," k2: ",k2," kc: ",kc," Iv_lower: ",Iv_lower[v+1]," Iv_big_lower: ",convert(Float64,Iv_big_lower[v+1]),
#    " Iv_raise: ",Iv_raise[v+1]," Iv_big_raise: ",convert(Float64,Iv_big_raise[v+1]),
#    " Jv_lower: ",Jv_lower[v+1]," Jv_big_lower: ",convert(Float64,Jv_big_lower[v+1]),
#    " Jv_raise: ",Jv_raise[v+1]," Jv_big_raise: ",convert(Float64,Jv_big_raise[v+1]))
##      " dIvdk_lower: ",dIvdk[v+1]," dIvdk_num: ",dIvdk2_num[v+1]*2*k," diff: ",dIvdk[v+1]-dIvdk2_num[v+1]*2*k," dJvdk: ",dJvdk[v+1]," dJvdk_num: ",dJvdk2_num[v+1]*2*k," diff: ",dJvdk[v+1]-dJvdk2_num[v+1]*2*k)
#  end
end
return Iv_lower,Jv_lower,Iv_raise,Jv_raise,dIvdk_lower,dJvdk_lower,dIvdk_raise,dJvdk_raise,convert(Array{Float64,1},Iv_big_lower),convert(Array{Float64,1},Jv_big_lower),convert(Array{Float64,1},Iv_big_raise),convert(Array{Float64,1},Jv_big_raise),convert(Array{Float64,1},dIvdk_big_lower),convert(Array{Float64,1},dJvdk_big_lower),convert(Array{Float64,1},dIvdk_big_raise),convert(Array{Float64,1},dJvdk_big_raise)
end

nk2 = 100
v_max = 20
Iv_lower = zeros(Float64,nk2,v_max+1); Jv_lower = zeros(Float64,nk2,v_max+1)
Iv_raise = zeros(Float64,nk2,v_max+1); Jv_raise = zeros(Float64,nk2,v_max+1)
dIvdk_lower = zeros(Float64,nk2,v_max+1); dJvdk_lower = zeros(Float64,nk2,v_max+1)
dIvdk_raise = zeros(Float64,nk2,v_max+1); dJvdk_raise = zeros(Float64,nk2,v_max+1)
Iv_big_lower = zeros(Float64,nk2,v_max+1); Jv_big_lower = zeros(Float64,nk2,v_max+1)
Iv_big_raise = zeros(Float64,nk2,v_max+1); Jv_big_raise = zeros(Float64,nk2,v_max+1)
dIvdk_big_lower = zeros(Float64,nk2,v_max+1); dJvdk_big_lower = zeros(Float64,nk2,v_max+1)
dIvdk_big_raise = zeros(Float64,nk2,v_max+1); dJvdk_big_raise = zeros(Float64,nk2,v_max+1)
k2 = logspace(0,2,nk2)
for i=1:nk2
  Iv_lower[i,:],Jv_lower[i,:],Iv_raise[i,:],Jv_raise[i,:],dIvdk_lower[i,:],dJvdk_lower,dIvdk_raise[i,:],dJvdk_raise[i,:],Iv_big_lower[i,:],Jv_big_lower[i,:],Iv_big_raise[i,:],Jv_big_raise[i,:],dIvdk_big_lower[i,:],dJvdk_big_lower[i,:],dIvdk_big_raise[i,:],dJvdk_big_raise = optimize_IJv_derivative(k2[i])
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)

# D, Δx = 1.0, 0.35

Γ = readmesh(joinpath(@__DIR__,"sphere2.in"))
X, Y = raviartthomas(Γ), buffachristiansen(Γ)

# Δt, Nt = 0.11, 200
Δt, Nt = 0.6, 200
δ = timebasisdelta(Δt, Nt)
h = timebasisc0d1(Δt, Nt)
p = timebasiscxd0(Δt, Nt)
V = X ⊗ p
W = Y ⊗ p

width, delay, scaling = 8.0, 12.0, 1.0
gaussian = creategaussian(width, delay, scaling)

direction, polarisation = ẑ, x̂
E = BEAST.planewave(polarisation, direction, gaussian, 1.0)
# E = BEAST.planewave(polarisation, direction, derive(gaussian), 1.0)
H = direction × E

@hilbertspace j; @hilbertspace m′
K, I, N = MWDoubleLayerTDIO(1.0, 1.0, 0), Identity(), NCross()
tdmfie = @discretise (0.5*(N⊗I) + 1.0*K)[m′,j] == -1H[m′]   j∈V  m′∈W
xmfie = solve(tdmfie)

Xmfie, Δω, ω0 = fouriertransform(xmfie, Δt, 0.0, 2)
ω = collect(ω0 .+ (0:Nt-1)*Δω)
_, i1 = findmin(abs.(ω .- 1.0))

ω1 = ω[i1]
um = Xmfie[:,i1] / fouriertransform(gaussian)(ω1)


text:
"""
    increment_grade(grades, bonuses, exam)

One way of apply a bonus. We simply increment the grade by one level.
But only if the `grade` is `passed => true`.
"""
function increment_grade(grades, bonuses, exam)
    all_possible_grades = get_grades(exam)

    grades_bonus_eligible = []
    for g in exam["grades"]
        if g["passed"]
            push!(grades_bonus_eligible, g["grade"])
        end
    end

    new_grades = []
    for (grade, bonus) in zip(grades, bonuses)
        ind = findfirst(g -> g == grade, all_possible_grades)
        if ind > 1 && bonus && grade ∈ grades_bonus_eligible
            push!(new_grades, all_possible_grades[ind-1])
        else
            push!(new_grades, grade)
        end
    end

    return new_grades
end


text:
using Test, BenchmarkTools, BesselK, SpecialFunctions, FiniteDifferences, ForwardDiff

const VGRID   = range(0.25, 10.0, length=100)
const XGRID   = range(0.0, 50.0,  length=201)[2:end]
const VX      = collect(Iterators.product(VGRID, XGRID))
const REF_FD1 = central_fdm(10,1)
const REF_FD2 = central_fdm(10,2)

atolfun(tru, est) = isnan(est) ? NaN : (isinf(tru) ? 0.0 : abs(tru-est))

besselkxv(v,x) = besselk(v,x)*(x^v)

fd_dbesselk_dv(v, x) = REF_FD1(_v->besselk(_v, x), v)
fd2_dbesselk_dv_dv(v, x) = REF_FD2(_v->besselk(_v, x), v)

fd_dbesselkxv_dv(v, x) = REF_FD1(_v->besselkxv(_v, x), v)
fd2_dbesselkxv_dv_dv(v, x) = REF_FD2(_v->besselkxv(_v, x), v)

ad_dbesselk_dv(v, x) = ForwardDiff.derivative(_v->BesselK.adbesselk(_v, x), v)
ad2_dbesselk_dv_dv(v, x) = ForwardDiff.derivative(_v->ad_dbesselk_dv(_v, x), v)

ad_dbesselkxv_dv(v, x) = ForwardDiff.derivative(_v->BesselK.adbesselkxv(_v, x), v)
ad2_dbesselkxv_dv_dv(v, x) = ForwardDiff.derivative(_v->ad_dbesselkxv_dv(_v, x), v)

# direct accuracy:
@testset "direct eval" begin
  println("\nDirect evaluations:")
  for (ref_fn, cand_fn, case) in ((besselk,   BesselK._besselk, :standard),
                                  (besselkxv, BesselK.adbesselkxv, :rescaled))
    amos_ref  = map(vx->ref_fn(vx[1], vx[2]), VX)
    candidate = map(vx->cand_fn(vx[1], vx[2]), VX)
    atols     = map(a_c->atolfun(a_c[1], a_c[2]), zip(amos_ref, candidate))
    ix        = findall(x-> x <= 1000.0, amos_ref)
    thresh    = case == :standard ? 5e-11 : 2e-12
    (maxerr, maxix) = findmax(abs, atols[ix])
    (maxerr_v, maxerr_x) = VX[ix][maxix]
    println("Case $case:")
    println("worst (v,x):  ($maxerr_v, $maxerr_x)")
    println("Ref value:     $(amos_ref[ix][maxix])")
    println("Est value:     $(candidate[ix][maxix])")
    println("Abs error:     $(round(maxerr, sigdigits=3))")
    @test maxerr < thresh
  end
  println()
end

# test derivative accuracy:
@testset "first derivative" begin
  println("\nFirst derivatives:")
  for (ref_fn, cand_fn, case) in ((fd_dbesselk_dv, ad_dbesselk_dv, :standard),
                                  (fd_dbesselkxv_dv, ad_dbesselkxv_dv, :rescaled))
    amos_ref  = map(vx->ref_fn(vx[1], vx[2]), VX)
    candidate = map(vx->cand_fn(vx[1], vx[2]), VX)
    atols     = map(a_c->atolfun(a_c[1], a_c[2]), zip(amos_ref, candidate))
    ix        = findall(x-> x <= 1000.0, amos_ref)
    thresh    = case == :standard ? 4e-9 : 8e-7
    (maxerr, maxix) = findmax(abs, atols[ix])
    (maxerr_v, maxerr_x) = VX[ix][maxix]
    println("Case $case:")
    println("worst (v,x):  ($maxerr_v, $maxerr_x)")
    println("Ref value:     $(amos_ref[ix][maxix])")
    println("Est value:     $(candidate[ix][maxix])")
    println("Abs error:     $(round(maxerr, sigdigits=3))")
    @test maxerr < thresh
  end
  println()
end

# test second derivative accuracy:
@testset "second derivative" begin
  println("\nSecond derivatives:")
  for (ref_fn, cand_fn, case) in ((fd2_dbesselk_dv_dv, ad2_dbesselk_dv_dv, :standard),
                                  (fd2_dbesselkxv_dv_dv, ad2_dbesselkxv_dv_dv, :rescaled))
    amos_ref  = map(vx->ref_fn(vx[1], vx[2]), VX)
    candidate = map(vx->cand_fn(vx[1], vx[2]), VX)
    atols     = map(a_c->atolfun(a_c[1], a_c[2]), zip(amos_ref, candidate))
    ix        = findall(x-> x <= 100.0, amos_ref)
    thresh    = case == :standard ? 5e-7 : 2e-6
    (maxerr, maxix) = findmax(abs, atols[ix])
    (maxerr_v, maxerr_x) = VX[ix][maxix]
    println("Case $case:")
    println("worst (v,x):  ($maxerr_v, $maxerr_x)")
    println("Ref value:     $(amos_ref[ix][maxix])")
    println("Est value:     $(candidate[ix][maxix])")
    println("Abs error:     $(round(maxerr, sigdigits=3))")
    @test maxerr < thresh
  end
  println()
end

# Testing the _xv versions really slows down the test script, and in general
# there are no no routines.
@testset "confirm no allocations" begin
  VGRID_ALLOC = (0.25, 1.0-1e-8, 1.0, 1.5, 2.1, 3.0, 3.5, 4.8)
  XGRID_ALLOC = range(0.0, 50.0, length=11)[2:end]
  VX_ALLOC    = collect(Iterators.product(VGRID_ALLOC, XGRID_ALLOC))
  ad_alloc_test(v,x)     = @ballocated ad_dbesselk_dv($v,$x) samples=5
  ad2_alloc_test(v,x)    = @ballocated ad2_dbesselk_dv_dv($v,$x) samples=5
  ad_alloc_test_xv(v,x)  = @ballocated ad_dbesselkxv_dv($v,$x) samples=5
  ad2_alloc_test_xv(v,x) = @ballocated ad2_dbesselkxv_dv_dv($v,$x) samples=5
  ad_allocs  = map(vx->ad_alloc_test(vx[1], vx[2]), VX_ALLOC)
  ad2_allocs = map(vx->ad2_alloc_test(vx[1], vx[2]), VX_ALLOC)
  ad_allocs_xv  = map(vx->ad_alloc_test_xv(vx[1], vx[2]), VX_ALLOC)
  ad2_allocs_xv = map(vx->ad2_alloc_test_xv(vx[1], vx[2]), VX_ALLOC)
  @test all(iszero, ad_allocs)
  @test all(iszero, ad2_allocs)
  @test all(iszero, ad_allocs_xv)
  @test all(iszero, ad2_allocs_xv)
end


text:
abstract type AbstractCoarsening end

abstract type AbstractAlgebraicCoarsening <: AbstractCoarsening end

abstract type AbstractStrength end


text:
#examples of using segy2jld2 functions to convert PASSCAL segy files to jld2 format.
using SeisIO, SeisConvert
# Saves individual segy file (in) to JLD2 file.
# The output file name will be the same as the input, with new extension *.jld2.

infile = "02.089.00.01.50.7460.1"
segy2jld2(infile)
#

# """
# Saves individual segy file (in) to JLD2 file (out), which is explicitly specified.
# """

segy2jld2(infile,"02.089.00.01.50.7460.1_younameit.jld2")

# """
# Saves segy files based on a filelist generated using ls() in Julia to
# JLD2 file (out), which is explicitly specified.
# """

filelist = ls("segy/02.089.00.01.50/02*")

segy2jld2(filelist,"2002.089T00.01.50","2002.089T00.01.50.jld2")

# """
# Searches for all segy files scanning throgh the segydirlist (Array{String,1}). The segy files
#     in each segy directory are assigned to one time stamp based on timestamplist (Array{String,1}).
#     The two Arrays must be equal in size. The list can be generated from a file list (using ls in shell),
#     using a utility function in this package: readfilelist(infilelist::String).
#
# All segy files will be saved into one single JLD2 file (out), which is explicitly specified as outfile.
# """
segydirlist = ls("segy/02.089*")
timestamplist=["2002.089T00.01.50","2002.089T00.02.10"]
segy2jld2(segydirlist,timestamplist,"testoutput_2groups.jld2")

"""
examples when we need to update site information in the segy files.
compile_stationinfo() is a function in utilities.jl
"""
siteinfodict=compile_stationinfo("station_list.txt",separator=',',header=true)
segy2jld2(infile,"02.089.00.01.50.7460.1_updatedsiteinfo.jld2",siteinfodict=siteinfodict)
#other methods have the similar usage with an additional argument: siteinfodict=your_siteinfodict


text:
@testset "sparse_constructor" begin
    A = sprand(10,10,0.1)
    s = SparseTensor(A)
    @test run(sess, s)≈A
    I = [1;2;4;2;3;5]
    J = [1;3;2;2;2;1]
    V = rand(6)
    A = sparse(I,J,V,6,5)
    s = SparseTensor(I,J,V,6,5)
    @test run(sess, s)≈A
    indices = [I J]
    s = SparseTensor(I,J,V,6,5)
    @test run(sess, s)≈A

    S = Array(s)
    @test run(sess, S)≈Array(A)

    @test size(s)==(6,5)
    @test size(s,1)==6
    @test size(s,2)==5
end

@testset "sparse_arithmetic" begin
    A1 = rand(10,5); B1 = sprand(10,5,0.3)
    A = constant(A1)
    B = SparseTensor(B1)
    @test run(sess, -B)≈-B1
    for op in [+,-]
        C = op(A, B)
        C1 = op(A1, B1)
        @test run(sess, C)≈C1
    end
    @test run(sess, B-A)≈B1-A1
end

@testset "sparse_adjoint" begin
    A = sprand(10,5,0.3)
    A1 = SparseTensor(A)
    @test run(sess, A1')≈sparse(A')
end

@testset "sparse_mul" begin
    A1 = rand(10,10); B1 = sprand(10,10,0.3)
    C1 = rand(10)
    A = constant(A1)
    B = SparseTensor(B1)
    C = constant(C1)
    @test run(sess, B*A) ≈ B1*A1
    @test run(sess, A*B) ≈ A1*B1
    @test run(sess, B*A1) ≈ B1*A1
    @test run(sess, A1*B) ≈ A1*B1
    @test run(sess, B*C) ≈ B1*C1
    @test run(sess, B*C1) ≈ B1*C1
end

@testset "sparse_vcat_hcat" begin
    B1 = sprand(10,3,0.3)
    B = SparseTensor(B1)
    @test run(sess, [B;B])≈[B1;B1]
    @test run(sess, [B B])≈[B1 B1]
end

@testset "sparse_indexing" begin
    B1 = sprand(10,10,0.3)
    B = SparseTensor(B1)
    @test run(sess, B[2:3,2:3])≈B1[2:3,2:3]
    @test run(sess, B[2:3,:])≈B1[2:3,:]
    @test run(sess, B[:,2:3])≈B1[:,2:3]
end

@testset "sparse_solve" begin
    A = sparse(I, 10,10) + sprand(10,10,0.1)
    b = rand(10)
    A1 = SparseTensor(A)
    b1 = constant(b)
    u = A1\b1
    @test run(sess, u) ≈ A\b
end

@testset "sparse_assembler" begin

    m = 20
    n = 100
    handle = SparseAssembler(100, m, 0.0)
    op = PyObject[]
    A = zeros(m, n)
    for i = 1:1
        ncol = rand(1:n, 10)
        row = rand(1:m)
        v = rand(10)
        for (k,val) in enumerate(v)
            @show k
            A[row, ncol[k]] += val
        end
        @show v
        push!(op, accumulate(handle, row, ncol, v))
    end
    op = vcat(op...)
    J = assemble(m, n, op)
    B = run(sess, J)
    @test norm(A-B)<1e-8

    handle = SparseAssembler(100, 5, 1.0)
    op1 = accumulate(handle, 1, [1;2;3], [2.0;0.5;0.5])
    J = assemble(5, 5, op1)
    B = run(sess, J)
    @test norm(B-[2.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0])<1e-8

    handle = SparseAssembler(100, 5, 0.0)
    op1 = accumulate(handle, 1, [1;1], [1.0;1.0])
    op2 = accumulate(handle, 1, [1;2], [1.0;1.0])

    J = assemble(5, 5, [op1;op2])
    B = run(sess, J)
    @test norm(B-[3.0  1.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0
        0.0  0.0  0.0  0.0  0.0])<1e-8

end

@testset "sparse_least_square" begin
    ii = Int32[1;1;2;2;3;3]
    jj = Int32[1;2;1;2;1;2]
    vv = Float64[1;2;3;4;5;6]
    ff = Float64[1;1;1]
    A = SparseTensor(ii, jj, vv, 3, 2)
    o = A\ff
    @test norm(run(sess, o)-[-1;1])<1e-6
end

@testset "sparse mat mul" begin
    A = sprand(10,5,0.3)
    B = sprand(5,20,0.3)
    C = A*B
    CC = SparseTensor(A)*SparseTensor(B)
    C_ = run(sess, CC)
    @test C_≈C

    A = spdiagm(0=>[1.;2.;3;4;5])
    B = sprand(5,20,0.3)
    C = A*B
    CC = SparseTensor(A)*SparseTensor(B)
    C_ = run(sess, CC)
    @test C_≈C

    A = sprand(10,5,0.5)
    B = spdiagm(0=>[1.;2.;3;4;5])
    C = A*B
    CC = SparseTensor(A)*SparseTensor(B)
    C_ = run(sess, CC)
    @test C_≈C
end

@testset "spdiag" begin
    p = rand(10)
    A = spdiagm(0=>p)
    B = spdiag(constant(p))
    C = spdiag(10)
    @test run(sess, B)≈A
    @test B._diag
    @test run(sess, C)≈spdiagm(0=>ones(10))
end

@testset "spzero" begin
    q = spzero(10)
    @test run(sess, q)≈sparse(zeros(10,10))
    q = spzero(10,20)
    @test run(sess, q)≈sparse(zeros(10,20))
end

@testset "sparse indexing" begin
    i1 = unique(rand(1:20,3))
    j1 = unique(rand(1:30,3))
    A = sprand(20,30,0.3)
    Ad = Array(A[i1, j1])
    B = SparseTensor(A)
    Bd = Array(B[i1, j1])
    Bd_ = run(sess, Bd)
    @test Ad≈Bd_
end

@testset "sum" begin
    s = sprand(10,20,0.2)
    S = SparseTensor(s)
    @test run(sess, sum(S)) ≈ sum(s)
    @test run(sess, sum(S,dims=1)) ≈ sum(Array(s),dims=1)[:]
    @test run(sess, sum(S,dims=2)) ≈ sum(Array(s),dims=2)[:]
end

@testset "dense_to_sparse" begin
    A = sprand(10,20,0.3)
    B = Array(A)
    @test run(sess, dense_to_sparse((B))) ≈ A
    @test run(sess, dense_to_sparse(constant(B))) ≈ A
end

@testset "spdiagm" begin
    a = rand(10)
    b = rand(9)
    A = spdiag(
        10,
        0=>a,
        -1=>b
    )
    B = diagm(
        0=>a,
        -1=>b
    )
    @test Array(run(sess, A)) ≈ B

    b = rand(7)
    A = spdiag(
        10,
        0=>a,
        -3=>b
    )
    B = diagm(
        0=>a,
        -3=>b
    )
    @test Array(run(sess, A)) ≈ B

    b = rand(7)
    A = spdiag(
        10,
        0=>a,
        -3=>b,
        3=>4b
    )
    B = diagm(
        0=>a,
        -3=>b,
        3=>4b
    )
    @test Array(run(sess, A)) ≈ B

    b = rand(7)
    A = spdiag(
        10,
        0=>a,
        -3=>b
    )
    B = diagm(
        0=>a,
        -3=>b
    )
    @test Array(run(sess, A)) ≈ B
end

@testset "hvcat" begin
    A = sprand(10,5,0.3)
    B = sprand(10,5,0.2)
    C = sprand(5,10,0.4)
    D = [A B;C]
    D_ = [SparseTensor(A) SparseTensor(B); SparseTensor(C)]
    @test run(sess, D_)≈D
end

@testset "find" begin
    A = sprand(10,10, 0.3)
    ii = Int64[]
    jj = Int64[]
    vv = Float64[]
    for i = 1:10
        for j = 1:10
            if A[i,j]!=0
                push!(ii, i)
                push!(jj, j)
                push!(vv, A[i,j])
            end
        end
    end
    a = SparseTensor(A)
    i, j, v = find(a)
    @test run(sess, i)≈ii
    @test run(sess, j)≈jj
    @test run(sess, v)≈vv

    @test run(sess, rows(a))≈ii
    @test run(sess, cols(a))≈jj
    @test run(sess, values(a))≈vv
end

@testset "sparse scatter update add" begin
    A = sprand(10,10,0.3)
    B = sprand(3,3,0.6)
    ii = [1;4;5]
    jj = [2;4;6]
    u = scatter_update(A, ii, jj, B)
    C = copy(A)
    C[ii,jj] = B
    @test run(sess, u)≈C

    u = scatter_add(A, ii, jj, B)
    C = copy(A)
    C[ii,jj] += B
    @test run(sess, u)≈C
end

@testset "constant sparse" begin
    A = sprand(10,10,0.3)
    B = constant(A)
    @test run(sess, B)≈A
end

@testset "get index" begin
    idof = [false;true]
    M = spdiag(constant(ones(2)))
    Md = M[idof, idof]
    @test run(sess, Md) ≈ sparse(reshape([1.0],1,1))
end

@testset "sparse_factorization_and_solve" begin
    A = sprand(10,10,0.7)
    rhs1 = rand(10)
    rhs2 = rand(10)
    Afac = factorize(constant(A))
    v1 = Afac\rhs1
    v2 = Afac\rhs2

    @test norm(run(sess, v1) - A\rhs1)<1e-10
    @test norm(run(sess, v2) - A\rhs2)<1e-10
end


text:
## Particle filter initialization ##
export pf_initialize

ParticleFilterState(trs::Vector{T}) where {T <: Trace} =
    ParticleFilterState{T}(trs, Vector{T}(undef, length(trs)),
                           fill(0., length(trs)), 0., collect(1:length(trs)))

ParticleFilterState(trs::Vector{T}, ws::Vector{Float64}) where {T <: Trace} =
    ParticleFilterState{T}(trs, Vector{T}(undef, length(trs)),
                           ws, 0., collect(1:length(trs)))

"""
    pf_initialize(model::GenerativeFunction, model_args::Tuple,
                  observations::ChoiceMap, n_particles::Int, dynamic::Bool=false)
    pf_initialize(model::GenerativeFunction, model_args::Tuple,
                  observations::ChoiceMap, proposal::GenerativeFunction,
                  proposal_args::Tuple, n_particles::Int, dynamic::Bool=false)

Initialize the state of a particle filter, generating traces (i.e. particles)
from the `model` constrained to the provided `observations`. A custom `proposal`
can optionally be used to propose choices for the initial set of traces.

If `dynamic` is `true`, the particle filter will not be specialized to a fixed
trace type, allowing for sequential Monte Carlo inference over a sequence
of models with differing structure, at the expense of more memory usage.
"""
function pf_initialize(model::GenerativeFunction{T,U}, model_args::Tuple,
                       observations::ChoiceMap, n_particles::Int,
                       dynamic::Bool=false) where {T,U}
    traces = Vector{Any}(undef, n_particles)
    log_weights = Vector{Float64}(undef, n_particles)
    for i=1:n_particles
        (traces[i], log_weights[i]) = generate(model, model_args, observations)
    end
    return dynamic ?
        ParticleFilterState{Trace}(traces, Vector{Trace}(undef, n_particles),
                               log_weights, 0., collect(1:n_particles)) :
        ParticleFilterState{U}(traces, Vector{U}(undef, n_particles),
                               log_weights, 0., collect(1:n_particles))
end

function pf_initialize(model::GenerativeFunction{T,U}, model_args::Tuple,
                       observations::ChoiceMap, proposal::GenerativeFunction,
                       proposal_args::Tuple, n_particles::Int,
                       dynamic::Bool=false) where {T,U}
    traces = Vector{Any}(undef, n_particles)
    log_weights = Vector{Float64}(undef, n_particles)
    for i=1:n_particles
        (prop_choices, prop_weight, _) = propose(proposal, proposal_args)
        (traces[i], model_weight) =
            generate(model, model_args, merge(observations, prop_choices))
        log_weights[i] = model_weight - prop_weight
    end
    return dynamic ?
        ParticleFilterState{Trace}(traces, Vector{Trace}(undef, n_particles),
                               log_weights, 0., collect(1:n_particles)) :
        ParticleFilterState{U}(traces, Vector{U}(undef, n_particles),
                               log_weights, 0., collect(1:n_particles))
end


text:
module Snappy

include("LibSnappy.jl")

export compress, uncompress

# snappy status
const SnappyOK             = Cint(0)
const SnappyInvalidInput   = Cint(1)
const SnappyBufferTooSmall = Cint(2)

function statuscheck(st, op::AbstractString="compress")
    st == LibSnappy.SNAPPY_OK && return nothing
    error("snappy failed to $op: " * if st == LibSnappy.SNAPPY_INVALID_INPUT
        "invalid input"
    elseif st == LibSnappy.SNAPPY_BUFFER_TOO_SMALL
        "buffer too small"
    else
        "unknown error code ($st)"
    end)
end

function compress(input::Vector{UInt8}; validate::Bool=false)
    ilen = length(input)
    maxlen = snappy_max_compressed_length(UInt(ilen))
    compressed = Array{UInt8}(undef, maxlen)
    olen, st = snappy_compress(input, compressed)
    statuscheck(st)
    resize!(compressed, olen)
    compressed
end
compress(input::AbstractVector{UInt8}) = compress(convert(Vector{UInt8}, input))

function uncompress(input::Array{UInt8})
    ilen = length(input)
    explen, st = snappy_uncompressed_length(input)
    if st != SnappyOK
        error("faield to guess the length of the uncompressed data (the compressed data may be broken?)")
    end
    uncompressed = Array{UInt8}(undef, explen)
    olen, st = snappy_uncompress(input, uncompressed)
    statuscheck(st, "uncompress")
    if explen ≠ olen
        error("snappy expected uncompressed length $explen but got $olen")
    end
    resize!(uncompressed, olen)
    uncompressed
end
uncompress(input::AbstractVector{UInt8}) = uncompress(convert(Vector{UInt8}, input))

# Low-level Interfaces

function snappy_compress(input::Vector{UInt8}, compressed::Vector{UInt8})
    olen = Ref{Csize_t}(length(compressed))
    st = LibSnappy.snappy_compress(input, length(input), compressed, olen)
    olen[], Int(st)
end

function snappy_uncompress(compressed::Vector{UInt8}, uncompressed::Vector{UInt8})
    olen = Ref{Csize_t}(length(uncompressed))
    st = LibSnappy.snappy_uncompress(compressed, length(compressed), uncompressed, olen)
    olen[], st
end

function snappy_max_compressed_length(source_length::Integer)
    LibSnappy.snappy_max_compressed_length(UInt(source_length))
end

function snappy_uncompressed_length(compressed::Vector{UInt8})
    result = Ref{Csize_t}(0)
    st = LibSnappy.snappy_uncompressed_length(compressed, length(compressed), result)
    result[], st
end

function snappy_validate_compressed_buffer(compressed::Vector{UInt8})
    LibSnappy.snappy_validate_compressed_buffer(compressed, length(compressed))
end

end # module


text:
module Flow

using DataStructures, HttpCommon, Morsel, JSON, Logging

app = Morsel.app()
state = {:console => {:input => "",
                      :output => ""},

         :tabs => {{:id => -1,
                    :name => "test"
                    :values => [{:name => "x",
                                 :type => :constant,
                                 :value => 8
                                 }]}},
         :next_tab_id => 0,
         :active_tab_id => -1,

         :dirty_queue => {}
         }

## utility
load_static(path) = open(readall, Pkg.dir("Flow") * "/" * path)

function get_in(m, ks)
    n = m
    for k = ks
        n = n[k]
    end
    n
end

function assoc_in(m, ks, v)
    merge(m, assoc_in(m, ks, v, 1))
end
function assoc_in(m, ks, v, i)
    if i < length(ks)
        if haskey(m, ks[i])
            merge(m, assoc_in(m[ks[i]], ks, v, i+1))
        else
            n = [ks[i] => assoc_in(Dict(), ks, v, i+1)]
            merge(m, n)
        end
    else
        merge(m, {ks[end] => v})
    end
end

function assoc_in!(m, ks, v)
    n = m
    for k = ks[1:end-1]
        if !haskey(n, k)
            n[k] = Dict()
        end
        n = n[k]
    end
    n[ks[end]] = v
end

## rest
get(app, "/eval") do req, res
    res.headers["Content-Type"] = mimetypes["json"]
    text = decodeURI(urlparam(req, :text))
    info("EVAL " * text)
    ptext = parse(text)
    json({:value => eval(ptext)})
end

put(app, "/store") do req, res
    ks = map(symbol, split(urlparam(req, :keys), ' '))
    v = decodeURI(urlparam(req, :val))
    if v == None
        println("none")
        v = ""
    end
    info(@sprintf("STORE %s => %s (%s)", ks, v, typeof(v)))
    store_silent!(ks, v)
    "ok"
end

get(app, "/fetch") do req, res
    res.headers["Content-Type"] = mimetypes["json"]
    ks = map(symbol, split(urlparam(req, :keys), ' '))
    info(@sprintf("FETCH %s", ks))
    json({:value => fetch(ks)})
end

get(app, "/dirty") do req, res
    res.headers["Content-Type"] = mimetypes["json"]

    if isempty(state[:dirty_queue])
        json({:status => :empty})
    else
        ks = shift!(state[:dirty_queue])
        res = json({:status => :ready,
              :keys => ks,
              :value => fetch(ks)})
        info(@sprintf("DIRTY %s", ks))
        debug(@sprintf("DIRTY => %s", res))
        res
    end
end

## process
function start(loglevel=INFO)
    Logging.configure(level=loglevel)
    Morsel.start(app, 8080)
end

function store_silent!(ks, v)
    debug(@sprintf("store_silent! %s => %s", ks, v))
    assoc_in!(state, ks, v)
end
function store!(ks, v)
    debug(@sprintf("store! %s => %s", ks, v))
    store_silent!(ks, v)
    mark_dirty!(ks)
end
fetch(ks) = get_in(state, ks)

function mark_dirty!(ks)
    push!(state[:dirty_queue], ks)
end

function sshow(val)
    ## http://thenewphalls.wordpress.com/2014/03/21/capturing-output-in-julia/
    originalSTDOUT = STDOUT
    (outRead, outWrite) = redirect_stdout()
    show(val)
    close(outWrite)

    data = readavailable(outRead)
    close(outRead)
    redirect_stdout(originalSTDOUT)

    data
end

eval_console() = store!([:console, :output],
                        sshow(eval(parse(fetch([:console, :input])))))

function new_tab(name)
    ts = fetch([:tabs])
    id = fetch([:next_tab_id])
    store_silent!([:next_tab_id], id + 1)

    info(@sprintf("new_tab: '%s' (%i)", name, id))

    push!(ts, {:id => id,
               :name => name,
               :active => true})

    store!([:tabs], ts)
    id
end

find_named(ms, n) = filter(m-> m[:name] == n, ms)

function tab_value(tab, value)
    tabs = find_named(state[:tabs], tab)
    if length(tabs) == 0
        throw("no such tab " * tab)
    elseif length(tabs) > 1
        throw("multiple tabs named " * tab)
    end

    vals = find_named(tabs[1], value)
    if length(vals) == 0
        throw("no such value " * value)
    elseif length(vals) > 1
        throw("multiple values named " * value)
    end

    vals[1]
end

end # module

max_stars_repo_name: farrellm/Flow

max_issues_repo_name: farrellm/Flow

max_forks_repo_name: farrellm/Flow


text:
@testset "incomplete beta" begin
    # Compared with Wolfram Engine's BetaRegularized
    @testset "a=0.5, b=0.5, x=$x" for (x, val) in zip(
        0.0:0.01:1.0,
        [0.e-25,0.06376856085851985,0.0903344706017331,0.1108246866044594,0.1281884336979499,
         0.1435662931287063,0.1575424254240443,0.170463428383199,0.1825548909938284,0.1939733680413566,
         0.2048327646991335,0.2152190255436425,0.2251989006184253,0.2348254690531804,0.2441417753461602,
         0.253183311106635,0.2619797608689092,0.2705562622683007,0.2789343361135709,0.2871325862574125,
         0.2951672353008665,0.3030525408934699,0.3108011236524005,0.3184242286319003,0.3259319361244763,
         0.3333333333333333,0.3406366554771171,0.3478494027630338,0.3549784381258894,0.3620300695016579,
         0.3690101195655454,0.375923985234415,0.3827766887550388,0.3895729218328954,0.3963170839725418,
         0.4030133159793217,0.4096655293982669,0.4162774325283954,0.42285255354095,0.4293942611422476,
         0.4359057831510251,0.4423902233029032,0.4488505765480783,0.4552897430704916,0.4617105412258534,
         0.4681157195707401,0.4745079681345231,0.48088992906936,0.4872642068002451,0.49363337778673,0.5,
         0.50636662221327,0.5127357931997549,0.51911007093064,0.5254920318654769,0.5318842804292599,
         0.5382894587741466,0.5447102569295084,0.5511494234519217,0.5576097766970968,0.5640942168489749,
         0.5706057388577524,0.57714744645905,0.5837225674716046,0.5903344706017331,0.5969866840206783,
         0.6036829160274582,0.6104270781671046,0.6172233112449612,0.624076014765585,0.6309898804344546,
         0.6379699304983421,0.6450215618741106,0.6521505972369662,0.6593633445228829,0.6666666666666667,
         0.6740680638755237,0.6815757713680997,0.6891988763475995,0.6969474591065301,0.7048327646991335,
         0.7128674137425875,0.7210656638864291,0.7294437377316993,0.7380202391310908,0.746816688893365,
         0.7558582246538398,0.7651745309468196,0.7748010993815747,0.7847809744563575,0.7951672353008665,
         0.8060266319586434,0.8174451090061716,0.829536571616801,0.8424575745759557,0.8564337068712937,
         0.8718115663020501,0.8891753133955406,0.9096655293982669,0.9362314391414802,1.])
        @test beta_inc(0.5, 0.5, x, 1.0 - x)[1] ≈ val rtol=1e-14
    end
    @testset "a=0.5, b=0.8, x=$x" for (x, val) in zip(
        0.0:0.01:1.0,
        [0.e-24,0.08704165448867688,0.1231783937242339,0.1509643905099393,0.1744376469592185,
         0.1951612417150757,0.2139363526723944,0.2312386644356416,0.247377962648835,0.2625693322193375,
         0.2769693157845956,0.2906960523171838,0.3038413001517716,0.3164780157848126,0.3286653449719896,
         0.3404520303383471,0.3518788093124823,0.3629801455598847,0.3737855072250243,0.3843203290282974,
         0.394606748829225,0.4046641800910133,0.4145097628356969,0.424158723204309,0.4336246632967352,
         0.4429197971428956,0.4520551445675527,0.4610406917926526,0.4698855255076289,0.4785979455868333,
         0.4871855604805254,0.4956553684394569,0.5040138270750072,0.5122669132520822,0.5204201749213043,
         0.5284787761920692,0.5364475367081124,0.5443309661970611,0.5521332949136762,0.5598585005745779,
         0.5675103322837338,0.5750923318679227,0.5826078529759856,0.5900600782419893,0.5974520347681598,
         0.6047866081468104,0.612066555210059,0.6192945156707836,0.6264730227971188,0.6336045132451182,
         0.6406913361594411,0.6477357616396057,0.6547399886591094,0.6617061525162714,0.6686363318887601,
         0.6755325555582559,0.6823968088674286,0.6892310399682792,0.6960371659188544,0.7028170786843475,
         0.7095726510986684,0.7163057428437263,0.7230182065060113,0.7297118937736915,0.7363886618425424,
         0.7430503801058253,0.7496989372120396,0.7563362485856976,0.7629642645204542,0.7695849789717461,
         0.7762004391984992,0.7828127564316306,0.7894241177826232,0.7960367996505097,0.8026531829430502,
         0.8092757705016214,0.8159072072147228,0.8225503034294581,0.8292080624343213,0.8358837130049195,
         0.8425807482983071,0.8493029727827216,0.8560545594443592,0.8628401202921954,0.869664794295266,
         0.8765343585075641,0.8834553705447895,0.8904353542429415,0.8974830460562918,0.9046087289748337,
         0.9118246961145984,0.9191459128120601,0.9265909945948355,0.9341837118759302,0.941955425134812,
         0.9499492886199519,0.9582281548957946,0.9668913290074657,0.9761171938969928,0.9863137214671616,1.])
        @test beta_inc(0.5, 0.8, x, 1.0 - x)[1] ≈ val rtol=1e-14
    end
    @testset "a=0.9, b=0.8, x=$x" for (x, val) in zip(
        0.0:0.01:1.0,
        [0.e-45,0.01287348375077395,0.0240457784957195,0.0346688685639659,0.04495791979284836,
         0.0550110492927223,0.06488448192684902,0.07461461051118363,0.08422688488805261,0.0937400884197129,
         0.1031686494864254,0.1125239993592263,0.1218154216032231,0.1310506102439909,0.1402360514251891,
         0.1493772930740891,0.1584791407604112,0.1675458033326384,0.1765810034326174,0.1855880628612014,
         0.1945699695603456,0.2035294309106093,0.2124689166764504,0.2213906940068579,0.2302968562599333,
         0.2391893469703827,0.2480699799571045,0.2569404563342797,0.2658023790171647,0.2746572651853306,
         0.2835065570691575,0.2923516313514671,0.3011938074192502,0.3100343546562188,0.3188744989322686,
         0.3277154284186073,0.3365582988355958,0.3454042382230052,0.3542543513084796,0.3631097235387736,
         0.3719714248292748,0.3808405130799951,0.3897180375002914,0.3986050417798133,0.4075025671393664,
         0.4164116552923807,0.4253333513453621,0.4342687066639933,0.4432187817303744,0.452184649016201,
         0.4611673958964416,0.4701681276282771,0.4791879704207046,0.4882280746212869,0.4972896180480853,
         0.5063738094968672,0.5154818924563003,0.5246151490670925,0.5337749043650042,0.5429625308524718,
         0.5521794534493672,0.5614271548803832,0.5707071815648901,0.5800211500851545,0.5893707543209051,
         0.5987577733528309,0.608184080255283,0.6176516519199679,0.6271625800787011,0.6367190837255744,
         0.6463235231787504,0.6559784160716227,0.6656864556250031,0.6754505316299672,0.6852737546699122,
         0.695159484236884,0.7051113615604331,0.7151333481797695,0.7252297715695021,0.7354053795009442,
         0.74566540532109,0.7560156470130347,0.7664625638438512,0.7770133957285906,0.7876763123275041,
         0.7984606016409664,0.8093769119481182,0.8204375671426882,0.8316569852132864,0.8430522452167303,
         0.8546438740859578,0.8664569696979065,0.8785228586162501,0.8908816447296843,0.9035863305210574,
         0.9167099250794181,0.9303587974356043,0.9447009522476523,0.9600378941965647,0.977058098680342,1.])
        @test beta_inc(0.9, 0.8, x, 1.0 - x)[1] ≈ val rtol=1e-14
    end
    @testset "a=80.9, b=0.8, x=$x" for (x, val) in zip(
        0.0:0.01:1.0,
        [0.e-4037,5.659850188208874e-163,1.279388616124076e-138,2.257693896707542e-124,
         2.897994605366003e-114,2.009164686416781e-106,5.124902812310911e-100,1.338436546592981e-94,
         6.592846323627736e-90,9.083404801612933e-86,4.581171164611029e-82,1.024712535445972e-78,
         1.171290435847841e-75,7.618998306182543e-73,3.066399996029232e-70,8.159340151736502e-68,
         1.514230271101736e-65,2.047520095589173e-63,2.091671353695881e-61,1.664040903808287e-59,
         1.057759838815721e-57,5.49132667443713e-56,2.372578688736847e-54,8.671833527840728e-53,
         2.719793317701389e-51,7.411825180712226e-50,1.774472557976592e-48,3.768914466846032e-47,
         7.163917757923168e-46,1.228169152371666e-44,1.912408289817481e-43,2.721825167444e-42,
         3.561050588818447e-41,4.305131971329379e-40,4.832096799800394e-39,5.057055939273883e-38,
         4.954351648811546e-37,4.560116807292394e-36,3.956512105211734e-35,3.245869269592197e-34,
         2.525020123177697e-33,1.867480042455247e-32,1.316327433960537e-31,8.862853181709876e-31,
         5.712179225668737e-30,3.531051129157696e-29,2.097389394183402e-28,1.199152837592699e-27,
         6.609838550976534e-27,3.517914083771492e-26,1.810402708594319e-25,9.020730533511399e-25,
         4.357446464530653e-24,2.042971920816074e-23,9.307222587505059e-23,4.124448596487034e-22,
         1.779654898083663e-21,7.4841565165423e-21,3.070285970414923e-20,1.22974643568586e-19,
         4.812895987534055e-19,1.841991715270897e-18,6.898889614425955e-18,2.530381730672556e-17,
         9.094925675449364e-17,3.205503996194959e-16,1.108518744159741e-15,3.763511868841169e-15,
         1.255134520239131e-14,4.114026783167532e-14,1.326012138884187e-13,4.204809280909274e-13,
         1.312411918693827e-12,4.033834438688928e-12,1.221467371969832e-11,3.645413830024844e-11,
         1.072734311570093e-10,3.113805924559142e-10,8.918947491417688e-10,2.521869395260185e-9,
         7.041721770169991e-9,1.942404070962948e-8,5.294920235920981e-8,1.426893004118986e-7,
         3.802663384681951e-7,1.002539998693497e-6,2.615699030509946e-6,6.756237677059152e-6,
         0.00001728297704388556,0.00004380303345697832,0.0001100404885846187,0.0002741427949192357,
         0.0006776797945177599,0.001663395040736671,0.004057633817143372,0.009848744173397283,
         0.02382855932322502,0.05763883015096727,0.1402012176136536,0.348133320555652,1.])
        @test beta_inc(80.9, 0.8, x, 1.0 - x)[1] ≈ val rtol=1e-13
    end
    @testset "a=1.7, b=10.5, x=$x" for (x, val) in zip(
        0.0:0.01:1.0,
        [0.e-83,0.01396204974701841,0.04273844739996358,0.0802355082453466,0.1233212707222279,
         0.169883883601507,0.2183876645881966,0.2676795794438451,0.3168824713253145,0.3653277197583937,
         0.4125089764419372,0.4580484843629353,0.5016715394987723,0.5431865527277862,0.5824691495850261,
         0.6194492913831362,0.654100724564198,0.6864322667088425,0.7164805688276856,0.7443040822000412,
         0.7699780198879489,0.7935901474806426,0.8152372703488115,0.8350223093424488,0.8530518758192943,
         0.8694342717267099,0.8842778522619746,0.8976896981615276,0.9097745524576825,0.9206339829823968,
         0.9303657372794607,0.9390632611240216,0.9468153557047174,0.9537059518240998,0.9598139823148841,
         0.9652133363295071,0.9699728812994702,0.974156540228006,0.9778234136145069,0.9810279367445075,
         0.9838200643417606,0.9862454756915709,0.9883457943259302,0.9901588172271333,0.9917187492710794,
         0.993056439306095,0.9941996148579756,0.9951731129758233,0.9959991051938595,0.9966973149884605,
         0.997285226463178,0.997778283302743,0.9981900773047402,0.9985325260289722,0.998816039303301,
         0.9990496744943506,0.9992412805949365,0.9993976313002439,0.9995245473440961,0.9996270084474357,
         0.9997092552954444,0.9997748820094408,0.9998269196165406,0.9998679110455995,0.9998999781936217,
         0.9999248816138942,0.9999440733768028,0.9999587436476631,0.9999698615139476,0.9999782105779107,
         0.9999844198105916,0.9999889901402758,0.9999923172233462,0.9999947108186769,0.9999964111588281,
         0.9999976026827751,0.9999984254661612,0.9999989846564765,0.9999993581924477,0.9999996030595673,
         0.9999997603073096,0.9999998590283993,0.9999999194766476,0.9999999554775051,0.9999999762646713,
         0.9999999878569335,0.9999999940719059,0.9999999972575227,0.9999999988080026,0.999999999518498,
         0.999999999821743,0.999999999940621,0.99999999998264,0.999999999995699,0.999999999999142,
         0.999999999999873,0.999999999999988,0.999999999999999,1.,1.,1.])
        @test beta_inc(1.7, 10.5, x, 1.0 - x)[1] ≈ val rtol=1e-14
    end
    @testset "a=100.5, b=100.5, x=$x" for (x, val) in zip(
        0.0:0.01:1.0,
        [0.e-4966,3.356014382599029e-143,2.21351837164392e-113,4.014552683228967e-96,
         5.21091312517625e-84,1.020389199029434e-74,3.267737964647878e-67,6.101573725809464e-61,
         1.418580292320946e-55,6.698584764694289e-51,8.971972641791594e-47,4.325143541605987e-43,
         8.94664782857427e-40,9.074256619313614e-37,5.005701167291761e-34,1.630336905229558e-31,
         3.349388657777907e-29,4.581338772370011e-27,4.362832768445309e-25,3.002816845923933e-23,
         1.541642496405598e-21,6.064729213175483e-20,1.870837362242701e-18,4.616584448897111e-17,
         9.27255359484341e-16,1.539075197361062e-14,2.139403432844863e-13,2.520059271944931e-12,
         2.541836676624415e-11,2.215842537778625e-10,1.683408890736215e-9,1.122874450435998e-8,
         6.620173165982372e-8,3.470763841511403e-7,1.626934118513085e-6,6.852626959254337e-6,
         0.00002605241980849797,0.00008977130313754974,0.0002814332703651759,0.0008055416902949821,
         0.002112019781014849,0.00508797075711991,0.01129547597461094,0.0231746884205069,
         0.04406489065169599,0.07787036008215303,0.1282703308247337,0.1975636232169813,0.2854862534447546,
         0.3885022505229126,0.5,0.6114977494770874,0.7145137465552454,0.8024363767830187,0.8717296691752663,
         0.922129639917847,0.955935109348304,0.9768253115794931,0.9887045240253891,0.9949120292428801,
         0.9978879802189852,0.999194458309705,0.9997185667296348,0.9999102286968625,0.9999739475801915,
         0.9999931473730407,0.9999983730658815,0.9999996529236158,0.9999999337982683,0.9999999887712555,
         0.9999999983165911,0.999999999778416,0.999999999974582,0.99999999999748,0.999999999999786,
         0.999999999999985,0.999999999999999,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,1.,
         1.,1.,1.])
        @test beta_inc(100.5, 100.5, x, 1.0 - x)[1] ≈ val rtol=1e-13
    end
    @testset "a=1.5, b=1e-20, x=$x" for (x, val) in zip(
        0.0:0.01:1.0,
        [0.e-95,6.706954621511613e-24,1.908573860989784e-23,3.527823559800854e-23,
         5.465108108164382e-23,7.685476511624814e-23,1.016654428420256e-22,1.289283267454361e-22,
         1.585315312526521e-22,1.903920840622343e-22,2.244476844084102e-22,2.606514996691386e-22,
         2.989687056036453e-22,3.393740920821702e-22,3.818503532521727e-22,4.26386835711828e-22,
         4.729786038720361e-22,5.21625731553983e-22,5.72332759176175e-22,6.251082749838933e-22,
         6.799645911929102e-22,7.369174942192775e-22,7.95986053851056e-22,8.571924801968895e-22,
         9.205620200883313e-22,9.861228866810969e-22,1.053906217532698e-21,1.12394605758909e-21,
         1.196279364399263e-21,1.270946033569195e-21,1.347988943018207e-21,1.427454015049574e-21,
         1.50939029562088e-21,1.593850050518841e-21,1.680888878423608e-21,1.77056584110116e-21,
         1.862943611198906e-21,1.958088638346952e-21,2.056071334492481e-21,2.156966279623312e-21,
         2.260852449274566e-21,2.367813465464866e-21,2.477937872981194e-21,2.591319443230187e-21,
         2.708057508204485e-21,2.828257327482878e-21,2.952030491600198e-21,3.079495365596476e-21,
         3.210777577095268e-21,3.346010553881064e-21,3.48533611665991e-21,3.628905133513816e-21,
         3.776878243519303e-21,3.929426658119426e-21,4.086733050147924e-21,4.248992541941733e-21,
         4.416413805790527e-21,4.589220292116253e-21,4.767651603322342e-21,4.951965034288911e-21,
         5.142437304125938e-21,5.339366508167657e-21,5.543074324470629e-21,5.753908515483088e-21,
         5.972245773362194e-21,6.19849496698491e-21,6.433100860486935e-21,6.67654838776915e-21,
         6.929367585613005e-21,7.192139310873532e-21,7.465501896037407e-21,7.75015893407293e-21,
         8.046888430431106e-21,8.356553620653897e-21,8.680115830952776e-21,9.018649862807561e-21,
         9.373362520215628e-21,9.745615082692113e-21,1.013695077726479e-20,1.05491286460553e-20,
         1.098416568357789e-20,1.14443897916644e-20,1.193250706559674e-20,1.245168833304345e-20,
         1.300568195962486e-20,1.359896310974523e-20,1.423693458276194e-20,1.492620220780857e-20,
         1.567496069508049e-20,1.649354777878116e-20,1.739526322363106e-20,1.839762082220938e-20,
         1.952434090962006e-20,2.080868037972059e-20,2.229934708861127e-20,2.407185551639959e-20,
         2.625271544895813e-20,2.907909078985462e-20,3.308342534794925e-20,3.996470818039522e-20,1.])
        @test beta_inc(1.5, 1e-20, x, 1.0 - x)[1] ≈ val rtol=1e-14
    end
    @testset "a=1e-20, b=0.99, x=$x" for (x, val) in zip(
        0.0:0.01:1.0,
        [1.0,4.621640571507914e-20,3.928392625456373e-20,3.522826232996826e-20,3.235042350168417e-20,
         3.011796455118398e-20,2.82937201372135e-20,2.675117900735502e-20,2.541482518522368e-20,
         2.423594928780016e-20,2.318129286271936e-20,2.222713398379447e-20,2.135595723375987e-20,
         2.05544611885328e-20,1.981230641875294e-20,1.912129648210598e-20,1.847482377915803e-20,
         1.786748370079641e-20,1.729479923207514e-20,1.675302011465343e-20,1.623897358450671e-20,
         1.574995156483639e-20,1.528362412553572e-20,1.483797219538297e-20,1.441123460536074e-20,
         1.400186594951488e-20,1.360850271543896e-20,1.322993581024732e-20,1.286508808543862e-20,
         1.251299580739944e-20,1.217279327043997e-20,1.184369993372369e-20,1.152500960103253e-20,
         1.121608126599496e-20,1.091633132431086e-20,1.062522691510309e-20,1.034228020045673e-20,
         1.006704342884576e-20,9.799104656964201e-21,9.538084027304864e-21,9.283630517029258e-21,
         9.035419088275798e-21,8.793148181840232e-21,8.556537505729393e-21,8.325326077896178e-21,
         8.099270488866332e-21,7.878143355244009e-21,7.661731939451283e-21,7.449836914688898e-21,
         7.242271257138218e-21,7.03885924996661e-21,6.839435585837694e-21,6.643844556434574e-21,
         6.451939319035191e-21,6.263581231480807e-21,6.078639247988984e-21,5.896989369212464e-21,
         5.718514140760424e-21,5.543102195099824e-21,5.370647832359387e-21,5.201050636081688e-21,
         5.034215120421851e-21,4.870050405684712e-21,4.708469919434276e-21,4.549391120707029e-21,
         4.392735245120113e-21,4.238427068891405e-21,4.086394689985318e-21,3.936569324769163e-21,
         3.788885118712737e-21,3.643278969790933e-21,3.499690363356917e-21,3.358061217343224e-21,
         3.21833573672026e-21,3.080460276196264e-21,2.944383210178798e-21,2.810054809033594e-21,
         2.677427120668968e-21,2.546453856438284e-21,2.41709028028191e-21,2.289293099913057e-21,
         2.163020358672793e-21,2.038231326414587e-21,1.914886387391528e-21,1.79294692255393e-21,
         1.67237518283336e-21,1.553134148749398e-21,1.435187369793683e-21,1.318498774123754e-21,
         1.20303243443583e-21,1.088752268197901e-21,9.756216372555088e-22,8.63602788212576e-22,
         7.526560302762232e-22,6.427384567004279e-22,5.338018164784628e-22,4.25788652801385e-22,
         3.186244273953916e-22,2.121983888512568e-22,1.063002842170579e-22,0])
        @test beta_inc(1e-20, 0.99, x, 1.0 - x)[2] ≈ val rtol=1e-14
    end

    @test beta_inc(1.5, 200.5,  0.07,0.93)[1] ≈ 0.99999790408564
    @test beta_inc(1e-20, 0.000001, 0.2)[2] ≈ 1.0000013862929421e-14

    # test promotions and return types
    for T in (Float64, Float32)
        a, b = randexp(T, 2)
        x = rand(T)
        @test beta_inc(a, b, x) isa Tuple{T,T}
        @test beta_inc(a, b, x, 1 - x) isa Tuple{T,T}
        @test beta_inc(a, b, x) == beta_inc(a, b, x, 1 - x)
    end
    a = randexp()
    b = randexp(Float32)
    x = rand(Float32)
    @test beta_inc(a, b, x) isa Tuple{Float64,Float64}
    @test beta_inc(a, b, x) == beta_inc(a, Float64(b), Float64(x))
    @test beta_inc(a, b, x, 1 - x) isa Tuple{Float64,Float64}
    @test beta_inc(a, b, x, 1 - x) == beta_inc(a, Float64(b), Float64(x), 1 - x)

    @test SpecialFunctions.loggammadiv(13.89, 21.0001) ≈ log(gamma(big(21.0001))/gamma(big(21.0001)+big(13.89)))
    @test SpecialFunctions.stirling_corr(11.99, 100.1) ≈ SpecialFunctions.stirling_error(11.99) + SpecialFunctions.stirling_error(100.1) - SpecialFunctions.stirling_error(11.99 + 100.1)
end

@testset "inverse of incomplete beta" begin
    f(a,b,p) = beta_inc_inv(a,b,p)[1]
    @test f(.5, .5, 0.6376856085851985E-01) ≈ 0.01
    @test f(.5, .5, 0.20483276469913355) ≈ 0.1
    @test f(.5, .5, 1.0000) ≈ 1.0000
    @test f(1.0, .5, 0.0) ≈ 0.00
    @test f(1.0, .5, 0.5012562893380045E-02) ≈ 0.01
    @test f(1.0, .5, 0.5131670194948620E-01) ≈ 0.1
    @test f(1.0, .5, 0.2928932188134525) ≈ 0.5
    @test f(1.0, 1.0, .5) ≈ 0.5
    @test f(2.0, 2.0, .028) ≈ 0.1
    @test f(2.0, 2.0, 0.104) ≈ 0.2
    @test f(2.0, 2.0, .216) ≈ 0.3
    @test f(2.0, 2.0, .352) ≈ 0.4
    @test f(2.0, 2.0, .5) ≈ 0.5
    @test f(2.0, 2.0, 0.648) ≈ 0.6
    @test f(2.0, 2.0, 0.784) ≈ 0.7
    @test f(2.0, 2.0, 0.896) ≈ 0.8
    @test f(2.0, 2.0, .972) ≈ 0.9
    @test f(5.5, 5.0, 0.4361908850559777) ≈ .5
    @test f(10.0, .5, 0.1516409096347099) ≈ 0.9
    @test f(10.0, 5.0, 0.8978271484375000E-01) ≈ 0.5
    @test f(10.0, 5.0, 1.00) ≈ 1.00
    @test f(10.0, 10.0, .5) ≈ .5
    @test f(20.0, 5.0, 0.4598773297575791) ≈ 0.8
    @test f(20.0, 10.0, 0.2146816102371739) ≈ 0.6
    @test f(20.0, 10.0, 0.9507364826957875) ≈ 0.8
    @test f(20.0, 20.0, .5) ≈ .5
    @test f(20.0, 20.0, 0.8979413687105918) ≈ 0.6
    @test f(30.0, 10.0, 0.2241297491808366) ≈ 0.7
    @test f(30.0, 10.0, 0.7586405487192086) ≈ 0.8
    @test f(40.0, 20.0, 0.7001783247477069) ≈ 0.7
    @test f(1.0, 0.5, 0.5131670194948620E-01) ≈ 0.1
    @test f(1.0, 0.5, 0.1055728090000841) ≈ 0.2
    @test f(1.0, 0.5, 0.1633399734659245) ≈ 0.3
    @test f(1.0, 0.5, 0.2254033307585166) ≈ 0.4
    @test f(1.0, 2.0, .36) ≈ 0.2
    @test f(1.0, 3.0, .488) ≈ 0.2
    @test f(1.0, 4.0, .5904) ≈ 0.2
    @test f(1.0, 5.0, .67232) ≈ 0.2
    @test f(2.0, 2.0, .216) ≈ 0.3
    @test f(3.0, 2.0, 0.837e-1) ≈ 0.3
    @test f(4.0, 2.0, 0.3078e-1) ≈ 0.3
    @test f(5.0, 2.0, 0.10935e-1) ≈ 0.3
    @test f(1.30625000, 11.75620000, 0.9188846846205182) ≈ 0.225609
    @test f(1.30625000, 11.75620000, 0.21053116418502513) ≈ 0.033557
    @test f(1.30625000, 11.75620000, 0.18241165418408148) ≈ 0.029522
    @test f(1000.0, 2.0, 9.0797754e-317) ≈ 0.48 # This one is a bit slow (but also hard)
end


text:
"""
A named tuple that contains the magic values that are used globally for
whatever purposes.
"""
const _constants = (
    default_stoich_show_size = 50_000,
    default_reaction_bound = 1e3,
    tolerance = 1e-6,
    sampling_keep_iters = 100,
    sampling_size = 1000,
    loopless_strict_inequality_tolerance = 1,
    extracellular_suffixes = ["_e", "[e]", "(e)"],
    exchange_prefixes = ["EX_", "Exch_", "Ex_", "R_EX_", "R_Ex", "R_Exch_"],
    biomass_strings = ["BIOMASS", "biomass", "Biomass"],
    keynames = (
        rxns = ["reactions", "rxns", "RXNS", "REACTIONS", "Reactions", "Rxns"],
        mets = ["metabolites", "mets", "METS", "METABOLITES", "Metabolites", "Mets"],
        genes = ["genes", "GENES", "Genes"],
        lbs = ["lbs", "lb", "lowerbounds", "lower_bounds"],
        ubs = ["ubs", "ub", "upperbounds", "upper_bounds"],
        stochiometry = ["S"],
        balance = ["b"],
        objective = ["c"],
        grrs = ["gene_reaction_rules", "grRules", "rules"],
        ids = ["id", "description"],
    ),
    gene_annotation_checks = (
        "ncbigene",
        "ncbigi",
        "refseq_locus_tag",
        "refseq_name",
        "refseq_synonym",
        "uniprot",
    ),
    reaction_annotation_checks = (
        "bigg.reaction",
        "biocyc",
        "ec-code",
        "kegg.reaction",
        "metanetx.reaction",
        "rhea",
        "sabiork",
        "seed.reaction",
    ),
    metabolite_annotation_checks = (
        "kegg.compound",
        "bigg.metabolite",
        "chebi",
        "inchi_key",
        "sabiork",
        "hmdb",
        "seed.compound",
        "metanetx.chemical",
        "reactome.compound",
        "biocyc",
    ),
    T = 298.15, # Kelvin
    R = 8.31446261815324e-3, # kJ/K/mol
)

const MAX_SENSE = MOI.MAX_SENSE
const MIN_SENSE = MOI.MIN_SENSE


text:
@testset "DataFrame integration" begin
    LearnBase.getobs(df::DataFrame, idx) = df[idx,:]
    StatsBase.nobs(df::DataFrame) = nrow(df)

    @testset "targets" begin
        y = [:a,:a,:b,:a,:b]
        df = DataFrame(x1 = rand(5), x2 = rand(5), y = y)
        @test y == targets(row->row.y, df)
    end

    @testset "eachtarget" begin
        y = [:a,:a,:b,:a,:b]
        df = DataFrame(x1 = rand(5), x2 = rand(5), y = y)
        iter = eachtarget(row->row.y, df)
        @test y == collect(iter)
    end

    @testset "gettarget" begin
        LearnBase.gettarget(col::Symbol, df::DataFrameRow) = df[col]
        y = [:a,:a,:b,:a,:b]
        df = DataFrame(x1 = rand(5), x2 = rand(5), y = y)
        @test y == targets(:y, df)
    end

    @testset "gettarget default" begin
        LearnBase.gettarget(df::DataFrameRow) = df[end]
        y = [:a,:a,:b,:a,:b]
        df = DataFrame(x1 = rand(5), x2 = rand(5), y = y)
        @test y == targets(df)
    end
end


text:
# Note: implement as generic, in order to also support views which are of type SubArray
function insertionsort!(a::AbstractVector, isless=isless)
    for i = 2:length(a)
        for j = i:-1:2
            isless(a[j], a[j-1]) || break
            swap!(a, j, j-1)
        end
    end
    a
end


text:
import StatsBase: mode

mutable struct NcDim
    name::String
    dimlen::Int
    vals::Union{AbstractArray,Nothing}
    atts::Dict
end

"""
    NcDim(name::AbstractString, vals::AbstractArray, atts::Dict = Dict())

Represents a NetCDF dimension of name `name` optionally holding the dimension values.
"""
NcDim(name::AbstractString, vals::AbstractArray, atts::Dict = Dict()) =
    NcDim(name, length(vals), vals, atts)

function Ipaper.names(dims::Vector{NcDim})
    map(x -> x.name, dims)
end

function Base.getindex(dims::Vector{NcDim}, name::AbstractString)
    names = Ipaper.names(dims)
    ind = indexin([name], names)[1]
    dims[ind]
end

# function Base.getindex(dims::Vector{NcDim}, name::Vector{AbstractString})
#     names = Ipaper.names(dims)
#     ind = indexin(name, names)
#     print(ind)
#     dims[ind]
# end
# Base.getindex(dims::Vector{NcDim}, name::AbstractString) = Base.getindex(dims, [name])
function nc_dim_shape(ds::NCdata, name = "time")
    n = ds.dim[name]
    vals = 1:n
    NcDim(name, n, vals, Dict())
    # NcDim(name, n, nothing, Dict()) # TODO
end

function nc_dim(ds::NCdata, name = "time")
    if haskey(ds, name)
        try
            x = ds[name]
            NcDim(name, x.var[:], Dict(x.attrib))
        catch
            nc_dim_shape(ds, name)
        end
    else
        nc_dim_shape(ds, name)
    end
end

function nc_dim(file::NCfiles, name = "time")
    nc_open(file) do ds
        nc_dim(ds, name)
    end
end

function nc_dims(ds::NCdata)
    items = keys(ds.dim) #|> reverse
    map(x -> nc_dim(ds, x), items)
end

function nc_dims(file::NCfiles)
    nc_open(file) do ds
        nc_dims(ds)
    end
end

function nc_dimsize(file::NCfiles)
    nc_open(file) do ds
        var = nc_bands(ds)[1]
        size(ds[var])
    end
end

"""
    nc_cellsize(ds::NCdata)
"""
function nc_cellsize(ds::NCdata)
    diflon = nc_dim(ds, "lon").vals |> diff
    diflat = nc_dim(ds, "lat").vals |> diff

    cell_x = mode(diflon)
    cell_y = mode(diflat)
    # regular or complex grid
    regular = length(unique(diflon)) == 1 && length(unique(diflat)) == 1
    cell_x, cell_y, regular
end

function nc_cellsize(file::AbstractString)
    nc_open(file) do ds
        nc_cellsize(ds)
    end
end

function nc_cellsize(files::Vector{<:AbstractString})
    n = length(files)
    cell_x = zeros(n)
    cell_y = zeros(n)
    regular = zeros(Bool, n)
    for i = 1:n
        cell_x[i], cell_y[i], regular[i] = nc_cellsize(files[i])
    end
    cell_x, cell_y, regular
end

function ncvar_dim(ds::NCdata, varname::Union{String, Nothing} = nothing; varid = 1)
    if varname === nothing
        varname = nc_bands(ds)[varid]
    end

    dims = nc_dims(ds)
    ids = ds[varname].var.dimids
    ids = collect(ids) .+ 1
    dims[ids]
end

function ncvar_dim(file::NCfiles, varname::Union{String,Nothing} = nothing; kwargs...)
    nc_open(file) do ds
        ncvar_dim(ds, varname; kwargs...)
    end
end


text:
# this file defines iterators used for looping over a grid
mutable struct CellData
    current_cellid::Int
    cell
end

struct CellIterator{dim,T}
    mesh::AbstractPolytopalMesh{dim,T}
    celldata::CellData
    function CellIterator{dim,T}(mesh::AbstractPolytopalMesh{dim,T}) where {dim,T}
        cellid = 0
        cell = nothing
        return new{dim,T}(mesh, CellData(cellid,cell))
    end
end

CellIterator(mesh::AbstractPolytopalMesh{dim,T}) where {dim,T} =
    CellIterator{dim,T}(mesh)

# iterator interface
function Base.iterate(ci::CellIterator, state = 1)
    if state > getncells(ci.mesh)
        return nothing
    else
        return (reinit!(ci, state), state+1)
    end
end
Base.length(ci::CellIterator)  = getncells(ci.mesh)

Base.IteratorSize(::Type{T})   where {T<:CellIterator} = Base.HasLength() # this is default in Base
Base.IteratorEltype(::Type{T}) where {T<:CellIterator} = Base.HasEltype() # this is default in Base
Base.eltype(::Type{T})         where {T<:CellIterator} = T

# utility
@inline cellid(ci::CellIterator) = ci.celldata.current_cellid

function reinit!(ci::CellIterator{dim}, i::Int) where {dim}
    ci.celldata.current_cellid = i
    ci.celldata.cell = getcell(ci.mesh,i)
    return ci
end

@inline reinit!(fs::AbstractFEMFunctionSpace{dim,T,FE}, ci::CellIterator{dim,T}) where {dim,T,FE} = reinit!(fs, getverticescoords(ci.mesh, ci.celldata.cell))
function_value(f::Function, fs::AbstractFEMFunctionSpace{dim,T,FE}, ci::CellIterator,q_point::Int) where {dim,T,FE} = function_value(f,fs, ci.celldata.current_cellid,q_point)
celldofs!(v::Vector, dh::DofHandler, ci::CellIterator) = celldofs!(v, dh, ci.celldata.current_cellid)

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
import LinearAlgebra
### matrix free S (slicing) matrix

function binfunc(ii::Int32,binFactor::Int32)
	ii -=1;
	ii = div(ii,binFactor);
	ii +=1;
	return ii;
end

export sampleSlices,sampleSlicesT
function sampleSlices(u::Array{Float64,3},cellVolume::Float64,samplingBinning::Int64=1)
n_tup = size(u);
n = collect(n_tup);
traceLength = div(n[3],samplingBinning);
s = zeros(eltype(u),traceLength);
for k = 1:traceLength
	@inbounds s[k] = sum(view(u,:,:,(samplingBinning*(k-1)+1):(samplingBinning*k)));
end
s.*=cellVolume;
return s;
end

function sampleSlicesT(s::Array,cellVolume::Float64,u::Array{Float64,3})
traceLength = length(s);
samplingBinning = div(size(u,3),traceLength);
for k = 1:traceLength
	@inbounds u[:,:,(samplingBinning*(k-1)+1):(samplingBinning*k)] = cellVolume*s[k];
end
end

export getData
function getData(m::Array,pFor::DirectParam,doClear::Bool=false)
d = 0;
Mesh = pFor.Mesh;
n = Mesh.n;
ndips = 1;
samplingBinning = pFor.samplingBinning;

if pFor.method == RBFBased || pFor.method == RBF10Based || pFor.method == RBF5Based
	if pFor.method == RBF10Based
		numParamOfRBF = 10;
	else
		numParamOfRBF = 5;
	end

	sigmaH = getDefaultHeavySide();
	u,I1,J1,V1 = ParamLevelSetModelFunc(Mesh,m;computeJacobian = 1,sigma = sigmaH,bf = 1,numParamOfRBF = numParamOfRBF);
	J1 = sparse(I1,J1,V1,prod(n),length(m));

	d = u;
	#d = reshape(d,div(length(d),ndips),ndips);
	println("getData jacobian:",size(J1))
	#println("Jrot size:",size(Jrot))
	pFor.Jacobian = J1;
# elseif pFor.method == RBFBasedSimple2 || pFor.method == RBF10BasedSimple2
	# if pFor.method == RBF10BasedSimple2
		# numParamOfRBF = 10;
	# else
		# numParamOfRBF = 5;
	# end
	# #mrot,Jrot = rotateAndMoveRBFsimple(m,Mesh,theta_phi,b;computeJacobian = 1,numParamOfRBF=numParamOfRBF);
	# d,JacT = getDipDataRBF(pFor, m,numParamOfRBF);
	# # multiply Jacobian with Jrot, and then take the transpose.
	# pFor.Jacobian = JacT;
# elseif pFor.method == RBFBased || pFor.method == RBF10Based
	# if pFor.method == RBF10Based
		# numParamOfRBF = 10;
	# else
		# numParamOfRBF = 5;
	# end
	# ## here theta_phi is not constant - so we ignore the theta_phi in pFor
	# nRBF 				= div(length(m)- 5*pFor.ndipsAll,numParamOfRBF) ;
	# (m1,theta_phi,b) 	= splitRBFparamAndRotationsTranslations(m,nRBF,pFor.ndipsAll,numParamOfRBF);
	# theta_phi 			= theta_phi[pFor.workerSubIdxs,:];
	# b 					= b[pFor.workerSubIdxs,:];
	# mrot,Jrot 			= rotateAndMoveRBF(m1,Mesh,theta_phi,b;computeJacobian = 1,numParamOfRBF=numParamOfRBF);
	# d,JacT 				= getDipDataRBF(pFor, mrot,theta_phi,b,numParamOfRBF);
	# # multiply Jacobian with Jrot, and then take the transpose.

	# JacobianT = convert(SparseMatrixCSC{Float64,Int32},spzeros(length(m),prod(size(d))));
	# IpIdxs = getIpIdxs(pFor.workerSubIdxs,nRBF,pFor.ndipsAll,numParamOfRBF);

	# JacobianT[IpIdxs,:] = Jrot'*JacT;
	# pFor.Jacobian = JacobianT;
	# # pFor.Jacobian 		= Jrot'*JacT;
end
return d,pFor;
end

function getDirectDataRBF(pFor, m,numParamOfRBF = 5)
samplingBinning = pFor.samplingBinning;
Mesh = pFor.Mesh;
n = pFor.Mesh.n;
h = pFor.Mesh.h;
traceLength = div(n[3],samplingBinning);
ndips = 1;
d = zeros(Float32,traceLength,ndips);
lengthRBFparams = size(m,1);
Ihuge = zeros(Int32,0);
I1 = zeros(Int32,0);
J1 = zeros(Int32,0);
V1 = zeros(Float64,0);
sigmaH = getDefaultHeavySide();
u = zeros(prod(n));
dsu = zeros(prod(n));
Xc = convert(Array{Float32,2},getCellCenteredGrid(Mesh));
binningFactor = convert(Int32,n[1]*n[2]*samplingBinning);
Jacobians = Array{SparseMatrixCSC{Float64,Int32}}(ndips);

iifunc = (ii)->binfunc(ii,binningFactor)
nz = 1;
volCell = prod(Mesh.h);
for ii = 1:ndips
	u = vec(u);
	u,I1,J1,V1,Ihuge = ParamLevelSetModelFunc(Mesh,m;computeJacobian = 1,sigma = sigmaH,
				Xc = Xc,u = u,dsu = dsu,Ihuge = Ihuge,Is = I1, Js = J1,Vs = V1,iifunc = iifunc,numParamOfRBF=numParamOfRBF);
	u = reshape(u,tuple(n...));
	d[:,ii] = sampleSlices(u,volCell,samplingBinning);
	Jacobians[ii] = sparse(J1,I1,V1,lengthRBFparams,traceLength);
	(Jacobians[ii].nzval).*=volCell;
end
JacT = blockdiag(Jacobians...);
# JacT = 0.0;

return d,JacT
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
type MongoStore

	graph::Any

	uri::String

	function MongoStore(g::Graph, uri::String)
		# Construct an edge from a pair of vertices with a set of property values

		ms = new()
		ms.graph = g
		ms.uri = uri

		return ms
	end

	MongoStore(g::Graph, host::String, port::Int) = MongoStore(g, "mongodb://$host:$port/")
    MongoStore(g::Graph) = MongoStore(g, "localhost", 27017)
end

function buildcollectionkey(id::UUID.Uuid, prefix::String, suffix::String = "")
	return "$UUID.asString(id)_$(prefix)_$(suffix)"
end

function getorbuildcollection(client::MongoClient, collectiondict::Dict{String, MongoCollection}, g::Graph, collectionkey::String, collectiontype::String)
	if haskey(collectiondict, collectionkey)
		return collectiondict[collectionkey]
	else
		dbname = "@pg_$(UUID.asString(g.id))"
		collection = MongoCollection(client, dbname, collectionkey)

		# ensure indexes
		#ensureindex(collection, Set{String}("@pg_id"), Dict{String, String}("unique"=>"true"))
		#ensureindex(collection, Set{String}("@pg_typelabel"), Dict{String, String}("unique"=>"false"))

		#if collectiontype == "edge"
		#	ensureindex(collection, Set{String}("@pg_head_id"), Dict{String, String}("unique"=>"false"))
		#
		#	ensureindex(collection, Set{String}("@pg_head_typelabel"), Dict{String, String}("unique"=>"false"))
		#
		#		ensureindex(collection, Set{String}("@pg_tail_id"), Dict{String, String}("unique"=>"false"))
		#
		#	ensureindex(collection, Set{String}("@pg_tail_typelabel"), Dict{String, String}("unique"=>"false"))
		#end

		collectiondict[collectionkey] = collection
		return collection
	end
end

function persistnewvertices(client::MongoClient, newvertices::Dict{UUID.Uuid, Vertex}, graph::Graph, collectiondict::Dict{String, MongoCollection})

	result = PersistenceOperationResult{Vertex}()

	try
		for (k,v) in newvertices
			result.lastitemattempted = v

			collectionkey = buildcollectionkey(graph.id, "vertex", v.typelabel)
			collection = getorbuildcollection(client, collectiondict, graph, collectionkey, "vertex")

			oid = BSONOID()
			props = Dict{Any,Any}("_id"=>oid,
									 "@pg_id"=>UUID.asString(v.id),
									 "@pg_typelabel"=>v.typelabel)

			merge!(props, v.attachedproperties)
			insert(collection, props)

			push!(result.successfulset, v)
		end
	catch e
		result.errorencountered = true

		if isdefined(result, :lastitemattempted)
			push!(result.failedset, result.lastitemattempted)
		end

		showerror(STDOUT, e)
	end

	result.issuccessful = !result.errorencountered

	return result
end

function persistnewedges(client::MongoClient, newedges::Dict{UUID.Uuid, Edge}, graph::Graph, collectiondict::Dict{String, MongoCollection})

	result = PersistenceOperationResult{Edge}()

	try
		for (k,v) in newedges

			result.lastitemattempted = v

			collectionkey = buildcollectionkey(graph.id, "edge", v.typelabel)
			collection = getorbuildcollection(client, collectiondict, graph, collectionkey,"edge")
			oid = BSONOID()
			props = Dict{Any,Any}("_id"=>oid,
									 "@pg_id"=>UUID.asString(v.id),
									 "@pg_typelabel"=>v.typelabel,
									 "@pg_head_id"=>UUID.asString(v.head.id),
									 "@pg_head_typelabel"=>v.head.typelabel,
									 "@pg_tail_id"=>UUID.asString(v.tail.id),
									 "@pg_tail_typelabel"=>v.tail.typelabel)
			merge!(props, v.attachedproperties)

			insert(collection, props)

			push!(result.successfulset, v)
		end
	catch
		result.errorencountered = true

		if isdefined(result, :lastitemattempted)
			push!(result.failedset, result.lastitemattempted)
		end
	end

	result.issuccessful = !result.errorencountered

	return result
end

function persistchangededges(client::MongoClient, changededges::Dict{UUID.Uuid, Edge}, graph::Graph, collectiondict::Dict{String, MongoCollection})

	result = PersistenceOperationResult{Edge}()

	try
		for (k,v) in changededges

			result.lastitemattempted = v

			collectionkey = buildcollectionkey(graph.id, "edge")
			collection = getorbuildcollection(client, collectiondict, graph, collectionkey,"edge")

			oid = BSONOID()
			props = Dict{Any,Any}("_id"=>oid,
									 "@pg_id"=>UUID.asString(v.id),
									 "@pg_typelabel"=>v.typelabel,
									 "@pg_head_id"=>UUID.asString(v.head.id),
									 "@pg_head_typelabel"=>v.head.typelabel,
									 "@pg_tail_id"=>UUID.asString(v.tail.id),
									 "@pg_tail_typelabel"=>v.tail.typelabel)
			merge!(props, v.attachedproperties)

			update(collection,
				   {"pg_id"=>UUID.asString(v.id)},
				   {"\$set"=>props}
				   )

			push!(result.successfulset, v)
		end
	catch
		result.errorencountered = true

		if isdefined(result, :lastitemattempted)
			push!(result.failedset, result.lastitemattempted)
		end
	end

	result.issuccessful = !result.errorencountered

	return result
end

function persistchangedvertices(client::MongoClient, changedvertices::Dict{UUID.Uuid, Vertex}, graph::Graph, collectiondict::Dict{String, MongoCollection})
	result = PersistenceOperationResult{Vertex}()

	try
		for (k,v) in changedvertices

			result.lastitemattempted = v

			collectionkey = buildcollectionkey(graph.id, "vertex", v.typelabel)
			collection = getorbuildcollection(client, collectiondict, graph, collectionkey,"vertex")
			oid = BSONOID()
			props = Dict{Any,Any}("_id"=>oid, "@pg_id"=>UUID.asString(v.id), "@pg_typelabel"=>v.typelabel)

			merge!(props, v.attachedproperties)

			update(collection,
				   {"pg_id"=>UUID.asString(v.id)},
				   {"\$set"=>props}
				   )

			push!(result.successfulset, v)
		end
	catch
		result.errorencountered = true

		if isdefined(result, :lastitemattempted)
			push!(result.failedset, result.lastitemattempted)
		end
	end

	result.issuccessful = !result.errorencountered

	return result
end

function processremovededges(client::MongoClient, removededges::Dict{UUID.Uuid, Edge}, graph::Graph, collectiondict::Dict{String, MongoCollection})
	result = PersistenceOperationResult{Edge}()

	try
		for (k,v) in removededges

			result.lastitemattempted = v
			collectionkey = buildcollectionkey(graph.id, "edge")
			collection = getorbuildcollection(client, collectiondict, graph, collectionkey,"vertex")
			delete(
				collection,
				{"@pg_id"=>UUID.asString(v.id)}
			)

			push!(result.successfulset, v)
		end
	catch
		result.errorencountered = true

		if isdefined(result, :lastitemattempted)
			push!(result.failedset, result.lastitemattempted)
		end
	end

	result.issuccessful = !result.errorencountered

	return result
end

function processremovedvertices(client::MongoClient, removedvertices::Dict{UUID.Uuid, Vertex}, graph::Graph, collectiondict::Dict{String, MongoCollection})
	result = PersistenceOperationResult{Vertex}()

	try
		for (k,v) in removedvertices
			result.lastitemattempted = v
			collectionkey = buildcollectionkey(graph.id, "vertex")
			collection = getorbuildcollection(client, collectiondict, graph, collectionkey,"vertex")
			delete(
				collection,
				{"@pg_id"=>UUID.asString(v.id)}
			)

			push!(result.successfulset, v)
		end
	catch
		result.errorencountered = true

		if isdefined(result, :lastitemattempted)
			push!(result.failedset, result.lastitemattempted)
		end
	end

	result.issuccessful = !result.errorencountered

	return result
end

function persist(store::MongoStore)

	result = PersistenceBatchResult()

	if !isdefined(store.graph, :tracker)
		throw(ChangeTrackerNotDefinedException())
	end

	if !store.graph.tracker.ischanged
		# if nothing changed, nothing to persist
		return false
	end

	#client = MongoClient(store.uri)
	client = MongoClient()

	# a dictionary used to store references to MongoCollections by key
	collectiondict = Dict{String, MongoCollection}()

	try
		# process new vertices
		result.insertverticesresult = persistnewvertices(client, store.graph.tracker.newvertices, store.graph, collectiondict)

		# process new edges
		result.insertedgesresult = persistnewedges(client, store.graph.tracker.newedges, store.graph, collectiondict)

		# process changed vertices
		result.updateverticesresult = persistchangedvertices(client, store.graph.tracker.changedvertices, store.graph, collectiondict)

		# process changed edges
		result.updateedgesresult = persistchangededges(client, store.graph.tracker.changededges, store.graph, collectiondict)

		# process removed edges
		result.removeedgesresult = processremovededges(client, store.graph.tracker.removededges, store.graph, collectiondict)

		# process removed vertices
		result.removeverticesresult = processremovedvertices(client, store.graph.tracker.removedvertices, store.graph, collectiondict)

	finally
		updatebatchstatus!(result)

		if result.issuccessful
			# mark the graph as unchanged as the change has been processed
			clearchanges(store.graph.tracker)
		else
			# more selectively remove record of changes
			for item in result.insertverticesresult.successfulset
				clearnewvertex(store.graph.tracker, item)
			end
			for item in result.insertedgesresult.successfulset
				clearnewedge(store.graph.tracker, item)
			end
			for item in result.updateverticesresult.successfulset
				clearchangedvertex(store.graph.tracker, item)
			end
			for item in result.updateedgesresult.successfulset
				clearchangededge(store.graph.tracker, item)
			end
			for item in result.removeverticesresult.successfulset
				clearremovedvertex(store.graph.tracker, item)
			end
			for item in result.removeedgesresult.successfulset
				clearremovededge(store.graph.tracker, item)
			end
		end
	end

	return result
end