text:
"Parameters for right-hand side function of fully implicit Runge-Kutta methods."
mutable struct ParametersFIRKimplicit{DT, TT, D, S, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    tab::Tableau{TT}
    Δt::TT

    t::TT
    q::Vector{DT}
    θ::Vector{DT}

    function ParametersFIRKimplicit{DT,D}(equs::ET, tab::Tableau{TT}, Δt::TT) where {DT, TT, D, ET <: NamedTuple}
        new{DT, TT, D, tab.s, ET}(equs, tab, Δt, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

@doc raw"""
Fully implicit Runge-Kutta integrator cache.

### Fields

* `q̃`: initial guess of solution
* `ṽ`: initial guess of vector field
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of solution
* `V`: internal stages of vector field
* `Θ`: implicit function of internal stages
* `F`: vector field of implicit function
"""
struct IntegratorCacheFIRKimplicit{DT,D,S} <: ODEIntegratorCache{DT,D}
    q::Vector{DT}
    v::Vector{DT}
    θ::Vector{DT}

    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    Θ::Vector{Vector{DT}}
    F::Vector{Vector{DT}}

    function IntegratorCacheFIRKimplicit{DT,D,S}() where {DT,D,S}
        Q = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        Θ = create_internal_stage_vector(DT, D, S)
        F = create_internal_stage_vector(DT, D, S)
        new(zeros(DT,D), zeros(DT,D), zeros(DT,D),
            Q, V, Θ, F)
    end
end

function IntegratorCache{ST}(params::ParametersFIRKimplicit{DT,TT,D,S}; kwargs...) where {ST,DT,TT,D,S}
    IntegratorCacheFIRKimplicit{ST,D,S}(; kwargs...)
end

@inline CacheType(ST, params::ParametersFIRKimplicit{DT,TT,D,S}) where {DT,TT,D,S} = IntegratorCacheFIRKimplicit{ST,D,S}

@doc raw"""
Fully implicit Runge-Kutta integrator solving the system
```math
\begin{aligned}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
\vartheta(Q_{n,i}) &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, f (Q_{n,j}, V_{n,j}) , &
\vartheta(q_{n+1}) &= q_{n} + h \sum \limits_{i=1}^{s} b_{i}  \, f (Q_{n,j}, V_{n,j}) .
\end{aligned}
```
"""
struct IntegratorFIRKimplicit{DT, TT, D, S, PT <: ParametersFIRKimplicit{DT,TT},
                                    ST <: NonlinearSolver{DT},
                                    IT <: InitialGuessODE{TT}} <: AbstractIntegratorRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorFIRKimplicit(params::ParametersFIRKimplicit{DT,TT,D,S}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,ST,IT}
        new{DT, TT, D, S, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorFIRKimplicit{DT,D}(equations::NamedTuple, tableau::Tableau{TT}, Δt::TT; exact_jacobian=true) where {DT,TT,D}
        # get number of stages
        S = tableau.s

        # check if tableau is fully implicit
        if get_config(:verbosity) ≥ 1
            if isexplicit(tableau)
                @warn "Initializing IntegratorFIRK with explicit tableau $(tableau.name).\nYou might want to use IntegratorERK instead."
            elseif isdiagnonallyimplicit(tableau)
                @warn "Initializing IntegratorFIRK with diagonally implicit tableau $(tableau.name).\nYou might want to use IntegratorDIRK instead."
            end
        end

        # create params
        params = ParametersFIRKimplicit{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, D*(S+1), params, caches)

        # create initial guess
        iguess = InitialGuessODE(get_config(:ig_extrapolation), equations[:v̄], Δt)

        # create integrator
        IntegratorFIRKimplicit(params, solver, iguess, caches)
    end

    # function IntegratorFIRKimplicit{DT,D}(v::Function, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
    #     IntegratorFIRKimplicit{DT,D}(NamedTuple{(:v,)}((v,)), tableau, Δt; kwargs...)
    # end

    # function IntegratorFIRKimplicit{DT,D}(v::Function, h::Function, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
    #     IntegratorFIRKimplicit{DT,D}(NamedTuple{(:v,:h)}((v,h)), tableau, Δt; kwargs...)
    # end

    function IntegratorFIRKimplicit(equation::Union{IODE{DT,TT}, LODE{DT,TT}}, tableau::Tableau{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorFIRKimplicit{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

@inline Base.ndims(::IntegratorFIRKimplicit{DT,TT,D,S}) where {DT,TT,D,S} = D

Solutions.AtomicSolution(integrator::IntegratorFIRKimplicit{DT,TT}) where {DT,TT} =
    AtomicSolutionPODE(DT, TT, ndims(integrator), get_internal_variables(integrator))

function initialize!(int::IntegratorFIRKimplicit, sol::AtomicSolutionODE)
    sol.t̄ = sol.t - timestep(int)

    equations(int)[:v̄](sol.t, sol.q, sol.v)

    initialize!(int.iguess, sol.t, sol.q, sol.v,
                            sol.t̄, sol.q̄, sol.v̄)
end

function update_params!(int::IntegratorFIRKimplicit, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
    equations(int)[:ϑ](sol.t, sol.q, sol.v, int.params.θ)
end

function initial_guess!(int::IntegratorFIRKimplicit{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                        cache::IntegratorCacheFIRKimplicit{DT}=int.caches[DT]) where {DT,TT}

    # compute initial guess for internal stages
    # for i in eachstage(int)
    #     evaluate!(int.iguess, sol.q̄, sol.v̄, sol.q, sol.v, cache.Q[i], cache.V[i], tableau(int).c[i])
    # end
    for i in eachstage(int)
        for k in eachdim(int)
            int.solver.x[ndims(int)*(i-1)+k] = 0#cache.V[i][k]
        end
    end

    # compute initial guess for solution
    # evaluate!(int.iguess, sol.q̄, sol.v̄, sol.q, sol.v, cache.q, cache.v, one(TT))
    for k in eachdim(int)
        int.solver.x[ndims(int)*nstages(int)+k] = sol.q[k]#cache.q[k]
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}},
                                        Θ::Vector{Vector{ST}}, F::Vector{Vector{ST}},
                                        q::Vector{ST}, v::Vector{ST}, θ::Vector{ST},
                                        params::ParametersFIRKimplicit{DT,TT,D,S}) where {ST,DT,TT,D,S}

    # temporary variables
    local y1::ST
    local y2::ST
    local tᵢ::TT

    # copy x to V
    for i in eachindex(V)
        for k in eachindex(V[i])
            V[i][k] = x[D*(i-1)+k]
        end
    end
    for k in eachindex(q)
        q[k] = x[D*S+k]
    end

    # compute Q = q + Δt A V, Θ = ϑ(Q), F = f(Q,V)
    for i in eachindex(Q,F,Θ)
        tᵢ = params.t + params.Δt * params.tab.c[i]
        for k in eachindex(Q[i],params.q)
            y1 = 0
            y2 = 0
            for j in eachindex(V)
                y1 += params.tab.a[i,j] * V[j][k]
                y2 += params.tab.â[i,j] * V[j][k]
            end
            Q[i][k] = params.q[k] + params.Δt * (y1 + y2)
        end
        params.equs[:ϑ](tᵢ, Q[i], V[i], Θ[i])
        params.equs[:f](tᵢ, Q[i], V[i], F[i])
    end

    # compute q̄ = q + Δt B V, Θ = ϑ(q̄)
    tᵢ = params.t + params.Δt
    params.equs[:ϑ](tᵢ, q, v, θ)
end

# Compute stages of fully implicit Runge-Kutta methods.
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersFIRKimplicit{DT,TT,D,S},
                          caches::CacheDict) where {ST,DT,TT,D,S}
    # temporary variables
    local y1::ST
    local y2::ST

    # get cache for internal stages
    cache = caches[ST]

    # compute stages from nonlinear solver solution x
    compute_stages!(x, cache.Q, cache.V, cache.Θ, cache.F, cache.q, cache.v, cache.θ, params)

    # compute b
    for i in eachindex(cache.Θ)
        for k in eachindex(cache.Θ[i], params.θ)
            y1 = 0
            y2 = 0
            for j in eachindex(cache.F)
                y1 += params.tab.a[i,j] * cache.F[j][k]
                y2 += params.tab.â[i,j] * cache.F[j][k]
            end
            b[D*(i-1)+k] = cache.Θ[i][k] - params.θ[k] - params.Δt * (y1 + y2)
        end
    end
    for k in eachindex(cache.θ, params.θ)
        y1 = 0
        y2 = 0
        for j in eachindex(cache.F)
            y1 += params.tab.b[j] * cache.F[j][k]
            y2 += params.tab.b̂[j] * cache.F[j][k]
        end
        b[D*S+k] = cache.θ[k] - params.θ[k] - params.Δt * (y1 + y2)
    end
end

function integrate_step!(int::IntegratorFIRKimplicit{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                         cache::IntegratorCacheFIRKimplicit{DT}=int.caches[DT]) where {DT,TT}

    # update nonlinear solver parameters from atomic solution
    update_params!(int, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset atomic solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector field at internal stages
    compute_stages!(int.solver.x, cache.Q, cache.V, cache.Θ, cache.F, cache.q, cache.v, cache.θ, int.params)

    # compute final update
    sol.q .= cache.q
    sol.p .= cache.θ

    # compute vector field for initial guess
    equations(int)[:v̄](sol.t, sol.q, sol.v)
    # update_vector_fields!(int.iguess, sol.t, sol.q, sol.v)
end


text:
using PlotlyJS, DataFrames, CSV

function maps1()
    marker = attr(size=[20, 30, 15, 10],
                  color=[10, 20, 40, 50],
                  cmin=0,
                  cmax=50,
                  colorscale="Greens",
                  colorbar=attr(title="Some rate",
                                ticksuffix="%",
                                showticksuffix="last"),
                  line_color="black")
    trace = scattergeo(;mode="markers", locations=["FRA", "DEU", "RUS", "ESP"],
                        marker=marker, name="Europe Data")
    layout = Layout(geo_scope="europe", geo_resolution=50, width=500, height=550,
                    margin=attr(l=0, r=0, t=10, b=0))
    plot(trace, layout)
end

function maps2()
    # read Data into dataframe
    nm = tempname()
    url = "https://raw.githubusercontent.com/plotly/datasets/master/2014_us_cities.csv"
    download(url, nm)
    df = CSV.read(nm)
    rm(nm)

    trace = scattergeo(;locationmode="USA-states",
                        lat=df[:lat],
                        lon=df[:lon],
                        hoverinfo="text",
                        text=[string(x[:name], " pop: ", x[:pop]) for x in eachrow(df)],
                        marker_size=df[:pop]/50_000,
                        marker_line_color="black", marker_line_width=2)
    geo = attr(scope="usa",
               projection_type="albers usa",
               showland=true,
               landcolor="rgb(217, 217, 217)",
               subunitwidth=1,
               countrywidth=1,
               subunitcolor="rgb(255,255,255)",
               countrycolor="rgb(255,255,255)")

    layout = Layout(;title="2014 US City Populations", showlegend=false, geo=geo)
    plot(trace, layout)
end

ext: jl


text:
function maximum_hourglass()
    # Defined 2d Matrix
    B = Matrix{Int64}(undef,6,6)

    for i in 1:1:6
        A = map((x) -> parse(Int64,x) , split(readline(), " "))
        for j in 1:1:6
            B[i,j] = A[j]
        end
    end
    max_val = -100

    for i in 2:5
        for j in 3:5
            max_val = max(max_val, B[i,j] + B[i-1,j] + B[i-1,j-1] + B[i-1,j+1] + B[i+1, j] + B[i+1, j-1] + B[i+1,j+1])
        end
    end
    println(max_val)

end

maximum_hourglass()

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 11.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
### A Pluto.jl notebook ###
# v0.19.0

using Markdown
using InteractiveUtils

# ╔═╡ bc9041de-0784-4542-b770-c3e0df68a04c
begin
	using GLMakie
end

# ╔═╡ 55429836-d535-11ec-0975-57ba934ea296
md"# m1 cairomakie demo"

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╠═55429836-d535-11ec-0975-57ba934ea296
# ╠═bc9041de-0784-4542-b770-c3e0df68a04c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_path: m1-cairomakie-1 backup 1.jl

max_stars_repo_name: phpsmarter/plutonote

max_issues_repo_path: m1-cairomakie-1 backup 1.jl

max_issues_repo_name: phpsmarter/plutonote

max_forks_repo_path: m1-cairomakie-1 backup 1.jl

max_forks_repo_name: phpsmarter/plutonote


text:
"""
This is a Julia port of the MATLAB routine ``BVN'' by Alan Genz.
It has been fairly comprehensively tested against the results of the original
in MATLAB, and runs currently close to 10x faster.

Translated by Alex Bird.
"""

function gaussCDF(z::T)::T where T <: AbstractFloat
    cdf(Normal(zero(T), one(T)), z)
end

function bvn(xl::Real, xu::Real, yl::Real, yu::Real, r::Real)
    xl, xu, yl, yu, r = promote(xl, xu, yl, yu, r)
    bvn(xl, xu, yl, yu, r)
end

function bvn(xl::T, xu::T, yl::T, yu::T, r::T) where T <: AbstractFloat
# BVN
#  A function for computing bivariate normal probabilities.
#  bvn calculates the probability that
#    xl < x < xu and yl < y < yu,
#  with correlation coefficient r.
#   p = bvn( xl, xu, yl, yu, r )
#
#
#   Author
#       Alan Genz, Department of Mathematics
#       Washington State University, Pullman, Wa 99164-3113
#       Email : alangenz@wsu.edu

    p = bvnu(xl,yl,r) - bvnu(xu,yl,r) - bvnu(xl,yu,r) + bvnu(xu,yu,r);
    return clamp(p, 0, 1)
end

function bvnu(dh::T, dk::T, r::T)::T where T <: AbstractFloat
# BVNU
#  A function for computing bivariate normal probabilities.
#  bvnu calculates the probability that x > dh and y > dk.
#    parameters
#      dh 1st lower integration limit
#      dk 2nd lower integration limit
#      r   correlation coefficient
#  Example: p = bvnu( -3, -1, .35 )
#  Note: to compute the probability that x < dh and y < dk,
#        use bvnu( -dh, -dk, r ).
#
#
#   Author
#       Alan Genz
#       Department of Mathematics
#       Washington State University
#       Pullman, Wa 99164-3113
#       Email : alangenz@wsu.edu
#
#    This function is based on the method described by
#        Drezner, Z and G.O. Wesolowsky, (1989),
#        On the computation of the bivariate normal inegral,
#        Journal of Statist. Comput. Simul. 35, pp. 101-107,
#    with major modifications for double precision, for |r| close to 1,
#    and for Matlab by Alan Genz. Minor bug modifications 7/98, 2/10.
#
#
# Copyright (C) 2013, Alan Genz,  All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided the following conditions are met:
#   1. Redistributions of source code must retain the above copyright
#      notice, this list of conditions and the following disclaimer.
#   2. Redistributions in binary form must reproduce the above copyright
#      notice, this list of conditions and the following disclaimer in
#      the documentation and/or other materials provided with the
#      distribution.
#   3. The contributor name(s) may not be used to endorse or promote
#      products derived from this software without specific prior
#      written permission.
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
# COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
# OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
# TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

    # deal with easy cases first (if one or other limit is ∞ or ⟂)
    (dh == Inf || dk == Inf) && return 0
    (dh == -Inf) && return (dk == -Inf) ? 1 : gaussCDF(-dk);
    (dk == -Inf) && return gaussCDF(-dh);
    (r == 0) && return gaussCDF(-dh)*gaussCDF(-dk);

    # now we have that the bounds are finite, and corrcoeff is !== 0.
    tp = 2*pi; h = dh; k = dk; hk = h*k; bvn = 0.0;

    # set up Gauss Legendre points / weights for integration
    if abs(r) < 0.3      # n =  6
        w = [0.1713244923791705, 0.3607615730481384, 0.4679139345726904];
        x = [0.9324695142031522, 0.6612093864662647, 0.2386191860831970];
    elseif abs(r) < 0.75 # n = 12
        w = [.04717533638651177, 0.1069393259953183, 0.1600783285433464,
             0.2031674267230659, 0.2334925365383547, 0.2491470458134029];
        x = [0.9815606342467191, 0.9041172563704750, 0.7699026741943050,
             0.5873179542866171, 0.3678314989981802, 0.1252334085114692];
    else                # n = 20
        w = [.01761400713915212, .04060142980038694, .06267204833410906,
             .08327674157670475, 0.1019301198172404, 0.1181945319615184,
             0.1316886384491766, 0.1420961093183821, 0.1491729864726037,
             0.1527533871307259];
        x = [0.9931285991850949, 0.9639719272779138, 0.9122344282513259,
             0.8391169718222188, 0.7463319064601508, 0.6360536807265150,
             0.5108670019508271, 0.3737060887154196, 0.2277858511416451,
             0.07652652113349733];
    end

    w, x = convert.(T, vcat(w, w)), convert.(T, vcat(1.0 .- x, 1.0 .+ x));

    # easier case when correlation is not too large.
    if abs(r) < 0.925
        hs  = (h^2 + k^2)/2;
        asr = asin(r)/2;
        sn  = sin.(asr*x);
        bvn = exp.((sn*hk .- hs)./(1.0 .- sn.^2))'*w;
        bvn = bvn * asr/tp + gaussCDF(-h)*gaussCDF(-k);
    else
        if r < 0
            k, hk = -k, -hk
        end
        if abs(r) < 1
            as = 1.0-r^2;
            a = sqrt(as);
            bs = (h-k)^2;
            asr = -( bs/as + hk )/2;
            c = (4.0-hk)/8.0;
            d = (12.0-hk)/80.0;
            if asr > -100;
                bvn = a*exp(asr)*(1-c*(bs-as)*(1-d*bs)/3+c*d*as^2);
            end
            if hk  > -100;
                b = sqrt(bs);
                sp = sqrt(tp)*gaussCDF(-b/a);
                bvn -= exp(-hk/2)*sp*b*( 1 - c*bs*(1-d*bs)/3 );
            end
            a = a/2;
            xs = (a*x).^2;
            asr = -( bs./xs .+ hk )/2;
            ix = findall( asr .> -100.0 );
            xs = xs[ix];
            sp = ( 1.0 .+ c*xs.*(1.0 .+ 5.0*d*xs) );
            rs = sqrt.(1.0 .- xs);
            ep = exp.( -(hk/2)*xs./(1.0 .+ rs).^2 )./rs;
            bvn = ( a*( (exp.(asr[ix]).*(sp-ep))'*w[ix] ) - bvn )/tp;
        end
        if r > 0
            bvn =  bvn + gaussCDF( -max( h, k ) );
        elseif h >= k
            bvn = -bvn;
        else
            if h < 0
                L = gaussCDF(k)-gaussCDF(h);
            else
                L = gaussCDF(-h)-gaussCDF(-k);
            end
            bvn =  L - bvn;
        end
    end
    p = clamp(bvn, 0.0, 1.0)
    return p
end

max_stars_repo_name: ornithos/ProbLineSearch


text:
@testset "autoencoder" begin

    # autoencoder cost: sum of squared error
    # function autoencoder_cost(We1, We2, Wd, b1, b2, x)
    #     firstLayer = logistic.(We1 * x .+ b1)
    #     encodedInput = logistic.(We2 * firstLayer .+ b2)
    #     reconstructedInput = logistic.(Wd * encodedInput)
    #     cost = sum((reconstructedInput .- x) .^ 2.0)
    #     return cost
    # end

    function autoencoder_cost(We1, We2, Wd, b1, b2, x)
        firstLayer = We1 * x .+ b1
        encodedInput = We2 * firstLayer .+ b2
        reconstructedInput = Wd * encodedInput
        cost = sum((reconstructedInput .- x) .^ 2.0)
        return cost
    end

    # @rdcmp autoencoder_cost We1=rand(4,5) We2=rand(3,4) Wd=rand(5,3) b1=rand(4) b2=rand(3) x=rand(5)
    test_compare(autoencoder_cost; We1=rand(4,5), We2=rand(3,4), Wd=rand(5,3),
                 b1=rand(4), b2=rand(3), x=rand(5))

    # slightly different version that uncovered a bug in previous version
    function autoencoder_cost2(We1, We2, Wd, b1, b2, x)
        firstLayer = logistic.(We1 * x .+ b1)
        encodedInput = logistic.(We2 * firstLayer .+ b2)
        reconstructedInput = logistic.(Wd * encodedInput)
        cost = sum(reconstructedInput .- x .^ 2.0)
        return cost
    end

    # @rdcmp autoencoder_cost We1=rand(4,5) We2=rand(3,4) Wd=rand(5,3) b1=rand(4) b2=rand(3) x=rand(5)
    test_compare(autoencoder_cost2; We1=rand(4,5), We2=rand(3,4), Wd=rand(5,3),
                 b1=rand(4), b2=rand(3), x=rand(5))

    # and the most simple autoencoder
    function autoencoder_cost3(We, Wd, x)
        y = We * x
        xd = Wd * y
        cost = sum(xd .- x)
        return cost
    end
    test_compare(autoencoder_cost3; We=rand(3,4), Wd=rand(4,3), x=rand(4,2))

end


text:
using Pkg
using ArgParse
using Dates

include("DataPacket.jl")
using .DataPacket

include("SensorsWebApp.jl")
using .SensorsWebApp

function parse_commandline()
    s = ArgParseSettings()

		s.prog = "Serial packet receiver"
		s.description = "Receives, prints and writes serial data."
		s.usage = "Usage: julia lorasensors.jl portname baudrate path"

    @add_arg_table! s begin
		"--path", "-p"
			help = "Path to save the data."
			arg_type = String
			required = false
			default = string(pwd(), "/", Dates.today(), "_serialdata.txt")

		"--baudrate", "-b"
			help = "Baudrate. Default = 115200."
			arg_type = Int
			required = false
			default = 115200

        "portname"
        	help = "Portname used for serial communication."
        	arg_type = String
        	required = true
    end

    return parse_args(ARGS, s)
end

function main()
	Pkg.activate(pwd())
	Pkg.instantiate()

	parsed_args = parse_commandline()
	path_to_datafile = parsed_args["path"]
	baudrate = parsed_args["baudrate"]
	portname = parsed_args["portname"]

	# Open a serial connection to the microcontroller
	mcu = open(portname, baudrate; mode=SP_MODE_READ)

	# Retrieve a first data packet to populate the datafile
	datapacket = retrieve_data(mcu)
	write_data(datapacket, path_to_datafile)
	############################################################################
	# Reference data written in datafile
	data = load_data(path_to_datafile)

	# Instantiate a Stipple's ReactiveModel:
	sensors_model = init_model(data)

	# Route to localhost:
	route("/") do
	  ui(sensors_model) |> html
	end

	# Launch the webserver and open a browser window:
	up()
	############################################################################

	while true
		# Get a new datapacket!
		datapacket = retrieve_data(mcu)
		# Feed the ReactiveModel
		update_data!(data, datapacket, sensors_model)
		# Save the data
		write_data(datapacket, path_to_datafile)
	end
end

main()

max_stars_repo_name: alflhomme/lorasensors

max_issues_repo_name: alflhomme/lorasensors

max_forks_repo_name: alflhomme/lorasensors


text:
using Base.Test
using QuantumOptics

distance(a::SparseOperator, b::SparseOperator) = tracedistance(full(a), full(b))

for spinnumber=1//2:1//2:5//2
    spinbasis = SpinBasis(spinnumber)
    I = operators.identity(spinbasis)
    Zero = SparseOperator(spinbasis)
    sx = sigmax(spinbasis)
    sy = sigmay(spinbasis)
    sz = sigmaz(spinbasis)
    sp = sigmap(spinbasis)
    sm = sigmam(spinbasis)

    # Test traces
    @test_approx_eq 0. trace(sx)
    @test_approx_eq 0. trace(sy)
    @test_approx_eq 0. trace(sz)

    # Test kommutation relations
    kommutator(x, y) = x*y - y*x

    @test_approx_eq_eps 0. distance(kommutator(sx, sx), Zero) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sx, sy), 2im*sz) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sx, sz), -2im*sy) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sy, sx), -2im*sz) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sy, sy), Zero) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sy, sz), 2im*sx) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sz, sx), 2im*sy) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sz, sy), -2im*sx) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sz, sz), Zero) 1e-12

    # Test creation and anihilation operators
    @test_approx_eq 0. distance(sp, 0.5*(sx + 1im*sy))
    @test_approx_eq 0. distance(sm, 0.5*(sx - 1im*sy))
    @test_approx_eq 0. distance(sx, (sp + sm))
    @test_approx_eq 0. distance(sy, -1im*(sp - sm))

    # Test commutation relations with creation and anihilation operators
    @test_approx_eq_eps 0. distance(kommutator(sp, sm), sz) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sz, sp), 2*sp) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sz, sm), -2*sm) 1e-12

    # Test v x (v x u) relation: [sa, [sa, sb]] = 4*(1-delta_{ab})*sb
    @test_approx_eq_eps 0. distance(kommutator(sx, kommutator(sx, sx)), Zero) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sx, kommutator(sx, sy)), 4*sy) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sx, kommutator(sx, sz)), 4*sz) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sy, kommutator(sy, sx)), 4*sx) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sy, kommutator(sy, sy)), Zero) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sy, kommutator(sy, sz)), 4*sz) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sz, kommutator(sz, sx)), 4*sx) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sz, kommutator(sz, sy)), 4*sy) 1e-12
    @test_approx_eq_eps 0. distance(kommutator(sz, kommutator(sz, sz)), Zero) 1e-12

    # Test spinup and spindown states
    @test_approx_eq_eps 1. norm(spinup(spinbasis)) 1e-11
    @test_approx_eq_eps 1. norm(spindown(spinbasis)) 1e-11
    @test_approx_eq_eps 0. norm(sp*spinup(spinbasis)) 1e-11
    @test_approx_eq_eps 0. norm(sm*spindown(spinbasis)) 1e-11
end

# Test special relations for spin 1/2

spinbasis = SpinBasis(1//2)
I = identity(spinbasis)
Zero = SparseOperator(spinbasis)
sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

# Test antikommutator
antikommutator(x, y) = x*y + y*x

@test_approx_eq 0. distance(antikommutator(sx, sx), 2*I)
@test_approx_eq 0. distance(antikommutator(sx, sy), Zero)
@test_approx_eq 0. distance(antikommutator(sx, sz), Zero)
@test_approx_eq 0. distance(antikommutator(sy, sx), Zero)
@test_approx_eq 0. distance(antikommutator(sy, sy), 2*I)
@test_approx_eq 0. distance(antikommutator(sy, sz), Zero)
@test_approx_eq 0. distance(antikommutator(sz, sx), Zero)
@test_approx_eq 0. distance(antikommutator(sz, sy), Zero)
@test_approx_eq 0. distance(antikommutator(sz, sz), 2*I)

# Test if involutory for spin 1/2
@test_approx_eq 0. distance(sx*sx, I)
@test_approx_eq 0. distance(sy*sy, I)
@test_approx_eq 0. distance(sz*sz, I)
@test_approx_eq 0. distance(-1im*sx*sy*sz, I)

# Test consistency of spin up and down with sigmap and sigmam
@test_approx_eq_eps 0. norm(sm*spinup(spinbasis) - spindown(spinbasis)) 1e-11
@test_approx_eq_eps 0. norm(sp*spindown(spinbasis) - spinup(spinbasis)) 1e-11


text:
struct ∇SemML_2 <: DiffFunction end

function (diff::∇SemML_2)(par, grad, model::Sem{O, I, L, D}) where
    {O <: SemObs, L <: Loss, I <: Imply, D <: SemAnalyticDiff}
    model.imply.imp_fun(model.imply.imp_cov, par)
    model.imply.gradient_fun(model.imply.∇Σ, par)
    a = cholesky(Symmetric(model.imply.imp_cov); check = false)
    if !isposdef(a)
        grad .= 0.0
    else
        Σ_inv = inv(a)
        grad .= ((vec(Σ_inv')'-vec((Σ_inv*model.observed.obs_cov*Σ_inv)')')*model.imply.∇Σ)'
    end
end


text:
using Flux
include("base.jl")
include("individual.jl")
include("population.jl")
include("lemonade.jl")

ind = Individual()
ind.genotype = [MConvBlock((40, 40, 4), (40, 40, 16)), MDepthwiseConvBlock((40, 40, 16), (40, 40, 32)), MSkipConnection((40, 40, 16), (40, 40, 48)), MConv((40, 40, 48), (40, 40, 128))]
println(ind.genotype)
construct_model!(ind, 4, LAYER_OUTPUT)
println(ind.model)
X = randn(Float32, 40, 40, 4, 1)
println(ind.model(X))

mutate!(ind)
println(ind.genotype)

max_stars_repo_name: mahmoud-al-najar/neuroevolution

max_issues_repo_name: mahmoud-al-najar/neuroevolution

max_forks_repo_name: mahmoud-al-najar/neuroevolution


text:
# ------------------------------------------------------------------------------------------
# # Modeling data 1
#
# Machine learning and data science is about modeling data. **Modeling** is the
# representation of an idea with some parameters and  a mathematical representation which we
# will encode in software. All machine learning methods are about training a computer to fit
# a model to some data. Even the fanciest neural networks are simply choices for models. In
# this notebook, we will begin to start building our first computational model of data.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Modeling data is hard!
#
# Let's pick up where we left off in notebook 1 with fruit. We were left with a riddle: when
# we load images of apples and bananas,
# ------------------------------------------------------------------------------------------

using Images

apple = load("data/10_100.jpg")

banana = load("data/104_100.jpg")

# ------------------------------------------------------------------------------------------
# and then compare their average value for the color red, we end up with something that is
# perhaps surprising:
# ------------------------------------------------------------------------------------------

apple_red_amount = mean(Float64.(red.(apple)))
banana_red_amount = mean(Float64.(red.(banana)));

"The average value of red in the apple is $apple_red_amount, " *
"while the average value of red in the banana is $banana_red_amount."

# ------------------------------------------------------------------------------------------
# We see that the banana's mean red value is higher than the apple's, even though the apple
# looks much redder. Can you guess why?
#
# There are actually two reasons. One of the reasons is the background: the image of the
# banana has a lot more white background than the apple, and that white background has a red
# value of 1! In our minds we ignore the background and say "the banana is bright yellow,
# the apple is dark red", but a computer just has a bundle of numbers and does not know
# where it should be looking.
#
# The other issue is that "bright yellow" isn't a color that exists in a computer. The
# computer has three colors: red, green, and blue. "Bright yellow" in a computer is a
# mixture of red and green, and it just so happens that to get this color yellow, it needs
# more red than the apple!
# ------------------------------------------------------------------------------------------

"The amount of red in the apple at (60, 60) is $(Float64(red(apple[60, 60]))), " *
"while the amount of red in the banana at (60, 60) is $(Float64(red(banana[60, 60])))."

apple[60,60]

banana[60,60]

# ------------------------------------------------------------------------------------------
# This is a clear example that modeling data is hard!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### A note on string interpolation
#
# In the last two input cells, we *interpolated a string*. This means that when we write the
# string using quotation marks (`"  "`), we insert a placeholder for some **value** we want
# the string to include. When the string is evaluated, the value we want the string to
# include replaces the placeholder. For example, in the following string,
#
# ```julia
# mystring = "The average value of red in the apple is $apple_red_amount"
# ```
#
# `$apple_red_amount` is a placeholder for the value stored in the variable
# `apple_red_amount`. Julia knows that we want to use the value bound to the variable
# `apple_red_amount` and *not* the word "apple_red_amount" because of the dollar sign, `$`,
# that comes before `apple_red_amount`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# Execute the following code to see what the dollar sign does:
#
# ```julia
# mypi = 3.14159
# println("I have a variable called mypi that has a value of $mypi.")
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 2
#
# Alter and execute the code that creates `mystring` below
#
# ```julia
# apple_blue_amount = mean(Float64.(blue.(apple)))
# mystring = "The average amount of blue in the apple is apple_blue_amount"
# ```
#
# so that `println(mystring)` prints a string that reports the mean value of blue coloration
# in our image of an apple.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Take some time to think about the data
#
# Apples and bananas are very different, but how could we use the array of RGB values (which
# is how the images are represented in the computer, as we saw in notebook 1) to tell the
# difference between the two? Here are some quick ideas:
#
# - We could use the shape of the object in the image. But how can we encode ideas about
# shape from an array?
# - We could use the size of the object in the image. But how do we calculate that size?
# - We could use another color, or combinations of colors, from the image. Which colors?
#
# Let's go with the last route. The banana is yellow, which is a combination of red and
# green, while the apple is red. This means that the color that clearly differentiates
# between the two is not red, but green!
# ------------------------------------------------------------------------------------------

apple_green_amount = mean(Float64.(green.(apple)))
banana_green_amount = mean(Float64.(green.(banana)));

"The average value of green in the apple is $apple_green_amount, " *
"while the average value of green in the banana is $banana_green_amount."

# ------------------------------------------------------------------------------------------
# The processes that we just went through are assigned fancy names: feature selection and
# data munging.
#
# **Feature selection** is the process of subsetting the data to a more relevant and
# informative set. We took the full image data and decided to select out the green channel.
#
# **Data munging** is transforming the data into a format more suitable for modeling. Here,
# instead of keeping the full green channel, we transformed it down to a single data point:
# the average amount of green.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Building a model
#
# We want to model the connection between "the average amount of green" and "is an apple or
# banana".
#
#
#
# This model is a mathematical function which takes in our data and spits out a number that
# we will interpret as "is an apple" or "is a banana".
#
#
#
#
# We will interpret the output of the function as "is an apple" if the output is close to 0,
# and "is a banana" if it's close to 1. Anything in the middle is something we are unsure
# about. Here we're using a mathematical function to perform a **classification**.
#
# Knowing how to declare and work with functions will allow us to model our data in the
# coming sections, so this is the subject of the next notebook!
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/03. ML - Representing data with models.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# https://docs.julialang.org/en/v1/base/parallel/

using Test

function learn_Tasks()
    ## Core.Task
    # Task from a function
    a() = sum([i for i in 1:1000])
    b = Task(a)

    ## Base.@task
    b = @task a()

    @test false == istaskstarted(b)
    schedule(b)
    yield()
    @test true == istaskstarted(b)
    @test true == istaskdone(b)

    x = fetch(b)
    @show x
end

learn_Tasks()

nothing


text:
"""
    MongoStreamFile(filepath;
        [flags=JL_O_RDONLY], [mode=0],
        [timeout_msec=DEFAULT_TIMEOUT_MSEC],
        [chunk_size=DEFAULT_CHUNK_SIZE])

Creates a stream from file located at `filepath`.

`flags` is the input to pass to `open`. Must be one of the constants
defined at `Base.Filesystem`: `Base.Filesystem.JL_O_RDONLY`,
`Base.Filesystem.JL_O_CREAT`, etc.

`mode` is an optional mode to pass to `open`.
"""
function MongoStreamFile(filepath::AbstractString;
    flags::Integer=Base.Filesystem.JL_O_RDONLY,
    mode::Integer=0,
    timeout_msec::Integer=DEFAULT_TIMEOUT_MSEC,
    chunk_size::Integer=DEFAULT_CHUNK_SIZE)

    handle = mongoc_stream_file_new_for_path(filepath, flags, mode)
    if handle == C_NULL
        error("Couldn't create stream for $filepath.")
    end
    return MongoStreamFile(nothing, handle, timeout_msec=timeout_msec, chunk_size=chunk_size)
end

function check_stream_error(stream::MongoIOStream)
    err_ref = Ref{BSONError}()

    if mongoc_gridfs_bucket_stream_error(stream.handle, err_ref)
        throw(err_ref[])
    end

    nothing
end

function check_stream_error(stream::AbstractMongoStream)
    nothing
end

"""
    abort_upload(io::MongoIOStream)

Aborts the upload of a GridFS upload stream.
"""
function abort_upload(io::MongoIOStream)
    ok = mongoc_gridfs_bucket_abort_upload(io.handle)
    if !ok
        check_stream_error(io)
        error("Error aborting upload.")
    end
    nothing
end

function Base.flush(stream::AbstractMongoStream)
    ok = mongoc_stream_flush(stream.handle)
    if ok == -1
        check_stream_error(stream)
        error("Error flushing stream.")
    end
    @assert ok == 0 "Unexpected result from mongoc_stream_flush: $(Int(ok))"
    nothing
end

function Base.close(stream::AbstractMongoStream)
    if stream.isopen
        ok = mongoc_stream_close(stream.handle)
        if ok == -1
            check_stream_error(stream)
            error("Error closing stream.")
        end
        @assert ok == 0 "Unexpected result from mongoc_stream_close: $(Int(ok))"
        stream.isopen = false
    end
    nothing
end

Base.isopen(stream::AbstractMongoStream) = stream.isopen

function Base.readbytes!(s::AbstractMongoStream, buffer::AbstractArray{UInt8}, nb=length(b))

    chunk = Vector{UInt8}(undef, s.chunk_size)
    total_nr = 0

    while true
        next_nb = min(s.chunk_size, nb - total_nr)
        nr = fillbuffer!(s, chunk, next_nb)

        if nr == 0
            # reached end of stream
            break
        end

        transfer_to_buffer!( buffer, total_nr, chunk, nr )

        total_nr += nr

        # should never happen
        @assert total_nr <= nb

        if total_nr == nb
            # read the number of bytes requested
            break
        end
    end

    return total_nr
end

# returns the number of available bytes in the buffer
function nb_free(buffer, buffer_pos::Integer)
    @assert length(buffer) >= buffer_pos
    return length(buffer) - buffer_pos
end

# returns the number of bytes that will be copied to the buffer
function nb_to_fill(buffer, buffer_pos::Integer, nr::Integer)
    fr = nb_free(buffer, buffer_pos)
    return min(nr, fr)
end

# returns the number of bytes that will be appended to the buffer
function nb_to_append(buffer, buffer_pos::Integer, nr::Integer)
    result = nr - nb_to_fill(buffer, buffer_pos, nr)
    @assert result >= 0
    return result
end

# transfers nr bytes from data to target.buffer
# writes data after target.pos, increasing buffer size if needed
function transfer_to_buffer!(buffer::T, buffer_pos::Integer, data::T, nr::Integer) where {T}

    # sanity checks
    @assert nr > 0
    @assert length(data) >= nr
    @assert buffer_pos <= length(buffer)

    resize!(buffer, length(buffer) + nb_to_append(buffer, buffer_pos, nr))

    for i in 1:nr
        buffer[buffer_pos + i] = data[i]
    end
end

# reads at most nb byter to buffer without resizing it
# returns the number of bytes read
function fillbuffer!(s::AbstractMongoStream, buffer::AbstractArray{UInt8}, nb::Integer)
    @assert nb <= length(buffer)
    nr = mongoc_stream_read(s.handle, pointer(buffer), nb, 0, s.timeout_msec)

    if nr == -1
        check_stream_error(stream)
        error("Error closing stream.")
    end

    return nr
end

function Base.unsafe_write(io::AbstractMongoStream, p::Ptr{UInt8}, n::UInt)
    written = mongoc_stream_write(io.handle, p, n, io.timeout_msec)

    if written == -1
        check_stream_error(io)
        error("Error writing to stream.")
    end

    @assert written >= 0
    return Int(written)
end


text:
"""
Simple management of distributed jobs for parallel work that doesn't block the
REPL.

# Example
```julia
Jobs.launch_workers(16)

@everywhere function my_calculation(arg1, arg2)
    sleep(arg1)
    return arg2
end

Jobs.run_on_workers(
    "run1", my_calculation,
    [((10-i), i) for i in 1:10]
)

# Fetch any completed results without blocking
Jobs.current_results("run1")
Jobs.current_times("run1")  # Job runtimes
```
"""
module Jobs

using Distributed
using OrderedCollections

export launch_workers, run_on_workers,
    current_results, current_times, current_time_info

capacity = 10000
const jobs = RemoteChannel(()->Channel{Tuple{Any, Function, Any}}(capacity))
const results = RemoteChannel(()->Channel{Pair{Tuple{Any, Any}, Tuple}}(
                    capacity))
const collected = OrderedDict()

function kill_workers(; waitfor=30)
    ws = workers()
    (length(ws) == 1 && 1 in ws) && return
    rmprocs(ws, waitfor=waitfor)
end

function clear_jobs()
    while isready(jobs)
        take!(jobs)
    end
    nothing
end

function take_results()
    while isready(results)
        (job_id, job_args), (
            result, t, bytes, gctime, memallocs) = take!(results)
        collected[job_id][job_args] = (result, t, bytes, gctime, memallocs)
    end
    nothing
end

function launch_workers(n::Integer)
    num_new = n - length(setdiff!(Set(workers()), [1]))
    # Create worker processes
    new_workers = addprocs(num_new)
    length(workers()) == n || error("didn't start the right number of workers")
    1 in workers() && error("main process is listed as a worker")
    nothing
end

function run_on_workers(job_id, job_f, job_args_list)
    job_id in keys(collected) && error("job_id <$job_id> has already been used")

    1 in workers() && error("run launch_workers(n) first")

    # Ensure do_work is defined
    @everywhere function _do_work(jobs, results)
        while isready(jobs)
            job_id, job_f, job_args = take!(jobs)
            # Run the job
            val, t, bytes, gctime, memallocs = @timed job_f(job_args...)
            put!(results,
                 (job_id, job_args) => (val, t, bytes, gctime, memallocs))
        end
        #print("Worker $(myid()) is done.")
    end

    # Make jobs
    collected[job_id] = OrderedDict()
    function make_jobs()
        for job_args in job_args_list
            put!(jobs, (job_id, job_f, job_args))
        end
    end
    #@async
    make_jobs()

    # Start worker tasks
    for pid in workers()
        remote_do(Main._do_work, pid, jobs, results)
    end
    nothing
end

function current_results(job_id)
    take_results()
    OrderedDict(
         k => val
         for (k, (val,)) in collected[job_id]
    )
end

function current_times(job_id)
    take_results()
    OrderedDict(
         k => time
         for (k, (val, time)) in collected[job_id]
    )
end

function current_time_info(job_id)
    take_results()
    OrderedDict(
         k => v[2:end]
         for (k, v) in collected[job_id]
    )
end

end  # module

max_stars_repo_name: yongshanding/VLQ

max_issues_repo_name: yongshanding/VLQ

max_forks_repo_name: yongshanding/VLQ


text:
_in_interval(v, x, ε) = x - ε <= v <= x + ε

for N in [Float64, Float32, Rational{Int}]
    ε = N(1e-3)

    b1 = BallInf(zeros(N, 2), N(1))
    b2 = BallInf(ones(N, 2), N(1))

    hd = hausdorff_distance(b1, b2; p=N(Inf), ε=ε)
    @test _in_interval(hd, N(1), ε)

    hd = hausdorff_distance(b1, b1; p=N(Inf), ε=ε)
    @test _in_interval(hd, N(0), ε)
end

if test_suite_polyhedra
    for N in [Float64]
        ε = N(1e-3)

        b1 = N[1 0; 0 1] * BallInf(zeros(N, 2), N(1))
        b2 = BallInf(ones(N, 2), N(1))

        hd = hausdorff_distance(b1, b2; p=N(Inf), ε=ε)
        @test _in_interval(hd, N(1), ε)
    end
end


text:
function parserect(line::String)
    a, b, c, d = split(line, ' ')
    x, y = parse.(Int, split(c[1:end-1], ',')) .+ 1
    w, h = parse.(Int, split(d, 'x')) .- 1
    return x, y, w, h
end

function foo()
    lines = readlines(open("data/3.txt", "r"))
    claims = zeros(Int, 1000, 1000)
    for line in lines
        x, y, w, h = parserect(line)
        claims[x:x+w, y:y+h] .+= 1
    end
    for line in lines
        a, = split(line, ' ')
        x, y, w, h = parserect(line)
        if all(claims[x:x+w, y:y+h] .== 1)
            return a[2:end]
        end
    end
end

println(foo())


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

mutable struct VirtualMachineScaleSetNetworkConfigurationDnsSettings <: SwaggerModel
    dnsServers::Any # spec type: Union{ Nothing, Vector{String} } # spec name: dnsServers

    function VirtualMachineScaleSetNetworkConfigurationDnsSettings(;dnsServers=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetNetworkConfigurationDnsSettings, Symbol("dnsServers"), dnsServers)
        setfield!(o, Symbol("dnsServers"), dnsServers)
        o
    end
end # type VirtualMachineScaleSetNetworkConfigurationDnsSettings

const _property_map_VirtualMachineScaleSetNetworkConfigurationDnsSettings = Dict{Symbol,Symbol}(Symbol("dnsServers")=>Symbol("dnsServers"))
const _property_types_VirtualMachineScaleSetNetworkConfigurationDnsSettings = Dict{Symbol,String}(Symbol("dnsServers")=>"Vector{String}")
Base.propertynames(::Type{ VirtualMachineScaleSetNetworkConfigurationDnsSettings }) = collect(keys(_property_map_VirtualMachineScaleSetNetworkConfigurationDnsSettings))
Swagger.property_type(::Type{ VirtualMachineScaleSetNetworkConfigurationDnsSettings }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetNetworkConfigurationDnsSettings[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetNetworkConfigurationDnsSettings }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetNetworkConfigurationDnsSettings[property_name]

function check_required(o::VirtualMachineScaleSetNetworkConfigurationDnsSettings)
    true
end

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


text:
module MetricsProject

using LinearAlgebra, Statistics, NLsolve

using Random, Distributions
using Optim
using Optim: converged, maximum, maximizer, minimizer, iterations

# Including the modules
include("Berry1994.jl")
include("sdc.jl")

# Exporting the functions of the modules
export Berry1994, sdc

end # module


text:
abstract type Sim end

struct POMDPSim <: Sim
    simulator::Simulator
    pomdp::POMDP
    policy::Policy
    updater::Updater
    initial_belief::Any
    initial_state::Any
    metadata::Dict{Symbol}
end

problem(sim::POMDPSim) = sim.pomdp

struct MDPSim <: Sim
    simulator::Simulator
    mdp::MDP
    policy::Policy
    initial_state::Any
    metadata::Dict{Symbol}
end

problem(sim::MDPSim) = sim.mdp

"""
    Sim(p::POMDP, policy::Policy, metadata=Dict(:note=>"a note"))
    Sim(p::POMDP, policy::Policy[, updater[, initial_belief[, initial_state]]]; kwargs...)

Create a `Sim` object that represents a POMDP simulation.
"""
function Sim(pomdp::POMDP,
             policy::Policy,
             up=updater(policy),
             initial_belief=initial_state_distribution(pomdp),
             initial_state=nothing;
             rng::AbstractRNG=Base.GLOBAL_RNG,
             max_steps::Int=typemax(Int),
             simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps),
             metadata::Dict{Symbol}=Dict{Symbol, Any}()
            )

    if initial_state == nothing && state_type(pomdp) != Void
        is = rand(rng, initial_belief)
    else
        is = initial_state
    end
    return POMDPSim(simulator, pomdp, policy, up, initial_belief, is, metadata)
end

"""
    Sim(p::MDP, policy::Policy, metadata=Dict(:note=>"a note"))
    Sim(p::MDP, policy::Policy[, initial_state]; kwargs...)

Create a `Sim` object that represents a MDP simulation.

A vector of `Sim` objects can be executed with `run` or `run_parallel`.

## Keyword Arguments
- `rng::AbstractRNG=Base.GLOBAL_RNG`
- `max_steps::Int=typemax(Int)`
- `simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps)`
- `metadata::Dict{Symbol}=Dict{Symbol, Any}()` a dictionary of metadata for the sim that will be recorded, e.g. `Dict(:solver_iterations=>500)`.
"""
function Sim(mdp::MDP,
             policy::Policy,
             initial_state=nothing;
             rng::AbstractRNG=Base.GLOBAL_RNG,
             max_steps::Int=typemax(Int),
             simulator::Simulator=HistoryRecorder(rng=rng, max_steps=max_steps),
             metadata::Dict{Symbol}=Dict{Symbol, Any}()
            )

    if initial_state == nothing && state_type(mdp) != Void
        is = POMDPs.initial_state(mdp, rng)
    else
        is = initial_state
    end
    return MDPSim(simulator, mdp, policy, is, metadata)
end

POMDPs.simulate(s::POMDPSim) = simulate(s.simulator, s.pomdp, s.policy, s.updater, s.initial_belief, s.initial_state)
POMDPs.simulate(s::MDPSim) = simulate(s.simulator, s.mdp, s.policy, s.initial_state)

default_process(s::Sim, r::Float64) = :reward=>r
default_process(s::Sim, hist::SimHistory) = default_process(s, discounted_reward(hist))

run_parallel(queue::AbstractVector; kwargs...) = run_parallel(default_process, queue; kwargs...)

"""
    run_parallel(queue::Vector{Sim})
    run_parallel(f::Function, queue::Vector{Sim})

Run `Sim` objects in `queue` in parallel and return results as a `DataFrame`.

By default, the `DataFrame` will contain the reward for each simulation and the metadata provided to the sim.

# Arguments
- `queue`: List of `Sim` objects to be executed
- `f`: Function to process the results of each simulation
This function should take two arguments, (1) the `Sim` that was executed and (2) the result of the simulation, by default a `SimHistory`. It should return a dictionary or vector of pairs of `Symbol`s and values that will appear in the dataframe. See Examples below.

## Keyword Arguments
- `progress`: a `ProgressMeter.Progress` for showing progress through the simulations; `progress=false` will suppress the progress meter

# Examples

```julia
run_parallel(queue) do sim, hist
    return [:n_steps=>n_steps(hist), :reward=>discounted_reward(hist)]
end
```
will return a dataframe with with the number of steps and the reward in it.
"""
function run_parallel(process::Function, queue::AbstractVector;
                      progress=Progress(length(queue), desc="Simulating..."),
                      proc_warn=true)

    #=
    frame_lines = pmap(progress, queue) do sim
        result = simulate(sim)
        return process(sim, result)
    end
    =#

    np = nprocs()
    if np == 1 && proc_warn
        warn("""
             run_parallel(...) was started with only 1 process, so simulations will be run in serial.

             To suppress this warning, use run_parallel(..., proc_warn=false).

             To use multiple processes, use addprocs() or the -p option (e.g. julia -p 4).
             """)
    end
    n = length(queue)
    i = 1
    prog = 0
    # based on the simple implementation of pmap here: https://docs.julialang.org/en/latest/manual/parallel-computing
    frame_lines = Vector{Any}(n)
    nextidx() = (idx=i; i+=1; idx)
    prog_lock = ReentrantLock()
    @sync begin
        for p in 1:np
            if np == 1 || p != myid()
                @async begin
                    while true
                        idx = nextidx()
                        if idx > n
                            break
                        end
                        frame_lines[idx] = remotecall_fetch(p, queue[idx]) do sim
                                result = simulate(sim)
                                output = process(sim, result)
                                append_metadata(output, sim.metadata)
                        end
                        if progress isa Progress
                            lock(prog_lock)
                            update!(progress, prog+=1)
                            unlock(prog_lock)
                        end
                    end
                end
            end
        end
    end
    if progress isa Progress
        lock(prog_lock)
        finish!(progress)
        unlock(prog_lock)
    end

    return create_dataframe(frame_lines)
end

Base.run(queue::AbstractVector) = run(default_process, queue)

"""
    run(queue::Vector{Sim})
    run(f::Function, queue::Vector{Sim})

Run the `Sim` objects in `queue` on a single process and return the results as a dataframe.

See `run_parallel` for more information.
"""
function Base.run(process::Function, queue::AbstractVector; show_progress=true)
    lines = []
    if show_progress
        @showprogress for sim in queue
            result = simulate(sim)
            output = process(sim, result)
            line = append_metadata(output, sim.metadata)
            push!(lines, line)
        end
    else
        for sim in queue
            try
                result = simulate(sim)
                output = process(sim, result)
                line = append_metadata(output, sim.metadata)
                push!(lines, line)
            catch ex
                if ex isa InterruptException
                    rethrow(ex)
                else
                    println("Caught $ex in Base.run()")
                end
            end
        end
    end
    return create_dataframe(lines)
end

append_metadata(single::Pair, metadata::Dict) = append!(Any[single], collect(metadata))
append_metadata(pairvec::AbstractVector, metadata::Dict) = vcat(pairvec, collect(metadata))
append_metadata(d::Dict, metadata::Dict) = merge!(d, metadata)

metadata_as_pairs(s::Sim) = convert(Array{Any}, collect(s.metadata))

function create_dataframe(lines::Vector)
    master = Dict{Symbol, AbstractVector}()
    for line in lines
        push_line!(master, line)
    end
    return DataFrame(master)
end

function _push_line!(d::Dict{Symbol, AbstractVector}, line)
    if isempty(d)
        len = 0
    else
        len = length(first(values(d)))
    end
    for (key, val) in line
        T = typeof(val)
        if !haskey(d, key)
            d[key] = fill!(Array{Union{T,Missing}}(len), missing)
        end
        data = d[key]
        if !isa(val,eltype(data))
            d[key] = convert(Array{Any,1}, data)
        end
        push!(d[key], val)
    end
    for da in values(d)
        if length(da) < len + 1
            push!(da, missing)
        end
    end
    return d
end
push_line!(d::Dict{Symbol, AbstractVector}, line::Dict) = _push_line!(d, line)
push_line!(d::Dict{Symbol, AbstractVector}, line::DataFrame) = _push_line!(d, n=>first(line[n]) for n in names(line))
push_line!(d::Dict{Symbol, AbstractVector}, line::AbstractVector) = _push_line!(d, line)
push_line!(d::Dict{Symbol, AbstractVector}, line::Tuple) = _push_line!(d, line)
push_line!(d::Dict{Symbol, AbstractVector}, line::Pair) = _push_line!(d, (line,))


text:
using AdversarialDriving
using AutomotiveSimulator
using AutomotiveVisualization
using POMDPSimulators
using POMDPPolicies
using POMDPs
using Test
using Distributions

## Test construction of a BV
bv1 = BlinkerVehicleAgent(up_left(id = 10), TIDM(Tint_TIDM_template))
@test bv1.get_initial_entity().state isa BlinkerState
@test bv1.model.idm == Tint_TIDM_template.idm
@test bv1.entity_dim == BLINKERVEHICLE_ENTITY_DIM
@test bv1.disturbance_dim == BLINKERVEHICLE_DISTURBANCE_DIM
@test length(bv1.entity_to_vec(bv1.get_initial_entity())) == BLINKERVEHICLE_ENTITY_DIM
@test bv1.entity_to_vec(bv1.get_initial_entity())  == bv1.entity_to_vec(bv1.vec_to_entity(bv1.entity_to_vec(bv1.get_initial_entity()), id(bv1), Tint_roadway, bv1.model))
@test bv1.disturbance_model.actions == get_bv_actions().actions
@test bv1.disturbance_model.probs == get_bv_actions().probs
@test id(bv1) == 10

## Test construction of a Adversarial Pedestrian
ped1 = NoisyPedestrianAgent(get_pedestrian(id=2, s=0., v=4.), AdversarialPedestrian())
@test ped1.get_initial_entity().state isa NoisyPedState
@test ped1.entity_dim == PEDESTRIAN_ENTITY_DIM
@test ped1.disturbance_dim == PEDESTRIAN_DISTURBANCE_DIM
@test length(ped1.entity_to_vec(ped1.get_initial_entity())) == PEDESTRIAN_ENTITY_DIM
@test ped1.entity_to_vec(ped1.get_initial_entity())  == ped1.entity_to_vec(ped1.vec_to_entity(ped1.entity_to_vec(ped1.get_initial_entity()), id(ped1), ped_roadway, ped1.model))
@test ped1.disturbance_model.actions == get_ped_actions().actions
@test ped1.disturbance_model.probs == get_ped_actions().probs
@test id(ped1) == 2

# #Test construction of a full MDP with BV in the Tintersection
bv2 = BlinkerVehicleAgent(left_straight(id=2), TIDM(Tint_TIDM_template))
bv3 = BlinkerVehicleAgent(right_turnleft(id=3), TIDM(Tint_TIDM_template))
bv4 = BlinkerVehicleAgent(left_turnright(id=4, s=40.), TIDM(Tint_TIDM_template))
mdp = AdversarialDrivingMDP(bv1, [bv2, bv3, bv4], Tint_roadway, 0.1)

@test agents(mdp) == [adversaries(mdp)..., sut(mdp)]
@test sutid(mdp) == 10
@test length(agents(mdp)) == 4
@test model(mdp, 10) == agents(mdp)[4].model
@test model(mdp, 2) == agents(mdp)[1].model
@test model(mdp, 3) == agents(mdp)[2].model
@test model(mdp, 4) == agents(mdp)[3].model
@test mdp.num_adversaries == 3
@test mdp.roadway == Tint_roadway
scene = rand(initialstate(mdp))
@test scene[1] == bv2.get_initial_entity()
@test scene[2] == bv3.get_initial_entity()
@test scene[3] == bv4.get_initial_entity()
@test scene[4] == bv1.get_initial_entity()
@test mdp.dt == 0.1
@test mdp.last_observation == Float64[]
@test length(convert_s(AbstractArray, scene, mdp)) == 16
@test mdp.last_observation == convert_s(AbstractArray, scene, mdp)
d = combine_discrete(collect(adversaries(mdp)))
@test mdp.disturbance_model.actions == d.actions
@test mdp.disturbance_model.probs == d.probs
# @test actionindex(mdp, acts[4]) == 4
acts = d.actions
@test exp(logpdf(mdp, rand(initialstate(mdp)), acts[6])) ≈ d.probs[6]
@test mdp.γ == discount(mdp)
@test mdp.γ == 1.
@test mdp.ast_reward == false
@test mdp.no_collision_penalty == 1e3
@test mdp.scale_reward == true
@test mdp.end_of_road == Inf

mdp_temp = AdversarialDrivingMDP(bv1, [bv2, bv3, bv4], Tint_roadway, 0.1, γ=0.9, ast_reward = true, no_collision_penalty = 1e7, scale_reward = false, end_of_road = 70.)
@test mdp_temp.γ == .9
@test mdp_temp.ast_reward == true
@test mdp_temp.no_collision_penalty == 1e7
@test mdp_temp.scale_reward == false
@test mdp_temp.end_of_road == 70.

## Test the update_adversary! function
bv5 = BlinkerVehicleAgent(left_straight(id=5), TIDM(Tint_TIDM_template))
s = rand(initialstate(mdp))
sbefore = deepcopy(s)

update_adversary!(bv5, actions(mdp)[2][1], s)
@test all([noise(s[i]) == noise(sbefore[i]) for i=1:4])

noise_action = BlinkerVehicleControl(noise = Noise((0.2, 0.3), 1.))
update_adversary!(agents(mdp)[1], noise_action, s).state.noise.pos
@test noise(s[1]).pos == noise_action.noise.pos
@test noise(s[1]).vel == noise_action.noise.vel
@test model(mdp, id(agents(mdp)[1])).next_action == noise_action

## Test reward, isterminal
empty_scene = Scene(Entity{BlinkerState, VehicleDef, Int64})
coll_scene = Scene([left_straight(id=1, s=50.)(), right_turnleft(id=2, s=50.)(), up_left(id=10)()])
ego_coll_scene = Scene([left_straight(id=1, s=50.)(), right_turnleft(id=2)(), up_left(id=10, s=50.)()])

@test !isterminal(mdp, s)
@test isterminal(mdp, empty_scene)
@test isterminal(mdp, coll_scene)
@test isterminal(mdp, ego_coll_scene)
@test reward(mdp, s, actions(mdp)[1], coll_scene) == 0.
@test reward(mdp, s, actions(mdp)[1], ego_coll_scene) == 1.0

# Test gen using action comparisons
sp, r = gen(mdp, s, actions(mdp)[1])
sp2, r2 = gen(mdp, s, actions(mdp)[2]) # Slowdown of first adversary
@test vel(get_by_id(sp, 2)) > vel(get_by_id(sp2, 2))

# Run full simulation with random policy
mdp.dt = 0.1
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, FunctionPolicy((s) -> actions(mdp)[1]))
# @test length(hist) == 52
# nticks = length(hist)
# timestep = mdp.dt
# scenes = state_hist(hist)
# using Reel
# animation = roll(fps=1.0/timestep, duration=nticks*timestep) do t, dt
#     i = Int(floor(t/dt)) + 1
#     render([Tint_roadway, crosswalk, scenes[i]], canvas_width=1200, canvas_height=800)
# end
# write("ped_roadway_animated.gif", animation)

## Test the reward function from rollouts (both with and without ast)
sut_agent = BlinkerVehicleAgent(up_left(id=1, s=25., v=15.), TIDM(Tint_TIDM_template, noisy_observations = true))
adv1 = BlinkerVehicleAgent(left_straight(id=2, s=20., v=15.0), TIDM(Tint_TIDM_template))
mdp = AdversarialDrivingMDP(sut_agent, [adv1], Tint_roadway, 0.15)
blinker_action = mdp.disturbance_model.actions[7]
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, FunctionPolicy((s) -> posf(get_by_id(s,1)).s <=25 ? blinker_action : actions(mdp)[1]))
@test undiscounted_reward(hist) == 1

mdp = AdversarialDrivingMDP(sut_agent, [adv1], Tint_roadway, 0.15, ast_reward = true)
hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, FunctionPolicy((s) -> posf(get_by_id(s,1)).s <=25 ? blinker_action : actions(mdp)[1]))
@test undiscounted_reward(hist) < 1
@test undiscounted_reward(hist) > 0.8

hist = POMDPSimulators.simulate(HistoryRecorder(), mdp, FunctionPolicy((s) -> actions(mdp)[1]))
@test undiscounted_reward(hist) < 0

## Test random initial states
sut_agent = BlinkerVehicleAgent(rand_up_left(id=1, s_dist=Uniform(15,35.), v_dist=Uniform(7,25.)), TIDM(Tint_TIDM_template, noisy_observations = true))
adv1 = BlinkerVehicleAgent(rand_left(id=2, s_dist=Uniform(15,35.), v_dist=Uniform(7,25.)), TIDM(Tint_TIDM_template))
mdp = AdversarialDrivingMDP(sut_agent, [adv1], Tint_roadway, 0.15)

s1 = rand(initialstate(mdp))
s2 = rand(initialstate(mdp))

@test posf(get_by_id(s1, 1)) != posf(get_by_id(s2, 1))
for i=1:100
    s1 = rand(initialstate(mdp))
    @test posf(get_by_id(s1, 1)).s >= 15. && posf(get_by_id(s1, 1)).s <= 35
    @test vel(get_by_id(s1, 1)) >= 7. && vel(get_by_id(s1, 1)) <= 25.
end

## Test decompose indices
sut_agent = BlinkerVehicleAgent(up_left(id=1, s=25., v=15.), TIDM(Tint_TIDM_template, noisy_observations = true))
adv1 = BlinkerVehicleAgent(left_straight(id=2, s=0., v=15.0), TIDM(Tint_TIDM_template))
adv2 = BlinkerVehicleAgent(left_turnright(id=3, s=10., v=15.0), TIDM(Tint_TIDM_template))
adv3 = BlinkerVehicleAgent(right_straight(id=4, s=40., v=20.0), TIDM(Tint_TIDM_template))
adv4 = BlinkerVehicleAgent(right_turnleft(id=5, s=30., v=20.0), TIDM(Tint_TIDM_template))

mdp = AdversarialDrivingMDP(sut_agent, [adv1, adv2, adv3, adv4], Tint_roadway, 0.1)

@test decompose_indices(mdp, ids = [2,3]) == [1,2,3,4,5,6,7,8]
@test decompose_indices(mdp, ids = [5,1]) == [13,14,15,16,17,18,19,20]
@test decompose_indices(mdp, ids = [1,4]) == [9,10,11,12,17,18,19,20]

## Create actions from dictionary
d = Dict(:da => rand(10),
         :toggle_blinker => rand(Bernoulli(0.1), 10),
         :toggle_goal => rand(Bernoulli(0.1), 10),
         :noise_s => rand(10),
         :noise_v => rand(10))

as = create_actions_1BV(d)

for i in 1:length(as)
    a = as[i][1]
    @test a.da == d[:da][i]
    @test a.toggle_goal == d[:toggle_goal][i]
    @test a.toggle_blinker == d[:toggle_blinker][i]
    @test a.noise.pos[1] == d[:noise_s][i]
    @test a.noise.vel == d[:noise_v][i]
end

sut_agent = BlinkerVehicleAgent(up_left(id=1, s=25., v=15.), TIDM(Tint_TIDM_template, noisy_observations = true))
adv1 = BlinkerVehicleAgent(left_straight(id=2, s=20., v=15.0), TIDM(Tint_TIDM_template))
mdp = AdversarialDrivingMDP(sut_agent, [adv1], Tint_roadway, 0.15)

d = Dict(:da => Uniform(),
         :toggle_blinker => Bernoulli(0.1),
         :toggle_goal => Bernoulli(0.1),
         :noise_s => Uniform(),
         :noise_v => Uniform())

h = POMDPSimulators.simulate(HistoryRecorder(),mdp, continuous_policy(d))

# Sample with blindspot
blindspot = Blindspot(rand([rand(Uniform(-π/3,-π/9)), rand(Uniform(π/9,π/3))]), rand(Uniform(π/20, π/9)))
adv_ped = NoisyPedestrianAgent(get_rand_pedestrian(id=1, s_dist=Uniform(9, 10), v_dist=Uniform(0, 0.1)), AdversarialPedestrian(ignore_idm = true))
sut_agent = BlinkerVehicleAgent(get_rand_ped_vehicle(id=2, s_dist=Uniform(0, 10.), v_dist=Uniform(6, 15)), TIDM(ped_TIDM_template, idm = IntelligentDriverModel(v_des = 15), blindspot = blindspot))
mdp = AdversarialDrivingMDP(sut_agent, [adv_ped], ped_roadway, 0.1, ast_reward=true, γ=0.95)

r = mean([POMDPSimulators.simulate(RolloutSimulator(max_steps = 100), mdp, FunctionPolicy((s) -> rand(actions(mdp)))) for i=1:100])
@test r < 0

max_stars_repo_name: ancorso/MAST

max_issues_repo_name: ancorso/MAST


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl
=#

insts = OrderedDict()
insts["minimal"] = [
    ((1, 1, 1, false, false, true),),
    ((1, 1, 1, false, true, false),),
    ((1, 1, 1, true, false, false),),
    ]
insts["fast"] = [
    ((3, 1, 3, false, false, true),),
    ((3, 1, 3, false, true, false),),
    ((3, 1, 3, true, false, false),),
    ((1, 2, 3, false, false, true),),
    ((1, 2, 3, false, true, false),),
    ((1, 2, 3, true, false, false),),
    ]
insts["various"] = [
    ((4, 2, 5, false, false, true),),
    ((4, 2, 5, false, true, false),),
    ((4, 2, 5, true, false, false),),
    ((4, 2, 10, false, true, false),),
    ((4, 2, 10, true, false, false),),
    ((4, 4, 5, false, true, false),),
    ((4, 3, 6, true, false, false),),
    ((4, 3, 6, false, true, false),),
    ]
return (NearestPolyMatJuMP, insts)


text:
using BinDeps
using Compat

@BinDeps.setup

libopus = library_dependency("libopus", aliases = ["libopus"])

@osx_only begin
  using Homebrew
  provides( Homebrew.HB, "opus", libopus, os = :Darwin )
end

provides( AptGet, "libopus0", libopus )

# Source build
provides( Sources,
          URI("http://downloads.xiph.org/releases/opus/opus-1.1.2.tar.gz"),
          SHA="0e290078e31211baa7b5886bcc8ab6bc048b9fc83882532da4a1a45e58e907fd",
          libopus )
provides( BuildProcess,
          Autotools(libtarget=".libs/libopus."*BinDeps.shlib_ext, configure_options = ["--libdir=$(BinDeps.libdir(libopus))"]),
          libopus )

@compat @BinDeps.install Dict(:libopus => :libopus)


text:
module TestSetMacro

module Clone
    import Accessors: setmacro, opticmacro, modifymacro
    macro optic(ex)
        opticmacro(identity, ex)
    end
    macro set(ex)
        setmacro(identity, ex)
    end
    macro modify(f, ex)
        modifymacro(identity, f, ex)
    end
end#module Clone

using Accessors: Accessors
using Test
using .Clone: Clone

using StaticArrays: @SMatrix
using StaticNumbers

@testset "setmacro, opticmacro isolation" begin

    # test that no symbols like `IndexLens` are needed:
    Clone.@optic(_                                         )
    Clone.@optic(_.a                                       )
    Clone.@optic(_[1]                                      )
    Clone.@optic(first(_)                                  )
    Clone.@optic(_[end]                                    )
    Clone.@optic(_[static(1)]                              )
    Clone.@optic(_.a[1][end, end-2].b[static(1), static(1)])

    @test Accessors.@optic(_.a) === Clone.@optic(_.a)
    @test Accessors.@optic(_.a.b) === Clone.@optic(_.a.b)
    @test Accessors.@optic(_.a.b[1,2]) === Clone.@optic(_.a.b[1,2])

    o = (a=1, b=2)
    @test Clone.@set(o.a = 2) === Accessors.@set(o.a = 2)
    @test Clone.@set(o.a += 2) === Accessors.@set(o.a += 2)

    @test Clone.@modify(x -> x+1, o.a) === Accessors.@modify(x -> x+1, o.a)

    m = @SMatrix [0 0; 0 0]
    m2 = Clone.@set m[end-1, end] = 1
    @test m2 === @SMatrix [0 1; 0 0]
    m3 = Clone.@set(first(m) = 1)
    @test m3 === @SMatrix[1 0; 0 0]
end

# inference

macro setfield156(expr)
    # Example of macro that caused inference issues,
    # see https://github.com/jw3126/Setfield.jl/pull/156
    quote
        function f($(esc(:x)))
            $(Accessors.setmacro(identity, expr, overwrite=true))
            $(Accessors.setmacro(identity, expr, overwrite=true))
            $(Accessors.setmacro(identity, expr, overwrite=true))
            $(Accessors.setmacro(identity, expr, overwrite=true))
            $(Accessors.setmacro(identity, expr, overwrite=true))
            return $(esc(:x))
        end
    end
end

function test_all_inferrable(f, argtypes)
    # hacky, maybe JETTest can help?
    # * JETTest.@test_nodispatch does not detect the problem
    typed = first(code_typed(f, argtypes))
    code = typed.first
    @test all(T -> !(T isa UnionAll || T === Any), code.slottypes)
end

@testset "setmacro multiple usage" begin
    let f = @setfield156(x[end] = 1)
        test_all_inferrable(f, (Vector{Float64}, ))
    end
end

end#module


text:
# # Building and analysing a small community model

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/@__NAME__.ipynb)

# Here we will use `COBREXA` to build and analyze a small community model
# consisting of three *E. coli* mutants using the `CoreModel`. We will use an
# objective function that enforces equal growth rates.

# We will first construct a community of only two mutants to illustrate the
# effect of the community biomass objective function. Then we will add a third
# member that has a lethal knockout. Due to the bounds on the exchange reactions
# these three models are incapable of sharing resources - hence the maximum
# growth rate will be zero. By changing the bounds we can allow resource sharing,
# "saving" the community.

# ## Load the base model

!isfile("e_coli_core.json") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json");

using COBREXA
using Tulip

# ## Load the models and inspect fba solutions

base_model = load_model(CoreModel, "e_coli_core.json") # base from from which the knockouts will be constructed

cytbd_knockout_model = remove_reaction(base_model, "CYTBD") # knockout the CYTBD (cytochrome oxidase) reaction
sol = flux_balance_analysis_dict(cytbd_knockout_model, Tulip.Optimizer)
sol["BIOMASS_Ecoli_core_w_GAM"] # Cytochrome oxidase knockout μ (growth rate)
#
atps4r_knockout_model = remove_reaction(base_model, "ATPS4r") # knockout the ATP synthase reaction
sol = flux_balance_analysis_dict(atps4r_knockout_model, Tulip.Optimizer)
sol["BIOMASS_Ecoli_core_w_GAM"] # ATP synthase knockout μ
#
eno_knockout_model = remove_reaction(base_model, "ENO") # knockout the enolase reaction
sol = flux_balance_analysis_dict(eno_knockout_model, Tulip.Optimizer)
sol["BIOMASS_Ecoli_core_w_GAM"] # Enolase knockout μ, cannot grow by itself

# ## Build a community model of the cytochrome oxidase knockout and the ATP synthase knockout models

ex_rxn_mets = Dict(
    ex_rxn => first(keys(reaction_stoichiometry(base_model, ex_rxn))) for
    ex_rxn in filter(looks_like_exchange_reaction, reactions(base_model))
) # identify exchange reactions heuristically
#
model_names = ["cytbd_ko", "atps4r_ko"]
community_model = join_with_exchanges(
    CoreModel,
    [cytbd_knockout_model, atps4r_knockout_model],
    ex_rxn_mets;
    biomass_ids = ["BIOMASS_Ecoli_core_w_GAM", "BIOMASS_Ecoli_core_w_GAM"],
    model_names = model_names,
)

# ## Set exchange reaction bounds of community model based on the bounds of the individual models

env_ex_rxn_idxs = indexin(keys(ex_rxn_mets), reactions(community_model)) # identify the global (environmental exchange reactions)
cytbd_ex_rxn_idxs = indexin(keys(ex_rxn_mets), reactions(cytbd_knockout_model)) # identify the indices of the corresponding exchange reactions in the original models
atps4r_ex_rxn_idxs = indexin(keys(ex_rxn_mets), reactions(atps4r_knockout_model))

# In case some exchange reactions are not present in both models, set
# environmental exchange bound to the sum of the individual exchange bounds
for (env_ex, m2_ex, m1_ex) in zip(env_ex_rxn_idxs, cytbd_ex_rxn_idxs, atps4r_ex_rxn_idxs)
    m2lb = isnothing(m2_ex) ? 0.0 : atps4r_knockout_model.xl[m2_ex]
    m2ub = isnothing(m2_ex) ? 0.0 : atps4r_knockout_model.xu[m2_ex]
    m1lb = isnothing(m1_ex) ? 0.0 : cytbd_knockout_model.xl[m1_ex]
    m1ub = isnothing(m1_ex) ? 0.0 : cytbd_knockout_model.xu[m1_ex]
    change_bounds!(community_model, [env_ex]; lower = [m1lb + m2lb], upper = [m1ub + m2ub])
end

# ## Add objective function to community model`

biomass_ids_weights = Dict(model_names .* "_BIOMASS_Ecoli_core_w_GAM" .=> 1.0)
update_community_objective!(community_model, "community_biomass", biomass_ids_weights)

# ## Perform community FBA

d = flux_balance_analysis_dict(
    community_model,
    Tulip.Optimizer;
    modifications = [change_optimizer_attribute("IPM_IterationsLimit", 1000)],
)
d["community_biomass"] # community μ
# Notice, the growth rate is limited to the slowest organism as per the objective function

# ## Add the enolase knockout to the community model

community_model = add_model_with_exchanges(
    community_model,
    eno_knockout_model,
    ex_rxn_mets;
    model_name = "eno_ko",
    biomass_id = "BIOMASS_Ecoli_core_w_GAM",
)

push!(model_names, "eno_ko")
biomass_ids_weights = Dict(model_names .* "_BIOMASS_Ecoli_core_w_GAM" .=> 1.0)
update_community_objective!(community_model, "community_biomass", biomass_ids_weights)

d = flux_balance_analysis_dict(
    community_model,
    Tulip.Optimizer;
    modifications = [change_optimizer_attribute("IPM_IterationsLimit", 1000)],
)
println("Community μ = ", d["community_biomass"])
# Notice that the high communal growth rate is 0, due to the enolase knockout.
# The reason for this behaviour: enolase is a central reaction in glycolysis - without
# it the organism cannot access the lower glycolysis pathways or the TCA cycle, hence
# the model predicts no growth for the knockout, and hence no growth for the system since
# they all have to have the same growth rate.

# ## Allow the mutants to rescue each other by sharing pyruvate

pyr_exs = model_names .* "_EX_pyr_e"
change_bounds!(community_model, pyr_exs; lower = fill(-1000.0, 3), upper = fill(1000.0, 3))

d = flux_balance_analysis_dict(
    community_model,
    Tulip.Optimizer;
    modifications = [change_optimizer_attribute("IPM_IterationsLimit", 1000)],
)
d["community_biomass"] # community μ
# Notice that the growth rate is now above 0! Nutrient sharing saved the day!


text:
using FractionalDiffEq
using Plots, LaTeXStrings

s="\$Bagley\\ Torvik\\ Equation\$"

T = 30;h = 0.05
tspan = collect(0:h:T)
result = bagleytorvik(1, 1, 1, T, h)

plot(tspan, result, title=s, legend=:bottomright)


text:
"""
    kfilter!(settings::KalmanSettings, status::KalmanStatus)

Kalman filter: a-priori prediction and a-posteriori update.

# Model
The state space model used below is,

``Y_{t} = B*X_{t} + e_{t}``

``X_{t} = C*X_{t-1} + v_{t}``

Where ``e_{t} ~ N(0, R)`` and ``v_{t} ~ N(0, V)``.

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
"""
function kfilter!(settings::KalmanSettings, status::KalmanStatus)

    # Update status.t
    status.t += 1;

    # A-priori prediction
    apriori!(typeof(status.X_prior), settings, status);

    # Handle missing observations
    ind_not_missings = find_observed_data(settings, status);

    # Ex-post update
    aposteriori!(settings, status, ind_not_missings);

    # Update history of *_prior and *_post
    if settings.store_history == true
        push!(status.history_X_prior, status.X_prior);
        push!(status.history_X_post, status.X_post);
        push!(status.history_P_prior, status.P_prior);
        push!(status.history_P_post, status.P_post);
    end
end

"""
    apriori!(::Type{Nothing}, settings::KalmanSettings, status::KalmanStatus)

Kalman filter a-priori prediction for t==1.

# Arguments
- `::Type{Nothing}`: first prediction
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct

    apriori!(::Type{FloatVector}, settings::KalmanSettings, status::KalmanStatus)

Kalman filter a-priori prediction.

# Arguments
- `::Type{FloatVector}`: standard prediction
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
"""
function apriori!(::Type{Nothing}, settings::KalmanSettings, status::KalmanStatus)

    status.X_prior = apriori(settings.X0, settings);
    status.P_prior = apriori(settings.P0, settings);

    if settings.compute_loglik == true
        status.loglik = 0.0;
    end

    if settings.store_history == true
        status.history_X_prior = Array{FloatVector,1}();
        status.history_X_post = Array{FloatVector,1}();
        status.history_P_prior = Array{SymMatrix,1}();
        status.history_P_post = Array{SymMatrix,1}();
    end
end

function apriori!(::Type{FloatVector}, settings::KalmanSettings, status::KalmanStatus)
    status.X_prior = apriori(status.X_post, settings);
    status.P_prior = apriori(status.P_post, settings);
end

"""
    apriori(X::FloatVector, settings::KalmanSettings)

Kalman filter a-priori prediction for X.

# Arguments
- `X`: Last expected value of the states
- `settings`: KalmanSettings struct

    apriori(P::SymMatrix, settings::KalmanSettings)

Kalman filter a-priori prediction for P.

# Arguments
- `P`: Last conditional covariance the states
- `settings`: KalmanSettings struct
"""
apriori(X::FloatVector, settings::KalmanSettings) = settings.C * X;
apriori(P::SymMatrix, settings::KalmanSettings) = Symmetric(settings.C * P * settings.C' + settings.V)::SymMatrix;

"""
    find_observed_data(settings::KalmanSettings, status::KalmanStatus)

Return position of the observed measurements at time t.

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
"""
function find_observed_data(settings::KalmanSettings, status::KalmanStatus)
    if status.t <= settings.T
        Y_t_all = @view settings.Y[:, status.t];
        ind_not_missings = findall(ismissing.(Y_t_all) .== false);
        if length(ind_not_missings) > 0
            return ind_not_missings;
        end
    end
end

"""
    aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::Array{Int64,1})

Kalman filter a-posteriori update. Measurements are observed (or partially observed) at time t.

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
- `ind_not_missings`: Position of the observed measurements

    aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::Nothing)

Kalman filter a-posteriori update. All measurements are not observed at time t.

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
- `ind_not_missings`: Empty array
"""
function aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::Array{Int64,1})

    Y_t = @view settings.Y[ind_not_missings, status.t];
    B_t = @view settings.B[ind_not_missings, :];
    R_t = @view settings.R[ind_not_missings, ind_not_missings];

    # Forecast error
    ε_t = Y_t - B_t*status.X_prior;
    Σ_t = Symmetric(B_t*status.P_prior*B_t' + R_t)::SymMatrix;

    # Kalman gain
    K_t = status.P_prior*B_t'*inv(Σ_t);

    # A posteriori estimates
    status.X_post = status.X_prior + K_t*ε_t;
    status.P_post = Symmetric(status.P_prior - K_t*B_t*status.P_prior)::SymMatrix;

    # Update log likelihood
    if settings.compute_loglik == true
        update_loglik!(status, ε_t, Σ_t);
    end
end

function aposteriori!(settings::KalmanSettings, status::KalmanStatus, ind_not_missings::Nothing)
    status.X_post = copy(status.X_prior);
    status.P_post = copy(status.P_prior);
end

"""
    update_loglik!(status::KalmanStatus, ε_t::FloatVector, Σ_t::SymMatrix)

Update status.loglik.

# Arguments
- `status`: KalmanStatus struct
- `ε_t`: Forecast error
- `Σ_t`: Forecast error covariance
"""
function update_loglik!(status::KalmanStatus, ε_t::FloatVector, Σ_t::SymMatrix)
    status.loglik -= 0.5*(logdet(Σ_t) + ε_t'*inv(Σ_t)*ε_t);
end

"""
    kforecast(settings::KalmanSettings, X::Union{FloatVector, Nothing}, h::Int64)

Forecast X up to h-steps ahead.

    kforecast(settings::KalmanSettings, X::Union{FloatVector, Nothing}, P::Union{SymMatrix, Nothing}, h::Int64)

Forecast X and P up to h-steps ahead.

# Model
The state space model used below is,

``Y_{t} = B*X_{t} + e_{t}``

``X_{t} = C*X_{t-1} + v_{t}``

Where ``e_{t} ~ N(0, R)`` and ``v_{t} ~ N(0, V)``.
"""
function kforecast(settings::KalmanSettings, Xt::Union{FloatVector, Nothing}, h::Int64)

    # Initialise forecast history
    history_X = Array{FloatVector,1}();

    X = copy(Xt);

    # Loop over forecast horizons
    for horizon=1:h
        X = apriori(X, settings);
        push!(history_X, X);
    end

    # Return output
    return history_X;
end

function kforecast(settings::KalmanSettings, Xt::Union{FloatVector, Nothing}, Pt::Union{SymMatrix, Nothing}, h::Int64)

    # Initialise forecast history
    history_X = Array{FloatVector,1}();
    history_P = Array{SymMatrix,1}();

    X = copy(Xt);
    P = copy(Pt);

    # Loop over forecast horizons
    for horizon=1:h
        X = apriori(X, settings);
        P = apriori(P, settings);
        push!(history_X, X);
        push!(history_P, P);
    end

    # Return output
    return history_X, history_P;
end

"""
    ksmoother(settings::KalmanSettings, status::KalmanStatus)

Kalman smoother: RTS smoother from the last evaluated time period in status to t==0.

# Model
The state space model used below is,

``Y_{t} = B*X_{t} + e_{t}``

``X_{t} = C*X_{t-1} + v_{t}``

Where ``e_{t} ~ N(0, R)`` and ``v_{t} ~ N(0, V)``.

# Arguments
- `settings`: KalmanSettings struct
- `status`: KalmanStatus struct
"""
function ksmoother(settings::KalmanSettings, status::KalmanStatus)

    # Initialise smoother history
    history_X = Array{FloatVector,1}();
    history_P = Array{SymMatrix,1}();
    push!(history_X, copy(status.X_post));
    push!(history_P, copy(status.P_post));

    # Loop over t
    for t=status.t:-1:2

        # Pointers
        Xp = status.history_X_prior[t];
        Pp = status.history_P_prior[t];
        Xf_lagged = status.history_X_post[t-1];
        Pf_lagged = status.history_P_post[t-1];
        Xs = history_X[1];
        Ps = history_P[1];

        # Smoothed estimates for t-1
        J1 = compute_J1(Pf_lagged, Pp, settings);
        pushfirst!(history_X, backwards_pass(Xf_lagged, J1, Xs, Xp));
        pushfirst!(history_P, backwards_pass(Pf_lagged, J1, Ps, Pp));
    end

    # Pointers
    Xp = status.history_X_prior[1];
    Pp = status.history_P_prior[1];
    Xs = history_X[1];
    Ps = history_P[1];

    # Compute smoothed estimates for t==0
    J1 = compute_J1(settings.P0, Pp, settings);
    X0 = backwards_pass(settings.X0, J1, Xs, Xp);
    P0 = backwards_pass(settings.P0, J1, Ps, Pp);

    # Return output
    return history_X, history_P, X0, P0;
end

"""
    compute_J1(Pf_lagged::SymMatrix, Pp::SymMatrix, settings::KalmanSettings)

Compute J_{t-1} as in Shumway and Stoffer (2011, chapter 6) to operate the RTS smoother.
"""
compute_J1(Pf_lagged::SymMatrix, Pp::SymMatrix, settings::KalmanSettings) = Pf_lagged*settings.C'*inv(Pp);

"""
    backwards_pass(Xf_lagged::FloatVector, J1::FloatArray, Xs::FloatVector, Xp::FloatVector)

Backwards pass for X to get the smoothed states at time t-1.

# Arguments
- `Xf_lagged`: Filtered states for time t-1
- `J1`: J_{t-1} as in Shumway and Stoffer (2011, chapter 6)
- `Xs`: Smoothed states for time t
- `Xp`: A-priori states for time t

    backwards_pass(Pf_lagged::SymMatrix, J1::FloatArray, Ps::SymMatrix, Pp::SymMatrix)

Backwards pass for P to get the smoothed conditional covariance at time t-1.

# Arguments
- `Pf_lagged`: Filtered conditional covariance for time t-1
- `J1`: J_{t-1} as in Shumway and Stoffer (2011, chapter 6)
- `Ps`: Smoothed conditional covariance for time t
- `Pp`: A-priori conditional covariance for time t
"""
backwards_pass(Xf_lagged::FloatVector, J1::FloatArray, Xs::FloatVector, Xp::FloatVector) = Xf_lagged + J1*(Xs-Xp);
backwards_pass(Pf_lagged::SymMatrix, J1::FloatArray, Ps::SymMatrix, Pp::SymMatrix) = Symmetric(Pf_lagged + J1*(Ps-Pp)*J1')::SymMatrix;


text:
using Gadfly, DataArrays, RDatasets

custom_theme = Theme(panel_fill=colorant"black",
                     default_color=colorant"pink",
                     background_color=colorant"red")

plot(dataset("datasets", "iris"), x=:SepalLength, y=:SepalWidth, Geom.point,
     custom_theme)


text:
import LightXML: free, parse_string
import SeisIO: KW, TimeSpec, check_for_gap!, dtconst, endtime, mk_xml!,
  parsetimewin, read_station_xml!, split_id, sxml_mergehdr!, t_win, trunc_x!
import SeisIO.Quake:event_xml!, new_qml!, write_qml!


text:
using YaoBlocks.AD
using SparseArrays, LuxurySparse
using YaoArrayRegister
using LinearAlgebra
using Test

@testset "outer prod" begin
    D = 8
    T = ComplexF64
    for (l, r) in [(randn(T, D), randn(T, D)), (randn(T, D, 3), randn(T, D, 3))]
        a = outerprod(l, r)
        A = Matrix(a)
        B = randn(T, D, D)
        V = randn(T, D)
        @test A * B ≈ a * B
        @test A * A ≈ a * a
        @test B * A ≈ B * a
        @test V' * A ≈ V' * a
        @test A * V ≈ a * V
        @test size(a) == size(A)
        @test size(a, 2) == size(A, 2)
        @test A[4] == a[4]
        @test A[4, 4] == a[4, 4]
        @test A' == a'
        @test transpose(A) == transpose(a)
        @test conj(A) == conj(a)
        @test conj!(A) == conj!(a)
        if l isa Vector
            outerprod(ArrayReg(l), ArrayReg(r)) == a
        else
            outerprod(BatchedArrayReg(l), BatchedArrayReg(r)) == a
        end

        a = outerprod(l, r)
        @test rmul!(copy(a), 0.3) ≈ a .* 0.3

        if ndims(a) == 2
            @test AD._sum_A_Bconj(a, B) ≈ AD._sum_A_Bconj(Matrix(a), B)
        end

        @test 3 * a ≈ Matrix(a) * 3
        @test a * 3 ≈ Matrix(a) * 3
    end
end

@testset "projection" begin
    T = ComplexF64
    D = 8
    for y in [pmrand(T, D), sprand(T, D, D, 0.5), Diagonal(randn(T, D))]
        @test projection(y, y) == y
        l, r = randn(T, D), randn(T, D)
        op = outerprod(l, r)
        @test projection(y, op) == projection(y, Matrix(op))
    end
end


text:
using StochasticAiry
using Test
using Random: seed!

@testset "StochasticAiry.jl" begin
    seed!(1)

    n_re_lambda = 1000
    n_im_lambda = 200

    re_grid = collect(range(-8, 15, length=n_re_lambda))
    im_grid = collect(range(-8, 8, length=n_im_lambda))

    lambda_grid = re_grid' .+ im*im_grid;

    ##

    t_grid_len = 3*10^3

    t = -20
    s = 10
    t_grid = range(s, t, length=t_grid_len);

    dt = (t - s)/t_grid_len;

    for beta in [1, Inf]
        phi, phi_dash= rand(SAi(beta), :space, 0, lambda_grid)

        @test size(phi) == size(lambda_grid)
        @test size(phi_dash) == size(lambda_grid)

        phi_re, phi_dash_re = rand(SAi(beta), :space, 0, -8:8, brownian_gridlen = 100)

        @test eltype(phi_re) <: Real
        @test eltype(phi_dash_re) <: Real

        sai, sai_dash = rand(SAi(beta), :time, 1 + im, t_grid)

        @test ndims(sai) == 1
        @test ndims(sai_dash) == 1

        sai_re, sai_dash_re = rand(SAi(beta), :time, 0, t_grid)

        @test eltype(sai_re) <: Real
        @test eltype(sai_dash_re) <: Real
    end

end


text:
# include("../SimpleTasks.jl")

module AWSScheduler

using SimpleTasks.Types
using SimpleTasks.Services.AWSQueue
using SimpleTasks.Services.CLIBucket
# using SimpleTasks.Services.AWSCLIProvider
using ContactAnalysis.CountEdgesTask

import AWS
import JSON
import SimpleTasks.Services.Bucket
import SimpleTasks.Services.Queue
import SimpleTasks.Tasks.BasicTask

function schedule_count_edges(args...; queue_name="")
    env = AWS.AWSEnv()
    queue = AWSQueueService(env, queue_name)
    task = CountEdgesTask.CountEdgesTaskDetails(args...)
    # {"segmentation_storage":"in","slices":[[1,2],[3,4],[5,6]],"scale_idx":1,"output_storage":"out"}
    Queue.push_message(queue; message_body = JSON.json(task));
end

end # module AWSScheduler


text:
# macro.jl
# Macro-exercising speed tests
using JuMP

function test_linear(N)
    m = Model()
    @variable(m, x[1:10N,1:5N])
    @variable(m, y[1:N,1:N,1:N])

    for z in 1:10
        @constraint(m,
            9*y[1,1,1] - 5*y[N,N,N] -
            2*sum{ z*x[j,i*N],                j=((z-1)*N+1):z*N, i=3:4} +
              sum{ i*(9*x[i,j] + 3*x[j,i]),   i=N:2N,            j=N:2N} +
            x[1,1] + x[10N,5N] + x[2N,1] +
            1*y[1,1,N] + 2*y[1,N,1] + 3*y[N,1,1] +
            y[N,N,N] - 2*y[N,N,N] + 3*y[N,N,N]
             <=
            sum{sum{sum{N*i*j*k*y[i,j,k] + x[i,j],k=1:N; i!=j && j!=k},j=1:N},i=1:N} +
            sum{sum{x[i,j], j=1:5N; j % i == 3}, i=1:10N; i <= N*z}
            )
    end
end

function test_quad(N)
    m = Model()
    @variable(m, x[1:10N,1:5N])
    @variable(m, y[1:N,1:N,1:N])

    for z in 1:10
        @constraint(m,
            9*y[1,1,1] - 5*y[N,N,N] -
            2*sum{ z*x[j,i*N],                j=((z-1)*N+1):z*N, i=3:4} +
              sum{ i*(9*x[i,j] + 3*x[j,i]),   i=N:2N,            j=N:2N} +
            x[1,1] + x[10N,5N] * x[2N,1] +
            1*y[1,1,N] * 2*y[1,N,1] + 3*y[N,1,1] +
            y[N,N,N] - 2*y[N,N,N] * 3*y[N,N,N]
             <=
            sum{sum{sum{N*i*j*k*y[i,j,k] * x[i,j],k=1:N; i!=j && j!=k},j=1:N},i=1:N} +
            sum{sum{x[i,j], j=1:5N; j % i == 3}, i=1:10N; i <= N*z}
            )
    end
end

# Warmup
println("Test 1")
test_linear(1)
test_quad(1)
for N in [20,50,100]
    println("  Running N=$(N)...")
    N1_times = Any[]
    N2_times = Any[]
    for iter in 1:10
        tic()
        test_linear(N)
        push!(N1_times, toq())
        tic()
        test_quad(N)
        push!(N2_times, toq())
    end
    println("    N=$(N) min $(minimum(N1_times))")
    println("    N=$(N) min $(minimum(N2_times))")
end


text:
"""
    The id of the master worker
"""
MASTERW = 1
_set_MASTERW(w) = (global MASTERW = w)
set_MASTERW(MW::Int, ws::Vector{Int} = workers()) =
    (remotecall_wait.([_set_MASTERW], ws, [MW]); MW)

print_inmw(io::IO, ss...) = (remotecall_wait(Base.print, MASTERW, io, ss...); nothing)
print_inmw(ss...) = (remotecall_wait(Base.print, MASTERW, ss...); nothing)
println_inmw(io::IO, ss...) = (remotecall_wait(Base.println, MASTERW, io, ss...); nothing)
println_inmw(ss...) = (remotecall_wait(Base.println, MASTERW, ss...); nothing)
print_ifmw(io::IO, ss...) = (myid() == MASTERW) ? Base.print(io, ss...) : nothing
print_ifmw(ss...) = (myid() == MASTERW) ? Base.print(ss...) : nothing
println_ifmw(io::IO, ss...) = (myid() == MASTERW) ? Base.println(io, ss...) : nothing
println_ifmw(ss...) = (myid() == MASTERW) ? Base.println(ss...) : nothing

function wtag(io::IO)
    ws = min(displaysize(io) |> last, 80)
    return lpad(
        string( " Worker ", myid(), " (", getpid(), ") [", Time(now()), "] "),
        ws, "-"
    )
end

_tagprint(printfun::Function, io::IO, ss...; tag::String = wtag(io)) =
    !isempty(tag) ? printfun(io, tag, "\n", ss...) : printfun(io, ss...)
_tagprint(printfun::Function, ss...; kwargs...) = _tagprint(printfun, stdout, ss...; kwargs...)

tagprint(ss...; kwargs...) = _tagprint(print, ss...; kwargs...)
tagprintln(ss...; kwargs...) = _tagprint(println, ss...; kwargs...)

tagprint_inmw(ss...; tag::String = wtag(stdout)) = _tagprint(print_inmw, ss...; tag)
tagprintln_inmw(ss...; tag::String = wtag(stdout)) = _tagprint(println_inmw, ss...; tag)

tagprint_ifmw(ss...; tag::String = wtag(stdout)) = _tagprint(print_ifmw, ss...; tag)
tagprintln_ifmw(ss...; tag::String = wtag(stdout)) = _tagprint(println_ifmw, ss...; tag)


text:
module Issue980

import MathOptInterface
const MOI = MathOptInterface

abstract type StaticArray end
(::Type{SA})(x...) where {SA <: StaticArray} = SA(x)

function crash()
    model = MOI.Utilities.Model{Float64}()
    bridged = MOI.Bridges.LazyBridgeOptimizer(MOI.Utilities.Model{Float64}())
    MOI.copy_to(bridged, model)
end

crash()

end


text:
using Documenter
using EnergyExpressions
using AtomicLevels
using Atoms

isdefined(Main, :NOPLOTS) && NOPLOTS || include("plots.jl")

DocMeta.setdocmeta!(Atoms, :DocTestSetup, :(using Atoms, AtomicLevels); recursive=true)
makedocs(
    modules = [Atoms],
    sitename = "Atoms",
    pages = [
        "Home" => "index.md",
        "Examples" => "examples.md",
        "Theory" => [
            "Density matrices" => "density_matrices.md"
        ],
        "Implementation" => [
            "Radial orbitals" => "radial_orbitals.md",
            "Atom types" => "atom_types.md",
            "One-body Hamiltonians" => "one_body.md",
            "Hydrogenic initialization" => "hydrogenic.md",
            "Orbital equations" => "orbital_equations.md",
            "Equation systems" => "equation_systems.md",
            "Observables" => "observables.md"
        ]
    ],
    format = Documenter.HTML(assets = ["assets/latex.js"])
)

deploydocs(repo = "github.com/JuliaAtoms/Atoms.jl.git")


text:
@dict_readable struct HoverClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
    contentFormat::Union{Vector{String},Missing}
end

@dict_readable struct HoverParams <: Outbound
    textDocument::TextDocumentIdentifier
    position::Position
    workDoneToken::Union{Int,String,Missing} # ProgressToken
end

struct HoverOptions <: Outbound
    workDoneProgress::Union{Bool,Missing}
end

struct HoverRegistrationOptions <: Outbound
    workDoneProgress::Union{Bool,Missing}
    documentSelector::Union{DocumentSelector,Nothing}
end

struct Hover <: Outbound
    contents::Union{MarkedString,Vector{MarkedString},MarkupContent}
    range::Union{Range,Missing}
end