text:
# Some methods for validated integration of ODEs

"""
    remainder_taylorstep!(f!, t, x, dx, xI, dxI, δI, δt, params)

Returns a remainder for the integration step for the dependent variables (`x`)
checking that the solution satisfies the criteria for existence and uniqueness.
"""
function remainder_taylorstep!(f!::Function, t::Taylor1{T},
        x::Vector{Taylor1{TaylorN{T}}}, dx::Vector{Taylor1{TaylorN{T}}},
        xI::Vector{Taylor1{Interval{T}}}, dxI::Vector{Taylor1{Interval{T}}},
        δI::IntervalBox{N,T}, δt::Interval{T}, params) where {N,T}

    orderT = get_order(dx[1])
    aux = δt^(orderT+1)
    Δx  = IntervalBox( [  xI[i][orderT+1] * aux for i in eachindex(xI)] )
    Δ0  = IntervalBox( [  dx[i][orderT](δI) * aux / (orderT+1) for i in eachindex(x)] )
    Δdx = IntervalBox( [ dxI[i][orderT+1] * aux for i in eachindex(xI)] )
    Δ = Δ0 + Δdx * δt
    Δxold = Δx

    # Checking existence and uniqueness
    iscontractive(Δ, Δx) && return Δx

    # If the check didn't work, compute new remainders. A new Δx is proposed,
    # and the corresponding Δdx is computed
    xxI  = Array{Taylor1{TaylorN{Interval{T}}}}(undef, N)
    dxxI = Array{Taylor1{TaylorN{Interval{T}}}}(undef, N)
    vv = Array{Interval{T}}(undef, N)
    for its = 1:10
        # Remainder of Picard iteration
        Δ = picard_remainder!(f!, t, x, dx, xxI, dxxI, δI, δt, Δx, Δ0, params)

        # Checking existence and uniqueness
        iscontractive(Δ, Δx) && return Δx
        # iscontractive(Δ, Δx) && return _contract_iteration!(f!, t, x, dx, xxI, dxxI, δI, δt, Δx, Δdx, Δ0, params)

        # Expand Δx in the directions needed
        Δxold = Δx
        if Δ ⊆ Δx
            @inbounds for ind in 1:N
                # Widen the directions where ⊂ does not hold
                vv[ind] = Δx[ind]
                if Δ[ind] == Δx[ind]
                    vv[ind] = widen.(Δ[ind])
                end
            end
            Δx = IntervalBox(vv)
            continue
        end
        Δx = Δ
    end

    # If it didn't work, throw an error
    @format full
    error("""
    Error: it cannot prove existence and unicity of the solution:
        t0 = $(t[0])
        δt = $(δt)
        Δ  = $(Δ)
        Δxo = $(Δxold)
        Δx = $(Δx)
        $(Δ .⊆ Δxold)
    """)
end

"""
    iscontractive(Δ, Δx)

Checks if `Δ .⊂ Δx` is satisfied. If ``Δ ⊆ Δx` is satisfied, it returns
`true` if all cases where `==` holds corresponds to the zero `Interval`.
"""
function iscontractive(Δ::IntervalBox{N,T}, Δx::IntervalBox{N,T}) where{N,T}
    zi = Interval{T}(0, 0)
    @inbounds for ind in 1:N
        Δ[ind] ⊂ Δx[ind] && continue
        Δ[ind] == Δx[ind] == zi && continue
        return false
    end
    return true
end

"""
    picard_remainder!(f!, t, x, dx, xxI, dxxI, δI, δt, Δx, Δ0, params)

Return the remainder of Picard operator
"""
function picard_remainder!(f!::Function, t::Taylor1{T},
    x::Vector{Taylor1{TaylorN{T}}}, dx::Vector{Taylor1{TaylorN{T}}},
    xxI::Vector{Taylor1{TaylorN{Interval{T}}}},
    dxxI::Vector{Taylor1{TaylorN{Interval{T}}}},
    δI::IntervalBox{N,T}, δt::Interval{T},
    Δx::IntervalBox{N,T}, Δ0::IntervalBox{N,T}, params) where {N,T}

    # Extend `x` and `dx` to have interval coefficients
    zI = zero(Interval{T})
    @inbounds for ind in eachindex(x)
        xxI[ind]  = x[ind] + Δx[ind]
        dxxI[ind] = dx[ind] + zI
    end

    # Compute `dxxI` from the equations of motion
    f!(dxxI, xxI, params, t)

    # Picard iteration, considering only the bound of `f` and the last coeff of f
    Δdx = IntervalBox( evaluate.( (dxxI - dx)(δt), δI... ) )
    Δ = Δ0 + Δdx * δt
    return Δ
end

# Picard iterations to contract further Δx, once Δ ⊂ Δx holds
# **Currently not used**
function _contract_iteration!(f!::Function, t::Taylor1{T},
        x::Vector{Taylor1{TaylorN{T}}}, dx::Vector{Taylor1{TaylorN{T}}},
        xxI::Vector{Taylor1{TaylorN{Interval{T}}}}, dxxI::Vector{Taylor1{TaylorN{Interval{T}}}},
        δI::IntervalBox{N,T}, δt::Interval{T},
        Δx::IntervalBox{N,T}, Δdx::IntervalBox{N,T}, Δ0::IntervalBox{N,T}, params) where {N,T}

    # Some abbreviations
    Δ = Δ0 + Δdx * δt
    Δxold = Δx

    # Picard contractions
    for its = 1:10
        # Remainder of Picard iteration
        Δ = picard_remainder!(f!, t, x, dx, xxI, dxxI, δI, δt, Δx, Δ0, params)

        # If contraction doesn't hold, return old bound
        iscontractive(Δ, Δx) || return Δxold

        # Contract estimate
        Δxold = Δx
        Δx = Δ
    end

    return Δxold
end

"""
    absorb_remainder(a::TaylorModelN)

Returns a TaylorModelN, equivalent to `a`, such that the remainder
is mostly absorbed in the constant and linear coefficients. The linear shift assumes
that `a` is normalized to the `IntervalBox(-1..1, Val(N))`.

Ref: Xin Chen, Erika Abraham, and Sriram Sankaranarayanan,
"Taylor Model Flowpipe Construction for Non-linear Hybrid
Systems", in Real Time Systems Symposium (RTSS), pp. 183-192 (2012),
IEEE Press.
"""
function absorb_remainder(a::TaylorModelN{N,T,T}) where {N,T}
    Δ = remainder(a)
    orderQ = get_order(a)
    δ = IntervalBox(Interval{T}(-1,1), Val(N))
    aux = diam(Δ)/(2N)
    rem = Interval{T}(0, 0)

    # Linear shift
    lin_shift = mid(Δ) + sum((aux*TaylorN(i, order=orderQ) for i in 1:N))
    bpol = a.pol + lin_shift

    # Compute the new remainder
    aI = a(δ)
    bI = bpol(δ)

    if bI ⊆ aI
        rem = Interval(aI.lo-bI.lo, aI.hi-bI.hi)
    elseif aI ⊆ bI
        rem = Interval(bI.lo-aI.lo, bI.hi-aI.hi)
    else
        r_lo = aI.lo-bI.lo
        r_hi = aI.hi-bI.hi
        if r_lo > 0
            rem = Interval(-r_lo, r_hi)
        else
            rem = Interval( r_lo, -r_hi)
        end
    end

    return TaylorModelN(bpol, rem, a.x0, a.dom)
end

"""
    shrink_wrapping!(xTMN::TaylorModelN)

Returns a modified inplace `xTMN`, which has absorbed the remainder
by the modified shrink-wrapping method of Florian Bünger.

Ref: Florian B\"unger, Shrink wrapping for Taylor models revisited,
Numer Algor 78:1001–1017 (2018), https://doi.org/10.1007/s11075-017-0410-1
"""
function shrink_wrapping!(xTMN::Vector{TaylorModelN{N,T,T}}) where {N,T}
    # Original domain of TaylorModelN
    B = IntervalBox(Interval{T}(-1,1), Val(N))
    zI = zero(Interval{T})

    # Vector of independent TaylorN variables
    order = get_order(xTMN[1])
    X = [TaylorN(T, i, order=order) for i in 1:N]

    # Remainder of original TaylorModelN and componentwise mag
    rem = remainder.(xTMN)
    r = mag.(rem)
    qB = r .* B

    # Shift to remove constant term
    xTN0 = constant_term.(xTMN)
    xTNcent = polynomial.(xTMN) .- xTN0

    # Jacobian (at zero) and its inverse
    jac = TaylorSeries.jacobian(xTNcent)
    ## If conditional number is too large, use absorb_remainder
    if cond(jac) > 1.0e4
        for i in eachindex(xTMN)
            # If remainder is still too big, do it again
            j = 0
            while (j < 10) && (mag(rem[i]) > 1.0e-12)
                j += 1
                xTMN[i] = absorb_remainder(xTMN[i])
                rem[i] = remainder(xTMN[i])
            end
        end
        return xTMN
    end
    invjac = inv(jac)

    # Componentwise bound
    r̃ = mag.(invjac * qB) # qB <-- r .* B
    qB´ = r̃ .* B
    @assert invjac * qB ⊆ qB´

    # Nonlinear part (linear part is close to identity)
    g = invjac*xTNcent .- X
    # g = invjac*(xTNcent .- linear_polynomial.(xTNcent))
    # ... and its jacobian matrix (full dependence!)
    jacmatrix_g = TaylorSeries.jacobian(g, X)

    # Step 7 of Bünger algorithm: obtain an estimate of `q`
    # Some constants/parameters
    q_tol = 1.0e-12
    q = 1.0 .+ r̃
    ff = 65/64
    q_max = ff .* q
    zs = zero(q)
    s = similar(q)
    q_old = similar(q)
    q_1 = similar(q)
    jaq_q1 = jacmatrix_g * r̃

    iter_max = 100
    improve = true
    iter = 0
    while improve && iter < iter_max
        qB .= q .* B
        s .= zs
        q_old .= q
        q_1 .= q_old .- 1.0
        mul!(jaq_q1, jacmatrix_g, q_1)
        @inbounds for i in eachindex(xTMN)
            s[i] = mag( jaq_q1[i](qB) )
            q[i] = 1.0 + r̃[i] + s[i]
            # If q is too large, return xTMN unchanged
            q[i] > q_max[i] && return xTMN
        end
        improve = any( ((q .- q_old)./q) .> q_tol )
        iter += 1
    end
    # (improve || q == one.(q)) && return xTMN

    # Compute final q and rescale X
    @. q = 1.0 + r̃ + ff * s
    @. X = q * X

    # Postverify and define Taylor models to be returned
    @inbounds for i in eachindex(xTMN)
        pol = polynomial(xTMN[i])
        ppol = fp_rpa(TaylorModelN(pol(X), zI, xTMN[i].x0, xTMN[i].dom ))
        if xTMN[i](B) ⊆ polynomial(ppol)(B) # assumes zero remainder
            xTMN[i] = TaylorModelN(pol, zI, xTMN[i].x0, xTMN[i].dom )
        else
            xTMN[i] = fp_rpa(TaylorModelN(xTMN[i](X), zI, xTMN[i].x0, xTMN[i].dom ))
        end
    end

    return xTMN
end

"""
    validated-step!
"""
function validated_step!(f!, t::Taylor1{T},
        x::Vector{Taylor1{TaylorN{T}}}, dx::Vector{Taylor1{TaylorN{T}}}, xaux::Vector{Taylor1{TaylorN{T}}},
        tI::Taylor1{T},
        xI::Vector{Taylor1{Interval{T}}}, dxI::Vector{Taylor1{Interval{T}}}, xauxI::Vector{Taylor1{Interval{T}}},
        t0::T, tmax::T, xTMN::Vector{TaylorModelN{N,T,T}},
        xv::Vector{IntervalBox{N,T}},
        rem::Vector{Interval{T}}, δq_norm::IntervalBox{N,T},
        q0::IntervalBox{N,T}, q0box::IntervalBox{N,T}, nsteps::Int,
        orderT::Int, abstol::T, params, parse_eqs::Bool,
        check_property::Function=(t, x)->true) where {N,T}

    # One step integration (non-validated)
    # TaylorIntegration.__jetcoeffs!(Val(parse_eqs), f!, t, x, dx, xaux, params)
    # δt = TaylorIntegration.stepsize(x, abstol)
    δt = TaylorIntegration.taylorstep!(f!, t, x, dx, xaux,
        t0, tmax, orderT, abstol, params, parse_eqs)
    # One step integration for the initial box
    # TaylorIntegration.__jetcoeffs!(Val(parse_eqs), f!, tI, xI, dxI, xauxI, params)
    # δtI = TaylorIntegration.stepsize(xI, abstol)
    δtI = TaylorIntegration.taylorstep!(f!, tI, xI, dxI, xauxI,
        t0, tmax, orderT+1, abstol, params, parse_eqs)

    # # Step size
    # δt = min(δt, tmax-t0, inf(δtI))
    δtI = Interval{T}(0.0, δt)

    # This updates the `dx[:][orderT]` and `dxI[:][orderT+1]`, which are currently zero
    f!(dx, x, params, t)
    f!(dxI, xI, params, tI)

    # Test if `check_property` is satisfied; if not, half the integration time.
    # If after 25 checks `check_property` is not satisfied, throw an error.
    nsteps += 1
    issatisfied = false
    rem_old = copy(rem)
    for nchecks = 1:25
        # Validate the solution: remainder consistent with Schauder thm
        Δ = remainder_taylorstep!(f!, t, x, dx, xI, dxI, δq_norm, δtI, params)
        rem .= rem_old .+ Δ

        # Create TaylorModelN to store remainders and evaluation
        @inbounds begin
            for i in eachindex(x)
                xTMN[i] = fp_rpa( TaylorModelN(x[i](0..δt), rem[i], q0, q0box) )

                # If remainder is still too big, do it again
                j = 0
                while (j < 10) && (mag(rem[i]) > 1.0e-10)
                    j += 1
                    xTMN[i] = absorb_remainder(xTMN[i])
                    rem[i] = remainder(xTMN[i])
                end
            end
            xv[nsteps] = evaluate(xTMN, δq_norm) # IntervalBox

            if !check_property(t0+δt, xv[nsteps])
                δt = δt/2
                continue
            end
        end # @inbounds

        issatisfied = true
        break
    end

    if !issatisfied
        error("""
            `check_property` is not satisfied:
            $t0 $nsteps $δt
            $(xv[nsteps])
            $(check_property(t0+δt, xv[nsteps]))""")
    end

    return δt
end

function validated_integ(f!, qq0::AbstractArray{T,1}, δq0::IntervalBox{N,T},
        t0::T, tmax::T, orderQ::Int, orderT::Int, abstol::T, params=nothing;
        maxsteps::Int=500, parse_eqs::Bool=true,
        check_property::Function=(t, x)->true) where {N, T<:Real}

    # Set proper parameters for jet transport
    @assert N == get_numvars()
    dof = N

    # Some variables
    zI = zero(Interval{T})
    q0 = IntervalBox(qq0)
    t   = t0 + Taylor1(orderT)
    tI  = t0 + Taylor1(orderT+1)
    δq_norm = IntervalBox(Interval{T}(-1, 1), Val(N))
    q0box = q0 .+ δq_norm

    # Allocation of vectors
    # Output
    tv    = Array{T}(undef, maxsteps+1)
    xv    = Array{IntervalBox{N,T}}(undef, maxsteps+1)
    xTM1v = Array{TaylorModel1{TaylorN{T},T}}(undef, dof, maxsteps+1)
    # Internals: jet transport integration
    x     = Array{Taylor1{TaylorN{T}}}(undef, dof)
    dx    = Array{Taylor1{TaylorN{T}}}(undef, dof)
    xaux  = Array{Taylor1{TaylorN{T}}}(undef, dof)
    xTMN  = Array{TaylorModelN{N,T,T}}(undef, dof)
    # Internals: Taylor1{Interval{T}} integration
    xI    = Array{Taylor1{Interval{T}}}(undef, dof)
    dxI   = Array{Taylor1{Interval{T}}}(undef, dof)
    xauxI = Array{Taylor1{Interval{T}}}(undef, dof)

    # Set initial conditions
    rem = Array{Interval{T}}(undef, dof)
    @inbounds for i in eachindex(x)
        qaux = normalize_taylor(qq0[i] + TaylorN(i, order=orderQ), δq0, true)
        x[i] = Taylor1( qaux, orderT)
        dx[i] = x[i]
        xTMN[i] = TaylorModelN(qaux, zI, q0, q0box)
        #
        xI[i] = Taylor1( q0[i]+δq0[i], orderT+1 )
        dxI[i] = xI[i]
        rem[i] = zI
        #
        xTM1v[i, 1] = TaylorModel1(deepcopy(x[i]), zI, zI, zI)
    end

    # Output vectors
    @inbounds tv[1] = t0
    @inbounds xv[1] = evaluate(xTMN, δq_norm)

    # Determine if specialized jetcoeffs! method exists (built by @taylorize)
    parse_eqs = parse_eqs && (length(methods(TaylorIntegration.jetcoeffs!)) > 2)
    if parse_eqs
        try
            TaylorIntegration.jetcoeffs!(Val(f!), t, x, dx, params)
        catch
            parse_eqs = false
        end
    end

    # Integration
    nsteps = 1
    while t0 < tmax

        # Validated step of the integration
        δt = validated_step!(f!, t, x, dx, xaux, tI, xI, dxI, xauxI,
            t0, tmax, xTMN, xv, rem, δq_norm, q0, q0box,
            nsteps, orderT, abstol, params, parse_eqs, check_property)

        # New initial conditions and time
        nsteps += 1
        t0 += δt
        @inbounds t[0] = t0
        @inbounds tI[0] = t0
        @inbounds tv[nsteps] = t0
        @inbounds for i in eachindex(x)
            xTM1v[i, nsteps] = TaylorModel1(deepcopy(x[i]), rem[i], zI, Interval{T}(0, δt))
            aux = x[i](δt)
            x[i]  = Taylor1( aux, orderT )
            dx[i] = Taylor1( zero(aux), orderT )
            auxI = xTMN[i](δq_norm)
            xI[i] = Taylor1( auxI, orderT+1 )
            dxI[i] = xI[i]
        end

        if nsteps > maxsteps
            @warn("""
            Maximum number of integration steps reached; exiting.
            """)
            break
        end

    end

    return view(tv,1:nsteps), view(xv,1:nsteps), view(xTM1v, :, 1:nsteps)
end


text:
function ichampf(t, fft_u, fft_v, A1, A2, Ntaumm, matr, conjmatr, Ktaubis, llambda, sigma)

    u = ifft(exp(1im * t * A1) * fft_u) * matr
    v = ifft(exp(1im * t * A1) * fft_v) * matr

    z = (u + conj(v)) / 2

    fz1 = abs(z) .^ (2 * sigma) * z
    champu = -1im * llambda * A2 * exp(-1im * t * A1) * fft(conjmatr * fz1)
    champv = -1im * llambda * A2 * exp(-1im * t * A1) * fft(conjmatr * conj(fz1))

    champu_fft = fft(champu, dims=1)
    champv_fft = fft(champv, dims=1)

    champumoy = champu_fft[0, :] / Ntaumm
    champvmoy = champv_fft[0, :] / Ntaumm

    champu_fft[0, :] = 0 * champu_fft[0, :]
    champv_fft[0, :] = 0 * champv_fft[0, :]
    champu = ifft(champu_fft / (1im * Ktaubis), dims=1)
    champv = ifft(champv_fft / (1im * Ktaubis), dims=1)

    return champu, champv, champumoy, champvmoy

end

ext: jl


text:
scenario(nx = 60, ny = 60, nz = 60; case = :torus2) = begin
    F = Float32
    vol = zeros(F, nx, ny, nz)

    sx, sy, sz = size(vol) ./ F(16)
    tx = F(nx) / 2sx
    ty = F(ny) / 2sy + F(1.5)
    tz = F(nz) / 2sz

    r = F(1.85)
    R = F(4)

    cushin(x, y, z) =
        z^2 * x^2 - z^4 - 2z * x^2 + 2z^3 + x^2 - z^2 - (x^2 - z)^2 - y^4 - 2x^2 * y^2 -
        y^2 * z^2 +
        2y^2 * z +
        y^2

    torus2(x, y, z) =
        ((x^2 + y^2 + z^2 + R^2 - r^2)^2 - 4R^2 * (x^2 + y^2)) *
        ((x^2 + (y + R)^2 + z^2 + R^2 - r^2)^2 - 4R^2 * ((y + R)^2 + z^2))

    sphere(x, y, z) = (
        ((x - 2)^2 + (y - 2)^2 + (z - 2)^2 - 1) *
        ((x + 2)^2 + (y - 2)^2 + (z - 2)^2 - 1) *
        ((x - 2)^2 + (y + 2)^2 + (z - 2)^2 - 1)
    )

    plane(x, y, z) = x + y + z - 3

    cassini(x, y, z) =
        (x^2 + y^2 + z^2 + F(0.45)^2)^2 - 16 * F(0.45)^2 * (x^2 + z^2) - F(0.5)^2

    blooby(x, y, z) = x^4 - 5x^2 + y^4 - 5y^2 + z^4 - 5z^2 + F(11.8)

    hyperboloid(x, y, z) = x^2 + y^2 - z^2 - 1

    callback = (;
        cushin = cushin,
        torus2 = torus2,
        sphere = sphere,
        plane = plane,
        cassini = cassini,
        blooby = blooby,
        hyperboloid = hyperboloid,
    )[case]

    for k ∈ 1:nz, j ∈ 1:ny, i ∈ 1:nx
        vol[i, j, k] = callback((i - 1) / sx - tx, (j - 1) / sy - ty, (k - 1) / sz - tz)
    end

    MC(vol, Int32)
end

output(PlyIO::Module, m::MC, fn::AbstractString = "test.ply") = begin
    ply = PlyIO.Ply()
    push!(
        ply,
        PlyIO.PlyElement(
            "vertex",
            PlyIO.ArrayProperty("x", Float32[v[1] for v ∈ m.vertices]),
            PlyIO.ArrayProperty("y", Float32[v[2] for v ∈ m.vertices]),
            PlyIO.ArrayProperty("z", Float32[v[3] for v ∈ m.vertices]),
            PlyIO.ArrayProperty("nx", Float32[n[1] for n ∈ m.normals]),
            PlyIO.ArrayProperty("ny", Float32[n[2] for n ∈ m.normals]),
            PlyIO.ArrayProperty("nz", Float32[n[3] for n ∈ m.normals]),
        ),
    )
    vertex_indices = PlyIO.ListProperty("vertex_indices", UInt8, Int32)
    for i ∈ eachindex(m.triangles)
        push!(vertex_indices, m.triangles[i] .- 1)
    end
    push!(ply, PlyIO.PlyElement("face", vertex_indices))

    PlyIO.save_ply(ply, fn, ascii = true)
    return
end


text:
@testset "Test the pivoted QR" begin

    ψ = randn(100, 20)
    u = randn(100)
    m, n = size(ψ)

    ctrue = ψ\u

    # Force a greedy approach over the entire set of indices
    # With pivoted QR
    idxpivot, cpivot, ϵpivot = pivotedqr(ψ, u; invert = true, verbose = true, ϵrel = eps())

    @test norm(ctrue[idxpivot]-cpivot)<1e-14

    # Residual error must decrease
    @test all(ϵpivot[2:end]-ϵpivot[1:end-1] .< 0) == true

    for k=1:n
        cverif = view(ψ,:, idxpivot[1:k])\u
        @test abs(ϵpivot[k+1] - norm(u - view(ψ,:, idxpivot[1:k])*cverif))<1e-14
    end
end

max_forks_repo_name: mleprovost/QROMP


text:
"""A `Behavior` object contains the mutable information in a filter bank. The
values of these fields may be changed *after* the construction of the filter
bank, without having to recompute the underlying architecture. Fields:

* `is_ϕ_applied`: true if and only if the lowpass filter `ϕ` (also known as scaling function) in addition to the bandpass filters `ψ`'s.

* `j_range`: range of wavelet octaves that are actually used in the convolutions. Default is all of them.

* `log2_oversampling`: base-2 logarithm of the oversampling factor with respect to the critical sampling rate. Must be positive. Default is 0, i.e. no oversampling.

* `max_log2_stride`: base-2 logarithm of the maximum distance between neighboring coefficients after subsampling (also known as *hop size* or *stride*). Must be positive. Default is the number of octaves minus one, which imposes no oversampling per se.

* `pathkey`: key of the variable over which the filter bank is applied.

* `weighting`: scalar weighting of each wavelet frequency band. Default is `EqualWeighting()`, i.e. all weights are one. Can alternatively be set to `LoudnessWeighting(samplerate)` to model the relative loudness perceived by the human ear, as defined by the international standard 61672:2003.

* `weights`: vector of floating-point weights corresponding to the frequencies of `ψ`'s. The weight assigned to `ϕ` is the same as the weight of the `ψ` with the lowest center frequency."""
type Behavior{T<:Real}
    ϕ_log2_sampling::Int
    ψ_log2_samplings::Vector{Int}
    is_ϕ_applied::Bool
    j_range::UnitRange
    log2_oversampling::Int
    max_log2_stride::Int
    pathkey::PathKey
    weighting::AbstractWeighting
    weights::Array{T, 3}
end

"""Given a lowpass filter `ϕ`, an array of wavelets `ψs`, and the corresponding
filter bank specification `spec`, the `Behavior` outer constructor computes
the critical sampling rates of all `ψs` and `ϕ`. Then, with the inputs
`log2_oversampling` and `max_log2_stride`, these sampling rates are bounded:
* from below, by the inverse of the maximum stride.
* from above, by 1, to avoid unnecessary upsampling with respect to the original
  sample rate.
NB: if `log2_oversampling` and `max_log2_stride` are left as defaults, all
sample rates are critical."""
function Behavior{
    T<:Number,
    D<:AbstractDomain,
    G<:AbstractPointGroup,
    W<:RedundantWaveletClass}(
        ϕ::AbstractFilter{T,D},
        ψs::AbstractArray{AbstractFilter{T,D},3},
        spec::AbstractSpec{T,D,G,W},
        is_ϕ_applied::Bool,
        j_range::UnitRange{Int},
        log2_oversampling::Int,
        max_log2_stride::Int,
        pathkey::PathKey,
        weighting::AbstractWeighting)
    ϕ_critical_log2_sampling = critical_log2_sampling(ϕ, spec)
    ϕ_log2_sampling =
        clamp(ϕ_critical_log2_sampling + log2_oversampling, -max_log2_stride, 0)
    ψ_critical_log2_samplings =
        Int[ critical_log2_sampling(ψ, spec) for ψ in ψs[1, end, :] ]
    ψ_log2_samplings = clamp(ψ_critical_log2_samplings + log2_oversampling,
        -max_log2_stride, 0)
    max_log2_stride = - min(ϕ_log2_sampling, minimum(ψ_log2_samplings))
    ξs = map(get_centerfrequency, spec.ψmetas)
    weights = map(T, weight_frequencies(weighting, ξs))
    Behavior(ϕ_log2_sampling, ψ_log2_samplings, is_ϕ_applied, j_range,
        log2_oversampling, max_log2_stride, pathkey, weighting, weights)
end


text:
using Test
using ComparativeAutograder
using ComparativeAutograder: @capture_stdstreams, smallrepr, truncatestring, errorstring

@testset "Test @capture_stdstreams" begin
    out, err = @capture_stdstreams begin
        println(stdout, "Hello, stdout!")
        println(stderr, "Hello, stderr!")
    end
    @test out == "Hello, stdout!\n"
    @test err == "Hello, stderr!\n"

    out, err = @capture_stdstreams begin
        @debug "Hello, world."
    end
    @test occursin("Debug:", err)
    @test occursin("Hello, world.", err)
    @test out == ""
end

@testset "Test smallrepr" begin
    # The `repr` of this has length 4001000 and the smallrepr has length 1587.
    onesrepr = smallrepr(ones(1000, 1000))
    @test length(onesrepr) <= 2000
    @test startswith(onesrepr, "1000×1000 Array{Float64,2}:\n")
    @test smallrepr("foo") == "\"foo\""
end

@testset "Test truncatestring" begin
    @test truncatestring("foo") == "foo"
    @test truncatestring("foo", 1) == "f...\n[ TRUNCATED ]"
    foos = repeat("foo", 12345)
    @test length(truncatestring(foos)) < length(foos)
end

function get_exception(f::Function)
    try
        f()
    catch exc
        return exc
    end
    return nothing
end

@testset "Test errorstring" begin
    @test errorstring("foo") == "foo"

    @testset "Test errorstring for builtin exception types." begin
        boundserror = errorstring(get_exception() do
            myarray = [1, 2, 3]
            value = myarray[4]
        end)
        @test occursin("BoundsError", boundserror)
        @test occursin("at index", boundserror)

        errorexception = errorstring(get_exception() do
            error("foo")
        end)
        @test occursin("ErrorException", errorexception)
        @test occursin("foo", errorexception)
    end

    # In Julia, anything can be thrown (including structs, primitives, etc.)
    @testset "Test errorstring for nonstandard exception types." begin
        interror = errorstring(get_exception(() -> throw(123)))
        @test occursin("123", interror)
    end
end


text:
# This file was generated, do not modify it. # hide
cap_sum = DataFrame(cap_sum);
cap_sum_ctry = DataFrame(cap_sum_ctry);
cap_share = leftjoin(cap_sum, cap_sum_ctry, on = :country, makeunique = true)
cap_share.capacity_mw_share = cap_share.capacity_mw_sum ./ cap_share.capacity_mw_sum_1;


text:
# Script to run the MIO formulation from Bertsimas and Mišić for n in [200,500], K in [10, 100, 200, 500, 100].
# This script will only solve the MIO directly, and will skip the n = 500, K = 1000 case (this is in MIOexpdata1_BMonly_n500_K1000.jl).
# This script produces part of the results needed for Table 3 of the paper.
#
# The script takes command line arguments to determine which instances from MIOexpdata1 to do. To run the first 20 instances
# (as per the paper), type
#
# > julia MIOexpdata1_BMonly.jl 1 20
#
# into the command line.

using JuMP, Gurobi, Stats, Distributions

include("../optimalPLD_code/optimalPLD_createPath.jl")
include("../optimalPLD_code/optimalPLD_solveSAO_constraints.jl")

srand(1204)
machine = 2;

n_vec = [200, 500]
numPerm_vec = [100, 200, 500, 1000]

numReps = 100;

counter = 0;

LR = parse(ARGS[1]);
UR = parse(ARGS[2]);

expdirpath = optimalPLD_createPath(machine, "exp_dir", "MIOexpdata1_BMonly_reps$(LR)to$(UR)");
mkdir(expdirpath)

b_vec = [5; 10];

MIP_timeLimit = 0 #1*3600;
isRelaxation = false;

data_prefix = "MIOexpdata1"

outcsvfilepath = string(expdirpath, "out.csv");
outcsvhandle = open(outcsvfilepath, "w")
row_string = string("n", ",",
					"numPermutations", ",",
					"rep", ",",
					"b_feas", ",",
					"sumxi", ",",
					"BM_expRevenue", ",",
					"BM_upperBound", ",",
					"BM_solutiontime", ",",
					"BM_MIPgap")
row_string = string(row_string, "\n");

print(outcsvhandle, row_string);
flush(outcsvhandle);

rep_range = LR:UR

counter = 0;
for nind = 1:length(n_vec)
	n = n_vec[nind]
	for nPind = 1:length(numPerm_vec)
		numPermutations = numPerm_vec[nPind];

		println("**** Reading from data directory...")
		data_name = string(data_prefix,"_neq",n,"_Keq",numPermutations)
		data_dir = optimalPLD_createPath(machine, "data_dir", data_name)
		revenues_mat = readcsv(string(data_dir, "revenues_mat.csv"))
		lambda_mat = readcsv(string(data_dir, "lambda_mat.csv"))
		orderings_mat = readcsv(string(data_dir, "orderings_mat.csv"))

		orderings_mat = convert(Array{Int64}, orderings_mat);

		# Skip the n = 500, K = 1000 case
		if ( (n == 500) & (numPermutations == 1000) )
			break;
		end

		for rep in rep_range

			revenues = revenues_mat[rep,:][:]
			lambda = lambda_mat[rep,:] #[1:numPermutations]
			orderings = orderings_mat[numPermutations*(rep-1)+1:rep*numPermutations,:]

			b_vec_actual = [b_vec; n]

			for b_feas in b_vec_actual

				A_feas = ones(1,n);

				println("**** Solving n = $n, nP = $numPermutations, rep = $rep ****")

				counter = counter + 1;
				isRelaxation = false;
				x_initial = [];

				if (b_feas >= n)
					println("**** Solving Bertsimas and Misic MIO (unconstrained) ****")
					x_val, y_val, BM_expRevenue, BM_upperBound, BM_solutiontime, BM_MIPgap = optimalPLD_solveSAO_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, [], [] )
				else
					println("**** Solving Bertsimas and Misic MIO (constrained) ****")
					x_val, y_val, BM_expRevenue, BM_upperBound, BM_solutiontime, BM_MIPgap = optimalPLD_solveSAO_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas)
				end

				sumxi = sum(x_val);

				row_string = string(n, ",",
						numPermutations, ",",
						rep, ",",
						b_feas, ",",
						sumxi, ",",
						BM_expRevenue, ",",
						BM_upperBound, ",",
						BM_solutiontime, ",",
						BM_MIPgap) #, ",",

				row_string = string(row_string, "\n");

				print(outcsvhandle, row_string);
				flush(outcsvhandle);
			end
		end
	end
end

close(outcsvhandle)

max_stars_repo_name: vvmisic/EFCPLO

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: vvmisic/EFCPLO

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: vvmisic/EFCPLO

max_forks_repo_licenses: MIT Unlicense


text:
function gen_night_valid(source_square::UInt64)
    target_squares = zeros(UInt64, 0)
    for cj in zip(NIGHT_CLEAR_FILES, NIGHT_JUMPS)
        candidate_square = (source_square & cj[1]) >> cj[2]
        if candidate_square != EMPTY
            push!(target_squares, candidate_square)
        end
    end
    return target_squares
end

function gen_all_night_valid_moves()
    night_moves = Dict{UInt64, Array{UInt64,1}}()
    for i in 1:64
        night_moves[INT2UINT[i]] = gen_night_valid(INT2UINT[i])
    end
    return night_moves
end
const NIGHT_MOVES = gen_all_night_valid_moves()


text:
include("SentinelMissings.jl")
import .SentinelMissings
import YAXArrays.DAT: DATConfig
import YAXArrays.YAXTools: PickAxisArray

struct CubeIterator{R,ART,ARTBC,LAX,ILAX,S}
    dc::DATConfig
    r::R
    inars::ART
    inarsBC::ARTBC
    loopaxes::LAX
end
Base.IteratorSize(::Type{<:CubeIterator}) = Base.HasLength()
Base.IteratorEltype(::Type{<:CubeIterator}) = Base.HasEltype()
Base.eltype(i::Type{<:CubeIterator{A,B,C,D,E,F}}) where {A,B,C,D,E,F} = F

tuplelen(::Type{<:NTuple{N,<:Any}}) where {N} = N

lift64(::Type{Float32}) = Float64
lift64(::Type{Int32}) = Int64
lift64(T) = T

defaultval(t::Type{<:AbstractFloat}) = convert(t, NaN)
defaultval(t::Type{<:Signed}) = typemin(t) + 1
defaultval(t::Type{<:Unsigned}) = typemax(t) - 1

function CubeIterator(
    dc,
    r;
    varnames::Tuple = ntuple(i -> Symbol("x$i"), length(dc.incubes)),
    include_loopvars = (),
)
    loopaxes = ntuple(i -> dc.LoopAxes[i], length(dc.LoopAxes))
    inars, _ = getCubeCache(dc)
    length(varnames) == length(dc.incubes) ||
        error("Supplied $(length(varnames)) varnames and $(length(dc.incubes)) cubes.")
    inarsbc = map(dc.incubes, inars) do ic, cache
        allax = getindsall(geticolon(ic), ic.loopinds, i -> true)
        if ic.colonperm === nothing
            pa = PickAxisArray(cache, allax)
        else
            pa = PickAxisArray(cache, allax, ic.colonperm)
        end
    end
    et = map(inarsbc) do ibc
        ti = eltype(ibc)
        if ti <: AbstractArray
            ti
        else
            ti = Base.nonmissingtype(ti)
            SentinelMissings.SentinelMissing{ti,defaultval(ti)}
        end
    end
    if include_loopvars == true
        include_loopvars = map(axname, (loopaxes...,))
    end
    if isa(include_loopvars, Tuple) && !isempty(include_loopvars)
        ilax = map(i -> findAxis(i, collect(loopaxes)), include_loopvars)
        any(isequal(nothing), ilax) && error("Axis not found in cubes")
        et = (et..., map(i -> eltype(loopaxes[i]), ilax)...)
    else
        ilax = ()
    end

    elt = NamedTuple{
        (map(Symbol, varnames)..., map(Symbol, include_loopvars)...),
        Tuple{et...},
    }
    CubeIterator{typeof(r),typeof(inars),typeof(inarsbc),typeof(loopaxes),ilax,elt}(
        dc,
        r,
        inars,
        inarsbc,
        loopaxes,
    )
end
iternames(::CubeIterator{<:Any,<:Any,<:Any,<:Any,<:Any,E}) where {E} = E
tuplenames(t::Type{<:NamedTuple{N}}) where {N} = string.(N)
function Base.show(io::IO, ci::CubeIterator{<:Any,<:Any,<:Any,<:Any,<:Any,E}) where {E}
    print(
        io,
        "Datacube iterator with ",
        length(ci),
        " elements with fields: ",
        tuplenames(E),
    )
end
Base.length(ci::CubeIterator) = prod(length.(ci.loopaxes))
function Base.iterate(ci::CubeIterator)
    rnow, blockstate = iterate(ci.r)
    updatears(ci.dc.incubes, rnow, :read, ci.inars)
    innerinds = CartesianIndices(length.(rnow))
    indnow, innerstate = iterate(innerinds)
    offs = map(i -> first(i) - 1, rnow)
    getrow(ci, ci.inarsBC, indnow, offs),
    (rnow = rnow, blockstate = blockstate, innerinds = innerinds, innerstate = innerstate)
end
function Base.iterate(ci::CubeIterator, s)
    t1 = iterate(s.innerinds, s.innerstate)
    N = tuplelen(eltype(ci.r))
    if t1 === nothing
        t2 = iterate(ci.r, s.blockstate)
        if t2 === nothing
            return nothing
        else
            rnow = t2[1]
            blockstate = t2[2]
            updatears(ci.dc.incubes, rnow, :read, ci.inars)
            innerinds = CartesianIndices(length.(rnow))
            indnow, innerstate = iterate(innerinds)

        end
    else
        rnow, blockstate = s.rnow, s.blockstate
        innerinds = s.innerinds
        indnow, innerstate = iterate(innerinds, s.innerstate)
    end
    offs = map(i -> first(i) - 1, rnow)
    getrow(ci, ci.inarsBC, indnow, offs),
    (
        rnow = rnow::NTuple{N,UnitRange{Int64}},
        blockstate = blockstate::Int64,
        innerinds = innerinds::CartesianIndices{N,NTuple{N,Base.OneTo{Int64}}},
        innerstate = innerstate::CartesianIndex{N},
    )
end
abstract type CubeRow end
abstract type CubeRowAx <: CubeRow end

function getrow(
    ci::CubeIterator{<:Any,<:Any,<:Any,<:Any,ILAX,S},
    inarsBC,
    indnow,
    offs,
) where {ILAX,S}
    axvalsall = map((ax, i, o) -> ax.values[i+o], ci.loopaxes, indnow.I, offs)
    axvals = map(i -> axvalsall[i], ILAX)
    cvals = map(i -> i[indnow], inarsBC)
    S((cvals..., axvals...))
end
function getrow(
    ci::CubeIterator{<:Any,<:Any,<:Any,<:Any,(),S},
    inarsBC,
    indnow,
    offs,
) where {S}
    cvals = map(i -> i[indnow], inarsBC)
    S(cvals)
end

function Base.show(io::IO, s::CubeRow)
    print(io, "Cube Row: ")
    for n in propertynames(s)
        print(io, string(n), "=", getproperty(s, n), " ")
    end
end
function Base.show(io::IO, s::CubeRowAx)
    print(io, "Cube Row: ")
    for n in propertynames(s)
        print(io, string(n), "=", getproperty(s, n), " ")
    end
end
function Base.show(io::IO, s::Type{<:CubeRow})
    foreach(fieldnames(s)) do fn
        print(io, fn, "::", fieldtype(s, fn), ", ")
    end
end
function Base.iterate(s::CubeRow, state = 1)
    allnames = propertynames(s)
    if state <= length(allnames)
        (getproperty(s, allnames[state]), state + 1)
    else
        nothing
    end
end

export CubeTable
"""
    CubeTable()

Function to turn a DataCube object into an iterable table. Takes a list of as arguments,
specified as a `name=cube` expression. For example
`CubeTable(data=cube1,country=cube2)` would generate a Table with the entries `data` and `country`,
where `data` contains the values of `cube1` and `country` the values of `cube2`. The cubes
are matched and broadcasted along their axes like in `mapCube`.

In addition, one can specify
`include_axes=(ax1,ax2...)` when one wants to include the values of certain axes in the table. For example
the command `(CubeTable(tair=cube1 axes=("lon","lat","time"))` would produce an iterator over a data structure
with entries `tair`, `lon`, `lat` and `time`.
"""
function CubeTable(; include_axes = (), expandaxes = (), cubes...)
    c = (map((k, v) -> v, keys(cubes), values(cubes))...,)
    all(i -> isa(i, Union{YAXArray,AbstractArray}), c) ||
        throw(ArgumentError("All inputs must be DataCubes"))
    varnames = map(string, keys(cubes))
    expandaxes = isa(expandaxes, Tuple) ? expandaxes : (expandaxes,)
    inaxnames = Set{String}()
    indims = if isempty(expandaxes)
        map(i -> InDims(), c)
    else
        map(c) do i
            axn = filter(collect(expandaxes)) do ax
                findAxis(ax, i) !== nothing
            end
            foreach(j -> push!(inaxnames, axname(getAxis(j, i))), axn)
            InDims(axn...)
        end
    end
    axnames = map(i -> axname.(caxes(i)), c)
    if isempty(intersect(axnames...))
        @warn "Input cubes to the table do not share a common axis, please check the axis names"
    end
    allvars = union(axnames...)
    allnums = collect(1:length(allvars))

    configiter =
        mapCube(identity, c, debug = true, indims = indims, outdims = (), ispar = false)
    # if inax !== nothing
    #   linax = length(inax)
    #   foreach(configiter.incubes) do ic1
    #     if !isempty(ic1.axesSmall)
    #       empty!(ic1.axesSmall)
    #       map!(i->i+1,ic1.loopinds,ic1.loopinds)
    #       pushfirst!(ic1.loopinds,1)
    #     else
    #       map!(i->i+1,ic1.loopinds,ic1.loopinds)
    #     end
    #   end
    # end
    r = collect(
        distributeLoopRanges(
            (configiter.loopcachesize...,),
            (map(length, configiter.LoopAxes)...,),
            getchunkoffsets(configiter),
        ),
    )
    ci = CubeIterator(configiter, r, include_loopvars = include_axes, varnames = varnames)
end

import Tables
Tables.istable(::Type{<:CubeIterator}) = true
Tables.rowaccess(::Type{<:CubeIterator}) = true
Tables.rows(x::CubeIterator) = x
Tables.schema(x::CubeIterator) = Tables.schema(typeof(x))
Tables.schema(x::Type{<:CubeIterator}) = Tables.Schema(
    fieldnames(eltype(x)),
    map(s -> fieldtype(eltype(x), s), fieldnames(eltype(x))),
)

Tables.istable(::Type{<:YAXArray}) = true
Tables.rowaccess(::Type{<:YAXArray}) = true
Tables.rows(x::YAXArray) = CubeTable(value = x, include_axes = true)
Tables.schema(x::YAXArray) = Tables.schema(Tables.rows(x))


text:
using OpenCV
using Base.Test

Base.include(@__MODULE__, "./test/jl/tests.jl"))


text:
module LazyCSV

import Mmap
import Parsers
import WeakRefStrings: WeakRefString

include("io.jl")
include("util.jl")

const SizeType = UInt64
const ZERO_UINT8 = UInt8(0)
const ZERO = zero(SizeType)
const ONE = one(SizeType)

include("buffered_vector.jl")
include("file.jl")
include("consumer.jl")

# the default value that determines whether the CSV parser should assume that the CSV file
# has a header not. Currently, by default it's assumed that there's no header in the CSV
# file if not otherwise stated.
const DEFAULT_HEADER_EXISTS = false

"""
    read_csv_line!(f::File)

read/parses a single line of input and returns it.

If no line is remained, `nothing` will be returned.
Otherwise, a pair of file position (`FilePos`) and line string is returned.

Please note that the returned line string is a `WeakRefString`.
"""
function read_csv_line!(f::File)::Union{Nothing, Tuple{FilePos, WeakRefString{UInt8}}}
    next_line_start_byte_pos_in_lines_buff = ZERO
    empty!(f.recs)
    resize!(f.recs, ONE)
    res = read_csv_line!(f.input, f.lines_buff, next_line_start_byte_pos_in_lines_buff, f.delim,
                   f.quotechar, f.escapechar, first_rec_fields(f), f.eager_parse_fields,
                   f.current_byte_pos)
    if res === nothing
        nothing
    else
        (completed_line, pos, line, _) = res
        @assert completed_line
        current_rec = first_rec(f)
        current_rec.pos = pos
        current_rec.line = line
        f.recs.size = ONE
        (pos, line)
    end
end

"""
A helper function to read/parse multiple lines of CSV files.
"""
function read_csv_lines!(f::File)
    # the byte index where the next line starts in the buffer
    next_line_start_byte_pos_in_lines_buff::SizeType = ZERO
    # number of lines read until now
    num_lines::SizeType = ZERO

    # clear out all the previous records
    empty!(f.recs)
    # the maximum number of lines to read/parse together
    max_lines = capacity(f.recs)

    # a flag that indicates whether we ran out of buffer in the middle of a line
    has_uncompleted_line = false

    while true
    	# if we read enough lines, we can stop for now
        num_lines >= max_lines && break

        # this will be the next line
        num_lines += ONE
        current_rec = f.recs[num_lines]
        res = read_csv_line!(f.input, f.lines_buff, next_line_start_byte_pos_in_lines_buff, f.delim,
                       f.quotechar, f.escapechar, current_rec.fields_buff, f.eager_parse_fields,
                       f.current_byte_pos)
        if res === nothing
        	# if there's no more line to read, we have to decrement the line counter that
        	# we have already incremented before.
            num_lines -= ONE
            break
        end
        (completed_line, pos, line, bytes_read) = res
        if completed_line
        	# if we were able to read a full line, then we can update the record information
        	# and prepare for reading the next line
            current_rec.pos = pos
            current_rec.line = line
            next_line_start_byte_pos_in_lines_buff += bytes_read
        else
        	# if we ran out of buffer in the middle of the run, we should stop for now and
        	# finish the line in the next batch
            has_uncompleted_line = true
            # this is done by reseting the input IO to the beginning of half-read line
            Parsers.fastseek!(f.input, pos.v)
            num_lines -= ONE
            break
        end
    end
    # fix the number of records
    f.recs.size = num_lines
    return has_uncompleted_line
end

@inline function private_add_field(num_fields_read, fields_len, fields_buff, buff_ptr, prev_field_index, num_bytes_read)
    num_fields_read += ONE
    if num_fields_read > fields_len
        # extend the fields buffer if required
        fields_len <<= ONE
        resize_vec!(fields_buff, fields_len)
    end
    @inbounds fields_buff[num_fields_read] = WeakRefString{UInt8}(buff_ptr+prev_field_index,
                                                                  num_bytes_read-prev_field_index-ONE)
    (num_fields_read, fields_len)
end

is_end_of_line_char(current_char::UInt8) = current_char == ASCII_NEWLINE ||
                                           current_char == ASCII_RETURN

"""
This function is at the heart of LazyCSV and it's actually an inlined parser-combinator
that knows how to parse a line of a CSV source.

It currently supports specifying:
  - `delim`: CSV delimiter character
  - `quotechar`: CSV quote character
  - `escapechar`: CSV escape character for escaping the quote and delimiter characters

In addition, by passing `eager_parse_fields=true`, then all the fields are also parsed
and stored as separate `WeakRefString`s. Otherwise, if `eager_parse_fields=false`, then
fields are not parsed and only a single line is extracted from the input (as a
`WeakRefString`).

This function returns a `Tuple{FilePos, WeakRefString{UInt8}}`. The first element of
this tuple contains the byte position of the starting point of this record (in the
source) and the second element is the record line content.

Note: the logic is a little bit complex, but being a core function, it's intended to
      be the most efficient code one can write. Any optimization is welcome!
"""
function read_csv_line!(s::IO, buff::Vector{UInt8}, buff_start::SizeType, delim::UInt8,
                        quotechar::UInt8, escapechar::UInt8,
                        fields::BufferedVector{WeakRefString{UInt8}},
                        eager_parse_fields::Bool, pos::MutFilePos
                       )::Union{Nothing, Tuple{Bool, FilePos, WeakRefString{UInt8}, SizeType}}
    prev_pos = pos.v
    current_pos = prev_pos
    # take the pointer to buffer only once
    buff_ptr::Ptr{UInt8} = pointer(buff) + buff_start
    # cache the buffer length
    buff_len::SizeType = SizeType(length(buff))
    buff_len_remained::SizeType = buff_len - buff_start

    # accumulates the number of bytes read for the current line
    # (it won't count the leading whitespaces)
    num_bytes_read::SizeType = ZERO

    # the buffer for storing fields (if `eager_parse_fields` is true)
    fields_buff = fields.buff
    # cache the length of fields buffer
    fields_len::SizeType = SizeType(length(fields_buff))
    # accumulates the number of fields read in this line (if `eager_parse_fields` is true)
    num_fields_read::SizeType = ZERO
    # stores the index of the separator for the previous field.
    prev_field_index::SizeType = ZERO

    # `inside_quote` determines whether the given fields is quoted and we are scanning
    # inside the quote this variable is initialized here, as there might be leading
    # whitespaces inside the quote
    inside_quote::Bool = false

    is_escaped::Bool = false

    same_quote_and_escape::Bool = quotechar == escapechar

    # this first while-loop is for bypassing the empty lines or the lines with only whitespaces
    while num_bytes_read == ZERO && !eof(s)
        # stores the next byte value
        # it's important that in the UTF-8 format, all non-ASCII characters start with 1-bit
        # on the leftmost (in all their parts) and it makes it easy to find special characters
        # (e.g., new-line or ASCII delimiters) using the byte value.
        current_char::UInt8 = ZERO_UINT8

        # we look for a new-line character or end-of-file
        while !is_end_of_line_char(current_char)
            if num_bytes_read >= buff_len_remained
                return (false, FilePos(prev_pos), WeakRefString{UInt8}(buff_ptr, num_bytes_read), num_bytes_read)
            end

            prev_char::UInt8 = current_char
            # use a new-line character if we've already reached end of the file
            current_char = eof(s) ? ASCII_NEWLINE : Parsers.readbyte(s)
            current_pos += 1

            # skip the leading whitespaces if not inside a quote
            num_bytes_read == ZERO && !inside_quote && iswhitespace(current_char) && continue

            # handles the quote character
            # there are two cases:
            if same_quote_and_escape
                # case 1: if both quotechar and escapechar are the same
                if current_char == quotechar
                    if inside_quote
                        # if we are already inside a quote, then this is either an end-quote or
                        # a escape character and we have to escape it anyways
                        inside_quote = false
                        continue
                    else
                        inside_quote = true
                        # if we were not in a quote, it could be because of a escape quote
                        # before this one and the following check does the job. If the prev.
                        # character was also a `quotechar`, then we have to write the character.
                        prev_char != quotechar && continue
                    end
                end
            else
                # case 2: if quotechar and escapechar are different
                if is_escaped
                    is_escaped = false
                    if !(current_char == quotechar ||  current_char == escapechar)
                        @warn("Malformed input string: escaped '\\$(Char(current_char))' is not supported.")
                    end
                elseif current_char == escapechar
                    is_escaped = true
                    continue
                elseif current_char == quotechar
                    inside_quote = !inside_quote
                    continue
                end
            end

            num_bytes_read += ONE
            if eager_parse_fields && !inside_quote && current_char == delim
                (num_fields_read, fields_len) = private_add_field(num_fields_read, fields_len,
                                                                  fields_buff, buff_ptr,
                                                                  prev_field_index, num_bytes_read)
                prev_field_index = num_bytes_read
            end
            unsafe_store!(buff_ptr + num_bytes_read - ONE, current_char)
        end
    end

    if num_bytes_read > ZERO
        # handle the last field if anything is read at all
        if eager_parse_fields

            (num_fields_read, fields_len) = private_add_field(num_fields_read, fields_len,
                                                              fields_buff, buff_ptr,
                                                              prev_field_index, num_bytes_read)
            # correctly set the number of fields read to `num_fields_read`
            fields.size = num_fields_read
        end
        pos.v = current_pos
        last_char = unsafe_load(buff_ptr, num_bytes_read - ONE)
        num_chars_in_line = !is_end_of_line_char(last_char) ? num_bytes_read - ONE : num_bytes_read
        (true, FilePos(prev_pos), WeakRefString{UInt8}(buff_ptr, num_chars_in_line), num_bytes_read)
    else
        # if fields were supposed to be parsed eagerly, then correctly set the size of buffer to zero
        eager_parse_fields && (fields.size = ZERO)
        # if nothing is read, then return `nothing`
        nothing
    end
end

"""
    csvread(input::IO, delim=','; <arguments>...)
Read CSV from `file`. Returns a tuple of 2 elements:
1. A tuple of columns each either a `Vector`, or `StringArray`
2. column names if `header_exists=true`, empty array otherwise
# Arguments:
- `input`: an input file path or an IO object
- `delim`: the delimiter character
- `lazy`: if `true`, then an instance of `LazyCSV.File` is returned.
          Otherwise, if `false`, the result gets materialized and then returned.
- `quotechar`: character used to quote strings, defaults to `"`
- `escapechar`: character used to escape quotechar in strings. (could be the same as quotechar)
- `header_exists`: boolean specifying whether CSV file contains a header
- `num_rows`: the number of rows to be loaded at once. These rows will be processed
              column-at-a-time.
- `lines_buff_len`: the buffer size (in bytes) used for storing lines to be parsed
"""
function csvread(input::Union{IO,AbstractString}; delim::Union{Char, Nothing}=nothing,
                 lazy=true, eager_parse_fields=DEFAULT_EAGER_PARSE_FIELDS,
                 quotechar::Char=DEFAULT_QUOTE, escapechar::Char=quotechar,
                 header_exists::Bool=DEFAULT_HEADER_EXISTS,
                 lines_buff_len::SizeType=DEFAULT_LINE_LEN,
                 num_rows::SizeType=DEFAULT_NUM_ROWS, kw...)
    input_io = read_mmap_data(input)

    if delim === nothing
        if isa(input, AbstractString)
            if endswith(input, ".tsv")
                delim = '\t'
            elseif endswith(input, ".wsv")
                delim = ' '
            elseif endswith(input, ".tbl")
                delim = '|'
            else
                delim = ','
            end
        else
            delim = ','
        end
    end

    file = File(input_io; delim=delim, eager_parse_fields=eager_parse_fields,
                quotechar=quotechar, escapechar=escapechar, header_exists=header_exists,
                fields_buff_len=DEFAULT_NUM_FIELDS,
                lines_buff_len=lines_buff_len, num_rows=num_rows)

    lazy ? file : materialize(file)
end

export csvread

end # module


text:
function demoBurglar()

    burglar,earthquake,alarm,radio=1,2,3,4 # Variable order is arbitary
    yes,no= 1,2 # define states, starting from 1

    # This defines the meaning of the variables and their states.
    # It is not necessary for the computation, but is useful for
    # printing the results in a meaningful form:
    DictVariable=Dict{Integer,DiscreteVariable}()
    DictVariable[burglar]=DiscreteVariable("burglar",["yes","no"])
    DictVariable[earthquake]=DiscreteVariable("earthquake",["yes","no"])
    DictVariable[alarm]=DiscreteVariable("alarm",["yes","no"])
    DictVariable[radio]=DiscreteVariable("radio",["yes","no"])

    # Define the discrete probability tables:
    pot=Array{PotArray}(4) # there are 4 discrete potentials (arrays)
    pot[burglar]=PotArray(burglar,[0.01 0.99])
    pot[earthquake]=PotArray(earthquake,[0.000001 1-0.000001])

    table=zeros(2,2)
    table[yes,yes]=1
    table[no,yes]=0
    table[yes,no]=0
    table[no,no]=1
    pot[radio]=PotArray([radio earthquake],table)

    table=zeros(2,2,2)
    table[yes,yes,yes]=0.9999
    table[yes,yes,no]=0.99
    table[yes,no,yes]=0.99
    table[yes,no,no]=0.0001
    table[no,:,:]=1-table[yes,:,:]
    pot[alarm]=PotArray([alarm burglar earthquake],table)

    jointpot=prod(pot)

    println("p(burglar|alarm=yes):")
    show(condpot(setpot(jointpot,alarm,yes),burglar),DictVariable)

    println("p(burglar|alarm=yes,radio=yes):")
    show(condpot(setpot(jointpot,[alarm radio],[yes yes]),burglar),DictVariable)

    L=dag(pot)
    #try
        # This is just for plotting the DAG:
        DictPlot=Dict{Integer,VariablePlot}()
        DictPlot[burglar]=VariablePlot("Burglar",x=0,y=1,nodesize=0.3)
        DictPlot[earthquake]=VariablePlot("Earthquake",x=3,y=1,nodesize=0.35)
        DictPlot[alarm]=VariablePlot("Alarm",x=0,y=0,nodesize=0.2)
        DictPlot[radio]=VariablePlot("Radio",x=3,y=0,nodesize=0.2)

        PlotGraph(DictPlot,L)

        #PlotGraph(PlaceVertices(DictPlot,L,scale=5),L,arrowsize=0.2) # automatic node placement

    #end

end


text:
"""
    write(C::ClimGrid, filename::String)

Write to disk ClimGrid C to netCDF file. Still experimental, some attributes are hardcoded.
"""
function write(C::ClimGrid, filename::String)

    # Test extension
    if extension(filename) .!= ".nc"
        filename = string(filename, ".nc")
    end

    # This creates a new NetCDF file
    # The mode "c" stands for creating a new file (clobber)
    ds = Dataset(filename, "c")

    latsymbol, lonsymbol = ClimateTools.getsymbols(C)
    x, y, timevec = ClimateTools.getdims(C)
    longrid, latgrid = ClimateTools.getgrids(C)

    # DIMENSIONS
    defDim(ds, string(lonsymbol),length(x))
    defDim(ds, string(latsymbol), length(y))
    defDim(ds, "time", length(timevec))

    nclon = defVar(ds, string(lonsymbol), Float64, (string(lonsymbol),))
    nclon.attrib["units"] = C.lonunits#"degrees_east"
    nclon.attrib["long_name"] = "longitude"
    nclon.attrib["standard_name"] = "longitude"
    nclon.attrib["actual_range"] = [minimum(x), maximum(x)]

    nclat = defVar(ds, string(latsymbol), Float64, (string(latsymbol),))
    nclat.attrib["units"] = C.latunits#"degrees_north"
    nclat.attrib["long_name"] = "latitude"
    nclat.attrib["standard_name"] = "latitude"
    nclat.attrib["actual_range"] = [minimum(y), maximum(y)]

    if latsymbol != :lat # a lon-lat grid is needed
        ncrlat = defVar(ds, "lat", Float64, (string(lonsymbol), string(latsymbol)))
        ncrlat.attrib["long_name"] = "latitude"
        ncrlat.attrib["units"] = "degrees"
        ncrlat.attrib["standard_name"] = "grid_latitude"
        ncrlat.attrib["axis"] = "Y"
        ncrlat.attrib["coordinate_defines"] = "point"
        ncrlat.attrib["actual_range"] = [minimum(latgrid), maximum(latgrid)]
    end

    if lonsymbol != :lon
        ncrlon = defVar(ds,"lon", Float64, (string(lonsymbol), string(latsymbol)))
        ncrlon.attrib["long_name"] = "longitude"
        ncrlon.attrib["units"] = "degrees"
        ncrlon.attrib["standard_name"] = "grid_longitude"
        ncrlon.attrib["axis"] = "X"
        ncrlon.attrib["coordinate_defines"] = "point"
        ncrlon.attrib["actual_range"] = [minimum(longrid), maximum(longrid)]
    end

    if ClimateTools.@isdefined C.grid_mapping["grid_mapping"]
        ncmapping = defVar(ds, C.grid_mapping["grid_mapping"], Char, ())
    elseif ClimateTools.@isdefined C.grid_mapping["grid_mapping_name"]
        ncmapping = defVar(ds, C.grid_mapping["grid_mapping_name"], Char, ())
    end
    for iattr in keys(C.grid_mapping)
        ncmapping.attrib[iattr] = C.grid_mapping[iattr]
    end

    # time
    timeout = Array{Float64, 1}(undef, length(timevec))
    for itime = 1:length(timevec)
        timeout[itime] = NCDatasets.timeencode([timevec[itime]], C.timeattrib["units"], C.timeattrib["calendar"])[1]
    end

    # daysfrom = string("days since ", string(timevec[1] - Dates.Day(1))[1:10], " ", string(timevec[1])[12:end])

    nctime = defVar(ds,"time", Float64, ("time",))
    nctime.attrib["long_name"] = "time"
    nctime.attrib["standard_name"] = "time"
    nctime.attrib["axis"] = "T"
    nctime.attrib["calendar"] = C.typeofcal#"365_day"
    nctime.attrib["units"] = C.timeattrib["units"]#daysfrom
    # nctime.attrib["bounds"] = "time_bnds"
    nctime.attrib["coordinate_defines"] = "point"

    # Define the variables contained in ClimGrid C
    v = defVar(ds, C.variable, Float32, (string(lonsymbol),string(latsymbol), "time"))

    # write attributes
    for iattr in keys(C.varattribs)
        v.attrib[iattr] = C.varattribs[iattr]
    end

    # Define global attributes
    for iattr in keys(C.globalattribs)
        ds.attrib[iattr] = C.globalattribs[iattr]
    end

    v[:] = C[1].data

    # Time vector
    nctime[:] = timeout

    # Dimensions
    nclon[:] = x
    nclat[:] = y

    # Dimensions
    if latsymbol != :lat
        ncrlon[:] = longrid
    end
    if latsymbol != :lat
        ncrlat[:] = latgrid
    end

    # Close file
    close(ds)

    return true

end


text:
# Encode MNIST images as compressed vectors that can later be decoded back into
# images.
using Flux, Flux.Data.MNIST
using Flux: @epochs, onehotbatch, mse, throttle
using Base.Iterators: partition
using Parameters: @with_kw
using CUDAapi
if has_cuda()
    @info "CUDA is on"
    import CuArrays
    CuArrays.allowscalar(false)
end

@with_kw mutable struct Args
    lr::Float64 = 1e-3		# Learning rate
    epochs::Int = 10		# Number of epochs
    N::Int = 32			# Size of the encoding
    batchsize::Int = 1000	# Batch size for training
    sample_len::Int = 20 	# Number of random digits in the sample image
    throttle::Int = 5		# Throttle timeout
end

function get_processed_data(args)
    # Loading Images
    imgs = MNIST.images()
    #Converting image of type RGB to float
    imgs = channelview.(imgs)
    # Partition into batches of size 1000
    train_data = [float(hcat(vec.(imgs)...)) for imgs in partition(imgs, args.batchsize)]

    train_data = gpu.(train_data)
    return train_data
end

function train(; kws...)
    args = Args(; kws...)

    train_data = get_processed_data(args)

    @info("Constructing model......")
    # You can try to make the encoder/decoder network larger
    # Also, the output of encoder is a coding of the given input.
    # In this case, the input dimension is 28^2 and the output dimension of
    # encoder is 32. This implies that the coding is a compressed representation.
    # We can make lossy compression via this `encoder`.
    encoder = Dense(28^2, args.N, leakyrelu) |> gpu
    decoder = Dense(args.N, 28^2, leakyrelu) |> gpu

    # Defining main model as a Chain of encoder and decoder models
    m = Chain(encoder, decoder)

    @info("Training model.....")
    loss(x) = mse(m(x), x)
    ## Training
    evalcb = throttle(() -> @show(loss(train_data[1])), args.throttle)
    opt = ADAM(args.lr)

    @epochs args.epochs Flux.train!(loss, params(m), zip(train_data), opt, cb = evalcb)

    return m, args
end

using Images

img(x::Vector) = Gray.(reshape(clamp.(x, 0, 1), 28, 28))

function sample(m, args)
    imgs = MNIST.images()
    #Converting image of type RGB to float
    imgs = channelview.(imgs)
    # `args.sample_len` random digits
    before = [imgs[i] for i in rand(1:length(imgs), args.sample_len)]
    # Before and after images
    after = img.(map(x -> cpu(m)(float(vec(x))), before))
    # Stack them all together
    hcat(vcat.(before, after)...)
end

cd(@__DIR__)
m, args= train()
# Sample output
@info("Saving image sample as sample_ae.png")
save("sample_ae.png", sample(m, args))

max_stars_repo_name: SubhadityaMukherjee/model-zoo

max_issues_repo_name: SubhadityaMukherjee/model-zoo

max_forks_repo_name: SubhadityaMukherjee/model-zoo


text:
include("galis_example.jl")

GAMMA_0 = [bet    0     0  0;
           1      sig   0  0;
           0      0     0  0;
           0      0     0  1]

GAMMA_1 = [1      -kappa  0  0;
           0       sig    1  0;
           -phi_pi  -phi_y  1 -1;
           0       0      0  rho_v]

PSI = [0; 0; 0; 1]

PI = [bet  0;
      1    sig;
      0    0;
      0    0]

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
using ExprRules, Statistics
using BenchmarkTools
using Random

protectedDiv(x, y) = iszero(y) ? 1.0 : x / y

const grammar = @cfgrammar begin
    R = x
    R = R + R
    R = R - R
    R = R * R
    R = protectedDiv(R, R)
    R = -R
    R = cos(R)
    R = sin(R)
    R = |(1.0:5.0)
end
const S = SymbolTable(grammar, Main)

gt(x) = x^4 + x^3 + x^2 + x
function loss1(tree::RuleNode, grammar::Grammar)
    ex = get_executable(tree, grammar)
    los = 0.0
    global x
    for y = -1.0:0.1:1.0
        x = y
        los += abs2(eval(ex) - gt(x))
    end
    mean(los)
end
function loss2(tree::RuleNode, grammar::Grammar)
    ex = get_executable(tree, grammar)
    los = 0.0
    for y = -1.0:0.1:1.0
        S[:x] = y
        los += abs2(Core.eval(S, ex) - gt(y))
    end
    mean(los)
end

function run_benchmark(seed::Int=1)
    Random.seed!(seed)
    tree = rand(RuleNode, grammar, :R, 10)
    ex = get_executable(tree, grammar)
    @show ex
    @btime loss1($tree, $grammar)
    @btime loss2($tree, $grammar)
end


text:
@auto_hash_equals type OPCSPState <: State
    i::Int
    open::IntSet # unvisited
    remaining::Float64
    d::Vector{Float64}
end
# ==(u::OPCSPState, v::OPCSPState) = u.i==v.i && u.v==v.v && u.d==v.d && u.remaining==v.remaining
# ==(u::OPCSPState, v::OPCSPState) = error("== is being used for a state") # want to see if this is being used
# hash(s::OPCSPState) = hash(s.i, hash(s.v, hash(s.d, hash(s.remaining))))
create_transition_distribution(problem::OPCSP) = OPCSPState(0, IntSet(), 0.0, Array(Float64, length(problem)))
create_state(op::OPCSP) = OPCSPState(0, IntSet(), 0.0, Array(Float64, length(op)))
get_d(s::OPCSPState) = s.d

@auto_hash_equals type OPCSPAction <: Action
    next::Int
    OPCSPAction() = new()
    OPCSPAction(i::Int) = new(i)
end
# ==(u::OPCSPAction, v::OPCSPAction) = u.next==v.next
# hash(a::OPCSPAction) = hash(a.next)
create_action(::OPCSP) = OPCSPAction()

@auto_hash_equals type OPCSPObs <: Observation
    d::Float64
    OPCSPObs() = new()
    OPCSPObs(d::Float64) = new(d)
end
create_observation_distribution(::OPCSP) = OPCSPObs()
create_observation(::OPCSP) = OPCSPObs()

type OPCSPDistribution <: Belief
    i::Int
    open::IntSet # unvisited
    remaining::Float64
    dist::MVN
end
function rand!(rng::AbstractRNG, s::OPCSPState, d::OPCSPDistribution)
    s.i = d.i
    copy!(s.open, d.open)
    rand!(rng, s.d, d.dist, robust=true)
    s.remaining = d.remaining
    return s
end
create_belief(problem::OPCSP) = OPCSPDistribution(
    0,
    IntSet(),
    0.0,
    MVN(Array(Float64,length(problem)),Array(Float64, length(problem), length(problem))))
type OPCSPUpdater <: BeliefUpdater
    problem::OPCSP
end
function update(u::OPCSPUpdater, bo::OPCSPDistribution, a::OPCSPAction, o::OPCSPObs, bn::OPCSPDistribution=create_belief(u))
    bn.i = a.next
    copy!(bn.open, bo.open)
    delete!(bn.open, a.next)
    bn.remaining = bo.remaining-u.problem.distances[bo.i, a.next]
    bn.dist = apply_measurement(bo.dist, a.next, o.d) #XXX creating a copy here seems bad maybe
    return bn
end
updater(p::OPCSP) = OPCSPUpdater(p)
create_belief(updater::OPCSPUpdater) = create_belief(updater.problem)
initial_belief(p::OPCSP) = OPCSPDistribution(
    p.start,
    delete!(IntSet(1:length(p)), p.start),
    p.distance_limit,
    MVN(zeros(length(p)), p.covariance)
    )
get_d(d::OPCSPDistribution) = d.dist.mean
Base.mean(d::OPCSPDistribution) = OPCSPState(

# this could be inefficient because we are copying too much
# note that state transitions are deterministic
function transition(problem::OPCSP, s::OPCSPState, a::OPCSPAction, sp::OPCSPState=create_transition_distribution(problem))
    sp.i = a.next
    copy!(sp.open, s.open) # if performance is an issue, then change this to = and move delete! to rand!
    delete!(sp.open, a.next)
    sp.d = s.d
    sp.remaining = s.remaining - problem.distances[s.i, a.next]
    return sp
end
function rand!(::AbstractRNG, s::OPCSPState, sd::OPCSPState)
    s.i = sd.i
    copy!(s.open, sd.open)
    s.d = sd.d
    s.remaining = sd.remaining
    return s
end

reward(op::OPCSP, s::OPCSPState, ::OPCSPAction, ::OPCSPState) = op.r[s.i]+s.d[s.i]
discount(op::OPCSP) = 1.0

function observation(op::OPCSP, s::OPCSPState, a::OPCSPAction, sp::OPCSPState, od::OPCSPObs=create_observation_distribution(op))
    od.d = sp.d[a.next]
    return od
end
function rand!(rng::AbstractRNG, sample::OPCSPObs, od::OPCSPObs)
    sample.d = od.d
    return sample
end
isterminal(op::OPCSP, s::OPCSPState) = op.stop == s.i
discount(op::OPCSP) = 1.0


text:
@testset "sauvola" begin
    @info "Test: Sauvola"

    @testset "API" begin
        img_gray = imresize(testimage("lena_gray_256"); ratio=0.25)
        img = copy(img_gray)

        # Sauvola
        @test Sauvola() == Sauvola(7, 0.2)
        @test Sauvola(7, 0.2) == Sauvola(window_size=7, bias=0.2)

        # window_size non-positive integer
        @test_throws ArgumentError Sauvola(window_size = -10)
        @test_throws TypeError Sauvola(window_size = 32.5)

        # binarize
        f = Sauvola(window_size=7, bias=0.2)
        binarized_img_1 = binarize(img, f)
        @test img == img_gray # img unchanged
        @test eltype(binarized_img_1) == Gray{N0f8}

        binarized_img_2 = binarize(Gray{Bool}, img, f)
        @test img == img_gray # img unchanged
        @test eltype(binarized_img_2) == Gray{Bool}

        binarized_img_3 = similar(img, Bool)
        binarize!(binarized_img_3, img, f)
        @test img == img_gray # img unchanged
        @test eltype(binarized_img_3) == Bool

        binarized_img_4 = copy(img_gray)
        binarize!(binarized_img_4, f)
        @test eltype(binarized_img_4) == Gray{N0f8}

        @test binarized_img_1 == binarized_img_2
        @test binarized_img_1 == binarized_img_3
        @test binarized_img_1 == binarized_img_4
    end

    @testset "Types" begin
        # Gray
        img_gray = imresize(testimage("lena_gray_256"); ratio=0.25)
        f = Sauvola(window_size=7, bias=0.2)

        type_list = generate_test_types([Float32, N0f8], [Gray])
        for T in type_list
            img = T.(img_gray)
            @test_reference "References/Sauvola_Gray.png" Gray.(binarize(img, f))
        end

        # Color3
        img_color = imresize(testimage("lena_color_256"); ratio=0.25)
        f = Sauvola(window_size=7, bias=0.2)

        type_list = generate_test_types([Float32, N0f8], [RGB, Lab])
        for T in type_list
            img = T.(img_gray)
            @test_reference "References/Sauvola_Color3.png" Gray.(binarize(img, f))
        end
    end

    @testset "Numerical" begin
        # Check that the image only has ones or zeros.
        img = imresize(testimage("lena_gray_256"); ratio=0.25)
        f = Sauvola(window_size=7, bias=0.2)
        img₀₁ = binarize(img, f)
        non_zeros = findall(x -> x != 0.0 && x != 1.0, img₀₁)
        @test length(non_zeros) == 0

        # Check that ones and zeros have been assigned to the correct side of the threshold.
        maxval, maxpos = findmax(Gray.(img))
        @test img₀₁[maxpos] == 1
        minval, minpos = findmin(Gray.(img))
        @test img₀₁[minpos] == 0

        for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
            img = T.([i <= 25 && j <= 25 ? 0.8 : 1.0 for i = 1:50, j = 1:50])
            target_row = target_col = 13
            img[target_row,target_col] = 0

            for i in 0:10:50, j in 0:10:50
                img₀ = circshift(img, (i,j))
                target_row₀ = (target_row + i) % 50
                target_col₀ = (target_col + j) % 50

                img_bin = binarize(img₀, Sauvola(window_size = 7, bias = 0.21))
                @test sum(img_bin .== 0) == 1
                @test img_bin[target_row₀, target_col₀] == 0
            end
        end
    end

end


text:
mutable struct CTM <: TopicModel
	K::Int
	M::Int
	V::Int
	N::Vector{Int}
	C::Vector{Int}
	corp::Corpus
	topics::VectorList{Int}
	mu::Vector{Float64}
	sigma::Matrix{Float64}
	beta::Matrix{Float64}
	lambda::VectorList{Float64}
	vsq::VectorList{Float64}
	lzeta::Float64
	phi::Matrix{Float64}
	invsigma::Matrix{Float64}
	newbeta::Matrix{Float64}
	elbo::Float64
	newelbo::Float64

	function CTM(corp::Corpus, K::Integer)
		@assert ispositive(K)
		@assert !isempty(corp)
		checkcorp(corp)

		M, V, U = size(corp)
		N = [length(doc) for doc in corp]
		C = [size(doc) for doc in corp]

		topics = [collect(1:V) for _ in 1:K]

		mu = zeros(K)
		sigma = eye(K)
		beta = rand(Dirichlet(V, 1.0), K)'
		lambda = [zeros(K) for _ in 1:M]
		vsq = [ones(K) for _ in 1:M]
		lzeta = 0.5
		phi = ones(K, N[1]) / K

		model = new(K, M, V, N, C, copy(corp), topics, mu, sigma, beta, lambda, vsq, lzeta, phi)
		fixmodel!(model, check=false)

		for d in 1:M
			model.phi = ones(K, N[d]) / K
			updateNewELBO!(model, d)
		end
		model.phi = ones(K, N[1]) / K
		updateELBO!(model)
		return model
	end
end

function Elogpeta(model::CTM, d::Int)
	x = 0.5 * (logdet(model.invsigma) - model.K * log(2pi) - dot(diag(model.invsigma), model.vsq[d]) - dot(model.lambda[d] - model.mu, model.invsigma * (model.lambda[d] - model.mu)))
	return x
end

function Elogpz(model::CTM, d::Int)
	counts = model.corp[d].counts
	x = dot(model.phi' * model.lambda[d], counts) + model.C[d] * model.lzeta
	return x
end

function Elogpw(model::CTM, d::Int)
	terms, counts = model.corp[d].terms, model.corp[d].counts
	x = sum(model.phi .* log.(@boink model.beta[:,terms]) * counts)
	return x
end

function Elogqeta(model::CTM, d::Int)
	x = -entropy(MvNormal(model.lambda[d], diagm(model.vsq[d])))
	return x
end

function Elogqz(model::CTM, d::Int)
	counts = model.corp[d].counts
	x = -sum([c * entropy(Categorical(model.phi[:,n])) for (n, c) in enumerate(counts)])
	return x
end

function updateELBO!(model::CTM)
	model.elbo = model.newelbo
	model.newelbo = 0
	return model.elbo
end

function updateNewELBO!(model::CTM, d::Int)
	model.newelbo += (Elogpeta(model, d)
					+ Elogpz(model, d)
					+ Elogpw(model, d)
					- Elogqeta(model, d)
					- Elogqz(model, d))
end

function updateMu!(model::CTM)
	model.mu = sum(model.lambda) / model.M
end

function updateSigma!(model::CTM)
	model.sigma = diagm(sum(model.vsq)) / model.M + Base.covm(hcat(model.lambda...), model.mu, 2, false)
	model.invsigma = inv(model.sigma)
end

function updateBeta!(model::CTM)
	model.beta = model.newbeta ./ sum(model.newbeta, 2)
	model.newbeta = zeros(model.K, model.V)
end

function updateNewBeta!(model::CTM, d::Int)
	terms, counts = model.corp[d].terms, model.corp[d].counts
	model.newbeta[:,terms] += model.phi .* counts'
end

function updateLambda!(model::CTM, d::Int, niter::Integer, ntol::Real)
	"Newton's method."

	counts = model.corp[d].counts
	for _ in 1:niter
		lambdaGrad = model.invsigma * (model.mu - model.lambda[d]) + model.phi * counts - model.C[d] * exp.(model.lambda[d] + 0.5 * model.vsq[d] - model.lzeta)
		lambdaInvHess = -inv(eye(model.K) + model.C[d] * model.sigma * diagm(exp.(model.lambda[d] + 0.5 * model.vsq[d] - model.lzeta))) * model.sigma
		model.lambda[d] -= lambdaInvHess * lambdaGrad
		if norm(lambdaGrad) < ntol
			break
		end
	end
end

function updateVsq!(model::CTM, d::Int, niter::Integer, ntol::Real)
	"Newton's method."

	for _ in 1:niter
		rho = 1.0
		vsqGrad = -0.5 * (diag(model.invsigma) + model.C[d] * exp.(model.lambda[d] + 0.5 * model.vsq[d] - model.lzeta) - 1 ./ model.vsq[d])
		vsqInvHess = -1 ./ (0.25 * model.C[d] * exp.(model.lambda[d] + 0.5 * model.vsq[d] - model.lzeta) + 0.5 ./ model.vsq[d].^2)
		p = vsqInvHess .* vsqGrad

		while minimum(model.vsq[d] - rho * p) <= 0
			rho *= 0.5
		end
		model.vsq[d] -= rho * p

		if norm(vsqGrad) < ntol
			break
		end
	end
	@bumper model.vsq[d]
end

function updateLzeta!(model::CTM, d::Int)
	model.lzeta = logsumexp(model.lambda[d] + 0.5 * model.vsq[d])
end

function updatePhi!(model::CTM, d::Int)
	terms = model.corp[d].terms
	model.phi = addlogistic(log.(model.beta[:,terms]) .+ model.lambda[d], 1)
end

function train!(model::CTM; iter::Integer=150, tol::Real=1.0, niter::Integer=1000, ntol::Real=1/model.K^2, viter::Integer=10, vtol::Real=1/model.K^2, chkelbo::Integer=1)
	@assert all(.!isnegative.([tol, ntol, vtol]))
	@assert all(ispositive.([iter, niter, viter, chkelbo]))

	for k in 1:iter
		chk = (k % chkelbo == 0)
		for d in 1:model.M
			for _ in 1:viter
				oldlambda = model.lambda[d]
				updatePhi!(model, d)
				updateLzeta!(model, d)
				updateLambda!(model, d, niter, ntol)
				updateVsq!(model, d, niter, ntol)
				if norm(oldlambda - model.lambda[d]) < vtol
					break
				end
			end
			updateNewBeta!(model, d)
			chk && updateNewELBO!(model, d)
		end
		updateMu!(model)
		updateSigma!(model)
		updateBeta!(model)
		if checkELBO!(model, k, chk, tol)
			break
		end
	end
	updatePhi!(model, 1)
	updateLzeta!(model, 1)
	model.topics = [reverse(sortperm(vec(model.beta[i,:]))) for i in 1:model.K]
	nothing
end


text:
using ParallelUtilities
using Base.Test

# write your own tests here
include("testmacros.jl")


text:
name = "libLLVM"
version = v"12.0.1+1"

# Include common tools.
include("../common.jl")
build_tarballs(ARGS, configure_extraction(ARGS, version, name; experimental_platforms=true, assert=true)...; skip_audit=true, julia_compat="1.7")

max_stars_repo_name: Alexander-Barth/Yggdrasil

max_issues_repo_name: Alexander-Barth/Yggdrasil

max_forks_repo_name: Alexander-Barth/Yggdrasil


text:
using HomotopyContinuation
using DynamicPolynomials
using BenchmarkTools

function triangle()
    a = 5
    b = 4
    c = 3;

    @polyvar sθ cθ
    f1 = cθ^2 + sθ^2 - (1.0 + 0im)
    f2 = (a*cθ - b)^2 + (1.0 + 0im) * (a*sθ)^2 - c^2
    F = PolySystem([f1, f2])

    # lets use a total degree homotopy
    G, start_solutions = totaldegree(F)

    # we will use a simple straight line homotopy
    H = StraightLineHomotopy(G, F)

    algorithm = PredictorCorrector.Spherical()
    #now we can solve
    results = solve(H, start_solutions, algorithm)
    results = map(x -> x.solution, filter(x -> x.returncode == :Success, results))
    @assert (length(results) == 2)

    @benchmark solve($H, $start_solutions, $algorithm)
end

println("\n\nRun polysystem benchmark:")
triangle_benchmark = triangle()
show(STDOUT, MIME"text/plain"(), triangle_benchmark)

println("\nDone.")


text:
using BinaryBuilder, Pkg

# abc
name = "ABC"
version = v"1.01.0"
ABC_ver = "1.01.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/berkeley-abc/abc.git", "645752f7d651b07a0ebd7cac827d6ce5dd5bd266")
]

# ABC is setup for native builds (the arch_flags program and running the
# shell to get other flags), so we need to manually add the arch flags
# here.  Also, the combo of BinaryBuilder and ABC's Makefile always seems
# to override rather than extend CFLAGS, so I've included the whole mess
# here.  Even though we use ABC_USE_STDINT_H=1 it doesn't work as
# expected, so we only use it here to prevent the Makefile from running
# the arch_flags executable, which will always fail when cross compiling.
# Note that neither -DLIN nor -DLIN64 actually appear to be used anymore but
# still set here anyway

# Bash recipe for building across all platforms
script = raw"""

EXFLGS="-Wall -Wno-unused-function -Wno-write-strings -Wno-sign-compare -DSIZEOF_INT=4 -DABC_USE_CUDD=1 -DABC_USE_PTHREADS -Wno-unused-but-set-variable -fPIC"

if [[ ${nbits} == 32 ]]; then
    EXFLGS="${EXFLGS} -DSIZEOF_VOID_P=4"
    EXFLGS="${EXFLGS} -DSIZEOF_LONG=4"
    EXFLGS="${EXFLGS} -DLIN"
else
    EXFLGS="${EXFLGS} -DSIZEOF_VOID_P=8"
    EXFLGS="${EXFLGS} -DSIZEOF_LONG=8"
    EXFLGS="${EXFLGS} -DLIN64"
fi

# musl doesn't like readline
if [[ ${target} != *"musl"* ]]; then
    EXFLGS="${EXFLGS} -DABC_READLINE"
fi

cd ${WORKSPACE}/srcdir/abc

# it appears necessary to have the file(s) for the FileProduct(s) live in the destdir
cp abc.rc ${prefix}/.

make -j${nproc} ABC_USE_STDINT_H=1 CFLAGS+="${EXFLGS}" libabc.so
# the abc Makefile always makes a .so  Fix that here.
mv libabc.so libabc.${dlext}
mkdir -p "${libdir}"
cp libabc.${dlext} ${libdir}

make -j${nproc} ABC_USE_STDINT_H=1 CFLAGS+="${EXFLGS}"
mkdir -p "${bindir}"
cp abc${exeext} ${bindir}
chmod +x ${bindir}/*
install_license copyright.txt
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(exclude= x -> (Sys.iswindows(x) || Sys.isfreebsd(x)))

# The products that we will ensure are always built
products = [
    ExecutableProduct("abc", :abc),
    LibraryProduct("libabc", :libabc),
    FileProduct("abc.rc", :abc_rc),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Readline_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function _funccall(func::PySliceType, args, ::Nothing)
    a, b, c = args
    PySlice(a, b, c)
end

function _funccall(func::PyRangeType, args, ::Nothing)
    a, b, c = args
    PyRange(a, b, c)
end

function _funccall(func::PyComplexType, args, ::Nothing)
    a, b = args
    PyComplex(a, b)
end

function _funccall(func::PyByteArrayType, args, ::Nothing)
    if length(args) == 0
        PyByteArray()
    elseif length(args) == 1
        a = args[1]::PyBytes
        PyByteArray(a.values)
    else
        error()
    end
end

function _funccall(func::PyArrayReconstructor, args, ::Nothing)
    # TODO: not sure what n is?
    # TODO: do i need to worry about endianness?
    cls, t, n, data = args
    (cls === PyArrayType() && t isa String && data isa PyBytes) || error()
    if t == "b"
        T = Cchar
    elseif t == "B"
        T = Cuchar
    elseif t == "u"
        T = Cwchar_t
    elseif t == "h"
        T = Cshort
    elseif t == "H"
        T = Cushort
    elseif t == "i"
        T = Cint
    elseif t == "I"
        T = Cuint
    elseif t == "l"
        T = Clong
    elseif t == "L"
        T = Clong
    elseif t == "q"
        T = Clonglong
    elseif t == "Q"
        T = Culonglong
    elseif t == "f"
        T = Cfloat
    elseif t == "d"
        T = Cdouble
    else
        error()
    end
    arr = reinterpret(T, data.values)
    PyArray{T}(arr)
end


text:
push!(LOAD_PATH, "src/", "../src/")
using Documenter
using ConstructiveGeometry
makedocs(
	pages = [
		"index.md",
		"primitives.md",
		"transformations.md",
		"operations.md",
		"meshing.md",
		"io.md",
	],
# 	modules = [ConstructiveGeometry],
	sitename="ConstructiveGeometry.jl",
	format=Documenter.HTML(prettyurls=get(ENV,"CI",nothing) !=nothing),
)

deploydocs(
	repo = "github.com/plut/ConstructiveGeometry.jl.git",
)


text:
isjll(name::AbstractString) = endswith(lowercase(strip(name)), lowercase(strip("_jll"))) # TODO check for Artifacts.toml?

function format_compat(v::VersionNumber, drop_patch::Bool=false)
    isempty(v.prerelease) || @warn "Ignoring prerelease suffix"
    isempty(v.build) || @warn "Ignoring build suffix"
    if v.patch == 0 || drop_patch
        if v.minor == 0
            return "$(v.major)"
        else # v.minor is nonzero, v.patch is 0
            return "$(v.major).$(v.minor)"
        end
    else # v.patch is nonzero
        return "$(v.major).$(v.minor).$(v.patch)"
    end
end

function format_compat(compat::AbstractString, drop_patch::Bool)
    compat = String(compat)
    try
        return format_compat(VersionNumber(compat), drop_patch)
    catch
        try
            spec = VersionSpec(semver_spec(compat)) # check to make sure valid
            @assert spec isa VersionSpec
            return String(strip(compat))
        catch
            throw(ArgumentError("not a valid compat entry: $compat"))
        end
    end
end

function format_compat(old, new, drop_patch)
    "$(format_compat(old, false)), $(format_compat(new, drop_patch))"
end

lowerbound(spec::VersionSpec) = minimum(r -> bound2ver(r.lower), spec.ranges)

format_lowerbound_compat(spec::VersionSpec) = ">= $(format_compat(lowerbound(spec)))"

bound2ver(bound::VersionBound) = VersionNumber(bound[1], bound[2], bound[3])

function with_sandbox(fn::Function, ctx::Context)
    project = deepcopy(ctx.env.project)
    manifest = deepcopy(ctx.env.manifest)
    abspath!(ctx, manifest)
    mktempdir() do tmp
        write_project(project, projectfile_path(tmp))
        write_manifest(manifest, manifestfile_path(tmp))
        with_temp_env(tmp) do
            Pkg.resolve()
            return fn(tmp)
        end
    end
end


text:
module Variational

import AdvancedVI
import Bijectors
import DistributionsAD
import DynamicPPL
import StatsBase
import StatsFuns

import Random

# Reexports
using AdvancedVI: vi, ADVI, ELBO, elbo, TruncatedADAGrad, DecayedADAGrad
export
    vi,
    ADVI,
    ELBO,
    elbo,
    TruncatedADAGrad,
    DecayedADAGrad

"""
    make_logjoint(model::Model; weight = 1.0)
Constructs the logjoint as a function of latent variables, i.e. the map z → p(x ∣ z) p(z).
The weight used to scale the likelihood, e.g. when doing stochastic gradient descent one needs to
use `DynamicPPL.MiniBatch` context to run the `Model` with a weight `num_total_obs / batch_size`.
## Notes
- For sake of efficiency, the returned function is closes over an instance of `VarInfo`. This means that you *might* run into some weird behaviour if you call this method sequentially using different types; if that's the case, just generate a new one for each type using `make_logjoint`.
"""
function make_logjoint(model::DynamicPPL.Model; weight = 1.0)
    # setup
    ctx = DynamicPPL.MiniBatchContext(
        DynamicPPL.DefaultContext(),
        weight
    )
    varinfo_init = DynamicPPL.VarInfo(model, ctx)

    function logπ(z)
        varinfo = DynamicPPL.VarInfo(varinfo_init, DynamicPPL.SampleFromUniform(), z)
        model(varinfo)

        return DynamicPPL.getlogp(varinfo)
    end

    return logπ
end

function logjoint(model::DynamicPPL.Model, varinfo, z)
    varinfo = DynamicPPL.VarInfo(varinfo, DynamicPPL.SampleFromUniform(), z)
    model(varinfo)

    return DynamicPPL.getlogp(varinfo)
end

# objectives
function (elbo::ELBO)(
    rng::Random.AbstractRNG,
    alg::AdvancedVI.VariationalInference,
    q,
    model::DynamicPPL.Model,
    num_samples;
    weight = 1.0,
    kwargs...
)
    return elbo(rng, alg, q, make_logjoint(model; weight = weight), num_samples; kwargs...)
end

# VI algorithms
include("advi.jl")

end


text:
module PowerMarketAdeq

include("dummy.jl")

end


text:
@testset "Disjunction2Reification: $(fct_type), $(array_size) items, $(T)" for fct_type in ["vector of variables", "vector affine function"], array_size in [2, 4], T in [Int, Float64]
    base_model = if T == Int
        IntReificationEqualToModel{T}()
    elseif T == Float64
        FloatReificationEqualToModel{T}()
    else
        @assert false
    end
    mock = MOIU.MockOptimizer(base_model)
    model = COIB.Disjunction2Reification{T}(mock)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Reification{MOI.EqualTo{T}},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.Disjunction{NTuple{array_size, MOI.LessThan{T}}},
    )

    if T == Int
        x_array, _ = MOI.add_constrained_variables(model, [MOI.Integer() for _ in 1:array_size])
    elseif T == Float64
        x_array = MOI.add_variables(model, array_size)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables(x_array)
    elseif fct_type == "vector affine function"
        MOIU.vectorize(x_array)
    else
        @assert false
    end
    c = MOI.add_constraint(model, fct,
        CP.Disjunction(tuple(collect(MOI.LessThan(one(T)) for i in 1:array_size)...))
    )

    for i in 1:array_size
        @test MOI.is_valid(model, x_array[i])
    end
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.Reification}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.Reification) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.ZeroOne,)]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.VectorAffineFunction{T}, CP.Reification),
            (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == array_size
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.ZeroOne}()) == array_size
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VectorAffineFunction{T}, CP.Reification}()) == array_size
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == 1

        @test MOI.get(bridge, MOI.ListOfVariableIndices()) == bridge.vars
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.ZeroOne}()) == bridge.vars_bin
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{T}, CP.Reification}()) == bridge.cons_reif
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}()) == [bridge.con_disjunction]
    end

    @testset "Set of variables" begin
        @test length(bridge.vars) == array_size
        @test length(bridge.vars_bin) == array_size
        for i in 1:array_size
            @test MOI.is_valid(model, bridge.vars[i])
            @test MOI.is_valid(model, bridge.vars_bin[i])
            @test MOI.get(model, MOI.ConstraintFunction(), bridge.vars_bin[i]) == bridge.vars[i]
            @test MOI.get(model, MOI.ConstraintSet(), bridge.vars_bin[i]) == MOI.ZeroOne()
        end
    end

    @testset "Reification" begin
        @test length(bridge.cons_reif) == array_size
        for i in 1:array_size
            @test MOI.is_valid(model, bridge.cons_reif[i])
            f = MOI.get(model, MOI.ConstraintFunction(), bridge.cons_reif[i])
            @test length(f.terms) == 2
            @test MOI.get(model, MOI.ConstraintSet(), bridge.cons_reif[i]) == CP.Reification(MOI.LessThan(one(T)))

            t1 = f.terms[1]
            @test t1.output_index == 1
            @test t1.scalar_term.coefficient === one(T)
            @test t1.scalar_term.variable == bridge.vars[i]

            t2 = f.terms[2]
            @test t2.output_index == 2
            @test t2.scalar_term.coefficient === one(T)
            @test t2.scalar_term.variable == x_array[i]
        end
    end

    @testset "Sum" begin
        @test MOI.is_valid(model, bridge.con_disjunction)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_disjunction)
        @test length(f.terms) == array_size
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_disjunction) == MOI.GreaterThan(one(T))

        for i in 1:array_size
            t = f.terms[i]
            @test t.coefficient === one(T)
            @test t.variable == bridge.vars[i]
        end
    end
end


text:
"""
    SearchProblem

A minimalistic, simple interface for search problems
`S` - the state type
`A` - the action type
"""
abstract type SearchProblem{S, A} end
"""
    start_state(p::SearchProblem)

Returns the start state of the problem.
"""
function start_state end
"""
    is_goal_state(p::SearchProblem, s::S)

Returns true if a given state is a goal state
"""
function is_goal_state end
"""
    successors(p::SearchProblem, s::S)

Returns a list of tuples `(sp, a, c)`, i.e. the next state, the action that
yields this state and the cost for this transition.
"""
function successors end
"""
    search_node_type

Returns the SearchNode parametrized according to the SearchProblem
"""
search_node_type(::SearchProblem{S, A}) where {S, A} = SearchNode{S, A}
"""
    state_type

Returns the type of the state the search is run over.
"""
state_type(::SearchProblem{S, A}) where {S, A} = S
"""
    action_type

Returns the type of the action the search is run over.
"""
action_type(::SearchProblem{S, A}) where {S, A} = A


text:
using HSA
include("helpers.jl")
using FactCheck

@with_agents facts("The ISA") do
    rt = Runtime()

    context("Can be retrieved from an agent") do
        agent = get_testagent()

        isa = HSA.agent_info_isa(agent)

        name = HSA.isa_info_name(isa)

        @fact name  --> not("")
    end
end


text:
# This is for rebuilding packages
using Pkg

my_registered = [
    "AbstractLattices"
    "BigCombinatorics"
    "Bijections"
    "ClosedIntervals"
    "Counters"
    "Clines"
    "HyperbolicPlane"
    "LatexPrint"
    "LinearFractionalTransformations"
    "Mods"
    "Multisets"
    "Permutations"
    "RiemannComplexNumbers"
    "ShowSet"
    "SimpleDrawing"
    "SimpleLife"
    "SimplePartitions"
    "SimpleRandom"
    "SimpleTools"
    "SimpleTropical"
    "TwentyFour"
    "Mazes"
    "DrawSimpleGraphs"
    "ChooseOptimizer"
    "SimpleGraphs"
    "FlexLinearAlgebra"
    "SimplePosets"
    "Sudoku"
    "DiscreteFunctions"
    "LatinSquares"
    "SimpleGraphAlgorithms"
    "SimplePosetAlgorithms"
    "SimplePolynomials"
    "LinearAlgebraX"
    "SimpleQuaternions"
    "RingLists"
    "SimplePadics"
    "Spirograph"
    "PlayingCards52"
    "ImplicitGraphs"
]

my_unregistered = [
    "Triangle"
    "SimplePlanarGraphs"
    "SimpleWorld"
    "IntPrint"
    "War"
    "SimpleGF2"
    "SimpleGraphRepresentations"
    "SimplePosetRepresentations"
]

other_packages = [
    "SimplePosetDrawings"
    "SemiIsomorphism"
    "BalancedIncompleteBlockDesigns"
    "HyperbolicDrawSimpleGraphs"
    "FreeCell"
    "Factorions"
    "Diodes"
    "HalfSine"
]

julia_packages = [
    "Cbc"
   # "Clp"
    "Combinatorics"
    "DataStructures"
    "Distributions"
    "FFTW"
    "GLPK"
    "GR"
    "Gurobi"
    "IJulia"
    "IterTools"
    "JuMP"
    "JuliaFormatter"
    "LaTeXStrings"
    "LightXML"
    "LinearAlgebra"
    "MathProgBase"
    "Memoize"
    "Optim"
    "Plots"
    "Polynomials"
    "Primes"
    "ProgressMeter"
    "PyPlot"
    "Random"
    "Nemo"
    "PkgTemplates"
    "Documenter"
    "BenchmarkTools"
]

function load_julia_favorites()
    for pkg in julia_packages
        @info "Adding Julia package: $pkg"
        try
            Pkg.add(pkg)
        catch
            @warn "Failed to add $pkg"
        end
    end
    nothing
end

function my_dev(pkg_name::String)
    pre = "https://github.com/scheinerman/"
    post = ".jl.git"
    full_name = pre * pkg_name * post
    try
        Pkg.develop(PackageSpec(url = full_name))
        @info "Added $pkg_name in development mode"
    catch
        @warn "Unable to add $pkg_name in development mode"
    end
    nothing
end

function my_add(pkg::String)
    my_dev(pkg)
    try
        Pkg.free(pkg)
        @info "$pkg freed"
    catch
        @warn "Unable to add/free $pkg"
    end
    nothing
end

function load_my_unregistered()
    for pkg in sort(my_unregistered)
        my_dev(pkg)
    end
    nothing
end

function load_my_registered()
    for pkg in sort(my_registered)
        my_add(pkg)
    end
    nothing
end

"""
`status_report()` is a version of `Pkg.status()`
organized for my stuff.
"""
function status_report()
    println("Status of my registered packages")
    for pkg in sort(my_registered)
        Pkg.status(pkg)
    end
    println("================================="^2)
    println("Status of my development packages")
    for pkg in sort(my_unregistered)
        Pkg.status(pkg)
    end
    nothing
end

"""
`load_all()`: This is the master command to including into
the `.julia` directory all the packages I like to use.
"""
function load_all()
    load_julia_favorites()
    load_my_registered()
    load_my_unregistered()
    nothing
end


text:
Base.string(x::XFloat16) = string(Float16(reinterpret(Float32, x)))
Base.string(x::XFloat32) = string(Float32(reinterpret(Float64, x)))

Base.show(io::IO, x::XFloat16) = print(io, Float16(reinterpret(Float32, x)))
Base.show(io::IO, x::XFloat32) = print(io, Float32(reinterpret(Float64, x)))

Base.tryparse(::Type{XFloat16}, x::String) = reinterpret(XFloat16, tryparse(Float32, x))
Base.tryparse(::Type{XFloat32}, x::String) = reinterpret(XFloat32, tryparse(Float64, x))

Printf.ini_dec(d::XFloat16, ndigits, digits) =
    Printf.ini_dec(reinterpret(Float32, d), ndigits, digits)

Printf.ini_dec(d::XFloat32, ndigits, digits) =
    Printf.ini_dec(reinterpret(Float64, d), ndigits, digits)

Base.hash(x::XFloat16) = hash(reinterpret(Float32, x))
Base.hash(x::XFloat32) = hash(reinterpret(Float64, x))
Base.hash(x::XFloat16, h::UInt) = hash(Float64(reinterpret(Float32,x)), h)
Base.hash(x::XFloat32, h::UInt) = hash(reinterpret(Float64,x), h)


text:
# anonymous functions:
(x, y) -> x^3 - y + x * y
f = (x, y) -> x^3 - y + x * y
println(f(3, 2) ) # 31

function (x)
    x + 2
end

(x) -> x + 2
x -> x + 2  # lambda syntax
() -> println("hello, Julia")
(x, y, z) -> 3x + 2y - z


text:
## ----------------------------------------------------------------------------
const RED_PIX = RGB{N0f8}(1.0,0.0,0.0)
const GREEN_PIX = RGB{N0f8}(0.0,1.0,0.0)
const BLUE_PIX = RGB{N0f8}(0.0,0.0,1.0)
const BLACK_PIX = RGB{N0f8}(0.0,0.0,0.0)
const WHITE_PIX = RGB{N0f8}(1.0,1.0,1.0)

## ----------------------------------------------------------------------------
function add_margin(a::Matrix, topm::Int, botm::Int, leftm::Int, rightm::Int, tofill)
    M, N = size(a)
    new_a = fill(tofill, M + topm + botm, N + leftm + rightm)
    M, N = size(new_a)
    new_a[(topm + 1):(M - botm), (leftm + 1):(N - rightm)] .= a
    new_a
end
add_margin(a::Matrix, m::Int, tofill) = add_margin(a, m, m, m, m, tofill)

## ----------------------------------------------------------------------------
function _auto_layout(d)
    w = round(Int, sqrt(d), RoundUp)
    h = round(Int, d / w, RoundUp)
    @assert w * h >= d
    (w, h)
end

## ----------------------------------------------------------------------------
function centered(img, d1_nsize, d2_nsize, tofill)
    d1_size, d2_size = size(img)
    any((d1_nsize, d2_nsize) .< (d1_size, d2_size)) &&
        error("size $((d1_size, d2_size)) do not fit in nsize $((d1_nsize, d2_nsize))")

    d1_marg0 = div(d1_nsize - d1_size, 2)
    d1_marg1 = d1_nsize - d1_size - d1_marg0

    d2_marg0 = div(d2_nsize - d2_size, 2)
    d2_marg1 = d2_nsize - d2_size - d2_marg0

    add_margin(img, d1_marg0, d1_marg1, d2_marg0, d2_marg1, tofill)
end
centered(img, gsize, tofill) = centered(img, gsize..., tofill)

## ----------------------------------------------------------------------------
function _compute_size(sizes, layout, margin)
    d1_img_count, d2_img_count = layout
    c = 1
    d1_img_sizes = Dict{Int, Int}()
    d2_img_sizes = Dict{Int, Int}()
    for d1_imgi in 1:d1_img_count
        for d2_imgi in 1:d2_img_count
            c > length(sizes) && break
            d1_size, d2_size = sizes[c] .+ 2 * margin
            d1_old_size = get!(d1_img_sizes, d1_imgi, 0)
            d1_img_sizes[d1_imgi] = max(d1_size, d1_old_size)

            d2_old_size = get!(d2_img_sizes, d2_imgi, 0)
            d2_img_sizes[d2_imgi] = max(d2_size, d2_old_size)
            c += 1
        end
    end
    return d1_img_sizes, d2_img_sizes
end

## ----------------------------------------------------------------------------
function _make_grid(arrs::Vector{Matrix{T}};
        layout = _auto_layout(length(arrs)),
        margin::Int = 5,
        tofill::T = one(T)
    ) where {T}

    d1_img_count, d2_img_count = layout

    img_sizes = size.(arrs)
    d1_img_sizes, d2_img_sizes = _compute_size(img_sizes, layout, margin)

    d1_tot_size, d2_tot_size = sum(values(d1_img_sizes)), sum(values(d2_img_sizes))
    grid::Matrix{T} = fill(tofill, d1_tot_size, d2_tot_size)

    c = 1
    d1_img_idx0 = 1
    for d1_imgi in 1:d1_img_count

        d1_img_size = d1_img_sizes[d1_imgi]
        d1_img_idx1 = d1_img_idx0 + d1_img_size - 1
        d1_range = d1_img_idx0:d1_img_idx1

        d2_img_idx0 = 1
        for d2_imgi in 1:d2_img_count
            c > length(arrs) && break

            d2_img_size = d2_img_sizes[d2_imgi]
            d2_img_idx1 = d2_img_idx0 + d2_img_size - 1
            d2_range = d2_img_idx0:d2_img_idx1

            centered_a = centered(arrs[c], d1_img_size, d2_img_size, tofill)

            grid[d1_range, d2_range] .= centered_a
            c += 1

            d2_img_idx0 = d2_img_idx1 + 1
        end

        d1_img_idx0 = d1_img_idx1 + 1
    end
    grid
end

## ----------------------------------------------------------------------------
function make_grid(imgs::Vector; kwargs...)

    img_untyped_ = map(imgs) do img
        (img isa AbstractPlot) ? plot_to_img(img) : img
    end

    T = eltype(first(img_untyped_))
    for img in img_untyped_
        T = promote_type(T, eltype(img))
    end

    imgs_typed_ = Vector{Matrix{T}}(undef, length(img_untyped_))
    for i in eachindex(img_untyped_)
        imgs_typed_[i] = img_untyped_[i]
    end
    img_untyped_ = nothing

    _make_grid(imgs_typed_; kwargs...)
end


text:
@testset "NfAbs" begin
  include("NfAbs/Conjugates.jl")
  include("NfAbs/NonSimple.jl")
  include("NfAbs/Hilbert.jl")
  include("NfAbs/Subfields.jl")
  include("NfAbs/Elem.jl")
  include("NfAbs/Simplify.jl")
  include("NfAbs/NormRelation.jl")
  include("NfAbs/MPolyGcd.jl")
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractMapper

A mapping strategy for mapping spatial data onto a domain.
"""
abstract type AbstractMapper end

"""
    map(spatialdata, domain, targetvars, mapper)

Map the `targetvars` in `spatialdata` to `domain` with `mapper`.
"""
function map end

# ----------------
# IMPLEMENTATIONS
# ----------------
include("mappers/nearest_mapper.jl")
include("mappers/copy_mapper.jl")