text:
using Test
using Mongoc
using Dates
using UUIDs
using Random

include("../../src/DecompositionDB/getqueries.jl")
include("../../src/DecompositionDB/setqueries.jl")

function construct_test_instances(collection::Mongoc.Collection)
    push_instance!(collection,
                   "case1354",
                   Mongoc.BSON("OPF_ctr" => "path/to/OPF/ctr",
                               "OPF_mat" => "path/to/OPF/mat",
                               "matpower" => "path/to/Matpower"),
                   Mongoc.BSON("OPF" => Mongoc.BSON("opf1" => 1, "opf2" => 2),
                               "graph" => Mongoc.BSON("graph1" => 1, "graph2" => 4))
                  )
    push_instance!(collection,
                   "case118",
                   Mongoc.BSON("OPF_ctr" => "path/to/OPF/ctr",
                               "OPF_mat" => "path/to/OPF/mat",
                               "matpower" => "path/to/Matpower"),
                   Mongoc.BSON("OPF" => Mongoc.BSON("opf1" => 1, "opf2" => 2),
                               "graph" => Mongoc.BSON("graph1" => 1, "graph2" => 4))
                  )
end

function construct_test_decompositions(collection::Mongoc.Collection)
    oids = Vector{String}(undef, 3)
    n, oids[1] = push_decomposition!(collection,
                                  "case1354",
                                  Dict("cliques" => [[1, 2, 3], [4, 5, 6, 7], [8, 7]],
                                       "cliquetree" => [[1, 2], [2, 3]],
                                       "nb_added_edges" => 100),
                                  Mongoc.BSON("clique" => Mongoc.BSON("nb_clique" => 1201, "size_var" => 0.05495)))
    n, oids[2] =  push_decomposition!(collection,
                                   "case1354",
                                   Dict("cliques" => [[1, 2], [2, 3, 4], [3, 4, 5], [5, 6, 7, 8]],
                                        "cliquetree" => [[1, 2], [2, 3], [3, 4]],
                                        "nb_added_edges" => 50),
                                   Mongoc.BSON("clique" => Mongoc.BSON("nb_clique" => 531, "size_var" => 1.4568),
                                               "solver" => Mongoc.BSON("solving_time" => 587, "nb_iter" => 37)))
    n, oids[3] = push_decomposition!(collection,
                                  "case118",
                                  Dict("cliques" => [[1, 2, 3], [4, 5, 6, 7], [8, 7]],
                                       "cliquetree" => [[1, 2], [2, 3]],
                                       "nb_added_edges" => 10),
                                  Mongoc.BSON("clique" => Mongoc.BSON("nb_clique" => 380, "size_var" => 0.1868),
                                              "solver" => Mongoc.BSON("solving_time" => 12.5, "nb_iter" => 58)))
    return oids
end

function get_test_collections()
    client = Mongoc.Client()
    Mongoc.ping(client)
    db = client["choleskytest"]
    Mongoc.delete_many(db["instances"], Mongoc.BSON())
    Mongoc.delete_many(db["decompositions"], Mongoc.BSON())
    return return db["instances"], db["decompositions"]
end

@testset "queries" begin
    instances, decompositions = get_test_collections()
    construct_test_instances(instances)
    oids = construct_test_decompositions(decompositions)
    oid = nothing
    @testset "setqueries.jl" begin
        @testset "push_instance!" begin
            push_instance!(instances, "case118mod")
            @test isinstance(instances, "case118mod")
        end
        @testset "push_decomposition!" begin
            instance_name, oid = push_decomposition!(decompositions,
                                                     "case118mod",
                                                     Dict("cliques" => [[7, 2], [2, 9], [1, 3]],
                                                          "cliquetree" => [[1, 2], [2, 3]],
                                                          "nb_added_edges" => 25))
            dec = getdecomposition(decompositions, instance_name, oid)
            @test dec["cliques"] == [[1, 3], [2, 7], [2, 9]]
            @test dec["cliquetree"] == [[2, 3], [3, 1]]
        end
        @testset "add_decomposition_in_instance!" begin
            add_decomposition_in_instance!(instances, "case118mod", oid)
            instance = getinstance(instances, "case118mod")
            @test instance["decompositions"][1] == oid
        end
        @testset "setfeatures_instance!" begin
            setfeatures_instance!(instances,
                                  "case118mod",
                                  Mongoc.BSON("graph" => Mongoc.BSON("testgraph1" => 1, "testgraph2" => 3),
                                              "OPF" => Mongoc.BSON("testOPF1" => "2", "testOPF2" => "4")))
            instance = getinstance(instances, "case118mod")
            @test instance["features"]["graph"]["testgraph1"] == 1
            @test instance["features"]["graph"]["testgraph2"] == 3
            @test instance["features"]["OPF"]["testOPF1"] == "2"
            @test instance["features"]["OPF"]["testOPF2"] == "4"
            setfeatures_instance!(instances,
                                  "case118mod",
                                  Mongoc.BSON("graph" => Mongoc.BSON("testgraph1" => 2, "testgraph2" => 4),
                                              "OPF" => Mongoc.BSON("testOPF1" => "1", "testOPF2" => "3")))
            instance = getinstance(instances, "case118mod")
            @test instance["features"]["graph"]["testgraph1"] == 2
            @test instance["features"]["graph"]["testgraph2"] == 4
            @test instance["features"]["OPF"]["testOPF1"] == "1"
            @test instance["features"]["OPF"]["testOPF2"] == "3"
        end
        @testset "setfeatures_decomposition!" begin
            setfeatures_decomposition!(decompositions,
                                        "case118mod",
                                        oid,
                                        Mongoc.BSON("clique" => Mongoc.BSON("testclique1" => "1", "testclique2" => "3"),
                                                    "solver" => Mongoc.BSON("solving_time" => 2, "nb_iter" => 4)))
            decomposition = getdecomposition(decompositions, "case118mod", oid)
            @test decomposition["features"]["clique"]["testclique1"] == "1"
            @test decomposition["features"]["clique"]["testclique2"] == "3"
            @test decomposition["features"]["solver"]["solving_time"] == 2
            @test decomposition["features"]["solver"]["nb_iter"] == 4

            setfeatures_decomposition!(decompositions,
                                       "case118mod",
                                       oid,
                                       Mongoc.BSON("clique" => Mongoc.BSON("testclique1" => "2", "testclique2" => "4"),
                                                   "solver" => Mongoc.BSON("solving_time" => 1, "nb_iter" => 3)))
            decomposition = getdecomposition(decompositions, "case118mod", oid)
            @test decomposition["features"]["clique"]["testclique1"] == "2"
            @test decomposition["features"]["clique"]["testclique2"] == "4"
            @test decomposition["features"]["solver"]["solving_time"] == 1
            @test decomposition["features"]["solver"]["nb_iter"] == 3
        end
        @testset "setMPIinfos" begin
            setMPIinfos(decompositions, "case1354", oids[2], Mongoc.BSON("ntasks" => 1))
            decomposition = getdecomposition(decompositions, "case1354", oids[2])
            @test decomposition["mpi"]["ntasks"] == 1
        end
        @testset "delete_decomposition" begin
            push!(decompositions, Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => "test", "oid" => "1234")))
            decomposition = Mongoc.find_one(decompositions, Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => "test", "oid" => "1234")))
            @test !isnothing(decomposition)
            delete_decomposition(decompositions, "test", "1234")
            decomposition = Mongoc.find_one(decompositions, Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => "test", "oid" => "1234")))
            @test isnothing(decomposition)
        end
        @testset "delete_similar_decomposition" begin
            instance_name1, oid1 = push_decomposition!(decompositions,
                                                       "test",
                                                        Dict("cliques" => [[1, 2], [3, 4]],
                                                             "cliquetree" => [[1, 2]],
                                                             "nb_added_edges" => 5)
                                                      )
            instance_name2, oid2 = push_decomposition!(decompositions,
                                                      "test",
                                                      Dict("cliques" => [[1, 2], [3, 4]],
                                                           "cliquetree" => [[1, 2]],
                                                           "nb_added_edges" => 5)
                                                      )
            selector = Mongoc.BSON("_id.instance_name" => "test", "cliques" => [[1, 2], [3, 4]])
            @test Mongoc.count_documents(decompositions, selector) == 2
            delete_similar_decomposition(decompositions)
            @test Mongoc.count_documents(decompositions, selector) == 1
            delete_decomposition(decompositions, instance_name1, oid1)

            instance_name1, oid1 = push_decomposition!(decompositions,
                                                       "test",
                                                       Dict("cliques" => [[1, 2], [3, 4]],
                                                            "cliquetree" => [[1, 2]],
                                                            "nb_added_edges" => 5)
                                                      )
            instance_name2, oid2 = push_decomposition!(decompositions,
                                                       "test",
                                                       Dict("cliques" => [[1, 2], [3, 4]],
                                                            "cliquetree" => [[1, 2]],
                                                            "nb_added_edges" => 5)
                                                      )
            setfeatures_decomposition!(decompositions, instance_name2, oid2, Mongoc.BSON("solver" => Mongoc.BSON("solving_time" => 5)))
            selector = Mongoc.BSON("_id.instance_name" => "test", "cliques" => [[1, 2], [3, 4]])
            @test Mongoc.count_documents(decompositions, selector) == 2
            delete_similar_decomposition(decompositions)
            @test Mongoc.count_documents(decompositions, selector) == 1
            selector = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name2, "oid" => oid2))
            @test Mongoc.count_documents(decompositions, selector) == 1
            selector = Mongoc.BSON("_id" => Mongoc.BSON("instance_name" => instance_name1, "oid" => oid1))
            @test Mongoc.count_documents(decompositions, selector) == 0
            delete_decomposition(decompositions, instance_name2, oid2)
        end
    end
    @testset "getqueries.jl" begin
        @testset "isinstance" begin
            @test isinstance(instances, "case1354")
            @test isinstance(instances, "case118")
            @test !(isinstance(instances, "case1354pegase"))
        end
        @testset "isdecomposition" begin
            @test isdecomposition(decompositions, "case1354", oids[1])
            @test isdecomposition(decompositions, "case1354", oids[2])
            @test isdecomposition(decompositions, "case118", oids[3])
            @test !isdecomposition(decompositions, "case118", "aaaaa")
            @test !isdecomposition(decompositions, "case1354", "bbbbbb")
        end
        @testset "issolved" begin
            @test issolved(decompositions, "case118mod", oid)
            @test !issolved(decompositions, "case1354", oids[1])
        end
        @testset "getOPFpath_all" begin
            paths = getOPFpath_all(instances)
            @test paths["case118"] == Dict("mat" => "path/to/OPF/mat","ctr" => "path/to/OPF/ctr")
            @test paths["case1354"] == Dict("mat" => "path/to/OPF/mat","ctr" => "path/to/OPF/ctr")
        end
        @testset "getOPFpath_one" begin
            paths = getOPFpath_one(instances, "case1354")
            @test paths == Dict("mat" => "path/to/OPF/mat","ctr" => "path/to/OPF/ctr")
            paths = getOPFpath_one(instances, "case118mod")
            @test isnothing(paths)
        end
        @testset "getcholesky" begin
            cholesky_decompositions = getcholesky(decompositions)
            @test isempty(cholesky_decompositions)
            push_decomposition!(decompositions,
                                "case118",
                                Dict("cliques" => [[1, 2], [5, 7], [6, 8]],
                                     "cliquetree" => [[1, 2], [2, 3]],
                                     "nb_added_edges" => 0,
                                    "dec_type" => "cholesky"),
                                Mongoc.BSON("clique" => Mongoc.BSON("nb_clique" => 380, "size_var" => 0.1868),
                                            "solver" => Mongoc.BSON("solving_time" => 12.5, "nb_iter" => 58),
                                            "dec_type" => "cholesky")
                               )
            cholesky_decompositions = getcholesky(decompositions)
            @test length(collect(cholesky_decompositions)) == 1
        end
        @testset "getmatpowerpath_all" begin
            paths = getmatpowerpath_all(instances)
            @test paths[1] == "path/to/Matpower"
            @test paths[2] == "path/to/Matpower"
        end
        @testset "getdecompositions_index" begin
            @test getdecompositions_index(decompositions) == [1]
            @test getdecompositions_index(decompositions; unsolved=[true, false]) == [1, 2, 3, 4]
            @test isempty(getdecompositions_index(decompositions; original=[true]))
            @test getdecompositions_index(decompositions; unsolved=[true, false], original=[true]) == [5]
        end
        @testset "getcholeskytime" begin
           @test getcholeskytime(decompositions, "case118") ==  12.5
           @test isnothing(getcholeskytime(decompositions, "case1354"))
        end
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
#
#  NLS: Non Linear Schrödinger
#

function NLSHam(u,p)

    N = 5

    @inbounds begin
    q = view(u,1:N)
    p = view(u,N+1:2*N)

    H1=zero(eltype(u))
    H2=zero(eltype(u))

    for i in 1:N
        H1+=1/4*(q[i]^2+p[i]^2)^2
    end

    for i in 2:N
        H2+=(p[i-1]^2*p[i]^2+q[i-1]^2*q[i]^2-q[i-1]^2*p[i]^2-p[i-1]^2*q[i]^2+
             4*p[i-1]*p[i]*q[i-1]*q[i])
    end

    return(H1-H2)
    end

end

function NLSODE!(du,u,p,t)

    N = 5

    q = view(u,1:N)
    p = view(u,N+1:2*N)

    du[1]=p[1]*(q[1]^2+p[1]^2)-(2*p[1]*p[2]^2-2*p[1]*q[2]^2+4*p[2]*q[1]*q[2])
    du[N+1]=-q[1]*(q[1]^2+p[1]^2)+(2*q[1]*q[2]^2-2*q[1]*p[2]^2+4*p[1]*p[2]*q[2])

    for i in 2:N-1

        du[i]=p[i]*(q[i]^2+p[i]^2)-(2*p[i-1]^2*p[i]-2*q[i-1]^2*p[i]+4*p[i-1]*q[i-1]*q[i])-(2*p[i]*p[i+1]^2-2*p[i]*q[i+1]^2+4*p[i+1]*q[i]*q[i+1])

        du[i+N]=-q[i]*(q[i]^2+p[i]^2)+(2*q[i-1]^2*q[i]-2*p[i-1]^2*q[i]+4*p[i-1]*p[i]*q[i-1])+(2*q[i]*q[i+1]^2-2*q[i]*p[i+1]^2+4*p[i]*p[i+1]*q[i+1])

   end

   du[N]=p[N]*(q[N]^2+p[N]^2)-(2*p[N-1]^2*p[N]-2*q[N-1]^2*p[N]+4*p[N-1]*q[N-1]*q[N])
   du[2*N]=-q[N]*(q[N]^2+p[N]^2)+(2*q[N-1]^2*q[N]-2*p[N-1]^2*q[N]+4*p[N-1]*p[N]*q[N-1])

end

ext: jl


text:
module timeevolution_mcwf

export mcwf, mcwf_h, mcwf_nh

using ...bases
using ...operators
using ...states
using ...ode_dopri

"""
Integrate a single Monte Carlo wave function trajectory.

Arguments
---------

dmcwf
    A function f(t, psi, dpsi) that calculates the time-derivative of psi at
    time t and stores the result in dpsi.
jumpfun
    A function f(rng, t, psi, dpsi) that uses the random number generator rng
    to determine if a jump is performed and stores the result in dpsi.
tspan
    Vector specifying the points of time for which output should be displayed.
psi0
    Initial state vector.
seed
    Seed used for the random number generator to make trajectories repeatable.

Keyword arguments
-----------------

fout (optional)
    If given this function fout(t, psi) is called every time an output should
    be displayed. To limit copying to a minimum the given state psi
    is further used and therefore must not be changed.

kwargs
    Further arguments are passed on to the ode solver.
"""
function integrate_mcwf(dmcwf::Function, jumpfun::Function, tspan, psi0::Ket, seed;
                fout=nothing,
                kwargs...)
    tmp = deepcopy(psi0)
    as_ket(x::Vector{Complex128}) = Ket(psi0.basis, x)
    as_vector(psi::Ket) = psi.data
    rng = MersenneTwister(convert(UInt, seed))
    jumpnorm = Float64[rand(rng)]
    djumpnorm(t, x::Vector{Complex128}) = norm(as_ket(x))^2 - (1-jumpnorm[1])
    function dojump(t, x::Vector{Complex128})
        jumpfun(rng, t, as_ket(x), tmp)
        for i=1:length(x)
            x[i] = tmp.data[i]
        end
        jumpnorm[1] = rand(rng)
        return ode_dopri.jump
    end
    if fout==nothing
        tout = Float64[]
        xout = Ket[]
        function fout_(t, x::Vector{Complex128})
            psi = deepcopy(as_ket(x))
            psi /= norm(psi)
            push!(tout, t)
            push!(xout, psi)
            nothing
        end
    else
        fout_(t, x::Vector{Complex128}) = fout(t, as_ket(x))
    end
    dmcwf_(t, x::Vector{Complex128}, dx::Vector{Complex128}) = dmcwf(t, as_ket(x), as_ket(dx))
    ode_event(dmcwf_, float(tspan), as_vector(psi0), fout_,
        djumpnorm, dojump;
        kwargs...)
    return fout==nothing ? (tout, xout) : nothing
end

"""
Default jump function.

Arguments
---------

rng
    Random number generator
t
    Point of time where the jump is performed.
psi
    State vector before the jump.
J
    List of jump operators.
psi_new
    Result of jump.
"""
function jump(rng, t::Float64, psi::Ket, J::Vector, psi_new::Ket)
    if length(J)==1
        operators.gemv!(complex(1.), J[1], psi, complex(0.), psi_new)
        N = norm(psi_new)
        for i=1:length(psi_new.data)
            psi_new.data[i] /= N
        end
    else
        probs = zeros(Float64, length(J))
        for i=1:length(J)
            operators.gemv!(complex(1.), J[i], psi, complex(0.), psi_new)
            #probs[i] = norm(psi_new)^2
            probs[i] = dagger(psi_new)*psi_new
        end
        cumprobs = cumsum(probs./sum(probs))
        r = rand(rng)
        i = findfirst(cumprobs.>r)
        operators.gemv!(complex(1.)/sqrt(probs[i]), J[i], psi, complex(0.), psi_new)
    end
    return nothing
end

"""
Evaluate non-hermitian Schroedinger equation.

The non-hermitian Hamiltonian is given in two parts - the hermitian part H and
the jump operators J.
"""
function dmcwf_h(psi::Ket, H::Operator,
                 J::Vector, Jdagger::Vector, dpsi::Ket, tmp::Ket)
    operators.gemv!(complex(0,-1.), H, psi, complex(0.), dpsi)
    for i=1:length(J)
        operators.gemv!(complex(1.), J[i], psi, complex(0.), tmp)
        operators.gemv!(-complex(0.5,0.), Jdagger[i], tmp, complex(1.), dpsi)
    end
    return dpsi
end

"""
Evaluate non-hermitian Schroedinger equation.

The given Hamiltonian is already the non-hermitian version.
"""
function dmcwf_nh(psi::Ket, Hnh::Operator, dpsi::Ket)
    operators.gemv!(complex(0,-1.), Hnh, psi, complex(0.), dpsi)
    return dpsi
end

"""
Integrate master equation using MCWF method with mcwf_h as derivative function.
"""
function mcwf_h(tspan, psi0::Ket, H::Operator, J::Vector;
                seed=rand(UInt), fout=nothing, Jdagger::Vector=map(dagger, J),
                tmp::Ket=deepcopy(psi0),
                display_beforeevent=false, display_afterevent=false,
                kwargs...)
    f(t, psi, dpsi) = dmcwf_h(psi, H, J, Jdagger, dpsi, tmp)
    j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new)
    return integrate_mcwf(f, j, tspan, psi0, seed; fout=fout,
                display_beforeevent=display_beforeevent,
                display_afterevent=display_afterevent,
                kwargs...)
end

"""
Integrate master equation using MCWF method with mcwf_nh as derivative function.
"""
function mcwf_nh(tspan, psi0::Ket, Hnh::Operator, J::Vector;
                seed=rand(UInt), fout=nothing,
                display_beforeevent=false, display_afterevent=false,
                kwargs...)
    f(t, psi, dpsi) = dmcwf_nh(psi, Hnh, dpsi)
    j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new)
    return integrate_mcwf(f, j, tspan, psi0, seed; fout=fout,
                display_beforeevent=display_beforeevent,
                display_afterevent=display_afterevent,
                kwargs...)
end

"""
Integrate the master equation using the MCWF method.

There are two implementations for integrating the non-hermitian
schroedinger equation:

* ``mcwf_h``: Usual formulation with Hamiltonian + jump operators separately.
* ``mcwf_nh``: Variant with non-hermitian Hamiltonian.

The ``mcwf`` function takes a normal Hamiltonian, calculates the
non-hermitian Hamiltonian and then calls mcwf_nh which is slightly faster.

Arguments
---------

tspan
    Vector specifying the points of time for which output should be displayed.
psi0
    Initial state vector.
H
    DenseOperator specifying the Hamiltonian.
J
    Vector containing all jump operators.

Keyword Arguments
-----------------

seed (optional)
    Seed used for the random number generator to make trajectories repeatable.
fout (optional)
    If given this function fout(t, psi) is called every time an output should
    be displayed. To limit copying to a minimum the given state psi
    is further used and therefore must not be changed.
Jdagger (optional)
    Vector containing the hermitian conjugates of the jump operators. If they
    are not given they are calculated automatically.
display_beforeevent [false]
    fout is called before every jump.
display_afterevent [false]
    fout is called after every jump.
kwargs
    Further arguments are passed on to the ode solver.
"""
function mcwf(tspan, psi0::Ket, H::Operator, J::Vector;
                seed=rand(UInt), fout=nothing, Jdagger::Vector=map(dagger, J),
                display_beforeevent=false, display_afterevent=false,
                kwargs...)
    Hnh = deepcopy(H)
    for i=1:length(J)
        Hnh -= 0.5im*Jdagger[i]*J[i]
    end
    f(t, psi, dpsi) = dmcwf_nh(psi, Hnh, dpsi)
    j(rng, t, psi, psi_new) = jump(rng, t, psi, J, psi_new)
    return integrate_mcwf(f, j, tspan, psi0, seed;
                fout=fout,
                display_beforeevent=display_beforeevent,
                display_afterevent=display_afterevent,
                kwargs...)
end

end #module


text:
using DataFrames
using VectorAutoregression
using Plots
pyplot()
default(show=true, reuse=false)

qdata = readtable("importable_data/qdata_processed.csv")
mdata = readtable("importable_data/mdata_processed.csv")

l = @layout([a b; c d])
plot(plot(mdata[:,2], legend=false, title="IP"),
     plot(mdata[:,3], legend=false, title="INFL"),
     plot(mdata[:,4], legend=false, title="UNEMP"),
     plot(qdata[:,2], legend=false, title="GDP"))

savefig("variable_plots.pdf")

max_stars_repo_name: lnsongxf/MF-VAR

max_issues_repo_name: lnsongxf/MF-VAR

max_forks_repo_name: lnsongxf/MF-VAR


text:
__precompile__()

"""
Provides an [`@private`](@ref) macro for hiding unexported symbols and scoped import macro
[`@local`](@ref).
"""
module PrivateModules

using Base.Meta, Compat

# Private module macro. #

export @private

"""
Make unexported symbols in a module private.

**Signature**

```julia
@private module ... end
```

**Example**

```julia
using PrivateModules

@private module M

export f

f(x) = g(x, 2x)

g(x, y) = x + y

end

using .M

f(1)      # works
M.g(1, 2) # fails
```
"""
macro private(x) private(x) end

function private(x)
    isexpr(x, :module) || error("not a valid module expression:\n\n$x")

    # Rename module.
    outer = x.args[2]
    inner = Symbol(outer, "#private")
    x.args[2] = inner

    # The inner module's `eval` should eval inside the inner module, not the outer one.
    unshift!(x.args[end].args, :(eval{T}(x::T) = Core.eval($inner, x)))

    # Build outer module.
    out = :(module $outer end)
    push!(out.args[end].args,
        x,                                   # Original module.
        :($(exports)($outer, $inner)),       # Exported symbols.
        macroexpand(:(Base.@__doc__ $outer)) # Documentation.
    )

    Expr(:toplevel, esc(out))
end

"""
**Signature**

```julia
exports(outer, inner)
```

Import all exported symbols from `inner` module into `outer` one and then re-export them.
"""
function exports(outer, inner)
    symbols = names(inner)
    imports = [Expr(:import, :., module_name(inner), s) for s in symbols]
    eval(outer, Expr(:toplevel, imports...))
    eval(outer, Expr(:toplevel, Expr(:export, symbols...)))
end

# Local imports macro. #

export @local

@eval macro $(:local)(x) localimports(x) end
"""
Local `import`, `importall` and `using` macro.

**Signature**

```julia
@local expression
```

**Examples**

```julia
function func(args...)
    @local using A, ..B, C.D
    # ...
end
```

Exported symbols from `A`, `..B`, and `C.D` are bound to local constants in `func`'s scope.

```julia
function func(args...)
    @local importall A, ..B, C.D
    # ...
end
```

All symbols from `A`, `..B`, and `C.D` are bound to local constants in `func`'s scope.

```julia
function func(args...)
    @local import A, ..B, C.D
    # ...
end
```

`A`, `B`, and `D` are bound to local constants in `func`'s scope.

**Notes**

- Macros cannot be imported using the [`@local`](@ref) macro.
- Modules listed in `@local` calls must be literals - not variables.
"""
:(@local)

function localimports(x)
    isexpr(x, :toplevel)  ? build(x)          :
    isexpr(x, :import)    ? build(x, nothing) :
    isexpr(x, :using)     ? build(x, false)   :
    isexpr(x, :importall) ? build(x, true)    :
    error("invalid expression given: '@local $x'.")
end

function build(x, flag)
    x.args[1] === :. && unshift!(x.args, :.) # Relative modules are off by one.
    build(newmodule(current_module(), x)..., flag)
end
build(x) = Expr(:block, [:(@local $arg) for arg in x.args]...)

function build(object :: Module, name, flag)
    symbols = flag ≡ nothing ? [name] : filter(Base.isidentifier, names(object, flag))
    esc(Expr(:block, [:(const $x = $object.$x) for x in symbols]...))
end
build(object, name, flag) = esc(:(const $name = $object))

function newmodule(parent, x)
    name = gensym("<local-module>")
    eval(parent, :(module $name end))
    anon = getfield(parent, name)
    eval(anon, Expr(:import, x.args...))
    getfield(anon, x.args[end]), x.args[end]
end

end # module


text:
using PWD
using Test

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


text:
const calllike = Set([:call, :foreigncall])

function extract_inner_call!(stmt, idx, once::Bool=false)
    isa(stmt, Expr) || return nothing
    (stmt.head == :toplevel || stmt.head == :thunk) && return nothing
    once |= stmt.head ∈ calllike
    for (i, a) in enumerate(stmt.args)
        isa(a, Expr) || continue
        # Make sure we don't "damage" special syntax that requires literals
        if i == 1 && stmt.head == :foreigncall
            continue
        end
        if i == 2 && stmt.head == :call && stmt.args[1] == :cglobal
            continue
        end
        ret = extract_inner_call!(a, idx, once) # doing this first extracts innermost calls
        ret !== nothing && return ret
        iscalllike = a.head ∈ calllike
        if once && iscalllike
            stmt.args[i] = NewSSAValue(idx)
            return a
        end
    end
    return nothing
end

function replace_ssa!(stmt, ssalookup)
    isa(stmt, Expr) || return nothing
    for (i, a) in enumerate(stmt.args)
        if isa(a, SSAValue)
            stmt.args[i] = SSAValue(ssalookup[a.id])
        elseif isa(a, NewSSAValue)
            stmt.args[i] = SSAValue(a.id)
        else
            replace_ssa!(a, ssalookup)
        end
    end
    return nothing
end

function renumber_ssa!(stmts::Vector{Any}, ssalookup)
    for (i, stmt) in enumerate(stmts)
        if isa(stmt, GotoNode)
            stmts[i] = GotoNode(ssalookup[stmt.label])
        elseif isa(stmt, SSAValue)
            stmts[i] = SSAValue(ssalookup[stmt.id])
        elseif isa(stmt, NewSSAValue)
            stmts[i] = SSAValue(stmt.id)
        elseif isa(stmt, Expr)
            replace_ssa!(stmt, ssalookup)
            if (stmt.head == :gotoifnot || stmt.head == :enter) && isa(stmt.args[end], Int)
                stmt.args[end] = ssalookup[stmt.args[end]]
            end
        end
    end
    return stmts
end

# Pre-frame-construction lookup
function lookup_stmt(stmts, arg)
    if isa(arg, SSAValue)
        arg = stmts[arg.id]
    end
    if isa(arg, QuoteNode)
        arg = arg.value
    end
    return arg
end

function smallest_ref(stmts, arg, idmin)
    if isa(arg, SSAValue)
        idmin = min(idmin, arg.id)
        return smallest_ref(stmts, stmts[arg.id], idmin)
    elseif isa(arg, Expr)
        for a in arg.args
            idmin = smallest_ref(stmts, a, idmin)
        end
    end
    return idmin
end

function lookup_global_refs!(ex::Expr)
    (ex.head == :isdefined || ex.head == :thunk || ex.head == :toplevel) && return nothing
    for (i, a) in enumerate(ex.args)
        ex.head == :(=) && i == 1 && continue # Don't look up globalrefs on the LHS of an assignment (issue #98)
        if isa(a, GlobalRef)
            if isdefined(a.mod, a.name) && isconst(a.mod, a.name)
                r = getfield(a.mod, a.name)
                ex.args[i] = QuoteNode(r)
            end
        elseif isa(a, Expr)
            lookup_global_refs!(a)
        end
    end
    return nothing
end

"""
    optimize!(code::CodeInfo, mod::Module)

Perform minor optimizations on the lowered AST in `code` to reduce execution time
of the interpreter.
Currently it looks up `GlobalRef`s (for which it needs `mod` to know the scope in
which this will run) and ensures that no statement includes nested `:call` expressions
(splitting them out into multiple SSA-form statements if needed).
"""
function optimize!(code::CodeInfo, mod::Module)
    code.inferred && error("optimization of inferred code not implemented")
    # TODO: because of builtins.jl, for CodeInfos like
    #   %1 = Core.apply_type
    #   %2 = (%1)(args...)
    # it would be best to *not* resolve the GlobalRef at %1

    ## Replace GlobalRefs with QuoteNodes
    for (i, stmt) in enumerate(code.code)
        if isa(stmt, GlobalRef)
            code.code[i] = QuoteNode(getfield(stmt.mod, stmt.name))
        elseif isa(stmt, Expr)
            if stmt.head == :call && stmt.args[1] == :cglobal  # cglobal requires literals
                continue
            else
                lookup_global_refs!(stmt)
            end
        end
    end

    ## Un-nest :call expressions (so that there will be only one :call per line)
    # This will allow us to re-use args-buffers rather than having to allocate new ones each time.
    old_code, old_codelocs = code.code, code.codelocs
    code.code = new_code = eltype(old_code)[]
    code.codelocs = new_codelocs = Int32[]
    ssainc = fill(1, length(old_code))
    for (i, stmt) in enumerate(old_code)
        loc = old_codelocs[i]
        inner = extract_inner_call!(stmt, length(new_code)+1)
        while inner !== nothing
            push!(new_code, inner)
            push!(new_codelocs, loc)
            ssainc[i] += 1
            inner = extract_inner_call!(stmt, length(new_code)+1)
        end
        push!(new_code, stmt)
        push!(new_codelocs, loc)
    end
    # Fix all the SSAValues and GotoNodes
    ssalookup = cumsum(ssainc)
    renumber_ssa!(new_code, ssalookup)
    code.ssavaluetypes = length(new_code)

    # Replace :llvmcall and :foreigncall with compiled variants. See
    # https://github.com/JuliaDebug/JuliaInterpreter.jl/issues/13#issuecomment-464880123
    methodtables = Vector{Union{Compiled,TypeMapEntry}}(undef, length(code.code))
    for (idx, stmt) in enumerate(code.code)
        if isexpr(stmt, :call)
            # Check for :llvmcall
            arg1 = stmt.args[1]
            if arg1 == :llvmcall || lookup_stmt(code.code, arg1) == Base.llvmcall
                uuid = uuid4()
                ustr = replace(string(uuid), '-'=>'_')
                methname = Symbol("llvmcall_", ustr)
                nargs = length(stmt.args)-4
                argnames = [Symbol("arg", string(i)) for i = 1:nargs]
                # Run a mini-interpreter to extract the types
                framecode = FrameCode(CompiledCalls, code; optimize=false)
                frame = Frame(framecode, prepare_framedata(framecode, []))
                idxstart = idx
                for i = 2:4
                    idxstart = smallest_ref(code.code, stmt.args[i], idxstart)
                end
                frame.pc = idxstart
                while true
                    pc = step_expr!(Compiled(), frame)
                    pc == idx && break
                    pc === nothing && error("this should never happen")
                end
                str, RetType, ArgType = @lookup(frame, stmt.args[2]), @lookup(frame, stmt.args[3]), @lookup(frame, stmt.args[4])
                def = quote
                    function $methname($(argnames...))
                        return Base.llvmcall($str, $RetType, $ArgType, $(argnames...))
                    end
                end
                f = Core.eval(CompiledCalls, def)
                stmt.args[1] = QuoteNode(f)
                deleteat!(stmt.args, 2:4)
                methodtables[idx] = Compiled()
            end
        end
    end

    return code, methodtables
end


text:
#needed within OrdinaryDiffEq somehow:
import Base: zero
zero(tmp::Array) = zero.(tmp)

"""
    dxyz_dt!(du,u,p::NamedTuple,tim)

Interpolate velocity from gridded fields (3D; with halos) to position `u`
(`x,y,z,fIndex`) to compute the derivative of position v time  `du_dt`.

```jldoctest
using IndividualDisplacements, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/worldwide/three_dimensional_ocean.jl"))
ref=[211. 34. -70.]
prod(isapprox.([mean(𝐼.🔴.lon) mean(𝐼.🔴.lat) mean(𝐼.🔴.z)],ref,atol=50.0))

# output

true
```
"""
function dxyz_dt!(du::Array{T,1},u::Array{T,1},𝑃::NamedTuple,tim) where T
    #compute positions in index units
    dt=(tim-𝑃.𝑇[1])/(𝑃.𝑇[2]-𝑃.𝑇[1])
    dt>1.0 ? error("dt>1.0") : nothing
    dt<0.0 ? error("dt>0.0") : nothing
    g=𝑃.u0.grid
    #
    while location_is_out(u,g)
        g.class=="PeriodicDomain" ? update_location_dpdo!(u,g) : nothing
        g.class=="CubeSphere" ? update_location_cs!(u,𝑃) : nothing
        g.class=="LatLonCap" ? update_location_llc!(u,𝑃) : nothing
    end

    x,y,z = u[1:3]
    fIndex = Int(u[4])
    nx,ny=g.fSize[fIndex]
    nz=size(𝑃.u0,2)
    #
    dx,dy,dz=[x - floor(x),y - floor(y),z - floor(z)]
    i_c,j_c = Int32.(floor.([x y])) .+ 2
    k_c = Int32.(floor.(z)) .+ 1
    #
    i_w,i_e=[i_c i_c+1]
    j_s,j_n=[j_c j_c+1]
    k_l,k_r=[k_c k_c+1]
    #
    k_c=min(k_c,nz)
    k_l=min(k_l,nz)
    k_r=min(k_r,nz)
    k_c=max(k_c,1)
    k_l=max(k_l,1)
    k_r=max(k_r,1)
    #interpolate u to position and time
    du[1]=(1.0-dx)*(1.0-dt)*𝑃.u0.f[fIndex,k_c][i_w,j_c]+
    dx*(1.0-dt)*𝑃.u0.f[fIndex,k_c][i_e,j_c]+
    (1.0-dx)*dt*𝑃.u1.f[fIndex,k_c][i_w,j_c]+
    dx*dt*𝑃.u1.f[fIndex,k_c][i_e,j_c]
    #interpolate v to position and time
    du[2]=(1.0-dy)*(1.0-dt)*𝑃.v0.f[fIndex,k_c][i_c,j_s]+
    dy*(1.0-dt)*𝑃.v0.f[fIndex,k_c][i_c,j_n]+
    (1.0-dy)*dt*𝑃.v1.f[fIndex,k_c][i_c,j_s]+
    dy*dt*𝑃.v1.f[fIndex,k_c][i_c,j_n]
    #interpolate w to position and time
    du[3]=(1.0-dz)*(1.0-dt)*𝑃.w0.f[fIndex,k_l][i_c,j_c]+
    dz*(1.0-dt)*𝑃.w0.f[fIndex,k_r][i_c,j_c]+
    (1.0-dz)*dt*𝑃.w1.f[fIndex,k_l][i_c,j_c]+
    dz*dt*𝑃.w1.f[fIndex,k_r][i_c,j_c]

    #leave face index unchanged
    du[4]=0.0
    #
    return du
end

function dxyz_dt!(du::Array{T,2},u::Array{T,2},𝑃::NamedTuple,tim) where T
    [dxyz_dt!(du[i],u[i],𝑃,tim) for i=1:size(u,2)]
end

"""
    dxy_dt!(du,u,p::NamedTuple,tim)

Interpolate velocity from gridded fields (2D; with halos) to position `u`
(`x,y,fIndex`) to compute the derivative of position v time  `du_dt`.

```jldoctest
using IndividualDisplacements, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/basics/random_flow_field.jl"))
ref=[4. 6.]
prod(isapprox.([mean(𝐼.🔴.x) mean(𝐼.🔴.y)],ref,atol=10.0))

# output

true
```

```jldoctest
using IndividualDisplacements, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/worldwide/global_ocean_circulation.jl"))
ref=[78. 88.]
prod(isapprox.([mean(𝐼.🔴.x) mean(𝐼.🔴.y)],ref,atol=10.0))

# output

true
```
"""
function dxy_dt!(du::Array{T,1},u::Array{T,1},𝑃::NamedTuple,tim) where T
    #compute positions in index units
    dt=(tim-𝑃.𝑇[1])/(𝑃.𝑇[2]-𝑃.𝑇[1])
    dt>1.0 ? error("dt>1.0") : nothing
    dt<0.0 ? error("dt>0.0") : nothing
    g=𝑃.u0.grid
    #
    while location_is_out(u,g)
        g.class=="PeriodicDomain" ? update_location_dpdo!(u,g) : nothing
        g.class=="CubeSphere" ? update_location_cs!(u,𝑃) : nothing
        g.class=="LatLonCap" ? update_location_llc!(u,𝑃) : nothing
    end

    x,y = u[1:2]
    fIndex = Int(u[3])
    nx,ny=g.fSize[fIndex]
    #
    dx,dy=[x - floor(x),y - floor(y)]
    i_c,j_c = Int32.(floor.([x y])) .+ 2
    #
    i_w,i_e=[i_c i_c+1]
    j_s,j_n=[j_c j_c+1]
    #interpolate u to position and time
    du[1]=(1.0-dx)*(1.0-dt)*𝑃.u0.f[fIndex][i_w,j_c]+
    dx*(1.0-dt)*𝑃.u0.f[fIndex][i_e,j_c]+
    (1.0-dx)*dt*𝑃.u1.f[fIndex][i_w,j_c]+
    dx*dt*𝑃.u1.f[fIndex][i_e,j_c]
    #interpolate v to position and time
    du[2]=(1.0-dy)*(1.0-dt)*𝑃.v0.f[fIndex][i_c,j_s]+
    dy*(1.0-dt)*𝑃.v0.f[fIndex][i_c,j_n]+
    (1.0-dy)*dt*𝑃.v1.f[fIndex][i_c,j_s]+
    dy*dt*𝑃.v1.f[fIndex][i_c,j_n]
    #leave face index unchanged
    du[3]=0.0
    #
    return du
end

function dxy_dt!(du::Array{T,2},u::Array{T,2},𝑃::NamedTuple,tim) where T
    [dxy_dt!(du[i],u[i],𝑃,tim) for i=1:size(u,2)]
end

"""
    dxyz_dt(du,u,𝑃::NamedTuple,tim)

Interpolate velocity from gridded fields (3D; NO halos) to position `u`
(`x,y,z`) to compute the derivative of position v time  `du_dt`.

```jldoctest
using IndividualDisplacements
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/basics/solid_body_rotation.jl"))
ref=[7.767441577479032 9.513402495574852 0.7065855989421701]
prod(isapprox.(𝐼.📌',ref,atol=1.0))

# output

true
```
"""
function dxyz_dt(du::Array{T,1},u::Array{T,1},𝑃::NamedTuple,tim) where T
    #compute positions in index units
    dt=(tim-𝑃.𝑇[1])/(𝑃.𝑇[2]-𝑃.𝑇[1])
    #
    x,y,z = u[1:3]
    nx,ny,nz = 𝑃.ioSize
    x,y=[mod(x,nx),mod(y,ny)]
    z=mod(z,nz)
    #
    dx,dy,dz=[x - floor(x),y - floor(y),z - floor(z)]
    i_c,j_c,k_c = Int32.(floor.([x y z])) .+ 1
    i_c==0 ? i_c=nx : nothing
    j_c==0 ? j_c=ny : nothing
    #k_c==0 ? k_c=nz : nothing
    #
    i_w,i_e=[i_c i_c+1]
    x>=nx-1 ? (i_w,i_e)=(nx,1) : nothing
    j_s,j_n=[j_c j_c+1]
    y>=ny-1 ? (j_s,j_n)=(ny,1) : nothing
    k_l,k_r=[k_c k_c+1]
    #z>=nz-1 ? (k_l,k_r)=(nz,1) : nothing

    #interpolate u to position and time
    du[1]=(1.0-dx)*(1.0-dt)*𝑃.u0.f[1,k_c][i_w,j_c]+
    dx*(1.0-dt)*𝑃.u0.f[1,k_c][i_e,j_c]+
    (1.0-dx)*dt*𝑃.u1.f[1,k_c][i_w,j_c]+
    dx*dt*𝑃.u1.f[1,k_c][i_e,j_c]
    #interpolate v to position and time
    du[2]=(1.0-dy)*(1.0-dt)*𝑃.v0.f[1,k_c][i_c,j_s]+
    dy*(1.0-dt)*𝑃.v0.f[1,k_c][i_c,j_n]+
    (1.0-dy)*dt*𝑃.v1.f[1,k_c][i_c,j_s]+
    dy*dt*𝑃.v1.f[1,k_c][i_c,j_n]
    #interpolate w to position and time
    du[3]=(1.0-dz)*(1.0-dt)*𝑃.w0.f[1,k_l][i_c,j_c]+
    dz*(1.0-dt)*𝑃.w0.f[1,k_r][i_c,j_c]+
    (1.0-dz)*dt*𝑃.w1.f[1,k_l][i_c,j_c]+
    dz*dt*𝑃.w1.f[1,k_r][i_c,j_c]
    #
    return du
end

function dxyz_dt(du::Array{T,2},u::Array{T,2},𝑃::NamedTuple,tim) where T
    [dxyz_dt(du[i],u[i],𝑃,tim) for i=1:size(u,2)]
end

"""
    dxy_dt(du,u,𝑃::NamedTuple,tim)

Interpolate velocity from gridded fields (2D; NO halos) to position `u`
(`x,y`) to compute the derivative of position v time  `du_dt`.

```jldoctest
using IndividualDisplacements, Statistics
p=dirname(pathof(IndividualDisplacements))
include(joinpath(p,"../examples/basics/particle_cloud.jl"))
ref=[29.381183342468674  19.890831699436823]
prod(isapprox.([mean(𝐼.🔴.x) mean(𝐼.🔴.y)],ref,atol=1.0))

# output

true
```
"""
function dxy_dt(du::Array{T,1},u::Array{T,1},𝑃::NamedTuple,tim) where T
    #compute positions in index units
    dt=(tim-𝑃.𝑇[1])/(𝑃.𝑇[2]-𝑃.𝑇[1])
    #
    x,y = u[1:2]
    nx,ny=𝑃.ioSize
    x,y=[mod(x,nx),mod(y,ny)]
    #
    dx,dy=[x - floor(x),y - floor(y)]
    i_c,j_c = Int32.(floor.([x y])) .+ 1
    i_c==0 ? i_c=nx : nothing
    j_c==0 ? j_c=ny : nothing
    #
    i_w,i_e=[i_c i_c+1]
    x>=nx-1 ? (i_w,i_e)=(nx,1) : nothing
    j_s,j_n=[j_c j_c+1]
    y>=ny-1 ? (j_s,j_n)=(ny,1) : nothing
    #interpolate u to position and time
    du[1]=(1.0-dx)*(1.0-dt)*𝑃.u0.f[1][i_w,j_c]+
    dx*(1.0-dt)*𝑃.u0.f[1][i_e,j_c]+
    (1.0-dx)*dt*𝑃.u1.f[1][i_w,j_c]+
    dx*dt*𝑃.u1.f[1][i_e,j_c]
    #interpolate v to position and time
    du[2]=(1.0-dy)*(1.0-dt)*𝑃.v0.f[1][i_c,j_s]+
    dy*(1.0-dt)*𝑃.v0.f[1][i_c,j_n]+
    (1.0-dy)*dt*𝑃.v1.f[1][i_c,j_s]+
    dy*dt*𝑃.v1.f[1][i_c,j_n]
    #
    return du
end

function dxy_dt(du::Array{T,2},u::Array{T,2},𝑃::NamedTuple,tim) where T
    [dxy_dt(du[i],u[i],𝑃,tim) for i=1:size(u,2)]
end

"""
    dxy_dt_CyclicArray(du,u,𝑃::NamedTuple,tim)

Nearest neighbor (?) velocity from gridded fields (2D; NO halos but
not needed when CyclicArrays is used to extend valid indice ranges).

_notes:_ spatial interpolation & temporal interpolation are lacking
"""
function dxy_dt_CyclicArray(du::Array{T,2},u::Array{T,2},𝑃::NamedTuple,tim) where T
    np=size(du,2)
    xi,yi=(u[1,:],u[2,:])

    #not needed:
    #xi=floor.(𝑃.xg[1,Int.(sign.(xi).*floor.(abs.(xi)))])+rem.(xi,1)
    #yi=floor.(𝑃.yg[Int.(sign.(yi).*floor.(abs.(yi))),1])+rem.(yi,1)

    i=Int.(floor.(𝑃.xg[1,Int.(floor.(xi))]))
    j=Int.(floor.(𝑃.yg[Int.(floor.(yi)),1]))
    du[1,:]=[𝑃.u[i[ii],j[ii]] for ii in 1:np]
    du[2,:]=[𝑃.v[i[ii],j[ii]] for ii in 1:np]
    return du
end

"""
    dict_to_nt(tmp::Dict)

Attempt to convert from Dict to NamedTuple
"""
dict_to_nt(tmp::Dict) = tmp=(; zip(Symbol.(keys(tmp)), values(tmp))...)

dxy_dt!(du,u,𝑃::Dict,tim) = dxy_dt!(du,u,dict_to_nt(𝑃),tim)
dxy_dt(du,u,𝑃::Dict,tim) = dxy_dt(du,u,dict_to_nt(𝑃),tim)
dxyz_dt(du,u,𝑃::Dict,tim) = dxyz_dt(du,u,dict_to_nt(𝑃),tim)

"""
    dxy_dt_replay(du,u,p::DataFrame,t)

Interpolate velocity from MITgcm float_trajectories output and return
position increment `du`.
"""
function dxy_dt_replay(du,u,p::DataFrame,t)
    tt=t/3600.0
    tt0=Int32(floor(tt))
    w=tt-tt0
    du[1]=(1.0-w)*p[tt0+1,:uVel]+w*p[tt0+2,:uVel]
    du[2]=(1.0-w)*p[tt0+1,:vVel]+w*p[tt0+2,:vVel]
end


text:
using Base: @__doc__

"""
    @cudnnDescriptor(XXX, setter=cudnnSetXXXDescriptor)

Defines a new type `cudnnXXXDescriptor` with a single field `ptr::cudnnXXXDescriptor_t` and
its constructor. The second optional argument is the function that sets the descriptor
fields and defaults to `cudnnSetXXXDescriptor`. The constructor is memoized, i.e. when
called with the same arguments it returns the same object rather than creating a new one.

The arguments of the constructor and thus the keys to the memoization cache depend on the
setter: If the setter has arguments `cudnnSetXXXDescriptor(ptr::cudnnXXXDescriptor_t,
args...)`, then the constructor has `cudnnXXXDescriptor(args...)`. The user can control
these arguments by defining a custom setter.
"""
macro cudnnDescriptor(x, set = Symbol("cudnnSet$(x)Descriptor"))
    sname = Symbol("cudnn$(x)Descriptor")
    tname = Symbol("cudnn$(x)Descriptor_t")
    cache = Symbol("cudnn$(x)DescriptorCache")
    cache_lock = Symbol("cudnn$(x)DescriptorCacheLock")
    create = Symbol("cudnnCreate$(x)Descriptor")
    destroy = Symbol("cudnnDestroy$(x)Descriptor")
    return quote
        @__doc__ mutable struct $sname                      # needs to be mutable for finalizer
            ptr::$tname
            $sname(p::$tname) = new(p)                      # prevent $sname(::Any) default constructor
        end
        Base.unsafe_convert(::Type{<:Ptr}, d::$sname)=d.ptr # needed for ccalls
        const $cache = Dict{Tuple,$sname}()                 # Dict is 3x faster than IdDict!
        const $cache_lock = ReentrantLock()
        function $sname(args...)
            lock($cache_lock) do
                get!($cache, args) do
                    ptr = $tname[C_NULL]
                    $create(ptr)
                    $set(ptr[1], args...)
                    d = $sname(ptr[1])
                    finalizer(x->$destroy(x.ptr), d)
                    return d
                end
            end
        end
    end |> esc
end

"""
    cudnnActivationDescriptor(mode::cudnnActivationMode_t,
                              reluNanOpt::cudnnNanPropagation_t,
                              coef::Cfloat)
"""
@cudnnDescriptor(Activation)

"""
    cudnnAttnDescriptor(attnMode::Cuint,
                        nHeads::Cint,
                        smScaler::Cdouble,
                        dataType::cudnnDataType_t,
                        computePrec::cudnnDataType_t,
                        mathType::cudnnMathType_t,
                        attnDropoutDesc::cudnnDropoutDescriptor_t,
                        postDropoutDesc::cudnnDropoutDescriptor_t,
                        qSize::Cint,
                        kSize::Cint,
                        vSize::Cint,
                        qProjSize::Cint,
                        kProjSize::Cint,
                        vProjSize::Cint,
                        oProjSize::Cint,
                        qoMaxSeqLength::Cint,
                        kvMaxSeqLength::Cint,
                        maxBatchSize::Cint,
                        maxBeamSize::Cint)
"""
@cudnnDescriptor(Attn)

"""
    cudnnCTCLossDescriptor(compType::cudnnDataType_t,
                           normMode::cudnnLossNormalizationMode_t,
                           gradMode::cudnnNanPropagation_t,
                           maxLabelLength::Cint)
"""
@cudnnDescriptor(CTCLoss, cudnnSetCTCLossDescriptor_v8)

"""
cudnnConvolutionDescriptor(pad::Vector{Cint},
                           stride::Vector{Cint},
                           dilation::Vector{Cint},
                           mode::cudnnConvolutionMode_t,
                           dataType::cudnnDataType_t,
                           groupCount::Cint,
                           mathType::cudnnMathType_t,
                           reorderType::cudnnReorderType_t)
"""
@cudnnDescriptor(Convolution, cudnnSetConvolutionDescriptor)

"""
    cudnnDropoutDescriptor(dropout::Real)
"""
@cudnnDescriptor(Dropout, cudnnSetDropoutDescriptorFromFloat)

"""
    cudnnFilterDescriptor(dataType::cudnnDataType_t,
                          format::cudnnTensorFormat_t,
                          nbDims::Cint,
                          filterDimA::Vector{Cint})
"""
@cudnnDescriptor(Filter, cudnnSetFilterNdDescriptor)

"""
    cudnnLRNDescriptor(lrnN::Cuint,
                       lrnAlpha::Cdouble,
                       lrnBeta::Cdouble,
                       lrnK::Cdouble)
"""
@cudnnDescriptor(LRN)

"""
    cudnnOpTensorDescriptor(opTensorOp::cudnnOpTensorOp_t,
                            opTensorCompType::cudnnDataType_t,
                            opTensorNanOpt::cudnnNanPropagation_t)
"""
@cudnnDescriptor(OpTensor)

"""
    cudnnPoolingDescriptor(mode::cudnnPoolingMode_t,
                           maxpoolingNanOpt::cudnnNanPropagation_t,
                           nbDims::Cint,
                           windowDimA::Vector{Cint},
                           paddingA::Vector{Cint},
                           strideA::Vector{Cint})
"""
@cudnnDescriptor(Pooling, cudnnSetPoolingNdDescriptor)

"""
    cudnnRNNDescriptor(algo::cudnnRNNAlgo_t,
                       cellMode::cudnnRNNMode_t,
                       biasMode::cudnnRNNBiasMode_t,
                       dirMode::cudnnDirectionMode_t,
                       inputMode::cudnnRNNInputMode_t,
                       dataType::cudnnDataType_t,
                       mathPrec::cudnnDataType_t,
                       mathType::cudnnMathType_t,
                       inputSize::Int32,
                       hiddenSize::Int32,
                       projSize::Int32,
                       numLayers::Int32,
                       dropoutDesc::cudnnDropoutDescriptor_t,
                       auxFlags::UInt32)
"""
@cudnnDescriptor(RNN, cudnnSetRNNDescriptor_v8)

"""
    cudnnRNNDataDescriptor(dataType::cudnnDataType_t,
                           layout::cudnnRNNDataLayout_t,
                           maxSeqLength::Cint,
                           batchSize::Cint,
                           vectorSize::Cint,
                           seqLengthArray::Vector{Cint},
                           paddingFill::Ptr{Cvoid})
"""
@cudnnDescriptor(RNNData)

"""
    cudnnReduceTensorDescriptor(reduceTensorOp::cudnnReduceTensorOp_t,
                                reduceTensorCompType::cudnnDataType_t,
                                reduceTensorNanOpt::cudnnNanPropagation_t,
                                reduceTensorIndices::cudnnReduceTensorIndices_t,
                                reduceTensorIndicesType::cudnnIndicesType_t)
"""
@cudnnDescriptor(ReduceTensor)

"""
    cudnnSeqDataDescriptor(dataType::cudnnDataType_t,
                           nbDims::Cint,
                           dimA::Vector{Cint},
                           axes::Vector{cudnnSeqDataAxis_t},
                           seqLengthArraySize::Csize_t,
                           seqLengthArray::Vector{Cint},
                           paddingFill::Ptr{Cvoid})
"""
@cudnnDescriptor(SeqData)

"""
    cudnnSpatialTransformerDescriptor(samplerType::cudnnSamplerType_t,
                                      dataType::cudnnDataType_t,
                                      nbDims::Cint,
                                      dimA::Vector{Cint})
"""
@cudnnDescriptor(SpatialTransformer, cudnnSetSpatialTransformerNdDescriptor)

"""
    cudnnTensorDescriptor(format::cudnnTensorFormat_t,
                          dataType::cudnnDataType_t,
                          nbDims::Cint,
                          dimA::Vector{Cint})
"""
@cudnnDescriptor(Tensor, cudnnSetTensorNdDescriptorEx)

"""
    cudnnTensorTransformDescriptor(nbDims::UInt32,
                                   destFormat::cudnnTensorFormat_t,
                                   padBeforeA::Vector{Int32},
                                   padAfterA::Vector{Int32},
                                   foldA::Vector{UInt32},
                                   direction::cudnnFoldingDirection_t)
"""
@cudnnDescriptor(TensorTransform)


text:
# # Temperature dependencies

## load packages
using Photosynthesis
using PlotPlants
FT = Float32;
#------------------------------------------------------------------------------

# ## Simple example
## define photosynthesis system and leaf (C3 and C4), and envir
c3_set = C3CLM(FT);
c4_set = C4CLM(FT);
leaf_3 = Leaf{FT}(APAR=1000);
leaf_4 = Leaf{FT}(APAR=1000);
envir  = AirLayer{FT}();

## define leaf temperature, total leaf conductance to CO₂, and internal CO₂
T   = FT(300);
glc = FT(0.1);
p_i = rand(FT) + 20;

## remember to update the temperature dependencies when temperature changes
println("initialize temperature dependencies");
leaf_temperature_dependence!(c3_set, leaf_3, envir, T);
leaf_temperature_dependence!(c4_set, leaf_4, envir, T);

println("calculate photosynthesis from known internal CO₂ partial pressure");
leaf_photosynthesis!(c3_set, leaf_3, envir, PCO₂Mode(), p_i);
leaf_photosynthesis!(c4_set, leaf_4, envir, PCO₂Mode(), p_i);
@show leaf_3.An;
@show leaf_4.An;

println("calculate photosynthesis from known leaf conductance to CO₂");
leaf_photosynthesis!(c3_set, leaf_3, envir, GCO₂Mode(), glc);
leaf_photosynthesis!(c4_set, leaf_4, envir, GCO₂Mode(), glc);
@show leaf_3.An;
@show leaf_4.An;
#------------------------------------------------------------------------------

# ## A-Ci curve
# Here we show an example of the A-Ci curves for C3 and C4 photosynthesis. As
#     stomatal conductance may differ when other environmental conditions
#     change, we leave the examples of photosynthesis responses to the
#     environment for `StomataModels` package.
## temperature not changing, no temperature dependencies update required
_p3 = collect(FT, 5:1:200);
_p4 = collect(FT, 0:0.1:15.01);
_a3 = similar(_p3);
_a4 = similar(_p4);
for i in eachindex(_p3)
    leaf_photosynthesis!(c3_set, leaf_3, envir, PCO₂Mode(), _p3[i]);
    _a3[i] = leaf_3.An;
end
for i in eachindex(_p4)
    leaf_photosynthesis!(c4_set, leaf_4, envir, PCO₂Mode(), _p4[i]);
    _a4[i] = leaf_4.An;
end

_fig,_axes = create_canvas("A-Ci curve"; ncol=2);
_ax1,_ax2 = _axes;
_ax1.plot(_p3, _a3, "k-", label="C3");
_ax2.plot(_p4, _a4, "k-", label="C4");
_ax1.set_xlabel("Leaf internal CO₂ (Pa)");
_ax1.set_ylabel("Anet (μmol m⁻² s⁻¹)");
_ax2.set_xlabel("Leaf internal CO₂ (Pa)");
_ax1.legend(loc="lower right");
_ax2.legend(loc="lower right");
_fig
#------------------------------------------------------------------------------


text:
const DEFAULT_TRAIT_SUPERTYPE = Any

# -----------------------------------------------------------------------------

"""
    @trait <name> [as <category>] [prefix <positive>,<negative>] [with <trait1,trait2,...>]

Create a new trait type for `name` called `\$(name)Trait`:
* If the `as` clause is provided, then `category` (an abstract type) will be used as the
  super type of the trait type.
* If the `prefix` clause is provided, then it allows the user to choose different prefixes
  than the default ones (`Can` and `Cannot`) e.g. `prefix Is,Not` or `prefix Has,Not`.
* If the `with` clause is provided, then it defines a composite trait from existing traits.
Note that you must specify at least 2 traits to make a composite trait.
"""
macro trait(name::Symbol, args...)
    category, prefixes, underlying_traits = parse_trait_args(args)
    pos, neg = prefixes.args
    mod = __module__

    trait_type = trait_type_name(name)
    this_can_type = Symbol("$(pos)$(name)")
    this_cannot_type = Symbol("$(neg)$(name)")
    lower_name = lowercase(String(name))

    # The default is "cannot".  But if it's a composite trait, then the
    # default is "can" only when all of it's underlying traits are also "can".
    default_trait_function = trait_func_name(name)
    default_expr = if underlying_traits !== nothing
        # Construct something like: flytrait(x) === CanFly() && swimtrait(x) === CanSwim()
        traits_func_names = [trait_func_name(mod, sym) for sym in underlying_traits.args]
        traits_can_types  = underlying_traits.args
        condition =
            Expr(:(&&),
                [Expr(:call, :(===), Expr(:call, f, :x), Expr(:call, g))
                        for (f,g) in zip(traits_func_names, traits_can_types)]...)

        # Construct expression like: [condition] ? CanFlySwim() : CannotFlySwim()
        Expr(:if, condition, Expr(:call, this_can_type), Expr(:call, this_cannot_type))
    else
        # The default is "cannot" for every type
        Expr(:call, this_cannot_type)
    end

    # If it's composite trait, then I want to maintain a mapping from the
    # can-type to the underlying's can-types.  It is needed for interface checks.
    composite_expr = if underlying_traits !== nothing
        traits_can_types = underlying_traits.args
        quote
        BinaryTraits.push_composite_map!($mod, $this_can_type, Set([$(traits_can_types...)]))
        end
    else
        nothing
    end

    prefixes = (pos, neg)
    name_node = QuoteNode(name)

    expr = quote
        abstract type $trait_type <: $category end
        struct $this_can_type <: $trait_type end
        struct $this_cannot_type <: $trait_type end
        $(default_trait_function)(x::Any) = $default_expr

        BinaryTraits.istrait(::Type{$trait_type}) = true

        $composite_expr
    end
    display_expanded_code(expr)
    return esc(expr)
end

"""
Parse arguments for the @trait macro.
"""
function parse_trait_args(args)

    category = Symbol(DEFAULT_TRAIT_SUPERTYPE)
    prefixes = Expr(:tuple, :Can, :Cannot)
    traits = nothing

    usage = "Invalid @trait usage. See doc string for details."
    length(args) % 2 === 0 || throw(SyntaxError(usage))

    if length(args) > 0 && args[1] == :as
        category = args[2]
        category isa Symbol || throw(SyntaxError("Not a symbol: $category. $usage"))
        args = args[3:end]
    end

    if length(args) > 0 && args[1] == :prefix
        prefixes =  args[2]
        is_tuple_of_symbols(prefixes; n = 2) || throw(SyntaxError(usage))
        args = args[3:end]
    end

    if length(args) > 0 && args[1] == :with
        traits =  args[2]
        (traits isa Symbol || is_tuple_of_symbols(traits; n = 2, op = >=)) ||
            throw(SyntaxError(usage))
        args = args[3:end]
    end

    return (category, prefixes, traits)
end


text:
function rate_plot(sepp::SEPP; step = nothing, from_idx = nothing, to_idx = nothing)

    ts = sepp.data
    isnothing(ts) && error("No data in the model, can't plot")

    times = ts.times
    isnothing(from_idx) && (from_idx = 1)
    isnothing(to_idx) && (to_idx = length(times))
    starttime = times[from_idx]
    endtime = times[to_idx]

    ((starttime isa TimeType) && !(step isa Union{TimeType, Nothing})) && (@warn "If times are TimeType step must be TimeType, ignoring step value."; step=nothing)

    if isnothing(step)
        if starttime isa TimeType
            anytimes = DateTime(starttime):Dates.Hour(1):DateTime(endtime)
        else
            anytimes = (starttime:oneunit(starttime-endtime):endtime)
        end
    else
        anytimes = (starttime:step:endtime)
    end

    μ = sepp.μ
    ϕ = sepp.ϕ
    γ = sepp.γ
    δ = sepp isa SEMPPExpKern ? sepp.δ : 0

    lamb = μ .+ ϕ .* volfunc(anytimes, ts, γ, δ)

    points = fill(0, (to_idx-from_idx + 1,))

    rate_layer = layer(x = anytimes, y = lamb, color = [color("black")], Geom.line)

    points_layer = layer(x = times[from_idx:to_idx], y = points, alpha = [0.7], Geom.point)

    plt = plot(rate_layer, points_layer, Guide.xlabel("time"), Guide.ylabel("intensity"), Theme(grid_line_width=0mm, highlight_width = 0mm))

    return plt
end

function marks_plot(mts::MarkedTimeSeries; from_idx::Union{Nothing, Integer} = nothing, to_idx::Union{Nothing, Integer} = nothing)
    marks = mts.marks
    times = mts.times

    isnothing(from_idx) && (from_idx = 1)
    isnothing(to_idx) && (to_idx = length(times))

    plt = plot(x=times[from_idx:to_idx], y=marks[from_idx:to_idx], Geom.hair, Geom.point, Guide.xlabel("time"), Guide.ylabel("marks"))

    return plt
end

function marks_plot(sempp::SEMPPExpKern)
    mts = sempp.data
    isnothing(mts) && error("No data in the model, can't plot")

    return marks_plot(mts)
end

function marked_rate_plot(sempp::SEMPPExpKern; step = nothing, from_idx = nothing, to_idx = nothing)

    mts = sempp.data
    isnothing(mts) && error("No data in the model, can't plot")

    rate_plt = rate_plot(sempp, step = step, from_idx = from_idx, to_idx = to_idx)

    marks_plt = marks_plot(mts, from_idx = from_idx, to_idx = to_idx)

    plt = vstack(rate_plt, marks_plt)

    return plt
end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb02nd(datfile, io=stdout)
    ZERO = 0.0e0
    NIN = 5
    NOUT = 6
    MMAX = 20
    NMAX = 20
    LMAX = 20
    LDC = max( MMAX, NMAX+LMAX )
    LDX = NMAX
    LENGQ = 2*min(MMAX,NMAX+LMAX)-1
    LIWORK = NMAX+2*LMAX
    LDWORK = max(2, max( MMAX, NMAX+LMAX ) + 2*min( MMAX, NMAX+LMAX ), min( MMAX, NMAX+LMAX ) + max( ( NMAX+LMAX )*( NMAX+LMAX-1 )÷2, MMAX*( NMAX+LMAX-( MMAX-1 )÷2 ) ) + max( 6*(NMAX+LMAX)-5, LMAX*LMAX + max( NMAX+LMAX, 3*LMAX ) ) )
    LBWORK = NMAX+LMAX
    C = Array{Float64,2}(undef, LDC,NMAX+LMAX)
    X = Array{Float64,2}(undef, LDX,LMAX)
    Q = Array{Float64,1}(undef, LENGQ)
    INUL = Array{BlasBool,1}(undef, NMAX+LMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    IWORK = Array{BlasInt,1}(undef, LIWORK)
    BWORK = Array{BlasBool,1}(undef, LBWORK)
f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    M = parse(BlasInt, vs[1])
    N = parse(BlasInt, vs[2])
    L = parse(BlasInt, vs[3])
    RANK = parse(BlasInt, vs[4])
    THETA = parse(Float64, replace(vs[5],'D'=>'E'))
    TOL = parse(Float64, replace(vs[6],'D'=>'E'))
    RELTOL = parse(Float64, replace(vs[7],'D'=>'E'))
if ( M<0 || M>MMAX )
elseif ( N<0 || N>NMAX )
elseif ( L<0 || L>LMAX )
elseif ( RANK>min( MMAX, NMAX ) )
elseif ( RANK<0 && THETA<ZERO )
else

    vs = String[]
    _isz,_jsz = (M,N+L)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       C[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    RANK1 = RANK
    THETA1 = THETA
# interp call 1

    RANK, THETA, INFO, IWARN = SLICOT.mb02nd!(M, N, L, RANK, THETA, C, X, Q, INUL, TOL, RELTOL)
    println(io, "RANK = $RANK")
    println(io, "THETA = $THETA")
    @test INFO == 0
    INFO == 0 || return
    println(io, "IWARN = $IWARN")

if ( INFO!=0 )
else
if ( IWARN!=0 )
else
end # if
    MINMNL = min( M, N+L )
    LOOP = MINMNL - 1
    println(io, "Q:")
    show(io, "text/plain", Bidiagonal(Q[1:MINMNL],Q[MINMNL+1:MINMNL+LOOP],'U'))
    println(io)
    println(io, "X:")
    show(io, "text/plain", X[1:N,1:L])
    println(io)

# interp output 1
    println(io, "C:")
    _nc = N+L
    _nr = max( M, N + L )
    show(io, "text/plain", C[1:_nr,1:_nc])
    println(io)

    inul = convert(Vector{Bool}, INUL[1:N+L])
    println(io, "INUL:")
    show(io, "text/plain", inul[1:N+L])
    println(io)
end # if
end # if
    close(f)
end # run_mb02nd()


text:
using MinFEM
using LinearAlgebra

function solve_gradient(file_name::String)

  mesh = import_mesh(file_name)

  G = asmGradient(mesh, qdim=2)

  f(x) = [x[1]+x[2]; x[1]-x[2]]
  s = evaluateMeshFunction(mesh, f, qdim=2)

  y = G*s;

  return y

end

@test norm(abs.(solve_gradient("test_v4.msh")).-1.0, Inf) < 1e-10


text:
module TestOnce

using ConcurrentUtils
using ConcurrentUtils.Internal: NOTSET
using Test

const O1 = Once{Vector{Int}}()
const O2 = Once{Vector{Int}}(() -> zeros(Int, 3))
const O3 = Once(() -> zeros(Int, 3))

function test_once_identity()
    @test O1[] === O1[]
    @test O2[] === O2[]
    @test O3[] === O3[]
end

const O1_NOUSE = Once{Vector{Int}}()
const O2_NOUSE = Once{Vector{Int}}(() -> zeros(Int, 3))
const O3_NOUSE = Once(() -> zeros(Int, 3))

function test_notset()
    @test O1_NOUSE.promise.value === NOTSET
    @test O2_NOUSE.promise.value === NOTSET
    @test O3_NOUSE.promise.value === NOTSET
end

end  # module


text:
function unwrap!(A::AbstractVector,
                 wrap_around::NTuple{1, Bool}=(false,),
                 seed::Int=-1)
    @inbounds previous_element = A[1]
    difference = zero(previous_element)
    periods = 0
    @inbounds for i in 2:length(A)
        difference = A[i] - previous_element
        if difference > π
            periods -= 1
        elseif difference < -π
            periods += 1
        end
        previous_element = A[i]
        A[i] = previous_element + 2 * one(previous_element) * π * periods
    end
    return A
end

max_forks_repo_name: platawiec/Unwrap


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# # Attempt to use ConcreteStructs.jl with Parameters.jl
#
# * Gen Kuroki
# * 2021-09-06
# * https://github.com/jonniedie/ConcreteStructs.jl/issues/4#issuecomment-913998030
#
# I have tried to make ConcreteStructs.jl and Parameters.jl work together well.
#
# __Conclusion:__ It is possible to do so by making the following two changes.
#
# * Change `@concrete` not creating the inner constructor, so that the `Foo{__T_a, __T_b, __T_c}(a, b, c)`-type default constructor will be defined.
# * Change `@with_kw` expanding macros in the argument, [like `Base.@kwdef`](https://github.com/JuliaLang/julia/blob/4931faa34a8a1c98b39fb52ed4eb277729120128/base/util.jl#L455).
#
# Then `@concrete` works well with `@with_kw` and more completely with `Base.@kwdef`.
#
# See below for details.

# %%
VERSION

# %%
using ConcreteStructs
using Parameters

macro macroexpand_rmln(code)
    :(macroexpand($__module__, $(QuoteNode(code)), recursive=true) |>
        Base.remove_linenums!)
end

# %% [markdown]
# ## Plain struct

# %%
struct Foo{A, B, C}
    a::A
    b::B
    c::C
end

# %%
methods(Foo)

# %%
methods(Foo{1,2,3})

# %% [markdown]
# The default constructor `Foo{A, B, C}(a, b, c)` is defined.

# %% [markdown]
# ## @concrete struct
#
# `@concrete` removes the `Foo_concrete{__T_a, __T_b, __T_c}(a, b, c)`-type default constructor.

# %%
@concrete struct Foo_concrete
    a
    b
    c
end

# %%
methods(Foo_concrete)

# %%
methods(Foo_concrete{1,2,3})

# %% [markdown]
# `Foo_concrete{__T_a, __T_b, __T_c}(a, b, c)` is not defined because only the inner constructor `Foo_concrete(a::__T_a, b::__T_b, c::__T_c)` is defined.

# %%
@macroexpand_rmln @concrete struct Foo_concrete
    a
    b
    c
end

# %% [markdown]
# ## Base.@kwdef concrete struct

# %%
Base.@kwdef @concrete struct Foo_kwdef_concrete
    a = 1
    b = 2.0
    c = "three"
end

# %%
Foo_kwdef_concrete()

# %%
Foo_kwdef_concrete{Int, Float64, String}(a = 4, b = 5.0, c = "six")

# %% [markdown]
# The reason for this error is that `Foo_kwdef_concrete{__T_a, __T_b, __T_c}(a, b, c)` is not defined.

# %%
methods(Foo_kwdef_concrete)

# %%
methods(Foo_kwdef_concrete{1,2,3})

# %%
@macroexpand_rmln Base.@kwdef @concrete struct Foo_kwdef_concrete
    a
    b
    c
end

# %% [markdown]
# ## @with_kw @concrete struct causes error

# %%
@with_kw @concrete struct Foo_with_kw_concrete
    a = 1
    b = 2.0
    c = "three"
end

# %% [markdown]
# [The first line of `Base.@kwdef`](https://github.com/JuliaLang/julia/blob/4931faa34a8a1c98b39fb52ed4eb277729120128/base/util.jl#L455) expands macros in the argument expression:
#
# ```julia
#     expr = macroexpand(__module__, expr) # to expand @static
# ```
#
# This is what makes `Base.@kwdef @concrete struct` possible.  Make a similar change to `Parameters.@with_kw`.

# %% [markdown]
# ## Change @with_kw expanding macros in the argument

# %%
@eval Parameters macro with_kw(typedef)
    typedef = macroexpand(__module__, typedef) # inserted
    return esc(with_kw(typedef, __module__, true))
end

# %%
@with_kw @concrete struct Foo_with_kw_concrete
    a = 1
    b = 2.0
    c = "three"
end

# %% [markdown]
# Okay, it seems to have worked.  But...

# %%
methods(Foo_with_kw_concrete)

# %%
methods(Foo_with_kw_concrete{1,2,3})

# %%
Foo_with_kw_concrete()

# %% [markdown]
# The reason for this error is that `Foo_with_kw_concrete{__T_a, __T_b, __T_c}(a, b, c)` is not defined again.

# %%
@macroexpand_rmln @with_kw @concrete struct Foo_with_kw_concrete
    a = 1
    b = 2.0
    c = "three"
end

# %% [markdown]
# ## Change `@concrete` not creating the inner constructor
#
# Change `ConcreteStructs._concretize(expr)` not creating the inner constructor.

# %%
@eval ConcreteStructs function _concretize(expr)
    expr isa Expr && expr.head == :struct || error("Invalid usage of @concrete")

    is_mutable = expr.args[1]
    struct_name, type_params, super = _parse_head(expr.args[2])
    line_tuples = _parse_line.(expr.args[3].args)
    lines = first.(line_tuples)
    type_params_full = (type_params..., filter(x -> x!==nothing, last.(line_tuples))...)

    struct_type = if length(type_params_full) == 0
        struct_name
    else
        Expr(:curly, struct_name, type_params_full...)
    end

    head = Expr(:(<:), struct_type, super)
    # constructor_expr = _make_constructor(struct_name, type_params, type_params_full, lines)
    # body = Expr(:block, lines..., constructor_expr)
    body = Expr(:block, lines...)
    struct_expr = Expr(:struct, is_mutable, head, body)

    return struct_expr, struct_name, type_params
end

# %%
@macroexpand_rmln @concrete struct Bar_concrete
    a = 1
    b = 2.0
    c = "three"
end

# %% [markdown]
# The inner constructor has been deleted.

# %% [markdown]
# ## @concrete works well with @with_kw

# %%
@with_kw @concrete struct Bar_with_kw_concrete
    a = 1
    b = 2.0
    c = "three"
end

# %%
methods(Bar_with_kw_concrete)

# %%
methods(Bar_with_kw_concrete{1,2,3})

# %% [markdown]
# The default constructor `Bar_with_kw_concrete{__T_a, __T_b, __T_c}(a, b, c)` is implicitly defined.

# %%
Bar_with_kw_concrete()

# %%
Bar_with_kw_concrete(c = '3')

# %%
Bar_with_kw_concrete{Int, Float64, String}(a = 4, b = 5.0, c = "six")

# %%
Bar_with_kw_concrete(4, 5.0, "six")

# %% [markdown]
# `ConcreteStructs.@concrete` works well with `Parameters.@with_kw`!

# %%
@macroexpand_rmln @with_kw @concrete struct Bar_with_kw_concrete
    a = 1
    b = 2.0
    c = "three"
end

# %% [markdown]
# ## @concrete works well more completely with Base.@kwdef

# %%
Base.@kwdef @concrete struct Bar_kwdef_concrete
    a = 1
    b = 2.0
    c = "three"
end

# %%
Bar_kwdef_concrete()

# %%
Bar_kwdef_concrete(c = '3')

# %%
Bar_kwdef_concrete{Int, Float64, String}(a = 4, b = 5.0, c = "six") # not error

# %%
Bar_kwdef_concrete(4, 5.0, "six")

# %% [markdown]
# `ConcreteStructs.@concrete` works more completely well with `Base.@kwdef`!

# %%
@macroexpand_rmln Base.@kwdef @concrete struct Bar_kwdef_concrete
    a = 1
    b = 2.0
    c = "three"
end

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
mutable struct Cvor
    color::Char
    d::Int
    pred::Int
    ime::String
    prezime::String
    pol::Char
end
mutable struct graph
   AdjMatrix::Array{Int, 2}
   V::Array{Cvor, 1}
end

function printGraph(G)
  n = size(G.AdjMatrix, 1)
  for i in 1:n
    for j in 1:n
      print(" $(G.AdjMatrix[i,j])")
    end
  println()
  end
  for i in 1:n
     println("$(G.V[i].color) $(G.V[i].d) $(G.V[i].pred) $(G.V[i].ime) $(G.V[i].prezime) $(G.V[i].pol) ")
  end
end

function emptyGraph(n)
   AdjMatrix = zeros(n, n)
   V = Array{Cvor}(undef, n)
   for i in 1:n
     V[i] = Cvor('W', 0, 0, "", "", ' ')
   end
   return graph(AdjMatrix, V)
end

function nonEmptyGraph(AdjMat)
   n = size(AdjMat, 1)
   V = Array{Cvor}(undef, n)
   for i in 1:n
     V[i] = Cvor('W', 0, 0, "", "", ' ')
   end
   return graph(copy(AdjMat), V)
end


text:
# TODO currently the program is dynamically allocated

_mode(instr::Symbol) = (:Implicit, instr, nothing)

function _mode(instr::Expr, line::Integer=-1)
    if @capture(instr, (op_,))
        return (:Implicit, op, nothing)
    elseif @capture(instr, (op_, arg_UInt8))
        return (:Immediate, op, arg)
    elseif @capture(instr, (op_, Π(arg_UInt8)))
        return (:ZeroPage, op, arg)
    elseif @capture(instr, (op_, Π(arg_UInt16)))
        return (:Absolute, op, arg)
    elseif @capture(instr, (op_, Π(arg_UInt8), X))
        return (:ZeroPageX, op, arg)
    elseif @capture(instr, (op_, Π(arg_UInt8), Y))
        return (:ZeroPageY, op, arg)
    elseif @capture(instr, (op_, Π(arg_UInt16), X))
        return (:AbsoluteX, op, arg)
    elseif @capture(instr, (op_, Π(arg_UInt16), Y))
        return (:AbsoluteY, op, arg)
    elseif @capture(instr, (op_, (Π(arg_UInt16),)))
        return (:Indirect, op, arg)
    elseif @capture(instr, (op_, (Π(arg_UInt8),X)))
        return (:IndirectX, op, arg)
    elseif @capture(instr, (op_, (Π(arg_UInt8),Y)))
        return (:IndirectY, op, arg)
    end
    throw(ArgumentError("Assembly Error: Unparsable instruction at line $line: `$instr`."))
end

_assembly_lookup(op::String, mode::Symbol) = ASSEMBLY_DICT[(op,mode)]

_assembly_args(::Nothing) = nothing
_assembly_args(x::UInt8) = x
_assembly_args(x::UInt16) = reinterpret(UInt8, [x])

function assemble!(program::Vector{UInt8}, instr::Union{Expr,Symbol})
    mode, op, arg = _mode(instr)
    op = lowercase(string(op))
    opcode = _assembly_lookup(op, mode)
    args = _assembly_args(arg)
    if args == nothing
        push!(program, opcode)
    else
        append!(program, UInt8[opcode; args])
    end
end

macro assemble(block::Expr)
    program = Vector{UInt8}(undef)  # for now we'll dynamically allocate
    if @capture(block, begin instrs__ end)
        for instr ∈ instrs
            assemble!(program, instr)
        end
    else
        throw(ArgumentError("Assembly Error: Malformed instruction block."))
    end
    :($program)
end


text:
convert_to_int_vector(line) = map(x -> parse(Int64, x), split(line, ","))

function read_data(lines)
    function read_interval(s)
        a, b = map(x -> parse(Int64, x), split(s, "-"))
        (a=a, b=b)
    end

    ranges = Dict{String, Set{NamedTuple{(:a, :b), Tuple{Int64, Int64}}}}()
    ix1, ix2 = findall(isempty, lines)
    for line in lines[1:ix1-1]
        label, rest = split(line, ": ")
        intervals = split(rest, " or ")
        ranges[label] = Set(map(read_interval, intervals))
    end

    tickets = Vector{Int64}[]
    for line in lines[ix2+2:end]
        push!(tickets, convert_to_int_vector(line))
    end

    (ranges, tickets)
end

function is_valid(x, ranges)
    item_in_intervals(item, intervals) = any(a <= item <= b for (a, b) in intervals)

    any(item_in_intervals(x, intervals) for intervals in values(ranges))
end

function invalid_nearby_tickets(tickets, ranges)
    invalid_values = Int64[]
    for ticket in tickets
        for x in ticket
            if !is_valid(x, ranges)
                push!(invalid_values, x)
            end
        end
    end
    sum(invalid_values)
end

function main()
    lines = readlines("../resources/input_16.txt")
    ranges, nearby_tickets = read_data(lines)
    result = invalid_nearby_tickets(nearby_tickets, ranges)
    println(result)
end

main()


text:
module Example

    greet() = println("Hi")

end # module


text:
module ANOVA
import GLM: ccdf, cholesky!, CoefTable, deviance, dof_residual, FDist, fit,
            ftest, LinearModel, nobs, drop_term, StatsModels.hasintercept, term,
            width, InteractionTerm

import Base: show, summary

"""
    Anova(mod, anovatype = 3)

Compute an ANOVA table for the given linear model.
"""
struct Anova
  Source::Vector{String}
  DF::Vector{Float64}
  SS::Vector{Float64}
  MSS::Vector{Float64}
  F::Vector{Float64}
  p::Vector{Float64}
  function Anova(mod; anovatype = 3)
    mf = mod.mf
    mm = mod.mm
    mod = mod.model

    eff = effects(mod) # get effects
    response = mf.f.lhs
    terms = drop_term(mf.f.rhs, term(1))

    ## calculate variables for residuals (this is the full modell) ##
    DFres = dof_residual(mod)           # degrees of freedom
    SSres = deviance(mod)               # calculate sum of squares of weighted residuals
    MSSres = SSres / DFres              # calculate mean sum of squares for residuals

    # check if there is an intercept
    v = hasintercept(mf.f)

    # get assigned numbers for factors
    factors_assig = unique(mm.assign)

    ## calculate some variables for computation ##
    n = nobs(mod)                                # sample size
    Nfactors = length(factors_assig) - v         # number of factors

    ## create arrays for ANOVA factors ##
    Source = fill("", Nfactors + 1)
    DF = zeros(Nfactors + 1)
    SS = zeros(Nfactors + 1)
    MSS = zeros(Nfactors + 1)
    RSS = zeros(Nfactors + 1)
    F = zeros(Nfactors + 1)
    p = zeros(Nfactors + 1)

    for i in 1:Nfactors
      Source[i] = string(terms.terms[i])
      mask = mm.assign .== factors_assig[i+v]
      if anovatype == 3 # if this is a type III ANOVA, we use a similar procedure to R's drop1() function
        fullmod = mod
        newmod = droptermbymask(fullmod, mask)
        RSS[i] = deviance(newmod)
        SS[i] = RSS[i] - deviance(fullmod)
        DF[i] = count(mask)
        MSS[i] = SS[i]/DF[i]
        F[i] = ftest(fullmod,newmod).fstat[1]
      elseif anovatype == 2 # if this is a type II ANOVA
        full_model_mask = falses(length(mm.assign))
        for j in 1:Nfactors
          term = terms.terms[j]
          if typeof(term) <: InteractionTerm
              if (!(typeof(terms.terms[i]) <: InteractionTerm) && # this is not an interaction term
                  terms.terms[i] ∈ term.terms) ||
                  (typeof(terms.terms[i]) <: InteractionTerm && # this is an interaction term
                  all(x -> x ∈ term.terms, terms.terms[i].terms) && term ≠ terms.terms[i])
                    full_model_mask = merge_bool_array(full_model_mask, mm.assign .== factors_assig[j+v])
              end
          end
        end
        fullmod = droptermbymask(mod, full_model_mask) # we have to drop all interaction terms containing the current term in the full modell
        newmod = droptermbymask(mod, merge_bool_array(full_model_mask, mask))
        RSS[i] = deviance(newmod)
        SS[i] = RSS[i] - deviance(fullmod)
        DF[i] = count(mask)
        MSS[i] = SS[i] / DF[i]
        F[i] = MSS[i] / MSSres
      else # type I ANOVA
        SS[i] = sum(abs2, eff[mask])
        RSS[i] = 0
        DF[i] = count(mask)
        MSS[i] = SS[i] / DF[i]
        F[i] = MSS[i] / MSSres
      end
      p[i] = ccdf(FDist(DF[i], DFres), F[i])
    end
    Source[end] = "Residuals"
    DF[end] = DFres
    SS[end] = SSres
    MSS[end] = MSSres
    F[end] = 0.
    p[end] = 0.
    new(Source, DF, SS, MSS, F, p)
  end
  Anova(Source::AbstractVector{<:AbstractString},
        DF::AbstractVector{<:Real},
        SS::AbstractVector{<:Real},
        F::AbstractVector{<:Real},
        p::AbstractVector{<:Real}) = new(Source, DF, SS, SS ./ DF, F, p)
end

## helper functions ##

merge_bool_array(a, b) = ifelse.(b, b, a)

# calculate effects for type I ANOVA
effects(mod) = (mod.pp.X / cholesky!(mod.pp).U)' * mod.rr.y

# this will drop variables from the matrix using the mask argument then fit a new linear model
droptermbymask(mod, mask) = fit(LinearModel, mod.pp.X[:,.!mask], mod.rr.y)

summary(io::IO, obj::Anova) =
  println(io, "ANOVA table for linear model")

show(io::IO, obj::Anova) =
  println(io,
          summary(obj),
          CoefTable(hcat(obj.DF, obj.SS, obj.MSS, obj.F, obj.p),
                    ["DF", "SS", "MSS", "F", "p"],
                    obj.Source,
                    5))

export Anova

end # module


text:
struct PSummary
    results
    table
end

function Base.display(x::PSummary)
     println("");println(x.table)
end

function A_mul_B!(C, A, B)
    @turbo for n ∈ indices((C,B), 2), m ∈ indices((C,A), 1)
        Cmn = zero(eltype(C))
        for k ∈ indices((A,B), (2,1))
            Cmn += A[m,k] * B[k,n]
        end
        C[m,n] = Cmn
    end
end
function sumdiag(x,n)
    s = 0
    @turbo for i in 1:n
    s += x[i,i]
    end
    return s
    end
function sumdiag(x,y,n)
    s = 0
    r = 0
    @turbo for i in 1:n
        s += x[i,i]
        r += y[i,i] - x[i,i]
    end
    return s,r
end

function doube_center(D)
    A = D .- mean(D,dims = 1)
    return A .- mean(A,dims = 2)
end

function QRfit(x,y)
    Q, R = LinearAlgebra.qr(x)
    R_inv_Q_trans = R^(-1)*Q'
    β = R^(-1)*Q' *y
    fitted =   x * β
    residuals = y .-fitted
    return sum(Diagonal(fitted) ), sum(Diagonal(residuals) ),R_inv_Q_trans
end

function get_fitted_residuals(R_inv_Q_trans,x,y,c,n)
    β = R_inv_Q_trans*y
    A_mul_B!(c, x, β)

    return sumdiag(c,y,n)
end
function get_fitted(R_inv_Q_trans,x,y,c,n)
    β = R_inv_Q_trans *y
    A_mul_B!(c, x, β)
    return sumdiag(c,n)
end

function term_fit(x,y)
    Q, R = LinearAlgebra.qr(x)
    R_inv_Q_trans = R^(-1)*Q'
    β = R_inv_Q_trans*y
    fitted =   x * β
    return sum(Diagonal(fitted)), R_inv_Q_trans
end

function unpack_formula(form)
    rhs = form.rhs
    lhs = form.lhs
    if typeof(rhs) == Term # if theres only a single rhs term
        return [form], [string(rhs)]
    else
        return     FormulaTerm.([lhs],rhs), string.(rhs)
    end
end

function permanova(data,D, formula = @formula(1~1), n_perm = 999;blocks = false)

    n = size(D)[1]
    formulae,terms = unpack_formula(formula) #list of model formulae for analysis and names for coef table
    try
        inds = vec(sum(ismissing.(Array(data[!,terms])),dims = 2) .==0)
        cc = sum(inds) # no. complete cases
        if cc < n
            data = data[(1:n)[inds],:]
            D = D[inds,inds]
            @warn "$(n-cc) data row(s) dropped due to missing values."
            n = cc
        end
    catch e
        @warn "Unable to check for missing values. Please make sure there are no missing values in relevant variables."
    end

    data = float.(data)

    G = Hermitian(-0.5 * doube_center(float.(D) .^2))
    fitted, residual = 0.0,0.0
    n_terms = length(terms)
    mod_mats = Vector{Matrix}(undef,n_terms)
    R_inv_Q_trans = Vector{Matrix}(undef,n_terms)
    Df = zeros(Int64,n_terms)
    sumsq = zeros(Float64,n_terms)
    term_matrix = ones(Int64,n,1) # just the intercept (not run, but added to later matrices)
    @inbounds for i in 1:n_terms
        mat_size = size(term_matrix)[2]
        formula = formulae[i]
        term_matrix = hcat(term_matrix,StatsModels.modelmatrix(formula,data))
        mod_mats[i] =term_matrix # store model matrices for permutation
        Df[i] = size(term_matrix)[2] -mat_size
        if i == n_terms # if full model
            fitted, residual,R_inv_Q_trans[i] = QRfit(term_matrix,G)
            sumsq[i] =fitted-sum(sumsq)
        else
        tmp,R_inv_Q_trans[i] = term_fit(term_matrix,G)
        sumsq[i] =tmp-sum(sumsq)
        end
    end
    Tot = fitted + residual
    Resid_Df = n - 1- sum(Df)
    f_terms = (sumsq ./Df) ./(residual /Resid_Df)
    r2 = sumsq ./Tot
    C = Vector{Array{Float64}}(undef,n_terms)
    for i in 1:n_terms
        C[i] = Array{Float64}(undef,size(mod_mats[i])[1],size(G)[2])
    end
    if blocks == false
        p = permute(G, n, n_terms, mod_mats,R_inv_Q_trans,n_perm,C)
    else
        block_mat = StatsModels.modelmatrix(blocks,data)
        p = permute(G, n, n_terms, mod_mats,R_inv_Q_trans,n_perm,C,block_mat)
    end

    regtab = get_output(terms,Df,sumsq,r2,f_terms,residual,Tot,p,n)
    return PSummary(regtab[1],regtab[2])

end

function permanova(data,M,metric ::DataType, formula= @formula(1~1), n_perm = 999;blocks = false)
    D = pairwise(metric(),M',M')
return  permanova(data,D, formula,n_perm , blocks = blocks)
end

hydra = permanova
function permute(G ::Hermitian, n, n_terms , mod_mats ,R_inv_Q_trans,n_perm , C)
    inds = collect(1:n)
    indscopy = copy(inds)
    perms = Array{Float64}(undef,n_terms,n_perm+1)
    fit = zeros(n_terms)
    f_terms = zeros(n_terms)
    res = 0.0
    g = Array(G)
    @inbounds for j in 1:n_perm +1
        g .= g[inds,inds]
       fit .= zeros(n_terms)
       @inbounds for i in 1:n_terms
            if i == n_terms
                prevsum = sum(fit)
                fit[i], res = get_fitted_residuals(R_inv_Q_trans[i],mod_mats[i],g,C[i],n)
                fit[i] -=prevsum
            else
                fit[i] = get_fitted(R_inv_Q_trans[i],mod_mats[i],g,C[i],n) - sum(fit)
            end
    end
    f_terms .= (fit) ./(res)
    perms[:,j] .= f_terms
    shuffle!(inds)
end
p = sum( perms[:,1] .<=perms[:,1:end],dims = 2) ./(n_perm+1)
return p
end

function permute(G ::Hermitian, n , n_terms , mod_mats ,R_inv_Q_trans,n_perm , C,block_mat)
    blockviews = []
    inds = collect(1:n)
    for col in eachcol(block_mat)
        push!(blockviews,view(inds,Bool.(collect(col))))
    end
    perms = Array{Float64}(undef,n_terms,n_perm+1)
    fit = zeros(n_terms)
    f_terms = zeros(n_terms)
    res = 0.0
    g = Array(G)
    @inbounds for j in 1:n_perm +1
        g .= g[inds,inds]
       fit .= zeros(n_terms)
       @inbounds for i in 1:n_terms
        if i == n_terms
            prevsum = sum(fit)
            fit[i], res = get_fitted_residuals(R_inv_Q_trans[i],mod_mats[i],g,C[i],n)
            fit[i] -=prevsum
        else
            fit[i] = get_fitted(R_inv_Q_trans[i],mod_mats[i],g,C[i],n) - sum(fit)
        end
    end
    f_terms .= (fit) ./(res)
    perms[:,j] .= f_terms
    shuffle!.(blockviews)

end
p = sum( perms[:,1] .<=perms[:,1:end],dims = 2) ./(n_perm +1)
return p
end


text:
import_all() = import_all(Main)

function import_all(m::Module)::Nothing
    package_list::Vector{String} = sort(unique(strip.(_package_list())))
    for p in package_list
        Base.eval(
            m,
            Base.Meta.parse(
                string(
                    "try ",
                    "import $(string(p)); ",
                    "catch e1 ",
                    "@debug(\"ignoring exception: \", e1,); ",
                    "try ",
                    "import Pkg; ",
                    "Pkg.add(\"$(string(p))\"); ",
                    "catch e2 ",
                    "@debug(\"ignoring exception: \", e2,); ",
                    "end ",
                    "end ",
                    ),
                ),
            )
        Base.eval(
            m,
            Base.Meta.parse(
                string(
                    "try ",
                    "import $(string(p)); ",
                    "@debug(\"imported $(string(p))\"); ",
                    "catch e3 ",
                    "@error(\"ignoring exception: \", e3,); ",
                    "end ",
                    ),
                ),
            )
    end
end


text:
using FactCheck

include("timearray.jl")
include("dataframe.jl")
include("auth_token.jl")

exitstatus()


text:
import Distributions
const D = Distributions
import Base: +, iterate

struct TimeInfo
    time::Float64
    bytes::Int64
    gctime::Float64
    memallocs::Base.GC_Diff
end

const notime = TimeInfo(0.0, 0, 0.0, Base.GC_Diff(0, 0, 0, 0, 0, 0, 0, 0, 0))

macro timeinfo(expr)
    quote
        val , time, bytes, gctime, memallocs = @timed $(esc(expr))
        val, TimeInfo(time, bytes, gctime, memallocs)
    end
end

maximumby(f, collection) = mapreduce(x -> (v=x, f=f(x)),
                                     (a, b) -> ifelse(isless(b.f, a.f), a, b),
                                     collection).v

randview(collection, n) = view(collection, rand(eachindex(collection), n))

# see: https://docs.python.org/3/library/itertools.html#itertools-recipes
repeatfunc(f, args...) = Iterators.map(((),) -> f(args...), Iterators.repeated(()))
repeatfunc(f, n, args...) = Iterators.map(((),) -> f(args...), Iterators.repeated((), n))

preparecache!(cache, n) = sizehint!(empty!(cache), n)

# see https://github.com/JuliaLang/julia/blob/master/doc/src/devdocs/ast.md
# using MacroTools
# macro parametrized(fundef)
#     fsplit = splitdef(fundef)
#     tname = gensym(fsplit[:name])
#     obj = gensym(:fp)
#     arg = gensym(:x)
#     funfield = gensym(:fun)
#     kwargs = [gensym(kw) for kw in fsplit[:kwargs]]

#     _struct = Expr(:struct, true, # mutable = true
#                    Expr(:(<:), tname, :Function),
#                    Expr(:block, funfield, kwargs...))
#     _call = Expr(:(=), Expr(:call, Expr(:(::), obj, tname), arg),
#                  Expr(:call,
#                       Expr(:(.), obj, Meta.quot(funfield)),
#                       arg,
#                       (Expr(:(.), obj, Meta.quot(kw)) for kw in kwargs)...))
#     Expr(:block, _struct, _call)
# end

# Will work if `D` is a location family.
struct Shifted{F<:D.VariateForm,
               S<:D.ValueSupport,
               T<:D.Distribution{F, S},
               E} <: D.Distribution{F, S}
    dist::T
    δ::E

    function Shifted(dist::T, δ) where {F<:D.VariateForm, S<:D.ValueSupport, T<:D.Distribution{F, S}}
        E = eltype(dist)
        new{F, S, T, E}(dist, convert(E, δ))
    end
end

D.rand(d::Shifted) = D.rand(d.dist) + d.δ
D.pdf(d::Shifted, x::Real) = D.pdf(d.dist, x - d.δ)
D.logpdf(d::Shifted, x::Real) = D.logpdf(d.dist, x - d.δ)
D.cdf(d::Shifted, x::Real) = D.cdf(d.dist, x - d.δ)
D.quantile(d::Shifted, q::Real) = D.quantile(d.dist, q) + d.δ
D.minimum(d::Shifted) = D.minimum(d.dist) + d.δ
D.maximum(d::Shifted) = D.maximum(d.dist) + d.δ
D.insupport(d::Shifted, x::Real) = D.insupport(d.dist, x - d.δ)
D.mean(d::Shifted) = D.mean(d.dist) + d.δ
D.var(d::Shifted) = D.var(d.dist)
D.modes(d::Shifted) = D.modes(d.dist) .+ d.δ
D.mode(d::Shifted) = D.mode(d.dist) + d.δ
# StatsBase.skewness(d::Shifted) = StatsBase.skewness(d.D)
# StatsBase.kurtosis(d::Distribution, correction::Bool) = StatsBase.kurtosis(d.D, correction)
# StatsBase.entropy(d::Shifted, b::Real) = StatsBase.entropy(d.D, b)
# mgf(d::Shifted, ::Any) = ???
# cf(d::Shifted, ::Any) = ???


text:
mutable struct DiscreteSystem
    probability::Vector{Vector{Int64}}
    table::Vector{Array{Float64, N} where N}
    incidence::SparseMatrixCSC{Float64,Int64}
    incidenceTranspose::SparseMatrixCSC{Float64,Int64}
end

mutable struct DiscreteTreeGraph
    Nvariable::Int64
    Nfactor::Int64
    Nlink::Int64
    root::Int64
    virtualMessage::Float64
    state::Array{Int64,1}
    messageDirect::Vector{Vector{Float64}}
    rowForward::Vector{Vector{Int64}}
    rowBackward::Vector{Vector{Int64}}
    colForward::Vector{Vector{Int64}}
    colBackward::Vector{Vector{Int64}}
    incomingToFactor::Vector{Vector{Int64}}
    incomingToVariable::Vector{Vector{Int64}}
    iterateFactor::Array{Int64,1}
    iterateVariable::Array{Int64,1}
    passFactorVariable::Int64
    passVariableFactor::Int64
    forward::Bool
    backward::Bool
end

mutable struct DiscreteInference
    fromFactor::Array{Int64,1}
    toVariable::Array{Int64,1}
    messageFactorVariable::Vector{Vector{Float64}}
    fromVariable::Array{Int64,1}
    toFactor::Array{Int64,1}
    messageVariableFactor::Vector{Vector{Float64}}
    marginal::Vector{Vector{Float64}}
end

struct DiscreteTreeModel
    graph::DiscreteTreeGraph
    inference::DiscreteInference
    system::DiscreteSystem
end

########## Form a graph with continuous variables and initialize messages ##########
function discreteTreeModel(
    args...;
    virtual::Float64 = 1.0,
    root::Int64 = 1)

    system = readDiscreteArguments(args)
    graph, inference = makeDiscreteTreeGraph(system, virtual, root)

    return DiscreteTreeModel(graph, inference, system)
end

########## Read discrete system model ##########
function readDiscreteArguments(args)
    if isa(args[1], Dict)
        probability = []; table = []
        @inbounds for (key, value) in args[1]
            if key == "probability"
                for k in value
                    push!(probability, k)
                end
            end
            if key == "table"
                for k in value
                    push!(table, k)
                end
            end
        end
    else
        probability = args[1]
        table = args[2]
    end

    incidence, incidenceTranspose = incidenceMatrix(probability)

    return DiscreteSystem(probability, table, incidence, incidenceTranspose)
end

########## Incidence matrix ##########
function incidenceMatrix(probabilities)
    row = Int64[]; col = Int64[]

    @inbounds for (k, probability) in enumerate(probabilities)
        for variable in probability
            push!(row, k)
            push!(col, variable)
        end
    end
    incidence = sparse(row, col, 1)
    incidenceTranspose = sparse(col, row, 1)

    return incidence, incidenceTranspose
end

########## Produce the graphical model ##########
function makeDiscreteTreeGraph(system, virtualMessage, root)
    ### Number of factor and variable nodes
    Nfactor, Nvariable = size(system.incidence)

    ### Pass through probabilities
    Nlink = 0
    states = fill(0, Nvariable)
    messageDirect = [Float64[] for i = 1:Nvariable]
    rowForward =  [Int64[] for i = 1:Nfactor]
    colForward = [Int64[] for i = 1:Nvariable]
    @inbounds for (k, probability) in enumerate(system.probability)
        for (kk, variable) in enumerate(probability)
            if states[variable] == 0
                states[variable] = size(system.table[k], kk)
                messageDirect[variable] = fill(virtualMessage, states[variable])
            elseif states[variable] != size(system.table[k], kk)
                error("The possible states of the variable $variable do not agree through tables.")
            end
            if length(probability) != 1
                push!(rowForward[k], variable)
                push!(colForward[variable], k)
                Nlink += 1
            else
                messageDirect[variable] .*= system.table[k]
            end
        end
    end
    rowBackward = deepcopy(rowForward)
    colBackward = deepcopy(colForward)

    ### Initialize the first forward step
    iterateVariable = fill(0, Nvariable)
    counter = 0
    @inbounds for col = 1:Nvariable
        if length(colForward[col]) == 1 && col != root
            counter += 1
            iterateVariable[counter] = col
        end
    end
    resize!(iterateVariable, counter)

    ### Initialize data
    fromVariable = fill(0, Nlink)
    toFactor = fill(0, Nlink)
    messageVariableFactor = [Float64[] for i = 1:Nlink]

    fromFactor = fill(0, Nlink)
    toVariable = fill(0, Nlink)
    messageFactorVariable = [Float64[] for i = 1:Nlink]

    marginal = [Float64[] for i = 1:Nvariable]
    incomingToVariable = [Int64[] for i = 1:Nvariable]
    incomingToFactor = [zeros(length(system.probability[i])) for i = 1:Nfactor]
    passFactorVariable = 0
    passVariableFactor = 0
    iterateFactor = Int64[]

    return DiscreteTreeGraph(Nvariable, Nfactor, Nlink, root, virtualMessage, states, messageDirect,
        rowForward, rowBackward, colForward, colBackward, incomingToFactor, incomingToVariable,
        iterateFactor, iterateVariable, passFactorVariable, passVariableFactor, true, true),
           DiscreteInference(fromFactor, toVariable, messageFactorVariable, fromVariable, toFactor, messageVariableFactor, marginal)
end


text:
using CloudEvents
using Test

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