text:
#=
examples
x => b => c  ==> b = m[1](x) ; c = m[2](b)
x => 3 ==> x => a => a => a ==> x = m[1](a); a = m[1](a); a = m[1](a)
(x, m) => a => b => c ==> a = m[1](x , m); b = m[2](b); c = m[3](b)
((x, m) => x) => 3 ==> (x = m[1](x, m)); (x = m[2](x, m)); (x = m[3](x, m))
(((x, m) => x:(x, m)) => 3) ==> (x = m[1](x,m)); (x = m[2](x,m)) ;(x = m[3](x,m))
=#

struct NNTopo{F}
    fs::String
    f::F
end

NNTopo(ex::Expr) = NNTopo(string(ex), tofunc(ex))
NNTopo(s::String) = NNTopo(s, tofunc(s))
(nt::NNTopo)(xs...) = nt.f(xs...)

macro nntopo_str(str)
    NNTopo(str)
end

macro nntopo(expr)
    NNTopo(string(expr), tofunc(expr))
end

tofunc(sf::String) = tofunc(Meta.parse(sf))
function tofunc(pattern)
    !islegal(pattern) && error("topo pattern illegal")

    m = gensym(:model)
    fname = gensym(:topo_func)

    code = to_code(pattern)

    if istuple(code.in)
        fname = Expr(:call, fname, m, code.in.args...)
    else
        fname = Expr(:call, fname, m, code.in)
    end

    fbody = Any[:block]
    for (i, l) ∈ enumerate(code.lines)
        push!(fbody, genline(l..., m, i))
    end

    push!(fbody, code.out)
    func = Expr(:function, fname, Expr(fbody...))

    eval(func)
end

function Base.show(io::IO, nt::NNTopo)
    println(io, "NNTopo{\"$(nt.fs)\"}")
    print_topo(io, nt)
    io
end

print_topo(nt::NNTopo; models=nothing) = print_topo(stdout, nt; models=models)
function print_topo(io::IO, nt::NNTopo; models=nothing)
    code = to_code(Meta.parse(nt.fs))
    farg = istuple(code.in) ? join(code.in.args, ", ") : string(code.in)
    println(io, "topo_func(model, $farg)")
    for (i, l) ∈ enumerate(code.lines)
        name = string(l[1])
        args = istuple(l[2]) ? string(l[2]) : "($(l[2]))"
        model = models === nothing ? "model[$i]" : string(models[i])
        println(io, "\t$name = $model$args")
    end
    println(io, "\t$(code.out)")
    println("end")
end

const legalsym = (:(=>),:→,:(:))

islegal(ex) = false
islegal(ex::Symbol) = true
islegal(ex::Int) = true
function islegal(ex::Expr)
    if istuple(ex)
        all(x->x isa Symbol, ex.args)
    elseif iscolon(ex)
        length(ex.args) == 3 && isdup(ex) ?
            getsym(ex) isa Symbol || istuple(getsym(ex)) :
            (ex.args[2] isa Symbol || istuple(ex.args[2])) &&
            (ex.args[3] isa Symbol || istuple(ex.args[3]))
    else
        ex.head == :call && ex.args[1] ∈ legalsym &&
            length(ex.args) == 3 && islegal(ex.args[2]) && islegal(ex.args[3])
    end
end

genline(name, arg::Symbol, m, i::Int) = Expr(:(=), name, Expr(:call, :($m[$i]), arg))
genline(name, args::Expr, m, i::Int) = Expr(:(=), name, Expr(:call, :($m[$i]), args.args...))

struct Code
    in
    out
    lines::Vector{Tuple}
end

Code(in, out) = Code(in, out, [(out, in)])

function Base.show(io::IO, code::Code)
    println(io, "Code: input: $(code.in), output: $(code.out)")
    println(io, "body: ")
    for l in code.lines
        println("\t$(l[1]) = $(l[2])")
    end
    io
end

out(code::Code) = code.out
out(c) = c

function add(code1::Code, code2::Code)
    in = code1.in
    out = code2.out
    lines = cat(code1.lines, code2.lines; dims=1)
    Code(in, out, lines)
end

function add(code::Code, newline)
    in = code.in
    out = newline
    lines = [code.lines..., (out, code.out)]
    Code(in, out, lines)
end

function add(left, right)
    in = left
    out = right
    Code(in, out)
end

iscolon(ex::Expr) = ex.head == :call && ex.args[1] == :(:)
iscolon(x) = false
istuple(ex::Expr) = ex.head == :tuple
istuple(ex) = false

isleaf(ex::Expr) = istuple(ex) || iscolon(ex)
isleaf(ex::Symbol) = true
isleaf(ex::Int) = true

leftmost(ex) = isleaf(ex) ? ex : leftmost(ex.args[2])

leftmostnode(ex) = isleaf(ex.args[2]) ? ex : leftmostnode(ex.args[2])

function duplicate(code::Code, n::Int)
    in = code.in
    out = code.out
    n <= 0 && error("n should > 0")

    tp_cls = [(first(code.lines)[1], out), code.lines[2:end]...]
    tp_cls = [tp_cls for i = 1:(n-1)]

    Code(in, out, cat(code.lines, tp_cls...; dims=1))
end

function duplicate(c, n::Int)
    in = c
    out = c
    lines = [(c, c) for i = 1:n]
    Code(in, out, lines)
end

isdup(ex) = false
isdup(ex::Int) = true
isdup(ex::Expr) = iscolon(ex) && any(x-> x isa Int, ex.args) &&
    any(x -> istuple(x) || x isa Symbol, ex.args)

getint(ex::Int) = ex
getint(ex::Expr) = getleft(ex) isa Int ? getleft(ex) : getright(ex)
getsym(ex::Expr) = getleft(ex) isa Int ? getright(ex) : getleft(ex)

getleft(ex::Expr) = ex.args[2]
getright(ex::Expr) = ex.args[3]

function _to_code(node)
    if !isleaf(getleft(node))
        pre_code = _to_code(getleft(node))
    else
        pre_code = getleft(node)
    end

    rL = leftmost(getright(node)) #Symbol(or tuple) or Int
    if isdup(rL)#rL isa Int
        rightnode = getright(node)
        if isdup(rightnode)#getright(node) isa Int
            node.args[3] = rL isa Int ? out(pre_code) : getsym(rL)
        else
            leftmostnode(rightnode).args[2] = rL isa Int ? out(pre_code) : getsym(rL)
        end
        if getint(rL) == 1 #don't need duplicate
            code = pre_code
        else
            code = duplicate(pre_code, getint(rL))
        end
    else
        code = add(pre_code, rL)
    end

    if !isleaf(getright(node))
        next_code = _to_code(getright(node))
        code = add(code, next_code)
    end

    code
end

getin(ex::Expr) = iscolon(ex) ? getleft(ex) : ex
getin(x::Symbol) = x
getout(ex::Expr) = iscolon(ex) ? getright(ex) : ex
getout(x::Symbol) = x

function _postcode(code::Code)
    in = getin(code.in)
    out = getin(code.out)
    lines = map(code.lines) do (vab, farg)
        (getin(vab), getout(farg))
    end
    Code(in, out, lines)
end

function to_code(ex::Expr)
    code = _to_code(ex)
    _postcode(code)
end


text:
using SparseArrays, LinearAlgebra, KrylovKit #Arpack
include("deps/ed.jl")

# -----------------------------------------------------
#              Real-space ED (free)
# -----------------------------------------------------
"""
    perform_ed(beta = 8., nmax::Integer = 22) -> obs

Perform ED for L=2 system in real-space.

Returns a couple of observable `obs` as a `Dict{String, Any}`.

`nmax` specifies how many low energy state we keep in the ED.
"""
function perform_ed(; beta::Float64 = 8., nmax::Integer = 22)
    println("Performing real-space ED..."); flush(stdout);
    ns = 8
    params = Dict("txh"=>1.0, "txv"=>0.5, "tyh"=>-0.5, "tyv"=>-1.0, "mu"=>-0.5,
        "r"=>1.0, "lambdax"=>0.0, "lambday"=>0.0, "lambdaz"=>0.0, "lambda0"=>0.0)

    cup, cdn = fermiops(ns);
    cs = vcat(cup, cdn);

    @assert check_fermiops(cup)
    @assert check_fermiops(cdn)

    H = generate_H_SDW(params, ns, cup, cdn);

    # ARPACK:
    # evals, evecs, nconverged = eigs(H, nev=nmax+1, which=:SR); # eigenstates are columns of evecs

    # KrylovKit:
    r = eigsolve(H, nmax+1, :SR)
    evecs = hcat(r[2]...)
    evals = r[1]

    idcs = findall(x -> abs(x) < 1e-15, evecs)
    evecs[idcs] .= 0.;

    obs = Dict{String, Any}()

    # calc ETGF
    g = GF(cs, evecs, evals, beta);
    obs["greens"] = real(g)

    # calc occupation
    n = spzeros(size(cs[1])...)
    for c in cs
        n = n + c' * c
    end
    obs["occupation"] = EV(n, evecs, evals, beta)

    # calc ETPC, ETCDC
    obs["etpc_minus"], obs["etpc_plus"] = etpc(cs, evecs, evals, beta)
    obs["etcdc_minus"], obs["etcdc_plus"] = etcdc(cs, evecs, evals, beta)

    #@assert isapprox(real(sum(1 .- diag(g))), nbar)

    println("Done with ED."); flush(stdout);

    return obs
end

@testset "Compare to real-space ED (free)" begin
    mc_ed = mc_from_inxml("parameters/free_L_2_B_8.in.xml")
    ed = perform_ed(beta = 8.0)

    L, N = 2, 4
    flv = 2

    # ETGF must be real
    @test isreal(mc_ed.s.greens)
    g = real(mc_ed.s.greens)

    # Compare ETGF (not exactly because real-space ED isn't exact)
    @test maximum(g - ed["greens"][1:8, 1:8]) < 1e-5

    @test abs(occupation(mc_ed)*N*flv - ed["occupation"]/2) < 1e-5

    # Compare ETPC
    allocate_etpc!(mc_ed)
    measure_etpc!(mc_ed, g)
    @test isapprox(mc_ed.s.meas.etpc_minus, ed["etpc_minus"], atol=1e-5)
    @test isapprox(mc_ed.s.meas.etpc_plus, ed["etpc_plus"], atol=1e-5)

    # Compare ETCDC
    allocate_etcdc!(mc_ed)
    measure_etcdc!(mc_ed, g)
    @test isapprox(mc_ed.s.meas.etcdc_minus, ed["etcdc_minus"], atol=1e-5)
    @test isapprox(mc_ed.s.meas.etcdc_plus, ed["etcdc_plus"], atol=1e-4)
end

# -----------------------------------------------------
#              Momentum-space ED (free)
# -----------------------------------------------------
function concat_greens_k_space(gkx, gky)
    Gkx = Diagonal(gkx[:])
    Gky = Diagonal(gky[:])
    # gk = cat(Gkx, Gky, Gkx, Gky, dims=(1,2)); # Full 4Nx4N momentum-space Green's function
    gk = cat(Gkx, Gky, dims=(1,2)); # 2Nx2N momentum-space Green's function
end

"""
    perform_ed_k_space(L = 4, beta = 8.) -> obs

Perform ED for L=2 system in momentum-space single mode basis.

Returns a couple of observable `obs` as a `Dict{String, Any}`.
"""
function perform_ed_k_space(; L::Integer = 8, beta::Float64 = 10.)
    N = L^2
    params_x = Dict("th" => 1.0, "tv" => 0.5, "mu" => -0.5)
    params_y = Dict("th" => -0.5, "tv" => -1.0, "mu" => -0.5)
    obs = Dict{String, Any}()

    # ETGF
    gkx = ifftshift(GF_k_space(params_x, L, beta));
    gky = ifftshift(GF_k_space(params_y, L, beta));

    gk = concat_greens_k_space(gkx, gky); # Full momentum-space Green's function

    g = ifft_greens(gk, flv=2)
    @assert maximum(imag(g)) < 1e-12
    obs["greens"] = real(g) # Full (4xN, 4xN) real-space Green's function

    # TDGF
    taus = collect(0.0:0.1:beta)[1:end-1]

    # Gt0
    gt0_k_x = [TDGF_k_space(params_x, L, beta, tau) for tau in taus]
    gt0_k_y = [TDGF_k_space(params_y, L, beta, tau) for tau in taus]
    gt0_k = concat_greens_k_space.(gt0_k_x, gt0_k_y)
    obs["Gt0"] = real.(ifft_greens.(gt0_k; ifftshift=true, flv=2));

    # G0t
    g0t_k_x = [TDGF_k_space(params_x, L, beta, tau; G0t=true) for tau in taus]
    g0t_k_y = [TDGF_k_space(params_y, L, beta, tau; G0t=true) for tau in taus]
    g0t_k = concat_greens_k_space.(g0t_k_x, g0t_k_y)
    obs["G0t"] = real.(ifft_greens.(g0t_k; ifftshift=true, flv=2));

    return obs
end

@testset "Compare to momentum-space ED (free)" begin
    mc_edk = mc_from_inxml("parameters/free_L_8_B_10.in.xml")
    edk = perform_ed_k_space(L = 8, beta = 10.)

    L = mc_edk.p.L
    N = mc_edk.l.sites
    M = mc_edk.p.slices

    # Calculate TDGFs: Gt0, G0t
    allocate_tdgfs!(mc_edk)
    measure_tdgfs!(mc_edk)
    Gt0 = mc_edk.s.meas.Gt0
    G0t = mc_edk.s.meas.G0t

    @testset "ETGF" begin
        # ETGF
        @test isapprox(mc_edk.s.greens, edk["greens"])
    end

    @testset "TDGF" begin
        # --------------- DQMC self-consistency checks -----------
        # TDGF == ETGF consistency
        @test isapprox(Gt0[1], mc_edk.s.greens)

        # TDGF: tau = 0 difference btw Gt0 and G0t
        @test isapprox(Gt0[1] - G0t[1], I)

        # TDGF: G(tau,0) == -G(0,beta-tau) (with tau=0 excluded)
        should_be_zeros = Vector{Float64}(undef, length(2:M))
        t = () -> begin
            for (i, tau) in enumerate(2:M)
                should_be_zero = maximum(Gt0[tau] + G0t[end-(tau-2)])
                should_be_zeros[i] = should_be_zero
                should_be_zero < 1e-10 || (return false)
            end
            return true
        end
        @test t()

        # TDGF: Check realness of FFTs
        @test maximum(maximum.(imag.(fft_greens.(Gt0, flv=2)))) < 1e-12
        @test maximum(maximum.(imag.(fft_greens.(G0t, flv=2)))) < 1e-12

        # --------------- DQMC vs ED -----------
        # Gt0
        f = () -> begin
            for i in 1:M
                isapprox(Gt0[i], edk["Gt0"][i]) || return false
            end
            return true
        end
        @test f()

        # G0t
        f = () -> begin
            for i in 1:M
                isapprox(G0t[i], edk["G0t"][i]) || return false
            end
            return true
        end
        @test f()
    end

    @testset "ETPC" begin
        allocate_etpc!(mc_edk)
        measure_etpc!(mc_edk, mc_edk.s.greens)
        Pm = mc_edk.s.meas.etpc_minus
        Pp = mc_edk.s.meas.etpc_plus

        @test isreal(Pm)
        @test isreal(Pp)
        @test is_reflection_symmetric(Pm)
        @test is_reflection_symmetric(Pp)
        @test Pm ≈ Pp
    end

    @testset "ETCDC" begin
        allocate_etcdc!(mc_edk)
        measure_etcdc!(mc_edk, mc_edk.s.greens)
        Cm = mc_edk.s.meas.etcdc_minus
        Cp = mc_edk.s.meas.etcdc_plus

        @test isreal(Cm)
        @test isreal(Cp)
        @test is_reflection_symmetric(Cm)
        @test is_reflection_symmetric(Cp)
    end

    @testset "ZFPC" begin
        allocate_zfpc!(mc_edk)
        measure_zfpc!(mc_edk, Gt0)
        Pm = mc_edk.s.meas.zfpc_minus
        Pp = mc_edk.s.meas.zfpc_plus

        @test isreal(Pm)
        @test isreal(Pp)
        @test is_reflection_symmetric(Pm)
        @test is_reflection_symmetric(Pp)
        @test Pm ≈ Pp
    end

    @testset "ZFCDC" begin
        allocate_zfcdc!(mc_edk)
        measure_zfcdc!(mc_edk, mc_edk.s.greens, Gt0, G0t)
        Cm = mc_edk.s.meas.zfcdc_minus
        Cp = mc_edk.s.meas.zfcdc_plus

        @test isreal(Cm)
        @test isreal(Cp)
        @test is_reflection_symmetric(Cm)
        @test is_reflection_symmetric(Cp, tol=1e-10) # TODO: tol necessary for travis. why?
    end

    @testset "ZFCCC + Superfluid density" begin
        allocate_zfccc!(mc_edk)
        measure_zfccc!(mc_edk, mc_edk.s.greens, Gt0, G0t)
        Λxx = mc_edk.s.meas.zfccc_xx
        Λyy = mc_edk.s.meas.zfccc_yy

        @test isreal(Λxx)
        @test is_reflection_symmetric(Λxx)
        @test isreal(Λyy)
        @test is_reflection_symmetric(Λyy)

        @test isapprox(measure_sfdensity(mc_edk, Λxx, Λyy), -0.5356976525573391) # This value is only self-consistent.
    end

end

nothing

ext: jl

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
using WRS
using Distributions
using Test

n = 10000
d = Normal()

# generate data
sample_a = rand(d, n)
sample_b = rand(d, n)

# run comparison
results = pb2gen(sample_a, sample_b)

@test length(results) == 5


text:
module DProfile

export @dprofile

const samples = Dict{Int, Vector{UInt64}}()
const dict = Profile.getdict(UInt64[])

"""
`@profile <expression>`

runs expression while taking periodic backtraces.
The expression may spawn sub tasks on workers processes.
These are appended to an internal buffer of backtraces.

`DProfile.print()` will print the collated profile
"""
macro dprofile(ex)
    quote
        try
            @everywhere begin
                status = Profile.start_timer()
                if status < 0
                    error(Base.Profile.error_codes[status])
                end
            end
            $(esc(ex))
        finally
            @everywhere Profile.stop_timer()

            ps = procs()

            results = fetch.(map(p->remotecall(p) do
                samples = Base.Profile.fetch()
                (samples, Base.Profile.getdict(samples))
            end, ps))

            for (p, res) in zip(ps, results)
                samples, lineinfo = res
                if !haskey(DProfile.samples, p)
                    DProfile.samples[p] = UInt64[]
                end

                append!(DProfile.samples[p], samples)
                dict = reduce(merge, last.(results))

                merge!(DProfile.dict, dict)
            end
        end
    end
end

"""
`DProfile.print(;pids=procs())`

Prints the results of distributed profiling. `pids` will let you
constrain samples to only selected `pids`. Any arguments to
`Base.Profile.print` are also supported.
"""
function print(args...;pids = keys(samples), kwargs...)
    samps = UInt64[]
    for p in pids
        append!(samps, samples[p])
    end
    Base.Profile.print(args..., samps, dict; kwargs...)
end

"""
Clear backtraces from profile buffers on all processes.
"""
function clear()
    @everywhere Profile.clear()
    empty!(samples)
    empty!(dict)
    nothing
end

"""
`DProfile.init(; n::Integer, delay::Float64)`

Runs `Profile.init` with the given args on all processes.
"""
function init(args...;kwargs...)
    if length(args) == 0 && length(kwargs) == 0
        fetch.(
            remotecall.((pid) -> pid => Profile.init(), procs(), procs())
        )
    else
        @everywhere Profile.init($(args...); $(kwargs...))
    end
end

function save(file)
    open(file, "w") do io
        serialize(io, DProfile.samples)
        serialize(io, DProfile.dict)
    end
end

function load(file)
    open(file, "r") do io
        if !isempty(samples)
            warn("Overwriting previous samples")
        end
        empty!(DProfile.samples)
        merge!(DProfile.samples, deserialize(io))
        empty!(DProfile.dict)
        merge!(DProfile.dict, deserialize(io))
        nothing
    end
end

end # module


text:
function Base.rand(::Type{RandomSelection}, presences::T) where {T <: SimpleSDMLayer}
    @assert SimpleSDMLayers._inner_type(presences) <: Bool
    n_occ = sum(presences)
    iszero(n_occ) && throw(ArgumentError("The presences layer is empty"))
    _msk = mask(presences, presences)
    acceptable_cells = setdiff(keys(presences), keys(_msk))
    if length(acceptable_cells) < n_occ
        @warn "There are fewer acceptable pseudo-absences than occurrences"
    end
    pa_places = sample(acceptable_cells, min(length(acceptable_cells), n_occ), replace=false)
    pa = similar(presences, Bool)
    pa[pa_places] = fill(true, length(pa_places))
    return pa
end


text:
using Plots

# Frequency domain analysis of a simple system with time delays

s = tf("s")

P = delay(0.2) * ss(1/(s+1))

K = 3; Ti = 0.3;
C = ss(K*(1 + 1/s))

ω = exp10.(LinRange(-2,2,500))

nyquistplot(P*C, ω, xlims=(-2,1), ylims=(-2,2), gaincircles=false)

G_yd = feedback(P, C)
bodeplot(G_yd, ω, plotphase=false, title="Transfer function from load disturbances to output.")


text:
using YaoBlocks, YaoBlocks.AD
using YaoArrayRegister
using YaoBase: NotImplementedError
using Random, Test

AA(i, j) = control(i, j => shift(2π / (1 << (i - j + 1))))
B(n, i) = chain(n, i == j ? put(i => H) : AA(j, i) for j = i:n)
qftcirc(n) = chain(B(n, i) for i = 1:n)

function state_numgrad(f, reg)
    map(1:length(reg.state)) do i
        reg.state[i] += 1e-5
        pos = f(reg)
        reg.state[i] -= 2e-5
        neg = f(reg)
        reg.state[i] += 1e-5 + 1e-5im
        ipos = f(reg)
        reg.state[i] -= 2e-5im
        ineg = f(reg)
        reg.state[i] += 1e-5im
        ((pos - neg) + im * (ipos - ineg)) / 2e-5 / 2
    end
end

@testset "expect grad" begin
    Random.seed!(2)
    nbit = 4
    c = qftcirc(nbit)
    H = repeat(nbit, X, 1:nbit)
    for reg in [rand_state(nbit), rand_state(nbit; nbatch = 10)]
        adjin, adjparams = expect'(H, reg => c)
        numgrad =
            YaoBlocks.AD.ng(x -> sum(expect(H, reg => dispatch!(c, x))), parameters(c))
        @test adjparams ≈ vec(numgrad)
        numgrad2 = state_numgrad(reg -> sum(expect(H, reg => c)), reg)
        @test isapprox(vec(adjin.state), numgrad2, atol = 1e-4)

        adjin2 = expect'(H, reg)
        numgrad3 = state_numgrad(reg -> sum(expect(H, reg)), reg)
        @test isapprox(vec(adjin2.state), numgrad3, atol = 1e-4)
    end
end

@testset "fidelity grad" begin
    nbit = 4
    Random.seed!(2)
    for nbatch in [1, 10]
        reg1 = rand_state(nbit; nbatch = nbatch)
        reg2 = rand_state(nbit; nbatch = nbatch)
        c1 = qftcirc(nbit)
        c2 = chain(put(nbit, 2 => Rx(0.5)), control(nbit, 1, 3 => Ry(0.5)))

        g1, g2 = fidelity'(reg1, reg2)
        @test isapprox(
            vec(g1.state),
            state_numgrad(reg1 -> sum(fidelity(reg1, reg2)), reg1),
            atol = 1e-4,
        )
        @test isapprox(
            vec(g2.state),
            state_numgrad(reg2 -> sum(fidelity(reg1, reg2)), reg2),
            atol = 1e-4,
        )

        (g1, pg1), (g2, pg2) = fidelity'(reg1 => c1, reg2 => c2)
        npg1 = YaoBlocks.AD.ng(
            x -> sum(fidelity(reg1 => dispatch!(c1, x), reg2 => c2)),
            parameters(c1),
        )
        npg2 = YaoBlocks.AD.ng(
            x -> sum(fidelity(reg1 => c1, reg2 => dispatch!(c2, x))),
            parameters(c2),
        )
        @test isapprox(pg1, vec(npg1), atol = 1e-5)
        @test isapprox(pg2, vec(npg2), atol = 1e-5)
        @test isapprox(
            vec(g1.state),
            state_numgrad(reg1 -> sum(fidelity(reg1 => c1, reg2 => c2)), reg1),
            atol = 1e-4,
        )
        @test isapprox(
            vec(g2.state),
            state_numgrad(reg2 -> sum(fidelity(reg1 => c1, reg2 => c2)), reg2),
            atol = 1e-4,
        )
    end

    nbatch = 1
    reg1 = rand_state(nbit; nbatch = nbatch) |> focus!(2, 1, 4)
    reg2 = rand_state(nbit; nbatch = nbatch) |> focus!(2, 1, 4)
    c1 = qftcirc(3)
    c2 = chain(put(3, 2 => Rx(0.5)), control(3, 1, 3 => Ry(0.5)))

    @test_throws ArgumentError fidelity'(reg1, reg2)
    #g1, g2 = fidelity'(reg1, reg2)
    #@test isapprox(vec(g1.state), state_numgrad(reg1->sum(fidelity(reg1, reg2)), reg1), atol=1e-5)
    #@test isapprox(vec(g2.state), state_numgrad(reg2->sum(fidelity(reg1, reg2)), reg2), atol=1e-5)

    #(g1,pg1), (g2,pg2) = fidelity'(reg1 => c1, reg2 => c2)
    #npg1 = YaoBlocks.AD.ng(x -> sum(fidelity(reg1 => dispatch!(c1, x), reg2=>c2)), parameters(c1))
    #npg2 = YaoBlocks.AD.ng(x -> sum(fidelity(reg1 => c1, reg2 => dispatch!(c2, x))), parameters(c2))
    #@test isapprox(pg1, vec(npg1),atol=1e-5)
    #@test isapprox(pg2, vec(npg2),atol=1e-5)
    #@test isapprox(vec(g1.state), state_numgrad(reg1->sum(fidelity(reg1=>c1, reg2=>c2)), reg1), atol=1e-5)
    #@test isapprox(vec(g2.state), state_numgrad(reg2->sum(fidelity(reg1=>c1, reg2=>c2)), reg2), atol=1e-5)
end

@testset "operator fideliy" begin
    nbit = 4
    c1 = qftcirc(nbit)
    c2 = chain(put(nbit, 2 => Rx(0.5)), control(nbit, 1, 3 => Ry(0.5)))

    g1, g2 = operator_fidelity'(c1, c2)
    npg1 = YaoBlocks.AD.ng(x -> operator_fidelity(dispatch!(c1, x), c2), parameters(c1))
    npg2 = YaoBlocks.AD.ng(x -> operator_fidelity(c1, dispatch!(c2, x)), parameters(c2))
    @test isapprox(g1, vec(npg1), atol = 1e-5)
    @test isapprox(g2, vec(npg2), atol = 1e-5)
end


text:
sir_model = @reaction_network rn begin
    0.1/1000, s + i --> 2i
    p1, i --> r
end p1
function sir_f(du,u,p,t)
    du[1] = -(0.1/1000)*u[1]*u[2]
    du[2] = (0.1/1000)*u[1]*u[2] - 0.01 * u[2]
    du[3] = p[1] * u[2]
end

for i = 1:100
    p1 = 0.01; p = [p1];
    u0 = [990. + 200 * rand(1)[1],10. + 2* rand(1)[1],0.]
    tspan = (0.,10. + 5 * rand(1)[1])

    prob1 = ODEProblem(sir_model,u0,tspan,p)
    sol1 = solve(prob1,Tsit5())
    prob2 = ODEProblem(sir_f,u0,tspan,p)
    sol2 = solve(prob2,Tsit5())
    @test sol1[end][1] == sol2[end][1]
end

model = @reaction_network rn begin
    (d,1000), X ↔ 0
end d

## For Julialang/julia#28356
i = 1.0
prob = SDEProblem(model,[1000.0+i],(0.,200.),[i])
sol = solve(prob, EM(), dt = 0.01)
@test sol[end][1] < 2000
##

for i in [1., 2., 3., 4., 5.]
    prob = SDEProblem(model,[1000.0+i],(0.,200.),[i])
    sol = solve(prob, EM(), dt = 0.01)
    @test sol[end][1] < 2000
end

#Tests time dependence
equi_model = @reaction_network rn begin
    100*(tanh(t-100)+2), 0 --> X
    1, X --> 0
end
prob1 = ODEProblem(equi_model,[100.],(0.,100.))
sol1 = solve(prob1,Tsit5())
prob2 = ODEProblem(equi_model,[100.],(0.,200.))
sol2 = solve(prob2,Tsit5())
@test 1.5*sol1[end][1] < sol2[end][1]


text:
using InteractiveUtils: code_llvm
using Test

include("../benchmark/gemm.jl")

"""
    llvm_ir(f, args) :: String

Get LLVM IR of `f(args...)` as a string.
"""
llvm_ir(f, args) = sprint(code_llvm, f, Base.typesof(args...))

nmatches(r, s) = count(_ -> true, eachmatch(r, s))

@testset "manmul!" begin
    M = ones(0, 0)
    @test nmatches(r"fmul contract <4 x double>", llvm_ir(manmul!, (M, M, M))) >= 4
end

@testset "ndrmul!" begin
    M = ones(0, 0)
    if VERSION < v"1.2-"
        CAB = _ndrmul_reducible(M, M, M)
        rf = _ndrmul_rf()
        ir = llvm_ir(transduce, (rf, nothing, CAB))
    else
        ir = llvm_ir(ndrmul!, (M, M, M))
    end
    @test nmatches(r"fmul contract <4 x double>", ir) >= 4
end


text:
using eFEM, Test

@testset "Laplace order 1" begin
  # define mesh through eFEMpart code
  mesh = squareMesh([-1.2,1.1,-1.2,1.1],12,1)

  # define variable name and operator type
  OperatorType = :Poisson2D

  # Define auxiliary information
  dNodes = Dirichlet(:left,:right)
  nNodes = Neumann(:top,:bottom)

  # Define dirichlet, neumann, and forcing functions
  dBCf = Dirichlet((x,y) -> (x==-1.2 ? -1.0 : y))
  nBCf = Neumann((x,y) -> 0.0)
  ff   = Forcing((x,y) -> 0.0)

  # required technicality for grouping auxiliary information
  Nodes = [dNodes,nNodes]; bcfun = [dBCf,nBCf,ff]

  # define problem variable
  prob = Problem(mesh,Nodes,bcfun,OperatorType)

  # solve problem
  sol = solve(prob,mesh)

  FEMstats(mesh)
  FEMstats(mesh,prob)

  # save solution
  vtkname = Path("test-data","laplace_soln")
  sd = ScalarData(sol.u,sol.u)
  sn = ScalarNames("testing_var","testing_var_2")
  vtksave(mesh,sd,sn,vtkname)

  rm("test-data/laplace_soln.vtk")
end

@testset "Advection-Diffusion order 2" begin
  mesh = squareMesh([-1.2,1.1,-1.2,1.1],12,2)

  OperatorType = :AdvDiff2D

  dNodes = Dirichlet(:left,:top)
  nNodes = Neumann(:right,:bottom)

  dBCf   = Dirichlet((x,y) -> (x==-1.2 ? -0.3 : 0.0))
  ff     = Forcing((x,y) -> 0.0)

  Pe = 7.3
  windX(x,y) =  -y;  windY(x,y) = x
  wx = [windX.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]
  wy = [windY.(mesh.xy[i].x,mesh.xy[i].y) for i=1:length(mesh.xy)]
  param  = AdvDiffParam(wx,wy,Pe)

  Nodes = [dNodes]
  bcfun = [dBCf,ff]

  prob = Problem(mesh,Nodes,bcfun,OperatorType)
  sol = solve(prob,mesh,param)

  vtkname = Path("test-data","advdiff_soln")
  sd = ScalarData(sol.u)
  sn = ScalarNames("testing_var")
  vd = VectorData([wx,wy],[wx,wy])
  vn = VectorNames("velocity","velocity_2")
  vtksave(mesh,sd,sn,vd,vn,vtkname)

  rm("test-data/advdiff_soln.vtk")
end

@testset "Stokes Example" begin
  mesh = squareMeshFluid([-1.2,1.1,-1.2,1.1],12)

  param = FluidParam(2.3)
  OperatorType = :Stokes2D

  dNodes = Dirichlet(:left,:bottom,:top)
  nNodes = Neumann(:right)
  Nodes = [dNodes,nNodes]

  dUBCf = Dirichlet((x,y) -> (x==-1.2 ? -(1.1-y)*(-1.2-y) : 0.0))
  dVBCf = Dirichlet((x,y) -> 0.0)
  bcfun = [dUBCf,dVBCf]

  prob = Problem(mesh,Nodes,bcfun,OperatorType)
  sol = solve(prob,mesh,param)

  # save solution
  vtkname = Path("test-data","stokes_soln")
  sd = ScalarData(sol.p)
  sn = ScalarNames("pressure")
  vd = VectorData([sol.u,sol.v])
  vn = VectorNames("velocity")
  vtksave(mesh,sd,sn,vd,vn,vtkname)

  rm("test-data/stokes_soln.vtk")
end

@testset "Stokes Axisymmetric Example" begin
  mesh = squareMeshFluid([-1.2,1.1,0.0,1.1],12)

  param = FluidParam(2.3)
  OperatorType = :StokesAS

  dUNodes = Dirichlet(:left,:top,:right)
  dVNodes = Dirichlet(:left,:top,:bottom,:right)
  Nodes = [dUNodes,dVNodes]

  dUBCf = Dirichlet((x,y) -> (x==-1.2 ? -(1.1-y)*(-1.2-y) : 0.0))
  dVBCf = Dirichlet((x,y) -> 0.0)
  bcfun = [dUBCf,dVBCf]

  prob = Problem(mesh,Nodes,bcfun,OperatorType)
  sol = solve(prob,mesh,param)

  # save solution
  vtkname = Path("test-data","stokes_as_soln")
  sd = ScalarData(sol.p)
  sn = ScalarNames("pressure")
  vd = VectorData([sol.u,sol.v])
  vn = VectorNames("velocity")
  vtksave(mesh,sd,sn,vd,vn,vtkname)

  rm("test-data/stokes_as_soln.vtk")
end

@testset "Multiphase Brinkman" begin
  x0 = 2.0
  mesh = squareMeshFluid([-x0,x0,-1,1],32)

  μ0 = 1.0
  a10 = 0.0
  a20 = 1.0

  μ  = μ0*ones(length(mesh.xy))
  a1  = a10*ones(length(mesh.xy))
  a2  = a20*ones(length(mesh.xy))
  G  = 1e0
  param = BrinkmanMPParam(μ,a1,a2)

  OperatorType = :BrinkmanMP2D

  dNodes = Dirichlet(:left,:top,:bottom)
  nNodes = Neumann(:right)
  Nodes = [dNodes,nNodes]

  dUBCf = Dirichlet((x,y) -> (x==-x0 ? 1.0-y^2 : 0.0))
  dVBCf = Dirichlet((x,y) -> 0.0)
  bcfun = [dUBCf,dVBCf]

  prob = Problem(mesh,Nodes,bcfun,OperatorType)
  sol = solve(prob,mesh,param)

  vtkname = Path("test-data","mpb_soln")
  vd = VectorData([sol.u,sol.v])
  vn = VectorNames("velocity")
  vtksave(mesh,vd,vn,vtkname)

  rm("test-data/mpb_soln.vtk")
end


text:
struct TDVPOrder{order,direction} end

TDVPOrder(order::Int, direction::Base.Ordering) = TDVPOrder{order,direction}()

orderings(::TDVPOrder) = error("Not implemented")
sub_time_steps(::TDVPOrder) = error("Not implemented")

function orderings(::TDVPOrder{1,direction}) where {direction}
  return [direction, Base.ReverseOrdering(direction)]
end
sub_time_steps(::TDVPOrder{1}) = [1.0, 0.0]

function orderings(::TDVPOrder{2,direction}) where {direction}
  return [direction, Base.ReverseOrdering(direction)]
end
sub_time_steps(::TDVPOrder{2}) = [1.0 / 2.0, 1.0 / 2.0]

function orderings(::TDVPOrder{4,direction}) where {direction}
  return [direction, Base.ReverseOrdering(direction)]
end
function sub_time_steps(::TDVPOrder{4})
  s = 1.0 / (2.0 - 2.0^(1.0 / 3.0))
  return [s / 2.0, s / 2.0, (1.0 - 2.0 * s) / 2.0, (1.0 - 2.0 * s) / 2.0, s / 2.0, s / 2.0]
end


text:
function fibonacci(n::Int)
    n==0 && return 0
    n==1 && return 1
    fibonacci(n-1) + fibonacci(n-2)
    end

print(fibonacci(40))


text:
module TestConsistency
    import Base.Test.@test
    import Benchmark.compare

    function f1()
        x = 0.0
        for i in 1:1_000_000
            x = rand()
        end
        return
    end

    comparisons = compare([f1, f1, f1], 100)
    @test maximum(comparisons[:Relative]) < 1.1
end


text:
mutable struct SubpassDescription
    flags::vk.VkSubpassDescriptionFlags
    pipelineBindPoint::vk.VkPipelineBindPoint
    resolveAttachments::Vector{vk.VkAttachmentReference}
    depthStencilAttachment::Ref{vk.VkAttachmentReference}
    preserveAttachments::Vector{UInt32}
    colorAttachments::Vector{vk.VkAttachmentReference}
    inputAttachments::Vector{vk.VkAttachmentReference}

    mHandleRef::Ref{vk.VkSubpassDescription}

    function SubpassDescription()
        this = new()
        this.flags = 0
        this.pipelineBindPoint = vk.VK_PIPELINE_BIND_POINT_GRAPHICS
        this.resolveAttachments = Vector{vk.VkAttachmentReference}()
        this.preserveAttachments = Vector{UInt32}()
        this.colorAttachments = Vector{vk.VkAttachmentReference}()
        this.inputAttachments = Vector{vk.VkAttachmentReference}()
        return this
    end
end

function depth(this::SubpassDescription, index::UInt32)
    this.depthStencilAttachment = Ref(vk.vkAttachmentReference(
                                            index, #attachment::UInt32
                                            vk.VK_IMAGE_LAYOUT_DEPTH_STENCIL_ATTACHMENT_OPTIMAL #layout::VkImageLayout
                                        ))
end

function colors(this::SubpassDescription, indices::Vector{UInt32})
    resize!(this.colorAttachments, length(indices))
    for idx in indices
        push!(this.colorAttachments, vk.VkAttachmentReference(
                                        idx, #attachment::UInt32
                                        vk.VK_IMAGE_LAYOUT_COLOR_ATTACHMENT_OPTIMAL #layout::VkImageLayout
                                    ))
    end
end

function inputsColor(this::SubpassDescription, indices::Vector{UInt32})
    oldSize = length(this.inputAttachments)
    resize!(this.inputAttachments, oldSize + length(indices))
    for idx in indices
        push!(this.inputAttachments, vk.VkAttachmentReference(
                                        idx, #attachment::UInt32
                                        vk.VK_IMAGE_LAYOUT_SHADER_READ_ONLY_OPTIMAL #layout::VkImageLayout
                                    ))
    end
end

function inputsDepth(this::SubpassDescription, indices::Vector{UInt32})
    oldSize = length(this.inputAttachments)
    resize!(this.inputAttachments, oldSize + length(indices))
    for idx in indices
        push!(this.inputAttachments, vk.VkAttachmentReference(
                                        idx, #attachment::UInt32
                                        vk.VK_IMAGE_LAYOUT_DEPTH_ATTACHMENT_STENCIL_READ_ONLY_OPTIMAL #layout::VkImageLayout
                                    ))
    end
end

function perViewAttributes(this::SubpassDescription, active::Bool = true, allComponents::Bool = false)
    if active
        this.flags = vk.VK_SUBPASS_DESCRIPTION_PER_VIEW_ATTRIBUTES_BIT_NVX
        if !allComponents
            this.flags |= vk.VK_SUBPASS_DESCRIPTION_PER_VIEW_POSITION_X_ONLY_BIT_NVX
        end
    else
        this.flags = 0
    end
end

function commit(this::SubpassDescription)
    colorAttCount = length(this.colorAttachments)
    resolveAttCount = length(this.resolveAttachments)
    this.mHandleRef = Ref(vk.VkSubpassDescription(
        this.flags, #flags::VkSubpassDescriptionFlags
        this.pipelineBindPoint, #pipelineBindPoint::VkPipelineBindPoint
        length(this.inputAttachments), #inputAttachmentCount::UInt32
        pointer(this.inputAttachments), #pInputAttachments::Ptr{VkAttachmentReference}
        colorAttCount, #colorAttachmentCount::UInt32
        pointer(this.colorAttachments), #pColorAttachments::Ptr{VkAttachmentReference}
        colorAttCount <= resolveAttCount ? pointer(this.resolveAttachments) : C_NULL, #pResolveAttachments::Ptr{VkAttachmentReference}
        Base.unsafe_convert(Ptr{vk.VkAttachmentReference}, this.depthStencilAttachment), #pDepthStencilAttachment::Ptr{VkAttachmentReference}
        length(this.preserveAttachments), #preserveAttachmentCount::UInt32
        pointer(this.preserveAttachments) #pPreserveAttachments::Ptr{UInt32}
    ))
end

function handleRef(this::SubpassDescription)::Ref{vk.VkSubpassDescription}
    return this.mHandleRef
end


text:
testdims(L::Integer) = ntuple(i -> 3 + i, Val(L))

randlike(x::Number) = rand(typeof(x))
randlike(A::AbstractArray{<:Number}) = rand(eltype(A), size(A)...)
randlike(A::AbstractArray{<:Number,0}) = rand!(zeros(eltype(A)))

function randlike(A::AbstractArray{<:AbstractArray{V,M},N}) where {V,M,N}
    B = Array{Array{V,M},N}(undef, size(A)...)
    for I in eachindex(B)
        B[I] = randlike(A[I])
    end
    return B
end

function randlike!(A::AbstractArray{<:AbstractArray})
    for I in eachindex(A)
        copyto!(A[I], randlike(A[I]))
    end
    return A
end

function randlike!(A::AbstractArray{<:Number})
    for I in eachindex(A)
        A[I] = randlike(A[I])
    end
    return A
end

macro qe(ex)
    Expr(:quote, MacroTools.striplines(esc(ex)))
end

macro test_inferred(ex)
    @qe begin
        $Test.@test (($Test.@inferred $ex); true)
    end
end

macro test_noalloc(ex)
    @qe begin
        local nbytes = $memory($BenchmarkTools.@benchmark $ex samples = 1 evals = 1)
        $iszero(nbytes) ? true : $error("Allocated $(BenchmarkTools.prettymemory(nbytes))")
    end
end

macro test_getindex(A, I)
    @qe begin
        if $I isa $Tuple
            @test_inferred $A[$I...]
            if $index_dimsum($I...) == ()
                $Test.@test $typeof($A[$I...]) == $eltype($A)
            else
                $Test.@test $size($A[$I...]) == $map($length, $index_shape($to_indices($A, $I)...))
                $Test.@test $eltype($A[$I...]) == $eltype($A)
            end
        else
            @test_inferred $A[$I]
            if $index_dimsum($I) == () && !($I isa AbstractArray{<:Any,0})
                $Test.@test $typeof($A[$I]) == $eltype($A)
            else
                $Test.@test $size($A[$I]) == $map($length, $index_shape($to_indices($A, ($I,))...))
                $Test.@test $eltype($A[$I]) == $eltype($A)
            end
        end
    end
end

macro test_setindex!(A, I)
    @qe begin
        if $I isa $Tuple
            local x = $randlike($A[$I...])
            @test_inferred $setindex!($A, x, $I...)
            $Test.@test $setindex!($A, x, $I...) === $A
            $Test.@test $A[$I...] == x
        else
            local x = $randlike($A[$I])
            @test_inferred $setindex!($A, x, $I)
            $Test.@test $setindex!($A, x, $I) === $A
            $Test.@test $A[$I] == x
        end
    end
end

macro test_getindices(A, Is)
    @qe begin
        $Test.@test $all($Is) do I
            if I isa $Tuple
                if $index_dimsum(I...) == ()
                    $typeof($A[I...]) == $eltype($A)
                else
                    $eltype($A[I...]) == $eltype($A) &&
                    $size($A[I...]) == $map($length, $index_shape($to_indices($A, I)...))
                end
            else
                #if $index_dimsum(I) == ()
                if $index_dimsum(I) == () && !(I isa AbstractArray{<:Any,0})
                    $typeof($A[I]) == $eltype($A)
                else
                    $eltype($A[I]) == $eltype($A) &&
                    $size($A[I]) == $map($length, $index_shape($to_indices($A, (I,))...))
                end
            end
        end
    end
end

macro test_setindices!(A, Is)
    @qe begin
        $Test.@test $all($Is) do I
            if I isa $Tuple
                x = $randlike($A[I...])
                $A[I...] != x && $setindex!($A, x, I...) === $A && $A[I...] == x
            else
                x = $randlike($A[I])
                $A[I] != x && $setindex!($A, x, I) === $A && $A[I] == x
            end
        end
    end
end

macro test_all_equal(A, B, Is)
    @qe begin
        $Test.@test all($Is) do I
            if I isa Tuple
                $A[I...] == $B[I...] #&& typeof($A[I...]) == typeof($B[I...])
            else
                $A[I] == $B[I] #&& typeof($A[I]) == typeof($B[I])
            end
        end
    end
end

macro test_copyto!(dest, src)
    @qe begin
        $randlike!($dest)
        @test_inferred $copyto!($dest, $src)
        $Test.@test $copyto!($dest, $src) === $dest
        $Test.@test $all(I -> $dest[I] == $src[I], $eachindex($dest, $src))
        $Test.@test $dest == $src

        if $length($dest) > 2
            $randlike!($dest)
            @test_inferred $copyto!($dest, 2, $src, 2, 1)
            $Test.@test $copyto!($dest, 2, $src, 2, 1) === $dest
            $Test.@test $dest[2] == $src[2] && $dest[1] != $src[1] && $dest[3:end] != $src[3:end]
        else
            $Test.@test_skip @test_inferred $copyto!($dest, 2, $src, 2, 1)
            $Test.@test_skip $copyto!($dest, 2, $src, 2, 1) === $dest
            $Test.@test_skip $dest[2] == $src[2] &&
                             $dest[1] != $src[1] &&
                             $dest[3:end] != $src[3:end]
        end
    end
end

macro test_array_attributes(A)
    @qe begin
        $Test.@test $eltype($A) === $typeof($first($A))
        $Test.@test $ndims($A) == $length($axes($A))
        $Test.@test $size($A) == $map($length, $axes($A))
        $Test.@test $length($A) == $prod($size($A))

        $Test.@test $all(1:$ndims($A)) do dim
            $size($A)[dim] == $size($A, dim)
        end
        $Test.@test $size($A, $ndims($A) + 1) == 1

        $Test.@test $all(1:$ndims($A)) do dim
            $axes($A)[dim] == $axes($A, dim)
        end
        $Test.@test $axes($A, $ndims($A) + 1) == $Base.OneTo(1)

        @test_inferred $eltype($A)
        @test_inferred $ndims($A)
        @test_inferred $axes($A)
        @test_inferred $size($A)
        @test_inferred $length($A)
    end
end

macro test_similar(A, T, dims)
    @qe begin
        @test_inferred $similar($A)
        $Test.@test $eltype($similar($A)) === $eltype($A)
        $Test.@test $size($similar($A)) == $size($A)

        @test_inferred $similar($A, $T)
        $Test.@test $eltype($similar($A, $T)) === $T
        $Test.@test $size($similar($A, $T)) == $size($A)

        @test_inferred $similar($A, $dims)
        $Test.@test eltype($similar($A, $dims)) === $eltype($A)
        $Test.@test size($similar($A, $dims)) == $dims

        @test_inferred $similar($A, $T, $dims)
        $Test.@test eltype($similar($A, $T, $dims)) === $T
        $Test.@test size($similar($A, $T, $dims)) == $dims
    end
end

function test_indexing_AB(AB::Function)
    @testset "LinearIndices" begin
        let (A, _) = AB()
            @test_getindex A first(LinearIndices(A))
            @test_setindex! A first(LinearIndices(A))
        end
        let (A, B) = AB()
            @test LinearIndices(A) === LinearIndices(B)
            @test_all_equal A B LinearIndices(A)
            @test_getindices A LinearIndices(A)
            @test_setindices! A LinearIndices(A)
        end
        let (A, B) = AB()
            @test A[LinearIndices(A)] == B[LinearIndices(A)]
            @test_getindex A LinearIndices(A)
            @test_setindex! A LinearIndices(A)
        end
        let (A, B) = AB()
            @test begin
                I = hcat(LinearIndices(A), LinearIndices(A))
                A[I] == B[I]
            end
            @test_getindex A hcat(LinearIndices(A), LinearIndices(A))
            #@test_setindex! A hcat(LinearIndices(A), LinearIndices(A))
        end
    end
    @testset "CartesianIndices" begin
        let (A, _) = AB()
            @test_getindex A first(CartesianIndices(A))
            @test_setindex! A first(CartesianIndices(A))
        end
        let (A, B) = AB()
            @test CartesianIndices(A) === CartesianIndices(B)
            @test_all_equal A B CartesianIndices(A)
            @test_getindices A CartesianIndices(A)
            @test_setindices! A CartesianIndices(A)
        end
        let (A, B) = AB()
            @test A[CartesianIndices(A)] == B[CartesianIndices(A)]
            @test_getindex A CartesianIndices(A)
            @test_setindex! A CartesianIndices(A)
        end
        let (A, B) = AB()
            @test begin
                I = hcat(CartesianIndices(A), CartesianIndices(A))
                A[I] == B[I]
            end
            @test_getindex A hcat(CartesianIndices(A), CartesianIndices(A))
            #@test_setindex! A CartesianIndices(A)
        end
    end
    @testset "trailing singleton" begin
        let (A, _) = AB()
            @test_getindex A (Tuple(first(CartesianIndices(A)))..., 1)
            @test_setindex! A (Tuple(first(CartesianIndices(A)))..., 1)
        end
        let (A, B) = AB()
            @test_all_equal A B [(Tuple(I)..., 1) for I in CartesianIndices(A)]
            @test_getindices A [(Tuple(I)..., 1) for I in CartesianIndices(A)]
            @test_setindices! A [(Tuple(I)..., 1) for I in CartesianIndices(A)]
        end
    end
    @testset "trailing colon" begin
        let (A, _) = AB()
            @test_getindex A (Tuple(first(CartesianIndices(A)))..., :)
            @test_setindex! A (Tuple(first(CartesianIndices(A)))..., :)
        end
        let (A, B) = AB()
            @test_all_equal A B [(Tuple(I)..., :) for I in CartesianIndices(A)]
            @test_getindices A [(Tuple(I)..., :) for I in CartesianIndices(A)]
            @test_setindices! A [(Tuple(I)..., :) for I in CartesianIndices(A)]
        end
    end
    @testset "dropped singleton" begin
        _, B = AB()
        if ndims(B) > 0 && size(B, ndims(B)) == 1
            let (A, _) = AB()
                @test_getindex A front(Tuple(first(CartesianIndices(A))))
                @test_setindex! A front(Tuple(first(CartesianIndices(A))))
            end
            let (A, B) = AB()
                @test_all_equal A B [front(Tuple(I)) for I in CartesianIndices(A)]
                @test_getindices A [front(Tuple(I)) for I in CartesianIndices(A)]
                @test_setindices! A [front(Tuple(I)) for I in CartesianIndices(A)]
            end
        end
    end
    @testset "logical" begin
        let (A, B) = AB()
            @test begin
                I = map(isodd, LinearIndices(A))
                A[I] == B[I]
            end
            @test_getindex A map(isodd, LinearIndices(A))
            @test_setindex! A map(isodd, LinearIndices(A))
        end
        let (A, B) = AB()
            @test begin
                I = Tuple(map(isodd, ax) for ax in axes(A))
                A[I...] == B[I...]
            end
            @test_getindex A Tuple(map(isodd, ax) for ax in axes(A))
            @test_setindex! A Tuple(map(isodd, ax) for ax in axes(A))
        end
    end
    @testset "single colon" begin
        A, B = AB()
        @test A[:] == B[:]
    end
end

function test_array_AB(AB::Function)
    @testset "array attributes" begin
        A, _ = AB()
        @test_array_attributes A

        @test_noalloc eltype($A)
        @test_noalloc ndims($A)
        @test_noalloc axes($A)
        @test_noalloc size($A)
        @test_noalloc length($A)
    end

    @testset "indexing" begin
        test_indexing_AB(AB)
    end

    @testset "copy/copyto!/equality" begin
        let (A1, _) = AB()
            A2 = copy(A1)
            @test A1 == A2
            randlike!(A1)
            @test A1 != A2
        end
        let (A, B) = AB()
            @test_copyto! A B
        end
        let (A, B) = AB()
            @test_copyto! B A
        end
        let A1 = first(AB()), A2 = first(AB())
            @test_copyto! A1 A2
        end
    end
end


text:
using KUnet
using MAT
include("../src/KUparser.jl")
include("../src/flist.jl")
isdefined(:d1mat) || (f=matopen("d1.mat");d1mat=read(f,"d1");close(f))

function readcorpus(m)
    d = KUparser.Sentence[]
    for ss in m
        s = KUparser.Sentence()
        s.wvec = ss["wvec"]
        s.head = isa(ss["head"], Number) ? [ss["head"]] : vec(ss["head"])
        push!(d, s)
    end
    return d
end

dev=readcorpus(d1mat["dev"])
trn=readcorpus(d1mat["trn"])
tst=readcorpus(d1mat["tst"])
s=dev[1]
p=KUparser.ArcHybrid(size(s.wvec,2))
v=KUparser.valid(p)
feats=Flist.fv021a
KUnet.gpu(false)
net = [KUnet.Layer("m11.h5"), KUnet.Layer("m12.h5")]
x=KUparser.features(p,s,feats)
y=KUnet.predict(net,x)


text:
module Environment

include("space/Space.jl")
using .Space
export BoxSpace, DiscreteSpace, MultiBinarySpace, MultiDiscreteSpace, gymspace2jlspace, sample

include("abstractenvironment.jl")
include("gymenvironment.jl")
export GymEnv, step!, reset!, render, close!, seed!, get_state, get_reward, is_end, action_space, observation_space

end # module


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

## dummy stub for https://github.com/JuliaBinaryWrappers/MPFR_jll.jl
baremodule MPFR_jll
using Base, Libdl, GMP_jll
Base.Experimental.@compiler_options compile=min optimize=0 infer=false

const PATH_list = String[]
const LIBPATH_list = String[]

export libmpfr

# These get calculated in __init__()
const PATH = Ref("")
const LIBPATH = Ref("")
artifact_dir = ""
libmpfr_handle = C_NULL
libmpfr_path = ""

if Sys.iswindows()
    const libmpfr = "libmpfr-6.dll"
elseif Sys.isapple()
    const libmpfr = "@rpath/libmpfr.6.dylib"
else
    const libmpfr = "libmpfr.so.6"
end

function __init__()
    global libmpfr_handle = dlopen(libmpfr)
    global libmpfr_path = dlpath(libmpfr_handle)
    global artifact_dir = dirname(Sys.BINDIR)
    LIBPATH[] = dirname(libmpfr_path)
    push!(LIBPATH_list, LIBPATH[])
end

# JLLWrappers API compatibility shims.  Note that not all of these will really make sense.
# For instance, `find_artifact_dir()` won't actually be the artifact directory, because
# there isn't one.  It instead returns the overall Julia prefix.
is_available() = true
find_artifact_dir() = artifact_dir
dev_jll() = error("stdlib JLLs cannot be dev'ed")
best_wrapper = nothing
get_libmpfr_path() = libmpfr_path

end  # module MPFR_jll

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
"""
Average forcing data.
"""
function average_forcings(base_folder, target_folder, soil_param, wsh_single, resolution)

    iboxes = getfield(wsh_single, resolution)
    gridcel = getfield(wsh_single, :gridcel)

    for ibox in unique(iboxes)

        isel = findin(soil_param[:gridcel], gridcel[iboxes .== ibox])

        ncells = length(isel)

        df_tmp = soil_param[isel, :]

        prec_final, tmin_final, tmax_final, wind_final = 0, 0, 0, 0

        for row in eachrow(df_tmp)

            lat = @sprintf("%0.5f", row[:lat][1])
            lon = @sprintf("%0.5f", row[:lon][1])

            file_src = joinpath(base_folder, "forcing/data_$(lat)_$(lon)")

            prec, tmin, tmax, wind = read_vic_forcing(file_src)

            prec_final += prec/ncells
            tmin_final += tmin/ncells
            tmax_final += tmax/ncells
            wind_final += wind/ncells

        end

        lat = @sprintf("%0.5f", mean(df_tmp[:lat]))
        lon = @sprintf("%0.5f", mean(df_tmp[:lon]))

        file_dst = joinpath(target_folder, "forcing/data_$(lat)_$(lon)")

        write_vic_forcing(file_dst, prec_final, tmin_final, tmax_final, wind_final)

    end

end

"""
Average soil parameters, stored in a dataframe, for a grid resolution
using metadata information stored in wsh_single.
"""
function average_soilparams(soil_param, wsh_single, resolution)

    param_all = convert(Array{Float64,2}, soil_param)

    iboxes = getfield(wsh_single, resolution)
    gridcel = getfield(wsh_single, :gridcel)

    param_tmp = zeros(length(unique(iboxes)), size(param_all, 2))

    for ibox in unique(iboxes)

        isel = findin(param_all[:, 2], gridcel[iboxes .== ibox])

        for icol in 1:size(param_all, 2)
            param_tmp[ibox, icol] = mean(param_all[isel, icol])
        end

    end

    param_tmp[:, 2] .= floor.(param_tmp[:, 2])

    df_final = DataFrame(param_tmp)
    names!(df_final, names(soil_param))

    sort!(df_final, :gridcel)

    return df_final

end

"""
Average vegetation parameters, stored in a dataframe, for a grid resolution
using metadata information stored in wsh_single.
"""
function average_vegparams(veg_param, wsh_single, resolution)

    param_all = convert(Array{Float64, 2}, veg_param)

    iboxes = getfield(wsh_single, resolution)
    gridcel = getfield(wsh_single, :gridcel)

    param_tmp = []

    for ibox in unique(iboxes)

        isel = findin(param_all[:, 1], gridcel[iboxes .== ibox])

        param_mean = average_vegparams(param_all, isel)

        if isempty(param_tmp)
            param_tmp = param_mean
        else
            param_tmp = [param_tmp; param_mean]
        end

    end

    df_final = DataFrame(param_tmp)

    names!(df_final, names(veg_param))

    sort!(df_final, :gridcel)

    return df_final

end

"""
Average vegetation parameters, stored in an array, for a set
of gridcells.
"""
function average_vegparams(param_all, isel)

    # Find parameters belonging to one grid box

    param_box = param_all[isel, :]

    gridcel_mean = floor(mean(param_box[:, 1]))

    # Find vegetation classes within the grid box

    veg_classes = unique(param_box[:, 3])

    Nveg = length(veg_classes)

    # Average parameters over vegetation classes within the box

    param_mean = zeros(Nveg, 20)

    irow = 1

    for veg_class in veg_classes

        # Get parameters for one vegetation class

        iveg = find(param_box[:, 3] .== veg_class)

        param_veg = param_box[iveg, :]

        # Compute an area  weighted average of the parameters

        weights = repmat(param_veg[:, 4], 1, 16)

        tmp_mean = sum(param_veg[:, 5:end] .* weights, 1) ./ sum(weights, 1)

        param_mean[irow, 1] = gridcel_mean
        param_mean[irow, 2] = Nveg
        param_mean[irow, 3] = veg_class
        param_mean[irow, 4] = sum(param_veg[:, 4])
        param_mean[irow, 5:end] = tmp_mean

        irow += 1

    end

    # Scale area fractions to one

    param_mean[:, 4] = param_mean[:, 4] / sum(param_mean[:, 4])

    return param_mean

end


text:
# [src/rkfixed.jl]

"""
Runge-Kutta fixed-step integration using the specified weights, nodes,
and Runge-Kutta matrix (or the Runge-Kutta 4th order "3/8" method by
default).

Implements numerical propagation of an ordinary differential equation
from some initial value over the desired range. This function is similar
to MATLAB's variable-step ODE propagators (e.g., ode45), but uses a
fixed step method. This is useful either when one knows an appropriate
step size or when a process is interrupted frequently (ode45 and the
similar functions in MATLAB will always make at least a certain number of
steps between ts(1) and ts(2), which may be very many more than are
necessary).

This function is generic for all fixed-step Runge-Kutta methods. That is,
any fixed-step Runge-Kutta propagator can be created by passing the
weightes, nodes, and Runge-Kutta matrix (together, the Butcher tableau)
into this function. See the example below.

    t, x = rkfixed(ode, ts, x0, dt);
    t, x = rkfixed(ode, ts, x0, dt, a, b, c);
    t, x = rkfixed(ode, ts, x0, options);
    t, x = rkfixed(ode, ts, x0, options, a, b, c);

Inputs:

    ode     Ordinary differential equation function
    ts      Time span, [t_start, t_end]
    x0      Initial state (column vector)
    options One can specify an options structure instead of dt
                so that this function is compatible with ode45 and its ilk. The
                only valid fields are MaxStep (the time step) and OutputFcn
    a       Runge-Kutta matrix
    b       Weights
    c       Nodes

Outputs:

    t       Time history
    x       State history, with each row containing the state corresponding to
                the time in the same row of t.

Example:

    # Simulate an undamped harmonic oscillator for 10s with a 0.1s time
    #    step, starting from an initial state of [1; 0] using RK 4th order
    #    integration (via the Butcher tableau specified by a, b, and c). This
    #    is exactly the same as the rk4 function

    a = [ 0    0    0  0;
          0.5  0    0  0;
          0    0.5  0  0;
          0    0    1  0];

    b = [ 1  2  2  1]/6;
    c = [ 0  0.5  0.5  1];

    t, x = rkfixed( (t,x) -> [-x[2]; x[1]], [0 10], [1; 0], 0.1, a, b, c);
    plot(t, x);

See "Runge-Kutta methods" on Wikipedia for discussion of the Butcher
tableau (a, b, and c).

See also: odehybrid, ode45, odeset, rk4.

Online doc: http://www.anuncommonlab.com/doc/odehybrid/rkfixed.html

Copyright 2014 An Uncommon Lab
"""
function rkfixed(ode, ts, x0, options::ODE_SET, a, b, c)

    dt = options.MaxStep
    if isempty(dt)
        @error "RKFixed: Specify the time step with the 'MaxStep' option."
    end

    # Time history
    t = ts[1]:dt:ts[end]

    # Go up to the end
    if abs(t[end] - ts[end]) > abs(dt * 0.0000001) # != ts[end]  <-  '≂' led to some numerical issues
        t = [t; ts[end]]
    end

    ns = maximum(size(t));            # Number of samples
    nx = length(x0);                  # Number of states
    x  = [x0[:]'; zeros(ns-1, nx)];   # State history
    xk = x0[:];                       # Current state

    s = maximum(size(b));             # Length of weights
    d = zeros(nx, s);                 # Matrix of derivatives

    # # If the user provided an OutputFcn, use it.
    if (options.OutputFcn != [])
        options.OutputFcn(ts, x0, "init")
    end

    # Propagate.
    for k = 1:ns-1

        # The last sample may be cut short.
        if k == ns-1
            dt = t[k+1] - t[k];
        end

        # Current time
        tk = t[k];

        # Calculate derivatives.
        d[:, 1] = dt * ode(tk, xk);
        for z = 2:s
            dxk = d[:, 1:z-1] * a[z, 1:z-1]
            d[:, z] = dt * ode(tk + c[z] * dt, xk + dxk);
        end

        # Update the state.
        for z = 1:s
            xk = xk + b[z] * d[:, z];
        end

        # Store.
        x[k+1, :] = xk[:]'

        # # If the user provided an OutputFcn, use it.
        if (options.OutputFcn != [])

            if (options.OutputFcn(t[k+1], xk, ""))
                t = t[1:k+1]
                x = x[1:k+1, :]
                break
            end
        end

    end

    # If the user provided an OutputFcn, use it.
    if (options.OutputFcn != [])
        options.OutputFcn([], [], "done")
    end

    return t, x
end

function rkfixed(ode, ts, x0, options::ODE_SET)
    """
        Wrapper function for when an RK tableau is not provided. Defaults to RK 3/8
    """
    # Default to Runge and Kutta's 3/8 formulation (4th order).
    a = [ 0    0 0 0;
          1/3  0 0 0;
         -1/3  1 0 0;
          1   -1 1 0];

    b = [1 3 3 1]/8;

    c = [0 1/3 2/3 1];

    return rkfixed(ode, ts, x0, options, a, b, c)
end

function rkfixed(ode, ts, x0, dt, a, b, c)
    """
        Wrapper function for when a time step (dt) is provided rather than an ODE SET.
            Generates the default ODE SET and updates the MaxStep to dt
    """

    options = ODE_SET()
    options.MaxStep = dt

    return rkfixed(ode, ts, x0, options, a, b, c)
end

function rkfixed(ode, ts, x0, dt)
    """
        Wrapper function for when no RK tableau is provided and a time step (dt) is provided rather than an ODE SET.
            Defaults to RK 3/8 and the default ODE SET with an updated MaxStep
    """

    # Default to Runge and Kutta's 3/8 formulation (4th order).
    a = [ 0    0 0 0;
          1/3  0 0 0;
         -1/3  1 0 0;
          1   -1 1 0];

    b = [1 3 3 1]/8;

    c = [0 1/3 2/3 1];

    return rkfixed(ode, ts, x0, dt, a, b, c)
end


text:
"""
```
@ode_def name begin
    differential equation
end parameters :: ODEFunction
```

## Definition of the Domain-Specific Language (DSL)

A helper macro is provided to make it easier to define a `ParameterizedFunction`,
and it will symbolically compute a bunch of extra functions to make the differential
equation solvers run faster. For example, to define the previous `LotkaVolterra`,
you can use the following command:

```julia
f = @ode_def LotkaVolterra begin
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
end a b c d
```

or you can define it anonymously:

```julia
f = @ode_def begin
    dx = a*x - b*x*y
    dy = -c*y + d*x*y
end a b c d
```

`@ode_def` uses ModelingToolkit.jl internally and returns an `ODEFunction` with the
extra definitions (Jacobian, parameter Jacobian, etc.) defined through the MTK
symbolic tools.
"""
macro ode_def(name,ex,params...)
    opts = Dict{Symbol,Bool}(
    :build_tgrad => true,
    :build_jac => true,
    :build_expjac => false,
    :build_invjac => false,
    :build_invW => false,
    :build_hes => false,
    :build_invhes => false,
    :build_dpfuncs => true)
    name isa Expr ? ode_def_opts(gensym(),opts,__module__,name,ex,params...) :
        ode_def_opts(name,opts,__module__,ex,params...)
end

"""
```
@ode_def_bare name begin
    differential equation
end parameters :: ODEFunction
```

Like `@ode_def` but the `opts` options are set so that no symbolic functions are generated.
See the `@ode_def` docstring for more details.
"""
macro ode_def_bare(name,ex,params...)
    opts = Dict{Symbol,Bool}(
    :build_tgrad => false,
    :build_jac => false,
    :build_expjac => false,
    :build_invjac => false,
    :build_invW => false,
    :build_hes => false,
    :build_invhes => false,
    :build_dpfuncs => false)
    name isa Expr ? ode_def_opts(gensym(),opts,__module__,name,ex,params...) :
        ode_def_opts(name,opts,__module__,ex,params...)
end

"""
```
@ode_def_all name begin
    differential equation
end parameters :: ODEFunction
```

Like `@ode_def` but the `opts` options are set so that all possible symbolic functions are generated.
See the `@ode_def` docstring for more details.
"""
macro ode_def_all(name,ex,params...)
    opts = Dict{Symbol,Bool}(
    :build_tgrad => true,
    :build_jac => true,
    :build_expjac => false,
    :build_invjac => false,
    :build_invW => true,
    :build_hes => false,
    :build_invhes => false,
    :build_dpfuncs => true)
    name isa Expr ? ode_def_opts(gensym(),opts,__module__,name,ex,params...) :
        ode_def_opts(name,opts,__module__,ex,params...)
end

max_stars_repo_name: JuliaDiffEq/ParameterizedFunctions

max_issues_repo_name: JuliaDiffEq/ParameterizedFunctions

max_forks_repo_name: JuliaDiffEq/ParameterizedFunctions


text:
using Pkg;
Pkg.activate("./")
Pkg.instantiate()
using TimeseriesSurrogates, BenchmarkTools, Plots, Distributions, DynamicalSystems, Random, StatsPlots;
gr()

# Here, we compute the mean time to generate a single surrogate for each of the methods
# implemented in TimeseriesSurrogates.jl at the time of submission to JOSS.
# We'll use the same test process as in Lancaster et al, and generate a 2000-pt long
# time series with additive noise from a normal distribution with zero mean and
# standard deviation 0.3.
f(t) = sin(2π*(t + 0.5*sin(2π*t/10))/3)
npts = 2000

# Add some noise to the signal.
rng = Random.MersenneTwister(1234)
t = collect(1:npts) .+ rand(rng, npts)
x = f.(t) .+ rand(rng, Normal(0, 0.3), npts)

# Write to file, to import in matlab
using DelimitedFiles
writedlm("./paper/timing/code/data.csv",  x, ',')

# Pre-initialize all different surrogate types
# Shuffle based surrogates
rs = surrogenerator(x, RandomShuffle(), rng)
rp = surrogenerator(x, RandomFourier(true), rng)
aaft = surrogenerator(x, AAFT(), rng)
iaaft = surrogenerator(x, IAAFT(), rng)

# Pseudoperiodic surrogates (PPS)
τ = estimate_delay(x, "ac_min")
d = findmin(delay_fnn(x, τ))[2]
ρ = noiseradius(x, d, τ, 0.025:0.025:0.5)
pps = surrogenerator(x, PseudoPeriodic(d, τ, ρ, true))

# Collect surrogate generators and run once to trigger precompilation
surrogens = [rs, rp, aaft, iaaft, pps];
surrogates = [s() for s in surrogens];

# Define labels for plot.
methods = ["RandomShuffle()", "RandomFourier(true)", "AAFT()", "IAAFT()", "PseudoPeriodic()"];
n_methods = length(methods)

function generate_surrs(surrogenerator, n)
    for i = 1:n
        surrogenerator()
    end
end

"""
    time_n(methods; n = 100)

Compute `n` surrogates using each pre-initialized `Surrogenerator` in `methods` and compute the
mean time to generate a surrogate.
"""
function time_n(methods; n = 100)
    benchmarks = zeros(length(methods))
    for (i, surrogen) in enumerate(methods)
        benchmarks[i] = @belapsed generate_surrs($surrogen, $n)
    end

    return benchmarks ./ n
end
time_n(surrogens, n = 1);

# Compute average time to generate surrogate over 30 realizations of each surrogate type
timings = time_n(surrogens, n = 100);

# [rp, ft, aaft, iaaft, pps]; these were generated in `surrogate.m`
matlab_timings = [2.0335e-04, 6.3073e-05, 1.4029e-04, 0.0048, 0.1482]

# Plot the results
x = repeat(methods, outer = 2)
gp = repeat(["TimeseriesSurrogates.jl", "MATLAB"], inner = 5)
c = repeat(["Black", "Red"], inner = 5)

groupedbar(
    [timings matlab_timings],
    group = gp,
    c = c,
    xrotation = 45,
    ylabel = "Mean runtime (s)",
    yaxis = :log10,
    yticks = [10.0^(i) for i = -6:1],
    xticks = (1:5, ["rp", "ft", "aaft", "iaaft", "pps"]),
    legend = :topleft,
)

savefig("timings.png")


text:
# ---
# title: 1363. Largest Multiple of Three
# id: problem1363
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Dynamic Programming
# link: <https://leetcode.com/problems/largest-multiple-of-three/description/>
# hidden: true
# ---
#
# Given an integer array of `digits`, return the largest multiple of **three**
# that can be formed by concatenating some of the given digits in any order.
#
# Since the answer may not fit in an integer data type, return the answer as a
# string.
#
# If there is no answer return an empty string.
#
#
#
# **Example 1:**
#
#
#
#     Input: digits = [8,1,9]
#     Output: "981"
#
#
# **Example 2:**
#
#
#
#     Input: digits = [8,6,7,1,0]
#     Output: "8760"
#
#
# **Example 3:**
#
#
#
#     Input: digits = [1]
#     Output: ""
#
#
# **Example 4:**
#
#
#
#     Input: digits = [0,0,0,0,0,0]
#     Output: "0"
#
#
#
#
# **Constraints:**
#
#   * `1 <= digits.length <= 10^4`
#   * `0 <= digits[i] <= 9`
#   * The returning answer must not contain unnecessary leading zeros.
#
#
## @lc code=start
using LeetCode

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


text:
### routiens that exploit matrix structure for computational advantages.

"""
Assumes diagonal entries of L are all positive.
logdettriangular(L::AbstractMatrix{T})::T
"""
function logdettriangular(L::AbstractMatrix{T})::T where T <: Real
    return sum(log(L[i,i]) for i = 1:size(L)[1])
end

### routines for matrix transforms and operations.

"""
evalinverse!( Y::Matrix{T}.
                inv_E::Matrix{T},
                f::Vector{T},
                h::T)

Computes Y := inv(M), where M := [E f; f' h], where E is DxD, f is Dx1, h is 1x1.
Uses the Matrix inversion lemma; see Theorem 4.3.2 from Murphy's book.

Test:
M = Utilities.generaterandomposdefmat(3)
E = M[1:2,1:2]
f = vec(M[1:2,3])
h = M[3,3]
inv_E = inv(E)

Y = zeros(3,3)
evalinverse!( Y,
                inv_E,
                f,
                h)
println("inv(M)-Y = ", inv(M)-Y)
"""
function evalinverse!(  Y::Matrix{T},
                        inv_E::Matrix{T},
                        f::Vector{T},
                        h::T) where T <: Real
    #
    D = length(f)
    @assert size(inv_E) == (D,D)
    @assert size(Y) == (D+1,D+1)

    # compute the Schur complement of M w.r.t. E.
    schur = h - dot(f,inv_E*f)

    b = inv_E*f

    Y[1:D,1:D] = inv_E + (b*b')./schur
    Y[1:D, D+1] = -b ./ schur
    Y[D+1, 1:D] = -b ./ schur
    Y[D+1, D+1] = one(T)/schur

    return nothing
end

function evalinverse(   inv_E::Matrix{T},
                        f::Vector{T},
                        h::T)::Matrix{T} where T <: Real
    #
    Y = Matrix{T}(undef,size(inv_E))
    evalinverse!(Y, inv_E, f, h)

    return Y
end

function forcesymmetric(A::Matrix{T})::Matrix{T} where T <: Real
    return (A+A')./2
end

# replace the upper triangular entires with the lower ones.
function replaceupperbylowertriangular!(A::Matrix{T})::Nothing where T <: Real
    N = size(A,1)
    @assert size(A,2) == size(A,1)
    for j = 1:N
        for i = j+1:N
            A[j,i] = A[i,j]
        end
    end

    return nothing
end

function naivesqrtpsdmatrix(A)
    @assert LinearAlgebra.isposdef(A)
    s,Q = LinearAlgebra.eigen(A)

    s2 = sqrt.(s)

    B = Q'*LinearAlgebra.Diagonal(s2)*Q
    B = forcesymmetric(B)

    return B
end

function array2matrix(X::Array{Vector{T},L})::Matrix{T} where {T,L}

    N = length(X)
    D = length(X[1])

    out = Matrix{T}(undef,D,N)
    for n = 1:N
        out[:,n] = X[n]
    end

    return out
end

"""
    makeblockdiagonalmatrix(A::Vector{Matrix{T}})::Matrix{T}
"""
function makeblockdiagonalmatrix(A::Vector{Matrix{T}})::Matrix{T} where T <: Real
    K = length(A)

    N_rows = 0
    N_cols = 0

    for i = 1:K
        N_rows += size(A[i], 1)
        N_cols += size(A[i], 2)
    end

    out = zeros(T, N_rows, N_cols)
    st_rows = 1
    st_cols = 1

    for k = 1:K
        M,N = size(A[k])

        out[st_rows:st_rows+M-1,st_cols:st_cols+N-1] = A[k]

        st_rows += M
        st_cols += N
    end
    return out
end

"""
    invertblockmatrix!( out::Vector{Matrix{T}},
                    src::Vector{Matrix{T}})::Nothing
"""
function invertblockmatrix!(out::Vector{Matrix{T}},
                            src::Vector{Matrix{T}})::Nothing where T <: Real
    @assert length(out) == length(src)
    for i = 1:length(src)
        out[i] = inv(src[i])
    end
    return nothing
end

"""
    blockmatrixmatrixproduct(   A::Matrix{T},
                                B_set::Vector{Matrix{T}})::Matrix{T}

does A*makeblockdiagonalmatrix(B_set).

# Example:
```
D = 30
p = 30
B2_set = collect(randn(D,D) for i=1:p)
B2 = makeblockdiagonalmatrix( B2_set )
B = randn(D*p,D*p)
norm(B*B2 - blockmatrixmatrixproduct(B,B2_set))
```
"""
function blockmatrixmatrixproduct(  A::Matrix{T},
                                    B_set::Vector{Matrix{T}})::Matrix{T} where T <: Real
    K = length(B_set)

    M,N = size(B_set[1])
    @assert K*N == size(A,2)

    st = 1
    out = Matrix{T}(undef, size(A,1), K*M)
    for i = 1:length(B_set)
        # parse.
        st = (i-1)*size(B_set[i],1) + 1
        W = view(A, :, st:st+N-1)

        # matrix-vector product.
        out[:,st:st+N-1] = W*B_set[i]
    end

    return out
end

"""
    blockmatrixvectorproduct(   B_set::Vector{Matrix{T}},
                                v::Vector{T})::Vector{T}

does makeblockdiagonalmatrix(B_set)*v.

# Example:
```
D = 30
p = 30
B2_set = collect(randn(D,D) for i=1:p)
B2 = makeblockdiagonalmatrix( B2_set )
b = randn(D*p)
norm(B2*b - blockmatrixvectorproduct(B2_set,b))
```
"""
function blockmatrixvectorproduct(  B_set::Vector{Matrix{T}},
                                    v::Vector{T})::Vector{T} where T <: Real
    K = length(B_set)

    M,N = size(B_set[1])
    @assert K*N == length(v)

    st = 1
    out = Array{T}(K*M)
    for i = 1:length(B_set)
        # parse.
        st = (i-1)*size(B_set[i],2) + 1
        w = view(v, st:st+N-1)

        # matrix-vector product.
        out[st:st+N-1] = B_set[i]*w
    end

    return out
end

#### symmetric matrix efficient querying.

"""
    readuppertriangular(i::Int, j::Int, N::Int, v::Vector{T})::T

Query the (i-j)-th entry of an upper triangular matrix that is stored in v.

```
A = randn(7,7)

v = packageuppertriangle(A)

function packagemat(v, N)
    B = zeros(Float64, N, N)

    for j = 1:N
        for i = 1:j
            B[i,j] = readuppertriangular(i, j, N, v)
        end
    end

    return B
end

C = packagemat(v, size(A,1))
display(A-C)
```
"""
function readuppertriangular(i::Int, j::Int, N::Int, v::Vector{T})::T where T

    linear_ind = div(j*(j-1), 2) + i

    #Printf.@printf("(i,j) = (%d,%d), linear index = %d, offset = %d\n", i, j, linear_ind, offset)
    return v[linear_ind]
end

"""
    readsymmetric(i::Int, j::Int, N::Int, v::Vector{T})::T

Queries a symmetric matrix stored in a vector format.
"""
function readsymmetric(i::Int, j::Int, N::Int, v::Vector{T})::T where T
    if i < j
        return readuppertriangular(i, j, N, v)
    end

    return readuppertriangular(j, i, N, v)
end

"""
    writeuppertriangular!(v::Vector{T}, i::Int, j::Int, N::Int, val::T)::Nothing

Writes a value to an upper triangular matrix storage.
"""
function writeuppertriangular!(v::Vector{T}, i::Int, j::Int, N::Int, val::T)::Nothing where T

    linear_ind = div(j*(j-1), 2) + i
    v[linear_ind] = val

    return nothing
end

"""

```
D = 3
v = Vector{Float64}(undef, 6)
@time Utilities.writesymmetric!(v, 2, 3, 3, 1.23)

function fillv!(v, D)
    c = 0.0
    for j = 1:D
        for i = 1:j
            Utilities.writesymmetric!(v, i, j, D, c)
            c += 1
        end
    end

end

fillv!(v, D)
display(v)
```
"""
function writesymmetric!(v::Vector{T}, i::Int, j::Int, N::Int, val::T)::Nothing where T
    if i < j
        writeuppertriangular!(v, i, j, N, val)

        return nothing
    end

    writeuppertriangular!(v, j, i, N, val)

    return nothing
end

"""
    packageuppertriangle(L::Matrix{T})::Vector{T}

Extracts the upper triangular portion of a matrix. Includes diagonal entries.
"""
function packageuppertriangle(L::Matrix{T})::Vector{T} where T
    N = size(L,1)

    a = Vector{T}(undef, div(N*(N+1),2) )

    p = 1
    for j = 1:N
        for i = 1:j
            a[p] = L[i,j]
            p += 1
        end
    end

    return a
end

max_stars_repo_name: RoyCCWang/RWJuliaUtilities

max_issues_repo_name: RoyCCWang/RWJuliaUtilities

max_forks_repo_name: RoyCCWang/RWJuliaUtilities


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 = "DGGRID7"
version = v"0.9.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/allixender/DGGRID.git", "0dabdfd3f8fd6ac27626064a94c483e610d578bb")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/DGGRID/src/
if [[ "${target}" == x86_64-linux-musl ]]; then
    # Remove libexpat to avoid it being picked up by mistake
    rm /usr/lib/libexpat.so*
fi
make -j${nproc} CCOMP="${CC}" CPPCOMP="${CXX}"
cp "apps/dggrid/dggrid${exeext}" "${bindir}/."
"""

# 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("dggrid", :dggrid)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="GDAL_jll", uuid="a7073274-a066-55f0-b90d-d619367d196c"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
#module Product_quantization
using Clustering

#for computing distances
using Distances

#for computing the error with k-means
include("../utils.jl")
include("../opq/kmeans.jl")

# Quantize using PQ codebooks
function quantize_pq{T <: AbstractFloat}(
  X::Matrix{T},         # d-by-n. Data to encode
  C::Vector{Matrix{T}}, # codebooks
  V::Bool=false)        # whether to print progress

  d, n = size( X )
  m    = length( C )
  h    = size( C[1], 2 )

  B       = zeros(Int, m, n) # codes
  subdims = splitarray( 1:d, m )

  # auxiliary variables for update_assignments! function
  costs     = zeros(Float32, n)
  counts    = zeros(Int, h)
  to_update = zeros(Bool, h)
  unused    = Int[]

  for i = 1:m
    if V print("Encoding on codebook $i / $m... ") end

    # Find distances from X to codebook
    dmat = pairwise( SqEuclidean(), C[i], X[subdims[i],:] )
    dmat = convert(Array{T}, dmat)
    update_assignments!( dmat, true, view(B,i,:), costs, counts, to_update, unused )

    if V println("done"); end
  end
  return convert(Matrix{Int16}, B) # return the codes
end

# Train product quantization
function train_pq{T <: AbstractFloat}(
  X::Matrix{T},  # d-by-n. Data to learn codebooks from
  m::Integer,    # number of codebooks
  h::Integer,    # number of entries per codebook
  V::Bool=false) # whether to print progress

  d, n = size( X )

  C = Vector{Matrix{T}}(m); # codebooks

  B        = zeros(Int16, m, n); # codes
  subdims  = splitarray(1:d, m); # subspaces

  for i = 1:m
    if V print("Working on codebook $i / $m... "); end
    cluster = kmeans( X[ subdims[i],: ], h, init=:kmpp);
    C[i], B[i,:] = cluster.centers, cluster.assignments;

    if V
      subdim_cost = cluster.totalcost ./ n;
      nits        = cluster.iterations;
      converged   = cluster.converged;

      println("done.");
      println("  Ran for $nits iterations");
      println("  Error in subspace is $subdim_cost");
      println("  Converged: $converged")
    end
  end

  error = qerror_pq( X, B, C )
  return C, B, error
end

#end # module PQ

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
# keywords
export
    ADD, ALL, ALTER, AND, AS, ASC, BEGIN, BETWEEN, BIGINT, BY, COMMIT, COLUMN, CONSTRAINT, CREATE, DATABASE, DEFAULT, DELETE, DESC, DISTINCT, DROP, EXCEPT, EXECUTE, EXISTS, FOREIGN, FROM, FULL, GROUP,
    HAVING, IN, INDEX, IF, INNER, INSERT, INTERSECT, INTO, IS, JOIN, KEY, LEFT, LIKE, LIMIT, NULL, OFF, OFFSET, ON, OR, ORDER, OUTER, OVER,
    PARTITION, PREPARE, PRIMARY, RECURSIVE, REFERENCES, RELEASE, RIGHT, ROLLBACK, SAVEPOINT, SELECT, SET, TABLE, TO, TRANSACTION, TRIGGER, UNION, UPDATE, USING, VALUES, WHERE, WITH

# aggregate functions
export AVG, COUNT, EVERY, MAX, MIN, NOT, SOME, SUM

# ranking functions
export DENSE_RANK, RANK, ROW_NUMBER

export Octo, Repo, Schema, Pretty
export PlaceHolder, Raw, Enclosed
export @sql_keywords, @sql_functions
export from, as
export to_sql


text:
"""
    LogNormal

* Parametrization
parametrized in \\mu and \\sigma^2

* Score

* Fisher Information

* `time_varying_params` map.

* Default link
"""
LogNormal

function score!(score_til::Matrix{T}, y::T, ::Type{LogNormal}, param::Matrix{T}, t::Int) where T
    score_til[t, 1] = (log(y) - param[t, 1])/param[t, 2]
    score_til[t, 2] = -(0.5/param[t, 2]) * (1 - ((log(y) - param[t, 1])^2)/param[t, 2])
    return
end

function fisher_information!(aux::AuxiliaryLinAlg{T}, ::Type{LogNormal}, param::Matrix{T}, t::Int) where T
    aux.fisher[1, 1] = 1/(param[t, 2])
    aux.fisher[2, 2] = 1/(2*(param[t, 2]^2))
    aux.fisher[2, 1] = 0
    aux.fisher[1, 2] = 0
    return
end

function log_likelihood(::Type{LogNormal}, y::Vector{T}, param::Matrix{T}, n::Int) where T
    loglik = -0.5*log(2*pi)*n
    for t in 1:n
        loglik -= log(y[t] * sqrt(param[t, 2])) + 0.5*(log(y[t]) - param[t, 1])^2/param[t, 2]
    end
    return -loglik
end

# Links
function link!(param_tilde::Matrix{T}, ::Type{LogNormal}, param::Matrix{T}, t::Int) where T
    param_tilde[t, 1] = link(IdentityLink, param[t, 1])
    param_tilde[t, 2] = link(LogLink, param[t, 2], zero(T))
    return
end
function unlink!(param::Matrix{T}, ::Type{LogNormal}, param_tilde::Matrix{T}, t::Int) where T
    param[t, 1] = unlink(IdentityLink, param_tilde[t, 1])
    param[t, 2] = unlink(LogLink, param_tilde[t, 2], zero(T))
    return
end
function jacobian_link!(aux::AuxiliaryLinAlg{T}, ::Type{LogNormal}, param::Matrix{T}, t::Int) where T
    aux.jac[1] = jacobian_link(IdentityLink, param[t, 1])
    aux.jac[2] = jacobian_link(LogLink, param[t, 2], zero(T))
    return
end

# utils
function update_dist(::Type{LogNormal}, param::Matrix{T}, t::Int) where T
    # lognormal here is parametrized as sigma^2
    return LogNormal(param[t, 1], sqrt(param[t, 2]))
end

function params_sdm(d::LogNormal)
    pars = Distributions.params(d)
    return (pars[1], pars[2]^2)
end

function num_params(::Type{LogNormal})
    return 2
end


text:
# Parser-lexer definitions

mutable struct Pos
    column
    line
    index
end

mutable struct LTable
    actions
    terminal_states
    tok2index
end

mutable struct Tok
    id
    terminal
    data
    start_pos
    end_pos
end

mutable struct Prod
    left
    right
    pos
    init
end

ltable_new(actions, terminal_states, tok2index) = LTable(actions, terminal_states, tok2index)

pos_new(; column = 1, line = 1, index = 1) = Pos(column, line, index)

Base.isequal(a::Pos, b::Pos) = a.column == b.column && a.line == b.line
Base.print(io::IO, p::Pos) = print(io, "(line: $(p.line), column: $(p.column))")

function tok_new(id; terminal = false, data = nothing, start_pos = nothing, end_pos = nothing)
    return Tok(id, terminal, data, start_pos, end_pos)
end

tok_eq(a, b) = a.id == b.id

# Special tokens
tok_eps() = tok_new("eps", terminal=true)

function tok_end(; start_pos = nothing, end_pos = nothing)
    return tok_new("\$", terminal=true, start_pos=start_pos, end_pos=end_pos)
end

# Only id matters for sets / dicts
Base.hash(a::Tok) = hash(a.id)
Base.isequal(a::Tok, b::Tok) = a.id == b.id
Base.print(io::IO, t::Tok) = print(io, t.id)

prod_new(left::Tok, right::Array{Tok}; pos = 1, init = false) = Prod(left, right, pos, init)

function Base.hash(a::Prod)
    h = hash(a.left)
    h ⊻= hash(a.pos)
    for (i, r) in enumerate(a.right)
        h ⊻= hash(r) + i
    end

    return h
end

Base.isequal(a::Prod, b::Prod) =
    isequal(a.left, b.left) &&
    all([isequal(u, v) for (u, v) in zip(a.right, b.right)]) &&
    isequal(a.pos, b.pos)

function Base.print(io::IO, p::Prod)
    r = [t.id for t in p.right]
    insert!(r, p.pos, ".")
    str = "$(p.left.id) ->" * reduce((a, b) -> "$a $b", r, init="")

    Base.print(io, str)
end

# Default character to table index mapping
ldefault_chartok2index = (tok) -> codepoint(tok)
ldefault_chartok2index_length = 128


text:
module TestValDate

using Test
using TSML
using DataFrames: DataFrame, nrow

function generateXY()
    Random.seed!(123)
    gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2016,1,1)
    gval = Array{Union{Missing,Float64}}(rand(length(gdate)))
    gmissing = 50000
    gndxmissing = Random.shuffle(1:length(gdate))[1:gmissing]
    X = DataFrame(Date=gdate,Value=gval)
    X.Value[gndxmissing] .= missing
    Y = rand(length(gdate))
    (X,Y)
end

const (XX,YY)=generateXY()
const (X1,Y1)=generateXY()

function test_datevalizer()
    dvzr1 = DateValizer(Dict(:dateinterval=>Dates.Hour(1)))
    dvzr2 = DateValizer(dvzr1.model)
    @test dvzr1.model[:dateinterval] == dvzr2.model[:dateinterval]
    res = fit_transform!(dvzr2,XX)
    @test sum(ismissing.(res.Value)) == 0
    @test sum(X1.Value .!== XX.Value) == 0
    @test sum(Y1 .!== YY) == 0
    @test abs(round(sum(res.Value),digits=2) - 8798.2) <= 100
    @test nrow(dvzr2.model[:medians]) == 24
    res = fit_transform(dvzr2,XX)
    @test sum(ismissing.(res.Value)) == 0
    @test sum(X1.Value .!== XX.Value) == 0
    @test sum(Y1 .!== YY) == 0
    @test abs(round(sum(res.Value),digits=2) - 8798.2) <= 100
    @test nrow(dvzr2.model[:medians]) == 24
end
@testset "DateValizer: Fill missings with medians" begin
    test_datevalizer()
end

function test_datevalgator()
    maxlim = 200
    dtvl = DateValgator(Dict(:dateinterval=>Dates.Hour(1)))
    res = fit_transform!(dtvl,XX)
    @test abs(sum(ismissing.(res.Value)) - 4466) <=200
    @test abs(round(sum(skipmissing(res.Value)),digits=2) - 6556.17) <= maxlim
    @test sum(X1.Value .!== XX.Value) == 0
    @test sum(Y1 .!== YY) == 0
    res = fit_transform(dtvl,XX)
    @test abs(sum(ismissing.(res.Value)) - 4466) <= maxlim
    @test abs(round(sum(skipmissing(res.Value)),digits=2) - 6556.17) <= maxlim
    @test sum(X1.Value .!== XX.Value) == 0
    @test sum(Y1 .!== YY) == 0

    dtvlmean = DateValgator(Dict(
	  :dateinterval=>Dates.Hour(1),
	  :aggregator => :mean))
    res = fit_transform!(dtvlmean,XX)
    @test abs(sum(ismissing.(res.Value)) - 4466) <= maxlim
    @test abs(round(sum(skipmissing(res.Value)),digits=2) - 6557.97) <= maxlim
    res = fit_transform(dtvlmean,XX)
    @test abs(sum(ismissing.(res.Value)) - 4466) <= maxlim
    @test abs(round(sum(skipmissing(res.Value)),digits=2) - 6557.97) <= maxlim

    dtvlmax = DateValgator(Dict(
	  :dateinterval=>Dates.Hour(1),
	  :aggregator => :maximum))
    res = fit_transform!(dtvlmax,XX)
    @test abs(sum(ismissing.(res.Value)) - 4466) <= maxlim
    @test abs(round(sum(skipmissing(res.Value)),digits=2) - 7599.95) <= maxlim
    res = fit_transform(dtvlmax,XX)
    @test abs(sum(ismissing.(res.Value)) - 4466) <= maxlim
    @test abs(round(sum(skipmissing(res.Value)),digits=2) - 7599.95) <= maxlim

    dtvlmin = DateValgator(Dict(
	  :dateinterval=>Dates.Hour(1),
	  :aggregator => :minimum))
    res = fit_transform!(dtvlmin,XX)
    @test abs(sum(ismissing.(res.Value)) - 4466) <= maxlim
    @test abs(round(sum(skipmissing(res.Value)),digits=2) - 5516.92) <= maxlim
    res = fit_transform(dtvlmin,XX)
    @test sum(ismissing.(res.Value)) >= 4000
    @test abs(round(sum(skipmissing(res.Value)),digits=2) - 5516.92) <= maxlim
end
@testset "DateValgator: aggregate by timeperiod without filling missings" begin
    test_datevalgator()
end

function test_datevalnner()
    dnnr = DateValNNer(Dict(
	  :dateinterval=>Dates.Hour(1),
	  :nnsize=>10,
	  :missdirection => :symmetric,
	  :strict=>true,
	  :aggregator => :mean))
    res=fit_transform!(dnnr,XX)
    @test sum(size(res) .== (17521,2)) == 2
    @test round(sum(res.Value),digits=2) > 8000.00
    res=fit_transform(dnnr,XX)
    @test sum(size(res) .== (17521,2)) == 2
    @test round(sum(res.Value),digits=2) > 8000.00

    dnnr = DateValNNer(Dict(
	  :dateinterval=>Dates.Hour(1),
	  :nnsize=>10,
	  :missdirection => :symmetric,
	  :strict=>true,
	  :aggregator => :maximum))
    res=fit_transform!(dnnr,XX)
    @test sum(size(res) .== (17521,2)) == 2
    @test round(sum(res.Value),digits=2) > 10000.00
    res=fit_transform(dnnr,XX)
    @test sum(size(res) .== (17521,2)) == 2
    @test round(sum(res.Value),digits=2) > 10000.00
    dnnr = DateValNNer(Dict(
	  :dateinterval=>Dates.Hour(1),
	  :nnsize=>10,
	  :missdirection => :symmetric,
	  :strict=>true,
	  :aggregator => :minimum))
    res=fit_transform!(dnnr,XX)
    @test sum(size(res) .== (17521,2)) == 2
    @test round(sum(res.Value),digits=2) > 7000.00
    res=fit_transform(dnnr,XX)
    @test sum(size(res) .== (17521,2)) == 2
    @test round(sum(res.Value),digits=2) > 7000.00

    dnnr = DateValNNer(Dict(
	  :dateinterval=>Dates.Hour(25),
	  :nnsize=>10,
	  :missdirection => :forward,
	  :strict=>true))
    res=fit_transform!(dnnr,XX)
    @test sum(size(res) .== (701,2)) == 2
    @test round(sum(res.Value),digits=2) > 300.00
    res=fit_transform(dnnr,XX)
    @test sum(size(res) .== (701,2)) == 2
    @test round(sum(res.Value),digits=2) > 300.00

    dnnr.model[:missdirection] = :reverse
    res=TSML.transform!(dnnr,XX)
    @test sum(size(res) .== (701,2)) == 2
    @test round(sum(res.Value),digits=2) > 300.00
    dnnr.model[:dateinterval]=Dates.Hour(1)
    #@test_throws ErrorException res=TSML.transform!(dnnr,XX)
    dnnr.model[:missdirection] = :reverse
    #@test_throws ErrorException res=TSML.transform!(dnnr,XX)
    dnnr.model[:missdirection] = :symmetric
    @test sum(size(TSML.transform!(dnnr,XX)) .== (17521,2)) == 2

    # testing boundaries
    Random.seed!(123)
    dlnr = DateValNNer(Dict(
	    :dateinterval=>Dates.Hour(1),
	    :nnsize=>2,
	    :missdirection => :forward,
	    :strict=>true))
    v1=DateTime(2014,1,1,1,0):Dates.Hour(1):DateTime(2014,1,3,1,0)
    val=Array{Union{Missing,Float64}}(rand(length(v1)))
    x=DataFrame(Date=v1,Value=val)
    x[45:end-1,:Value] .= missing
    fit!(dlnr,x,[])
    @test_throws ErrorException TSML.transform!(dlnr,x)
    @test_throws ErrorException fit_transform(dlnr,x)
    dlnr.model[:missdirection] = :reverse
    res = TSML.transform!(dlnr,x)
    @test sum((size(res)) .== (49,2)) == 2
    @test round(sum(res.Value),digits=2) > 10.00
    dlnr.model[:missdirection] = :forward
    dlnr.model[:strict] = false
    @test sum(ismissing.(TSML.transform!(dlnr,x)[!,:Value])) == 2
    dlnr.model[:missdirection] = :symmetric
    @test sum(ismissing.(TSML.transform!(dlnr,x)[!,:Value])) == 0
    defdnr=DateValNNer(Dict(:strict=>false))
    fit!(defdnr,XX,YY)
    @test sum((size(TSML.transform!(defdnr,XX))) .== (17521,2)) == 2
    m=fit(defdnr,XX,YY)
    @test sum((size(TSML.transform(m,XX))) .== (17521,2)) == 2

    # test with controlled locations of missings
    dlnr = DateValNNer(Dict(:dateinterval=>Dates.Hour(1),
			    :nnsize=>10,:strict=>false,
			    :missdirection => :forward))
    v1=DateTime(2014,1,1,1,0):Dates.Hour(1):DateTime(2014,1,3,1,0)
    val=Array{Union{Missing,Float64}}(collect(1:(length(v1))))
    x=DataFrame(Date=v1,Value=val)
    x[45:end,:Value] .= missing
    x[1:10,:Value] .= missing
    x[20:30,:Value] .= missing
    fit!(dlnr,x,[])
    res = TSML.transform!(dlnr,x)
    @test sum(ismissing.(TSML.transform!(dlnr,x)[!,:Value])) == 6
    dlnr.model[:missdirection] = :reverse
    @test sum(ismissing.(TSML.transform!(dlnr,x)[!,:Value])) == 5
    dlnr.model[:missdirection] = :symmetric
    @test sum(ismissing.(TSML.transform!(dlnr,x)[!,:Value])) == 0
end
@testset "DateValNNer: replace missings with nearest neighbors" begin
    test_datevalnner()
end

function test_dateifier()
  dtr = Dateifier(Dict())
  lower = DateTime(2017,1,1)
  upper = DateTime(2018,1,31)
  dat=lower:Dates.Day(1):upper |> collect
  vals = rand(length(dat))
  x=DataFrame(Date=dat,Value=vals)
  fit!(dtr,x)
  res = TSML.transform!(dtr,x)
  @test sum(size(res) .== (389,8)) == 2
  m=fit(dtr,x)
  res = TSML.transform(m,x)
  @test sum(size(res) .== (389,8)) == 2
  dtr.model[:stride]=2
  res = TSML.transform!(dtr,x)
  @test sum(size(res) .== (194,8)) == 2
end
@testset "Dateifier: extract sliding windows date features" begin
    test_dateifier()
end

function test_matrifier()
  mtr = Matrifier(Dict(:ahead=>24,:size=>24,:stride=>5))
  sz = mtr.model[:size]
  lower = DateTime(2017,1,1)
  upper = DateTime(2017,1,5)
  dat=lower:Dates.Hour(1):upper |> collect
  vals = 1:length(dat)
  x = DataFrame(Date=dat,Value=vals)
  res = fit_transform!(mtr,x)
  @test sum(size(res) .== (10,25)) == 2
  res = fit_transform(mtr,x)
  @test sum(size(res) .== (10,25)) == 2
  mtr.model = Dict(:ahead=>24,:size=>24,:stride=>1)
  res = TSML.transform!(mtr,x)
  @test sum(size(res) .== (50,25)) == 2
  mtr.model = Dict(:ahead=>1,:size=>24,:stride=>12)
  res = TSML.transform!(mtr,x)
  @test sum(size(res) .== (6,25)) == 2
  dtr = Matrifier()
  fit!(dtr,x)
  res = TSML.transform!(dtr,x)
  @test sum(size(res) .== (90,8)) == 2
  m=fit(dtr,x)
  res = TSML.transform(m,x)
  @test sum(size(res) .== (90,8)) == 2
  dtr.model = Dict(:ahead=>-1,:size=>24,:stride=>12)
  @test_throws AssertionError TSML.transform!(dtr,x)
end
@testset "Dateifier: extract sliding windows date features" begin
    test_matrifier()
end

function test_csvreaderwriter()
  inputfile =joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
  outputfile = joinpath(dirname(pathof(TSML)),"../data/testdata_output.csv")
  csvreader = CSVDateValReader(Dict(:filename=>inputfile,:dateformat=>"d/m/y H:M"))
  csvwtr = CSVDateValWriter(Dict(:filename=>outputfile,:dateformat=>"d/m/y H:M"))
  filter1 = DateValgator()
  filter2 = DateValNNer(Dict(:nnsize=>1,:strict=>:false))
  mypipeline = csvreader |> filter1 |> filter2
  res=fit_transform!(mypipeline)
  @test nrow(res) == 8761
  @test ncol(res) == 2
  @test sum(ismissing.(res.Value)) == 0
  @test floor(sum(res.Value)) == 93080.0
  res=fit_transform(mypipeline)
  @test nrow(res) == 8761
  @test ncol(res) == 2
  @test sum(ismissing.(res.Value)) == 0
  @test floor(sum(res.Value)) == 93080.0
  dat = fit_transform!(csvreader)
  res1=fit_transform!(filter1,dat)
  res2=fit_transform!(filter2,res1)
  @test nrow(res2) == 8761
  @test ncol(res2) == 2
  @test sum(ismissing.(res2.Value)) == 0
  @test floor(sum(res2.Value)) == 93080.0
  dat = fit_transform(csvreader)
  res1=fit_transform(filter1,dat)
  res2=fit_transform(filter2,res1)
  @test nrow(res2) == 8761
  @test ncol(res2) == 2
  @test sum(ismissing.(res2.Value)) == 0
  @test floor(sum(res2.Value)) == 93080.0
  mypipeline = csvreader |> filter1 |> filter2 |> csvwtr
  res=fit_transform!(mypipeline)
  @test nrow(res2) == 8761
  @test ncol(res2) == 2
  @test sum(ismissing.(res2.Value)) == 0
  @test floor(sum(res2.Value)) == 93080.0
  @test filesize(csvwtr.model[:filename]) > 209220
  rm(outputfile,force=true)
end
@testset "CSVDateValReaderWriter: reading csv with Date,Value columns" begin
  test_csvreaderwriter()
end

#function test_bzcsvreaderwriter()
#  inputfile =joinpath(dirname(pathof(TSML)),"../data/testdata.csv.bz2")
#  csvreader = BzCSVDateValReader(Dict(:filename=>inputfile,:dateformat=>"d/m/y H:M"))
#  filter1 = DateValgator()
#  filter2 = DateValNNer(Dict(:nnsize=>1))
#  mypipeline = csvreader |> filter1 |> filter2
#  res=fit_transform!(mypipeline)
#  @test nrow(res) == 8761
#  @test ncol(res) == 2
#end
#@testset "BzCSVDateValReader: reading bzipped csv with Date,Value columns" begin
#  test_bzcsvreaderwriter()
#end

function test_statoutputwriter()
  inputfile =joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
  statoutputfile = joinpath(dirname(pathof(TSML)),"../data/testdata_stat.csv")
  csvreader = CSVDateValReader(Dict(:filename=>inputfile,:dateformat=>"d/m/y H:M"))
  csvstatwtr = CSVDateValWriter(Dict(:filename=>statoutputfile))
  statfier =  Statifier(Dict(:processmissing=>true))
  filter1 = DateValgator(Dict(:nnsize=>1))
  filter2 = DateValNNer(Dict(:nnsize=>1,:strict=>false))

  mypipeline = csvreader |> filter1  |> filter2 |> statfier  |> csvstatwtr
  res=fit_transform!(mypipeline)
  @test nrow(res) == 1
  @test ncol(res) == 26
  @test filesize(csvstatwtr.model[:filename]) > 400
  rm(statoutputfile,force=true)

end
@testset "CSVDateValReaderWriter: writing stat output" begin
  test_statoutputwriter()
end

function test_datevalmultinner()
  Random.seed!(123)
  gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2016,1,1)
  gval1 = Array{Union{Missing,Float64}}(rand(length(gdate)))
  gval2 = Array{Union{Missing,Float64}}(rand(length(gdate)))
  gval3 = Array{Union{Missing,Float64}}(rand(length(gdate)))
  gmissing = 50000
  gndxmissing1 = Random.shuffle(1:length(gdate))[1:gmissing]
  gndxmissing2 = Random.shuffle(1:length(gdate))[1:gmissing]
  gndxmissing3 = Random.shuffle(1:length(gdate))[1:gmissing]
  X = DataFrame(Date=gdate,Temperature=gval1,Humidity=gval2,Ozone=gval3)
  X.Temperature[gndxmissing1] .= missing
  X.Humidity[gndxmissing2] .= missing
  X.Ozone[gndxmissing3] .= missing
  dnnr = DateValMultiNNer(Dict(:type=>:linear))
  res = fit_transform!(dnnr,X)
  mysum(x)=sum(skipmissing(x))
  @test mysum.(eachcol(res[:,2:end])) |> sum > 26000.00
  res = fit_transform(dnnr,X)
  @test mysum.(eachcol(res[:,2:end])) |> sum > 26000.00
  dnnr = DateValMultiNNer(Dict(:type=>:knn))
  res = fit_transform!(dnnr,X)
  @test mysum.(eachcol(res[:,2:end]))  |> sum > 26000.00
  res = fit_transform(dnnr,X)
  @test mysum.(eachcol(res[:,2:end]))  |> sum > 26000.00
  dnnr = DateValMultiNNer()
  @test_throws ArgumentError fit!(dnnr,DataFrame(Date=X.Date,Value1=X.Date,Value2=X.Temperature))
  @test_throws ArgumentError fit(dnnr,DataFrame(Date=X.Date,Value1=X.Date,Value2=X.Temperature))
end
@testset "DateValMultiNNer: multicolumn imputation" begin
  test_datevalmultinner()
end

function test_datevallinearimputer()
    Random.seed!(123)
    gdate = DateTime(2014,1,1):Dates.Minute(15):DateTime(2016,1,1)
    gval = Array{Union{Missing,Float64}}(rand(length(gdate)))
    gmissing = 50000
    gndxmissing = Random.shuffle(1:length(gdate))[1:gmissing]
    X = DataFrame(Date=gdate,Value=gval)
    X.Value[gndxmissing] .= missing
    dnnr = DateValLinearImputer()
    @test fit_transform!(dnnr,X) |> x-> sum(x.Value) > 8000.00
    @test fit_transform(dnnr,X) |> x-> sum(x.Value) > 8000.00
end
@testset "DateValLinearImputer: linear imputation" begin
  test_datevallinearimputer()
end

end

max_stars_repo_name: IBM/tsml

max_issues_repo_name: IBM/tsml

max_forks_repo_name: IBM/tsml


text:
module SeismicBSS

include("CICA.jl")
include("Decon.jl")

end # module


text:
export subpopPlot

function subpopPlot( mpSim::MPsim, timeGrid::Vector{T},
    subpops::Subpopulation...; showPlot::Bool=true, savePlot::Bool=false,
    filename::AbstractString="", timeFactor::Real=12.0 ) where T <: Real

    savePlot |= filename != ""

    # Don't do anything if the plot mustn't be shown or saved.
    if !( showPlot || savePlot )
        return
    end  # if !( showPlot || savePlot )

    # Generate the report.
    subpopReport = subpopulationPopReport( mpSim, timeGrid, subpops... )

    if isempty( subpopReport )
        return
    end  # if isempty( subpopReport )

    # Generate the plot.
    plotTitle = string( "Subpopulation plot" )
    plotData = Matrix( subpopReport[:, 2:end] )
    labels = hcat( string.( names( subpopReport )[2:end] )... )
    ymax = max( maximum( plotData ), 20 )

    plt = plot( subpopReport[:, :timePoint] / timeFactor, plotData,
        size=(960, 540), title=plotTitle, labels=labels, lw=2, ylim=[0, ymax] )

    # Show plot if needed.
    if showPlot
        gui( plt )
    end  # if showPlot

    # Save plot if needed.
    if savePlot
        # If the user doesn't specify a filename, a random one is generated.
        if filename == ""
            dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
            filename = string( "Subpopulation plot  ", dateStr )
        end  # if filename == ""

        if !ispath( dirname( filename ) )
            mkpath( dirname( filename ) )
        end  # if !ispath( dirname( filename ) )

        extension = split( filename, "." )
        extension = length( extension ) == 1 ? "" : extension[end]
        filename = string( filename, extension ∈ extensions ? "" : ".svg" )
        savefig( plt, filename )
    end  # if savePlot

end  # subpopPlot( mpSim, timeGrid, subpops, showPlot, savePlot, fileName,
     #   timeFactor )

subpopPlot( mpSim::MPsim, timeRes::Real,
    subpops::Subpopulation...; breakdownType::Symbol=:none, showPlot::Bool=true,
    savePlot::Bool=false, filename::AbstractString="", timeFactor::Real=12.0 ) =
    subpopPlot( mpSim, MP.generateTimeGrid( mpSim, timeRes ), subpops...,
    showPlot=showPlot, savePlot=savePlot, filename=filename,
    timeFactor=timeFactor )