text:
abstract type Sampler end

Base.@kwdef struct NormalSampler <: Sampler
    # (Public) Data
    rng::MersenneTwister
    y0::Vector{Float64}
    X0::Matrix{Float64}
    y1::Vector{Float64}
    X1::Matrix{Float64}
    # (Private) Data
    N0::Int = size(X0, 1)
    N1::Int = size(X1, 1)
    D0::Int = size(X0, 2)
    D1::Int = size(X1, 2)
    # (Public) HyperParameters
    m0β::Vector{Float64} = zeros(D0)
    Σ0β::Matrix{Float64} = 3 * I(D0)
    m0μ::Float64 = 0.0
    c0μ::Float64 = 2.0
    a0τ::Float64 = 1.0
    b0τ::Float64 = 1.0
    # (Private) Parameters
    ϕ::Vector{Float64} = ones(N0) / 2
    τ::Vector{Float64} = ones(N0)
    μ::Vector{Float64} = zeros(N0)
    n::Vector{Int} = [N0; zeros(Int, N0 - 1)]
    r::Vector{Int} = ones(Int, N0)
    d::Vector{Int} = ones(Int, N0)
    f::Vector{Float64} = zeros(N1)
    rmodel::BayesNegativeBinomial.Sampler =
        BayesNegativeBinomial.Sampler(ones(Int, N0), X0)
    β::Vector{Float64} = rmodel.β
    s::Vector{Int} = rmodel.s
    ȳ::Vector{Float64} = zeros(N0)
    v::Vector{Float64} = zeros(N0)
end

Base.@kwdef struct ErlangSampler <: Sampler
    # (Public) Data
    rng::MersenneTwister
    ỹ0::Vector{Float64}
    y0::Vector{Float64} = deepcopy(ỹ0)
    X0::Matrix{Float64}
    y1::Vector{Float64}
    X1::Matrix{Float64}
    event::Vector{Bool} = ones(length(y0))
    # (Private) Data
    N0::Int = size(X0, 1)
    N1::Int = size(X1, 1)
    D0::Int = size(X0, 2)
    D1::Int = size(X1, 2)
    # (Public) HyperParameters
    m0β::Vector{Float64} = zeros(D0)
    Σ0β::Matrix{Float64} = 3 * I(D0)
    a0φ::Float64 = 2.0
    b0φ::Float64 = 0.1
    a0λ::Float64 = 0.1
    b0λ::Float64 = 0.1
    # (Private) Parameters
    ϕ::Vector{Float64} = ones(N0) / 2
    φ::Vector{Float64} = rand(rng, Gamma(a0φ, 1.0 / b0φ), N0)
    λ::Vector{Float64} = rand(rng, Gamma(a0λ, 1.0 / b0λ), 1)
    n::Vector{Int} = [N0; zeros(Int, N0 - 1)]
    r::Vector{Int} = ones(Int, N0)
    d::Vector{Int} = ones(Int, N0)
    f::Vector{Float64} = zeros(N1)
    rmodel::BayesNegativeBinomial.Sampler =
        BayesNegativeBinomial.Sampler(ones(Int, N0), X0)
    β::Vector{Float64} = rmodel.β
    s::Vector{Int} = rmodel.s
    sumy::Vector{Float64} = [0.0]
    sumlogy::Vector{Float64} = zeros(N0)
end

function y(sampler::Sampler)
    return sampler.y0
end

function m(sampler::Sampler)
    return maximum(sampler.r)
end

function f0(sampler::NormalSampler, yi, j)
    kernel = Normal(sampler.μ[j], 1 / √sampler.τ[j])
    return pdf(kernel, yi)
end

function f0(sampler::ErlangSampler, yi, j)
    kernel = Erlang(ceil(Int, sampler.φ[j]), 1.0 / sampler.λ[])
    return pdf(kernel, yi)
end

function update_suffstats!(sampler::NormalSampler)
    (; N0, y0, d, ȳ, v, n) = sampler
    jmax = m(sampler)
    for j in 1:jmax
        ȳ[j] = 0.0
        v[j] = 0.0
        n[j] = 0
    end
    for i in 1:N0
        n[d[i]] += 1
        ȳ[d[i]] += y0[i]
    end
    for j in 1:jmax
        iszero(n[j]) && continue
        ȳ[j] /= n[j]
    end
    for i in 1:N0
        v[d[i]] += (y0[i] - ȳ[d[i]])^2
    end
end

function update_suffstats!(sampler::ErlangSampler)
    (; y0, sumlogy, sumy, d, n, N0) = sampler
    jmax = m(sampler)
    for j in 1:jmax
        n[j] = 0
        sumlogy[j] = 0.0
    end
    sumy[] = 0.0
    for i in 1:N0
        n[d[i]] += 1
        sumy[] += y0[i]
        sumlogy[d[i]] += log(y0[i])
    end
    return nothing
end

function logpφ(sampler::ErlangSampler, φ0, j)
    (; λ, a0φ, b0φ, sumlogy, n) = sampler
    δ0 = ceil(Int, φ0)
    return (
        n[j] * δ0 * log(λ[]) +
        (δ0 - 1) * sumlogy[j] -
        n[j] * logfactorial(δ0 - 1) +
        (a0φ - 1) * log(φ0) -
        b0φ * φ0
    )
end

function update_φ!(rng::AbstractRNG, sampler::ErlangSampler)
    (; φ, a0φ, b0φ, n) = sampler
    pφ = Gamma(a0φ, 1.0 / b0φ)
    for j in 1:m(sampler)
        if n[j] == 0
            φ[j] = rand(rng, pφ)
        else
            φ0 = φ[j]
            d0 = truncated(Normal(φ0, 0.1), 0, Inf)
            φ1 = rand(rng, d0)
            d1 = truncated(Normal(φ1, 0.1), 0, Inf)
            log_ar = 0.0
            log_ar += (logpφ(sampler, φ1, j) + logpdf(d1, φ0))
            log_ar -= (logpφ(sampler, φ0, j) + logpdf(d0, φ1))
            if rand(rng) < exp(log_ar)
                (φ[j] = φ1)
            end
        end
    end
    return nothing
end

function update_λ!(rng::AbstractRNG, sampler::ErlangSampler)
    (; λ, φ, a0λ, b0λ, d, sumy) = sampler
    a1λ = a0λ
    b1λ = b0λ + sumy[]
     for di in d
        a1λ += ceil(φ[di])
    end
    λ[] = rand(rng, Gamma(a1λ, 1.0 / b1λ))
    return nothing
end

function update_χ!(rng::AbstractRNG, sampler::ErlangSampler)
    update_suffstats!(sampler)
    update_φ!(rng, sampler)
    update_λ!(rng, sampler)
end

function update_χ!(rng::AbstractRNG, sampler::NormalSampler)
    (; μ, τ, n, ȳ, v, m0μ, c0μ, a0τ, b0τ) = sampler
    update_suffstats!(sampler)
    for j in 1:m(sampler)
        mμ1 = (c0μ * n[j] * ȳ[j] + m0μ) / (c0μ * n[j] + 1)
        cμ1 = c0μ / (c0μ * n[j] + 1)
        aτ1 = a0τ + n[j] / 2
        bτ1 = b0τ + n[j] * (ȳ[j] - m0μ)^2 / 2 / (c0μ * n[j] + 1) + v[j] / 2
        τ[j] = rand(rng, Gamma(aτ1, 1 / bτ1))
        μ[j] = mμ1 + randn(rng) * √(cμ1 / τ[j])
    end
    update_y!(rng, mdl)
end

function update_d!(rng::AbstractRNG, sampler::Sampler)
    (; N0, y0, d, r) = sampler
    p = zeros(m(sampler))
    for i in 1:N0
        resize!(p, r[i])
        for j in 1:r[i]
            p[j] = f0(sampler, y0[i], j)
        end
        p ./= sum(p)
        d[i] = rand(rng, Categorical(p))
    end
    return nothing
end

function update_n!(sampler::Sampler)
    (; N0, d, n) = sampler
    jmax = m(sampler)
    for j in 1:jmax
        n[j] = 0
    end
    for i in 1:N0
        n[d[i]] += 1
    end
    return nothing
end

function update_ϕ!(sampler::Sampler)
    (; X0, β, ϕ) = sampler
    mul!(ϕ, X0, β)
    @. ϕ = 1.0 / (1.0 + exp(ϕ))
end

function update_r!(rng::AbstractRNG, sampler::Sampler)
    (; N0, d, r, ϕ, s) = sampler
    for i in 1:N0
        zi = rand(rng, Beta(r[i] - d[i] + 1, d[i]))
        vi = rand(rng, Gamma(r[i] + s[] - 1, 1))
        r[i] = rand(rng, Poisson((1 - ϕ[i]) * zi * vi)) + d[i]
    end
end

function get_w(θ0, s0, j)
    exp(
        - log(j) +
        logabsbinomial(j + s0 - 2, j - 1)[1] +
        s0 * log(θ0) +
        (j - 1) * log(1 - θ0) +
        log(_₂F₁(BigFloat(j + s0 - 1), 1, j + 1, 1 - θ0))
    )
end

function update_β!(rng::AbstractRNG, sampler::Sampler)
    (; r, rmodel) = sampler
    @. rmodel.y = r - 1
    BayesNegativeBinomial.step!(rng, rmodel)
end

function update_f!(sampler::Sampler)
    (; N1, y1, X1, β, f, s) = sampler
    for i in 1:N1
        f[i] = 0.0
        ϕ1 = 1 / (1 + exp(X1[i, :] ⋅ β))
        for j in 1:m(sampler)
            wj = get_w(ϕ1, s[], j)
            f[i] += wj * f0(sampler, y1[i], j)
        end
    end
end

function update_y!(rng::AbstractRNG, mdl::ErlangSampler)
    (; N0, y0, ỹ0, event) = mdl
    for i in 1:N0
        if !event[i]
            pdf = Erlang(ceil(Int, sampler.φ[d[i]]), 1.0 / sampler.λ[])
            tpdf = Truncated(pdf, ỹ0[i], Inf)
            y0[i] = rand(rng, tpdf)
        end
    end
end

function step!(rng::AbstractRNG, mdl::Sampler)
    update_χ!(rng, mdl)
    update_d!(rng, mdl)
    # update_n!(mdl)
    update_f!(mdl)
    update_ϕ!(mdl)
    update_r!(rng, mdl)
    update_β!(rng, mdl)
end

function sample(rng::AbstractRNG, smpl::Sampler; mcmcsize = 4000, burnin = 2000)
    (; N1, D1, f, β) = smpl
    chainf = [zeros(N1) for _ in 1:(mcmcsize - burnin)]
    chainβ = [zeros(D1) for _ in 1:(mcmcsize - burnin)]
    for iter in 1:mcmcsize
        step!(rng, smpl)
        if iter > burnin
            chainf[iter - burnin] .= f
            chainβ[iter - burnin] .= β
        end
    end
    return chainf, chainβ
end


text:
function fizzbuzz()
    for x in 1 : 100
        if x % 15 == 0
            print("FizzBuzz")
        elseif x % 3 == 0
            print("Fizz")
        elseif x % 5 == 0
            print("Buzz")
        else
            print(x)
        end
        print(" ")
    end
end

fizzbuzz()

ext: jl


text:
import Pkg; Pkg.add("Documenter")

using Documenter, PhaseSpaceTools

makedocs(sitename="PhaseSpaceTools.jl",
pages = [
"Overview" => "index.md",
"Introduction to Phase Space" => "phasespace.md",
"Quantum States" => "states.md",
"Distributions" => "distributions.md",
"Solving SDES" => "sdes.md"
],
format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true")
)

deploydocs(
    repo = "github.com/AshtonSBradley/PhaseSpaceTools.jl.git",
)


text:
module TestConstructors
    using Base.Test
    using CategoricalArrays
    using CategoricalArrays: DefaultRefType
    using Compat

    typealias String Compat.ASCIIString

    pool = CategoricalPool{String}()

    @test isa(pool, CategoricalPool{String})

    @test isa(pool.index, Vector{String})
    @test length(pool.index) == 0

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 0

    pool = CategoricalPool{Int, UInt8}()

    @test isa(pool, CategoricalPool{Int, UInt8})

    @test isa(pool.index, Vector{Int})
    @test length(pool.index) == 0

    @test isa(pool.invindex, Dict{Int, UInt8})
    @test length(pool.invindex) == 0

    pool = CategoricalPool(["a", "b", "c"])

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector{String})
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)

    pool = CategoricalPool{String, UInt8}(["a", "b", "c"])

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector{String})
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict{String, UInt8})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === UInt8(1)
    @test pool.invindex["b"] === UInt8(2)
    @test pool.invindex["c"] === UInt8(3)

    pool = CategoricalPool(
        Dict(
            "a" => DefaultRefType(1),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(3),
        )
    )

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector{String})
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)

    # TODO: Make sure that invindex input is exhaustive
    # Raise an error if map misses any entries
    pool = CategoricalPool(
        Dict(
            "a" => 1,
            "b" => 2,
            "c" => 3,
        )
    )

    @test isa(pool, CategoricalPool)

    @test isa(pool.index, Vector{String})
    @test length(pool.index) == 3
    @test pool.index[1] == "a"
    @test pool.index[2] == "b"
    @test pool.index[3] == "c"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["a"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["c"] === DefaultRefType(3)

    pool = CategoricalPool(["c", "b", "a"])

    @test isa(pool, CategoricalPool)

    @test length(pool.index) == 3
    @test pool.index[1] == "c"
    @test pool.index[2] == "b"
    @test pool.index[3] == "a"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["c"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["a"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(1)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(3)

    pool = CategoricalPool(
        Dict(
            "a" => DefaultRefType(3),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(1),
        )
    )

    @test isa(pool, CategoricalPool)

    @test length(pool.index) == 3
    @test pool.index[1] == "c"
    @test pool.index[2] == "b"
    @test pool.index[3] == "a"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["c"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["a"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(1)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(3)

    pool = CategoricalPool(["c", "b", "a"], ["c", "b", "a"])

    @test isa(pool, CategoricalPool)

    @test length(pool.index) == 3
    @test pool.index[1] == "c"
    @test pool.index[2] == "b"
    @test pool.index[3] == "a"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["c"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["a"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(1)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(3)

    pool = CategoricalPool(
        Dict(
            "a" => DefaultRefType(3),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(1),
        ),
        ["c", "b", "a"]
    )

    @test isa(pool, CategoricalPool)

    @test length(pool.index) == 3
    @test pool.index[1] == "c"
    @test pool.index[2] == "b"
    @test pool.index[3] == "a"

    @test isa(pool.invindex, Dict{String, DefaultRefType})
    @test length(pool.invindex) == 3
    @test pool.invindex["c"] === DefaultRefType(1)
    @test pool.invindex["b"] === DefaultRefType(2)
    @test pool.invindex["a"] === DefaultRefType(3)

    @test isa(pool.order, Vector{DefaultRefType})
    @test length(pool.order) == 3
    @test pool.order[1] === DefaultRefType(1)
    @test pool.order[2] === DefaultRefType(2)
    @test pool.order[3] === DefaultRefType(3)
end


text:
"""
    SerendipityRefFE(::Type{T},p::Polytope,order::Int) where T
    SerendipityRefFE(::Type{T},p::Polytope,orders::Tuple) where T

Returns an instance of `LagrangianRefFE`, whose underlying approximation space
is the serendipity space of order `order`. Implemented for order from 1 to 4.
The type of the polytope `p` has to implement all the queries detailed in the
constructor [`LagrangianRefFE(::Type{T},p::Polytope{D},orders) where {T,D}`](@ref).

# Examples

```jldoctest
using Gridap.ReferenceFEs

order = 2
reffe = SerendipityRefFE(Float64,QUAD,order)

println( num_dofs(reffe) )

# output
8

```
"""
function SerendipityRefFE(::Type{T},p::Polytope,order::Int) where T
  @assert is_n_cube(p) "Polytope not compatible with serendipity elements"
  if order > 0
    sp = SerendipityPolytope(p)
  else
    sp = p
  end
  LagrangianRefFE(T,sp,order)
end

function SerendipityRefFE(::Type{T},p::Polytope,orders::Tuple) where T
  order = first(orders)
  @assert all( orders .== order ) "Anisotropic serentopity FEs not allowed"
  SerendipityRefFE(T,p,order)
end

# Helper private type
struct SerendipityPolytope{D,P} <: Polytope{D}
  hex::P
  SerendipityPolytope(p::Polytope{D}) where D = new{D,typeof(p)}(p)
end

# Implemented polytope interface

function get_faces(p::SerendipityPolytope)
  get_faces(p.hex)
end

function get_dimranges(p::SerendipityPolytope)
  get_dimranges(p.hex)
end

function Polytope{N}(p::SerendipityPolytope,Nfaceid::Integer) where N
  face_hex = Polytope{N}(p.hex, Nfaceid)
  SerendipityPolytope(face_hex)
end

function Polytope{D}(p::SerendipityPolytope{D},Dfaceid::Integer) where D
  @assert Dfaceid == 1 "Only one D-face"
  p
end

function get_vertex_coordinates(p::SerendipityPolytope)
  get_vertex_coordinates(p.hex)
end

function (==)(a::SerendipityPolytope{D},b::SerendipityPolytope{D}) where D
  a.hex == b.hex
end

function get_vertex_permutations(p::SerendipityPolytope)
  get_vertex_permutations(p.hex)
end

LagrangianRefFE(p::SerendipityPolytope) = LagrangianRefFE(p.hex)

is_simplex(p::SerendipityPolytope) = false

is_n_cube(p::SerendipityPolytope) = true

get_extrusion(p::SerendipityPolytope{D}) where D = Point(tfill(HEX_AXIS,Val{D}()))

# Implemented polytope interface for LagrangianRefFEs

function _s_filter(e,order)
  sum( [ i for i in e if i>1 ] ) <= order
end

function compute_monomial_basis(::Type{T},p::SerendipityPolytope{D},orders) where {T,D}
  MonomialBasis{D}(T,orders,_s_filter)
end

function compute_own_nodes(p::SerendipityPolytope{0},orders)
  compute_own_nodes(p.hex,orders)
end

function compute_own_nodes(p::SerendipityPolytope{1},orders)
  compute_own_nodes(p.hex,orders)
end

function compute_own_nodes(p::SerendipityPolytope{2},orders)
  order, = orders
  if order == 4
    o = (2,2)
  elseif order in (0,1,2,3)
    o=(1,1)
  else
    @unreachable "Serendipity elements only up to order 4"
  end
  compute_own_nodes(p.hex,o)
end

function compute_own_nodes(p::SerendipityPolytope{3},orders)
  Point{3,Float64}[]
end

function compute_own_nodes(p::SerendipityPolytope,orders)
  @unreachable "Serendipity elements only up to 3d"
end

function compute_face_orders(
  p::SerendipityPolytope{D},face::SerendipityPolytope{N},iface::Int,orders) where {D,N}
  compute_face_orders(p.hex,face.hex,iface,orders)
end


text:
"""
    `build_blp_model(bp::BilevelLP, solver; [comp_method])`

Build a JuMP model (constraints and objective) based
on the data from bp and with a solver.

Returns a tuple `(m, x, y, λ, s)`.

An optional keyword `comp_method` can be passed for
handling complementarity constraints, default is `SOS1Complementarity`.
"""
function build_blp_model(bp::BilevelLP, optimizer; comp_method = SOS1Complementarity())
    m = JuMP.Model(optimizer)
    @variable(m, bp.xl[j] <= x[j=1:bp.nu] <= bp.xu[j])
    @variable(m, y[j=1:bp.nl])
    @constraint(m, uppercons[i=1:bp.mu],
        sum(bp.G[i,j]*x[j] for j in 1:bp.nu) +
        sum(bp.H[i,j]*y[j] for j in 1:bp.nl) <= bp.q[i])
    @objective(m, Min, sum(bp.cx .* x) + dot(bp.cy, y))
    for j in bp.Jx
        JuMP.set_integer(x[j])
    end
    # adding SOS1 constraints
    return build_blp_model(m, bp, x, y, comp_method = comp_method)
end

"""
    `build_blp_model(m::JuMP.Model, bp::BilevelLP, x, y; [comp_method])`

Adds the lower-level constraints and optimality conditions to
an existing JuMP model. This assumes the upper-level feasibility
constraints and objective have already been set.

Returns a tuple `(m, x, y, λ, s)`.

An optional keyword `comp_method` can be passed for
handling complementarity constraints, default is `SOS1Complementarity`.
"""
function build_blp_model(m::JuMP.Model, bp::BilevelLP, x, y; comp_method = SOS1Complementarity())
    @variable(m, s[1:bp.ml] >= 0) # lower-level slack variables
    @constraint(m, lowercons[i=1:bp.ml],
        sum(bp.A[i,j]*x[j] for j in 1:bp.nu) +
        sum(bp.B[i,j]*y[j] for j in 1:bp.nl) + s[i] == bp.b[i]
    )
    @variable(m, λ[1:bp.ml] >= 0)
    @variable(m, σ[1:bp.nl] >= 0) # dual of lower-level lower bound
    for j in Base.OneTo(bp.nl)
        if !bp.yl[j]
            # sigma at 0 if free variable
            @constraint(m, σ[j] == 0)
        else
            JuMP.set_lower_bound(y[j], 0)
        end
    end
    @constraint(m, [j in 1:bp.nl], bp.d[j] + dot(bp.F[:,j], x) + dot(bp.B[:,j], λ)  - σ[j] == 0)
    add_complementarity_constraint(m, comp_method, s, λ, y, σ)
    return (m, x, y, λ, s)
end

"""
    `build_blp_model(m, B, d, x, y, s, F; [comp_method])`

Build the bilevel JuMP model from the data without
grouping everything into a `BilevelLP`.
An optional keyword `comp_method` can be passed for
handling complementarity constraints, default is `SOS1Complementarity`.
"""
function build_blp_model(m::JuMP.Model, B::M, d, x, y, s, F; comp_method = SOS1Complementarity()) where {M<:MT}
    @variable(m, λ[1:bp.ml] >= 0)
    @variable(m, σ[1:bp.nl] >= 0.) # dual of lower-level lower bound
    for j in Base.OneTo(bp.nl)
        if !bp.yl[j]
            # sigma at 0 if free variable
            @constraint(m, σ[j] == 0.)
        else
            JuMP.set_lower_bound(y[j], 0.)
        end
    end
    @constraint(m, d .+ F' * x .+ B' * λ .- σ .== 0)
    add_complementarity_constraint(m, comp_method, s, λ, y, σ)
    return (m, λ, s)
end

function build_blp_model(m::JuMP.Model, B::M, d, s; comp_method = SOS1Complementarity()) where {M<:MT}
    (ml, nl) = size(B)
    @variable(m, λ[1:ml] >= 0)
    @constraint(m, kkt[j=1:nl], d[j] + dot(B[:,j], λ) == 0)
    add_complementarity_constraint(m, comp_method, s, λ, [], [])
    return (m, λ, s)
end


text:
using Flux
using Random

mutable struct GVF_DQNAgent{M, TN, O, LU, AP<:AbstractValuePolicy, Φ, ERP} <: AbstractAgent
    model::M
    target_network::TN
    opt::O
    lu::LU
    ap::AP
    er::ExperienceReplay{ERP}
    γ::Float32
    batch_size::Int64
    tn_counter_init::Int64
    target_network_counter::Int64
    action::Int64
    prev_s::Φ
end

GVF_DQNAgent(model, target_network, opt, lu, ap, size_buffer, γ, batch_size, tn_counter_init, s) =
    GVF_DQNAgent(model,
                 target_network,
                 opt,
                 lu,
                 ap,
                 ExperienceReplay(100000,
                                  (Array{Float32, 1}, Int64, Array{Float32, 1}, Float32, Bool),
                                  (:s, :a, :sp, :r, :t)),
                 γ,
                 batch_size,
                 tn_counter_init,
                 tn_counter_init,
                 0,
                 s)

function RLCore.start!(agent::DQNAgent, env_s_tp1, rng::AbstractRNG; kwargs...)
    # Start an Episode
    agent.action = sample(agent.ap,
                          agent.model(env_s_tp1),
                          rng)
    return agent.action
end

function RLCore.step!(agent::DQNAgent, env_s_tp1, r, terminal, rng::AbstractRNG; kwargs...)

    add!(agent.er, (Float32.(agent.prev_s), agent.action, copy(Float32.(env_s_tp1)), r, terminal))

    if size(agent.er)[1] > 1000
        e = sample(agent.er, agent.batch_size; rng=rng)
        update_params!(agent, e)
    end

    agent.prev_s .= env_s_tp1
    agent.action = sample(agent.ap,
                          agent.model(agent.prev_s),
                          rng)

    return agent.action
end

function update_params!(agent::DQNAgent, e)

    if agent.tn_counter_init > 0
        update!(agent.model, agent.lu, agent.opt, e.s, e.a, e.sp, e.r, e.t, agent.target_network)

        if agent.target_network_counter == 1
            agent.target_network_counter = agent.tn_counter_init
            agent.target_network = mapleaves(
                Flux.Tracker.data,
                deepcopy(agent.model))
        else
            agent.target_network_counter -= 1
        end
    else
        update!(agent.model, agent.lu, agent.opt, e.s, e.a, e.sp, e.r, e.t)
    end
    return nothing

end


text:
# Learn basic Julia stuff
using Test

module MyJulia
using Test
# Learn operator overloading
mutable struct Node
    id::Int64

    Node() = ( # incomplete initialization
        me = new();
        me
    )
    Node(id) = (
        me = new();
        me.id = id;
        me
    )
end
# https://discourse.julialang.org/t/proper-way-to-overload-operators/19872

# One can use function name or symbol to implement operator overloading

Base.:(<)(x::Node, y::Node) = x.id < y.id # Base.isless(x::Node, y::Node) = x.id < y.id
function Base.:(==)(x::Node,y::Node)
    return x.id == y.id
end
# used by sort()
Base.isless(x::Node, y::Node) = error("Unsupported operation")

# isequal => ===
#Base.isequal(x::Node, y::Node) = x.id == y.id

function learn_Node()
    @test Node(7) == Node(7)
    @test Node(7) != Node(13)

    @test Node(7) < Node(13)
    @test Node(7) <= Node(13)
    @test Node(13) <= Node(13)

    @test Node(13) > Node(7)
    @test Node(13) >= Node(7)
    @test Node(13) >= Node(13)
end

end # module MyJulia

mj = MyJulia

#mj.learn_Node()

function learn_resize!()
    x = [1,2,3]
    y = x
    @test y === x
    @test 3 == length(y)

    resize!(x, 1)
    @test y === x
    @test 1 == length(y)
    @test 1 == y[1]

    resize!(x, 0)
    @test y === x
    @test isempty(y)

end

function learn_searchsortedfirst()
    v = Float64[1, 2, 3, 4, 5]

    i = searchsortedfirst(v, 1)
    @test 1 == i

    i = searchsortedfirst(v, 1.5)
    @test 2 == i

    i = searchsortedfirst(v, 2)
    @test 2 == i
end

"""

Learn for-loop syntax
```
for i in 1:2, j in 3:4
end
```
"""
function learn_for_comma()
    ijs = []

    for i in 1:2, j in 3:4
        push!(ijs, (i,j))
    end

    @test [
        (1,3),
        (1,4),
        (2,3),
        (2,4)
    ] == ijs
end

"""
    learn_comprehension()

There are two ways to do two for-loops using comprehension.
"""
function learn_comprehension()
    # List of pairs
    x1 = [(i,j) for i in 1:2 for j in 3:4]
    @test [
        (1, 3),
        (1, 4),
        (2, 3),
        (2, 4)
    ] == x1

    # Matrix of pairs
    x2 = [(i,j) for i in 1:2,j in 3:4]
    @test [
        (1, 3) (1, 4);
        (2, 3) (2, 4)
    ] == x2

    # Flatten
    # Notice the order is different from x1
    @test [(1, 3), (2, 3), (1, 4), (2, 4)] == x2[:]
    @test [
        (1, 3),
        (2, 3),
        (1, 4),
        (2, 4)
    ] == vec(x2)
end
@testset "Node" begin
    MyJulia.learn_Node()
end
@testset "Basic" begin
    learn_resize!()
    learn_searchsortedfirst()
end
@testset "For-Loop" begin
    learn_for_comma()
    learn_comprehension()
end

nothing


text:
using ImageFiltering: Fill, KernelFactors.gaussian, centered, kernelfactors, imfilter!, padarray, Pad
using LinearAlgebra: qr
using PyPlot: Figure
using PyPlot, BenchmarkTools
using Images: assess_psnr

#TODO: add see section about timing in the docs to the api
"""
    pflap(args; kwargs)

Find a transformation flow (complex displacment field), that transforms image `source` to image `target`.
Returns the transformation a [`Flow`](@ref) and the registered `source` image.

# Arguments
- `target::Image`: the image we want `source` to look like.
- `source::Image`: warped image we want to transform into `target`.

# Keyword Arguments
- `filter_count::Integer=3`: the number of basis filters used in `single_lap` calls (so far only =3 implemented).
- `max_repeats::Integer=1`: the maximum number of times an iteration of one filter size can be repeated.
- `display::Bool=true`: use verbose prints and return an array of figures.
- `timer::TimerOutput=TimerOutput("pflap")`: provide a timer which times certain blocks in the function.
- `match_source_histogram::Bool=true`: choose to match the `source` image histogram to the `target` image histogram.
- `rescale_intensities::Bool=false`: choose to rescale! intensities of the images to the whole `[0, 1]` spectrum.

# Outputs
- `flow::Flow`: is the complex vector field that transforms `source` closer to `target`.
- `source_reg::Image`: is the image `source` transformed by `flow`.
- [`figs::Matrix{Figure}`: is a 2D array of `PyPlot` Figures which shows the work of the algorithm at each iteration. For each iteration there are 3 Figures in this order: 1) current `u_est`, 2) newest addition to `u_est` `Δ_u`, 3) current `source_reg`.]

# Describtion
Implements the basic concept of `Algorithm 2` from the [paper](http://www.ee.cuhk.edu.hk/~tblu/monsite/pdfs/gilliam1701.pdf) without some features.
It uses `single_lap` iteratively; in each iteration using the transformation estimated by `single_lap` to warp the `source` image closer to
the `target` image and then using this warped closer image as the source image in the next iteration, while using progressively smaller
`filter_half_sizes` to estimate even small and faster varying displacements.

## Note:
If `display=true` is chosen the function will also return an array of figures demonstrating the
work of the registration.

See also: [`single_lap`](@ref), [`imgshow`](@ref),[`imgshowflow`](@ref), [`warp_imgshowflow`](@ref), [`Flow`](@ref)
"""
function pflap(target::Image,
               source::Image;
               filter_count::Integer=3,
               max_repeats::Integer=3,
               display::Bool=true,
               timer::TimerOutput=TimerOutput("pflap"),
               prefilter::Bool=false,
               match_source_histogram::Bool=true,
               rescale_intensities::Bool=false,
               warp_border_strat::Symbol=:fill_img)

    @timeit_debug timer "setup" begin
        # convert images to floats
        source = Float64.(source)
        target = Float64.(target)

        # rescale images to have the whole [0, 1] spectrum.
        if rescale_intensities
            target, source = rescale!(target, source)
        end

        # match histogram
        if match_source_histogram
            @timeit_debug timer "hist match" begin
                source = adjust_histogram(source, Matching(targetimg = target)) # takes 8ms on average
            end
        end

        # pad with zeros if sizes difer.
        target, source = pad_images(target, source)

        image_size = size(target)

        # set number of layers in the filter pyramid.
        level_count = floor(Int64, log2(minimum(size(target))/8)+1)+1

        @assert ((2^(level_count)+1) <= minimum(image_size)) "level number results in a filter larger than the size of the input images."

        # displacement init.
        u_est = zeros(ComplexF64, size(target))

        # filter half sizes array eg. [16, 8, 4, 2, 1]
        half_size_pyramid = Int64.(2 .^ range(level_count-1, stop=0, length=level_count))

        # threshold above which we repeat the iteration.
        psnr_threshold=0.3

        if display
            num_plots = 4
            figs = Array{Figure}(undef, level_count, max_repeats, num_plots)
        end

        psnr_before = 0

        source_reg = source
    end

    for level in 1:level_count

        @timeit_debug timer "single filter pyramid level" begin
            if display
                println("###################")
                println("ITERATION: ", level)
                println("filter_half_size: ", half_size_pyramid[level])
            end

            filter_half_size = half_size_pyramid[level]::Int
            window_size = Int64.(2 .* filter_half_size .* (1, 1) .+ 1)
            window_half_size = Int64.((window_size .- 1) ./ 2)

            if prefilter
                @timeit_debug timer "prefiltering" begin
                    target_inner = highpass_image(target, filter_half_size)
                end
            else
                target_inner = target
            end

            for iter_repeat in 1:max_repeats

                if display
                    println("inner loop iter: $iter_repeat")
                end

                if prefilter
                    @timeit_debug timer "prefiltering" begin
                        source_reg_inner = highpass_image(source_reg, filter_half_size)
                    end
                else
                    source_reg_inner = source_reg
                end

                @timeit_debug timer "lap" begin
                    Δ_u = single_lap(target_inner, source_reg_inner, filter_half_size, window_size, filter_count, timer=timer)
                end

                # USE INPAINTING TO CORRECT U_EST:
                # 1) replicate borders
                @timeit_debug timer "inpainting" begin
                    @timeit_debug timer "replicating borders" begin
                        middle_vals = Δ_u[window_half_size[1]+1:end-window_half_size[1],
                                          window_half_size[2]+1:end-window_half_size[2]]
                        Δ_u = parent(padarray(middle_vals, Pad(:replicate, window_half_size...)))
                    end # "replication borders"
                        # 2) inpaint middle missing values
                    if all(isnan.(real(Δ_u))) # if all are NaNs
                        Δ_u = zeros(ComplexF64, image_size)
                    elseif any(isnan.(Δ_u)) # if some are NaNs
                        if display
                            println("NaN count: $(count(isnan, Δ_u))")
                        end
                        inpaint_nans!(Δ_u)
                    end
                end # "inpainting"

                # SMOOTH U_EST WITH A GAUSSIAN FILTER:
                @timeit_debug timer "smoothing" begin
                    Δ_u = smooth_with_gaussian!(Δ_u, window_half_size)
                end

                # add the Δ_u to my u_est
                u_est_adept = u_est + Δ_u

                @timeit_debug timer "image interpolation" begin
                    # linear interpolation
                    if warp_border_strat != :fill_img
                        source_reg = warp_img(source, -real(u_est_adept), -imag(u_est_adept), border_strat=warp_border_strat)
                    else
                        source_reg = warp_img(source, -real(u_est_adept), -imag(u_est_adept), target)
                    end
                end # "interpolation"

                if prefilter
                    @timeit_debug timer "prefiltering" begin
                        source_reg_inner = highpass_image(source_reg, filter_half_size)
                    end
                else
                    source_reg_inner = source_reg
                end

                # calculate the psnr after warping the source by u_est_adept
                psnr_after = assess_psnr(source_reg_inner, target_inner)

                # compare psnr before and after to see if there was an improvement larger than psnr_threshold
                if psnr_after - psnr_before > psnr_threshold
                    # nice, the improvement is large -> try to run again with the same filter size.
                    u_est = u_est_adept
                    psnr_before = psnr_after
                    if display
                        println("\t status: success, repeating filter size")
                    end
                else
                    # the improvement isn't larger than psnr_threshold -> next filter size
                    # check to see if it was even an improvement:
                    if psnr_after - psnr_before > 0
                        # it was an improvement, save the variables:
                        u_est = u_est_adept
                        psnr_before = psnr_after
                        if display
                            println("\t status: success, next filter size")
                            add_figs_pflap(figs, level, iter_repeat, u_est, Δ_u, source_reg, level_count, max_repeats, filter_half_size)
                        end
                    else
                        # it wasn't an improvement, skip
                        if display
                            println("\t status: failure, trashing")
                        end
                    end
                    break # go for next filter size.
                end

                if display
                    add_figs_pflap(figs, level, iter_repeat, u_est, Δ_u, source_reg, level_count, max_repeats, filter_half_size)
                end
            end
        end # "single filter pyramid level"
    end

    if display; return u_est, source_reg, figs; end
    return u_est, source_reg
end

"""
    add_figs_pflap(figs, level, iter_repeat, u_est, Δ_u, source_reg, level_count, max_repeats)

Fill the figs array specifically for the `pflap` method.
"""
function add_figs_pflap(figs, level, iter_repeat, u_est, Δ_u, source_reg, level_count, max_repeats, fhs)
    title_string = "(Level: $(string(level))/$(string(level_count))), (Iter: $(string(iter_repeat))/$(string(max_repeats))), (fhs: $fhs)"
    figs[level, iter_repeat, 1] = showflow(u_est, figtitle="U_EST " * title_string)
    figs[level, iter_repeat, 2] = showflow(Δ_u, figtitle="Δ_U " * title_string)
    figs[level, iter_repeat, 3] = imgshow(source_reg, figtitle="SOURCE_REG " * title_string)
end

"""
    sparse_pflap(args; kwargs)

Find a transformation flow (complex displacment field), that transforms image `source` to closer to image `target`.
Returns the transformation a [`Flow`](@ref) and the registered `source` image.

# Arguments:
- `target::Image`: target/fixed grayscale image.
- `source::Image`: source/moving grayscale image.

# Keyword Arguments:
- `filter_count::Integer=3`: the number of basis filters used (so far only =3 implemented).
- `timer::TimerOutput=TimerOutput("sparse pflap")`: provide a timer which times certain blocks in the function.
- `display::Bool=false`: verbose and debug prints.
- `point_count::Int=500`: the number of points attempted to be found at the edges of the target image.
- `spacing::Int=10`: the minimal distance between two points.
- `match_source_histogram::Bool=true`: choose to match the `source` image histogram to the `target` image histogram.

## Note:
If `display=true` is chosen the function will also return an array of figures demonstrating the
work of the registration.

See also: [`sparse_lap`](@ref), [`pflap`](@ref), [`showflow`](@ref), [`single_lap`](@ref), [`sparse_pflap`](@ref).
"""
function sparse_pflap(target::Image,
                      source::Image;
                      filter_count::Integer=3,
                      max_repeats::Integer=1,
                      display::Bool=false,
                      point_count::Int=500,
                      spacing::Int=13,
                      timer::TimerOutput=TimerOutput("sparse pflap"),
                      match_source_histogram::Bool=true,
                      rescale_intensities::Bool=false)

    @timeit_debug timer "setup" begin
        # convert images to floats
        source = Float64.(source)
        target = Float64.(target)

        # rescale images to have the whole [0, 1] spectrum.
        if rescale_intensities
            target, source = rescale!(target, source)
        end

        # adjust the source histogram to be like target
        if match_source_histogram
            @timeit_debug timer "hist match" begin
                source = adjust_histogram(source, Matching(targetimg = target)) # takes 8ms on average
            end
        end

        # pad with zeros if sizes difer.
        target, source, mask = pad_images_mask(target, source)

        image_size = size(target)

        # set number of layers in the filter pyramid.
        level_count = floor(Int64, log2(minimum(size(target))/8)+1)+1

        @assert ((2^(level_count)+1) <= minimum(image_size)) "level number results in a filter larger than the size of the input images."

        # displacement init.
        u_est = zeros(Complex{Float64}, image_size...)

        # filter half sizes array eg. [16, 8, 4, 2, 1]
        half_size_pyramid = Int64.(2 .^ range(level_count-1, stop=0, length=level_count))

        if display
            num_plots = 4
            figs = Array{Figure}(undef, level_count, max_repeats, num_plots)
        end

        source_reg = source

        #TODO edit
        # get edge points
        # mask = falses(size(target))
        # fhs = 3
        # mask[fhs+1:end-fhs, fhs+1:end-fhs] .= true
        @timeit_debug timer "find edge points" begin
            inds = find_edge_points(target, spacing=spacing, point_count=point_count, mask=mask)
        end
        if display
            println("ind count: ", length(inds))
        end

    end # "setup"

    for level in 1:level_count

        @timeit_debug timer "single filter pyramid level" begin
            if display
                println("###################")
                println("ITERATION: ", level)
                println("filter_half_size: ", half_size_pyramid[level])
            end

            fhs = half_size_pyramid[level]
            window_size = Int64.(2 .* fhs .+ (1, 1))

            for iter_repeat in 1:max_repeats

                @timeit_debug timer "filter inds" begin
                    current_inds = filter_border_inds(inds, size(target), fhs)
                end
                if display
                    println("inner loop iter: $iter_repeat")
                    println("\tcurrent inds: ", length(current_inds))
                end
                @timeit_debug timer "single lap at points" begin
                    new_estim_at_inds, current_inds = single_lap_at_points(target, source_reg, fhs, window_size, current_inds, timer=timer, display=display)
                end

                # interpolate flow
                if all(isnan, new_estim_at_inds)
                    @timeit_debug timer "flow interpolation" begin
                        if display
                            println("\tIMPORTANT: all new estim vectors are NaN.")
                        end
                    end
                else
                    @timeit_debug timer "flow interpolation" begin
                        # extract u_est_vecs from the previous estimate and add them to the newly estimated ones
                        # then interpolate a new flow.
                        if display
                            println("\tnon NaN new estim flow vector count: ", length(new_estim_at_inds))
                        end
                        u_est_at_inds = map(ind -> u_est[ind], current_inds) .+ new_estim_at_inds
                        @assert any(.!isnan.(u_est_at_inds)) "new estim nans", count(isnan, u_est_at_inds), length(u_est_at_inds)
                        u_est = interpolate_flow(u_est_at_inds, current_inds, size(u_est))
                    end
                end

                if any(isnan.(u_est))
                    if display
                        println("\tIMPORTNT: interpolation is full of NaNs, skip")
                    end
                    continue;
                end
                # @assert any(.!isnan.(u_est)) count(isnan, u_est), length(u_est)

                # linear interpolation
                @timeit_debug timer "image interpolation" begin
                    # source_reg = warp_img(source, -real(u_est), -imag(u_est))
                    source_reg = warp_img(source, -real(u_est), -imag(u_est), target)
                end

                if display
                    add_figs_sparse_pflap(figs, level, iter_repeat, u_est, new_estim_at_inds, source_reg, level_count, max_repeats, current_inds, fhs)
                end
            end
        end # "single filter pyramid level"
    end

    if display; return u_est, source_reg, figs; end
    return u_est, source_reg
end

function add_figs_sparse_pflap(figs, level, iter_repeat, u_est, new_estim_at_inds, source_reg, level_count, max_repeats, current_inds, fhs)
    level_iter_fhs = "(Level: $(string(level))/$(string(level_count))), (Iter: $(string(iter_repeat))/$(string(max_repeats))), (fhs: $fhs)"
    figs[level, iter_repeat, 1] = showflow(u_est, figtitle="U_EST " * level_iter_fhs)
    figs[level, iter_repeat, 2] = imgshow(source_reg, figtitle="SOURCE_REG " * level_iter_fhs, origin_bot=true)
    figs[level, iter_repeat, 3] = showflow(create_sparse_flow_from_sparse(new_estim_at_inds, current_inds, size(u_est)), disp_type=:sparse, figtitle="SPARSE Δ_U " * level_iter_fhs)
    figs[level, iter_repeat, 4] = showflow(create_sparse_flow_from_full(u_est, current_inds), figtitle="u_est_at_points " * level_iter_fhs)
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module Objective

import MathOptInterface

const MOI = MathOptInterface

include("bridge.jl")
include("map.jl")
include("single_bridge_optimizer.jl")

include("bridges/functionize.jl")
include("bridges/slack.jl")

"""
    add_all_bridges(model, ::Type{T}) where {T}

Add all bridges defined in the `Bridges.Objective` submodule to `model`.

The coefficient type used is `T`.
"""
function add_all_bridges(model, ::Type{T}) where {T}
    MOI.Bridges.add_bridge(model, FunctionizeBridge{T})
    MOI.Bridges.add_bridge(model, SlackBridge{T})
    return
end

end


text:
# ---
# title: 395. Longest Substring with At Least K Repeating Characters
# id: problem395
# author: Indigo
# date: 2021-02-27
# difficulty: Medium
# categories: Divide and Conquer, Recursion, Sliding Window
# link: <https://leetcode.com/problems/longest-substring-with-at-least-k-repeating-characters/description/>
# hidden: true
# ---
#
# Given a string `s` and an integer `k`, return _the length of the longest
# substring of_ `s` _such that the frequency of each character in this substring
# is greater than or equal to_ `k`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aaabb", k = 3
#     Output: 3
#     Explanation: The longest substring is "aaa", as 'a' is repeated 3 times.
#
#
# **Example 2:**
#
#
#
#     Input: s = "ababbc", k = 2
#     Output: 5
#     Explanation: The longest substring is "ababb", as 'a' is repeated 2 times and 'b' is repeated 3 times.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 104`
#   * `s` consists of only lowercase English letters.
#   * `1 <= k <= 105`
#
#
## @lc code=start
using LeetCode

function longest_substring_395(s::String, k::Int)
    function dfs(l::Int, r::Int)
        l > r && return 0
        cnt = counter(SubString(s, l, r))
        split_set = Int[]
        for i in l:r
            (cnt[s[i]] < k) && push!(split_set, i)
        end
        isempty(split_set) && return r - l + 1
        lft, res = l, 0
        for rt in split_set
            res = max(res, dfs(lft, rt - 1))
            lft = rt + 1
        end
        res
    end
    dfs(1, length(s))
end
## @lc code=end


text:
SRC_DIR = joinpath(Base.source_dir(), "..", "..", "src")
OUT_DIR = joinpath(Base.source_dir(), "..", "..", "output")
mkpath(OUT_DIR)

#include(joinpath(SRC_DIR, "BridgeSDEInference.jl"))
#using Main.BridgeSDEInference
include(joinpath(SRC_DIR, "BridgeSDEInference_for_tests.jl"))

using StaticArrays
using Distributions # to define priors
using Random        # to seed the random number generator

# Let's generate the data
# -----------------------
using Bridge
include(joinpath(SRC_DIR, "auxiliary", "data_simulation_fns.jl"))
Random.seed!(4)
# values taken from table 1 of Golightly and Wilkinson
θ_init = [0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1]
K = 10.0
Pˣ = Prokaryote(θ_init..., K)

x0, dt, T = ℝ{4}(7.0, 10.0, 4.0, 6.0), 1/5000, 20.0
tt = 0.0:dt:T
XX, _ = simulateSegment(ℝ{4}(0.0, 0.0, 0.0, 0.0), x0, Pˣ, tt)
skip = 2500
obs_time, obs_vals = XX.tt[1:skip:end], XX.yy[1:skip:end]

auxFlag = Val{(true, true, true, true)}()
P̃ = [ProkaryoteAux(θ_init..., K, t₀, u, T, v, auxFlag) for (t₀, T, u, v)
     in zip(obs_time[1:end-1], obs_time[2:end], obs[1:end-1], obs[2:end])]

Σdiagel = 10^-4
Σ = SMatrix{4,4}(1.0I)*Σdiagel
L = SMatrix{4,4}(1.0I)

setup = MCMCSetup(P˟, P̃, PartObs())
set_observations!(setup, [L for _ in P̃], [Σ for _ in P̃], obs, obs_time)
set_imputation_grid!(setup, 1/1000)
set_transition_kernels!(setup,
                        [RandomWalk([0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2],
                                    collect(1:8))],
                        0.7, false, 1,
                        (MetropolisHastingsUpdt(),))
set_priors!(setup,
            Priors((ImproperPrior(), ImproperPrior(), ImproperPrior(),
                    ImproperPrior(), ImproperPrior(), ImproperPrior(),
                    ImproperPrior(), ImproperPrior())),
            KnownStartingPt(x0)
            )
set_mcmc_params!(setup, 1*10^3, 1*10^0, 10^1, 10^0, 0)
set_blocking!(setup)
set_solver!(setup, Vern7(), NoChangePt())
initialise!(eltype(x0), setup)

Random.seed!(4)
out, elapsed = @timeit mcmc(setup)
display(out.accpt_tracker)

include(joinpath(SRC_DIR, DIR, "plotting_fns.jl"))
plot_chains(out; truth=[0.1, 0.7, 0.35, 0.2, 0.1, 0.9, 0.3, 0.1])
plot_paths(out; obs=(times=obs_time[2:end], vals=obs_vals[2:end], indices=1:4))


text:
using CompScienceMeshes, BEAST

h = 2π / 51; Γ = meshcircle(1.0, h)
X, Y = lagrangecxd0(Γ), lagrangec0d1(Γ)

κ = 1.0
S, Dᵀ = SingleLayer(κ), DoubleLayerTransposed(κ)
D, N  = DoubleLayer(κ), HyperSingular(κ)
I  = Identity()

# TM scattering
E = PlaneWaveDirichlet(κ, point(1.0,0.0))
H = PlaneWaveNeumann(  κ, point(1.0,0.0))

@hilbertspace u; @hilbertspace u′
@hilbertspace t; @hilbertspace t′

tm_efie = @discretise           S[u′,u] == E[u′] u∈X u′∈X
tm_mfie = @discretise (0.5I + Dᵀ)[t′,u] == H[t′] u∈X t′∈X

u1, u2 = solve(tm_efie), solve(tm_mfie)

# TE scattering
E = PlaneWaveNeumann(  κ, point(1.0, 0.0))
H = PlaneWaveDirichlet(κ, point(1.0, 0.0))

te_efie = @discretise          N[t,t′] == E[t′] t∈Y t′∈Y
te_mfie = @discretise (0.5I - D)[t,u′] == H[u′] t∈Y u′∈Y

t1, t2 = solve(te_efie), solve(te_mfie)

using Plots
plotlyjs()

nx = numfunctions(X);
Δα = 2π/nx; α = (collect(1:nx) - 0.5) * Δα
plt1 = Plots.plot(α, real(u1), c=:blue, label="TM-EFIE")
Plots.scatter!(α, real(u2), c=:red, m=:circle, label="TM-MFIE")
Plots.title!("current vs. angle")

ny = numfunctions(Y);
Δα = 2π/ny; α = collect(0:ny-1) * Δα
plt2 = Plots.plot(α, real(t1), c=:blue, label="TE-EFIE")
Plots.scatter!(α, real(t2), c=:red, m=:circle, label="TE-MFIE")
Plots.title!("current vs. angle")

Plots.plot(plt1,plt2)


text:
using DFControl

#This example shows how to setup and interact with the user_defaults.jl file which will get loaded each time you do `using DFControl`
setdefault_server("blabla@123.312.124")
getdefault_pseudodirs()
setdefault_pseudodir(:sssp, "~/pseudos/sssp")
setdefault_pseudodir(:pbesol, "~/pseudos/pbesol/PSEUDOPOTENTIALS")
setdefault_pseudodir(:pbesolrel, "~/pseudos/pbesolrel/PSEUDOPOTENTIALS")
configuredefault_pseudos()
#look at the user_defaults.jl file inside ~/.julia/config/DFControl

# stored default pseudos can be accessed through e.g.
getdefault_pseudo(:Si, :sssp)


text:
# kuramoto_8-18-19.jl
# Geez, have *I* been gone for a while. :P

#----------------------------------------------------------#

# Reading through
#
# https://docs.julialang.org/en/v1/manual/functions/
#
# now, so this'll be some practice with what's in there.

function f(x, y)
    x = 2*x
    y = 2*y
    x + y
end


text:
gear = [KPDSupply("map", 9, 150, 1),
        KPDSupply("compass", 13, 35, 1),
        KPDSupply("water", 153, 200, 2),
        KPDSupply("sandwich", 50, 60, 2),
        KPDSupply("glucose", 15, 60, 2),
        KPDSupply("tin", 68, 45, 3),
        KPDSupply("banana", 27, 60, 3),
        KPDSupply("apple", 39, 40, 3),
        KPDSupply("cheese", 23, 30, 1),
        KPDSupply("beer", 52, 10, 3),
        KPDSupply("suntan cream", 11, 70, 1),
        KPDSupply("camera", 32, 30, 1),
        KPDSupply("T-shirt", 24, 15, 2),
        KPDSupply("trousers", 48, 10, 2),
        KPDSupply("umbrella", 73, 40, 1),
        KPDSupply("waterproof trousers", 42, 70, 1),
        KPDSupply("waterproof overclothes", 43, 75, 1),
        KPDSupply("note-case", 22, 80, 1),
        KPDSupply("sunglasses", 7, 20, 1),
        KPDSupply("towel", 18, 12, 2),
        KPDSupply("socks", 4, 50, 1),
        KPDSupply("book", 30, 10, 2)]

pack = solve(gear, 400)
println("The hiker should pack: \n - ", join(pack, "\n - "))
println("\nPacked weight: ", sum(getfield.(pack, :weight)), " kg")
println("Packed value: ", sum(getfield.(pack, :value)), " €")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
show(io::IO, c::Colorant) = get(io, :compact, false) ? _showcompact(io, c) : _show(io, c)
show(io::IO, c::ColorantNormed) = show_normed(io, c)

# Nonparametric types
show_normed(io::IO, c::Gray24) = print(io, "Gray24(", gray(c), ')')
show_normed(io::IO, c::RGB24)  = print(io, "RGB24(", red(c), ',', green(c), ',', blue(c), ')')
show_normed(io::IO, c::ARGB32) = print(io, "ARGB32(", red(c), ',', green(c), ',', blue(c), ',', alpha(c), ')')

for N = 1:4
    component = N >= 3 ? (:comp1, :comp2, :comp3, :alpha) : (:comp1, :alpha)
    printargs = Array{Any}(undef, 2, N)
    for i = 1:N
        printargs[1,i] = :(show(io, $(component[i])(c)))
        chr = i < N ? ',' : ')'
        printargs[2,i] = :(print(io, $chr))
    end
    @eval begin
        function _show(io::IO, c::Colorant{T,$N}) where T
            print(io, colorant_string(typeof(c)), "{", T, "}(")
            $(printargs[:]...)
        end
    end
    for i = 1:N
        printargs[1,i] = :(show(IOContext(io, :compact => true), $(component[i])(c)))
    end
    @eval begin
        function _showcompact(io::IO, c::Colorant{T,$N}) where T
            print(io, colorant_string(typeof(c)), "{", T, "}(")
            $(printargs[:]...)
        end
        # Special handling for Normed types: don't print the giant type name
        function show_normed(io::IO, c::Colorant{FixedPointNumbers.Normed{T,f},$N}) where {T,f}
            print(io, colorant_string(typeof(c)), '{')
            FixedPointNumbers.showtype(io, eltype(typeof(c)))
            print(io, "}(")
            $(printargs[:]...)
        end
    end
end

function Base.showarg(io::IO, a::Array{C}, toplevel) where C<:Colorant
    toplevel || print(io, "::")
    print(io, "Array{")
    colorant_string_with_eltype(io, C)
    print(io, ",$(ndims(a))}")
    toplevel && print(io, " with eltype ", C)
end


text:
include("../src/read/read_datasets.jl")
include("../src/utils.jl")
include("../src/pq/PQ.jl");
include("../src/lsq_sparse/LSQ_SPGL1.jl");
include("../src/linscan/Linscan.jl");

function demo_lsq_sparse(
  dataset_name="SIFT1M",
  nread::Integer=Int(1e4))

  # === Hyperparams ===
  m       = 7 # In LSQ we use m-1 codebooks
  h       = 256
  verbose = true
  nquery  = Int(1e4)
  knn     = Int(1e3) # Compute recall up to
  b       = Int( log2(h) * m )
  niter   = 10

  # === PQ initialization ===
  x_train           = read_dataset(dataset_name, nread )
  d, _              = size( x_train )
  C, B, train_error = train_pq(x_train, m, h, verbose)
  @printf("Error after PQ is %e\n", train_error)

  # === LSQ sparse train ===
  ilsiter = 8
  icmiter = 4
  randord = true
  npert   = 4
  randord = true
  S       = d*h # SLSQ1 in the paper. Use S = d*h + (d.^2) for SLSQ2
  tau     = 0.7 # 0.7 for SLSQ1. Use 0.9 for SLSQ2

  # Multiply tau times the l1 norm of the PQ solution
  taus    = zeros( d )
  subdims = splitarray( 1:d, m )
  for i = 1:m
    taus[ subdims[i] ] += sum( abs(C[i]), 2 ) .* tau
  end
  tau = sum( taus )

  spgl1_path = joinpath( pwd(), "matlab/spgl1")
  C, B, R, train_error, cbnorms, objs =
    train_lsq_sparse(x_train, m, h, niter, ilsiter, icmiter, randord, npert, S, tau,
    B, C, eye(Float32, d), spgl1_path)
  cbnorms = vec( cbnorms[:] )

  # === Encode the base set ===
  nread_base   = Int(1e6)
  x_base       = read_dataset(dataset_name * "_base", nread_base )
  B_base       = randinit(nread_base, m, h) # initialize B at random

  ilsiter_base = 16 # LSQ-16 in the paper
  for i = 1:ilsiter_base
    @printf("Iteration %02d / %02d\n", i, ilsiter_base)
    @time B_base = encoding_icm( x_base, B_base, C, icmiter, randord, npert, verbose )
  end
  base_error = qerror( x_base, B_base, C )
  @printf("Error in base is %e\n", base_error)

  # Compute and quantize the database norms
  B_base_norms = quantize_norms( B_base, C, cbnorms )
  db_norms     = vec( cbnorms[ B_base_norms ] )

  # === Compute recall ===
  x_query = read_dataset( dataset_name * "_query", nquery, verbose )
  gt      = read_dataset( dataset_name * "_groundtruth", nquery, verbose )
  if dataset_name == "SIFT1M" || dataset_name == "GIST1M"
    gt = gt + 1;
  end
  gt           = convert( Vector{UInt32}, gt[1,1:nquery] )
  B_base       = convert( Matrix{UInt8}, B_base-1 )
  B_base_norms = convert( Vector{UInt8}, B_base_norms-1 )

  print("Querying m=$m ... ")
  @time dists, idx = linscan_lsq( B_base, x_query, C, db_norms, eye(Float32, d), knn )
  println("done")

  idx = convert( Matrix{UInt32}, idx );
  rec = eval_recall( gt, idx, knn )

end

# train
demo_lsq_sparse()

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
using Test, LinearAlgebra
using OrdinaryDiffEq
using DiffEqSensitivity, StochasticDiffEq, DiffEqBase
using ForwardDiff, Calculus, ReverseDiff
using Random
import Tracker, Zygote

@info "SDE Adjoints"

seed = 100
Random.seed!(seed)
abstol = 1e-4
reltol = 1e-4

u₀ = [0.5]
tstart = 0.0
tend = 0.1
dt = 0.005
trange = (tstart, tend)
t = tstart:dt:tend
tarray = collect(t)

function g(u,p,t)
  sum(u.^2.0/2.0)
end

function dg!(out,u,p,t,i)
  (out.=-u)
end

p2 = [1.01,0.87]

@testset "SDE oop Tests" begin
  f_oop_linear(u,p,t) = p[1]*u
  σ_oop_linear(u,p,t) = p[2]*u

  function f_oop_linear(u::Tracker.TrackedArray,p,t)
    dx = p[1]*u[1]
    Tracker.collect([dx])
  end

  function σ_oop_linear(u::Tracker.TrackedArray,p,t)
    dx = p[2]*u[1]
    Tracker.collect([dx])
  end

  p = [1.01,0.0]

  # generate ODE adjoint results

  prob_oop_ode = ODEProblem(f_oop_linear,u₀,(tstart,tend),p)
  sol_oop_ode = solve(prob_oop_ode,Tsit5(),saveat=t,abstol=abstol,reltol=reltol)
  res_ode_u0, res_ode_p = adjoint_sensitivities(sol_oop_ode,Tsit5(),dg!,t
    ,abstol=abstol,reltol=reltol,sensealg=BacksolveAdjoint())

  function G(p)
    tmp_prob = remake(prob_oop_ode,u0=eltype(p).(prob_oop_ode.u0),p=p,
                    tspan=eltype(p).(prob_oop_ode.tspan),abstol=abstol, reltol=reltol)
    sol = solve(tmp_prob,Tsit5(),saveat=Array(t),abstol=abstol, reltol=reltol)
    res = g(sol,p,nothing)
  end
  res_ode_forward = ForwardDiff.gradient(G,p)
  #res_ode_reverse = ReverseDiff.gradient(G,p)

  res_ode_trackeru0, res_ode_trackerp = Zygote.gradient((u0,p)->sum(solve(prob_oop_ode,Tsit5();u0=u0,p=p,abstol=abstol,reltol=reltol,saveat=Array(t),sensealg=TrackerAdjoint()).^2.0/2.0),u₀,p)

  @test isapprox(res_ode_forward[1], sum(@. u₀^2*exp(2*p[1]*t)*t), rtol = 1e-4)
  #@test isapprox(res_ode_reverse[1], sum(@. u₀^2*exp(2*p[1]*t)*t), rtol = 1e-4)
  @test isapprox(res_ode_p'[1], sum(@. u₀^2*exp(2*p[1]*t)*t), rtol = 1e-4)
  @test isapprox(res_ode_p', res_ode_trackerp, rtol = 1e-4)

  # SDE adjoint results (with noise == 0, so should agree with above)

  Random.seed!(seed)
  prob_oop_sde = SDEProblem(f_oop_linear,σ_oop_linear,u₀,trange,p)
  sol_oop_sde = solve(prob_oop_sde,RKMil(interpretation=:Stratonovich),dt=1e-4,adaptive=false,save_noise=true)
  res_sde_u0, res_sde_p = adjoint_sensitivities(sol_oop_sde,
    EulerHeun(),dg!,t,dt=1e-2,sensealg=BacksolveAdjoint())

  @info res_sde_p

  function GSDE1(p)
    Random.seed!(seed)
    tmp_prob = remake(prob_oop_sde,u0=eltype(p).(prob_oop_sde.u0),p=p,
                    tspan=eltype(p).(prob_oop_sde.tspan))
    sol = solve(tmp_prob,RKMil(interpretation=:Stratonovich),dt=tend/10000,adaptive=false,sensealg=DiffEqBase.SensitivityADPassThrough(),saveat=t)
    A = convert(Array,sol)
    res = g(A,p,nothing)
  end
  res_sde_forward = ForwardDiff.gradient(GSDE1,p)
  res_sde_reverse = ReverseDiff.gradient(GSDE1,p)

  res_sde_trackeru0, res_sde_trackerp = Zygote.gradient((u0,p)->sum(solve(prob_oop_sde,RKMil(interpretation=:Stratonovich),dt=tend/1400,adaptive=false,u0=u0,p=p,saveat=Array(t),sensealg=TrackerAdjoint()).^2.0/2.0),u₀,p)

  noise = vec((@. sol_oop_sde.W(tarray)))
  Wfix = [W[1][1] for W in noise]
  @test isapprox(res_sde_forward[1], sum(@. u₀^2*exp(2*p[1]*t)*t), rtol = 1e-4)
  @test isapprox(res_sde_reverse[1], sum(@. u₀^2*exp(2*p[1]*t)*t), rtol = 1e-4)
  @test isapprox(res_sde_p'[1], sum(@. u₀^2*exp(2*p[1]*t)*t), rtol = 1e-4)
  @test isapprox(res_sde_p'[2], sum(@. (Wfix)*u₀^2*exp(2*(p[1])*tarray+2*p[2]*Wfix)), rtol = 1e-4)
  @test isapprox(res_sde_p'[1], res_sde_trackerp[1], rtol = 1e-4)

  # SDE adjoint results (with noise != 0)

  Random.seed!(seed)
  prob_oop_sde2 = SDEProblem(f_oop_linear,σ_oop_linear,u₀,trange,p2)
  sol_oop_sde2 = solve(prob_oop_sde2,RKMil(interpretation=:Stratonovich),dt=tend/1e6,adaptive=false,save_noise=true)

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol_oop_sde2,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint())

  @info res_sde_p2

  # test consitency for different switches for the noise Jacobian
  res_sde_u02a, res_sde_p2a = adjoint_sensitivities(sol_oop_sde2,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=false))

  @test isapprox(res_sde_u02, res_sde_u02a, rtol = 1e-6)
  @test isapprox(res_sde_p2, res_sde_p2a, rtol = 1e-6)

  @info res_sde_p2a

  res_sde_u02a, res_sde_p2a = adjoint_sensitivities(sol_oop_sde2,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

  @test isapprox(res_sde_u02, res_sde_u02a, rtol = 1e-6)
  @test isapprox(res_sde_p2, res_sde_p2a, rtol = 1e-6)

  @info res_sde_p2a

  res_sde_u02a, res_sde_p2a = adjoint_sensitivities(sol_oop_sde2,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise()))

  @test isapprox(res_sde_u02, res_sde_u02a, rtol = 1e-6)
  @test isapprox(res_sde_p2, res_sde_p2a, rtol = 1e-6)

  @info res_sde_p2a

  function GSDE2(p)
    Random.seed!(seed)
    tmp_prob = remake(prob_oop_sde2,u0=eltype(p).(prob_oop_sde2.u0),p=p,
                      tspan=eltype(p).(prob_oop_sde2.tspan)
                      #,abstol=abstol, reltol=reltol
                      )
    sol = solve(tmp_prob,RKMil(interpretation=:Stratonovich),dt=tend/1e6,adaptive=false,sensealg=DiffEqBase.SensitivityADPassThrough(),saveat=Array(t))
    A = convert(Array,sol)
    res = g(A,p,nothing)
  end
  res_sde_forward2 = ForwardDiff.gradient(GSDE2,p2)
  res_sde_reverse2 = ReverseDiff.gradient(GSDE2,p2)

  Random.seed!(seed)
  res_sde_trackeru02, res_sde_trackerp2 = Zygote.gradient((u0,p)->sum(solve(prob_oop_sde2,RKMil(interpretation=:Stratonovich),dt=tend/1e3,adaptive=false,u0=u0,p=p,saveat=Array(t),sensealg=TrackerAdjoint()).^2.0/2.0),u₀,p2)

  Wfix = [sol_oop_sde2.W(t)[1][1] for t in tarray]
  resp1 = sum(@. tarray*u₀^2*exp(2*(p2[1])*tarray+2*p2[2]*Wfix))
  resp2 = sum(@. (Wfix)*u₀^2*exp(2*(p2[1])*tarray+2*p2[2]*Wfix))
  resp = [resp1, resp2]

  @test isapprox(res_sde_forward2, resp, rtol = 2e-6)
  @test isapprox(res_sde_reverse2, resp, rtol = 2e-6)
  @test isapprox(res_sde_trackerp2, resp, rtol = 4e-1)

  @test isapprox(res_sde_p2', res_sde_forward2, rtol = 1e-6)
  @test isapprox(res_sde_p2', resp, rtol = 2e-6)

  @info "ForwardDiff" res_sde_forward2
  @info "ReverseDiff" res_sde_reverse2
  @info "Exact" resp
  @info "Adjoint SDE" res_sde_p2

  # consistency check with respect to tracker
  Random.seed!(seed)
  prob_oop_sde2 = SDEProblem(f_oop_linear,σ_oop_linear,u₀,trange,p2)
  sol_oop_sde2 = solve(prob_oop_sde2,RKMil(interpretation=:Stratonovich),dt=tend/1e3,adaptive=false,save_noise=true)
  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol_oop_sde2,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e3,adaptive=false,sensealg=BacksolveAdjoint())

  @test isapprox(res_sde_p2', res_sde_trackerp2, rtol = 3e-4)

  # Free memory to help Travis

  Wfix = nothing
  res_sde_forward2 = nothing
  res_sde_reverse2 = nothing
  resp = nothing
  res_sde_trackerp2 = nothing
  res_sde_u02 = nothing
  sol_oop_sde2 = nothing
  res_sde_u02a = nothing
  res_sde_p2a = nothing
  res_sde_p2 = nothing
  sol_oop_sde = nothing
  GC.gc()

  # SDE adjoint results with diagonal noise

  Random.seed!(seed)
  prob_oop_sde2 = SDEProblem(f_oop_linear,σ_oop_linear,[u₀;u₀;u₀],trange,p2)
  sol_oop_sde2 = solve(prob_oop_sde2,EulerHeun(),dt=tend/1e6,adaptive=false,save_noise=true)

  @info "Diagonal Adjoint"

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol_oop_sde2,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint())

  res_sde_u03, res_sde_p3 = adjoint_sensitivities(sol_oop_sde2,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=false))

  res_sde_u04, res_sde_p4 = adjoint_sensitivities(sol_oop_sde2,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

  res_sde_u05, res_sde_p5 = adjoint_sensitivities(sol_oop_sde2,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise()))

  @info  res_sde_p2

  sol_oop_sde2 = nothing
  GC.gc()

  @test isapprox(res_sde_p2, res_sde_p3, rtol = 1e-7)
  @test isapprox(res_sde_p2, res_sde_p4, rtol = 1e-7)
  @test isapprox(res_sde_p2, res_sde_p5, rtol = 1e-7)

  @test isapprox(res_sde_u02, res_sde_u03, rtol = 1e-7)
  @test isapprox(res_sde_u02, res_sde_u04, rtol = 1e-7)
  @test isapprox(res_sde_u02, res_sde_u05, rtol = 1e-7)

  @info "Diagonal ForwardDiff"
  res_sde_forward2 = ForwardDiff.gradient(GSDE2,p2)
  #@info "Diagonal ReverseDiff"
  #res_sde_reverse2 = ReverseDiff.gradient(GSDE2,p2)

  #@test isapprox(res_sde_forward2, res_sde_reverse2, rtol = 1e-6)
  @test isapprox(res_sde_p2', res_sde_forward2, rtol = 1e-3)
  #@test isapprox(res_sde_p2', res_sde_reverse2, rtol = 1e-3)

  # u0
  function GSDE3(u)
    Random.seed!(seed)
    tmp_prob = remake(prob_oop_sde2,u0=u,p=eltype(p).(prob_oop_sde2.p),
                      tspan=eltype(p).(prob_oop_sde2.tspan)
                      #,abstol=abstol, reltol=reltol
                      )
    sol = solve(tmp_prob,RKMil(interpretation=:Stratonovich),dt=tend/1e6,adaptive=false,saveat=Array(t))
    A = convert(Array,sol)
    res = g(A,p,nothing)
  end

  @info "ForwardDiff u0"
  res_sde_forward2 = ForwardDiff.gradient(GSDE3,[u₀;u₀;u₀])

  @test isapprox(res_sde_u02, res_sde_forward2, rtol = 1e-4)

end

##
## Inplace
##
@testset "SDE inplace Tests" begin
  f!(du,u,p,t) = du.=p[1]*u
  σ!(du,u,p,t) = du.=p[2]*u

  Random.seed!(seed)
  prob_sde = SDEProblem(f!,σ!,u₀,trange,p2)
  sol_sde = solve(prob_sde,RKMil(interpretation=:Stratonovich),dt=tend/1e6,adaptive=false, save_noise=true)

  # Wfix = [sol_sde.W(t)[1][1] for t in tarray]
  # resp1 = sum(@. tarray*u₀^2*exp(2*(p2[1])*tarray+2*p2[2]*Wfix))
  # resp2 = sum(@. (Wfix)*u₀^2*exp(2*(p2[1])*tarray+2*p2[2]*Wfix))
  # resp = [resp1, resp2]
  # resu0 = sum(@. u₀*exp(2*(p2[1])*tarray+2*p2[2]*Wfix))

  function GSDE(p)
    Random.seed!(seed)
    tmp_prob = remake(prob_sde,u0=eltype(p).(prob_sde.u0),p=p,
                      tspan=eltype(p).(prob_sde.tspan))
    sol = solve(tmp_prob,RKMil(interpretation=:Stratonovich),dt=tend/1e6,adaptive=false,saveat=Array(t))
    A = convert(Array,sol)
    res = g(A,p,nothing)
  end

  res_sde_forward = ForwardDiff.gradient(GSDE,p2)

  res_sde_u0, res_sde_p = adjoint_sensitivities(sol_sde,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint())

  @test isapprox(res_sde_p', res_sde_forward, rtol = 1e-5)

  @info res_sde_p

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol_sde,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=false))

  @info res_sde_p2

  @test isapprox(res_sde_p, res_sde_p2, rtol = 1e-5)
  @test isapprox(res_sde_u0, res_sde_u02, rtol = 1e-5)

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol_sde,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

  @info res_sde_p2

  @test isapprox(res_sde_p, res_sde_p2, rtol = 1e-5) # not broken here because it just uses the vjps
  @test isapprox(res_sde_u0 ,res_sde_u02, rtol = 1e-5)

  res_sde_u02, res_sde_p2 = adjoint_sensitivities(sol_sde,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise()))

  @info res_sde_p2

  @test isapprox(res_sde_p, res_sde_p2, rtol = 1e-10)
  @test isapprox(res_sde_u0 ,res_sde_u02, rtol = 1e-10)

  # diagonal noise

  #compare with oop version
  f_oop_linear(u,p,t) = p[1]*u
  σ_oop_linear(u,p,t) = p[2]*u
  Random.seed!(seed)
  prob_oop_sde = SDEProblem(f_oop_linear,σ_oop_linear,[u₀;u₀;u₀],trange,p2)
  sol_oop_sde = solve(prob_oop_sde,EulerHeun(),dt=tend/1e6,adaptive=false,save_noise=true)
  res_oop_u0, res_oop_p = adjoint_sensitivities(sol_oop_sde,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint())

  @info res_oop_p

  Random.seed!(seed)
  prob_sde = SDEProblem(f!,σ!,[u₀;u₀;u₀],trange,p2)
  sol_sde = solve(prob_sde,EulerHeun(),dt=tend/1e6,adaptive=false,save_noise=true)

  res_sde_u0, res_sde_p = adjoint_sensitivities(sol_sde,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint())

  @test isapprox(res_sde_p, res_oop_p, rtol = 1e-6)
  @test isapprox(res_sde_u0 ,res_oop_u0, rtol = 1e-6)

  @info res_sde_p

  res_sde_u0, res_sde_p = adjoint_sensitivities(sol_sde,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=false))

  @test isapprox(res_sde_p, res_oop_p, rtol = 1e-6)
  @test isapprox(res_sde_u0 ,res_oop_u0, rtol = 1e-6)

  @info res_sde_p

  res_sde_u0, res_sde_p = adjoint_sensitivities(sol_sde,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ZygoteNoise()))

  @test_broken isapprox(res_sde_p, res_oop_p, rtol = 1e-6)
  @test isapprox(res_sde_u0 ,res_oop_u0, rtol = 1e-6)

  @info res_sde_p

  res_sde_u0, res_sde_p = adjoint_sensitivities(sol_sde,EulerHeun(),dg!,Array(t)
      ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint(noise=DiffEqSensitivity.ReverseDiffNoise()))

  @test isapprox(res_sde_p, res_oop_p, rtol = 1e-6)
  @test isapprox(res_sde_u0 ,res_oop_u0, rtol = 1e-6)

  @info res_sde_p
end

# scalar noise
@testset "SDE scalar noise tests" begin
  using DiffEqNoiseProcess

  f!(du,u,p,t) = (du .= p[1]*u)
  σ!(du,u,p,t) = (du .= p[2]*u)

  @info "scalar SDE"

  Random.seed!(seed)
  W = WienerProcess(0.0,0.0,0.0)
  u0 = rand(2)

  linear_analytic_strat(u0,p,t,W) = @.(u0*exp(p[1]*t+p[2]*W))

  prob = SDEProblem(SDEFunction(f!,σ!,analytic=linear_analytic_strat),σ!,u0,trange,p2,
    noise=W
    )
  sol = solve(prob,EulerHeun(), dt=tend/1e6, save_noise=true)

  @test isapprox(sol.u_analytic,sol.u, atol=2e-5)

  res_sde_u0, res_sde_p = adjoint_sensitivities(sol,EulerHeun(),dg!,Array(t)
    ,dt=tend/1e6,adaptive=false,sensealg=BacksolveAdjoint())

  function compute_grads(sol, scale=1.0)
    xdis = sol(tarray)
    helpu1 = [u[1] for u in xdis.u]
    tmp1 = sum((@. xdis.t*helpu1*helpu1))

    Wtmp = [sol.W(t)[1][1] for t in tarray]
    tmp2 = sum((@. Wtmp*helpu1*helpu1))

    tmp3 = sum((@. helpu1*helpu1))/helpu1[1]

    return [tmp3, scale*tmp3], [tmp1*(1.0+scale^2), tmp2*(1.0+scale^2)]
  end

  @test isapprox(compute_grads(sol, u0[2]/u0[1])[2], res_sde_p', atol=1e-6)
  @test isapprox(compute_grads(sol, u0[2]/u0[1])[1], res_sde_u0, atol=1e-6)
end


text:
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Imports
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

using Main.JuliaSub: countTextualConstr
using Main.JuliaSub: subtc, suptc
using Main.JuliaSub: TxtConstrStat, LBValsFreq, LBStat, FileLBInfo
using Main.JuliaSub: extractLowerBound, extractLowerBounds
using Main.JuliaSub: nonVacuous, lbStatInfo, lbFileInfo
using Main.JuliaSub: FilesLBInfo, PackageStat
using Main.JuliaSub: processPkg, processPkgsDir

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Aux values and functions
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Tests
#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

# Analysis relies on the fact that there are only few patterns
# where a lower bound can appear.
# If language changes and there are more, the analysis needs
# to be extended
@testset "lb-analysis :: type bounds format" begin
    @test (Meta.parse("f(x::T) where T = 0"); true)
    @test (Meta.parse("f(x::T) where T<:Number = 0"); true)
    @test (Meta.parse("f(x::T) where T>:Int = 0"); true)
    @test (Meta.parse("f(x::T) where Int<:T<:Number = 0"); true)
    @test (Meta.parse("Vector{>:Int}"); true)

    @test_throws Exception (eval(Meta.parse("f(x::T) where Number>:T>:Int = 0")); true)
    @test_throws UndefVarError (eval(Meta.parse("f(x::T) where Int<:T = 0")); true)
end

@testset "lb-analysis :: sub/sup symbols" begin
    @test countTextualConstr(subtc, "")  == 0
    @test countTextualConstr(subtc, ":") == 0
    @test countTextualConstr(subtc, "<") == 0
    @test countTextualConstr(subtc, ">:")  == 0
    @test countTextualConstr(suptc, "  ab ") == 0
    @test countTextualConstr(suptc, "> :") == 0
    @test countTextualConstr(suptc, "<:")  == 0

    @test countTextualConstr(subtc, "<:")   == 1
    @test countTextualConstr(subtc, "<:<")  == 1
    @test countTextualConstr(suptc, ">:")   == 1
    @test countTextualConstr(suptc, "{>:}") == 1
    @test countTextualConstr(suptc, " >::") == 1

    @test countTextualConstr(subtc, "where T <: Int") == 1
    @test countTextualConstr(subtc, "where Int<:T<:Number") == 2
    @test countTextualConstr(suptc, "where T >:   Int ") == 1

    @test countTextualConstr("<: T>:Int Nothing<:S<:Number") ==
        TxtConstrStat(3, 1)
    @test countTextualConstr("< : T>:Int >:>: S>:Number </:") ==
        TxtConstrStat(0, 4)
end

@testset "lb-analysis :: capture direct lower bound" begin
    @test extractLowerBound(:(<:Int)) == nothing
    @test extractLowerBound(:(T <: Number)) == nothing
    @test extractLowerBound(:(Number>:T>:Int)) == nothing

    @test extractLowerBound(:(>:Int)) == :Int
    @test extractLowerBound(:(T>:Int)) == :Int
    @test extractLowerBound(:( T >: Nothing )) == :Nothing
    @test extractLowerBound(:(Int <: T <: Bool)) == :Int
end

@testset "lb-analysis :: capture lower bounds" begin
    @test extractLowerBounds(:(g(xs::Vector{<:Real}) = 0)) ==
        Multiset()
    @test extractLowerBounds(:(f(x::T) where T = 0)) ==
        Multiset()

    @test extractLowerBounds(:(g(xs::Vector{>:Int}) = 0)) ==
        Multiset(:Int)
    @test extractLowerBounds(:(g(xs::Dict{>:String, >:Int}) = 0)) ==
        Multiset(:String, :Int)

    @test extractLowerBounds(:(f(x::T) where T>:Int)) ==
        Multiset(:Int)
    @test extractLowerBounds(:(f(x::T) where {T>:Int, Nothing<:S<:Number} = 0)) ==
        Multiset(:Int, :Nothing)
    @test extractLowerBounds(:(f(x::T) where Int<:T where S<:Number = 0)) ==
        Multiset()
    @test extractLowerBounds(:(f(x::T) where T>:Int where S>:Int = 0)) ==
        Multiset(:Int, :Int)
    @test extractLowerBounds(:(f(x::T) where Int<:T<:S where S<:Number = 0)) ==
        Multiset(:Int)
    @test extractLowerBounds(:(f(x::T) where Int<:T<:S where S<:Number where U>:Bool = 0)) ==
        Multiset(:Int, :Bool)
    @test extractLowerBounds(:(f(x::T) where Int<:T<:S where {S<:Number, Int<:Q<:Nothing} where U>:Bool = 0)) ==
        Multiset(:Int, :Int, :Bool)
    @test extractLowerBounds(quote
            f(x) = 5
            f(x, y::T) where T>:S where S = x
            const X = 7
            g(x::T) where {Int<:T<:Number, S<:Bool} = x
            h(d::Dict{K,>:Int}) where K<:Vector{>:Bool} :: Vector{>:Missing} = d
        end) ==
        Multiset(:Int, :S, :Int, :Bool, :Missing)
end

@testset "lb-analysis :: nonVacuous" begin
    @test nonVacuous(TxtConstrStat(1, 0))
    @test nonVacuous(TxtConstrStat(0, 6))
    @test nonVacuous(TxtConstrStat(4, 1))
    @test !nonVacuous(TxtConstrStat(0, 0))
end

@testset "lb-analysis :: lower-bounds statistics" begin
    @test lbStatInfo("f(x::T) where T>:Int = 0") ==
        LBStat(1, 1, LBValsFreq(:Int))
    @test lbStatInfo("f(x::T) where {T>:Int, Int<:S<:Number} = 0") ==
        LBStat(2, 1, LBValsFreq(:Int, :Int))
    @test lbStatInfo("f(x::T) where {T>:Int, Int<:S<:Number, Q>:Missing} = 0") ==
        LBStat(3, 2, LBValsFreq(:Int, :Int, :Missing))
end

@testset "lb-analysis :: lower-bounds file statistics" begin
    @test lbFileInfo("")  == FileLBInfo(TxtConstrStat())
    @test lbFileInfo("<") == FileLBInfo(TxtConstrStat())
    @test lbFileInfo(">:") ==
        FileLBInfo(TxtConstrStat(0,1), LBStat(LBValsFreq()))
    @test lbFileInfo("f(x::T) where {T>:Int, Int<:S<:Number, Q>:Missing} = 0") ==
        FileLBInfo(TxtConstrStat(2,2), LBStat(LBValsFreq(:Int, :Int, :Missing)))

    @test lbFileInfo(testFilePath("fRet10.jl"); isPath=true) ==
        FileLBInfo(TxtConstrStat(3,4), LBStat(LBValsFreq(:Int, :Bool, :Nothing, :String, :T)))

    failedFileLBInfo = lbFileInfo("f(x::T) wher T>:Int = 0")
    @test failedFileLBInfo.txtStat == TxtConstrStat(0, 1)
    @test isa(failedFileLBInfo.err, Base.Meta.ParseError)
    @test failedFileLBInfo.lbStat == nothing
end

@testset "lb-analysis :: lower-bounds package statistics" begin
    pkgBad = "FAKE-bad.jl"
    pkgStatBad = PackageStat(pkgBad, false)
    @test processPkg(testFilePath(joinpath("pkgs", pkgBad)), pkgBad) ==
        pkgStatBad

    pkgGood = "FAKE-fRet10.jl"
    pkgStatGood = PackageStat(
        pkgGood, true, 3, 0, 2,
        FilesLBInfo(
            joinpath("src", "fRet10.jl") =>
                FileLBInfo(TxtConstrStat(3,4), LBStat(LBValsFreq(:Int, :Bool, :Nothing, :String, :T))),
            joinpath("src", "id.jl") =>
                FileLBInfo(TxtConstrStat(0,2), LBStat(LBValsFreq(:Any, :Int)))),
        LBStat(LBValsFreq(:Int,:Int, :Bool, :Nothing, :Any, :String, :T)))
    @test processPkg(testFilePath(joinpath("pkgs", pkgGood)), pkgGood) ==
        pkgStatGood

    (badPkgs, goodPkgs, totalStat) = processPkgsDir(testFilePath("pkgs"))
    getName(pkg :: PackageStat) = pkg.name
    @test badPkgs == [pkgStatBad]
    @test sort(goodPkgs, by=getName) ==
        sort([pkgStatGood, PackageStat("FAKE-empty.jl", true)], by=getName)
    @test totalStat == LBStat(7, 6, LBValsFreq(:Int,:Int, :Bool, :Nothing, :Any, :String, :T))
end

max_stars_repo_name: prl-julia/julia-sub

max_issues_repo_name: prl-julia/julia-sub

max_forks_repo_name: prl-julia/julia-sub


text:
# 2D FWI on Overthrust model with SPG using minConf library
# Author: Philipp Witte, pwitte@eoas.ubc.ca
# Date: December 2017
#

using Statistics, Random, LinearAlgebra, PyPlot
using JUDI, SlimOptim, HDF5, SegyIO

# Load starting model
n,d,o,m0 = read(h5open("../../data/overthrust_model.h5","r"), "n", "d", "o", "m0")
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)

# Bound constraints
v0 = sqrt.(1f0 ./ model0.m)
vmin = ones(Float32,model0.n) .* 1.3f0
vmax = ones(Float32,model0.n) .* 6.5f0
vmin[:,1:21] .= v0[:,1:21]   # keep water column fixed
vmax[:,1:21] .= v0[:,1:21]

# Slowness squared [s^2/km^2]
mmin = vec((1f0 ./ vmax).^2)
mmax = vec((1f0 ./ vmin).^2)

# Load data
block = segy_read("../../data/overthrust_shot_records.segy")
d_obs = judiVector(block)

# Set up wavelet
src_geometry = Geometry(block; key="source")
wavelet = ricker_wavelet(src_geometry.t[1],src_geometry.dt[1],0.008f0)  # 8 Hz wavelet
q = judiVector(src_geometry,wavelet)

############################### FWI ###########################################

# Optimization parameters
fevals = 16
batchsize = 8

# Objective function for minConf library
count = 0
function objective_function(x)
    model0.m .= reshape(x,model0.n);

    # fwi function value and gradient
    i = randperm(d_obs.nsrc)[1:batchsize]
    fval, grad = fwi_objective(model0, q[i], d_obs[i])
    grad = .125f0*grad/maximum(abs.(grad))  # scale for line search

    global count; count+= 1
    return fval, vec(grad.data)
end

# Bound projection
ProjBound(x) = median([mmin x mmax]; dims=2)

# FWI with SPG
options = spg_options(verbose=3, maxIter=fevals, memory=3)
sol = spg(objective_function, vec(model0.m), ProjBound, options)

# Plot result
imshow(reshape(sqrt.(1f0 ./ sol.sol), model0.n)', extent=[0, 10, 3, 0])
xlabel("Lateral position [km]")
ylabel("Depth [km]")

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
# Julia wrapper for source: julia/src/llvm-api.cpp

# initialization functions

function LLVMInitializeAllTargetInfos()
    ccall(:LLVMExtraInitializeAllTargetInfos,Cvoid,())
end

function LLVMInitializeAllTargets()
    ccall(:LLVMExtraInitializeAllTargets,Cvoid,())
end

function LLVMInitializeAllTargetMCs()
    ccall(:LLVMExtraInitializeAllTargetMCs,Cvoid,())
end

function LLVMInitializeAllAsmPrinters()
    ccall(:LLVMExtraInitializeAllAsmPrinters,Cvoid,())
end

function LLVMInitializeAllAsmParsers()
    ccall(:LLVMExtraInitializeAllAsmParsers,Cvoid,())
end

function LLVMInitializeAllDisassemblers()
    ccall(:LLVMExtraInitializeAllDisassemblers,Cvoid,())
end

function LLVMInitializeNativeTarget()
    ccall(:LLVMExtraInitializeNativeTarget,LLVMBool,())
end

function LLVMInitializeNativeAsmPrinter()
    ccall(:LLVMExtraInitializeNativeAsmPrinter,LLVMBool,())
end

function LLVMInitializeNativeAsmParser()
    ccall(:LLVMExtraInitializeNativeAsmParser,LLVMBool,())
end

function LLVMInitializeNativeDisassembler()
    ccall(:LLVMExtraInitializeNativeDisassembler,LLVMBool,())
end

# infrastructure for writing LLVM passes in Julia

mutable struct LLVMOpaquePass
end

const LLVMPassRef = Ptr{LLVMOpaquePass}

function LLVMAddPass(PM, P)
    ccall(:LLVMExtraAddPass,Cvoid,
        (LLVMPassManagerRef, LLVMPassRef),
        PM, P)
end

function LLVMCreateModulePass(Name, Callback)
    ccall(:LLVMExtraCreateModulePass,LLVMPassRef,
        (Cstring, Any),
        Name, Callback)
end

function LLVMCreateFunctionPass(Name, Callback)
    ccall(:LLVMExtraCreateFunctionPass,LLVMPassRef,
        (Cstring, Any),
        Name, Callback)
end

function LLVMCreateBasicBlockPass(Name, Callback)
    ccall(:LLVMExtraCreateBasicBlockPass,LLVMPassRef,
        (Cstring, Any),
        Name, Callback)
end

function LLVMCreateModulePass2(Name, Callback, Data)
    ccall(:LLVMExtraCreateModulePass2,LLVMPassRef,
        (Cstring, Ptr{Cvoid}, Ptr{Cvoid}),
        Name, Callback, Data)
end

function LLVMCreateFunctionPass2(Name, Callback, Data)
    ccall(:LLVMExtraCreateFunctionPass2,LLVMPassRef,
        (Cstring, Ptr{Cvoid}, Ptr{Cvoid}),
        Name, Callback, Data)
end

# various missing functions

function LLVMAddInternalizePassWithExportList(PM, ExportList, Length)
    ccall(:LLVMExtraAddInternalizePassWithExportList,Cvoid,(LLVMPassManagerRef,Ptr{Cstring},Csize_t), PM, ExportList, Length)
end

function LLVMAddTargetLibraryInfoByTriple(Triple, PM)
    ccall(:LLVMExtraAddTargetLibraryInfoByTiple,Cvoid,(Cstring, LLVMPassManagerRef), Triple, PM)
end

if VERSION < v"1.2.0-DEV.531"
function LLVMAddNVVMReflectPass(PM, smversion)
    ccall(:LLVMExtraAddMVVMReflectPass,Cvoid,(LLVMPassManagerRef,), PM)
end
else

if version() < v"8.0"
    function LLVMAddNVVMReflectPass(PM, smversion)
        ccall(:LLVMExtraAddNVVMReflectPass,Cvoid,(LLVMPassManagerRef,), PM)
    end
else
    function LLVMAddNVVMReflectPass(PM, smversion)
        ccall(:LLVMExtraAddNVVMReflectFunctionPass,Cvoid,(LLVMPassManagerRef, Cuint), PM, smversion)
    end
end

function LLVMAddAllocOptPass(PM)
    ccall(:LLVMExtraAddAllocOptPass,Cvoid,(LLVMPassManagerRef,), PM)
end

function LLVMAddBarrierNoopPass(PM)
    ccall(:LLVMExtraAddBarrierNoopPass,Cvoid,(LLVMPassManagerRef,), PM)
end

function LLVMAddGCInvariantVerifierPass(PM, Strong)
    ccall(:LLVMExtraAddGCInvariantVerifierPass,Cvoid,(LLVMPassManagerRef,LLVMBool), PM, Strong)
end

function LLVMAddLowerExcHandlersPass(PM)
    ccall(:LLVMExtraAddLowerExcHandlersPass,Cvoid,(LLVMPassManagerRef,), PM)
end

function LLVMAddCombineMulAddPass(PM)
    ccall(:LLVMExtraAddCombineMulAddPass,Cvoid,(LLVMPassManagerRef,), PM)
end

function LLVMAddMultiVersioningPass(PM)
    ccall(:LLVMExtraAddMultiVersioningPass,Cvoid,(LLVMPassManagerRef,), PM)
end

function LLVMAddPropagateJuliaAddrspaces(PM)
    ccall(:LLVMExtraAddPropagateJuliaAddrspaces,Cvoid,(LLVMPassManagerRef,), PM)
end

function LLVMAddLowerPTLSPass(PM, imaging_mode)
    ccall(:LLVMExtraAddLowerPTLSPass,Cvoid,(LLVMPassManagerRef,LLVMBool), PM, imaging_mode)
end

function LLVMAddLowerSimdLoopPass(PM)
    ccall(:LLVMExtraAddLowerSimdLoopPass,Cvoid,(LLVMPassManagerRef,), PM)
end

function LLVMAddLateLowerGCFramePass(PM)
    ccall(:LLVMExtraAddLateLowerGCFramePass,Cvoid,(LLVMPassManagerRef,), PM)
end

end

if VERSION >= v"1.3.0-DEV.95"
function LLVMAddFinalLowerGCPass(PM)
    ccall(:LLVMExtraAddFinalLowerGCPass,Cvoid,(LLVMPassManagerRef,), PM)
end
end

if VERSION >= v"1.5.0-DEV.802"
function LLVMAddRemoveJuliaAddrspacesPass(PM)
    ccall(:LLVMExtraAddRemoveJuliaAddrspacesPass,Cvoid,(LLVMPassManagerRef,), PM)
end
end

function LLVMGetValueContext(V)
    ccall(:LLVMExtraGetValueContext,LLVMContextRef,(LLVMValueRef,),V)
end

if VERSION >= v"0.7.0-alpha.37"
function LLVMGetSourceLocation(V, index, Name, Filename, Line, Column)
    ccall(:LLVMExtraGetSourceLocation,Cint,(LLVMValueRef,Cint,Ptr{Cstring},Ptr{Cstring},Ptr{Cuint},Ptr{Cuint}), V, index, Name, Filename, Line, Column)
end
end

if VERSION >= v"1.5" && !(v"1.6-" <= VERSION < v"1.6.0-DEV.90")
function LLVMExtraAppendToUsed(Mod, Values, Count)
    ccall(:LLVMExtraAppendToUsed,Cvoid,(LLVMModuleRef,Ptr{LLVMValueRef},Csize_t), Mod, Values, Count)
end

function LLVMExtraAppendToCompilerUsed(Mod, Values, Count)
    ccall(:LLVMExtraAppendToCompilerUsed,Cvoid,(LLVMModuleRef,Ptr{LLVMValueRef},Csize_t), Mod, Values, Count)
end

function LLVMExtraAddGenericAnalysisPasses(PM)
    ccall(:LLVMExtraAddGenericAnalysisPasses, Cvoid, (LLVMPassManagerRef,), PM)
end
end

if version() >= v"8.0"
    @cenum(LLVMDebugEmissionKind,
        LLVMDebugEmissionKindNoDebug = 0,
        LLVMDebugEmissionKindFullDebug = 1,
        LLVMDebugEmissionKindLineTablesOnly = 2,
        LLVMDebugEmissionKindDebugDirectivesOnly = 3,
    )
end


text:
using Documenter, RAFF

push!(LOAD_PATH, "../test/scripts/")

makedocs(
    assets = ["assets/favicon.ico"],
    sitename = "RAFF- Robust Algebraic Fitting Function",
    pages = ["Overview" => "index.md",
	     "Tutorial"=> "tutorial.md",
	     "Examples" => "examples.md",
             "API" => "api.md",
             "Advanced" => "advanced.md",
             "Random generation" => "randomgeneration.md"],
    #html_prettyurls = false
    #format = Documenter.HTML(prettyurls = false),
    modules = [RAFF]
	)

deploydocs(
	repo = "github.com/fsobral/RAFF.jl.git"
	)


text:
/* Fibonacci. */

int main () {
  int lo,hi,mx ;
  lo = 1 ;
  hi = lo ;
  mx = 5000000 ;
  printInt(lo) ;
  while (hi < mx) {
    printInt(hi) ;
    hi = lo + hi ;
    lo = hi - lo ;
  }
  return 0 ;

}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
# convenient way to get library information
using BinDeps

# export avcodec_dir, avformat_dir, avutil_dir, swscale_dir, avdevice_dir, avfilter_dir, avresample_dir, swresample_dir,
#        avcodec_version, avformat_version, avutil_version, swscale_version, avdevice_version, avfilter_version, avresample_version, swresample_version,
#        libavcodec, libavformat, libavutil, libswscale, libavdevice, libavfilter, libavresample, libswresample,

if !isdefined(:libavcodec)
    include("../deps/deps.jl")
end

INSTALL_ROOT = realpath(joinpath(splitdir(libavcodec)[1], ".."))

if !isdefined(:ffmpeg_or_libav)
    include("version.jl")
end

av_load_path = joinpath(dirname(@__FILE__), ffmpeg_or_libav)
!(av_load_path in LOAD_PATH) && unshift!(LOAD_PATH, av_load_path)

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
module RealDot

using LinearAlgebra: LinearAlgebra

export realdot

"""
    realdot(x, y)

Compute `real(dot(x, y))` while avoiding computing the imaginary part if possible.

This function can be useful if you work with derivatives of functions on complex
numbers. In particular, this computation shows up in pullbacks for non-holomorphic
functions.
"""
@inline realdot(x, y) = real(LinearAlgebra.dot(x, y))
@inline realdot(x::Complex, y::Complex) = muladd(real(x), real(y), imag(x) * imag(y))
@inline realdot(x::Real, y::Number) = x * real(y)
@inline realdot(x::Number, y::Real) = real(x) * y
@inline realdot(x::Real, y::Real) = x * y

end


text:
@testset "Standard synapse tests" begin
    test_net="""
spike_duration: 0.1
refractory_duration: 0.2
delay: 0.3
plateau_duration: 50.0

# Comments should work, tool
inputs:
    - id: i1 # even inline comments
    - id: i2
    - id: i3

outputs:
    - id: o1
    - id: o2

random_unsupported_network_stuff1: nothing_of_relevance
random_unsupported_network_stuff2: nothing_of_relevance

neurons:
    - id: n1
      random_unsupported_neuron_stuff: nothing_of_relevance
    - id: n2
      refractory_duration: 0.4
      θ_syn: 2
      θ_seg: 2
      branches:
        - id: s1
          θ_syn: 10
        - id: s2
          branches:
            - id: s21
              plateau_duration: 20.0
              random_unsupported_segment_stuff: nothing_of_relevance
              θ_seg: 0
            - id: s22

synapses:
    - {id: syn1, source: i1, target: n1}
    - {id: syn2, source: i2, target: s22}
    - {id: syn3, source: i3, target:  s2, spike_duration: 3.0, random_unsupported_synapse_stuff: nothing_of_relevance}
    - {id: syn4, source: n1, target:  o1, delay: 10.0, weight: 2}
    """

    (net,obj) = (
        @test_logs (
                :warn, "The following network parameters could not be parsed: random_unsupported_network_stuff1, random_unsupported_network_stuff2"
            ) (
                :warn, "The following neuron parameters could not be parsed: random_unsupported_neuron_stuff"
            ) (
                :warn, "The following segment parameters could not be parsed: random_unsupported_segment_stuff"
            ) (
                :warn, "The following synapse parameters could not be parsed: random_unsupported_synapse_stuff"
            ) load_network(YAML_source=test_net)
        )

    @testset "Network structure tests" begin
        @test net.neurons[1] == obj[:n1]
        @test net.neurons[2] == obj[:n2]
        @test net.neurons[2].next_upstream[1] == obj[:s1]
        @test net.neurons[2].next_upstream[2].next_upstream[1] == obj[:s21]
        @test net.inputs[1] == obj[:i1]
        @test net.inputs[1].synapses[1] == obj[:syn1]
        @test net.neurons[1].synapses[1] == obj[:syn4]
        @test net.outputs[1] == obj[:o1]
    end

    @testset "Neuron parameters tests" begin
        @test obj[:n1].id == :n1
        @test obj[:n1].refractory_duration == 0.2
        @test obj[:n1].θ_syn == 1
        @test obj[:n1].θ_seg == 1
        @test obj[:n2].refractory_duration == 0.4
        @test obj[:n2].θ_syn == 2
        @test obj[:n2].θ_seg == 2
    end

    @testset "Segment parameters tests" begin
        @test obj[:s1].id == :s1
        @test obj[:s1].plateau_duration == 50.0
        @test obj[:s1].θ_syn == 10
        @test obj[:s1].θ_seg == 1
        @test obj[:s21].plateau_duration == 20.0
        @test obj[:s21].θ_syn == 1
        @test obj[:s21].θ_seg == 0
    end

    @testset "Synapse parameters tests" begin
        @test obj[:syn3].id == :syn3
        @test obj[:syn3].target == obj[:s2]
        @test obj[:syn3].delay == 0.3
        @test obj[:syn3].spike_duration == 3.0
        @test obj[:syn3].weight == 1
        @test obj[:syn4].target == obj[:o1]
        @test obj[:syn4].delay == 10.0
        @test obj[:syn4].weight == 2
    end

    @testset "Input parameters tests" begin
        @test obj[:i1].id == :i1
    end

    @testset "Output parameters tests" begin
        @test obj[:o1].id == :o1
    end

    @testset "Saving tests" begin
        tmp_str = save_network(net)
        (net2,obj2) = @test_logs  load_network(YAML_source=tmp_str)
        @test all(pairs(obj2)) do (key,value)
            key ∈ keys(obj)
        end
    end
end

@testset "Bernoulli synapse tests" begin
    test_net="""
    spike_duration: 0.1
    delay: 0.3
    plateau_duration: 50.0

    inputs:
        - id: i1
        - id: i2

    neurons:
        - id: n1

    synapses:
        - {id: syn1, source: i1, target: n1}
        - {id: syn2, source: i2, target: n1, weight: 1.2345, spike_duration: 3.0, random_unsupported_synapse_stuff: nothing_of_relevance}
        - {id: syn3, source: n1, target: n1, delay: 10.0, weight: [1.2345, 0.75]}
    """

    (net,obj) = (
        @test_logs (
                :warn, "The following synapse parameters could not be parsed: random_unsupported_synapse_stuff"
            ) load_network(YAML_source=test_net, weight_type=BernoulliSynapseWeight{Float64})
        )

    @testset "Synapse parameters tests" begin
        @test obj[:syn1].id == :syn1
        @test obj[:syn1].target == obj[:n1]
        @test obj[:syn1].delay == 0.3
        @test obj[:syn1].spike_duration == 0.1
        @test isa(obj[:syn1].weight, BernoulliSynapseWeight{Float64}) && obj[:syn1].weight.magnitude ≈ 1.0 && obj[:syn1].weight.probability ≈ 1.0
        @test obj[:syn2].id == :syn2
        @test obj[:syn2].target == obj[:n1]
        @test obj[:syn2].delay == 0.3
        @test obj[:syn2].spike_duration == 3.0
        @test isa(obj[:syn2].weight, BernoulliSynapseWeight{Float64}) && obj[:syn2].weight.magnitude ≈ 1.2345 && obj[:syn2].weight.probability ≈ 1.0
        @test obj[:syn3].id == :syn3
        @test obj[:syn3].target == obj[:n1]
        @test obj[:syn3].delay == 10.0
        @test obj[:syn3].spike_duration == 0.1
        @test isa(obj[:syn3].weight, BernoulliSynapseWeight{Float64}) && obj[:syn3].weight.magnitude ≈ 1.2345 && obj[:syn3].weight.probability ≈ 0.75
    end

    @testset "Saving tests" begin
        tmp_str = save_network(net)
        (net2,obj2) = @test_logs  load_network(YAML_source=tmp_str, weight_type=BernoulliSynapseWeight{Float64})
        @test all(pairs(obj2)) do (key,value)
            key ∈ keys(obj)
        end
    end
end


text:
# Base category on Julia types

abstract type Arrow{A, B}; end

source(::Arrow{A,B}) where {A, B} = A
source(::Type{Arrow{A,B}}) where {A, B} = A
target(::Arrow{A,B}) where {A, B} = B
target(::Type{Arrow{A,B}}) where {A, B} = B

export source, target, @arrow, @alias, typestr

"Compute the source and target of the composed morphism g o f, or error if they are not composable."
function composition_obj(g::Arrow, f::Arrow)
    A = source(f)
    Bf = target(f)
    Bg = source(g)
    C = target(g)

    if Bf <: Bg
        return (A, C) # Composition will be an Arrow{A, C}
    else
        error("Cannot compose ($A~>$Bf) to ($Bg~>$C)")
    end
end

parent_category(a::Type{<:Arrow}) = __module__

function compose
end

function lift
end

function typestr(a::Arrow)
    "$(typeof(a).name) :: $(source(a)) ~> $(target(a))"
end

"
```
@arrow C Name{T1, T2} :: A --> B mutable struct
    p1::T1
    p2::T2
end
```

or

```
@arrow C Adder{T<:Number} :: T --> T
```
"
macro arrow(category, type_expr, structure=:())
    @capture(type_expr, name_ :: src_ --> tgt_)
    # TODO: Walk src, tgt and convert, e.g., (A, B) to Tuple{A, B}
    if isexpr(src, :tuple)
        src = :(Tuple{$(src.args...)})
    end
    has_type_params = @capture(name, sname_{targs__})
    new_expr = has_type_params ? :(new{$(targs...)}()) : :(new())
    where_expr = has_type_params ? :(where{$(targs...)}) : :()
    # TODO: Parse out structure fields and inner constructors
    #       then auto-composify inner constructors
    esc(quote
        mutable struct $name <: $category.Arrow{$src, $tgt}
           $(structure.args...)
           $name(args...) = compose($new_expr, args...)
        end
        end)
end

"Aliases an arrow constructor to a function"
macro alias(category, arrow_type, alias, arity)
    # TODO: automatically determine arity
    exprs = []
    args = convert(Vector{Any}, [Symbol("arg$k") for k in 1:arity])
    aargs = [:($a::$category.Arrow) for a in args]

    for q=1:(2^arity-1)
        pattern = q
        vargs = copy(args)
        k = 1
        while pattern > 0
            if pattern & 1 == 1
                vargs[k] = aargs[k]
            end
            pattern >>= 1
            k += 1
        end
        push!(exprs, :($alias($(vargs...)) = $arrow_type($(args...))))
    end

    esc(quote
            $(exprs...)
        end)
end

"Defines composition rules making m1 and m2 inverses (an isomorphism pair)."
macro inverse(category, m1, m2)
    # TODO: auto-infer category
    esc(quote
        function compose(g::$m1, f::$m2)
            A, B = composition_obj(g, f)
            @assert A == B
            $category.Identity{A}()
        end
        compose(g::$m2, f::$m1) = Cat.compose(f, g)

        function compose(g::$m1, f::$category.Composed{<:$m2})
            A, B = composition_obj(g, f.g)
            @assert A == B
            f.f
        end

        function compose(g::$m2, f::$category.Composed{<:$m1})
            A, B = composition_obj(g, f.g)
            @assert A == B
            f.f
        end
    end)
end

"Define a new category, which manifests itself as a module with some populated types."
macro category(arrow_type, flags...)
    product_code = include("$(@__DIR__)/product.jl")
    e =
        quote
            module $arrow_type
            export Arrow, Identity, Terminal, Composed, Constant
            using Cat
            abstract type Arrow{A, B} <: Cat.Arrow{A, B}; end
            struct Identity{A} <: Arrow{A, A}; end
            struct Terminal{A} <: Arrow{A, Nothing}; end

            Cat.parent_category(a::Type{<:Arrow}) = __module__

            Cat.compose(a::Arrow{A, B}, b::Identity{A}) where {A,B} = a
            Cat.compose(a::Identity{B}, b::Arrow{A, B}) where {A,B} = b

            "Free composition"
            struct Composed{T2, T1, S2, S1} <: Arrow{S1, S2}
                g::T2
                f::T1
            end

            $(product_code.args...)
            # Force associativity
            Cat.compose(b::Composed, a::Arrow) = compose(b.g, compose(b.f, a))

            # Handle tuples automatically
            function Cat.compose(g::T2, fs...) where {T2 <: Arrow}
                if isnothing(fs) || isempty(fs)
                    g
                else
                    fs = lift.(fs)
                    Cat.compose(g, Product(fs...))
                end
            end

            # Handle terminals
            function Cat.compose(g::Terminal{A}, f::Arrow{B,A}) where {A, B}
                Terminal{B}()
            end

            # General composition
            function Cat.compose(g::T2, f::T1) where {T1 <: Arrow, T2 <: Arrow}
                A, C = composition_obj(g, f)
                Composed{T2, T1, C, A}(g, f)
            end

            struct Constant{A} <: Arrow{Nothing, A}
                val::A
            end

            "Lifts values of type A to arrows Nothing~>A while leaving arrows alone"
            lift(a::Arrow) = a
            lift(a) = Constant(a)
           end
    end
    # This is necessary to allow defining modules
    esc(Expr(:toplevel, e.args...))
end