text:
#=
# DG test
# 1D Poisson using SIPG formulation for DG
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("test");

@useLog("dgtestlog")

@domain(1)
@solver(DG)
@functionSpace(LEGENDRE, 2)

@mesh(LINEMESH, 10)

@variable(u)

@testSymbol(v)

@boundary(u, 1, DIRICHLET, 0)

@coefficient(f, "-pi*pi*sin(pi*x)")
@coefficient(beta, 1)
@weakForm(u, "dot(grad(u),grad(v)) + f*v - surface( ave_normdotgrad(u) * jump(v)) - surface(jump(u) * ave_normdotgrad(v)) + surface(beta*jump(u)*jump(v))")
@weakForm(u, "dot(grad(u),grad(v)) + f*v - surface( ave_normdotgrad(u) * jump(v)) - surface(jump(u) * ave_normdotgrad(v)) + surface(beta*jump(u)*jump(v))")

solve(u);

# solution is stored in the variable's "values"
using Plots
pyplot();
display(plot(Femshop.grid_data.allnodes[1,:], u.values[1,:], reuse=false))

@finalize()


text:
"""
Returns an array of tuples of directions.

```julia
    get_directions(dim::Integer; include_zero::Bool = false) -> Array{NTuple{N, Integer}, 1}

julia> get_directions(3)
27-element Array{Tuple{Int64,Int64,Int64},1}:
 (-1, -1, -1)
 (0, -1, -1)
 (1, -1, -1)
 (-1, 0, -1)
 (0, 0, -1)
 (1, 0, -1)
 (-1, 1, -1)
 (0, 1, -1)
 (1, 1, -1)
 (-1, -1, 0)
 (0, -1, 0)
 (1, -1, 0)
 (-1, 0, 0)
 (0, 0, 0)
 (1, 0, 0)
 (-1, 1, 0)
 (0, 1, 0)
 (1, 1, 0)
 (-1, -1, 1)
 (0, -1, 1)
 (1, -1, 1)
 (-1, 0, 1)
 (0, 0, 1)
 (1, 0, 1)
 (-1, 1, 1)
 (0, 1, 1)
 (1, 1, 1)
```
"""
function get_directions(dim::T; include_zero::Bool = false) where T <: Integer
    D = reshape(collect(Iterators.product([-one(T):one(T) for _ in 1:dim]...)), :)
    return include_zero ? D : filter(d -> d ≠ 𝟘(dim), D)
end
get_directions(M::AbstractArray{T, N}; include_zero::Bool = false) where {T, N} = get_directions(ndims(M), include_zero = include_zero)
# get_directions(M::AbstractArray{T, N}; include_zero::Bool = false) where {T, N} = get_directions(N, include_zero = include_zero)

function _expand_as_required(M::AbstractArray{T, N}, expand_by::T, inds::AbstractIndexOrIndices{N}...) where {T, N}
    indices = Vector{Union{T, Nothing}}()
    ind_extrema = extrema_indices(inds...)

    for d in 1:ndims(M)
        ith_extrema = ind_extrema[d]
        if ! all(map(m -> m ≤ size(M, d), ith_extrema))
            for invalid_idx in filter(i -> i > size(M, i) || i < 1, ith_extrema)
                difference = invalid_idx - size(M, d)
                M = sign(difference) == 1 ? append_n_times(M, abs(difference), expand_by, dims = d) : append_n_times_backwards(M, abs(difference), expand_by, dims = d)
            end
        end
    end

    return M
end

"""
    expand_as_required(M::AbstractArray{T, N}, expand_by::T, inds::AbstractIndexOrIndices{N})

Given indices, `expand_as_required` will fill in a matrix with `expand_by` where needed if such indices are not currently accessible.  If no `expand_by` is given, and the element type of the matrix is a number, fills with zero.

```julia
julia> A = rand(Int8, 2, 2)
2×2 Array{Int8,2}:
 -96   83
 -94  -39

julia> expand_as_required(A, zero(Int8), [(1, 1), (2, 2), (3, 3)])
3×3 Array{Int8,2}:
 -96   83  0
 -94  -39  0
   0    0  0
```
"""
expand_as_required(M::AbstractArray{T, N}, expand_by::T, inds::AbstractIndices{N}) where {T, N} = _expand_as_required(M, expand_by, inds)
expand_as_required(M::AbstractArray{T, N}, expand_by::T, inds::AbstractIndex{N}...) where {T, N} = _expand_as_required(M, expand_by, inds)
expand_as_required(M::AbstractArray{T, N}, inds::AbstractIndices{N}) where {T <: Number, N} = _expand_as_required(M, zero(T), inds)
expand_as_required(M::AbstractArray{T, N}, inds::AbstractIndex{N}...) where {T <: Number, N} = _expand_as_required(M, zero(T), inds)

"""
    adjacencies(M::AbstractArray{T, N}, idx::AbstractIndex{N})

`adjacencies` will get all adjacent element to an index in a given matrix.
"""
function adjacencies(M::AbstractArray{T, N}, idx::AbstractIndex{N}) where {T, N}
    direction_multipliers = get_directions(M, include_zero = false)
    adjacent_indices = AbstractIndex{N}[idx .+ j for j in direction_multipliers]
    return T[k for k in tryindex(M, adjacent_indices) if ! isnothing(k)]
end

"""
    expanded_adjacencies(M::Array{T, N}, expand_by::T, idx::AbstractIndex{N})

`expanded_adjacencies` will get all adjacencies of an index in the matrix M, given that the matrix is infinitely expanding by a single element beyond the specified M.  See also `adjacencies` and `expand_as_required`.  If no `expand_by` is given, and the element type of the matrix is a number, fills with zero.
"""
function expanded_adjacencies(M::Array{T, N}, expand_by::T, idx::AbstractIndex{N}) where {T, N}
    𝟎 = ntuple(_ -> zero(Int), N)
    idx = Tuple(idx)
    direction_multipliers = get_directions(M, include_zero = false)
    adjacent_indices = AbstractIndex{N}[idx .+ j for j in direction_multipliers]
    M = expand_as_required(M, expand_by, adjacent_indices)
    idx_shift = abs.(size(M) .- idx) # in case we have had to shift the array for non-positive indices
    return T[M[(i .+ idx_shift .- 1)...] for i in adjacent_indices]
end
expanded_adjacencies(M::Array{T, N}, idx::AbstractIndex{N}) where {T <: Number, N} = expanded_adjacencies(M, zero(T), idx)

"""
    n_adjacent_to(M::AbstractArray{T, N}, idx::AbstractIndex{N}, adj_elem::T)

Given a matrix M, counts the number of adjacent elements to index `idx` that are exactly the `adj_elem`.
"""
function n_adjacent_to(M::AbstractArray{T, N}, idx::AbstractIndex{N}, adj_elem::T) where {N, T}
    n = 0

    for i in adjacencies(M, idx)
        if i == adj_elem
            n += 1
        end
    end

    return n
end

"""
    expanded_n_adjacent_to(M::AbstractArray{T, N}, expand_by::T, idx::AbstractIndex{N}, adj_elem::T)

Given a matrix M, counts the number of adjacent elements to index `idx` that are exactly the `adj_elem`, expanding the matrix if needed (see `expanded_adjacencies`).  If no `expand_by` is given, and the element type of the matrix is a number, fills with zero.
"""
function expanded_n_adjacent_to(M::AbstractArray{T, N}, expand_by::T, idx::AbstractIndex{N}, adj_elem::T) where {N, T}
    n = 0

    for i in expanded_adjacencies(M, expand_by, idx)
        if i == adj_elem
            n += 1
        end
    end

    return n
end
expanded_n_adjacent_to(M::AbstractArray{T, N}, idx::AbstractIndex{N}, adj_elem::T) where {N, T <: Number} = expanded_n_adjacent_to(M, zero(T), idx, adj_elem)

"""
    global_adjacencies(M::AbstractArray{T, N}, idx::AbstractIndex{N}, ignored_elem::T)

Returns the "globally" adjacent indices in arbitrary positions in the cardinal directions from a specified index in matrix `M`, ignoring certain adjacent elements (i.e., skipping over them).
"""
function global_adjacencies_indices(M::AbstractArray{T, N}, idx::AbstractIndex{N}, ignored_elem::T) where {T, N}
    adjacent_indices = AbstractIndex{N}[]
    directional_shifts = get_directions(M, include_zero = false)
    n_adjacent, adjacent_count = n_adjacencies(M), 0

    while adjacent_count < n_adjacent
        for directional_shift in directional_shifts
            while true
                adj_index = idx .+ directional_shift

                if nothing ∈ tryindex(M, adj_index)
                    n_adjacent -= 1
                    break
                end

                if M[adj_index...] ≠ ignored_elem
                    adjacent_count += 1
                    push!(adjacent_indices, adj_index)
                    break
                else
                    directional_shift = (abs.(directional_shift) .+ 1) .* sign.(directional_shift)
                end
            end
        end
    end

    return adjacent_indices
end

"""
    global_adjacencies(M::AbstractArray{T}, idx::AbstractIndex{N}, ignored_elem::T)

Using `global_adjacencies_indices`, returns the elements of each globally adjacent index.  This function is mainly used for `global_n_adjacent_to`.
"""
function global_adjacencies(M::AbstractArray{T, N}, idx::AbstractIndex{N}, ignored_elem::T) where {N, T}
    return T[M[i...] for i in global_adjacencies_indices(M, idx, ignored_elem)]
end

"""
    global_n_adjacent_to(M::AbstractArray{T, N},idx::AbstractIndex{N}, adj_elem::T, ignored_elem::T) where {N, T}

Given a matrix M, counts the number of adjacent elements to index `idx` that are exactly the `adj_elem`, skipping over certain elements if needed (see `expanded_adjacencies`).
"""
function global_n_adjacent_to(M::AbstractArray{T, N},idx::AbstractIndex{N}, adj_elem::T, ignored_elem::T) where {N, T}
    n = 0

    for i in global_adjacencies(M, idx, ignored_elem)
        if i == adj_elem
            n += 1
        end
    end

    return n
end

ext: jl


text:
module TestApp
  include("../../CellDataTests.jl")
  include("../../FESpacesTests.jl")
  include("../../GeometryTests.jl")
  include("../../MultiFieldTests.jl")
  include("../../PLaplacianTests.jl")
  include("../../PoissonTests.jl")
  include("../../PeriodicBCsTests.jl")
end


text:
# a more efficient checkbounds implementation to minimize the overhead
_checkbounds(::Type{Bool}, A::AbstractArray, I::CartesianIndex) = _in(I, CartesianIndices(A))
function _checkbounds(::Type{Bool}, A::AbstractArray, p::CartesianIndex, q::CartesianIndex)
    _checkbounds(Bool, A, CartesianIndex(p.I..., q.I...))
end

# # short-circuiting gives better performance
_in(i::CartesianIndex, r::CartesianIndices) = _alltsc(map(in, i.I, r.indices))
@inline function _alltsc(t::Tuple{Bool,Vararg{Bool}})
    first(t) || return false
    return _alltsc(Base.tail(t))
end
_alltsc(::Tuple{}) = true


text:
using PowerModels
using Ipopt
using DelimitedFiles
using Memento
using ArgParse
include("qc_relaxation_methods.jl")

function parse_cli()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--log", "-l"
            help = "Define log level. Possible values are: [debug | info |
            notice | warn | error]."
            arg_type = String
            default = "info"
        "--vol", "-v"
            help = "Save A,b polytope data. Define the frequency of saving files
            e.g.: -v 1000 means that A,b data saved in every 1000th iteration."
            arg_type = Int
            default = 0
        "--bit"
            help = "Bound tightening ITerations. Max number of iterations in OBBT."
            arg_type = Int
            default = 5
        "--bto"
            help = "Bound tightening Time Out. Max time spent in OBBT."
            arg_type = Int
            default = 120
        "--tightened", "-t"
            help = "Tightened output case file name (relative or full path with
            extension). Supported file formats:\n\t- matpower v2.0"
            arg_type = String
        "--out_file_name_A", "-A"
            help = "Output file of the convex polytope's A matrix. File fullpath
            should be defined including the extension (.csv)."
            arg_type = String
        "--out_file_name_B", "-b"
            help = "Output file of the convex polytope's b vector. File fullpath
            should be defined including the extension (.csv)."
        "--out_file_name_X", "-x"
            help = "Output file of optimal set points. The first row contains the
            variable names."
            arg_type = String
        "case_file_name"
            help = "Case file name (relative or full path with extension).
            Supported file formats:\n\t- matpower v2.0"
            required = true
        "number_of_iterations"
            help = "Number of polytopes added."
            arg_type = Int
            required = true
        "out_samples_file_name"
            help = "Output file containing N2 number of uniform samples from the
            convex polytope. File fullpath should be defined including the
            extension (.csv)."
            arg_type = String
            required = true
        "n2_number_of_samples"
            help = "N2 number of samples from the convex polytope."
            arg_type = Int
            required = true
    end
    # to debug: parse_args([raw"C:\Users\Timon\OneDrive - Danmarks Tekniske Universitet\Denmark\DTU\2019_20_II\software\case_files\case9.m","10",raw"C:\Users\Timon\OneDrive - Danmarks Tekniske Universitet\Denmark\DTU\2019_20_II\software\.data\case9_QCRM_1.csv","10"],s)
    return parse_args(s)
end

function main()
    """ Run N number of iterations that executes a modified QC relaxed ACOPF to created
        separating hyper-planes and event reduce."""
    parsed_inputs = parse_cli()
    # Logger settings
    # Silence all: PowerModels.silence()
    logger = Memento.config!(parsed_inputs["log"]; fmt="[{level} | {name}]: {msg}")
    Memento.setlevel!(Memento.getlogger(PowerModels.InfrastructureModels), "error")
    if parsed_inputs["log"] == "debug"
        # enable debug:
        ENV["JULIA_DEBUG"] = "all"
    end

    # Input data: number of iterations , model name
    debug(logger, string("ARGS: ", parsed_inputs))

    start = time();
    # Load network data from file
    Memento.setlevel!(Memento.getlogger(PowerModels), "error")
    network_data = PowerModels.parse_file(parsed_inputs["case_file_name"])
    Memento.setlevel!(Memento.getlogger(PowerModels), "info")
    # Run bound tightening
    info(logger, string("Run obbt for ", parsed_inputs["case_file_name"],"."))

    to = convert(AbstractFloat, parsed_inputs["bto"])
    # suppress Ipopt message
    local network_data_tight, stats
    @suppress_out network_data_tight, stats = PowerModels.run_obbt_opf!(network_data,
        optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0);
        model_constructor=ACPPowerModel, max_iter=parsed_inputs["bit"],
        time_limit=to);

    if  parsed_inputs["bit"] < 2
        warn(logger,"obbt opf max_iter should be set to higher!")
    end
    # MATPOWER export
    # update original data with solution

    if parsed_inputs["tightened"] != nothing
        info(logger, string("Writing tightened case file: ", parsed_inputs["tightened"],"."))
        update_data!(network_data, network_data_tight) #network_data_tight["solution"]
        f = open(parsed_inputs["tightened"],"w")
        export_matpower(f,network_data)
        close(f)
    end

    # Initialize abstract power system model
    #TODO: I think ACPPowerModel is the acopf formulation and PowerModels.QCRMPowerModel is the quadratic relaxation
    # ref: https://lanl-ansi.github.io/PowerModels.jl/stable/formulation-details/#Formulation-Details-1
    power_model = instantiate_model(network_data_tight, QCRMPowerModel, PowerModels.build_opf)

    info(logger,string("Run QC relaxation for, ",  parsed_inputs["case_file_name"],"."))
    A_out,b_out,x_opt,header,scaling_factor = run_qc_relax(power_model, parsed_inputs["number_of_iterations"],parsed_inputs["vol"])
    debug(logger,string("scaling_factor",scaling_factor))

    if parsed_inputs["out_file_name_A"] != nothing
        info(logger,string("Write file:\n\t- ", parsed_inputs["out_file_name_A"]))
        writedlm(parsed_inputs["out_file_name_A"], A_out, ',')
    end
    if parsed_inputs["out_file_name_B"] != nothing
        info(logger,string("Write file:\n\t- ", parsed_inputs["out_file_name_B"]))
        writedlm(parsed_inputs["out_file_name_B"], b_out, ',')
    end
    if parsed_inputs["out_file_name_X"] != nothing
        writedlm(parsed_inputs["out_file_name_X"], header, ',')
        open(parsed_inputs["out_file_name_X"],"a") do io
            writedlm(io, x_opt, ',')
        end
    end

    info(logger,string("Sample ",  parsed_inputs["n2_number_of_samples"]," points."))
    N2_out = sample_polyhedron(parsed_inputs["n2_number_of_samples"])
    N2_scaled = tt(scaling_factor,N2_out)
    writedlm(parsed_inputs["out_samples_file_name"], header, ',')
    open(parsed_inputs["out_samples_file_name"],"a") do io
        writedlm(io, N2_scaled, ',')
    end

    elapsed = time() - start
    info(logger,"Elapsed time: $elapsed sec .")
end

logger = Memento.config!("info"; fmt="[{level} | {name}]: {msg}")
main()

max_stars_repo_name: timonviola/PSDC

max_issues_repo_name: timonviola/PSDC

max_forks_repo_name: timonviola/PSDC


text:
############### in file hailstone.jl ###############
module Hailstone

function hailstone(n)
    ret = [n]
    while n > 1
        if n & 1 > 0
            n = 3n + 1
        else
            n = Int(n//2)
        end
        append!(ret, n)
    end
    return ret
end

export hailstone

end

if PROGRAM_FILE == "hailstone.jl"
    using Hailstone
    h = hailstone(27)
    n = length(h)
    println("The sequence of hailstone(27) is:\n $h.\nThis sequence is of length $n. It starts with $(h[1:4]) and ends with $(h[n-3:end]).")
end
############ in file moduletest.jl ####################
include("hailstone.jl")
using Hailstone
function countstones(mi, mx)
    lengths2occurences = Dict()
    mostfreq = mi
    maxcount = 1
    for i in mi:mx
        h = hailstone(i)
        n = length(h)
        if haskey(lengths2occurences, n)
            newoccurences = lengths2occurences[n] + 1
            if newoccurences > maxcount
                maxcount = newoccurences
                mostfreq = n
            end
            lengths2occurences[n] = newoccurences
        else
            lengths2occurences[n] = 1
        end
    end
    mostfreq, maxcount
end

nlen, cnt = countstones(1,99999)

print("The most common hailstone sequence length for hailstone(n) for 1 <= n < 100000 is $nlen, which occurs $cnt times.")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function ImGui_ImplGlfw_ErrorCallback(code::Cint, description::Ptr{Cchar})::Cvoid
    @error "GLFW ERROR: code $code msg: $(unsafe_string(description))"
    return nothing
end

function ImGui_ImplGlfw_MouseButtonCallback(window::Ptr{GLFWwindow}, button::Cint, action::Cint, mods::Cint)::Cvoid
    ctx = unsafe_pointer_to_objref(glfwGetWindowUserPointer(window))
    if ctx.PrevUserCallbackMousebutton != C_NULL
        ccall(ctx.PrevUserCallbackMousebutton, Cvoid, (Ptr{GLFWwindow}, Cint, Cint, Cint), window, button, action, mods)
    end

    if action == GLFW_PRESS && button ≥ 0 && button < length(ctx.MouseJustPressed)
        ctx.MouseJustPressed[button+1] = true
    end

    return nothing
end

function ImGui_ImplGlfw_ScrollCallback(window::Ptr{GLFWwindow}, xoffset::Cdouble, yoffset::Cdouble)::Cvoid
    ctx = unsafe_pointer_to_objref(glfwGetWindowUserPointer(window))
    if ctx.PrevUserCallbackScroll != C_NULL
        ccall(ctx.PrevUserCallbackScroll, Cvoid, (Ptr{GLFWwindow}, Cdouble, Cdouble), window, xoffset, yoffset)
    end

    io::Ptr{ImGuiIO} = igGetIO()
    io.MouseWheelH = unsafe_load(io.MouseWheelH) + Cfloat(xoffset)
    io.MouseWheel = unsafe_load(io.MouseWheel) + Cfloat(yoffset)

    return nothing
end

function ImGui_ImplGlfw_KeyCallback(window::Ptr{GLFWwindow}, key::Cint, scancode::Cint, action::Cint, mods::Cint)::Cvoid
    ctx = unsafe_pointer_to_objref(glfwGetWindowUserPointer(window))
    if ctx.PrevUserCallbackKey != C_NULL
        ccall(ctx.PrevUserCallbackKey, Cvoid, (Ptr{GLFWwindow}, Cint, Cint, Cint, Cint), window, key, scancode, action, mods)
    end

    io::Ptr{ImGuiIO} = igGetIO()
    if key ≥ 0 && key < length(unsafe_load(io.KeysDown))
        if action == GLFW_PRESS
            c_set!(io.KeysDown, key, true)
            ctx.KeyOwnerWindows[key+1] = window
        end
        if action == GLFW_RELEASE
            c_set!(io.KeysDown, key, false)
            ctx.KeyOwnerWindows[key+1] = C_NULL
        end
    end

    # modifiers are not reliable across systems
    io.KeyCtrl = c_get(io.KeysDown, GLFW_KEY_LEFT_CONTROL) || c_get(io.KeysDown, GLFW_KEY_RIGHT_CONTROL)
    io.KeyShift = c_get(io.KeysDown, GLFW_KEY_LEFT_SHIFT) || c_get(io.KeysDown, GLFW_KEY_RIGHT_SHIFT)
    io.KeyAlt = c_get(io.KeysDown, GLFW_KEY_LEFT_ALT) || c_get(io.KeysDown, GLFW_KEY_RIGHT_ALT)
    if Sys.iswindows()
        io.KeySuper = false
    else
        io.KeySuper = c_get(io.KeysDown, GLFW_KEY_LEFT_SUPER) || c_get(io.KeysDown, GLFW_KEY_RIGHT_SUPER)
    end

    return nothing
end

function ImGui_ImplGlfw_CharCallback(window::Ptr{GLFWwindow}, x::Cuint)::Cvoid
    ctx = unsafe_pointer_to_objref(glfwGetWindowUserPointer(window))
    if ctx.PrevUserCallbackChar != C_NULL
        ccall(ctx.PrevUserCallbackChar, Cvoid, (Ptr{GLFWwindow}, Cuint), window, x)
    end

    0 < x < 0x10000 && ImGuiIO_AddInputCharacter(igGetIO(), x)

    return nothing
end

function ImGui_ImplGlfw_MonitorCallback(monitor::Ptr{GLFWmonitor}, x::Cint)::Cvoid
    io::Ptr{ImGuiIO} = igGetIO()
    ctx = unsafe_pointer_to_objref(unsafe_load(io.UserData))
    ctx.WantUpdateMonitors = true
    return nothing
end

function ImGui_ImplGlfw_WindowCloseCallback(window::Ptr{GLFWwindow})::Cvoid
    viewport::Ptr{ImGuiViewport} = igFindViewportByPlatformHandle(window)
    if viewport != C_NULL
        viewport.PlatformRequestClose = true
    end
    return nothing
end

function ImGui_ImplGlfw_WindowPosCallback(window::Ptr{GLFWwindow}, x::Cint, y::Cint)::Cvoid
    viewport::Ptr{ImGuiViewport} = igFindViewportByPlatformHandle(window)
    if viewport != C_NULL
        data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
        if data != C_NULL
            ignore_event = igGetFrameCount() ≤ (unsafe_load(data.IgnoreWindowPosEventFrame) + 1)
            ignore_event && return nothing
        end
        viewport.PlatformRequestMove = true
    end
    return nothing
end

function ImGui_ImplGlfw_WindowSizeCallback(window::Ptr{GLFWwindow}, x::Cint, y::Cint)::Cvoid
    viewport::Ptr{ImGuiViewport} = igFindViewportByPlatformHandle(window)
    if viewport != C_NULL
        data::Ptr{ImGuiViewportDataGlfw} = unsafe_load(viewport.PlatformUserData)
        if data != C_NULL
            ignore_event = igGetFrameCount() ≤ (unsafe_load(data.IgnoreWindowSizeEventFrame) + 1)
            ignore_event && return nothing
        end
        viewport.PlatformRequestResize = true
    end
    return nothing
end


text:
using Documenter, ManifoldHMC

makedocs(
    modules = [ManifoldHMC],
    format = Documenter.HTML(prettyurls = false),
    checkdocs = :exports,
    sitename = "ManifoldHMC",
    pages = [
        "Overview" => "index.md",
        "(Co)geodesic Flow" => "flow.md"
    ]
)

deploydocs(
    repo = "github.com/sdaxen/ManifoldHMC.jl.git",
)


text:
# Utilities to filter occurrences of UMLS semantic type in pubmed articles
# Date: September 7, 2016
# Authors: Isabel Restrepo
# BCBI - Brown University
# Version: Julia 0.4.5

using BioMedQuery.PubMed
using DataFrames
using SparseArrays

"""
    umls_semantic_occurrences(db, umls_semantic_type)

Return a sparse matrix indicating the presence of MESH descriptors associated
with a given umls semantic type in all articles of the input database

## Output

* `des_ind_dict`: Dictionary matching row number to descriptor names
* `disease_occurances` : Sparse matrix. The columns correspond to a feature vector, where each row is a MESH descriptor. There are as many columns as articles. The occurance/abscense of a descriptor is labeled as 1/0
"""
function umls_semantic_occurrences(db, umls_concepts...)

    filtered_mesh = Set(filter_mesh_by_concept(db, umls_concepts...))

    #create a map of filtered descriptor name to index to guarantee order
    des_ind_dict = Dict{String, Int}()

    for (i, fm) in enumerate(filtered_mesh)
        des_ind_dict[fm]= i
    end

    articles = all_pmids(db)

    #create the data-matrix
    disease_occurances = spzeros(length(filtered_mesh), length(articles))

    #Can this process be more efficient using database join/select?
    narticle = 0

    for (i, pmid) in enumerate(articles)

        #not all mesh are of the desired semantic type
        article_filtered_mesh = PubMed.get_article_mesh_by_concept(db, pmid, umls_concepts...)

        #skip if empty
        if isempty(article_filtered_mesh)
            continue
        end

        #otherwise form feature vector for this article
        indices = []
        for d in article_filtered_mesh
            push!(indices, des_ind_dict[d])
        end

        #TO DO: Not sure about the type. Should we choose bool to save space
        # or float to support opperations
        article_dis_feature  = zeros(Int, (length(filtered_mesh),1))
        article_dis_feature[indices] .= 1

        #append to data matrix
        disease_occurances[:, i] = article_dis_feature
        narticle+=1
    end

    println("-------------------------------------------------------------")
    println("Found ", narticle, " articles with valid descriptors")
    println("-------------------------------------------------------------")
    return des_ind_dict, disease_occurances

end

# Retrieve all mesh descriptors associated with the given umls_concept
function filter_mesh_by_concept(db, umls_concepts...)

    concept_set_string = "'" * join(umls_concepts,"','") * "'"

    query = string("SELECT mesh FROM mesh2umls WHERE umls IN ($(concept_set_string))")
    println("Filter mesh query string : $(query)")

    sel  = db_query(db, query)
    return sel[1]
end

"""
    umls_semantic_occurrences(dfs, mesh2umls_df, umls_semantic_type)

Return a sparse matrix indicating the presence of MESH descriptors associated
with a given umls semantic type in all articles of the input database

## Output

* `des_ind_dict`: Dictionary matching row number to descriptor names
* `disease_occurances` : Sparse matrix. The columns correspond to a feature vector, where each row is a MESH descriptor. There are as many columns as articles. The occurance/abscense of a descriptor is labeled as 1/0
"""
function umls_semantic_occurrences(dfs::Dict{String,DataFrame}, mesh2umls_df::DataFrame, umls_concepts...)

    filtered_mesh = Set(filter_mesh_by_concept(mesh2umls_df, umls_concepts...))

    #create a map of filtered descriptor name to index to guarantee order
    des_ind_dict = Dict{String, Int}()

    for (i, fm) in enumerate(filtered_mesh)
        des_ind_dict[fm]= i
    end

    articles = dfs["basic"][:pmid]

    #create the data-matrix
    disease_occurances = spzeros(length(filtered_mesh), length(articles))

    #Can this process be more efficient using database join/select?
    narticle = 0

    for (i, pmid) in enumerate(articles)

        #not all mesh are of the desired semantic type
        article_filtered_mesh = get_article_mesh_by_concept(dfs, filtered_mesh, pmid)

        #skip if empty
        if isempty(article_filtered_mesh)
            continue
        end

        #otherwise form feature vector for this article
        indices = []
        for d in article_filtered_mesh
            push!(indices, des_ind_dict[d])
        end

        #TO DO: Not sure about the type. Should we choose bool to save space
        # or float to support opperations
        article_dis_feature  = zeros(Int, (length(filtered_mesh),1))
        article_dis_feature[indices] .= 1

        #append to data matrix
        disease_occurances[:, i] = article_dis_feature
        narticle+=1
    end

    println("-------------------------------------------------------------")
    println("Found ", narticle, " articles with valid descriptors")
    println("-------------------------------------------------------------")
    return des_ind_dict, disease_occurances

end

# Retrieve all mesh descriptors associated with the given umls_concept
function filter_mesh_by_concept(mesh2umls_df::DataFrame, umls_concepts...)

    res = filter(row -> row[:concept] in umls_concepts, mesh2umls_df)

    return res[:descriptor]
end

#convert sparse matrix from semantic occurances function into a dataframe ready to be plugged into the apriori function
function occurances_to_itemsets(des_ind_dict, disease_occurances)
    name_dict = sort(collect(des_ind_dict), by=x->x[2])
    col_names = DataArray(String, length(des_ind_dict))
    for i in 1:length(des_ind_dict)
        name = name_dict[i][1]
        col_names[i] = name
    end
    itemsets = DataFrame(Matrix(convert(Array{Int64}, disease_occurances')))
    names!(itemsets, [symbol(col_names[i]) for i in 1:length(col_names)])
    return itemsets
end

# dataframe version of SQL function in pubmed_sql_utils
function get_article_mesh_by_concept(dfs::Dict{String,DataFrame}, filtered_mesh, article::Int)

    filtered_by_pmid = filter(row -> row[:pmid] == article , dfs["mesh_heading"])

    get_uids = filter(row -> row[:uid] in filtered_by_pmid[:desc_uid], dfs["mesh_desc"])

    filtered_by_mesh = filter(row -> row[:name] in filtered_mesh, get_uids)

    return filtered_by_mesh[:name]

end


text:
import TurbulenceConvection
const TC = TurbulenceConvection

const tc_dir = dirname(dirname(pathof(TurbulenceConvection)))
include(joinpath(tc_dir, "driver", "generate_namelist.jl"))
include(joinpath(tc_dir, "driver", "Cases.jl"))
include(joinpath(tc_dir, "driver", "parameter_set.jl"))
include(joinpath(tc_dir, "driver", "dycore.jl"))
include(joinpath(tc_dir, "driver", "main.jl"))
import .NameList

TurbulenceConvection.initialize_io(sim::Simulation1d) = nothing
TurbulenceConvection.io(sim::Simulation1d) = nothing

update_n(sim, tendencies, N::Int) = update_n(sim, tendencies, Val(N))

function update_n(sim, tendencies, ::Val{N}) where {N}
    grid = sim.grid
    TS = sim.TS
    prog = sim.state.prog
    aux = sim.state.aux
    params = (; edmf = sim.edmf, grid = grid, gm = sim.gm, case = sim.Case, TS = TS, aux = aux)
    for i in 1:N
        TC.∑tendencies!(tendencies, prog, params, TS.t)
    end
    return nothing
end

function init_sim(case_name)
    @info "Initializing $case_name for single timestep, with no IO."
    @info "call update_n(sim, tendencies, n) to run update n-times"
    namelist = NameList.default_namelist(case_name)
    namelist["time_stepping"]["t_max"] = namelist["time_stepping"]["dt_max"]
    namelist["stats_io"]["frequency"] = 10.0e10
    namelist["stats_io"]["skip"] = true
    namelist["meta"]["uuid"] = "01"
    sim = Simulation1d(namelist)

    Cases.initialize_profiles(sim.Case, sim.grid, sim.gm, sim.state)
    satadjust(sim.gm, sim.grid, sim.state)

    Cases.initialize_surface(sim.Case, sim.grid, sim.state, sim.param_set)
    Cases.initialize_forcing(sim.Case, sim.grid, sim.state, sim.gm, sim.param_set)
    Cases.initialize_radiation(sim.Case, sim.grid, sim.state, sim.gm, sim.param_set)

    initialize_edmf(sim.edmf, sim.grid, sim.state, sim.Case, sim.gm, sim.TS)
    return sim
end


text:
function area(r::Robot)
    to_Sud =  moves!(r,Sud)
    to_West = moves!(r,West)

    side = Ost
    putmarkers!(r,side)
    while isborder(r,Nord)==false
        move!(r,Nord)
        side=HorizonSide(mod(Int(side)+2, 4))
        putmarkers!(r,side)
    end

    to_sup_side(r,Sud)
    to_sup_side(r,West)

    moves!(r,Ost,to_West)
    moves!(r,Nord,to_Sud)
end
function putmarkers!(r::Robot, side::HorizonSide)
    putmarker!(r)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end
#функция при которой вся линии помечается маркерами(включая первоначальную точку)
function to_sup_side(r::Robot,side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
    end
end
function moves!(r::Robot, side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end
#функция,при которой считаются шаги до границы

function moves!(r::Robot,side::HorizonSide,num_steps::Int)
    for _ in 1:num_steps
        move!(r,side)
    end
end
#функция,при которой робот идет до исходной точки

max_stars_repo_name: ilya-design/mirea-progs

max_issues_repo_name: ilya-design/mirea-progs

max_forks_repo_name: ilya-design/mirea-progs


text:
"""
Mark a coordinate.
img is should be RBG format. If not, try `colorImg = convert.(RGB, img)`
It will mark a given x and y position with a red dot.
change red dot size with 'r'.
"""
function markImg!(img, y0, x0, r)
  x1 = max(x0 - r, 1)
  y1 = max(y0 - r, 1)
  x2 = min(x0 + r, size(img, 1))
  y2 = min(y0 + r, size(img, 2))
  for x in x1:x2
    for y in y1:y2
      if r^2 > (x0 - x)^2 + (y0 - y)^2
        img[x, y] = typeof(img[1])(1, 0, 0)
      end
    end
  end
  img
end

function markImg!(img, y0, x0, z, r) # Maybe mapped array can be used
  x1 = max(x0 - r, 1)
  y1 = max(y0 - r, 1)
  x2 = min(x0 + r, size(img, 1))
  y2 = min(y0 + r, size(img, 2))
  for x in x1:x2
    for y in y1:y2
      if r^2 > (x0 - x)^2 + (y0 - y)^2
        img[x, y, z] = typeof(img[1])(1, 0, 0)
      end
    end
  end
  img
end

function markImg!(img, pos::AbstractVector, rad)
for i	 in 1:length(pos)
markImg!(img, pos[i].location[2], pos[i].location[1	], pos[i].location[3], rad)
  end
end


text:
using JLD
using RealNeuralNetworks
using RealNeuralNetworks.SWCs

data = load("/tmp/point_clouds.jld")
point_cloud = data["point_cloud"]
dbf = data["dbf"]

swc = nodeNetize(point_cloud; dbf=dbf)
SWCs.save(swc, "/tmp/76880.swc")


text:
module Capsule

export make

"""
    Instance(name, value)

Declares a value auto-wireing for the functions of the encapsulated module.
"""
struct Instance{T}
    name::Symbol
    value::T
end

"""
    make(instance_name, target, instances, [privates])

For use inside of `macro`, generate an instance module named `instance_name` with proxies for all
functions defined within a `target` module. Functions that defines arguments that match the `instances`
provided are auto-wired and removed from the arguments in the instance module.

Optionally, a list of `privates` functions from the `target` module can be specified. These functions
will not have proxy generated for them. They will effectively be unaccessible from the instance module.

# Examples
```julia
module MyApi
macro make(instance_name::Symbol, baseurl::String)
    Capsule.make(
        instance_name,
        MyApi,
        [Capsule.Instance(:ctx, ApiContext(baseurl))],
        [MyApi.my_util_fn]
    )
end
end
...

MyApi.@make ma "https://my.service.io"
ma.get_resource(7)
```
"""
function make(
    instance_name::Symbol,
    target::Module,
    instances::Vector{Instance{T}} where T,
    privates::Vector{T} where T <: Function = Vector{Function}()
)
    module_block = quote
        target_name = nameof($target) |> Symbol
        quote
            using $target_name
        end |> eval

        # Setup instance parameters.
        for i ∈ $instances
            quote
                $(i.name) = $(i.value)
            end |> eval
        end

        $import_methods($target, $instances, $privates) |> eval
    end

    return Expr(:toplevel, Expr(:module, true, esc(instance_name), esc(module_block)))
end

function import_methods(
    target::Module,
    instances::Vector{Instance{T}} where T,
    privates::Vector{T} where T <: Function = []
)
    instance_filter = instances .|> pw -> (string(pw.name), string(typeof(pw.value)))

    quote
        target = $target

        for symbol ∈ names(target; all=true)
            if match(r"^[a-z]\w+$", string(symbol)) !== nothing
                for method ∈ eval(Meta.parse("$target.$symbol")) |> methods
                    decl_parts = Base.arg_decl_parts(method)[2]
                    m_name = Symbol(decl_parts[1][2])

                    # We do not import these methods.
                    if m_name ∈ [:eval, :include] ∪ nameof.($privates)
                        continue
                    end

                    # Regular arguments assignment.
                    m_args = decl_parts[2:end] |>
                        argdef_list -> filter(arg -> arg ∉ $instance_filter, argdef_list) .|>
                        argdef -> argdef[2] === "" ? Symbol(argdef[1]) : Meta.parse("$(argdef[1])::$(argdef[2])")
                    m_args_fwd = decl_parts[2:end] .|> p -> Symbol(p[1])

                    # Key-word argument assignment.
                    m_has_kwargs = length(Base.kwarg_decl(method)) > 0
                    m_kwargs = Base.kwarg_decl(method)
                    m_kwargs_fwd = Base.kwarg_decl(method) .|> kw -> Meta.parse("$kw = $kw")

                    # Generating the signature.
                    signature = m_has_kwargs ? quote # With kwargs
                        $m_name($(m_args...); $(m_kwargs...)) = $target.$m_name($(m_args_fwd...); $(m_kwargs_fwd...))
                    end : quote  # Without kwargs
                        $m_name($(m_args...)) = $target.$m_name($(m_args_fwd...))
                    end

                    # Add it to the module.
                    eval(signature)
                end
            end
        end
    end
end

end


text:
using BandedMatrices, LinearAlgebra, LazyArrays, Test

import BandedMatrices: rowstart, rowstop, colstart, colstop,
                       rowlength, collength, diaglength, BandedColumns
import LazyArrays: MemoryLayout, DenseColumnMajor

@testset "Indexing" begin
    @testset "BandedMatrix Indexing" begin
        let
            for n in (10,50), m in (12,50), Al in (0,1,2,30), Au in (0,1,2,30)
                A=brand(n,m,Al,Au)
                kr,jr=3:10,5:12
                @test Matrix(A[kr,jr]) ≈ Matrix(A)[kr,jr]
            end
        end
    end

    @testset "rowstart/rowstop" begin
        A = BandedMatrix(Ones(7, 5), (1, 2))
        # 1.0  1.0  1.0  0.0  0.0
        # 1.0  1.0  1.0  1.0  0.0
        # 0.0  1.0  1.0  1.0  1.0
        # 0.0  0.0  1.0  1.0  1.0
        # 0.0  0.0  0.0  1.0  1.0
        # 0.0  0.0  0.0  0.0  1.0
        # 0.0  0.0  0.0  0.0  0.0
        @test (rowstart(A, 1), rowstop(A, 1), rowlength(A, 1)) == (1, 3, 3)
        @test (rowstart(A, 2), rowstop(A, 2), rowlength(A, 2)) == (1, 4, 4)
        @test (rowstart(A, 3), rowstop(A, 3), rowlength(A, 3)) == (2, 5, 4)
        @test (rowstart(A, 4), rowstop(A, 4), rowlength(A, 4)) == (3, 5, 3)
        @test (rowstart(A, 5), rowstop(A, 5), rowlength(A, 5)) == (4, 5, 2)
        @test (rowstart(A, 6), rowstop(A, 6), rowlength(A, 6)) == (5, 5, 1)
        @test (rowstart(A, 7), rowstop(A, 7), rowlength(A, 7)) == (6, 5, 0) # zero length
        @test (colstart(A, 1), colstop(A, 1), collength(A, 1)) == (1, 2, 2)
        @test (colstart(A, 2), colstop(A, 2), collength(A, 2)) == (1, 3, 3)
        @test (colstart(A, 3), colstop(A, 3), collength(A, 3)) == (1, 4, 4)
        @test (colstart(A, 4), colstop(A, 4), collength(A, 4)) == (2, 5, 4)
        @test (colstart(A, 5), colstop(A, 5), collength(A, 5)) == (3, 6, 4)

        A = BandedMatrix(Ones(3, 6), (1, 2))
        # 1.0  1.0  1.0  0.0  0.0  0.0
        # 1.0  1.0  1.0  1.0  0.0  0.0
        # 0.0  1.0  1.0  1.0  1.0  0.0
        @test (rowstart(A, 1), rowstop(A, 1), rowlength(A, 1)) == (1, 3, 3)
        @test (rowstart(A, 2), rowstop(A, 2), rowlength(A, 2)) == (1, 4, 4)
        @test (rowstart(A, 3), rowstop(A, 3), rowlength(A, 3)) == (2, 5, 4)
        @test (colstart(A, 1), colstop(A, 1), collength(A, 1)) == (1, 2, 2)
        @test (colstart(A, 2), colstop(A, 2), collength(A, 2)) == (1, 3, 3)
        @test (colstart(A, 3), colstop(A, 3), collength(A, 3)) == (1, 3, 3)
        @test (colstart(A, 4), colstop(A, 4), collength(A, 4)) == (2, 3, 2)
        @test (colstart(A, 5), colstop(A, 5), collength(A, 5)) == (3, 3, 1)
        @test (colstart(A, 6), colstop(A, 6), collength(A, 6)) == (4, 3, 0) # zero length

        A = BandedMatrix(Ones(3, 4), (-1, 2))
        # 0.0  1.0  1.0  0.0
        # 0.0  0.0  1.0  1.0
        # 0.0  0.0  0.0  1.0
        @test rowrange(A, 1) == 2:3
        @test rowrange(A, 2) == 3:4
        @test rowrange(A, 3) == 4:4
        @test colrange(A, 1) == 1:0
        @test colrange(A, 2) == 1:1
        @test colrange(A, 3) == 1:2
        @test colrange(A, 4) == 2:3
    end

    @testset "test length of diagonal" begin
        A = BandedMatrix(Ones(4, 6), (2, 3))
        # 4x6 BandedMatrices.BandedMatrix{Float64}:
        #  1.0  1.0  1.0  1.0
        #  1.0  1.0  1.0  1.0  1.0
        #  1.0  1.0  1.0  1.0  1.0  1.0
        #       1.0  1.0  1.0  1.0  1.0
        @test diaglength(A, -2) == 2
        @test diaglength(A, -1) == 3
        @test diaglength(A,  0) == 4
        @test diaglength(A,  1) == 4
        @test diaglength(A,  2) == 4
        @test diaglength(A,  3) == 3

        A = BandedMatrix(Ones(6, 5), (2, 1))
        # 6x5 BandedMatrices.BandedMatrix{Float64}:
        #  1.0  1.0
        #  1.0  1.0  1.0
        #  1.0  1.0  1.0  1.0
        #       1.0  1.0  1.0  1.0
        #            1.0  1.0  1.0
        #                 1.0  1.0
        @test diaglength(A, -2) == 4
        @test diaglength(A, -1) == 5
        @test diaglength(A,  0) == 5
        @test diaglength(A,  1) == 4

        A = BandedMatrix(Ones(4, 4), (1, 1))
        # 4x4 BandedMatrices.BandedMatrix{Float64}:
        #  1.0  1.0
        #  1.0  1.0  1.0
        #       1.0  1.0  1.0
        #            1.0  1.0
        @test diaglength(A, -1) == 3
        @test diaglength(A,  0) == 4
        @test diaglength(A,  1) == 3

        # test with Band type
        @test diaglength(A,  band(1)) == 3
    end

    @testset "_firstdiagrow/_firstdiagcol" begin
        A = BandedMatrix(Ones(6, 5), (2, 1))
        # 6x5 BandedMatrices.BandedMatrix{Float64}:
        #  1.0  1.0
        #  1.0  1.0  1.0
        #  1.0  1.0  1.0  1.0
        #       1.0  1.0  1.0  1.0
        #            1.0  1.0  1.0
        #                 1.0  1.0
        @test BandedMatrices._firstdiagrow(A, 1) == 2
        @test BandedMatrices._firstdiagrow(A, 2) == 2
        @test BandedMatrices._firstdiagrow(A, 3) == 2
        @test BandedMatrices._firstdiagrow(A, 4) == -3
        @test BandedMatrices._firstdiagrow(A, 5) == -3
        @test BandedMatrices._firstdiagrow(A, 6) == -3

        @test BandedMatrices._firstdiagcol(A, 1) == -3
        @test BandedMatrices._firstdiagcol(A, 2) == -3
        @test BandedMatrices._firstdiagcol(A, 3) == 2
        @test BandedMatrices._firstdiagcol(A, 4) == 2
        @test BandedMatrices._firstdiagcol(A, 5) == 2
    end

    @testset "scalar - integer - integer" begin
        a = BandedMatrix(Ones(5, 5), (1, 1))
        # 1.0  1.0  0.0  0.0  0.0
        # 1.0  1.0  1.0  0.0  0.0
        # 0.0  1.0  1.0  1.0  0.0
        # 0.0  0.0  1.0  1.0  1.0

        # in band
        a[1, 1] = 2
        @test a[1, 1] == 2

        # out of band
        @test_throws BandError a[1, 3] = 1
        @test_throws BandError a[3, 1] = 1

        # out of range
        @test_throws BoundsError a[0, 0] = 1
        @test_throws BoundsError a[0, 1] = 1
        @test_throws BoundsError a[1, 0] = 1
        @test_throws BoundsError a[5, 6] = 1
        @test_throws BoundsError a[6, 5] = 1
        @test_throws BoundsError a[6, 6] = 1

        a = BandedMatrix(Ones(5, 5), (-1, 1))
        a[1, 2] = 2
        @test a[1, 2] == 2
        @test_throws BandError a[1, 1] = 1
    end

    @testset "indexing along a column" begin
        @testset "scalar - BandRange/Colon - integer" begin
            a = BandedMatrix(Ones(5, 5), (1, 1))
            # 1.0  1.0  0.0  0.0  0.0
            # 1.0  1.0  1.0  0.0  0.0
            # 0.0  1.0  1.0  1.0  0.0
            # 0.0  0.0  1.0  1.0  1.0
            # 0.0  0.0  0.0  1.0  1.0

            # in band
            a[BandRange, 1] .= 2
            a[BandRange, 2] .= 3
            a[BandRange, 3] .= 4
            a[BandRange, 4] .= 5
            a[BandRange, 5] .= 6
            @test a ==  [2  3  0  0  0;
                         2  3  4  0  0;
                         0  3  4  5  0;
                         0  0  4  5  6;
                         0  0  0  5  6]

            @test a[BandRange, 1] == [2, 2]
            @test a[BandRange, 2] == [3, 3, 3]
            @test a[BandRange, 3] == [4, 4, 4]
            @test a[BandRange, 4] == [5, 5, 5]
            @test a[BandRange, 5] == [6, 6]

            @test_throws BoundsError a[:, 0] .= 1
            @test_throws BandError a[:, 1] .= 1
            @test_throws BoundsError a[BandRange, 0] .= 1
            @test_throws BoundsError a[BandRange, 6] .= 1

            a = BandedMatrix(Ones(3, 5), (-1, 2))
            @test isempty(a[BandRange,1])
            a[BandRange,2] = [1]
            a[BandRange,3] = [2, 2]
            a[BandRange,4] = [3, 3]
            a[BandRange,5] = [4]

            @test a ==   [0 1 2 0 0;
                          0 0 2 3 0;
                          0 0 0 3 4]
        end

        @testset "vector - BandRange/Colon - integer" begin
            a = BandedMatrix(Ones{Int}(5, 7), (2, 1))
            # 5x7 BandedMatrices.BandedMatrix{Float64}:
            #  1.0  1.0    0    0    0    0   0   0
            #  1.0  1.0  1.0    0    0    0   0   0
            #  1.0  1.0  1.0  1.0    0    0   0   0
            #    0  1.0  1.0  1.0  1.0    0   0   0
            #    0    0  1.0  1.0  1.0  1.0   0   0

            # in band
            a[BandRange, 1] = [1,   2,  3]
            a[BandRange, 2] = [4,   5,  6,  7]
            a[BandRange, 3] = [8,   9, 10, 11]
            a[BandRange, 4] = [12, 13, 14]
            a[BandRange, 5] = [15, 16]
            a[BandRange, 6] = [17]

            @test a == [ 1  4   0   0   0   0   0;
                         2  5   8   0   0   0   0;
                         3  6   9  12   0   0   0;
                         0  7  10  13  15   0   0;
                         0  0  11  14  16  17   0]

            @test a[BandRange, 1] == [1,   2,  3]
            @test a[BandRange, 2] == [4,   5,  6,  7]
            @test a[BandRange, 3] == [8,   9, 10, 11]
            @test a[BandRange, 4] == [12, 13, 14]
            @test a[BandRange, 5] == [15, 16]
            @test a[BandRange, 6] == [17]

            @test_throws BoundsError a[:, 0] = [1, 2, 3]
            @test_throws DimensionMismatch a[:, 1] = [1, 2, 3]
            @test_throws BoundsError a[BandRange, 0] = [1, 2, 3]
            @test_throws BoundsError a[BandRange, 8] = [1, 2, 3]
            @test_throws DimensionMismatch a[BandRange, 1] = [1, 2]
        end

        # scalar - range - integer
        let
            a = BandedMatrix(Ones(3, 4), (2, 1))
            # 1.0  1.0  0.0  0.0
            # 1.0  1.0  1.0  0.0
            # 1.0  1.0  1.0  1.0

            # Matrix column span
            a[1:3, 1] .= 1
            a[1:3, 2] .= 2
            a[2:3, 3] .= 3
            a[3:3, 4] .= 4
            @test a == [ 1  2  0  0;
                         1  2  3  0;
                         1  2  3  4]

            # partial span
            a[1:1, 1] .= 3
            a[2:3, 2] .= 4
            a[3:3, 3] .= 5
            @test a == [ 3  2  0  0;
                         1  4  3  0;
                         1  4  5  4]

             a[:, 1] = [1, 1, 1]
             a[:, 2] = [2, 2, 2]
             a[:, 3] = [0, 3, 3]
             a[:, 4] = [0, 0, 4]

             @test a == [ 1  2  0  0;
                          1  2  3  0;
                          1  2  3  4]

            # wrong range input
            @test_throws BoundsError   a[0:1, 1] .= 3
            @test_throws BoundsError   a[1:4, 1] .= 3
            @test_throws BandError a[2:3, 4] .= 3
            @test_throws BoundsError   a[3:4, 4] .= 3
        end

        # vector - range - integer
        let
            a = BandedMatrix(Ones(5, 4), (1, 2))
            # 1.0  1.0  1.0  0.0
            # 1.0  1.0  1.0  1.0
            # 0.0  1.0  1.0  1.0
            # 0.0  0.0  1.0  1.0
            # 0.0  0.0  0.0  1.0
            # Matrix column span
            a[1:2, 1] = 1:2
            a[1:3, 2] = 3:5
            a[1:4, 3] = 6:9
            a[2:5, 4] = 10:13
            @test a == [1  3  6   0;
                        2  4  7  10;
                        0  5  8  11;
                        0  0  9  12;
                        0  0  0  13]

            # partial span
            a[1:1, 1] = [2]
            a[2:3, 2] = 1:2
            a[3:4, 3] = 3:4
            a[4:5, 4] = 3:4
            @test a == [2  3  6   0;
                        2  1  7  10;
                        0  2  3  11;
                        0  0  4   3;
                        0  0  0   4]

            # wrong range input
            @test_throws BoundsError a[1:2, 0] = 1:2
            @test_throws BoundsError a[0:1, 1] = 1:2
            @test_throws BoundsError a[1:8, 1] = 1:8
            @test_throws BoundsError a[2:6, 4] = 2:6
            @test_throws BoundsError a[3:4, 5] = 3:4

            a[:, 1] = [1,1,0,0,0]
            a[:, 2] = [2,2,2,0,0]
            a[:, 3] = [3,3,3,3,0]
            a[:, 4] = [0,4,4,4,4]

            @test a == [1  2  3   0;
                        1  2  3   4;
                        0  2  3   4;
                        0  0  3   4;
                        0  0  0   4]
        end
    end

    @testset "indexing along a row" begin
        # scalar - integer - BandRange/colon
        let
            a = BandedMatrix(Ones(5, 5), (1, 1))
            # 1.0  1.0  0.0  0.0  0.0
            # 1.0  1.0  1.0  0.0  0.0
            # 0.0  1.0  1.0  1.0  0.0
            # 0.0  0.0  1.0  1.0  1.0
            # 0.0  0.0  0.0  1.0  1.0

            # in band
            a[1, BandRange] .= 2
            a[2, BandRange] .= 3
            a[3, BandRange] .= 4
            a[4, BandRange] .= 5
            a[5, BandRange] .= 6

            @test a == [2  2  0  0  0;
                        3  3  3  0  0;
                        0  4  4  4  0;
                        0  0  5  5  5;
                        0  0  0  6  6]

            @test_throws BoundsError a[0, :] .= 1
            @test_throws BandError a[1, :] .= 1
            @test_throws BoundsError a[0, BandRange] .= 1
            @test_throws BoundsError a[6, BandRange] .= 1

            a[1, :] .= 0
            @test vec(a[1, :]) == [0, 0, 0, 0, 0]
        end

        # vector - integer - BandRange/colon
        let
            a = BandedMatrix(Ones(7, 5), (1, 2))
            # 7x5 BandedMatrices.BandedMatrix{Float64}:
            # 1.0  1.0  1.0  0.0  0.0
            # 1.0  1.0  1.0  1.0  0.0
            # 0.0  1.0  1.0  1.0  1.0
            # 0.0  0.0  1.0  1.0  1.0
            # 0.0  0.0  0.0  1.0  1.0
            # 0.0  0.0  0.0  0.0  1.0
            # 0.0  0.0  0.0  0.0  0.0

            # in band
            a[1, BandRange] = [1,   2,  3]
            a[2, BandRange] = [4,   5,  6,  7]
            a[3, BandRange] = [8,   9, 10, 11]
            a[4, BandRange] = [12, 13, 14]
            a[5, BandRange] = [15, 16]
            a[6, BandRange] = [17]

            @test a ==  [1  2   3   0   0;
                         4  5   6   7   0;
                         0  8   9  10  11;
                         0  0  12  13  14;
                         0  0   0  15  16;
                         0  0   0   0  17;
                         0  0   0   0   0]

            @test_throws BoundsError a[0, :] = [1, 2, 3]
            @test_throws DimensionMismatch a[1, :] = [1, 2, 3]
            @test_throws BoundsError a[0, BandRange] = [1, 2, 3]
            @test_throws BoundsError a[8, BandRange] = [1, 2, 3]
            @test_throws DimensionMismatch a[1, BandRange] = [1, 2]
            @test_throws DimensionMismatch a[7, BandRange] = [1, 2, 3]

            a[7, BandRange] .= 1
            @test a ==  [1  2   3   0   0;
                         4  5   6   7   0;
                         0  8   9  10  11;
                         0  0  12  13  14;
                         0  0   0  15  16;
                         0  0   0   0  17;
                         0  0   0   0   0]
        end

        # scalar - integer - range
        let
            a = BandedMatrix(Ones(7, 5), (1, 2))
            # 7x5 BandedMatrices.BandedMatrix{Float64}:
            # 1.0  1.0  1.0  0.0  0.0
            # 1.0  1.0  1.0  1.0  0.0
            # 0.0  1.0  1.0  1.0  1.0
            # 0.0  0.0  1.0  1.0  1.0
            # 0.0  0.0  0.0  1.0  1.0
            # 0.0  0.0  0.0  0.0  1.0
            # 0.0  0.0  0.0  0.0  0.0

            # in band
            a[1, 1:3] .= 2
            a[2, 1:4] .= 3
            a[3, 2:5] .= 4
            a[4, 3:5] .= 5
            a[5, 4:5] .= 6
            a[6, 5:5] .= 7

            @test a ==  [2  2  2  0  0;
                         3  3  3  3  0;
                         0  4  4  4  4;
                         0  0  5  5  5;
                         0  0  0  6  6;
                         0  0  0  0  7;
                         0  0  0  0  0]

            a[1, 1:4] .= 0

            @test a ==  [0  0  0  0  0;
                         3  3  3  3  0;
                         0  4  4  4  4;
                         0  0  5  5  5;
                         0  0  0  6  6;
                         0  0  0  0  7;
                         0  0  0  0  0]

            @test_throws BoundsError a[0, 1:3] .= 1
            @test_throws BoundsError a[8, 2:3] .= 1
            @test_throws BoundsError a[1, 0:3] .= 1
            @test_throws BoundsError a[4, 4:6] .= 1

            @test_throws BandError a[1, 1:4] .= 1
        end

        # vector - integer - range
        let
            a = BandedMatrix(Ones(7, 5), (1, 2))
            # 7x5 BandedMatrices.BandedMatrix{Float64}:
            # 1.0  1.0  1.0  0.0  0.0
            # 1.0  1.0  1.0  1.0  0.0
            # 0.0  1.0  1.0  1.0  1.0
            # 0.0  0.0  1.0  1.0  1.0
            # 0.0  0.0  0.0  1.0  1.0
            # 0.0  0.0  0.0  0.0  1.0
            # 0.0  0.0  0.0  0.0  0.0

            # in band
            a[1, 1:3] = [1, 2, 3]
            a[2, 1:4] = [1, 2, 3, 4]
            a[3, 2:5] = [1, 2, 3, 4]
            a[4, 3:5] = [1, 2, 3]
            a[5, 4:5] = [1, 2]
            a[6, 5:5] = [1]

            @test a ==  [1  2  3  0  0;
                         1  2  3  4  0;
                         0  1  2  3  4;
                         0  0  1  2  3;
                         0  0  0  1  2;
                         0  0  0  0  1;
                         0  0  0  0  0]

            @test_throws BoundsError a[0, 1:3] = [1, 2, 3]
            @test_throws BoundsError a[8, 2:3] = [1, 2]
            @test_throws BoundsError a[1, 0:3] = [1, 2, 3, 4]
            @test_throws BoundsError a[4, 4:6] = [2, 3]

            @test_throws BandError a[1, 1:4] = [1, 2, 3, 4]
            @test_throws BoundsError a[7, 6:7] = [1, 2]
            @test_throws DimensionMismatch a[1, 1:3] = [1, 2]
        end
    end

    @testset "indexing along a band" begin
        let
            a = BandedMatrix(Zeros(5, 4), (2, 1))
            # 5x4 BandedMatrices.BandedMatrix{Float64}:
            #  0.0  0.0
            #  0.0  0.0  0.0
            #  0.0  0.0  0.0  0.0
            #       0.0  0.0  0.0
            #            0.0  0.0
            a[band(-2)] .= 5
            a[band(-1)] .= 1
            a[band( 0)] .= 2
            a[band( 1)] .= 3

            @test Matrix(a) == [ 2  3  0  0;
                               1  2  3  0;
                               5  1  2  3;
                               0  5  1  2;
                               0  0  5  1]

            @test_throws BandError a[band(-3)] .= 1

            a[band(-2)] = [4, 4, 4]
            a[band(-1)] = [1, 2, 3, 4]
            a[band( 0)] = [5, 6, 7, 8]
            a[band( 1)] = [9, 10, 11]

            @test Matrix(a) == [ 5  9   0   0;
                               1  6  10   0;
                               4  2   7  11;
                               0  4   3   8;
                               0  0   4   4]

            @test a[band(-2)] == [4, 4, 4]
            @test a[band(-1)] == [1, 2, 3, 4]
            @test a[band( 0)] == [5, 6, 7, 8]
            @test a[band( 1)] == [9, 10, 11]

            @test_throws BandError a[band(-3)] = [1, 2, 3]
        end

        let
            a = BandedMatrix(Zeros(5, 4), (2, 2))
            # 5x4 BandedMatrices.BandedMatrix{Float64}:
            #  0.0  0.0  0.0
            #  0.0  0.0  0.0  0.0
            #  0.0  0.0  0.0  0.0
            #       0.0  0.0  0.0
            #            0.0  0.0
            a[band(-2)] .= 5
            a[band(-1)] .= 1
            a[band( 0)] .= 2
            a[band( 1)] .= 3
            a[band( 2)] .= 4

            @test Matrix(a) == [ 2  3  4  0;
                               1  2  3  4;
                               5  1  2  3;
                               0  5  1  2;
                               0  0  5  1]

            @test_throws BandError a[band(-3)] .= 1

            a[band(-2)] = [4,  4,  4]
            a[band(-1)] = [1,  2,  3, 4]
            a[band( 0)] = [5,  6,  7, 8]
            a[band( 1)] = [9,  10, 11]
            a[band( 2)] = [12, 13]

            @test Matrix(a) == [ 5  9   12  0;
                               1  6  10  13;
                               4  2   7  11;
                               0  4   3   8;
                               0  0   4   4]

            @test a[band(-2)] == [4, 4, 4]
            @test a[band(-1)] == [1, 2, 3, 4]
            @test a[band( 0)] == [5, 6, 7, 8]
            @test a[band( 1)] == [9, 10, 11]
            @test a[band( 2)] == [12, 13]

            @test_throws BandError a[band(-3)] = [1, 2, 3]
        end
    end

    @testset "other special indexing" begin
        @testset "all elements" begin
            a = BandedMatrix(Ones(3, 3), (1, 1))
            a[:] .= 0
            @test a == [0 0 0;
                        0 0 0;
                        0 0 0]

            # all rows/cols
            a.data .= 2
            @test a == [2 2 0;
                        2 2 2;
                        0 2 2]

            a[:, :] = [3 3 0;
                        3 3 3;
                        0 3 3]
            @test a == [3 3 0;
                        3 3 3;
                        0 3 3]

            @test_throws BandError a[:] .= 1
            @test_throws BandError a[:,:] .= 1
            @test_throws BandError a[:,:] = ones(3,3)
        end

        @testset "replace a block in the band" begin
            a = BandedMatrix(Zeros(5, 4), (2, 1))
            # 5x4 BandedMatrices.BandedMatrix{Float64}:
            #  0.0  0.0
            #  0.0  0.0  0.0
            #  0.0  0.0  0.0  0.0
            #       0.0  0.0  0.0
            #            0.0  0.0
            a[1:3, 1:2] .= 2
            a[3:5, 3:4] = [1 2;
                           3 4;
                           5 6]
            @test Matrix(a) == [2  2  0 0;
                              2  2  0 0;
                              2  2  1 2;
                              0  0  3 4;
                              0  0  5 6]

            @test_throws BoundsError a[0:3, 1:2] .= 2
            @test_throws BoundsError a[3:6, 3:4] .= 2
            @test_throws BoundsError a[1:3, 0:2] .= 2
            @test_throws BoundsError a[3:5, 3:5] .= 2
            @test_throws BoundsError a[0:3, 1:2] = [1 2; 3 4]
            @test_throws BoundsError a[3:6, 3:4] = [1 2; 3 4]
            @test_throws BoundsError a[1:3, 0:2] = [1 2; 3 4]
            @test_throws BoundsError a[3:5, 3:5] = [1 2; 3 4]

            @test_throws BandError a[1:3, 1:3] .= 2
            @test_throws BandError a[1:3, 1:3] = rand(3, 3)

            @test_throws DimensionMismatch a[1:3, 1:2] = rand(3, 3)

            s = BandedMatrix(view(a,1:3,1:2))
            @test isa(s,BandedMatrix)
            @test Matrix(s) == a[1:3,1:2]
        end

        @testset "tests bug" begin
            a = BandedMatrix(Zeros(1,1), (3,-1))
            @test a[band(-2)] == Float64[]
        end

        @testset "band views" begin
            for A in (rand(11,10), brand(11,10,2,3), brand(Float32, 11,10,2,3),
                         brand(ComplexF64, 11,10,2,3))
                for k = -5:5
                    V = view(A, band(k))
                    bs = BandedMatrices.parentindices(V)[1] # a bandslice
                    @test bs.indices == diagind(A, k)
                    @test bs.band == Band(k)
                    @test collect(bs) == collect(diagind(A, k))
                    @test Vector{eltype(A)}(V) == collect(V) == A[diagind(A,k)] == A[band(k)]
                    @test Vector{ComplexF64}(V) == Vector{ComplexF64}(A[diagind(A,k)]) ==
                            convert(AbstractVector{ComplexF64}, V) == convert(AbstractArray{ComplexF64}, V)
                    @test V ≡ convert(AbstractArray, V) ≡ convert(AbstractArray{eltype(A)}, V) ≡
                            convert(AbstractArray, V) ≡ convert(AbstractVector, V)
                end
            end
        end
    end

    @testset "Sub-banded views" begin
        A = brand(10,10,2,1)
        V = view(A,Base.OneTo(5),Base.OneTo(6))
        @test MemoryLayout(V) == BandedColumns(DenseColumnMajor())
        @test A[Base.OneTo(5),Base.OneTo(6)] isa BandedMatrix
        @test A[2:5,Base.OneTo(6)] isa BandedMatrix
        @test A[2:2:5,Base.OneTo(6)] isa Matrix
    end
end


text:
module GlassBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Glass Block" => Ahorn.EntityPlacement(
        Maple.GlassBlock,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::Maple.GlassBlock) = 8, 8
Ahorn.resizable(entity::Maple.GlassBlock) = true, true

Ahorn.selection(entity::Maple.GlassBlock) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.GlassBlock, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (1.0, 1.0, 1.0, 0.5), (1.0, 1.0, 1.0, 0.5))
end

end


text:
#Must define isweak, weaken, strengthen and @lift_gi the datastructure.
import Base.getindex

#Don't weaken arrays or dicts
isweak(a) = false
weaken(a) = nothing
strengthen(a) = nothing
@lift_gi(Vector, 1)
@lift_gi(Matrix, 2)
@lift_gi(Dict, 1)

#If
export If, @If
type If
    f_true::Function
    f_false::Function
    val_true
    val_false
    If(f1::Function, f2::Function) = new(f1, f2, uneval, uneval)
end

type Uneval
end
const uneval = Uneval()

getindex(i::If, cond::Int) = getindex(i, bool(cond))
getindex(i::If, cond::Bool) =
    if cond
        if isa(i.val_true, Uneval)
            i.val_true = i.f_true()
        end
        i.val_true
    else
        if isa(i.val_false, Uneval)
            i.val_false = i.f_false()
        end
        i.val_false
    end

macro If(cond, val_true, val_false)
    esc(:(If(() -> $val_true, () -> $val_false)[$cond]))
end

isweak(i::If) = isa(i.val_true, WeakRef)
weaken(i::If) = begin
    i.val_true = WeakRef(i.val_true)
    i.val_false = WeakRef(i.val_false)
end
strengthen(wr::WeakRef) = (wr == WeakRef()) ? uneval : wr.value
strengthen(i::If) = begin
    i.val_true = strengthen(i.val_true)
    i.val_false = strengthen(i.val_false)
end
@lift_gi(If, 1)

export Mem#, InfiniteVector

#type InfiniteVector{T}
#    vec::Vector{T}
#    defined::Vector{Bool}
#end
#InfiniteVector{T}() = InfiniteVector(Array(T, 0), Array(Bool, 0))
#InfiniteVector() = InfiniteVector(Array(Any, 0), Array(Bool, 0))
#getindex(iv::InfiniteVector, i::Int) = begin
#    if length(iv.vec) < i || !iv.defined[i]
#        error("Index not defined in InfiniteVector")
#    end
#    iv.vec[i]
#end
#setindex!(iv::InfiniteVector, val, i::Int) = begin
#    while length(iv.vec) < i
#        old_vec = iv.vec
#        iv.vec = Array(eltype(iv.vec), length(old_vec))
#        append!(iv.vec, Array(eltype(iv.vec), i - length(iv.vec)))
#        append!(iv.defined, fill(false, i - length(iv.defined)))
#    end
#    iv.defined[i] = true
#    iv.vec[i] = val
#end
#import Base.haskey
#haskey(iv::InfiniteVector, i::Int) = begin
#    @assert length(iv.vec) == length(iv.defined)
#    length(iv.vec) >= i && iv.defined[i]
#end

type Mem
    f::Function
    dict::Dict{Any, WRV}
end
Mem(f::Function) = Mem(f, Dict{Any, WRV}())
getindex(m::Mem, args...) = begin
    if !haskey(m.dict, args)
        m.dict[args] = m.f(args...)
    end
    m.dict[args]
end
#Need to do something here!
#isweak(m::Mem) = isa(first(m.dict)[2], WeakRef)
isweak(m::Mem) = begin
    isa(first(m.dict)[2], WeakRef)
end
weaken(m::Mem) = begin
    for key in keys(m.dict)
        m.dict[key] = WeakRef(m.dict[key])
    end
    nothing
end
strengthen(m::Mem) = begin
    for key in collect(keys(m.dict))
        if m.dict[key] == WeakRef()
            delete!(m.dict, key)
        else
            @assert isa(m.dict[key], WeakRef)
            m.dict[key] = m.dict[key].value
        end
    end
    nothing
end
@lift_gi(Mem, 1)


text:
module NLExample

greet() = print("Hello World!")

include("haircut_file.jl")
# this is a file to contain the haircut code

export therightsum
end # module


text:
"Domain Specific Language for SEMs"
module SEMLang

# export @SEM, interpret, SEMSyntaxError
export SEM
using ..CausalCore: ExogenousVariable, EndogenousVariable

struct SEMSyntaxError <: Exception
  msg
end

SEMSyntaxError() = SEMSyntaxError("")

"Parse exogenous variable `line`"
function parseexo(line)
  new_var = line.args[2]
  dist = line.args[3]
  :($new_var = ExogenousVariable( $(Meta.quot(new_var)), $(dist)))
end

"Parse endogenous variable"
function parseendo(line)
    line
end

"Structural Equation Model"
function SEM(sem)
  if sem.head != :block
    throw(SEMSyntaxError("@SEM expects a block expression as input"))
  end
  semlines = Expr[]
  for line in sem.args
    if typeof(line) == Expr
      if line.head == :(=)
        expr = parseendo(line)
      elseif line.head == :call
        expr = parseexo(line)
      else
        throw(SEMSyntaxError())
      end
      push!(semlines, expr)
    end
  end
  Expr(:block, semlines...)
end

end


text: Dict{String,Any}()

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
if basename(pwd()) == "aoc"
    cd("2016/7")
end

regex = r"\[([^\]]+)\]"

function supportstls(line::AbstractString)
    all(!hasabba(m.match) for m in eachmatch(regex, line)) && hasabba(line)
end

function hasabba(str::AbstractString)
    m = match(r"(.)(.)\2\1", str)
    !isnothing(m) && m[1] != m[2]
end

part1(filename::AbstractString) = count(supportstls, eachline(filename))
@assert part1("example.txt") == 2
part1("input.txt")

function supportssl(line::AbstractString)
    str = replace(line, regex => "")
    any(
        occursin(m[2] * m[1] * m[2], str)
        for section in eachmatch(regex, line)
        for m in eachmatch(r"(.)(.)\1", first(section), overlap = true)
        if m[1] != m[2]
    )
end

part2(filename::AbstractString) = count(supportssl, eachline(filename))
@assert part2("example2.txt") == 3

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
using FluxRM
using FluxUtils

using FluxRM
const flux = Flux()

const kvs = FluxRM.KVS(flux)
FluxRM.transaction(kvs) do txn
    FluxRM.mkdir!(txn, "run_snap")
end

const __job_counter = Ref{Int64}(0)
function job_counter()
    counter = __job_counter[]
    __job_counter[] = counter + 1
    counter
end

function md_job(params)
    counter = job_counter()

    path = "run_snap.job_$counter"

    FluxRM.transaction(kvs) do txn
        FluxRM.mkdir!(txn, path)
        put_raw!(txn, path * ".input", Dict(:data => params))
    end

    script = """
    const kvs = FluxRM.KVS(Flux())
    input = FluxRM.lookup(kvs, $path * ".input")
    params = input[:data]

    out = run_md(params)

    FluxRM.transaction(kvs) do txn
        put_raw!(txn, $path * ".output", Dict(:data => out))
    end
    """

    return FluxUtils.Job(".", "utils.jl", script)
end

tasks = []
for (i, β) in enumerate(β_samples)
    md_params = MDParams(
        vcat([βi for βi in β]...),
        params, starting_file, save_dir*"$(i)/", Tend,
        dt, i, Temp, dT)

    job = md_job(md_params)
    push!(tasks, run(job))
end

# forall(tasks) do task
#     # wait(task)
#     if success(task)
#         out = fetch(tasks)
#         #  something productive here
#     else
#         @warn "task failed" task
#     end
# end

# function get_params!(server)
# end

# function update_params!(server, gradient)
# end

# for iter in interations
#     params = get_params!(server)
#     gradient = Zygote.get_gradients(run_md, params)
#     update_params!(server, gradient)
# end


text:
immutable Point{T} <: Domain{T,0}
    x::T
end

==(a::Point,b::Point)=a.x==b.x

for op in (:*,:+,:-,:.*,:.+,:.-,:.^)
    @eval begin
        $op(c::Number,d::Point)=Point($op(c,d.x))
        $op(d::Point,c::Number)=Point($op(d.x,c))
    end
end

for op in (:/,:./)
    @eval $op(d::Point,c::Number)=Point($op(d.x,c))
end

for op in (:+,:-,:.+,:.-)
    @eval $op(a::Point,b::Point)=Point($op(a.x,b.x))
end

for op in (:*,:+)
    @eval begin
        $op(a::Point,v::Vector)=map(y->$op(a,y),v)
        $op(v::Vector,a::Point)=map(y->$op(y,a),v)
    end
end

isambiguous(d::Point) = isnan(d.x)
Base.convert{PT<:Point}(::Type{PT},::AnyDomain) = PT(NaN)

Base.norm(p::Point)=norm(p.x)

Base.getindex(p::Point,k...)=p.x[k...]

Base.in(x,d::Point) = isapprox(x,d.x)

Base.isnan(d::Point) = false

Base.issubset(a::Point,d::UnionDomain)=a.x in d
Base.issubset(a::Point,b::Domain)=a.x in b

Base.intersect(a::Point,b::Point)=b.x in a?b:EmptyDomain()
Base.intersect(a::UnionDomain,b::Point)=b.x in a?b:EmptyDomain()
Base.intersect(a::Domain,b::Point)=b.x in a?b:EmptyDomain()
Base.intersect(b::Point,a::UnionDomain)=b.x in a?b:EmptyDomain()
Base.intersect(b::Point,a::Domain)=b.x in a?b:EmptyDomain()

Base.setdiff(a::Point,b::Point)=a==b?EmptyDomain():a
Base.reverse(a::Point)=a

canonicaldomain(a::Point)=Point(0.)


text:
"""
[Ceva's theorem](https://en.wikipedia.org/wiki/Ceva%27s_theorem) is a theorem about triangles in
plane geometry. Given a triangle `ABC`, let the lines `AO`, `BO` and `CO` be drawn from the vertices
to a common point `O` (not on one of the sides of `ABC`), to meet opposite sides at `D`, `E` and `F`
respectively.  (The segments `AD`, `BE`, and `CF` are known as
[cevians](https://en.wikipedia.org/wiki/Cevian).) Then, using signed lengths of segments,
```math
BD \\times CE \\times AF = DC \\times EA \\times FB.
```
"""
module Ceva

using Plots
using SymPy
using PlaneGeometry
using PlaneGeometry.GeoPlots

export
    ceva_proof, ceva_draw, ceva_rand, ceva, ceva_check

"""
    ceva_draw(📐, O)

Verify Ceva's Theorem for the triangle `📐`.
"""
function ceva_draw(tri, O)

    plist, elist = ceva(tri, O)
    A, B, C, D, E, F = plist
    push!(plist, O)

    plt = plot(tri, fill=(0, :green), aspect_ratio=:equal, fillalpha= 0.2, leg=false)

    plt = scatter!(shape(plist),
            leg=false,
            aspect_ratio=:equal,
            color=:red,
            series_annotations = text.(['A':'F'..., 'O'], :bottom))

    AD, BE, CF = [Edge(epair...) for epair in zip([A, B, C], [D, E, F])]

    for e in [AD, BE, CF]
        plot!(e, color=:orange, leg=false)
    end

    plt, ceva_check(elist...)
end

"""
    ceva_rand()

Verify Ceva's Theorem for a random triangle.
"""
function ceva_rand()
    pts = rand(0:1//10:1,6);
    A, B, C = [Point(pts[i], pts[i+3]) for i in 1:3]

    # Choose O within ABC so the picture looks nicer
    ptx = [pt.x for pt in [A, B, C]]
    pty = [pt.y for pt in [A, B, C]]
    xmin, xmax = N.([f(ptx...) for f in [min, max]])
    ymin, ymax = N.([f(pty...) for f in [min, max]])

    tri = Triangle(A, B, C)
    AB, BC, CA = edges(tri)

    # Use rejections sampling until the point is in ABC
    while true
        x = rand(xmin:1//100:xmax)
        y = rand(ymin:1//100:ymax)
        O = Point(x, y)

        plist, elist = ceva(tri, O)

        # We are unlucky
        if plist == nothing
            continue
        end

        BD, DC, CE, EA, AF, FB = elist

        if BD > BC || DC > BC || CE > CA || EA > CA || AF > AB || FB > AB
            # This mean O is not in ABC
            continue
        end

        return ceva_draw(Triangle(A, B, C), O)
    end
end

"""
    ceva_proof()

Prove Ceva's Theorem.
"""
function ceva_proof()
    @vars by cx positive=true
    @vars cy ox oy

    A = Point(0, 0); B = Point(0, by); C = Point(cx, cy); O = Point(ox, oy)
    tri = Triangle(A, B, C)

    plist, elist = ceva(tri, O)
    ceva_check(elist...)
end

"""
    ceva(Triangle(A, B, C), O)

Returns `([A, B, C, D, E, F], [BD, DC, CE, EA, AF, FB])`.
"""
function ceva(tri, O)
    A, B, C = vertices(tri)
    AB, BC, CA = edges(tri)
    AO, BO, CO = [Edge(pt, O) for pt in vertices(tri)]
    D, E, F = [concurrent(epair...) for epair in zip([BC, CA, AB], [AO, BO, CO])]

    # This mean no such points can be found. (This is possible)
    if D == nothing || E == nothing || F == nothing
        return nothing, nothing
    end

    elist = [Edge(pts...) for pts in zip([B, D, C, E, A, F], [D, C, E, A, F, B])]
    ([A, B, C, D, E, F], elist)
end

"""
    ceva_check(BD, DC, CE, EA, AF, FB)

Check if the line segments satisfies Ceva's theorem.
"""
function ceva_check(BD, DC, CE, EA, AF, FB)
    eq = Eq((BD*CE*AF)^2, (DC*EA*FB)^2)
    N(simplify(eq))
end

end


text:
using ArgParse

using Random
Random.seed!(0)

using Transformers.Datasets
using Transformers.Datasets: WMT, IWSLT

function parse_commandline()
  s = ArgParseSettings()

  @add_arg_table s begin
    "--gpu", "-g"
    help = "use gpu"
    action = :store_true

    "task"
    help = "task name"
    required = true
    range_tester = x-> x ∈ ["wmt14", "iwslt2016", "copy"]
  end

  return parse_args(ARGS, s)
end

const args = parse_commandline()

enable_gpu(args["gpu"])

if args["task"] == "copy"
  const N = 2
  const V = 10
  const Smooth = 1e-6
  const Batch = 32
  const lr = 1e-4

  startsym = "11"
  endsym = "12"
  unksym = "0"
  labels = [unksym, startsym, endsym, collect(map(string, 1:V))...]

  function gen_data()
    global V
    d = map(string, rand(1:V, 10))
    (d,d)
  end

  function preprocess(data)
      x, t = data
      x = mkline.(x)
      t = mkline.(t)
      x_mask = getmask(x)
      t_mask = getmask(t)
      x, t = vocab(x, t)
     todevice(x,t,x_mask,t_mask)
  end

  function train!()
    global Batch
    println("start training")
    model = (embed=embed, encoder=encoder, decoder=decoder)
    i = 1
    for i = 1:320*7
      data = batched([gen_data() for i = 1:Batch])
      x, t, x_mask, t_mask = preprocess(data)
      grad = gradient(ps) do
          l = loss(model, x, t, x_mask, t_mask)
          l
      end
        i%8 == 0 && @show loss(model, x, t, x_mask, t_mask)
      update!(opt, ps, grad)
    end
  end

  mkline(x) = [startsym, x..., endsym]
elseif args["task"] == "wmt14" || args["task"] == "iwslt2016"
  const N = 6
  const Smooth = 0.4
  const Batch = 8
  const lr = 1e-6
  const MaxLen = 100

  const task = args["task"]

  if task == "wmt14"
    wmt14 = WMT.GoogleWMT()

    datas = dataset(Train, wmt14)
    vocab = get_vocab(wmt14)
  else
    iwslt2016 = IWSLT.IWSLT2016(:en, :de)

    datas = dataset(Train, iwslt2016)
    vocab = get_vocab(iwslt2016)
  end

  startsym = "<s>"
  endsym = "</s>"
  unksym = "</unk>"
  labels = [unksym, startsym, endsym, collect(keys(vocab))...]

  function preprocess(batch)
      x = mkline.(batch[1])
      t = mkline.(batch[2])
      x_mask = getmask(x)
      t_mask = getmask(t)
      x, t = vocab(x, t)
      todevice(x,t,x_mask,t_mask)
  end

  function train!()
    global Batch
    println("start training")
    i = 1
    model = (embed=embed, encoder=encoder, decoder=decoder)
    while (batch = get_batch(datas, Batch)) != []
      x, t, x_mask, t_mask = preprocess(batch)
      grad = gradient(ps) do
          loss(model, x, t, x_mask, t_mask)
      end
      i+=1
      i%8 == 0 && @show loss(model, x, t, x_mask, t_mask)

      @time update!(opt, ps, grad)
    end
  end

  if task == "wmt14"
    function mkline(x)
      global MaxLen
      xi = split(x)
      if length(xi) > MaxLen
        xi = xi[1:100]
      end

      [startsym, xi..., endsym]
    end
  else
    function mkline(x)
      global MaxLen
      xi = tokenize(x)
      if length(xi) > MaxLen
        xi = xi[1:100]
      end

      [startsym, xi..., endsym]
    end
  end
else
  error("task not define")
end


text:
# Automatically generated using Clang.jl

const helics_input = Ptr{Cvoid}
const helics_publication = Ptr{Cvoid}
const helics_endpoint = Ptr{Cvoid}
const helics_filter = Ptr{Cvoid}
const helics_core = Ptr{Cvoid}
const helics_broker = Ptr{Cvoid}
const helics_federate = Ptr{Cvoid}
const helics_federate_info = Ptr{Cvoid}
const helics_query = Ptr{Cvoid}
const helics_message_object = Ptr{Cvoid}
const helics_time = Cdouble
const helics_bool = Cint

@cenum helics_iteration_request::UInt32 begin
    helics_iteration_request_no_iteration = 0
    helics_iteration_request_force_iteration = 1
    helics_iteration_request_iterate_if_needed = 2
end

@cenum helics_iteration_result::UInt32 begin
    helics_iteration_result_next_step = 0
    helics_iteration_result_error = 1
    helics_iteration_result_halted = 2
    helics_iteration_result_iterating = 3
end

@cenum helics_federate_state::UInt32 begin
    helics_state_startup = 0
    helics_state_initialization = 1
    helics_state_execution = 2
    helics_state_finalize = 3
    helics_state_error = 4
    helics_state_pending_init = 5
    helics_state_pending_exec = 6
    helics_state_pending_time = 7
    helics_state_pending_iterative_time = 8
    helics_state_pending_finalize = 9
end

struct helics_complex
    real::Cdouble
    imag::Cdouble
end

mutable struct helics_message
    time::helics_time
    data::Cstring
    length::Int64
    messageID::Int32
    flags::Int16
    original_source::Cstring
    source::Cstring
    dest::Cstring
    original_dest::Cstring
end

struct helics_error
    error_code::Int32
    message::Cstring
end

# Skipping MacroDefinition: HELICS_EXPORT __attribute__ ( ( visibility ( "default" ) ) )
# Skipping MacroDefinition: HELICS_NO_EXPORT __attribute__ ( ( visibility ( "hidden" ) ) )
# Skipping MacroDefinition: HELICS_DEPRECATED __attribute__ ( ( __deprecated__ ) )

const HELICS_DEPRECATED_EXPORT = HELICS_EXPORT
const HELICS_DEPRECATED_NO_EXPORT = HELICS_NO_EXPORT

@cenum helics_data_type::UInt32 begin
    helics_data_type_string = 0
    helics_data_type_double = 1
    helics_data_type_int = 2
    helics_data_type_complex = 3
    helics_data_type_vector = 4
    helics_data_type_complex_vector = 5
    helics_data_type_named_point = 6
    helics_data_type_boolean = 7
    helics_data_type_time = 8
    helics_data_type_raw = 25
    helics_data_type_multi = 33
    helics_data_type_any = 25262
end

const helics_data_type_char = helics_data_type_string

@cenum helics_core_type::UInt32 begin
    helics_core_type_default = 0
    helics_core_type_zmq = 1
    helics_core_type_mpi = 2
    helics_core_type_test = 3
    helics_core_type_interprocess = 4
    helics_core_type_ipc = 5
    helics_core_type_tcp = 6
    helics_core_type_udp = 7
    helics_core_type_zmq_test = 10
    helics_core_type_nng = 9
    helics_core_type_tcp_ss = 11
    helics_core_type_http = 12
    helics_core_type_websocket = 14
    helics_core_type_inproc = 18
    helics_core_type_null = 66
end

@cenum helics_federate_flags::UInt32 begin
    helics_flag_observer = 0
    helics_flag_uninterruptible = 1
    helics_flag_interruptible = 2
    helics_flag_source_only = 4
    helics_flag_only_transmit_on_change = 6
    helics_flag_only_update_on_change = 8
    helics_flag_wait_for_current_time_update = 10
    helics_flag_restrictive_time_policy = 11
    helics_flag_rollback = 12
    helics_flag_forward_compute = 14
    helics_flag_realtime = 16
    helics_flag_single_thread_federate = 27
    helics_flag_slow_responding = 29
    helics_flag_delay_init_entry = 45
    helics_flag_enable_init_entry = 47
    helics_flag_ignore_time_mismatch_warnings = 67
    helics_flag_terminate_on_error = 72
end

@cenum helics_log_levels::Int32 begin
    helics_log_level_no_print = -1
    helics_log_level_error = 0
    helics_log_level_warning = 1
    helics_log_level_summary = 2
    helics_log_level_connections = 3
    helics_log_level_interfaces = 4
    helics_log_level_timing = 5
    helics_log_level_data = 6
    helics_log_level_trace = 7
end

@cenum helics_error_types::Int32 begin
    helics_error_fatal = -404
    helics_error_external_type = -203
    helics_error_other = -101
    helics_error_insufficient_space = -18
    helics_error_execution_failure = -14
    helics_error_invalid_function_call = -10
    helics_error_invalid_state_transition = -9
    helics_warning = -8
    helics_error_system_failure = -6
    helics_error_discard = -5
    helics_error_invalid_argument = -4
    helics_error_invalid_object = -3
    helics_error_connection_failure = -2
    helics_error_registration_failure = -1
    helics_ok = 0
end

@cenum helics_properties::UInt32 begin
    helics_property_time_delta = 137
    helics_property_time_period = 140
    helics_property_time_offset = 141
    helics_property_time_rt_lag = 143
    helics_property_time_rt_lead = 144
    helics_property_time_rt_tolerance = 145
    helics_property_time_input_delay = 148
    helics_property_time_output_delay = 150
    helics_property_int_max_iterations = 259
    helics_property_int_log_level = 271
    helics_property_int_file_log_level = 272
    helics_property_int_console_log_level = 274
end

@cenum helics_multi_input_mode::UInt32 begin
    helics_multi_input_no_op = 0
    helics_multi_input_vectorize_operation = 1
    helics_multi_input_and_operation = 2
    helics_multi_input_or_operation = 3
    helics_multi_input_sum_operation = 4
    helics_multi_input_diff_operation = 5
    helics_multi_input_max_operation = 6
    helics_multi_input_min_operation = 7
    helics_multi_input_average_operation = 8
end

@cenum helics_handle_options::UInt32 begin
    helics_handle_option_connection_required = 397
    helics_handle_option_connection_optional = 402
    helics_handle_option_single_connection_only = 407
    helics_handle_option_multiple_connections_allowed = 409
    helics_handle_option_buffer_data = 411
    helics_handle_option_strict_type_checking = 414
    helics_handle_option_ignore_unit_mismatch = 447
    helics_handle_option_only_transmit_on_change = 452
    helics_handle_option_only_update_on_change = 454
    helics_handle_option_ignore_interrupts = 475
    helics_handle_option_multi_input_handling_method = 507
    helics_handle_option_input_priority_location = 510
    helics_handle_option_clear_priority_list = 512
    helics_handle_option_connections = 522
end

@cenum helics_filter_type::UInt32 begin
    helics_filter_type_custom = 0
    helics_filter_type_delay = 1
    helics_filter_type_random_delay = 2
    helics_filter_type_random_drop = 3
    helics_filter_type_reroute = 4
    helics_filter_type_clone = 5
    helics_filter_type_firewall = 6
end


text:
## (Conditional) Probability Queries
## =================================
function prob(Y::SymbolicRandVar{Bool};
              partition_alg::Type{BFSPartition} = BFSPartition,
              args...)
  # init_box = unit_box(LazyBox{Float64}, dims(Y))
  partition = pre_partition(Y, partition_alg; args...)
  res = measure(partition)
  AbstractDomains.Interval(max(0.0, res.l), min(1.0, res.u))
end

function prob_conj(XY::SymbolicRandVar{Bool},
              Y::RandVar{Bool};
              partition_alg::Type{BFSPartition} = BFSPartition,
              args...)
  # init_box = unit_box(LazyBox{Float64}, dims(Y))
  Y_partition = pre_partition(Y, partition_alg; args...)
  Y_measure = measure(Y_partition)
  isequal(Y_measure, zero(Y_measure)) && error("Cannot condition on unsatisfiable events")

  XY_partition = pre_partition(XY, partition_alg; args...)
  res = measure(XY_partition) / Y_measure
  AbstractDomains.Interval(max(0.0, res.l), min(1.0, res.u))
end

# "Lower and upper bounds for marginal probability that 'Y' is true"
# function prob(Y::SymbolicRandVar{Bool}; RandVarType = default_randvar(), args...)
#   executable_Y = convert(RandVarType{Bool}, Y)
#   prob(executable_Y; args...)
# end

"Lower/upper bounds for conditional probability that `X` is true given `Y` is true"
function prob(X::SymbolicRandVar{Bool},
              Y::SymbolicRandVar{Bool};
              args...)
  # executable_Y = convert(RandVarType{Bool}, Y)
  # # We do & here because after conversion & may not be defined on
  # # compiled RandVar types
  # executable_XY = convert(RandVarType{Bool}, X & Y)
  prob_conj(X & Y, Y; args...)
end


text:
using Test, YaoBlocks, YaoArrayRegister
using YaoBlocks.ConstGate
using StatsBase: mean
using LinearAlgebra: tr

@testset "blockfilter expect" begin
    ghz = (ArrayReg(bit"0000") + ArrayReg(bit"1111")) |> normalize!
    obs1 = kron(4, 2=>Z)
    obs2 = kron(4, 2=>X)
    obs3 = repeat(4, X)
    @test expect(obs1, ghz) ≈ 0
    @test expect(obs2, ghz) ≈ 0
    @test expect(obs3, ghz) ≈ 1

    @test expect(obs1+obs2+obs3, ghz) ≈ 1
    @test expect(obs1+obs2+obs3, repeat(ghz, 3)) ≈ [1,1,1]
    @test expect(2*obs3, ghz) ≈ 2
    @test expect(2*obs3, repeat(ghz, 3)) ≈ [2,2,2]

    @test blockfilter(ishermitian, chain(2, kron(2, X, P0), repeat(2, Rx(0), (1,2)), kron(2, 2=>Rz(0.3)))) == Any[X, P0, kron(2, X, P0), Rx(0), repeat(2, Rx(0), (1,2))]
    @test blockfilter(b->ishermitian(b) && b isa PrimitiveBlock, chain(2, kron(2, X, P0), repeat(2, Rx(0), (1,2)), kron(2, 2=>Rz(0.3)))) == [X, P0, Rx(0)]
end

@testset "density matrix" begin
    reg = rand_state(4)
    dm = reg |> density_matrix
    op = put(4, 3=>X)
    @test expect(op, dm) ≈ expect(op, reg)

    reg = rand_state(4; nbatch=3)
    dm = reg |> density_matrix
    op = put(4, 3=>X)
    @test expect(op, dm) ≈ expect(op, reg)
end

@testset "insert_qubits!" begin
    reg = rand_state(5; nbatch=10)
    insert_qubits!(reg, 3; nqubits=2)
    @test reg |> nqubits == 7
    @test expect(put(7, 3=>Z), reg) .|> tr |> mean ≈ 1
    @test expect(put(7, 4=>Z), reg) .|> tr |> mean ≈ 1
end

@testset "expect" begin
    reg = rand_state(3,nbatch=10)
    e1 = expect(put(2, 2=>X), reg |> copy |> focus!(1,2) |> ρ)
    e2 = expect(put(2, 2=>X), reg |> copy |> focus!(1,2))
    e3 = expect(put(3, 2=>X), reg |> ρ)
    e4 = expect(put(3, 2=>X), reg)
    e5 = expect(put(3, 2=>-X), reg)
    @test e1 ≈ e2
    @test e1 ≈ e3
    @test e1 ≈ e4
    @test e5 ≈ -e4
end

@testset "viewbatch apply" begin
    reg = zero_state(1, nbatch=2)
    viewbatch(reg, 1) |> X
    @test reg.state[:,1] == [0, 1]
    @test reg.state[:,2] == [1, 0]
end


text:
using Cropbox
using LeafGasExchange
using Test

"Kim et al. (2007), Kim et al. (2006)"
ge_maize1 = :C4 => (
    Vpm25 = 70, Vcm25 = 50, Jm25 = 300,
    Rd25 = 2
)

"In von Cammerer (2000)"
ge_maize2 = :C4 => (
    Vpm25 = 120, Vcm25 = 60, Jm25 = 400,
)

"In Kim et al.(2006), under elevated CO2, YY"
ge_maize3 = :C4 => (
    Vpm25 = 91.9, Vcm25 = 71.6, Jm25 = 354.2,
    Rd25 = 2, # Values in Kim (2006) are for 31C, and the values here are normalized for 25C. SK
)

"""
switchgrass params from Albaugha et al. (2014)
https://doi.org/10.1016/j.agrformet.2014.02.013
"""
ge_switchgrass1 = :C4 => (Vpm25 = 52, Vcm25 = 26, Jm25 = 145)

"""
switchgrass Vcmax from Le et al. (2010),
others multiplied from Vcmax (x2, x5.5)
"""
ge_switchgrass2 = :C4 => (Vpm25 = 96, Vcm25 = 48, Jm25 = 264)

ge_switchgrass3 = :C4 => (Vpm25 = 100, Vcm25 = 50, Jm25 = 200)
ge_switchgrass4 = :C4 => (Vpm25 = 70, Vcm25 = 50, Jm25 = 180.8)

"switchgrass params from Albaugha et al. (2014)"
ge_switchgrass_base = :C4 => (
    Rd25 = 3.6, # not sure if it was normalized to 25 C
    θ = 0.79,
)

"In Sinclair and Horie, Crop Sciences, 1989"
ge_ndep1 = :NitrogenDependence => (s = 4, N0 = 0.2)

"In J Vos et al. Field Crop Research, 2005"
ge_ndep2 = :NitrogenDependence => (s = 2.9, N0 = 0.25)

"In Lindquist, Weed Science, 2001"
ge_ndep3 = :NitrogenDependence => (s = 3.689, N0 = 0.5)

"""
in P. J. Sellers, et al.Science 275, 502 (1997)
g0 is b, of which the value for c4 plant is 0.04
and g1 is m, of which the value for c4 plant is about 4 YY
"""
ge_stomata1 = :StomataBallBerry => (g0 = 0.04, g1 = 4.0)

"""
Ball-Berry model parameters from Miner and Bauerle 2017,
used to be 0.04 and 4.0, respectively (2018-09-04: KDY)
"""
ge_stomata2 = :StomataBallBerry => (g0 = 0.017, g1 = 4.53)

"calibrated above for our switchgrass dataset"
ge_stomata3 = :StomataBallBerry => (g0 = 0.04, g1 = 1.89)
ge_stomata4 = :StomataBallBerry => (g0 = 0.02, g1 = 2.0)

"parameters from Le et. al (2010)"
ge_stomata5 = :StomataBallBerry => (g0 = 0.008, g1 = 8.0)

"for garlic"
ge_stomata6 = :StomataBallBerry => (g0 = 0.0096, g1 = 6.824)

ge_water1 = :StomataTuzet => (
    sf = 2.3, # sensitivity parameter Tuzet et al. 2003 Yang
    ϕf = -1.2, # reference potential Tuzet et al. 2003 Yang
)

"switchgrass params from Le et al. (2010)"
ge_water2 = :StomataTuzet => (
    #? = -1.68, # minimum sustainable leaf water potential (Albaugha 2014)
    sf = 6.5,
    ϕf = -1.3,
)

"""
August-Roche-Magnus formula gives slightly different parameters
https://en.wikipedia.org/wiki/Clausius–Clapeyron_relation
"""
ge_vaporpressure1 = :VaporPressure => (
    a = 0.61094, # kPa
    b = 17.625, # C
    c = 243.04, # C
)

ge_weather = :Weather => (
    PFD = 1500,
    CO2 = 400,
    RH = 60,
    T_air = 30,
    wind = 2.0,
)

ge_spad = :Nitrogen => (
    _a = 0.0004,
    _b = 0.0120,
    _c = 0,
    SPAD = 60,
)

ge_water = :StomataTuzet => (
    #WP_leaf = 0,
    sf = 2.3,
    Ψf = -1.2,
)

ge_base = (ge_weather, ge_spad, ge_water)

#HACK: zero CO2 prevents convergence of bisection method
ge_step_c = :Weather => :CO2 => 10:10:1500
ge_step_q = :Weather => :PFD => 0:20:2000
ge_step_t = :Weather => :T_air => -10:1:50

@testset "gasexchange" begin
    @testset "C3" begin
        @testset "A-Ci" begin
            Cropbox.visualize(LeafGasExchange.ModelC3MD, :Ci, [:A_net, :Ac, :Aj, :Ap]; config=ge_base, xstep=ge_step_c) |> println
        end

        @testset "A-Q" begin
            Cropbox.visualize(LeafGasExchange.ModelC3MD, :PFD, [:A_net, :Ac, :Aj, :Ap]; config=ge_base, xstep=ge_step_q) |> println
        end

        @testset "A-T" begin
            Cropbox.visualize(LeafGasExchange.ModelC3MD, :T_air, [:A_net, :Ac, :Aj, :Ap]; config=ge_base, xstep=ge_step_t) |> println
        end
    end

    @testset "C4" begin
        @testset "A-Ci" begin
            Cropbox.visualize(LeafGasExchange.ModelC4MD, :Ci, [:A_net, :Ac, :Aj]; config=ge_base, xstep=ge_step_c) |> println
        end

        @testset "A-Q" begin
            Cropbox.visualize(LeafGasExchange.ModelC4MD, :PFD, [:A_net, :Ac, :Aj]; config=ge_base, xstep=ge_step_q) |> println
        end

        @testset "A-T" begin
            Cropbox.visualize(LeafGasExchange.ModelC4MD, :T_air, [:A_net, :Ac, :Aj]; config=ge_base, xstep=ge_step_t) |> println
        end
    end

    @testset "N vs Ψv" begin
        Cropbox.visualize(LeafGasExchange.ModelC4MD, :N, :Ψv, :A_net;
            config=ge_base,
            kind=:heatmap,
            xstep=:Nitrogen=>:N=>0:0.05:2,
            ystep=:StomataTuzet=>:WP_leaf=>-2:0.05:0,
        ) |> println
    end
end


text:
using TimerOutputs
using LinearAlgebra
include("settings.jl")
#This file defines buffer structs for passing theta-dependent values, as well as funtions to initialize them

"""
Buffer of θ-dependent parameters
"""
struct θ_params{O<:Array{Float64, 2}, C<:Cholesky{Float64,Array{Float64, 2}}}
    Eθ::O
    Σθ::O
    Bθ::O
    Dθ::O
    Hθ::O
    Cθ::O
    Σθ_inv_X::O
    choleskyΣθ::C
    choleskyXΣX::C
end

"""
Buffer of θ-dependent parameters, including derivatives
"""
struct θ_param_derivs{O<:Array{Float64, 2}, C<:Cholesky{Float64,Array{Float64, 2}}}
    Eθ::O
    Eθ_prime::O
    Eθ_prime2::O
    Σθ::O
    Σθ_prime::O
    Σθ_prime2::O
    Bθ::O
    Bθ_prime::O
    Bθ_prime2::O
    Dθ::O
    Dθ_prime::O
    Dθ_prime2::O
    Hθ::O
    Hθ_prime::O
    Hθ_prime2::O
    Cθ::O
    Cθ_prime::O
    Cθ_prime2::O
    Σθ_inv_X::O
    tripleΣ::O #X'Sigma^-1 dSigma Sigma^-1 X
    choleskyΣθ::C
    choleskyXΣX::C
end

"""
Compute theta-dependent quantities
Return a struct of type θ_param_derivs if type is \"Turan\" and a struct of type θ_params if type is \"Gaussian\"
"""
function funcθ(θ::Float64, train::trainingData{A, B}, test::testingData{A}, type = "Gaussian") where A<:Array{Float64, 2} where B<:Array{Float64, 1} #::Union{θ_params{W, C}, θ_param_derivs{W, C}} where {W<:Array{Float64, 2}, C<:Cholesky{Float64,Array{Float64, 2}}}
    s = train.s
    s0 = test.s0
    X = train.X
    X0 = test.X0
    z = train.z

    #zeroth order expressions and easily computed higher derivatives (building blocks)
    #@timeit "zero order" begin
        Eθ = fastK(s0, s0, θ, rbf_single)
        Σθ =  fastK(s, s, θ, rbf_single)
        Bθ =  fastK(s0, s, θ, rbf_single)
        @time choleskyΣθ = cholesky(Σθ)
        choleskyXΣX = cholesky(Hermitian(X'*(choleskyΣθ\X)))
        Dθ = Eθ - Bθ*(choleskyΣθ\Bθ')
        Hθ = X0 - Bθ*(choleskyΣθ\X)
        Cθ = Dθ + Hθ*(choleskyXΣX\Hθ')
        Σθ_inv_X = choleskyΣθ\X

    if type == "Turan"
        Eθ_prime = K(s0, s0, θ, rbf_prime)
        Eθ_prime2 = K(s0, s0, θ, rbf_prime2)
        Σθ_prime = fastK(s, s, θ, rbf_prime_single)
        Σθ_prime2 = fastK(s, s, θ, rbf_prime2_single)
        Bθ_prime = K(s0, s, θ, rbf_prime)
        Bθ_prime2 = K(s0, s, θ, rbf_prime2)

    #abstractions used to compute higher derivatives
    dQ = Y -> (choleskyΣθ\((Σθ_prime2*(choleskyΣθ\Y)))) - 2*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\(Y))))))
    Q = Y -> (choleskyΣθ\(Σθ_prime*(choleskyΣθ\(Y))))
    dPinv = Y -> choleskyXΣX\(X'*(choleskyΣθ\(Σθ_prime*(choleskyΣθ\(X*(choleskyXΣX\Y))))))
    XQX = X'*Q(X)
    dP = Y -> -XQX*Y
    d2P = Y -> -X'*(dQ(X)*Y)
    d2Pinv = Y -> -dPinv(dP(choleskyXΣX\Y)) - (choleskyXΣX\(d2P(choleskyXΣX\Y))) - (choleskyXΣX\(dP(dPinv(Y))))
    tripleΣ = Σθ_inv_X' * Σθ_prime * Σθ_inv_X

    #higher derivatives
    @timeit "higher derivs" begin
        Hθ_prime = compute_Hθ_prime(Bθ_prime, Σθ_inv_X, choleskyΣθ, Σθ_prime, Bθ, X)
        Dθ_prime = compute_Dθ_prime(choleskyΣθ, Bθ, Eθ_prime, Σθ_prime, Bθ_prime)
        Cθ_prime = compute_Cθ_prime(Dθ_prime,Hθ,  Hθ_prime, choleskyXΣX, Σθ_inv_X, Σθ_prime)
        Hθ_prime2 = compute_Hθ_prime2(Bθ, Bθ_prime, Bθ_prime2, Σθ_inv_X, choleskyΣθ, X, Q, dQ)
        Dθ_prime2 = compute_Dθ_prime2(choleskyΣθ, Bθ, Bθ_prime, Bθ_prime2, Eθ_prime2, Q, dQ)
        Cθ_prime2 = compute_Cθ_prime2(Dθ_prime2, Hθ, Hθ_prime, Hθ_prime2, choleskyXΣX, dPinv, d2Pinv)
    end
        return θ_param_derivs(Eθ, Eθ_prime, Eθ_prime2, Σθ, Σθ_prime, Σθ_prime2, Bθ, Bθ_prime, Bθ_prime2, Dθ, Dθ_prime, Dθ_prime2, Hθ, Hθ_prime, Hθ_prime2, Cθ, Cθ_prime, Cθ_prime2, Σθ_inv_X, tripleΣ, choleskyΣθ, choleskyXΣX)

    elseif type == "Gaussian"
        return θ_params(Eθ, Σθ, Bθ, Dθ, Hθ, Cθ, Σθ_inv_X, choleskyΣθ, choleskyXΣX)
    else
        throw(ArgumentError("Quadrature type undefined. Please use \"Gaussian\" or \"Turan\" for last arg"))
    end
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
"""
Pretty-prints a `:lambda` expression.
"""
function print_lambda(expr)
    @assert expr.head == :lambda
    arg_symbols = expr.args[1]
    local_vars = expr.args[2][1]
    env_vars = expr.args[2][2]
    gen_syms = expr.args[2][3]
    param_symbols = expr.args[2][4]
    tree = expr.args[3]

    println("Lambda function:")
    println("\t arguments: $arg_symbols")
    println("\t local variables: ")
    println(join(map(o -> "\t\t$((o...,))", local_vars), "\n"))
    println("\t closure variables: ")
    println(join(map(o -> "\t\t$((o...,))", env_vars), "\n"))
    println("\t SSA types: $gen_syms")
    println("\t parameters: $param_symbols")
    println(strip_box(strip_lineno(tree)))
end

"""
Returns a typed AST for function `f`.
"""
function get_typed_ast(f::Function, argtype::Type{Tuple})
    (ast, ret_type) = Core.Inference.typeinf(f.code, argtype, f.env)
    ast_typed = ccall(:jl_uncompress_ast, Any, (Any,Any), f.code, ast)
end

function get_ast(linfo::LambdaStaticData, ast=linfo.ast)
    if !isa(ast, Expr)
        ast = ccall(:jl_uncompress_ast, Any, (Any,Any), linfo, ast)
    end
    ast::Expr
end

"""
Performs type inference on a lambda AST.
"""
function typeinf(linfo::LambdaStaticData, argtypes)
    (ast, ret_type) = Core.Inference.typeinf_uncached(linfo, Tuple{argtypes...}, Base.svec())
    get_ast(linfo, ast)
end

"""
Sets the types of closed vars to those found in the given Dict.
"""
function set_closure_types!(linfo::LambdaStaticData,
                            vartypes::Dict{Symbol, DataType})
    # Make sure to uncompress the AST:
    linfo.ast = get_ast(linfo)

    # Set the types of outer variables:
    closures = linfo.ast.args[2][2]
    linfo.ast.args[2][2] = Any[Any[symbol, get(vartypes, symbol, typ), kind]
                               for (symbol, typ, kind) in closures]
    linfo
end

"""
Calls `func()` on every node in `expr`.
"""
function walk(func, expr)
    func(expr)
    if isa(expr, Expr)
        func(expr.head)
        for o in expr.args
            walk(func, o)
        end
    end
end

"""
Traverses `expr` and replaces nodes according to the
substitutions in dict/function `subst`.
"""

function substitute(subst::Dict, obj::Any)
    get(subst, obj, obj)
end

function substitute(subst::Dict, expr::Expr)
    if haskey(subst, expr)
        return subst[expr]
    else
        head = substitute(subst, expr.head)
        args = [substitute(subst, a) for a in expr.args]
        return Expr(head, args...)
    end
end

function substitute(subst::Function, obj::Any)
    result = subst(obj)
    if result != nothing
        return result
    else
        return obj
    end
end

function substitute(subst::Function, obj::LambdaStaticData)
    new_obj = deepcopy(obj)
    new_obj.ast = substitute(subst, get_ast(new_obj))
    new_obj
end

function substitute(subst::Function, expr::Expr)
    result = subst(expr)
    if result != nothing
        return result
    else
        head = substitute(subst, expr.head)
        args = [substitute(subst, a) for a in expr.args]
        return Expr(head, args...)
    end
end

function sexpr(obj)
    io = IOBuffer()
    Base.Meta.show_sexpr(io, obj)
    takebuf_string(io)
end

function rename_fdef(func)
    funcname = func.args[1].args[1]
    argspec = func.args[1].args[2:end]
    body = func.args[2]
    if isa(funcname, Expr) && funcname.head == :curly
        newsym = gensym(funcname.args[1])
        newname = Expr(:curly, newsym, funcname.args[2:end]...)
    elseif isa(funcname, Symbol)
        newsym = gensym(funcname)
        newname = newsym
    else
        error("Couldn't parse function name: $funcname")
    end
    stageddef = Expr(func.head, Expr(:call, newname, argspec...), body)
#    stagedcall = Expr(:call, newsym, argspec...)
    return newsym, stageddef
end

function get_argname(e::Expr)
    if e.head == :(::)
        return e.args[1]
    elseif e.head == :kw
        return e.args[1]
    end
    error("argname($(sexpr(e))")
end
get_argname(s::Symbol) = s


text: Dict{String, Any}("oct1" => Dict{String, Any}("value" => "342391", "type" => "integer"), "oct2" => Dict{String, Any}("value" => "493", "type" => "integer"), "hex3" => Dict{String, Any}("value" => "3735928559", "type" => "integer"), "hex4" => Dict{String, Any}("value" => "2439", "type" => "integer"), "bin1" => Dict{String, Any}("value" => "214", "type" => "integer"), "bin2" => Dict{String, Any}("value" => "5", "type" => "integer"), "hex2" => Dict{String, Any}("value" => "3735928559", "type" => "integer"), "hex1" => Dict{String, Any}("value" => "3735928559", "type" => "integer"), "oct3" => Dict{String, Any}("value" => "501", "type" => "integer"))


text:
module NLSolversTools

using NLSolversBase
include("objective_types.jl")

end # module