text:
module Direct

using Winston

using DtPlot
using DtProtos.Pdfs

using DtProtos

function buildA(f, fo)
    info("Building A")

end

"""
maximise ... blah

"""
function maximise()
    # blah
    info("Maximise")
    niter = -1
    nmax = 20
    while niter < nmax
        n += 1

    end
end

#-----------------------------
# From here down move to separate file

function create_candidate(mu, sigma, weight)
    n = GaussianPdf(mu, sigma, weight)
end

"""
Create historical forecasts.

"""

function create_historic_obs(d1, dn)
    # 0,100
    dates = [d1:dn]
    obs = sin(dates)
    return obs
end

function create_historic_forecasts(obs::Array{Float64},
                                   errmu::Float64, errsig::Float64)
    fs = SimplePdf[]
    for o in obs
        mu = o + errmu*(rand() - 0.5)
        sig = errsig*rand()
        p = GaussianPdf(mu, sig, 1)
        push!(fs, p) # list of simplePdfs
    end
    return fs
end

function expected_score(obs::Array{Float64}, fs::Array{SimplePdf}, score)
    #as = Array{Float64}
    as = Float64[]
    n = length(obs)
    for i in [1:n]
        o = obs[i]
        f = fs[i]
        a = score(f,o)
        #println(a)
        push!(as, a)
    end
    A = mean(as)
end

function ign(p::SimplePdf, o::Float64)
    if isinf(-log(pdf(p,o)))
        println(o, ", ", p.mu, ", ", p.sigma, ", ", pdf(p,o))
    end
    return -log(pdf(p,o))
end

function buildF(fo, ao, lambda, mu) # B1T, B2
    info("Building F")
    F = zeros(nX)
    for i = [1:1:nX]
        g = (1/tau) * fo[i] .* log(fo[i])
        # for j = [1:1:ncAND]
        #     B1Tyk = a .*
        # end
    end

end

function buildA()
    info("Building A")
    A = speye(nX)
    A = A ./ tau
end

function solve(A, F)
    info("Solve")
    blas_set_num_threads(16)
    #v = Sparsesolve(A,F)
    f = \(A, F)
    return f

end

function updateLambda(lambda, f, a)
    for i in [1:length(lambda)]
        lambda[i] += sum(f[i] * a[1])
    end
    return lambda
end

function updateMu(mu, f)
    intf = sum(f^2) #* dx
    mu += omega * (intf - 1)
    return mu
end

function run_direct()
    # Setup data (historic)
    #
    obs = create_historic_obs(0, 20)

    ncandidates = 2
    c1 = create_historic_forecasts(obs, 0.25, 2.)
    c2 = create_historic_forecasts(obs, 0.50, 1.)

    # do leaf plot here
    # date; obs, c1, c2
    fp = FramedPlot()
    fp = leaf2(fp, c1, "red")
    fp = leaf2(fp, c2, "blue")

    obsp = Points([1:length(obs)], obs, color="black")
    style(obsp, kind="circle", size=0.3)
    obsc = Curve([1:length(obs)], obs, color="grey")
    add(fp, obsp)
    add(fp, obsc)

    A1 = expected_score(obs, c1, ign)
    A1p = Slope(0, (0,A1), kind="dotted", color="red")
    setattr(A1p, label="A1")
    add(fp, A1p)
    println("------",A1)

    A2 = expected_score(obs, c2, ign)
    A2p = Slope(0, (0,A2), kind="dotted", color="blue")
    setattr(A2p, label="A2")
    add(fp, A2p)
    println("------",A2)

    l = Legend(.1, .9, {A1p, A2p})
    add(fp,l)

    display(fp)

    # Setup numerics
    #
    # X: observation domain (and space over which integration done)
    dobs = maximum(obs) - minimum(obs)
    obsrange = [minimum(obs)-dobs : maximum(obs)+dobs] # ie X
    X = obsrange
    nX = length(X)

    # # Main loop
    # #
    # f1 = zeros(X)
    # f0 = f1
    # mu = 0.
    # lambda = zeros(ncadidates)
    # tau = 1.0
    # ncount = -1
    # while ncount < maxcount
    #     ncount += 1
    #     F  = buildF(f, fo, mu)
    #     A  = buildA(f, fo)
    #     fo = f
    #     f  = solve(A, F)

    #     f[1] = f[2]
    #     f[end-1] = f[end-2]

    #     mu = updateMu(mu, f)
    #     info("mu = ", mu)

    #     lambda = updateLambda(lambda, f)
    #     info("lambda = ", lambda)
    # end
end

end # module


text:
# Test various interconnections using MATLAB

using MATLAB, LinearAlgebra
include("StateSpace.jl");
include("Interconnections.jl");

# Test Parallel(...)
println("\nTesting Parallel(...)");
for i=1:10
    A1 = rand(3,3);
    B1 = rand(3,3);
    C1 = rand(3,3);
    D1 = rand(3,3);

    A2 = rand(3,3);
    B2 = rand(3,3);
    C2 = rand(3,3);
    D2 = rand(3,3);

    s1 = StateSpace(A1,B1,C1,D1);
    s2 = StateSpace(A2,B2,C2,D2);
    s3 = Parallel(s1,s2,Vector(1:2),Vector(1:2),Vector(1:3),Vector(1:3));

    # Get MATLAB's implementation.
    mat"s1 = ss($A1,$B1,$C1,$D1)";
    mat"s2 = ss($A2,$B2,$C2,$D2)";
    mat"s3 = parallel(s1,s2,[1:2],[1:2],[1:3],[1:3])";

    mat"Ap=s3.A; Bp=s3.B;Cp=s3.C;Dp=s3.D"
    @mget Ap Bp Cp Dp

    eA = norm(s3.A-Ap);
    eB = norm(s3.B-Bp)
    eC = norm(s3.C-Cp)
    eD = norm(s3.D-Dp)

    println("Norm of error: ", norm([eA,eB,eC,eD]));
end

# Test Feedback(...)
println("\n\nTesting Feedback(...)");
for i=1:10
    A1 = rand(3,3);
    B1 = rand(3,3);
    C1 = rand(3,3);
    D1 = rand(3,3);

    A2 = rand(3,3);
    B2 = rand(3,2);
    C2 = rand(2,3);
    D2 = rand(2,2);

    s1 = StateSpace(A1,B1,C1,D1);
    s2 = StateSpace(A2,B2,C2,D2);
    s3 = Feedback(s1,s2,Vector(2:3),Vector(2:3));

    # Get MATLAB's implementation.
    mat"s1 = ss($A1,$B1,$C1,$D1)";
    mat"s2 = ss($A2,$B2,$C2,$D2)";
    mat"s3 = feedback(s1,s2,[2:3],[2:3])";
    mat"Ap=s3.A; Bp=s3.B;Cp=s3.C;Dp=s3.D"
    @mget Ap Bp Cp Dp

    eA = norm(s3.A-Ap);
    eB = norm(s3.B-Bp)
    eC = norm(s3.C-Cp)
    eD = norm(s3.D-Dp)

    println("Norm of error: ", norm([eA,eB,eC,eD]));
end

# Test Series(...)
println("\n\nTesting Series(...)");
for i=1:10
    A1 = rand(3,3);
    B1 = rand(3,3);
    C1 = rand(3,3);
    D1 = rand(3,3);

    A2 = rand(3,3);
    B2 = rand(3,3);
    C2 = rand(2,3);
    D2 = rand(2,3);

    s1 = StateSpace(A1,B1,C1,D1);
    s2 = StateSpace(A2,B2,C2,D2);
    s3 = Series(s1,s2,Vector(2:3),Vector(1:2));
    s33 = Subsystem(s3,Vector(2:4),Vector(1:2)); # Matlab returns this.

    # Get MATLAB's implementation.
    mat"s1 = ss($A1,$B1,$C1,$D1)";
    mat"s2 = ss($A2,$B2,$C2,$D2)";
    mat"s3 = series(s1,s2,[2:3],[1:2])";
    mat"Ap=s3.A; Bp=s3.B;Cp=s3.C;Dp=s3.D"
    @mget Ap Bp Cp Dp

    # The ordering of states is different in Matlab.
    # We check sorted elements of matrices
    eA = sort(s33.A[:])-sort(Ap[:]);
    eB = sort(s33.B[:])-sort(Bp[:]);
    eC = sort(s33.C[:])-sort(Cp[:]);
    eD = sort(s33.D[:])-sort(Dp[:]);

    println("Norm of error: ",norm([eA;eB;eC;eD]))
end

ext: jl

max_stars_repo_name: isrlab/RobustControl

max_issues_repo_name: isrlab/RobustControl

max_forks_repo_name: isrlab/RobustControl


text:
# function mpi_collapse_solu()# (comms_ctx, horizontal_mesh, quad, z_max, z_elem, z_stretch, Y) # TODO: it's using the global scope here
if is_distributed # replace sol.u on the root processor with the global sol.u

    if ClimaComms.iamroot(comms_ctx)
        global_h_space = make_horizontal_space(horizontal_mesh, quad, nothing)
        global_center_space, global_face_space = make_hybrid_spaces(global_h_space, z_max, z_elem, z_stretch)
        global_Y_c_type = Fields.Field{typeof(Fields.field_values(Y.c)), typeof(global_center_space)}
        global_Y_f_type = Fields.Field{typeof(Fields.field_values(Y.f)), typeof(global_face_space)}
        global_Y_type = Fields.FieldVector{FT, NamedTuple{(:c, :f), Tuple{global_Y_c_type, global_Y_f_type}}}
        global_Y_slab_type =
            Fields.Field{typeof(Fields.field_values(slab_sim.integrator.u.T_sfc)), typeof(global_h_space)}
        global_sol_u_atmos = similar(sol_atm.u, global_Y_type)
        global_sol_u_slab = similar(sol_slab.u, global_Y_slab_type)
    end
    for i in 1:length(sol_atm.u)
        global_Y_c = DataLayouts.gather(comms_ctx, Fields.field_values(sol_atm.u[i].c))
        global_Y_f = DataLayouts.gather(comms_ctx, Fields.field_values(sol_atm.u[i].f))
        global_Y_slab = DataLayouts.gather(comms_ctx, Fields.field_values(sol_slab.u[i].T_sfc))
        if ClimaComms.iamroot(comms_ctx)
            global_sol_u_atmos[i] = Fields.FieldVector(
                c = Fields.Field(global_Y_c, global_center_space),
                f = Fields.Field(global_Y_f, global_face_space),
            )
            global_sol_u_slab[i] = Fields.Field(global_Y_slab, global_h_space)
        end
    end
    if ClimaComms.iamroot(comms_ctx)
        solu_atm = global_sol_u_atmos
        solu_slab = global_sol_u_slab
        return solu_atm, solu_slab
    end
else
    solu_atm = sol_atm.u
    h_space = make_horizontal_space(horizontal_mesh, quad, nothing) #TODO move this to the beginning (once same the instance error sorted)
    solu_slab = Fields.FieldVector(T_sfc = [Fields.Field(Fields.field_values(u.T_sfc), h_space) for u in sol_slab.u])
end

max_stars_repo_name: CliMA/ClimaCoupler


text:
export init_random_pars!

"""
    build_rng_generator_T(T::abstractArray, seed)

builds an RNG generator for the type of T.
    Defaults to MersenneTWister.
"""
function build_rng_generator_T(arrT::Array, seed)
    return MersenneTwister(seed)
end

"""
    init_random_pars!([rng=GLOBAL_RNG], net; sigma=0.01 )

Initializes all weights of the neural network to a random distribution with
variance `sigma`.
"""
init_random_pars!(net::NeuralNetwork, args...; kwargs...) = init_random_pars!(GLOBAL_RNG, net, args...; kwargs...)
function init_random_pars!(rng::AbstractRNG, net; sigma=0.01)
    for f=trainable(net)
        init_random_pars!(rng, f; sigma=sigma)
    end
    return net
end

function init_random_pars!(rng::AbstractRNG, f::AbstractArray; sigma=0.01)
    randn!(rng, f)
    f .*= sqrt(sigma)

    return f
end


text:
module YORUBA

import ..LOCATIONS
import ..LABELS
import ..NON_SPEECH_UNITS

using Dates
using NaturalSort

const URL_FEMALE_AUDIO = "https://www.openslr.org/resources/86/yo_ng_female.zip"
const URL_FEMALE_TRANS = "https://www.openslr.org/resources/86/line_index_female.tsv"
const URL_MALE_AUDIO = "https://www.openslr.org/resources/86/yo_ng_male.zip"
const URL_MALE_TRANS = "https://www.openslr.org/resources/86/line_index_male.tsv"
const URL_ALI = "https://raw.githubusercontent.com/beer-asr/beer/master/recipes/aud/local/google_lr/ali"

function prepare(datadir)

    ###################################################################
    # Download the data

    rawdata = mkpath(joinpath(datadir, LOCATIONS[:local], "data"))
    orig_wavs_f = mkpath(joinpath(rawdata, "orig_wavs_female"))
    orig_wavs_m = mkpath(joinpath(rawdata, "orig_wavs_male"))
    if ! ispath(joinpath(rawdata, ".done"))
        @info "Dowloading the data..."
        run(`wget -P $rawdata $URL_FEMALE_AUDIO`)
        run(`wget -P $rawdata $URL_FEMALE_TRANS`)
        run(`wget -P $rawdata $URL_MALE_AUDIO`)
        run(`wget -P $rawdata $URL_MALE_TRANS`)

        run(`unzip -d $orig_wavs_f $(joinpath(rawdata, "yo_ng_female.zip"))`)

        run(`unzip -d $orig_wavs_m $(joinpath(rawdata, "yo_ng_male.zip"))`)

        open(joinpath(rawdata, ".done"), "w") do f print(f, now()) end
    else
        @info "Data already dowloaded."
    end

    ###################################################################
    # Fix the timing and down-sample to 16kHz

    wavdir = mkpath(joinpath(datadir, LOCATIONS[:local], "wavs"))

    if ! ispath(wavdir, ".done_f")
        @info "Downsampling the data (female) to 16kHz"
        for fname in readdir(orig_wavs_f)
            if endswith(fname, "wav")
                out = joinpath(wavdir, fname)
                run(`sox -G $(joinpath(orig_wavs_f, fname)) -t wav $out rate 16k`)
            end
        end
        open(joinpath(wavdir, ".done_f"), "w") do f print(f, now()) end
    else
        @info "Data (female) already downsampled to 16kHz."
    end

    if ! ispath(wavdir, ".done_m")
        @info "Downsampling the data (male) to 16kHz"
        for fname in readdir(orig_wavs_m)
            if endswith(fname, "wav")
                out = joinpath(wavdir, fname)
                run(`sox -G $(joinpath(orig_wavs_m, fname)) -t wav $out rate 16k`)
            end
        end
        open(joinpath(wavdir, ".done_m"), "w") do f print(f, now()) end
    else
        @info "Data (male) already downsampled to 16kHz."
    end

    ###################################################################
    # Prepare the full set.

    allwavs = filter(x -> endswith(x, "wav"), readdir(wavdir))

    set = "full"
    setdir = mkpath(joinpath(datadir, set))
    if ! ispath(joinpath(setdir, ".done"))
        @info "Preparing $set set..."

        uttids = []
        speakers = Set{String}()
        utt2spk = Dict()
        utt2wav = Dict()
        for fname in allwavs
            uttid = replace(fname, ".wav" => "")
            tokens = split(uttid, "_")
            speaker = tokens[2]
            push!(uttids, uttid)
            push!(speakers, speaker)
            utt2spk[uttid] = speaker
            utt2wav[uttid] = abspath(fname)
        end
        uttids = sort(uttids)
        speakers = sort([s for s in speakers])

        # uttids
        open(joinpath(setdir, LOCATIONS[:uttids]), "w") do f
            for uttid in uttids println(f, uttid) end
        end

        # wav.scp
        open(joinpath(setdir, LOCATIONS[:wavs]), "w") do f
            for uttid in uttids println(f, uttid, "\t", utt2wav[uttid]) end
        end

        # uttids_speakers
        open(joinpath(setdir, LOCATIONS[:utt2spk]), "w") do f
            for uttid in uttids println(f, uttid, "\t", utt2spk[uttid]) end
        end

        # speakers
        open(joinpath(setdir, LOCATIONS[:speakers]), "w") do f
            for s in speakers println(f, s) end
        end

        alifile = joinpath(setdir, LOCATIONS[:ali])
        newalifile = joinpath(setdir, "ali_new")
        run(`rm -f $alifile $newalifile`)
        run(`wget -P $(abspath(setdir)) $URL_ALI`)
        write(newalifile, read(pipeline(`cat $alifile`, `sed s/sil/\<sil\>/g`)))
        run(`mv $newalifile $alifile`)

        # trans.wrd
        write(joinpath(rawdata, "line_index.tsv"),
              read(`cat $rawdata/line_index_male.tsv $rawdata/line_index_female.tsv`))
        cmd = pipeline(`cat $rawdata/line_index.tsv`, `sed "s/$(Char('\t'))/ /g"`)
        write(joinpath(setdir, "trans.wrd"), read(cmd))

        open(joinpath(setdir, ".done"), "w") do f print(f, now()) end
    else
        @info "$set already prepared"
    end

    ###################################################################
    # Prepare the lang directory.

    langdir = mkpath(joinpath(datadir, LOCATIONS[:lang]))
    @info "preparing $langdir..."

    phones = Set()
    open(joinpath(datadir, set, LOCATIONS[:ali]), "r") do f
        for line in eachline(f)
            tokens = split(line)
            for phone in tokens[2:end]
                push!(phones, phone)
            end
        end
    end

    open(joinpath(langdir, LOCATIONS[:units]), "w") do f
        for p in sort(collect(phones), lt = natural)
            type = p ∈ NON_SPEECH_UNITS ? LABELS[:nonspeechunit] : LABELS[:speechunit]
            println(f, p, "\t", type)
        end
    end

end

end # module

max_issues_repo_name: lucasondel/SpeechDatasets

max_forks_repo_name: lucasondel/SpeechDatasets


text:
import LinearAlgebra: norm, I
@testset "Laplacian in two-d" begin
    q = 5
    A, coarse_domains, scales, basis_functions, multicolor_ordering, fine_domains = CompressingSolvers.FD_Laplacian_subdivision_2d(q)
end


text:
# Autogenerated wrapper script for SDL2_ttf_jll for i686-w64-mingw32
export libsdl2_ttf

using Bzip2_jll
using FreeType2_jll
using Glib_jll
using Graphite2_jll
using libpng_jll
using PCRE_jll
using SDL2_jll
JLLWrappers.@generate_wrapper_header("SDL2_ttf")
JLLWrappers.@declare_library_product(libsdl2_ttf, "SDL2_ttf.dll")
function __init__()
    JLLWrappers.@generate_init_header(Bzip2_jll, FreeType2_jll, Glib_jll, Graphite2_jll, libpng_jll, PCRE_jll, SDL2_jll)
    JLLWrappers.@init_library_product(
        libsdl2_ttf,
        "bin\\SDL2_ttf.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export @checked, @debug_ccall

"""
    @checked function foo(...)
        rv = ...
        return rv
    end

Macro for wrapping a function definition returning a status code. Two versions of the
function will be generated: `foo`, with the function body wrapped by an invocation of the
`@check` macro (to be implemented by the caller of this macro), and `unsafe_foo` where no
such invocation is present and the status code is returned to the caller.
"""
macro checked(ex)
    # parse the function definition
    @assert Meta.isexpr(ex, :function)
    sig = ex.args[1]
    @assert Meta.isexpr(sig, :call)
    body = ex.args[2]
    @assert Meta.isexpr(body, :block)

    # generate a "safe" version that performs a check
    safe_body = quote
        @check $body
    end
    safe_sig = Expr(:call, sig.args[1], sig.args[2:end]...)
    safe_def = Expr(:function, safe_sig, safe_body)

    # generate a "unsafe" version that returns the error code instead
    unsafe_sig = Expr(:call, Symbol("unsafe_", sig.args[1]), sig.args[2:end]...)
    unsafe_def = Expr(:function, unsafe_sig, body)

    return esc(:($safe_def, $unsafe_def))
end

macro debug_ccall(ex)
    @assert Meta.isexpr(ex, :(::))
    call, ret = ex.args
    @assert Meta.isexpr(call, :call)
    target, argexprs... = call.args
    args = map(argexprs) do argexpr
        @assert Meta.isexpr(argexpr, :(::))
        argexpr.args[1]
    end

    ex = Expr(:macrocall, Symbol("@ccall"), __source__, ex)

    quote
        print($(string(target)), '(')
        for (i, arg) in enumerate(($(map(esc, args)...),))
            i > 1 && print(", ")
            render_arg(stdout, arg)
        end
        print(')')

        rv = $(esc(ex))

        println(" = ", rv)
        for (i, arg) in enumerate(($(map(esc, args)...),))
            if arg isa Base.RefValue
                println(" $i: ", arg[])
            end
        end
        rv
    end
end

render_arg(io, arg) = print(io, arg)
render_arg(io, arg::Union{<:Base.RefValue, AbstractArray}) = summary(io, arg)


text:
function benchmarksamplefunction(f, pc, params, percents=[99,90,60,30,10,1]; seed=1234)
    octree = RANSAC.buildoctree(pc)
    for percent in percents
        enablepoints!(pc, percent, seed)
        drawN = params.iteration.drawN
        sd = Vector{Int}(undef, drawN)
        Random.seed!(seed)
        countvalidsampling(f, pc, octree, params, 10000)
        Random.seed!(seed)
        result = @benchmark $f($sd, $pc, $octree, $params)
        display(result)
    end
end

function enablepoints!(pc, percent, seed = nothing)
    if ! (seed === nothing)
        Random.seed!(seed)
    end
    for i in eachindex(pc.isenabled)
        pc.isenabled[i] = rand() < percent/100 ? true : false
    end
    println("$(100*count(pc.isenabled)/pc.size) % is enabled")
    return pc
end

function countvalidsampling(f, pc, octree, params, n)
    drawN = params.iteration.drawN
    sd = Vector{Int}(undef, drawN)
    res = [f(sd, pc, octree, params) for i in 1:n]
    # valids
    ns = count(i->i[1], res)
    # false because simiar
    sims = count(i->((! i[1]) && i[2]==1), res)
    println("$(100*ns/n) % is valid; $(100*sims/(n-ns)) % of invalid is similar")
end


text:
## =============================================================================
# Hidden markov model with zero inflated likelihood

struct RandomWalkProblem{Y,N,I,L,IL,P} <: SurveyProblem
    y::Y
    n::N
    indices::I
    lnk::L
    invlnk::IL
    prior::P
end

function RandomWalkProblem(y,n,indices, lnk, invlnk)
    prior = PriorRandomWalkProblem(n, invlnk)
    RandomWalkProblem(y,n,indices, lnk, invlnk, prior)
end

struct PriorRandomWalkProblem{P,L,F,S,N,I}
    prior_p::P
    prior_μ0::L
    prior_ϕ::F
    prior_σ::S
    n::N
    invlnk::I
end

function PriorRandomWalkProblem(n, invlnk)
    prior_p  = Beta()
    prior_μ0 = truncated(Normal(10,10),0,100)
    prior_ϕ  = truncated(GammaMeanCv(50,.5),1.,Inf)
    prior_σ  = truncated(Normal(0.2, .1), 0., Inf)
    PriorRandomWalkProblem( prior_p, prior_μ0, prior_ϕ, prior_σ, n, invlnk)
end

function Base.rand(priors::PriorRandomWalkProblem)
    @unpack prior_p, prior_μ0, prior_ϕ, prior_σ, n, invlnk = priors
    p        = rand(prior_p)
    μ0       = rand(prior_μ0)
    ϕ        = rand(prior_ϕ)
    σ        = rand(prior_σ)
    latent_μ = rand( CovidSurvey.RandomWalk(n-1, σ, invlnk(μ0)) )
    return (; p, μ0, ϕ, σ, latent_μ)
end

function Distributions.logpdf(prior::PriorRandomWalkProblem, θ)
    @unpack p, μ0, ϕ, σ, latent_μ = θ
    @unpack prior_p, prior_μ0, prior_ϕ, prior_σ, n, invlnk = prior
    ℓ  = logpdf(prior_p, p)
    ℓ += logpdf(prior_μ0, μ0)
    ℓ += logpdf(prior_ϕ, ϕ)
    ℓ += logpdf(prior_σ, σ)
    ℓ += logpdf(CovidSurvey.RandomWalk(n-1, σ, invlnk(μ0)), latent_μ)
    return ℓ
end

likelihood_dist(p, μ, ϕ) =
    ZeroInflatedDistribution.(
        Ref(Bernoulli(p)),
        NegativeBinomial2.(μ, Ref(ϕ)
    )
)

function lnk!(μ, problem::RandomWalkProblem, θ)
    @unpack μ0, latent_μ = θ
    @unpack lnk = problem

    μ[1] = μ0
    for (i, x) in enumerate(latent_μ)
        μ[i+1] = lnk(x)
    end
    return nothing
end

function (problem::RandomWalkProblem)(θ)
    @unpack p, μ0, ϕ, σ, latent_μ  = θ               # extract the parameters
    @unpack y, n, indices, lnk, invlnk, prior = problem
    !(ϕ > 0) && (@error "not ϕ>0"; return -Inf)
    !(0 <= p <= 1) && (@error "not 0 <= p <= 1"; return -Inf)

    ℓ  = logpdf(prior, θ)

    μ = Vector{eltype(latent_μ)}(undef, n)
    lnk!(μ, problem, θ)
    !all(0 .<= μ ) && (@error "not all(0 .<= μ )"; return -Inf)

    d = likelihood_dist(p, μ, ϕ)
    ℓ += logpdf(arraydist(d[indices]),y)
    return ℓ
end

function problem_transformation(p::RandomWalkProblem)
    as((p = as𝕀, μ0 = as(Real, 0., 100.), ϕ = asℝ₊, σ = asℝ₊, latent_μ = as(Array, p.n-1)))
end

function RandomWalkProblem(data, col)
    y = data[:,col]
    indices = date2index( data.date )
    n = last(indices)
    lnk = KLogistic(100.)
    invlnk = KLogit(100.)
    return RandomWalkProblem(y, n, indices, lnk, invlnk)
end

function generate_quantities(problem::RandomWalkProblem, θ::NamedTuple)
    @unpack n = problem

    μ = zeros(n)
    lnk!(μ, problem, θ)
    return μ
end

function StatsBase.predict(problem::RandomWalkProblem, θ)
    @unpack p, ϕ = θ               # extract the parameters
    @unpack y, n, indices = problem

    μ = Vector{typeof(ϕ)}(undef, n)
    lnk!(μ, problem, θ)

    d = likelihood_dist(p, μ, ϕ)
    y = rand(arraydist(d[indices]))
    return y
end

max_stars_repo_name: andreaskoher/CovidSurvey

max_issues_repo_name: andreaskoher/CovidSurvey

max_forks_repo_name: andreaskoher/CovidSurvey


text:
using ScrewDimensions
using Test

@testset "ScrewDimensions.jl" begin
    # Write your tests here.
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 5934500e-fbe3-11ea-3ea2-e7f2c1cf9359
begin
	using Plots; gr()
	using ForwardDiff
	using DifferentialEquations
	using LinearAlgebra
	using Interpolations
end

# ╔═╡ 80a89f0e-fc02-11ea-3924-a11494022a97
md"""
# Ray Tracing
This notebook is completely self-contained. It solves the eikonal equation for tracing rays in a cylindrically-symmetric ocean environment, as a 2D slice of ocean (ocn). It includes the bathymetry (bty), the altimetry (ati) and variable ocean sound speed, all range and depth dependent.

The eikonal equation has been expressed in terms of path-length as a system of first order ODEs. This enables simple utilization of the `DifferentialEquations.jl` package.

For $r$ the range, $z$ the depth, $c(r, z)$ the sound speed field, $s$ the path length along the ray, and $(\xi, \zeta)$ a tangent to the ray, the equations are

```math
\newcommand{\Dif}[2]{\frac{d{#1}}{d{#2}}}
\newcommand{\Part}[2]{\frac{\partial{#1}}{\partial{#2}}}
\begin{aligned}
\Dif{r}{s} &= c(r, z) \xi(s) &
\Dif{\xi}{s} &= \frac{-1}{c(r, z)^2} \Part{c}{r} \\
\Dif{z}{s} &= c(r, z) \zeta(s) &
\Dif{\zeta}{s} &= \frac{-1}{c(r, z)^2} \Part{c}{z}
\end{aligned}
```

The functionality demonstrated here has been successfully run without `struct`s. This notebook is my attempt to port that working code to a function that takes in `struct`s.
"""

# ╔═╡ a47ae3c0-fc03-11ea-2322-a5e2ddd51990
md"## Augmenting Code"

# ╔═╡ 735054c0-fc02-11ea-10fe-0b7a22aa8dca
md"### Preamble"

# ╔═╡ e9036770-fc02-11ea-3e27-379927eec193
md"### Generic Variables and Augmenting Functions"

# ╔═╡ 9cc6fab2-fbe5-11ea-2b09-cfb7ec9c627f
begin
	rPeak = 5e3
	rMax = 10e3
	zMax = 1e3
	zMin = 8e2
	zAtiMin = -10
	zAtiMax = 50
	r = range(0, rMax, length = 100)
	z = range(zAtiMin, zMax, length = 51)
end

# ╔═╡ e6961940-fcb8-11ea-2bfa-f32a7d60f480
begin
	function InterpolatingFunction(rng, val)
		Itp = LinearInterpolation(rng, val, extrapolation_bc = Flat())
		return ItpFcn(r) = Itp(r)
	end

	function InterpolatingFunction(rng, dpt, val)
		display(rng)
		display(dpt)
		Itp = LinearInterpolation((dpt, rng), val, extrapolation_bc = Flat())
		return ItpFcn(r, z) = Itp(z, r)
	end
end

# ╔═╡ f43d8cb0-fc02-11ea-2271-0dca2ac46977
md"### Ocean Environment: Functions and Structs"

# ╔═╡ 8e0bc300-fbfb-11ea-374f-9bcf6e06ba54
"""
	t_rfl::Vector = boundary_reflection(t_inc::Vector, t_bnd::Vector)

`t_inc`    tangent vector of incident ray
`t_bnd`    tangent vector of boundary
`t_rfl`    tangent vector of reflected ray
"""
function boundary_reflection(t_inc::Vector, t_bnd::Vector)
	n_bnd = [-t_bnd[2], t_bnd[1]]
# 	t_rfl = t_inc - 2(t_inc ⋅ n_bnd)*n_bnd
	t_rfl = t_inc - 2LinearAlgebra.dot(t_inc, n_bnd)*n_bnd
	println(t_inc)
	println(t_bnd)
	println(t_rfl)
	return t_rfl
end

# ╔═╡ 17994f60-fbe4-11ea-1bd3-3925f7a09733
begin
	"""
		Boundary(z)

	`z(r)`    depth of boundary in terms of range `r`
	"""
	struct Boundary
		z::Function
		dzdr::Function
		condition::Function
		affect!::Function
	end
	function Boundary(z::Function)
		dzdr(r) = ForwardDiff.derivative(z, r)
		condition(u, t, ray) = z(u[1]) - u[2]
		affect!(ray) = ray.u[3], ray.u[4] = boundary_reflection([ray.u[3], ray.u[4]], [1, dzdr(ray.u[1])])
		return Boundary(z, dzdr, condition, affect!)
	end
	function Boundary(r::AbstractArray)
		r_ = [rng for rng ∈ r[:, 1]]
		z_ = [dpt for dpt ∈ r[:, 2]]
		zFcn = InterpolatingFunction(r_, z_)
		return Boundary(zFcn)
	end
end

# ╔═╡ e4ad3740-fc03-11ea-3cec-27b6ed6a7704
md"""
#### Bathymetry
A Gaussian seabed shape is chosen.
"""

# ╔═╡ 82b8dd10-fbe4-11ea-292c-710024b9f9eb
begin
	Aᵣ = (2rPeak/3)^2/log((9zMax - 11zMin)/(10(zMax - zMin)))
	zBty(r) = zMax - (zMax - zMin)*exp(-(r - rPeak)^2/4e5)
	bty = Boundary(zBty)
end

# ╔═╡ 9b62f350-fbe4-11ea-2bd8-a99db7360093
begin
	pBty = plot(r, bty.z,
		legend = false,
		yaxis = (:flip, (0, zMax), raw"$z$"),
		title = "Bathymetry")
	pdBty = plot(r, bty.dzdr,
		legend = false,
		yaxis = (:flip, raw"$\partial{z}/\partial{r}$"),
		xaxis = "Range (m)")
	lBty = @layout [a; b]
	plot(pBty, pdBty, layout = lBty)
end

# ╔═╡ f2713bb0-fc03-11ea-1213-0f56f7df8683
md"""
#### Altimetry
A sinusoidal sea surface is chosen.
"""

# ╔═╡ 46a62e80-fbe5-11ea-3dd8-014cec82f60f
begin
	zAti(r) = zAtiMin + (zAtiMax - zAtiMin)*(sin(r/1e3) + 1.)/2
	ati = Boundary(zAti)
end

# ╔═╡ 97152790-fbe5-11ea-1ca6-a5018f9e1228
begin
	pAti = plot(r, ati.z,
		legend = false,
		yaxis = (:flip, (zAtiMin, zMax), raw"$z$"),
		title = "Altimetry")
	pdAti = plot(r, ati.dzdr,
		legend = false,
		yaxis = (:flip, raw"$\partial{z}/\partial{r}$"),
		xaxis = raw"$r$")
	lAti = @layout [a; b]
	plot(pAti, pdAti, layout = lAti)
end

# ╔═╡ 0c53d1f0-fc04-11ea-1996-c50e9f3050f4
md"#### Ocean Sound Speed"

# ╔═╡ cb5709b0-fbdb-11ea-1d9b-375f555a2a8f
begin
	"""
		Medium(c, R)

	`c(r, z)`    sound speed at range `r` and depth `z`
	`R`          maximum range of slice
	"""
	struct Medium
		c::Function
		∂c∂r::Function
		∂c∂z::Function
		R::Real
	end
	function Medium(c::Function, R::Real)
		c_(x) = c(x[1], x[2])
		∇c_(x) = ForwardDiff.gradient(c_, x)
		∇c(r, z) = ∇c_([r, z])
		∂c∂r(r, z) = ∇c(r, z)[1]
		∂c∂z(r, z) = ∇c(r, z)[2]
		return Medium(c, ∂c∂r, ∂c∂z, R)
	end
	function Medium(c::AbstractArray, R::Real)
		r_ = [rc for rc ∈ c[1, 2:end]]
		z_ = [zc for zc ∈ c[2:end, 1]]
		c_ = c[2:end, 2:end]

		cFcn = InterpolatingFunction(r_, z_, c_)
		return Medium(cFcn, R)
	end
end

# ╔═╡ 0ff87620-fbe7-11ea-0ef0-6b827f36c4a5
begin
	cMin = 1500
	cMax = 1600
	C(r) = [1 zAti(r) zAti(r)^2
		1 (zAti(r) + zBty(r))/2 ((zAti(r) + zBty(r))/2)^2
		1 zBty(r) zBty(r)^2]
	C_(r) = C(r)\[cMax, cMin, cMax]
	c₀(r) = C_(r)[1]
	c₁(r) = C_(r)[2]
	c₂(r) = C_(r)[3]
	c(r, z) = c₀(r) + c₁(r)*z + c₂(r)*z^2
	ocn = Medium(c, rMax)
end

# ╔═╡ 41932180-fbe7-11ea-05be-15c01a179539
begin
	pc = heatmap(r, z, ocn.c,
		yaxis = (:flip, raw"$z$"),
		title = "Celerity")
	plot!(r, zBty, color = :red, label = "")
	plot!(r, zAti, color = :white, label = "")
	pdcr = heatmap(r, z, ocn.∂c∂r,
		yaxis = (:flip, raw"$z$"))
	plot!(r, zBty, color = :red, label = "")
	plot!(r, zAti, color = :white, label = "")
	pdcz = heatmap(r, z, ocn.∂c∂z,
		yaxis = (:flip, raw"$z$"),
		xaxis = raw"$r$")
	plot!(r, zBty, color = :red, label = "Bathymetry")
	plot!(r, zAti, color = :white, label = "Altimetry")

	lc = @layout [a; b; c]
	plot(pc, pdcr, pdcz, layout = lc)
end

# ╔═╡ 694b88f0-fc07-11ea-2f5b-93fcfa7c1463
md"I've designed the sound speed to be parabolic but varying with range such that the maximum value of $cMax is at the boundaries, and $cMin is halfway between."

# ╔═╡ fcc629b0-fc06-11ea-3eaa-adf50eadbce0
begin
	pcAtAti = plot(r, r -> c(r, zAti(r)),
		yaxis = raw"$c\left(r, z_\textrm{ati}\right)$",
		legend = false)
	pcAtHalf = plot(r, r -> c(r, (zBty(r) + zAti(r))/2),
		yaxis = raw"$c\left(r, \frac{z_\textrm{ati} + z_\textrm{bty}}{2}\right)$",
		legend = false)
	pcAtBty = plot(r, r -> c(r, zBty(r)),
		yaxis = raw"$c\left(r, z_\textrm{bty}\right)$",
		xaxis = raw"$r$",
		legend = false)
	lcAt = @layout [a; b; c]
	plot(pcAtAti, pcAtHalf, pcAtBty, layout = lcAt)
end

# ╔═╡ 19e1f950-fc04-11ea-1287-5b8d3307de12
md"#### Sound Source"

# ╔═╡ bcdf3810-fbfb-11ea-015c-998ab7c3d380
"""
	Entity(r, z, θ)

`r`    range of entity
`z`    depth of entity
`θ`    angle of ray
"""
struct Entity
	r::Real
	z::Real
end

# ╔═╡ 84f01650-fc0e-11ea-2861-079e697543db
begin
	r₀ = 0.0
	z₀ = (zBty(r₀) + zAti(r₀))/2
	src = Entity.(r₀, z₀)
end

# ╔═╡ 0ca4cd90-fc03-11ea-3d70-03cf319100f2
md"""
## Acoustic Propagation
The solution implementation is similar to the general Julia differential equation solution form of first defining the problem, then computing the solution.
"""

# ╔═╡ 4d44d700-fc03-11ea-1914-e30845160f20
md"### Problem"

# ╔═╡ 3d33d420-fbe9-11ea-3bc2-bd9475b53a2e
function acoustic_propagation_problem(
		θ₀::Real,
		Src::Entity,
		Ocn::Medium,
		Bty::Boundary,
		Ati::Boundary)

	function eikonal!(du, u, p, s)
		r = u[1]
		z = u[2]
		ξ = u[3]
		ζ = u[4]
		τ = u[5]
		du[1] = dr = Ocn.c(r, z)*ξ
		du[2] = dz = Ocn.c(r, z)*ζ
		du[3] = dξ = -Ocn.∂c∂r(r, z)/Ocn.c(r, z)^2
		du[4] = dζ = -Ocn.∂c∂z(r, z)/Ocn.c(r, z)^2
		du[5] = dτ = 1/Ocn.c(r, z)
	end

	rng_condition(u, t, ray) = Ocn.R/2 - abs(u[1] - Ocn.R/2)
	rng_affect!(ray) = terminate!(ray)
	CbRng = ContinuousCallback(rng_condition, rng_affect!)
	CbBty = ContinuousCallback(Bty.condition, Bty.affect!)
	CbAti = ContinuousCallback(Ati.condition, Ati.affect!)
	CbBnd = CallbackSet(CbRng, CbBty, CbAti)

	r₀ = Src.r
	z₀ = Src.z
	ξ₀ = cos(θ₀)/Ocn.c(r₀, z₀)
	ζ₀ = sin(θ₀)/Ocn.c(r₀, z₀)
	τ₀ = 0.
	u₀ = [r₀, z₀, ξ₀, ζ₀, τ₀]

	TLmax = 100
	S = 10^(TLmax/10)
	sSpan = (0., S)

	prob_eikonal = ODEProblem(eikonal!, u₀, sSpan)

	return prob_eikonal, CbBnd
end

# ╔═╡ 3e7f8f30-fc03-11ea-0f07-bd4c7b8e7ce2
md"### Solution"

# ╔═╡ 78232bd0-fbfd-11ea-0de5-7fb83321c1f0
function solve_acoustic_propagation(prob_eikonal, CbBnd)
	@time RaySol = solve(prob_eikonal, callback = CbBnd)
	return RaySol
end

# ╔═╡ d268f650-fc0d-11ea-151e-fd59c8065050
struct Ray
	θ₀
	Sol
	function Ray(θ₀, Src::Entity, Ocn::Medium, Bty::Boundary, Ati::Boundary)
		Prob, CbBnd = acoustic_propagation_problem(θ₀, Src, Ocn, Bty, Ati)
		Sol = solve_acoustic_propagation(Prob, CbBnd)
		return new(θ₀, Sol)
	end
end

# ╔═╡ 6bf42d00-fc02-11ea-0cc1-27333ef247b9
md"""
### Single Ray
Here we demonstrate the working code for a single ray.
"""

# ╔═╡ a9618660-fbfd-11ea-0c25-3bdaf8e181d2
begin
	θ₀ = acos(c(r₀, z₀)/cMax)
	ray = Ray(θ₀, src, ocn, bty, ati)
end

# ╔═╡ bbc9d5d0-fc0e-11ea-0187-a3a914b19d12
begin
	plot(r, zAti, label = "Altimetry")
	plot!(r, zBty, label = "Bathymetry")
	plot!(ray.Sol, vars = (1, 2), yaxis = :flip, label = "")
end

# ╔═╡ 5c903430-fc02-11ea-28fd-df2119a5b104
md"""
### Multirays
Next we demonstrate the working code for multiple rays. Hopefully this can work via broadcasting.
"""

# ╔═╡ 5ac1f0f0-fc0f-11ea-0d9f-d17ef57c699c
begin
	Base.broadcastable(m::Entity) = Ref(m)
	Base.broadcastable(m::Medium) = Ref(m)
	Base.broadcastable(m::Boundary) = Ref(m)
end

# ╔═╡ 08d817b0-fc0f-11ea-375d-4154a74e15ef
rays = Ray.(θ₀.*(-1.5:0.125:1.5), src, ocn, bty, ati)

# ╔═╡ 72debb50-fc0f-11ea-2192-77bcc5c839f0
begin
	pt = plot(yaxis = :flip)
	plot!(r, zAti, label = "Altimetry")
	plot!(r, zBty, label = "Bathymetry")
	for nRay = 1:length(rays)
		plot!(rays[nRay].Sol, vars = (1, 2), label = "")
	end
	pt
end

# ╔═╡ cb685e02-fcbb-11ea-2e8a-43da1616798e
md"""
### Interpolation
Now we take in vectors of values for the range and depth dependent ocean elements.
"""

# ╔═╡ 26609b80-fcbf-11ea-13f7-71bbfc3d08b3
begin
	rVec = [0, 10, 20, 50, 100, 200, 500, 1e3, 2e3, 5e3, 1e4]
	zVec = [0, 10, 20, 100, 500, 1e3]
end

# ╔═╡ 36f265f2-fcbf-11ea-1d21-87106ecdf10e
md"#### Ocean Sound Speed"

# ╔═╡ 28d91cf0-fcbc-11ea-0b90-41e002050d51
begin
	cMat_ = [1550 + sin(r/1e2) + sin(z/1e2) for z ∈ zVec, r ∈ rVec]
	cMat = vcat(hcat(0, transpose(rVec)), hcat(zVec, cMat_))
end

# ╔═╡ 989d5e10-fcbd-11ea-2e0e-5f918774e04b
heatmap(rVec, zVec, cMat_, yaxis = :flip)

# ╔═╡ 3bf2a3b0-fcbc-11ea-0f1b-73c25779bc06
ocnMat = Medium(cMat, rMax)

# ╔═╡ a586af60-fcbc-11ea-110b-7de105f36ff7
begin
	rNew = range(0, rVec[end], length = 201)
	zNew = range(0, zVec[end], length = 101)
	pcMat = heatmap(rNew, zNew, ocnMat.c, yaxis = :flip)
end

# ╔═╡ d41640d0-fcc0-11ea-20f8-6fd11922f72b
heatmap(rNew, zNew, ocnMat.∂c∂r)

# ╔═╡ 4644ab80-fcbf-11ea-3140-f1b695b3b20d
md"#### Ocean Boundaries"

# ╔═╡ 4fd36aae-fcbf-11ea-0017-ad3e19160af2
begin
	zAtiVec_ = [10(1 - sin(r/2e3)) for r ∈ rVec]
	zAtiVec = hcat(rVec, zAtiVec_)
end

# ╔═╡ a805e35e-fcc0-11ea-0efd-cf82b90f857b
atiVec = Boundary(zAtiVec)

# ╔═╡ d0d4f880-fcc0-11ea-2ca2-5f761b9af8d9
begin
	scatter(rVec, zAtiVec_)
	plot!(rNew, atiVec.z)
end

# ╔═╡ 3b145f60-fcc1-11ea-25d0-332f3ad41f3a
begin
	zBtyVec_ = [zVec[end] - 10sin(r/1e2)^2 for r ∈ rVec]
	zBtyVec = hcat(rVec, zBtyVec_)
end

# ╔═╡ 9defd650-fcc1-11ea-17f0-9752875085d2
btyVec = Boundary(zBtyVec)

# ╔═╡ b282b6a0-fcc1-11ea-24be-09f82bb1cbc4
begin
	scatter(rVec, zBtyVec_)
	plot!(rNew, btyVec.z)
end

# ╔═╡ d720a530-fcc1-11ea-0afb-0df93dd0d9a4
begin
	srcVec = Entity(0.0, 500.)
	θ₀Vec = π/5 .*(-1.5:0.125:1.5)
# 	θ₀Vec = π/5*(-0.125) # issue
end

# ╔═╡ 1c594990-fcc2-11ea-02fb-c36fcef2e847
raysVec = Ray.(θ₀Vec, srcVec, ocnMat, btyVec, atiVec)

# ╔═╡ 58cc1bb0-fcd0-11ea-0c35-7d16b2ce2599
md"Notice a few incorrect rays penetrating the boundaries, but reflecting back."

# ╔═╡ c0b2c252-fccc-11ea-3a60-352b68b6ccc7
begin
	ptVec = plot(yaxis = :flip)
	plot!(rVec, zAtiVec_, label = "Altimetry")
	plot!(rVec, zBtyVec_, label = "Bathymetry")
	for nRay = 1:length(raysVec)
		plot!(raysVec[nRay].Sol, vars = (1, 2), label = "")
	end
# 	plot!(raysVec.Sol, vars = (1, 2), label = "")
	ptVec
end

# ╔═╡ Cell order:
# ╟─80a89f0e-fc02-11ea-3924-a11494022a97
# ╟─a47ae3c0-fc03-11ea-2322-a5e2ddd51990
# ╟─735054c0-fc02-11ea-10fe-0b7a22aa8dca
# ╠═5934500e-fbe3-11ea-3ea2-e7f2c1cf9359
# ╟─e9036770-fc02-11ea-3e27-379927eec193
# ╠═9cc6fab2-fbe5-11ea-2b09-cfb7ec9c627f
# ╠═e6961940-fcb8-11ea-2bfa-f32a7d60f480
# ╟─f43d8cb0-fc02-11ea-2271-0dca2ac46977
# ╠═8e0bc300-fbfb-11ea-374f-9bcf6e06ba54
# ╠═17994f60-fbe4-11ea-1bd3-3925f7a09733
# ╟─e4ad3740-fc03-11ea-3cec-27b6ed6a7704
# ╠═82b8dd10-fbe4-11ea-292c-710024b9f9eb
# ╠═9b62f350-fbe4-11ea-2bd8-a99db7360093
# ╟─f2713bb0-fc03-11ea-1213-0f56f7df8683
# ╠═46a62e80-fbe5-11ea-3dd8-014cec82f60f
# ╠═97152790-fbe5-11ea-1ca6-a5018f9e1228
# ╟─0c53d1f0-fc04-11ea-1996-c50e9f3050f4
# ╠═cb5709b0-fbdb-11ea-1d9b-375f555a2a8f
# ╠═0ff87620-fbe7-11ea-0ef0-6b827f36c4a5
# ╠═41932180-fbe7-11ea-05be-15c01a179539
# ╟─694b88f0-fc07-11ea-2f5b-93fcfa7c1463
# ╠═fcc629b0-fc06-11ea-3eaa-adf50eadbce0
# ╟─19e1f950-fc04-11ea-1287-5b8d3307de12
# ╠═bcdf3810-fbfb-11ea-015c-998ab7c3d380
# ╠═84f01650-fc0e-11ea-2861-079e697543db
# ╟─0ca4cd90-fc03-11ea-3d70-03cf319100f2
# ╟─4d44d700-fc03-11ea-1914-e30845160f20
# ╠═3d33d420-fbe9-11ea-3bc2-bd9475b53a2e
# ╟─3e7f8f30-fc03-11ea-0f07-bd4c7b8e7ce2
# ╠═78232bd0-fbfd-11ea-0de5-7fb83321c1f0
# ╠═d268f650-fc0d-11ea-151e-fd59c8065050
# ╟─6bf42d00-fc02-11ea-0cc1-27333ef247b9
# ╠═a9618660-fbfd-11ea-0c25-3bdaf8e181d2
# ╠═bbc9d5d0-fc0e-11ea-0187-a3a914b19d12
# ╟─5c903430-fc02-11ea-28fd-df2119a5b104
# ╠═5ac1f0f0-fc0f-11ea-0d9f-d17ef57c699c
# ╠═08d817b0-fc0f-11ea-375d-4154a74e15ef
# ╠═72debb50-fc0f-11ea-2192-77bcc5c839f0
# ╟─cb685e02-fcbb-11ea-2e8a-43da1616798e
# ╠═26609b80-fcbf-11ea-13f7-71bbfc3d08b3
# ╟─36f265f2-fcbf-11ea-1d21-87106ecdf10e
# ╠═28d91cf0-fcbc-11ea-0b90-41e002050d51
# ╠═989d5e10-fcbd-11ea-2e0e-5f918774e04b
# ╠═3bf2a3b0-fcbc-11ea-0f1b-73c25779bc06
# ╠═a586af60-fcbc-11ea-110b-7de105f36ff7
# ╠═d41640d0-fcc0-11ea-20f8-6fd11922f72b
# ╟─4644ab80-fcbf-11ea-3140-f1b695b3b20d
# ╠═4fd36aae-fcbf-11ea-0017-ad3e19160af2
# ╠═a805e35e-fcc0-11ea-0efd-cf82b90f857b
# ╠═d0d4f880-fcc0-11ea-2ca2-5f761b9af8d9
# ╠═3b145f60-fcc1-11ea-25d0-332f3ad41f3a
# ╠═9defd650-fcc1-11ea-17f0-9752875085d2
# ╠═b282b6a0-fcc1-11ea-24be-09f82bb1cbc4
# ╠═d720a530-fcc1-11ea-0afb-0df93dd0d9a4
# ╠═1c594990-fcc2-11ea-02fb-c36fcef2e847
# ╟─58cc1bb0-fcd0-11ea-0c35-7d16b2ce2599
# ╠═c0b2c252-fccc-11ea-3a60-352b68b6ccc7


text:
module RulesBetaOutTest

using Test
using ReactiveMP
using Random
import ReactiveMP: @test_rules

@testset "rules:Beta:out" begin

     @testset "Belief Propagation: (m_a::PointMass, m_b::PointMass)" begin

        @test_rules [ with_float_conversions = true ] Beta(:out, Marginalisation) [
            (input = (m_a = PointMass(1.0), m_b = PointMass(2.0)), output = Beta(1.0, 2.0)),
            (input = (m_a = PointMass(2.0), m_b = PointMass(2.0)), output = Beta(2.0, 2.0)),
            (input = (m_a = PointMass(3.0), m_b = PointMass(3.0)), output = Beta(3.0, 3.0)),
        ]

    end

    @testset "Variational Message Passing: (q_a::PointMass, q_b::PointMass)" begin

        @test_rules [ with_float_conversions = true ] Beta(:out, Marginalisation) [
            (input = (q_a = PointMass(1.0), q_b = PointMass(2.0)), output = Beta(1.0, 2.0)),
            (input = (q_a = PointMass(2.0), q_b = PointMass(2.0)), output = Beta(2.0, 2.0)),
            (input = (q_a = PointMass(3.0), q_b = PointMass(3.0)), output = Beta(3.0, 3.0)),
        ]

    end
end
end


text:
"""
    test_constvel(; kwargs...)

Advects a gaussian concentration c0(x, y, t) with a constant velocity flow
u(x, y) = uvel and v(x, y) = vvel and compares the final state with
cfinal = c0(x-uvel*tfinal, y-vvel*tfinal)
"""
function test_constvel(stepper, dt, nsteps)

  nx, Lx = 128, 2π
  uvel, vvel = 0.2, 0.1
  u(x, y) = uvel
  v(x, y) = vvel

  prob = TracerAdvDiff.Problem(; nx=nx, Lx=Lx, kap=0.0, u=u, v=v, dt=dt, stepper=stepper, steadyflow=true)
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid

  x, y = gridpoints(g)

  σ = 0.1
  c0func(x, y) = 0.1*exp.(-(x.^2+y.^2)/(2σ^2))

  c0 = c0func.(x, y)
  tfinal = nsteps*dt
  cfinal = @. c0func(x - uvel*tfinal, y - vvel*tfinal)

  TracerAdvDiff.set_c!(prob, c0)

  stepforward!(prob, nsteps)
  TracerAdvDiff.updatevars!(prob)

  isapprox(cfinal, v.c, rtol=g.nx*g.ny*nsteps*1e-12)
end

"""
    test_timedependenttvel(; kwargs...)

Advects a gaussian concentration c0(x, y, t) with a time-varying velocity flow
u(x, y, t) = uvel and v(x, y, t) = vvel*sign(-t+tfinal/2) and compares the final
state with cfinal = c0(x-uvel*tfinal, y)
"""
function test_timedependentvel(stepper, dt, tfinal)

  nx, Lx = 128, 2π
  nsteps = round(Int, tfinal/dt)

  if !isapprox(tfinal, nsteps*dt, rtol=1e-12)
      error("tfinal is not multiple of dt")
  end

  uvel, vvel = 0.2, 0.1
  u(x, y, t) = uvel
  v(x, y, t) = t <= tfinal/2 ? vvel : -vvel

  prob = TracerAdvDiff.Problem(; nx=nx, Lx=Lx, kap=0.0, u=u, v=v, dt=dt, stepper=stepper)
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
  x, y = gridpoints(g)

  σ = 0.1
  c0func(x, y) = @. 0.1*exp(-(x^2+y^2)/(2σ^2))

  c0 = @. c0func(x, y)
  tfinal = nsteps*dt
  cfinal = @. c0func(x - uvel*tfinal, y)

  TracerAdvDiff.set_c!(prob, c0func)

  stepforward!(prob, nsteps)
  TracerAdvDiff.updatevars!(prob)

  isapprox(cfinal, v.c, rtol=g.nx*g.ny*nsteps*1e-12)
end

"""
    test_diffusion(; kwargs...)

Diffuses a gaussian concentration c0(x, y, t) and compares the final state with
the analytic solution of the heat equation, cfinal
"""
function test_diffusion(stepper, dt, tfinal; steadyflow = true)

  nx  = 128
  Lx  = 2π
  kap = 0.01
  nsteps = round(Int, tfinal/dt)

  if !isapprox(tfinal, nsteps*dt, rtol=1e-12)
      error("tfinal is not multiple of dt")
  end

  grid = TwoDGrid(nx, Lx)
  prob = TracerAdvDiff.Problem(; steadyflow=steadyflow, grid=grid, nx=nx,
    Lx=Lx, kap=kap, dt=dt, stepper=stepper)
  sol, cl, v, p, g = prob.sol, prob.clock, prob.vars, prob.params, prob.grid
  x, y = gridpoints(g)

  c0ampl, σ = 0.1, 0.1
  c0func(x, y) = @. c0ampl*exp(-(x^2+y^2)/(2σ^2))

  c0 = @. c0func.(x, y)
  tfinal = nsteps*dt
  σt = sqrt(2*kap*tfinal + σ^2)
  cfinal = @. c0ampl*σ^2/σt^2 * exp(-(x^2+y^2)/(2*σt^2))

  TracerAdvDiff.set_c!(prob, c0)

  stepforward!(prob, nsteps)
  TracerAdvDiff.updatevars!(prob)

  isapprox(cfinal, v.c, rtol=g.nx*g.ny*nsteps*1e-12)
end

"""
    test_hyperdiffusion(; kwargs...)

Diffuses a gaussian concentration c0(x, y, t) using hyperdiffusivity and
compares the final state with the analytic solution of the heat equation, cfinal
"""
function test_hyperdiffusion(stepper, dt, tfinal; steadyflow = true)

    nx  = 128
    Lx  = 2π
    kap = 0.0   # no diffusivity
    eta = kap   # no diffusivity
   kaph = 0.01  # hyperdiffusivity coeff
  nkaph = 1     # nkaph=1 converts hyperdiffusivity to plain diffusivity
                # so we can compare with the analytic solution of heat equation

  nsteps = round(Int, tfinal/dt)

  if !isapprox(tfinal, nsteps*dt, rtol=1e-12)
      error("tfinal is not multiple of dt")
  end

   g = TwoDGrid(nx, Lx)
  x, y = gridpoints(g)

  u, v = zero(x), zero(x) #0*x, 0*x

  vs = TracerAdvDiff.Vars(g)
  pr = TracerAdvDiff.ConstDiffSteadyFlowParams(eta, kap, kaph, nkaph, u, v)
  eq = TracerAdvDiff.Equation(pr, g)
  prob = FourierFlows.Problem(eq, stepper, dt, g, vs, pr)

  c0ampl, σ = 0.1, 0.1
  c0func(x, y) = @. c0ampl*exp(-(x^2+y^2)/(2σ^2))

  c0 = @. c0func(x, y)
  tfinal = nsteps*dt
  σt = sqrt(2*kaph*tfinal + σ^2)
  cfinal = @. c0ampl*σ^2/σt^2 * exp(-(x^2+y^2)/(2*σt^2))

  TracerAdvDiff.set_c!(prob, c0)

  stepforward!(prob, nsteps)
  TracerAdvDiff.updatevars!(prob)

  isapprox(cfinal, vs.c, rtol=g.nx*g.ny*nsteps*1e-12)
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

# All the types supported by mshadow. See `mshadow/base.h`
const DType = Union{Float32, Float64, Float16, UInt8, Int32, Int8, Int64}
@enum TypeFlag kFloat32 kFloat64 kFloat16 kUint8 kInt32 kInt8 kInt64
const DEFAULT_DTYPE = Float32  # MSHADOW_DEFAULT_DTYPE

function toTypeFlag(T::Type{<:DType})
  if T == Float32
    return kFloat32
  elseif T == Float64
    return kFloat64
  elseif T == Float16
    return kFloat16
  elseif T == UInt8
    return kUint8
  elseif T == Int32
    return kInt32
  elseif T == Int8
    return kInt8
  elseif T == Int64
    return kInt64
  else
    throw(ArgumentError("Can't convert $T to DType."))
  end
end

function fromTypeFlag(T::TypeFlag)
  if T == kFloat32
    return Float32
  elseif T == kFloat64
    return Float64
  elseif T == kFloat16
    return Float16
  elseif T == kUint8
    return UInt8
  elseif T == kInt32
    return Int32
  elseif T == kInt8
    return Int8
  elseif T == kInt64
    return Int64
  else
    throw(ArgumentError("Can't convert DType $T."))
  end
end

# create a NDArray handle of specific shape
function _ndarray_alloc(shape::NTuple{N,Int}, ctx::Context, delay_alloc::Bool) where N
  h_ref  = Ref{MX_handle}(0)
  shape  = collect(reverse(MX_uint.(shape)))
  @mxcall(:MXNDArrayCreate, (Ptr{MX_uint}, MX_uint, Cint, Cint, Cint, Ref{MX_handle}),
      shape, N, ctx.device_type, ctx.device_id, delay_alloc, h_ref)
  handle = MX_NDArrayHandle(h_ref[])
  return handle
end

# create a NDArray handle of specific shape type
function _ndarray_alloc(::Type{T}, shape::NTuple{N,Int}, ctx::Context, delay_alloc::Bool) where {T<:DType,N}
  h_ref  = Ref{MX_handle}(0)
  shape  = collect(reverse(MX_uint.(shape)))
  dtype  = toTypeFlag(T)
  @mxcall(:MXNDArrayCreateEx, (Ptr{MX_uint}, MX_uint, Cint, Cint, Cint, Cint, Ref{MX_handle}),
      shape, N, ctx.device_type, ctx.device_id, delay_alloc, dtype, h_ref)
  handle = MX_NDArrayHandle(h_ref[])
  return handle
end

# create a handle to an empty NDArray, this handle can be used to hold
# results returned by libmx API calls
function _ndarray_alloc()
  h_ref = Ref{MX_handle}(0)
  @mxcall(:MXNDArrayCreateNone, (Ref{MX_handle},), h_ref)
  return MX_NDArrayHandle(h_ref[])
end

################################################################################
# NDArray Type
################################################################################
"""
    NDArray{T,N}

Wrapper of the `NDArray` type in `libmxnet`. This is the basic building block
of tensor-based computation.

!!! note
      since C/C++ use row-major ordering for arrays while Julia follows a
      column-major ordering. To keep things consistent, we keep the underlying data
      in their original layout, but use *language-native* convention when we talk
      about shapes. For example, a mini-batch of 100 MNIST images is a tensor of
      C/C++/Python shape (100,1,28,28), while in Julia, the same piece of memory
      have shape (28,28,1,100).
"""
mutable struct NDArray{T,N}
  handle   :: MX_NDArrayHandle
  writable :: Bool

  NDArray{T,N}(handle::MX_NDArrayHandle, writable::Bool = true) where {T,N} =
    new(handle, writable)
end

# UndefInitializer constructors
NDArray{T,N}(::UndefInitializer, dims::NTuple{N,Integer};
             writable = true, ctx::Context = current_context()) where {T,N} =
  NDArray{T,N}(_ndarray_alloc(T, dims, ctx, false), writable)
NDArray{T,N}(::UndefInitializer, dims::Vararg{Integer,N}; kw...) where {T,N} =
  NDArray{T,N}(undef, dims; kw...)

NDArray{T}(::UndefInitializer, dims::NTuple{N,Integer}; kw...) where {T,N} =
  NDArray{T,N}(undef, dims; kw...)
NDArray{T}(::UndefInitializer, dims::Vararg{Integer,N}; kw...) where {T,N} =
  NDArray{T,N}(undef, dims; kw...)

NDArray(::UndefInitializer, dims::NTuple{N,Integer}; kw...) where {N} =
  NDArray{DEFAULT_DTYPE,N}(undef, dims; kw...)
NDArray(::UndefInitializer, dims::Vararg{Integer,N}; kw...) where {N} =
  NDArray{DEFAULT_DTYPE,N}(undef, dims; kw...)

NDArray(x::AbstractArray{<:DType}) = copy(collect(x), cpu())
NDArray(x::Array{<:DType})         = copy(x, cpu())

NDArray(::Type{T}, x::AbstractArray) where {T<:DType} =
  copy(convert(AbstractArray{T}, x), cpu())

NDArray(handle, writable = true) =
  NDArray{eltype(handle), ndims(handle)}(handle, writable)

# type aliases
const NDArrayOrReal = Union{NDArray,Real}
const VecOfNDArray = AbstractVector{<:NDArray}

Base.unsafe_convert(::Type{MX_handle}, x::NDArray) =
  Base.unsafe_convert(MX_handle, x.handle)
Base.convert(T::Type{MX_handle}, x::NDArray) = Base.unsafe_convert(T, x)
Base.cconvert(T::Type{MX_handle}, x::NDArray) = Base.unsafe_convert(T, x)

MX_handle(x::NDArray) = Base.convert(MX_handle, x)


text:
####
#### Parametric gates
####

"""
    Rphi(ϕ)

The phase shift gate ``R_ϕ``.
This is equivalent to qiskit's `u1` gate.
"""
function Rphi(ϕ)
    return [1 0; 0 cis(ϕ)]
end

# Could use complex(reverse(sincospi(z))...)
"""
    _exp_ipi(z)

Return `exp(im * pi * z)`. This uses the accurate `cospi` and `sinpi` functions.
"""
_exp_ipi(z) = cispi(z)

"""
    Rphipi(z)

Return the ``R_ϕ`` gate for ``ϕ = zπ``. This is more accurate
than `Rphi(pi * z)`.
"""
function Rphipi(z)
    return [1 0; 0 _exp_ipi(z)]
end

"""
    U2(ϕ, λ)

u2 gate. Need to put math definitions in here.
"""
function U2(ϕ, λ)
    return [1 -cis(λ); cis(ϕ) cis(ϕ + λ)] / sqrt(2)
end

# This does not appar any faster than computing each separately
# function RXalt(θ)
#     (s, c) = sincos(θ/2)
#     si = - im * s
#     return [c si; si c]
# end

function RX(θ)
    c = cos(θ/2)
    s = - im * sin(θ/2)
    return [c s; s c]
end

function RXpi(z)
    c = cospi(z/2)
    s = -im * sinpi(z/2)
    return [c s; s c]
end

function RY(θ)
    c = cos(θ/2)
    s = sin(θ/2)
    return [c -s; s c]
end

function RYpi(z)
    c = cospi(z/2)
    s = sinpi(z/2)
    return [c -s; s c]
end

function RZ(θ)
    return [cis(-θ/2) 0; 0 cis(θ/2)]
end

function RZpi(z)
    return [_exp_ipi(-z/2) 0; 0 _exp_ipi(z/2)]
end

"""
    Ualt(θ, ϕ, λ)

Matrix from SU(2).
Alternative parameterization.
"""
function Ualt(θ, ϕ, λ)
    c = cos(θ/2)
    s = sin(θ/2)
    fpl =  (ϕ + λ) / 2
    cfpl = cos(fpl)
    sfpl = sin(fpl)
    fml =  (ϕ - λ) / 2
    cfml = cos(fml)
    sfml = sin(fml)
    f00 = complex(cfpl, -sfpl) * c
    f01 = complex(-cfml, sfml) * s
    f10 = complex(cfml, sfml) * s
    f11 = complex(cfpl, sfpl) * c
    return [f00 f01; f10 f11]
end

"""
    U(θ, ϕ, λ)

Matrix from SU(2). This is
the same as qiskit's U or U3.
"""
function U(θ, ϕ, λ)
    c = cos(θ/2)
    s = sin(θ/2)
    fpl =  (ϕ + λ)
    cfpl = cos(fpl)
    sfpl = sin(fpl)
    f00 = c
    f01 = -s * complex(cos(λ), sin(λ))
    f10 = s * complex(cos(ϕ), sin(ϕ))
    f11 = c * complex(cfpl, sfpl)
    return [f00 f01; f10 f11]
end

"""
    Ualtpi(θ, ϕ, λ)

Matrix from SU(2), with `θ`, `ϕ`, and `λ` given as multiples of `π`.
This is more accurate than `Ualt`.
"""
function Ualtpi(θ, ϕ, λ)
    c = cospi(θ/2)
    s = sinpi(θ/2)
    fpl =  (ϕ + λ) / 2
    fml =  (ϕ - λ) / 2
    cfpl = cospi(fpl)
    sfpl = sinpi(fpl)
    cfml = cospi(fml)
    sfml = sinpi(fml)
    f00 = complex(cfpl, -sfpl) * c
    f01 = complex(-cfml, sfml) * s
    f10 = complex(cfml, sfml) * s
    f11 = complex(cfpl, sfpl) * c
    return [f00 f01; f10 f11]
end

@doc raw"""
    R(θ, ϕ)

The R gate.
```math
R(\theta, \phi) = e^{-i \frac{\theta}{2} (\cos{\phi} x + \sin{\phi} y)}
```
"""
function R(θ, ϕ)
    c = cos(θ/2)
    s = sin(θ/2)
    return [c -im*exp(-im*ϕ)*s;
            -im*exp(im*ϕ)*s c]
end

@doc raw"""
    Rpi(θ, ϕ)

The R gate with θ and ϕ reduced by π. This is may be more accurate,
for example with integral and half-integral multiples of π.
```math
Rpi(\theta, \phi) = e^{-i \frac{\pi\theta}{2} (\cos{\pi\phi} x + \sin{\pi\phi} y)}
```
"""
function Rpi(θ, ϕ)
    c = cospi(θ/2)
    s = sinpi(θ/2)
    return [c -im*cispi(-ϕ)*s;
            -im*cispi(ϕ)*s c]
end

"""
    RXXYY(θ)

The R_XXYY gate, `exp(-im θ / 2 (X⊗X + Y⊗Y))`.
"""
function RXXYY(θ)
    _RXXYY(θ, cos, sin)
end

"""
    RXXYYpi(θ)

The R_XXYY gate with `θ` scaled by `π`. `exp(-im * pi * θ / 2 * (X⊗X + Y⊗Y))`
"""
function RXXYYpi(θ)
    _RXXYY(θ, cospi, sinpi)
end

function _RXXYY(θ, cfunc, sfunc)
    II = I2 ⊗ I2 / 2
    XX = X ⊗ X / 2
    YY = Y ⊗ Y / 2
    ZZ = Z ⊗ Z / 2
    return II + ZZ + cfunc(θ) * (II - ZZ) -im * sfunc(θ) * (XX + YY)
end


text:
using Test,FermionicOperatorStringAnalyzer

@testset "FOSA.jl" begin
  @testset "$filename" for filename in [
    "FO.jl",
    "FOS.jl",
    "simplify.jl"
  ]
    println("Running $filename")
    include(filename)
  end
end


text:
# Test Inertia

using ModelingToolkit, OrdinaryDiffEq, CardiovascularModels

@variables t
D = Differential(t)

@named vessel1 = DampedElasticVessel(Ees = 5e6, Vd = 0.0, R_damp = 1e6)
@named vessel2 = DampedElasticVessel(Ees = 1e7, Vd = 0.0, R_damp = 1e6)

@named inertia = InertiaCompartment(Area = 0.4e-3, V = 0.2e-3)

eqs_con_inertia = [
    connect(vessel1.out, inertia.in),
    connect(inertia.out, vessel2.in),
    vessel1.in.Q ~ 0,
    vessel2.out.Q ~ 0
    ]

time_span = (0.0, 5.0)
@named connect_sys_inertia = ODESystem(eqs_con_inertia, t)
@named test_inertia = compose(connect_sys_inertia, vessel1, vessel2, inertia)
problem_inertia = ODEProblem(structural_simplify(test_inertia), volume_start, time_span, [])

sol_inertia = solve(problem_inertia, Rodas5(), dtmax = 0.01, reltol = 1e-6)

# Plots.plot(sol_inertia, vars = [vessel1.V, vessel2.V])


text:
"""
T: number of obvervations
M: number of joint prob
K: number of hidden states
"""
function likelihoods!(L::AbstractArray, hmm::AbstractHMM{Univariate}, observations)
    T, M, K = size(observations, 2), size(observations, 1), size(hmm, 1)
    @argcheck size(L) == (T, M, K)
    @inbounds for k in OneTo(K), t in OneTo(T), m in OneTo(M)
        L[t,m,k] = pdf(hmm.B[m, k], observations[m, t])
    end
end

function loglikelihoods!(LL::AbstractArray, hmm::AbstractHMM{Univariate}, observations)
    T, M, K = size(observations, 2), size(observations, 1), size(hmm, 1)
    @argcheck size(LL) == (T, M, K)
    @inbounds for k in OneTo(K), t in OneTo(T), m in OneTo(M)
        LL[t,m,k] = logpdf(hmm.B[m, k], observations[m, t])
    end
end

function likelihoods!(L::AbstractArray, hmm::AbstractHMM{Multivariate}, observations)
    T, M, K = size(observations, 2), size(observations, 1), size(hmm, 1)
    @argcheck size(L) == (T, M, K)
    @inbounds for k in OneTo(K), t in OneTo(T), m in OneTo(M)
        L[t,m,k] = pdf(hmm.B[m, k], view(observations, m, t))
    end
end

function loglikelihoods!(LL::AbstractArray, hmm::AbstractHMM{Multivariate}, observations)
    T, M, K = size(observations, 2), size(observations, 1), size(hmm, 1)
    @argcheck size(LL) == (T, M, K)
    @inbounds for k in OneTo(K), t in OneTo(T), m in OneTo(M)
        LL[t,m,k] = logpdf(hmm.B[m, k], view(observations, m, t))
    end
end

# function likelihoods(hmm::AbstractHMM, observations; logl = false, robust = false)
#     T, M, K = size(observations, 2), size(observations, 1), size(hmm, 1)
#     L = Array{Float64}(undef,T,M,K)#Matrix{Float64}(undef, T, M, K)
#
#     if logl
#         loglikelihoods!(L, hmm, observations)
#         robust && replace!(L, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))
#     else
#         likelihoods!(L, hmm, observations)
#         robust && replace!(L, -Inf => nextfloat(-Inf), Inf => prevfloat(Inf))
#     end
#
#     L
# end
# hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1) Normal(10,1); Normal(10,1) Normal(0,1)])
# y = rand(hmm, 1000)
# likelihoods(hmm, y)
# hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1) Normal(10,1); Normal(10,1) Normal(0,1)])
# obs = [rand(hmm,10);rand(hmm,10)];
# obs = reshape(obs,2,10);
# size(obs);
# likelihoods(hmm, obs)

#
# T, M, K = size(obs, 2), size(obs, 1), size(hmm, 1)
# observations = obs
# L = Array{Float64}(undef,T,M,K)
# observations[1, 2]
# hmm.B[1,2]
#
# L[1,1,1] = pdf(hmm.B[1, 1], observations[1, 1])
#
# for k in OneTo(K), t in OneTo(T), m in OneTo(M)
#     L[t,m,k] = pdf(hmm.B[m, k], observations[m, t])
# end

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase


text:
module MeasureEvalMethods

# Import the necessary packages.
import FastGaussQuadrature
import JuMP
import Distributions
const JuMPC = JuMP.Containers
using ..InfiniteOpt

# include jl files here
include("methods.jl")

end # end module


text:
module test_40_Warnings

using  ModiaPlot

include("test_base/test_40_Warnings.jl")

end


text:
using Statistics

include("update.jl")

function randomwalk(num::Integer, nsteps::Integer, prob::Real=0.5)
    means, vars = zeros(nsteps+1), zeros(nsteps+1)

    xs = zeros(num)
    means[1] = mean(xs)
    vars[1] = var(xs)
    for t in 1:nsteps
        map!(x -> update(x, prob), xs, xs)
        means[t+1] = mean(xs)
        vars[t+1] = var(xs)
    end
    return means, vars
end

function main(num, nsteps, prob, filename)
    if isempty(filename)
        io = stdout
    else
        io = open(filename, "w")
    end

    means, vars = randomwalk(num, nsteps, prob)

    for (i, (m, v)) in enumerate(zip(means, vars))
        println(io, "$(i-1) $m $v")
    end

    if io != stdout
        close(io)
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    num = length(ARGS)>0 ? parse(Int, ARGS[1]) : 1000
    nsteps = length(ARGS)>1 ? parse(Int, ARGS[2]) : 1000
    prob = length(ARGS)>2 ? parse(Float64, ARGS[3]) : 0.5
    filename = length(ARGS)>3 ? ARGS[4] : ""

    main(num, nsteps, prob, filename)
end

max_stars_repo_name: JuliaTagBot/JuliaBook-Samples

max_issues_repo_name: JuliaTagBot/JuliaBook-Samples

max_forks_repo_name: JuliaTagBot/JuliaBook-Samples


text:
function dedup_names(ns)
    count = Dict{Symbol,Int}()
    for n in ns
        if haskey(count, n)
            count[n] += 1
        else
            count[n] = 1
        end
    end

    repeated = filter(((k,v),) -> v > 1, count)
    for k in keys(repeated)
        repeated[k] = 0
    end
    [haskey(repeated, n) ? Symbol(n, "_", repeated[n]+=1) : n for n in ns]
end

function mapslices(f, x::NDSparse, dims; name = nothing)
    iterdims = setdiff([1:ndims(x);], map(d->fieldindex(columns(x.index),d), dims))
    idx = Any[Colon() for v in columns(x.index)]

    iter = Columns(astuple(columns(x.index))[[iterdims...]])
    if !isempty(dims) || !issorted(iter)
        iter = sort(iter)
    end

    for j in 1:length(iterdims)
        d = iterdims[j]
        idx[d] = iter[1][j]
    end

    if isempty(dims)
        idx[end] = vcat(idx[end])
    end

    y = f(x[idx...]) # Apply on first slice

    if isa(y, NDSparse)
        # this means we need to concatenate outputs into a big NDSparse
        ns = vcat(collect(dimlabels(x)[iterdims]), collect(dimlabels(y)))
        if !all(x->isa(x, Symbol), ns)
            ns = nothing
        else
            ns = dedup_names(ns)
        end
        n = length(y)
        index_first = similar(iter, n)
        for j=1:n
            @inbounds index_first[j] = iter[1]
        end
        index = Columns((columns(index_first)..., astuple(columns(copy(y.index)))...); names=ns)
        data = copy(y.data)
        output = NDSparse(index, data)
        if isempty(dims)
            _mapslices_itable_singleton!(f, output, x, 2)
        else
            _mapslices_itable!(f, output, x, iter, iterdims, 2)
        end
    else
        ns = dimlabels(x)[iterdims]
        if !all(x->isa(x, Symbol), ns)
            ns = nothing
        end
        index = Columns(Tuple(columns(iter[1:1])); names=ns)
        if isa(y, Tup)
            vec = convert(Columns, [y])
        else
            vec = [y]
        end
        if name === nothing
            output = NDSparse(index, vec)
        else
            output = NDSparse(index, Columns(Tuple(columns(vec)), names=[name]))
        end
        if isempty(dims)
            error("calling mapslices with no dimensions and scalar return value -- use map instead")
        else
            _mapslices_scalar!(f, output, x, iter, iterdims, 2, name!==nothing ? x->(x,) : identity)
        end
    end
end

function _mapslices_scalar!(f, output, x, iter, iterdims, start, coerce)
    idx = Any[Colon() for v in columns(x.index)]

    for i = start:length(iter)
        if i != 1 && roweq(iter, i-1, i) # We've already visited this slice
            continue
        end
        for j in 1:length(iterdims)
            d = iterdims[j]
            idx[d] = iter[i][j]
        end
        if length(idx) == length(iterdims)
            idx[end] = vcat(idx[end])
        end
        y = f(x[idx...])
        push!(output.index, iter[i])
        push!(output.data, coerce(y))
    end
    output
end

function _mapslices_itable_singleton!(f, output, x, start)
    I = output.index
    D = output.data

    I1 = Columns(columns(I)[1:ndims(x)])
    I2 = Columns(columns(I)[ndims(x)+1:end])
    i = start
    for i in start:length(x)
        k = x.index[i]
        y = f(NDSparse(x.index[i:i], x.data[i:i]))
        n = length(y)

        foreach((x,y)->append_n!(x,y,n), columns(I1), k)
        append!(I2, y.index)
        append!(D, y.data)
    end
    NDSparse(I,D)
end

function _mapslices_itable!(f, output, x, iter, iterdims, start)
    idx = Any[Colon() for v in columns(x.index)]
    I = output.index
    D = output.data
    initdims = length(iterdims)

    I1 = Columns(getsubfields(columns(I), 1:initdims)) # filled from existing table
    I2 = Columns(getsubfields(columns(I), initdims+1:fieldcount(typeof(columns(I))))) # filled from output tables

    for i = start:length(iter)
        if i != 1 && roweq(iter, i-1, i) # We've already visited this slice
            continue
        end
        for j in 1:length(iterdims)
            d = iterdims[j]
            idx[d] = iter[i][j]
        end
        if length(idx) == length(iterdims)
            idx[end] = vcat(idx[end])
        end
        subtable = x[idx...]
        y = f(subtable)
        n = length(y)

        foreach((x,y)->append_n!(x,y,n), columns(I1), iter[i])
        append!(I2, y.index)
        append!(D, y.data)
    end
    NDSparse(I,D)
end

function _flatten!(others, vecvec, out_others, out_vecvec)
    for i in 1:length(others)
        vec = vecvec[i]
        try
            for x in vec
                push!(out_vecvec, x)
                pushrow!(out_others, others, i)
            end
        catch err
            # FIXME: this is slow!
            if err isa MethodError && err.f === iterate
                push!(out_vecvec, vec)
                pushrow!(out_others, others, i)
            end
        end
    end
end

"""
    flatten(t::Table, col=length(columns(t)))

Flatten `col` column which may contain a vector of vectors while repeating the other fields.
If column argument is not provided, default to last column.

# Examples:

    x = table([1,2], [[3,4], [5,6]], names=[:x, :y])
    flatten(x, 2)

    t1 = table([3,4],[5,6], names=[:a,:b])
    t2 = table([7,8], [9,10], names=[:a,:b])
    x = table([1,2], [t1, t2], names=[:x, :y]);
    flatten(x, :y)
"""
function flatten(t::IndexedTable, col=length(columns(t)); pkey=nothing)
    vecvec = rows(t, col)
    hasmethod(iterate, (eltype(vecvec),)) || return t
    everythingbut = excludecols(t, col)

    order_others = Int[colindex(t, everythingbut)...]
    order_vecvec = Int[colindex(t, col)...]

    others = rows(t, everythingbut)
    out_others = similar(others, 0)
    out_vecvec = similar(arrayof(eltype(eltype(vecvec))), 0)

    _flatten!(others, vecvec, out_others, out_vecvec)

    cols = Any[columns(out_others)...]
    cs = columns(out_vecvec)
    newcols = isa(cs, Tup) ? Any[cs...] : Any[cs]
    ns = colnames(out_vecvec)
    i = colindex(t, col)
    cns = convert(Array{Any}, collect(colnames(t)))
    if length(ns) == 1 && !(ns[1] isa Symbol)
        ns = [colname(t, col)]
    end
    deleteat!(cns, i)
    for (n,c) in zip(reverse(ns), reverse(newcols))
        insert!(cns, i, n)
        insert!(cols, i, c)
    end
    if pkey === nothing
        if all(p -> p in order_others, t.pkey)
            pkey = t.pkey
        else
            pkey = []
        end
    end
    table(cols...; names=cns, pkey=pkey)
end


text:
module Peaks

using Compat

export argmaxima, argminima, maxima, minima, findmaxima, findminima, findnextmaxima,
    findnextminima, peakproms, peakproms!, peakwidths, peakwidths!

include("minmax.jl")
include("peakprom.jl")
include("peakwidth.jl")

end # module Peaks

max_stars_repo_name: halleysfifthinc/Peaks

max_issues_repo_name: halleysfifthinc/Peaks

max_forks_repo_name: halleysfifthinc/Peaks


text:
using Plots
using DelimitedFiles
using LinearAlgebra
using ProgressMeter
using LsqFit
using KernelDensity
using Statistics
using LaTeXStrings
using StatsBase
using Measures

image_type = "png"

using ArgParse
function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--pdf"
            action = :store_true
        "--svg"
            action = :store_true
    end
    return parse_args(s)
end

args = parse_commandline()
@info args

if args["pdf"]
    image_type="pdf"
elseif args["svg"]
    image_type="svg"
elseif args["pdf"] && args["svg"]
    @info "Multiple filetypes, defaulting to png"
end

@info "Saving plots as .$(image_type)"

beetle_length=13.0
beetle_width=7.0
dt = 1.0/30.0

include("AlphaShapes/AlphaShapes.jl")
include("analysis-utils.jl")

# using CSV rather than Julia specific JLD format
X50 = readdlm("../Data/50-X.csv",',')
X50 = zeros(size(X50)...,3)
X50[:,:,1] = readdlm("../Data/50-X.csv",',')
X50[:,:,2] = readdlm("../Data/50-Y.csv",',')
X50[:,:,3] = readdlm("../Data/50-theta.csv",',')

X100 = readdlm("../Data/100-X.csv",',')
X100 = zeros(size(X100)...,3)
X100[:,:,1] = readdlm("../Data/100-X.csv",',')
X100[:,:,2] = readdlm("../Data/100-Y.csv",',')
X100[:,:,3] = readdlm("../Data/100-theta.csv",',')

X200 = readdlm("../Data/200-X.csv",',')
X200 = zeros(size(X200)...,3)
X200[:,:,1] = readdlm("../Data/200-X.csv",',')
X200[:,:,2] = readdlm("../Data/200-Y.csv",',')
X200[:,:,3] = readdlm("../Data/200-theta.csv",',');

# collect speed and density
@info "collecting speed and density distributions"
D = []
V = []
@showprogress for X in [X50,X100,X200]
    d = zeros(size(X,1),size(X,2))
    v = zeros(size(X,1)-1,size(X,2))
    for i in 1:size(X,1)
        d[i,:] = Density(X[i,:,1:2] ./ beetle_length)
        if 1 < i < size(X,1)-1
            for j in 1:size(X,2)
                # central difference about i
                v[i,j] = norm(X[i+1,j,1:2] .- X[i-1,j,1:2])./(2.0*beetle_length*dt)
            end
        end
    end
    push!(D,d)
    push!(V,v)
end

names = ["50","100","200"]
linestyles = [:dash,:solid,:dot]
plot(layout=3,dpi=600,size=(900,900))
Xs = [X50,X100,X200]
for i in 1:size(Xs,1)
    X = Xs[i]
    bx = Float64.(logspace(1.0e-4,3.0,100))
    # bin across time and find avg and std for density
    # pdf across time bins
    pdfs = TimeBinDensity(bx,100,X,length_scale=beetle_length)
    pdf = mean(pdfs,dims=1)[1,:]
    pdf_e = std(pdfs,dims=1)[1,:]
    plot!(bx,pdf,ribbon=pdf_e,label="N=$(names[i])",linestyle=linestyles[i],linewidth=3,subplot=i)
    # effectively 0.0 at this point
    xlims!(0.0,3.0)
    #ylims!(0.0,5.0)
end
plot!()
ylabel!("PDF")
xlabel!("Local Density (square body lengths)")
savefig("density-distributions-time-errors.$(image_type)")

names = ["50","100","200"]
linestyles = [:dash,:solid,:dot]
plot(dpi=600)
for i in 1:size(D,1)
    d = D[i][isnan.(D[i]).==false];
    pdf = kde(d)
    plot!(pdf.x,pdf.density,label="N=$(names[i])",linestyle=linestyles[i],linewidth=3)
end
plot!()
# effectively 0.0 at this point
xlims!(0,3.0)
ylabel!("Kernel Density Estimate of PDF")
xlabel!("Local Density (square body lengths)")
savefig("density-distributions.$(image_type)")

names = ["50","100","200"]
linestyles = [:dash,:solid,:dot]
plot(layout=(1,3),dpi=150,size=(1200,300))
for i in 1:size(V,1)
    v = V[i][isnan.(V[i]).==false]
    speed_bins=linspace(0.0,60.0,100)
    pdfs = TimeBinDistribution(speed_bins,v;bins=100)
    pdf = mean(pdfs,dims=1)[1,:]
    pdf_e = std(pdfs,dims=1)[1,:]
    plot!(speed_bins,pdf,ribbon=pdf_e,label="N=$(names[i])",linestyle=linestyles[i],linewidth=3,subplot=i,bottommargin=5mm,leftmargin=5mm)
    ylims!(0.0,maximum(pdf)+0.5,subplot=i)
end
plot!()
ylabel!("PDF")
xlabel!("Speed (body lengths per second)")
savefig("speed-distributions.$(image_type)")

# speed and density relationship data
@info "Computing speed and density relationship"

DV50 = []
DV100 = []
DV200 = []

DV = [DV50,DV100,DV200]
@showprogress for (i,X) in enumerate([X50,X100,X200])
    for t in 1:size(X,1)-1
        for b in 1:size(X,2)
            push!(DV[i],[D[i][t,b],V[i][t,b],X[t,b,3]])
        end
    end
end

@. model(x,p) = p[1].*x.^(-p[2]).+p[3]
n = ["50","100","200"]
Data = []
fits = []
for (j,DV) in enumerate([DV50,DV100,DV200])
    d = []
    v = []
    o = []
    for i in 1:length(DV)
        push!(d,DV[i][1])
        push!(v,DV[i][2])
        push!(o,DV[i][3])
    end

    ind = isinf.(d).==false
    d = d[ind]
    v = v[ind]
    o = o[ind]

    ind = (d) .< Inf
    # scale both back to pixels for plot function
    d = d[ind] ./ beetle_length^2.0
    v = v[ind] .* (2.0*dt*beetle_length)
    o = o[ind]

    #println(maximum(d))

    # generate a binned plot with errorbars
    # 1.0/(2dt) because of the central difference!
    local D,V,SE=DVBinPlot(D=d,V=v,n_bins=25,beetlescale=true,beetle_length=beetle_length,xscale=:log10,FPS=1.0/(2.0*dt),errors=true)
    inds = isnan.(D) .== false
    D = D[inds]
    V = V[inds]
    inds = isnan.(V) .== false
    D = D[inds]
    V = V[inds]

    inds = isnan.(o) .== false
    d = d[inds]
    o = o[inds]
    d_,o_,seo_ = DVBinPlot(D=d,V=o,n_bins=25,beetlescale=false,errors=true)
    push!(Data,(D,V,SE,d_,o_,seo_))
    fit = curve_fit(model, D, V, [1.,1.,1.])
    @info "power law fit converged: ", fit.converged, " N = $(n[j])"
    push!(fits,fit.param)
end

replace!(Data[2][3],NaN=>0.0)
plot(dpi=600)
labels=[L"N=50",L"N=100",L"N=200"]
colours = ["red","green","blue"]
shapes = [:square, :diamond, :circle]
for i in 1:length(Data)
    scatter!(Data[i][1],Data[i][2],yerr=Data[i][3],label=labels[i],c=colours[i],mark=shapes[i])
end
plot!(xaxis=:log10,yaxis=:log10)
xlabel!(L"\rho")
ylabel!(L"v(\rho)")

ρ = logspace(1e-2,1,100)
plot!(ρ,10*ρ.^-0.4,label=L"\propto\rho^{-0.4}",c="black")
savefig("speed-density-relationship.$(image_type)")

# Diffusion
@info "Diffusion statistics"
T = copy(X200)

# unwravel phase
for i in 1:size(T,2)
    T[:,i,3] = unwrap(T[:,i,3])
end

MSAD,t = msdθ(T,size(T,1),800)
MSAD = MSAD[:,1]
@. model(t,p) = p[1]*t
ind = findfirst(x->x.>2.5,t*dt)
fit_ = curve_fit(model,t[ind:end].*dt,MSAD[ind:end],[1.0])
fit_.converged

k = 15

x = t[1:k:end]*dt
y = MSAD[1:k:end]
xt = t[ind:end]*dt
yt = model(t[ind:end]*dt,fit_.param)

xt = xt[1:k:end]
yt = yt[1:k:end]

scatter(x,y,label="",legend=:topleft,title="")
plot!(xt,yt,label="Fit to 2.0*Dr*t, Dr = $(round(fit_.param[1]/2.0,digits=3))",linewidth=2)
xlabel!("Seconds")
ylabel!("Mean Square Angular Displacement (rad^2 s)")

T = copy(X200)
Dr = 2.341
taur = 1.0 / Dr

MSAD,t = MSD(T.*(5e-3/beetle_length),size(T,1),30*1)
MSAD = MSAD[:,1];

@. model(t,p) = (4.0*p[1]+(2.0*p[2]^2.0)*taur)*t + (2.0*p[2]^2.0)*(taur^2.0)*(exp(-t/taur)-1)
ind = 1
fit_ = curve_fit(model,t[ind:end].*dt,MSAD[ind:end],[1.0,1.0],lower=[0.0,5e-3], upper=[Inf,Inf])
fit_.converged

scatter(t*dt,MSAD,label="",legend=:topleft,title="",dpi=600)
plot!(t[ind:end]*dt,model(t[ind:end]*dt,fit_.param),label="Dt = $(round(fit_.param[1],digits=5)), v0 = $(round(fit_.param[2],digits=5))")
xlabel!("Seconds")
ylabel!("Mean Square Displacement (m^2 s)")
savefig("translational-diffusion.$(image_type)")

@info "Example of delay between velocity and orientation"
b = 3
t_max = 20
K = 5
X = X50
plot(aspect_ratio=:equal,grid=false,dpi=600)
plot!(xaxis=nothing,yaxis=nothing)
for t in 2:t_max
    H = VelocityHeading(X,t,b)
    o = [cos(X[t,b,3]),sin(X[t,b,3])]
    Arrow(H,X[t,b,1:2],K,"black")
    Arrow(o,X[t,b,1:2],K,"red")
    E = Ellipse(1000,a=13.0/2.0,b=7.0/2.0,θ=X[t,b,3],center=X[t,b,1:2])
    s = Shape(E[:,1],E[:,2])
    plot!(s,aspect_ratio=:equal,c="blue",alpha=0.25,label="")
end
plot!([],[],label="Body Orientation",c="red")
plot!([],[],label="Velocity Direction",c="Black")
plot!([],[],label="Template Indicative of Beetle Outline",c="blue",alpha=0.25)
savefig("delay-example.$(image_type)")

### CollisonFreeAnalysis

@info "Collision Free Statistics"

data = [X200]
C = Collisions(beetle_length,data) # collision filtered
T,L,O = CollisionFreeTrajectories(C[1],X200)

histogram(L,label="",dpi=600)
savefig("collision-free-trajectory-length-distribution.$(image_type)")

function thresh(x;l=minimum(x),u=maximum(x),dx=0.01)
    X = collect(l:dx:u)
    c = zeros(length(X))
    for (i,xx) in enumerate(X)
        c[i] = sum(x.>X[i])
    end
    return X,c
end

x,y = thresh(L.*1.0/30.0,dx=0.1)
plot(x,y,label="")
xlabel!("Minimum Collision Free Trajectory Length (seconds)")
ylabel!("Number of Trajectories")
savefig("collision-free-trajectory-count.$(image_type)")

@info "Gathering orientation Velocity Correlation"

C,t = NVCorr(O,1.0/30.0,4,2);
c = length.(C)
e = std.(C)
Cm = mean.(C)
t = t[2:end]
Cm = Cm[2:end]

@. model(t,p) = p[1].*exp.(-0.5.*((t.-p[2])^2.0)./p[3]^2.0)
fit_ = curve_fit(model,t,Cm,[1.0,1.0,1.0])
@info  fit_.converged,fit_.param
@info "Gaussian fit peak location ", fit_.param[2]
@info "MSE ", mean((Cm.-model(t,fit_.param)).^2.0)
tt = collect(minimum(t):0.001:maximum(t))

p1 = scatter(t,Cm,yerr=e./sqrt.(c),label="")
plot!(tt,model(tt,fit_.param),label="Gaussian Fit",linewidth=3,xtickfontsize=18,ytickfontsize=18)
savefig("velocity-orientation-correlation-short-times.$(image_type)")

C,t = NVCorr(O,1.0/30.0,30,2);
c = length.(C)
e = std.(C)
Cm = mean.(C)
t = t[2:end]
Cm = Cm[2:end]

p2 = scatter(t,Cm,yerr=e./sqrt.(c),label="")
plot!(tt,model(tt,fit_.param),label="Gaussian Fit",linewidth=3)
savefig("velocity-orientation-correlation.$(image_type)")

data = [X50,X100,X200]
C = Collisions(beetle_length,data)
for j in 1:3
    local t,l,o = CollisionFreeTrajectories(C[j],data[j])
    local S = []
    @showprogress for i in 1:length(o)
        d = o[i]
        if length(d) >= 3
            for t in 2:size(d,1)-1
                v = d[t+1][1:2] .- d[t-1][1:2]
                if isnan(norm(v)) == false
                    push!(S,norm(v)/(2*beetle_length*dt))
                end
            end
        end
    end
    @info "N = $(names[j]): Mean Collision Free Speed", mean(S), " body lengths per second"
end

max_stars_repo_name: harveydevereux/CUDA-Whirligigs

max_issues_repo_name: harveydevereux/CUDA-Whirligigs

max_forks_repo_name: harveydevereux/CUDA-Whirligigs


text:
function gettypes(
    ::Type{T}, # number type
    ::Type{Val{:Static}}, # matrix type
    ::Type{Val{Kesize}}, # matrix size
) where {T,Kesize}
    return SMatrix{Kesize,Kesize,T,Kesize^2}, SVector{Kesize,T}
end
function gettypes(
    ::Type{T}, # number type
    ::Type{Val{:SMatrix}}, # matrix type
    ::Type{Val{Kesize}}, # matrix size
) where {T,Kesize}
    return SMatrix{Kesize,Kesize,T,Kesize^2}, SVector{Kesize,T}
end
function gettypes(
    ::Type{T}, # number type
    ::Type{Val{:MMatrix}}, # matrix type
    ::Type{Val{Kesize}}, # matrix size
) where {T,Kesize}
    return MMatrix{Kesize,Kesize,T,Kesize^2}, MVector{Kesize,T}
end
function gettypes(
    ::Type{BigFloat}, # number type
    ::Type{Val{:Static}}, # matrix type
    ::Type{Val{Kesize}}, # matrix size
) where {Kesize}
    return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat}
end
function gettypes(
    ::Type{BigFloat}, # number type
    ::Type{Val{:SMatrix}}, # matrix type
    ::Type{Val{Kesize}}, # matrix size
) where {Kesize}
    return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat}
end
function gettypes(
    ::Type{BigFloat}, # number type
    ::Type{Val{:MMatrix}}, # matrix type
    ::Type{Val{Kesize}}, # matrix size
) where {Kesize}
    return SizedMatrix{Kesize,Kesize,BigFloat,Kesize^2}, SizedVector{Kesize,BigFloat}
end
function gettypes(
    ::Type{T}, # number type
    ::Any, # matrix type
    ::Any, # matrix size
) where {T}
    return Matrix{T}, Vector{T}
end

initialize_K(sp::StiffnessTopOptProblem) = Symmetric(create_sparsity_pattern(sp.ch.dh))

initialize_f(sp::StiffnessTopOptProblem{dim,T}) where {dim,T} = zeros(T, ndofs(sp.ch.dh))

function make_Kes_and_fes(problem, quad_order=2)
    return make_Kes_and_fes(problem, quad_order, Val{:Static})
end
function make_Kes_and_fes(problem, ::Type{Val{mat_type}}) where {mat_type}
    return make_Kes_and_fes(problem, 2, Val{mat_type})
end
function make_Kes_and_fes(problem, quad_order, ::Type{Val{mat_type}}) where {mat_type}
    T = floattype(problem)
    dim = getdim(problem)
    geom_order = getgeomorder(problem)
    dh = getdh(problem)
    E = getE(problem)
    ν = getν(problem)
    ρ = getdensity(problem)

    refshape = Ferrite.getrefshape(dh.field_interpolations[1])

    λ = E * ν / ((1 + ν) * (1 - 2 * ν))
    μ = E / (2 * (1 + ν))
    δ(i, j) = i == j ? T(1) : T(0)
    g(i, j, k, l) = λ * δ(i, j) * δ(k, l) + μ * (δ(i, k) * δ(j, l) + δ(i, l) * δ(j, k))
    C = SymmetricTensor{4,dim}(g)

    # Shape functions and quadrature rule
    interpolation_space = Lagrange{dim,refshape,geom_order}()
    quadrature_rule = QuadratureRule{dim,refshape}(quad_order)
    cellvalues = CellScalarValues(quadrature_rule, interpolation_space)
    facevalues = FaceScalarValues(
        QuadratureRule{dim - 1,refshape}(quad_order), interpolation_space
    )

    # Calculate element stiffness matrices
    n_basefuncs = getnbasefunctions(cellvalues)

    Kesize = dim * n_basefuncs
    MatrixType, VectorType = gettypes(T, Val{mat_type}, Val{Kesize})
    Kes, weights = _make_Kes_and_weights(
        dh,
        Tuple{MatrixType,VectorType},
        Val{n_basefuncs},
        Val{dim * n_basefuncs},
        C,
        ρ,
        quadrature_rule,
        cellvalues,
    )
    dloads = _make_dloads(weights, problem, facevalues)

    return Kes, weights, dloads, cellvalues, facevalues
end

const g = [0.0, 9.81, 0.0] # N/kg or m/s^2

# Element stiffness matrices are StaticArrays
# `weights` : a vector of `xdim` vectors, element_id => self-weight load vector
function _make_Kes_and_weights(
    dh::DofHandler{dim,N,T},
    ::Type{Tuple{MatrixType,VectorType}},
    ::Type{Val{n_basefuncs}},
    ::Type{Val{Kesize}},
    C,
    ρ,
    quadrature_rule,
    cellvalues,
) where {dim,N,T,MatrixType<:StaticArray,VectorType,n_basefuncs,Kesize}
    # Calculate element stiffness matrices
    nel = getncells(dh.grid)
    body_force = ρ .* g # Force per unit volume
    Kes = Symmetric{T,MatrixType}[]
    sizehint!(Kes, nel)
    weights = [zeros(VectorType) for i in 1:nel]
    Ke_e = zeros(T, dim, dim)
    fe = zeros(T, Kesize)
    Ke_0 = Matrix{T}(undef, Kesize, Kesize)
    celliterator = CellIterator(dh)
    for (k, cell) in enumerate(celliterator)
        Ke_0 .= 0
        reinit!(cellvalues, cell)
        fe = weights[k]
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for b in 1:n_basefuncs
                ∇ϕb = shape_gradient(cellvalues, q_point, b)
                ϕb = shape_value(cellvalues, q_point, b)
                for d2 in 1:dim
                    fe = @set fe[(b - 1) * dim + d2] += ϕb * body_force[d2] * dΩ
                    for a in 1:n_basefuncs
                        ∇ϕa = shape_gradient(cellvalues, q_point, a)
                        Ke_e .= dotdot(∇ϕa, C, ∇ϕb) * dΩ
                        for d1 in 1:dim
                            #if dim*(b-1) + d2 >= dim*(a-1) + d1
                            Ke_0[dim * (a - 1) + d1, dim * (b - 1) + d2] += Ke_e[d1, d2]
                            #end
                        end
                    end
                end
            end
        end
        weights[k] = fe
        if MatrixType <: SizedMatrix # Work around because full constructor errors
            push!(Kes, Symmetric(SizedMatrix{Kesize,Kesize,T}(Ke_0)))
        else
            push!(Kes, Symmetric(MatrixType(Ke_0)))
        end
    end
    return Kes, weights
end
# Fallback
function _make_Kes_and_weights(
    dh::DofHandler{dim,N,T},
    ::Type{Tuple{MatrixType,VectorType}},
    ::Type{Val{n_basefuncs}},
    ::Type{Val{Kesize}},
    C,
    ρ,
    quadrature_rule,
    cellvalues,
) where {dim,N,T,MatrixType,VectorType,n_basefuncs,Kesize}
    # Calculate element stiffness matrices
    nel = getncells(dh.grid)
    body_force = ρ .* g # Force per unit volume
    Kes = let Kesize = Kesize, nel = nel
        [Symmetric(zeros(T, Kesize, Kesize), :U) for i in 1:nel]
    end
    weights = let Kesize = Kesize, nel = nel
        [zeros(T, Kesize) for i in 1:nel]
    end
    Ke_e = zeros(T, dim, dim)
    celliterator = CellIterator(dh)
    for (k, cell) in enumerate(celliterator)
        reinit!(cellvalues, cell)
        fe = weights[k]
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for b in 1:n_basefuncs
                ∇ϕb = shape_gradient(cellvalues, q_point, b)
                ϕb = shape_value(cellvalues, q_point, b)
                for d2 in 1:dim
                    fe[(b - 1) * dim + d2] += ϕb * body_force[d2] * dΩ
                    for a in 1:n_basefuncs
                        ∇ϕa = shape_gradient(cellvalues, q_point, a)
                        Ke_e .= dotdot(∇ϕa, C, ∇ϕb) * dΩ
                        for d1 in 1:dim
                            #if dim*(b-1) + d2 >= dim*(a-1) + d1
                            Kes[k].data[dim * (a - 1) + d1, dim * (b - 1) + d2] += Ke_e[
                                d1, d2
                            ]
                            #end
                        end
                    end
                end
            end
        end
    end
    return Kes, weights
end

"""
    _make_dload(problem)

Assemble a sparse vector for boundary (face) distributed loads
"""
function _make_dloads(fes, problem, facevalues)
    dim = getdim(problem)
    N = nnodespercell(problem)
    T = floattype(problem)
    dloads = deepcopy(fes)
    for i in 1:length(dloads)
        if eltype(dloads) <: SArray
            dloads[i] = zero(eltype(dloads))
        else
            dloads[i] .= 0
        end
    end
    pressuredict = getpressuredict(problem)
    dh = getdh(problem)
    grid = dh.grid
    boundary_matrix = grid.boundary_matrix
    cell_coords = zeros(Ferrite.Vec{dim,T}, N)
    n_basefuncs = getnbasefunctions(facevalues)
    for k in keys(pressuredict)
        t = -pressuredict[k] # traction = negative the pressure
        faceset = getfacesets(problem)[k]
        for (cellid, faceid) in faceset
            boundary_matrix[faceid, cellid] ||
                throw("Face $((cellid, faceid)) not on boundary.")
            fe = dloads[cellid]
            getcoordinates!(cell_coords, grid, cellid)
            reinit!(facevalues, cell_coords, faceid)
            for q_point in 1:getnquadpoints(facevalues)
                dΓ = getdetJdV(facevalues, q_point) # Face area
                normal = getnormal(facevalues, q_point) # Nomral vector at quad point
                for i in 1:n_basefuncs
                    ϕ = shape_value(facevalues, q_point, i) # Shape function value
                    for d in 1:dim
                        if fe isa SArray
                            fe = @set fe[(i - 1) * dim + d] += ϕ * t * normal[d] * dΓ
                        else
                            fe[(i - 1) * dim + d] += ϕ * t * normal[d] * dΓ
                        end
                    end
                end
            end
            dloads[cellid] = fe
        end
    end

    return dloads
end

"""
    make_cload(problem)

Assemble a sparse vector for concentrated loads
"""
function make_cload(problem)
    T = floattype(problem)
    dim = getdim(problem)
    cloads = getcloaddict(problem)
    dh = getdh(problem)
    metadata = getmetadata(problem)
    node_dofs = metadata.node_dofs
    inds = Int[]
    vals = T[]
    for nodeidx in keys(cloads)
        for (dofidx, force) in enumerate(cloads[nodeidx])
            if force != 0
                dof = node_dofs[(nodeidx - 1) * dim + dofidx]
                push!(inds, dof)
                push!(vals, force)
            end
        end
    end
    return sparsevec(inds, vals, ndofs(dh))
end


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

name = "LEAN_Community"
version = v"3.15.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/leanprover-community/lean.git", "56f8877f1efa22215aca0b82f1c0ce2ff975b9c3"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/lean

for patch_file in ../patches/*.patch; do
    atomic_patch -p1 "${patch_file}"
done

mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DBUILD_TESTING=OFF \
    ../src
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"))
]

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


text:
#!/usr/bin/env julia

# This attempts to implement call/cc feature in julia

type Continuation_box{name}
    value
end

macro call_cc(name::Symbol, content)
    quote
        u = Symbol(randstring())
        $name = v -> throw(Continuation_thrown{u}(v))
        try
            $content
        catch e
            if e isa Continuation_box{u}
                e.value
            end
        end
    end
end

max_stars_repo_name: sclereid/collection

max_issues_repo_name: sclereid/collection

max_forks_repo_name: sclereid/collection


text:
using Test
using LinearAlgebra
using QCMaterial

import PyCall: pyimport
import Random

@testset "vqe.solve_gs" begin
    Random.seed!(100)

    nsite = 2
    n_qubit = 2*nsite
    U = 1.0
    t = -0.01
    μ = 0

    ham = FermionOperator()
    for i in 1:nsite
        up = up_index(i)
        down = down_index(i)
        ham += FermionOperator("$(up)^ $(down)^ $(up) $(down)", -U)
    end

    for i in 1:nsite-1
        ham += FermionOperator("$(up_index(i+1))^ $(up_index(i))", t)
        ham += FermionOperator("$(up_index(i))^ $(up_index(i+1))", t)
        ham += FermionOperator("$(down_index(i+1))^ $(down_index(i))", t)
        ham += FermionOperator("$(down_index(i))^ $(down_index(i+1))", t)
    end

    for i in 1:nsite
        up = up_index(i)
        down = down_index(i)
        ham += FermionOperator("$(up)^  $(up) ", -μ)
        ham += FermionOperator("$(down)^ $(down)", -μ)
    end

    # Compute exact ground-state energy
    n_electron = 2
    @assert mod(n_electron, 2) == 0
    sparse_mat = get_number_preserving_sparse_operator(ham, n_qubit, n_electron);
    EigVal_min = minimum(eigvals(sparse_mat.toarray()));

    # Ansatz
    state0 = create_hf_state(n_qubit, n_electron)
    circuit = uccgsd(n_qubit, orbital_rot=true)
    theta_init = rand(num_theta(circuit))

    # VQE
    cost_history, thetas_opt =
       QCMaterial.solve_gs(jordan_wigner(ham), circuit, state0, theta_init=theta_init, verbose=true,
           comm=QCMaterial.MPI_COMM_WORLD
       )

    @test abs(EigVal_min-cost_history[end]) < 1e-6
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
"""
    GraphsOptim

A package for graph-related optimization algorithms that rely on Linear Programming.

The package JuMP.jl and one of its supported solvers are required.
"""
module GraphsOptim

using Graphs
using JuMP
using SimpleTraits
using SparseArrays

include("utils.jl")
include("minimum_cost_flow.jl")
include("maximum_weight_matching.jl")

export minimum_cost_flow
export maximum_weight_matching, maximum_weight_maximal_matching

end


text:
println("Stats Tests")
println("...........")
using Lathe.stats: NormalDist, T_Dist, mean, std
# Distributions
@testset "Distributions" begin
x = randn(200)
dist = NormalDist(x)
mu = mean(x)
zero_std = dist.apply(mu)
@test Real(zero_std[1]) == 0
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_draggabledashboardresponsive

"""
    ''_draggabledashboardresponsive(;kwargs...)
    ''_draggabledashboardresponsive(children::Any;kwargs...)
    ''_draggabledashboardresponsive(children_maker::Function;kwargs...)

A DraggableDashboardResponsive component.
DraggableDashboard is a component for building
dashboards with draggable and resizable items.
It takes a list of children and display them in
div elements that can be moved around the page.
The initial size of each element can either be
defined with the layout argument or by wrapping
each element with the DashboardItem component.
By default, DraggableDashboard will saved the
position of the elements on client side, when
moved on the page. But you can also save it
on server side by defining a callback with :
`Input("<my-id>", "layout")`.
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional): Children is a list of the items (dash Components and/or
DashboardItem) to diplay on the layout.
By default all the items can be dragged and resized.
- `id` (String; optional): (string) The ID used to identify this component in Dash callbacks.
The id is also used to automatically save the layout on client side.
- `layouts` (Dict; optional): Layout is a list(python)/vector(R) of dictionnary(Python)/list(R) with the format:
{x: number, y: number, w: number, h: number}
The index into the layout must match the id used on each item component with DashboardItem.
If you choose to use custom keys, you can specify that key in the layout
array objects like so:
{i: string, x: number, y: number, w: number, h: number}
The ID used to identify this component in Dash callbacks.
The id is also used to automatically save the layout on client side.
- `breakpoints` (Dict; optional): ({breakpoint: number}) The breakpoints for the responsive layout.
For each screen size (breakpoint) we can define a different layout.
(see also 'layouts' and 'gridCols' arguments)
Default value is {lg: 1200, md: 996, sm: 768, xs: 480, xxs: 0}
- `gridCols` (Dict; optional): ({breakpoint: number}) the number of columns in the grid layout.
Default value is {lg: 12, md: 10, sm: 6, xs: 4, xxs: 2}.
- `save` (Bool; optional): (bool) If True, then the layout is automatically saved on client browser.
Default value is True
- `clearSavedLayout` (Bool; optional): (bool) If set to true, then the layout saved in the client browser
will be cleared in the next page load.
- `ncols` (Real; optional): ({breakpoint: number}) the default number of columns by item.
Default value is {lg: 6, md: 5, sm: 3, xs: 4, xxs: 2}.
- `nrows` (Real; optional): (number) the default number of row by item.
Default value is 8.
- `height` (Real; optional): (number) height of a row (in px).
Default value is 30.
- `className` (String; optional): (string) class passed to the react-grid-layout component
- `style` (Dict; optional): (dict) css style passed to the react-grid-layout component
"""
function ''_draggabledashboardresponsive(; kwargs...)
        available_props = Symbol[:children, :id, :layouts, :breakpoints, :gridCols, :save, :clearSavedLayout, :ncols, :nrows, :height, :className, :style]
        wild_props = Symbol[]
        return Component("''_draggabledashboardresponsive", "DraggableDashboardResponsive", "dash_draggable", available_props, wild_props; kwargs...)
end

''_draggabledashboardresponsive(children::Any; kwargs...) = ''_draggabledashboardresponsive(;kwargs..., children = children)
''_draggabledashboardresponsive(children_maker::Function; kwargs...) = ''_draggabledashboardresponsive(children_maker(); kwargs...)

max_stars_repo_name: tcgvn/dash-draggable

max_issues_repo_name: tcgvn/dash-draggable

max_forks_repo_name: tcgvn/dash-draggable


text:
mutable struct simulation_command_t <: LCMType
    timestamp::Int64
    command_type::Int8
    string_data::String
    num_float_data::Int32
    float_data::Vector{Float32}
end

@lcmtypesetup(simulation_command_t,
    float_data => (num_float_data,),
)


text:
function qsort(array)
    if  (isempty(array))
        return array
    else
        pivot = array[1]
        less = [i for i in array[2:end] if i < pivot]
        greater = [i for i in array[2:end] if i >= pivot]
        return vcat(qsort(less), pivot, qsort(greater))
    end
end

v = [144 89 4 9 95 12 86 25]
println("unsorted array = ",v)

sorted_v = qsort(v)

println("selection sorted array = ", sorted_v)

max_stars_repo_name: landrumb/julia

max_issues_repo_name: landrumb/julia

max_forks_repo_name: landrumb/julia


text:
eval(quote
    using .Symbolics: gradient,
                      simplify,
                      Num,  # variable like, e.g. x[1]
                      Term, # term like, eg. x[1] + x[2] == 1
                      Symbolic,
                      operation,
                      arguments

    """
        _vec(vars)

    Transform a tuple of operations into one vector of operations.

    ### Input

    - `vars` -- tuple where each element is either variable-like (`Num`) or a
                vector of variables (`Vector{Num}`)

    ### Output

    A vector of `Operation` obtained by concatenating each tuple component.

    ## Examples

    ```jldoctest
    julia> using Symbolics

    julia> vars = @variables x[1:2] y
    2-element Vector{Any}:
      x[1:2]
     y

    julia> LazySets._vec(vars)
    3-element Vector{Num}:
     x[1]
     x[2]
        y
    ```
    """
    function _vec end

    # reduce for several variables e.g. when vars = @variables x[1:3] t
    _vec(vars::Vector{Any}) = reduce(vcat, vars)
    _vec(vars::Vector{Num}) = vars
    _vec(vars::Vector{Symbolics.Arr{Num, 1}}) = reduce(vcat, vars)
    _vec(vars::Vector{Vector{Num}}) = reduce(vcat, vars)
    _vec(vars::Vector{Real}) = reduce(vcat, vars)

    _get_variables(expr::Num) = convert(Vector{Num}, Symbolics.get_variables(expr))
    _get_variables(expr::Vector{<:Num}) = unique(reduce(vcat, _get_variables(ex) for ex in expr))
end)

eval(load_symbolics_hyperplane())
eval(load_symbolics_halfspace())
eval(load_symbolics_hpolyhedron())
eval(load_symbolics_hpolytope())