text:
using IterTools
using KernelDensity

function integrate_on_grid(f, ranges, density)
    rs = []
    for r in ranges
        step = (r[2] - r[1]) / density
        rmin = r[1] + step / 2
        rmax = r[2] - step / 2
        push!(rs, range(rmin, rmax, length=density))
    end

    logsum = -1e300
    for v in Iterators.product(rs...)
        logsum = log(exp(logsum) + f(v))
    end
    logsum -= length(ranges) * log(density)

    return logsum
end

function integrate_on_grid(f, ranges)
    density = 100
    logsum_old = -Inf
    while true
        logsum = integrate_on_grid(f, ranges, density)
        if abs(logsum - logsum_old) < 0.001
            return logsum
        end
        logsum_old = logsum
        density *= 2
    end
end

## Contrib from Firefly.jl
function findpeaks(samples::AbstractVector)
    k = kde(samples)
    # the sign of the difference will tell use whether we are increasing or decreasing
    # using rle gives us the points at which the sign switches (local extreema)
    runs = rle(sign.(diff(k.density)))
    # if we start going up, first extreme will be maximum, else minimum
    start = runs[1][1] == 1 ? 1 : 2
    # find the peak indices at the local minima
    peak_idx = cumsum(runs[2])[start:2:end]
    sorted_idx = sortperm(k.density[peak_idx], rev=true)
    return k.x[peak_idx[sorted_idx]]
end


text:
# return +-

function bias_adjusted_gaussian_ci(se; maxbias=0.0, level=0.9)
    rel_bias = maxbias/se
    zz = fzero( z-> cdf(Normal(), rel_bias-z) + cdf(Normal(), -rel_bias-z) +  level -1,
        0, rel_bias - quantile(Normal(),(1- level)/2.1))
    zz*se
end

ext: jl


text:
function StopOnEvent()
    return StopOnEvent(())
end

function event_occurred(obj::StopOnEvent, arg0::SpacecraftState, arg1::EventDetector, arg2::jboolean)
    return jcall(obj, "eventOccurred", Action, (SpacecraftState, EventDetector, jboolean), arg0, arg1, arg2)
end


text:
using PatModules
@auxmodule! two ("utils",) begin
    export say_bye

    function say_bye()
        say("bye")
    end
end


text:
using AstroTime
using Test
using ERFA
using SPICE: furnsh, et2utc, utc2et

furnsh(joinpath("data", "naif0012.tls"))
AstroTime.load_test_eop()

const speed_of_light = 299792458.0 # m/s
const astronomical_unit = 149597870700.0 # m

@timescale GMT UTC
AstroTime.Epochs.getoffset(::GMTScale, ::CoordinatedUniversalTime, _, _) = 0.0
AstroTime.Epochs.getoffset(::CoordinatedUniversalTime, ::GMTScale, _, _) = 0.0

@timescale SCET UTC
function AstroTime.Epochs.getoffset(::SCETScale, ::CoordinatedUniversalTime,
                                    _, _, distance)
    return distance / speed_of_light
end
function AstroTime.Epochs.getoffset(::CoordinatedUniversalTime, ::SCETScale,
                                    _, _, distance)
    return -distance / speed_of_light
end

@timescale Dummy TDB

@timescale Lonely
@timescale Together Lonely
function AstroTime.Epochs.getoffset(::LonelyScale, ::TogetherScale, _, _)
    return 5.0
end
function AstroTime.Epochs.getoffset(::TogetherScale, ::LonelyScale, _, _)
    return -5.0
end

@timescale OneWay TDB true

@testset "AstroTime" begin
    include("periods.jl")
    include("astrodates.jl")
    include("epochs.jl")
    @testset "Custom Time Scales" begin
        utc = UTCEpoch(2000, 1, 1, 0, 0, 0.1)
        gmt = GMTEpoch(2000, 1, 1, 0, 0, 0.1)
        @test utc == UTCEpoch(gmt)
        @test gmt == GMTEpoch(utc)

        @test string(GMT) == "GMT"
        @test typeof(GMT) == GMTScale
        @test string(gmt) == "2000-01-01T00:00:00.100 GMT"
        @test find_path(GMT, UTC) == [GMT, UTC]

        @test string(SCET) == "SCET"
        @test typeof(SCET) == SCETScale
        @test find_path(SCET, UTC) == [SCET, UTC]
        scet = SCETEpoch(2000, 1, 1, 0, 0, 0.1)
        utc_exp = UTCEpoch(2000, 1, 1, 0, 8, 19.10478383615643)
        @test UTCEpoch(scet, astronomical_unit) ≈ utc_exp
        @test SCETEpoch(utc_exp, astronomical_unit) ≈ scet

        @test string(Dummy) == "Dummy"
        @test typeof(Dummy) == DummyScale
        @test find_path(Dummy, UTC) == [Dummy, TDB, TT, TAI, UTC]
        dummy = DummyEpoch(2000, 1, 1)
        @test_throws NoOffsetError UTCEpoch(dummy)
        @test_throws NoOffsetError TDBEpoch(dummy)

        lonely = LonelyEpoch(2000, 1, 1)
        together = TogetherEpoch(lonely)
        @test together == TogetherEpoch(2000, 1, 1, 0, 0, 5.0)
        @test lonely == LonelyEpoch(together)
        @test find_path(Lonely, TDB) == []

        @test find_path(TDB, OneWay) == [TDB, OneWay]
        @test find_path(OneWay, TDB) == []
        @test_throws NoPathError TDBEpoch(OneWayEpoch(2000, 1, 1))
    end
end


text:
using NaturalES
using Test
using LinearAlgebra

f(x) = sum(cosh.(x.-(1:length(x)))) - length(x)
sphere(x) = sum(abs2,x)
function rosenbrock2d(x::AbstractVector{T}) where T
    s=(1.0 - x[1])^2
    for i in 1:(length(x)-1)
        s+=100.0 * (x[i+1] - x[i]^2)^2
    end
    return s
end

@testset "Optimizers" begin
    for method in [sNES,xNES], x in [0.5,1.0,4.0], fun in [f,sphere,rosenbrock2d]
        ix=[x,-x]
        @info "method: $method, testing init condition $ix, function: $fun"
        @test optimize(fun,ix,1.0,method).cost < 1e-7
        @test optimize(fun,ix,1.0,method,samples=10).cost < 1e-7
        @test optimize(fun,ix,2.0,method,ησ=0.01,atol=1e-8,iterations=10^4).cost < 1e-7
        @test optimize(fun,ix,2.0,method,ησ=0.01,atol=1e-8,iterations=10).cost > 1e-7
    end
end


text:
using DelimitedFiles
crabdata = readdlm("crabpos.txt",',',Int)
#@show crabdata

function fuelguage(crabdata)
    crabmax = maximum(crabdata)
    crabfuel = zeros(Int,1,crabmax)
    for crabfuelx in 1:crabmax
        for crabx in crabdata
            crabfuel[crabfuelx] += abs((crabfuelx-1) - crabx)
        end
    end

    #@show crabfuel
    crabfuel
end

@show minimum(fuelguage(crabdata))


text:
using Stheno: project, pairwise, lhsfinite, unbox, Xt_A_X, Xt_A_Y

@testset "project" begin

let
    rng, N, N′, NZ, NZ′ = MersenneTwister(123456), 10, 11, 5, 3
    X, X′, Z, Z′ = randn(rng, N), randn(rng, N′), randn(rng, NZ), randn(rng, NZ′)

    # Construct projection with non-block matrices.
    f = GP(ConstantMean(3.0), EQ(), GPC())
    ϕ = lhsfinite(EQ(), Z)
    f_pr = project(ϕ, f(Z), sin)

    @test mean_vec(f_pr(X)) == pairwise(EQ(), X, Z) * mean_vec(f(Z)) + sin.(X)
    @test cov(f_pr(X)) ≈ pairwise(EQ(), X, Z) * cov(f(Z)) * pairwise(EQ(), X, Z)'
    @test maximum(abs.(xcov(f_pr(X), f(X′)) - pairwise(EQ(), X, Z) * xcov(f(Z), f(X′)))) < eps()
    @test maximum(abs.(xcov(f(X), f_pr(X′)) - xcov(f(X), f(Z)) * pairwise(EQ(), Z, X′))) < eps()

    # Construct projection with BlockGP as the thing over which we're projecting.
    Zb = BlockData([Z, Z′])
    ϕb = lhsfinite(BlockCrossKernel(reshape([EQ(), EQ()], :, 1)), Zb)
    f_pr_b = project(ϕb, f(Zb), cos)

    @test mean_vec(f_pr_b(X)) == pairwise(EQ(), Zb, X)' * mean_vec(f(Zb)) .+ cos.(X)
    ϕZX, ϕZX′, ϕZ = pairwise(EQ(), Zb, X), pairwise(EQ(), Zb, X′), cov(f(Zb))
    @test cov(f_pr_b(X)) == Xt_A_X(ϕZ, ϕZX)
    @test xcov(f_pr_b(X), f_pr_b(X′)) == Xt_A_Y(ϕZX, ϕZ, ϕZX′)
    @test xcov(f_pr_b(X), f(X′)) == ϕZX' * xcov(f(Zb), f(X′))
    @test xcov(f(X), f_pr_b(X′)) == xcov(f(X), f(Zb)) * ϕZX′

    # Check that Blocked version is project is consistent with dense.
    Zb2 = BlockData([Z[1:3], Z[4:5]])
    ϕb2 = lhsfinite(BlockCrossKernel(reshape([EQ(), EQ()], :, 1)), Zb2)
    f_pr_b2 = project(ϕb2, f(Zb2), sin)

    @test mean_vec(f_pr_b2(X)) ≈ mean_vec(f_pr(X))
    @test cov(f_pr_b2(X)) ≈ cov(f_pr(X))
    @test xcov(f_pr_b2(X), f_pr_b2(X′)) ≈ xcov(f_pr(X), f_pr(X′))
    @test xcov(f_pr_b2(X), f(X′)) ≈ xcov(f_pr(X), f(X′))
    @test xcov(f(X), f_pr_b2(X′)) ≈ xcov(f(X), f_pr(X′))
end

end


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

#####################
# structs/constants #
#####################

# The type of a value might be constant
struct Const
    val
    actual::Bool  # if true, we obtained `val` by actually calling a @pure function
    Const(@nospecialize(v)) = new(v, false)
    Const(@nospecialize(v), a::Bool) = new(v, a)
end

# The type of this value might be Bool.
# However, to enable a limited amount of back-propagagation,
# we also keep some information about how this Bool value was created.
# In particular, if you branch on this value, then may assume that in
# the true branch, the type of `var` will be limited by `vtype` and in
# the false branch, it will be limited by `elsetype`. Example:
# ```
# cond = isa(x::Union{Int, Float}, Int)::Conditional(x, Int, Float)
# if cond
#    # May assume x is `Int` now
# else
#    # May assume x is `Float` now
# end
# ```
struct Conditional
    var::Slot
    vtype
    elsetype
    function Conditional(
                var,
                @nospecialize(vtype),
                @nospecialize(nottype))
        return new(var, vtype, nottype)
    end
end

struct PartialTypeVar
    tv::TypeVar
    # N.B.: Currently unused, but would allow turning something back
    # into Const, if the bounds are pulled out of this TypeVar
    lb_certain::Bool
    ub_certain::Bool
    PartialTypeVar(tv::TypeVar, lb_certain::Bool, ub_certain::Bool) = new(tv, lb_certain, ub_certain)
end

# Wraps a type and represents that the value may also be undef at this point.
struct MaybeUndef
    typ
end

# The type of a variable load is either a value or an UndefVarError
struct VarState
    typ
    undef::Bool
    VarState(@nospecialize(typ), undef::Bool) = new(typ, undef)
end

const VarTable = Array{Any,1}

struct StateUpdate
    var::Union{Slot,SSAValue}
    vtype::VarState
    state::VarTable
end

struct NotFound end

const NOT_FOUND = NotFound()

#####################
# lattice utilities #
#####################

function rewrap(@nospecialize(t), @nospecialize(u))
    isa(t, Const) && return t
    isa(t, Conditional) && return t
    return rewrap_unionall(t, u)
end

_typename(a) = Union{}
_typename(a::Vararg) = Any
_typename(a::TypeVar) = Any
function _typename(a::Union)
    ta = _typename(a.a)
    tb = _typename(a.b)
    ta === tb ? tb : (ta === Any || tb === Any) ? Any : Union{}
end
_typename(union::UnionAll) = _typename(union.body)

_typename(a::DataType) = Const(a.name)

# N.B.: typename maps type equivalence classes to a single value
typename_static(@nospecialize(t)) = isType(t) ? _typename(t.parameters[1]) : Any
typename_static(t::Const) = _typename(t.val)

#################
# lattice logic #
#################

function issubconditional(a::Conditional, b::Conditional)
    avar = a.var
    bvar = b.var
    if (isa(avar, Slot) && isa(bvar, Slot) && slot_id(avar) === slot_id(bvar)) ||
       (isa(avar, SSAValue) && isa(bvar, SSAValue) && avar === bvar)
        if a.vtype ⊑ b.vtype
            if a.elsetype ⊑ b.elsetype
                return true
            end
        end
    end
    return false
end

maybe_extract_const_bool(c::Const) = isa(c.val, Bool) ? c.val : nothing
function maybe_extract_const_bool(c::Conditional)
    (c.vtype === Bottom && !(c.elsetype === Bottom)) && return false
    (c.elsetype === Bottom && !(c.vtype === Bottom)) && return true
    nothing
end
maybe_extract_const_bool(c) = nothing

function ⊑(@nospecialize(a), @nospecialize(b))
    if isa(a, MaybeUndef) && !isa(b, MaybeUndef)
        return false
    end
    isa(a, MaybeUndef) && (a = a.typ)
    isa(b, MaybeUndef) && (b = b.typ)
    (a === NOT_FOUND || b === Any) && return true
    (a === Any || b === NOT_FOUND) && return false
    a === Union{} && return true
    b === Union{} && return false
    if isa(a, Conditional)
        if isa(b, Conditional)
            return issubconditional(a, b)
        elseif isa(b, Const) && isa(b.val, Bool)
            return maybe_extract_const_bool(a) === b.val
        end
        a = Bool
    elseif isa(b, Conditional)
        return a === Bottom
    end
    if isa(a, Const)
        if isa(b, Const)
            return a.val === b.val
        end
        return isa(a.val, widenconst(b))
    elseif isa(b, Const)
        return a === Bottom
    elseif !(isa(a, Type) || isa(a, TypeVar)) ||
           !(isa(b, Type) || isa(b, TypeVar))
        return a === b
    else
        return a <: b
    end
end

widenconst(c::Conditional) = Bool
function widenconst(c::Const)
    if isa(c.val, Type)
        if isvarargtype(c.val)
            return Type
        end
        return Type{c.val}
    else
        return typeof(c.val)
    end
end
widenconst(m::MaybeUndef) = widenconst(m.typ)
widenconst(c::PartialTypeVar) = TypeVar
widenconst(@nospecialize(t)) = t

issubstate(a::VarState, b::VarState) = (a.typ ⊑ b.typ && a.undef <= b.undef)

function tmerge(@nospecialize(typea), @nospecialize(typeb))
    typea ⊑ typeb && return typeb
    typeb ⊑ typea && return typea
    if isa(typea, MaybeUndef) || isa(typeb, MaybeUndef)
        return MaybeUndef(tmerge(
            isa(typea, MaybeUndef) ? typea.typ : typea,
            isa(typeb, MaybeUndef) ? typeb.typ : typeb))
    end
    if isa(typea, Conditional) && isa(typeb, Conditional)
        if typea.var === typeb.var
            vtype = tmerge(typea.vtype, typeb.vtype)
            elsetype = tmerge(typea.elsetype, typeb.elsetype)
            if vtype != elsetype
                return Conditional(typea.var, vtype, elsetype)
            end
        end
        return Bool
    end
    typea, typeb = widenconst(typea), widenconst(typeb)
    typea === typeb && return typea
    if !(isa(typea,Type) || isa(typea,TypeVar)) || !(isa(typeb,Type) || isa(typeb,TypeVar))
        return Any
    end
    if (typea <: Tuple) && (typeb <: Tuple)
        if isa(typea, DataType) && isa(typeb, DataType) && length(typea.parameters) == length(typeb.parameters) && !isvatuple(typea) && !isvatuple(typeb)
            return typejoin(typea, typeb)
        end
        if isa(typea, Union) || isa(typeb, Union) || (isa(typea,DataType) && length(typea.parameters)>3) ||
            (isa(typeb,DataType) && length(typeb.parameters)>3)
            # widen tuples faster (see #6704), but not too much, to make sure we can infer
            # e.g. (t::Union{Tuple{Bool},Tuple{Bool,Int}})[1]
            return Tuple
        end
    end
    u = Union{typea, typeb}
    if unionlen(u) > MAX_TYPEUNION_LEN || type_too_complex(u, MAX_TYPE_DEPTH)
        # don't let type unions get too big
        # TODO: something smarter, like a common supertype
        return Any
    end
    return u
end

function smerge(sa::Union{NotFound,VarState}, sb::Union{NotFound,VarState})
    sa === sb && return sa
    sa === NOT_FOUND && return sb
    sb === NOT_FOUND && return sa
    issubstate(sa, sb) && return sb
    issubstate(sb, sa) && return sa
    return VarState(tmerge(sa.typ, sb.typ), sa.undef | sb.undef)
end

@inline tchanged(@nospecialize(n), @nospecialize(o)) = o === NOT_FOUND || (n !== NOT_FOUND && !(n ⊑ o))
@inline schanged(@nospecialize(n), @nospecialize(o)) = (n !== o) && (o === NOT_FOUND || (n !== NOT_FOUND && !issubstate(n, o)))

function widenconditional(typ::Conditional)
    if typ.vtype == Union{}
        return Const(false)
    elseif typ.elsetype == Union{}
        return Const(true)
    else
        return Bool
    end
end

function stupdate!(state::Tuple{}, changes::StateUpdate)
    newst = copy(changes.state)
    if isa(changes.var, Slot)
        changeid = slot_id(changes.var::Slot)
        newst[changeid] = changes.vtype
        # remove any Conditional for this Slot from the vtable
        for i = 1:length(newst)
            newtype = newst[i]
            if isa(newtype, VarState)
                newtypetyp = newtype.typ
                if isa(newtypetyp, Conditional) && slot_id(newtypetyp.var) == changeid
                    newst[i] = VarState(widenconditional(newtypetyp), newtype.undef)
                end
            end
        end
    end
    return newst
end

function stupdate!(state::VarTable, changes::StateUpdate)
    if !isa(changes.var, Slot)
        return stupdate!(state, changes.state)
    end
    newstate = false
    changeid = slot_id(changes.var::Slot)
    for i = 1:length(state)
        if i == changeid
            newtype = changes.vtype
        else
            newtype = changes.state[i]
        end
        oldtype = state[i]
        # remove any Conditional for this Slot from the vtable
        if isa(newtype, VarState)
            newtypetyp = newtype.typ
            if isa(newtypetyp, Conditional) && slot_id(newtypetyp.var) == changeid
                newtype = VarState(widenconditional(newtypetyp), newtype.undef)
            end
        end
        if schanged(newtype, oldtype)
            newstate = state
            state[i] = smerge(oldtype, newtype)
        end
    end
    return newstate
end

function stupdate!(state::VarTable, changes::VarTable)
    newstate = false
    for i = 1:length(state)
        newtype = changes[i]
        oldtype = state[i]
        if schanged(newtype, oldtype)
            newstate = state
            state[i] = smerge(oldtype, newtype)
        end
    end
    return newstate
end

stupdate!(state::Tuple{}, changes::VarTable) = copy(changes)

stupdate!(state::Tuple{}, changes::Tuple{}) = false

function stupdate1!(state::VarTable, change::StateUpdate)
    if !isa(change.var, Slot)
        return false
    end
    changeid = slot_id(change.var::Slot)
    # remove any Conditional for this Slot from the catch block vtable
    for i = 1:length(state)
        oldtype = state[i]
        if isa(oldtype, VarState)
            oldtypetyp = oldtype.typ
            if isa(oldtypetyp, Conditional) && slot_id(oldtypetyp.var) == changeid
                state[i] = VarState(widenconditional(oldtypetyp), oldtype.undef)
            end
        end
    end
    # and update the type of it
    newtype = change.vtype
    oldtype = state[changeid]
    if schanged(newtype, oldtype)
        state[changeid] = smerge(oldtype, newtype)
        return true
    end
    return false
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
using BenchmarkTools, FinancialMonteCarlo, ForwardDiff

S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 100000;
Nstep = 30;
sigma = 0.2
mc = MonteCarloConfiguration(Nsim, Nstep);
toll = 1e-3;

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));

f__1(x::Vector) = pricer(BlackScholesProcess(x[1], Underlying(x[2], x[4])), ZeroRate(x[3]), mc, EuropeanOption(x[5], K));
x = Float64[sigma, S0, r, d, T]
g__1 = x -> ForwardDiff.gradient(f__1, x);
y0 = g__1(x);
@btime f__1(x);
@btime g__1(x);

f_(x::Vector) = pricer(BlackScholesProcess(x[1], Underlying(x[2], x[4])), ZeroRate(x[3]), mc, AmericanOption(x[5], K));
g_ = x -> ForwardDiff.gradient(f_, x);
y0 = g_(x);
@btime f_(x);
@btime g_(x);


text:
# issue 462

using Gadfly

if VERSION < v"0.4-dev"
    using Dates
else
    using Base.Dates
end

a = [unix2datetime(100)]
b = [10]

plot(x=a, y=b, Geom.point)


text:
using RoadRunner
#println("in the callAPI file: ", rrlib)

RoadRunner.setConfigInt("LOADSBMLOPTIONS_CONSERVED_MOIETIES", 1)
#rrlib = Libdl.dlopen("C:/vs_rebuild/install/roadrunner/bin/roadrunner_c_api.dll")
ant_str = """
    const Xo, X1
     Xo -> S1; k1*Xo - k2*S1
     S1 -> S2; k3*S2
     S2 -> X1; k4*S2

     Xo = 1;   X1 = 0
     S1 = 0;   S2 = 0
     k1 = 0.1; k2 = 0.56
     k3 = 1.2; k4 = 0.9
"""
rr = loada(ant_str)
println("conserved moieties: ", RoadRunner.getConfigInt("LOADSBMLOPTIONS_CONSERVED_MOIETIES"))
data = simulate(rr, 0, 50, 51)
println(data)
ss = steadyState(rr)
println("this is the steady state value: ", ss)
ssValues = RoadRunner.computeSteadyStateValues(rr)
# str = RoadRunner.getListOfConfigKeys()
ids = getFloatingSpeciesIds(rr)
println(ids)

sbmlFile = "C:/Users/lukez/OneDrive/Desktop/Network Generator/Network-Generator/minimal.xml"
f = open(sbmlFile)
sbmlStr = read(f,String)
close(f)

rr = createRRInstance()         # Start up roadRunner
RoadRunner.loadSBML(rr, sbmlStr)
RoadRunner.addCompartment(rr, "compartment", 1.0, false)
RoadRunner.addSpecies(rr, "S1", "compartment", 0.1, "concentration", false)
RoadRunner.addSpecies(rr, "S2", "compartment", 0.1, "concentration", false)
RoadRunner.addParameter(rr, "k1", 0.1, false)
RoadRunner.addSpecies(rr, "inter-S1_S2", "compartment", 0.0, "concentration", false)
RoadRunner.addReaction(rr, "J1", ["S1"], ["S2"], "k1 * S1", true)
println("Here are the floating species: ", RoadRunner.getFloatingSpeciesIds(rr))


text:
"""
Check if should skip data point i, for skip arg which may be `nothing`.
"""
shouldskip(skip, i) = !isnothing(skip) && skip(i)

# Filter array by predicate on indices
skipindices!(skip, a) = deleteat!(a, findall(skip, axes(a, 1)))
skipindices!(skip::Nothing, a) = a

skipindices(skip, a) = [x for (i, x) in enumerate(a) if !skip(i)]
skipindices(skip::Nothing, a) = a

# Filter array by predicate on values
skipvalues(skip, a) = filter(i -> !skip(i), a)
skipvalues(skip::Nothing, a) = a

skipvalues!(skip, a) = filter!(i -> !skip(i), a)
skipvalues!(skip::Nothing, a) = a

"""
Hamming distance between two bit strings (as integers)
"""
function hamming(a::Integer, b::Integer)
    count_ones(xor(a, b))
end

"""
Euclidean distance between two vectors.
"""
euclidean(x, y) = sqrt(sum((x .- y).^2))

"""
Brute-force implementation of find().
"""
function find_bruteforce(data, metric, query, radius, skip=nothing)
    return [i for (i, d) in enumerate(data) if !shouldskip(skip, i) && metric(query, d) <= radius]
end

"""
Brute-force implementation of find().

Note - correct answer is only unambiguous if all distances are distinct
(otherwise there may be a tie in choosing k-th closest data point)."""
function find_nearest_bruteforce(data, metric, query, k, skip=nothing)
    dists = [metric(query, d) for d in data]
    idxs = skipindices!(skip, sortperm(dists))
    return idxs[1:min(k, end)]
end

"""
Brute-force implementation of find(), but return ties for k-th nearest-neighbor
as separate set of results.
"""
function find_nearest_bruteforce_ties(data, metric, query, k, skip=nothing)
    dists = [metric(query, d) for d in data]
    idxs = skipvalues!(skip, sortperm(dists))
    k = min(k, length(idxs))

    maxd = dists[idxs[k]]

    kmin = k
    while kmin > 1 && dists[idxs[kmin - 1]] == maxd
        kmin -= 1
    end

    kmax = k
    while kmax < length(idxs) && dists[idxs[kmax + 1]] == maxd
        kmax += 1
    end

    return idxs[1:(kmin-1)], idxs[kmin:kmax]
end

"""
Test that find() gives identical results to find_bruteforce().
"""
function test_find(vptree, query, radius, skip=nothing)
    result = find(vptree, query, radius, skip)
    bfresult = find_bruteforce(vptree.data, vptree.metric, query, radius, skip)
    @test issetequal(result, bfresult)
end

"""
Test that find_nearest() gives identical results to find_nearest_bruteforce().

This accounts for cases in which there is a tie for the k-th closest point.
"""
function test_find_nearest(vptree, query, k, skip=nothing)
    result = find_nearest(vptree, query, k, skip)

    nearest, ties = find_nearest_bruteforce_ties(vptree.data, vptree.metric, query, k, skip)

    @test length(result) == min(k, length(vptree.data))
    @test nearest ⊆ result
    @test setdiff(result, nearest) ⊆ ties
end

# Skip every 3rd data point
skip3(i) = i % 3 == 0

@testset "Tie for kth-nearest-neighbor" begin
    Random.seed!(1)

    metric = hamming
    query = 0b0001
    data = 0b0000:0b1111
    data = [
        # 4 @ D=1
        0b0000,
        0b0011,
        0b0101,
        0b1001,
        # 6 @ D=2
        0b0010,
        0b0100,
        0b1000,
        0b1101,
        0b1011,
        0b0111,
        # 4 @ D=3
        0b1100,
        0b1010,
        0b0110,
        0b1111,
        # 1 @ D=4
        0b1110,
    ]
    vptree = VPTree(data, metric)

    # KNN can include any two data points @ D=2
    k = 6

    nearest, ties = find_nearest_bruteforce_ties(data, metric, query, k)
    @test issetequal(nearest, 1:4)
    @test issetequal(ties, 5:10)

    test_find_nearest(vptree, query, k)

    # With skip arg
    # This skips 1 of the D=1 and two of the D=2, so D=2 is still the tied value
    # with k=6.
    # So, sets should be the same as before just with skipped indices filtered out.
    skip = skip3
    nearest_s, ties_s = find_nearest_bruteforce_ties(data, metric, query, k, skip)
    @test issetequal(nearest_s, skipvalues(skip, nearest))
    # @test issetequal(ties_s, skipindices(skip, ties))
    @test Set(ties_s) == Set(skipvalues(skip, ties))
    test_find_nearest(vptree, query, k, skip)
end

@testset "Hamming distance" begin
    Random.seed!(1)

    metric = hamming
    data = collect(0x0000:0xFFFF)

    queries = [0x0000, 0x0001, 0x4eaf, 0xa44a, 0xFFFF]
    k = 50
    r = 5

    for threaded in [false, true]
        vptree = VPTree(data, metric, threaded=threaded)

        for query in queries
            for skip in [nothing, skip3]
                test_find(vptree, query, r, skip)
                test_find_nearest(vptree, query, k, skip)
            end
        end
    end
end

@testset "Euclidean distance" begin
    Random.seed!(1)

    metric = euclidean
    data = [rand(2) for _ in 1:10000]

    queries = [rand(2) for _ in 1:10]
    k = 50
    r = .25

    for threaded in [false, true]
        vptree = VPTree(data, metric)

        for query in queries
            for skip in [nothing, skip3]
                test_find(vptree, query, r, skip)
                test_find_nearest(vptree, query, k, skip)
            end
        end
    end
end

@testset "Levenshtein distance" begin
    Random.seed!(1)

    queries = ["bla", "blub", "asdf", ":assd", "ast", "baube"]
    alphabet = union(queries..., "xyz")

    metric = Levenshtein()
    data = unique(randstring(alphabet, rand(3:5)) for _ in 1:10000)

    k = 10
    r = 3

    for threaded in [false, true]
        vptree = VPTree(data, metric)

        for query in queries
            for skip in [nothing, skip3]
                test_find(vptree, query, r, skip)
                test_find_nearest(vptree, query, k, skip)
            end
        end
    end
end

@testset "Construct invalid" begin

    # Empty data set
    @test_throws AssertionError VPTree(Vector{Float64}[], euclidean)

    # Non-callable
    @test_throws AssertionError VPTree([1.], 1)

    # No method for metric with data type
    @test_throws MethodError VPTree(["foo"], euclidean)
    @test_throws MethodError VPTree(["foo"], hamming)
end


text:
module AdvancedHMC

const DEBUG = Bool(parse(Int, get(ENV, "DEBUG_AHMC", "0")))

using Statistics: mean, var, middle
using LinearAlgebra: Symmetric, UpperTriangular, mul!, ldiv!, dot, I, diag, cholesky
using StatsFuns: logaddexp
using LazyArrays: BroadcastArray
using Random: GLOBAL_RNG, AbstractRNG
using ProgressMeter
using Parameters: @unpack, reconstruct
using ArgCheck: @argcheck

import StatsBase: sample

# Support of passing a vector of RNGs
Base.rand(rng::AbstractVector{<:AbstractRNG}) = rand.(rng)
Base.randn(rng::AbstractVector{<:AbstractRNG}) = randn.(rng)
function Base.randn(
    rng::AbstractVector{<:AbstractRNG},
    T,
    sz::Int...)
    return cat(randn.(rng, T, first(sz))...; dims=2)
end

# Notations
# ℓπ: log density of the target distribution
# ∂ℓπ∂θ: gradient of the log density of the target distribution
# θ: position variables / model parameters
# r: momentum variables

include("adaptation/Adaptation.jl")
export UnitEuclideanMetric, DiagEuclideanMetric, DenseEuclideanMetric
export NesterovDualAveraging, Preconditioner, NaiveHMCAdaptor, StanHMCAdaptor
using .Adaptation

using .Adaptation: AbstractScalarOrVec
import .Adaptation: adapt!, NesterovDualAveraging

include("hamiltonian.jl")
export Hamiltonian
include("integrator.jl")
export Leapfrog, JitteredLeapfrog, TemperedLeapfrog
include("trajectory.jl")
export StaticTrajectory, HMCDA, NUTS, EndPointTS, SliceTS, MultinomialTS, ClassicNoUTurn, GeneralisedNoUTurn, find_good_eps

include("diagnosis.jl")
include("sampler.jl")
export sample

### Init

using Requires

function __init__()
    include(joinpath(@__DIR__, "contrib/diffeq.jl"))
    include(joinpath(@__DIR__, "contrib/ad.jl"))
end

end # module


text:
using Pkg

function main()
    deps = [
        "Random",
        "Statistics",
        "Distributions",
        "LinearAlgebra",
        "InvertedIndices",
        "IntervalSets",
        "Parameters",
        "Roots",
        "Reexport",
    ]

    for pkg in deps
        Pkg.add(pkg)
    end
end


text:
#Julia Regular Expressions

ismatch(r"^\([0-9]{3}\)[0-9]{3}-[0-9]{4}$", "(781)244-1212")
ismatch(r"^\([0-9]{3}\)[0-9]{3}-[0-9]{4}$", "-781-244-1212")


text:
# TODO: random transposition; once we have this, monte carlo over
#  permutations is easy to implement


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ProtoBuf.google.protobuf

mutable struct SessionStartInfo_HparamsEntry <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SessionStartInfo_HparamsEntry(; kwargs...)
        obj = new(meta(SessionStartInfo_HparamsEntry), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SessionStartInfo_HparamsEntry (mapentry)
const __meta_SessionStartInfo_HparamsEntry = Ref{ProtoMeta}()
function meta(::Type{SessionStartInfo_HparamsEntry})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SessionStartInfo_HparamsEntry)
            __meta_SessionStartInfo_HparamsEntry[] = target = ProtoMeta(SessionStartInfo_HparamsEntry)
            allflds = Pair{Symbol,Union{Type,String}}[:key => AbstractString, :value => ProtoBuf.google.protobuf.Value]
            meta(target, SessionStartInfo_HparamsEntry, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SessionStartInfo_HparamsEntry[]
    end
end
function Base.getproperty(obj::SessionStartInfo_HparamsEntry, name::Symbol)
    if name === :key
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :value
        return (obj.__protobuf_jl_internal_values[name])::ProtoBuf.google.protobuf.Value
    else
        getfield(obj, name)
    end
end

mutable struct SessionStartInfo <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SessionStartInfo(; kwargs...)
        obj = new(meta(SessionStartInfo), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SessionStartInfo
const __meta_SessionStartInfo = Ref{ProtoMeta}()
function meta(::Type{SessionStartInfo})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SessionStartInfo)
            __meta_SessionStartInfo[] = target = ProtoMeta(SessionStartInfo)
            allflds = Pair{Symbol,Union{Type,String}}[:hparams => Base.Dict{AbstractString,ProtoBuf.google.protobuf.Value}, :model_uri => AbstractString, :monitor_url => AbstractString, :group_name => AbstractString, :start_time_secs => Float64]
            meta(target, SessionStartInfo, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SessionStartInfo[]
    end
end
function Base.getproperty(obj::SessionStartInfo, name::Symbol)
    if name === :hparams
        return (obj.__protobuf_jl_internal_values[name])::Base.Dict{AbstractString,ProtoBuf.google.protobuf.Value}
    elseif name === :model_uri
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :monitor_url
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :group_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :start_time_secs
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct SessionEndInfo <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SessionEndInfo(; kwargs...)
        obj = new(meta(SessionEndInfo), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SessionEndInfo
const __meta_SessionEndInfo = Ref{ProtoMeta}()
function meta(::Type{SessionEndInfo})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SessionEndInfo)
            __meta_SessionEndInfo[] = target = ProtoMeta(SessionEndInfo)
            allflds = Pair{Symbol,Union{Type,String}}[:status => Int32, :end_time_secs => Float64]
            meta(target, SessionEndInfo, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SessionEndInfo[]
    end
end
function Base.getproperty(obj::SessionEndInfo, name::Symbol)
    if name === :status
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :end_time_secs
        return (obj.__protobuf_jl_internal_values[name])::Float64
    else
        getfield(obj, name)
    end
end

mutable struct HParamsPluginData <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function HParamsPluginData(; kwargs...)
        obj = new(meta(HParamsPluginData), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct HParamsPluginData
const __meta_HParamsPluginData = Ref{ProtoMeta}()
function meta(::Type{HParamsPluginData})
    ProtoBuf.metalock() do
        if !isassigned(__meta_HParamsPluginData)
            __meta_HParamsPluginData[] = target = ProtoMeta(HParamsPluginData)
            allflds = Pair{Symbol,Union{Type,String}}[:version => Int32, :experiment => Experiment, :session_start_info => SessionStartInfo, :session_end_info => SessionEndInfo]
            oneofs = Int[0,1,1,1]
            oneof_names = Symbol[Symbol("data")]
            meta(target, HParamsPluginData, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, oneofs, oneof_names)
        end
        __meta_HParamsPluginData[]
    end
end
function Base.getproperty(obj::HParamsPluginData, name::Symbol)
    if name === :version
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :experiment
        return (obj.__protobuf_jl_internal_values[name])::Experiment
    elseif name === :session_start_info
        return (obj.__protobuf_jl_internal_values[name])::SessionStartInfo
    elseif name === :session_end_info
        return (obj.__protobuf_jl_internal_values[name])::SessionEndInfo
    else
        getfield(obj, name)
    end
end

export HParamsPluginData, SessionStartInfo_HparamsEntry, SessionStartInfo, SessionEndInfo
# mapentries: "SessionStartInfo_HparamsEntry" => ("AbstractString", "ProtoBuf.google.protobuf.Value")


text:
# strings of parameter names that are not fixed
function get_parameter_names(samples::DensitySampleVector)
    parameter_names = collect(string.(keys(BAT.varshape(samples))))
    parameter_names = [parameter_names[i] for i in 1:length(parameter_names) if BAT.varshape(samples)._accessors[i].len != 0]
    return parameter_names
end

@recipe function f(
    maybe_shaped_samples::DensitySampleVector,
    p::Float64;
    parameter_names = get_parameter_names(maybe_shaped_samples),
    y_positions = collect(1:length(parameter_names))*-1,
    y_offset = 0,
    bins = 200,
    atol = 0,
)
    samples = BAT.unshaped.(maybe_shaped_samples)

    for i in 1:length(y_positions)

        lower, upper =  get_smallest_interval_edges(samples, i, p, bins=bins, atol=atol)
        w = 0.5*(upper-lower)
        middle = lower .+ w

        ypos = fill(y_positions[i], length(middle)).-y_offset

        if i == 1
            intern_label = get(plotattributes, :label, "")
        else
            intern_label = ""
        end

        internal_color = get(plotattributes, :markerstrokecolor , :steelblue)

        @series begin
            seriestype --> :path
            markerstrokealpha := 0
            linewidth := 0
            markerstrokewidth --> 1
            markersize --> 5
            markerstrokecolor --> :steelblue
            linecolor := internal_color
            xerror := w
            yticks --> (y_positions, parameter_names)
            label := intern_label
            xguide --> "Ci"

            (middle, ypos)
        end
    end
end


text:
import FilePathsBase: /
using Base.Libc

"""
Has a field fn which takes stream as it's first argument

 fn(stream, ...)
"""
abstract type AbstractHandler end

(handler::AbstractHandler)(args...) = handler.fn(args...)

struct Middleware  <: AbstractHandler
    fn
end

struct HttpHandler  <: AbstractHandler
	fn
end

struct Folder{T <: AbstractPath}  <: AbstractHandler
    path::T
    fn

    function Folder(path::T) where T <: AbstractPath

         function fn(io, file; strict=true)

            if file isa AbstractString
                file = T(file)
            end

            # collecting and then splating handles joining paths starting with /
            file = normalize(joinpath(path, collect(file)...))
            hasprefix = startswith(string(path))
            file_str = string(file)

            if strict && !hasprefix(file_str)
                error("File above folder")
            end

            try

                if !isfile(file)
                    @warn "File not found $file"
                        # the http server errors if we write a empty string
                    Bonsai.write(io, "File not found", ResponseCodes.NotFound())
                    return
                end

                body = Base.read(file)

                # disable caching util we can implement it more robustly

                # body = get!(folder.lru, file_str) do
                #     read(file)
                # end

                Base.write(io, body)

                if io isa HTTP.Stream
                    HTTP.setheader(io, "Content-Type" => mime_type(file))
                end
            catch e
                # I'm unsure if the error codes are the same on windows
                # https://www.thegeekstuff.com/2010/10/linux-error-codes/
                if e isa SystemError && e.errnum == Libc.ENOENT && io isa HTTP.Stream
                    @warn e
                else
                    rethrow(e)
                end
            end
        end

        new{T}(
            normalize(path),
            fn
        )

    end
end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
# Original implementations by @nantonel
# https://github.com/maxmouchet/HMMBase.jl/pull/6

# In-place forward pass, where α and c are allocated beforehand.
function forwardlog!(
    α::AbstractMatrix,
    c::AbstractVector,
    a::AbstractVector,
    A::AbstractMatrix,
    LL::AbstractMatrix,
)
    @argcheck size(α, 1) == size(LL, 1) == size(c, 1)
    @argcheck size(α, 2) == size(LL, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    T, K = size(LL)
    (T == 0) && return

    fill!(α, 0.0)
    fill!(c, 0.0)

    m = vec_maximum(view(LL, 1, :))

    for j in OneTo(K)
        α[1, j] = a[j] * exp(LL[1, j] - m)
        c[1] += α[1, j]
    end

    for j in OneTo(K)
        α[1, j] /= c[1]
    end

    c[1] = log(c[1]) + m

    @inbounds for t = 2:T
        m = vec_maximum(view(LL, t, :))

        for j in OneTo(K)
            for i in OneTo(K)
                α[t, j] += α[t-1, i] * A[i, j]
            end
            α[t, j] *= exp(LL[t, j] - m)
            c[t] += α[t, j]
        end

        for j in OneTo(K)
            α[t, j] /= c[t]
        end

        c[t] = log(c[t]) + m
    end
end

# In-place backward pass, where β and c are allocated beforehand.
function backwardlog!(
    β::AbstractMatrix,
    c::AbstractVector,
    a::AbstractVector,
    A::AbstractMatrix,
    LL::AbstractMatrix,
)
    @argcheck size(β, 1) == size(LL, 1) == size(c, 1)
    @argcheck size(β, 2) == size(LL, 2) == size(a, 1) == size(A, 1) == size(A, 2)

    T, K = size(LL)
    L = zeros(K)
    (T == 0) && return

    fill!(β, 0.0)
    fill!(c, 0.0)

    for j in OneTo(K)
        β[end, j] = 1.0
    end

    @inbounds for t = T-1:-1:1
        m = vec_maximum(view(LL, t + 1, :))

        for i in OneTo(K)
            L[i] = exp(LL[t+1, i] - m)
        end

        for j in OneTo(K)
            for i in OneTo(K)
                β[t, j] += β[t+1, i] * A[j, i] * L[i]
            end
            c[t+1] += β[t, j]
        end

        for j in OneTo(K)
            β[t, j] /= c[t+1]
        end

        c[t+1] = log(c[t+1]) + m
    end

    m = vec_maximum(view(LL, 1, :))

    for j in OneTo(K)
        c[1] += a[j] * exp(LL[1, j] - m) * β[1, j]
    end

    c[1] = log(c[1]) + m
end

# In-place posterior computation, where γ is allocated beforehand.
function posteriors!(γ::AbstractMatrix, α::AbstractMatrix, β::AbstractMatrix)
    @argcheck size(γ) == size(α) == size(β)
    T, K = size(α)
    for t in OneTo(T)
        c = 0.0
        for i in OneTo(K)
            γ[t, i] = α[t, i] * β[t, i]
            c += γ[t, i]
        end

        for i in OneTo(K)
            γ[t, i] /= c
        end
    end
end

"""
    forward(a, A, LL) -> (Vector, Float)

Compute forward probabilities using samples likelihoods.
See [Forward-backward algorithm](https://en.wikipedia.org/wiki/Forward–backward_algorithm).

**Output**
- `Vector{Float64}`: forward probabilities.
- `Float64`: log-likelihood of the observed sequence.
"""
function forward(a::AbstractVector, A::AbstractMatrix, LL::AbstractMatrix; logl = nothing)
    (logl !== nothing) && deprecate_kwargs("logl")
    m = Matrix{Float64}(undef, size(LL))
    c = Vector{Float64}(undef, size(LL, 1))
    forwardlog!(m, c, a, A, LL)
    m, sum(c)
end

"""
    backward(a, A, LL) -> (Vector, Float)

Compute backward probabilities using samples likelihoods.
See [Forward-backward algorithm](https://en.wikipedia.org/wiki/Forward–backward_algorithm).

**Output**
- `Vector{Float64}`: backward probabilities.
- `Float64`: log-likelihood of the observed sequence.
"""
function backward(a::AbstractVector, A::AbstractMatrix, LL::AbstractMatrix; logl = nothing)
    (logl !== nothing) && deprecate_kwargs("logl")
    m = Matrix{Float64}(undef, size(LL))
    c = Vector{Float64}(undef, size(LL, 1))
    backwardlog!(m, c, a, A, LL)
    m, sum(c)
end

"""
    forward(hmm, observations; robust) -> (Vector, Float)

Compute forward probabilities of the `observations` given the `hmm` model.

**Output**
- `Vector{Float64}`: forward probabilities.
- `Float64`: log-likelihood of the observed sequence.

**Example**
```julia
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
y = rand(hmm, 1000)
probs, tot = forward(hmm, y)
```
"""
function forward(hmm::AbstractHMM, observations; logl = nothing, robust = false)
    (logl !== nothing) && deprecate_kwargs("logl")
    LL = loglikelihoods(hmm, observations; robust = robust)
    forward(hmm.a, hmm.A, LL)
end

"""
    backward(hmm, observations; robust) -> (Vector, Float)

Compute forward probabilities of the `observations` given the `hmm` model.

**Output**
- `Vector{Float64}`: backward probabilities.
- `Float64`: log-likelihood of the observed sequence.

**Example**
```julia
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
y = rand(hmm, 1000)
probs, tot = forward(hmm, y)
```
"""
function backward(hmm::AbstractHMM, observations; logl = nothing, robust = false)
    (logl !== nothing) && deprecate_kwargs("logl")
    LL = loglikelihoods(hmm, observations; robust = robust)
    backward(hmm.a, hmm.A, LL)
end

"""
    posteriors(α, β) -> Vector

Compute posterior probabilities from `α` and `β`.

**Arguments**
- `α::AbstractVector`: forward probabilities.
- `β::AbstractVector`: backward probabilities.
"""
function posteriors(α::AbstractMatrix, β::AbstractMatrix)
    γ = Matrix{Float64}(undef, size(α))
    posteriors!(γ, α, β)
    γ
end

"""
    posteriors(a, A, LL; kwargs...) -> Vector

Compute posterior probabilities using samples likelihoods.
"""
function posteriors(a::AbstractVector, A::AbstractMatrix, LL::AbstractMatrix; kwargs...)
    α, _ = forward(a, A, LL; kwargs...)
    β, _ = backward(a, A, LL; kwargs...)
    posteriors(α, β)
end

"""
    posteriors(hmm, observations; robust) -> Vector

Compute posterior probabilities using samples likelihoods.

**Example**
```julia
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
y = rand(hmm, 1000)
γ = posteriors(hmm, y)
```
"""
function posteriors(hmm::AbstractHMM, observations; logl = nothing, robust = false)
    (logl !== nothing) && deprecate_kwargs("logl")
    LL = loglikelihoods(hmm, observations; robust = robust)
    posteriors(hmm.a, hmm.A, LL)
end

"""
    loglikelihood(hmm, observations; robust) -> Float64

Compute the log-likelihood of the observations under the model.
This is defined as the sum of the log of the normalization coefficients in the forward filter.

**Output**
- `Float64`: log-likelihood of the observations sequence under the model.

**Example**
```jldoctest
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
loglikelihood(hmm, [0.15, 0.10, 1.35])
# output
-4.588183811489616
```
"""
function loglikelihood(hmm::AbstractHMM, observations; logl = nothing, robust = false)
    (logl !== nothing) && deprecate_kwargs("logl")
    forward(hmm, observations, robust = robust)[2]
end


text:
using Test
using ExtensibleEffects
using DataTypesBasic
using TypeClasses
using Documenter

splitln(str) = split(strip(str), "\n")
import ExtensibleEffects: autorun

@test isempty(detect_ambiguities(ExtensibleEffects))

if  v"1.6" <= VERSION < v"1.7"
  # doctests are super instable, hence we only do it for a specific Julia Version
  doctest(ExtensibleEffects)
end

@testset "Identity/Const" begin
  include("OptionTryEither.jl")
end

@testset "Vector/Iterable" begin
  include("VectorIterable.jl")
end

@testset "Task/Future" begin
  include("TaskFuture.jl")
end

@testset "ContextManager" begin
  include("ContextManager.jl")
end

@testset "Writer" begin
  include("Writer.jl")
end

@testset "Callable" begin
  include("Callable.jl")
end

@testset "State" begin
  include("State.jl")
end

@testset "NoEffect" begin
  include("NoEffect.jl")
end

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
export bisect_left, bisect_right, bisect

"""
    bisect_left(a, x, lo = 1, hi = length(a) + 1)

Return the index where to insert item `x` in array `a`, assuming `a` is in an
non-decreasing order.

The return value `i` is such that all `e` in `a[:(i - 1)]` have `e < x`, and all `e` in
`a[i:]` have `e >= x`.  So if `x` already appears in the array, `insert!(a, i, x)` will
insert just before the leftmost `x` already there.

# Arguments
Optional args `lo` (default `1`) and `hi` (default `length(a) + 1`) bound the
slice of `a` to be searched.

# Examples

```jldoctest
julia> bisect_left([1, 2, 3, 4, 5], 3.5)
4

julia> bisect_left([1, 2, 3, 4, 5], 2)
2

julia> bisect_left([1, 2, 3, 3, 3, 5], 3)
3
```
"""
function bisect_left(a, x, lo = 1, hi = nothing)
    if lo < 1
        throw(BoundsError(a, lo))
    end
    if hi === nothing
        hi = length(a) + 1  # It's not `length(a)`!
    end
    while lo < hi
        mid = (lo + hi) ÷ 2
        a[mid] < x ? lo = mid + 1 : hi = mid
    end
    return lo
end

"""
    bisect_right(a, x, lo = 1, hi = length(a) + 1)

Return the index where to insert item `x` in array `a`, assuming `a` is in an
non-decreasing order.

The return value `i` is such that all `e` in `a[:(i - 1)]` have `e <= x`, and all `e` in
`a[i:]` have `e > x`.  So if `x` already appears in the array, `insert!(a, i, x)` will
insert just after the rightmost `x` already there.

# Arguments
Optional args `lo` (default `1`) and `hi` (default `length(a) + 1`) bound the
slice of `a` to be searched.

# Examples

```jldoctest
julia> bisect_right([1, 2, 3, 4, 5], 3.5)
4

julia> bisect_right([1, 2, 3, 4, 5], 2)
3

julia> bisect_right([1, 2, 3, 3, 3, 5], 3)
6
```
"""
function bisect_right(a, x, lo = 1, hi = nothing)
    if lo < 1
        throw(BoundsError(a, lo))
    end
    if hi === nothing
        hi = length(a) + 1  # It's not `length(a)`!
    end
    while lo < hi
        mid = (lo + hi) ÷ 2
        x < a[mid] ? hi = mid : lo = mid + 1
    end
    return lo
end

const bisect = bisect_right


text:
using Pkg

Pkg.add("PackageCompiler")
Pkg.build("PackageCompiler")
Pkg.add("DiffEqBase")
using PackageCompiler

function installed()
    deps = Pkg.dependencies()
    installs = Dict{String, VersionNumber}()
    for (uuid, dep) in deps
        dep.is_direct_dep || continue
        dep.version === nothing && continue
        installs[dep.name] = dep.version
    end
    return installs
end

function adduninstalled(pkgs)
    for i in pkgs
        i_str = String(i)
        if i_str ∉ keys(installed())
            Pkg.add(String(i))
        end
    end
end

# Add your packages here
conservative = [
    :ForwardDiff
    :Plots
    :StaticArrays
    :BenchmarkTools
    :Parameters
    :NLsolve
    :FiniteDiff
    :LabelledArrays
    :DifferentialEquations
]

# Add packages if not exists
adduninstalled(conservative)

diffeqbase = [
    :DiffEqBase
]

PackageCompiler.create_sysimage(vcat(conservative,diffeqbase);replace_default=true)

max_stars_repo_name: stevengogogo/JuliaPackageImage

max_issues_repo_name: stevengogogo/JuliaPackageImage

max_forks_repo_name: stevengogogo/JuliaPackageImage


text:
# CheckerboardFalse
const DQMC_CBFalse = DQMC{M, CheckerboardFalse} where M

"""
    slice_matrix(mc::DQMC_CBFalse, m::Model, slice::Int, power::Float64=1.)

Direct calculation of effective slice matrix, i.e. no checkerboard.
Calculates `Beff(slice) = exp(−1/2∆tauT)exp(−1/2∆tauT)exp(−∆tauV(slice))`.
"""
@bm function slice_matrix(
        mc::DQMC_CBFalse, m::Model, slice::Int, power::Float64 = 1.0,
        conf::AbstractArray = mc.conf
    )
    eT2 = mc.stack.hopping_matrix_exp_squared
    eTinv2 = mc.stack.hopping_matrix_exp_inv_squared
    eV = mc.stack.eV

    interaction_matrix_exp!(mc, m, eV, conf, slice, power)

    if power > 0
        return eT2 * eV
    else
        return eV * eTinv2
    end
end
@bm function slice_matrix!(
        mc::DQMC_CBFalse, m::Model, slice::Int, power::Float64 = 1.0,
        result::AbstractArray = mc.stack.tmp2, conf::AbstractArray = mc.conf
    )
    eT2 = mc.stack.hopping_matrix_exp_squared
    eTinv2 = mc.stack.hopping_matrix_exp_inv_squared
    eV = mc.stack.eV

    interaction_matrix_exp!(mc, m, eV, conf, slice, power)

    if power > 0
        # eT * (eT * eV)
        vmul!(result, eT2, eV)
    else
        # ev * (eTinv * eTinv)
        vmul!(result, eV, eTinv2)
    end
    return result
end

@bm function multiply_slice_matrix_left!(
        mc::DQMC_CBFalse, m::Model, slice::Int, M::AbstractMatrix,
        conf::AbstractArray = mc.conf
    )
    slice_matrix!(mc, m, slice, 1.0, mc.stack.tmp2, conf)
    vmul!(mc.stack.tmp1, mc.stack.tmp2, M)
    M .= mc.stack.tmp1
    nothing
end
@bm function multiply_slice_matrix_right!(
        mc::DQMC_CBFalse, m::Model, slice::Int, M::AbstractMatrix,
        conf::AbstractArray = mc.conf
    )
    slice_matrix!(mc, m, slice, 1.0, mc.stack.tmp2, conf)
    vmul!(mc.stack.tmp1, M, mc.stack.tmp2)
    M .= mc.stack.tmp1
    nothing
end
@bm function multiply_slice_matrix_inv_right!(
        mc::DQMC_CBFalse, m::Model, slice::Int, M::AbstractMatrix,
        conf::AbstractArray = mc.conf
    )
    slice_matrix!(mc, m, slice, -1.0, mc.stack.tmp2, conf)
    vmul!(mc.stack.tmp1, M, mc.stack.tmp2)
    M .= mc.stack.tmp1
    nothing
end
@bm function multiply_slice_matrix_inv_left!(
        mc::DQMC_CBFalse, m::Model, slice::Int, M::AbstractMatrix,
        conf::AbstractArray = mc.conf
    )
    slice_matrix!(mc, m, slice, -1.0, mc.stack.tmp2, conf)
    vmul!(mc.stack.tmp1, mc.stack.tmp2, M)
    M .= mc.stack.tmp1
    nothing
end
@bm function multiply_daggered_slice_matrix_left!(
        mc::DQMC_CBFalse, m::Model, slice::Int, M::AbstractMatrix,
        conf::AbstractArray = mc.conf
    )
    slice_matrix!(mc, m, slice, 1.0, mc.stack.tmp2, conf)
    vmul!(mc.stack.tmp1, adjoint(mc.stack.tmp2), M)
    M .= mc.stack.tmp1
    nothing
end

################################################################################
### CheckerboardTrue
################################################################################

# TODO:
# allow passing of conf for global updates

const DQMC_CBTrue = DQMC{M, CheckerboardTrue} where M

@bm function slice_matrix(mc::DQMC_CBTrue, m::Model, slice::Int,
                    power::Float64=1.)
    N = length(lattice(m)) * nflavors(m)
    M = Matrix{heltype(mc)}(I, N, N)
    if power > 0
        multiply_slice_matrix_left!(mc, m, slice, M)
    else
        multiply_slice_matrix_inv_left!(mc, m, slice, M)
    end
    return M
end
@bm function slice_matrix!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    power::Float64=1., M = mc.stack.U)
    copyto!(M, I)
    if power > 0
        multiply_slice_matrix_left!(mc, m, slice, M)
    else
        multiply_slice_matrix_inv_left!(mc, m, slice, M)
    end
    return M
end

@bm function multiply_slice_matrix_left!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    M::AbstractMatrix{T}) where T<:Number
    s = mc.stack
    interaction_matrix_exp!(mc, m, s.eV, mc.conf, slice, 1.)

    mul!(s.tmp1, s.eV, M)
    M .= s.tmp1
    mul!(s.tmp1, s.chkr_mu, M)
    M .= s.tmp1

    @inbounds begin
        for i in reverse(2:s.n_groups)
            mul!(s.tmp1, s.chkr_hop_half[i], M)
            M .= s.tmp1
        end
        mul!(s.tmp1, s.chkr_hop[1], M)
        M .= s.tmp1
        for i in 2:s.n_groups
            mul!(s.tmp1, s.chkr_hop_half[i], M)
            M .= s.tmp1
        end
    end
    nothing
end
@bm function multiply_slice_matrix_right!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    M::AbstractMatrix{T}) where T<:Number
    s = mc.stack
    @inbounds begin
        for i in reverse(2:s.n_groups)
            mul!(s.tmp1, M, s.chkr_hop_half[i])
            M .= s.tmp1
        end
        mul!(s.tmp1, M, s.chkr_hop[1])
        M .= s.tmp1
        for i in 2:s.n_groups
            mul!(s.tmp1, M, s.chkr_hop_half[i])
            M .= s.tmp1
        end
    end

    interaction_matrix_exp!(mc, m, s.eV, mc.conf, slice, 1.)
    mul!(s.tmp1, M, s.chkr_mu)
    M .= s.tmp1
    mul!(s.tmp1, M, s.eV)
    M .= s.tmp1
    nothing
end
@bm function multiply_slice_matrix_inv_left!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    M::AbstractMatrix{T}) where T<:Number
    s = mc.stack
    @inbounds begin
        for i in reverse(2:s.n_groups)
            mul!(s.tmp1, s.chkr_hop_half_inv[i], M)
            M .= s.tmp1
        end
        mul!(s.tmp1, s.chkr_hop_inv[1], M)
        M .= s.tmp1
        for i in 2:s.n_groups
            mul!(s.tmp1, s.chkr_hop_half_inv[i], M)
            M .= s.tmp1
        end
    end

    interaction_matrix_exp!(mc, m, s.eV, mc.conf, slice, -1.)
    mul!(s.tmp1, s.chkr_mu_inv, M)
    M .= s.tmp1
    mul!(s.tmp1, s.eV, M)
    M .= s.tmp1
    nothing
end
@bm function multiply_slice_matrix_inv_right!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    M::AbstractMatrix{T}) where T<:Number
    s = mc.stack
    interaction_matrix_exp!(mc, m, s.eV, mc.conf, slice, -1.)
    mul!(s.tmp1, M, s.eV)
    M .= s.tmp1
    mul!(s.tmp1, M, s.chkr_mu_inv)
    M .= s.tmp1

    @inbounds begin
        for i in reverse(2:s.n_groups)
            mul!(s.tmp1, M, s.chkr_hop_half_inv[i])
            M .= s.tmp1
        end
        mul!(s.tmp1, M, s.chkr_hop_inv[1])
        M .= s.tmp1
        for i in 2:s.n_groups
            mul!(s.tmp1, M, s.chkr_hop_half_inv[i])
            M .= s.tmp1
        end
    end
    nothing
end
@bm function multiply_daggered_slice_matrix_left!(mc::DQMC_CBTrue, m::Model, slice::Int,
                    M::AbstractMatrix{T}) where T<:Number
    s = mc.stack

    @inbounds begin
        for i in reverse(2:s.n_groups)
            mul!(s.tmp1, s.chkr_hop_half_dagger[i], M)
            M .= s.tmp1
        end
        mul!(s.tmp1, s.chkr_hop_dagger[1], M)
        M .= s.tmp1
        for i in 2:s.n_groups
            mul!(s.tmp1, s.chkr_hop_half_dagger[i], M)
            M .= s.tmp1
        end
    end

    interaction_matrix_exp!(mc, m, s.eV, mc.conf, slice, 1.)
    # s.eV == adjoint(s.eV) and s.chkr_mu == adjoint(s.chkr_mu)
    mul!(s.tmp1, s.chkr_mu, M)
    M .= s.tmp1
    mul!(s.tmp1, s.eV, M)
    M .= s.tmp1
    nothing
end


text:
# Autogenerated wrapper script for Readline_jll for i686-w64-mingw32
export libhistory, libreadline

using Ncurses_jll
JLLWrappers.@generate_wrapper_header("Readline")
JLLWrappers.@declare_library_product(libhistory, "libhistory8.dll")
JLLWrappers.@declare_library_product(libreadline, "libreadline8.dll")
function __init__()
    JLLWrappers.@generate_init_header(Ncurses_jll)
    JLLWrappers.@init_library_product(
        libhistory,
        "bin\\libhistory8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libreadline,
        "bin\\libreadline8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using RegressionBenchmarks, Gurobi
using Base.Test, Distributions

srand(1)
n = [20]
d = 10
sparsity = 5
Xdist = MvNormal
Xcorr = MatrixCorrelation(0.1)
Xdata = XData(Xdist, Xcorr, d)
wdist = BinChoice()
w = getw(d, sparsity, wdist)
X = getX(n[1], d, Xdata)
Y = X * w

@testset "Types" begin
    # We are getting Y = X * w
    @test isapprox(Y[1], -1.72637, atol=1e-4)

    # We don't allow to specify X with covariance matrix of different dimensions
    # to the number of features
    @test_throws ErrorException getdata(Xdata = Xdata,
        wdist = BinChoice(),
        noisedist = NoNoise(),
        n = 100, nfeatures = 100, sparsity = 10)

    # We should also get Y = X * w with this interface
    srand(1)
    rd = getdata(Xdata = Xdata, wdist = BinChoice(), n = n[1], nfeatures = d, sparsity = sparsity)
    @test isapprox(rd.Y[1], -1.72637, atol=1e-4)

    # ... and we should also get Y = X * w with this interface
    srand(1)
    bd = BenchmarkData(Xdata = Xdata, wdist = BinChoice(),
        noisedist = NoNoise(), SNR = 0.0, n = n, nfeatures = d, sparsity = sparsity)
    rd = getdata(bd, 1)
    @test isapprox(rd.Y[1], -1.72637, atol=1e-4)

    # We don't allow stupid values for the SNR (currently 0) if user wants noise
    bd.noisedist = Normal()
    @test_throws ErrorException getdata(bd, 1)

end
@testset "Utils" begin
    @testset "Accuracy and False Alarm" begin
        pred = [1; 2; 5]
        truth = [1; 2; 3]
        a = accuracy(pred, truth)
        f = falsepositive(pred, truth)
        @test isapprox(a, 2 / 3, atol = 1e-6)
        @test isapprox(f, 1 / 3, atol = 1e-6)
    end
    @testset "Saving results" begin
        bd = BenchmarkData(Xdata = Xdata, wdist = BinChoice(),
            noisedist = NoNoise(), SNR = 0.0, n = n, nfeatures = d,
            sparsity = sparsity)
        @test data2str(bd) ==
            "x_normal_corr_rho_0.1_w_binchoice_noise_nonoise_snr_0.0_d_$(d)_k_$(sparsity)"
        m = RelaxDualSubgradient(5.0, PolyakStepping(2.0, 30), 100)
        @test method2str(m) ==
            "relax_dual_subgradient_polyak_2.0_maxiter_100"
        m = ExactPrimalCuttingPlane(0.1, 30.0, GurobiSolver, NoWarmStart())
        @test method2str(m) ==
            "exact_primal_tlimit_30.0_gurobi"
        m = ExactPrimalCuttingPlane(0.1, 30.0, GurobiSolver, RelaxDualSubgradient())
        @test method2str(m) ==
            "exact_primal_tlimit_30.0_gurobi_ws"
    end
end

@testset "getslope!" begin
    ∇ = zeros(n[1])
    α = 0.1 * ones(n[1])
    ax = X[:, 1:sparsity]' * α
    inds = collect(1:sparsity)
    γ = 5.0
    c = SubsetSelection.Cache(n[1], d)
    RegressionBenchmarks.getslope!(∇, X, Y, α, inds, γ, sparsity, ax)
    g = SubsetSelection.grad_dual(SubsetSelection.OLS(), Y, X, -α, inds, sparsity, γ, c)
    @test ∇ ≈ -g
    db = RegressionBenchmarks.dual_bound(SubsetSelection.OLS(), Y, X, -α, inds, sparsity, γ, c)
    lb = RegressionBenchmarks.getlb(Y, α, ax, γ)
    @test db ≈ lb
end

@testset "Methods" begin
    srand(1)
    bd = BenchmarkData(Xdata = Xdata, wdist = BinChoice(),
        noisedist = NoNoise(), SNR = 0.0, n = n, nfeatures = d, sparsity = sparsity)
    rd = getdata(bd, 1)
    @testset "Subgradient" begin
        @testset "Ordinary" begin
            m = RelaxDualSubgradient(5.0, ConstantStepping(1e-3), 100)
            indices, weights = solve_problem(m, rd.X, rd.Y, sparsity)
            @test indices == find(abs.(w) .> 1e-6)
            @test isapprox(weights, w[abs.(w) .> 1e-6], atol = 1e-1)
        end
        @testset "Varying stepsizes" begin
            m = RelaxDualSubgradient(50.0, PolyakStepping(0.5, 30), 100)
            indices, weights = solve_problem(m, rd.X, rd.Y, sparsity)
            @test indices == find(abs.(w) .> 1e-6)
            @test isapprox(weights, w[abs.(w) .> 1e-6], atol = 1e-1)
        end
    end
    @testset "Dual cutting planes" begin
        lpsolver = GurobiSolver(OutputFlag = 0)
        m = RelaxDualCutting(5.0, lpsolver)
        indices, weights = solve_problem(m, rd.X, rd.Y, sparsity)
        @test indices == find(abs.(w) .> 1e-6)
        @test isapprox(weights, w[abs.(w) .> 1e-6], atol = 1e-1)
    end
    @testset "Warm starts" begin
        warm_start = RelaxDualSubgradient(50.0, PolyakStepping(0.5, 30), 100)
        m = ExactPrimalCuttingPlane(5.0, 10.0, GurobiSolver, warm_start)
        indices, weights = solve_problem(m, rd.X, rd.Y, sparsity)
        @test indices == find(abs.(w) .> 1e-6)
        @test isapprox(weights, w[abs.(w) .> 1e-6], atol = 1e-1)
    end
end

@testset "Validation" begin
    srand(1)
    bd = BenchmarkData(Xdata = Xdata, wdist = BinChoice(),
        noisedist = NoNoise(), SNR = 0.0, n = n, nfeatures = d, sparsity = sparsity)
    rd = getdata(bd, 1)
    for method in [ExactPrimalCuttingPlane(GurobiSolver);
                   PrimalWithHeuristics(GurobiSolver);
                   RelaxDualSubgradient();
                   RelaxDualSubgradient(50.0, PolyakStepping(), 100)
                  ]
        v_results = validate_params!(rd.X, rd.Y, bd.sparsity, method)
        @test any(isapprox.(v_results.valid_scores, 1.0, atol=1e-2))
        if isa(method, Union{ExactPrimalCuttingPlane, PrimalWithHeuristics})
            @test method.time_limit == 30.0
        end
    end
end


text:
using Printf

parse_line(line) = [c == '.' ? false : true for c in line]

function load(path)
  # NOTE: I do not know why vcat doesn't work here, but this is just the input loader, so I don't
  # care too much
  hcat(map(parse_line, readlines(path))...)'
end

function check_slopes(data)
  y = 1
  # NOTE: We could do this in a more general way (maybe?) but I'm feeling lazy
  trees :: Array{Int} = zeros(5)
  rows = size(data, 1)
  cols = size(data, 2)
  while y <= rows
    for (i, h) in enumerate((1, 3, 5, 7))
      x = (((y - 1) * h) % cols) + 1
      if data[y, x]
        trees[i] += 1
      end
    end

    if (y - 1) % 2 == 0 && data[y, (div(y, 2) % cols) + 1]
      trees[5] += 1
    end

    y += 1
  end

  trees, reduce(*, trees)
end

function solve(data)
  println(@time check_slopes(data))
end

data = load("input.txt")
solve(data)

max_stars_repo_name: wbthomason/advent

max_issues_repo_name: wbthomason/advent

max_forks_repo_name: wbthomason/advent


text:
# Postcompose with scaling and sum

export Postcompose

"""
**Postcomposition with an affine transformation**

    Postcompose(f, a=1.0, b=0.0)

Returns the function
```math
g(x) = a\\cdot f(x) + b.
```
"""
struct Postcompose{T <: ProximableFunction, R <: Real} <: ProximableFunction
    f::T
    a::R
    b::R
    function Postcompose{T,R}(f::T, a::R, b::R) where {T <: ProximableFunction, R <: Real}
        if a <= 0.0
            error("parameter `a` must be positive")
        else
            new(f, a, b)
        end
    end
end

is_prox_accurate(f::Postcompose) = is_prox_accurate(f.f)
is_separable(f::Postcompose) = is_separable(f.f)
is_convex(f::Postcompose) = is_convex(f.f)
is_set(f::Postcompose) = is_set(f.f)
is_singleton(f::Postcompose) = is_singleton(f.f)
is_cone(f::Postcompose) = is_cone(f.f)
is_affine(f::Postcompose) = is_affine(f.f)
is_smooth(f::Postcompose) = is_smooth(f.f)
is_quadratic(f::Postcompose) = is_quadratic(f.f)
is_generalized_quadratic(f::Postcompose) = is_generalized_quadratic(f.f)
is_strongly_convex(f::Postcompose) = is_strongly_convex(f.f)

Postcompose(f::T, a::R=R(1), b::R=R(0)) where {T <: ProximableFunction, R <: Real} = Postcompose{T, R}(f, a, b)

Postcompose(f::Postcompose{T, R}, a::R=R(1), b::R=R(0)) where {T <: ProximableFunction, R <: Real} = Postcompose{T, R}(f.f, a*f.a, b+a*f.b)

function (g::Postcompose)(x::AbstractArray{T}) where T <: RealOrComplex
    return g.a*g.f(x) + g.b
end

function gradient!(y::AbstractArray{T}, g::Postcompose, x::AbstractArray{T}) where T <: RealOrComplex
    v = gradient!(y, g.f, x)
    y .*= g.a
    return g.a*v + g.b
end

function prox!(y::AbstractArray{T}, g::Postcompose, x::AbstractArray{T}, gamma=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    v = prox!(y, g.f, x, g.a*gamma)
    return g.a*v + g.b
end

function prox_naive(g::Postcompose, x::AbstractArray{T}, gamma=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    y, v = prox_naive(g.f, x, g.a*gamma)
    return y, g.a*v + g.b
end

fun_name(f::Postcompose) = string("Postcomposition of ", fun_name(f.f))
fun_dom(f::Postcompose) = fun_dom(f.f)
fun_expr(f::Postcompose) = "x ↦ a*f(x)+b"
fun_params(f::Postcompose) = string("f(x) = ", fun_expr(f.f), ", a = $(f.a), b = $(f.b)")


text:
using Colors, TestImages, ImageIO, ImageShow, FileIO
using LabView0mqJl, Serialization

# img = testimage("lena_color_256")
# bd = LabView0mqJl.BinDescr(5, 25, [3,4], "mytype", "kvn", "mummarrays")
# sers = (;img, bd)

pth = raw"C:\_LabView_projects\ZMQ\LabView0mqJl.jl\src\Julia\this-and-that\sers.jsr"
      # joinpath(@__DIR__, "sers.jsr")

open(pth, "w") do io
    @show typeof(LabView0mqJl.gbd) typeof(LabView0mqJl.gbdd)
    serialize(io, LabView0mqJl.gbd)
    serialize(io, LabView0mqJl.gbdd)
end;

# sers2=nothing
# open(pth, "r") do io
#     global sers2=deserialize(io)
#     # @show typeof(sers2)
# end;

nbd = nbdd = nothing
open(pth, "r") do io
    global nbd=deserialize(io)
    global nbdd=deserialize(io)

end;
@show typeof(nbd)
@show typeof(nbdd)


text:
@testset "Partial MD" begin
    st = raw"""
        \newcommand{\com}{HH}
        \newcommand{\comb}[1]{HH#1HH}
        A list
        * \com and \comb{blah}
        * $f$ is a function
        * a last element
        """

    steps = explore_md_steps(st)
    lxdefs, tokens, braces, blocks, lxcoms = steps[:latex]

    @test length(braces) == 1
    @test F.content(braces[1]) == "blah"

    @test length(blocks) == 1
    @test blocks[1].name == :MATH_A
    @test F.content(blocks[1]) == "f"

    b2insert, = steps[:b2insert]

    inter_md, mblocks = F.form_inter_md(st, b2insert, lxdefs)
    @test inter_md == "\n\nA list\n*  ##FDINSERT##  and  ##FDINSERT## \n*  ##FDINSERT##  is a function\n* a last element\n"
    inter_html = F.md2html(inter_md)
    @test inter_html == "
A list
\n
\n

##FDINSERT##  and  ##FDINSERT##
\n
\n

##FDINSERT##  is a function
\n
\n

a last element
\n
\n
\n"
end

# index arithmetic over a string is a bit trickier when using all symbols
# we can use `prevind` and `nextind` to make sure it works properly
@testset "Inter Md 2" begin
    st = raw"""
        ~~~
        this⊙ then ⊙ ⊙ and
        ~~~
        finally ⊙⊙𝛴⊙ and
        ~~~
        escape ∀⊙∀
        ~~~
        done
        """

    inter_md, = explore_md_steps(st)[:inter_md]
    @test inter_md == " ##FDINSERT## \nfinally ⊙⊙𝛴⊙ and\n ##FDINSERT## \ndone\n"
end

@testset "Latex eqa" begin
    st = raw"""
        a\newcommand{\eqa}[1]{\begin{eqnarray}#1\end{eqnarray}}b@@d .@@
        \eqa{\sin^2(x)+\cos^2(x) &=& 1}
        """

    steps = explore_md_steps(st)
    lxdefs, tokens, braces, blocks, lxcoms = steps[:latex]
    b2insert, = steps[:b2insert]
    inter_md, mblocks = steps[:inter_md]
    @test inter_md == "ab ##FDINSERT## \n ##FDINSERT## \n"

    inter_html, = steps[:inter_html]

    @test F.convert_block(b2insert[1], lxdefs) == "
.
"
    @test isapproxstr(F.convert_block(b2insert[2], lxdefs), "\\[\\begin{array}{c} \\sin^2(x)+\\cos^2(x) &=& 1\\end{array}\\]")
    hstring = F.convert_inter_html(inter_html, b2insert, lxdefs)
    @test isapproxstr(hstring, raw"""

                          ab
.

                          \[\begin{array}{c}
                            \sin^2(x)+\cos^2(x) &=& 1
                          \end{array}\]

""")
end

@testset "MD>HTML" begin
    st = raw"""
        text A1 \newcommand{\com}{blah}text A2 \com and
        ~~~
        escape B1
        ~~~
        \newcommand{\comb}[ 1]{\mathrm{#1}} text C1 $\comb{b}$ text C2
        \newcommand{\comc}[ 2]{part1:#1 and part2:#2} then \comc{AA}{BB}.
        """

    steps = explore_md_steps(st)
    lxdefs, tokens, braces, blocks, lxcoms = steps[:latex]
    b2insert, = steps[:b2insert]
    inter_md, mblocks = steps[:inter_md]
    inter_html, = steps[:inter_html]

    @test isapproxstr(inter_md, """
                                text A1 text A2  ##FDINSERT##  and
                                ##FDINSERT##
                                text C1  ##FDINSERT##  text C2
                                then  ##FDINSERT## .""")

    @test isapproxstr(inter_html, """
text A1 text A2  ##FDINSERT##  and  ##FDINSERT##   text C1  ##FDINSERT##  text C2  then  ##FDINSERT## .
""")

    hstring = F.convert_inter_html(inter_html, b2insert, lxdefs)
    @test isapproxstr(hstring, """

text A1 text A2 blah and
                                escape B1
                                text C1 \\(\\mathrm{ b}\\) text C2
                                then part1: AA and part2: BB.
""")
end

@testset "headers" begin
    set_curpath("index.md")
    h = """
        # Title
        and then
        ## Subtitle cool!
        done
        """ |> seval
    @test isapproxstr(h, """

Title

                        and then

Subtitle cool!

                        done
                        """)
end


text:
"""
A [`Plugin`](@ref) that only adds a [`Badge`](@ref) to the [`Readme`](@ref) file.

Concrete subtypes only need to implement a [`badges`](@ref) method.
"""
abstract type BadgePlugin <: Plugin end

"""
    BlueStyleBadge()

Adds a [`BlueStyle`](https://github.com/invenia/BlueStyle) badge to the [`Readme`](@ref) file.
"""
struct BlueStyleBadge <: BadgePlugin end

function badges(::BlueStyleBadge)
    return Badge(
        "Code Style: Blue",
        "https://img.shields.io/badge/code%20style-blue-4495d1.svg",
        "https://github.com/invenia/BlueStyle",
    )
end

"""
    ColPracBadge()

Adds a [`ColPrac`](https://github.com/SciML/ColPrac) badge to the [`Readme`](@ref) file.
"""
struct ColPracBadge <: BadgePlugin end

function badges(::ColPracBadge)
    return Badge(
        "ColPrac: Contributor's Guide on Collaborative Practices for Community Packages",
        "https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet",
        "https://github.com/SciML/ColPrac",
    )
end


text:
using DataFrames
using CSV
using HTTP
using Optim
using GLM
using LinearAlgebra
using Random
using Statistics
using DataFrames
using FreqTables

 #:::::::::::::::::::::::::::::::::::::::::::::::::::
 # question 1
 #:::::::::::::::::::::::::::::::::::::::::::::::::::

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2021/master/ProblemSets/PS3-gev/nlsw88w.csv"
df = CSV.read(HTTP.get(url).body,DataFrame)
X = [df.age df.white df.collgrad]
Z = hcat(df.elnwage1, df.elnwage2, df.elnwage3, df.elnwage4,
df.elnwage5, df.elnwage6, df.elnwage7, df.elnwage8)
y = df.occupation

function mlogit(alpha, X, Z, y)
    J = size(unique(y),1)
    beta = reshape(alpha[1:end-1],size(X,2),J-1)
    gamma = alpha[end]
    Z_diff = Z[:,1:end-1].-Z[:,end]
    loglike = 0
    for i = 1:size(X,1)
        term1 = sum([convert(Float64,y[i]==j)*(sum(beta[:,j].*X[i,:])+sum(gamma.*Z_diff[i,j])) for j =1:J-1])
        term2 = log(1 + sum([exp((sum(beta[:,j].*X[i,:])+sum(gamma.*Z_diff[i,j]))) for j =1:J-1]))
        loglike += term1 - term2
    end
    return -loglike
end

J = size(unique(y),1)
start_beta = rand((J-1)*size(X,2))
start_gamma = rand(1)
alpha0 = vcat(start_beta,start_gamma)
obj = optimize(alpha->mlogit(alpha, X, Z, y), alpha0, LBFGS(), Optim.Options(g_tol=1e-5, iterations=100_000, show_trace=true, show_every=50))
beta = obj.minimizer[1:end-1]
gamma = obj.minimizer[end]
println("beta is ",beta)
println("gamma is ",gamma)

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 2
#:::::::::::::::::::::::::::::::::::::::::::::::::::

#lamdba is the coefficient of marginal effect of wage difference.

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 3
#:::::::::::::::::::::::::::::::::::::::::::::::::::

function nestedlogit(alpha, X, Z, y)
    K1= size(X,2)
    K2= size(Z,2)
    #WC:1 Professional/Technical; 2 Managers/Administrators; 3 Sales
    J_WC = 3
    s0 = 1
    s1 = K1*J_WC
    beta_wc = reshape(alpha[s0:s1],K1,J_WC)
    #BC: 4 Clerical/Unskilled; 5 Craftsmen;6 Operatives; 7 Transport
    J_BC = 4
    s2 = s1 + 1
    s3 = s1 + K1 * J_BC
    beta_bc = reshape(alpha[s2:s3],K1,J_BC)

    lambda_wc = alpha[s3+1]
    lambda_bc = alpha[s3+2]
    gamma = alpha[s3+3]
    # Other occupations
    J_O = 1
    beta_O = zeros(K1)

    beta = [beta_wc beta_bc beta_O]

    N = length(y)
    WC = zeros(N,J_WC)
    BC = zeros(N,J_BC)
    O = zeros(N,J_O)
    J = length(unique(y))
    Z_diff = Z[:,1:end].-Z[:,end]

    bigY = zeros(N,J)
    for j=1:J
        bigY[:,j] = y.==j
    end

    num = zeros(N,J)
    num1 = zeros(N,J)
    num2 = zeros(N,3)
    dem = zeros(N)
    for j=1:J
        if j <= 2
            num1[:,j] = exp.((X*beta[:,j]+gamma*Z_diff[:,j])/lambda_wc)
            num2[:,1] .+= num1[:,j]
        elseif j == 3
            num1[:,j] = exp.((X*beta[:,j]+gamma*Z_diff[:,j])/lambda_wc)
            num2[:,1] .+= num1[:,j]

            dem .+= num2[:,1].^lambda_wc
            for m =1:3
                num[:,m]=num1[:,m].*num2[:,1].^(lambda_wc-1)
            end

        elseif j >3 && j<7
            num1[:,j] = exp.((X*beta[:,j]+gamma*Z_diff[:,j])/lambda_bc)
            num2[:,2] .+= num1[:,j]

        elseif j == 7
            num1[:,j] = exp.((X*beta[:,j]+gamma*Z_diff[:,j])/lambda_bc)
            num2[:,2] .+= num1[:,j]

            dem .+= num2[:,2].^lambda_bc
            for m =4:7
                num[:,m]=num1[:,m].*num2[:,2].^(lambda_bc-1)
            end

        else
            num1[:,j] = exp.(X*beta[:,j]+gamma*Z_diff[:,j])
            num2[:,3] .+= num1[:,j]

            num[:,8] = num2[:,3]
            dem .+= num2[:,3]
        end

    end

    P = num./repeat(dem,1,J)
    loglike = -sum( bigY.*log.(P))

    return loglike
end

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# question 4
#:::::::::::::::::::::::::::::::::::::::::::::::::::
K1= size(X,2)
J_WC = 3
J_BC = 4
J = length(unique(y))
alpha0 = rand(K1*J+3)
alpha_hat_optim = optimize(alpha -> nestedlogit(alpha, X, Z, y), alpha0, LBFGS(), Optim.Options(g_tol = 1e-5, iterations=1000, show_trace=true, show_every=50))
alpha_hat_mle = alpha_hat_optim.minimizer
#println(alpha_hat_mle)
println("beta_wc is ",alpha_hat_mle[1:K1*J_WC])
println("beta_bc is ",alpha_hat_mle[K1*J_WC+1:K1*J_WC+K1*J_BC])
println("lambda_wc is ",alpha_hat_mle[K1*J_WC+K1*J_BC+1])
println("lambda_bc is ",alpha_hat_mle[K1*J_WC+K1*J_BC+2])
println("gamma is ",alpha_hat_mle[K1*J_WC+K1*J_BC+3])


text:
function sparsification(m₀,m₁,μ,ξ,Ω,graph)

    sg = size(graph,1)
    rows = ones(length(m₁))
    cols = m₁ .+ 1
    adjacency = sparse(rows,cols,true,sg,sg)
    graph .⊻= adjacency .| adjacency'

    edges = [[i+1 j+1] for i ∈ union(m₀,m₁),j ∈ m₁ if i ≠ j]
    edges = vcat(edges...)
    adjacency = sparse(edges[:,1],edges[:,2],true,sg,sg)
    graph .|= adjacency .| adjacency'

    ind₀ = [[2i+2,2i+3] for i ∈ sort(m₀)]
    ind₀ = vcat(ind₀...)
    ind₁ = [[2i+2,2i+3] for i ∈ sort(m₁)]
    ind₁ = vcat(ind₁...)

    Ω₂ = Ω[[1:3;ind₀;ind₁],ind₁]
    l₁ = Ω₂ / Symmetric(Ω[ind₁,ind₁]) * Ω₂'

    Ω₃ = Ω[[1:3;ind₀;ind₁],[1:3;ind₁]]
    l₂ = Ω₃ / Symmetric(Ω[[1:3;ind₁],[1:3;ind₁]]) * Ω₃'

    Ω₄ = Ω[:,1:3]
    l₃ = Ω₄ / Symmetric(Ω[1:3,1:3]) * Ω₄'

    Ω₁ = copy(Ω)
    Ω[[1:3;ind₀;ind₁],[1:3;ind₀;ind₁]] .-= l₁
    Ω[[1:3;ind₀;ind₁],[1:3;ind₀;ind₁]] .+= l₂
    Ω .-= l₃

    ξ .+= (Ω-Ω₁)*μ
end


text:
"""
**QuadratureVI**

Variational Inference solver by approximating gradients via numerical integration via Quadrature

```julia
QuadratureVI(ϵ::T=1e-5,nGaussHermite::Integer=20,optimizer::Optimizer=Momentum(η=0.0001))
```

**Keyword arguments**

    - `ϵ::T` : convergence criteria
    - `nGaussHermite::Int` : Number of points for the integral estimation
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl](https://github.com/jacobcvt12/GradDescent.jl) package. Default is `Momentum(η=0.0001)`
"""
mutable struct QuadratureVI{T<:Real} <: NumericalVI{T}
    ϵ::T #Convergence criteria
    nIter::Integer #Number of steps performed
    optimizer::LatentArray{Optimizer} #Learning rate for stochastic updates
    nPoints::Int64 #Number of points for the quadrature
    nodes::Vector{T}
    weights::Vector{T}
    clipping::T
    Stochastic::Bool #Use of mini-batches
    nSamples::Int64 #Number of samples of the data
    nSamplesUsed::Int64 #Size of mini-batches
    MBIndices::Vector #Indices of the minibatch
    ρ::T #Stochastic Coefficient
    HyperParametersUpdated::Bool #To know if the inverse kernel matrix must updated
    ∇η₁::LatentArray{Vector{T}}
    ∇η₂::LatentArray{Symmetric{T,Matrix{T}}}
    ν::LatentArray{Vector{T}} #Derivative -<dv/dx>_qn
    λ::LatentArray{Vector{T}} #Derivative  <d²V/dx²>_qm
    x::SubArray{T,2,Matrix{T}}#,Tuple{UnitRange{Int64},Base.Slice{Base.OneTo{Int64}}},false}#SubArray{T,2,Matrix{T},Tuple{Base.Slice{Base.OneTo{Int64}},Base.Slice{Base.OneTo{Int64}}},true}
    y::LatentArray{SubArray}
    function QuadratureVI{T}(ϵ::T,nPoints::Integer,nIter::Integer,optimizer::Optimizer,Stochastic::Bool,clipping::Real,nSamplesUsed::Integer=1) where T
        gh = gausshermite(nPoints)
        return new{T}(ϵ,nIter,[optimizer],nPoints,gh[1],gh[2]./sqrtπ,clipping,Stochastic,1,nSamplesUsed)
    end
end

function QuadratureVI(;ϵ::T=1e-5,nGaussHermite::Integer=1000,optimizer::Optimizer=Momentum(η=1e-5),clipping::Real=0.0) where {T<:Real}
    QuadratureVI{T}(ϵ,nGaussHermite,0,optimizer,false,clipping)
end

"""
**QuadratureSVI**

Stochastic Variational Inference solver by approximating gradients via numerical integration via Quadrature

```julia
QuadratureSVI(nMinibatch::Integer;ϵ::T=1e-5,nGaussHermite::Integer=20,optimizer::Optimizer=Adam(α=0.1))
```
    -`nMinibatch::Integer` : Number of samples per mini-batches

**Keyword arguments**

    - `ϵ::T` : convergence criteria, which can be user defined
    - `nGaussHermite::Int` : Number of points for the integral estimation (for the QuadratureVI)
    - `optimizer::Optimizer` : Optimizer used for the variational updates. Should be an Optimizer object from the [GradDescent.jl](https://github.com/jacobcvt12/GradDescent.jl) package. Default is `Momentum(η=0.001)`
"""
function QuadratureSVI(nMinibatch::Integer;ϵ::T=1e-5,nGaussHermite::Integer=20,optimizer::Optimizer=Momentum(η=1e-5),clipping::Real=0.0) where {T<:Real}
    QuadratureVI{T}(ϵ,nGaussHermite,0,optimizer,true,clipping,nMinibatch)
end

function expecLogLikelihood(model::VGP{T,L,<:QuadratureVI}) where {T,L}
    tot = 0.0
    for k in 1:model.nLatent
        for i in 1:model.nSample
            nodes = model.inference.nodes*sqrt2*sqrt(model.Σ[k][i,i]) .+ model.μ[k][i]
            tot += dot(model.inference.weights,logpdf.(model.likelihood,model.y[k][i],nodes))
        end
    end
    return tot
end

function expecLogLikelihood(model::SVGP{T,L,<:QuadratureVI}) where {T,L}
    tot = 0.0
    for k in 1:model.nLatent
        k_correct = model.nLatent == 1 ? 1 : k
        μ = model.κ[k_correct]*model.μ[k]
        Σ = opt_diag(model.κ[k_correct]*model.Σ[k],model.κ[k_correct])
        for i in 1:model.inference.nSamplesUsed
            nodes = model.inference.nodes*sqrt2*sqrt(Σ[i]) .+ μ[i]
            tot += dot(model.inference.weights,logpdf.(model.likelihood,model.inference.y[k][i],nodes))
        end
    end
    return tot
end

function compute_grad_expectations!(model::VGP{T,L,<:QuadratureVI}) where {T,L}
    for k in 1:model.nLatent
        for i in 1:model.nSample
            model.inference.ν[k][i], model.inference.λ[k][i] = grad_quad(model.likelihood, model.y[k][i], model.μ[k][i], model.Σ[k][i,i],model.inference)
        end
    end
end

function compute_grad_expectations!(model::SVGP{T,L,<:QuadratureVI}) where {T,L}
    for k in 1:model.nLatent
        k_correct = model.nLatent == 1 ? 1 : k
        μ = model.κ[k_correct]*model.μ[k]
        Σ = opt_diag(model.κ[k_correct]*model.Σ[k],model.κ[k_correct])
        for i in 1:model.inference.nSamplesUsed
            model.inference.ν[k][i], model.inference.λ[k][i] = grad_quad(model.likelihood, model.inference.y[k][i], μ[i], Σ[i], model.inference)
        end
    end
end

#Compute the first and second derivative of the log-likelihood using the quadrature nodes
function grad_quad(likelihood::Likelihood{T},y::Real,μ::Real,σ²::Real,inference::Inference) where {T<:Real}
    nodes = inference.nodes*sqrt2*sqrt(σ²) .+ μ
    Edlogpdf = dot(inference.weights,grad_log_pdf.(likelihood,y,nodes))
    Ed²logpdf = dot(inference.weights,hessian_log_pdf.(likelihood,y,nodes))
    if inference.clipping != 0
        return (abs(Edlogpdf) > inference.clipping ? sign(Edlogpdf)*inference.clipping : -Edlogpdf::T,
                abs(Ed²logpdf) > inference.clipping ? sign(Ed²logpdf)*inference.clipping : -Ed²logpdf::T)
    else
        return -Edlogpdf::T, Ed²logpdf::T
    end
end


text:
using JuMP

include(joinpath("src", "EMP-model", "model.jl"))
import .GlnModel

prev_sol=Dict{Symbol, Float64}() # a variable to store the previous solution

nadh_nad_ratio = 0.1

gln_model = GlnModel.gln_model(
    glc_ext = 50e-3,
    lac_ext = 1e-4,
    nh4_ext = 10e-3,
    ac_ext = 1e-16,
    etoh_ext = 1e-4,
    atp_adp_ratio = 10,
    nadh_nad_ratio = nadh_nad_ratio,
    num_ms = 10,
    silence=false,
)

solve_1 = GlnModel.max_mu!(gln_model, prev_sol)

solve_2 = GlnModel.find_nearest!(gln_model, prev_sol) # find a solution close to the previous one, smoothen things

# update old previous solution
GlnModel.update!(gln_model, prev_sol)


text:
module SPICE

include("../deps/deps.jl")

using SIUnits
using DataFrames
using DataArrays

export read_ascii_table

# read_ascii_table(filename, colnames)
#
# Reads SPICE data written to an ASCII file using the spice command
# `set file-type=ascii`. Unfortunately this data format does not include
# metadata but just the raw data dump, so you'll have to manually specify
# the column names. This function will also apply units based on the
# first character of the specified column names Volts for "V", Amperes for
# "I" and Seconds for "t". In `ngspice` this format is produced by the
# `wrdata` command
#

function read_ascii_table(f,cnames; sweepcol = "in")
    a=Array(Float64,0)
    c=0
    open(f) do fh
        for l in EachLine(fh)
            if c== 0
                c=length(split(l,' ',-1,false))-1
            end
            for x in split(l,' ',-1,false)
                x=="\n" && continue
                push!(a,parsefloat(x))
            end
        end
    end
    a = reshape(a,(c,itrunc(length(a)/c)))'
    d = DataFrame(a)
    issweep = true
    to_delete = ASCIIString[]
    for col in Base.names(d)
        if issweep
            push!(to_delete,col)
            issweep = false
        else
            #Keep this column, next one is sweep
            issweep = true
        end
    end
    sweep = d["x1"]
    delete!(d,to_delete)
    rename!(d,Base.names(d),[cnames])
    d[sweepcol] = sweep
    d = unitify(d,[c[1] == 'I' ? Ampere : c[1] == 't' ? Second  : Volt for c in Base.names(d)]);
end

function unitify(d,units)
    names = Base.names(d)
    ret = DataFrame()
    for i = 1:length(names)
        ret[names[i]] = d[names[i]]*units[i]
    end
    ret
end

function beginswith_ci(line,query)
    if endof(query) > endof(line)
        return false
    end
    lowercase(line[1:(prevind(line,endof(query)+1))]) == query
end

immutable TokIterator
    s::String
    offset::Int
end

import Base: start, next, done

start(t::TokIterator) = t.offset
function next(t::TokIterator,pos)
    e = endof(t.s)
    paren = 0
    orgpos = pos
    while done(t.s,pos)
        lastpos = pos
        c, pos = next(t.s,pos)
        if c == '('
            paren += 1
        elseif c == ')'
            paren -= 1
        elseif c == ',' && paren < 1
            return t.s[orgpos:lastpos], pos
        elseif isspace(c)
            if orgpos == lastpos #leading whitespace
                orgpos = pos
            else
                return t.s[orgpos:lastpos], skipwspace(s,lastpos)
            end
        end
    end
    t.s[orgpos:pos], pos
end
done(t::TokIterator,pos) = done(t.s,pos)

function skipwspace(s,offset)
    while !done(s,offset)
        c, n = next(s,offset)
        if !isspace(c)
            break
        end
        offset = n
    end
    offset
end

function getnum(s,pos)
    orgpos = pos
    while !done(s,pos)
        c, n = next(s,pos)
        if !isdigit(c)
            break
        end
        pos = n
    end
    offset, parseint(s[orgpos:pos])
end

function atodims(s,offset)
    offset = skipwspace(s,offset)
    c, offset = next(s,offset)
    needbracket = false
    state = 0
    sep = '\0'
    data = Array(Uint,0)
    if c == '['
        offset = skipwspace(s,offset)
        needbracket = true
    end
    while state != 3
        if state == 0
            if !isdigit(c)
                push!(data,0)
            else
                offset, d = getnum(s,offset)
                push!(data,d)
            end
            state = 1
        elseif state == 1 # p after a number, looking for ',' or ']'
            if sep == '\0'
                sep = c

            if c == ']' && c == sep
                state = 2
            elseif c == ',' && c == sep
                state = 0
            end
            # skip other characters, this is the behavior of SPICE 3. I believe it
            # to be a bug in the original implementation that went unnoticed because
            # this is an error case. Either way, I will not change it for compatibility
            # reasons
        elseif state == 2 #after a ']', either at the end or looking for '['
            if c == '['
                state = 0
            else
                state = 3
            end
        end
        offset = skipwspace(s,offset)
        if done(s,offset)
            break
        end
    end
    if state == 3
        err = !needbracket
    else
        err = needbracket
    end
    if err
        error("Incomplete input!")
    end
    data
end

immutable Variable{T}
    name
    length
    data::Array{T}
end

using SIUnits

function vtype_to_unit(vtype)
    if vtype == "voltage"
        return Volts
    elseif vtype == "current"
        return Ampere
    end
end

function read_raw(io::IO)
    title = "default title"
    date = "unknown date"
    plotname = "unnamed plot"
    flags = 0
    raw_padded = true
    nvars = 0
    npoints = 0
    dimensions = (0,)
    variables = Array(Variable,0)
    while !eof(io)
        line = readline(io)
        if beginswith_ci(line,"title:")
            title = strip(line[7:end])
        elseif beginswith_ci(line,"date:")
            date = strip(line[6:end])
        elseif beginswith_ci(line,"plotname:")
            plotname = strip(line[10:end])
        elseif beginswith_ci(line,"flags:")
            for tok in TokIterator(line,7)
                if lowercase(tok) == "real"
                    flags |= VF_REAL
                elseif lowercase(tok) == "complex"
                    flags |= VF_COMPLEX
                elseif lowercase(tok) == "unpadded"
                    raw_padded = false
                elseif lowercase(tok) == "padded"
                    raw_padded = true
                else
                    error("Unknown flag $tok")
                end
            end
        elseif beginswith_ci(line,"no. variables:")
            nvars = parseint(strip(line[14:end]))
        elseif beginswith_ci(line,"no. points:")
            npoints = parseint(strip(line[11:end]))
        elseif beginswith_ci(line,"dimensions:")
            npoints == 0 && error("Mispaced Dimensions: line")
            dimensions = tuple(atodims(line[11:end])...)
            prod(dimensions) != npoints && error("dimension mismatch")
        elseif beginswith_ci(line,"command:")
            error("Unimplemented")
        elseif beginswith_ci(line,"option:")
            error("Unimplemented")
        elseif beginswith_ci(line,"variables:")
            for i = 1:nvars
                # Here we assume one line per variable. I don't know if that's always the case. If not,
                # This needs to be changed
                s = readline(io)
                it = TokIterator(line,0)
                state = start(it)
                _, state = next(it,state)
                name, state = next(it,state)
                vartype, state = next(it,state)

                if isdigit(name)
                    error("Unimplemented")
                end

                while !done(s,state)
                    tok, state
                    if beginswith_ci(tok,"min=")
                        vflags |= VF_MINGIVEN
                        parseint(tok[4:end])
                    elseif beginswith_ci(tok,"max=")
                        vflags |= VF_MAXGIVEN
                        parseint(tok[4:end])
                    elseif beginswith_ci(tok,"color=")
                        vcolor = tok[6:end]
                    elseif beginswith_ci(tok,"scale=")
                        v_scale = tok[6:end]
                    elseif beginswith_ci(tok,"grid=")
                        vgrid = tok[5:end]
                    elseif beginswith_ci(tok,"plot=")
                        vplot = tok[5:end]
                    elseif beginswith_ci(tok,"dims=")

                    else
                        error("Bad parameter")
                    end
                end
                push!(variables,Variable(name,npoints,Array(quantity(flags & VF_REAL > 0 ? Float64 : Complex128,vtype_to_unit(v_type)),0)))
            end
        elseif beginswith_ci(line,"binary:") || beginswith_ci("values:")
            c = read(io,Uint8)
            is_ascii = c == uint8('v') || c == uint8('V')

            for i=1:npoints
                if is_ascii
                    error("Unimplemented")
                else
                    for v in variables
                        if i <= v.length
                            v.data[i] = read(io,eltype(v.data))
                        elseif raw_padded
                            read(io,eltype(v.data))
                        end
                    end
                end
            end
            break
        else
            error("Malformed file")
        end
    end
end

end # module


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible ideal GLM-MHD equations

gamma = 5/3
equations = IdealGlmMhdEquations2D(gamma)

initial_condition = initial_condition_convergence_test

# Get the DG approximation space
volume_flux = (flux_central, flux_nonconservative_powell)
solver = DGSEM(polydeg=4, surface_flux=(flux_hll, flux_nonconservative_powell),
               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = (0.0      , 0.0      )
coordinates_max = (sqrt(2.0), sqrt(2.0))

trees_per_dimension = (8, 8)
mesh = P4estMesh(trees_per_dimension,
                 polydeg=3, initial_refinement_level=0,
                 coordinates_min=coordinates_min, coordinates_max=coordinates_max,
                 periodicity=true)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

cfl = 0.9
stepsize_callback = StepsizeCallback(cfl=cfl)

glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        stepsize_callback,
                        glm_speed_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary