text:
"""
    find_graph_component(gfa_result::BioGraph.GFAResult)

Return graph components of `GFAResult`:

- Simple Graphs
- Lone Cycles
- Lone Nodes
"""
function find_graph_component(gfa_result::GFAResult)
    source_nodes::Array{Int64} = Int64[]
    sink_nodes::Array{Int64} = Int64[]
    start_nodes::Array{Int64} = Int64[]
    end_nodes::Array{Int64} = Int64[]

    g = gfa_result.g
    w = gfa_result.w
    l = gfa_result.l
    e = gfa_result.e
    p = gfa_result.p

    l_dict, l_s_dict, l_i_dict = _label_map_array_to_dict(l)
    w_dict = _weight_array_to_dict(w)
    e_o_dict, e_i_dict = _edge_label_array_to_dict(e)

    for e::AbstractEdge in edges(g)
        push!(start_nodes, e.src)
        push!(end_nodes, e.dst)
    end

    sink_nodes = collect(intersect(setdiff!(Set(1:nv(g)), Set(start_nodes)), Set(end_nodes)))
    source_nodes = collect(intersect(setdiff!(Set(1:nv(g)), Set(end_nodes)), Set(start_nodes)))

    g_com = weakly_connected_components(g)

    graphs = GraphComponent()
    graphs.lone_node = []
    graphs.lone_cycle = []
    graphs.graph = []
    for com in g_com
        if length(com) == 1
            g_dummy = GraphResult()
            g_dummy.weight_inherited = gfa_result.weight_inherited
            sg, vmap = induced_subgraph(g, com)
            g_dummy.edge_label = []
            g_dummy.g = sg
            g_dummy.path = []
            if l == []
                g_dummy.node_label = [NodeLabel(1, string(vmap[1]), "", "")]
                g_dummy.source_node = [NodeLabel(1, string(vmap[1]), "", "")]
                g_dummy.sink_node = [NodeLabel(1, string(vmap[1]), "", "")]
            else
                g_dummy.node_label = [NodeLabel(1, l_dict[vmap[1]], l_s_dict[vmap[1]], l_i_dict[vmap[1]])]
                g_dummy.source_node = [NodeLabel(1, l_dict[vmap[1]], l_s_dict[vmap[1]], l_i_dict[vmap[1]])]
                g_dummy.sink_node = [NodeLabel(1, l_dict[vmap[1]], l_s_dict[vmap[1]], l_i_dict[vmap[1]])]
            end
            if w == []
                g_dummy.weight = [Weight(1, 1)]
            else
                g_dummy.weight = [Weight(1, w_dict[vmap[1]])]
            end
            g_dummy.sg_map = [SubGraphMap(1, vmap[1])]
            push!(graphs.lone_node, g_dummy)
        else
            sg_source = collect(intersect(Set(source_nodes), Set(com)))
            sg_sink = collect(intersect(Set(sink_nodes), Set(com)))
            g_dummy = GraphResult()
            g_dummy.weight_inherited = gfa_result.weight_inherited
            sg, vmap = induced_subgraph(g, com)
            g_dummy.g = sg
            g_dummy.sg_map = []
            for i in 1:nv(sg)
                push!(g_dummy.sg_map, SubGraphMap(i, vmap[i]))
            end
            sg_dict = _sg_map_array_to_dict(g_dummy.sg_map)
            reverse_sg_dict = Dict(value => key for (key, value) in sg_dict)

            g_dummy.node_label = []
            g_dummy.source_node = []
            g_dummy.sink_node = []
            g_dummy.weight = []
            g_dummy.edge_label = []
            g_dummy.path = []

            if p != []
                for p_dummy in p
                    if p_dummy.path[1] in vmap
                        p_arr = []
                        for p_index in p_dummy.path
                            push!(p_arr, reverse_sg_dict[p_index])
                        end
                        push!(g_dummy.path, Path(p_dummy.name,p_arr))
                    end
                end
            end

            if l == []
                for i in 1:nv(sg)
                    push!(g_dummy.node_label, NodeLabel(i, string(vmap[i]), "", ""))
                end
                for source_node in sg_source
                    push!(g_dummy.source_node, NodeLabel(reverse_sg_dict[source_node], string(vmap[reverse_sg_dict[source_node]]), "", ""))
                end
                for sink_node in sg_sink
                    push!(g_dummy.sink_node, NodeLabel(reverse_sg_dict[sink_node], string(vmap[reverse_sg_dict[sink_node]]), "", ""))
                end
            else
                for i in 1:nv(sg)
                    push!(g_dummy.node_label, NodeLabel(i, l_dict[vmap[i]], l_s_dict[vmap[i]], l_i_dict[vmap[i]]))
                end
                for source_node in sg_source
                    push!(g_dummy.source_node, NodeLabel(reverse_sg_dict[source_node], l_dict[vmap[reverse_sg_dict[source_node]]], l_s_dict[vmap[reverse_sg_dict[source_node]]], l_i_dict[vmap[reverse_sg_dict[source_node]]]))
                end
                for sink_node in sg_sink
                    push!(g_dummy.sink_node, NodeLabel(reverse_sg_dict[sink_node], l_dict[vmap[reverse_sg_dict[sink_node]]], l_s_dict[vmap[reverse_sg_dict[sink_node]]], l_i_dict[vmap[reverse_sg_dict[sink_node]]]))
                end
            end
            if w == []
                for i in 1:nv(sg)
                    push!(g_dummy.weight, Weight(i, 1))
                end
            else
                for i in 1:nv(sg)
                    push!(g_dummy.weight, Weight(i, w_dict[vmap[i]]))
                end
            end
            if e == []
                for edge in edges(sg)
                    push!(g_dummy.edge_label, EdgeLabel(string(sg_dict[edge.src]), string(sg_dict[edge.dst]), "", ""))
                end
            else
                for edge in edges(sg)
                    e_name = l_dict[sg_dict[edge.src]] * "::" * l_dict[sg_dict[edge.dst]]
                    push!(g_dummy.edge_label, EdgeLabel(l_dict[sg_dict[edge.src]], l_dict[sg_dict[edge.dst]], e_o_dict[e_name], e_i_dict[e_name]))
                end
            end
            if sg_sink == [] || sg_source == []
                push!(graphs.lone_cycle, g_dummy)
            else
                push!(graphs.graph, g_dummy)
            end
        end
    end

    return graphs
end

"""
    get_summary(g_coms::BioGraph.GraphComponent)

Summary of `GraphComponent`.
"""
function get_summary(g_coms::GraphComponent)
    println("No of simple graphs: " * string(length(g_coms.graph)))
    println("No of lone cycles: " * string(length(g_coms.lone_cycle)))
    println("No of lone nodes: " * string(length(g_coms.lone_node)))
end

"""
    get_summary(g_coms::BioGraph.GraphResult)

Summary of `GraphResult`:

- No of vertices
- No of edges
- No of source nodes
- No of end nodes
- No of path
"""
function get_summary(g_result::GraphResult)
    println("No of vertices: " * string(nv(g_result.g)))
    println("No of edges: " * string(ne(g_result.g)))
    println("No of source_node: " * string(length(g_result.source_node)))
    println("No of sink_node: " * string(length(g_result.sink_node)))
    println("No of path: " * string(length(g_result.path)))
end

"""
    get_terminus(g_result::BioGraph.GraphResult; outfile::String)

Get all source and sink nodes of `GraphResult`. Write to CSV, including data from GFA.
"""
function get_terminus(g_result::GraphResult, outfile::String="")
    if outfile != ""
        sink_name = outfile * "_sink_node.csv"
        source_name = outfile * "_source_node.csv"
        open(sink_name, "w") do io
            write(io, "Label,Direction,Sequence,Information" * "\n")
            for node in g_result.sink_node
                s_dummy = split(node.label, "\t")
                write(io, s_dummy[1] * "," * s_dummy[2] * "," * node.sequence * "," * node.info * "\n")
            end
        end
        open(source_name, "w") do io
            write(io, "Label,Direction,Sequence,Information" * "\n")
            for node in g_result.source_node
                s_dummy = split(node.label, "\t")
                write(io, s_dummy[1] * "," * s_dummy[2] * "," * node.sequence * "," * node.info * "\n")
            end
        end
    end

    println("Source node: ")
    for node in g_result.source_node
        println("Node: " * string(node.node) * "\tLabel: " * node.label)
    end
    println("Sink node: ")
    for node in g_result.sink_node
        println("Node: " * string(node.node) * "\tLabel: " * node.label)
    end
    return g_result.source_node, g_result.sink_node
end

"""
    get_gfa(g_result::BioGraph.GraphResult; outfile::String)

Get GFA output of `GraphResult`.
"""
function get_gfa(g_result::GraphResult, outfile::String)
    open(outfile, "w") do io
        write(io, "H\tVN:Z:1.0" * "\n")
        for node in g_result.node_label
            s_dummy = split(node.label, "\t")
            write(io, "S\t" * s_dummy[1] * "\t" * node.sequence * "\t" * node.info * "\n")
        end
        for edge in g_result.edge_label
            write(io, "L\t" * edge.src * "\t" * edge.dst * "\t" * edge.overlap * "\t" * edge.info * "\n")
        end
    end
end

ext: jl


text:
julia> [1 1 ; 1 0]^10
2x2 Array{Int64,2}:
 89  55
 55  34

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("NfRel/NfRel.jl")
include("NfRel/NfRel_ns.jl")
# include("NfRel/NfRelTypes.jl")
include("NfRel/NfRelOrd.jl")
include("NfRel/NfRelOrdElem.jl")
include("NfRel/NfRelOrdIdl.jl")
include("NfRel/NfRelOrdFracIdl.jl")
include("Map/NfRelOrd.jl")


text:
using DFTShims
using Base.Test
using Unitful
using AxisArrays

const DD = DFTShims.Dispatch.Dimensions
const DH = DFTShims.Dispatch.Hartree
const UH = DFTShims.UnitfulHartree

@testset "Scalars" begin
    @test 1u"ρ" isa DD.Scalars.ρ
    @test 1u"nm^-3" isa DD.Scalars.ρ
    @test 1u"ρ" isa DH.Scalars.ρ
    @test !(1u"nm^-3" isa DH.Scalars.ρ)

    @test DH.Scalars.ρ <: DD.Scalars.ρ
    @test DH.Scalars.ρ{Float32} <: DD.Scalars.ρ
    @test DH.Scalars.ρ{Float32} <: DD.Scalars.ρ{Float32}
    @test DH.Scalars.ρ{Float32} === DD.Scalars.ρ{Float32, typeof(UH.ρ)}
    @test DH.Scalars.ρ{Float32} <: DD.Scalars.All{Float32}
    @test !(DH.Scalars.ρ{Int32} <: DD.Scalars.All{Float32})
    @test !(DH.Scalars.ρ <: DD.Scalars.All{Float32})

    f = (::DH.Scalars.ρ{<: Integer}) -> :hartree_with_integer
    (::typeof(f))(::DH.Scalars.ρ) = :hartree_with_any_number
    (::typeof(f))(::DD.Scalars.ρ{<: Integer}) = :density_with_integer
    (::typeof(f))(::DD.Scalars.ρ) = :density_with_any_number
    @test f(1u"ρ") == :hartree_with_integer
    @test f(1.0u"ρ") == :hartree_with_any_number
    @test f(1u"nm^-3") == :density_with_integer
    @test f(1.0u"nm^-3") == :density_with_any_number
end

@testset "Arrays" begin
    @test [1u"ρ", 2.0u"ρ"] isa DFTShims.DFTArray
    @test [1u"ρ", 2.0u"ρ"] isa DD.DensityArray
    @test [1u"ρ", 2.0u"ρ"] isa DD.Arrays.ρ
    @test [1u"ρ", 2.0u"ρ"] isa DD.Arrays.ρ{Float64}
    @test [1u"ρ", 2.0u"ρ"] isa DH.DensityDenseArray{Float64, 1}
    @test [1u"ρ", 2.0u"ρ"] isa DH.DenseArrays.ρ{Float64, 1}
    @test [1u"nm^-3", 2.0u"nm^-3"] isa DD.DenseArrays.ρ{Float64, 1}
    @test !([1u"nm^-3", 2.0u"nm^-3"] isa DH.DenseArrays.ρ{Float64, 1})

    @test DH.DensityDenseArray{Int8, 2} === DD.DensityDenseArray{Int8, 2, typeof(u"ρ")}
    @test DH.Arrays.ρ{Float32} <: DD.Arrays.All{Float32}
    @test !(DH.Arrays.ρ{Int32} <: DD.Arrays.All{Float32})
    @test !(DH.Arrays.ρ <: DD.Arrays.All{Float32})
    @test DH.Arrays.ρ{Float32, 2} <: DD.Arrays.All{Float32, 2}
    @test [1u"nm^-3", 2.0u"nm^-3"] isa DD.Arrays.ρ{Float64, 1}
    @test !([1u"nm^-3", 2.0u"nm^-3"] isa DD.Arrays.ρ{Float64, 2})

    f = (::DH.Arrays.ρ{<: Integer}) -> :hartree_with_integer
    (::typeof(f))(::DH.Arrays.ρ) = :hartree_with_any_number
    (::typeof(f))(::DD.Arrays.ρ{<: Integer}) = :density_with_integer
    (::typeof(f))(::DD.Arrays.ρ) = :density_with_any_number
    (::typeof(f))(::DD.AxisArrays.ρ) = :axis_density_with_any_number
    @test f([1u"ρ"]) == :hartree_with_integer
    @test f([1.0u"ρ"]) == :hartree_with_any_number
    @test f([1u"nm^-3"]) == :density_with_integer
    @test f([1.0u"nm^-3"]) == :density_with_any_number
    @test f(AxisArray([1.0u"nm^-3"])) == :axis_density_with_any_number
end

max_stars_repo_name: mdavezac/DFTShims

max_issues_repo_name: mdavezac/DFTShims

max_forks_repo_name: mdavezac/DFTShims


text:
using RigidBodySim, Documenter, RigidBodyDynamics

makedocs(
    doctest = true,
    modules = [RigidBodySim],
    checkdocs = :none, # because of the reexported functionality
    format = :html,
    sitename ="RigidBodySim.jl",
    linkcheck = true,
    authors = "Twan Koolen and contributors.",
    pages = [
        "Home" => "index.md",
        "Quick start guide" => "quickstart.md",
        "Details" => "details.md"
    ]
)

deploydocs(
    deps = nothing,
    repo = "github.com/JuliaRobotics/RigidBodySim.jl.git",
    target = "build",
    make = nothing,
    julia = "1.0",
    osname = "linux"
)


text:
using LazySets, Reachability, Polyhedra, Optim
import LazySets.HalfSpace
include("DriveTrain.jl");

# compatibility of julia versions
if VERSION >= v"0.7"
    using LinearAlgebra, SparseArrays
end

N = Float64

function get_projection(sol::AbstractSolution, system_dimension::Int64) ::AbstractSolution
    sol_processed =  Reachability.ReachSolution([Reachability.ReachSet{CartesianProductArray{N}, N}(
                CartesianProductArray{N, HPolytope{N}}([LazySets.Approximations.overapproximate(rs.X, LazySets.Approximations.OctDirections(system_dimension))]),
                rs.t_start, rs.t_end) for rs in sol.Xk], sol.options)

    sol_proj = Reachability.ReachSolution(Reachability.project_reach(
        sol_processed.Xk, [1,2], system_dimension, sol.options), sol.options);

    return sol_proj;
end

function get_solution(drivetrain, θ, opD, t, max_jumps)
    # initial condition in mode 1
    c = vcat([-0.0432, -11, 0, 30, 0, 30, 360], repeat([-0.0013, 30], θ))
    g = zeros(2*θ + 7,2*θ + 7)
    g[1,1], g[2,2], g[3,3], g[4,4], g[5,5], g[6,6], g[7,7] = 0.0056, 4.67, 0, 10, 0, 10, 120
    i = 1;
    while (i <= θ)
        g[2*i+6,2*i+6], g[2*i+7,2*i+7] = 0.0006, 10
        i+=1
    end
    X0 = Zonotope(c, g)

    system = InitialValueProblem(drivetrain, [(1, X0)]);
    plot_vars = [1, 3]
    system_dimension = 2*θ+7
    options = Options(:mode=>"reach",:vars=>1:system_dimension, :T=>t, :δ=>0.01,
                      :max_jumps=>max_jumps, :plot_vars=>plot_vars,
                      :ε_proj=>0.001, :verbosity=>0, :project_reachset=>false);

    # default algorithm
    @time begin
        sol = solve(system, options, Reachability.BFFPSV18(), opD);
    end
end

opDs = [
        (Reachability.ReachSets.LazyDiscretePost(),  "L", 64)
        (Reachability.ReachSets.ApproximatingDiscretePost(), "A", 64)
       ];

#warmup run for each opD for low dimension
println("Warmup run just after restart REPL")
for (opD, name, upper_bound) in opDs
    θ = 1;
    while (θ <= 2)
        HS = drivetrain(θ)
        sol = get_solution(HS, θ, opD, 2., 4);
        θ += 1;
    end
end
println("End of warmup run")

project_and_store = false
if project_and_store
    println("Note: using output projection and storage")
else
    println("Note: skipping output projection and storage")
end

results = Vector{Tuple{AbstractSolution, String}}()
θ = 1
while θ <= 47
    for (opD, name, upper_bound) in opDs
        println("**********************************")
        if θ > upper_bound
            continue
        end
        println("$name $(n0)")
        T = 2.
        max_jumps = 4
        end
        HS = drivetrain(θ)
        sol = get_solution(HS, θ, opD, T, max_jumps);
        if project_and_store
            sol_proj = get_projection(sol, 2*θ+7);
            push!(results, (sol_proj, name));
        end
    end
    θ += 1;
end

return results

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
using Nettle

function parse_file(config_file)
    local x
    open( config_file ) do io
        x = JSON.parse(io)
    end
    for key in keys(x)
        setindex!(configuration, getindex(x, key), key )
    end
    return nothing
end # function

function write_file(config_file, x)
    open( config_file, "w" ) do io
        JSON.print( io, x )
    end # do
    return nothing
end # function

function register( path )
    if path ∈ configuration["paths"]
        return
    end
    push!(configuration["paths"], path)
    return nothing
end # function

function unregister( path )
    path_ind = findfirst(path, configuration["paths"])
    if path_ind === nothing
        @info "Path $path not found.\nRegistered paths are:\n$(configuration["paths"])"
        return
    end
    deleteat!( configuration["paths"], path_ind )
    return nothing
end # function

function reset_paths()
    configuration["paths"] = String[]
    return nothing
end # function

function add_token( pat; key = key[] )
    configuration["token"] = String( encrypt(
        Encryptor("AES256", key),
        add_padding_PKCS5(Vector{UInt8}(pat),16) ) )
    return nothing
end # function


text:
using Distributions, StatsBase, LinearAlgebra, Random
Random.seed!(123)
function histVals(μ,σ; npts=1000, izq = -8, der = 6, paso = 0.15)
    dist = rand(Normal(μ,σ), npts)
    h = fit(Histogram, dist, izq:paso:der) #nbins = 50
    h = normalize(h)
    edg = collect(h.edges[1])
    centros = (edg[2:end] .+ edg[1:end-1])/2
    centros, h.weights
end
@gp histVals(0.25,0.5)... "w boxes t '0.25,0.5' lc '#E69F00'"
@gp :- "set key title 'μ,σ' box 3" xlab = "x" ylab = "P(x)"
@gp :- "set style fill transparent solid 0.45 noborder"
@gp :- histVals(-1,2)... "w histeps t '-1,2' lc 'black'"
@gp :- histVals(2,1)... "w fillsteps t '2,1' lc '#009E73' fs transparent solid 0.2 noborder"
save(term="pngcairo size 600,400", output="./code/plt_ex16.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
# Text provided under a Creative Commons Attribution license, CC-BY.  All code is made available under the FSF-approved BSD-3 license.  (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784.
# [@LorenaABarba](https://twitter.com/LorenaABarba)
# Translation to Julia by [@miguelraz_](https://twitter.com/miguelraz_)

## 12 steps to Navier–Stokes
# =====
# ***

# This lesson complements the first interactive module of the online [CFD Python](https://github.com/barbagroup/CFDPython) class, by Prof. Lorena A. Barba, called **12 Steps to Navier–Stokes.** It was written with BU graduate student Gilbert Forsyth.

## Array Operations with NumPy
# ----------------

# For more computationally intensive programs, the use of built-in Numpy functions can provide an  increase in execution speed many-times over.  As a simple example, consider the following equation:

# $$u^{n+1}_i = u^n_i-u^n_{i-1}$$

# Now, given a vector $u^n = [0, 1, 2, 3, 4, 5]\ \ $   we can calculate the values of $u^{n+1}$ by iterating over the values of $u^n$ with a for loop.

u = [0, 1, 2, 3, 4, 5]

for i in 1:length(u)
    println(u[i] - u[i-1])
end

# This is the expected result and the execution time was nearly instantaneous.
# We can calculate each operation with one function:

diff(u)

# What this command says is subtract the 1st, 2nd, 3rd, 4th and 5th elements of $u$ from the 2nd, 3rd, 4th, 5th and 6th elements of $u$.

### Speed Increases
# In NumPy, we would have to learn a bucketful about how to restructure our program to exploit "vectorization".
# This unfortunately makes code very brittle and diminishes readability, but it's require by NumPy for speed.
# We have no such obstacles in Julia - if you write simple for loops, they will start out being "fast enough",
# and then we can move on to talk about worthy optimizations later.
# Note: There

nₓ = 81
ny = 81 # note we can't call this n\_y because of the Unicode standard, not because of Julia
nₜ = 100
c = 1
Δx = 2 / (nx - 1)
Δy = 2 / (ny - 1)
σ = 0.2
Δt = σ * Δx

rangex = range(start = 0, stop = 2, length = nₓ)
rangey = range(start = 0, stop = 2, length = ny)

### TODO
u = ones(ny, nₓ) ##create a 1xn vector of 1's
uₙ = ones(ny, nₓ)

###Assign initial conditions
#We can try the following to start the initial conditions:
# TODO - cartesian indexes
rangey = Int(.5/Δy):Int(1/Δy)
rangex = Int(.5/Δx):Int(1/Δx)
for (_, j) in rangey
    for (_, i) in rangex
        u[i, j] = 2.0
    end
end
u

# But we could also try doing this with `CartesianIndexes` - they will subset a rectangular region of our array (for any dimension!)
# and loop efficiently through it. You can find more about them if you press `?CartesianIndex` in the REPL
#u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

u[CartesianIndices((rangex,rangey))] .= 2

# With our initial conditions all set up, let's first try running our original nested loop code, making use of the
# BenchmarkTools.jl `@btime` macro, which will help us evaluate the performance of our code.

# **Note**: The `@btime` macro will run the code several times and then give an average execution time as a result.
# If you have any figures being plotted within a code cell where you run `@btime`,
# it will plot those figures repeatedly which can be a bit messy and slow.

# The execution times below will vary from machine to machine.  Don't expect your times to match these times,
# but you _should_ expect to see the same general trend in decreasing execution time as we switch to array operations.

using BenchmarkTools

# TODO
@btime u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2 :setup=(u = ones(ny, nx))

@btime begin
    for n in 1:nₜ ##loop across number of time steps
    uₙ = copy(u)
    row, col = size(u)
    for j in 1:row
        for i in 1:col
            u[j, i] = (uₙ[j, i] - (c * Δt / Δx *
                                  (uₙ[j, i] - uₙ[j, i - 1])) -
                                  (c * Δt / Δy *
                                   (uₙ[j, i] - uₙ[j - 1, i])))
            u[1, :] = 1.0
            u[end, :] = 1.0
            u[:, 1] = 1.0
            u[:, end] = 1.0
        end
    end
    end
end :setup=(u = ones(ny, nx))

## Python comparison
# With the "raw" Python code above, the mean execution time achieved was 3.07 seconds (on a MacBook Pro Mid 2012).
# Keep in mind that with these three nested loops, that the statements inside the **j** loop are being evaluated more than 650,000 times.
# Let's compare that with the performance of the same code implemented with array operations:
# I've left the vectorized version so that people can compare the code and **appreciate** how substantial the difference from Julia is.

# ```python
# %%timeit
# u = numpy.ones((ny, nx))
# u[int(.5 / dy): int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2

# for n in range(nt + 1): ##loop across number of time steps
#     un = u.copy()
#     u[1:, 1:] = (un[1:, 1:] - (c * dt / dx * (un[1:, 1:] - un[1:, 0:-1])) -
#                               (c * dt / dy * (un[1:, 1:] - un[0:-1, 1:])))
#     u[0, :] = 1
#     u[-1, :] = 1
#     u[:, 0] = 1
#     u[:, -1] = 1
# ```

#     7.38 ms ± 105 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# As you can see, the speed increase is substantial in Python.  The same calculation goes from 3.07 seconds to 7.38 milliseconds.
# 3 seconds isn't a huge amount of time to wait, but these speed gains will increase exponentially
# with the size and complexity of the problem being evaluated.

## Julia timings:
# We're going to be doing some mild modifications to the code, which will involve
# - avoid global variables (even though Julia 1.8 will not penalize this), as it's good practice
# - diminishing allocations (because [views in Julia copy by default](TODO perf manual and views))
# - putting code inside of a function
# There's very little point to write efficient array code if it's not within a function in Julia.
# This is because Julia needs to JIT compile your code to get all the speedups.
# The JIT compiler does not need you to use type hints or annotations to work properly, it's only necessary that
# [TODO type stability]() function output types be predictable on input types only.
#
function time1(u, params)
    u = ones(ny, nₓ)
    (;Δt, Δx, c, nₜ) = params

    for n in 1:nₜ ##loop across number of time steps
    uₙ = copy(u)
            @views u[:,:] .= (uₙ[2:end, 2:end] - (c * Δt / Δx *
                                  (uₙ[2:end, 2:end] - uₙ[2:end, 1:end])) -
                                  (c * Δt / Δy *
                                   (uₙ[2:end, 2:end] - uₙ[1:end, 2:end])))
            @views u[1, :] .= 1.0
            @views u[end, :] .= 1.0
            @views u[:, 1] .= 1.0
            @views u[:, end] .= 1.0
        end
    end
    end
end
params = (;Δt, Δx, c, nₜ)
@btime time1(u, params)

max_stars_repo_name: mkitti/CFDJulia

max_issues_repo_name: mkitti/CFDJulia

max_forks_repo_name: mkitti/CFDJulia


text:
using Dates
using PrettyPrinting: PrettyPrinting, @isexpr

import Base: isempty, parse

macro unpack(ex)
    if @isexpr ex Expr(:struct, mut::Bool, decl, Expr(:block, args...))
        if @isexpr decl Expr(:(<:), T, _)
        else
            T = decl
        end
        struct_slots = Any[]
        ctr_slots = Any[]
        new_slots = Any[]
        unpack_slots = Any[]
        quoteof_slots = Any[]
        for arg in args
            if arg isa LineNumberNode
                push!(struct_slots, arg)
                continue
            end
            if @isexpr arg Expr(:(=), arg′, default)
                arg = arg′
                has_default = true
            else
                has_default = false
            end
            if @isexpr arg Expr(:call, :(=>), key, arg′)
                arg = arg′
                has_key = true
            else
                has_key = false
            end
            if @isexpr arg Expr(:(::), name, FT)
            else
                error("expected field declaration; got $(repr(arg))")
            end
            @isexpr FT Expr(:curly, :Union, FT, _)
            push!(struct_slots, arg)
            push!(ctr_slots, has_default ? Expr(:kw, name, default) : name)
            push!(new_slots, name)
            unpack_slot = Expr(:call, :unpack!, FT, :data)
            if has_key
                push!(unpack_slot.args, key)
                if has_default
                    push!(unpack_slot.args, default)
                end
            end
            unpack_slot = Expr(:kw, name, unpack_slot)
            push!(unpack_slots, unpack_slot)
            quoteof_slot = :(push!(ex.args, Expr(:kw, $(QuoteNode(name)), obj.$name)))
            if has_default
                quoteof_slot = :(obj.$name == $default || $quoteof_slot)
            end
            push!(quoteof_slots, quoteof_slot)
        end
        return quote
            struct $decl
                $(struct_slots...)
                $T(; $(ctr_slots...)) = new($(new_slots...))
            end
            unpack!(::Type{$T}, data::Dict) = $T($(unpack_slots...))
            function PrettyPrinting.quoteof(obj::$T)
                ex = Expr(:call, nameof($T))
                $(quoteof_slots...)
                ex
            end
        end |> esc
    else
        error("expected a struct; got $(repr(ex))")
    end
end

unpack!(T::Type, data::Dict, key::String, default) =
    haskey(data, key) ?
        something(unpack!(T, data, key), default) :
        default

function unpack!(T::Type, data::Dict, key::String)
    bucket = data[key]
    retval = unpack!(T, bucket)
    if !(bucket isa Union{Dict,Vector}) || isempty(bucket)
        delete!(data, key)
    end
    return retval
end

unpack!(::Type{String}, data) =
    data

function unpack!(T::Type{<:Union{Date, Number, Enum}}, data)
    if data isa String
        return parse(T, data)
    end
    return T(data)
end

function unpack!(::Type{Vector{T}}, items::Vector{Any}) where {T}
    retval = T[]
    for item in items
        push!(retval, unpack!(T, item))
    end
    filter!(item -> !isempty(item), items)
    return retval
end

@enum RangeOp GT GTE LT LTE EQ NEQ BT NBT
Base.parse(::Type{RangeOp}, s::String) =
    s == "gt" ? GT :
    s == "gte" ? GTE :
    s == "lt" ? LT :
    s == "lte" ? LTE :
    s == "eq" ? EQ :
    s == "!eq" ? NEQ :
    s == "bt" ? BT :
    s == "!bt" ? NBT :
    throw(DomainError(s, "Unknown Range Operation"))

@enum TextOp CONTAINS NCONTAINS STARTSWITH NSTARTSWITH ENDSWITH NENDSWITH
Base.parse(::Type{TextOp}, s::String) =
    s == "contains" ? CONTAINS :
    s == "!contains" ? NCONTAINS :
    s == "startsWith" ? STARTSWITH :
    s == "!startsWith" ? NSTARTSWITH :
    s == "endsWith" ? ENDSWITH :
    s == "!endsWith" ? NENDSWITH :
    throw(DomainError(s, "Unknown Text Operation"))

@unpack struct DateRange
    "Value" => value::Date
    "Op" => op::RangeOp
    "Extent" => extent::Union{Date, Nothing} = nothing
end

@unpack struct TextFilter
    "Text" => text::String = "null"
    "Op" => op::TextOp
end

@unpack struct NumericRange
    "Value" => value::Number
    "Op" => op::RangeOp
    "Extent" => extent::Union{Number, Nothing} = nothing
end

@enum InvalidReasonFlag UNKNOWN_REASON VALID INVALID
InvalidReasonFlag(::Nothing) = UNKNOWN_REASON
Base.parse(::Type{InvalidReasonFlag}, s::Union{String, Nothing}) =
    s == "V" ? VALID :
    s == "D" ? INVALID :
    s == "U" ? INVALID :
    isnothing(s) ? UNKNOWN_REASON :
         throw(DomainError(s, "Unknown Invalid Reason Flag"))

@enum StandardConceptFlag UNKNOWN_STANDARD STANDARD NON_STANDARD CLASSIFICATION
StandardConceptFlag(::Nothing) = UNKNOWN_STANDARD
Base.parse(::Type{StandardConceptFlag}, s::Union{String, Nothing}) =
    s == "N" ? NON_STANDARD :
    s == "S" ? STANDARD :
    s == "C" ? CLASSIFICATION :
    isnothing(s) ? UNKNOWN_STANDARD :
         throw(DomainError(s, "Unknown Standard Concept Flag"))

@unpack struct Concept
    "CONCEPT_CLASS_ID" => concept_class_id::String = ""
    "CONCEPT_CODE" => concept_code::String
    "CONCEPT_ID" => concept_id::Int
    "CONCEPT_NAME" => concept_name::String
    "DOMAIN_ID" => domain_id::String
    "INVALID_REASON" => invalid_reason::InvalidReasonFlag = UNKNOWN_REASON
    "INVALID_REASON_CAPTION" => invalid_reason_caption::String
    "STANDARD_CONCEPT" => standard_concept::StandardConceptFlag = UNKNOWN_STANDARD
    "STANDARD_CONCEPT_CAPTION" => standard_concept_caption::String
    "VOCABULARY_ID" => vocabulary_id::String
end

abstract type Criteria end

function Base.getproperty(obj::Criteria, prop::Symbol)
    if prop in fieldnames(BaseCriteria)
        return getfield(obj.base, prop)
    else
        return getfield(obj, prop)
    end
end

@unpack struct Endpoint
    "Days" => days::Union{Int, Nothing} = nothing
    "Coeff" => coeff::Int
end

@unpack struct Window
    "Start" => start::Endpoint
    "End" => end_::Endpoint
    "UseIndexEnd" => use_index_end::Union{Bool, Nothing} = nothing
    "UseEventEnd" => use_event_end::Union{Bool, Nothing} = nothing
end

@enum OccurrenceType EXACTLY=0 AT_MOST=1 AT_LEAST=2
Base.parse(::Type{OccurrenceType}, s::String) =
    s == "0" ? EXACTLY :
    s == "1" ? AT_MOST :
    s == "2" ? AT_LEAST :
         throw(DomainError(s, "Unknown Occurrence Type"))

@unpack struct Occurrence
    "Type" => type::OccurrenceType
    "Count" => count::Int
    "IsDistinct" => is_distinct::Bool = false
    "CountColumn" => count_column::Union{String, Nothing} = nothing
end

@unpack struct CorrelatedCriteria
    "Criteria" => criteria::Union{Criteria, Nothing} = nothing
    "EndWindow" => end_window::Union{Window, Nothing} = nothing
    "IgnoreObservationPeriod" => ignore_observation_period::Bool = false
    "Occurrence" => occurrence::Union{Occurrence, Nothing} = nothing
    "RestrictVisit" => restrict_visit::Bool = false
    "StartWindow" => start_window::Window
end

@unpack struct DemographicCriteria
    "Age" => age::Union{NumericRange, Nothing} = nothing
    "Ethnicity" => ethnicity::Vector{Concept} = Concept[]
    "Gender" => gender::Vector{Concept} = Concept[]
    "OccurrenceEndDate" => occurrence_end_date::Union{DateRange, Nothing} = nothing
    "OccurrenceStartDate" => occurrence_start_date::Union{DateRange, Nothing} = nothing
    "Race" => race::Vector{Concept} = Concept[]
end

@enum CriteriaGroupType ALL_CRITERIA ANY_CRITERIA AT_LEAST_CRITERIA AT_MOST_CRITERIA
Base.parse(::Type{CriteriaGroupType}, s::String) =
    s == "ALL" ? ALL_CRITERIA :
    s == "ANY" ? ANY_CRITERIA :
    s == "AT_LEAST" ? AT_LEAST_CRITERIA :
    s == "AT_MOST" ? AT_MOST_CRITERIA :
    throw(DomainError(s, "Unknown Criteria Group Type"))

@unpack struct CriteriaGroup
    "Count" => count::Union{Int, Nothing} = nothing
    "CriteriaList" => correlated_criteria::Vector{CorrelatedCriteria} = CorrelatedCriteria[]
    "DemographicCriteriaList" => demographic_criteria::Vector{DemographicCriteria} = DemographicCriteria[]
    "Groups" => groups::Vector{CriteriaGroup} = CriteriaGroup[]
    "Type" => type::CriteriaGroupType
end

isempty(g::CriteriaGroup) =
    isempty(d.correlated_criteria) &&
    isempty(d.demographic_criteria) &&
    isempty(d.groups)

@enum CollapseType UNKNOWN_COLLAPSE ERA
CollapseType(::Nothing) = UNKNOWN_COLLAPSE
Base.parse(::Type{CollapseType}, s::Union{String, Nothing}) =
    s == "ERA" ? ERA :
    isnothing(s) ? UNKNOWN_COLLAPSE :
         throw(DomainError(s, "Unknown Collapse Type"))

@unpack struct CollapseSettings
    "CollapseType" => collapse_type::CollapseType
    "EraPad" => era_pad::Int = 0
end

@unpack struct Period
    "StartDate" => start_date::Union{Date, Nothing} = nothing
    "EndDate" => end_date::Union{Date, Nothing} = nothing
end

@unpack struct ConceptSetItem
    "concept" => concept::Concept
    "isExcluded" => is_excluded::Bool = false
    "includeDescendants" => include_descendants::Bool = false
    "includeMapped" => include_mapped::Bool = false
end

function unpack!(T::Type{Vector{ConceptSetItem}}, data::Dict)
    items = data["items"]
    retval = unpack!(T, items)
    if isempty(items)
        delete!(data, "items")
    end
    return retval
end

@unpack struct ConceptSet
    "id" => id::Int
    "name" => name::String
    "expression" => items::Vector{ConceptSetItem} = ConceptSetItem[]
end

abstract type EndStrategy end

@unpack struct CustomEraStrategy <: EndStrategy
    "DrugCodesetId" => drug_codeset_id::Union{Int, Nothing} = nothing
    "GapDays" => gap_days::Int = 0
    "Offset" => offset::Int = 0
    "DaysSupplyOverride" => days_supply_override::Union{Int, Nothing} = nothing
end

@enum DateField START_DATE END_DATE
Base.parse(::Type{DateField}, s::String) =
    s == "StartDate" ? START_DATE :
    s == "EndDate" ? END_DATE :
    throw(DomainError(s, "Unknown Date Field"))

@unpack struct DateOffsetStrategy <: EndStrategy
    "Offset" => offset::Integer
    "DateField" => date_field::DateField
end

function unpack!(::Type{EndStrategy}, data::Dict)
    if haskey(data, "DateOffset")
        (key, type) = ("DateOffset", DateOffsetStrategy)
    else
        (key, type) = ("CustomEra", CustomEraStrategy)
    end
    subdata = data[key]
    retval = unpack!(type, subdata)
    if isempty(subdata)
        delete!(data, key)
    end
    return retval
end

@unpack struct InclusionRule
    "name" => name::String
    "description" => description::String = ""
    "expression" => expression::CriteriaGroup
end

@unpack struct ObservationFilter
    "PriorDays" => prior_days::Int = 0
    "PostDays" => post_days::Int = 0
end

@enum ResultLimitType FIRST LAST ALL
Base.parse(::Type{ResultLimitType}, s::Union{String, Nothing}) =
    s == "First" ? FIRST :
    s == "Last" ? LAST :
    s == "All" ? ALL :
    isnothing(s) ? FIRST :
        throw(DomainError(s, "Unknown Result Limit Type"))

@unpack struct ResultLimit
    "Type" => type::ResultLimitType = FIRST
end

@unpack struct PrimaryCriteria
    "CriteriaList" => criteria_list::Vector{Criteria}
    "ObservationWindow" => observation_window::ObservationFilter
    "PrimaryCriteriaLimit" => primary_limit::ResultLimit
end

@unpack struct BaseCriteria
    "Age" => age::Union{NumericRange, Nothing} = nothing
    "CodesetId" => codeset_id::Union{Int, Nothing} = nothing
    "CorrelatedCriteria" => correlated_criteria::Union{CriteriaGroup, Nothing} = nothing
    "First" => first::Bool = false
    "Gender" => gender::Vector{Concept} = Concept[]
    "OccurrenceEndDate" => occurrence_end_date::Union{DateRange, Nothing} = nothing
    "OccurrenceStartDate" => occurrence_start_date::Union{DateRange, Nothing} = nothing
    "ProviderSpecialty" => provider_specialty::Vector{Concept} = Concept[]
    "VisitType" => visit_type::Vector{Concept} = Concept[]
end

struct UnknownCriteria <: Criteria
end

unpack!(::Type{UnknownCriteria}, data::Dict) = UnknownCriteria()

PrettyPrinting.quoteof(obj::UnknownCriteria) =
    Expr(:call, nameof(UnknownCriteria))

@unpack struct ConditionEra <: Criteria
    # like DrugEra, but missing gap_length?
    base::BaseCriteria
    "EraEndDate" => era_end_date::Union{DateRange, Nothing} = nothing
    "EraStartDate" => era_start_date::Union{DateRange, Nothing} = nothing
    "EraLength" => era_length::Union{NumericRange, Nothing} = nothing
    "OccurrenceCount" => occurrence_count::Union{NumericRange, Nothing} = nothing
    "AgeAtStart" => age_at_start::Union{NumericRange, Nothing} = nothing
    "AgeAtEnd" => age_at_end::Union{NumericRange, Nothing} = nothing
end

@unpack struct ConditionOccurrence <: Criteria
    base::BaseCriteria
    "ConditionSourceConcept" => condition_source_concept::Union{Int, Nothing} = nothing
    "ConditionStatus" => condition_status::Vector{Concept} = Concept[]
    "ConditionType" => condition_type::Vector{Concept} = Concept[]
    "ConditionTypeExclude" => condition_type_exclude::Bool = false
    "StopReason" => stop_reason::Union{TextFilter, Nothing} = nothing
end

@unpack struct Death <: Criteria
    base::BaseCriteria
    "DeathSourceConcept" => death_source_concept::Union{Int, Nothing} = nothing
    "DeathType" => death_type::Vector{Concept} = Concept[]
    "DeathTypeExclude" => death_type_exclude::Bool = false
end

@unpack struct DeviceExposure <: Criteria
    base::BaseCriteria
    "DeviceSourceConcept" => device_source_concept::Union{Int, Nothing} = nothing
    "DeviceType" => device_type::Vector{Concept} = Concept[]
    "DeviceTypeExclude" => device_type_exclude::Bool = false
    "Quantity" => quantity::Union{NumericRange, Nothing} = nothing
    "UniqueDeviceId" => unique_device_id::Union{TextFilter, Nothing} = nothing
end

@unpack struct DrugEra <: Criteria
    base::BaseCriteria
    "EraEndDate" => era_end_date::Union{DateRange, Nothing} = nothing
    "EraStartDate" => era_start_date::Union{DateRange, Nothing} = nothing
    "EraLength" => era_length::Union{NumericRange, Nothing} = nothing
    "OccurrenceCount" => occurrence_count::Union{NumericRange, Nothing} = nothing
    "GapDays" => gap_days::Union{NumericRange, Nothing} = nothing
    "AgeAtStart" => age_at_start::Union{NumericRange, Nothing} = nothing
    "AgeAtEnd" => age_at_end::Union{NumericRange, Nothing} = nothing
end

@unpack struct DrugExposure <: Criteria
    base::BaseCriteria
    "DrugSourceConcept" => drug_source_concept::Union{Int, Nothing} = nothing
    "DrugType" => drug_type::Vector{Concept} = Concept[]
    "DrugTypeExclude" => drug_type_exclude::Bool = false
    "Refills" => refills::Union{NumericRange, Nothing} = nothing
    "Quantity" => quantity::Union{NumericRange, Nothing} = nothing
    "DaysSupply" => days_supply::Union{NumericRange, Nothing} = nothing
    "RouteConcept" => route_concept::Vector{Concept} = Concept[]
    "EffectiveDrugDose" => effective_drug_dose::Union{NumericRange, Nothing} = nothing
    "DoseUnit" => dose_unit::Vector{Concept} = Concept[]
    "LotNumber" => lot_number::Union{TextFilter, Nothing} = nothing
    "StopReason" => stop_reason::Union{TextFilter, Nothing} = nothing
end

@unpack struct DoseEra <: Criteria
    base::BaseCriteria
    "DoseValue" => dose_value::Union{NumericRange, Nothing} = nothing
    "EraEndDate" => era_end_date::Union{DateRange, Nothing} = nothing
    "EraStartDate" => era_start_date::Union{DateRange, Nothing} = nothing
    "EraLength" => era_length::Union{NumericRange, Nothing} = nothing
    "AgeAtStart" => age_at_start::Union{NumericRange, Nothing} = nothing
    "AgeAtEnd" => age_at_end::Union{NumericRange, Nothing} = nothing
    "Unit" => unit::Vector{Concept} = Concept[]
end

@unpack struct LocationRegion <: Criteria
    "CodesetId" => codeset_id::Union{Int, Nothing} = nothing
    "StartDate" => start_date::Union{DateRange, Nothing} = nothing
    "EndDate" => end_date::Union{DateRange, Nothing} = nothing
end

Base.getproperty(obj::LocationRegion, prop::Symbol) =
    getfield(obj, prop)

@unpack struct Measurement <: Criteria
    base::BaseCriteria
    "MeasurementSourceConcept" => measurement_source_concept::Union{Int, Nothing} = nothing
    "MeasurementType" => measurement_type::Vector{Concept} = Concept[]
    "MeasurementTypeExclude" => measurement_type_exclude::Bool = false
    "Abnormal" => abnormal::Union{Bool, Nothing} = nothing
    "RangeLow" => range_low::Union{NumericRange, Nothing} = nothing
    "RangeHigh" => range_high::Union{NumericRange, Nothing} = nothing
    "RangeLowRatio" => range_low_ratio::Union{NumericRange, Nothing} = nothing
    "RangeHighRatio" => range_high_ratio::Union{NumericRange, Nothing} = nothing
    "ValueAsNumber" => value_as_number::Union{NumericRange, Nothing} = nothing
    "ValueAsConcept" => value_as_concept::Vector{Concept} = Concept[]
    "Operator" => operator::Vector{Concept} = Concept[]
    "Unit" => unit::Vector{Concept} = Concept[]
end

@unpack struct Observation <: Criteria
    base::BaseCriteria
    "ObservationSourceConcept" => observation_source_concept::Union{Int, Nothing} = nothing
    "ObservationType" => observation_type::Vector{Concept} = Concept[]
    "ObservationTypeExclude" => observation_type_exclude::Bool = false
    "ValueAsString" => value_as_string::Union{TextFilter, Nothing} = nothing
    "ValueAsNumber" => value_as_number::Union{NumericRange, Nothing} = nothing
    "ValueAsConcept" => value_as_concept::Vector{Concept} = Concept[]
    "Qualifier" => qualifier::Vector{Concept} = Concept[]
    "Unit" => unit::Vector{Concept} = Concept[]
end

@unpack struct ObservationPeriod <: Criteria
    base::BaseCriteria
    "PeriodType" => period_type::Vector{Concept} = Concept[]
    "PeriodTypeExclude" => period_type_exclude::Bool = false
    "PeriodStartDate" => period_start_date::Union{DateRange, Nothing} = nothing
    "PeriodEndDate" => period_end_date::Union{DateRange, Nothing} = nothing
    "PeriodLength" => period_length::Union{NumericRange, Nothing} = nothing
    "AgeAtStart" => age_at_start::Union{NumericRange, Nothing} = nothing
    "AgeAtEnd" => age_at_end::Union{NumericRange, Nothing} = nothing
    "UserDefinedPeriod" => user_defined_period::Union{Period, Nothing} = nothing
end

@unpack struct PayerPlanPeriod <: Criteria
    base::BaseCriteria
    "PeriodType" => period_type::Vector{Concept} = Concept[]
    "PeriodTypeExclude" => period_type_exclude::Bool = false
    "PeriodStartDate" => period_start_date::Union{DateRange, Nothing} = nothing
    "PeriodEndDate" => period_end_date::Union{DateRange, Nothing} = nothing
    "PeriodLength" => period_length::Union{NumericRange, Nothing} = nothing
    "AgeAtStart" => age_at_start::Union{NumericRange, Nothing} = nothing
    "AgeAtEnd" => age_at_end::Union{NumericRange, Nothing} = nothing
    "PayerConcept" => payer_concept::Union{Int, Nothing} = nothing
    "PlanConcept" => plan_concept::Union{Int, Nothing} = nothing
    "SponsorConcept" => sponsor_concept::Union{Int, Nothing} = nothing
    "StopReasonConcept" => stop_reason_concept::Union{Int, Nothing} = nothing
    "StopReasonSourceConcept" => stop_reason_source_concept::Union{Int, Nothing} = nothing
    "PayerSourceConcept" => payer_source_concept::Union{Int, Nothing} = nothing
    "PlanSourceConcept" => plan_source_concept::Union{Int, Nothing} = nothing
    "SponsorSourceConcept" => sponsor_source_concept::Union{Int, Nothing} = nothing
    "UserDefinedPeriod" => user_defined_period::Union{Period, Nothing} = nothing
end

@unpack struct ProcedureOccurrence <: Criteria
    base::BaseCriteria
    "ProcedureSourceConcept" => procedure_source_concept::Union{Int, Nothing} = nothing
    "ProcedureType" => procedure_type::Vector{Concept} = Concept[]
    "ProcedureTypeExclude" => procedure_type_exclude::Bool = false
    "Modifier" => modifier::Vector{Concept} = Concept[]
    "Quantity" => quantity::Union{NumericRange, Nothing} = nothing
end

@unpack struct Specimen <: Criteria
    base::BaseCriteria
    "SpecimenSourceConcept" => specimen_source_concept::Union{Int, Nothing} = nothing
    "SpecimenType" => specimen_type::Vector{Concept} = Concept[]
    "SpecimenTypeExclude" => specimen_type_exclude::Bool = false
    "Quantity" => quantity::Union{NumericRange, Nothing} = nothing
    "Unit" => unit::Vector{Concept} = Concept[]
    "AnatomicSite" => anatomic_site::Vector{Concept} = Concept[]
    "DiseaseStatus" => disease_status::Vector{Concept} = Concept[]
    "SourceId" => source_id::Union{TextFilter, Nothing} = nothing
end

@unpack struct VisitOccurrence <: Criteria
    base::BaseCriteria
    "PlaceOfService" => place_of_service::Vector{Concept} = Concept[]
    "PlaceOfServiceLocation" => place_of_service_location::Union{Int, Nothing} = nothing
    "VisitSourceConcept" => visit_source_concept::Union{Int, Nothing} = nothing
    "VisitLength" => visit_length::Union{NumericRange, Nothing} = nothing
    "VisitTypeExclude" => visit_type_exclude::Bool = false
end

function unpack!(::Type{Criteria}, data::Dict)
    for type in (ConditionEra, ConditionOccurrence, Death,
                 DeviceExposure, DoseEra, DrugEra, DrugExposure,
                 LocationRegion, Measurement, Observation,
                 ObservationPeriod, PayerPlanPeriod,
                 ProcedureOccurrence, Specimen, VisitOccurrence)
        key = string(nameof(type))
        if haskey(data, key)
            subdata = data[key]
            retval = unpack!(type, subdata)
            if isempty(subdata)
                delete!(data, key)
            end
            return retval
        end
    end
    return unpack!(UnknownCriteria, data)
end

@unpack struct CohortExpression
    "AdditionalCriteria" => additional_criteria::Union{CriteriaGroup, Nothing} = nothing
    "CensorWindow" => censor_window::Union{Period, Nothing} = nothing
    "CensoringCriteria" => censoring_criteria::Vector{Criteria} = Criteria[]
    "CollapseSettings" => collapse_settings::CollapseSettings
    "ConceptSets" => concept_sets::Vector{ConceptSet} = ConceptSet[]
    "EndStrategy" => end_strategy::Union{EndStrategy, Nothing} = nothing
    "ExpressionLimit" => expression_limit::ResultLimit
    "InclusionRules" => inclusion_rules::Vector{InclusionRule} = InclusionRule[]
    "PrimaryCriteria" => primary_criteria::PrimaryCriteria
    "QualifiedLimit" => qualified_limit::ResultLimit
    "Title" => title::Union{String, Nothing} = nothing
    "cdmVersionRange" => version_range::Union{String, Nothing} = nothing
end

unpack!(data) = unpack!(CohortExpression, data)


text:
module Expression

using CSV
using DataFrames

function getExpression(filename)
    headerCount = open(filename) do file
        lineCounter = 1
        for line in eachline(file)
             if string(line[1]) == "#"
                lineCounter += 1
             else
                break
             end
        end

        lineCounter
    end

    return CSV.read(filename, delim="\t", header=headerCount)
end

function getTissue(df, tissue)
    # Where the Gene Name is the HUGO name
    tissueValues = Dict()
    for row in eachrow(df)
        value = row[tissue]
        if isnothing(value) == false
            tissueValues[row[Symbol("Gene Name")]] = value
        end
    end

    return tissueValues
end

end

max_stars_repo_name: OICR/mp-biopath

max_issues_repo_name: OICR/mp-biopath

max_forks_repo_name: OICR/mp-biopath


text:
import JSExpr: @js, @js_str, JSString, jsstring, jsexpr, @var, @new
using Lazy: @d

export js, @js, @js_, @var, @new

# include("jsexprs.jl")

include("callbacks.jl")

mutable struct JSError <: Exception
    name::String
    msg::String
end

Base.showerror(io::IO, e::JSError) =
    print(io, "Javascript error\t$(e.name): $(e.msg)")

# RPC API

export js, js_, @js, @js_, @var, @new

"""
    js(win, expr::JSString; callback=false)

Execute the javscript in `expr`, inside `win`.

If `callback==true`, returns the result of evaluating `expr`.
"""
function js end

"""
    JSString(str)

A wrapper around a string indicating the string contains javascript code.
"""
function JSString end

msg(o, m) = error("$(typeof(o)) object doesn't support JS messages")

function js(o, js::JSString; callback = true)
  cmd = @d(:type => :eval,
           :code => js.s)
  if callback
    id, cond = callback!()
    cmd[:callback] = id
  end
  msg(o, cmd)

  if callback
      val = wait(cond)
      if isa(val, AbstractDict) && get(val, "type", "") == "error"
          err = JSError(get(val, "name", "unknown"), get(val, "message", "blank"))
          throw(err)
      end
      return val
  else
      return o
  end
end

js(o, j; callback=true) = js(o, JSString(string(jsstring(j)...)); callback=callback)

"""
    @js win expr

Execute `expr`, converted to javascript, inside `win`, and return the result.

`expr` will be parsed as julia code, and then converted directly to the
equivalent javascript. Language keywords that don't exist in julia can be
represented with their macro equivalents, `@var`, `@new`, etc.

See also: `@js_`, the asynchronous version.

# Examples
```julia-repl
julia> @js win x = 5
5
julia> @js_ win for i in 1:x console.log(i) end
```
"""
macro js(o, ex)
    :(js($(esc(o)), $(Expr(:quote, ex))))
end

"""
    @js_ win expr

Execute `expr`, converted to javascript, asynchronously inside `win`, and return
immediately.

`expr` will be parsed as julia code, and then converted directly to the
equivalent javascript. Language keywords that don't exist in julia can be
represented with their macro equivalents, `@var`, `@new`, etc.

See also: `@js`, the synchronous version that returns its result.

# Examples
```julia-repl
julia> @js win x = 5
5
julia> @js_ win for i in 1:x console.log(i) end
```
"""
macro js_(o, ex)
    :(js($(esc(o)), $(Expr(:quote, ex)), callback=false))
end


text:
using PkgSRA
using Plots, Random, OrdinaryDiffEq, Statistics
using Distributions, Turing
using StatsPlots, DSP
using Flux
using BSON: @save
pyplot()
Random.seed!(11)

#####
##### Define the controlled dynamical system
#####
# Load helper functions
# include("../../utils/sindy_statistics_utils.jl")
include("../../utils/sindy_turing_utils.jl")
include("../scripts/paper_settings.jl")
# Load example problem
include(EXAMPLE_FOLDERNAME*"example_lorenz.jl")

# Define the multivariate forcing function
num_ctr = 50;
    U_starts = rand(3, num_ctr) .* tspan[2]
    U_widths = [0.1, 0.1, 0.0];
    amplitudes = [150.0, 100.0, 0.0]
my_U_func_time2(t, u) = U_func_time_multivariate(t, u,
                        U_widths, U_starts,
                        F_dim=[1, 2, 3],
                        amplitudes=amplitudes)

#####
##### Produce data
#####
sol = solve_lorenz_system(my_U_func_time2)
dat = Array(sol)
# plot(dat[1,:],dat[2,:],dat[3,:],label="Data")
# Derivatives
numerical_grad = numerical_derivative(dat, ts)

## Also get baseline true/ideal cases
# Uncontrolled
dat_raw = Array(solve_lorenz_system())
numerical_grad_raw = numerical_derivative(dat_raw, ts)
true_grad_raw = zeros(size(dat))
for i in 1:size(dat,2)
    true_grad_raw[:,i] = lorenz_system(dat_raw[:,i], Float64.(p), [0])
end

# Controlled; true variables
true_grad = zeros(size(dat))
for i in 1:size(dat,2)
    true_grad[:,i] = lorenz_system(true_grad[:,i], dat[:,i], p, [0])
end
let d = dat
    plot(d[1,:],d[2,:],d[3,:],label="Data")
end

# True control signal
U_true = zeros(size(dat))
for (i, t) in enumerate(ts)
    U_true[:,i] = my_U_func_time2(t, dat[:,i])
end
plot(U_true')

#####
##### Calculate distribution of residuals
#####
# Turing.setadbackend(:forward_diff)
# CHEATING: known structure
#   Calculate posterior distribution of parameters
# chain, train_ind = generate_chain(dat,
#                         numerical_grad, lorenz_grad_residual)
sindy_library = Dict("cross_terms"=>2,"constant"=>nothing);
Random.seed!(13)

# Upgrade: Use AIC
# val_list = calc_permutations(5,3)
# val_list = combinations(1:5, 3)
val_list = Iterators.product(1:3,1:3,1:3) # DEBUG
(sindy_unctr,best_criterion,all_criteria,all_models) =
    sindyc_ensemble(dat, numerical_grad, sindy_library, val_list,
                    selection_criterion=my_aicc,
                    sparsification_mode="num_terms",
                    selection_dist=Normal(0.0,20),
                    use_clustering_minimization=true)
print_equations(sindy_unctr)
sindy_grad1 = sindy_unctr(dat, 0)
# scatter(sum.(val_list), all_criteria)

# Generate a full trajectory from SINDy model
condition(u,t,integrator) = any(abs.(u).>2e2)
cb = DiscreteCallback(condition, terminate!)
prob_unctr = ODEProblem(sindy_unctr, u0, tspan, [0], callback=cb)
sindy_dat_unctr = Array(solve(prob_unctr, Tsit5(), saveat=ts));

plot(dat[1,:],dat[2,:],dat[3,:],label="Data")
    plot!(sindy_dat_unctr[1,:],sindy_dat_unctr[2,:],sindy_dat_unctr[3,:],label="Turing model with enforced zeros")
    title!("Integrated SINDy model")

# grad_unctr = sindy_unctr(dat)
# let d = grad_unctr, d0 = numerical_grad
#     plot(d[1,:],d[2,:],d[3,:],label="Model")
#     plot!(d0[1,:],d0[2,:],d0[3,:],label="Data")
#     title!("Gradients (L2 SINDy model)")
# end
turing_unctr = convert_sindy_to_turing_enforce_zeros(sindy_unctr;
                                dat_noise_prior=Normal(0.0, 20.0),
                                coef_noise_std=2.0)
chain = generate_chain(dat, numerical_grad, turing_unctr,
                            iterations=200,
                            num_training_pts=100, start_ind=101)[1]
turing_unctr_sample = sindy_from_chain(sindy_unctr, chain,
                                        enforced_zeros=true)
# plot(chain["all_coef[3]"])
# Generate a single full sample trajectory (FROM TURING)
# prob_unctr = ODEProblem(sindy_unctr, u0, tspan, [0], callback=cb)
prob_unctr = ODEProblem(turing_unctr_sample, u0, tspan, [0], callback=cb)
sol_unctr = solve(prob_unctr, Tsit5(), saveat=ts);
dat_unctr = Array(sol_unctr)

# plot(dat[1,:],dat[2,:],dat[3,:],label="Data")
plot(dat_raw[1,:],dat_raw[2,:],dat_raw[3,:],label="Data")
    plot!(dat_unctr[1,:],dat_unctr[2,:],dat_unctr[3,:],label="Turing model with enforced zeros")
    title!("Integrated Turing model")

# Get the full posterior distribution
sample_ind = 1:length(ts)
sample_gradients3d, sample_noise = sample_sindy_posterior_grad(chain,
            dat, sample_ind, sindy_unctr)
# let y = sample_gradients3d[1,:,:], y2 = sample_gradients3d[2,:,:], n=numerical_grad
#     plot3d(y[:,1], y[:,2], y[:,3], label="Naive Sample 1")
#     plot3d!(y2[:,1], y2[:,2], y2[:,3], label="Naive Sample 2")
#     plot3d!(n[1,:], n[2,:], n[3,:], label="Data")
#     title!("Sampled Turing Gradients")
# end

# Calculate the residuals
sample_gradients = transpose(drop_all_1dims(mean(sample_gradients3d, dims=1)))
    sample_gradients_std = transpose(drop_all_1dims(std(sample_gradients3d, dims=1)))
    sample_trajectory_noise = mean(sample_noise)
dat_grad = numerical_grad[:, sample_ind]

residual = dat_grad .- sample_gradients
ctr_guess = process_residual(residual, sample_trajectory_noise)

# Visualization
ind = 101:1000
# ind = sample_ind
i = 2
    plot(dat_grad[i,ind], label="Data gradient")
    plot!(sample_gradients[i,ind], label="Sample gradient")
i = 2
    plot(residual[i,ind], ribbon=sample_trajectory_noise)
    plot!(U_true[i,ind], label="True")
    title!("Residual and true control")
i = 1
    plot(sample_gradients_std[i,ind])
    plot!(U_true[i,ind], label="True")
    title!("Residual and true control")
i = 2
    plot(residual[i,ind],ribbon=
            sample_trajectory_noise.+sample_gradients_std[i,ind])
    plot!(U_true[i,ind], label="True")
    title!("Residual and true control (with std)")
# i = 1
#     plot(ctr_guess[i,ind], label="Control guess")
#     plot!(U_true[i,ind], label="True", show=true)

#####
##### INSTEAD: subsample the data based on the above residual
#####

# Get "confident" indices
accepted_ind = subsample_using_residual(residual,
            sample_trajectory_noise, min_length=4)

num_pts = 400
start_ind = 101
subsample_ind = accepted_ind[start_ind:num_pts+start_ind-1]
dat_sub = dat[:,subsample_ind]
grad_sub = numerical_grad[:,subsample_ind]

# Any control signal in the subset?
U_sub = U_true[:,subsample_ind]
plot(U_sub')
    # plot!(residual[:,subsample_ind]', ribbon=sample_trajectory_noise)
    title!("Control signals in the subsampled dataset")

# SINDY SETUP
# val_list = calc_permutations(5,3)
# val_list = combinations(1:5, 3)
val_list = Iterators.product(1:3,1:3,1:3) # DEBUG
(sindy_sub,best_criterion,all_criteria,all_models,best_index) =
    sindyc_ensemble(dat_sub, grad_sub, sindy_library, val_list,
                    selection_criterion=my_aicc,
                    sparsification_mode="num_terms",
                    selection_dist=Normal(0.0,sample_trajectory_noise),
                    use_clustering_minimization=true)
println("Best index is $best_index")
print_equations(sindy_sub)
scatter(sum.(val_list), all_criteria)
    title!("AIC for various sparsities")
    xlabel!("Number of nonzero terms")

# Generate a single trajectory
prob_ctr = ODEProblem(sindy_sub, u0, tspan, [0], callback=cb)
sindy_dat_ctr = Array(solve(prob_ctr, Tsit5(), saveat=ts));
sindy_grad_ctr = sindy_sub(dat, 0)

plot(dat[1,:],dat[2,:],dat[3,:],label="Data")
    plot!(sindy_dat_ctr[1,:],sindy_dat_ctr[2,:],sindy_dat_ctr[3,:],label="Turing model with enforced zeros")
    title!("Integrated SINDy model (after control)")

i = 1
    plot(dat_grad[i,ind], label="Data Gradient")
    plot!(sindy_grad1[i,ind], label="Old SINDy gradient")
    plot!(sindy_grad_ctr[i,ind], label="New SINDy Gradient")

# Just get the residual from the SINDy model, not Turing
residual_ctr = dat_grad .- sindy_grad_ctr

i = 2
    plot(residual[i,ind], ribbon=sample_trajectory_noise, label="Original residual")
    plot!(residual_ctr[i,ind], label="Residual2")
    plot!(U_true[i,ind], label="True")
    title!("Iteration 2: Residual and true control for i=$i")

# TURING ANALYSIS
turing_sub = convert_sindy_to_turing_enforce_zeros(sindy_sub;
                                dat_noise_prior=Normal(0.0, 1.0),
                                coef_noise_std=0.1)
chain_sub = generate_chain(dat, numerical_grad, turing_sub,
                            train_ind=subsample_ind,
                            iterations=300)[1]
turing_sub_sample = sindy_from_chain(sindy_sub, chain_sub,
                                        enforced_zeros=true)
println("Turing sample equations")
print_equations(turing_sub_sample)
println("SINDy equations")
print_equations(sindy_sub)
# Generate test trajectories from the CORRECTED posterior
sample_gradients3d_ctr, sample_noise_ctr =
            sample_sindy_posterior_grad(chain_sub,
                    dat, sample_ind, sindy_sub)

density(chain_sub)
# Calculate the residuals
sample_gradients_ctr = transpose(drop_all_1dims(mean(sample_gradients3d_ctr, dims=1)))
    sample_trajectory_noise_ctr = mean(sample_noise_ctr)
# dat_grad2 = numerical_grad[:, sample_ind]

residual2 = dat_grad .- sample_gradients_ctr
ctr_guess2 = process_residual(residual2, sample_trajectory_noise_ctr)

# Generate a single full sample trajectory (FROM TURING)
prob_ctr = ODEProblem(turing_sub_sample, u0, tspan, [0], callback=cb)
turing_dat_ctr = Array(solve(prob_ctr, Tsit5(), saveat=ts));

# Plot!
plot(dat[1,:],dat[2,:],dat[3,:],label="Data")
    plot!(dat_ctr[1,:],dat_ctr[2,:],dat_ctr[3,:],label="Turing model with enforced zeros")
    title!("Integrated Turing model (after control)")
i = 2
    plot(residual[i,ind], ribbon=sample_trajectory_noise, label="Original residual")
    plot!(residual2[i,ind], ribbon=sample_trajectory_noise_ctr, label="Residual2")
    plot!(U_true[i,ind], label="True")
    title!("Iteration 2: Residual and true control for i=$i")
i = 2
    plot(dat_grad[i,ind], label="Data gradient")
    plot!(sample_gradients_ctr[i,ind], label="Sample gradient2")
    plot!(sample_gradients[i,ind], label="Sample gradient1")
    title!("Gradient predictions")

println("SINDy of subset")
print_equations(sindy_sub)
println("Turing-SINDy of subset")
print_equations(turing_sub_sample)
println("True")
print_equations(core_dyn_true)

#####
##### Save everything; plotted in a different script
#####
this_dat_name = DAT_FOLDERNAME*"TMP_dat_flowchart_inkscape_"

# Intermediate data
fname = this_dat_name*"intermediate.bson";
@save fname chain chain_sub sol sample_gradients sample_gradients_ctr sample_noise sample_noise_ctr sample_ind numerical_grad

# Uncontrolled data
fname = this_dat_name*"uncontrolled_lorenz.bson"
@save fname dat_raw numerical_grad_raw

# Original ODE and controller variables
fname = this_dat_name*"ode_vars.bson";
@save fname dat dat_grad true_grad U_true

# SINDy variables 1: before control
fname = this_dat_name*"naive_vars_sindy.bson";
unctr_nz_terms = get_nonzero_terms(sindy_unctr)
unctr_nz_names = get_nonzero_term_names(sindy_unctr)
# NOTE: saving the SINDy model directly currently breaks BSON
@save fname sindy_dat_unctr unctr_nz_terms unctr_nz_names #sindy_unctr

# Bayesian variables 1: before control
fname = this_dat_name*"naive_vars_bayes.bson";
@save fname dat_unctr sample_trajectory_noise residual accepted_ind ctr_guess

# SINDY variables 2: after control
fname = this_dat_name*"ctr_vars_sindy.bson";
sub_nz_terms = get_nonzero_terms(sindy_sub)
sub_nz_names = get_nonzero_term_names(sindy_sub)
@save fname sindy_dat_ctr sindy_grad_ctr sub_nz_terms sub_nz_names #sindy_sub

# Bayesian variables 2: after control
fname = this_dat_name*"ctr_vars_turing.bson";
@save fname dat_ctr sample_trajectory_noise_ctr sample_gradients_ctr turing_dat_ctr ctr_guess2

# TODO: NN variables


text:
function dtsol(n, L::TronSparseMatrixCSC, r)
    zero = 0.0

    # Solve L'*x = r and store the result in r.

    r[n] = r[n]/L.diag_vals[n]
    @inbounds for j=n-1:-1:1
        temp = zero
        for k=L.colptr[j]:L.colptr[j+1]-1
            temp = temp + L.tril_vals[k]*r[L.rowval[k]]
        end
        r[j] = (r[j] - temp)/L.diag_vals[j]
    end

    return
end

function dtsol(n, L::TronDenseMatrix, r)
    zero = 0.0

    # Solve L'*x = r and store the result in r.

    r[n] = r[n]/L.vals[n,n]
    @inbounds for j=n-1:-1:1
        temp = zero
        for k=j+1:n
            temp = temp + L.vals[k,j]*r[k]
        end
        r[j] = (r[j] - temp)/L.vals[j,j]
    end

    return
end

@inline function dtsol(n::Int, L::CuDeviceArray{Float64,2},
                       r::CuDeviceArray{Float64,1})
    # Solve L'*x = r and store the result in r.

    tx = threadIdx().x
    ty = threadIdx().y

    @inbounds for j=n:-1:1
        if tx == 1 && ty == 1
            r[j] = r[j] / L[j,j]
        end
        CUDA.sync_threads()

        if tx < j && ty == 1
            r[tx] = r[tx] - L[tx,j]*r[j]
        end
        CUDA.sync_threads()
    end

    return
end


text:
"""
    circshift(v::AbstractArray, n)

Return a `CircShiftedArray` object, with underlying data `v`. The second argument gives the amount
to circularly shift in each dimension. If it is an integer, it is assumed to refer to the
first dimension.

# Examples

```jldoctest circshift
julia> v = [1, 3, 5, 4];

julia> ShiftedArrays.circshift(v, 1)
4-element ShiftedArrays.CircShiftedArray{Int64,1,Array{Int64,1}}:
 4
 1
 3
 5

julia> w = reshape(1:16, 4, 4);

julia> ShiftedArrays.circshift(w, (1, -1))
4×4 ShiftedArrays.CircShiftedArray{Int64,2,Base.ReshapedArray{Int64,2,UnitRange{Int64},Tuple{}}}:
 8  12  16  4
 5   9  13  1
 6  10  14  2
 7  11  15  3
```
"""
circshift(v::AbstractArray, n) = CircShiftedArray(v, n)


text:
"""
    covcor_axes(x, dim) -> NTuple{2}

Returns appropriate axes for a `cov` or `var` method on array `x`.

## Examples
```jldoctest
julia> covcor_axes((Index{:a}(1:4), Index{:b}(1:4)), 2)
(Index{a}(1:4 => Base.OneTo(4)), Index{a}(1:4 => Base.OneTo(4)))

julia> covcor_axes((Index{:a}(1:4), Index{:b}(1:4)), :b)
(Index{a}(1:4 => Base.OneTo(4)), Index{a}(1:4 => Base.OneTo(4)))

julia> covcor_axes((Index{:a}(1:4), Index{:b}(1:4)), 1)
(Index{b}(1:4 => Base.OneTo(4)), Index{b}(1:4 => Base.OneTo(4)))

julia> covcor_axes((Index{:a}(1:4), Index{:b}(1:4)), :a)
(Index{b}(1:4 => Base.OneTo(4)), Index{b}(1:4 => Base.OneTo(4)))
```
"""
covcor_axes(x::AbstractMatrix, dim) = _covcor_axes(axes(x), to_dims(x, dim))
covcor_axes(x::NTuple{2,Any}, dim) = _covcor_axes(x, to_dims(x, dim))
_covcor_axes(x::NTuple{2,Any}, dim::Int) = dim === 1 ? (x[2], x[2]) : (x[1], x[1])


text:
"""
    setinvalid_qualityflag!(df;
      vars=["LE","H","NEE","Tair","VPD","wind"],
      qc_suffix="_qc",
      good_quality_threshold = 1.0,
      missing_qc_as_bad = true,
      setvalmissing = true,
    )

Set records with quality flags indicating problems to false in :valid column.

# Arguments
- df: DataFrame with column :GPP
optional
- `vars=["LE","H","NEE","Tair","VPD","wind"]`: columns to theck for quality
- `qc_suffix="_qc"`: naming of the corresponding quality-flag column
- `good_quality_threshold = 1.0`: threshold in quality flag up to which
  data is considered good quality
- `missing_qc_as_bad = true`: set to false to not mark records with missing
  quality flag as invalid
- `setvalmissing = true`: set to false to prevent replacing values in value column
  corresponding to problematic quality flag to missing.

# Value
df with modified :valid and value columns.

# Example
```jldoctest; output = false
using DataFrames
df = DataFrame(
  NEE = 1:3, NEE_qc = [1,1,2],
  GPP = 10:10:30, GPP_qc = [1,missing,1])
setinvalid_qualityflag!(df; vars = ["NEE", "GPP"])
df.valid == [true, false, false]
ismissing(df.GPP[2]) && ismissing(df.NEE[3])
# output
true
```
"""
function setinvalid_qualityflag!(df; setvalmissing = true, kwargs...)
  if setvalmissing
    set_badquality_missing!(df; kwargs...)
  end
  _setinvalid_qualityflag!(df; kwargs...)
end

function set_badquality_missing!(df;
  vars=["LE","H","NEE","Tair","VPD","wind"],
  qc_suffix="_qc",
  good_quality_threshold = 1.0,
  missing_qc_as_bad = true
  )
  fqc(x,xqc) = if missing_qc_as_bad
     ifelse(!ismissing(xqc) && xqc <= good_quality_threshold, x, missing)
  else
    ifelse(ismissing(xqc) || xqc <= good_quality_threshold, x, missing)
  end
  tmp = map(vars) do var
      qcvar = var * qc_suffix
      [var, qcvar] => ByRow(fqc) => var
  end
  transform!(df, tmp...)
end

function _setinvalid_qualityflag!(df;
  vars=["LE","H","NEE","Tair","VPD","wind"],
  qc_suffix="_qc",
  good_quality_threshold = 1.0,
  missing_qc_as_bad=true
  )
  vars_qc = vars .* qc_suffix
  fvalid_var = missing_qc_as_bad ? valid_nonmissingvar_ : valid_missingvar_
  function fsel(valid, x...)
     nvar = length(x) ÷ 2
     vs = x[1:nvar]
     vs_qc = x[(nvar+1):end]
     #valid = repeat(@MVector([true]), length(vs[1]))
     for i = 1:nvar
         # TODO change when 1.7 is out which support .&&: @. valid = valid && fvalid_var(vs[i], vs_qc[i], good_quality_threshold)
         @. valid = valid & fvalid_var(vs[i], vs_qc[i], good_quality_threshold)
     end
     valid
  end
  if !hasproperty(df, :valid) df[!,:valid] .= true; end
  select!(df, :, Cols(vcat("valid", vars, vars_qc)) => fsel => :valid)
end

function valid_nonmissingvar_(x, xqc, good_quality_threshold)
  !ismissing(x) && !ismissing(xqc) && xqc <= good_quality_threshold
end
function valid_missingvar_(x, xqc, good_quality_threshold)
 !ismissing(x) && (ismissing(xqc) || xqc <= good_quality_threshold)
end

"""
    setinvalid_range!(df, var_ranges...; setvalmissing = true, ...)

Set records with values outside specified ranges to false in :valid column.

If their is no limit towards small or
large values, supply `-Inf` or `Inf` as the minimum or maximum respectively.
If there were false values in the :value column before, they are kept false.
In addition, set values outside ranges to missing.

# Arguments
- `df`: DataFrame with column :GPP
- `var_ranges`: Pair `Varname_symbol => (min,max)`: closed valid interval for
  respective column
optional
- setvalmissing: set to false to prevent replacing values in value column outside ranges to missing.

# Value
df with modified :valid and value columns.
```jldoctest; output = false
using DataFrames
df = DataFrame(NEE = [missing, 2,3], GPP = 10:10:30)
setinvalid_range!(df, :NEE => (-2.0,4.0), :GPP => (8.0,28.0))
df.valid == [false, true, false]
ismissing(df.GPP[3])
# output
true
```
"""
function setinvalid_range!(df, var_ranges::Vararg{Pair,N}; setvalmissing = true, kwargs...) where N
  if setvalmissing
    set_badrange_missing!(df, var_ranges...; kwargs...)
  end
  _setinvalid_range!(df, var_ranges...; kwargs...)
end

function set_badrange_missing!(df, var_ranges::Vararg{Pair,N}) where N
  tmp = map(var_ranges) do p
      var, (min, max) = p
      var => ByRow(x -> (ifelse(!ismissing(x) && min <= x <= max, x, missing))) => var
  end
  transform!(df, tmp...)
end

function _setinvalid_range!(df, var_ranges::Vararg{Pair,N}) where N
  function fval(valid, x...)
    for (p,xj) in zip(var_ranges, x)
      var, (min, max) = p
      #TODO simplify once 1.7 with .&& is out: @. valid = valid && !ismissing(xi) && (min <= xi <= max)
      for i in eachindex(xj)
        valid[i] = valid[i] && !ismissing(xj[i]) && (min <= xj[i] <= max)
      end
    end
    valid
  end
  vars = map(x -> x.first, var_ranges)
  if !hasproperty(df, :valid) df[!,:valid] .= true; end
  select!(df, :, SA[:valid, vars...] => fval => :valid)
end

"""
    setinvalid_nongrowingseason!(df, tGPP; kwargs...)

Set non-growseason to false in :valid column.

# Arguments
- `df`: DataFrame with columns `:GPP` and `:datetime`
- `tGPP`: scalar threshold of daily GPP (see [`get_growingseason`](@ref))
optional:
- `update_GPPd`: set to true additionally update `:GPPd_smoothed` column to
  results from [`get_growingseason`](@ref)
- and others passed to [`get_growingseason`](@ref)

# Value
df with modified columns :valid and if  `:GPPd_smoothed`,
where all non-growing season records are set to false.
"""
function setinvalid_nongrowingseason!(df, tGPP; update_GPPd_smoothed = false, kwargs...)
  if !hasproperty(df, :valid) df[!,:valid] .= true; end
  # non-copying dataframe where we can add grouping column __day
  dft = transform(df, :datetime => ByRow(Date) => :__day, copycols = false)
  function mean_nonmissing(x)
      mx = mean(skipmissing(x))
      # if there is no non-missing record in a group, mean returns nan
      ifelse(isnan(mx), missing, mx)
  end
  dfday = combine(groupby(dft, :__day), :GPP => mean_nonmissing => :GPPd, nrow)
  transform!(dfday,
    :GPPd => (x -> get_growingseason(x, tGPP; kwargs...)) => SA[:valid, :GPPd_smoothed])
  # rep(dfd.valid, each = df.nrow)
  #TODO simplify once 1.7 is out: df[!,:valid] .= df.valid .&& vcat(fill.(dfday.valid, dfday.nrow)...)
  df[!,:valid] .= df.valid .& vcat(fill.(dfday.valid, dfday.nrow)...)
  if update_GPPd_smoothed
      df[!,:GPPd_smoothed] .= vcat(fill.(dfday.GPPd_smoothed, dfday.nrow)...)
  end
  df
end

"""
    get_growingseason(GPPd, tGPP; ws=15, min_int=5, warngap=true)

Filters annual time series for growing season based on smoothed daily GPP data.

# Arguments
- `GPPd`:    daily GPP (any unit)
- `tGPP`:    GPP threshold (fraction of 95th percentile of the GPP time series).
               Takes values between 0 and 1.
optional
- `ws`:      window size used for GPP time series smoothing
- `min_int`: minimum time interval in days for a given state of growing season
- `warngap`: set to false to suppress warning on too few non-missing data

# Details
The basic idea behind the growing season filter is that vegetation is
considered to be active when its carbon uptake (GPP) is above a specified
threshold, which is defined relative to the peak GPP (95th percentile)
observed in the year.
The GPP-threshold is calculated as:

``GPP_{threshold} = quantile(GPPd,0.95)*tGPP``

GPPd time series are smoothed with a moving average to avoid fluctuations
in the delineation of the growing season. The window size defaults to 15
days, but depending on the ecosystem, other values can be appropriate.

The argument `min_int` serves to avoid short fluctuations in the
status growing season vs. no growing season by defining a minimum length
of the status. If a time interval shorter than `min_int` is labeled
as growing season or non-growing season, it is changed to the status of
the neighboring values, i.e its opposite.

# Value
A NamedTuple with entries
- `is_growingseason`: a `BitVector` of the same length as the input GPPd in which `false`
  indicate no growing season (dormant season) and `true` indicate growing season.
- `GPPd_smoothed`: smoothed GPPd
"""
function get_growingseason(GPPd,tGPP;ws=15,min_int=5,warngap=true)
  nday = length(GPPd)
  if sum(.!ismissing.(GPPd)) < 0.5*nday
    error("Expected number of available GPPd data " *
    "to be at least half the total number of days ($nday). " *
    "But was only ($(sum(.!ismissing.(GPPd)))).")
  end
  GPP_threshold = quantile(skipmissing(GPPd), 0.95)*tGPP
  #@show GPP_threshold
  # smooth GPP
  #fromR: GPPd_smoothed = filter(GPPd,method="convolution",filter=rep(1/ws,ws))
  GPPd_smoothed = moving_average(GPPd, ws)
  # set values at the beginning and end of the time series to the mean of original values
  wsd = floor(Int, ws/2)
  GPPd_smoothed[1:wsd] .= mean(skipmissing(GPPd[1:(2*wsd)]))
  GPPd_smoothed[(length(GPPd)-(wsd-1)):length(GPPd)] .=
    mean(skipmissing(GPPd[(length(GPPd)-(2*wsd-1)):length(GPPd)]))
  # check for occurence of missing values and set to mean of the values surrounding them
  imissing = findall(ismissing, GPPd_smoothed)
  if length(imissing) > 0
    warngap && length(imissing) > 10 &&
      @warn "Attention, there is a gap in 'GPPd' of length n = $(length(imissing))"
    #TODO check and correct for several gaps, see Impute package
    replace_val = mean(skipmissing(GPPd_smoothed[max(1,imissing[1] - 4):min(
      (imissing[length(imissing)] + 4),length(GPPd_smoothed))]))
    GPPd_smoothed = coalesce.(GPPd_smoothed, replace_val)
  end
  # filter daily GPP
  growseas = GPPd_smoothed .>= GPP_threshold
  # change short intervals to the surrounding values to avoid 'wrong' fluctuations
  # switch shortest interval successively and recompute interval lengths
  intervals = rle(growseas)[2]
  imin = argmin(intervals)
  while intervals[imin] < min_int
      end_int = cumsum(intervals)
      start_int = vcat(0, end_int[1:end-1]) .+ 1
      growseas[start_int[imin]:end_int[imin]] .= !growseas[start_int[imin]]
      intervals = rle(growseas)[2]
      imin = argmin(intervals)
  end
  (is_growingseason = growseas, GPPd_smoothed = GPPd_smoothed)
end

"""
    setinvalid_afterprecip!(df; min_precip = 0.01, hours_after = 24.0)

Set records after precipitation to false in `:valid` column.

# Arguments
- df: DataFrame with columns `:datetime` and `:precip` sorted by `:datetime`
  in increasing order.
optional:
- `min_precip` (in mm per timestep): minimum precip to be considered effective
  precipitation.
- `hours_after`: time after the precipitation event to be considered invalid

# Value
`df` with modified column `:valid`.
"""
function setinvalid_afterprecip!(df; min_precip = 0.01, hours_after = 24.0)
  dfo = @pipe copy(df) |>
    select!(_, :datetime, :precip) |>
    subset!(_, :precip => ByRow(>=(min_precip)); skipmissing = true) |>
    select!(_,
      :datetime => :p_start,
      :datetime => (x -> x .+ frac_hour(Second, hours_after)) => :p_end)
  dfno = get_nonoverlapping_periods(dfo)
  #
  if !hasproperty(df, :valid) df[!,:valid] .= true; end
  #ip = 1
  # since date.time is ordered we need to search for the next index only from last position
  i_start = i_end = 1
  for ip in 1:nrow(dfno)
      p_start, p_end = dfno.p_start[ip], dfno.p_end[ip]
      i_start = i_end - 1 + findfirst(>=(p_start), df.datetime[i_end:end])
      i_end = i_start - 1 + findlast(<=(p_end), df.datetime[i_start:end])
      df.valid[i_start:i_end] .= false
  end
  df
end


text:
"""
# module utils

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-06-07

Functions for plotting.
# Examples

```jldoctest
julia>
```
"""
module Utils
using PyPlot, Images, Statistics
using DrWatson
@quickactivate "MEngProject"

# saving plots
location = plotsdir()

# pretty layer names
layers = [
    "L6 (\$x\$)",
    "L6(\$x\$)",
    "L4 excit (\$y\$)",
    "L4 excit (\$y\$)",
    "L4 inhib (\$m\$)",
    "L4 inhib (\$m\$)",
    "L2/3 excit (\$z\$)",
    "L2/3 excit (\$z\$)",
    "L2/3 inhib (\$s\$)",
    "L2/3 inhib (\$s\$)",
    "LGN (\$v\$)",
    "LGN (\$v\$)",
]

# bare layer names
layers_1 = [
    "\$x_1\$",
    "\$x_2\$",
    "\$y_1\$",
    "\$y_2\$",
    "\$m_1\$",
    "\$m_2\$",
    "\$z_1\$",
    "\$z_2\$",
    "\$s_1\$",
    "\$s_2\$",
    "\$v^+\$",
    "\$v^-\$",
]

# short layer names
la = ["x", "x", "y", "y", "m", "m", "z", "z", "s", "s", "v", "v"]

# line types for matplotlib
lines = [
    "b--",
    "g--",
    "r-.",
    "c-.",
    "r:",
    "c:",
    "m-",
    "y-",
    "m:",
    "y:",
    "b-.",
    "g-.",
]

# colours for matplotlib
colours = [
    "tab:blue",
    "tab:orange",
    "tab:green",
    "tab:red",
    "tab:purple",
    "tab:brown",
    "tab:pink",
    "tab:gray",
    "tab:olive",
    "tab:cyan",
]

"""
Creates image plot with blue-red cmap.
"""
function plot_rb(
    img::AbstractArray;
    name = "img",
    save = false,
    axMin = -1,
    axMax = 1,
    clbar = false,
    loc = location,
    filetype = ".png",
)
    findmax(img)[1] > axMax && throw(ArgumentError(string(
        "Image has max ",
        findmax(img)[1],
        ",outside range",
    )))
    findmin(img)[1] < axMin && throw(ArgumentError(string(
        "Image has min ",
        findmin(img)[1],
        ",outside range",
    )))
    fig, ax = plt.subplots()

    im = ax.imshow(img, cmap = matplotlib.cm.RdBu_r, vmax = axMax, vmin = axMin)
    if clbar
        cbar = fig.colorbar(im, shrink = 0.9, ax = ax)
    end

    plt.axis("off")
    fig.tight_layout()
    plt.show()
    if save
        plt.savefig(string(loc, name, filetype))
    end
end

"""
Creates image plot with two images overlayed.
"""
function plot_two(
    img::AbstractArray,
    k,
    t;
    name = "img",
    save = false,
    axMin = -1,
    axMax = 1,
    clbar = false,
    loc = location,
    filetype = ".png",
)
    findmax(img)[1] > axMax && throw(ArgumentError(string(
        "Image has max ",
        findmax(img)[1],
        ",outside range",
    )))
    findmin(img)[1] < axMin && throw(ArgumentError(string(
        "Image has min ",
        findmin(img)[1],
        ",outside range",
    )))
    fig, ax = plt.subplots()

    im = ax.imshow(
        img[:, :, k, 1, t],
        cmap = matplotlib.cm.BrBG,
        vmax = axMax,
        vmin = axMin,
    )
    ax.imshow(
        img[:, :, k+1, 1, t],
        cmap = matplotlib.cm.BrBG,
        vmax = axMax,
        vmin = axMin,
    )
    if clbar
        cbar = fig.colorbar(im, shrink = 0.9, ax = ax)
    end

    plt.axis("off")
    fig.tight_layout()
    plt.show()
    if save
        plt.savefig(string(loc, name, filetype))
    end
end

"""
Creates image plot with grey cmap.
"""
function plot_gs(
    img::AbstractArray;
    name = "img",
    save = false,
    axMin = 0,
    axMax = 2,
    clbar = false,
    loc = location,
    filetype = ".png",
)
    findmax(img)[1] > axMax && throw(ArgumentError(string(
        "Image has max ",
        findmax(img)[1],
        ",outside range",
    )))
    findmin(img)[1] < axMin && throw(ArgumentError(string(
        "Image has min ",
        findmin(img)[1],
        ",outside range",
    )))
    fig, ax = plt.subplots()

    im = ax.imshow(img, cmap = matplotlib.cm.gray, vmax = axMax, vmin = axMin)
    if clbar
        cbar = fig.colorbar(im, shrink = 0.9, ax = ax)
    end

    plt.axis("off")
    fig.tight_layout()
    plt.show()
    if save
        plt.savefig(string(loc, name, filetype))
    end
end

"""
Plots all orientations with blue-red cmap.
"""
function save_orientations_rb(
    A::Array,
    name::String,
    axMin = -1,
    axMax = 1,
    clbar = false,
    loc = location,
    filetype = ".png",
)
    for k = 1:size(A)[3]
        fn = string(name, "_", k)
        plot_rb(A[:, :, k], fn, true, axMin, axMax, clbar, loc, filetype)
    end
end

"""
Plots all orientations with gray cmap.
"""
function save_orientations_gs(
    A::Array,
    name::String,
    axMin = 0,
    axMax = 1,
    clbar = false,
    loc = location,
    filetype = ".png",
)
    for k = 1:size(A)[3]
        fn = string(name, "_", k)
        plot_gs(A[:, :, k], fn, true, axMin, axMax, clbar, loc, filetype)
    end
end

# function save_2d_list_rb(
#     A::Array,
#     name::String,
#     loc = location,
#     filetype = ".png",
#     axMin = -1,
#     axMax = 1,
#     clbar = false,
# )
#     n = 0
#     for img in A
#         n += 1
#         plot_101_rb(img, name, true, axMin, axMax, clbar, loc, filetype)
#     end
# end
#
#
# function save_orientations_gs_(
#     A::Array,
#     name::String,
#     bright = 1,
#     loc = location,
#     filetype = ".png",
# )
#     for k = 1:size(A)[3]
#         fn = string(loc, name, "_", k, filetype)
#         save(fn, Gray.(bright .* (A[:, :, k])))
#     end
# end
#
#
# function save_2d_gs(
#     A::Array,
#     name::String,
#     bright = 1,
#     loc = location,
#     filetype = ".png",
# )
#     fn = string(loc, name, filetype)
#     save(fn, Gray.(bright .* (A)))
# end
#
# function save_2d_list_gs(
#     A::Array,
#     name::String,
#     bright = 1,
#     loc = location,
#     filetype = ".png",
# )
#     n = 0
#     for img in A
#         n += 1
#         save_2d(img, string(name, n), bright, loc, filetype)
#     end
# end

"""
Plots activation vs time with all layers and orientations.
Uses highest value pixel at end k=1 orientation for each layer.
"""
function plot_t_act(sol, name, batch, file; save = true)
    fig, axs = plt.subplots()
    @inbounds begin
        for k ∈ 1:12
            v1 = @view sol[:, :, k, 1, end]
            v2 = @view sol[findmax(v1)[2][1], findmax(v1)[2][2], k, 1, :]
            layer = Utils.layers_1[k]
            axs.plot(sol.t, v2, Utils.lines[k], label = "$layer", alpha = 0.8)
        end
        axs.set_xlabel("Time")
        axs.set_ylabel("Activation")
        plt.legend()
        fig.tight_layout()
        if save
            plt.savefig(plotsdir(
                string(name, batch),
                string(file, "_time.png"),
            ))
        end
        close("all")
    end
    return nothing
end

"""
Plots mean activation vs time with all layers and orientations.
Uses mean for each layer/orientation.
Currently very slow.
"""
function plot_t_act_mean(sol, name, batch, file; save = true)
    fig, axs = plt.subplots()
    @inbounds begin
        for k ∈ 1:12
            # v1 = @view sol[:, :, k, 1, end]
            # v2 = @view sol[findmax(v1)[2][1], findmax(v1)[2][2], k, 1, :]
			v2 = Array{eltype(float32)}(undef, size(sol, 5))
			for s in 1:size(sol, 5)
				v2[s] = mean(@view sol[:,:, k, 1, s])
			end
			layer = Utils.layers_1[k]
			axs.plot(sol.t, v2, Utils.lines[k], label = "$layer", alpha = 0.8)
        end
        axs.set_xlabel("Time")
        axs.set_ylabel("Activation")
        plt.legend()
        fig.tight_layout()
        if save
            plt.savefig(plotsdir(
                string(name, batch),
                string(file, "_time.png"),
            ))
        end
        close("all")
    end
    return nothing
end

"""
Plots activation vs time with all layers and orientations.
Uses highest value pixel at end of specified layer/orient.
Default is z, k=1
"""
function plot_t_act_px(sol, name, batch, file; d=7, save = true)
    fig, axs = plt.subplots()
    @inbounds begin
        v1 = @view sol[:, :, d, 1, end]
        for k ∈ 1:12
            v2 = @view sol[findmax(v1)[2][1], findmax(v1)[2][2], k, 1, :]
            layer = Utils.layers_1[k]
            axs.plot(sol.t, v2, Utils.lines[k], label = "$layer", alpha = 0.8)
        end
        axs.set_xlabel("Time")
        axs.set_ylabel("Activation")
        plt.legend()
        fig.tight_layout()
        if save
            plt.savefig(plotsdir(
                string(name, batch),
                string(file, "_time.png"),
            ))
        end
        close("all")
    end
    return nothing
end

"""
Plots activation vs time with all layers and orientations.
Takes in pixel address.
Default: (50,50)
Default is z, k=1
"""
function plot_t_act_spec(sol, name, batch, file; px=(50,50), save = true)
    fig, axs = plt.subplots()
    @inbounds begin
#         v1 = @view sol[:, :, d, 1, end]
        for k ∈ 1:12
            v2 = @view sol[px[1], px[2], k, 1, :]
            layer = Utils.layers_1[k]
            axs.plot(sol.t, v2, Utils.lines[k], label = "$layer", alpha = 0.8)
        end
        axs.set_xlabel("Time")
        axs.set_ylabel("Activation")
        plt.legend()
        fig.tight_layout()
        if save
            plt.savefig(plotsdir(
                string(name, batch),
                string(file, "_time.png"),
            ))
        end
        close("all")
    end
    return nothing
end

"""
Plots two orientantions together for all layers.
v^+ and v^- plotted together
Cbar alpha for lower image adjusted.
"""
function plot_k2(sol, t, name, batch, file; save = true)

    @inbounds begin

        v0 = @view sol(t)[:, :, :, 1]
        axMax = findmax(v0)[1]

        # plot x, y, z, m, s
        for k ∈ 1:2:10
            k2 = k + 1
            fig, ax = plt.subplots()

            v1 = @view sol(t)[:, :, k, 1]
            v2 = @view sol(t)[:, :, k+1, 1]
            im = ax.imshow(
                v1,
                cmap = matplotlib.cm.PRGn,
                vmax = axMax,
                vmin = -axMax,
            )
            im2 = ax.imshow(
                v2,
                cmap = matplotlib.cm.RdBu_r,
                vmax = axMax,
                vmin = -axMax,
                alpha = 0.5,
            )

            cbar = fig.colorbar(im2, shrink = 0.9, ax = ax)
            cbar.ax.set_xlabel("\$k=2\$")
            cbar = fig.colorbar(im, shrink = 0.9, ax = ax)
            cbar.set_alpha(0.5)
            cbar.draw_all()
            cbar.ax.set_xlabel("\$k=1\$")
            layer = Utils.layers[k]
            plt.title("Layer: $layer, \$t=$t\$")
            plt.axis("off")
            fig.tight_layout()
            if save
                plt.savefig(plotsdir(
                    string(name, batch),
                    string(file, "_", t, "_", Utils.la[k], ".png"),
                ))
            end
            close("all")
        end

        # plot lgn
        k = 11
        fig, ax = plt.subplots()
        v1 = @view sol(t)[:, :, k, 1]
        v2 = @view sol(t)[:, :, k+1, 1]
        im = ax.imshow(
            v1,
            cmap = matplotlib.cm.PRGn,
            vmax = axMax,
            vmin = -axMax,
        )
        im2 = ax.imshow(
            v2,
            cmap = matplotlib.cm.RdBu_r,
            vmax = axMax,
            vmin = -axMax,
            alpha = 0.5,
        )

        cbar = fig.colorbar(im2, shrink = 0.9, ax = ax)
        cbar.ax.set_xlabel("\$v^-\$")
        cbar = fig.colorbar(im, shrink = 0.9, ax = ax)
        cbar.ax.set_xlabel("\$v^+\$")
        cbar.set_alpha(0.5)
        cbar.draw_all()

        layer = Utils.layers[k]
        plt.title("Layer: $layer, \$t=$t\$")
        plt.axis("off")
        fig.tight_layout()
        if save
            plt.savefig(plotsdir(
                string(name, batch),
                string(file, "_", t, "_", Utils.la[k], ".png"),
            ))
        end
        close("all")
    end
    return nothing
end

end


text:
wrap_return(ex, type, jtype, next_types = nothing) = @match t = type begin
    :VkResult => :(@check($ex))
    :Cstring => :(unsafe_string($ex))
    GuardBy(is_opaque_pointer) => ex

    # Call handle constructor.
    GuardBy(in(spec_handles.name)) => :($(remove_vk_prefix(t))($ex))
    # Don't change enumeration variables since they won't be wrapped under a new name
    GuardBy(in(spec_enums.name)) => ex
    # Vulkan and idiomatic Julia types are the same (up to aliases).
    if is_fn_ptr(type) || follow_constant(type) == jtype || innermost_type(type) ∈ spec_flags.name
    end => ex
    # Fall back to the from_vk function for conversion.
    _ => isnothing(next_types) ? :(from_vk($jtype, $ex)) : :(from_vk($jtype, $ex, $next_types))
end

_wrap_implicit_return(params::AbstractVector{SpecFuncParam}; with_func_ptr = false) =
    length(params) == 1 ? _wrap_implicit_return(first(params); with_func_ptr) :
    Expr(:tuple, _wrap_implicit_return.(params; with_func_ptr)...)

"""
Build a return expression from an implicit return parameter.
Implicit return parameters are pointers that are mutated by the API, rather than returned directly.
API functions with implicit return parameters return either nothing or a return code, which is
automatically checked and not returned by the wrapper.
Such implicit return parameters are `Ref`s or `Vector`s holding either a base type or a core struct Vk*.
They need to be converted by the wrapper to their wrapping type.
"""
function _wrap_implicit_return(return_param::SpecFuncParam, next_types = nothing; with_func_ptr = false)
    p = return_param
    @assert is_ptr(p.type) "Invalid core type for an implicit return. Expected $(p.type) <: Ptr"
    pt = follow_alias(ptr_type(p.type))
    ex = @match p begin

        # array pointer
        GuardBy(is_arr) => @match ex = wrap_return(p.name, pt, innermost_type((idiomatic_julia_type(p))), next_types) begin
            ::Symbol => ex
            ::Expr => broadcast_ex(ex) # broadcast result
        end

        # pointer to a unique object
        GuardBy(is_data_with_retrievable_size) => begin
            @assert has_length(p)
            size = len(p)
            @assert is_size(size)
            :($(size.name)[], $(p.name))
        end
        _ => wrap_return(:($(p.name)[]), pt, innermost_type((idiomatic_julia_type(p))), next_types) # call return_expr on the dereferenced pointer
    end

    @match p begin
        if pt ∈ spec_handles.name
        end => wrap_implicit_handle_return(parent_spec(return_param), handle_by_name(pt), ex, with_func_ptr)
        _ => ex
    end
end

function wrap_implicit_return(spec::SpecFunc, args...; kwargs...)
    ex = _wrap_implicit_return(args...; kwargs...)
    if spec.name == :vkCreateInstance || spec.name == :vkCreateDevice
        # Insert an expression for the automatic dispatch.
        ex = :(@fill_dispatch_table $ex)
    end
    must_return_success_code(spec) ? :(($ex, _return_code)) : ex
end

function wrap_implicit_handle_return(handle::SpecHandle, ex::Expr, parent_handle::SpecHandle, parent_ex, destroy::Bool, with_func_ptr)
    @match ex begin
        :($f($v[])) => :($f($v[], $(!destroy ? :identity : destructor(handle, with_func_ptr)), $parent_ex))
        :($f.($v)) => :($f.($v, $(!destroy ? :identity : destructor(handle, with_func_ptr)), $parent_ex))
    end
end

function wrap_implicit_handle_return(handle::SpecHandle, ex::Expr, destroy::Bool, with_func_ptr)
    @match ex begin
        :($f($v[])) => :($f($v[], $(!destroy ? :identity : destructor(handle, with_func_ptr))))
        :($f.($v)) => :($f.($v, $(!destroy ? :identity : destructor(handle, with_func_ptr))))
    end
end

function wrap_implicit_handle_return(spec::SpecFunc, handle::SpecHandle, ex::Expr, with_func_ptr)
    destroy = spec.type ≠ FTYPE_QUERY
    args = @match parent_spec(handle) begin
        ::Nothing => (handle, ex)
        p::SpecHandle => @match spec.type begin
            &FTYPE_CREATE || &FTYPE_ALLOCATE => (handle, ex, p, retrieve_parent_ex(p, create_func(spec)))
            _ => (handle, ex, p, retrieve_parent_ex(p, spec)::Symbol)
        end
    end

    wrap_implicit_handle_return(args..., destroy, with_func_ptr)
end

function wrap_return_type(spec::SpecFunc, ret_type)
    if must_return_success_code(spec) && has_implicit_return_parameters(spec)
        ret_type = :(Tuple{$ret_type,$(enum_type(:VkResult))})
    end

    @match spec.return_type begin
        :VkResult => :(ResultTypes.Result{$ret_type,VulkanError})
        _ => ret_type
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
import Pkg;

Pkg.add("CSV") # install CSV
Pkg.add("DataFrames") # install Dataframes
Pkg.add("TimeSeries")
Pkg.add("Pandas")
Pkg.add("Plots") # plotting dataframe
Pkg.add("PyCall")
Pkg.add("RollingFunctions") # rolling functions

# import libraries
using CSV;
using Dates;
using DataFrames;
using Statistics;
using Plots;
using PyCall;
using RollingFunctions;
using Printf;

@pyimport matplotlib.pyplot as plt

# load data
df = CSV.File("../../database/microeconomic_data/hkex_ticks_day/hkex_0001.csv") |> DataFrame

ticker = "0001.HK"

# moving average crossover
signals = df[:,[:Date, :Close]]

dates = Array(convert(Matrix, select(df, :Date))) # get dates
close = convert(Matrix, select(df, :Close)) # get closing price

# short MA
short_window = 40
short_mavg = runmean(vec(close), short_window)
insertcols!(signals, 1, :short_mavg => short_mavg)

# long MA
long_window = 100
long_mavg = runmean(vec(close), long_window)
insertcols!(signals, 1, :long_mavg => long_mavg)

# Create signals
signal = Float64[]

for i in 1:length(short_mavg)
    if short_mavg[i] > long_mavg[i]
        x = 1.0
    else
        x = 0.0
    end
    push!(signal, x)
end

insertcols!(signals, 1, :signal => signal)

# Generate positions
function gen_pos(signal)
    positions = zeros(length(signal))
    positions[1] = 0
    for i in 2:length(signal)
        positions[i] = signal[i] - signal[i-1]
    end
    return positions
end

positions = gen_pos(signal)

insertcols!(signals, 1, :positions => positions)

# Generate tmp arrays to plot buy signals

buy_signals = DataFrame()
buy_dates = []
buy_prices = []

for i in 1:length(positions)
    if (positions[i] == 1.0)
       push!(buy_dates, dates[i])
       push!(buy_prices, close[i])
    end
end

insertcols!(buy_signals, 1, :Date => buy_dates)
insertcols!(buy_signals, 1, :Price => buy_prices)
#print(first(buy_signals,10))

# Generate tmp arrays to plot sell signals

sell_signals = DataFrame()
sell_dates = []
sell_prices = []

for i in 1:length(positions)
    if (positions[i] == -1.0)
       push!(sell_dates, dates[i])
       push!(sell_prices, close[i])
    end
end

insertcols!(sell_signals, 1, :Date => sell_dates)
insertcols!(sell_signals, 1, :Price => sell_prices)
#print(first(sell_signals,10)

# Plot signals and MAs

fig = plt.figure() # Initialise the plot figure
ax1 = fig.add_subplot(111, ylabel="Price in \$") # Add a subplot and label for y-axis

plt.plot(signals.Date, signals.short_mavg, color="blue", linewidth=1.0, label="Short MA")
plt.plot(signals.Date, signals.long_mavg, color="orange", linewidth=1.0, label="Long MA")
#plt.plot(signals.Date, signals.Close, color="red", linewidth=0.8, label="Closing price")
plt.plot(buy_signals.Date, buy_signals.Price, marker=10, markersize=7, color="m", linestyle="None", label="Buy signal")
plt.plot(sell_signals.Date, sell_signals.Price, marker=11, markersize=7, color="k", linestyle="None", label="Sell signal")

plt.title("MA crossover signals")
plt.show()

# save fig
fig.savefig("./figures/moving-average-crossover_signals", dpi=100)

# ============== Backtest ============== #

initial_capital = 100000.0

# Initialise the portfolio with value owned
portfolio = signals[:,[:Date, :Close, :positions]]
portfolio[:trade] = signals[:Close] .* (100 .* signals[:positions])

# Add `holdings` to portfolio
portfolio[:quantity] = cumsum(100 .* signals[:positions])
portfolio[:holdings] = portfolio[:Close] .* portfolio[:quantity]

# Add `cash` to portfolio
portfolio[:cash] = initial_capital .- cumsum(portfolio[:trade])

# Add `total` to portfolio
portfolio[:total] = portfolio[:cash] .+ portfolio[:holdings]

portfolio_total = Array(portfolio[:total])

# Generate returns
function gen_returns(portfolio_total)
    returns = zeros(length(portfolio_total))
    returns[1] = 0
    for i in 2:length(portfolio_total)
        returns[i] = (portfolio_total[i] - portfolio_total[i-1]) / portfolio_total[i-1]
    end
    return returns
end

returns = gen_returns(portfolio_total)
insertcols!(portfolio, 1, :returns => returns)

#print(first(portfolio, 100))
@printf("Final total value: %f\n", portfolio.total[size(portfolio,1)])

total_return = (portfolio.total[size(portfolio,1)] - portfolio.total[1]) / portfolio.total[1]
@printf("Total return: %f\n", total_return)

# ============== Evaluation metrics ============== #

function portfolio_return(portfolio)
    fig = plt.figure() # Initialise the plot figure
    ax1 = fig.add_subplot(111, ylabel="Total in \$") # Add a subplot and label for y-axis

    plt.plot(portfolio.Date, portfolio.returns, color="blue", linewidth=1.0, label="Returns")

    plt.title("MA crossover portfolio return")
    plt.show()

    # save fig
    fig.savefig("./figures/moving-average-crossover_returns", dpi=100)
end

function sharpe_ratio(portfolio)
    # Annualised Sharpe ratio
    sharpe_ratio = sqrt(252) * (mean(returns) / std(returns))

    return sharpe_ratio
end

function CAGR(portfolio)
    # Get the number of days in df
    format = DateFormat("y-m-d")
    days = portfolio.Date[size(portfolio,1)] - portfolio.Date[1]

    # Calculate the CAGR
    cagr = ^((portfolio.total[size(portfolio,1)] / portfolio.total[1]), (252.0 / Dates.value(days))) - 1

    return cagr
end

portfolio_return(portfolio)
sharpe = sharpe_ratio(portfolio)
cagr = CAGR(portfolio)

@printf("Sharpe ratio: %f\n", sharpe)
@printf("CAGR: %f\n", cagr)


text:
module VertexFinder

import Base: isempty

"""
    Queue{T}

Simple queue type, having several shortcomings but it's working in this case
well enough.
"""
mutable struct Queue{T}
    data :: T
    head :: Int
    tail :: Int
end

function Queue(T, n::Int)
    return Queue(zeros(T, n), 1, 1)
end

function isempty(q::Queue)
    return q.head == q.tail
end

function enqueue!(q::Queue, item)
    q.data[q.tail] = item
    q.tail += 1
    if q.tail > length(q.data)
        resize!(q.data, length(q.data)*2)
    end
end

function dequeue!(q::Queue)
    if isempty(q)
        error("Queue is empty.")
    end
    item = q.data[q.head]
    q.head += 1
    return item
end

"""
    find(G, v)

Find a better estimate of a pseudoperipheral vertex `w` of graph `G`, starting
from a vertex `v`. George-Liu algorithm is used to find pseudoperipheral vertex.
"""
function find(G, v)
    n = length(G)

    queue = Queue(Int, n)
    visited = zeros(Bool, n)
    dist = zeros(Int, n)

    enqueue!(queue, v)
    visited[v] = true

    n = 0
    while !isempty(queue)
        n += 1
        v = dequeue!(queue)
        adjacent_vertices = getindex(G, v)
        for w in adjacent_vertices
            visited[w] && continue
            visited[w] = true
            enqueue!(queue, w)
            dist[w] = dist[v] + 1
        end
    end

    maxdist = maximum(dist)
    mindeg = 2^32
    w = -1

    for i=1:n
        dist[i] != maxdist && continue
        deg = length(G[i])
        if deg < mindeg
            mindeg = deg
            w = i
        end
    end

    return w
end

end # module