text:
# This file defines some struct types used in the SMLMData package.
"""
    SMLD

 SMLD is the highest level abstract type in SMLMData.
"""
abstract type SMLD end

"""
    SMLD2D

Structure containing 2D localization data

    # Fields
    - connectID: integer ID associating localizations
    - x: x position
    - y: y position
    - σ_x: standard error uncertainty in x
    - σ_y: standard error uncertainty in y
    - photons: total photons
    - σ_photons: standard error uncertainty in σ_photons
    - bg: fit fluorescence background (photons/pixel)
    - σ_bg: standard error uncertainty in bg
    - framenum: integer framenumber within a dataset
    - datasetnum: integer dataset number
    - datasize::Vector{Int}: size of image area
    - nframes: integer frames per dataset
    - ndatasets: number of dataasets
    - datafields: bookeeping - do not modify
"""
mutable struct SMLD2D <: SMLD
    connectID::Vector{Int}
    x::Vector{Float64}
    y::Vector{Float64}
    σ_x::Vector{Float64}
    σ_y::Vector{Float64}
    photons::Vector{Float64}
    σ_photons::Vector{Float64}
    bg::Vector{Float64}
    σ_bg::Vector{Float64}
    framenum::Vector{Int}
    datasetnum::Vector{Int}
    datasize::Vector{Int}
    nframes::Int
    ndatasets::Int
    datafields::NTuple{11, Symbol}
    SMLD2D() = new()
end

"""
    SMLD2D(nlocalizations::Int)

Constructor to generate an empty `smld` with a specific size.

# Description
This is a constructor for the SMLD2D struct which allows you to populate the
structure with undefined values for a predefined number of localizations.
"""
function SMLD2D(nlocalizations::Int)
    smld = SMLD2D()
    smld.x = zeros(Float64, nlocalizations)
    smld.y = zeros(Float64, nlocalizations)
    smld.σ_x = zeros(Float64, nlocalizations)
    smld.σ_y = zeros(Float64, nlocalizations)
    smld.bg = zeros(Float64, nlocalizations)
    smld.σ_bg = zeros(Float64, nlocalizations)
    smld.photons = zeros(Float64, nlocalizations)
    smld.σ_photons = zeros(Float64, nlocalizations)
    smld.connectID = collect(1:nlocalizations)
    smld.framenum = zeros(Int, nlocalizations)
    smld.datasetnum = zeros(Int, nlocalizations)
    smld.nframes = 0
    smld.ndatasets = 0
    smld.datasize = [0; 0]
    smld.datafields = (:connectID, :x, :y, :σ_x, :σ_y,
        :photons, :σ_photons, :bg, :σ_bg, :framenum, :datasetnum)

    return smld
end

"""
    SMLD2D(data::DataFrames.DataFrame)

Constructor to generate an `smld` from a data frame.

# Description
This is a constructor for the SMLD2D struct which allows you to populate the
structure with data defined in the dataframe `data`. The intention is that a
.csv table can be organized with localizations on each row as
`[datasetnum, framenum, x, y, σ_x, σ_y]`, loaded using the CSV package, and
placed into a dataframe with the DataFrames package.
"""
function SMLD2D(data::DataFrames.DataFrame)
    smld = SMLD2D()
    smld.connectID = collect(1:size(data)[1])
    smld.datasetnum = Int.(data[:, 1])
    smld.framenum = Int.(data[:, 2])
    smld.x = Float64.(data[:, 3])
    smld.y = Float64.(data[:, 4])
    smld.σ_x = Float64.(data[:, 5])
    smld.σ_y = Float64.(data[:, 6])
    smld.photons = Float64.(data[:, 7])
    smld.σ_photons = Float64.(data[:, 8])
    smld.bg = Float64.(data[:, 9])
    smld.σ_bg = Float64.(data[:, 10])
    smld.nframes = Int(maximum(smld.framenum))
    smld.ndatasets = Int(Base.length(unique(smld.datasetnum)))
    smld.datasize = [ceil(maximum(data[:, 3]) - 0.5);
                     ceil(maximum(data[:, 4]) - 0.5)]
    smld.datafields = (:connectID, :x, :y, :σ_x, :σ_y,
        :photons, :σ_photons, :bg, :σ_bg, :framenum, :datasetnum)

    return smld
end

"""
    Base.length(smld::SMLMData.SMLD2D)

Returns the number of localizations in `smld`.
"""
function Base.length(smld::SMLMData.SMLD2D)
    # Determine the length of `smld` by counting the number of frames.
    return Base.length(smld.framenum)
end

"""
    Base.getindex(smld::SMLMData.SMLD2D, ind)

Returns the requested localization from `smld`.
"""
function Base.getindex(smld::SMLMData.SMLD2D, ind)
    return SMLMData.isolatesmld(smld, ind)
end


text:
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using JuliaFormatter

main_paths = ["./docs/src"]
for main_path in main_paths
    for folder in readdir(main_path)
        @show folder_path = joinpath(main_path, folder)
        if isfile(folder_path)
            !occursin(".md", folder_path) && continue
        end
        format(
            folder_path;
            format_markdown=true,
            whitespace_ops_in_indices=true,
            remove_extra_newlines=true,
            verbose=true,
            always_for_in=true,
            whitespace_typedefs=true,
            whitespace_in_kwargs=false,
            # always_use_return = true # removed since it has false positives.
        )
    end
end

ext: jl


text:
function test_primes_iscoprime()
    @test iscoprime(1, 3) == true
    @test iscoprime(2, 3) == true
    @test iscoprime(11*13*17, 13*17*19) == false
    @test iscoprime(2*3, 5*7) == true
end

function test_primes_primes()
    @test primes(10) == [2, 3, 5, 7]
    @test primes(11) == [2, 3, 5, 7, 11]
    @test primes(1, 1) == primes(1:1) == Int[]
    @test primes(10, 10) == primes(10:10) == Int[]
    @test primes(10, 20) == primes(10:20) == [11, 13, 17, 19]
    @test primes(11, 19) == primes(11:19) == [11, 13, 17, 19]
    @test primes(1000, 1100) == primes(1000:1100) == [1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097]
    @test eltype(primes(10)) == typeof(2)
    @test eltype(primes(10, 20)) == eltype(primes(10:20)) == typeof(11)
end

function test_primes_primesmask()
    @test primesmask(100000) == primesmask(1:100000) == [isprime(i) for i in 1:100000]
    @test primesmask(2:100000) == [isprime(i) for i in 2:100000]
    @test primesmask(1000:100000) == [isprime(i) for i in 1000:100000]
    @test primesmask(1009:100000) == [isprime(i) for i in 1009:100000]
    @test typeof(primesmask(100000)) == typeof(primesmask(1:100000)) == BitVector
end

function test_primes_primepi()
    @test primepi(10) == primepi(1:10) == primepi(2:10) == 4
    @test primepi(11) == primepi(1:11) == primepi(2:11) == 5
    @test primepi(10000) == primepi(1:10000) == primepi(2:10000) == 1229
    @test primepi(10, 20) == primepi(10:20) == 4
    @test primepi(11, 19) == primepi(11:19) == 4
end

function test_primes_nthprime()
    p = primes(10000)
    @test [nthprime(i) for i in 1:1000] == [p[i] for i in 1:1000]
    @test [nthprime(a, a + 100) for a in 1:1000] ==  [nthprime(a:(a+100)) for a in 1:1000] == [p[a:(a+100)] for a in 1:1000]
    @test typeof(nthprime(100)) == eltype(nthprime(1, 100)) == eltype(nthprime(1:100)) == typeof(100)
end

function test_primes_nextprime()
    p = primes(2000)
    @test [nextprime(i) for i in -1:1000] == [first(dropwhile(p -> p <= i, p)) for i in -1:1000]
    @test typeof(nextprime(3)) == typeof(3)
    @test typeof(nextprime(big(3))) == typeof(big(3))
end

function test_primes_prevprime()
    @test prevprime(0) == 0
    @test prevprime(2) == 0
    p = reverse(primes(2000))
    @test [prevprime(i) for i in 3:1000] == [first(dropwhile(p -> p >= i, p)) for i in 3:1000]
    @test typeof(prevprime(3)) == typeof(3)
    @test typeof(prevprime(big(3))) == typeof(big(3))
end

function test_primes_twinprimes()
    @test twinprimes(1000, 1050) == twinprimes(1000:1050) == [1019 1021; 1031 1033]
    @test length(twinprimes(20, 30)) == 0
end

function test_primes_all()
    label = "Math.NumberTheory.Primes..."
    print(rpad(label, 50, ' '))

    test_primes_iscoprime()
    test_primes_primes()
    test_primes_primesmask()
    test_primes_primepi()
    test_primes_nthprime()
    test_primes_nextprime()
    test_primes_prevprime()
    test_primes_twinprimes()

    println("PASS")
end


text:
"Parameters for right-hand side function of variational partitioned Runge-Kutta methods."
mutable struct ParametersVPRKpMidpoint{DT, TT, ET <: IODE{DT,TT}, D, S} <: AbstractParametersVPRK{DT,TT,ET,D,S}
    equ::ET
    tab::TableauVPRK{TT}
    Δt::TT

    R::Vector{TT}

    t̅::TT
    q̅::Vector{DT}
    p̅::Vector{DT}
end

function ParametersVPRKpMidpoint(equ::ET, tab::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: IODE{DT,TT}}
    R = convert(Vector{TT}, [1, tab.R∞])
    q̅ = zeros(DT, ndims(equ))
    p̅ = zeros(DT, ndims(equ))

    ParametersVPRKpMidpoint{DT, TT, ET, ndims(equ), tab.s}(equ, tab, Δt, R, zero(TT), q̅, p̅)
end

@generated function compute_projection_vprk!(x::Vector{ST},
                q::SolutionVector{ST}, p::SolutionVector{ST}, λ::SolutionVector{ST},
                V::Vector{Vector{ST}}, U::Vector{Vector{ST}}, G::Vector{Vector{ST}},
                params::ParametersVPRKpMidpoint{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    # create temporary vectors
    q̃  = zeros(ST,D)
    # qm = zeros(ST,D)

    quote
        local t₀::TT = params.t̅
        local t₁::TT = params.t̅ + params.Δt
        local tₘ::TT = (t₀+t₁)/2
        local y::ST

        # copy x to λ and q̅
        for k in 1:D
            q[k] = x[D*(S+0)+k]
            λ[k] = x[D*(S+1)+k]
        end

        # compute U=λ
        U[1] .= λ
        U[2] .= λ

        # compute G=g(q,λ)
        for k in 1:D
            y = 0
            for j in 1:S
                y += params.tab.q.b[j] * V[j][k]
            end
            $q̃[k] = params.q̅[k] + 0.5 * params.Δt * y + params.Δt * params.R[1] * U[1][k]
            # $qm[k] = params.q̅[k] + 0.5 * params.Δt * y + 0.5 * params.Δt * params.R[1] * U[k,1] + 0.5 * params.Δt * params.R[2] * U[k,2]
        end

        # println("q̃mid = ", $q̃)
        # println("qmid = ", $qm)

        params.equ.g(tₘ, $q̃, λ, G[1])
        G[2] .= G[1]

        # compute p=ϑ(q)
        params.equ.ϑ(t₁, q, λ, p)
    end
end

"Compute stages of variational partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersVPRKpMidpoint{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    fcache = NonlinearFunctionCacheVPRK{ST}(D,S)
    pcache = NonlinearFunctionCacheVPRKprojection{ST}(D,S)

    function_stages = quote
        compute_stages!(x, $pcache.q, $pcache.p, $pcache.λ, $fcache.Q, $fcache.V, $pcache.U, $fcache.P, $fcache.F, $pcache.G, params)

        # compute b = - [P-AF-U]
        compute_rhs_vprk!(b, $fcache.P, $fcache.F, $pcache.G, params)

        # compute b = - [p-bF-G]
        compute_rhs_vprk_projection_p!(b, $pcache.p, $fcache.F, $pcache.G, D*(S+0), params)

        # compute b = - [q-bV-U]
        compute_rhs_vprk_projection_q!(b, $pcache.q, $fcache.V, $pcache.U, D*(S+1), params)

        compute_rhs_vprk_correction!(b, $fcache.V, params)
    end

    return function_stages
end

"Variational partitioned Runge-Kutta integrator."
struct IntegratorVPRKpMidpoint{DT, TT, PT <: ParametersVPRKpMidpoint{DT,TT},
                                       ST <: NonlinearSolver{DT},
                                       IT <: InitialGuessPODE{DT,TT}} <: AbstractIntegratorVPRKwProjection{DT,TT}

    params::PT
    solver::ST
    iguess::IT
end

function IntegratorVPRKpMidpoint(equation::ET, tableau::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: IODE{DT,TT}}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create params
    params = ParametersVPRKpMidpoint(equation, tableau, Δt)

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

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create integrator
    IntegratorVPRKpMidpoint{DT, TT, typeof(params), typeof(solver), typeof(iguess)}(params, solver, iguess)
end

equation(integrator::IntegratorVPRKpMidpoint) = integrator.params.equ
timestep(integrator::IntegratorVPRKpMidpoint) = integrator.params.Δt
tableau(integrator::IntegratorVPRKpMidpoint) = integrator.params.tab
nstages(integrator::IntegratorVPRKpMidpoint) = integrator.params.tab.s

function initial_guess!(int::IntegratorVPRKpMidpoint{DT,TT}, cache::IntegratorCacheVPRK{DT,TT}) where {DT,TT}
    for i in 1:nstages(int)
        evaluate!(int.iguess, cache.q, cache.p, cache.v, cache.f,
                              cache.q̅, cache.p̅, cache.v̅, cache.f̅,
                              cache.q̃, cache.ṽ,
                              tableau(int).q.c[i])

        for k in 1:ndims(int)
            int.solver.x[ndims(int)*(i-1)+k] = cache.ṽ[k]
        end
    end

    evaluate!(int.iguess, cache.q, cache.p, cache.v, cache.f,
                          cache.q̅, cache.p̅, cache.v̅, cache.f̅,
                          cache.q̃, one(TT))
    for k in 1:ndims(int)
        int.solver.x[ndims(int)*(nstages(int)+0)+k] = cache.q̃[k]
    end
    for k in 1:ndims(int)
        int.solver.x[ndims(int)*(nstages(int)+1)+k] = 0
    end
end

"Integrate ODE with variational partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorVPRKpMidpoint{DT,TT}, cache::IntegratorCacheVPRK{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, cache)

    # compute initial guess
    initial_guess!(int, cache)

    # reset cache
    reset!(cache, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

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

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

    # compute vector fields at internal stages and projection vector fields
    compute_stages!(int.solver.x,
                    cache.q̃, cache.p̃, cache.λ,
                    cache.Q, cache.V, cache.U,
                    cache.P, cache.F, cache.G, int.params)

    # compute unprojected solution
    update_solution!(int, cache)

    # add projection to solution
    project_solution!(int, cache, int.params.R)

    # copy solution to initial guess
    update!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f)

    # take care of periodic solutions
    cut_periodic_solution!(cache, equation(int).periodicity)
end


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

# Text / HTML objects

import Base: print, show, ==, hash

export HTML, @html_str

export HTML, Text

"""
`HTML(s)`: Create an object that renders `s` as html.

    HTML("
foo
")

You can also use a stream for large amounts of data:

    HTML() do io
      println(io, "
foo
")
    end
"""
mutable struct HTML{T}
    content::T
end

function HTML(xs...)
    HTML() do io
        for x in xs
            show(io, MIME"text/html"(), x)
        end
    end
end

show(io::IO, ::MIME"text/html", h::HTML) = print(io, h.content)
show(io::IO, ::MIME"text/html", h::HTML{<:Function}) = h.content(io)

"""
    @html_str -> Docs.HTML

Create an `HTML` object from a literal string.
"""
macro html_str(s)
    :(HTML($s))
end

function catdoc(xs::HTML...)
    HTML() do io
        for x in xs
            show(io, MIME"text/html"(), x)
        end
    end
end

export Text, @text_str

"""
`Text(s)`: Create an object that renders `s` as plain text.

    Text("foo")

You can also use a stream for large amounts of data:

    Text() do io
      println(io, "foo")
    end
"""
mutable struct Text{T}
    content::T
end

print(io::IO, t::Text) = print(io, t.content)
print(io::IO, t::Text{<:Function}) = t.content(io)
show(io::IO, t::Text) = print(io, t)

==(t1::T, t2::T) where {T<:Union{HTML,Text}} = t1.content == t2.content
hash(t::T, h::UInt) where {T<:Union{HTML,Text}} = hash(T, hash(t.content, h))

"""
    @text_str -> Docs.Text

Create a `Text` object from a literal string.
"""
macro text_str(s)
    :(Text($s))
end

function catdoc(xs::Text...)
    Text() do io
        for x in xs
            show(io, MIME"text/plain"(), x)
        end
    end
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
@testset "heuristic models" begin
    @test true == true
end


text:
using DSGE, Test

# Initialize model object
m = Model990(testing = true)
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:use_population_forecast, true)

# Forecast under historical rule
params = [α.value for α in m.parameters]
df = load_data(m, verbose = :none)

output_vars = [:histpseudo, :forecastpseudo, :irfpseudo]
out_hist = DSGE.forecast_one_draw(m, :mode, :none, output_vars, params, df, verbose = :none)

# Assign alternative policy
m <= Setting(:alternative_policy, DSGE.taylor93())

# Compare equilibrium conditions under historical and alternative rules
Γ0_hist, Γ1_hist, _ = eqcond(m)
Γ0_alt,  Γ1_alt,  _ = alternative_policy(m).eqcond(m)

eq_mp    = m.equilibrium_conditions[:eq_mp]
eq_other = setdiff(1:n_equilibrium_conditions(m), eq_mp)

@testset "Compare AltPolicy eqconds under hist and alt rules" begin
    @test Γ0_hist[eq_other, :] ≈ Γ0_alt[eq_other, :]
    @test Γ1_hist[eq_other, :] ≈ Γ1_alt[eq_other, :]
    @test !(Γ0_hist[eq_mp, :]  ≈ Γ0_alt[eq_mp, :])
    @test !(Γ1_hist[eq_mp, :]  ≈ Γ1_alt[eq_mp, :])
end

# Test error thrown if trying to run shockdec products
@testset "Check AltPolicy error thrown with shockdec products" begin
    @test_throws ErrorException forecast_one(m, :mode, :none, [:shockdecobs])
    @test_throws ErrorException forecast_one(m, :mode, :none, [:dettrendobs])
    @test_throws ErrorException forecast_one(m, :mode, :none, [:trendobs])
end

# Forecast under Taylor 93
out_alt = DSGE.forecast_one_draw(m, :mode, :none, output_vars, params, df, verbose = :none)

@testset "Check forecast under Taylor 93" begin
    @test out_hist[:histpseudo] ≈ out_alt[:histpseudo]
    @test !(out_hist[:forecastpseudo] ≈ out_alt[:forecastpseudo])
end


text:
## NamedArrays.jl
## (c) 2013--2018 David A. van Leeuwen

## Julia type that implements a drop-in replacement of Array with named dimensions.

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

module NamedArrays

using Requires
import InvertedIndices.Not
using Random
using Statistics
using SparseArrays
using DelimitedFiles
using DataStructures
using LinearAlgebra

export NamedArray, NamedVector, NamedMatrix, Name, Not

## type definition
include("namedarraytypes.jl")

export names, dimnames, defaultnames, setnames!, setdimnames!, array, enamerate

include("constructors.jl")
include("arithmetic.jl")
include("linearalgebra.jl")
include("base.jl")
include("changingnames.jl")
include("index.jl")
include("keepnames.jl")
include("names.jl")
include("rearrange.jl")
include("show.jl")
include("convert.jl")
include("iterators.jl")

function __init__()
    @require KahanSummation="8e2b3108-d4c1-50be-a7a2-16352aec75c3" begin

        # NOTE: KahanSummation does not support Julia 0.7 dims keyword argument at the moment

        ## rename a dimension
        function KahanSummation.cumsum_kbn(a::NamedArray, dims::Integer)
            NamedArrays.fan(KahanSummation.cumsum_kbn, "cumsum_kbn", a, dims)
        end
    end
end

end


text:
using Comonicon, PlutoUtils
Comonicon.install(PlutoUtils)


text:
using LazyArtifacts
using Lerche: Lerche, Lark, Transformer, @rule, @inline_rule
using CodecZlib: GzipDecompressorStream, GzipCompressorStream

export FileFormat, GzipFormat

# version of model zoo to grab files from
const zoo_version = "/Circuit-Model-Zoo-0.1.4"

#  by default don't transform tokens in parser
abstract type JuiceTransformer <: Transformer end
Lerche.visit_tokens(t::JuiceTransformer) = false

# file formats supported by this package
abstract type FileFormat end

struct GzipFormat <: FileFormat
    inner_format::FileFormat
end

# usual comment format for DIMACS-based file formats
const dimacs_comments = raw"""
    COMMENT : ("c" | "cc") (_WS /[^\n]/*)? (_NL | /$/)
    %ignore COMMENT
"""

# individual io functionalities

include("vtree_io.jl")
include("fnf_io.jl")
include("jlc_io.jl")
include("sdd_io.jl")
include("nnf_io.jl")
include("data_load.jl")
include("plot.jl")

# if no logic circuit file format is given on read, infer file format from extension

function file2logicformat(file)
    if endswith(file,".gz")
        file_inner, _ = splitext(file)
        format_inner = file2logicformat(file_inner)
        GzipFormat(format_inner)
    elseif endswith(file,".jlc")
        JlcFormat()
    elseif endswith(file,".sdd")
        SddFormat()
    elseif endswith(file,".nnf")
        NnfFormat()
    elseif endswith(file,".cnf") || endswith(file,".dnf")
        FnfFormat()
    else
        throw("Unknown file extension in $file: provide a file format argument")
    end
end

"""
    Base.read(file::AbstractString, ::Type{C}) where C <: LogicCircuit

Reads circuit from file; uses extension to detect format type, for example ".sdd" for SDDs.
"""
Base.read(file::AbstractString, ::Type{C}) where C <: LogicCircuit =
    read(file, C, file2logicformat(file))

Base.read(files::Tuple{AbstractString,AbstractString},
          ::Type{C}) where C <: StructLogicCircuit =
    read(files, C, (file2logicformat(files[1]), VtreeFormat()))

"""
    Base.write(file::AbstractString, circuit::LogicCircuit)

Writes circuit to file; uses file name extention to detect file format.
"""
Base.write(file::AbstractString, circuit::LogicCircuit) =
    write(file, circuit, file2logicformat(file))

"""
    Base.write(files::Tuple{AbstractString,AbstractString}, circuit::StructLogicCircuit)

Saves circuit and vtree to file.
"""
Base.write(files::Tuple{AbstractString,AbstractString},
           circuit::StructLogicCircuit) =
    write(files, circuit, (file2logicformat(files[1]), VtreeFormat()))

#  when asked to parse/read as `LogicCircuit`, default to `PlainLogicCircuit`

Base.parse(::Type{LogicCircuit}, args...) =
    parse(PlainLogicCircuit, args...)

Base.read(io::IO, ::Type{LogicCircuit},  args...) =
    read(io, PlainLogicCircuit,  args...)

Base.read(io::IO, ::Type{LogicCircuit}, f::GzipFormat) =
    # avoid method ambiguity
    read(io, PlainLogicCircuit, f)

#  when asked to parse/read as `StructLogicCircuit`, default to `PlainStructLogicCircuit`

Base.parse(::Type{StructLogicCircuit}, args...) =
    parse(PlainStructLogicCircuit, args...)

Base.read(io::IO, ::Type{StructLogicCircuit},  args...) =
    read(io, PlainStructLogicCircuit,  args...)

Base.read(io::IO, ::Type{StructLogicCircuit}, f::GzipFormat) =
    # avoid method ambiguity
    read(io, PlainStructLogicCircuit, f)

Base.read(ios::Tuple{IO,IO}, ::Type{StructLogicCircuit},  args...) =
    read(ios, PlainStructLogicCircuit,  args...)

Base.read(ios::Tuple{IO,IO}, ::Type{StructLogicCircuit}, f::Tuple{GzipFormat,VtreeFormat}) =
    # avoid method ambiguity
    read(ios, PlainStructLogicCircuit, f)

# copy read/write API for tuples of files

function Base.read(files::Tuple{AbstractString, AbstractString}, ::Type{C}, args...) where C <: StructLogicCircuit
    open(files[1]) do io1
        open(files[2]) do io2
            read((io1, io2), C, args...)
        end
    end
end

function Base.write(files::Tuple{AbstractString,AbstractString},
                    circuit::StructLogicCircuit, args...)
    open(files[1], "w") do io1
        open(files[2], "w") do io2
            write((io1, io2), circuit, args...)
        end
    end
end

# (de)compress Gzip streams

Base.read(io::IO, circuit_type, f::GzipFormat) =
    read(GzipDecompressorStream(io), circuit_type, f.inner_format)

Base.write(io::IO, circuit, f::GzipFormat) = begin
    iogz = GzipCompressorStream(io)
    write(iogz, circuit, f.inner_format)
    close(iogz)
end

Base.read(ios::Tuple{IO,IO}, circuit_type, f::Tuple{GzipFormat,VtreeFormat}) = begin
    ios = (GzipDecompressorStream(ios[1]), ios[2])
    formats = (f[1].inner_format, f[2])
    read(ios, circuit_type, formats)
end

Base.write(io::Tuple{IO,IO}, circuit, f::Tuple{GzipFormat,VtreeFormat}) = begin
    iosgz = (GzipCompressorStream(io[1]), io[2])
    formats = (f[1].inner_format, f[2])
    write(iosgz, circuit, formats)
    close(iosgz[1])
end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
module Delaunator

using OffsetArrays

const EPSILON = 2^-52 # eps()
const EDGE_STACK = OffsetVector{UInt32}(undef,0:511)

mutable struct DelaunatorData
    _triangles::OffsetVector{UInt32}
    _halfedges::OffsetVector{UInt32}
end

function delaunator!(points)
    n = length(points)
    coords = OffsetVector{Float64}(undef, 0:(n*2))

    @show n
    @show coords
    for i = 0:n-1
        p = points[i+1]
        @show p
        coords[2 * i] = defaultGetX(p)
        coords[2 * i + 1] = defaultGetY(p)
    end
    @show "1"
    #return new Delaunator(coords);

    n = div(n,2) # coords.length >> 1
    #if (n > 0 && typeof coords[0] !== 'number') throw new Error('Expected coords to contain numbers.')

    # arrays that will store the triangulation graph
    maxTriangles = max(2 * n - 5, 0)
    @show typeof(maxTriangles),  maxTriangles
    _triangles =  OffsetVector{UInt32}(undef,0:(maxTriangles * 3)-1)
    _halfedges = OffsetVector{UInt32}(undef,0:(maxTriangles * 3)-1)

    # temporary arrays for tracking the edges of the advancing convex hull
    _hashSize = ceil(Int,sqrt(n))
    hullPrev = OffsetVector{UInt32}(undef,0:n-1) # edge to prev edge
    hullNext = OffsetVector{UInt32}(undef,0:n-1)# edge to next edge
    hullTri = OffsetVector{UInt32}(undef,0:n-1) # edge to adjacent triangle
    hullHash = fill!(OffsetVector{Int32}(undef,0:n-1), -1) # angular edge hash

    # temporary arrays for sorting points
    _ids = OffsetVector{UInt32}(undef,0:n-1)
    _dists = OffsetVector{Float64}(undef,0:n-1)

    #this.update()

    #n = n/2 # coords.length >> 1 # n/2

    # populate an array of point indices; calculate input data bbox
    minX = Inf
    minY = Inf
    maxX = -Inf
    maxY = -Inf

    for i = 0:n-1
        x = coords[2 * i]
        y = coords[2 * i + 1]
        x < minX && (minX = x)
        y < minY && (minY = y)
        x > maxX && (maxX = x)
        y > maxY && (maxY = y)
        _ids[i] = i
    end
    cx = (minX + maxX) / 2
    cy = (minY + maxY) / 2
    @show cx, cy

    minDist = Inf
    i0 = 0
    i1 = 0
    i2 = 0

    # pick a seed point close to the center
    for i = 0:n-1
        d = dist(cx, cy, coords[2 * i], coords[2 * i + 1])
        @show d
        if d < minDist
            i0 = i
            @show minDist
            minDist = d
        end
    end
    i0x = coords[2 * i0]
    i0y = coords[2 * i0 + 1]
    @show i0x,i0y
    minDist = Inf

    # find the point closest to the seed
    for i = 0:n-1
        i == i0 && continue
        @show minDist
        d = dist(i0x, i0y, coords[2 * i], coords[2 * i + 1])
        if d < minDist && d > 0
            i1 = i
            minDist = d
        end
    end
    i1x = coords[2 * i1]
    i1y = coords[2 * i1 + 1]

    @show i1x, i1y

    minRadius = Inf

    # find the third point which forms the smallest circumcircle with the first two
    for i = 1:n-1
        i == i0 || i == i1 && continue
        r = circumradius(i0x, i0y, i1x, i1y, coords[2 * i], coords[2 * i + 1])
        if r < minRadius
            i2 = i
            minRadius = r
        end
    end
    i2x = coords[2 * i2]
    i2y = coords[2 * i2 + 1]

    if minRadius == Inf
        # order collinear points by dx (or dy if all x are identical)
        # and return the list as a hull
        for i = 0:n-1
            cxv = (coords[2 * i] - coords[0])
            _dists[i] = !iszero(cxv) ? cxv : (coords[2 * i + 1] - coords[1])
        end
        quicksort(_ids, _dists, 0, n - 1)
        hull = OffsetVector{UInt32}(undef,0:n-1)
        j = 0
        d0 = -Inf
        for i = 0:n-1
            id = _ids[i]
            if _dists[id] > d0
                hull[j] = id
                j += 1
                d0 = _dists[id]
            end
        end
        hull = resize!(hull, j+1)
        triangles = OffsetVector{UInt32}(undef, 0:0)
        halfedges = OffsetVector{UInt32}(undef, 0:0)
        return hull # TODO
    end

    # swap the order of the seed points for counter-clockwise orientation
    if orient(i0x, i0y, i1x, i1y, i2x, i2y)
        i = i1
        x = i1x
        y = i1y
        i1 = i2
        i1x = i2x
        i1y = i2y
        i2 = i
        i2x = x
        i2y = y
    end

    center = circumcenter(i0x, i0y, i1x, i1y, i2x, i2y);
    _cx = center[1]
    _cy = center[2]

    for i = 0:n-1
        _dists[i] = dist(coords[2 * i], coords[2 * i + 1], _cx, _cy)
    end

    # sort the points by distance from the seed triangle circumcenter
    quicksort(_ids, _dists, 0, n - 1)

    # set up the seed triangle as the starting hull
    _hullStart = i0
    hullSize = 3

    hullNext[i0] = hullPrev[i2] = i1
    hullNext[i1] = hullPrev[i0] = i2
    hullNext[i2] = hullPrev[i1] = i0

    hullTri[i0] = 0
    hullTri[i1] = 1
    hullTri[i2] = 2

    fill!(hullHash,-1)
    hullHash[_hashKey(i0x, i0y,_cx,_cy,_hashSize)] = i0
    hullHash[_hashKey(i1x, i1y,_cx,_cy,_hashSize)] = i1
    hullHash[_hashKey(i2x, i2y,_cx,_cy,_hashSize)] = i2

    trianglesLen = _addTriangle(_triangles, _halfedges, 0, i0, i1, i2, -1, -1, -1)

    xp = 0.0
    yp = 0.0
    for k = 0:length(_ids)-1
        i = _ids[k]
        x = coords[2 * i]
        y = coords[2 * i + 1]

        # skip near-duplicate points
        if k > 0 && abs(x - xp) <= EPSILON && abs(y - yp) <= EPSILON
            continue
        end
        xp = x
        yp = y

        # skip seed triangle points
        if i === i0 || i === i1 || i === i2
            continue
        end

        # find a visible edge on the convex hull using edge hash
        start = 0
        key = _hashKey(x, y, _cx, _cy, _hashSize)
        @show key
        for j = 0:_hashSize-1
            start = hullHash[(key + j) % _hashSize]
            start != -1 && start != hullNext[start] && break
        end

        start = hullPrev[start]
        e = start
        q = hullNext[e]
        while !orient(x, y, coords[2 * e], coords[2 * e + 1], coords[2 * q], coords[2 * q + 1])
            e = q
            if e == start
                e = -1
                break
            end
            q = hullNext[e]
        end

        e == -1 && continue # likely a near-duplicate point; skip it

        # add the first triangle from the point
        t = _addTriangle(e, i, hullNext[e], -1, -1, hullTri[e])

        # recursively flip triangles from the point until they satisfy the Delaunay condition
        hullTri[i] = _legalize(t + 2)
        hullTri[e] = t # keep track of boundary triangles on the hull
        hullSize += 1

        # walk forward through the hull, adding more triangles and flipping recursively
        n = hullNext[e]
        q = hullNext[n]
        while orient(x, y, coords[2 * n], coords[2 * n + 1], coords[2 * q], coords[2 * q + 1])
            t = _addTriangle(n, i, q, hullTri[i], -1, hullTri[n])
            hullTri[i] = _legalize(t + 2)
            hullNext[n] = n # mark as removed
            hullSize -= 1
            n = q
            q = hullNext[n]
        end

        # walk backward from the other side, adding more triangles and flipping
        if e == start
            q = hullPrev[e]
            while orient(x, y, coords[2 * q], coords[2 * q + 1], coords[2 * e], coords[2 * e + 1])
                t = _addTriangle(q, i, e, -1, hullTri[e], hullTri[q])
                _legalize(t + 2)
                hullTri[q] = t
                hullNext[e] = e # mark as removed
                hullSize -= 1
                e = q
                q = hullPrev[e]
            end
        end

        # update the hull indices
        _hullStart = hullPrev[i] = e
        hullNext[e] = hullPrev[n] = i
        hullNext[i] = n

        # save the two new edges in the hash table
        hullHash[_hashKey(x, y,_cx,_cy,_hashSize)] = i
        hullHash[_hashKey(coords[2 * e], coords[2 * e + 1],_cx,_cy,_hashSize)] = e
    end

    hull = OffsetVector{UInt32}(undef,0:hullSize-1)
    for i = 0:hullSize-1
        e = _hullStart
        hull[i] = e
        e = hullNext[e]
    end

    # trim typed triangle mesh arrays
    return (resize!(_triangles, trianglesLen), resize!(_halfedges, trianglesLen))
end

function _hashKey(x, y, cx, cy, _hashSize)
    return floor(Int, pseudoAngle(x - cx, y - cy) * _hashSize) % _hashSize
end

function _legalize(a, triangles, halfedges, coords)

    i = 0
    ar = 0

    #recursion eliminated with a fixed-size stack
    while true
        b = halfedges[a]

         # if the pair of triangles doesn't satisfy the Delaunay condition
         # (p1 is inside the circumcircle of [p0, pl, pr]), flip them,
         # then do the same check/flip recursively for the new pair of triangles
         #
         #           pl                    pl
         #          /||\                  /  \
         #       al/ || \bl            al/    \a
         #        /  ||  \              /      \
         #       /  a||b  \    flip    /___ar___\
         #     p0\   ||   /p1   =>   p0\---bl---/p1
         #        \  ||  /              \      /
         #       ar\ || /br             b\    /br
         #          \||/                  \  /
         #           pr                    pr
         #
        a0 = a - a % 3
        ar = a0 + (a + 2) % 3

        if b == -1 #convex hull edge
            i == 0 && break
            i -= 1
            a = EDGE_STACK[i]
            continue
        end

        b0 = b - b % 3
        al = a0 + (a + 1) % 3
        bl = b0 + (b + 2) % 3

        p0 = triangles[ar]
        pr = triangles[a]
        pl = triangles[al]
        p1 = triangles[bl]

        illegal = inCircle(coords[2 * p0], coords[2 * p0 + 1],
                           coords[2 * pr], coords[2 * pr + 1],
                           coords[2 * pl], coords[2 * pl + 1],
                           coords[2 * p1], coords[2 * p1 + 1])

        if illegal
            triangles[a] = p1
            triangles[b] = p0

            hbl = halfedges[bl]

            # edge swapped on the other side of the hull (rare); fix the halfedge reference
            if hbl == -1
                e = _hullStart
                # CAUTION do-while
                if hullTri[e] == bl
                    hullTri[e] = a
                else
                    e = hullPrev[e]
                    while e != _hullStart
                        if this.hullTri[e] == bl
                            this.hullTri[e] = a
                            break
                        end
                        e = hullPrev[e]
                    end
                end
            end
            _link(halfedges, a, hbl)
            _link(halfedges, b, halfedges[ar])
            _link(halfedges, ar, bl)

            br = b0 + (b + 1) % 3

            # don't worry about hitting the cap: it can only happen on extremely degenerate input
            if i < EDGE_STACK.length
                EDGE_STACK[i] = br
                i += 1
            end
        else
            i == 0 && break
            i -= 1
            a = EDGE_STACK[i]
        end
    end

    return ar
end

function _link(_halfedges, a, b)
    @show a, b
    _halfedges[a] = b == -1 ? typemax(UInt32) : b
    b != -1 && (_halfedges[b] = a)
end

# add a new triangle given vertex indices and adjacent half-edge ids
function _addTriangle(_triangles, _halfedges, t, i0, i1, i2, a, b, c)

    _triangles[t] = i0
    _triangles[t + 1] = i1
    _triangles[t + 2] = i2

    _link(_halfedges, t, a)
    _link(_halfedges, t + 1, b)
    _link(_halfedges, t + 2, c)

    return t + 3
end

# monotonically increases with real angle, but doesn't need expensive trigonometry
function pseudoAngle(dx, dy)
    p = dx / (abs(dx) + abs(dy))
    return (dy > 0 ? 3 - p : 1 + p) / 4 # [0..1]
end

function dist(ax, ay, bx, by)
    dx = ax - bx
    dy = ay - by
    return dx * dx + dy * dy
end

# return 2d orientation sign if we're confident in it through J. Shewchuk's error bound check
function orientIfSure(px, py, rx, ry, qx, qy)
    l = (ry - py) * (qx - px)
    r = (rx - px) * (qy - py)
    return abs(l - r) >= 3.3306690738754716e-16 * abs(l + r) ? l - r : 0
end

# a more robust orientation test that's stable in a given triangle (to fix robustness issues)
function orient(rx, ry, qx, qy, px, py)
    #return (orientIfSure(px, py, rx, ry, qx, qy) ||
    #    orientIfSure(rx, ry, qx, qy, px, py) ||
    #    orientIfSure(qx, qy, px, py, rx, ry)) < 0;
    return orientIfSure(px, py, rx, ry, qx, qy) < 0 ||
           orientIfSure(rx, ry, qx, qy, px, py) < 0 ||
           orientIfSure(qx, qy, px, py, rx, ry) < 0
end

function inCircle(ax, ay, bx, by, cx, cy, px, py)
    dx = ax - px
    dy = ay - py
    ex = bx - px
    ey = by - py
    fx = cx - px
    fy = cy - py

    ap = dx * dx + dy * dy
    bp = ex * ex + ey * ey
    cp = fx * fx + fy * fy

    return dx * (ey * cp - bp * fy) -
           dy * (ex * cp - bp * fx) +
           ap * (ex * fy - ey * fx) < 0
end

function circumradius(ax, ay, bx, by, cx, cy)
    dx = bx - ax
    dy = by - ay
    ex = cx - ax
    ey = cy - ay

    bl = dx * dx + dy * dy
    cl = ex * ex + ey * ey
    d = 0.5 / (dx * ey - dy * ex)

    x = (ey * bl - dy * cl) * d
    y = (dx * cl - ex * bl) * d

    return x * x + y * y
end

function circumcenter(ax, ay, bx, by, cx, cy)
    dx = bx - ax
    dy = by - ay
    ex = cx - ax
    ey = cy - ay

    bl = dx * dx + dy * dy
    cl = ex * ex + ey * ey
    d = 0.5 / (dx * ey - dy * ex)

    x = ax + (ey * bl - dy * cl) * d
    y = ay + (dx * cl - ex * bl) * d

    return (x, y)
end

function quicksort(ids, dists, left, right)
    if right - left <= 20
        i = i = left + 1
        while i <= right
            temp = ids[i]
            tempDist = dists[temp]
            j = i - 1
            while j >= left && dists[ids[j]] > tempDist
                ids[j + 1] = ids[j]
                j -= 1
            end
            ids[j + 1] = temp
            i += 1
        end
    else
        median = (left + right) >> 1
        i = left + 1
        j = right
        swap(ids, median, i)
        dists[ids[left]] > dists[ids[right]] && swap(ids, left, right)
        dists[ids[i]] > dists[ids[right]] && swap(ids, i, right)
        dists[ids[left]] > dists[ids[i]] && swap(ids, left, i)

        temp = ids[i]
        tempDist = dists[temp]
        while true
            i += 1
            while dists[ids[i]] < tempDist
                i += 1
            end
            j -= 1
            while dists[ids[j]] > tempDist
                j -= 1
            end
            j < i && break
            swap(ids, i, j)
        end
        ids[left + 1] = ids[j]
        ids[j] = temp

        if (right - i + 1 >= j - left)
            quicksort(ids, dists, i, right)
            quicksort(ids, dists, left, j - 1)
        else
            quicksort(ids, dists, left, j - 1)
            quicksort(ids, dists, i, right)
        end
    end
end

function swap(arr, i, j)
    tmp = arr[i]
    arr[i] = arr[j]
    arr[j] = tmp
end

function defaultGetX(p)
    return p[1]
end

function defaultGetY(p)
    return p[2]
end

end # module


text:
using SafeTestsets

@safetestset "Inference Tests" begin include("inferencetests.jl") end
@safetestset "Descriptive Statistics Tests" begin include("descriptivestatisticstests.jl") end
@safetestset "Regression Tests" begin include("regressiontests.jl") end
@safetestset "Distributions Tests" begin include("distributionstests.jl") end

max_issues_repo_name: persson-io/Stochastics

max_forks_repo_name: persson-io/Stochastics


text:
using Plots
using DelimitedFiles

a = readdlm("vqe.txt")

using GR
heatmap(a, xlog=true)


text:
function solve_exp_ode(
    M::AbstractManifold,
    p,
    X;
    basis::AbstractBasis=DefaultOrthonormalBasis(),
    solver=AutoVern9(Rodas5()),
    backend=default_differential_backend(),
    kwargs...,
)
    n = length(p)
    iv = SVector{n}(1:n)
    ix = SVector{n}((n + 1):(2n))
    u0 = allocate(p, 2n)
    u0[iv] .= X
    u0[ix] .= p

    function exp_problem(u, params, t)
        M = params[1]
        dx = u[iv]
        p = u[ix]
        ddx = allocate(u, Size(n))
        du = allocate(u)
        Γ = christoffel_symbols_second(M, p, basis; backend=backend)
        @einsum ddx[k] = -Γ[k, i, j] * dx[i] * dx[j]
        du[iv] .= ddx
        du[ix] .= dx
        return Base.convert(typeof(u), du)
    end

    params = (M,)
    prob = ODEProblem(exp_problem, u0, (0.0, 1.0), params)
    sol = solve(prob, solver; kwargs...)
    q = sol.u[1][(n + 1):(2 * n)]
    return q
end


text:
@testset "390.elimination-game.jl" begin
    for (i, j) in
        zip([9, 246, 12, 634, 24, 32, 357, 4734], [6, 120, 6, 384, 14, 22, 118, 1920])
        @test last_remaining(i) == j
    end
end


text:
########################
#   Type Definition    #
########################

abstract type AbstractGKGᵀ end

mutable struct symGKGᵀ <: AbstractGKGᵀ

    support::LinRange{Float64}

    ℓ::Float64

    g::PhysicalTransferFunction

    X::Array{Array{Float64, 1}}

    mass::Float64

    accretion::Float64

    wavelength::Float64

    f::AbstractExtrapolation

    parallelevaluation::Bool

end

##############
#    Info    #
##############

function Base.show(io::IO, k::symGKGᵀ)

    print(io, "Symmetric GKGᵀ product (0.1.52)\n",
              @sprintf("parallel evaluation set to %s\n", k.parallelevaluation),
              @sprintf("mass is %e\n", k.mass),
              @sprintf("accretion is %2.2f\n", k.accretion),
              @sprintf("wavelength is %4.1f", k.wavelength))

end

###############
# Constructor #
###############

function symGKGᵀ(mass::Float64=1e8, accretion::Float64=0.5, wavelength::Float64=5000.0, parallelevaluation::Bool=false)

    gₖ = PhysicalTransferFunction(mass=mass, accretion=accretion, wavelength=wavelength)

    numberofsegments = 31 # crucial parameter of approximation

    Xₖ = pw(collect(0.0:0.01:width(gₖ)), gₖ.(collect(0.0:0.01:width(gₖ))), numberofsegments)

    support = LinRange(-30.0, 30.0, 271) # this must also be looked at again

    y = pmap(t -> convkernel(; Xₖ = Xₖ, Xₗ = Xₖ, tᵢ = t, tⱼ = 0.0, ℓ = 1.0), support, distributed = parallelevaluation)

    symGKGᵀ(support, 1.0, gₖ, Xₖ, mass, accretion, wavelength, CubicSplineInterpolation(support, y, extrapolation_bc=0.0), parallelevaluation)

end

##############
# Evaluation #
##############

function (k::symGKGᵀ)(tᵢ, tⱼ, ℓ)

    if ℓ !== k.ℓ # in this case recompute

        k.ℓ = ℓ

        y = pmap(t -> convkernel(; Xₖ = k.X, Xₗ = k.X, tᵢ = t, tⱼ = 0.0, ℓ = k.ℓ), k.support, distributed = k.parallelevaluation)

        k.f = CubicSplineInterpolation(k.support, y, extrapolation_bc=0.0)

    end

    k.f(tᵢ - tⱼ)

end


text:
using IcanProj
using Base.Test

##########################################################################

# Code for preparing data to vic and running simulations

# 1. Prepare metadata for the experiment catchments

# include("metadata_script.jl")

# 2. Prepare vic input files using the metadata

# include("vic_setup_data.jl")

# 3. Run vic simulations

# include("run_vic.jl")

##########################################################################

# Code for tranforming vic raw output files into usable format

# 1. Average vic output and store in dataframes

# include("../src/process_vic_output.jl")

##########################################################################

# Plotting scripts

# include("plots_script.jl")


text:
export batched_tr, batched_tr!

"""
    batched_tr(A::AbstractArray{T, 3}) where T

Batched version of trace.
"""
batched_tr(A::AbstractArray{T, 3}) where T = batched_tr!(A, fill!(similar(A, (size(A, 3), )), 0))

function batched_tr!(A::AbstractArray{T, 3}, B::AbstractVector{T}) where T
    @boundscheck size(A, 1) == size(A, 2) || error("Expect a square matrix")
    @boundscheck size(A, 3) == length(B) || error("Batch size mismatch")
    @inbounds for k in 1:size(A, 3), i in 1:size(A, 1)
        B[k] += A[i, i, k]
    end
    return B
end


text:
module ADChain

using ChainRulesCore
using ChainRules
using GeneralizedGenerated
using Flux

export gradient

include("lib.jl")
include("closure_pullback.jl")
include("pass.jl")

end # module


text:
@test GeneralizedMetropolisHastings._numsamplesinchain(Val{:burnin},10,20,5) == (50,50)
@test GeneralizedMetropolisHastings._numsamplesinchain(Val{:main},10,20,5) == (100,0)
@test GeneralizedMetropolisHastings._numsamplesinchain(Val{:all},10,20,5) == (151,51)

mhpolicy1 = policy(:mh,1)
mhpolicy2 = policy(:mh,1,store=:all)
mhpolicy3 = policy(:mh,1,store=:burnin,initialize=:default)

@test GeneralizedMetropolisHastings._storeduring(:burnin,mhpolicy1) == false
@test GeneralizedMetropolisHastings._storeduring(:burnin,mhpolicy2) == true
@test GeneralizedMetropolisHastings._storeduring(:burnin,mhpolicy3) == true
@test GeneralizedMetropolisHastings._storeduring(:main,mhpolicy1) == true
@test GeneralizedMetropolisHastings._storeduring(:main,mhpolicy2) == true
@test GeneralizedMetropolisHastings._storeduring(:main,mhpolicy3) == false

type TestMHRunner <: AbstractMHRunner
    numburnin::Int
    numiterations::Int
    policy::MHRuntimePolicy
end

#create samplerstates
mhsamplerstate1 = samplerstate(rsampler1,smhnprops1,Float64,Float64)
mhsamplerstate2 = samplerstate(rsampler1,gmhnprops1,Float64,Float64)

#create all required objects
mhrunner1 = TestMHRunner(rnburnin1,rniter1,mhpolicy1)
mhrunner2 = TestMHRunner(rnburnin1,rniter1,mhpolicy2)
mhrunner3 = TestMHRunner(rnburnin1,rniter1,mhpolicy3)
mhindicator1,mhtunerstate1,mhchain1 = GeneralizedMetropolisHastings.createcommon(mhrunner1,rtuner1,rnparas1,smhnprops1,smhnprops1)
mhindicator2,mhtunerstate2,mhchain2 = GeneralizedMetropolisHastings.createcommon(mhrunner2,rtuner2,rnparas1,gmhnprops1,gmhnprops1)

@test numproposals(mhindicator1) == smhnprops1 && numsamples(mhindicator1) == smhnprops1
@test numproposals(mhindicator2) == gmhnprops1 && numsamples(mhindicator2) == gmhnprops1
@test indicatorsamples(mhindicator1) == [0,0]
@test indicatorsamples(mhindicator2) == [0,0,0]
@test numtunesteps(mhtunerstate1) == 2
@test numtunesteps(mhtunerstate2) == 1
@test numparas(mhchain1) == rnparas1 && numsamples(mhchain1) == rniter1
@test numparas(mhchain2) == rnparas1 && numsamples(mhchain2) == (rnburnin1 + rniter1)*gmhnprops1 + 1

#test the initialize function
#prepare some prior values
srand(rinitseed1) ; rpriorvals1 = initvalues(trait(:initialize,:prior),rparas1,mhpolicy1.sampletype)
srand(rinitseed1+1) ; rpriorvals2 = initvalues(trait(:initialize,:prior),rparas1,mhpolicy1.sampletype)
#copy (for testing) the first values of the chaing
mcv1 = [mhchain1.values[1]]
mcll1 = mhchain1.loglikelihood[1]

#initialilze but without storing the first value in the chain
srand(rinitseed1) ; initialize!(mhrunner1,rmodel1,mhsamplerstate1,mhchain1,false)
mhll1 = [loglikelihood(rmodel1,evaluate!(rmodel1,rpriorvals1))]
testproposals(mhsamplerstate1,rpriorvals1,mhll1) #the value is firt generated in the proposals
testfrom(mhsamplerstate1,rpriorvals1,mhll1) #then copied to the from field
@test_approx_eq mean(mhsamplerstate1.density.distribution) rpriorvals1 #and the normal distribution of the mh sampler should also have been conditioned
teststore(mhchain1,mcv1,mcll1,1,0,0) #finally, the value has not been stored inside the chain

srand(rinitseed1+1) ; initialize!(mhrunner1,rmodel1,mhsamplerstate1,mhchain1,true)
mhll2 = [loglikelihood(rmodel1,evaluate!(rmodel1,rpriorvals2))]
testproposals(mhsamplerstate1,rpriorvals2,mhll2) #same tests as above
testfrom(mhsamplerstate1,rpriorvals2,mhll2)
@test_approx_eq mean(mhsamplerstate1.density.distribution) rpriorvals2 #and the normal distribution of the mh sampler should also have been conditioned
teststore(mhchain1,mhsamplerstate1.from.values,mhsamplerstate1.from.loglikelihood[1],1,0,1) #this time, the value should have been copied to the chain


text:
using Screws

m = 0.2
Iz = Matrix3(Diagonal([0.5, 0.5, 0.0375]))
r_A = Vector3(0.75,0.0,0.0)
cg = Vector3(0.45,0.0,0.0)
n = Wrench(ĵ, r_A)
v = Twist(-5.0*k̂, ô)
I = spi(m,Iz,cg)
e = 0
s = Twist(k̂, ô)
T = I*s/dot(s,I*s)
Λ = s*s'/dot(s,I*s)
J = (1+e)*dot(n,v)/dot(n,Λ*n)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ngspice_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ngspice")
JLLWrappers.@generate_main_file("ngspice", UUID("57e080ce-e004-5ece-a144-c95753423b2a"))
end  # module ngspice_jll


text:
empty!(Base.DEPOT_PATH)
push!(Base.DEPOT_PATH, mktempdir(; cleanup = true))

import Pkg
import Logging
import TOML

Pkg.add(; name = "Coverage", uuid = "a2441757-f6aa-5fb2-8edb-039e3f45d037", version = "1")
Pkg.precompile()

import Coverage

function process_folders()
    # `Coverage.process_folder` will have a LOT of `@info` statements that will make the log
    # way too long. So before we run `Coverage.process_folder`, we disable logging for `@info`
    # statements. After we run `Coverage.process_folder`, we re-enable logging for `@info`
    # statements.
    Logging.disable_logging(Logging.Info)
    fcs_base   = Coverage.process_folder("base");
    fcs_stdlib = Coverage.process_folder("stdlib");
    Logging.disable_logging(Logging.Debug)

    fcs = Coverage.merge_coverage_counts(
        fcs_base,
        fcs_stdlib,
    );

    return fcs
end

function get_external_stdlib_names(stdlib_dir::AbstractString)
    filename_list = filter(x -> isfile(joinpath(stdlib_dir, x)), readdir(stdlib_dir))
    # find all of the files like `Pkg.version`, `Statistics.version`, etc.
    regex_matches_or_nothing = match.(Ref(r"^([\w].*?)\.version$"), filename_list)
    regex_matches = filter(x -> x !== nothing, regex_matches_or_nothing)
    # get the names of the external stdlibs, like `Pkg`, `Statistics`, etc.
    external_stdlib_names = only.(regex_matches)
    unique!(external_stdlib_names)
    sort!(external_stdlib_names)
    @info "# Begin list of external stdlibs"
    for (i, x) in enumerate(external_stdlib_names)
        @info "$(i). $(x)"
    end
    @info "# End list of external stdlibs"
    return external_stdlib_names
end

function get_external_stdlib_prefixes(stdlib_dir::AbstractString)
    external_stdlib_names = get_external_stdlib_names(stdlib_dir)
    prefixes_1 = joinpath.(Ref(stdlib_dir), external_stdlib_names, Ref(""))
    prefixes_2 = joinpath.(Ref(stdlib_dir), string.(external_stdlib_names, Ref("-")))
    prefixes = vcat(prefixes_1, prefixes_2)
    unique!(prefixes)
    sort!(prefixes)
    # example of what `prefixes` might look like:
    # 4-element Vector{String}:
    # "stdlib/Pkg-"
    # "stdlib/Pkg/"
    # "stdlib/Statistics-"
    # "stdlib/Statistics/"
    return prefixes
end

function print_coverage_summary(fc::Coverage.FileCoverage)
    cov_lines, tot_lines = Coverage.get_summary(fc)
    if cov_lines == tot_lines == 0
        cov_pct = 0
    else
        cov_pct = floor(Int, cov_lines/tot_lines * 100)
    end
    pad_1 = 71
    pad_2 = 15
    pad_3 = 15
    col_1 = rpad(fc.filename, pad_1)
    col_2 = rpad(string(cov_pct, " %"), pad_2)
    col_3 = string(
        rpad(string(cov_lines), pad_3),
        string(tot_lines),
    )
    @info "$(col_1) $(col_2) $(col_3)"
    return nothing
end

function print_coverage_summary(
        fcs::Vector{Coverage.FileCoverage}, description::AbstractString,
    )
    cov_lines, tot_lines = Coverage.get_summary(fcs)
    if cov_lines == tot_lines == 0
        cov_pct = 0
    else
        cov_pct = floor(Int, cov_lines/tot_lines * 100)
    end
    @info "$(description): $(cov_pct)% ($(cov_lines)/$(tot_lines))"
    return (; cov_pct)
end

function buildkite_env(name::String)
    value = String(strip(ENV[name]))
    if isempty(value)
        throw(ErrorException("environment variable $(name) is empty"))
    end
    return value
end

function buildkite_env(name_1::String, name_2::String, default::String)
    value_1 = String(strip(ENV[name_1]))
    value_2 = String(strip(ENV[name_2]))
    !isempty(value_1) && return value_1
    !isempty(value_2) && return value_2
    return default
end

function buildkite_branch_and_commit()
    branch = buildkite_env("BUILDKITE_BRANCH")
    commit = buildkite_env("BUILDKITE_COMMIT")
    head_rev_parse = String(strip(read(`git rev-parse HEAD`, String)))
    if strip(commit) == "HEAD"
        commit = head_rev_parse
    end
    if commit !== head_rev_parse
        msg = "mismatch"
        @error msg commit head_rev_parse
        throw(ErrorException(msg))
    end
    if !occursin(r"^[a-f0-9]{40}$", commit)
        msg = "BUILDKITE_COMMIT does not look like a long commit SHA"
        @error msg commit
        throw(ErrorException(msg))
    end
    return (; branch, commit)
end

function codecov_buildkite_add_local_to_kwargs()
    branch, commit = buildkite_branch_and_commit()
    kwargs = Coverage.Codecov.set_defaults(
        Dict();
        branch,
        commit,
    )
    return kwargs
end

function coveralls_buildkite_query_git_info()
    branch, commit = buildkite_branch_and_commit()
    remote_name  = "origin"
    remote       = buildkite_env("BUILDKITE_REPO")
    message      = buildkite_env("BUILDKITE_MESSAGE")
    author_name  = buildkite_env(
        "BUILDKITE_BUILD_AUTHOR",
        "BUILDKITE_BUILD_CREATOR",
        "",
    )
    author_email = buildkite_env(
        "BUILDKITE_BUILD_AUTHOR_EMAIL",
        "BUILDKITE_BUILD_CREATOR_EMAIL",
        "",
    )
    remotes = [
        Dict(
            "name"  => remote_name,
            "url"   => remote,
        )
    ]
    head = Dict(
        "id"                => commit,
        "author_name"       => author_name,
        "author_email"      => author_email,
        "committer_name"    => author_name,
        "committer_email"   => author_email,
        "message"           => message,
    )
    git_info = Dict(
        "branch"  => branch,
        "remotes" => remotes,
        "head"    => head,
    )
    return git_info
end

const fcs = process_folders()

# Only include source code files. Exclude test files, benchmarking files, etc.
filter!(fcs) do fc
    occursin(r"^base\/", fc.filename) || occursin("/src/", fc.filename)
end;

# Exclude all external stdlibs (stdlibs that live in external repos).
const external_stdlib_prefixes = get_external_stdlib_prefixes("stdlib")
filter!(fcs) do fc
    all(x -> !startswith(fc.filename, x), external_stdlib_prefixes)
end;

# Exclude all stdlib JLLs (stdlibs of the form `stdlib/*_jll/`).
filter!(fcs) do fc
    !occursin(r"^stdlib\/[A-Za-z0-9]*?_jll\/", fc.filename)
end;

sort!(fcs; by = fc -> fc.filename);

print_coverage_summary.(fcs);
const total_cov_pct = print_coverage_summary(fcs, "Total").cov_pct

let
    git_info = coveralls_buildkite_query_git_info()
    @info "" git_info
    @info "" git_info["branch"]
    @info "" git_info["head"]

    # In order to upload to Coveralls, you need to have the `COVERALLS_TOKEN` environment variable defined.
    Coverage.Coveralls.submit_local(fcs, git_info)
end

let
    kwargs = codecov_buildkite_add_local_to_kwargs()
    @info "" kwargs

    # In order to upload to Codecov, you need to have the `CODECOV_TOKEN` environment variable defined.
    Coverage.Codecov.submit_generic(fcs, kwargs)
end

if total_cov_pct < 50
    msg = string(
        "The total coverage is less than 50%. This should never happen, ",
        "so it means that something has probably gone wrong with the code coverage job.",
    )
    @error msg total_cov_pct
    throw(ErrorException(msg))
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
#numerical 1d integration

#integrate func below
f(x) =sin(x)
#over
xmin = -π
xmax = π
#with number of divisions (linear)
n=1000

#Trapeziod rule
Δ = (xmax - xmin)/n
xseq = range(xmin,xmax,length = n+1)
s_trapeziod = Δ/2 *(sum(f.(xseq))+sum(f.(xseq[2:end-1])))

#error of trapeziod estimation for constinuous second differentiable f
#e_trapeziod = (xmax - xmin)^3/12*n^2 * max(abs(f))
# see this slide for more info
# https://slideplayer.com/slide/6419579/

#Mid-point rule
s_Midp = sum(f.(xseq[2:end] - xseq[1:end-1]))*Δ

# n has to be even, t.s., n+1 is odd, there are n+1 terms to play with
iseven(n)

s_Simpson = Δ/3 * sum([[1]; repeat([4; 2; 4; 1], Int(n/4))] .* f.(xseq))

max_stars_repo_path: 1D Num Integration.jl

max_issues_repo_path: 1D Num Integration.jl

max_forks_repo_path: 1D Num Integration.jl


text:
mutable struct MiniMaxSampleSize <: Objective
    λ::Real
    prior::Prior
end

function (objective::MiniMaxSampleSize)(design::AbstractDesign)

    ess = SampleSize(objective.prior)
    return (1 - objective.λ)*maximum(n(design)) + λ*ess(design)
end

function update(obj::MiniMaxSampleSize, prior::Prior)
    obj = deepcopy(obj)
    obj.prior = prior
    return obj
end

function add!(
        JuMP_model_and_indicator_variables::Tuple,
        objective::MiniMaxSampleSize,
        problem::Problem;
        partial::Tuple{TI,TI} = (0, 0)
    ) where {TI<:Integer}

    m, ind = JuMP_model_and_indicator_variables
    prior  = objective.prior

    @variable(m, nmax)
    for n1_ in problem.n1values
        for x1_ in partial[1]:n1_
            @constraint(m, nmax >= sum(
                (n1_ + n2) * ind[(n1_, x1_, n2, c2)]
                    for (x1, n2, c2) in grid(problem, n1_)
                    if  (x1 == x1_)
                )
            )
        end
    end
    @objective(m, Min,
        (1 - objective.λ)*nmax +
        objective.λ*sum(
            (n1 + n2) * pmf(x1, n1, objective.prior; partial = partial) * ind[(n1, x1, n2, c2)]
                for (n1, x1, n2, c2) in grid(problem)
        )
    )
end


text:
using OmniSci, Test, Dates, Random, DataFrames, DecFP, GeoInterface, LibGEOS, Decimals

#defaults for OmniSci CPU Docker image
host="localhost"
user="admin"
passwd="HyperInteractive"
port=6274
dbname="omnisci"

#global connection to database, used for most tests
conn = connect(host, port , user, passwd, dbname)

@testset "connect/disconnect" begin
   #set up connection just to disconnect
   conn_d = connect(host, port , user, passwd, dbname)
   @test typeof(conn_d) == OmniSci.OmniSciConnection

   disc = disconnect(conn_d)
   @test typeof(disc) == Nothing
end

@testset "get_status" begin
   #get properties of server: version, rendering enabled, start time, etc.
   sstatus = get_status(conn)
   @test typeof(sstatus) == Vector{OmniSci.TServerStatus}
   @test sstatus[1].start_time <= Dates.datetime2unix(now(UTC)) #test sstatus has values...assumes server started before test happened
end

@testset "get_users" begin
   #get users for the database specified in connect()
   users = get_users(conn)
   @test typeof(users) == DataFrame

   users_nodf = get_users(conn, as_df = false)
   @test typeof(users_nodf) == Vector{String}
end

@testset "get_databases" begin
   #get databases user specified during connect() can see
   databases = get_databases(conn)
   @test typeof(databases) == DataFrame

   databases_nodf = get_databases(conn, as_df = false)
   @test typeof(databases_nodf) == Vector{OmniSci.TDBInfo}
end

@testset "get_tables_meta" begin
   #gets tables/views and their properties
   met = get_tables_meta(conn)
   @test typeof(met) == DataFrame

   met_nodf = get_tables_meta(conn, as_df = false)
   @test typeof(met_nodf) == Vector{OmniSci.TTableMeta}
end

@testset "sql_execute df" begin
   #execute query, return result via Thrift interface
   se = sql_execute(conn, "select * from omnisci_counties")
   @test typeof(se) == DataFrame
   @test size(se) == (3250, 6)
   @test eltype.(eachcol(se)) ==  [Union{Missing, String},
                          Union{Missing, String},
                          Union{Missing, String},
                          Union{Missing, String},
                          Union{Missing, String},
                          GeoInterface.MultiPolygon
                          ]
end

@testset "create/get_dashboards and grantees" begin
   #cdash represents id of dashboard created, use for later tests
   cdash = OmniSci.create_dashboard(conn, randstring(10), "state", "image", "metadata")
   @test cdash > 0

   getdbs = get_dashboards(conn)
   @test typeof(getdbs) == DataFrame

   getdbs_nodf = get_dashboards(conn, as_df = false)
   @test typeof(getdbs_nodf) == Vector{OmniSci.TDashboard}

   getdashgrant = get_dashboard_grantees(conn, cdash)
   @test typeof(getdashgrant) == Vector{OmniSci.TDashboardGrantees}
end

@testset "get_roles" begin
   gr = get_roles(conn)
   @test typeof(gr) == DataFrame

   gr_nodf = get_roles(conn, as_df = false)
   @test typeof(gr_nodf) == Vector{String}
end

@testset "get_table_details" begin
   tbl_details = get_table_details(conn, "omnisci_counties")
   @test typeof(tbl_details) == DataFrame
   @test size(tbl_details) == (6,21)
   @test names(tbl_details) == [:col_name, :col_type, :comp_param, :encoding, :is_array, :is_physical,
                               :is_reserved_keyword, :is_system, :nullable, :precision, :scale,
                               :size, :src_name, :fragment_size, :page_size, :max_rows, :view_sql,
                               :shard_count, :key_metainfo, :is_temporary, :partition_detail]

   tbl_details_nodf = get_table_details(conn, "omnisci_counties", as_df = false)
   @test typeof(tbl_details_nodf) == OmniSci.TTableDetails
end

@testset "get_all_roles_for_user" begin
   #drop for convenience so later test will pass, just in case tests run multiple times
   "testuser" in get_roles(conn, as_df=false) ? sql_execute(conn, "drop role testuser") : nothing
   "mapd2" in get_users(conn, as_df=false)  ? sql_execute(conn, "drop user mapd2") : nothing

   #get roles assigned to user
   sql_execute(conn, "create role testuser")
   sql_execute(conn, "create user mapd2 (password = 'OmniSciRocks!', is_super = 'true')")
   sql_execute(conn, "grant testuser to mapd2")
   roleuser = get_all_roles_for_user(conn, "mapd2")
   @test typeof(roleuser) == DataFrame

   roleuser_nodf = get_all_roles_for_user(conn, "mapd2", as_df = false)
   @test typeof(roleuser_nodf) == Vector{String}
end

@testset "create and load: atomic types" begin

   tinyintcol = Union{Int8,Missing}[missing,3,2,1]
   smallintcol = Union{Int16,Missing}[4,missing,2,1]
   intcol = Union{Int32,Missing}[4,3,missing,1]
   bigintcol = Union{Int64,Missing}[4,3,2,missing]
   floatcol = Union{Float32,Missing}[missing, 4.1, 2.69, 3.8]
   doublecol = [3//2, missing, 9//72, 90/112]
   boolcol = [missing, false, true, true]
   textcol = ["hello", "world", "omnisci", "gpu"]
   datecol = [Date(2013,7,1), missing, Date(2013,7,1), Date(2013,7,1)]
   timecol = [Time(4), Time(5), missing, Time(7)]
   tscol = [DateTime(2013,7,1,12,30,59), DateTime(2013,7,1,12,30,59), DateTime(2013,7,1,12,30,59), missing]

   df = DataFrame(x1 = tinyintcol,
                  x2 = smallintcol,
                  x3 = intcol,
                  x4 = bigintcol,
                  x5 = floatcol,
                  x6 = doublecol,
                  x7 = boolcol,
                  x8 = textcol,
                  x9 = datecol,
                  x10 = timecol,
                  x11 = tscol)

   #drop table if it exists
   tables = get_tables_meta(conn)
   "test_int_float" in tables[!, :table_name] ? sql_execute(conn, "drop table test_int_float") : nothing

   @test create_table(conn, "test_int_float", df) == nothing

   #load data rowwise from dataframe
   @test load_table(conn, "test_int_float", df) == nothing

   #load data colwise from dataframe
   @test load_table_binary_columnar(conn, "test_int_float", df) == nothing

   #test roundtrip of data
   #NOTE: These tests implicitly rely on having the missing values for numeric
   #Column types almost always returned from OmniSci with a Missing type component
   #So isequal will likely fail if missings removed, as column types won't match
   tbldb = sql_execute(conn, "select * from test_int_float")
   @test size(tbldb) == (8,11)
   @test isequal(df, tbldb[1:4, :])
   @test isequal(vcat(df, df), tbldb)
end

@testset "create and load: geospatial" begin

   pointcol = ["POINT (30 10)",
               "POINT (-30.18764587 12.2)",
               "POINT (30 -10.437878634)",
               "POINT (-78 -25)"]

   linecol = ["LINESTRING (30 10, 10 30, 40 40)",
              "LINESTRING (30 10, 10 30, 40 40)",
              "LINESTRING (30 10, 10 30, 40 40)",
              "LINESTRING (30 10, 10 30, 40 40)"]

   polycol = ["POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))",
              "POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))",
              "POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))",
              "POLYGON ((30 10, 40 40, 20 40, 10 20, 30 10))"]

   mpolycol = ["MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)),((15 5, 40 10, 10 20, 5 10, 15 5)))",
               "MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)),((15 5, 40 10, 10 20, 5 10, 15 5)))",
               "MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)),((15 5, 40 10, 10 20, 5 10, 15 5)))",
               "MULTIPOLYGON (((30 20, 45 40, 10 40, 30 20)),((15 5, 40 10, 10 20, 5 10, 15 5)))"]

   pointcol_native = GeoInterface.Point.(readgeom.(pointcol))
   linecol_native = GeoInterface.LineString.(readgeom.(linecol))
   polycol_native = GeoInterface.Polygon.(readgeom.(polycol))
   mpolycol_native = GeoInterface.MultiPolygon.(readgeom.(mpolycol))

   pointcol_geos = readgeom.(pointcol)
   linecol_geos = readgeom.(linecol)
   polycol_geos = readgeom.(polycol)
   mpolycol_geos = readgeom.(mpolycol)

   df = DataFrame(x3 = pointcol_native,
                  x4 = linecol_native,
                  x5 = polycol_native,
                  x6 = mpolycol_native
                  )

   df2 = DataFrame(x3 = pointcol,
                  x4 = linecol,
                  x5 = polycol,
                  x6 = mpolycol
                  )

   df3 = DataFrame(x3 = pointcol_geos,
                  x4 = linecol_geos,
                  x5 = polycol_geos,
                  x6 = mpolycol_geos
                  )

   #drop table if it exists
   tables = get_tables_meta(conn)
   "test_geo_native" in tables[!, :table_name] ? sql_execute(conn, "drop table test_geo_native") : nothing

   @test create_table(conn, "test_geo_native", df) == nothing

   #load data rowwise from dataframe: GeoInterface
   @test load_table(conn, "test_geo_native", df) == nothing

   #load data rowwise from dataframe: WKT
   @test load_table(conn, "test_geo_native", df2) == nothing

   #load data rowwise from dataframe: LibGEOS
   @test load_table(conn, "test_geo_native", df3) == nothing

   #load data colwise from dataframe: GeoInterface
   @test load_table_binary_columnar(conn, "test_geo_native", df) == nothing

   #load data colwise from dataframe: WKT
   @test load_table_binary_columnar(conn, "test_geo_native", df2) == nothing

   #load data colwise from dataframe: LibGEOS
   @test load_table_binary_columnar(conn, "test_geo_native", df3) == nothing

   tbldb = sql_execute(conn, "select * from test_geo_native")
   @test size(tbldb) == (24,4)

   # test roundtrip of data, isequal on dataframe doesn't seem to work
   # WKT and GEOS test a bit hacky, since comparison is converted GeoInterface from WKT against what OmniSci returns
   # this test checks that tbldb matches the original df part on a coordinate level by cellwise indexing
   # each load table uploads 4 rows, so the comparison is the original df against the rows uploaded and indexed in tbldb
   for i in 1:4
      for j in 1:4
         @test isequal(df[i,j].coordinates, tbldb[i,j].coordinates)
         @test isequal(df[i,j].coordinates, tbldb[i + 4,j].coordinates)
         @test isequal(df[i,j].coordinates, tbldb[i + 8,j].coordinates)
         @test isequal(df[i,j].coordinates, tbldb[i + 12,j].coordinates)
         @test isequal(df[i,j].coordinates, tbldb[i + 16,j].coordinates)
         @test isequal(df[i,j].coordinates, tbldb[i + 20,j].coordinates)
      end
   end

end

@testset "create and load: arrays" begin
   #Missing in array not supported in OmniSci 4.4
   #TODO: when missing in array supported, add test
   tinyintcol = Int8[4,3,2,1]
   smallintcol = Int16[4,3,2,1]
   intcol = Int32[4,3,2,1]
   bigintcol = Int64[4,3,2,1]
   floatcol = Float32[0.0, 4.1, 2.69, 3.8]
   doublecol = [3//2, 1/2, 9//72, 90/112]
   boolcol = [true, false, true, true]
   textcol = ["hello", "world", "omnisci", "gpu"]
   datecol = [Date(2013,7,1), Date(2013,7,1), Date(2013,7,1), Date(2013,7,1)]
   timecol = [Time(4), Time(5), Time(6), Time(7)]
   tscol = [DateTime(2013,7,1,12,30,59), DateTime(2013,7,1,12,30,59), DateTime(2013,7,1,12,30,59), DateTime(2013,7,1,12,30,59)]

   tinyarrarr = [tinyintcol, tinyintcol, tinyintcol, tinyintcol]
   smallarrarr = [smallintcol, smallintcol, smallintcol, smallintcol]
   intarrarr = [intcol, intcol, intcol, intcol]
   bigintarrarr = [bigintcol, bigintcol, bigintcol, bigintcol]
   textarrarr = [textcol, textcol, textcol, textcol]
   floatarrarr = [floatcol, floatcol, floatcol, floatcol]
   dblarrarr = [doublecol, doublecol, doublecol, doublecol]
   boolarrarr = [boolcol, boolcol, boolcol, boolcol]
   datearrarr = [datecol, datecol, datecol, datecol]
   timearrarr = [timecol, timecol, timecol, timecol]
   tsarrarr = [tscol, tscol, tscol, tscol]

   df = DataFrame(c1 = tinyarrarr,
                  c2 = smallarrarr,
                  c3 = intarrarr,
                  c4 = bigintarrarr,
                  c5 = textarrarr,
                  c6 = floatarrarr,
                  c7 = dblarrarr,
                  c8 = boolarrarr,
                  c9 = datearrarr,
                  c10 = timearrarr,
                  c11 = tsarrarr)

   #drop table if it exists
   tables = get_tables_meta(conn)
   "test_array" in tables[!, :table_name] ? sql_execute(conn, "drop table test_array") : nothing

   @test create_table(conn, "test_array", df) == nothing

   #load data rowwise from dataframe
   @test load_table(conn, "test_array", df) == nothing

   #load arrays columnar
   @test load_table_binary_columnar(conn, "test_array", df) == nothing

   #Tests for arrays indirect due to issue querying. Test based on known properties
   #Write proper tests once https://github.com/omnisci/OmniSci.jl/issues/53 solved
   @test sql_execute(conn, "select count(*) as records from test_array")[!, :records][1] == 8

end

@testset "sql_execute with decimal column" begin

   #This tests that sql_execute returns columns defined in a table
   #load_table and actual decimal types to come later

   #drop table if it exists
   tables = get_tables_meta(conn)
   "test_decimals" in tables[!, :table_name] ? sql_execute(conn, "drop table test_decimals") : nothing

   #hand-create table here since test is that anything decimal-like will be loaded correctly
   #if table already exists
   sql = """create table test_decimals(
            x1 DECIMAL(7,3),
            x2 DECIMAL(7,3),
            x3 DECIMAL(7,3),
            x4 DECIMAL(7,3)
            )
         """

   @test sql_execute(conn, sql) == nothing

   str = String["1234.567", "345.999", "16.876", "9.987"]
   flo = Float32[1234.567, 345.999, 16.876, 9.987]
   doub = Float64[1234.567, 345.999, 16.876, 9.987]
   rat = Rational[1234567//1000, 345999//1000, 16876//1000, 9987//1000]

   df = DataFrame(x1 = str, x2 = flo, x3 = doub, x4 = rat)

   @test load_table(conn, "test_decimals", df) == nothing

   tbldb = sql_execute(conn, "select * from test_decimals")
   @test size(tbldb) == (4,4)

   #Validate roundtrip within reason
   #Because inputs aren't Dec64 (the default return type), test approx equal
   @test isequal(Float32.(tbldb[!, :x1]), [parse(Float32, x) for x in str])
   @test isequal(Float64.(tbldb[!, :x1]), [parse(Float64, x) for x in str])
   @test isequal(Float32.(tbldb[!, :x2]), flo)
   @test isequal(Float64.(tbldb[!, :x3]), doub)
   @test isequal(Float32.(tbldb[!, :x4]), Float32.(rat))
   @test isequal(Float64.(tbldb[!, :x4]), Float64.(rat))

end

@testset "create/load_table with decimal column" begin

   #drop table if it exists
   tables = get_tables_meta(conn)
   "test_decimals_createload" in tables[!, :table_name] ? sql_execute(conn, "drop table test_decimals_createload") : nothing

   d32 = DecFP.Dec32.(["123.456", "1333.758", "23.3", "234.22"])
   d64 = DecFP.Dec64.(["123.456", "1333.758", "23.3", "234.22"])
   d128 = DecFP.Dec128.(["123.456", "1333.758", "23.3", "234.22"])
   dd = Decimals.decimal.(["123.456", "1333.758", "23.3", "234.22"])

   df = DataFrame(x1 = d32,
                  x2 = d64,
                  x3 = d128,
                  x4 = dd)

   #Test that assertion thrown if decimal columns without specifying precision
   @test_throws AssertionError create_table(conn, "test_decimals_createload", df)

   @test create_table(conn, "test_decimals_createload", df, precision=(7,3)) == nothing

   #load data rowwise from dataframe
   @test load_table(conn, "test_decimals_createload", df) == nothing

   tbldb = sql_execute(conn, "select * from test_decimals_createload")

   @test size(tbldb) == (4,4)

   #Since data are returned as Dec64,
   @test isequal(tbldb[!, :x1], df[!, :x2])

end

@testset "get_hardware_info" begin
   #TODO: create a show method and/or return as dataframe
   hware = get_hardware_info(conn)
   @test typeof(hware) == OmniSci.TClusterHardwareInfo
end

@testset "sql_execute_df" begin
   #TODO: Figure out IPC
   cpu_arrow = sql_execute_df(conn,  "select id from omnisci_counties limit 100", 0, 0)
   @test typeof(cpu_arrow) == OmniSci.TDataFrame
end

######################################## not exported (essentially, OmniSci internal)

@testset "clear_cpu_memory" begin
   clear_cpu = OmniSci.clear_cpu_memory(conn)
   @test typeof(clear_cpu) == Nothing
end

@testset "set_execution_mode" begin
   execmode = OmniSci.set_execution_mode(conn, TExecuteMode.CPU)
   @test typeof(execmode) == Nothing
end

@testset "get_memory" begin
   mem = OmniSci.get_memory(conn, "cpu")
   @test typeof(mem) == Vector{OmniSci.TNodeMemoryInfo}
end


text:
using StaticArrays, LinearAlgebra, KiteUtils
using KiteModels, KitePodModels

if ! @isdefined kcu
    const kcu = KCU()
end
if ! @isdefined kps
    const kps = KPS4(kcu)
end

const dt = 0.05

clear!(kps)
KiteModels.set_depower_steering!(kps, kps.set.depower_offset/100.0, 0.0)
kps.stiffness_factor = 0.04

@time KiteModels.find_steady_state!(kps, prn=true)

println("\nlift, drag    [N]  : $(KiteModels.lift_drag(kps))")
# println("\nSpring forces:")
# spring_forces(kps)


text:
module Neovim

using Compat
using MsgPack
import MsgPack: pack, unpack
import Sockets: connect

export NvimClient, nvim_connect, nvim_env, nvim_spawn, nvim_child, start_host
export reply_result, reply_error

# types that have api methods (nvim itself + api defined types)
abstract type NvimObject end

include("client.jl")

# Only export API when not being checked by JuliaRegistries Automerge CI since
# `nvim` binary isn't available there,
# see https://github.com/JuliaRegistries/RegistryCI.jl/issues/283#issuecomment-715362834
if !(get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", "") == "true")
    export Buffer, Tabpage, Window
    include("api_gen.jl")
    include("interface.jl")
    include("plugin_host.jl")
end

# too inconvenient api to supply handler here?
function nvim_connect(path::String, args...)
    s = connect(path)
    NvimClient(s, s, args...)
end

nvim_env(args...) = nvim_connect(ENV["NVIM_LISTEN_ADDRESS"], args...)

function nvim_spawn(args...; cmd=`nvim --embed`)
    # TODO(smolck): Maybe just use the Process handle instead of p.in and p.out
    p = open(cmd, "r+")
    (NvimClient(p.in, p.out, args...), p)
end

function nvim_child(args...)
    # make stdio private. Reversed since from nvim's perspective
    input, output = stdout, stdin
    if haskey(ENV, "NEOVIM_JL_DEBUG")
        debug = open(ENV["NEOVIM_JL_DEBUG"], "w")
        redirect_stdout(debug)
        redirect_stderr(debug)
    else
        redirect_stdout()
        redirect_stderr()
    end
    redirect_stdin()

    NvimClient(input, output, args...)
end

end


text:
include("initialization.jl")
include("crossover.jl")
include("mutation.jl")

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using Test
using IntervalSets
using LinearAlgebra
using Literate

#include("test_finufft.jl")
include("test_splinepp.jl")
include("test_splinenn.jl")
include("test_geometry.jl")
include("linear_pendulum.jl")
include("test_meshes.jl")
include("test_cubic_splines.jl")
include("test_poisson_1d.jl")
include("test_poisson_2d.jl")
include("test_bsl.jl")
include("test_rotation_2d.jl")

examples = ["vlasov-ampere",
            "vlasov-poisson",
            "bump_on_tail",
            "rotation2d_bsl",
            "vlasov-hmf"]

output = joinpath(@__DIR__, "examples")

for example in examples
    filename = joinpath(@__DIR__, "..", "examples", string(example,".jl"))
    Literate.script(filename, output)
    include(joinpath(@__DIR__, "examples", string(example,".jl")))
end


text:
using ProbNumDiffEq
using Test
using Plots
using LinearAlgebra
using OrdinaryDiffEq
using Statistics
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems;
importodeproblems();
import DiffEqProblemLibrary.ODEProblemLibrary:
    prob_ode_linear, prob_ode_2Dlinear, prob_ode_lotkavoltera, prob_ode_fitzhughnagumo

@testset "Solution" begin
    prob = ProbNumDiffEq.remake_prob_with_jac(prob_ode_lotkavoltera)
    sol = solve(prob, EK1())

    @test length(sol) > 2
    @test length(sol.t) == length(sol.u)
    @test length(prob.u0) == length(sol.u[end])

    # Destats
    @testset "DEStats" begin
        @test length(sol.t) == sol.destats.naccept + 1
        @test sol.destats.naccept <= sol.destats.nf
    end

    @testset "Hit the provided tspan" begin
        @test sol.t[1] == prob.tspan[1]
        @test sol.t[end] == prob.tspan[2]
    end

    @testset "Prob and non-prob u" begin
        @test sol.u == sol.pu.μ
    end

    @testset "Call on known t" begin
        @test sol(sol.t).u == sol.pu
    end

    @testset "Correct initial values" begin
        @test sol.pu[1].μ == prob.u0
        @test iszero(sol.pu[1].Σ)
    end

    # Interpolation
    @testset "Dense Solution" begin
        t0 = prob.tspan[1]
        t1, t2 = t0 + 1e-2, t0 + 2e-2

        u0, u1, u2 = sol(t0), sol(t1), sol(t2)
        @test norm.(u0.μ - u1.μ) < norm.(u0.μ - u2.μ)

        @test all(diag(u1.Σ) .< diag(u2.Σ))

        @test sol.(t0:1e-3:t1) isa Array{Gaussian{T,S}} where {T,S}
        @test sol(t0:1e-3:t1).u isa StructArray{Gaussian{T,S}} where {T,S}
    end

    # Sampling
    @testset "Solution Sampling" begin
        n_samples = 10

        samples = ProbNumDiffEq.sample(sol, n_samples)
        dsamples, dts = ProbNumDiffEq.dense_sample(sol, n_samples)

        @test samples isa Array

        m, n, o = size(samples)
        @test m == length(sol)
        @test n == length(sol.u[1])
        @test o == n_samples

        u = ProbNumDiffEq.stack(sol.u)
        stds = sqrt.(ProbNumDiffEq.stack(diag.(sol.pu.Σ)))
        outlier_count = sum(abs.(u .- samples) .> 3stds)
        @test outlier_count < 0.05 * m * n * o

        # Dense sampling
        dense_samples, dense_times = ProbNumDiffEq.dense_sample(sol, n_samples)
        m, n, o = size(dense_samples)
        @test m == length(dense_times)
        @test n == length(sol.u[1])
        @test o == n_samples
    end

    @testset "Sampling states from the solution" begin
        n_samples = 10

        samples = ProbNumDiffEq.sample_states(sol, n_samples)
        dsamples, dts = ProbNumDiffEq.dense_sample_states(sol, n_samples)

        @test samples isa Array

        m, n, o = size(samples)
        @test m == length(sol)
        @test n == length(sol.u[1]) * (sol.interp.q + 1)
        @test o == n_samples

        x = ProbNumDiffEq.stack(mean(sol.x_smooth))
        stds = sqrt.(ProbNumDiffEq.stack(diag.(sol.x_smooth.Σ)))
        outlier_count = sum(abs.(x .- samples) .> 3stds)
        @test outlier_count < 0.05 * m * n * o

        # Dense sampling
        dense_samples, dense_times = ProbNumDiffEq.dense_sample_states(sol, n_samples)
        m, n, o = size(dense_samples)
        @test m == length(dense_times)
        @test n == length(sol.u[1]) * (sol.interp.q + 1)
        @test o == n_samples
    end

    @testset "Plotting" begin
        @test plot(sol) isa AbstractPlot
        @test plot(sol, denseplot=false) isa AbstractPlot
        @test plot(sol, vars=(1, 2)) isa AbstractPlot
    end

    @testset "Mean Solution" begin
        @test mean(sol) isa DiffEqBase.AbstractODESolution
        @test plot(mean(sol)) isa AbstractPlot
    end
end


text:
using Documenter
using OneDmin

makedocs(
    sitename = "OneDmin.jl",
    format = Documenter.HTML(assets = ["assets/style.css"], prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [OneDmin],
    pages = ["Home" => "index.md"]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(repo = "github.com/vepiteski/OneDmin.jl")
#https://juliadocs.github.io/Documenter.jl/stable/man/hosting/ ?


text:
# infrastructur.jl

include("../api/api.jl")

# process(training::Training) = begin
process() = begin
    # Dummy data for testing purposes
    training = create_training("LS", DateTime(2019, 8, 30), 2, "Learn Smiling", 1000.0)

    org1 = create_org(
        "Scrooge Investment Bank",
        "1180 Seven Seas Dr",
        "FL 32830",
        "Lake Buena Vista",
        "USA")
    org2 = create_org(
        "Duck City Chronicals",
        "1185 Seven Seas Dr",
        "FL 32830",
        "Lake Buena Vista",
        "USA")
    org3 = create_org(
        "Donalds Hardware Store",
        "1190 Seven Seas Dr",
        "FL 32830",
        "Lake Buena Vista",
        "USA")

    order1 = create_order(
        org1,
        training,
        "PO-456",
        "Scrooge McDuck",
        "scrooge@duckcity.com",
        ["Scrooge McDuck"])
    order2 = create_order(
        org2,
        training,
        "DD-001",
        "Mickey Mouse",
        "mickey@duckcity.com",
        ["Mini Mouse", "Goofy"])
    order3 = create_order(
        org3,
        training,
        "",
        "Donald Duck",
        "donald@duckcity.com",
        ["Daisy Duck"])

    orders = [order1, order2, order3]

    return orders
end


text:
function simple_toy_MetNet()
    net = Dict()
    net[:S] =
    # rxns: gt    ferm  resp  ldh   lt   biom    atpm  # mets
        [   1.0  -1.0   0.0   0.0   0.0   0.0    0.0;  #  G
            0.0   2.0  18.0   0.0   0.0  -55.0  -5.0;  #  E
            0.0   2.0  -1.0  -1.0   0.0   0.0    0.0;  #  P
            0.0   0.0   0.0   1.0   1.0   0.0    0.0;  #  L
        ]

    net[:mets] = ["G", "E", "P", "L"]
    net[:b] =    [0.0, 0.0, 0.0, 0.0] # demand

    net[:metNames] = ["Glucose", "Energy", "Intermediate Product" , "Lactate"];

    net[:rxns] = ["gt"  ,"ferm" ,"resp" , "ldh" ,  "lt" , "biom", "atpm"];
    net[:lb] =   [0.0   , 0.0   , 0.0   ,  0.0  , -100.0,   0.0,     0.5];
    net[:ub] =   [10.0 , 100.0 , 100.0 , 100.0 ,    0.0, 100.0,    100.0];
    net[:c] =    [ 0.0 ,   0.0 ,   0.0 ,   0.0 ,    0.0,   0.0,      0.0];
    net[:rxnNames] = ["Glucose transport", "Fermentation", "Respiration",
        "Lactate DH", "Lactate transport", "Biomass production rate", "atp demand"];

    return MetNet(;net...)
end
const TOY_MODEL_BIOMASS_IDER = "biom"
const TOY_MODEL_COST_IDER = "tot_cost"

function  toy_model(;resp_cost::Real = -0.1, E_demand::Real = 5.0)
    @assert resp_cost <= 0.0

    model = simple_toy_MetNet()

    # Atpm
    lb!(model, "atpm", E_demand)

    # cost
    # Add a new metabolite simulating the cost penalazing
    # reaction fluxes.
    # A new balance equations is then added:
    #     Σ(rᵢ*costᵢ) + tot_cost = 0
    # Because the cost coefficients (costᵢ) < 0, the system must allocate
    # the fluxes (rᵢ) so that Σ(rᵢ*costᵢ) = tot_cost, and tot_cost
    # are usually bounded [0.0, 1.0]

    # We will add 1 met and 1 rxn
    M, N = size(model)
    model = expanded_model(model, M + 1, N + 1)
    # Add cost met
    set_met!(model, M + 1, Met("cost"; rxns = ["resp"], S = [resp_cost]))
    set_rxn!(model, N + 1, Rxn("tot_cost"; mets = ["cost"], S = [1.0], lb = 0.0, ub = 1.0))
    return model
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c7952ee6-45b5-11eb-1158-5bb2ff274ce9
begin
	using DrWatson

	function ingredients(path::String)
		# this is from the Julia source code (evalfile in base/loading.jl)
		# but with the modification that it returns the module instead of the last object
		name = Symbol(basename(path))
		m = Module(name)
		Core.eval(m,
			Expr(:toplevel,
				 :(eval(x) = $(Expr(:core, :eval))($name, x)),
				 :(include(x) = $(Expr(:top, :include))($name, x)),
				 :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
				 :(include($path))))
		m
	end

	mtg_dir = projectdir() * "/games/MtG"
end;

# ╔═╡ 0c145632-3927-11eb-19b9-877e05c1bcdc
begin
	@quickactivate

	using JSON
	using HTTP
	using Plots
	using Images
	using PlutoUI
	using Serialization
	using ImageTransformations: imresize

	plotly()

	GS = ingredients("$mtg_dir/MtG.jl/notebooks/game_settings.jl")

	md"""
	*game settings loaded!*

	## MtG EDH Kraum and Tymna by ???
	"""
end

# ╔═╡ 621b08a4-384e-11eb-0109-61e9b9ecf125
deck = Dict{Symbol,Any}(
    :commanders => [
        "Kraum, Ludevic's Opus",
        "Tymna the Weaver"
        ],
    :cards => [
        "Ajani Vengeant",
        "Arid Mesa",
        "Aven Mindcensor",
        "Badlands",
        "Baleful Strix",
        "Blood Crypt",
        "Bloodstained Mire",
        "Brain Maggot",
        "Brainstorm",
        "Brazen Borrower",
        "Breya, Etherium Shaper",
        "Celestial Colonnade",
        "Chain Lightning",
        "City of Brass",
        "Command Tower",
        "Counterspell",
        "Creeping Tar Pit",
        "Dark Confidant",
        "Dardslick Shores",
        "Daze",
        "Deadly Rollick",
        "Demonic Tutor",
        "Duress",
        "Evasive Action",
        "Fatal Push",
        "Fiery Islet",
        "Fire // Ice",
        "Flawless Maneuver",
        "Flooded Strand",
        "Force of Negation",
        "Force of Will",
        "Force Spike",
        "Chop Down // Giant Killer",
        "Gitaxian Probe",
        "Giver of Runes",
        "Godless Shrine",
        "Grim Lavamancer",
        "Hallowed Fountain",
        "Inquisition of Kozilek",
        "Island",
        "Jace, the Mind Sculptor",
        "Kitesail Freebooter",
        "Kraum, Ludevic's Opus",
        "Lightning Bolt",
        "Lightning Helix",
        "Mana Confluence",
        "Mana Leak",
        "Mana Tithe",
        "Marsh Flats",
        "Mausoleum_Wanderer",
        "Memory Lapse",
        "Mental Misstep",
        "Mesmeric Fiend",
        "Miscalculation",
        "Misty Rainforest",
        "Mother of Runes",
        "Plains",
        "Plateau",
        "Polluted Delta",
        "Ponder",
        "Preordain",
        "Prismatic Vista",
        "Reflecting Pool",
        "Remorseful Cleric",
        "Sacred Foundry",
        "Scalding Tarn",
        "Scrubland",
        "Seachrome Coast",
        "Selfless Spirit",
        "Silent Clearing",
        "Smuggler's Copter",
        "Snapcaster Mage",
        "Spell Pierce",
        "Spell Queller",
        "Spell Snare",
        "Spirebluff Canal",
        "Steam Vents",
        "Stoneorge Mystic",
        "Sunbaked Canyon",
        "Swamp",
        "Sword of Feast and Famine",
        "Swords to Plowshares",
        "Tainted Pact",
        "Thoughtseize",
        "Tidehollow Sculler",
        "Tollegeist",
        "Tribal Flames",
        "True-Name Nemesis",
        "Tundra",
        "Umezawa's Jitte",
        "Underground Sea",
        "Vendillion Clique",
        "Verdant Catacombs",
        "Vindicate",
        "Volcanic Island",
        "Watery Grave",
        "Windswept Heath",
        "Wooded Foothills",
        ]
    )

# ╔═╡ fb61d01c-458d-11eb-2c2a-f711dc7ab7f4
md"""
Found $(length(deck[:cards]) + length(deck[:commanders])) cards!
"""

# ╔═╡ c61fa79e-4583-11eb-3b71-2d334aca843d
begin
	mtg_cards = JSON.parsefile("$mtg_dir/MtG.jl/json/oracle-cards-20201224220555.json")

	md"""
	###### MtG database loaded. Found $(length(mtg_cards)) unique cards (by name).
	Most recent data available here: https://scryfall.com/docs/api/bulk-data
	"""
end

# ╔═╡ 7420cf10-45cc-11eb-2780-4f320bd8a2cf
begin  # note: this func only downloads the first card with a matching name and then moves to the next.
	deck_cards = []

	for n in vcat(deck[:cards], deck[:commanders])

		for c in mtg_cards

			if n == c["name"]
				push!(deck_cards, c)
				break
			end
		end
	end

	deck_cards
end

# ╔═╡ 2ec9a2fe-4660-11eb-2b73-7127019757b5
deck[:cardnames] = vcat(deck[:cards], deck[:commanders])

# ╔═╡ 9f64efdc-4660-11eb-2333-45a14af374c6
missing_cards = [ c for c in deck[:cardnames] if !(c in [ n for n in deck[:cardnames] ]) ]

# ╔═╡ 489f3da8-4681-11eb-26af-f75d8ecc552e
@bind i Slider(1:length(deck_cards), show_value=true)

# ╔═╡ 4666005c-4d62-11eb-1215-e7e010c3125c
md"""
#### Look good? These images will be displayed in-game.
TODO: write support for previews of double-sided cards
"""

# ╔═╡ fb0a03ea-4ed5-11eb-1280-e70793eac293
@bind ratio Slider(0.1:0.05:1, show_value=true)

# ╔═╡ 5f7ebd78-3db7-11eb-0690-1b8ee4ebe7db
begin
	CARD_BACK_PATH = "$mtg_dir/MtG.jl/ui/cards/card_back.png"
	im = imresize(load(CARD_BACK_PATH), ratio=ratio)

	md"""
	`CARD_BACK_IMG`:

	$(CARD_BACK_IMG = im)
	"""
end

# ╔═╡ 457f0ee2-4d6f-11eb-310f-cf4784a06469
html"""

"""

# ╔═╡ 9ded258e-468d-11eb-3428-d915bfe9e13e
function get_card_img(img_uri::String)
	img_resp = HTTP.get(img_uri)
	card_img = img_resp.body |> IOBuffer |> load
end

# ╔═╡ e67591a6-4ed5-11eb-215e-e3497b7ef30e
deck_card_img = imresize(get_card_img(deck_cards[i]["image_uris"]["border_crop"]), ratio=ratio
)

# ╔═╡ a26b32f4-468d-11eb-1ddd-7958d61b4ac6
function search_cards_by_keyword(q::String, mtg_cards::Array)
	[ n for n in [ c["name"] for c in mtg_cards ] if occursin(q, n) ]
end

# ╔═╡ cb8e3b6c-4661-11eb-152e-3f2ce56cdd9d
search_cards_by_keyword("Lotus", mtg_cards)

# ╔═╡ d491b248-4da5-11eb-115b-193dfb0e1979
begin
	deck[:CARD_BACK_PATH] = CARD_BACK_PATH
	deck[:CARD_RATIO] = ratio
end

# ╔═╡ Cell order:
# ╟─c7952ee6-45b5-11eb-1158-5bb2ff274ce9
# ╟─0c145632-3927-11eb-19b9-877e05c1bcdc
# ╟─621b08a4-384e-11eb-0109-61e9b9ecf125
# ╟─fb61d01c-458d-11eb-2c2a-f711dc7ab7f4
# ╟─5f7ebd78-3db7-11eb-0690-1b8ee4ebe7db
# ╟─c61fa79e-4583-11eb-3b71-2d334aca843d
# ╟─7420cf10-45cc-11eb-2780-4f320bd8a2cf
# ╟─2ec9a2fe-4660-11eb-2b73-7127019757b5
# ╟─9f64efdc-4660-11eb-2333-45a14af374c6
# ╠═cb8e3b6c-4661-11eb-152e-3f2ce56cdd9d
# ╟─e67591a6-4ed5-11eb-215e-e3497b7ef30e
# ╟─489f3da8-4681-11eb-26af-f75d8ecc552e
# ╟─4666005c-4d62-11eb-1215-e7e010c3125c
# ╠═fb0a03ea-4ed5-11eb-1280-e70793eac293
# ╟─457f0ee2-4d6f-11eb-310f-cf4784a06469
# ╟─9ded258e-468d-11eb-3428-d915bfe9e13e
# ╟─a26b32f4-468d-11eb-1ddd-7958d61b4ac6
# ╠═d491b248-4da5-11eb-115b-193dfb0e1979

max_stars_repo_path: EDH/decks/Kraum and Tymna/Kraum and Tymna-backup.jl


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "quadruped", "visuals.jl"))
vis = Visualizer()
open(vis)
const ContactImplicitMPC = Main

s = get_simulation("quadruped", "flat_2D_lc", "flat")
model = s.model
env = s.env

function run_policy(s::Simulation; H_sim::Int = 2000, verbose = false,
		qinit::AbstractVector = zeros(s.model.dim.q))
	model = s.model
	env = s.env

	ref_traj = deepcopy(ContactImplicitMPC.get_trajectory(model, env,
	    joinpath(module_dir(), "src/dynamics/quadruped/gaits/gait2.jld2"),
	    load_type = :split_traj_alt))

	# time
	H = ref_traj.H
	h = ref_traj.h
	N_sample = 5
	H_mpc = 10
	h_sim = h / N_sample
	H_sim = H_sim

	# barrier parameter
	κ_mpc = 1.0e-4

	obj = TrackingObjective(model, env, H_mpc,
	    q = [Diagonal(1e-2 * [1.0; 0.02; 0.25; 0.75 * ones(model.dim.q-3)]) for t = 1:H_mpc],
	    u = [Diagonal(3e-2 * ones(model.dim.u)) for t = 1:H_mpc],
		γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
	    b = [Diagonal(1.0e-100 * ones(model.dim.c * friction_dim(env))) for t = 1:H_mpc])

	p = linearized_mpc_policy(ref_traj, s, obj,
	    H_mpc = H_mpc,
	    N_sample = N_sample,
	    κ_mpc = κ_mpc,
		mode = :configuration,
		ip_type = :mehrotra,
	    n_opts = NewtonOptions(
			r_tol = 3e-4,
			max_iter = 5,
			max_time = ref_traj.h, # HARD REAL TIME
			),
	    mpc_opts = LinearizedMPCOptions(
	        ),
		ip_opts = InteriorPointOptions(
			max_iter = 100,
			verbose = false,
			r_tol = 1.0e-4,
			κ_tol = 1.0e-4,
			diff_sol = true,
			# κ_reg = 1e-3,
			# γ_reg = 1e-1,
			solver = :empty_solver,
			),
	    )

	q1_ref = copy(qinit)
	q0_ref = copy(qinit)
	q1_sim = SVector{model.dim.q}(q1_ref)
	q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
	@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

	sim = ContactImplicitMPC.simulator(s, q0_sim, q1_sim, h_sim, H_sim,
	    p = p,
	    ip_opts = ContactImplicitMPC.InteriorPointOptions(
	        r_tol = 1.0e-8,
	        κ_init = 1.0e-6,
	        κ_tol = 2.0e-6,
	        diff_sol = true),
	    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

	status = ContactImplicitMPC.simulate!(sim, verbose = false)
	return deepcopy(sim.traj)
end

function collect_runs(s::Simulation; n::Int = 1, H_sim::Int = 2000, verbose::Bool = false)
	nq = s.model.dim.q
	trajs = []
	Random.seed!(100)
	conf_max = [0.05, 0.8, 0.8, 0.8, +0.2,  0.10]
	conf_min = [0.00, 0.6, 0.6, 0.6, -0.2, -0.30]
	conf_Δ = conf_max .- conf_min
	for i = 1:n
		verbose && println("sample = $i/$n")
		conf = conf_min .+ conf_Δ .* rand(length(conf_min))
		conf[end] = max(conf[end], 0.0)
		qinit = initial_configuration(s.model, conf...)
		traj = run_policy(s, H_sim = H_sim, qinit = qinit, verbose = verbose)
		push!(trajs, traj)
	end
	return trajs
end

function initial_configuration(model::Quadruped, θ0, θ1, θ2, θ3, x, Δz)
    q = zeros(model.dim.q)
	q[1] = x
    q[3] = pi / 2.0
    q[4] = -θ1
    q[5] = θ2

    q[8] = -θ1
    q[9] = θ2

    q[2] = model.l_thigh1 * cos(q[4]) + model.l_calf1 * cos(q[5])

    q[10] = -θ3
    q[11] = acos((q[2] - model.l_thigh2 * cos(q[10])) / model.l_calf2)

    q[6] = -θ3
    q[7] = acos((q[2] - model.l_thigh2 * cos(q[6])) / model.l_calf2)

	q[2] += Δz
	q[3] += θ0
    return q
end

function visualize_runs!(vis::Visualizer, model::ContactModel, trajs::AbstractVector;
		sample=max(1, Int(floor(trajs[1].H / 100))), h=trajs[1].h*sample,  α=1.0,
		anim::MeshCat.Animation=MeshCat.Animation(Int(floor(1/h))),
		name::Symbol=model_name(model))

	for (i, traj) in enumerate(trajs)
		name_i = Symbol(string(name) * string(i))
		anim = visualize_meshrobot!(vis, model, traj, sample = sample, h = h, α = α, anim = anim, name = name_i)
	end
	return anim
end

# monte carlo
H_sim = 1000
trajs = collect_runs(s, n = 100, H_sim = H_sim, verbose = true)
anim = visualize_runs!(vis, s.model, trajs, α=0.3, sample=5)
rep_ref_traj = repeat_ref_traj(ref_traj, Int(ceil(H_sim/N_sample/H)), idx_shift=[1])
anim = visualize_meshrobot!(vis, model, rep_ref_traj, name=:ref, anim=anim, sample=1, α=1.0)
settransform!(vis[:ref], Translation(0.0, -0.4, 0.0))
plot_surface!(vis, s.env, xlims=[-0.5,4], ylims=[-0.4,0.4])

# high drop
conf = [-0.05, 0.6, 0.6, 0.6, 0.0, 0.30]
qinit = initial_configuration(s.model, conf...)
traj = run_policy(s, H_sim = 300, qinit=qinit)
anim = visualize_meshrobot!(vis, model, traj, name=:drop, anim=anim, sample=5, α=1.0)

# filename = "quadruped_monte_carlo_clean"
# MeshCat.convert_frames_to_video(
#     "/home/simon/Downloads/$filename.tar",
#     "/home/simon/Documents/$filename.mp4", overwrite=true)
#
# convert_video_to_gif(
#     "/home/simon/Documents/$filename.mp4",
#     "/home/simon/Documents/$filename.gif", overwrite=true)