text:
using Lasso
using GLM, Distributions, Random, LinearAlgebra, SparseArrays
using GLMNet: glmnet

datapath = joinpath(dirname(@__FILE__), "data")

testpath(T::DataType, d::Normal, l::GLM.Link, nsamples::Int, nfeatures::Int) =
    joinpath(datapath, "$(T)_$(typeof(d).name.name)_$(typeof(l).name.name)_$(nsamples)_$(nfeatures).tsv")

function makeX(ρ, nsamples, nfeatures, sparse)
    Σ = fill(ρ, nfeatures, nfeatures)
    Σ[diagind(Σ)] .= 1
    X = permutedims(rand(testrng, MvNormal(Σ), nsamples))
    sparse && (X[randperm(testrng, length(X))[1:round(Int, length(X)*0.95)]] .= 0)
    β = [(-1)^j*exp(-2*(j-1)/20) for j = 1:nfeatures]
    (X, β)
end

randdist(::Normal, x) = rand(testrng, Normal(x))
randdist(::Binomial, x) = rand(testrng, Bernoulli(x))
randdist(::Poisson, x) = rand(testrng, Poisson(x))
function genrand(T::DataType, d::Distribution, l::GLM.Link, nsamples::Int, nfeatures::Int, sparse::Bool)
    X, coef = makeX(0.0, nsamples, nfeatures, sparse)
    y = X*coef
    for i = 1:length(y)
        y[i] = randdist(d, GLM.linkinv(l, y[i]))
    end
    (X, y)
end

function gen_penalty_factors(X,nonone_penalty_factors;frac1=0.7,frac0=0.05)
    if nonone_penalty_factors
        penalty_factor = ones(size(X,2))
        nzeros = Int(floor(size(X,2)*frac0))
        nonone = Int(floor(size(X,2)*(1-frac1)))
        Random.seed!(testrng, 7337)
        penalty_factor[1:nzeros] = zeros(nzeros)
        penalty_factor[end-nonone+1:end] = rand(testrng, Float64,nonone)
        penalty_factor_glmnet = penalty_factor
    else
        penalty_factor = nothing
        penalty_factor_glmnet = ones(size(X,2))
    end
    penalty_factor, penalty_factor_glmnet
end

# Test against GLMNet
@testset "LassoPath" begin
    @testset "$(typeof(dist).name.name) $(typeof(link).name.name)" for (dist, link) in ((Normal(), IdentityLink()), (Binomial(), LogitLink()), (Poisson(), LogLink()))
        @testset "sparse = $sp" for sp in (false, true)
            Random.seed!(testrng, 465)
            (X, y) = genrand(Float64, dist, link, 1000, 10, sp)
            yoff = randn(testrng, length(y))
            @testset "$(intercept ? "w/" : "w/o") intercept" for intercept = (false, true)
                @testset "unfitted LassoPath (dofit=false)" begin
                    @testset "$(spfit ? "as SparseMatrixCSC" : "as Matrix")" for spfit in (true,false)
                        l = fit(LassoPath, spfit ? sparse(X) : X, y, dist, link,
                            intercept=intercept, dofit=false)

                        p = size(X,2)
                        if intercept
                            p += 1
                        end
                        @test size(l) == (p,1)
                        @test coef(l) == zeros(eltype(X),p,1)
                    end
                end

                @testset "alpha = $alpha" for alpha = (1, 0.5)
                    @testset "$(nonone_penalty_factors ? "non-one" : "all-one") penalty factors" for nonone_penalty_factors in (false,true)
                        penalty_factor, penalty_factor_glmnet = gen_penalty_factors(X,nonone_penalty_factors)
                        @testset "$(isempty(offset) ? "w/o" : "w/") offset" for offset = Vector{Float64}[Float64[], yoff]
                            let y=y
                                # First fit with GLMNet
                                if isa(dist, Normal)
                                    yp = isempty(offset) ? y : y + offset
                                    ypstd = std(yp, corrected=false)
                                    # glmnet does this on entry, which changes λ mappings, but not
                                    # coefficients. Should we?
                                    yp = yp ./ ypstd
                                    !isempty(offset) && (offset = offset ./ ypstd)
                                    y = y ./ ypstd
                                    g = glmnet(X, yp, dist, intercept=intercept, alpha=alpha, tol=10*eps(); penalty_factor=penalty_factor_glmnet)
                                elseif isa(dist, Binomial)
                                    yp = zeros(size(y, 1), 2)
                                    yp[:, 1] = y .== 0
                                    yp[:, 2] = y .== 1
                                    g = glmnet(X, yp, dist, intercept=intercept, alpha=alpha, tol=10*eps(),
                                               offsets=isempty(offset) ? zeros(length(y)) : offset; penalty_factor=penalty_factor_glmnet)
                                else
                                    g = glmnet(X, y, dist, intercept=intercept, alpha=alpha, tol=10*eps(),
                                               offsets=isempty(offset) ? zeros(length(y)) : offset; penalty_factor=penalty_factor_glmnet)
                                end
                                gbeta = convert(Matrix{Float64}, g.betas)

                                @testset "$(randomize ? "random" : "sequential")" for randomize = [false, true]
                                    niter = 0
                                    @testset "$(algorithm == NaiveCoordinateDescent ? "naive" : "covariance")" for algorithm = (NaiveCoordinateDescent, CovarianceCoordinateDescent)
                                        @testset "$(spfit ? "as SparseMatrixCSC" : "as Matrix")" for spfit in (true,false)
                                            criterion = :coef
                                            #  for criterion in (:coef,:obj) # takes too long for travis
                                            #      @testset "criterion = $criterion" begin
                                            if criterion == :obj
                                                irls_tol = 100*eps()
                                                cd_tol = 100*eps()
                                            else
                                                irls_tol = 10*eps()
                                                cd_tol = 10*eps()
                                            end
                                            # Now fit with Lasso
                                            l = fit(LassoPath, spfit ? sparse(X) : X, y, dist, link,
                                                    λ=g.lambda, algorithm=algorithm, intercept=intercept,
                                                    cd_tol=cd_tol, irls_tol=irls_tol, criterion=criterion, randomize=randomize,
                                                    rng=StableRNG(1337),
                                                    α=alpha, offset=offset, penalty_factor=penalty_factor)
                                            rd = (l.coefs - gbeta)./gbeta
                                            rd[.!isfinite.(rd)] .= 0
                                            println("         coefs adiff = $(maximum(abs, l.coefs - gbeta)) rdiff = $(maximum(abs, rd))")
                                            rd = (l.b0 - g.a0)./g.a0
                                            rd[.!isfinite.(rd)] .= 0
                                            println("         b0    adiff = $(maximum(abs, l.b0 - g.a0)) rdiff = $(maximum(abs, rd))")
                                            if criterion==:obj
                                                # nothing to compare results against at this point, we just make sure the code runs
                                            else
                                                # @test l.λ ≈ g.lambda rtol=5e-7
                                                @test l.coefs ≈ gbeta rtol=5e-7
                                                @test l.b0 ≈ g.a0 rtol=2e-5

                                                # Ensure same number of iterations with all algorithms
                                                if niter == 0
                                                    niter = l.niter
                                                else
                                                    @test abs(niter - l.niter) <= 10
                                                end
                                            end
                                            #     end
                                            # end
                                            # With a function to generate λ
                                            λfunc(λmax) = range(λmax, stop=g.lambda[end], length=10)
                                            lf = fit(LassoPath, spfit ? sparse(X) : X, y, dist, link,
                                                     λ=λfunc, algorithm=algorithm, intercept=intercept,
                                                     cd_tol=cd_tol, irls_tol=irls_tol, criterion=criterion, randomize=randomize,
                                                     rng=StableRNG(1337),
                                                     α=alpha, offset=offset, penalty_factor=penalty_factor)
                                            @test length(lf.λ) == 10
                                            @test_skip all(iszero, lf.coefs[:,1])
                                            @test lf.coefs[:,end] ≈ l.coefs[:,end]
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

# Test for sparse matrices

# @testset "LassoPath Zero in" begin
#     for (dist, link) in ((Normal(), IdentityLink()), (Binomial(), LogitLink()), (Poisson(), LogLink()))[2:2]
#         @testset "$(typeof(dist).name.name) $(typeof(link).name.name)" begin
#             for sp in (false, true)[2:2]
#                 Random.seed!(testrng, 1337)
#                 @testset sp ? "sparse" : "dense" begin
#                     (X, y) = genrand(Float64, dist, link, 1000, 10, sp)
#                     yoff = randn(testrng, length(y))
#                     for intercept = (false, true)[2:2]
#                         @testset "$(intercept ? "w/" : "w/o") intercept" begin
#                             @testset "unfitted LassoPath (dofit=false)" begin
#                                 for spfit in (true,false)
#                                     @testset spfit ? "as SparseMatrixCSC" : "as Matrix" begin
#                                         l = fit(LassoPath, spfit ? sparse(X) : X, y, dist, link,
#                                             intercept=intercept, dofit=false)
#
#                                         p = size(X,2)
#                                         if intercept
#                                             p += 1
#                                         end
#                                         @test size(l) == (p,1)
#                                         @test coef(l) == zeros(eltype(X),p,1)
#                                     end
#                                 end
#                             end
#                             for alpha = [1, 0.5]
#                                 @testset "alpha = $alpha" begin
#                                     for nonone_penalty_factors in (false,true)
#                                         @testset "$(nonone_penalty_factors ? "non-one" : "all-one") penalty factors" begin
#                                             penalty_factor, penalty_factor_glmnet = gen_penalty_factors(X,nonone_penalty_factors)
#                                             for offset = Vector{Float64}[Float64[], yoff]
#                                                 @testset "$(isempty(offset) ? "w/o" : "w/") offset" begin
#                                                     let y=y
#                                                         # First fit with GLMNet
#                                                         if isa(dist, Normal)
#                                                             yp = isempty(offset) ? y : y + offset
#                                                             ypstd = std(yp, corrected=false)
#                                                             # glmnet does this on entry, which changes λ mappings, but not
#                                                             # coefficients. Should we?
#                                                             yp = yp ./ ypstd
#                                                             !isempty(offset) && (offset = offset ./ ypstd)
#                                                             y = y ./ ypstd
#                                                             g = glmnet(X, yp, dist, intercept=intercept, alpha=alpha, tol=10*eps(); penalty_factor=penalty_factor_glmnet)
#                                                         elseif isa(dist, Binomial)
#                                                             yp = zeros(size(y, 1), 2)
#                                                             yp[:, 1] = y .== 0
#                                                             yp[:, 2] = y .== 1
#                                                             g = glmnet(X, yp, dist, intercept=intercept, alpha=alpha, tol=10*eps(),
#                                                                        offsets=isempty(offset) ? zeros(length(y)) : offset; penalty_factor=penalty_factor_glmnet)
#                                                         else
#                                                             g = glmnet(X, y, dist, intercept=intercept, alpha=alpha, tol=10*eps(),
#                                                                        offsets=isempty(offset) ? zeros(length(y)) : offset; penalty_factor=penalty_factor_glmnet)
#                                                         end
#                                                         gbeta = convert(Matrix{Float64}, g.betas)
#
#                                                         for randomize = [false, true]
#                                                             @testset randomize ? "random" : "sequential" begin
#                                                                 niter = 0
#                                                                 for algorithm = [NaiveCoordinateDescent, CovarianceCoordinateDescent]
#                                                                      @testset algorithm == NaiveCoordinateDescent ? "naive" : "covariance" begin
#                                                                          for spfit in (true,false)
#                                                                              @testset spfit ? "as SparseMatrixCSC" : "as Matrix" begin
#                                                                                criterion = :coef
#                                                                                #  for criterion in (:coef,:obj) # takes too long for travis
#                                                                                #      @testset "criterion = $criterion" begin
#                                                                                if criterion == :obj
#                                                                                irls_tol = 100*eps()
#                                                                                cd_tol = 100*eps()
#                                                                                else
#                                                                                irls_tol = 10*eps()
#                                                                                cd_tol = 10*eps()
#                                                                                end
#                                                                                # Now fit with Lasso
#                                                                                l = fit(LassoPath, spfit ? sparse(X) : X, y, dist, link,
#                                                                                λ=g.lambda, algorithm=algorithm, intercept=intercept,
#                                                                                cd_tol=cd_tol, irls_tol=irls_tol, criterion=criterion, randomize=randomize,
#                                                                                α=alpha, offset=offset)#, penalty_factor=penalty_factor)
#                                                                                rd = (l.coefs - gbeta)./gbeta
#                                                                                rd[!isfinite(rd)] = 0
#                                                                                # println("         coefs adiff = $(maximum(abs,l.coefs - gbeta)) rdiff = $(maximum(abs,rd))")
#                                                                                rd = (l.b0 - g.a0)./g.a0
#                                                                                rd[!isfinite(rd)] = 0
#                                                                                # println("         b0    adiff = $(maximum(abs,l.b0 - g.a0)) rdiff = $(maximum(abs,rd))")
#                                                                                if criterion==:obj
#                                                                                # nothing to compare results against at this point, we just make sure the code runs
#                                                                                else
#                                                                                # @test l.λ ≈ g.lambda rtol=5e-7
#                                                                                @test l.coefs ≈ gbeta rtol=5e-7
#                                                                                @test l.b0 ≈ g.a0 rtol=2e-5
#
#                                                                                # Ensure same number of iterations with all algorithms
#                                                                                if niter == 0
#                                                                                niter = l.niter
#                                                                                else
#                                                                                @test abs(niter - l.niter) <= 10
#                                                                                end
#                                                                                end
#                                                                                #     end
#                                                                                # end
#                                                                             end
#                                                                         end
#                                                                     end
#                                                                 end
#                                                             end
#                                                         end
#                                                     end
#                                                 end
#                                             end
#                                         end
#                                     end
#                                 end
#                             end
#                         end
#                     end
#                 end
#             end
#         end
#     end
# end


text:
function Identity()
    return Identity(())
end

function value(obj::Identity, arg0::DerivativeStructure)
    return jcall(obj, "value", DerivativeStructure, (DerivativeStructure,), arg0)
end

function value(obj::Identity, arg0::jdouble)
    return jcall(obj, "value", jdouble, (jdouble,), arg0)
end

ext: jl


text:
using JuMP, FactCheck, AmplNLWriter

# Example testing basic use of NLExpr with AmplNLWriter.jl

if !isdefined(:solver); solver = IpoptNLSolver(); end

m = Model(solver=solver)

n = 30
l = -ones(n); l[1] = 0
u = ones(n)
@variable(m, l[i] <= x[i=1:n] <= u[i])
@NLexpression(m, f1, x[1])
@NLexpression(m, g, 1 + 9 * sum{x[j] ^ 2, j = 2:n} / (n - 1))
@NLexpression(m, h, 1 - (f1 / g) ^ 2)
@NLexpression(m, f2, g * h)

setvalue(x[1], 1)
setvalue(x[2:n], zeros(n - 1))
@NLobjective(m, :Min, f2)

context("example: jump_nlexpr") do
    @fact solve(m) --> :Optimal
    @fact getvalue(x[1]) --> roughly(1.0, 1e-5)
    @fact getvalue(x[2:end]) --> roughly(zeros(n - 1), 1e-5)
    @fact getvalue(f1) --> roughly(1.0, 1e-5)
    @fact getvalue(f2) --> roughly(0.0, 1e-5)
    @fact getvalue(g) --> roughly(1.0, 1e-5)
    @fact getvalue(h) --> roughly(0.0, 1e-5)
    @fact getobjectivevalue(m) --> roughly(0.0, 1e-5)
end


text:
#=
# 3D CG reaction diffusion with tet mesh
# div(K*grad(u)) -C*u = f
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("RDtet");

@useLog("RDtetlog")

@domain(3, SQUARE, UNSTRUCTURED)
@functionSpace(LEGENDRE, 1)

@mesh("cubefine.msh")

@variable(u)
@testSymbol(v)

@boundary(u, 1, DIRICHLET, "sin(pi*0.5*x)*sin(pi*y)*sin(pi*z)")

@coefficient(K, "x*x + 4")
@coefficient(C, 10)
@coefficient(f, "x*pi*cos(pi*0.5*x)*sin(pi*y)*sin(pi*z) + (-10-(x*x+4)*pi*pi*2.25)*sin(pi*0.5*x)*sin(pi*y)*sin(pi*z)")

# PDE: div(K*grad(u)) -C*u = f
@weakForm(u, "K*dot(grad(u), grad(v)) + C*u*v + f*v")

solve(u);

# exact solution is sin(pi*0.5*x)*sin(pi*y)*sin(pi*z)
# check error
maxerr = 0;
exact(x,y,z) = sin(pi*0.5*x)*sin(pi*y)*sin(pi*z);

for i=1:size(Femshop.grid_data.allnodes,2)
    x = Femshop.grid_data.allnodes[1,i];
    y = Femshop.grid_data.allnodes[2,i];
    z = Femshop.grid_data.allnodes[3,i];
    err = abs(u.values[i] - exact(x,y,z));
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

@finalize()


text:
#=

   N-queens problem in WebPPL

   This model use the standard formulation of the n-queens problem.

   For n = 4
   Distributions of variable (num:0)
   [3, 1, 4, 2]	=>	5129 (0.5129)
   [2, 4, 1, 3]	=>	4871 (0.4871)

   But for larger n's it don't give all the correct solutions, and sometimes
   incorrect solutions.

  Cf ~/webppl/queens.wppl

=#
using Turing, StatsPlots
include("jl_utils.jl")

@model function queens(n=4)

    queen ~ filldist(DiscreteUniform(1,n),n)

    # All different
    for i in 1:n
        for j in 1:n
            if i != j
                # All different rows
                true ~ Dirac(queen[i] != queen[j])
                # Different diagonal 1
                true ~ Dirac(queen[i]+i != queen[j]+j)
                # Different diagonal 2
                true ~ Dirac(queen[i]-i != queen[j]-j)
            end
        end
    end

    return queen
end

n = 4
model = queens(n)
# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# display(chns)
# display(plot(chns))

chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
genq = generated_quantities(model, chains_params)
show_var_dist_pct(genq)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Revise
using Dice

# Convert (value, weight) tuples to input for dice discrete distribution
# Example: unpack_sparse([(0, 0.2), (3, 0.8)]) => [0.2, 0, 0, 0.8]
function unpack_sparse(sparse::Vector{Tuple{Int,Float64}})
    dice_disc = zeros(maximum(val for (val, _) in sparse) + 1)
    for (val, weight) in sparse
        dice_disc[val + 1] = weight
    end
    return dice_disc
end

# Naively calculate PMF of gcd(X, Y) given distributions of X and Y
function gcd_enumeration(
    sparse1::Vector{Tuple{Int,Float64}},
    sparse2::Vector{Tuple{Int,Float64}}
)
    gcd_dist = Dict()
    for (val1, weight1) in sparse1
        for (val2, weight2) in sparse2
            res = gcd(val1, val2)
            gcd_dist[res] = get(gcd_dist, res, 0) + weight1 * weight2
        end
    end
    return gcd_dist
end

# Generate Dice code for P[gcd(X, Y) = res] given distributions of X and Y
function code_gcd(
    sparse1::Vector{Tuple{Int,Float64}},
    sparse2::Vector{Tuple{Int,Float64}},
    res::Int
)
    discrete_dist1 = unpack_sparse(sparse1)
    discrete_dist2 = unpack_sparse(sparse2)
    @dice begin
        function discrete(p::Vector{Float64})
            mb = length(p)
            v = Vector(undef, mb)
            sum = 1
            for i=1:mb
                v[i] = p[i]/sum
                sum = sum - p[i]
            end

            ans = DistInt(dicecontext(), mb-1)
            for i=mb-1:-1:1
                ans = if flip(v[i]) DistInt(dicecontext(), i-1) else ans end
            end
            return ans
        end

        # https://en.wikipedia.org/wiki/Euclidean_algorithm#Implementations
        function dice_gcd(a::DistInt, b::DistInt)
            a = a, flip(0)
            b = b, flip(0)
            for _ = 1 : 1 + max_bits(b[1]) ÷ log2(MathConstants.golden)
                b_zero = prob_equals(b[1], 0)
                t = b
                amb = (a[1] % b[1])
                b = Dice.ifelse(b_zero, b[1], amb[1]), Dice.ifelse(b_zero, b[2], amb[2] | a[2] | b[2])
                a = Dice.ifelse(b_zero, a[1], t[1]), Dice.ifelse(b_zero, a[2], t[2])
            end
            return a[1], a[2]
        end

        a = discrete(discrete_dist1)
        b = discrete(discrete_dist2)
        g = dice_gcd(a, b)[1]
    end
end

function test(sparse::Vector{Tuple{Int,Float64}})
    gcd_dist = gcd_enumeration(sparse, sparse)
    code = code_gcd(sparse, sparse, 4)
    bdd = compile(code)
    dice_p = infer(bdd)
    for res in Set(Iterators.flatten((keys(gcd_dist), 0:1)))
        @assert dice_p[res + 1] ≈ get(gcd_dist, res, 0)
    end
end

test([(4, 0.3), (6, 0.2), (9, 0.5)])  # 4 has probability 0.3, etc...
test([(0, 0.2), (4, 0.3), (6, 0.5)])
test([(0, 0.1), (2, 0.2), (3, 0.7)])
test([(0, 0.1), (8, 0.2), (13, 0.7)])
test([(0, 0.1), (34, 0.2), (55, 0.7)])
test([(144, 0.1), (233, 0.2), (610, 0.3), (987, 0.4)])


text:
using YaoBlocksQobj, DocThemeIndigo
using Documenter
indigo = DocThemeIndigo.install(YaoBlocksQobj)

makedocs(;
    modules = [YaoBlocksQobj],
    authors = "Arsh Sharma <sharmarsh15@gmail.com> and contributors",
    repo = "https://github.com/QuantumBFS/YaoBlocksQobj.jl/blob/{commit}{path}#{line}",
    sitename = "YaoBlocksQobj.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://QuantumBFS.github.io/YaoBlocksQobj.jl",
        assets = String[indigo, "assets/default.css"],
    ),
    pages = ["Quickstart" => "index.md", "API References" => "refs.md"],
)

deploydocs(; repo = "github.com/QuantumBFS/YaoBlocksQobj.jl", devbranch = "main")


text:
"""
build_system(
    category::Type{<:SystemCategory},
    name::String,
    print_stat::Bool = false;
    kwargs...,
)

Accepted Key Words:
- `add_forecasts::Bool`: Default is `true`
- `add_reserves::Bool`: Default is `false`
- `force_build::Bool`: `true` runs entire build process, `false` (Default) uses deserializiation if possible
- `skip_serialization::Bool`: Default is `false`
- `system_catelog::SystemCatalog`
- `assign_new_uuids::Bool`: Assign new UUIDs to the system and all components if
   deserialization is used. Default is `true`.
"""
function build_system(
    category::Type{<:SystemCategory},
    name::String,
    print_stat::Bool = false;
    assign_new_uuids = true,
    kwargs...,
)
    add_forecasts = get(kwargs, :add_forecasts, true)
    add_reserves = get(kwargs, :add_reserves, false)
    force_build = get(kwargs, :force_build, false)
    skip_serialization = get(kwargs, :skip_serialization, false)
    system_catelog = get(kwargs, :system_catelog, SystemCatalog(SYSTEM_CATALOG))
    sys_descriptor = get_system_descriptor(category, system_catelog, name)
    if !is_serialized(name, add_forecasts, add_reserves) || force_build
        check_serialized_storage()
        download_function = get_download_function(sys_descriptor)
        if !isnothing(download_function)
            filepath = download_function(; kwargs...)
            set_raw_data!(sys_descriptor, filepath)
        end
        @info "Building new system $(sys_descriptor.name) from raw data" sys_descriptor.raw_data
        build_func = get_build_function(sys_descriptor)
        start = time()
        sys = build_func(; raw_data = sys_descriptor.raw_data, kwargs...)
        construct_time = time() - start
        serialized_filepath = get_serialized_filepath(name, add_forecasts, add_reserves)
        start = time()
        if !skip_serialization
            PSY.to_json(sys, serialized_filepath; force = true)
            serialize_time = time() - start
        end
        # set_stats!(sys_descriptor, SystemBuildStats(construct_time, serialize_time))
    else
        @debug "Deserialize system from file" sys_descriptor.name
        start = time()
        # time_series_in_memory = get(kwargs, :time_series_in_memory, false)
        sys_kwargs = filter_kwargs(; kwargs...)
        file_path = get_serialized_filepath(name, add_forecasts, add_reserves)
        sys = PSY.System(file_path; assign_new_uuids = assign_new_uuids, sys_kwargs...)
        PSY.get_runchecks(sys)
        # update_stats!(sys_descriptor, time() - start)
    end
    # print_stat ? print_stats(sys_descriptor) : nothing
    return sys
end


text:
using Revise

using Pedestrians
using Plots

# create model
k = 500
model = build_model()

@time anim = @animate for i in 1:k
    title = string("Blind velocity: ", lpad(i, 3, " "))
    makeplot(model; title, addview = true)
    simulation_step!(model)
end;
gif(anim, "./assets/blind_velocity.gif", fps = 10)


text:
# Import modules
@everywhere using Pkg
@everywhere Pkg.activate("../..")
@everywhere using Distributions  # probability distributions and associated functions
@everywhere using StatsBase
@everywhere using LinearAlgebra
# Import Calibrate-Emulate-Sample modules
@everywhere using EnsembleKalmanProcesses.EnsembleKalmanProcessModule
@everywhere using EnsembleKalmanProcesses.Observations
@everywhere using EnsembleKalmanProcesses.ParameterDistributionStorage
@everywhere include(joinpath(@__DIR__, "helper_funcs.jl"))
using JLD
@everywhere using NPZ

###
###  Define the parameters and their priors
###

# Define the parameters that we want to learn

@everywhere param_names = ["entrainment_factor", "detrainment_factor", "sorting_power",
	"tke_ed_coeff", "tke_diss_coeff", "pressure_normalmode_adv_coeff",
         "pressure_normalmode_buoy_coeff1", "pressure_normalmode_drag_coeff", "static_stab_coeff"]
@everywhere n_param = length(param_names)

# Assume lognormal priors for all parameters
# Note: For the EDMF model to run, all parameters need to be nonnegative.
# The ekp update can result in violations of
# these constraints - therefore, we perform CES in log space, i.e.,
# (the parameters can then simply be obtained by exponentiating the final results).

# Prior: Log-normal in original space defined by mean and std
logmeans = zeros(n_param)
log_stds = zeros(n_param)
logmeans[1], log_stds[1] = logmean_and_logstd(0.3, 0.2)
logmeans[2], log_stds[2] = logmean_and_logstd(0.3, 0.2)
logmeans[3], log_stds[3] = logmean_and_logstd(2.0, 1.0)
logmeans[4], log_stds[4] = logmean_and_logstd(0.4, 0.2)
logmeans[5], log_stds[5] = logmean_and_logstd(0.2, 0.2)
logmeans[6], log_stds[6] = logmean_and_logstd(0.2, 0.2)
logmeans[7], log_stds[7] = logmean_and_logstd(0.2, 0.2)
logmeans[8], log_stds[8] = logmean_and_logstd(8.0, 1.0)
logmeans[9], log_stds[9] = logmean_and_logstd(0.4, 0.2)
prior_dist = [Parameterized(Normal(logmeans[1], log_stds[1])),
                        Parameterized(Normal(logmeans[2], log_stds[2])),
                        Parameterized(Normal(logmeans[3], log_stds[3])),
                        Parameterized(Normal(logmeans[4], log_stds[4])),
                        Parameterized(Normal(logmeans[5], log_stds[5])),
                        Parameterized(Normal(logmeans[6], log_stds[6])),
                        Parameterized(Normal(logmeans[7], log_stds[7])),
                        Parameterized(Normal(logmeans[8], log_stds[8])),
                        Parameterized(Normal(logmeans[9], log_stds[9]))]
@everywhere prior_dist = $prior_dist

###
###  Retrieve true LES samples from PyCLES data
###
# This is the true value of the observables (e.g. LES horizontal mean)
@everywhere data_dir = "/groups/esm/ilopezgo/padeops_data/"
@everywhere padeops_t = npzread( string(data_dir,"time_padeops.npy") )*3600.0 # s
@everywhere padeops_z = npzread( string(data_dir,"zCell_padeops.npy") )*1000.0 # m
@everywhere padeops_theta = npzread( string(data_dir,"potT_padeops.npy") )
@everywhere padeops_uh = npzread( string(data_dir,"wind_speed_padeops.npy") )

# Times on which to interpolate
@everywhere t_fig3 = ([4, 8, 12, 16, 20])*3600.0
# @everywhere t_fig3 = ([4, 6, 8, 10, 12, 15.5, 17, 18.5, 20, 21.5])*3600.0
# @everywhere t_fig5b = ([6, 7, 8, 9, 10, 11, 12, 13])*3600.0
# Get SCM vertical grid
@everywhere sim_names = ["Kumar_dc_init_m1"]
sim_dir = string("Output.", sim_names[1],".00000")
z_scm = get_profile(sim_dir, ["z_half"])
# Initialize objectives
@everywhere yt = zeros(0)
yt_var_list = []

for t_ in t_fig3
    yt_, yt_var_ = padeops_m_σ2(padeops_theta, padeops_z, padeops_t, z_scm, t_)
    append!(yt, yt_)
    push!(yt_var_list, yt_var_)
    yt_, yt_var_ = padeops_m_σ2(padeops_uh, padeops_z, padeops_t, z_scm, t_)
    append!(yt, yt_)
    push!(yt_var_list, yt_var_)
end

@everywhere yt = $yt
yt_var = zeros(length(yt), length(yt))
vars_num = 1
for sim_covmat in yt_var_list
    vars = length(sim_covmat[1,:])
    yt_var[vars_num:vars_num+vars-1, vars_num:vars_num+vars-1] = sim_covmat
    global vars_num = vars_num+vars
end
@everywhere yt_var = $yt_var
@everywhere n_observables = length(yt)
padeops_names = Array{String, 1}[]
push!(padeops_names, ["theta_fig3", "uh"])
@everywhere padeops_names=$padeops_names

# This is how many samples of the true data we have
n_samples = 1
samples = zeros(n_samples, length(yt))
samples[1,:] = yt
# Noise level of the samples, which scales the time variance of each output.
noise_level = 0.005
Γy = yt_var + Matrix(noise_level*I, size(yt_var)[1], size(yt_var)[2]) #Uncertainty inflation (nugget)
@everywhere Γy = $Γy
μ_noise = zeros(length(yt))
# We construct the observations object with the samples and the cov.
truth = Obs(Array(samples'), Γy, padeops_names[1])
@everywhere truth = $truth

###
###  Calibrate: Ensemble Kalman Inversion
###

@everywhere N_ens = 50 # number of ensemble members
@everywhere N_iter = 20 # number of ekp iterations.
@everywhere N_yt = length(yt) # Length of data array

@everywhere constraints = [[no_constraint()], [no_constraint()],
                [no_constraint()], [no_constraint()],
                [no_constraint()], [no_constraint()],
                [no_constraint()], [no_constraint()],[no_constraint()]]

@everywhere priors = ParameterDistribution(prior_dist, constraints, param_names)
@everywhere initial_params = construct_initial_ensemble(priors, N_ens)
@everywhere y_names = ["thetal_mean", "horizontal_vel"]
precondition_ensemble!(initial_params, priors, param_names, y_names, t_fig3)
@everywhere initial_params = $initial_params

# Check dimensionality
@everywhere N_par_, N_ens_ = size(initial_params)
@assert N_par_ == length(constraints)
@assert N_ens == N_ens_
@everywhere ekobj = EnsembleKalmanProcess(initial_params, yt, Γy, Inversion())

g_ens = zeros(N_ens, n_observables)
@everywhere scm_dir = "/home/ilopezgo/SCAMPy/"
@everywhere g_(x::Array{Float64,1}) = run_SCAMPy(x, param_names, y_names, scm_dir, t_fig3)

outdir_path = string("results_kumar_p", n_param,"_n", noise_level,"_e",
    N_ens, "_i", N_iter, "_d", N_yt)
command = `mkdir $outdir_path`
try
    run(command)
catch e
    println("Output directory already exists. Output may be overwritten.")
end

# ekp iterations
@everywhere Δt = 1.0
for i in 1:N_iter
    # Note that the parameters are exp-transformed when used as input
    # to SCAMPy
    @everywhere params_i = deepcopy(exp_transform(get_u_final(ekobj)))
    @everywhere params = [row[:] for row in eachrow(params_i')]
    g_ens_arr = pmap(g_, params)
    println(string("\n\nekp evaluation ",i," finished. Updating ensemble ...\n"))
    for j in 1:N_ens
      g_ens[j, :] = g_ens_arr[j]
    end
    update_ensemble!(ekobj, Array(g_ens') )
    println("\nEnsemble updated.\n")
    # Save ekp information to file
    save( string(outdir_path,"/ekp.jld"), "ekp_u", get_u(ekobj), "ekp_g", get_g(ekobj),
        "truth_mean", ekobj.obs_mean, "truth_cov", ekobj.obs_noise_cov, "ekp_err", ekobj.err)
end

# ekp results: Has the ensemble collapsed toward the truth? Store and analyze.
println("\nekp ensemble mean at last stage (original space):")
println(mean(deepcopy(exp_transform(get_u_final(ekobj) )), dims=1))

println("\nEnsemble covariance det. 1st iteration, transformed space.")
println(det(cov(deepcopy((get_u(ekobj, 1) )), dims=1)))

println("\nEnsemble covariance det. last iteration, transformed space.")
println(det(cov(deepcopy((get_u_final(ekobj) )), dims=1)))


text:
"""
    Statistics

"""
mutable struct Statistics{T<:AbstractFloat}
    flag::Int   # -1 mean 'not yet computed'
    rank::Int
    rank_spqr::Int
    tol::T    # to be copied from opts
    normest_A::T  # to be copied from options
    tol_alt::T
    est_sval_upper_bounds::Vector{T}
    est_sval_lower_bounds::Vector{T}

    est_svals::Vector{T}
    est_svals_of_R::Vector{T}
    est_error_bounds::Vector{T}
    iters::Int
    time::Int64
    time_initialize::Int64
    time_iters::Int64
    time_est_error_bounds::Int64
    time_svd::Int64
    time_basis::Int64
    norm_R_times_N::T
    norm_R_transpose_times_NT::T
    est_norm_A_times_N::T
    est_err_bound_norm_A_times_N::T
    est_norm_A_transpose_times_NT::T
    est_err_bound_norm_A_transpose_times_NT::T
    nsvals_large_found::Int
    final_blocksize::Int
    ssi_min_block_used::Int
    ssi_max_block_used::Int
    sval_numbers_for_bounds::AbstractUnitRange{Int}
    opts_used::NamedTuple
    info_spqr::NamedTuple
    stats_ssi::Union{Missing,Statistics{T}}
    stats_ssp_N::Union{Missing,Statistics{T}}
    stats_ssp_NT::Union{Missing,Statistics{T}}
end

function Statistics(::Type{T}; tol::T=0.0, normest_A::T=0.0) where T
    Statistics{T}(-1,
                  -1, -1,
                  tol, normest_A, -1.0,
                  T[], T[], T[], T[], T[],
                  0,
                  0, 0, 0, 0, 0, 0,
                  0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
                  0, -1, -1, -1, 0:-1,
                  NamedTuple(), NamedTuple(), missing, missing, missing)
end


text:
export read_con

"""
Use:   read_con(con::SeisCon;
                blocks::Array{Int,1} = Array(1:length(con)),
                prealloc_traces::Int = 50000)

Read 'blocks' out of 'con' into a preallocated array of size (ns x prealloc_traces).

If preallocated memory fills, it will be expanded again by 'prealloc_traces'.
"""
function read_con(con::SeisCon, blocks::Array{Int,1};
                                prealloc_traces::Int = 50000)
    nblocks = length(blocks)
    if maximum(blocks)>length(con) @error "Call for block $(maximum(blocks)) in a container with $(length(con)) blocks." end

    # Check dsf
    datatype = Float32
    if con.dsf == 1
        datatype = IBMFloat32
    elseif con.dsf != 5
        @error "Data type not supported ($(fh.bfh.DataSampleFormat))"
    end

    # Pre-allocate
    data = Array{datatype,2}(undef, con.ns, prealloc_traces)
    headers = Array{BinaryTraceHeader,1}(undef, prealloc_traces)
    fh = FileHeader(); set_fileheader!(fh.bfh, :ns, con.ns)
    set_fileheader!(fh.bfh, :DataSampleFormat, con.dsf)

    trace_count = 0
    # Read all blocks
    for block in blocks

        # Check size of next block and pass view to pre-alloc
        brange = con.blocks[block].endbyte - con.blocks[block].startbyte
        ntraces = Int((brange)/(240 + con.ns*4))

        # Check if there is room in pre-alloc'd mem
        isroom = (trace_count + ntraces) <= length(headers)
        if ~isroom
            println("Expanding preallocated memory")
            prealloc_traces *= 2
            data = hcat(data, Array{datatype,2}(undef, con.ns, ntraces+prealloc_traces))
            append!(headers, Array{BinaryTraceHeader,1}(undef, ntraces+prealloc_traces))
        end
        tmp_data = view(data, :,(trace_count+1):(trace_count+ntraces))
        tmp_headers = view(headers, (trace_count+1):(trace_count+ntraces))

        # Read the next block
        read_block!(con.blocks[block], con.ns, con.dsf, tmp_data, tmp_headers)
        trace_count += ntraces
    end

    return SeisBlock{datatype}(fh, headers[1:trace_count], data[:,1:trace_count])

end

function read_con(con::SeisCon, keys::Array{String,1}, blocks::Array{Int,1};
                                prealloc_traces::Int = 50000)
    nblocks = length(blocks)

    # Check dsf
    datatype = Float32
    if con.dsf == 1
        datatype = IBMFloat32
    elseif con.dsf != 5
        @error "Data type not supported ($(fh.bfh.DataSampleFormat))"
    end

    # Check for RecSrcScalar
    in("RecSourceScalar", keys) ? nothing : push!(keys, "RecSourceScalar")

    # Pre-allocate
    data = Array{datatype,2}(undef, con.ns, prealloc_traces)
    headers = Array{BinaryTraceHeader,1}(undef, prealloc_traces)
    fh = FileHeader(); set_fileheader!(fh.bfh, :ns, con.ns)
    set_fileheader!(fh.bfh, :DataSampleFormat, con.dsf)

    trace_count = 0
    # Read all blocks
    for block in blocks

        # Check size of next block and pass view to pre-alloc
        brange = con.blocks[block].endbyte - con.blocks[block].startbyte
        ntraces = Int((brange)/(240 + con.ns*4))

        # Check if there is room in pre-alloc'd mem
        isroom = (trace_count + ntraces) <= length(headers)
        if ~isroom
            println("Expanding preallocated memory")
            data = hcat(data, Array{datatype,2}(undef, con.ns, ntraces+prealloc_traces))
            append!(headers, Array{BinaryTraceHeader,1}(undef, ntraces+prealloc_traces))
            prealloc_traces *= 2
        end
        tmp_data = view(data, :,(trace_count+1):(trace_count+ntraces))
        tmp_headers = view(headers, (trace_count+1):(trace_count+ntraces))

        # Read the next block
        read_block!(con.blocks[block], keys, con.ns, con.dsf, tmp_data, tmp_headers)
        trace_count += ntraces
    end

    return SeisBlock{datatype}(fh, headers[1:trace_count], data[:,1:trace_count])

end

# RANGES & INT
function read_con(con::SeisCon, blocks::TR;
                  prealloc_traces::Int = 50000) where {TR<:AbstractRange}
    read_con(con, Array(blocks), prealloc_traces = prealloc_traces)
end
function read_con(con::SeisCon, blocks::Integer;
                  prealloc_traces::Int = 50000)
    read_con(con, [blocks], prealloc_traces = prealloc_traces)
end
function read_con(con::SeisCon, keys::Array{String,1}, blocks::TR;
                  prealloc_traces::Int = 50000) where {TR<:AbstractRange}
    read_con(con, keys, Array(blocks), prealloc_traces = prealloc_traces)
end
function read_con(con::SeisCon, keys::Array{String,1}, blocks::Integer;
                  prealloc_traces::Int = 50000)
    read_con(con, keys, [blocks], prealloc_traces = prealloc_traces)
end


text:
function BilinearInterpolator()
    return BilinearInterpolator(())
end

function interpolate(obj::BilinearInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::Vector{Vector{jdouble}})
    return jcall(obj, "interpolate", BilinearInterpolatingFunction, (Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1, arg2)
end


text:
@testset "Virtual Brownian Tree tests" begin
  using DiffEqNoiseProcess, DiffEqBase, Test, Random
  import RandomNumbers, Random123

  W = VirtualBrownianTree(0.0,0.0; tree_depth=3, search_depth=5)
  @test isinplace(W) == false
  @test length(W.seeds) == 2^(W.tree_depth)
  @test length(W.W) == 2^(W.tree_depth)+1
  @test length(W.t) == 2^(W.tree_depth)+1
  @test length(unique(W.seeds)) == length(W.seeds)

  # step to cached value
  dt = W.t[2]-W.t[1]
  calculate_step!(W,dt,nothing,nothing)
  @test W.dW == W.W[2]-W.W[1]

  for i in 1:length(W.t)-1
    accept_step!(W,dt,nothing,nothing)
  end
  @test_throws ErrorException accept_step!(W,dt,nothing,nothing)

  # test interpolation/binary search
  rng = RandomNumbers.Random123.Threefry4x() # instantiate PRNG
  rngcopy = copy(rng)
  W = VirtualBrownianTree(0.0,0.0; tree_depth=1, search_depth=5, rng=rng)
  t = (W.t[2]-W.t[1])/8
  # VBT but with deeper tree such that no binary search is required
  Wcopy = VirtualBrownianTree(0.0,0.0; tree_depth=4, search_depth=5, rng=rngcopy)
  @test W(t) == Wcopy(t)

  W = VirtualBrownianTree(0.0,0.0; tree_depth=1, rng=rng)
  for i in 1:10
    accept_step!(W,0.1,nothing,nothing)
  end

  # test NoiseProblem
  W = VirtualBrownianTree(0.0,0.0; tree_depth=1, rng=rng)
  prob = NoiseProblem(W,(0.0,1.0))
  sol = solve(prob;dt=0.05)

  # test rational numbers
  W = VirtualBrownianTree(0//1,0.0; tree_depth=1, rng=rng)
  prob = NoiseProblem(W,(0,1))
  sol = solve(prob;dt=1//2)

  W = VirtualBrownianTree!(0.0,zeros(8); tree_depth=3, search_depth=5)
  @test isinplace(W) == true
  dt = W.t[2]-W.t[1]
  calculate_step!(W,dt,nothing,nothing)
  @test W.dW == W.W[2]-W.W[1]

  for i in 1:length(W.t)-1
    accept_step!(W,dt,nothing,nothing)
  end

  # test memory size
  w0 = zeros(1_000_000)
  W = VirtualBrownianTree(0.0,w0; tree_depth=0)
  # size = 4*W0 to store W0, Wend, curW, and dW
  @test isapprox(Base.summarysize(W)/Base.summarysize(w0), 4, rtol=1e-4)

  prob = NoiseProblem(W,(0.0,1.0))
  sol = solve(prob;dt=1/10)

  @test isapprox(Base.summarysize(sol)/Base.summarysize(w0), 4, rtol=1e-4)

end


text:
module TestDiffView

using Test
using LeftChildRightSiblingTrees
using ProfileVega
using VegaLite
using ProfileVega: islastsibling, lcs

eq(x::Node{Int}, y::Node{Int}) = x.data == y.data

@testset "lcs basics" begin
    root1 = Node(1)
    root2 = Node(1)

    res = lcs(root1, root2, eq)
    @test res.l == 1
    @test length(res.v) == 1
    @test res.v[1] == (root1, root2)

    root1 = Node(1)
    root2 = Node(2)

    res = lcs(root1, root2, eq)
    @test res.l == 0
    @test isempty(res.v)
end

@testset "lcs simple sequences" begin
    x1 = Node(1)
    x2 = addsibling(x1, 2)

    y1 = Node(1)
    y2 = addsibling(y1, 2)
    res = lcs(x1, y1, eq)
    @test res.l == 2

    x1 = Node(1)
    x2 = addsibling(x1, 2)
    x3 = addsibling(x2, 3)

    y1 = Node(1)
    y2 = addsibling(y1, 3)

    res = lcs(x1, y1, eq)
    @test res.l == 2
end

@testset "lcs simple wildcard sequences" begin
    wildcard(x) = x.data == 0

    x1 = Node(1)
    x2 = addsibling(x1, 0)
    x3 = addsibling(x2, 2)
    x4 = addsibling(x3, 3)
    x2_1 = addchild(x2, 1)
    x2_2 = addsibling(x2_1, 2)
    x2_3 = addsibling(x2_2, 3)

    y1 = Node(1)
    y2 = addsibling(y1, 2)
    y3 = addsibling(y2, 3)
    y4 = addsibling(y3, 0)
    y4_1 = addchild(y4, 1)
    y4_2 = addsibling(y4_1, 2)
    y4_3 = addsibling(y4_2, 3)

    res = lcs(x1, y1, eq, wildcard)
    @test length(res.v) == 2
    @test res.l == 5
    @test res.v[1] == (x1, y1)
    @test res.v[2] == (x2, y4)
end

function f1(n)
    res = 0
    for _ in 1:n
        res += sum(rand(10_000))
    end

    for _ in 1:n
        res -= sum(rand(10_000))
    end
    res
end

function f2(n)
    res = 0
    for _ in 1:n
        res += sum(rand(20_000))
    end

    for _ in 1:n
        A = randn(10, 10, 10)
        res += maximum(A)
    end

    for _ in 1:n
        res -= sum(rand(5_000))
    end
    res
end

@testset "ProfileVega diffview" begin
    baseline = @profbase f1(2)
    p1 = @profdiffview baseline f2(2)

    @test p1 isa VegaLite.VLSpec
end

end # module

max_stars_repo_name: davidanthoff/ProfileVega

max_issues_repo_name: davidanthoff/ProfileVega

max_forks_repo_name: davidanthoff/ProfileVega


text:
using Quadrant_Data
using Test

@testset "Quadrant_Data.jl" begin
    # Write your tests here.
    xtraj = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    ytraj = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,]
    arraylist = ["quad 1", "quad 2", "quad 2", "quad 3", "quad 4", "quad 2", "quad 2", "quad 1", "quad 2", "quad 2", "quad 3", "quad 3", "quad 3", "quad 4", "quad 2", "quad 2", "quad 2", "quad 1", "quad 1", "quad 2", "quad 2", "quad 2", "quad 2", "quad 4", "quad 1", "quad 3", "quad 4", "quad 2", "quad 2", "quad 3", "quad 1", "quad 2", "quad 2", "quad 3", "quad 2", "quad 1", "quad 2", "quad 3", "quad 2", "quad 4", "quad 2", "quad 4", "quad 2", "quad 2", "quad 2", "quad 2", "quad 2", "quad 2", "quad 1", "quad 3", "quad 4", "quad 2", "quad 2", "quad 2", "quad 2", "quad 2", "quad 1", "quad 2", "quad 4"]
    quadrantdf = DataFrame(x = xtraj, y = ytraj, array = arraylist)
    @test calculatepvalue(quadrantdf, "quad 2") <= 0.05
end


text:
"""

    SuSiEGLMM

Sum of Single Effects model based on Generalized Linear Mixed Model for binary traits.

"""
module SuSiEGLMM

using Statistics, LinearAlgebra, Random, StatsBase, Distributions, Distributed,DelimitedFiles

### include function files
# include("Utilities.jl")

# include("VEM.jl")

include("GLMM.jl")

### call modules and export particular functions

end


text:
# In this example we are going to ingest and solve a real dataset from a Brookstone rover.

# 0. We need additional libraries to ingest the raw LCM datafile. Minimal installs, definitely worth the data :)
#  - Note that a full Caesar install will give all these to you (later examples), but in case you don't have them, let's go through the steps.
# ] then add LCMCore#v0.5.1 #Unsafe arrays don't work yet.
# ] then dev CaesarLCMTypes

using LCMCore, CaesarLCMTypes
using JSON
using GraffSDK

# 1. Import the initialization code.
cd(joinpath(dirname(pathof(GraffSDK)), "..", "examples"))
# 1a. Constants
graffConfig = loadGraffConfig("synchronyConfig.json")
graffConfig.sessionId = "Hackathon809"
graffConfig.robotId = "Brookstone"

# 2. Confirm that the robot already exists, create if it doesn't.
println(" - Creating or retrieving robot '$(graffConfig.robotId)'...")
robot = nothing
if(GraffSDK.isRobotExisting())
    println(" -- Robot '$(graffConfig.robotId)' already exists, retrieving it...")
    robot = getRobot()
else
    # Create a new one
    println(" -- Robot '$(graffConfig.robotId)' doesn't exist, creating it...")
    newRobot = RobotRequest(graffConfig.robotId, "My New Bot", "Description of my neat robot", "Active")
    robot = addRobot(graffConfig, newRobot)
end

# 3. Create or retrieve the session.
# Get sessions, if it already exists, add to it.
println(" - Creating or retrieving data session '$(graffConfig.sessionId)' for robot...")
session = nothing
if(GraffSDK.isSessionExisting())
    println(" -- Session '$(graffConfig.sessionId)' already exists for robot '$(graffConfig.robotId)', retrieving it...")
    session = getSession()
else
    # Create a new one
    println(" -- Session '$(graffConfig.sessionId)' doesn't exist for robot '$(graffConfig.robotId)', creating it...")
    newSessionRequest = SessionDetailsRequest(graffConfig.sessionId, "Brookstone test dataset - driving in a loop around the Work DoneGin.", "Pose2")
    session = addSession(newSessionRequest)
end
println(session)

# Adding the data!

# 5. Set up the LCM callbacks
function lcm_NewOdoAvailable(channel, msg::brookstone_supertype_t, graffConfig::GraffConfig)
    # @show variable = msg.newvariable
    factor = msg.newfactor
    img = msg.img
    # Factor
    dataFrame = JSON.parse(String(factor.data))
    measurement = Float64.(dataFrame["meas"])
    @show pOdo = Float64.(reshape(dataFrame["podo"], 3, 3))
    newOdometryMeasurement = AddOdometryRequest(measurement, pOdo)
    @time @show response = addOdometryMeasurement(newOdometryMeasurement)

    # Image
    # TODO: Finding the ID for the node - this shouldn't be necessary! Fix once we have response.variable.id
    # nodeLabel = response.variable.label
    # newNode = getNode(graffConfig, robotId, sessionId, nodeLabel)
    # nodeId = newNode.id
    # # Make a request payload
    # dataElementRequest = BigDataElementRequest("CamImage", "Mongo", "Brookstone camera data for timestamp $(factor.utime)", base64encode(img.data), "image/jpeg")
    #   # add the data to the database
    # addOrUpdateDataElement(graffConfig, robotId, sessionId, nodeId, dataElementRequest)
end

lcmLog = LCMLog("data/brookstone_3.lcmlog")
subscribe(lcmLog, "BROOKSTONE_ROVER", (c,m)->lcm_NewOdoAvailable(c,m, graffConfig), brookstone_supertype_t)
# subscribe(lcm, "BROOKSTONE_NEW_FACTOR", (c,m) -> newFactor_Callback(c, m, runtimeInfo), generic_factor_t)
# subscribe(lcm, "BROOKSTONE_CAM_IMAGE", (c,m)->keyframe_Callback(c,m, graffConfig, nodehist), image_t)

# Run while there is data
while handle(lcmLog)
end

# KAMEHAMEHA :D
# NOW PLAY THIS: https://www.youtube.com/watch?v=UT9w0PGykZ0
# How awesome is that?

# Let's take a look at the data:
visualizeSession(graffConfig, robotId, sessionId, "CamImage")


text:
# extract_spectrograms.jl
# This script takes audio recordings and dumps the corresponding spectrograms
# into trainpath, defined in the config file.
# author: Timothy Evans
# date: 30/07/14

push!(LOAD_PATH, pwd())
using DSP, WAV, DataFrames, Images
using config, dataprocessing, ProgressMeter

#load paths from config.jl
specpath,labelledpath,trainpath,modelpath = loadpaths()

#= loop through the levels of recording directories. i.e. BSR01/140107/concat/
 NOTE: this is as per the file structure received from Eric Watson with the
 Ark in the Park. If structure is different the following nested loop needs to
 be adapted. =#

fss = Uint[]
p = Progress(length(readdir(specpath)),1) # initialise progress bar

for recorder in readdir(wavpath)
  recdatepath = string(wavpath,recorder,"/")
  for day in readdir(recdatepath)
    recpath = string(recdatepath,day,"/concat/")
    for file in readdir(recpath)
      filepath = string(recpath,file)
      rec,fs = wavread(filepath)
      frec = filterwav(rec,fs)
      #spectrogram segments will be spaced as per defined speclength
      segments = 1:fs*speclength:length(frec)
      specname = string(recorder,"_",day,"_",file[1],"_")
      dumpspectrograms(frec,segments,trainpath,specname,specwidth,p)
    end
  end
  push!(fss,fs)
end

recorderdata = DataFrame(recorders=readdir(wavpath),fsample=fss)
writetable(string(modelpath,"recorderdata.csv"),recorderdata)

max_stars_repo_name: kokakodetect/kokako-detection

max_issues_repo_name: kokakodetect/kokako-detection

max_forks_repo_name: kokakodetect/kokako-detection


text:
"""
    DQMCParameters([p::DQMCParameters]; [kwargs...])

Specifies the parameters of the determinant quantum Monte-Carlo (DQMC)
simulation.

If an old set of parameters `p` is passed a copy will be generated. Any
parameter passed via kwargs will replace the respective parameter in the copy.

## Keyword Arguments:
- `thermalization = 100`: Number of thermalization sweeps. (First MC stage)
- `sweeps = 100`: Number of measurement sweeps. (Second MC stage)
- `silent = false`: Should error printing be supressed? (They are still recorded
 in DQMCAnalysis.)
- `check_sign_problem = true`: Should check for sign problems happen? (This fully
enables/disables sign problem checks, including those for DQMCAnalysis!)
- `check_propagation_error = true`: Should checks for propagation errors happen?
(This fully enables/disables propgation error checks, includign those for
DQMCAnalysis!)
- `safe_mult = 10`: Number of safe matrix multiplications in `propagate`. If you
see frequent/large propagation errors you may want to lower this number. Note
that the number of time slices must be divisable by `safe_mult`.
- `delta_tau`: The imaginary time discretization of the path integral.
- `beta`: The inverse temperature.
- `slices`: The total number of time slices.
- `measure_rate = 10`: Perform measurements every `measure_rate` sweeps.
- `print_rate = 10`: Print general information every `print_rate` sweeps.

If you construct `DQMCParameters` from scratch you need to specify two of
`delta_tau`, `beta` and `slices`. The following formula must hold:
`beta = delta_tau * slices`.
"""
struct DQMCParameters
    thermalization::Int
    sweeps::Int

    silent::Bool
    check_sign_problem::Bool
    check_propagation_error::Bool

    safe_mult::Int
    delta_tau::Float64
    beta::Float64
    slices::Int

    measure_rate::Int
    print_rate::Int
end

function DQMCParameters(
        p::DQMCParameters;
        thermalization = p.thermalization,
        sweeps = p.sweeps,
        silent = p.silent,
        check_sign_problem = p.check_sign_problem,
        check_propagation_error = p.check_propagation_error,
        safe_mult = p.safe_mult,
        delta_tau = p.delta_tau,
        measure_rate = p.measure_rate,
        print_rate = p.print_rate,
        kwargs...
    )
    if (!haskey(kwargs, :beta)) && (!haskey(kwargs, :slices))
        kwargs = Dict{Symbol, Any}(kwargs)
        push!(kwargs, :beta => p.beta)
    end
    DQMCParameters(
        thermalization = thermalization,
        sweeps = sweeps,
        silent = silent,
        check_sign_problem = check_sign_problem,
        check_propagation_error = check_propagation_error,
        safe_mult = safe_mult,
        delta_tau = delta_tau,
        measure_rate = measure_rate,
        print_rate = print_rate;
        kwargs...
    )
end

function DQMCParameters(;
        thermalization::Int = 100,
        sweeps::Int         = 100,
        silent::Bool        = false,
        check_sign_problem::Bool = true,
        check_propagation_error::Bool = true,
        safe_mult::Int      = 10,
        measure_rate::Int   = 10,
        warn_round::Bool    = true,
        print_rate::Int     = max(10, div(thermalization + sweeps, 100)),
        kwargs...
    )
    nt = (;kwargs...)
    keys(nt) == (:beta,) && (nt = (;beta=nt.beta, delta_tau=0.1))
    @assert length(nt) >= 2 "Invalid keyword arguments to DQMCParameters: $nt"
    if (Set ∘ keys)(nt) == Set([:delta_tau, :beta, :slices])
        delta_tau, beta = nt.delta_tau, nt.beta
        slices = round(Int, beta/delta_tau)
        if slices != nt.slices
            error(
                "Given slices ($(nt.slices)) does not match calculated slices" *
                " beta/delta_tau ≈ $(slices)"
            )
        end
    elseif (Set ∘ keys)(nt) == Set([:beta, :slices])
        beta, slices = nt.beta, nt.slices
        delta_tau = beta / slices
    elseif (Set ∘ keys)(nt) == Set([:delta_tau, :slices])
        delta_tau, slices = nt.delta_tau, nt.slices
        beta = delta_tau * slices
    elseif (Set ∘ keys)(nt) == Set([:delta_tau, :beta])
        delta_tau, beta = nt.delta_tau, nt.beta
        slices = round(beta/delta_tau)
        if warn_round && !(slices ≈ beta/delta_tau)
            @warn "beta/delta_tau = $(beta/delta_tau) not an integer. Rounded to $slices"
        end
    else
        error("Invalid keyword arguments to DQMCParameters $nt")
    end

    DQMCParameters(
        thermalization,
        sweeps,
        silent,
        check_sign_problem,
        check_propagation_error,
        safe_mult,
        delta_tau,
        beta,
        slices,
        measure_rate,
        print_rate
    )
end


text:
using KCenters, KNearestCenters, TextSearch, TextClassification, SearchModels, UnicodePlots, StatsBase, Statistics
using LIBLINEAR, JLD2

function score_by_name(score, y, ypred)
    if score === :macrorecall
        recall_score(y, ypred, weight=:macro)
    elseif score === :macrof1
        f1_score(y, ypred, weight=:macro)
    elseif score === :accuracy
        accuracy_score(y, ypred)
    elseif score isa NamedTuple
        if score.name == :classf1
            classification_scores(y, ypred).classf1[score.label]
        elseif score.name == :classrecall
            classification_scores(y, ypred).classrecall[score.label]
        else
            error("Unknown score $score")
        end
    else
        error("Unknown score $score")
    end
end

function search_models_for_task(D, paramsfile, space, search_kwargs; verbose=true)
    cv = pop!(search_kwargs, :cv)

    if cv.kind == :folds
        partitions = stratified_folds(D.labels; folds=cv.folds, shuffle=true)
    elseif cv.kind == :montecarlo
        partitions = [stratified_split(D.labels; p=cv.ratio, shuffle=true) for i in 1:cv.repeat]
    elseif cv.kind == :holdout
        partitions = [stratified_split(D.labels; p=cv.ratio, shuffle=true)]
    else
        error("Unknown $cv")
    end

    score = pop!(search_kwargs, :score, :recall)
    initialpopulation = pop!(search_kwargs, :initialpopulation, 32)
    function error_function(config::MicroTC_Config)
        S = Float64[]
        for (itrain, itest) in partitions
            tc = MicroTC(config, D.corpus[itrain], D.labels[itrain]; verbose=true)
            ypred = predict_corpus(tc, D.corpus[itest])
            push!(S, score_by_name(score, D.labels[itest], ypred))
        end
        s = mean(S)
        verbose && println(stderr, "score: $s, $(typeof(config)), config: $(JSON3.write(config))")
        (1.0 - s, S)
    end

    iter = 0

    function inspect_population(space, params, population)
        iter += 1
        println(stderr, "===== inspecting population iter: $iter, popsize: $(length(population))")
        filename = paramsfile * ".iter=$iter"
        println(stderr, "saving iteration to $filename")
        jldsave(filename; population)
        println(lineplot([first(last(p)) for p in population], title="Error", xlabel="configurations", ylabel="error"))
    end

    search_kwargs[:inspect_population] = inspect_population
    best_list = search_models(space, error_function, initialpopulation; by=first, search_kwargs...)

    for (i, b) in enumerate(best_list)
        @info "-- perf best_lists[$i]:", b[1] => b[2]
    end

    best_list
end

function run_params(D, paramsfile, space, search_kwargs)
    if !isfile(paramsfile)
        population = search_models_for_task(D, paramsfile, space, search_kwargs)
        jldsave(paramsfile; population)
        population
    else
        load(paramsfile, "population")
    end
end

function run_train(D, config, modelfile; tok=Tokenizer(config.textconfig, invmap=nothing))
    if !isfile(modelfile)
        model = MicroTC(config, D.corpus, D.labels; tok)
        jldsave(modelfile; model)
        model
    else
        load(modelfile, "model")
    end
end

function run_predict(T, model, predictedfile; textkey="text", labelkey="klass", levels_=levels(T.labels))
    ypred = predict_corpus(model, T.corpus)
    open(predictedfile, "w") do f
        for (message, ilabel) in zip(T.corpus, ypred)
            println(f, JSON3.write(Dict(textkey => message, labelkey => levels_[ilabel])))
        end
    end

    ypred
end

max_stars_repo_name: sadit/text-classification-scripts

max_issues_repo_name: sadit/text-classification-scripts

max_forks_repo_name: sadit/text-classification-scripts


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'dta'_dashtextannotate

"""
    'dta'_dashtextannotate(;kwargs...)

A DashTextAnnotate component.
Learnings:
Don't try to print (<p>..</p>) lists or dictionaries directly
Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `text` (String; required): The text that should be annotated
- `entities` (Array; optional): The entities that are currently available
- `tag` (String; required): The currently selected tag (e.g. DATE) with which the selections will be tagged. This must be managed by the Dash App.
- `tag_colors` (Dict; required): The currently selected tag color (e.g. GREEN) with which the selected tag will be colored. This must be managed by the Dash App.
"""
function 'dta'_dashtextannotate(; kwargs...)
        available_props = Symbol[:id, :text, :entities, :tag, :tag_colors]
        wild_props = Symbol[]
        return Component("'dta'_dashtextannotate", "DashTextAnnotate", "dash_text_annotate", available_props, wild_props; kwargs...)
end


text:
function each_line_matching(callback, filepath, regexp)
  f = open(filepath, "r")
  i = 0
  for(line in eachline(f))
    i += 1
    if ismatch(regexp, line)
      callback(i, line, line[length(match(regexp, line).match)+1:end])
    end
  end
  close(f)
end

function levenshtein(s::ASCIIString, t::ASCIIString)
  m = length(s)
  n = length(t)
  if n == 0
    return(m)
  elseif m == 0
    return(m)
  end

  d = Array(Int64, m+1, n+1)

  for(i in 1:(m+1))
    d[i, 1] = i
  end
  for(j in 1:(n+1))
    d[1, j] = j
  end
  for(j in 2:n)
    for(i in 2:m)
      if s[i-1] == t[j-1]
        d[i, j] = d[i-1, j-1]
      else
        d[i, j] = minimum([d[i-1,j]+1, d[i, j-1]+1, d[i-1, j-1]+1])
      end
    end
  end
  d[m,n]
end

# Find lines in file
function find_macro_invocations_similar_to(file, macroname, exstr)
  regex = Regex("^\s*@" * macroname)

  # For each line which invokes the macro take what comes after the macroname
  # and skip spaces etc. Then check which line is closest to exstr.
  lines_with_distance = Any[]
  cb(linenum, line, lineaftermacro) = begin
    lam = replace(lineaftermacro, " ", "")
    push!(lines_with_distance, (line, levenshtein(lam, exstr), linenum))
  end
  each_line_matching(cb, file, regex)

  # Now hope we have the right one... It might not be if there are several same/similar.
  sort(lines_with_distance, by=(t)->t[2])[1]
end

macro mytestmacro(x)
  if !isa(x.args[1], Expr) || x.args[1].head != :line
    file = Base.source_path()
    # Ok, we need to disambiguate..
    line, dist, linenum = find_macro_invocations_similar_to(file, "mytestmacro", string(x)[3:end-1])
    println("line: $(linenum) file: $(file)")
  else
    println("line: $(x.args[1].args[1]) file: $(x.args[1].args[2])")
  end
end

@mytestmacro begin mytestfun() end
@mytestmacro 1==2
@mytestmacro 3<2


text:
# \sum_{i, j} A_{i, j} min(w_i, w_j)
function minplus_quad(A::AbstractMatrix{T},
                      w::AbstractVector{T}) where T
    length(w) == size(A, 1) == size(A, 2) ||
        throw(DimensionMismatch("A $(size(A)) and w ($(length(w))) size mismatch"))
    isempty(w) && return zero(T)
    res = zero(T)
    @inbounds for i in eachindex(w)
        wi = w[i]
        ioffset = (i-1)*size(A, 1)
        resi = zero(T)
        if zero(T) < wi < one(T)
            for j in eachindex(w)
                resi += A[ioffset+j]*min(wi, w[j])
            end
        elseif wi == one(T) # wi >= maxw
            for j in eachindex(w)
                resi += A[ioffset+j] * w[j] # FIXME dotprod() faster?
            end
        end
        res += resi
    end
    return res
end

take2(_, v) = v

# calculate res_i = fold(res_i, \sum_{j = 1} A_{i, j} min(u_i, v_j)), i = 1..N
function minplus_bilinear!(foldl::Function,
                           res::AbstractVector{T},
                           A::AbstractMatrix{T},
                           u::AbstractVector{T},
                           v::AbstractVector{T}) where T
    length(res) == length(u) == size(A, 2) ||
        throw(DimensionMismatch("res ($(length(res))), u ($(length(u))) or A $(size(A)) size mismatch"))
    length(v) == size(A, 1) ||
        throw(DimensionMismatch("v ($(length(v))) and A $(size(A)) size mismatch"))
    n = length(res)
    @inbounds for i in axes(A, 2)
        ui = u[i]
        ioffset = (i-1)*size(A, 1)
        x = zero(T)
        if zero(T) < ui < one(T) # if ui == 0, nothing is done
            for j in eachindex(v)
                x += A[ioffset+j] * min(ui, v[j])
            end
        elseif ui == one(T)
            for j in eachindex(v)
                x += A[ioffset+j] * v[j] # FIXME dotprod() faster?
            end
        end
        res[i] = foldl(res[i], x)
    end
    return res
end

function maxplus_linear(v::AbstractVector, A::SparseMatrixCSC)
    length(v) == size(A, 1) ||
        throw(DimensionMismatch("v ($(length(v))) and A $(size(A)) size mismatch"))
    Arows = rowvals(A)
    Avals = nonzeros(A)
    res = 0.0
    @inbounds for i in axes(A, 2)
        wmax = 0.0
        for j in nzrange(A, i)
            jr = Arows[j]
            wmax = max(wmax, v[jr] * Avals[j])
        end
        res += wmax
    end
    return res
end

function maxplus_linear!(res::AbstractVector, A::SparseMatrixCSC, u::AbstractVector)
    length(res) == size(A, 1) ||
        throw(DimensionMismatch("result ($(length(res))) and A $(size(A)) size mismatch"))
    length(u) == size(A, 2) ||
        throw(DimensionMismatch("u ($(length(u))) and A $(size(A)) size mismatch"))
    Arows = rowvals(A)
    Avals = nonzeros(A)
    fill!(res, 0.0)
    @inbounds for i in eachindex(u)
        ((ui = u[i]) > 0) || continue
        for j in nzrange(A, i)
            jr = Arows[j]
            res[jr] = max(res[jr], ui * Avals[j])
        end
    end
    return res
end


text:
using GPUAcceleratedTracking, DrWatson, ProgressMeter
@quickactivate "GPUAcceleratedTracking"

allparams = Dict(
    "num_samples" => 2 .^ (11:18),
    "algorithm" => ["gmem", "textmem"]
)

dicts = dict_list(allparams)

@showprogress 0.5 "Benchmarking reduction algorithms (pure vs cplx vs cplx_multi)" for (_, d) in enumerate(dicts)
    benchmark_results = run_replica_benchmark(d)
    @tagsave(
        datadir("benchmarks/codereplica", savename("CodeReplica", d, "jld2")),
        benchmark_results
    )
end

max_stars_repo_name: mfkiwl/GPUAcceleratedTracking

max_issues_repo_name: mfkiwl/GPUAcceleratedTracking

max_forks_repo_name: mfkiwl/GPUAcceleratedTracking


text:
module VideoIO

using Libdl                                             #0.227492 seconds
using Requires, Dates, ProgressMeter                    #0.001749 seconds
using ImageCore: permuteddimsview, channelview, rawview #0.723065 seconds
using ColorTypes: RGB, Gray, N0f8, YCbCr                #0.529263 seconds
using ImageTransformations: restrict                    #3.156594 seconds!!

include("init.jl")
include("util.jl")
include(joinpath(av_load_path, "AVUtil", "src", "AVUtil.jl"))
include(joinpath(av_load_path, "AVCodecs", "src", "AVCodecs.jl"))
include(joinpath(av_load_path, "AVFormat", "src", "AVFormat.jl"))
include(joinpath(av_load_path, "AVDevice", "src", "AVDevice.jl"))
include(joinpath(av_load_path, "SWScale", "src", "SWScale.jl"))

using .AVUtil
using .AVCodecs
using .AVFormat
using .SWScale

if FFMPEG.have_avdevice()
    import .AVDevice
end

include("info.jl")
include("avdictionary.jl")
include("avio.jl")
include("encoding.jl")
include("testvideos.jl")
using .TestVideos

if Sys.islinux()
    import Glob
    function init_camera_devices()
        append!(CAMERA_DEVICES, Glob.glob("video*", "/dev"))
        DEFAULT_CAMERA_FORMAT[] = AVFormat.av_find_input_format("video4linux2")
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        DEFAULT_CAMERA_DEVICE[] = isempty(CAMERA_DEVICES) ? "" : CAMERA_DEVICES[1]
    end
end

if Sys.iswindows()
    function init_camera_devices()
        append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "dshow", "dummy"))
        DEFAULT_CAMERA_FORMAT[] = AVFormat.av_find_input_format("dshow")
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        DEFAULT_CAMERA_DEVICE[] = string(
            "video=",
            isempty(CAMERA_DEVICES) ? "0" : CAMERA_DEVICES[1]
        )
    end
end

if Sys.isapple()
    function init_camera_devices()
        try
            append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "avfoundation", "\"\""))
            DEFAULT_CAMERA_FORMAT[] = AVFormat.av_find_input_format("avfoundation")
        catch
            try
                append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "qtkit", "\"\""))
                DEFAULT_CAMERA_FORMAT[] = AVFormat.av_find_input_format("qtkit")
            catch
            end
        end
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        # Note: "Integrated" is another possible default value
        DEFAULT_CAMERA_OPTIONS["pixel_format"] = "uyvy422"
        DEFAULT_CAMERA_DEVICE[] = isempty(CAMERA_DEVICES) ? "0" : CAMERA_DEVICES[1]
    end
elseif Sys.isbsd()
    # copied loosely from apple above - needs figuring out
    function init_camera_devices()
        append!(CAMERA_DEVICES, get_camera_devices(ffmpeg, "avfoundation", "\"\""))
        DEFAULT_CAMERA_FORMAT[] = AVFormat.av_find_input_format("avfoundation")
    end
    function init_camera_settings()
        DEFAULT_CAMERA_OPTIONS["framerate"] = 30
        DEFAULT_CAMERA_OPTIONS["pixel_format"] = "uyvy422"
        DEFAULT_CAMERA_DEVICE[] = isempty(CAMERA_DEVICES) ? "0" : CAMERA_DEVICES[1]
    end
end

#Helper functions to explain about Makie load order requirement
function play(f; flipx=false, flipy=false)
    error("Makie must be loaded before VideoIO to provide video playback functionality. Try a new session with `using Makie, VideoIO`")
end
function playvideo(video;flipx=false,flipy=false)
    error("Makie must be loaded before VideoIO to provide video playback functionality. Try a new session with `using Makie, VideoIO`")
end
function viewcam(device=DEFAULT_CAMERA_DEVICE, format=DEFAULT_CAMERA_FORMAT)
    if FFMPEG.have_avdevice()
        error("Makie must be loaded before VideoIO to provide camera playback functionality. Try a new session with `using Makie, VideoIO`")
    else
        error("No AV device detected")
    end
end

function __init__()
    # Always check your dependencies from `deps.jl`
    # TODO remove uncessary ENV["LD_LIBRARY_PATH"] from check_deps, so that
    # it doesn't mess with LD_LIBRARY_PATH, which was causing CI download issues due to issues with julia's curl
    # since check_deps is optional, I hope this is ok for now

    #check_deps()

    loglevel!(AVUtil.AV_LOG_FATAL)
    # @info "VideoIO: Low-level FFMPEG reporting set to minimal (AV_LOG_FATAL). See `? VideoIO.loglevel!` for options"

    read_packet[] = @cfunction(_read_packet, Cint, (Ptr{AVInput}, Ptr{UInt8}, Cint))

    av_register_all()

    if FFMPEG.have_avdevice()
        AVDevice.avdevice_register_all()
        init_camera_devices()
        init_camera_settings()
    end

    @require Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" begin
        # Define read and retrieve for Images
        function play(f; flipx=false, flipy=false, pixelaspectratio=nothing)
            if pixelaspectratio ≡ nothing # if user did not specify the aspect ratio we'll try to use the one stored in the video file
                pixelaspectratio = aspect_ratio(f)
            end
            h = f.height
            w = round(typeof(h), f.width*pixelaspectratio) # has to be an integer
            scene = Makie.Scene(resolution = (w, h))
            buf = read(f)
            makieimg = Makie.image!(scene, 1:h, 1:w, buf, show_axis = false, scale_plot = false)[end]
            Makie.rotate!(scene, -0.5pi)
            if flipx && flipy
                Makie.scale!(scene, -1, -1, 1)
            else
                flipx && Makie.scale!(scene, -1, 1, 1)
                flipy && Makie.scale!(scene, 1, -1, 1)
            end
            display(scene)
            while !eof(f) && isopen(scene)
                read!(f, buf)
                makieimg[3] = buf
                sleep(1 / f.framerate)
            end

        end

        function playvideo(video;flipx=false,flipy=false,pixelaspectratio=nothing)
            f = VideoIO.openvideo(video)
            play(f,flipx=flipx,flipy=flipy,pixelaspectratio=pixelaspectratio)
        end

        if FFMPEG.have_avdevice()
            function viewcam(device=DEFAULT_CAMERA_DEVICE, format=DEFAULT_CAMERA_FORMAT, pixelaspectratio=nothing)
                init_camera_settings()
                camera = opencamera(device[], format[])
                play(camera, flipx=true, pixelaspectratio=pixelaspectratio)
                close(camera)
            end
        end
    end
end

end # VideoIO


text:
# Starting with the number 1 and moving to the right in a clockwise direction a
# 5 by 5 spiral is formed as follows:
#
# 21 22 23 24 25 26
# 20  7  8  9 10 27
# 19  6  1  2 11 28
# 18  5  4  3 12 29
# 17 16 15 14 13 30
#
# It can be verified that the sum of the numbers on the diagonals is 101.
#
# What is the sum of the numbers on the diagonals in a 1001 by 1001 spiral
# formed in the same way?

using ProjectEulerSolutions

# Create and sum the sequence of corner values, then increment.
function p028solution_increment(n::Integer=5)::Integer
    val = 1
    startindex = 3
    for i = 2:2:n
        val += sum(collect(startindex .+ (0:i:i*3)))
        startindex = startindex + 4*i + 2
    end
    return val
end

# Closed form solution by noticing that the number in the top right corner is
# n^2 and the other corners as: n^2-n+1, n^2-2n+2, and n^2-3n+3.  Summing gives
# 4n^2-6n+6
function p028solution_broadcast(n::Integer=5)::Integer
    nset = 3:2:n
    vals = 4 .* nset .^ 2 .- 6 .* nset .+ 6
    return sum(vals) + 1
end

# Even more closed form solution solving for sum of series, which is:
# 2n^3/3 + n^2/2 + 4n/3 - 5/2.  Somehow it is not faster than either previous
# solution.
function p028solution_closedform(n::Integer=5)::Integer
    return fld(4 * n^3 + 3 * n^2 + 8 * n - 9, 6)
end

p028 = Problems.Problem(Dict("Incremental" => p028solution_increment,
                             "Broadcast"   => p028solution_broadcast,
                             "Closed form" => p028solution_closedform))

Problems.benchmark(p028, 1001)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
"""
	readgpcpd(filein)
Read DAILY GPCP data

**Input**
* filein: daily (GPCP)[https://precip.gsfc.nasa.gov/gpcp_daily_comb.html] binary file

**Output**
* Array{Float63}(lon,lat,time) with unchanged units (mm/day). All negative values (flag=-99999) are set to NaN

**Example**
```
# Load data for whole month and select last day of the month
# should result in approx.: https://climatedataguide.ucar.edu/climate-data/gpcp-daily-global-precipitation-climatology-project
filein = "/test/input/gpcpd_data";
datagpcp = readgpcpd(filein);
data_june30 = datagpcp(:,:,end);
timegpcp = readgpcpd_time(filein);
lon,lat = readgpcpd_lonlat(filein);
head_string = readgpcpd_head(filein);
```
"""
function readgpcpd(filein::String)
	# get header info
	headout = readgpcpd_head(filein);
	lon,lat = readgpcpd_lonlat(headout,from_file=false);
	timeout = readgpcpd_time(headout,from_file=false);
	# declare output
	dataout = Array{Float32}(undef,length(lon)*length(lat)*length(timeout));
	open(filein,"r") do fid
		# remove header lines
		seek(fid,length(headout)-3);
		# read data
		read!(fid,dataout);
	end
	#dataout[dataout.<0.] == NaN;
	return reshape(convert(Vector{Float64},dataout),length(lon),length(lat),length(timeout))
end

"""
	readgpcpd_head(filein;headersize)
Read DAILY GPCP header line

**Input**
* filein: daily (GPCP)[https://precip.gsfc.nasa.gov/gpcp_daily_comb.html] binary file
* optional header size (1440 for daily data)

**Output**
* whole header string

**Example**
```
filein = "/test/input/gpcpd_data";
head_string = readgpcpd_head(filein);
```
"""
function readgpcpd_head(filein::String;headsize::Int=1440)
	headout = ' ';
	open(filein,"r") do fid
		for i in 1:headsize
			headout = string(headout,read(fid,Char)); # concatenate characters
		end
	end
	return headout[2:end] # first value = ' '
end

"""
	readgpcpd_lonlat(filein;from_file=true)
Read DAILY GPCP longitude and latitude

**Input**
* filein: daily (GPCP)[https://precip.gsfc.nasa.gov/gpcp_daily_comb.html] binary file
* read from file (true) or parse from string (false)

**Output**
* longitude (degree east), latitude (degree north to east)

**Example**
```
filein = "/test/input/gpcpd_data";
lon,lat = readgpcpd_lonlat(filein);
```
"""
function readgpcpd_lonlat(headout::String;from_file=true)
	if from_file==true
		headout = readgpcpd_head(headout);
	end
	string_size = match(r"[0-9]{3}x[0-9]{3}",headout);
	x,y = split(string_size.match,"x") |> x -> (Meta.parse(x[1]),Meta.parse(x[2]))
	return collect(0.5:1:(x-0.5)),collect(89.5:-1:(-y/2+0.5))
end

"""
	readgpcpd_time(filein)
Read DAILY GPCP time information

**Input**
* filein: daily (GPCP)[https://precip.gsfc.nasa.gov/gpcp_daily_comb.html] binary file
* read from file (true) or parse from string (false)

**Output**
* time vector in DateTime format

**Example**
```
filein = "/test/input/gpcpd_data";
timegpcp = readgpcpd_time(filein);
```
"""
function readgpcpd_time(headout::String;from_file=true)
	if from_file==true
		headout = readgpcpd_head(headout);
	end
	string_yyyy = match(r"year=[0-9]{4}",headout);
	string_mm = match(r"month=[0-9]{2}",headout);
	string_dd = match(r"days=.-[0-9]{2}",headout) |> x->split(x.match,"=") |> x->x[2]
	time_start = DateTime(Meta.parse(split(string_yyyy.match,"=")[2]),
						Meta.parse(split(string_mm.match,"=")[2]),
						Meta.parse(split(string_dd,"-")[1]));
	time_stop = DateTime(Meta.parse(split(string_yyyy.match,"=")[2]),
						Meta.parse(split(string_mm.match,"=")[2]),
						Meta.parse(split(string_dd,"-")[2]));
	return collect(time_start:Dates.Day(1):time_stop)
end

"""
	coorShift_lon(lon,datain,to)
Convert 0 to 360 longitude coordinate (e.g. GPCP) system to -180 to 180

**Input**
* lon: longitude vector
* datain: data matrix or 3D Array{T}. If 3D matrix => first 2 dimensions are lon*lat or lat*lon
* to: switch between systems (will be converted to this system): "-180to180" | "0to360"

**Output**
* lonout: new longitude vector
* dataout: transformed data in the same dimensions as input 'datain'

**Example**
```
# convert to 0-360 system
lon = collect(-179.5:1:179.5);
datain = repeat(-179.5:1:179.5,1,180);
lonout,dataout = coorShift_lon(lon,datain,to="0to360")

# convert to -180 to 180 system
lon = collect(0.125:0.25:359.875);
datain = Array{Float64}(720,length(lon),2);
datain[:,:,1] = transpose(repeat(0.125:0.25:359.875,1,720));
datain[:,:,2] = transpose(repeat(0.125:0.25:359.875,1,720));
lonout,dataout = coorShift_lon(lon,datain,to="-180to180");
```
"""
function coorShift_lon(lon,datain;to::String="-180to180")
	lonin = copy(lon);
	# prepare output longitude
	if to == "-180to180"
		lonout = sort(lon .- 180);
		lonin[lon.>180] = lon[lon.>180] .- 360.;
	elseif to == "0to360"
		lonout = sort(lon .+ 180);
		lonin[lon.<0] = lon[lon.<0] .+ 360.;
	end
	# find indices
	i_lon = zeros(Int,length(lon))
	for (i,v) in enumerate(lonout)
		i_lon[i] = findall(v.==lonin)[1];
	end
	# Get the correct dimension and transform the data
	if size(datain,1) == length(lon)
		if size(datain,3) == 1
			return lonout,datain[i_lon,:];
		else
			return lonout,datain[i_lon,:,:];
		end
	else
		if size(datain,3) == 1
			return lonout,datain[:,i_lon];
		else
			return lonout,datain[:,i_lon,:];
		end
	end
end


text:
# # Symmetry reduction

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Symmetry/symmetry_reduction.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Symmetry/symmetry_reduction.ipynb)
# **Adapted from**: [SymbolicWedderburn example](https://github.com/kalmarek/SymbolicWedderburn.jl/blob/tw/ex_sos/examples/ex_C4.jl)

import MutableArithmetics
const MA = MutableArithmetics
using MultivariatePolynomials
const MP = MultivariatePolynomials
using MultivariateBases
const MB = MultivariateBases

using Test #src
using DynamicPolynomials
@polyvar x[1:4]

# We would like to find the minimum value of the polynomial

poly = sum(x) + sum(x.^2)

# As we can decouple the problem for each `x[i]` for which `x[i] + x[i]^2` has
# minimum value 0.25, we would expect to get `-1` as answer.
# Can this decoupling be exploited by SumOfSquares as well ?
# For this, we need to use a certificate that can exploit the permutation symmetry of the polynomial.

using SumOfSquares

# We define the symmetry group as a permutation group in the variables.
# In order to do that, we define the action of a permutation on a monomial
# as the monomial obtained after permuting the variables.

using PermutationGroups
G = PermGroup([perm"(1,2,3,4)"])

# We can use this certificate as follows:

import CSDP
solver = CSDP.Optimizer
model = Model(solver)
@variable(model, t)
@objective(model, Max, t)
pattern = Symmetry.Pattern(G, Symmetry.VariablePermutation())
con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)
optimize!(model)
@test value(t) ≈ -1 rtol=1e-6 #src
value(t)

# We indeed find `-1`, let's verify that symmetry was exploited:

g = gram_matrix(con_ref).sub_gram_matrices    #src
@test length(g) == 4                          #src
@test length(g[1].basis.polynomials) == 2     #src
@test g[1].basis.polynomials[1] == 1.0        #src
@test g[1].basis.polynomials[2] ≈ -0.5 * sum(x) #src
@test size(g[1].Q) == (2, 2)                  #src
@test g[1].Q[1, 1] ≈ 1.0 atol=1e-6            #src
@test g[1].Q[1, 2] ≈ -1.0 atol=1e-6            #src
@test g[1].Q[2, 2] ≈ 1.0 atol=1e-6            #src
@test length(g[2].basis.polynomials) == 1     #src
@test g[2].basis.polynomials[1] ≈ (x[3] - x[1]) / √2 #src
@test size(g[2].Q) == (1, 1)                  #src
@test g[2].Q[1, 1] ≈ 1.0 atol=1e-6            #src
@test length(g[3].basis.polynomials) == 1     #src
@test g[3].basis.polynomials[1] ≈ (x[4] - x[2]) / √2 #src
@test size(g[3].Q) == (1, 1)                  #src
@test g[3].Q[1, 1] ≈ 1.0 atol=1e-6            #src
@test length(g[4].basis.polynomials) == 1     #src
@test g[4].basis.polynomials[1] ≈ (-x[1] + x[2] - x[3] + x[4]) / 2 #src
@test size(g[4].Q) == (1, 1)                  #src
@test g[4].Q[1, 1] ≈ 1.0 atol=1e-6            #src
for g in gram_matrix(con_ref).sub_gram_matrices
    println(g.basis.polynomials)
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
"""
    drawpreimageR2(f [; iterations])

Return the drawing of the preimage of the \$n\$-iterative of a function
    \$f:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$,
    in a rectangular region in \$\\mathbb{R}^2\$.

Using the a color map function \$c:\\mathbb{R}^2\\rightarrow Colors\$,
is drawn \$c\\circ f^n(x,y)\$ for \$(x,y)\$ in the rectangular region.

#### Arguments
- `f::Function`: A function \$f:\\mathbb{R}^2\\rightarrow\\mathbb{R}^2\$.
- `iterations::Integer`: Number of iterations.
"""
function drawpreimageR2(f::Function; iterations::Int=1)

    # Veryfying if graphics backend supports functions
    SDDGraphics.supported(:drawpixel)

    # Verifying functions
    @assert typeof(f(1., 1.)) <: Tuple{Real,Real}

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(iterations)

    #fn = @iterativeR2 f iterations
    fn = foriterativeR2(f,iterations)

    @sweeprectregion SDDGraphics.xlims() SDDGraphics.ylims() SDDGraphics.canvassize() begin
        SDDGraphics.color(fn(x,y)...)
        SDDGraphics.drawpixel(i,j)
    end # Implemented algorithm

    SDDGraphics.drawing()
end

"""
    drawpreimageC(f [; iterations])

Return the drawing of the preimage of the \$n\$-iterative of a function
    \$f:\\mathbb{C}\\rightarrow\\mathbb{C}\$,
    in a rectangular region in \$\\mathbb{C}\$.

Using the a color map function \$c:\\mathbb{C}\\rightarrow Colors\$,
is drawn \$c\\circ f^n(z)\$ for \$z\$ in the rectangulr region.

#### Arguments
- `f::Function`: A function \$f:\\mathbb{C}\\rightarrow\\mathbb{C}\$.
- `iterations::Integer`: Number of iterations.
"""
function drawpreimageC(f::Function; iterations::Int=1)

    # Veryfying if graphics backend supports functions
    SDDGraphics.supported(:drawpixel)

    # Verifying functions
    @assert typeof(f(1.0im)) <: Number

    SDDGraphics.newdrawing()
    SDDGraphics.updatecolorarray(iterations)

    #fn = @iterative f iterations
    fn = foriterative(f,iterations)

    @sweeprectregion SDDGraphics.xlims() SDDGraphics.ylims() SDDGraphics.canvassize() begin
        z = complex(x,y)
        SDDGraphics.color(fn(z))
        SDDGraphics.drawpixel(i,j)
    end # Implemented algorithm

    SDDGraphics.drawing()
end


text:
@sumProductRule(:node_type     => Exponential,
                :outbound_type => Message{LogNormal},
                :inbound_types => (Nothing, Message{Gaussian}),
                :name          => SPExponentialOutNG)

@sumProductRule(:node_type     => Exponential,
                :outbound_type => Message{GaussianMeanVariance},
                :inbound_types => (Message{LogNormal}, Nothing),
                :name          => SPExponentialIn1LN)

@sumProductRule(:node_type     => Exponential,
                :outbound_type => Message{PointMass},
                :inbound_types => (Nothing, Message{PointMass}),
                :name          => SPExponentialOutNP)

@sumProductRule(:node_type     => Exponential,
                :outbound_type => Message{PointMass},
                :inbound_types => (Message{PointMass}, Nothing),
                :name          => SPExponentialIn1PN)


text:
function ACQRS_ssi1(y,n,i)
	#
	# ssi algorithm for the identification of state space models
	# -- CVA version --
	#
	# Reference:
	# Subspace Identification for Linear Systems
	# Theory - Implementation - Applications
	# Peter Van Overschee / Bart De Moor
	# Kluwer Academic Publishers, 1996
	#
	# javier.cara@upm.es, 2015-08
	# --------------------------------------------------

	# Data by rows
	y,m,N = byrow(y)

	# Hankel matrix
	# --------------------------------------------------
	# number of Hankel matrix columns
	j = N - 2*i + 1
	if (j < m*2*i)
		# rows(H) has to be > columns(H)
		error("Not enough data for building the Hankel matrix")
	end
	Y = hankel_yij(y/sqrt(j),2*i,j)

	# LQ factorization
	# --------------------------------------------------
	F = qrfact(Y')
	L = F[:R]'
	L21 = L[m*i+1:2*m*i,1:m*i]

	# singular values
	# --------------------------------------------------
	F = svdfact(L21)
	U1 = F[:U][:,1:n]
	S1 = F[:S][1:n]

	# Matrices gam and gam1
	# --------------------------------------------------
	gam  = U1*diagm(sqrt(S1))
	gam1 = U1[1:m*(i-1),:]*diagm(sqrt(S1))
	# and pseudo-inverses
	gam_inv  = pinv(gam)
	gam1_inv = pinv(gam1)

	# Determine the states Xi and Xi1
	# --------------------------------------------------
	Xi  = gam_inv  * L21
	Xi1 = gam1_inv * L[m*(i+1)+1:2*m*i,1:m*(i+1)]

	# Computing the state matrices A and C
	# --------------------------------------------------
	Rhs = [ Xi zeros(n,m) ]	# Right hand side
	Lhs = [ Xi1; L[m*i+1:m*(i+1),1:m*(i+1)] ] # Left hand side

	# Least squares
	sol = Lhs/Rhs

	A = sol[1:n,1:n]
	C = sol[n+1:n+m,1:n]

	# Computing the covariance matrices Q, R and S
	# -------------------------------------------------
	# Residuals
	res = Lhs - sol*Rhs
	cov = res*res'

	Q = cov[1:n,1:n]
	S = cov[1:n,n+1:n+m]
	R = cov[n+1:n+m,n+1:n+m]

	return A,C,Q,R,S

end

##########################################################
function hankel_yij(y,i,j)
	# Hankel matrix with i rows and j columns
	# from data [y_1 y_2 ... y_N],  y_t in R^{m x 1}

	y,m,N = byrow(y)

	# Hankel matrix
	H = zeros(m*i,j)
	for k in 1:i
		H[(k-1)*m+1:k*m,:] = y[:,k:k+j-1]
	end

	return H

end


text:
include("arithmetic.jl")
include("exponential.jl")
include("extrema.jl")
include("hyperbolic.jl")
include("inverse_hyperbolic.jl")
include("inverse_trig.jl")
include("trig.jl")


text:
export iauPap
"""
Position-angle from two p-vectors.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  vector/matrix support function.

Given:
   a      double[3]  direction of reference point
   b      double[3]  direction of point whose PA is required

Returned (function value):
         double     position angle of b with respect to a (radians)

Notes:

   1. The result is the position angle, in radians, of direction b with
   respect to direction a.  It is in the range -pi to +pi.  The
   sense is such that if b is a small distance "north" of a the
   position angle is approximately zero, and if b is a small
   distance "east" of a the position angle is approximately +pi/2.

   2. The vectors a and b need not be of unit length.

   3. Zero is returned if the two directions are the same or if either
   vector is null.

   4. If vector a is at a pole, the result is ill-defined.

Called:
   iauPn        decompose p-vector into modulus and direction
   iauPm        modulus of p-vector
   iauPxp       vector product of two p-vectors
   iauPmp       p-vector minus p-vector
   iauPdp       scalar product of two p-vectors

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# double iauPap(double a[3], double b[3])
function iauPap(a::AbstractVector{<:Real}, b::AbstractVector{<:Real})
   return ccall((:iauPap, libsofa_c), Cdouble,
       (Ptr{Cdouble}, Ptr{Cdouble}),
       convert(Array{Float64, 1}, a),
       convert(Array{Float64, 1}, b))
end


text:
# Tools for dealing with nested dictionaries
"""
    nested_default_dict_get(dict,k,keylist...;default=nothing)

Query into an arbitrarily deep nested dict until finding a non-dict value or
returning `default` when the current key does not point to a dict.
Example
@julia```
dict = Dict()
nested_default_dict_set!(dict,1,2,"val")
nested_default_dict_set!(dict,1,2,3,"other val")
nested_default_dict_get(dict,1,2) == "val"
nested_default_dict_get(dict,1,2,3) == "val"
```
"""
function nested_default_dict_get(dict,k,keylist...;default=nothing,_default_key=:Default)
    if haskey(dict,k)
        val = dict[k]
        if isa(val,Dict)
            return nested_default_dict_get(val,keylist...;default=default,_default_key=_default_key)
        else
            # return value early
            return val
        end
    end
    if !(default === nothing)
        return default
    else
        return get(dict,_default_key,default)
    end
    return get(dict,_default_key,default)
end
function nested_default_dict_get(dict,k;default=nothing,_default_key=:Default)
    if haskey(dict,k)
        val = dict[k]
        if isa(val,Dict)
            if haskey(val,_default_key)
                return val[_default_key]
            end
        else
            return val
        end
    end
    if !(default === nothing)
        return default
    else
        return get(dict,_default_key,default)
    end
    return get(dict,_default_key,default)
end
"""
    nested_default_dict_set!(dict,keylist...,[val])

Sets the value of a nested dictionary, indexed by a list of keys, to `val`. If
the full `keylist` does not yet exist, create the intermediate dicts along the
way.
Example:
@julia```
dict = Dict()
nested_default_dict_set!(dict,1,2,3,"my val")
```
"""
function nested_default_dict_set!(dict,k,k2,args...;_default_key=:Default)
    subdict = get!(dict,k,Dict())
    if !isa(subdict,Dict)
        subdict = Dict{Any,Any}(_default_key=>subdict)
        dict[k] = subdict
    end
    nested_default_dict_set!(subdict,k2,args...;_default_key=_default_key)
end
function nested_default_dict_set!(dict,k,val;_default_key=:Default)
    if haskey(dict,k)
        old_val = dict[k]
        if isa(old_val,Dict)
            old_val[_default_key] = val
            return val
        end
    end
    dict[k] = val
end


text:
# module Bukdu

abstract ApplicationController

# verbs: get, post, delete, patch, put
const HTTP_VERBS = [:get, :post, :delete, :patch, :put]

import Base: get, show, edit, getindex

for verb in HTTP_VERBS
    @eval $verb{AC<:ApplicationController}(path::String, ::Type{AC}, action::Function; kw...) =
        Routing.match($verb, path, AC, action, Dict(kw))
end

function getindex{AC<:ApplicationController}(C::AC, sym::Symbol)
    if sym in [:query_params, :params, :host, :action, :private, :assigns]
        task = current_task()
        if haskey(Routing.task_storage, task)
            branch = Routing.task_storage[task]
            return getfield(branch, sym)
        else
            throw(ErrorException("no $task"))
        end
    end
    throw(KeyError(sym))
end

# actions: index, edit, new, show, create, update, delete
function index
end

function edit
end

function new
end

function show
end

function create
end

function update
end

function delete
end


text:
using NetMSA
using Documenter

makedocs(;
    modules=[NetMSA],
    authors="Aadam <aadimator@gmail.com> and contributors",
    repo="https://github.com/aadimator/NetMSA.jl/blob/{commit}{path}#L{line}",
    sitename="NetMSA.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://aadimator.github.io/NetMSA.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/aadimator/NetMSA.jl",
)


text:
# ---
# title: 1455. Check If a Word Occurs As a Prefix of Any Word in a Sentence
# id: problem1455
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/check-if-a-word-occurs-as-a-prefix-of-any-word-in-a-sentence/description/>
# hidden: true
# ---
#
# Given a `sentence` that consists of some words separated by a  **single
# space** , and a `searchWord`.
#
# You have to check if `searchWord` is a prefix of any word in `sentence`.
#
# Return _the index of the word_ in `sentence` where `searchWord` is a prefix of
# this word ( **1-indexed** ).
#
# If `searchWord` is a prefix of more than one word, return the index of the
# first word **(minimum index)**. If there is no such word return **-1**.
#
# A **prefix** of a string `S` is any leading contiguous substring of `S`.
#
#
#
# **Example 1:**
#
#
#
#     Input: sentence = "i love eating burger", searchWord = "burg"
#     Output: 4
#     Explanation: "burg" is prefix of "burger" which is the 4th word in the sentence.
#
#
# **Example 2:**
#
#
#
#     Input: sentence = "this problem is an easy problem", searchWord = "pro"
#     Output: 2
#     Explanation: "pro" is prefix of "problem" which is the 2nd and the 6th word in the sentence, but we return 2 as it's the minimal index.
#
#
# **Example 3:**
#
#
#
#     Input: sentence = "i am tired", searchWord = "you"
#     Output: -1
#     Explanation: "you" is not a prefix of any word in the sentence.
#
#
# **Example 4:**
#
#
#
#     Input: sentence = "i use triple pillow", searchWord = "pill"
#     Output: 4
#
#
# **Example 5:**
#
#
#
#     Input: sentence = "hello from the other side", searchWord = "they"
#     Output: -1
#
#
#
#
# **Constraints:**
#
#   * `1 <= sentence.length <= 100`
#   * `1 <= searchWord.length <= 10`
#   * `sentence` consists of lowercase English letters and spaces.
#   * `searchWord` consists of lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end