text:
using DashBootstrapComponents, DashHtmlComponents

dropdown = dbc_dropdownmenu(
    [
        dbc_dropdownmenuitem("Header", header=true),
        dbc_dropdownmenuitem("An item"),
        dbc_dropdownmenuitem(divider=true),
        dbc_dropdownmenuitem("Active and disabled", header=true),
        dbc_dropdownmenuitem("Active item", active=true),
        dbc_dropdownmenuitem("Disabled item", disabled=true),
        dbc_dropdownmenuitem(divider=true),
        html_p(
            "You can have other content here like text if you like.",
            className="text-muted px-4 mt-4",
        ),
    ],
    label="Menu",
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
@testset "$TEST $G" begin

#tests for the concrete types Graph and DiGraph
g = DiGraph(10, 20)
@test  nv(g) == 10
@test  ne(g) == 20
@test typeof(g) == DiGraph{Int}
@test vertextype(g) == Int
@test edgetype(g) == Edge{Int}

g = Graph(10, 20)
@test  nv(g) == 10
@test  ne(g) == 20
@test typeof(g) == Graph{Int}
@test edgetype(g) == Edge{Int}

w = sprand(10,10, 0.1)
@test Graph(w) == Graph{Int}(w)
@test DiGraph(w) == DiGraph{Int}(w)

@test Graph(0) == Graph() == Graph{Int}() == Graph{Int}(0)
@test DiGraph(0) == DiGraph() == DiGraph{Int}() == DiGraph{Int}(0)

g = G(10, 20)
@test G(g) isa G
@test g == G(g)
@test g == g |> Network |> G

g = DG(10, 20)
@test DG(g) isa DG
@test g == DG(g)
@test g == g |> DiNetwork |> DG

end # testset

ext: jl


text:
asarray(x::AbstractArray) = x
asarray(x::Number) = fill(x, ())

"""
    padzeros(A, N; dim = 1)

Pad array `A` with zeros to length `N` along dimension `dim`.
"""
function padzeros(A, N; dim = 1)
    shape = size(A)
    Nold = shape[dim]
    shape = Base.setindex(shape, N - Nold, dim)
    zs = zeros(eltype(A), shape)
    return cat(A, zs; dims = dim)
end

"""
    dot!(out, M, V; dim=1)

Matrix-vector multiplication along specific dimension of array `V`, storing result in `out`.

This is equivalent to iterating over all dimensions of `V` other than `dim` and applying the
matrix-vector multiplication `M * v[..., :, ...]`, but works by reshaping the array if
necessary to take advantage of faster matrix-matrix multiplication. If `dim==1`, `dot!` is
fastest and allocation-free.
"""
function dot!(out, M, V::AbstractVector; dim=1)
    size(V, dim) == size(M, 1) || throw(DomainError(
        "Size of V along dim must be same as size of M"))
    size(out) == size(V) || throw(DomainError("Input and output arrays must have same size"))
    dim == 1 || throw(DomainError("Cannot multiply along dimension $dim of 1-D array"))
    mul!(out, M, V)
end

function dot!(out, M, V::AbstractArray{T, 2}; dim=1) where T
    size(V, dim) == size(M, 1) || throw(DomainError(
        "Size of V along dim must be same as size of M"))
    size(out) == size(V) || throw(DomainError("Input and output arrays must have same size"))
    dim <= 2 || throw(DomainError("Cannot multiply along dimension $dim of 2-D array"))
    if dim == 1
        mul!(out, M, V)
    else
        Vtmp = permutedims(V, (2, 1))
        tmp = M * Vtmp
        permutedims!(out, tmp, (2, 1))
    end
end

function dot!(out, M, V::AbstractArray; dim=1) where T
    size(V, dim) == size(M, 1) || throw(DomainError(
        "Size of V along dim must be same as size of M"))
    size(out) == size(V) || throw(DomainError("Input and output arrays must have same size"))
    dim <= ndims(V) || throw(DomainError(
        "Cannot multiply along dimension $dim of $(ndims(V))-D array"))
    if dim == 1
        idxhi = CartesianIndices(size(V)[3:end])
        _dot!(out, M, V, idxhi)
    else
        dims = collect(1:ndims(V)) # [1, 2, ..., ndims(V)]
        otherdims = filter(d -> d ≠ dim, dims) # dims but without working dimension
        #= sort dimensions by their size so largest other dimension is part of matrix-matrix
            multiplication, for speed. other dimensions are iterated over in _dot! =#
        sidcs = sortperm(collect(size(V)[otherdims]), rev=true)
        perm = (dim, otherdims[sidcs]...)
        iperm = invperm(perm) # permutation to get back to original size
        Vtmp = permutedims(V, perm)
        tmp = similar(Vtmp)
        idxhi = CartesianIndices(size(Vtmp)[3:end])
        _dot!(tmp, M, Vtmp, idxhi)
        permutedims!(out, tmp, iperm)
    end
end

function _dot!(out, M, V, idxhi)
    for hi in idxhi
        mul!(view(out, :, :, hi), M, view(V, :, :, hi))
    end
end

"""
    dimdot(v, A; dim=1)

Calculate the dot product between vector `v` and one dimension of array `A`, iterating over
all other dimensions.
"""
function dimdot(v, A; dim=1)
    dims = collect(size(A))
    dims[dim] = 1
    out = Array{eltype(A)}(undef, Tuple(dims))
    dimdot!(out, v, A; dim=dim)
    return out
end

function dimdot(v, A::Vector; dim=1)
    return dot(v, A)
end

function dimdot!(out, v, A; dim=1)
    idxlo = CartesianIndices(size(A)[1:dim-1])
    idxhi = CartesianIndices(size(A)[dim+1:end])
    _dimdot!(out, v, A, idxlo, idxhi)
end

function _dimdot!(out, v, A, idxlo, idxhi)
    for lo in idxlo
        for hi in idxhi
            out[lo, 1, hi] = dot(v, view(A, lo, :, hi))
        end
    end
end

"""
    sphbesselj_scale(n)

Return the normalization factor for the (hyper)spherical Bessel function
([`sphbesselj`](@ref)) with spherical dimension ``n``, given as
``c_n = \\sqrt{\\frac{π}{2}}`` for even ``n`` and 1 otherwise.

After:

[1] J. S. Avery, J. E. Avery. Hyperspherical Harmonics And Their Physical Applications.
    Singapore: World Scientific Publishing Company, 2017.
"""
sphbesselj_scale(n) = isodd(n) ? 1.0 : √(π / 2)

"""
    sphbesselj(p, n, x)

(Hyper)spherical Bessel function of order ``p`` and spherical dimension ``n``.
The hyperspherical Bessel function generalizes the cylindrical and spherical Bessel
functions to the ``n``-sphere (embedded in ``ℝ^{n+1}``). It is given as

```math
j_p^{n}(x) = c_n x^{-(n-1)/2} J_{p + (n-1)/2}(x),
```

where ``c_n`` is a normalization factor defined by [`sphbesselj_scale`](@ref). Note that
``n`` is not an exponent here.

It has as its special cases:
- Cylindrical Bessel function (``n=1``): ``j_p^{1}(x) = J_p(x)``
- Spherical Bessel function (``n=2``): ``j_p^{2}(x) = j_p(x)``

After:

[1] J. S. Avery, J. E. Avery. Hyperspherical Harmonics And Their Physical Applications.
    Singapore: World Scientific Publishing Company, 2017.
"""
function sphbesselj(p, n, x)
    n == 1 && return besselj(p, x)
    α = (n - 1) / 2
    Jppa = besselj(p + α, x)
    cn = sphbesselj_scale(n)
    if abs(x) ≤ sqrt(eps(real(zero(Jppa))))
        if p == 0
            J0pa = cn / gamma(α + 1) / 2^α
            return oftype(Jppa, J0pa)
        else
            return zero(Jppa)
        end
    end
    jp = cn * Jppa / x^α
    return oftype(Jppa, jp)
end

"""
    besselj_zero_init(nu, n)

Get an initial guess of the ``n``th zero of the Bessel function of order ``\\nu``.
"""
function besselj_zero_init(nu, n)
    T = float(eltype(nu))
    iszero(nu) && return T(GSL.sf_bessel_zero_J0(n))
    isone(nu) && return T(GSL.sf_bessel_zero_J1(n))
    return T(GSL.sf_bessel_zero_Jnu(nu, n))
end

# Adapted from FunctionZeros.besselj_zero, with different initialization
"""
    besselj_zero(nu, n)

Get the ``n``th zero of the Bessel function of order ``\\nu``.
"""
function besselj_zero(nu, n)
    return Roots.find_zero(x -> besselj(nu, x), besselj_zero_init(nu, n), Roots.Order2())
end

"""
    sphbesselj_zero(p, n, m)

Get the ``m``th zero of the (hyper)spherical Bessel function of order ``p`` and spherical
dimension ``n``.

See [`sphbesselj`](@ref).
"""
sphbesselj_zero(p, n, m) = besselj_zero(p + (n - 1) / 2, m)

max_stars_repo_name: chrisbrahms/Hankel

max_issues_repo_name: chrisbrahms/Hankel

max_forks_repo_name: chrisbrahms/Hankel


text:
export ricker, dampedSinusoid, sinc, spectSynth

"""

    ricker(;f::Float64=20.0, n::Int64=100, dt::Float64=0.001, nargout::Bool=false)

Generate a Ricker wavelet.

# Arguments
- `f::Float64,`    : Peak frequnecy
- `n::Int64`    : Number of points
- `dt::Float64`    : Sampling interval

# Output
- `s`::Array{Float64,1}    : Computed Ricker wavelet
- `t`::Array{Float64,1}    : Time axis

"""
function ricker(;f::Float64=20.0, n::Int64=41, dt::Float64=0.05, tpoint::Float64=100.0)
    # generate time axis
    t = timeAxis(dt, 0., n)
    # create the wavelet
    s = (1.0 .- 2π^2 * f^2 .* (t.-tpoint).^2).*exp.(-π^2 * f^2 .* (t.-tpoint).^2)

    return (s, t)
end

"""

    dampedSinusoid(;A::Float64=1.0, ω::Float64=1.0, ϕ::Float64=0.0, n::Int64=100, dt::Float64=0.001, λ::Float64=2.0)

Generate a damped sinusoid.

# Arguments
- `A::Float64`    : Amplitude of the damped sinusoid
- `n::Int64`    : Number of points
- `dt::Float64`    : Sampling interval
- `ω::Float64`    : Frequency
- `λ::Float64`    : Decay constant
- `ϕ::Float64`    : Phase angle at t=0
- `η::Float64`    : Shift in function maximum

# Output
- `s::Array{Float64,1}`    : Computed damped sinusoid
- `t::Array{Float64,1}`    : Time axis

"""
function dampedSinusoid(;A::Float64=1.0, ω::Float64=1.0, ϕ::Float64=0.0, η::Float64=0.0, n::Int64=100, dt::Float64=0.001, t0::Float64=0.0, λ::Float64=2.0)
    # generate time axis
    t = timeAxis(dt, t0, n)
    # damped sinusoidal function
    s = A*exp.(-λ*abs.(t.-η)) .* cos.(ω*t .+ ϕ)

    return (s, t)
end

"""

    sinc(;A::Float64=1.0, ω::Float64=1.0, ϕ::Float64=0.0, n::Int64=100, dt::Float64=0.001, t0=0.0)

Generate a sinc (cardinal-sin) function

# Arguments
- `A::Float64`    : Amplitude of the damped sinusoid
- `n::Int64`    : Number of points
- `dt::Float64`    : Sampling interval
- `ω::Float64`    : Frequency
- `ϕ::Float64`    : Phase angle at t=0

# Output
- `s::Array{Float64,1}`    : Computed sinc
- `t::Array{Float64,1}`    : Time axis

"""
function sinc(;A::Float64=1.0, ω::Float64=1.0, ϕ::Float64=0.0, n::Int64=100, dt::Float64=0.001, t0::Float64=0.0)
    # generate time axis
    t = timeAxis(dt, t0, n)
    # sinc function
    s = A.*Base.sinc.(ω.*t.+ϕ)

    return (s, t)
end

"""

    chirp(;c::Array{Float64,1}=collect(range(0.15, stop=15.0, length=100)), tp::Array{Float64,1}=collect(range(0.04, stop=4.0, length=100)), mintp::Float64=0.04, maxtp::Float64=4.0, dist::Float64=1000.0, n::Int64=500, dt::Float64=0.002, t0::Float64=0.0)

# Arguments
- `c::Array{Float64,1}`    : phase velocity array
- `tp::Array{Float64,1}`    : period array
- `mintp::Float64`    : minimum period
- `maxtp::Float64`    : maximum period
- `dist::Float64`    : distance
- `n::Int64`    : number of samples in arrays
- `dt::Float64`    : sampling interval
- `t0::Float64`    : start time of time vector

# Output
- `s::Array{Float64,1}`    : Computed signal
- `t::Array{Float64,1}`    : Time axis

"""
function chirp(;c::Array{Float64,1}=collect(range(0.15, stop=15.0, length=100)), tp::Array{Float64,1}=collect(range(0.04, stop=4.0, length=100)), mintp::Float64=0.04, maxtp::Float64=4.0, dist::Float64=1000.0, n::Int64=500, dt::Float64=0.002, t0::Float64=0.0)
    # generate time axis
    t = timeAxis(dt, t0, n)
    # generate a chirp
    tp_ind = findall(x->(x≤maxtp && x≥mintp), tp)
    s = [sum(cos.((2π ./ tp[tp_ind]) .* (j .- dist ./ c[tp_ind]))) for j in time]

    return (s, t)
end

"""

    spectSynth(;A::Float64=1.0, n::Int64=100, dt::Float64=0.001, t0=0.0)

Generate a synthetic by definition in the frequency domain and inverse FFT

# Arguments
- `A::Float64`    : Amplitude of the damped sinusoid
- `n::Int64`    : Number of points
- `dt::Float64`    : Sampling interval

# Output
- `s::Array{Float64,1}`    : Computed signal
- `t::Array{Float64,1}`    : Time axis

"""
function spectSynth(;A::Float64=1.0, n::Int64=100, dt::Float64=0.001, t0::Float64=0.0)
    # generate time axis
    t = timeAxis(dt, t0, n)
    # generate a signal with flat amplitude spectrum and random phase in [-π, π]
    # this is of length (n/2)+1 to get a final signal of length n
    N = convert(Int64, floor(length(t)/2)+2)
    spectrum = A.+(((rand(N).-0.5).*2π).*im)
    # inverse FFT
    s = real(irfft(spectrum, n+1))

    return (s[2:end], t)
end

"""

    timeAxis(dt::Float64, t0::Float64, n::Float64)

Create a vector of times given sample frequency, number of points, and start time.

# Arguments
- `dt::Float64`    : Sample interval
- `t0::Int64`    : Start time
- `n::Int64`    : Number of points in the time vector

# Output
- `t::Array{Float64,1}`    : Time axis

"""
function timeAxis(dt::Float64, t0::Float64, n::Int64)
    T = dt * (n-1) # end time
    # create time axis from 0 to end time at dt spacing
    t = collect(0:dt:T) .+ t0

    return t
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
@polly function kernel_ludcmp(A, b, x, y)
    n = size(A,1)

    for i = 1:n
        for j = 1:(i-1)
            w = A[i,j]
            for k = 1:(j-1)
                w -= A[i,k] * A[k,j]
            end
            A[i,j] = w / A[j,j]
        end
        for j = i:n
            w = A[i,j]
            for k = 1:(i-1)
                w -= A[i,k] * A[k,j]
            end
            A[i,j] = w
        end
    end

    for i = 1:n
        w = b[i]
        for j = 1:(i-1)
            w -= A[i,j] * y[j]
        end
        y[i] = w
    end

    for i = n:-1:1
        w = y[i]
        for j = (i+1):n
            w -= A[i,j] * x[j]
        end
        x[i] = w / A[i,i]
    end
end

let
    n = 2000

    A = zeros(Float32, n, n)
    b = zeros(Float32, n)
    x = zeros(Float32, n)
    y = zeros(Float32, n)

    fn = Float32(n)

    for i = 1:n
        x[i] = 0
        y[i] = 0
        b[i] = (i+1)/fn/2.0 + 4
    end

    for i = 1:n
        for j = 1:i
            A[i,j] = (-j % n) / n + 1
        end
        for j = i+1:n
            A[i,j] = 0
        end
        A[i,i] = 1
    end

    # Make the matrix positive semi-definite.
    # not necessary for LU, but using same code as cholesky
    B = zeros(Float32, n, n)
    for r = 1:n, s = 1:n
        B[r,s] = zero(eltype(B))
    end
    for t = 1:n
        for r = 1:n, s = 1:n
            B[r,s] += A[r,t] * A[s,t]
        end
        for r = 1:n, s = 1:n
            A[r,s] = B[r,s]
        end
    end

    SUITE["ludcmp"] = @benchmarkable kernel_ludcmp(A, b, x, y) setup = (A = copy($A); b = copy($b); x = copy($x); y = copy($y))
end


text:
@testset "@js_str macro" begin
    @testset "@js_str character escaping" begin
        @test js"\"foo\"" == JSString("\"foo\"")
        @test js"\"foo\"" == JSString("\"foo\"")
        @test js"this.\$refs" == JSString(raw"this.$refs")

        @test js"foo\$" == JSString(raw"foo$")
        @test js"\$('div.my-id')" == JSString(raw"$('div.my-id')")

        @test js"'foo\\bar'" == JSString("'foo\\bar'")

        foo = "foo"
        @test js"\\$foo" == JSString("\\\"foo\"")

        # See note about Julia and weirdness with escaping when quotes are involved.
        @test js"""console.log("\\")""" == JSString("""console.log(\"\\\")""")
        #@test js"""console.log(\"\\\")""" == JSString("""console.log(\"\\\")""")
        @test js"""foo = '\\\\'""" == JSString(raw"""foo = '\\'""")
    end

    @testset "@js_str interpolation" begin
        mystr = "foo"
        @test js"const foo = $mystr" == JSString("const foo = \"foo\"")
        @test js"const foo = $(mystr)" == JSString("const foo = \"foo\"")

        mydict = Dict(:foo => "bar")
        @test js"const myDict = $mydict" == JSString("""const myDict = {"foo":"bar"}""")
    end

    @testset "@js_str interpolation overloading" begin
        person_struct_name = gensym()
        @eval begin
            struct $(person_struct_name)
                name::String
                age::Int
            end
            JSExpr.tojs(person::$(person_struct_name)) = "$(person.name)/$(person.age)"
        end
        Person = @eval $(person_struct_name)
        travis = Person("Travis", 22)

        @test js"const travis = $travis" == JSString("const travis = \"Travis/22\"")
    end

    @testset "JSString interpolation into normal strings" begin
        myjs = js"""console.log("foo");"""
        @test "<script>$(myjs)</script>" == """<script>console.log("foo");</script>"""
    end
end


text:
# ## Scaling

function scaleprojection(ratios::NTuple{N}, T = Float32) where N
    a = zeros(Float32, N, N)
    a[I(N)] = SVector{N}(ratios)
    return LinearMap(SArray{Tuple{N, N}}(a))
end

struct ScaleRatio{N} <: ProjectiveTransform
    ratios::NTuple{N}
end

function getprojection(scale::ScaleRatio, bounds; randstate = nothing)
    return scaleprojection(scale.ratios)
end

"""
    ScaleKeepAspect(minlengths) <: ProjectiveTransform

Scales the shortest side of `item` to `minlengths`, keeping the
original aspect ratio.

## Examples

{cell=ScaleKeepAspect}
```julia
using DataAugmentation, TestImages
image = testimage("lighthouse")
tfm = ScaleKeepAspect((200, 200))
apply(tfm, Image(image)) |> showitems
```
"""
struct ScaleKeepAspect{N} <: ProjectiveTransform
    minlengths::NTuple{N, Int}
end

function getprojection(scale::ScaleKeepAspect{N}, bounds; randstate = nothing) where N
    # If no scaling needs to be done, return a noop transform
    scale.minlengths == length.(bounds.rs) && return IdentityTransformation()

    # Offset `minlengths` by 1 to avoid black border on one side
    ratio = maximum((scale.minlengths .+ 1) ./ length.(bounds.rs))
    upperleft = SVector{N, Float32}(minimum.(bounds.rs)) .- 0.5
    P = scaleprojection(Tuple(ratio for _ in 1:N))
    if upperleft != SVector(0, 0)
        P = P ∘ Translation((Float32.(P(upperleft)) .+ 0.5f0))
    end
    return P
end

function projectionbounds(tfm::ScaleKeepAspect{N}, P, bounds::Bounds{N}; randstate = nothing) where N
    origsz = length.(bounds.rs)
    ratio = maximum((tfm.minlengths) ./ origsz)
    sz = floor.(Int,ratio .* origsz)
    bounds_ = transformbounds(bounds, P)
    bs_ = offsetcropbounds(sz, bounds_, ntuple(_ -> 0.5, N))
    return bs_
end

"""
    ScaleFixed(sizes)

Projective transformation that scales sides to `sizes`, disregarding
aspect ratio.

See also [`ScaleKeepAspect`](#).
"""
struct ScaleFixed{N} <: ProjectiveTransform
    sizes::NTuple{N, Int}
end

function getprojection(scale::ScaleFixed, bounds; randstate = nothing)
    ratios = (scale.sizes .+ 1) ./ length.(bounds.rs)
    upperleft = SVector{2, Float32}(minimum.(bounds.rs)) .- 1
    P = scaleprojection(ratios)
    if upperleft != SVector(0, 0)
        P = P  ∘ Translation(-upperleft)
    end
    return P
end

function projectionbounds(tfm::ScaleFixed{N}, P, bounds::Bounds{N}; randstate = nothing) where N
    bounds_ = transformbounds(bounds, P)
    return offsetcropbounds(tfm.sizes, bounds_, (1., 1.))
end

"""
    Zoom(scales = (1, 1.2)) <: ProjectiveTransform
    Zoom(distribution)

Zoom into an item by a factor chosen from the interval `scales`
or `distribution`.
"""
struct Zoom{D<:Sampleable} <: ProjectiveTransform
    dist::D
end

Zoom(scales::NTuple{2, T} = (1., 1.2)) where T = Zoom(Uniform(scales[1], scales[2]))

getrandstate(tfm::Zoom) = rand(tfm.dist)

function getprojection(tfm::Zoom, bounds::Bounds{N}; randstate = getrandstate(tfm)) where N
    ratio = randstate
    return scaleprojection(ntuple(_ -> ratio, N))
end

"""
    Rotate(γ)
    Rotate(γs)

Rotate 2D spatial data around the center by an angle chosen at
uniformly from [-γ, γ], an angle given in degrees.

You can also pass any `Distributions.Sampleable` from which the
angle is selected.

## Examples

```julia
tfm = Rotate(10)
```

"""
struct Rotate{S<:Sampleable} <: ProjectiveTransform
    dist::S
end
Rotate(γ) = Rotate(Uniform(-abs(γ), abs(γ)))

getrandstate(tfm::Rotate) = rand(tfm.dist)

function getprojection(
        tfm::Rotate,
        bounds::Bounds{2};
        randstate = getrandstate(tfm))
    γ = randstate
    middlepoint = SVector{2, Float32}(mean.(bounds.rs))
    r = γ / 360 * 2pi
    return recenter(RotMatrix(convert(Float32, r)), middlepoint)
end

"""
    Reflect(γ)
    Reflect(distribution)

Reflect 2D spatial data around the center by an angle chosen at
uniformly from [-γ, γ], an angle given in degrees.

You can also pass any `Distributions.Sampleable` from which the
angle is selected.

## Examples

```julia
tfm = Reflect(10)
```
"""
struct Reflect <: ProjectiveTransform
    γ
end

function getprojection(tfm::Reflect, bounds; randstate = getrandstate(tfm))
    r = tfm.γ / 360 * 2pi
    return centered(LinearMap(reflectionmatrix(r)), bounds)
end

"""
    centered(P, bounds)

Transform `P` so that is applied around the center of `bounds`
instead of the origin
"""
function centered(P, bounds::Bounds{2})
    upperleft = minimum.(bounds.rs)
    bottomright = maximum.(bounds.rs)

    midpoint = SVector{2, Float32}((bottomright .- upperleft) ./ 2) .+ SVector{2, Float32}(.5, .5)
    return recenter(P, midpoint)
end

FlipX() = Reflect(180)
FlipY() = Reflect(90)

function reflectionmatrix(r)
    A = SMatrix{2, 2, Float32}(cos(2r), sin(2r), sin(2r), -cos(2r))
    return round.(A; digits = 12)
end

"""
    PinOrigin()

Projective transformation that translates the data so that
the upper left bounding corner is at the origin `(0, 0)` (or
the multidimensional equivalent).

Projective transformations on images return `OffsetArray`s,
but not on keypoints. Hardware like GPUs do not support OffsetArrays,
so they will be unwrapped and no longer match up with the keypoints.

Pinning the data to the origin makes sure that the resulting
`OffsetArray` has the same indices as a regular array, starting
at one.
"""
struct PinOrigin <: ProjectiveTransform end

function getprojection(::PinOrigin, bounds; randstate = nothing)
    p = (-SVector{2, Float32}(minimum.(bounds.rs))) .+ 1
    P = Translation(p)
    return P
end

function apply(::PinOrigin, item::Union{<:Image, <:MaskMulti, <:MaskBinary}; randstate = nothing)
    item = @set item.data = parent(itemdata(item))
    item = @set item.bounds = Bounds(size(itemdata(item)))
    return item
end

function apply!(buf::AbstractItem, ::PinOrigin, item::Union{Image, MaskMulti, MaskBinary}; randstate = nothing)
    item = @set item.data = parent(itemdata(item))
    copyitemdata!(buf, item)
    return buf
end

# `PinOrigin` should not compose with a cropped transform otherwise the pinning won't work.
# This overwrites the default composition.

compose(cropped::CroppedProjectiveTransform, pin::PinOrigin) = Sequence(cropped, pin)
compose(cropped::ComposedProjectiveTransform, pin::PinOrigin) = Sequence(cropped, pin)
compose(cropped::ProjectiveTransform, pin::PinOrigin) = Sequence(cropped, pin)

# ## Resize crops

RandomResizeCrop(sz) = ScaleKeepAspect(sz) |> RandomCrop(sz) |> PinOrigin()
CenterResizeCrop(sz) = ScaleKeepAspect(sz) |> CenterCrop(sz) |> PinOrigin()

ResizePadDivisible(sz, by) = ScaleKeepAspect(sz) |> PadDivisible(by) |> PinOrigin()


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -2.3024137037198598 <= q <= 2.756640266423038)

                     add_NL_constraint(m, :(sigmoid(0.8745011162287222 + -0.4747792679508138*sigmoid(0.20979439694353186 + 0.2523963619940135*$(x[1]) + 0.6039440929888729*$(x[2])) + 0.06822380173369114*sigmoid(-0.192203947089018 + 0.35618645260012016*$(x[1]) + -0.9380253880003773*$(x[2]))) + sigmoid(-0.5284568968024108 + 0.8782860084347535*sigmoid(0.20979439694353186 + 0.2523963619940135*$(x[1]) + 0.6039440929888729*$(x[2])) + 0.990986450920027*sigmoid(-0.192203947089018 + 0.35618645260012016*$(x[1]) + -0.9380253880003773*$(x[2]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
GQL_GETSTATUSMESSAGES = """
query sdk_ls_statusmessages(\$id: ID!) {
    statusMessages(id: \$id) {
        requestId,
        action,
        state,
        timestamp,
        client {
            userId,
            robotId,
            sessionId
        }
    }
}
"""

GQL_GETSTATUSLATEST = """
query sdk_get_statuslatest(\$id: ID!) {
  statusLatest(id: \$id) {
    requestId,
    action,
    state,
    timestamp,
    client {
      userId,
      robotId,
      sessionId
    }
  }
}
"""


text:
module GPUInspector

# stdlibs
using Statistics
using Random
using LinearAlgebra
using Libdl
using Distributed: addprocs, rmprocs, @everywhere, workers
using Base: UUID
using Pkg: Pkg
using Logging

# external
using Requires
using Reexport
@reexport using CUDA
@reexport using ThreadPinning
using CpuId: cachesize
using HDF5: h5open
using UnicodePlots: UnicodePlots
using Glob: glob

# export BFloat16 for convenience
const BFloat16 = CUDA.BFloat16
export BFloat16

include("UnitPrefixedBytes.jl")
export UnitPrefixedBytes, B, KB, MB, GB, TB, KiB, MiB, GiB, TiB
export bytes, simplify, change_base, value
include("cuda_wrappers.jl")
export get_temperatures, get_power_usages, get_gpu_utilizations
include("utility.jl")
include("utility_unroll.jl")
export clear_gpu_memory,
    clear_all_gpus_memory,
    cublasGemmEx_wrapper!,
    cublasGemmEx_wrapper_wrapper!,
    toggle_tensorcoremath,
    hastensorcores, MultiLogger, multi_log
export get_cpusocket_temperatures, get_cpu_utilizations, get_cpu_utilization
export logspace
include("monitoring.jl")
export MonitoringResults,
    monitoring_start,
    monitoring_stop,
    plot_monitoring_results,
    savefig_monitoring_results,
    livemonitor_temperature,
    livemonitor_powerusage
include("workers.jl")
export @worker, @worker_create, @worker_killall

include("gpuinfo.jl")
export gpuinfo, gpuinfo_p2p_access, gpus
include("p2p_bandwidth.jl")
export p2p_bandwidth,
    p2p_bandwidth_all, p2p_bandwidth_bidirectional, p2p_bandwidth_bidirectional_all
include("host2device_bandwidth.jl")
export host2device_bandwidth
include("stresstest_tests.jl")
include("stresstest.jl")
include("stresstest_cpu.jl")
export stresstest, stresstest_cpu
include("peakflops_gpu.jl")
include("peakflops_gpu_matmul.jl")
include("peakflops_gpu_fmas.jl")
include("peakflops_gpu_wmmas.jl")
export peakflops_gpu,
    peakflops_gpu_fmas,
    peakflops_gpu_wmmas,
    peakflops_gpu_matmul,
    peakflops_gpu_matmul_graphs,
    peakflops_gpu_matmul_scaling
export theoretical_peakflops_gpu, theoretical_peakflops_gpu_tensorcores
include("memory_bandwidth.jl")
include("memory_bandwidth_saxpy.jl")
export memory_bandwidth,
    memory_bandwidth_saxpy,
    memory_bandwidth_scaling,
    memory_bandwidth_saxpy_scaling,
    theoretical_memory_bandwidth

include("hdf5.jl")
export save_monitoring_results, load_monitoring_results

function __init__()
    @require CairoMakie="13f3f980-e62b-5c42-98c6-ff1f3baf88f0" include("requires/cairomakie.jl")

    if CUDA.functional()
        toggle_tensorcoremath(true; verbose=false) # by default, use CUDA.FAST_MATH
    end
end

end


text:
"""
`sample_field` - sample the random potential V(x); this potential is log-normally distributed with mean m(x) and covariance

σ^2(-d^2/dx^2)^(-s/2),

where m, σ, and s are user-specified parameters.
"""
function sample_field(problem::RandomODEProblem)
    N = problem.nₓ + 1; σ = problem.σ; s = problem.s; m = problem.m

    ξ = randn(N)
    û = zeros(ComplexF64, 2*N)

    k = 1:N
    λ = @. σ^2/(π*k)^s

    @. û[2:N + 1] = 2 * N * √(λ) * √(2) * ξ

    u = exp.(m + imag.(ifft(û))[N + 2:end])
    return u
end


text:
include("prelude.jl")

@testset "functions" begin
    interact(`cat`, 1) do proc
        # Check that all reading function emit an ExpectTimeout exception
        @test_throws ExpectTimeout read(proc, UInt8)
        @test_throws ExpectTimeout readbytes!(proc, Vector{UInt8}(undef, 1))
        @test_throws ExpectTimeout readuntil(proc, '\n')
        @test_throws ExpectTimeout Base.wait_readnb(proc, 1)
        if VERSION ≤ v"1.3.0-rc1.0"
            @test_throws ExpectTimeout Base.wait_readbyte(proc, UInt8('\n'))
        else
            @test_throws ExpectTimeout Base.wait_close(proc)
        end
        @test_throws ExpectTimeout read(proc, String)
        @test_throws ExpectTimeout readline(proc)
        @test process_running(proc)
    end
end

@testset "with_timeout!" begin
    interact(`sh -c 'sleep 5 && echo test'`, 1) do proc
        @test_throws ExpectTimeout expect!(proc, "test")
        with_timeout!(proc, 5) do
            @test expect!(proc, "test") == ""
        end
    end
end

@testset "inf_timeout" begin
    interact(`sh -c 'sleep 5 && echo test'`, Inf) do proc
        @test expect!(proc, "test") == ""
    end
end


text:
export sod_shock_tube_builder

function sod_shock_tube_builder(;
    number_of_particles::Integer = 200,
    γ::Number = 1.4,
)
    exterior_boundary_positions = Bounds(0.0, 1.0)
    interior_boundary_positions = [0.5]
    mass_densities = [1.0, 0.125]
    pressures = [1.0, 0.1]

    eos = AdiabaticIdealGasEOS(γ, NORMALIZED)
    specific_energies =
        specific_energy_from.(Ref(eos), mass_densities, pressures)

    mass_distribution = PiecewiseUniformMassDistribution(
        exterior_boundary_positions,
        interior_boundary_positions,
        mass_densities,
        number_of_particles,
    )
    energy_distribution = PiecewiseUniformEnergyDistribution(
        interior_boundary_positions,
        specific_energies,
    )

    SimulationBuilder(
        mass_distribution,
        energy_distribution,
        eos = eos,
        velocity_distribution = StaticVelocityDistribution{1}(),
    )
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
# idea: use a type parameterized on functions as the stack element
# try to only allocate once
using DataStructures
import Base: convert, promote_rule, +, -, *, /, >, <, log

type RDiff <: Real
    val::Float64
    adj::Float64
    RDiff(val, adj) = new(val, adj)
end
typealias Converts Union{Float64, Int, Irrational, Rational}
RDiff(x::Converts) = RDiff(x, 0.)

convert(::Type{RDiff}, x::Converts) = RDiff(x)
promote_rule(::Type{Bool}, ::Type{RDiff}) = RDiff
promote_rule(::Type{Float64}, ::Type{RDiff}) = RDiff
promote_rule{T<:Integer}(::Type{T}, ::Type{RDiff}) = RDiff

value(x::RDiff) = x.val
adjoint(x::RDiff) = x.adj
>(x::RDiff, y::RDiff) = x.val > y.val
<(x::RDiff, y::RDiff) = x.val < y.val

type Node{F<:Function, T1<:Tuple, T2}
    op::F
    inputs::T1
    outputs::T2
end

global varstack = Stack(Node)

function grad(f, N)
    y = Array{RDiff}(N)
    function ∇f(x, storage)
        copy!(y, x)
        if isempty(varstack)
            f(y)
        end
        for node in reverse_iter(varstack)
            forward(node)
        end
        for node in varstack
            backprop(node)
        end
        for i in eachindex(storage)
            storage[i] = y[i].val
        end
    end

    return ∇f
end

function clear!(s::Stack)
    while !isempty(s)
        pop!(s)
    end
    nothing
end

backprop(n::Node) = backprop(n.op, n.inputs..., n.outputs)
forward(n::Node) = forward(n.op, n.inputs..., n.outputs)

macro diff_rule(sig, code)
    fn = _fn_of_vals(sig)
    op = sig.args[1]
    op_args = sig.args[2:end]
    argsyms = Symbol[s.args[1] for s in sig.args[2:end]]
    ex = quote
        function $op($(op_args...))
            res = map(RDiff, $fn)
            n = Node($op, tuple($(argsyms...)), res)
            push!(varstack, n)
            res
        end

        function forward(::typeof($op), $(op_args...), res)
            for (r, n) in zip(res, $fn)
                r.val = n
                r.adj = 0.
            end
            nothing
        end

        function backprop(::typeof($op), $(op_args...), res)
            $code
            nothing
        end
    end
    esc(ex)
end

_fn_of_vals(x) = x
function _fn_of_vals(ex::Expr)
    if ex.head == :(::)
        out = :(value($(ex.args[1])))
    else
        out = copy(ex)
        for i in eachindex(out.args)
            out.args[i] = _fn_of_vals(out.args[i])
        end
    end
    out
end

@diff_rule log(x::RDiff) begin
    x.adj += res.adj / x.val
end

@diff_rule +(x::RDiff, y::RDiff) begin
    x.adj += res.adj
    y.adj += res.adj
end

@diff_rule -(x::RDiff, y::RDiff) begin
    x.adj += res.adj
    y.adj += -res.adj
end

@diff_rule *(x::RDiff, y::RDiff) begin
    x.adj += res.adj * y.val
    y.adj += res.adj * x.val
end

@diff_rule /(x::RDiff, y::RDiff) begin
    x.adj += res.adj / y.val
    y.adj += -res.adj * x.val / y.val^2
end

@diff_rule -(x::RDiff) begin
    x.adj += -res.adj
end


text:
using FractionalSystems

G = fotf([1], [0], [1, 10, 35, 50, 24], [3.5, 2.8, 2.1, 1.4, 0.7])
rlocus(G)


text:
include("viral_load_infectivity_testpos.jl")

const scen_names = ["(b) Status Quo","(c1) Fortnightly concurrent PCR","(c2) Fortnightly random PCR", "(d) 3 LFDs per week","(e) 2 LFDs per week","(f) Daily LFDs","(g) Daily LFDs + PCR","(h) 3 LFDs + PCR",
    "(a) No testing"]

function scenario_1_setup(Ndays::Int)   #2 LFDs + 1 concurrent PCR
    #assume day0 is random for each person
    test_day0 = rand(1:7)
    LFD_test_days1 = collect((test_day0+3):7:Ndays)
    LFD_test_days2 = collect(test_day0:7:Ndays)
    PCR_test_days = collect(test_day0:7:Ndays)
    PCR_delays = draw_PCR_delays(length(PCR_test_days))
    PCR_result_days = PCR_test_days .+ PCR_delays
    if test_day0 > 4
        push!(LFD_test_days1,test_day0-4)
    end
    TestDays = vcat(LFD_test_days1,LFD_test_days2,PCR_test_days)
    ResultDays = vcat(LFD_test_days1,LFD_test_days2,PCR_result_days)
    TestTypes = vcat(ones(length(LFD_test_days1) + length(LFD_test_days2)),
                             zeros(length(PCR_test_days)))
    itd = sortperm(ResultDays)
    return TestDays[itd], TestTypes[itd], ResultDays[itd]
end

function scenario_2_setup(Ndays::Int)   #2 LFDs + 1 concurrent fortnightly PCR
    #assume day0 is random for each person
    test_day0 = rand(1:14)
    LFD_test_days1 = collect((test_day0+3):7:Ndays)
    LFD_test_days2 = collect(test_day0:7:Ndays)
    PCR_test_days = collect(test_day0:14:Ndays)
    PCR_delays = draw_PCR_delays(length(PCR_test_days))
    PCR_result_days = PCR_test_days .+ PCR_delays
    if test_day0 > 4
        push!(LFD_test_days1,test_day0-4)
        if test_day0 > 7
            push!(LFD_test_days1,test_day0-7)
            if test_day0 > 11
                push!(LFD_test_days1,test_day0-11)
            end
        end
    end
    TestDays = vcat(LFD_test_days1,LFD_test_days2,PCR_test_days)
    ResultDays = vcat(LFD_test_days1,LFD_test_days2,PCR_result_days)
    TestTypes = vcat(ones(length(LFD_test_days1) + length(LFD_test_days2)),
                             zeros(length(PCR_test_days)))
    itd = sortperm(ResultDays)
    return TestDays[itd], TestTypes[itd], ResultDays[itd]
end

function scenario_3_setup(Ndays::Int)  #2 LFDs + 1 random fortnightly PCR
    LFDtest_day0 = rand(1:7)
    PCRtest_day0 = rand(1:14)
    LFD_test_days1 = collect((LFDtest_day0+3):7:Ndays)
    LFD_test_days2 = collect(LFDtest_day0:7:Ndays)
    PCR_test_days = collect(PCRtest_day0:14:Ndays)
    PCR_delays = draw_PCR_delays(length(PCR_test_days))
    PCR_result_days = PCR_test_days .+ PCR_delays
    if LFDtest_day0 > 4
        push!(LFD_test_days1,LFDtest_day0-4)
    end
    TestDays = vcat(LFD_test_days1,LFD_test_days2,PCR_test_days)
    ResultDays = vcat(LFD_test_days1,LFD_test_days2,PCR_result_days)
    TestTypes = vcat(ones(length(LFD_test_days1) + length(LFD_test_days2)),
                             zeros(length(PCR_test_days)))
    itd = sortperm(ResultDays)
    return TestDays[itd], TestTypes[itd], ResultDays[itd]
end

function scenario_4_setup(Ndays::Int)  #3 LFDs
    test_day0 = rand(1:7)
    LFD_test_days1 = collect((test_day0+4):7:Ndays)
    LFD_test_days2 = collect((test_day0+2):7:Ndays)
    LFD_test_days3 = collect(test_day0:7:Ndays)
    if test_day0 > 3
        push!(LFD_test_days2,test_day0-3)
        if test_day0 > 5
            push!(LFD_test_days1,test_day0-5)
        end
    end
    TestDays = vcat(LFD_test_days1,LFD_test_days2,LFD_test_days3)
    TestTypes = ones(length(TestDays))
    itd = sortperm(TestDays)
    return TestDays[itd], TestTypes[itd], TestDays[itd] #result days same as test days
end

function scenario_5_setup(Ndays::Int)  #2 LFDs
    test_day0 = rand(1:7)
    LFD_test_days1 = collect((test_day0+3):7:Ndays)
    LFD_test_days2 = collect(test_day0:7:Ndays)
    if test_day0 > 4
        push!(LFD_test_days1,test_day0-4)
    end
    TestDays = vcat(LFD_test_days1,LFD_test_days2)
    TestTypes = ones(length(TestDays))
    itd = sortperm(TestDays)
    return TestDays[itd], TestTypes[itd], TestDays[itd] #result days same as test days
end

function scenario_6_setup(Ndays::Int)  #7 LFDs
    TestDays = collect(1:Ndays)
    TestTypes = ones(length(TestDays))
    return TestDays, TestTypes, TestDays #result days same as test days
end

function scenario_7_setup(Ndays::Int)   #Daily LFDs + 1 concurrent PCR
    #assume day0 is random for each person
    test_day0 = rand(1:7)
    LFD_test_days = collect(1:Ndays)
    PCR_test_days = collect(test_day0:7:Ndays)
    PCR_delays = draw_PCR_delays(length(PCR_test_days))
    PCR_result_days = PCR_test_days .+ PCR_delays
    TestDays = vcat(LFD_test_days,PCR_test_days)
    ResultDays = vcat(LFD_test_days,PCR_result_days)
    TestTypes = vcat(ones(length(LFD_test_days)), zeros(length(PCR_test_days)))
    itd = sortperm(ResultDays)
    return TestDays[itd], TestTypes[itd], ResultDays[itd]
end

function scenario_8_setup(Ndays::Int)   #3 LFDs + 1 concurrent PCR
    #assume day0 is random for each person
    test_day0 = rand(1:7)

    LFD_test_days1 = collect((test_day0+4):7:Ndays)
    LFD_test_days2 = collect((test_day0+2):7:Ndays)
    LFD_test_days3 = collect(test_day0:7:Ndays)
    if test_day0 > 3
        push!(LFD_test_days2,test_day0-3)
        if test_day0 > 5
            push!(LFD_test_days1,test_day0-5)
        end
    end
    LFD_test_days = vcat(LFD_test_days1,LFD_test_days2,LFD_test_days3)

    PCR_test_days = collect(test_day0:7:Ndays)
    PCR_delays = draw_PCR_delays(length(PCR_test_days))
    PCR_result_days = PCR_test_days .+ PCR_delays

    TestDays = vcat(LFD_test_days,PCR_test_days)
    ResultDays = vcat(LFD_test_days,PCR_result_days)
    TestTypes = vcat(ones(length(LFD_test_days)), zeros(length(PCR_test_days)))
    itd = sortperm(ResultDays)
    return TestDays[itd], TestTypes[itd], ResultDays[itd]
end

"""
    init_testing_random!(sim::Dict, testing_params::Dict, i_day::Int, Ndays::Int)

Function to initialise test positive profiles and test isolation probabilities

## Arguments:
`sim` = Dict generated by `init_VL_and_infectiousness(Ntot::Int, Pisol::Float64)`

`testing_params` = Dict containing testing options, must have:

    "scenario" => `String` = Choose from '2_LFD_+_concPCR', '3_LFD', '2_LFD_+_Fort_concPCR'
                                          '2_LFD_+_Fort_randPCR', '2_LFD'

    "comply_prob" => `Float64` = Chance of performing next LFD test (random)

`i_day` = Day to start simulation on (counting from 1)

`Ndays` = Total number of days to simulate (counting from 1)

## Returns:
`Array{Int64,1}` = Days when testing will occur

`Int64` = Next testing day index (from `i_day`)

## See also:
`init_VL_and_infectiousness(Ntot::Int, Pisol::Float64)`
"""
function init_testing_random!(sim::Dict, testing_params::Dict, Conf_PCR::Bool)
    #add test positivity profiles to simulation Dict, i_day is start day, Ndays is total length of sim
    #testing params contains "tperiod" (days between test)
    #testing params contains "protocol" (LFD_mass_protocol or PCR_mass_protocol)
    #optional: sens_rel: relative sensitivity of LFD as VL -> infinity
    Ndays = length.(sim["VL_profiles"])
    sim["will_isolate_with_test"] = ones(Bool,sim["Ntot"])
    sim["testing_paused"] = zeros(Bool,sim["Ntot"])
    sim["resume_testing"] = -ones(Int64,sim["Ntot"])
    sim["test_days"] = Array{Array{Int64,1},1}(undef,sim["Ntot"])
    sim["test_types"] = Array{Array{Int64,1},1}(undef,sim["Ntot"])
    sim["test_result_days"] = Array{Array{Int64,1},1}(undef,sim["Ntot"])
    sim["conf_PCR"] = Array{Array{Bool,1},1}(undef,sim["Ntot"])
    TestOutput = []
    if testing_params["scenario"] == scen_names[1]
        TestOutput = scenario_1_setup.(Ndays)
    elseif testing_params["scenario"] == scen_names[2]
        TestOutput = scenario_2_setup.(Ndays)
    elseif testing_params["scenario"] == scen_names[3]
        TestOutput = scenario_3_setup.(Ndays)
    elseif testing_params["scenario"] == scen_names[4]
        TestOutput = scenario_4_setup.(Ndays)
    elseif testing_params["scenario"] == scen_names[5]
        TestOutput = scenario_5_setup.(Ndays)
    elseif testing_params["scenario"] == scen_names[6]
        TestOutput = scenario_6_setup.(Ndays)
    elseif testing_params["scenario"] == scen_names[7]
        TestOutput = scenario_7_setup.(Ndays)
    elseif testing_params["scenario"] == scen_names[8]
        TestOutput = scenario_8_setup.(Ndays)
    end
    for i in 1:sim["Ntot"]
        sim["test_days"][i] = TestOutput[i][1]
        sim["test_types"][i] = TestOutput[i][2]
        sim["test_result_days"][i] = TestOutput[i][3]
        sim["conf_PCR"][i] = zeros(Bool,length(sim["test_days"][i]))
        if Conf_PCR
            sim["conf_PCR"][i][sim["test_types"][i] .== 1] .= true
        end
    end

    sim["test_pos_prob"] = get_test_probs.(sim["VL_profiles"], sim["test_days"],
        sim["test_types"]; Do_LFD_test_prob = testing_params["comply_prob"],
        all_or_nothing = testing_params["aon_compliance"])
    #assuming random compliance
end

"""
    run_testing_scenarios_impact(Ntot::Int, Pisol::Float64, LFD_comply::Float64,
                                      Conf_PCR::Bool; )

## Arguments:
`Ntot` = Number of profiles to generate
`LFD_comply` = probability of doing each LFD test
`Conf_PCR` = Do confirmatory PCR after positive LFD (if negative stop isolating early)

## Returns:
`sim_scens` = Array of Dicts, one for each scenario simulated
`sim_names` = Names of scenarios simulated
 """
function run_testing_scenarios_impact(Ntot::Int, Pisol::Float64, LFD_comply::Float64,
                            Conf_PCR::Bool; Day5release_bool::Bool = false,
                            LFD_AllorNone::Bool = false)
    sim_baseline = init_VL_and_infectiousness(Ntot, Pisol)
    Nscens = length(scen_names)
    sim_scens = Array{Dict,1}(undef,Nscens)
    for i in 1:(Nscens-1)
        sim_scens[i] = copy(sim_baseline)
        init_testing_random!(sim_scens[i], Dict("scenario"=>scen_names[i],
                    "comply_prob"=>LFD_comply, "aon_compliance"=>LFD_AllorNone), Conf_PCR)

        sim_scens[i]["inf_profile_isolation"] = copy.(sim_scens[i]["infection_profiles"])
        sim_scens[i]["isol_days"] = run_testing_scenario!.(sim_scens[i]["inf_profile_isolation"],
            sim_scens[i]["infection_profiles"],  sim_scens[i]["test_pos_prob"], sim_scens[i]["test_result_days"],
            sim_scens[i]["symp_day"] .+ 1,  sim_scens[i]["will_isolate"], sim_scens[i]["VL_profiles"],
            sim_scens[i]["conf_PCR"]; Day5release=Day5release_bool)
    end
    sim_baseline["test_pos_prob"] = fill(zeros(0),sim_baseline["Ntot"])
    sim_baseline["test_result_days"]  = fill(zeros(Int64,0),sim_baseline["Ntot"])
    Conf_PCR_h = fill(zeros(Bool,0),sim_baseline["Ntot"])
    sim_baseline["inf_profile_isolation"] = copy.(sim_baseline["infection_profiles"])
    sim_baseline["isol_days"] = run_testing_scenario!.(sim_baseline["inf_profile_isolation"],
        sim_baseline["infection_profiles"], sim_baseline["test_pos_prob"],
        sim_baseline["test_result_days"], sim_baseline["symp_day"] .+ 1,
        sim_baseline["will_isolate"], sim_baseline["VL_profiles"],
        Conf_PCR_h; Day5release=Day5release_bool)
    sim_scens[Nscens] = copy(sim_baseline)
    for i in 1:(Nscens)
        sim_scens[i]["inf_days"] = zeros(Int64, Ntot)
        for j in 1:Ntot
            sim_scens[i]["inf_days"][j] = sum((sim_scens[i]["inf_profile_isolation"][j] .> 0))
        end
    end

    return sim_scens, scen_names
end

function get_no_testing_scenario(Ntot::Int, Pisol::Float64;
        Day5release_bool::Bool = false)
    sim_baseline = init_VL_and_infectiousness(Ntot, Pisol)
    sim_baseline["test_pos_prob"] = fill(zeros(0),sim_baseline["Ntot"])
    sim_baseline["test_result_days"]  = fill(zeros(Int64,0),sim_baseline["Ntot"])
    Conf_PCR_h = fill(zeros(Bool,0),sim_baseline["Ntot"])
    sim_baseline["inf_profile_isolation"] = copy.(sim_baseline["infection_profiles"])
    sim_baseline["isol_days"] = run_testing_scenario!.(sim_baseline["inf_profile_isolation"],
        sim_baseline["infection_profiles"], sim_baseline["test_pos_prob"],
        sim_baseline["test_result_days"], sim_baseline["symp_day"] .+ 1,
        sim_baseline["will_isolate"], sim_baseline["VL_profiles"],
        Conf_PCR_h; Day5release=Day5release_bool)

    return sim_baseline
end

function run_testing_scenarios_vs_baseline(sim_baseline::Dict, LFD_comply::Float64,
                                      Conf_PCR::Bool; Day5release_bool::Bool = false,
                                      LFD_AllorNone::Bool = false)
    Ntot = sim_baseline["Ntot"]
    Nscens = length(scen_names)
    sim_scens = Array{Dict,1}(undef,Nscens)
    for i in 1:(Nscens-1)
        sim_scens[i] = copy(sim_baseline)
        init_testing_random!(sim_scens[i], Dict("scenario"=>scen_names[i],
                 "comply_prob"=>LFD_comply, "aon_compliance"=>LFD_AllorNone), Conf_PCR)

        sim_scens[i]["inf_profile_isolation"] = copy.(sim_scens[i]["infection_profiles"])
        sim_scens[i]["isol_days"] = run_testing_scenario!.(sim_scens[i]["inf_profile_isolation"],
            sim_scens[i]["infection_profiles"],  sim_scens[i]["test_pos_prob"], sim_scens[i]["test_result_days"],
            sim_scens[i]["symp_day"] .+ 1,  sim_scens[i]["will_isolate"], sim_scens[i]["VL_profiles"],
            sim_scens[i]["conf_PCR"]; Day5release=Day5release_bool)
    end
    sim_scens[Nscens] = copy(sim_baseline)
    for i in 1:(Nscens)
        sim_scens[i]["inf_days"] = zeros(Int64, Ntot)
        for j in 1:Ntot
            sim_scens[i]["inf_days"][j] = sum((sim_scens[i]["inf_profile_isolation"][j] .> 0))
        end
    end

    return sim_scens, scen_names
end


text:
include(joinpath(pwd(), "src/models/double_integrator.jl"))

# Horizon
T = 11

# Bounds
tf = 1.0
h0 = tf / (T-1) # timestep

# Initial and final states
x1 = [1.0; 1.0]
xT = [0.0; 0.0]

xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# Objective
obj = quadratic_tracking_objective(
        [t < T ? Diagonal(ones(model.n)) : Diagonal(ones(model.n)) for t = 1:T],
        [Diagonal(1.0e-1 * ones(model.m)) for t = 1:T-1],
        [xT for t = 1:T], [zeros(model.m) for t = 1:T])

# Problem
prob = trajectory_optimization_problem(model,
			   obj,
			   T,
               xl = xl,
               xu = xu,
               )

# Initialization
X0 = linear_interp(x1, xT, T) # linear interpolation for states
U0 = [ones(model.m) for t = 1:T-1]

# Pack trajectories into vector
Z0 = pack(X0, U0, prob)

# Solve
@time Z̄ = solve(prob, copy(Z0))

# Visualize
using Plots
X̄, Ū = unpack(Z̄, prob)
plot(hcat(X̄...)', width = 2.0)
plot(hcat(Ū...)', width = 2.0, linetype = :steppost)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
function lintpkg(pkg::AbstractString)
    p = joinpath(Pkg.dir(pkg), "src", basename(pkg) * ".jl")
    if !ispath(p)
        throw("cannot find path: " * p)
    end
    LintResult(lintpkgforfile(p))
end

"""
Lint the package for a file.
If file is in src, lint Package.jl and its includes.
If file is in test, lint runtests.jl and its includes.
If file is in base lint all files in base dir.
"""
function lintpkgforfile(path::AbstractString, ctx::LintContext=LintContext())
    path = abspath(path)
    if ispath(ctx.path)
        if is_windows()
            len = count(x -> x == '\\', path)
        else
            len = count(x -> x == '/', path) - 1
        end
        for _ = 1:len
            path, folder = splitdir(path)
            if folder == "src"
                file = joinpath(path, folder, basename(path) * ".jl")
                if ispath(file)
                    lintinclude(ctx, file)
                end
                break
            elseif folder == "test"
                file = joinpath(path, folder, "runtests.jl")
                if ispath(file)
                    lintinclude(ctx, file)
                end
                break
            elseif folder == "base"
                lintdir(joinpath(path, folder), ctx)
                break
            end
        end
    end
    finish(ctx)
    ctx.messages
end

function lintfile(f::AbstractString)
    if !ispath(f)
        throw("no such file exists")
    end
    str = open(readstring, f)
    lintfile(f, str)
end

function lintfile(f::AbstractString, code::AbstractString)
    ctx = LintContext(f)

    _lintstr(code, ctx)
    finish(ctx)
    msgs = ctx.messages

    # If we have an undeclared symbol, lint the package to try and resolve
    for message in msgs
        if message.code == :E321 # undeclared symbol
            ctx = LintContext(f)
            message = lintpkgforfile(file(ctx), ctx)
            break
        end
    end

    filter!(msg -> file(msg) == f, msgs)
    sort!(msgs)

    LintResult(msgs)
end

function _lintstr(str::AbstractString, ctx::LintContext, lineoffset = 0)
    linecharc = cumsum(map(x->endof(x)+1, split(str, "\n", keep=true)))
    numlines = length(linecharc)
    i = start(str)
    while !done(str,i)
        problem = false
        ex = nothing
        linerange = searchsorted(linecharc, i)
        if linerange.start > numlines # why is it not donw?
            break
        else
            linebreakloc = linecharc[linerange.start]
        end
        if linebreakloc == i || isempty(strip(str[i:(linebreakloc-1)]))# empty line
            i = linebreakloc + 1
            continue
        end
        ctx.line = ctx.lineabs = linerange.start + lineoffset
        try
            (ex, i) = parse(str,i)
        catch y
            if typeof(y) != ParseError || y.msg != "end of input"
                msg(ctx, :E111, string(y))
            end
            break
        end
        lintexpr(ex, ctx)
    end
end

"""
    lintstr(s::AbstractString)

Unlike other lint functions, like `lintpkg` and `lintfile`, `lintstr` does not
ignore any messages by default.
"""
function lintstr(str::AbstractString,
                 ctx::LintContext = (c = LintContext(); c.ignore = []; c))
    _lintstr(str, ctx)
    finish(ctx)
    ctx.messages
end

"""
Lint all .jl ending files at a given directory.
Will ignore LintContext file and already included files.
"""
function lintdir(dir::AbstractString, ctx::LintContext=LintContext())
    for file in readdir(dir)
        if endswith(file, ".jl")
            file = joinpath(dir, file)
            if !isdir(file)
                lintinclude(ctx, file)
            end
        end
    end
    finish(ctx)
    ctx.messages
end


text:
function JuMP.add_variable(
    model::RobustModel,
    v::RobustVariable,
    name::String,
    )

    model.nextvaridx+=1
    print("Index of new Variable will be ",model.nextvaridx ,"\n")

    if typeof(v) == UncertainVariable
        vref = UncertainVariableRef(model, model.nextvaridx, length(model.uncertainVariables)+1,'U')
        model.uncertainVariables[vref.type_idx] = v
        model.var_index_to_type_index[model.nextvaridx]=length(model.uncertainVariables)
        model.type_index_to_var_index[length(model.uncertainVariables)]=model.nextvaridx

    elseif typeof(v) == DecisionVariable
        vref = DecisionVariableRef(model, model.nextvaridx, length(model.decisionVariables)+1,'D')
        model.decisionVariables[vref.type_idx] = v
        model.var_index_to_type_index[model.nextvaridx]=length(model.decisionVariables)
        model.type_index_to_var_index[length(model.decisionVariables)]=model.nextvaridx
    end
    print("Now Size of uncertainVariables ",length(model.uncertainVariables),"\n")
    print("Now Size of decisionVariables ",length(model.decisionVariables),"\n")
    JuMP.set_name(vref, name)
    return vref
end

function JuMP.build_variable(
    _error::Function,
    info::JuMP.VariableInfo,
    var_type::RobustVariableType;
    extra_kwargs...
    )
    # check for unneeded keywords
    for (kwarg, _) in extra_kwargs
        _error("Keyword argument $kwarg is not for use with infinite variables.")
    end
    # get the parameter_refs
    thisStage = var_type.stage

    #ToDO
    if typeof(var_type) == Uncertain
        return UncertainVariable(info, thisStage)
    elseif typeof(var_type) == Decision
        thisPolicy = var_type.policy
        thisDependencies = var_type.depends_on
        print(thisDependencies)
        return DecisionVariable(info, thisStage,thisPolicy,thisDependencies)
    else
        _error("Weird Variable Type!")
    end
end

function Base.:(==)(v::RobustVariableRef, w::RobustVariableRef)
    return v.model === w.model && v.idx == w.idx
end

Base.broadcastable(v::RobustVariableRef) = Ref(v)
JuMP.isequal_canonical(v::RobustVariableRef, w::RobustVariableRef) = v == w

function JuMP.delete(model::RobustModel, vref::RobustVariableRef)
    @assert JuMP.is_valid(model, vref)
    if typeof(vref) == UncertainVariableRef
        delete!(model.uncertainVariables, vref.type_idx)
    elseif typeof(vref) == DecisionVariableRef
        delete!(model.decisionVariables, vref.type_idx)
    end
    return delete!(model.var_to_name, vref.idx)
end

function JuMP.delete(model::RobustModel, vrefs::Vector{RobustVariableRef})
    return JuMP.delete.(model, vrefs)
end
function JuMP.is_valid(model::RobustModel, vref::RobustVariableRef)
    if typeof(vref) == UncertainVariableRef
        return (model === vref.model && vref.type_idx in keys(model.uncertainVariables))
    else
        return (model === vref.model && vref.type_idx in keys(model.decisionVariables))
    end
end

variable_info(vref::UncertainVariableRef) = vref.model.uncertainVariables[vref.type_idx].info
variable_info(vref::DecisionVariableRef) = vref.model.decisionVariables[vref.type_idx].info

function update_variable_info(vref::RobustVariableRef, info::JuMP.VariableInfo)
    if typeof(vref) == UncertainVariableRef
        return vref.model.uncertainVariables[vref.type_idx] = JuMP.ScalarVariable(info)
    else
        return vref.model.decisionVariables[vref.type_idx] = JuMP.ScalarVariable(info)
    end
end

JuMP.has_lower_bound(vref::RobustVariableRef) = variable_info(vref).has_lb

function JuMP.lower_bound(vref::RobustVariableRef)::Float64
    @assert !JuMP.is_fixed(vref)
    return variable_info(vref).lower_bound
end

function JuMP.set_lower_bound(vref::RobustVariableRef, lower)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            true,
            lower,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end

function JuMP.delete_lower_bound(vref::RobustVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            false,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end

JuMP.has_upper_bound(vref::RobustVariableRef) = variable_info(vref).has_ub
function JuMP.upper_bound(vref::RobustVariableRef)::Float64
    @assert !JuMP.is_fixed(vref)
    return variable_info(vref).upper_bound
end
function JuMP.set_upper_bound(vref::RobustVariableRef, upper)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            true,
            upper,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
function JuMP.delete_upper_bound(vref::RobustVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            false,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end

JuMP.is_fixed(vref::RobustVariableRef) = variable_info(vref).has_fix
function JuMP.fix_value(vref::RobustVariableRef)::Float64
    return variable_info(vref).fixed_value
end
function JuMP.fix(vref::RobustVariableRef, value; force::Bool = false)
    info = variable_info(vref)
    if !force && (info.has_lb || info.has_ub)
        error(
            "Unable to fix $(vref) to $(value) because it has existing bounds.",
        )
    end
    update_variable_info(
        vref,
        JuMP.VariableInfo(
            false,
            info.lower_bound,
            false,
            info.upper_bound,
            true,
            value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
    return
end
function JuMP.unfix(vref::RobustVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            false,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            info.integer,
        ),
    )
end
function JuMP.start_value(vref::RobustVariableRef)::Union{Nothing,Float64}
    return variable_info(vref).start
end
function JuMP.set_start_value(vref::RobustVariableRef, start)
    info = variable_info(vref)
    # checking ot of bounds
    if info.has_ub && start > info.upper_bound
        error(
            "Unable to set start_value of $(vref) to $(start) because it is out of upper bounds",
        )
    end
    if info.has_lb && start < info.lower_bound
        error(
            "Unable to set start_value of $(vref) to $(start) because it is out of lower bounds",
        )
    end
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            true,
            start,
            info.binary,
            info.integer,
        ),
    )
end
JuMP.is_binary(vref::RobustVariableRef) = variable_info(vref).binary
function JuMP.set_binary(vref::RobustVariableRef)
    @assert !JuMP.is_integer(vref)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            true,
            info.integer,
        ),
    )
end
function JuMP.unset_binary(vref::RobustVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            false,
            info.integer,
        ),
    )
end
JuMP.is_integer(vref::RobustVariableRef) = variable_info(vref).integer
function JuMP.set_integer(vref::RobustVariableRef)
    @assert !JuMP.is_binary(vref)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            true,
        ),
    )
end
function JuMP.unset_integer(vref::RobustVariableRef)
    info = variable_info(vref)
    return update_variable_info(
        vref,
        JuMP.VariableInfo(
            info.has_lb,
            info.lower_bound,
            info.has_ub,
            info.upper_bound,
            info.has_fix,
            info.fixed_value,
            info.has_start,
            info.start,
            info.binary,
            false,
        ),
    )
end

function has_depends()

end


text:
#using SparseArrays
using MathProgBase
using NLPModels, NLPModelsJuMP

function test_compute_indicies()
    @testset "test_compute_indicies" begin
        J = SparseArrays.sparse([[1.0 1.0 1.0];
                    [-1.0 -1.0 -1.0];
                    [1.0 0.0 0.0];
                    [0.0 4.0 3.0];
                    [0.0 2.0 1.5];
                    [3.0 2.0 1.0];
                    [-2.0 -2.0 -2.0]])

        J_T = SparseArrays.sparse(J')

        sorted_cols = OnePhase.sorted_col_list(J_T)
        @test sorted_cols == [4, 5, 3, 6, 1, 2, 7]
        break_points = OnePhase.compute_breakpoints(J_T,sorted_cols)
        @test break_points == [1, 3, 4, 5]

        u = 1.0 ./ collect(1:size(J,1))
        no_para_indicies, para_row_info = OnePhase.compute_indicies(J,u)
        @test no_para_indicies == [1, 3, 4, 6]
        for i = 1:length(no_para_indicies)
            @test para_row_info[i].first == no_para_indicies[i]
        end
        @test para_row_info[1].ls[2].ind == 2
        @test para_row_info[1].ls[2].ratio == -1.0
        @test para_row_info[1].ls[3].ind == 7
        @test para_row_info[1].ls[3].ratio == -2.0
        @test para_row_info[3].ls[2].ind == 5
        @test para_row_info[3].ls[2].ratio == 0.5

        # rows with no duplicates
        @test para_row_info[2].u == u[3]
        @test para_row_info[2].ls[1].g == 1.0
        @test para_row_info[4].u == u[6]
        @test para_row_info[4].ls[1].g == 1.0
    end
end

function test_compute_indicies_speed()
    J = sprand(5000,5000,0.01)
    u = 1.0 ./ collect(1:size(J,1))
    tic()
    no_para_indicies, para_row_info = OnePhase.compute_indicies(J,u)
    toc()
end
#test_compute_indicies_speed()
########################################
# TODO
# test algorithm with various options
# e.g., KKT_solvers
# - schur
# - symmetric
#
#
#####

function test_kkt_solver(jump_model,pars)
    nlp_raw = OnePhase.MathOptNLPModel(jump_model)
    nlp = OnePhase.Class_CUTEst(nlp_raw)
    timer = OnePhase.class_advanced_timer()
    OnePhase.start_advanced_timer(timer)

    iter = OnePhase.mehrotra_init(nlp, pars, timer);
    OnePhase.update!(iter, timer, pars) # is this necessary ????

    #@show iter.point.x
    #@show iter.point.y
    #@show iter.point.mu
    #@show iter.point.s

    kkt_solver = OnePhase.pick_KKT_solver(pars);
    OnePhase.initialize!(kkt_solver, iter)
    OnePhase.form_system!(kkt_solver, iter, timer)
    OnePhase.factor!(kkt_solver, 1e-8, timer)
    reduct_factors = OnePhase.Reduct_affine()
    OnePhase.kkt_associate_rhs!(kkt_solver, iter, reduct_factors, timer)
    OnePhase.compute_direction!(kkt_solver, timer)

    #@show kkt_solver.dir.x
    #@show kkt_solver.dir.y
    #@show kkt_solver.dir.s

    OnePhase.pause_advanced_timer(timer)
    return kkt_solver.dir
end

function test_kkt_solvers(jump_model)
    pars = OnePhase.Class_parameters()
    pars.output_level = 0

#Original Code
    dir_schur = test_kkt_solver(jump_model,pars)
    if OnePhase.USE_HSL
        begin
            pars.kkt.kkt_solver_type=:symmetric
            pars.kkt.linear_solver_type=:HSL

            dir_sym = test_kkt_solver(jump_model,pars)

            @test LinearAlgebra.norm(dir_schur.x - dir_sym.x,2)<1e-6
            @test LinearAlgebra.norm(dir_schur.y - dir_sym.y,2)<1e-6
            @test LinearAlgebra.norm(dir_schur.s - dir_sym.s,2)<1e-6
        end
    end

#Updated Code

    begin
        pars.kkt.kkt_solver_type=:symmetric
        pars.kkt.linear_solver_type=:julia

        dir_sym = test_kkt_solver(jump_model,pars)

        @test LinearAlgebra.norm(dir_schur.x - dir_sym.x,2)<1e-6
        @test LinearAlgebra.norm(dir_schur.y - dir_sym.y,2)<1e-6
        @test LinearAlgebra.norm(dir_schur.s - dir_sym.s,2)<1e-6
    end

#Original Code
    if OnePhase.USE_HSL
        begin
            pars.kkt.kkt_solver_type=:clever_symmetric
            pars.kkt.linear_solver_type=:HSL

            dir_clever_sym =test_kkt_solver(jump_model,pars)

            @test LinearAlgebra.norm(dir_clever_sym.x - dir_sym.x,2)<1e-6
            @test LinearAlgebra.norm(dir_clever_sym.y - dir_sym.y,2)<1e-6
            @test LinearAlgebra.norm(dir_clever_sym.s - dir_sym.s,2)<1e-6
        end
    end

#Updated Code

    begin
        pars.kkt.kkt_solver_type=:clever_symmetric
        pars.kkt.linear_solver_type=:julia

        dir_clever_sym =test_kkt_solver(jump_model,pars)
        @test LinearAlgebra.norm(dir_schur.x - dir_clever_sym.x,2)<1e-6
        @test LinearAlgebra.norm(dir_clever_sym.y - dir_schur.y,2)<1e-6
        @test LinearAlgebra.norm(dir_clever_sym.s - dir_schur.s,2)<1e-6
    end

end

function test_kkt_solvers()
    @testset "test_kkt_solvers" begin
	jump_model = toy_lp0()
        test_kkt_solvers(jump_model)

        jump_model = toy_lp1()
        test_kkt_solvers(jump_model)

        jump_model = toy_lp2()
        test_kkt_solvers(jump_model)

        jump_model = toy_lp3()
        test_kkt_solvers(jump_model)

        jump_model = toy_lp4()
        test_kkt_solvers(jump_model)

        jump_model = toy_lp5()
        test_kkt_solvers(jump_model)

        jump_model = toy_lp6()
        test_kkt_solvers(jump_model)

        jump_model = toy_lp7()
        test_kkt_solvers(jump_model)

        jump_model = toy_lp8()
        test_kkt_solvers(jump_model)

    end
end

function test_compare_columns()
    @testset "test_compare_columns" begin
        A = SparseArrays.sparse(SparseArrays.sparse([[0 10 0.0]; [1 0 1]; [1 0 0]; [2 0 0];])')
        @test OnePhase.compare_columns(A,1,2) === true
        @test OnePhase.compare_columns(A,2,1) === false
        @test OnePhase.compare_columns(A,2,3) === false
        @test OnePhase.compare_columns(A,3,2) === true
        @test OnePhase.compare_columns(A,3,1) === false
        @test OnePhase.compare_columns(A,1,3) === true
    end
end


text:
using BenchmarkTools

using Yao
using LuxurySparse

Id = IMatrix{1<<16}()
Pm = pmrand(ComplexF64, 1<<16)
Dv = Diagonal(randn(ComplexF64, 1<<16))

bench = BenchmarkGroup()
bg = bench["To Sparse"] = BenchmarkGroup()
bg["permmatrix"] = @benchmarkable SparseMatrixCSC(Pm)
bg["imatrix"] = @benchmarkable SparseMatrixCSC(Id)
bg["diagonal"] = @benchmarkable SparseMatrixCSC(Dv)

showall(run(bench, verbose=true))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "GCCBootstrap"
version = v"9.4.0"

# Collection of sources required to complete build
sources = [
    # crosstool-ng will provide the build script
    ArchiveSource("http://crosstool-ng.org/download/crosstool-ng/crosstool-ng-1.25.0_rc1.tar.bz2",
                  "8a839df71bea7b2c411447b41b917f2142482df171733ddb92d615cec2f0f43a"),

    # We provide some configs for crostool-ng
    DirectorySource("./bundled"),

    # crosstool-ng can download the files, but we'd rather download them ourselves
    FileSource("http://mirrors.kernel.org/gnu/gcc/gcc-9.4.0/gcc-9.4.0.tar.xz",
               "c95da32f440378d7751dd95533186f7fc05ceb4fb65eb5b85234e6299eb9838e"),
    FileSource("https://mirrors.kernel.org/gnu/mpfr/mpfr-4.1.0.tar.xz",
               "0c98a3f1732ff6ca4ea690552079da9c597872d30e96ec28414ee23c95558a7f"),
    FileSource("https://mirrors.kernel.org/gnu/mpc/mpc-1.2.1.tar.gz",
               "17503d2c395dfcf106b622dc142683c1199431d095367c6aacba6eec30340459"),
    FileSource("https://gcc.gnu.org/pub/gcc/infrastructure/isl-0.24.tar.bz2",
               "fcf78dd9656c10eb8cf9fbd5f59a0b6b01386205fe1934b3b287a0a1898145c0"),
    FileSource("https://mirrors.kernel.org/gnu/gmp/gmp-6.2.1.tar.xz",
               "fd4829912cddd12f84181c3451cc752be224643e87fac497b69edddadc49b4f2"),
    FileSource("http://mirrors.kernel.org/pub/linux/kernel/v4.x/linux-4.1.49.tar.xz",
               "ff2e0ea5c536650aef64447c3aaa49c1a25e8f1db4ec4f7da700d3176f512ba8"),
    FileSource("https://mirrors.kernel.org/gnu/glibc/glibc-2.19.tar.xz",
               "2d3997f588401ea095a0b27227b1d50cdfdd416236f6567b564549d3b46ea2a2"),
    FileSource("https://musl.libc.org/releases/musl-1.2.2.tar.gz",
               "9b969322012d796dc23dda27a35866034fa67d8fb67e0e2c45c913c3d43219dd"),
    FileSource("https://sourceforge.net/projects/mingw-w64/files/mingw-w64/mingw-w64-release/mingw-w64-v9.0.0.tar.bz2",
               "1929b94b402f5ff4d7d37a9fe88daa9cc55515a6134805c104d1794ae22a4181"),
    FileSource("https://github.com/madler/zlib/archive/refs/tags/v1.2.11.tar.gz",
               "629380c90a77b964d896ed37163f5c3a34f6e6d897311f1df2a7016355c45eff",
               filename="zlib-1.2.11.tar.gz"),
    FileSource("http://mirrors.kernel.org/gnu/ncurses/ncurses-6.2.tar.gz",
               "30306e0c76e0f9f1f0de987cf1c82a5c21e1ce6568b9227f7da5b71cbea86c9d"),
    FileSource("http://mirrors.kernel.org/gnu/libiconv/libiconv-1.16.tar.gz",
               "e6a1b1b589654277ee790cce3734f07876ac4ccfaecbee8afa0b649cf529cc04"),
    FileSource("http://mirrors.kernel.org/gnu/gettext/gettext-0.21.tar.xz",
               "d20fcbb537e02dcf1383197ba05bd0734ef7bf5db06bdb241eb69b7d16b73192"),
    FileSource("http://mirrors.kernel.org/gnu/binutils/binutils-2.29.1.tar.xz",
               "e7010a46969f9d3e53b650a518663f98a5dde3c3ae21b7d71e5e6803bc36b577"),
    FileSource("http://mirrors.kernel.org/gnu/binutils/binutils-2.38.tar.xz",
               "e316477a914f567eccc34d5d29785b8b0f5a10208d36bbacedcc39048ecfe024"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/crosstool-ng*/

# These tools will help us to bootstrap
apk add build-base texinfo help2man ncurses-dev

# Copy in our extra patches for all packages
for package in ${WORKSPACE}/srcdir/patches/*; do
    package="$(basename "${package}")"
    for version in ${WORKSPACE}/srcdir/patches/${package}/*; do
        version="$(basename "${version}")"
        if [ ! -d packages/${package}/${version} ]; then
            continue
        fi
        cp -v ${WORKSPACE}/srcdir/patches/${package}/${version}/* packages/${package}/${version}/
    done
done

# Disable some checks that ct-ng performs
export CT_ALLOW_BUILD_AS_ROOT_SURE=1

# Unset some things that BB automatically inserts into the environment,
# but which crosstool-ng rightfully complains about.
unset LD_LIBRARY_PATH
for TOOL in CC CXX LD AS AR FC OBJCOPY OBJDUMP RANLIB STRIP LIPO MESON NM READELF; do
    unset "${TOOL}" "BUILD_${TOOL}" "${TOOL}_BUILD" "${TOOL}_FOR_BUILD" "HOST${TOOL}"
done
PATH="$(printf "%s" "$(echo -n "${PATH}" | tr ':' '\n' | grep -v "/opt")" | tr '\n' ':')"

# Build crosstool-ng for the current host
./configure --enable-local
make -j${nproc}

# Generate the appropriate crosstool-ng config file for our current target
${WORKSPACE}/srcdir/gen_config.sh > .config
cat .config

# This takes our stripped-down config and fills out all the other options
./ct-ng upgradeconfig

# Do the actual build!
./ct-ng build

# Fix case-insensitivity problems in netfilter headers
if [[ "${target}" == *linux* ]]; then
    NF="${prefix}/${target}/sysroot/usr/include/linux/netfilter"
    for NAME in CONNMARK DSCP MARK RATEEST TCPMSS; do
        mv "${NF}/xt_${NAME}.h" "${NF}/xt_${NAME}_.h"
    done
    for NAME in ECN TTL; do
        mv "${NF}_ipv4/ipt_${NAME}.h" "${NF}_ipv4/ipt_${NAME}_.h"
    done
    mv "${NF}_ipv6/ip6t_HL.h" "${NF}_ipv6/ip6t_HL_.h"
fi

# Move licenses to the right spot
mkdir -p /tmp/GCCBootstrap
mv ${prefix}/share/licenses/* /tmp/GCCBootstrap
mv /tmp/GCCBootstrap ${prefix}/share/licenses/

[[ -f "${bindir}/${target}-gcc" ]]
"""

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

# The products that we will ensure are always built
products = Product[
    #ExecutableProduct("gcc", :gcc)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

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

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
function linear_map_stable_understandable( seq::Vector{NTuple{N,T}}, poly::Vector{NTuple{N,I}} ) where {N,T,I <: Integer}
    n = length(poly)
    np = length(seq)
    Np1 = N+1
    X = Matrix{T}(undef, np*Np1, n)
    Xm = fill(zero(T), np*Np1, n)
    for i ∈ 1:np
        polies = gen_polies(poly[i], Val{1}(), T(3)^( -prod(poly[i]) ) )
        for j ∈ 1:np
            node = seq[j]
            l = np + N*(j-1)
            for (k,(c,ex)) ∈ enumerate(polies)
                X[l+k,i] = eval_poly_term(c, node, ex)
            end
            X[j,i] = eval_poly_term(one(T), node, poly[i])
        end
        Xm[i,i] = one(T)
    end
    for i ∈ np+1:n
        polies = gen_polies(poly[i], Val{1}(), T(3)^( -prod(poly[i]) ) )
        for j ∈ 1:np
            node = seq[j]
            l = np + N*(j-1)
            for (k,(c,ex)) ∈ enumerate(polies)
                Xm[l+k,i] = X[l+k,i] = eval_poly_term(c, node, ex)
            end
            X[j,i] = eval_poly_term(one(T), node, poly[i])
        end
    end
    XtX = Xm' * X
    poly, X, Xm, XtX, inv(XtX)
end
@generated function gen_polies(c::Float32, x::NTuple{N,T}) where {N,T,T2}
    Np1 = N+1
    quote
        out = Vector{Tuple{T2,NTuple{$N,T}}}($Np1)
        @nexprs $N (i -> out[i] = (T2(c*x[i]), ( @ntuple $N j -> i == j ? x[j] - one(T) : x[j] )) )
        out[$Np1] = ( T2(c), x )
        out
    end
end
@generated function gen_polies(x::NTuple{N,T}, ::Val{k}, c::T2 = one(Float32)) where {N,T,T2,k}
    N_level = binomial(N+k-1,N-1)
    quote
        buffer = Vector{T}(undef, $N)
        for i ∈ 1:$N
            buffer[i] = x[i]
        end
        N_0 = $N
        out = Vector{Tuple{T2,NTuple{$N,T}}}(undef, $N_level)
        itert = 0
        @nloops $k i l -> begin
            1:N_{$k-l}
        end l -> begin
            N_{$k-l+1} = i_l
            c_l = buffer[i_l]
            buffer[i_l] = max(zero(T), c_l - one(T))
        end l -> begin
        buffer[i_l] = c_l
        end begin
            itert += 1
            c_temp = c
            @nexprs $k j -> begin
                c_temp *= c_j
            end
            out[itert] = ( c_temp, ( @ntuple $N j -> buffer[j] ))
        end
        out
    end
end

@generated function eval_poly_term(coef::T, node::NTuple{N,T}, exponent) where {N,T}
    quote
        iszero(coef) && return zero(T)
        @nexprs $N i -> coef *= node[i] ^ exponent[i]
        if isinf(coef)
            throw("Inf error:\nnode:$node\nexponent:$exponent.")
        end
        coef
    end
end

function create_linear_map( seq::Vector{NTuple{N,T}}, poly::Vector{NTuple{N,I}}, coefs::AbstractArray  ) where {N,T,I <: Integer}
    n = length(poly)
    np = length(seq)
    Np1 = N+1
    X = Matrix{T}(n,n)
    for i ∈ 1:n
        polies = gen_polies(coefs[i], poly[i], T)
        for j ∈ 1:np
            node = seq[j]
            l = Np1*(j-1)
            for (k,(c,ex)) ∈ enumerate(polies)
                X[l+k,i] = eval_poly_term(c, node, ex)
            end
        end
    end
    poly, X, inv(X)
end
function create_linear_map( seq::Vector{NTuple{N,T}}, poly::Vector{NTuple{N,I}}  ) where {N,T,I <: Integer}
    n = length(poly)
    np = length(seq)
    Np1 = N+1
    @assert np * Np1 == n
    X = Matrix{T}(n,n)
    for i ∈ 1:n
        polies = gen_polies(one(T), poly[i], T)
        for j ∈ 1:np
            node = seq[j]
            l = Np1*(j-1)
            for (k,(c,ex)) ∈ enumerate(polies)
                X[l+k,i] = eval_poly_term(c, node, ex)
            end
        end
    end
    poly, X, inv(X)
end

function linear_map_stable_understandable( seq::Vector{NTuple{N,T}}, poly::Vector{NTuple{N,I}} ) where {N,T,I <: Integer}
    n = length(poly)
    np = length(seq)
    Np1 = N+1
    X = Matrix{T}(undef, np*Np1, n)
    Xm = fill(zero(T), np*Np1, n)
    for i ∈ 1:np
        polies = gen_polies(poly[i], Val{1}(), T(3)^( -prod(poly[i]) ) )
        for j ∈ 1:np
            node = seq[j]
            l = np + N*(j-1)
            for (k,(c,ex)) ∈ enumerate(polies)
                X[l+k,i] = eval_poly_term(c, node, ex)
            end
            X[j,i] = eval_poly_term(one(T), node, poly[i])
        end
        Xm[i,i] = one(T)
    end
    for i ∈ np+1:n
        polies = gen_polies(poly[i], Val{1}(), T(3)^( -prod(poly[i]) ) )
        for j ∈ 1:np
            node = seq[j]
            l = np + N*(j-1)
            for (k,(c,ex)) ∈ enumerate(polies)
                Xm[l+k,i] = X[l+k,i] = eval_poly_term(c, node, ex)
            end
            X[j,i] = eval_poly_term(one(T), node, poly[i])
        end
    end
    XtX = Xm' * X
    poly, X, Xm, XtX, inv(XtX)
end

@generated function gen_polies(c::Float32, x::NTuple{N,T}) where {N,T,T2}
    Np1 = N+1
    quote
        out = Vector{Tuple{T2,NTuple{$N,T}}}($Np1)
        @nexprs $N (i -> out[i] = (T2(c*x[i]), ( @ntuple $N j -> i == j ? x[j] - one(T) : x[j] )) )
        out[$Np1] = ( T2(c), x )
        out
    end
end
@generated function gen_polies(x::NTuple{N,T}, ::Val{k}, c::T2 = one(Float32)) where {N,T,T2,k}
    N_level = binomial(N+k-1,N-1)
    quote
        buffer = Vector{T}(undef, $N)
        for i ∈ 1:$N
            buffer[i] = x[i]
        end
        N_0 = $N
        out = Vector{Tuple{T2,NTuple{$N,T}}}(undef, $N_level)
        itert = 0
        @nloops $k i l -> begin
            1:N_{$k-l}
        end l -> begin
            N_{$k-l+1} = i_l
            c_l = buffer[i_l]
            buffer[i_l] = max(zero(T), c_l - one(T))
        end l -> begin
        buffer[i_l] = c_l
        end begin
            itert += 1
            c_temp = c
            @nexprs $k j -> begin
                c_temp *= c_j
            end
            out[itert] = ( c_temp, ( @ntuple $N j -> buffer[j] ))
        end
        out
    end
end

@generated function eval_poly_term(coef::T, node::NTuple{N,T}, exponent) where {N,T}
    quote
        iszero(coef) && return zero(T)
        @nexprs $N i -> coef *= node[i] ^ exponent[i]
        coef
    end
end

function create_linear_map( seq::Vector{NTuple{N,T}}, poly::Vector{NTuple{N,I}}, coefs::AbstractArray  ) where {N,T,I <: Integer}
    n = length(poly)
    np = length(seq)
    Np1 = N+1
    X = Matrix{T}(n,n)
    for i ∈ 1:n
        polies = gen_polies(coefs[i], poly[i], T)
        for j ∈ 1:np
            node = seq[j]
            l = Np1*(j-1)
            for (k,(c,ex)) ∈ enumerate(polies)
                X[l+k,i] = eval_poly_term(c, node, ex)
            end
        end
    end
    poly, X, inv(X)
end
function create_linear_map( seq::Vector{NTuple{N,T}}, poly::Vector{NTuple{N,I}}  ) where {N,T,I <: Integer}
    n = length(poly)
    np = length(seq)
    Np1 = N+1
    @assert np * Np1 == n
    X = Matrix{T}(n,n)
    for i ∈ 1:n
        polies = gen_polies(one(T), poly[i], T)
        for j ∈ 1:np
            node = seq[j]
            l = Np1*(j-1)
            for (k,(c,ex)) ∈ enumerate(polies)
                X[l+k,i] = eval_poly_term(c, node, ex)
            end
        end
    end
    poly, X, inv(X)
end

"""
Sequence and polynomial terms do not have to be of equal length.
"""
function linear_map_stable( seq::Vector{NTuple{N,T}}, poly::Vector{NTuple{N,I}} ) where {N,T,I <: Integer}
    n = length(poly)
    np = length(seq)
    Np1 = N+1
    X = Matrix{T}(undef, np*Np1, n)
    Xm = fill(zero(T), np*Np1, n)
    for i ∈ 1:np
        polies = gen_polies(poly[i], Val{1}(), T(3)^( -prod(poly[i]) ) )
        for j ∈ 1:np
            node = seq[j]
            l = np + N*(j-1)
            for (k,(c,ex)) ∈ enumerate(polies)
                X[l+k,i] = eval_poly_term(c, node, ex)
            end
            X[j,i] = eval_poly_term(one(T), node, poly[i])
        end
        Xm[i,i] = one(T)
    end
    for i ∈ np+1:n
        polies = gen_polies(poly[i], Val{1}(), T(3)^( -prod(poly[i]) ) )
        for j ∈ 1:np
            node = seq[j]
            l = np + N*(j-1)
            for (k,(c,ex)) ∈ enumerate(polies)
                Xm[l+k,i] = X[l+k,i] = eval_poly_term(c, node, ex)
            end
            Xm[j,i] = X[j,i] = eval_poly_term(one(T), node, poly[i])
        end
    end
    XtX = Xm' * X
    poly, X, Xm, XtX, inv(XtX)
end

@generated function pcvec2(::Val{N}, seq, total_size = (N+1)*64, ::Type{I} = Int) where {N,I<:Integer}
    quote
        l = length(seq) - 1
        j_0 = l
        l_0 = 1
        s_0 = 0
        out = Vector{NTuple{$N,I}}(undef, 0)
        sizehint!(out, total_size)
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N k p -> begin
                if p == $N
                    0:seq[i_p+1]
                else
                    1+seq[i_p]:seq[i_p+1]
                end
            end begin
                push!(out, (@ntuple $N p -> I(k_p)) )
            end
        end
        out
    end
end

@generated function max_degree_poly(::Val{N}, l::Int, T = Int) where N
    quote
        j_0 = l
        s_0 = 0
        ind = 0
        out = Vector{NTuple{N,T}}(undef, binomial($N+l,$N))
        @nloops $N i p -> begin
            0:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p
            j_{$N-p+1} = l - s_{$N-p+1}
        end begin
            ind += 1
            out[ind] = @ntuple $N j -> T(i_{j})
        end
        out
    end
end

using Compat
using Compat.LinearAlgebra, Sobol, Base.Cartesian, GenericSVD, Compat.Random, SpecialFunctions
rt = pcvec2(Val(3), [-1,0,1,3,4,5,6,7,8]); length(rt) #196
sbf3 = sobol_vec(Val(3), 100, BigFloat);
res3_step = linear_map_stable_understandable(sbf3, rt);

p3 = max_degree_poly(Val(3), 8); length(p3)
res3 = linear_map_stable_understandable(sbf3, p3);

size.(res3)
cond(res3[2])
cond(res3[4])

sbf8 = sobol_vec(Val(8), 500, BigFloat);
p8 = max_degree_poly(Val(8), 6); length(p8)
res8 = linear_map_stable_understandable(sbf8, p8);

size.(res3)
cond(res3[2])
cond(res3[4])


text:
distr = SymbolDistr(Dict{Int64, Float64}({0 => 0.9, 1 => 0.1}))
distr = SymbolDistr({0 => 0.9, 1 => 0.1})

# Check that only valid symbols are created
@repeatedly begin
  s = sample(distr)
  @test in(s, [0, 1])
end

# Check that an empirical distribution can be created based on samples.
ed1 = empirical_distribution([:a, :b])
@test ed1[:a] == 0.5
@test ed1[:b] == 0.5

ed2 = empirical_distribution([:a, :b, :a, :a])
@test ed2[:a] == 0.75
@test ed2[:b] == 0.25

ed3 = empirical_distribution([1, 0, 1, 1])
@test ed3[1] == 0.75
@test ed3[0] == 0.25

ed4 = empirical_distribution([1, "2", 0, 1, 1, 0, "2", 1, 1, 0])
@test ed4[0] == (3/10)
@test ed4[1] == (5/10)
@test ed4["2"] == (2/10)

# Check that we get roughly the same distribution when sampling.
N = 2000
# Should be within 0.05 in probability...
dm = JSDivergence()
expected_bound = max(evaluate(dm, [0.9, 0.1], [0.85, 0.15]), evaluate(dm, [0.9, 0.1], [0.95, 0.05]))
@repeatedly begin
  samples = [sample(distr) for i in 1:N]
  edistr = SymbolDistr(empirical_distribution(samples))
  @test distance(edistr, distr) < expected_bound
end

# A white noise stream from a binary alphabet.
wn = WhiteNoiseStream([0, 1])
@test typeof(wn) <: SymbolStream

# Sampling from a white noise stream gives an empirical distribution close to
# the uniform one.
expected_bound = evaluate(dm, [0.5, 0.5], [0.45, 0.55])
@repeatedly begin
  samples = [sample(wn) for i in 1:N]
  ps = collect(values(empirical_distribution(samples)))
  @test evaluate(dm, [0.5, 0.5], ps) < expected_bound
end

# A white noise stream from a symbol distribution.
wn2 = WhiteNoiseStream(distr)
@test typeof(wn2) <: SymbolStream

# A white noise stream from another symbol stream.
wn3 = WhiteNoiseStream(wn)
@test typeof(wn3) <: SymbolStream

# A fixed stream returns a fixed sequence of values and then returns nothing.
fs1 = FixedStream([0,0,1])
@test typeof(fs1) <: SymbolStream
@test next(fs1) == 0
@test next(fs1) == 0
@test next(fs1) == 1
@test next(fs1) == nothing
@test next(fs1) == nothing
@test next(fs1) == nothing

# An independent stream copy of a fixed stream can never return a sequence
# which is longer than the stream it is a copy of.


text:
# Ex 1.1
println("Hello World")

# Ex 1.2
for i in 1:10
    println("Hello World")
end

# Ex 1.3
for i in 1:10
    println(i)
end

# Ex 1.4
println("What is your name?")
x = chomp(readline())
if x == "Frederik"
    println("That is a good name, ", x, ". Welcome!")
else
    println("Welcome ", x, "!")
end

# Ex 1.5
println("Enter a number")
x = parse(Float64, chomp(readline()))
println("Enter another number")
y = parse(Float64, chomp(readline()))
println(x, " + ", y, " = ", x + y)

# Ex 1.6
println("Enter a number")
x = parse(Float64, chomp(readline()))
println("Enter another number")
y = parse(Float64, chomp(readline()))

println("Enter an operator (+, -, *, /)")
a = chomp(readline())
f = getfield(Main, Symbol(a))
println(x, " ", a, " ", y, " = ", f(x,y))

# Ex 1.7
println("Enter a number")
x = parse(Int64, chomp(readline()))
println("Enter another number")
y = parse(Int64, chomp(readline()))

for i in x:y
    println(i)
end

# Ex 1.8
println("Enter a number")
x = parse(Int64, chomp(readline()))
println("Enter another number")
y = parse(Int64, chomp(readline()))

if x > y
    x, y = y, x
end

for i in x:y
    println(i)
end

# Ex 1.9
function gt_previous(val, list)
    for elem in list
        if val < elem
            return false
        end
    end
    true
end

println("Enter a number")
x = parse(Int64, chomp(readline()))
num_list = []

while gt_previous(x, num_list)
    push!(num_list, x)
    println("Enter a number")
    global x = parse(Int64, chomp(readline()))
end

# Ex 1.10
println("Enter a non-negative number")
x = parse(Int64, chomp(readline()))

function fact(n::Int64)::Int64
    if n < 0
        error("Negative values are not allowed for factorial.")
    end
    res = 1
    for i in 1:n
        res *= i
    end
    return res
end
println(fact(x))

println("Enter a non-negative number")
x = parse(Int64, chomp(readline()))

fact(n) = (n > 0) ? *(1:n...) :
            (n < 0) ? error("Only non-negative vals allowed.") : 1

println(fact(x))

# Ex 1.11
println("Enter a number")
x = parse(Int64, chomp(readline()))

#if x > 0
#    r = 1:x
#else
#    r = x:-1
#end

#sum = +(r...)

sum = +(1:abs(x)...) * sign(x)

println(sum)


text:
abstract type BridgeAlgorithm <: DEAlgorithm end
struct BridgeEuler <: BridgeAlgorithm end
struct BridgeHeun <: BridgeAlgorithm end
struct BridgeSRK <: BridgeAlgorithm end
struct BridgeR3 <: BridgeAlgorithm end
struct BridgeBS3 <: BridgeAlgorithm end


text:
function heatmap(;
                 x::AbstractVector = Float64[],
                 y::AbstractVector = Float64[],
                 color::AbstractVector = Int[])

    v = VegaVisualization(width = 800, height = 500)
    add_data!(v, x = x, y = y, group = color)

    #Get unique table name
    table = v.data[1].name

    #Use defaults as shortcut, update properties as necessary
    default_scales!(v)
    v.scales[1]._type = "ordinal"
    v.scales[2]._type = "ordinal"

    v.scales[3]._type = "linear"
    v.scales[3].range = ["#313695", "#4575b4", "#74add1", "#abd9e9", "#e0f3f8", "#ffffbf", "#fee090", "#fdae61", "#f46d43", "#d73027", "#a50026"]
    cmin,cmax = extrema(color)
    v.scales[3].domain = collect(linspace(cmin,cmax,length(v.scales[3].range)))

    #Use defaults for axes
    default_axes!(v)

    #Use default legend
    legend!(v)
    v.legends[1].values = [cmin, cmax]
    v.legends[1].title = ""

    #Marks
    v.marks = Array{VegaMark}(1)
    v.marks[1] = VegaMark(_type = "rect",
                          from = VegaMarkFrom(data = table),
                          properties = VegaMarkProperties(enter = VegaMarkPropertySet(x = VegaValueRef(scale = "x", field = "x"),
                                                                                width = VegaValueRef(scale = "x", band = true),
                                                                                y = VegaValueRef(scale = "y", field = "y"),
                                                                                height = VegaValueRef(scale = "y", band = true),
                                                                                fill = VegaValueRef(scale = "group", field = "group")
                                                                                )
                                                          )
                                               )

    return v
end


text:
using Test
using CausalDiscovery.Autumn

@testset begin
  prog = au"""
    (program
    (external (: z Int))
    (: x Int)
    (= x 3)
    (: y Int)
    (= y (initnext (+ 1 2) (/ 3 this)))
    (: map (-> (-> a b) (List a) (List b)))
    (: f (-> Int Int))
    (= f (fn (x) (+ x z)))
    (= ys (map f xs))
    (= zs (map (--> a (a + 1)) xs))
    (type alias Particle ((: x Int) (: y Int)))
    (= field (.. particle x))
    (: g (-> (Particle a) (Int)))
    (= g (fn (particle)
            (case particle
                  (=> (Particle a_) 4)
                  (=> (Dog a_ b_) 5))))
    (: o Int)
    (= o (let ((= q 3) (= d 12)) (+ q d)))
    (if (x == 3) then (y == 4) else (y == 5))
  )
  """
end

## Should parse to
  """
  external z : Int
  x : Int
  x = 3
  y : Int
  y = init 1 + 2 next 3 / this
  map : (a -> b) -> List a -> List b
  f : (Int -> Int)
  f = (x) -> (x + z)
  ys = map f xs
  zs = map a --> a + 1 xs
  type alias Particle = { x : Int, y : Int }
  field = particle.x
  g : (Particle a -> Int)
  g = (particle) -> (
    case particle of
      (Particle a_) => 4
      (Dog a_ b_) => 5)
  o : Int
  o = let
    q = 3
    d = 12
   in
     q + d
  if x == 3 then y == 4 else y == 5
  """


text:
import StaticArrays.SVector
import LinearAlgebra.norm

abstract type AbstractMeasurement end
abstract type AbstractPointMeasurement{d}<:AbstractMeasurement end

function get_coordinate(m::AbstractPointMeasurement)
    return m.coordinate
end

struct PointMeasurement{d}<:AbstractPointMeasurement{d}
    coordinate::SVector{d,Float64}
end

struct ΔδPointMeasurement{Tv,d}<:AbstractPointMeasurement{d}
    coordinate::SVector{d,Float64}
    weight_Δ::Tv
    weight_δ::Tv
end

struct Δ∇δPointMeasurement{Tv,d}<:AbstractPointMeasurement{d}
    coordinate::SVector{d,Float64}
    weight_Δ::Tv
    weight_∇::SVector{d,Float64}
    weight_δ::Tv
end

function Δ∇δPointMeasurement(in::PointMeasurement{d}) where d
    return Δ∇δPointMeasurement{Float64,d}(in.coordinate, zero(Float64), SVector{d,Float64}(zeros(Float64,d)), one(Float64))
end

function ΔδPointMeasurement(in::PointMeasurement{d}) where d
    return ΔδPointMeasurement{Float64,d}(in.coordinate, zero(Float64), one(Float64))
end

function point_measurements(x::Matrix; dims=1)
    if dims == 2
        x = x'
    elseif dims !=1
        error("keyword argumend \"dims\" should be 1 or 2")
    end
    d = size(x, 1)
    return [PointMeasurement{d}(SVector{d,Float64}(x[:, k])) for k = 1 : size(x, 2)]
end


text:
"""

```
using Omega
using Distributions
using SoftPredicates
x = 1 ~ StdNormal{Float64}()
y = 2 ~ StdNormal{Float64}()
z = pw(==ₛ, x, y)
joint_post = @joint(x,y) |ᶜ z
loglikelihood = logerr ∘ softconstraints(joint_post)
# logprior = logenergyexo(x)
logposterior(ω) = logenergyexo(ω) + loglikelihood(ω)
randsample(joint_post, 10; alg = MH, logenergy = logposterior)
```
"""
function OmegaCore.randsample(rng,
                              ΩT,
                              x,
                              n,
                              ::MHAlg;
                              state_init = auto_init(rng, x, ΩT),
                              propose_and_logratio = defpropose_and_logratio(x),
                              logenergy = auto_logenergy(x),
                              kwargs...)
  ωsamples = mh(rng, logenergy, n, state_init, propose_and_logratio; kwargs...)
  map(x ∘ OmegaCore.tagignorecondition, ωsamples)
end

function OmegaCore.omegarandsample(rng,
                                   ΩT,
                                   logenergy,
                                   n,
                                   ::MHAlg;
                                   state_init = auto_init(rng, logenergy, ΩT),
                                   propose_and_logratio = defpropose_and_logratio(logenergy),
                                   kwargs...)
  ωsamples = mh(rng, logenergy, n, state_init, propose_and_logratio; kwargs...)
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using Test

test_modules = [
    "test_bracket",
    "test_golden_section",
    "test_ccs",
    "test_hooke_jeeves",
    "test_uninformed_search",
    # "test_airport",
    "test_exhaustive_search",
    "test_random_hill"
]

for mod in test_modules
    @testset "$(mod)" begin
        include("$mod.jl")
    end
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation