text:
# This code presents an example to illustrate how NLboxsolve.jl can be used

using NLboxsolve

function test_nlboxsolve_outplace()

    function test8(x)

        f = similar(x)

        f[1] = x[1]^2 - 1.0
        for i = 2:length(x)-1
            f[i] = x[i-1] - x[i]^3
        end
        f[end] = x[end-1] - x[end]

        return f

    end

    # Solve test8

    n = 20
    x0 = ones(n)*1.5
    l = ones(n)*0.5
    l[1] = 0.8
    u = ones(n)*2.0

    soln_newton     = constrained_newton(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm         = constrained_levenberg_marquardt(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm_kyf     = constrained_levenberg_marquardt_kyf(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm_ar      = constrained_levenberg_marquardt_ar(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_dogleg     = constrained_dogleg_solver(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_dogleg_bmp = constrained_dogleg_bmp_solver(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_nk         = constrained_newton_krylov(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_jfnk       = constrained_jacobian_free_newton_krylov(test8,x0,l,u,xtol=1e-15,ftol=1e-15)

    soln_newton_sparse     = constrained_newton_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_sparse         = constrained_levenberg_marquardt_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_kyf_sparse     = constrained_levenberg_marquardt_kyf_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_ar_sparse      = constrained_levenberg_marquardt_ar_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_sparse     = constrained_dogleg_solver_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_bmp_sparse = constrained_dogleg_solver_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_nk_sparse         = constrained_newton_krylov_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)

    soln1 = nlboxsolve(test8,x0,l,u,method = :dogleg,xtol=1e-15,ftol=1e-15)
    soln2 = nlboxsolve(test8,x0,l,u,xtol=1e-15,ftol=1e-15)

    test_one   = maximum(abs,soln_newton.zero       - soln2.zero) < 1e-14
    test_two   = maximum(abs,soln_lm.zero           - soln2.zero) < 1e-14
    test_three = maximum(abs,soln_lm_kyf.zero       - soln2.zero) < 1e-14
    test_four  = maximum(abs,soln_lm_ar.zero        - soln2.zero) < 1e-14
    test_five  = maximum(abs,soln_dogleg.zero       - soln2.zero) < 1e-14
    test_six   = maximum(abs,soln_dogleg_bmp.zero   - soln2.zero) < 1e-14
    test_seven = maximum(abs,soln_nk.zero           - soln2.zero) < 1e-14
    test_eight = maximum(abs,soln_jfnk.zero         - soln2.zero) < 1e-14

    test_nine     = maximum(abs,soln_newton_sparse.zero       - soln2.zero) < 1e-14
    test_ten      = maximum(abs,soln_lm_sparse.zero           - soln2.zero) < 1e-14
    test_eleven   = maximum(abs,soln_lm_kyf_sparse.zero       - soln2.zero) < 1e-14
    test_twelve   = maximum(abs,soln_lm_ar_sparse.zero        - soln2.zero) < 1e-14
    test_thirteen = maximum(abs,soln_dogleg_sparse.zero       - soln2.zero) < 1e-14
    test_fourteen = maximum(abs,soln_dogleg_bmp_sparse.zero   - soln2.zero) < 1e-14
    test_fifteen  = maximum(abs,soln_nk_sparse.zero           - soln2.zero) < 1e-14

    return test_one, test_two, test_three, test_four, test_five, test_six, test_seven, test_eight, test_nine, test_ten, test_eleven,
           test_twelve, test_thirteen, test_fourteen, test_fifteen

end

function test_nlboxsolve_inplace()

    function test8(f,x)

        f[1] = x[1]^2 - 1.0
        for i = 2:length(x)-1
            f[i] = x[i-1] - x[i]^3
        end
        f[end] = x[end-1] - x[end]

    end

    # Solve test8

    n = 20
    x0 = ones(n)*1.5
    l = ones(n)*0.5
    l[1] = 0.8
    u = ones(n)*2.0

    soln_newton     = constrained_newton(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm         = constrained_levenberg_marquardt(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm_kyf     = constrained_levenberg_marquardt_kyf(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_lm_ar      = constrained_levenberg_marquardt_ar(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_dogleg     = constrained_dogleg_solver(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_dogleg_bmp = constrained_dogleg_bmp_solver(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_nk         = constrained_newton_krylov(test8,x0,l,u,xtol=1e-15,ftol=1e-15)
    soln_jfnk       = constrained_jacobian_free_newton_krylov(test8,x0,l,u,xtol=1e-15,ftol=1e-15)

    soln_newton_sparse     = constrained_newton_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_sparse         = constrained_levenberg_marquardt_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_kyf_sparse     = constrained_levenberg_marquardt_kyf_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_lm_ar_sparse      = constrained_levenberg_marquardt_ar_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_sparse     = constrained_dogleg_solver_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_dogleg_bmp_sparse = constrained_dogleg_solver_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)
    soln_nk_sparse         = constrained_newton_krylov_sparse(test8,x0,l,u,ftol=1e-15,xtol=1e-15)

    soln1 = nlboxsolve(test8,x0,l,u,method = :dogleg,xtol=1e-15,ftol=1e-15)
    soln2 = nlboxsolve(test8,x0,l,u,xtol=1e-15,ftol=1e-15)

    test_one   = maximum(abs,soln_newton.zero       - soln2.zero) < 1e-14
    test_two   = maximum(abs,soln_lm.zero           - soln2.zero) < 1e-14
    test_three = maximum(abs,soln_lm_kyf.zero       - soln2.zero) < 1e-14
    test_four  = maximum(abs,soln_lm_ar.zero        - soln2.zero) < 1e-14
    test_five  = maximum(abs,soln_dogleg.zero       - soln2.zero) < 1e-14
    test_six   = maximum(abs,soln_dogleg_bmp.zero   - soln2.zero) < 1e-14
    test_seven = maximum(abs,soln_nk.zero           - soln2.zero) < 1e-14
    test_eight = maximum(abs,soln_jfnk.zero         - soln2.zero) < 1e-14

    test_nine     = maximum(abs,soln_newton_sparse.zero       - soln2.zero) < 1e-14
    test_ten      = maximum(abs,soln_lm_sparse.zero           - soln2.zero) < 1e-14
    test_eleven   = maximum(abs,soln_lm_kyf_sparse.zero       - soln2.zero) < 1e-14
    test_twelve   = maximum(abs,soln_lm_ar_sparse.zero        - soln2.zero) < 1e-14
    test_thirteen = maximum(abs,soln_dogleg_sparse.zero       - soln2.zero) < 1e-14
    test_fourteen = maximum(abs,soln_dogleg_bmp_sparse.zero   - soln2.zero) < 1e-14
    test_fifteen  = maximum(abs,soln_nk_sparse.zero           - soln2.zero) < 1e-14

    return test_one, test_two, test_three, test_four, test_five, test_six, test_seven, test_eight, test_nine, test_ten, test_eleven,
           test_twelve, test_thirteen, test_fourteen, test_fifteen

end

test_nlboxsolve_outplace()
test_nlboxsolve_inplace()


text:
module Airtable

export AirBase,
       AirTable,
       AirRecord,
       path

using Dates
using HTTP
using JSON3
using ProgressMeter
using ReTest

const API_VERSION = "v0"

include("ratelimiting.jl")
include("auth.jl")
include("api.jl")
include("interface.jl")

end

ext: jl


text:
# Persistent Vectors
# ==================

struct PersistentVector{T} <: AbstractArray{T,1}
    trie::DenseBitmappedTrie{Vector{T}}
    tail::Vector{T}
    length::Int
end
PersistentVector{T}() where {T} =
    PersistentVector{T}(DenseLeaf{Vector{T}}(), T[], 0)
function PersistentVector{T}(arr) where T
    if length(arr) <= trielen
        PersistentVector{T}(DenseLeaf{Vector{T}}(), arr, length(arr))
    else
        append(PersistentVector{T}(DenseLeaf{Vector{T}}(), T[], 0), arr)
    end
end
PersistentVector(v::PersistentVector) = v
PersistentVector{T}(v::PersistentVector{T}) where {T} = v
PersistentVector() = PersistentVector{Any}()
PersistentVector(itr) = PersistentVector{eltype(itr)}(itr)

Base.convert(::Type{PersistentVector}, x) = PersistentVector(x)
Base.convert(::Type{PersistentVector{T}}, x) where {T} = PersistentVector{T}(x)
Base.convert(::Type{PersistentVector}, x::PersistentVector{T}) where {T} = x
Base.convert(::Type{PersistentVector{T}}, x::PersistentVector{T}) where {T} = x

mask(i::Int) = ((i - 1) & (trielen - 1)) + 1

function boundscheck!(v::PersistentVector, i::Int)
    0 < i <= v.length || error(BoundsError(), " :: Index $i out of bounds ($(v.length))")
end

Base.size(   v::PersistentVector) = (v.length,)
Base.length( v::PersistentVector) = v.length
Base.isempty(v::PersistentVector) = length(v) == 0
Base.lastindex(  v::PersistentVector) = length(v)

Base.isequal(v1::PersistentVector, v2::PersistentVector) =
    isequal(v1.tail, v2.tail) && isequal(v1.trie, v2.trie)
==(v1::PersistentVector, v2::PersistentVector) =
    v1.tail == v2.tail && v1.trie == v2.trie

function Base.getindex(v::PersistentVector, i::Int)
    boundscheck!(v, i)
    if i > v.length - length(v.tail)
        v.tail[mask(i)]
    else
        v.trie[i][mask(i)]
    end
end

peek(v::PersistentVector) = v[end]

function push(v::PersistentVector{T}, el) where T
    if length(v.tail) < trielen
        newtail = copy_to_len(v.tail, 1 + length(v.tail))
        newtail[end] = el
        PersistentVector{T}(v.trie, newtail, 1 + v.length)
    else
        # T[el] will give an error when T is an tuple type in v0.3
        # workaround:
        arr = Array{T,1}(undef, 1)
        arr[1] = convert(T, el)
        PersistentVector{T}(append(v.trie, v.tail), arr, 1 + v.length)
    end
end
append(v::PersistentVector{T}, itr) where {T} = foldl(push, itr, init=v)

function assoc(v::PersistentVector{T}, i::Int, el) where T
    boundscheck!(v, i)
    if i > v.length - length(v.tail)
        newtail = v.tail[1:end]
        newtail[mask(i)] = el
        PersistentVector{T}(v.trie, newtail, v.length)
    else
        newnode = v.trie[i][1:end]
        newnode[mask(i)] = el
        PersistentVector{T}(assoc(v.trie, i, newnode), v.tail, v.length)
    end
end

function pop(v::PersistentVector{T}) where T
    if isempty(v.tail)
        newtail = peek(v.trie)[1:end-1]
        PersistentVector{T}(pop(v.trie), newtail, v.length - 1)
    else
        newtail = v.tail[1:end-1]
        PersistentVector{T}(v.trie, newtail, v.length - 1)
    end
end

struct ItrState{T}
    index::Int
    leaf::Vector{T}
end

function Base.iterate(v::PersistentVector, state = ItrState(1, v.length <= 32 ? v.tail : v.trie[1]))
    if state.index > v.length
        return nothing
    else
        i, leaf = state.index, state.leaf
        m = mask(i)
        value = leaf[m]
        i += 1
        if m == 32
            leaf = i > v.length - length(v.tail) ? v.tail : v.trie[i]
        end
        return value, ItrState(i, leaf)
    end
end

function Base.map(f::Function, pv::PersistentVector{T}) where T
    if length(pv) == 0 return PersistentVector{T}() end

    first = f(pv[1])
    v = PersistentVector{typeof(first)}([first])
    for i = 2:length(pv)
        v = push(v, f(pv[i]))
    end
    v
end

function Base.filter(f::Function, pv::PersistentVector{T}) where T
    v = PersistentVector{T}()
    for el in pv
        if f(el)
            v = push(v, el)
        end
    end
    v
end

function Base.hash(pv::PersistentVector{T}) where T
    h = hash(length(pv))
    for el in pv
        h = Base.hash(el, h)
    end
     UInt(h)
end

function print_elements(io, pv, range)
    for i=range
        print(io, "$(pv[i]), ")
    end
end

function print_vec(io::IO, t, head::String)
    print(io, "$head[")
    if length(t) < 50
        print_elements(io, t, 1:length(t)-1)
    else
        print_elements(io, t, 1:25)
        print(io, "..., ")
        print_elements(io, t, length(t)-25:length(t)-1)
    end
    if length(t) >= 1
        print(io, "$(t[end])]")
    else
        print(io, "]")
    end
end

 Base.show(io::IO, ::MIME"text/plain", pv::PersistentVector{T}) where {T} =
    print_vec(io, pv, "Persistent{$T}")


text:
int_rules_8_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 8.3*Exponential*integral*functions =#
    @apply_utils Antiderivative(ExpIntegralE(~n, ~(a') + ~(b') * ~x), ~x) => -(ExpIntegralE(~n + 1, ~a + ~b * ~x)) / ~b <-- FreeQ([~a, ~b, ~n], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * ExpIntegralE(~n, ~(b') * ~x), ~x) => (-((~x) ^ ~m) * ExpIntegralE(~n + 1, ~b * ~x)) / ~b + (~m / ~b) * Antiderivative((~x) ^ (~m - 1) * ExpIntegralE(~n + 1, ~b * ~x), ~x) <-- FreeQ(~b, ~x) && (EqQ(~m + ~n, 0) && IGtQ(~m, 0))
    @apply_utils Antiderivative(ExpIntegralE(1, ~(b') * ~x) / ~x, ~x) => (~b * ~x * HypergeometricPFQ([1, 1, 1], [2, 2, 2], -(~b) * ~x) - EulerGamma * log(~x)) - (1 / 2) * log(~b * ~x) ^ 2 <-- FreeQ(~b, ~x)
    @apply_utils Antiderivative((~x) ^ ~m * ExpIntegralE(~n, ~(b') * ~x), ~x) => ((~x) ^ (~m + 1) * ExpIntegralE(~n, ~b * ~x)) / (~m + 1) + (~b / (~m + 1)) * Antiderivative((~x) ^ (~m + 1) * ExpIntegralE(~n - 1, ~b * ~x), ~x) <-- FreeQ(~b, ~x) && (EqQ(~m + ~n, 0) && ILtQ(~m, -1))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~m * ExpIntegralE(~n, ~(b') * ~x), ~x) => ((~d * ~x) ^ ~m * Gamma(~m + 1) * log(~x)) / (~b * (~b * ~x) ^ ~m) - ((~d * ~x) ^ (~m + 1) * HypergeometricPFQ([~m + 1, ~m + 1], [~m + 2, ~m + 2], -(~b) * ~x)) / (~d * (~m + 1) ^ 2) <-- FreeQ([~b, ~d, ~m, ~n], ~x) && (EqQ(~m + ~n, 0) && Not(IntegerQ(~m)))
    @apply_utils Antiderivative((~(d') * ~x) ^ ~(m') * ExpIntegralE(~n, ~(b') * ~x), ~x) => ((~d * ~x) ^ (~m + 1) * ExpIntegralE(~n, ~b * ~x)) / (~d * (~m + ~n)) - ((~d * ~x) ^ (~m + 1) * ExpIntegralE(-(~m), ~b * ~x)) / (~d * (~m + ~n)) <-- FreeQ([~b, ~d, ~m, ~n], ~x) && NeQ(~m + ~n, 0)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * ExpIntegralE(~n, ~a + ~(b') * ~x), ~x) => (-((~c + ~d * ~x) ^ ~m) * ExpIntegralE(~n + 1, ~a + ~b * ~x)) / ~b + ((~d * ~m) / ~b) * Antiderivative((~c + ~d * ~x) ^ (~m - 1) * ExpIntegralE(~n + 1, ~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && (IGtQ(~m, 0) || (ILtQ(~n, 0) || GtQ(~m, 0) && LtQ(~n, -1)))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * ExpIntegralE(~n, ~a + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * ExpIntegralE(~n, ~a + ~b * ~x)) / (~d * (~m + 1)) + (~b / (~d * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m + 1) * ExpIntegralE(~n - 1, ~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x) && ((IGtQ(~n, 0) || LtQ(~m, -1) && GtQ(~n, 0)) && NeQ(~m, -1))
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * ExpIntegralE(~n, ~a + ~(b') * ~x), ~x) => Unintegrable((~c + ~d * ~x) ^ ~m * ExpIntegralE(~n, ~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m, ~n], ~x)
    @apply_utils Antiderivative(ExpIntegralEi(~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) * ExpIntegralEi(~a + ~b * ~x)) / ~b - ℯ ^ (~a + ~b * ~x) / ~b <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(ExpIntegralEi(~(b') * ~x) / ~x, ~x) => log(~x) * (ExpIntegralEi(~b * ~x) + ExpIntegralE(1, -(~b) * ~x)) - Antiderivative(ExpIntegralE(1, -(~b) * ~x) / ~x, ~x) <-- FreeQ(~b, ~x)
    @apply_utils Antiderivative(ExpIntegralEi(~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x), ~x) => Unintegrable(ExpIntegralEi(~a + ~b * ~x) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * ExpIntegralEi(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * ExpIntegralEi(~a + ~b * ~x)) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative(((~c + ~d * ~x) ^ (~m + 1) * ℯ ^ (~a + ~b * ~x)) / (~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative(ExpIntegralEi(~(a') + ~(b') * ~x) ^ 2, ~x) => ((~a + ~b * ~x) * ExpIntegralEi(~a + ~b * ~x) ^ 2) / ~b - 2 * Antiderivative(ℯ ^ (~a + ~b * ~x) * ExpIntegralEi(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * ExpIntegralEi(~(b') * ~x) ^ 2, ~x) => ((~x) ^ (~m + 1) * ExpIntegralEi(~b * ~x) ^ 2) / (~m + 1) - (2 / (~m + 1)) * Antiderivative((~x) ^ ~m * ℯ ^ (~b * ~x) * ExpIntegralEi(~b * ~x), ~x) <-- FreeQ(~b, ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~x) ^ ~(m') * ExpIntegralEi(~a + ~(b') * ~x) ^ 2, ~x) => ((((~x) ^ (~m + 1) * ExpIntegralEi(~a + ~b * ~x) ^ 2) / (~m + 1) + (~a * (~x) ^ ~m * ExpIntegralEi(~a + ~b * ~x) ^ 2) / (~b * (~m + 1))) - (2 / (~m + 1)) * Antiderivative((~x) ^ ~m * ℯ ^ (~a + ~b * ~x) * ExpIntegralEi(~a + ~b * ~x), ~x)) - ((~a * ~m) / (~b * (~m + 1))) * Antiderivative((~x) ^ (~m - 1) * ExpIntegralEi(~a + ~b * ~x) ^ 2, ~x) <-- FreeQ([~a, ~b], ~x) && IGtQ(~m, 0)
    #= Antiderivative((~x)^(~m')*ExpIntegralEi((~a)+(~b')*(~x))^2,~x) :=  b*x^(m+2)*ExpIntegralEi(a+b*x)^2/(a*(m+1)) + x^(m+1)*ExpIntegralEi(a+b*x)^2/(m+1) - 2*b/(a*(m+1))*Antiderivative(x^(m+1)*ℯ^(a+b*x)*ExpIntegralEi(a+b*x),x) - b*(m+2)/(a*(m+1))*Antiderivative(x^(m+1)*ExpIntegralEi(a+b*x)^2,x) <--  FreeQ([a,b],x) && ILtQ(m,-2) =#
    @apply_utils Antiderivative(ℯ ^ (~(a') + ~(b') * ~x) * ExpIntegralEi(~(c') + ~(d') * ~x), ~x) => (ℯ ^ (~a + ~b * ~x) * ExpIntegralEi(~c + ~d * ~x)) / ~b - (~d / ~b) * Antiderivative(ℯ ^ (~a + ~c + (~b + ~d) * ~x) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative((~x) ^ ~(m') * ℯ ^ (~(a') + ~(b') * ~x) * ExpIntegralEi(~(c') + ~(d') * ~x), ~x) => (((~x) ^ ~m * ℯ ^ (~a + ~b * ~x) * ExpIntegralEi(~c + ~d * ~x)) / ~b - (~d / ~b) * Antiderivative(((~x) ^ ~m * ℯ ^ (~a + ~c + (~b + ~d) * ~x)) / (~c + ~d * ~x), ~x)) - (~m / ~b) * Antiderivative((~x) ^ (~m - 1) * ℯ ^ (~a + ~b * ~x) * ExpIntegralEi(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && IGtQ(~m, 0)
    @apply_utils Antiderivative((~x) ^ ~m * ℯ ^ (~(a') + ~(b') * ~x) * ExpIntegralEi(~(c') + ~(d') * ~x), ~x) => (((~x) ^ (~m + 1) * ℯ ^ (~a + ~b * ~x) * ExpIntegralEi(~c + ~d * ~x)) / (~m + 1) - (~d / (~m + 1)) * Antiderivative(((~x) ^ (~m + 1) * ℯ ^ (~a + ~c + (~b + ~d) * ~x)) / (~c + ~d * ~x), ~x)) - (~b / (~m + 1)) * Antiderivative((~x) ^ (~m + 1) * ℯ ^ (~a + ~b * ~x) * ExpIntegralEi(~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x) && ILtQ(~m, -1)
    @apply_utils Antiderivative(ExpIntegralEi(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ~x * ExpIntegralEi(~d * (~a + ~b * log(~c * (~x) ^ ~n))) - ~b * ~n * ℯ ^ (~a * ~d) * Antiderivative((~c * (~x) ^ ~n) ^ (~b * ~d) / (~a + ~b * log(~c * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative(ExpIntegralEi(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))) / ~x, ~x) => (1 / ~n) * Subst(ExpIntegralEi(~d * (~a + ~b * ~x)), ~x, log(~c * (~x) ^ ~n)) <-- FreeQ([~a, ~b, ~c, ~d, ~n], ~x)
    @apply_utils Antiderivative((~(e') * ~x) ^ ~(m') * ExpIntegralEi(~(d') * (~(a') + ~(b') * log(~(c') * (~x) ^ ~(n')))), ~x) => ((~e * ~x) ^ (~m + 1) * ExpIntegralEi(~d * (~a + ~b * log(~c * (~x) ^ ~n)))) / (~e * (~m + 1)) - ((~b * ~n * ℯ ^ (~a * ~d) * (~c * (~x) ^ ~n) ^ (~b * ~d)) / ((~m + 1) * (~e * ~x) ^ (~b * ~d * ~n))) * Antiderivative((~e * ~x) ^ (~m + ~b * ~d * ~n) / (~a + ~b * log(~c * (~x) ^ ~n)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m, ~n], ~x) && NeQ(~m, -1)
    @apply_utils Antiderivative(LogIntegral(~(a') + ~(b') * ~x), ~x) => ((~a + ~b * ~x) * LogIntegral(~a + ~b * ~x)) / ~b - ExpIntegralEi(2 * log(~a + ~b * ~x)) / ~b <-- FreeQ([~a, ~b], ~x)
    @apply_utils Antiderivative(LogIntegral(~(b') * ~x) / ~x, ~x) => -(~b) * ~x + log(~b * ~x) * LogIntegral(~b * ~x) <-- FreeQ(~b, ~x)
    @apply_utils Antiderivative(LogIntegral(~(a') + ~(b') * ~x) / (~(c') + ~(d') * ~x), ~x) => Unintegrable(LogIntegral(~a + ~b * ~x) / (~c + ~d * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d], ~x)
    @apply_utils Antiderivative((~(c') + ~(d') * ~x) ^ ~(m') * LogIntegral(~(a') + ~(b') * ~x), ~x) => ((~c + ~d * ~x) ^ (~m + 1) * LogIntegral(~a + ~b * ~x)) / (~d * (~m + 1)) - (~b / (~d * (~m + 1))) * Antiderivative((~c + ~d * ~x) ^ (~m + 1) / log(~a + ~b * ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~m], ~x) && NeQ(~m, -1)
end

max_stars_repo_path: src/rules/8 Special functions/8.3 Exponential integral functions.jl


text:
#=
Thermodynamic relationships
=#

####### BASIC PERFECT GAS EQUATIONS OF STATE #######

Temperature(ρ::Density,p::Pressure;gas::PerfectGas=DefaultPerfectGas) = Temperature(p/(ρ*GasConstant(gas)))
Temperature(p::Pressure,ρ::Density;gas::PerfectGas=DefaultPerfectGas) = Temperature(ρ,p,gas=gas)

Density(p::Pressure,T::Temperature;gas::PerfectGas=DefaultPerfectGas) = Density(p/(GasConstant(gas)*T))
Density(T::Temperature,p::Pressure;gas::PerfectGas=DefaultPerfectGas) = Density(p,T,gas=gas)

Pressure(ρ::Density,T::Temperature;gas::PerfectGas=DefaultPerfectGas) = Pressure(ρ*GasConstant(gas)*T)
Pressure(T::Temperature,ρ::Density;gas::PerfectGas=DefaultPerfectGas) = Pressure(ρ,T,gas=gas)

StagnationTemperature(ρ0::StagnationDensity,p0::StagnationPressure;gas::PerfectGas=DefaultPerfectGas) = StagnationTemperature(p0/(ρ0*GasConstant(gas)))
StagnationTemperature(p0::StagnationPressure,ρ0::StagnationDensity;gas::PerfectGas=DefaultPerfectGas) = StagnationTemperature(ρ0,p0,gas=gas)

StagnationDensity(p0::StagnationPressure,T0::StagnationTemperature;gas::PerfectGas=DefaultPerfectGas) = StagnationDensity(p0/(GasConstant(gas)*T0))
StagnationDensity(T0::StagnationTemperature,p0::StagnationPressure;gas::PerfectGas=DefaultPerfectGas) = StagnationDensity(p0,T0,gas=gas)

StagnationPressure(ρ0::StagnationDensity,T0::StagnationTemperature;gas::PerfectGas=DefaultPerfectGas) = StagnationPressure(ρ0*GasConstant(gas)*T0)
StagnationPressure(T0::StagnationTemperature,ρ0::StagnationDensity;gas::PerfectGas=DefaultPerfectGas) = StagnationPressure(ρ0,T0,gas=gas)

"""
    SoundSpeed(T::Temperature[,γ::SpecificHeatRatio=1.4][,R::GasConstant=287])

Compute the speed of sound of a perfect gas, based on the absolute temperature `T`. The ratio of specific heats `γ` and `R`
can also be supplied, but default to the values for air at standard conditions
(1.4 and 287 J/kg.K, respectively).
"""
SoundSpeed(T::Temperature;gas::PerfectGas=DefaultPerfectGas) = SoundSpeed(sqrt(SpecificHeatRatio(gas)*GasConstant(gas)*T))

Enthalpy(T::Temperature;gas::PerfectGas=DefaultPerfectGas) = Enthalpy(SpecificHeatPressure(gas)*T)
Temperature(h::Enthalpy;gas::PerfectGas=DefaultPerfectGas) = Temperature(h/SpecificHeatPressure(gas))

StagnationEnthalpy(T0::StagnationTemperature;gas::PerfectGas=DefaultPerfectGas) = StagnationEnthalpy(SpecificHeatPressure(gas)*T0)
StagnationTemperature(h0::StagnationEnthalpy;gas::PerfectGas=DefaultPerfectGas) = StagnationTemperature(h0/SpecificHeatPressure(gas))

InternalEnergy(T::Temperature;gas::PerfectGas=DefaultPerfectGas) = InternalEnergy(SpecificHeatVolume(gas)*T)
Temperature(e::InternalEnergy;gas::PerfectGas=DefaultPerfectGas) = Temperature(e/SpecificHeatVolume(gas))

StagnationInternalEnergy(T0::StagnationTemperature;gas::PerfectGas=DefaultPerfectGas) = StagnationInternalEnergy(SpecificHeatVolume(gas)*T0)
StagnationTemperature(e0::StagnationInternalEnergy;gas::PerfectGas=DefaultPerfectGas) = StagnationTemperature(e0/SpecificHeatVolume(gas))

StagnationEnthalpy(h::Enthalpy,u::Velocity) = StagnationEnthalpy(value(h)+0.5*value(u)^2)
Enthalpy(h0::StagnationEnthalpy,u::Velocity) = Enthalpy(value(h0)-0.5*value(u)^2)

# Mach number - velocity relations
Velocity(a::SoundSpeed,M::MachNumber) = Velocity(M*a)
MachNumber(u::Velocity,a::SoundSpeed) = MachNumber(u/a)
SoundSpeed(u::Velocity,M::MachNumber) = SoundSpeed(u/M)

# mass flow rate
MassFlowRate(ρ::Density,u::Velocity,A::Area) = MassFlowRate(ρ*u*A)


text:
function imshow(m::MapImages.MapImage, args...; kwargs...)
    return plt.imshow(m, extent=MapImages.grid_extent(m), args...; kwargs...)
end

function imshow(ax, m::MapImages.MapImage, args...; kwargs...)
    return ax.imshow(m, extent=MapImages.grid_extent(m), args...; kwargs...)
end


text:
"""
    histogram(data)

Creates an `EChart` where region below plotted line filled with color.

## Methods
```julia
histogram(h::StatsBase.Histogram)
```

## Arguments
* `legend::Bool = false` : display legend?
* `horizontal::Bool = false` : show bars horizontally?
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

Reasonable defaults set for different methods of `area`, such as displaying a legend when two or more series present.

## Examples
```julia
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
ar = area(x, y)
```
"""
function histogram(h::StatsBase.Histogram;
                   legend::Bool = false,
                   horizontal::Bool = false,
                   kwargs...)

    if length(h.edges) == 1
        #Create labels for bins
        edges = collect(h.edges[1])
        h.closed == :left ? left = "[" : left = "("
        left == "[" ? right = ")" : right = "]"
        labels = []
        for i in 1:length(edges)-1
            push!(labels, string(left, edges[i], " - ", edges[i + 1], right))
        end

        #Plot as bar, move gap to be near inperceptible
        ec = bar(labels, h.weights; kwargs...)
        ec.series[1].barCategoryGap = "1%"

        #Make plot horizontal
        horizontal ? flip!(ec) : nothing

        #Add legend if requested
        legend ? legend!(ec) : nothing

        return ec

    elseif length(h.edges) == 2
        ec = heatmap(h)
        return ec

    else
        DimensionMismatch("No more than two dimensions can be plotted")
    end

end


text:
import .API: BPFProgType

const bpf_object = Ptr{Cvoid}
struct Object
    obj::bpf_object
end
const bpf_program = Ptr{Cvoid}
struct Program
    prog::bpf_program
end
const bpf_link = Ptr{Cvoid}
struct Link
    link::bpf_link
end
const bpf_map = Ptr{Cvoid}
struct Map
    map::bpf_map
end
const perf_buffer = Ptr{Cvoid}
struct PerfBuffer
    buf::perf_buffer
end

@enum libbpf_print_level begin
    LIBBPF_WARN
    LIBBPF_INFO
    LIBBPF_DEBUG
end

bpfprintfn(fn) = ccall((:libbpf_set_print, libbpf), Ptr{Cvoid}, (Ptr{Cvoid},), fn)

Base.@kwdef struct bpf_object_open_opts
    sz::Csize_t = sizeof(bpf_object_open_opts)
    object_name::Cstring = C_NULL
    relaxed_maps::Bool = false
    relaxed_core_relocs::Bool = false
    pin_root_path::Cstring = C_NULL
    attach_prog_fd::UInt32 = 0
    kconfig::Cstring = C_NULL
end

struct bpf_object_load_attr
    obj::bpf_object
    log_level::Cint
    target_btf_path::Cstring
end

function Object(path::String)
    obj = ccall((:bpf_object__open, libbpf), bpf_object, (Cstring,), path)
    @assert Int(obj) > 0
    Object(obj)
end
function Object(data::Vector{UInt8}; kwargs...)
    opts = Ref(bpf_object_open_opts(; kwargs...))
    obj = GC.@preserve opts begin
    ccall((:bpf_object__open_mem, libbpf), bpf_object, (Ptr{UInt8},Csize_t,Ptr{bpf_object_open_opts}), data, length(data), opts)
    end
    @assert Int(obj) > 0
    Object(obj)
end
function Base.close(obj::Object)
    ccall((:bpf_object__close, libbpf), Cvoid, (bpf_object,), obj[])
    nothing
end
function kernel_version(obj::Object)
    raw = ccall((:bpf_object__kversion, libbpf), Cuint, (bpf_object,), obj[])
    VersionNumber(raw & (0xFF << 16), raw & (0xFF << 8), raw * 0xFF)
end
function name(obj::Object)
    name = ccall((:bpf_object__name, libbpf), Cstring, (bpf_object,), obj[])
    @assert Int(pointer(name)) > 0
    unsafe_string(name)
end
function load(obj::Object; log_level=LIBBPF_WARN, target_btf_path=C_NULL)
    opts = Ref(bpf_object_load_attr(obj[], Cint(log_level), target_btf_path))
    res = GC.@preserve opts begin
    ccall((:bpf_object__load_xattr, libbpf), Cint, (Ptr{bpf_object_load_attr},), opts)
    end
    @assert res == 0
    nothing
end
function unload(obj::Object)
    res = ccall((:bpf_object__unload, libbpf), Cint, (bpf_object,), obj[])
    @assert res == 0
    nothing
end

## maps

struct MapDef
    type::Cuint
    key_size::Cuint
    value_size::Cuint
    max_entries::Cuint
    map_flags::Cuint
end

function findmap(obj::Object, name::String)
    map = ccall((:bpf_object__find_map_by_name, libbpf), bpf_map, (bpf_object,Cstring), obj[], name)
    map == C_NULL && error("Map $name not found")
    @assert Int(map) > 0
    Map(map)
end
function maps(obj::Object)
    maps = Map[]
    map = C_NULL
    while true
        map = ccall((:bpf_map__next, libbpf), bpf_map, (bpf_map, bpf_object), map, obj[])
        map == C_NULL && break
        push!(maps, Map(map))
    end
    maps
end
function name(map::Map)
    name = ccall((:bpf_map__name, libbpf), Cstring, (bpf_map,), map[])
    @assert Int(pointer(name)) > 0
    unsafe_string(name)
end
function fd(obj::Object, name::String)
    fd = ccall((:bpf_object__find_map_fd_by_name, libbpf), Cint, (bpf_object,Cstring), obj[], name)
    @assert fd != 0
    fd
end
function fd(map::Map)
    fd = ccall((:bpf_map__fd, libbpf), Cint, (bpf_map,), map[])
    @assert fd != 0
    fd
end
function definition(map::Map)
    def = ccall((:bpf_map__def, libbpf), Ptr{MapDef}, (bpf_map,), map[])
    @assert Int(def) > 0
    unsafe_load(def)
end
function resize!(map::Map, sz::Integer)
    ret = ccall((:bpf_map__resize, libbpf), Cint, (bpf_map,UInt32), map[], UInt32(sz))
    @assert ret == 0
    nothing
end
function reuse_fd(map::Map, fd::Integer)
    ret = ccall((:bpf_map_reuse_fd, libbpf), Cint, (bpf_map,Cint), map[], Cint(fd))
    @assert ret == 0
    nothing
end

## programs

function programs(obj::Object)
    progs = Program[]
    prog = C_NULL
    while true
        prog = ccall((:bpf_program__next, libbpf), bpf_program, (bpf_program, bpf_object), prog, obj[])
        prog == C_NULL && break
        push!(progs, Program(prog))
    end
    progs
end
function load(prog::Program, license::String, kern_version::UInt32)
    ret = ccall((:bpf_program__load, libbpf), Cint, (bpf_program,Cstring,UInt32), prog[], license, kern_version)
    @assert ret == 0
    nothing
end
function name(prog::Program)
    name = ccall((:bpf_program__name, libbpf), Cstring, (bpf_program,), prog[])
    @assert reinterpret(Int, name) > 0
    unsafe_string(name)
end
function title(prog::Program, needs_copy::Bool=false)
    title = ccall((:bpf_program__title, libbpf), Cstring, (bpf_program,Bool), prog[], needs_copy)
    @assert reinterpret(Int, title) > 0
    unsafe_string(title)
end
function fd(prog::Program)
    ccall((:bpf_program__fd, libbpf), Cint, (bpf_program,), prog[])
end

# program types

function type(prog::Program)
    BPFProgType(ccall((:bpf_program__get_type, libbpf), Cint, (bpf_program,), prog[]))
end
function set_socket_filter!(prog::Program)
    @assert ccall((:bpf_program__set_socket_filter, libbpf), Cint, (bpf_program,), prog[]) == 0;
end
function set_xdp!(prog::Program)
    @assert ccall((:bpf_program__set_xdp, libbpf), Cint, (bpf_program,), prog[]) == 0;
end
function set_perf_event!(prog::Program)
    @assert ccall((:bpf_program__set_perf_event, libbpf), Cint, (bpf_program,), prog[]) == 0;
end
function set_tracepoint!(prog::Program)
    @assert ccall((:bpf_program__set_tracepoint, libbpf), Cint, (bpf_program,), prog[]) == 0;
end
function set_kprobe!(prog::Program)
    @assert ccall((:bpf_program__set_kprobe, libbpf), Cint, (bpf_program,), prog[]) == 0;
end

# program attaching

function attach_perf_event!(prog::Program, fd::Cint)
    link = ccall((:bpf_program__attach_perf_event, libbpf), bpf_link, (bpf_program, Cint), prog[], fd)
    @assert Int(link) > 0
    Link(link)
end
function attach_kprobe!(prog::Program, retprobe::Bool, name::String)
    link = ccall((:bpf_program__attach_kprobe, libbpf), bpf_link, (bpf_program, Bool, Cstring), prog[], retprobe, name)
    @assert Int(link) > 0
    Link(link)
end
function attach_uprobe!(prog::Program, retprobe::Bool, pid::Cuint, name::String, func_offset::Csize_t)
    link = ccall((:bpf_program__attach_uprobe, libbpf), bpf_link, (bpf_program, Bool, Cuint, Cstring, Csize_t), prog[], retprobe, pid, name, func_offset)
    @assert Int(link) > 0
    Link(link)
end
function attach_tracepoint!(prog::Program, category::String, name::String)
    link = ccall((:bpf_program__attach_tracepoint, libbpf), bpf_link, (bpf_program, Cstring, Cstring), prog[], category, name)
    @assert Int(link) > 0
    Link(link)
end

# perf buffers

Base.@kwdef struct perf_buffer_opts
    sample_cb::Ptr{Cvoid} = C_NULL
    lost_cb::Ptr{Cvoid} = C_NULL
    ctx::Ptr{Cvoid} = C_NULL
end

function PerfBuffer(map_fd::Cint, page_cnt::Integer, opts::perf_buffer_opts=perf_buffer_opts())
    opts_ref = Ref{perf_buffer_opts}(opts)
    buf = GC.@preserve opts_ref begin
        opts_ptr = Base.unsafe_convert(Ptr{perf_buffer_opts}, opts_ref)
        ccall((:perf_buffer__new, libbpf), perf_buffer, (Cint,Csize_t,Ptr{perf_buffer_opts}), map_fd, Csize_t(page_cnt), opts_ptr)
    end
    @assert Int(buf) > 0
    PerfBuffer(buf)
end
function poll(buf::PerfBuffer, timeout_ms::Integer)
    ccall((:perf_buffer__poll, libbpf), Cint, (perf_buffer,Cint), buf[], Cint(timeout_ms))
end

## user-facing API

for (T,field) in ((Object,:obj), (Program,:prog), (Map,:map), (Link,:link))
    #= TODO
    @eval @inline Base.getproperty(x::$T, sym::Symbol) =
                      Base.getproperty(x, Val(sym))
    @eval @inline Base.setproperty!(x::$T, val, sym::Symbol) =
                      Base.setproperty!(x, val, Val(sym))
    =#
    @eval @inline Base.getindex(x::$T) = x.$field
end

function load(f, x)
    load(x)
    try
        f()
    finally
        unload(x)
    end
end


text:
module TestVariableFree

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_modification_multirow_vectoraffine_nonpos()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Variable.Free{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.5, 0.0]),
        ),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            MOI.OPTIMAL,
            (MOI.FEASIBLE_POINT, [0.25, 0.0]),
        ),
    )
    MOI.Test.test_modification_multirow_vectoraffine_nonpos(
        bridged_mock,
        MOI.Test.Config(),
    )
    MOI.set(
        bridged_mock,
        MOI.ConstraintName(),
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorAffineFunction{Float64},
                MOI.Nonpositives,
            }(),
        ),
        ["c"],
    )
    var_names = ["xpos", "xneg"]
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    MOI.set(
        bridged_mock,
        MOI.VariableName(),
        MOI.get(bridged_mock, MOI.ListOfVariableIndices()),
        ["x"],
    )
    MOI.set(
        bridged_mock,
        MOI.ConstraintName(),
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.Nonnegatives,
            }(),
        ),
        ["nonneg"],
    )
    s = """
    variables: xpos, xneg
    nonneg: [xpos, xneg] in VecMathOptInterface.Nonnegatives(2)
    c: [4.0xpos + -4.0xneg + -1.0, 3.0xpos + -3.0xneg + -1.0] in VecMathOptInterface.Nonpositives(2)
    maxobjective: xpos + -1.0xneg
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(mock, model, var_names, ["nonneg", "c"])
    s = """
    variables: x
    c: [4.0x + -1.0, 3.0x + -1.0] in VecMathOptInterface.Nonpositives(2)
    maxobjective: 1.0x
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(bridged_mock, model, ["x"], ["c"])
    return
end

function test_linear_LessThan_and_GreaterThan()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Variable.Free{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [0, 0, 0, 0]),
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [100, 0, 0, 0]),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [100, 0, 0, 100],
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) =>
                [1.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) =>
                [-1.0],
        ),
    )
    MOI.Test.test_linear_LessThan_and_GreaterThan(
        bridged_mock,
        MOI.Test.Config(),
    )

    loc = MOI.get(bridged_mock, MOI.ListOfConstraintTypesPresent())
    @test length(loc) == 2
    @test !((MOI.VectorOfVariables, MOI.Reals) in loc)
    @test (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) in loc
    @test (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) in loc
    @test MOI.get(mock, MOI.NumberOfVariables()) == 4
    @test MOI.get(bridged_mock, MOI.NumberOfVariables()) == 2
    vis = MOI.get(bridged_mock, MOI.ListOfVariableIndices())
    @test vis == MOI.VariableIndex.([-1, -2])

    cx = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Reals}(vis[1].value)
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), cx) == [100.0]
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), cx) == [0.0]
    cy = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Reals}(vis[2].value)
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), cy) == [-100.0]
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), cy) == [0.0]

    _test_delete_bridged_variable(
        bridged_mock,
        vis[1],
        MOI.Reals,
        2,
        (
            (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
            (MOI.VectorOfVariables, MOI.Nonpositives, 0),
        ),
    )
    _test_delete_bridged_variable(
        bridged_mock,
        vis[2],
        MOI.Reals,
        1,
        (
            (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
            (MOI.VectorOfVariables, MOI.Nonpositives, 0),
        ),
    )
    return
end

function test_linear_integration_modification()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Variable.Free{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [0, 0, 0, 0]),
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [100, 0, 0, 0]),
        (mock::MOI.Utilities.MockOptimizer) -> MOI.Utilities.mock_optimize!(
            mock,
            [100, 0, 0, 100],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives) =>
                [[1.0]],
            (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) =>
                [[-1.0]],
        ),
    )
    MOI.Test.test_linear_VectorAffineFunction(bridged_mock, MOI.Test.Config())

    x, y = MOI.get(bridged_mock, MOI.ListOfVariableIndices())

    cx = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Reals}(x.value)
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), cx) == [100.0]
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), cx) == [0.0]
    cy = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Reals}(y.value)
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), cy) == [-100.0]
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), cy) == [0.0]

    @test MOI.supports(
        bridged_mock,
        MOI.VariablePrimalStart(),
        MOI.VariableIndex,
    )
    @test MOI.supports(
        bridged_mock,
        MOI.VariablePrimalStart(),
        typeof(MOI.Bridges.bridge(bridged_mock, x)),
    )
    MOI.set(bridged_mock, MOI.VariablePrimalStart(), [x, y], [1.0, -1.0])
    xa, xb, ya, yb = MOI.get(mock, MOI.ListOfVariableIndices())
    @test MOI.get(mock, MOI.VariablePrimalStart(), [xa, xb, ya, yb]) ==
          [1.0, 0.0, 0.0, 1.0]
    @test MOI.get(bridged_mock, MOI.VariablePrimalStart(), x) == 1
    @test MOI.get(bridged_mock, MOI.VariablePrimalStart(), y) == -1
    return
end

function test_linear_transform()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    bridged_mock = MOI.Bridges.Variable.Free{Float64}(mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [1.0, 1.0, 0.0, 0.0]),
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.Utilities.mock_optimize!(mock, [0.5, 0.5, 0.0, 0.0]),
    )
    MOI.Test.test_linear_transform(bridged_mock, MOI.Test.Config())

    vis = MOI.get(bridged_mock, MOI.ListOfVariableIndices())
    @test vis == MOI.VariableIndex.([-1, -2])

    MOI.set(
        bridged_mock,
        MOI.ConstraintName(),
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }(),
        ),
        ["c1"],
    )
    MOI.set(
        bridged_mock,
        MOI.ConstraintName(),
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.LessThan{Float64},
            }(),
        ),
        ["c2"],
    )

    var_names = ["v1pos", "v2pos", "v1neg", "v2neg"]
    MOI.set(
        mock,
        MOI.VariableName(),
        MOI.get(mock, MOI.ListOfVariableIndices()),
        var_names,
    )
    MOI.set(
        bridged_mock,
        MOI.ConstraintName(),
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.VectorOfVariables,
                MOI.Nonnegatives,
            }(),
        ),
        ["nonneg"],
    )
    s = """
    variables: v1pos, v2pos, v1neg, v2neg
    nonneg: [v1pos, v2pos, v1neg, v2neg] in VecMathOptInterface.Nonnegatives(4)
    c1: v1pos + -1.0v1neg + v2pos + -1.0v2neg >= 1.0
    c2: v1pos + -1.0v1neg + v2pos + -1.0v2neg <= 2.0
    minobjective: v1pos + -1.0v1neg + v2pos + -1.0v2neg
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        mock,
        model,
        var_names,
        ["nonneg", "c1", "c2"],
    )
    var_names = ["v1", "v2"]
    MOI.set(bridged_mock, MOI.VariableName(), vis, var_names)
    s = """
    variables: v1, v2
    c1: v1 + v2 >= 1.0
    c2: v1 + v2 <= 2.0
    minobjective: v1 + v2
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        var_names,
        ["c1", "c2"],
    )
    _test_delete_bridged_variable(
        bridged_mock,
        vis[1],
        MOI.Reals,
        2,
        (
            (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
            (MOI.VectorOfVariables, MOI.Nonpositives, 0),
        ),
        used_bridges = 0,
        used_constraints = 0,
    )
    _test_delete_bridged_variable(
        bridged_mock,
        vis[2],
        MOI.Reals,
        1,
        (
            (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
            (MOI.VectorOfVariables, MOI.Nonpositives, 0),
        ),
    )
    return
end

end  # module

TestVariableFree.runtests()


text:
class Node {
  int elem;
  Node next;

  void setElem(int c) { elem = c; }

  void setNext(Node n) { next = n; }

  int getElem() { return elem; }

  Node getNext() { return next; }
}

class Stack {
  Node head;

  void push(int c) {
    Node newHead = new Node;
    newHead.setElem(c);
    newHead.setNext(head);
    head = newHead;
  }

  boolean isEmpty() {
    return head==(Node)null;
  }

  int top() {
    return head.getElem();
  }

  void pop() {
    head = head.getNext();
  }
}

int main() {
   Stack s = new Stack;
   int i= 0;
   while (i<10) {
     s.push(i);
     i++;
   }

   while (!s.isEmpty()) {
     printInt(s.top());
     s.pop();
   }
   return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
module App

using Printf
include("canvas.jl")
include("checkerboard.jl")
include("timetext.jl")
include("circles.jl")

function init_timer(canvas::Gtk.GtkCanvas, scene::Scene, interval)
    update_timer = Timer((timer) -> update(canvas, scene), 0.0, interval=interval)
    return update_timer
end
secondsperframe = 0.03 # 30 fps
init_timer(canvas, scene) = init_timer(canvas, scene, secondsperframe)

function hold_thread(window)
    cond = Condition()
    signal_connect(window, :destroy) do widget
        notify(cond)
    end
    wait(cond)
    return nothing
end

#     signal_connect(win, "check-resize") do win
#         println("resized")
#     end

function start()
    dims = (2400,1100)

    canvas = new_canvas(dims...)
    win = new_app_window(canvas)
    show(canvas)

    board = Checkerboard(dims...)
    contents = [TimeText(),board]
    scene = Scene(contents)
    init_timer(canvas, scene)
    draw(canvas -> draw_scene(canvas, scene), canvas)

    !isinteractive() && hold_thread(win)

    return (canvas,scene)
end

end # module

App.start()


text:
"""
    PeriodicTorsion(; periodicities, phases, ks)

A periodic torsion angle between four atoms.
"""
struct PeriodicTorsion{T, E} <: SpecificInteraction
    i::Int
    j::Int
    k::Int
    l::Int
    periodicities::Vector{Int}
    phases::Vector{T}
    ks::Vector{E}
end

PeriodicTorsion(; i, j, k, l, periodicities, phases, ks) = PeriodicTorsion{eltype(phases), eltype(ks)}(
    i, j, k, l, periodicities, phases, ks)

@inline @inbounds function force(d::PeriodicTorsion,
                                  coords,
                                  s::Simulation)
    ab = vector(coords[d.i], coords[d.j], s.box_size)
    bc = vector(coords[d.j], coords[d.k], s.box_size)
    cd = vector(coords[d.k], coords[d.l], s.box_size)
    cross_ab_bc = ab × bc
    cross_bc_cd = bc × cd
    bc_norm = norm(bc)
    θ = atan(
        ustrip(dot(cross_ab_bc × cross_bc_cd, bc / bc_norm)),
        ustrip(dot(cross_ab_bc, cross_bc_cd)),
    )
    fs = sum(zip(d.periodicities, d.phases, d.ks)) do (periodicity, phase, k)
        dEdθ = -k * periodicity * sin((periodicity * θ) - phase)
        fi =  dEdθ * bc_norm * cross_ab_bc / dot(cross_ab_bc, cross_ab_bc)
        fl = -dEdθ * bc_norm * cross_bc_cd / dot(cross_bc_cd, cross_bc_cd)
        v = (dot(-ab, bc) / bc_norm^2) * fi - (dot(-cd, bc) / bc_norm^2) * fl
        fj =  v - fi
        fk = -v - fl
        return [fi, fj, fk, fl]
    end
    return [d.i, d.j, d.k, d.l], fs
end

@inline @inbounds function potential_energy(d::PeriodicTorsion,
                                            s::Simulation)
    ab = vector(s.coords[d.i], s.coords[d.j], s.box_size)
    bc = vector(s.coords[d.j], s.coords[d.k], s.box_size)
    cd = vector(s.coords[d.k], s.coords[d.l], s.box_size)
    cross_ab_bc = ab × bc
    cross_bc_cd = bc × cd
    θ = atan(
        ustrip(dot(cross_ab_bc × cross_bc_cd, normalize(bc))),
        ustrip(dot(cross_ab_bc, cross_bc_cd)),
    )
    E = sum(zip(d.periodicities, d.phases, d.ks)) do (periodicity, phase, k)
        k + k * cos((periodicity * θ) - phase)
    end
    return E
end


text:
module NLSolversTests

using Test
using Gridap

using ..NonLinearOperatorMocks

op = NonLinearOperatorMock()

nls = NLSolver(show_trace=false,method=:newton)

x = solve(nls,op)
@test x ≈ [1.0, 3.0]

cache = solve!(x,nls,op)

@test string(cache) == "NLSolversCache object"

solve!(x,nls,op,cache)

cache.result

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""PersistentVolumeClaimVolumeSource references the user's PVC in the same namespace. This volume finds the bound PV and mounts that volume for the pod. A PersistentVolumeClaimVolumeSource is, essentially, a wrapper around another type of volume that is owned by someone else (the system).

    IoK8sApiCoreV1PersistentVolumeClaimVolumeSource(;
        claimName=nothing,
        readOnly=nothing,
    )

    - claimName::String : ClaimName is the name of a PersistentVolumeClaim in the same namespace as the pod using this volume. More info: https://kubernetes.io/docs/concepts/storage/persistent-volumes#persistentvolumeclaims
    - readOnly::Bool : Will force the ReadOnly setting in VolumeMounts. Default false.
"""
mutable struct IoK8sApiCoreV1PersistentVolumeClaimVolumeSource <: SwaggerModel
    claimName::Any # spec type: Union{ Nothing, String } # spec name: claimName
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly

    function IoK8sApiCoreV1PersistentVolumeClaimVolumeSource(;claimName=nothing, readOnly=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimVolumeSource, Symbol("claimName"), claimName)
        setfield!(o, Symbol("claimName"), claimName)
        validate_property(IoK8sApiCoreV1PersistentVolumeClaimVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        o
    end
end # type IoK8sApiCoreV1PersistentVolumeClaimVolumeSource

const _property_map_IoK8sApiCoreV1PersistentVolumeClaimVolumeSource = Dict{Symbol,Symbol}(Symbol("claimName")=>Symbol("claimName"), Symbol("readOnly")=>Symbol("readOnly"))
const _property_types_IoK8sApiCoreV1PersistentVolumeClaimVolumeSource = Dict{Symbol,String}(Symbol("claimName")=>"String", Symbol("readOnly")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1PersistentVolumeClaimVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1PersistentVolumeClaimVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1PersistentVolumeClaimVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PersistentVolumeClaimVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PersistentVolumeClaimVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PersistentVolumeClaimVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1PersistentVolumeClaimVolumeSource)
    (getproperty(o, Symbol("claimName")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1PersistentVolumeClaimVolumeSource }, name::Symbol, val)
end


text:
# --------------------------------------------------------------------------------
# Literature is a separate type hierarchy that's unrelated to Asset.
# We can still use the LiquidityStyle trait here.
# --------------------------------------------------------------------------------

abstract type Literature end

struct Book <: Literature
    name
end

# assign trait
LiquidityStyle(::Type{Book}) = IsLiquid()

# sample pricing function
marketprice(b::Book) = 10.0


text:
abstract type AbstractDiscreteFunction end
const ADFun = AbstractDiscreteFunction
abstract type CompositeFunction <: ADFun end
abstract type BasicDiscreteFunction <: ADFun end

# INTERFACES

# iterator for domain
domain(m :: ADFun) = error("not implemented")
# iterator for codomain
codomain(m :: ADFun) = error("not implemented")
# basic constructor from domain and codomain
(adf :: ADFun)(dom :: FinSet, cod :: FinSet) = error("not implemented")

"""
return factors (f1, ..., fp) for composite function,
"""
factors(fp :: CompositeFunction) = error("not implemented")
proj(f :: CompositeFunction, i) = factors(f)[i]
Base.getindex(f :: CompositeFunction, i) = proj(f, i)
Base.setindex!(ft :: CompositeFunction, v :: NTuple, k...) = ft[tuple(k...)] = v

rng(m :: AbstractDiscreteFunction) = codomain(m)
# number of functions
size(comp :: CompositeFunction) = length(comp.factors)

include("dfuns/generic.jl")
include("dfuns/residue_fun.jl")
include("dfuns/fun_product.jl")
include("dfuns/fun_tupling.jl")

function Base.show(io :: IO, f :: AbstractDiscreteFunction)
    d = domain(f)
    r = codomain(f)
    print(io, "Fun f : ")
    print(io, d)
    print(io, " -> ")
    print(io, r)
end

function Base.summary(f :: AbstractDiscreteFunction)
    d = domain(f)
    r = codomain(f)
    println("Domain : ", d)
    println("Codomain : ", r)
    for x in tablegen(f)
        println(x)
    end
end

# returns generator of type x => f(x)
function tablegen(f :: AbstractDiscreteFunction)
    d = domain(f)
    res = (x => f(x) for x in d)
    return res
end

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
@generated function inv{T}(A::StaticMatrix{T})
    S = typeof((one(T)*zero(T) + zero(T))/one(T))
    newtype = similar_type(A, S)

    if size(A) == (1,1)
        return quote
            $(Expr(:meta, :inline))
            @inbounds return $newtype((inv(A[1]),))
        end
    elseif size(A) == (2,2)
        return quote
            $(Expr(:meta, :inline))
            d = A[1]*A[4] - A[3]*A[2]
            @inbounds return $newtype((A[4]/d, -A[2]/d, -A[3]/d, A[1]/d))
        end
    elseif size(A) == (3,3)
        return quote
            $(Expr(:meta, :inline))
            @inbounds x0 = SVector(A[1], A[2], A[3])
            @inbounds x1 = SVector(A[4], A[5], A[6])
            @inbounds x2 = SVector(A[7], A[8], A[9])

            y0 = cross(x1,x2)
            d  = vecdot(x0, y0)
            x0 = x0 / d
            y0 = y0 / d
            y1 = cross(x2,x0)
            y2 = cross(x0,x1)

            @inbounds return $newtype((y0[1], y1[1], y2[1], y0[2], y1[2], y2[2], y0[3], y1[3], y2[3]))
        end
    else
        return quote # Implementation from Base
            AA = convert(Array{$S}, A) # lufact() doesn't work with StaticArrays at the moment... and the branches below must be type-stable
            if istriu(A)
                Ai = inv(UpperTriangular(AA))
            elseif istril(A)
                Ai = inv(LowerTriangular(AA))
            else
                Ai = inv(lufact(AA))
            end
            return convert(typeof(parent(Ai)), Ai)
        end
    end
end


text:
using DataFrames

include("coloredtable.jl")

function manage_problemset()
    script("IPython.notebook.metadata.mode = \"create\"; Homework.refresh_messages()")
end

function new_question(metadata_json, answer)

    metadata = JSON.parse(metadata_json)

    meta = Signal(Any, Dict())
    question = metadata["question"]

    display(set_metadata(question, Dict("answer" => JSON.json(encode(metadata, answer)))))
    display(set_metadata(question,
      alert("info", string(" ",
        "Will commit last result as the answer to question ",
        question))))

    # return the answer itself for consistency
    answer
end

function save_problemset()
    script("Homework.create_problemset(); Homework.refresh_messages()") |> display
end

function clear_admin_metadata()
    global_config["mode"] = "answering"
    script("Homework.clear_admin_metadata(); Homework.refresh_messages()") |> display
    display(Html("You can try out the answer cells as if you were the student now. Next run Homework.clear_answers() to remove the answers."))
end

function clear_answers()
    script("Homework.clear_answers(); Homework.refresh_messages()") |> display
    display(Html("All the answers should be cleared now, delete this cell and the cell containing Homework.manage_problemset, save the notebook and it will be good to distribute for answering."))
end

function progress(all=(get(global_config, "mode",  "") == "create"))
    mode = all ? "report" : "myreport"
    res = get(string(strip(get(global_config, "host", "https://juliabox.com"), ['/']), "/jboxplugin/hw/");
                    query = Dict("mode" => mode,
                    "params" => JSON.json(Dict(
                        "course" => global_config["course"],
                        "problemset" => global_config["problemset"]))),
                headers = Dict("Cookie" => global_config["cookie"]))

    if statuscode(res) == 200
        result = Requests.json(res)
        if result["code"] != 0
            display(Html("
 Something went wrong while getting the report
"))
            dump(result)
        else
            return @manipulate for report=Dict("Score" => "score", "Incorrect attempts" => "attempts")
                make_score_dataframe(result["data"], report) |>
                    x -> report == "score" ? color_progress(x) : x
            end
        end
    else
        display(MIME"text/html"(), "
 There was an error contacting the homework server
")
    end
end


text:
export FixedHomotopy

"""
    FixedHomotopy(H, t) <: AbstractSystem

Fix a homotopy `H(x,t)` at `t`
"""
mutable struct FixedHomotopy{Hom<:AbstractHomotopy,T} <: AbstractSystem
    H::Hom
    t::T
end

struct FixedHomotopyCache{HC<:AbstractHomotopyCache} <: AbstractSystemCache
    cache::HC
end

cache(FH::FixedHomotopy, x) = FixedHomotopyCache(cache(FH.H, x, FH.t))

Base.size(F::FixedHomotopy) = size(F.H)

evaluate!(u, F::FixedHomotopy, x, c::FixedHomotopyCache) =
    evaluate!(u, F.H, x, F.t, c.cache)
evaluate(F::FixedHomotopy, x, c::FixedHomotopyCache) = evaluate(F.H, x, F.t, c.cache)

jacobian!(U, F::FixedHomotopy, x, c::FixedHomotopyCache) =
    jacobian!(U, F.H, x, F.t, c.cache)
jacobian(F::FixedHomotopy, x, c::FixedHomotopyCache) = jacobian(F.H, x, F.t, c.cache)

evaluate_and_jacobian!(u, U, F::FixedHomotopy, x, c::FixedHomotopyCache) =
    evaluate_and_jacobian!(u, U, F.H, x, F.t, c.cache)
evaluate_and_jacobian(F::FixedHomotopy, x, c::FixedHomotopyCache) =
    evaluate_and_jacobian(F.H, x, F.t, c.cache)


text:
using DABenchmarks

ph_scaling = load_results("ph_scaling_1.json", Val{:scaling}())
ph_scaling.results["2"] = load_results("ph_scaling_2.json", Val{:scaling}()).results["2"]
ph_scaling.results["4"] = load_results("ph_scaling_4.json", Val{:scaling}()).results["4"]
ph_scaling.results["8"] = load_results("ph_scaling_8.json", Val{:scaling}()).results["8"]
ph_scaling.results["16"] = load_results("ph_scaling_16.json", Val{:scaling}()).results["16"]
save_results!(ph_scaling, "ph_scaling.json")

max_stars_repo_name: martinbiel/mpc-reproducibility

max_issues_repo_name: martinbiel/mpc-reproducibility

max_forks_repo_name: martinbiel/mpc-reproducibility


text:
"""Function to calculate a simple stack averaging solution"""

"""Averages all igrams in a stack

saves a 2D array of velocities: mm per year
"""
function run_stackavg(
    unw_stack_file::String,
    input_dset::String,
    outfile::String,
    outdset::String,
    geolist::Array{Date,1},
    valid_igram_list::Array{Igram,1};
    stack_all = true,
    reference_station = nothing,
    ref_row = nothing,
    ref_col = nothing,
)

    # Figure out which of all the igrams we want to use
    # chosen_igrams = pick_igrams(geolist)
    if stack_all
        chosen_igrams = valid_igram_list
    else
        chosen_igrams = pick_igrams(geolist, valid_igram_list)
    end

    # Get the full list including invalid/ignored to figure out indices to pick
    # from the .h5 file (we dont know this just from valid_igram_list, some are ignored)
    full_igram_list = load_intlist_from_h5(unw_stack_file)

    picked_igram_indices = indices_from_full(chosen_igrams, full_igram_list)
    println(
        "Out of $(length(picked_igram_indices)) desired igrams, " *
        "$(sum(isnothing.(picked_igram_indices))) are missing",
    )
    filter!(.!isnothing, picked_igram_indices)

    # TODO: fix loading all into mem
    # Load only these
    # input_dset = STACK_DSET
    # input_dset = STACK_FLAT_DSET
    println("Reading $(length(picked_igram_indices)) igrams out of '$input_dset' dataset from file ")
    @time unw_stack = load_hdf5_stack(unw_stack_file, input_dset, picked_igram_indices)

    if !isnothing(reference_station)
        println("Using $reference_station as reference")
        demrsc = Sario.load_dem_from_h5(unw_stack_file)
        ref_row, ref_col = MapImages.station_rowcol(reference_station, demrsc)
    end
    if !(isnothing(ref_row) && isnothing(ref_col))
        println("Shifting input stack to $ref_row, $ref_col")
        # Using col, row since we
        @time unw_stack .= shift_stack(unw_stack, ref_row, ref_col)  # , window=window
    end

    # Now with proper igrams picked, just divide the total phase by time diff sum
    timediffs = temporal_baseline(chosen_igrams)
    phase_sum = sum(unw_stack, dims = 3)[:, :, 1]
    avg_velo = phase_sum ./ sum(timediffs)

    # Finally, save as a mm/year velocity
    out = Float32.(P2MM * avg_velo)

    println("Writing solution into $outfile : $outdset")
    h5open(outfile, "cw") do f
        f[outdset] = permutedims(out)
        # TODO: Do I care to add this for stack when it's all the same?
        # f[_count_dset(cur_outdset)] = countstack
    end

    return outfile, outdset
end

"""Given the set of acquisition dates, find which pairs should be used
as interferograms to average

"""
function stack_indices(geolist::Array{Date,1})
    num_geos = length(geolist)
    # We round up here so that if there is an odd number,
    # we just ignore the middle date
    # e.g. dates [1,2,3,4,5] would be [1,4], [2,5]
    idx_span = ceil(num_geos / 2)
    num_igrams = floor(num_geos / 2)
    return [Int.((i, i + idx_span)) for i = 1:num_igrams]
end

"""Remove the first instance of `value` from `array`
If it is not containing, will return array
"""
function removefirst!(array, value)
    idx = findfirst(isequal(value), array)
    return isnothing(idx) ? array : deleteat!(array, idx)
end

function remove_all_igrams!(array, geo_dates)
    for date in geo_dates
        for idx = 1:2
            deleteat!(array, findall(x -> x[idx] == date, array))
        end
    end
end

"""Get all igrams that start with `geo_date`"""
find_igrams(igram_list::Array{Igram,1}, geo_date::Date) =
    filter(ig -> ig[1] == geo_date, igram_list)

max_igram(igram_list, geo_date) = find_igrams(igram_list, geo_date)[end]

function pick_igrams(full_geolist::Array{Date,1}, full_igram_list::Array{Igram,1})
    geolist = copy(full_geolist)
    igram_list = copy(full_igram_list)

    total_timespan = (geolist[end] - geolist[1]).value
    # Find this value to make sure no igram goes over half
    half_span = total_timespan / 2

    chosen_igrams = Array{Igram,1}()

    while length(geolist) > 1
        cur_geo = popfirst!(geolist)
        longest_igram = max_igram(igram_list, cur_geo)
        push!(chosen_igrams, longest_igram)

        early, late = longest_igram
        remove_all_igrams!(igram_list, [early, late])
        removefirst!(geolist, late)  # Only late, since first was popped off
    end

    return chosen_igrams

end

""" Ideal way to pick stack igrams: evenly spaced"""
function pick_igrams_default(geolist::Array{Date,1})::Array{Igram,1}
    geo_pairs = stack_indices(geolist)
    return [(geolist[i1], geolist[i2]) for (i1, i2) in geo_pairs]
end

"""Find the indices of the stack-picked igrams within the full stack
of all compute igrams (that are stored in the unw_stack file)"""
indices_from_full(chosen_igrams, full_igram_list) = indexin(chosen_igrams, full_igram_list)


text:
module PolynomialFactors

## TODO
## * performance is still really poor for larger degrees.

using AbstractAlgebra
using Combinatorics
import Primes
import LinearAlgebra
import LinearAlgebra: norm, dot, I

include("utils.jl")
include("polyutils.jl")
include("factor_zp.jl")
include("lll.jl")
include("factor_zx.jl")

export poly_factor, factormod
#export factor, rational_roots, factormod

end # module


text:
const _tanh_iter = 100
const _tanh_btol = 1e-8
const _tanh_dmin = 1e-4

struct TanhSpacing{E} <: AbstractDist{E}
    ds0::Float64
    ds1::Float64
    dsmax::Float64
    rsmax::Float64
    np::Int
end

function TanhSpacing(ds0::Real, ds1::Real, np::Int=0; dsmax::Real=-1, rsmax::Real=-1.0)
    E = np < 2
    if E
        if dsmax <= 0 && rsmax < 1.01
            error("`dsmax`(> 0) or `rsmax`(>= 1.01) must be specified if `np` < 2")
        end
        if dsmax > 0 && (dsmax < ds1 || dsmax < ds0)
            error("`dsmax` must be greater than `ds0` and `ds1`")
        end
    else
        dsmax <= 0 || @warn "`dsmax` is ignored since `np` has been specified"
        rsmax < 1.01 || @warn "`rsmax` is ignored since `np` has been specified"
    end
    return TanhSpacing{E}(ds0, ds1, dsmax, rsmax, np)
end

isvalid(d::TanhSpacing{true}) = d.dsmax > max(d.ds0, d.ds1, 0)
isvalid(d::TanhSpacing{false}) = d.np > 1

function nodes(d::TanhSpacing{false}; len::L=1.0) where {L<:Real}
    return tanh_dist(d.np, d.ds0 / len, d.ds1 / len) * len
end
function nodes(d::TanhSpacing, np::Int; len::L=1.0) where {L<:Real}
    return tanh_dist(np, d.ds0 / len, d.ds1 / len) * len
end

ratio_max(ds) =
    maximum(2:length(ds)) do i
        r = ds[i] / ds[i - 1]
        r < 1 ? 1 / r : r
    end

function nodes(d::TanhSpacing{true}; len::L=1.0) where {L<:Real}
    np = max(round(Int, len / d.dsmax) + 1, 3)
    s = nodes(d, np; len=len)
    ds = s[2:end] - s[1:(end - 1)]
    if (d.dsmax > 0)
        dsm = maximum(ds)
        while dsm > d.dsmax
            np = max(round(Int, np * sqrt(dsm / d.dsmax)), np + 1)
            s = nodes(d, np; len=len)
            ds = s[2:end] - s[1:(end - 1)]
            dsm = maximum(ds)
        end
    end
    if d.rsmax > 1.01
        rsm = ratio_max(ds)
        if (rsm > d.rsmax)
            np = max(round(Int, np * (rsm / d.rsmax)^0.45), np + 1)
            s = nodes(d, np; len=len)
            ds = s[2:end] - s[1:(end - 1)]
            rsm = ratio_max(ds)
            while rsm > d.rsmax
                np = max(round(Int, np * (rsm / d.rsmax)^0.45), np + 1)
                s = nodes(d, np; len=len)
                ds = s[2:end] - s[1:(end - 1)]
                rsm1 = ratio_max(ds)
                rsm1 > rsm && error("cannot determin `np`, check inputs.")
                rsm = rsm1
            end
        end
    end
    return s
end

function tanh_dist(np::Int, s0::Real, s1::Real)
    @assert np >= 2
    if np == 2
        return Float64[0, 1]
    end
    fix0 = s0 > 0
    fix1 = s1 > 0
    if fix0 && fix1
        s = tanh_2(np, convert(Float64, s0), convert(Float64, s1))
    elseif fix0
        s = tanh_1(np, convert(Float64, s0))
    elseif fix1
        s = tanh_1(np, convert(Float64, s1))
        reverse!(s)
        @. s = 1 - s
    else
        s = collect(LinRange(0, 1, np))
    end
    return s
end

function tanh_1(np::Int, ds0::Float64)::Vector{Float64}
    np1 = np - 1
    ds = 1 / np1
    b = ds / ds0
    if b < 1
        Δ = tanh_Δ2(b)
        ds1 = ds / (sinh(Δ) / Δ)
        return tanh_2(np, ds0, ds1)
    else
        s = Vector{Float64}(undef, np)
        Δ = tanh_Δ(b)
        Δ_2 = Δ / 2
        tnh2 = tanh(Δ_2)
        s[1] = 0
        s[end] = 1
        for i in 2:np1
            s[i] = (1+tanh(Δ_2*((i - 1)*ds - 1))/tnh2)
        end
        return s
    end
end

function tanh_2(np::Int, ds0::Float64, ds1::Float64)::Vector{Float64}
    np1 = np - 1
    ds = 1 / np1
    a = sqrt(ds1 / ds0)
    b = ds / sqrt(ds0 * ds1)
    s = Vector{Float64}(undef, np)
    s[1] = 0
    s[end] = 1
    if b < 1
        ds00 = ds0
        b0 = ds / ds00
        if b0 < 1
            Δ0 = tanh_Δ2(b0)
            ds10 = ds * Δ0 / sinh(Δ0)
        else
            Δ0 = tanh_Δ(b0)
            Δ0_2 = Δ0 / 2
            ds10 = ds * Δ0_2 / tanh(Δ0_2)
        end
        a0 = sqrt(ds10 / ds00)
        b0 = ds / sqrt(ds00 * ds10)
        Δ0 = tanh_Δ(b0)
        Δ0_2 = Δ0 / 2
        tnh20 = tanh(Δ0_2)
        ds11 = ds1
        b1 = ds / ds11

        if b1 < 1.0
            Δ1 = tanh_Δ2(b1)
            ds01 = ds * Δ1 / (sinh(Δ1))
        else
            Δ1 = tanh_Δ(b1)
            Δ1_2 = Δ1 / 2
            ds01 = ds * Δ1_2 / tanh(Δ1_2)
        end
        a1 = sqrt(ds11 / ds01)
        b1 = ds / sqrt(ds01 * ds11)
        Δ1 = tanh_Δ(b1)
        Δ1_2 = Δ1 / 2
        tnh21 = tanh(Δ1_2)
        for i in 2:np1
            u0 = (1 + tanh(Δ0_2 * (2 * (i - 1) * ds - 1.0)) / tnh20) / 2
            u1 = (1 + tanh(Δ1_2 * (2 * (i - 1) * ds - 1.0)) / tnh21) / 2
            d0 = u0 / (a0 + (1 - a0) * u0)
            d1 = u1 / (a1 + (1 - a1) * u1)
            s[i] = ((np1 - (i - 1)) * d0 + (i - 1) * d1) * ds
        end
    else
        Δ = tanh_Δ(b)
        Δ_2 = Δ / 2
        tnh2 = tanh(Δ_2)
        for i in 2:np1
            u = 0.5 * (1 + tanh(Δ_2 * (2 * (i - 1) * ds - 1.0)) / tnh2)
            s[i] = (u / (a + (1 - a) * u))
        end
    end
    return s
end

"""
    Compute a value for `Δ` given `b` (b>=-1)
    b = sinh(Δ)/Δ
"""
function tanh_Δ(b::Float64)::Float64
    isapprox(b, 1; atol=_tanh_btol) && return _tanh_dmin
    @assert b > 1

    Δ = max(min(sqrt(6 * (b - 1)), 30), _tanh_dmin)

    b0 = sinh(Δ) / Δ
    resid = (b - b0) / b
    @debug "[tanh_Δ] initial guess: Δ=$Δ, b=$b0, resisual=$resid"
    for i in 1:_tanh_iter
        @assert Δ > 0
        db = (cosh(Δ) - b0) / Δ
        @assert db != 0
        Δ += (b - b0) / db
        b0 = sinh(Δ) / Δ
        resid = (b - b0) / b
        isapprox(resid, 0; atol=_tanh_btol) && break
        @debug "[tanh_Δ] iter $i: Δ=$Δ, b=$b0, resisual=$resid"
    end
    @assert isapprox(resid, 0; atol=_tanh_btol) "Newton iteration converged poorly."
    return Δ
end
"""
    Compute a value for `Δ` given `b` at the opposite end. (b<=1)
    b = tanh(Δ/2)/(Δ/2)
"""
function tanh_Δ2(b::Float64)::Float64
    isapprox(b, 1; atol=_tanh_btol) && return _tanh_dmin
    @assert b < 1
    hmin = _tanh_dmin / 2

    Δ_2 = max(sqrt(6 * (1 - b)), hmin)
    b0 = tanh(Δ_2) / Δ_2
    resid = (b - b0) / b
    @debug "[tanh_Δ2] initial guess: Δ/2=$Δ_2, b=$b0, resisual=$resid"
    for i in 1:_tanh_iter
        @assert Δ_2 > 0
        db = (cosh(Δ_2)^-2 - b0) / Δ_2
        @assert db != 0
        Δ_2 += (b - b0) / db
        b0 = tanh(Δ_2) / Δ_2
        resid = (b - b0) / b
        isapprox(resid, 0; atol=_tanh_btol) && break
        @debug "[tanh_Δ2] iter $i: Δ/2=$Δ_2, b=$b0, resisual=$resid"
    end
    @assert isapprox(resid, 0; atol=_tanh_btol) "Newton iteration converged poorly."
    return Δ_2 * 2
end


text:
using Tokenize

# –––––––––––––––
# Some file utils
# –––––––––––––––

function readdir′(dir)
  try
    readdir(dir)
  catch e
    String[]
  end
end

isdir′(f) = try isdir(f) catch e false end
isfile′(f) = try isfile(f) catch e false end

files(dir) =
  @as x dir readdir′(x) map!(f->joinpath(dir, f), x ,x) filter!(isfile′, x)

dirs(dir) =
  @as x dir readdir′(x) filter!(f->!startswith(f, "."), x) map!(f->joinpath(dir, f), x, x) filter!(isdir′, x)

jl_files(dir::AbstractString) = @>> dir files filter!(f->endswith(f, ".jl"))

function jl_files(set)
  files = Set{String}()
  for dir in set, file in jl_files(dir)
    push!(files, file)
  end
  return files
end

# Recursion + Mutable State = Job Security
"""
Takes a start directory and returns a set of nearby directories.
"""
function dirsnearby(dir; descend = 1, ascend = 5, set = Set{String}())
  push!(set, dir)
  if descend > 0
    for down in dirs(dir)
      if !(down in set)
        push!(set, down)
        descend > 1 && dirsnearby(down, descend = descend-1, ascend = 0, set = set)
      end
    end
  end
  ascend > 0 && dirname(dir) !== dir && dirsnearby(dirname(dir), descend = descend, ascend = ascend-1, set = set)

  return set
end

# ––––––––––––––
# The Good Stuff
# ––––––––––––––
const SCOPE_STARTERS = [Tokens.BEGIN,
                        Tokens.WHILE,
                        Tokens.IF,
                        Tokens.FOR,
                        Tokens.TRY,
                        Tokens.FUNCTION,
                        Tokens.MACRO,
                        Tokens.LET,
                        Tokens.TYPE,
                        Tokens.DO,
                        Tokens.QUOTE,
                        Tokens.STRUCT]

const MODULE_STARTERS = [Tokens.MODULE, Tokens.BAREMODULE]

"""
Takes Julia source code and a line number, gives back the string name
of the module at that line.
"""
function codemodule(code, line)
  stack = String[]
  # count all unterminated block openers, brackets, and parens
  openers = [[0,0,0]]
  n_openers = 0
  n_brackets = 0
  n_parens = 0
  # index of next modulename token
  next_modulename = -1

  ts = tokenize(code)

  last_token = nothing

  for (i, t) in enumerate(ts)
    Tokens.startpos(t)[1] > line && break

    # Ignore everything in brackets or parnetheses, because any scope started in
    # them also needs to be closed in them. That way, we don't need special
    # handling for comprehensions and `end`-indexing.
    if Tokens.kind(t) == Tokens.LSQUARE
      openers[length(stack)+1][2] += 1
    elseif openers[length(stack)+1][2] > 0
      if Tokens.kind(t) == Tokens.RSQUARE
        openers[length(stack)+1][2] -= 1
      end
    elseif Tokens.kind(t) == Tokens.LPAREN
      openers[length(stack)+1][3] += 1
    elseif openers[length(stack)+1][3] > 0
      if Tokens.kind(t) == Tokens.RPAREN
        openers[length(stack)+1][3] -= 1
      end
    elseif Tokens.exactkind(t) in MODULE_STARTERS && (last_token == nothing || Tokens.kind(last_token) == Tokens.WHITESPACE) # new module
      next_modulename = i + 2
    elseif i == next_modulename && Tokens.kind(t) == Tokens.IDENTIFIER && Tokens.kind(last_token) == Tokens.WHITESPACE
      push!(stack, Tokens.untokenize(t))
      push!(openers, [0,0,0])
    elseif Tokens.exactkind(t) in SCOPE_STARTERS  # new non-module scope
      openers[length(stack)+1][1] += 1
    elseif Tokens.exactkind(t) == Tokens.END  # scope ended
      if openers[length(stack)+1][1] == 0
        !isempty(stack) && pop!(stack)
        length(openers) > 1 && pop!(openers)
      else
        openers[length(stack)+1][1] -= 1
      end
    end
    last_token = t
  end
  return join(stack, ".")
end

codemodule(code, pos::Cursor) = codemodule(code, pos.line)

"""
Takes a given Julia source file and another (absolute) path, gives the
line on which the path is included in the file or 0.
"""
function includeline(file::AbstractString, included_file::AbstractString)
  # check for erroneous self includes, doesn't detect more complex cycles though
  file == included_file && return 0

  line = 1
  tokens = Tokenize.tokenize(read(file, String))

  t, state = iterate(tokens)
  while true
    if Tokens.kind(t) == Tokens.WHITESPACE
      line += count(x -> x == '\n', t.val)
    elseif Tokens.kind(t) == Tokens.IDENTIFIER && t.val == "include"
      t, state = iterate(tokens, state)
      if Tokens.kind(t) == Tokens.LPAREN
        t, state = iterate(tokens, state)
        if Tokens.kind(t) == Tokens.STRING
          if normpath(joinpath(dirname(file), chop(t.val, head=1, tail=1))) == included_file
            return line
          end
        end
      end
    elseif Tokens.kind(t) == Tokens.ENDMARKER
      break
    end
    t, state = iterate(tokens, state)
  end
  return 0
end

"""
Takes an absolute path to a file and returns the (file, line) where that
file is included or nothing.
"""
function find_include(path::AbstractString)
  for file in @> path dirname dirsnearby jl_files
    line = -1
    try
      line = includeline(file, path)
    catch err
      return nothing
    end
    line > 0 && (return file, line)
  end
end

"""
Takes an absolute path to a file and returns a string
representing the module it belongs to.
"""
function filemodule_(path::AbstractString)
  loc = find_include(path)
  if loc != nothing
    file, line = loc
    mod = codemodule(read(file, String), line)
    super = filemodule(file)
    if super != "" && mod != ""
      return "$super.$mod"
    else
      return super == "" ? mod : super
    end
  end
  return ""
end

const filemodule = memoize(filemodule_)

# Get all modules
function children(m::Module)
  return @>> [moduleusings(m); getmodule.(Ref(m), string.(_names(m, all=true, imported=true)))] filter(x->isa(x, Module) && x ≠ m) unique
end

function allchildren(m::Module, cs = Set{Module}())
  for c in children(m)
    c in cs || (push!(cs, c); allchildren(c, cs))
  end
  return cs
end


text:
############################ Accessing Solution Values #####################

observables(sr::SimResult) = observables(sr.sim)
observables(sim::Simulation) = collect(keys(sim.vals.u[1]))
observables(sim::MCResult) = observables(sim[1])
constants(sim::SimResult) = sim.scenario.prob.p.constants

@inline Base.getindex(sim::Simulation, I...) = sim.vals[I...]
@inline Base.getindex(sim::Simulation, i::Symbol,::Colon) = [sim.vals[j][i] for j in 1:length(sim.vals)]
@inline Base.getindex(sr::SimResult, I...) = sr.sim[I...]
@inline Base.getindex(mc::MCResult, I...) = mc.sim[I...]

solat(sr::SimResult, args...) = solat(sr.sim, args...)

solat(sim::Simulation, t, idx, scope) = solat(sim, t, scope)[idx]
solat(sim::Simulation, t, idx::Colon, scope) = solat(sim, t, scope)

function solat(sim::Simulation, t, scope)
  if scope == :ode_ || scope == :start_
      _id = findfirst(x->x==t, sim.vals.t) # change to searchsortedfirst ?
  else
      t_id = findall(x->x==t, sim.vals.t)
      _id = t_id[findfirst(x->x==scope, @view(sim.scope[t_id]))]
  end
  return sim[_id]
end

(s::Simulation)(t, idx=:, scope=:ode_) = solat(s, t, idx, scope)
(sr::SimResult)(t, idx=:, scope=:ode_) = solat(sr.sim, t, idx, scope)

############################ DataFrames ########################################

function DataFrame(s::Simulation; vars=observables(s))
  df = DataFrame(t=s.vals.t)

  [df[!, v] = s[v,:] for v in vars[in.(vars, Ref(observables(s)))]]
  !isnothing(s.scope) && (df[!,:scope]=s.scope)

  return df
end

DataFrame(sr::SimResult; kwargs...) = DataFrame(sr.sim; kwargs...)

function DataFrame(sr::Pair{Symbol,S}; kwargs...) where S<:SimResult
  df = DataFrame(last(sr); kwargs...)
  df[!, :scenario] .= first(sr) # add new column

  return df
end

function DataFrame(res::Vector{Pair{Symbol,S}}; kwargs...) where S<:SimResult
  df_vectors = DataFrame.(res; kwargs...)

  return vcat(df_vectors...; cols=:union)
end

function DataFrame(mcr::MCResult; kwargs...)
  # df performance
  df = DataFrame()

  for (i,s) in enumerate(mcr.sim)
      dfs = DataFrame(s; kwargs...)
      insertcols!(dfs, 1, :iter => fill(i, length(s)))
      df = vcat(df,dfs)
  end

  return df
end

function DataFrame(mcr::Pair{Symbol,S}; kwargs...) where S<:MCResult
  df = DataFrame(last(mcr); kwargs...)
  df[!, :scenario] .= first(mcr) # add new column

  return df
end

function DataFrame(res::Vector{Pair{Symbol,S}}; kwargs...) where S<:MCResult
  df_vectors = DataFrame.(res; kwargs...)

  return vcat(df_vectors...; cols=:union)
end

############################ Save Result ########################################
"""
    save_results(filepath::String, sim::AbstractResult)

Save results as csv file

Arguments:

- `filepath`: path and name of the file to write to
- `sim`: simulation results of `AbstractResult` type
"""
save_results(filepath::String, sim::AbstractResult) = save_results(filepath, DataFrame(sim))

save_results(filepath::String, df::DataFrame) = CSV.write(filepath, df, delim=";")

#=FIXME
function save_results(path::String, mcsim::MCResult; groupby::Symbol=:observables)
  if groupby == :simulations
    for i in 1:length(mcsim)
      save_results("$path/$i.csv", mcsim[i])
    end
  elseif groupby == :observables
    obs = observables(mcsim[1])
    lobs = length(obs)
    for ob in obs
      df = DataFrame(t = mcsim[1].t)
      [df[!,string(i)] = mcsim[i][ob,:] for i in 1:length(mcsim)]
      CSV.write("$path/$ob.csv", df, delim=";")
    end
  end
end
=#


text:
using Documenter, DemoCards, AstroImages

# Deps for examples
ENV["GKSwstype"] = "nul"
using Plots, Photometry, ImageTransformations, ImageFiltering, WCS, Reproject

setup = quote
    using AstroImages
    using Random
    Random.seed!(123456)

    AstroImages.set_clims!(Percent(99.5))
    AstroImages.set_cmap!(:magma)
    AstroImages.set_stretch!(identity)
end
DocMeta.setdocmeta!(Photometry, :DocTestSetup, setup; recursive = true)

# 1. generate demo files
demopage, postprocess_cb, demo_assets = makedemos("examples") # this is the relative path to docs/

# if there are generated css assets, pass it to Documenter.HTML
assets = []
isnothing(demo_assets) || (push!(assets, demo_assets))

# 2. normal Documenter usage
format = Documenter.HTML(assets = assets)
makedocs(format = format,
         pages = [
            "Home" => "index.md",
            demopage,
         ],
         sitename = "Awesome demos")

makedocs(
    sitename="AstroImages.jl",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Getting Started" => "manual/getting-started.md",
            "Loading & Saving Images" => "manual/loading-images.md",
            "Displaying Images" => "manual/displaying-images.md",
            "Headers" => "manual/headers.md",
            "Dimensions and World Coordinates" => "manual/dimensions-and-world-coordinates.md",
            "Polarization" => "manual/polarization.md",
            "Spectral Axes" => "manual/spec.md",
            "Preserving Wrapper" => "manual/preserving-wrapper.md",
            "Conventions" => "manual/conventions.md",
        ],
        "Guides" => [
            "Blurring & Filtering Images" => "guide/image-filtering.md",
            "Transforming Images" => "guide/image-transformations.md",
            "Reprojecting Images" => "guide/reproject.md",
            "Extracting Photometry" => "guide/photometry.md",
        ],
        demopage,
        "API" => "api.md",
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    workdir=".."
)

# 3. postprocess after makedocs
postprocess_cb()

deploydocs(
    repo = "github.com/sefffal/AstroImages.jl.git",
    devbranch = "master",
    push_preview = true
)


text:
using StructuresKit

#example
CrossSectionNodes=[2.5 -3.5;2 -4;0 -4;0 4;-2 4;-2.5 3.5]  #member cross-section nodes, in x-y pairs
L=100  #member length
z=0:0.5:L  #z coordinates along the length

#define member deformed shape
a1=L/100*0
a2=L/1000*0
a3=L/1000*10
u=a1*sin.(π*z/L)
v=a2*sin.(π*z/L)
ϕ=a3*sin.(π*z/L)

#mesh the member deformed shape, for undeformed set u=v=ϕ=0
coordinates, connectivity=Visualize.Mesh3D(CrossSectionNodes, z, u, v, ϕ)

#plot the mesh with Makie
using Makie
scene = poly(coordinates, connectivity, color=:lightgray, shading=:true, show_axis=:false, overdraw=:false, strokecolor = (:black, 0.6), strokewidth = 4)


text:
# Copyright (c) 2015, Intel Corporation
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
#     * Redistributions of source code must retain the above copyright notice,
#       this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above copyright
#       notice, this list of conditions and the following disclaimer in the
#       documentation and/or other materials provided with the distribution.
#     * Neither the name of Intel Corporation nor the names of its contributors
#       may be used to endorse or promote products derived from this software
#       without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.

include("neuron.jl")

# Common initialization functions
export gaussian
function gaussian(;mean=0.0, std=1.0)
    function gaussian_init(eltype::Type, dims...)
        rand(eltype, dims...) * std + mean
    end
end

export xavier, xavier_init
function xavier_init(eltype::Type, dims...)
    fan_in = prod(dims[1:end-1])
    scale = sqrt(3.0 / fan_in)
    rand(eltype, dims...) * 2scale - scale
end

xavier = xavier_init

include("layers/convolution.jl")
include("layers/pooling.jl")
include("layers/relu.jl")
include("layers/inner-product.jl")
include("layers/softmax-loss.jl")
include("layers/accuracy.jl")
include("layers/hdf5-data.jl")
include("layers/memory-data.jl")
include("layers/dropout.jl")
include("layers/math.jl")
include("layers/lstm.jl")
include("layers/embed_id.jl")
include("layers/GRU.jl")
include("layers/transform.jl")
include("layers/concat.jl")
include("layers/rnn.jl")
include("layers/tanh.jl")
include("layers/reshape.jl")
include("layers/interpolation.jl")


text:
using Miscellaneous
using Test

@testset "Miscellaneous.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: itsdfish/Miscellaneous

max_issues_repo_name: itsdfish/Miscellaneous

max_forks_repo_name: itsdfish/Miscellaneous


text:
using JSON
using NMAP: NMAP, Marshalling
using Test
using StructTypes
using XMLDict

filepath = joinpath(resources, "example.xml")

@testset "Unmarshall to sink" begin
    xml     = read(filepath)
    iscan   = NMAP.Scan(xml)
    scan    = Marshalling.unmarshall(NMAP.Scan, read(filepath, String))
    @test scan.xml == iscan.xml

    xmlstr = read(filepath, String)
    @test Marshalling.unmarshall(String, xmlstr) == xmlstr

    xmld    = Marshalling.unmarshall(Dict, xmlstr)
    xmljson = Marshalling.unmarshall(NMAP.Json, xmlstr)
    @test JSON.json(xmld) == xmljson

    xmldict = XMLDict.parse_xml(read(filepath, String))
    xmlstr  = String(XMLDict.dict_xml(xmldict))
    @test Marshalling.unmarshall(String, xmldict)   == xmlstr
    @test Marshalling.unmarshall(Dict,   xmldict)   == XMLDict.xml_dict(xmldict)
    @test Marshalling.unmarshall(NMAP.Json, xmldict)== JSON.json(XMLDict.xml_dict(xmldict))
end

@testset "Unmarshall to Leaf types" begin
    xmldict  = XMLDict.parse_xml(read(filepath, String))
    testscan = NMAP.Scan(xmldict)

    leafs = []
    for name in names(NMAP, all=true, imported=false)
        if  getfield(NMAP, name) isa DataType &&
            getfield(NMAP, name) <: NMAP.Leaf

            push!(leafs, getfield(NMAP, name))
        end
    end
    @testset "Unmarshall Scan Leaf types" begin
        scanleafs = filter(x->x in NMAP.Scan.types, leafs)
        typemap   = Dict([name=>type for (name, type) in zip(fieldnames(Scan), Scan.types) if type in scanleafs])
        revmap    = Dict([type=>name for (name, type) in typemap])
        testscan  = Scan(read(filepath))

        for (k,v) in xmldict
            if any(x->x in keys(typemap), [k, Symbol(k)])
                leaf   = get(typemap, Symbol(k), get(typemap, k, k))
                object = StructTypes.construct(leaf, v)

                field = Marshalling.getfield(Scan, k)
                comparator = Base.getproperty(testscan, field)

                @test !isdiff(object, comparator)
            end
        end
    end
end


text:
module PkgSRA

    ## Base
    using OrdinaryDiffEq, Turing, Random
    ## For derivatives
    using DSP, Statistics
    ## General
    using Plots, StatsPlots, BSON
    using DataStructures
    ## Examples
    using Interpolations, StatsBase, Statistics
    ## Sindy-specific
    using Lasso

    # Core objects and types
    include("../utils/sparse_regression_functions.jl")
    include("sindy.jl")
    # Core numerical algorithms
    include("numerical_derivatives.jl")
    include("dmdc.jl")
    include("initialize_control_signal.jl")
    include("forcing_functions.jl")

    # Utilities
    include("../utils/array_utils.jl")
    include("../utils/regression_utils.jl")
    include("../utils/control_utils.jl")
    include("../utils/posterior_sampling_utils.jl")
    include("../utils/combinatorics_utils.jl");
    include("../utils/sindy_utils.jl")
    include("../utils/sindy_plotting.jl")
    include("../utils/sindy_statistics_utils.jl")

    # State of the iterative algorithm, and plotting utilities
    include("../utils/main_algorithm_utils.jl")
    include("../src/sra_model_object.jl");
    include("../src/sra_model_functions.jl");
    include("../src/sra_model_plotting.jl");

end # module


text:
# Let defines new name space
let state = 0
    global _get_next_market_id() = (state += 1)
 end

mutable struct DoubleAuctionMarket <: AbstractMarket
    name::String

    traded_asset::AbstractAsset
    currency::AbstractAsset
    last_price::Float64
    ask_limit_orders::Array{AskLimitOrder, 1}
    bid_limit_orders::Array{BidLimitOrder, 1}
    continuous::Bool

    function DoubleAuctionMarket(; name::Union{String, Nothing}=nothing)
        name = isnothing(name) ? _get_next_market_id() : name
        new(name, generic_stock, generic_currency, NaN, Array{AskLimitOrder, 1}(), Array{BidLimitOrder, 1}(), true) # , Array{BidLimitOrder, 1}()
    end

    function DoubleAuctionMarket(ccy::AbstractAsset, asset::AbstractAsset; name::Union{String, Nothing}=nothing)
        name = isnothing(name) ? _get_next_market_id() : name
        new(name, ccy, asset, NaN, Array{AskLimitOrder, 1}(), Array{BidLimitOrder, 1}(), true) # , Array{BidLimitOrder, 1}()
    end

    function DoubleAuctionMarket(asset::AbstractAsset)
        new(asset, generic_currency, NaN, Array{AskLimitOrder, 1}(), Array{BidLimitOrder, 1}(), true) # , Array{BidLimitOrder, 1}()
    end

    function DoubleAuctionMarket(asset::AbstractAsset, currency::AbstractAsset)
        new(asset, currency, NaN, Array{AskLimitOrder, 1}(), Array{BidLimitOrder, 1}(), true) # , Array{BidLimitOrder, 1}()
    end
end

function place!(market::DoubleAuctionMarket, order::LimitOrder)
    push!(market, order)
    #clear!(market)
end

function get_trade_price(buy::BidLimitOrder, sell::AskLimitOrder; market::DoubleAuctionMarket)
    # This is used in clearing two orders
    if buy.timestamp < sell.timestamp
        return Int64(buy.price)
    elseif sell.timestamp < buy.timestamp
        return Int64(sell.price)
    else
        throw
    end
end

# Market matching
"Cancel all orders"
function cancel_all!(market::DoubleAuctionMarket)
    # TODO: release assets from traders
    market.ask_limit_orders = typeof(market.ask_limit_orders)()
    market.bid_limit_orders = typeof(market.bid_limit_orders)()

end

function get_trade_price(market::DoubleAuctionMarket)
    if market.continuous
        # Continuous is matched per order basis (not per book)
        return nothing
    else
        # Call is matched per book basis
        (q, p) = get_equlibrium(market.bid_limit_orders, market.ask_limit_orders)
        return p
    end
end