text:
using SBP_GPU
using CUDA
using Test
using SafeTestsets

@testset "SBP_GPU.jl" begin
    # Write your tests here.
    @test π ≈ 3.14 atol=0.01
    Nx = Ny = 10
    y_in = randn(Nx*Ny)
    y_in_GPU = CuArray(y_in)
    y_out_GPU = similar(y_in_GPU)
    y_out_GPU = D2x_GPU(y_in_GPU,y_out_GPU,Nx,Ny)
    @test Array(y_out_GPU) ≈ D2x(y_in,Nx,Ny,1/Nx)
    @test Array(Dx_GPU(y_in_GPU,y_out_GPU,Nx,Ny)) ≈ Dx(y_in,Nx,Ny,1/Nx)

    Nx = Ny = 100
    y_in = randn(Nx*Ny)
    y_in_GPU = CuArray(y_in)
    y_out_GPU = similar(y_in_GPU)
    y_out_GPU = D2x_GPU(y_in_GPU,y_out_GPU,Nx,Ny)
    @test Array(y_out_GPU) ≈ D2x(y_in,Nx,Ny,1/Nx)
    @test Array(Dx_GPU(y_in_GPU,y_out_GPU,Nx,Ny)) ≈ Dx(y_in,Nx,Ny,1/Nx)

    Nx = Ny = 1000
    y_in = randn(Nx*Ny)
    y_in_GPU = CuArray(y_in)
    y_out_GPU = similar(y_in_GPU)
    y_out_GPU = D2x_GPU(y_in_GPU,y_out_GPU,Nx,Ny)
    @test Array(y_out_GPU) ≈ D2x(y_in,Nx,Ny,1/Nx)
    @test Array(Dx_GPU(y_in_GPU,y_out_GPU,Nx,Ny)) ≈ Dx(y_in,Nx,Ny,1/Nx)
end


text:
using Knet, Optim, CUDArt, Base.LinAlg

# v is our estimated unit gradient vector, last successful step.
# "step" is a step-size/scale that we adjust.  The mean step is going to be step*v.
# we grow and shrink s using the factor "grow".
# r is a random unit vector for exploration.
# The scale for r should be proportional to the scale for v.
# Call that fixed proportionality constant rstep, i.e. step*(v+rstep*r).
# We also have rmin, i.e. when s shrinks too much, we want the exploration vector to keep a minimum size.
# step*v + max(step*rstep,rmin)*r

# We assume fmin=0 here and stop when f(x) < ftol

function dhc2{T}(f::Function, x::CudaVector{T}; grow=2.0, rmin=1e-3, rstep=1.0, ftol=1e-4)
    v = randn!(similar(x))
    r = randn!(similar(x))
    x0 = copy(x)
    f0 = f(x)
    nf = np = nfail = 1
    avgstep = 0
    while f0 > ftol
        axpy!(1, v, x)
        # we want vecnorm(r)/vecnorm(v) = rstep
        # vecnorm(randn!(r)) = sqrt(length(r))
        # scale!(rstep*vecnorm(v)/sqrt(length(r)), randn!(r))
        rscale = max(rmin, rstep*vecnorm(v))/sqrt(length(r))
        axpy!(rscale, randn!(r), x)
        f1 = f(x); nf += 1
        if f1 < f0
            nfail = 0
            copy!(x0, x)
            f0 = f1
            axpy!(rscale, r, v)
            avgstep = (avgstep==0 ? vecnorm(v) : 0.99*avgstep + 0.01*vecnorm(v))
            scale!(grow, v)
        else
            nfail += 1
            copy!(x, x0)
            scale!(1/grow, v)
        end
        nf >= np && (println((nf,f0,avgstep,x[1],test02(mnist_f, mnist_tst, zeroone) )); np *= 2)
    end
    (nf,f0,avgstep,x)
end

# Set up mnist as an objective function

using Knet
using Knet: Net, params, axpy!
if !isdefined(:MNIST)
    include(Pkg.dir("Knet/examples/mnist.jl"))
end
batchsize = 1000
mnist_trn = minibatch(MNIST.xtrn, MNIST.ytrn, 60000)
mnist_tst = minibatch(MNIST.xtst, MNIST.ytst, 10000)

@knet function model02(x; winit=Gaussian(0,.1), hidden=32, f1=:relu, o...)
    h    = wf(x; out=hidden, f=f1, winit=winit)
    return wf(h; out=10, f=:soft, winit=winit)
end

# We'll try to create the two weight matrices from a single array

function initmodel{T}(w::CudaArray{T}; X=784, H=32, Y=10, B=100)
    @assert length(w) == X*H+H*Y
    f = compile(:model02; hidden=H)
    forw(f, zeros(T, X, B))
    f.reg[2].out0 = CudaArray(w.ptr, (H,X), w.dev)
    f.reg[5].out0 = CudaArray(w.ptr+H*X*sizeof(T), (Y,H), w.dev)
    return f
end

mnist_w = fill!(CudaArray(Float32,794*32),0)
mnist_f = initmodel(mnist_w)

function mnist_loss1(w)
    global mnist_f, mnist_w
    w === mnist_w || error("w")
    (x,y) = mnist_trn[1]
    softloss(forw(mnist_f, x), y)
end

function test02(f, data, loss)
    sumloss = numloss = 0
    for (x,ygold) in data
        ypred = forw(f, x)
        sumloss += loss(ypred, ygold)
        numloss += 1
    end
    sumloss / numloss
end

# dhc2(mnist_loss1, mnist_w)
# test02(mnist_f, mnist_tst, softloss)

ext: jl


text:
function create_integrator_cache(int::AbstractIntegratorSPARK{DT,TT}) where {DT,TT}
    IntegratorCacheSPARK{DT, TT, ndims(int), nstages(int), pstages(int)}()
end

function initialize!(int::AbstractIntegratorSPARK, cache::IntegratorCacheSPARK)
    cache.t̅ = cache.t - timestep(int)

    equation(int).v(cache.t, cache.q, cache.p, cache.v)
    equation(int).f(cache.t, cache.q, cache.p, cache.f)

    initialize!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f,
                            cache.t̅, cache.q̅, cache.p̅, cache.v̅, cache.f̅)
end

function update_solution!(int::AbstractIntegratorSPARK{DT,TT}, sol::AtomicSolutionPDAE{DT,TT}) where {DT,TT}
    # compute final update
    update_solution!(sol.q, sol.q̃, int.cache.Vi, int.params.tab.q.b, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.Fi, int.params.tab.p.b, timestep(int))

    # compute projection
    update_solution!(sol.q, sol.q̃, int.cache.Up, int.params.tab.q.β, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.Gp, int.params.tab.p.β, timestep(int))
    # TODO # update_multiplier!(sol.λ, int.cache.Λp, int.params.tab.λ.b)
end

"Integrate an implicit DAE with a specialised partitioned additive Runge-Kutta integrator."
function integrate_step!(int::AbstractIntegratorSPARK{DT,TT}, sol::AtomicSolutionPDAE{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, sol)

    # compute initial guess
    initial_guess!(int, sol)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # check_jacobian(int.solver)
    # print_jacobian(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, int.cache, int.params)

    # compute final update
    update_solution!(int, sol)

    # copy solution to initial guess
    update!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end


text:
# ------------------------------------------------------------------------------
#
#  Gmsh Python tutorial 10
#
#  Mesh size fields
#
# ------------------------------------------------------------------------------

# In addition to specifying target mesh sizes at the points of the geometry (see
# `t1.jl') you can use general mesh
# size "Fields".
import gmsh

gmsh.initialize(append!(["gmsh"], ARGS))

gmsh.model.add("t10")

# Let's create a simple rectangular geometry:
lc = .15
gmsh.model.geo.addPoint(0.0, 0.0, 0, lc, 1)
gmsh.model.geo.addPoint(1, 0.0, 0, lc, 2)
gmsh.model.geo.addPoint(1, 1, 0, lc, 3)
gmsh.model.geo.addPoint(0, 1, 0, lc, 4)
gmsh.model.geo.addPoint(0.2, .5, 0, lc, 5)

gmsh.model.geo.addLine(1, 2, 1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4, 3)
gmsh.model.geo.addLine(4, 1, 4)

gmsh.model.geo.addCurveLoop([1, 2, 3, 4], 5)
gmsh.model.geo.addPlaneSurface([5], 6)

gmsh.model.geo.synchronize()

# Say we would like to obtain mesh elements with size lc/30 near curve 2 and
# point 5, and size lc elsewhere. To achieve this, we can use two fields:
# "Distance", and "Threshold". We first define a Distance field (`Field[1]') on
# points 5 and on curve 2. This field returns the distance to point 5 and to
# (100 equidistant points on) curve 2.
gmsh.model.mesh.field.add("Distance", 1)
gmsh.model.mesh.field.setNumbers(1, "NodesList", [5])
gmsh.model.mesh.field.setNumber(1, "NNodesByEdge", 100)
gmsh.model.mesh.field.setNumbers(1, "EdgesList", [2])

# We then define a `Threshold' field, which uses the return value of the
# `Distance' field 1 in order to define a simple change in element size
# depending on the computed distances
#
# SizeMax -                     /------------------
#                              /
#                             /
#                            /
# SizeMin -o----------------/
#          |                |    |
#        Point         DistMin  DistMax
gmsh.model.mesh.field.add("Threshold", 2)
gmsh.model.mesh.field.setNumber(2, "IField", 1)
gmsh.model.mesh.field.setNumber(2, "LcMin", lc / 30)
gmsh.model.mesh.field.setNumber(2, "LcMax", lc)
gmsh.model.mesh.field.setNumber(2, "DistMin", 0.15)
gmsh.model.mesh.field.setNumber(2, "DistMax", 0.5)

# Say we want to modulate the mesh element sizes using a mathematical function
# of the spatial coordinates. We can do this with the MathEval field:
gmsh.model.mesh.field.add("MathEval", 3)
gmsh.model.mesh.field.setString(3, "F",
                                "Cos(4*3.14*x) * Sin(4*3.14*y) / 10 + 0.101")

# We could also combine MathEval with values coming from other fields. For
# example, let's define a `Distance' field around point 1
gmsh.model.mesh.field.add("Distance", 4)
gmsh.model.mesh.field.setNumbers(4, "NodesList", [1])

# We can then create a `MathEval' field with a function that depends on the
# return value of the `Distance' field 4, i.e., depending on the distance to
# point 1 (here using a cubic law, with minimum element size = lc / 100)
gmsh.model.mesh.field.add("MathEval", 5)
gmsh.model.mesh.field.setString(5, "F", "F4^3 + " * string(lc / 100))

# We could also use a `Box' field to impose a step change in element sizes
# inside a box
gmsh.model.mesh.field.add("Box", 6)
gmsh.model.mesh.field.setNumber(6, "VIn", lc / 15)
gmsh.model.mesh.field.setNumber(6, "VOut", lc)
gmsh.model.mesh.field.setNumber(6, "XMin", 0.3)
gmsh.model.mesh.field.setNumber(6, "XMax", 0.6)
gmsh.model.mesh.field.setNumber(6, "YMin", 0.3)
gmsh.model.mesh.field.setNumber(6, "YMax", 0.6)

# Many other types of fields are available: see the reference manual for a
# complete list. You can also create fields directly in the graphical user
# interface by selecting `Define->Size fields' in the `Mesh' module.

# Finally, let's use the minimum of all the fields as the background mesh field:
gmsh.model.mesh.field.add("Min", 7)
gmsh.model.mesh.field.setNumbers(7, "FieldsList", [2, 3, 5, 6])

gmsh.model.mesh.field.setAsBackgroundMesh(7)

# The API also allows to set a global mesh size callback, which is called each
# time the mesh size is queried
function meshSizeCallback(dim, tag, x, y, z)
    return 0.02 * x + 0.01
end
gmsh.model.mesh.setSizeCallback(meshSizeCallback)

# To determine the size of mesh elements, Gmsh locally computes the minimum of
#
# 1) the size of the model bounding box;
# 2) if `Mesh.MeshSizeFromPoints' is set, the mesh size specified at geometrical
#    points;
# 3) if `Mesh.MeshSizeFromCurvature' is set, the mesh size based on the
#    curvature and `Mesh.MinimumElementsPerTwoPi';
# 4) the background mesh field;
# 5) any per-entity mesh size constraint;
# 6) the mesh size returned by the mesh size callback, if any.
#
# This value is then constrained in the interval [`Mesh.MeshSizeMin',
# `Mesh.MeshSizeMax'] and multiplied by `Mesh.MeshSizeFactor'.  In addition,
# boundary mesh sizes (on curves or surfaces) are interpolated inside the
# enclosed entity (surface or volume, respectively) if the option
# `Mesh.MeshSizeExtendFromBoundary' is set (which is the case by default).
#
# When the element size is fully specified by a background mesh (as it is in
# this example), it is thus often desirable to set
gmsh.option.setNumber("Mesh.MeshSizeExtendFromBoundary", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromPoints", 0)
gmsh.option.setNumber("Mesh.MeshSizeFromCurvature", 0)

# This will prevent over-refinement due to small mesh sizes on the boundary.

gmsh.model.mesh.generate(2)
gmsh.write("t10.msh")

gmsh.finalize()


text:
Base.@kwdef mutable struct MiniProgressBar
    max::Int = 1.0
    header::String = ""
    color::Symbol = :white
    width::Int = 40
    current::Int = 0.0
    prev::Int = 0.0
    has_shown::Bool = false
    time_shown::Float64 = 0.0
    percentage::Bool = true
    always_reprint::Bool = false
    indent::Int = 4
end

const NONINTERACTIVE_TIME_GRANULARITY = Ref(2.0)
const PROGRESS_BAR_PERCENTAGE_GRANULARITY = Ref(0.1)

function showprogress(io::IO, p::MiniProgressBar)
    if p.max == 0
        perc = 0.0
        prev_perc = 0.0
    else
        perc = p.current / p.max * 100
        prev_perc = p.prev / p.max * 100
    end
    # Bail early if we are not updating the progress bar,
    # Saves printing to the terminal
    if !p.always_reprint && p.has_shown && !((perc - prev_perc) > PROGRESS_BAR_PERCENTAGE_GRANULARITY[])
        return
    end
    if !isinteractive()
        t = time()
        if p.has_shown && (t - p.time_shown) < NONINTERACTIVE_TIME_GRANULARITY[]
            return
        end
        p.time_shown = t
    end
    p.prev = p.current
    p.has_shown = true
    n_filled = ceil(Int, p.width * perc / 100)
    n_left = p.width - n_filled
    print(io, " "^p.indent)
    printstyled(io, p.header, color=p.color, bold=true)
    print(io, " [")
    print(io, "="^n_filled, ">")
    print(io, " "^n_left, "]  ", )
    if p.percentage
        @printf io "%2.1f %%" perc
    else
        print(io, p.current, "/",  p.max)
    end
    print(io, "\r")
end


text:
#!/usr/bin/env julia
using BrainFlow

include("julia_tests.jl")
include("brainflow_get_data.jl")
include("brainflow_get_data_twice.jl")
include("signal_filtering.jl")
include("denoising.jl")
include("downsampling.jl")
include("csp.jl")
include("serialization.jl")
include("transforms.jl")
include("eeg_metrics.jl")
include("markers.jl")


text:
"""
## Extract epoch data from SSR

#### Arguments
* `a`: A SSR object
* `valid_triggers`: Trigger numbers that are considered valid ([1,2])
* `remove_first`: Remove the first n triggers (0)
* `remove_last`: Remove the last n triggers (0)

#### Example

```julia
epochs = extract_epochs(SSR, valid_triggers=[1,2])
```
"""
function extract_epochs(a::SSR; valid_triggers::Union{AbstractArray, Int}=[1,2], remove_first::Int=0, remove_last::Int=0, kwargs...)

    merge!(a.processing, Dict("epochs" => extract_epochs(a.data, a.triggers, valid_triggers, remove_first, remove_last)))

    return a
end

function epoch_rejection(a::SSR; retain_percentage::Number=0.95, kwargs...)

    a.processing["epochs"] = epoch_rejection(a.processing["epochs"], retain_percentage; kwargs...)

    return a
end

function create_sweeps(a::SSR; epochsPerSweep::Int=64, kwargs...)

    if epochsPerSweep > size(a.processing["epochs"], 2)
        error("Sweep length is longer than number of epochs will allow")
    end

    merge!(a.processing, Dict("sweeps" => create_sweeps(a.processing["epochs"], epochsPerSweep)))

    return a
end

#######################################
#
# Add triggers for more epochs
#
#######################################

function add_triggers(a::SSR; kwargs...)

    @debug("Adding triggers to reduce SSR. Using SSR modulation frequency")

    add_triggers(a, modulationrate(a); kwargs...)
end

function add_triggers(a::SSR, mod_freq::Number; kwargs...)

    @debug("Adding triggers to reduce SSR. Using $(mod_freq)Hz")

    epochIndex = DataFrame(Code = a.triggers["Code"], Index = [Int(i) for i in a.triggers["Index"]]);
    epochIndex[:Code] = epochIndex[:Code] - 252

    add_triggers(a, mod_freq, epochIndex; kwargs...)
end

function add_triggers(a::SSR, mod_freq::Number, epochIndex; cycle_per_epoch::Int=1, kwargs...)

    @info("Adding triggers to reduce SSR. Reducing $(mod_freq)Hz to $cycle_per_epoch cycle(s).")

    # Existing epochs
    existing_epoch_length   = median(diff(epochIndex[:Index]))     # samples
    existing_epoch_length_s = existing_epoch_length / samplingrate(a)
    @debug("Existing epoch length: $(existing_epoch_length_s)s")

    # New epochs
    new_epoch_length_s = cycle_per_epoch / mod_freq
    new_epochs_num     = round.(existing_epoch_length_s / new_epoch_length_s) - 2
    new_epoch_times    = collect(1:new_epochs_num)*new_epoch_length_s
    new_epoch_indx     = [0; round.(new_epoch_times * samplingrate(a))]
    @debug("New epoch length = $new_epoch_length_s")
    @debug("New # epochs     = $new_epochs_num")

    # Place new epoch indices
    @debug("Was $(length(epochIndex[:Index])) indices")
    new_indx = epochIndex[:Index][1:end-1] .+ new_epoch_indx'
    new_indx = reshape(new_indx', length(new_indx), 1)[1:end-1]
    @debug("Now $(length(new_indx)) indices")

    # Place in dict
    new_code = round.(Int, ones(1, length(new_indx))) .+ 252
    a.triggers = Dict("Index" => vec((new_indx)'), "Code" => vec(new_code), "Duration" => ones(length(new_code), 1)')
    #TODO Possible the trigger duration of one is not long enough

    return a
end

#######################################
#
# Rejection channels
#
#######################################

function channel_rejection(a::SSR; threshold_abs::Number=1000, threshold_var::Number=2, kwargs...)

    # Run on epochs if available, else run on the raw data
    if haskey(a.processing, "epochs")
        data = reshape(a.processing["epochs"],
                size(a.processing["epochs"], 1) * size(a.processing["epochs"], 2), size(a.processing["epochs"],3))
    else
        data = a.data
    end

    valid = vec(channel_rejection(data, threshold_abs, threshold_var))

    @info("Rejected $(sum(.!valid)) channels $(join(channelnames(a)[find(.!valid)], " "))")

    remove_channel!(a, channelnames(a)[.!valid])

    return a
end


text:
"""
    function simu(; nid = 100, nlc = 1000, nql = 100)
---
## Simulate:
- 1000 loci
- 100  ID
- 100  QTL
- Breeding values
- VR-I G
"""
function simu(; nid = 100, nlc = 1000, nql = 100)
    Z   = rand([0., 1., 2.], nid, nlc) # sample genotypes
    qlc = sort(rand(1:nlc, nql))       # sample QTL
    eql = randn(nql) ./ sqrt(nql)      # QTL effects
    ID  = String[]                     # random ID names
    for _ in 1:nid
        push!(ID, randstring(6))
    end
    bv = Z[:, qlc] * eql      # breeding values

    twop = mean(Z, dims=1)    # for G matrix
    Z  .-= twop
    s2pq = (1 .- .5twop) * twop'
    r2pq = 1. / s2pq[1]         # mul faster than div
    G    = Z * Z' .* r2pq + 0.0001I

    ID, bv, G
end

"""
    function homo_R(bv; h2 = 0.8)
---
Simulate phenotypes with homogeneous residuals.
"""
function homo_R(bv; h2 = 0.8)
    nid = length(bv)
    vg  = var(bv)
    ve  = vg/h2*(1-h2)
    se  = sqrt(ve)
    pp  = bv + randn(nid) .* se
end

"""
    function hete_R(bv; h2 = 0.8, mo = 20)
---
Suppose ID may have 1:mo records
"""
function hete_R(bv, h2 = 0.8, mo = 20)
    nid = length(bv)
    nob = rand(1:mo, nid)
    vg  = var(bv)
    ve  = vg/h2*(1-h2)
    se  = sqrt(ve)
    wt  = 1 ./ sqrt.(nob)
    pp  = bv + randn(nid) .* wt
    nob, pp
end


text:
function TOML_print_conversion(x)
    x isa VersionNumber && return "$x"
    error("TOML unhandled type $(typeof(x)).")
end

function tag(
        package_name::AbstractString,
        version::Union{Symbol,VersionNumber,Nothing}=nothing;
        kwargs...
        )

    ctx = Pkg.Types.Context()
    haskey(ctx.env.project.deps, package_name) || error("Unknown package $package_name.")
    pkg_uuid = ctx.env.project.deps[package_name]
    pkg_path = ctx.env.manifest[pkg_uuid].path
    pkg_path===nothing && error("Package must be deved to be tagged.")

    tag_internal(package_name, pkg_uuid, pkg_path, version; kwargs...)
end

function tag(
        package_path::AbstractPath,
        version::Union{Symbol,VersionNumber,Nothing}=nothing;
        kwargs...)

    project_toml = isfile(joinpath(package_path, "Project.toml")) ? joinpath(package_path, "Project.toml") : isfile(joinpath(package_path, "JuliaProject.toml")) ? joinpath(package_path, "JuliaProject.toml") : nothing
    project_toml === nothing && error("Could not find a 'Project.toml' at $package_path.")

    project_content = TOML.parsefile(string(project_toml))

    package_name = get(project_content, "name", nothing)
    package_name===nothing && error("The project toml for the package doesn't contain a name.")

    pkg_uuid = get(project_content, "uuid", nothing)
    pkg_uuid===nothing && error("The project toml for the package doesn't contain a uuid.")

    tag_internal(package_name, pkg_uuid, string(package_path), version; kwargs...)
end

# This is adapted over from LocalRegistry.jl
function our_collect_registries()
    registries = []
    for depot in Pkg.depots()
        isdir(depot) || continue
        reg_dir = joinpath(depot, "registries")
        isdir(reg_dir) || continue
        for name in readdir(reg_dir)
            file = joinpath(reg_dir, name, "Registry.toml")
            if !isfile(file)
                # Packed registry in Julia 1.7+.
                file = joinpath(reg_dir, "$(name).toml")
            end

            if isfile(file)
                content =TOML.parsefile(file)

                spec = (
                    name = content["name"]::String,
                    uuid = UUID(content["uuid"]::String),
                    url = get(content, "repo", nothing)::Union{String,Nothing},
                    path = file
                )

                push!(registries, spec)
            end
        end
    end
    return registries
end

function tag_internal(
        package_name::AbstractString,
        pkg_uuid, pkg_path::AbstractString,
        version::Union{Symbol,VersionNumber,Nothing}=nothing;
        registry::Union{AbstractString,Nothing}=nothing,
        release_notes::Union{AbstractString,Nothing}=nothing,
        credentials::Union{AbstractString, Nothing}=nothing,
        github_username::Union{AbstractString, Nothing} = nothing)

    general_reg_url = "https://github.com/JuliaRegistries/General"

    if github_username===nothing
        github_username = LibGit2.getconfig("github.user", "")
        github_username == "" && error("You need to configure the github.user setting.")
    end

    isdir(pkg_path) || error("Path for package does not exist on disc.")

    all_registries = our_collect_registries()

    private_reg_url, private_reg_uuid = if registry===nothing
        registries_that_contain_the_package = []
        for reg_spec in all_registries
            reg_data = TOML.parsefile(reg_spec.path)

            if haskey(reg_data["packages"], string(pkg_uuid))
                push!(registries_that_contain_the_package, reg_spec)
            end
        end

        if length(registries_that_contain_the_package)==0
            (nothing, nothing)
        elseif length(registries_that_contain_the_package)>1
            error("Package is registered in more than one registry, please specify in which you want to register the tag.")
        else
            (registries_that_contain_the_package[1].url, registries_that_contain_the_package[1].uuid)
        end
    else
        relevant_registry = findfirst(i->i.name==registry, all_registries)

        relevant_registry===nothing && error("The registry $registry does not exist.")

        (all_registries[relevant_registry].url, all_registries[relevant_registry].uuid)
    end

    if private_reg_uuid == UUID("23338594-aafe-5451-b93e-139f81909106")
        private_reg_url = nothing
        private_reg_uuid = nothing
    end

    if credentials===nothing
        creds = LibGit2.GitCredential(GitConfig(), "https://github.com")

        creds.password===nothing && error("Did not find credentials for github.com in the git credential manager.")

        credentials = read(creds.password, String)
        Base.shred!(creds.password)
    end

    myauth = GitHub.authenticate(credentials)

    registry_github_owner_repo_name = private_reg_url===nothing ? "JuliaRegistries/General" : get_repo_onwer_from_url(private_reg_url)

    gh_registry_repo = GitHub.repo(registry_github_owner_repo_name, auth=myauth)
    gh_forks = GitHub.forks(gh_registry_repo, auth=myauth)
    fork_index = findfirst(i->i.owner.login==github_username, gh_forks[1])
    fork_index===nothing && error("You need to have a fork of the registry in your github account.")
    registry_fork_url = string(gh_forks[1][fork_index].html_url)

    pkg_repo = GitRepo(pkg_path)

    try

        LibGit2.isdirty(pkg_repo) && error("The repo for the package cannot be dirty.")

        pkg_remote = LibGit2.lookup_remote(pkg_repo, "origin")

        pkg_remote===nothing && error("The package must have a remote called origin.")

        pkg_url = LibGit2.url(pkg_remote)

        pkg_project_toml_path = isfile(joinpath(pkg_path, "JuliaProject.toml")) ? joinpath(pkg_path, "JuliaProject.toml") : isfile(joinpath(pkg_path, "Project.toml")) ? joinpath(pkg_path, "Project.toml") : error("Couldn't find Project.toml.")

        pkg_toml_content = TOML.parsefile(pkg_project_toml_path)

        haskey(pkg_toml_content, "version") || error("Project.toml must have a version field.")

        current_version_in_pkg = VersionNumber(pkg_toml_content["version"])

        version_to_be_tagged = if version===nothing
            current_version_in_pkg.prerelease!=("DEV",) && current_version_in_pkg.build!=() && error("Version in Project.toml must have format of x.y.z-DEV.")
            VersionNumber(current_version_in_pkg.major, current_version_in_pkg.minor, current_version_in_pkg.patch)
        elseif version isa VersionNumber
            version
        elseif version==:major
            VersionNumber(current_version_in_pkg.major+1, 0, 0)
        elseif version==:minor
            VersionNumber(current_version_in_pkg.major, current_version_in_pkg.minor+1, 0)
        elseif version==:patch && current_version_in_pkg.prerelease==("DEV",)
            VersionNumber(current_version_in_pkg.major, current_version_in_pkg.minor, current_version_in_pkg.patch)
        elseif version==:patch
            VersionNumber(current_version_in_pkg.major, current_version_in_pkg.minor, current_version_in_pkg.patch+1)
        else
            error("Invalid argument for version, must be nothing, a VersionNumber, :major, :minor or :patch.")
        end

        next_version = VersionNumber(version_to_be_tagged.major, version_to_be_tagged.minor, version_to_be_tagged.patch+1, ("DEV",))

        # TODO Check whether that version already exists, and if so, error.

        name_of_release_branch = "release-$version_to_be_tagged"

        LibGit2.lookup_branch(pkg_repo, name_of_release_branch)===nothing || error("A branch named $name_of_release_branch already exists in the package.")

        name_of_old_branch_in_pkg = LibGit2.headname(pkg_repo)

        LibGit2.branch!(pkg_repo, name_of_release_branch, force=true)

        # Now update the version field in the Project.toml

        pkg_toml_content["version"] = version_to_be_tagged

        open(pkg_project_toml_path, "w") do f
            TOML.print(TOML_print_conversion, f, pkg_toml_content)
        end

        project_as_it_should_be_tagged = Pkg.Types.read_project(pkg_project_toml_path)

        LibGit2.add!(pkg_repo, splitdir(pkg_project_toml_path)[2])
        hash_of_commit_to_be_tagged = LibGit2.commit(pkg_repo, "Set version to v$version_to_be_tagged")

        tree_hash_of_commit_to_be_tagged = LibGit2.GitHash(LibGit2.peel(LibGit2.GitTree, LibGit2.GitCommit(pkg_repo, hash_of_commit_to_be_tagged)))

        # Now update the version field in the Project.toml

        pkg_toml_content = TOML.parsefile(pkg_project_toml_path)

        pkg_toml_content["version"] = next_version

        open(pkg_project_toml_path, "w") do f
            TOML.print(TOML_print_conversion, f, pkg_toml_content)
        end

        LibGit2.add!(pkg_repo, splitdir(pkg_project_toml_path)[2])
        LibGit2.commit(pkg_repo, "Set version to v$next_version")

        run(Cmd(`git push origin refs/heads/$name_of_release_branch`, dir=pkg_path))

        LibGit2.branch!(pkg_repo, name_of_old_branch_in_pkg)

        LibGit2.delete_branch(LibGit2.lookup_branch(pkg_repo, name_of_release_branch))

        pkg_owner_repo_name = get_repo_onwer_from_url(pkg_url)

        gh_pkg_repo = GitHub.repo(pkg_owner_repo_name, auth=myauth)

        GitHub.create_pull_request(gh_pkg_repo, auth=myauth, params=Dict(:title=>"New version: v$version_to_be_tagged", :head=>name_of_release_branch, :base=>name_of_old_branch_in_pkg, :body=>""))

        if private_reg_url===nothing
            GitHub.create_comment(gh_pkg_repo, string(hash_of_commit_to_be_tagged), :commit, params = Dict("body"=>"@JuliaRegistrator register()"), auth=myauth)
        else
            mktempdir() do tmp_path
                cd(tmp_path) do
                    folder_for_registry = nothing
                    regbranch = if private_reg_url===nothing
                        folder_for_registry = joinpath(tmp_path, "registries", "23338594-aafe-5451-b93e-139f81909106")
                        RegistryTools.register(pkg_url, project_as_it_should_be_tagged, string(tree_hash_of_commit_to_be_tagged); registry=general_reg_url, push=false)
                    else
                        folder_for_registry = joinpath(tmp_path, "registries", string(private_reg_uuid))
                        RegistryTools.register(pkg_url, project_as_it_should_be_tagged, string(tree_hash_of_commit_to_be_tagged); registry=private_reg_url, registry_deps=[general_reg_url], push=false)
                    end

                    @info regbranch.metadata

                    registry_repo = GitRepo(folder_for_registry)
                    try
                        run(Cmd(`git push $registry_fork_url refs/heads/$(regbranch.branch)`, dir=folder_for_registry))
                    finally
                        close(registry_repo)
                    end

                    body = ""
                    if release_notes !== nothing
                        # Prepend every line with '> ' to quote it (this format is expected by TagBot).
                        notes = join(map(line -> "> $line", split(release_notes, "\n")), "\n")
                        body *= """

                            Release notes:
                            <!-- BEGIN RELEASE NOTES -->
                            $notes
                            <!-- END RELEASE NOTES -->
                            """
                    end

                    GitHub.create_pull_request(gh_registry_repo, auth=myauth, params=Dict(:title=>"New version: $package_name v$version_to_be_tagged", :head=>"$github_username:$(regbranch.branch)", :base=>"master", :body=>strip(body)))
                end
            end
        end
    finally
        close(pkg_repo)
    end

    return nothing
end


text:
const FORCE_COLOR = Ref(false)
const FORCE_256_COLORS = Ref(false)

force_color(b::Bool) = FORCE_COLOR[] = b
force_256_colors(b::Bool) = FORCE_256_COLORS[] = b

_force_color() = FORCE_COLOR[] || haskey(ENV, "FORCE_COLOR")
_force_256_colors() = FORCE_256_COLORS[] || haskey(ENV, "FORCE_256_COLORS")

const CSI = "\e["
const ESCAPED_CSI = "\\e["
const END_ANSI = "m"

# Add 30 to get fg ANSI
# Add 40 to get bg ANSI
const COLORS = Dict(
    :black         => 0,
    :red           => 1,
    :green         => 2,
    :yellow        => 3,
    :blue          => 4,
    :magenta       => 5,
    :cyan          => 6,
    :light_gray    => 7,
    :default       => 9,
    :dark_gray     => 60,
    :light_red     => 61,
    :light_green   => 62,
    :light_yellow  => 63,
    :light_blue    => 64,
    :light_magenta => 65,
    :light_cyan    => 66,
    :white         => 67
)

@enum(ColorMode,
RESET,
COLORS_16,
COLORS_256,
COLORS_24BIT)

struct ANSIColor
    r::Int # [0-9, 60-69], for 16 colors, 0-255 for 256 colors and 24 bit
    g::Int # [0-255] Only used for 24 bit colors
    b::Int # [0-255] Only used for 24 bit colors
    style::ColorMode
    active::Bool
    function ANSIColor(r::Int, g::Int = 0, b::Int = 0, style::ColorMode = COLORS_16, active = true)
        for v in (r, g, b)
            !(0 <= v <= 255) && throw(ArgumentError("RGB color component has to be between 0 and 255"))
        end
        return new(r, g, b, style, active)
    end
end

ANSIColor() = ANSIColor(0, 0, 0, COLORS_16, false)
ANSIColor(val::Int, style::ColorMode, active::Bool = true) = ANSIColor(val, 0, 0, style, active)

red(x::ANSIColor) = x.r
green(x::ANSIColor) = x.g
blue(x::ANSIColor) = x.b
val(x::ANSIColor) = x.r

# The inverse sets the color to default.
# No point making active if color already is default
Base.inv(x::ANSIColor) = ANSIColor(9, 0, 0, COLORS_16, x.active && !(x.style == COLORS_16 && x.r == 9))

struct ANSIStyle
    on::Bool
    active::Bool
end

ANSIStyle() = ANSIStyle(false, false)
ANSIStyle(v::Bool) = ANSIStyle(v, true)

# The inverse always sets the thing to false
# No point in setting active if the style is off.
Base.inv(x::ANSIStyle) = ANSIStyle(false, x.active && x.on)

struct Crayon
    fg::ANSIColor
    bg::ANSIColor

    reset::ANSIStyle
    bold::ANSIStyle
    faint::ANSIStyle
    italics::ANSIStyle
    underline::ANSIStyle
    blink::ANSIStyle
    negative::ANSIStyle
    conceal::ANSIStyle
    strikethrough::ANSIStyle
end

anyactive(x::Crayon) = ((x.reset.active && x.reset.on) ||
                        x.fg.active    || x.bg.active       || x.bold.active      ||
                        x.faint.active || x.italics.active  || x.underline.active ||
                        x.blink.active || x.negative.active || x.conceal.active   || x.strikethrough.active)

Base.inv(c::Crayon) = Crayon(inv(c.fg), inv(c.bg), ANSIStyle(), # no point taking inverse of reset,
                             inv(c.bold), inv(c.faint), inv(c.italics), inv(c.underline),
                             inv(c.blink), inv(c.negative), inv(c.conceal), inv(c.strikethrough))

function Base.print(io::IO, x::Crayon)
    io_has_color = get(io, :color, false)
    if anyactive(x) && (io_has_color || _force_color())
        print(io, CSI)
        if (x.fg.style == COLORS_24BIT || x.bg.style == COLORS_24BIT) && _force_256_colors()
            x = _to256(x)
        end
        _print(io, x)
        print(io, END_ANSI)
    end
end

function Base.show(io::IO, x::Crayon)
    io_has_color = get(io, :color, false)
    if anyactive(x) && io_has_color
        print(io, x)
        print(io, ESCAPED_CSI)
        _print(io, x)
        print(io, END_ANSI, CSI, "0", END_ANSI)
    end
end

_ishex(c::Char) = isdigit(c) || ('a' <= c <= 'f') || ('A' <= c <= 'F')
_torgb(hex::UInt32) = Int(hex << 8 >> 24), Int(hex << 16 >> 24), Int(hex << 24 >> 24)

function _parse_color(c::Union{Integer,Symbol,NTuple{3,Integer},UInt32})
    ansicol = ANSIColor()
    if c != :nothing
        if isa(c, Symbol)
            ansicol = ANSIColor(COLORS[c], COLORS_16)
        elseif isa(c, UInt32)
            r, g, b = _torgb(c)
            ansicol = ANSIColor(r, g, b, COLORS_24BIT)
        elseif isa(c, Integer)
            ansicol = ANSIColor(c, COLORS_256)
        elseif isa(c, NTuple{3,Integer})
            ansicol = ANSIColor(c[1], c[2], c[3], COLORS_24BIT)
        else
            error("should not happen")
        end
    end
    return ansicol
end

function Crayon(;foreground::Union{Int,Symbol,NTuple{3,Integer},UInt32} = :nothing,
                 background::Union{Int,Symbol,NTuple{3,Integer},UInt32} = :nothing,
                 reset = :nothing,
                 bold = :nothing,
                 faint = :nothing,
                 italics = :nothing,
                 underline = :nothing,
                 blink = :nothing,
                 negative = :nothing,
                 conceal = :nothing,
                 strikethrough = :nothing)

    fgcol = _parse_color(foreground)
    bgcol = _parse_color(background)

    _reset         = ANSIStyle()
    _bold          = ANSIStyle()
    _faint         = ANSIStyle()
    _italics       = ANSIStyle()
    _underline     = ANSIStyle()
    _blink         = ANSIStyle()
    _negative      = ANSIStyle()
    _conceal       = ANSIStyle()
    _strikethrough = ANSIStyle()

    reset         != :nothing && (_reset         = ANSIStyle(reset))
    bold          != :nothing && (_bold          = ANSIStyle(bold))
    faint         != :nothing && (_faint         = ANSIStyle(faint))
    italics       != :nothing && (_italics       = ANSIStyle(italics))
    underline     != :nothing && (_underline     = ANSIStyle(underline))
    blink         != :nothing && (_blink         = ANSIStyle(blink))
    negative      != :nothing && (_negative      = ANSIStyle(negative))
    conceal       != :nothing && (_conceal       = ANSIStyle(conceal))
    strikethrough != :nothing && (_strikethrough = ANSIStyle(strikethrough))

    return Crayon(fgcol,
                  bgcol,
                  _reset,
                  _bold,
                  _faint,
                  _italics,
                  _underline,
                  _blink,
                  _negative,
                  _conceal,
                  _strikethrough)
end

# Prints the crayon without the inital and terminating ansi escape sequences
function _print(io::IO, c::Crayon)
    first_active = true
    if c.reset.active && c.reset.on
        first_active = false
        print(io, "0")
    end

    for (col, num) in ((c.fg, 30),
                       (c.bg, 40))
        if col.active
            !first_active && print(io, ";")
            first_active = false

            col.style == COLORS_16    && print(io, val(col) + num)
            col.style == COLORS_256   && print(io, num + 8, ";5;", val(col))
            col.style == COLORS_24BIT && print(io, num + 8, ";2;", red(col), ";", green(col), ";", blue(col))
        end
    end

    for (style, val) in ((c.bold, 1),
                         (c.faint, 2),
                         (c.italics, 3),
                         (c.underline, 4),
                         (c.blink, 5),
                         (c.negative, 7),
                         (c.conceal, 8),
                         (c.strikethrough, 9))

        if style.active
            !first_active && print(io, ";")
            first_active = false

            style.on && print(io, val)
            # Bold off is actually 22 so special case for val == 1
            !style.on && print(io, val == 1 ? val + 21 : val + 20)
        end
    end
    return nothing
end

function Base.merge(a::Crayon, b::Crayon)
    fg            = b.fg.active            ? b.fg            : a.fg
    bg            = b.bg.active            ? b.bg            : a.bg
    reset         = b.reset.active         ? b.reset         : a.reset
    bold          = b.bold.active          ? b.bold          : a.bold
    faint         = b.faint.active         ? b.faint         : a.faint
    italics       = b.italics.active       ? b.italics       : a.italics
    underline     = b.underline.active     ? b.underline     : a.underline
    blink         = b.blink.active         ? b.blink         : a.blink
    negative      = b.negative.active      ? b.negative      : a.negative
    conceal       = b.conceal.active       ? b.conceal       : a.conceal
    strikethrough = b.strikethrough.active ? b.strikethrough : a.strikethrough

    return Crayon(fg,
                  bg,
                  reset,
                  bold,
                  faint,
                  italics,
                  underline,
                  blink,
                  negative,
                  conceal,
                  strikethrough)
end

Base.:*(a::Crayon, b::Crayon) = merge(a, b)

function Base.merge(toks::Crayon...)
    if length(toks) == 0
        return Crayon()
    end
    tok = toks[1]
    for i in 2:length(toks)
        tok = merge(tok, toks[i])
    end
    return tok
end

# 24bit -> 256 colors
function _to256(crayon::Crayon)
    fg = crayon.fg
    bg = crayon.bg
    crayon.fg.style == COLORS_24BIT && (fg = _to256(crayon.fg))
    crayon.bg.style == COLORS_24BIT && (bg = _to256(crayon.bg))
    return Crayon(
        fg,
        bg,
        crayon.reset,
        crayon.bold,
        crayon.faint,
        crayon.italics,
        crayon.underline,
        crayon.blink,
        crayon.negative,
        crayon.conceal,
        crayon.strikethrough,
    )
end

function _to256(color::ANSIColor)
    @assert color.style == COLORS_24BIT
    r, g, b = color.r, color.g, color.b
    r24, g24, b24 = map(c->round(Int, c * 23 / 256), (r, g, b))
    if r24 == g24 == b24
        return ANSIColor(232 + r24, COLORS_256, color.active)
    else
        r6, g6, b6 = map(c->round(Int, c * 5  / 256), (r, g, b))
        return ANSIColor(16 + 36 * r6 + 6 * g6 + b6, COLORS_256, color.active)
    end
end


text:
using Test
using Dispatcher
using Memento
using TranscodingStreams
using CodecBzip2
using CodecZlib
using JSON
using DispatcherCache

# Set Dispatcher logging level to warning
setlevel!(getlogger("Dispatcher"), "warn")

# Run tests
include("compression.jl")
include("hash.jl")
include("core.jl")


text: Pkg.build("Ipopt")


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Testing

@testset "test tet10 stiffness matrix" begin
    el = Element(Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
    el["youngs modulus"] = 480.0
    el["poissons ratio"] = 1/3
    x1 = [2.0, 3.0, 4.0]
    x2 = [6.0, 3.0, 2.0]
    x3 = [2.0, 5.0, 1.0]
    x4 = [4.0, 3.0, 6.0]
    x5 = 0.5*(x1+x2)
    x6 = 0.5*(x2+x3)
    x7 = 0.5*(x3+x1)
    x8 = 0.5*(x1+x4)
    x9 = 0.5*(x2+x4)
    x10 = 0.5*(x3+x4)
    X = Dict{Int64, Vector{Float64}}(
        1 => x1, 2 => x2, 3 => x3, 4 => x4, 5 => x5,
        6 => x6, 7 => x7, 8 => x8, 9 => x9, 10 => x10)
    u = Dict{Int64, Vector{Float64}}()
    for i=1:10
        u[i] = [0.0, 0.0, 0.0]
    end
    update!(el, "geometry", X)
    update!(el, "displacement", u)
    pr = Problem(Elasticity, "tet10", 3)
    ass = Assembly()
    assemble!(ass, pr, el, 0.0)
    Kt = full(ass.K)
    eigs = real(eigvals(Kt))
    eigs_expected = [8809.45, 4936.01, 2880.56, 2491.66, 2004.85,
                     1632.49, 1264.32, 1212.42, 817.905,
                     745.755, 651.034, 517.441, 255.1, 210.955,
                     195.832, 104.008, 72.7562, 64.4376, 53.8515,
                     23.8417, 16.6354, 9.54682, 6.93361, 2.22099,
                     0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
    @test isapprox(eigs, eigs_expected; atol=1.0e-2)
end


text:
"""
    bimanual block
        particle with contacts at each corner
"""
struct BimanipulationBlockV2{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int

    mass_block
	mass_pusher

    J
    μ_surface
	μ_pusher
    g

    contact_corner_offset
	block_dim
	block_rnd

    nq
    nu
    nc
    nf
    nb

    idx_u
    idx_λ
    idx_b
    idx_ψ
    idx_η
    idx_s
end

# signed distance for a box
function sd_box(p, dim)
	q = abs.(p) - dim
	norm(max.(q, 1.0e-32)) + min(maximum(q), 0.0)
end

function sd_2d_box(p, pose, dim, rnd)
	x, y, θ = pose
	R = rotation_matrix(-θ)
	p_rot = R * (p - pose[1:2])

	return sd_box(p_rot, dim .- rnd) - rnd
end

# Dimensions
nq = 7 # configuration dimension
nu = 4 # control dimension
nc = 8 # number of contact points
nf = 2 # number of faces for friction cone pyramid
nb = nc * nf
ns = 1

n = 2 * nq
m = nu + nc + nb + nc + nb + 1
d = 0

idx_u = (1:nu)
idx_λ = nu .+ (1:nc)
idx_b = nu + nc .+ (1:nb)
idx_ψ = nu + nc + nb .+ (1:nc)
idx_η = nu + nc + nb + nc .+ (1:nb)
idx_s = nu + nc + nb + nc + nb .+ (1:ns)

# Kinematics
r = 0.1

# contact corner
cc1 = @SVector [r, r]
cc2 = @SVector [-r, r]
cc3 = @SVector [r, -r]
cc4 = @SVector [-r, -r]

# control location
# cu1 = @SVector [r, r]
# cu2 = @SVector [-r, r]
# cu3 = @SVector [-r, -r]
# cu4 = @SVector [r, -r]
# cu5 = @SVector [r, 0.0]
# cu6 = @SVector [0.0, r]
# cu7 = @SVector [-r, 0.0]
# cu8 = @SVector [0.0, -r]

contact_corner_offset = @SVector [cc1, cc2, cc3, cc4]
# control_input_offset = @SVector [cu1, cu2, cu3, cu4, cu5, cu6, cu7, cu8]
# control_input_offset = @SVector [cu5, cu6, cu7, cu8]

# Parameters
μ_surface = 0.1  # coefficient of friction
μ_pusher = 0.5
g = 9.81
mass_block = 1.0   # mass
mass_pusher = 10.0
J = 1.0 / 12.0 * mass_block * ((2.0 * r)^2 + (2.0 * r)^2)

rnd = 0.01
dim = [r, r]

# Methods
M_func(model::BimanipulationBlockV2, q) = Diagonal(@SVector [model.mass_block, model.mass_block,
	model.J, model.mass_pusher, model.mass_pusher, model.mass_pusher, model.mass_pusher])

function C_func(model::BimanipulationBlockV2, q, q̇)
	SVector{7}([0.0, model.mass_block * model.g, 0.0, 0.0, 0.0, 0.0, 0.0])
end

function rotation_matrix(x)
	SMatrix{2,2}([cos(x) -sin(x); sin(x) cos(x)])
end

function ϕ_func(model::BimanipulationBlockV2, q)
    p_block = view(q, 1:3)
	p_pusher1 = view(q, 4:5)
	p_pusher2 = view(q, 6:7)

	θ = q[3]
    R = rotation_matrix(θ)

	@SVector [(p_block[1:2] + R * model.contact_corner_offset[1])[2],
              (p_block[1:2] + R * model.contact_corner_offset[2])[2],
              (p_block[1:2] + R * model.contact_corner_offset[3])[2],
              (p_block[1:2] + R * model.contact_corner_offset[4])[2],
			  q[5],
			  q[7],
			  sd_2d_box(p_pusher1, p_block, model.block_dim, model.block_rnd),
			  sd_2d_box(p_pusher2, p_block, model.block_dim, model.block_rnd)]
end

function B_func(model::BimanipulationBlockV2, q)
	SMatrix{7,4}([0.0 0.0 0.0 0.0;
				  0.0 0.0 0.0 0.0;
				  0.0 0.0 0.0 0.0;
				  1.0 0.0 0.0 0.0;
				  0.0 1.0 0.0 0.0;
				  0.0 0.0 1.0 0.0;
				  0.0 0.0 0.0 1.0;])
end

function N_func(model::BimanipulationBlockV2, q)
    tmp(z) = ϕ_func(model, z)
    ForwardDiff.jacobian(tmp, q)
end

function P_func(model::BimanipulationBlockV2, q)
	map1 = [1.0;
	        -1.0]

    # map2 = [1. 0.;
    #        0. 1.;
    #        -1. 0.;
    #        0. -1.]

    function p(x)
        pos = view(x, 1:2)
		θ = x[3]
        R = rotation_matrix(θ)

        [map1 * (pos + R * model.contact_corner_offset[1])[1];
         map1 * (pos + R * model.contact_corner_offset[2])[1];
         map1 * (pos + R * model.contact_corner_offset[3])[1];
         map1 * (pos + R * model.contact_corner_offset[4])[1]]
    end

    P_block = ForwardDiff.jacobian(p, q)

	# pusher block
	p_block = view(q, 1:3)
	p_pusher1 = view(q, 4:5)
	p_pusher2 = view(q, 6:7)

	sd_p(x) = sd_2d_box(x, p_block, model.block_dim, model.block_rnd)

	Np1 = ForwardDiff.gradient(sd_p, p_pusher1)

	n_dir1 = Np1[1:2] ./ norm(Np1[1:2])
	t_dir1 = rotation_matrix(0.5 * π) * n_dir1

	r1 = p_pusher1 - p_block[1:2]
	m1 = cross([r1; 0.0], [t_dir1; 0.0])[3]

	P_pusherblock1 = map1 * [t_dir1[1]; t_dir1[2]; m1; -t_dir1[1]; -t_dir1[2]; 0.0; 0.0]'

	Np2 = ForwardDiff.gradient(sd_p, p_pusher2)

	n_dir2 = Np2[1:2] ./ norm(Np2[1:2])
	t_dir2 = rotation_matrix(0.5 * π) * n_dir2

	r2 = p_pusher2 - p_block[1:2]
	m2 = cross([r2; 0.0], [t_dir2; 0.0])[3]

	P_pusherblock2 = map1 * [t_dir2[1]; t_dir2[2]; m2; 0.0; 0.0; -t_dir2[1]; -t_dir2[2]]'

	return [P_block; map1 * [0.0; 0.0; 0.0; 1.0; 0.0; 0.0; 0.0]'; map1 * [0.0; 0.0; 0.0; 0.0; 0.0; 1.0; 0.0]'; P_pusherblock1; P_pusherblock2]
end

function friction_cone(model::BimanipulationBlockV2,u)
    λ = u[model.idx_λ]
    b = u[model.idx_b]

    @SVector [model.μ_surface[1] * λ[1] - sum(b[1:2]),
              model.μ_surface[2] * λ[2] - sum(b[3:4]),
              model.μ_surface[3] * λ[3] - sum(b[5:6]),
              model.μ_surface[4] * λ[4] - sum(b[7:8]),
			  model.μ_surface[5] * λ[5] - sum(b[9:10]),
			  model.μ_surface[6] * λ[6] - sum(b[11:12]),
			  model.μ_pusher * λ[7] - sum(b[13:14]),
			  model.μ_pusher * λ[8] - sum(b[15:16])]
end

function maximum_dissipation(model::BimanipulationBlockV2, x⁺, u, h)
    q3 = x⁺[model.nq .+ (1:model.nq)]
	q2 = x⁺[1:model.nq]

    ψ = u[model.idx_ψ]
    ψ_stack = [ψ[1] * ones(2);
               ψ[2] * ones(2);
               ψ[3] * ones(2);
               ψ[4] * ones(2);
			   ψ[5] * ones(2);
			   ψ[6] * ones(2);
			   ψ[7] * ones(2);
			   ψ[8] * ones(2)]

    η = u[model.idx_η]

    P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function lagrangian_derivatives(model, q, v)
	D1L = -1.0 * C_func(model, q, v)
    D2L = M_func(model, q) * v
	return D1L, D2L
end

function fd(model::BimanipulationBlockV2{Discrete, FixedTime}, x⁺, x, u, w, h, t)
	q3 = view(x⁺, model.nq .+ (1:model.nq))
	q2⁺ = view(x⁺, 1:model.nq)
	q2⁻ = view(x, model.nq .+ (1:model.nq))
	q1 = view(x, 1:model.nq)
	u_ctrl = view(u, model.idx_u)
	λ = view(u, model.idx_λ)
	b = view(u, model.idx_b)

	qm1 = 0.5 * (q1 + q2⁺)
    vm1 = (q2⁺ - q1) / h
    qm2 = 0.5 * (q2⁺ + q3)
    vm2 = (q3 - q2⁺) / h

	D1L1, D2L1 = lagrangian_derivatives(model, qm1, vm1)
	D1L2, D2L2 = lagrangian_derivatives(model, qm2, vm2)

	[q2⁺ - q2⁻;
     (0.5 * h * D1L1 + D2L1 + 0.5 * h * D1L2 - D2L2
     + B_func(model, qm2) * SVector{4}(u_ctrl[1:4])
     + transpose(N_func(model, q3)) * SVector{8}(λ)
     + transpose(P_func(model, q3)) * SVector{16}(b))]
end

model = BimanipulationBlockV2{Discrete, FixedTime}(n, m, d,
			mass_block, mass_pusher, J,
			[μ_surface for i = 1:nc], μ_pusher,
			g,
			contact_corner_offset,
			dim,
			rnd,
            nq, nu, nc, nf, nb,
            idx_u,
            idx_λ,
            idx_b,
            idx_ψ,
            idx_η,
            idx_s)

function visualize!(vis, model::BimanipulationBlockV2, q, u; r = r,
        Δt = 0.1)

	default_background!(vis)

    setobject!(vis["box"], GeometryBasics.Rect(Vec(-1.0 * r,
		-1.0 * r,
		-1.0 * r),
		Vec(2.0 * r, 2.0 * r, 2.0 * r)),
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    setobject!(vis["pusher1"], Sphere(Point3f0(0),
        convert(Float32, 0.01)),
        MeshPhongMaterial(color = RGBA(1.0, 1.0, 0, 1.0)))

    setobject!(vis["pusher2"], Sphere(Point3f0(0),
        convert(Float32, 0.01)),
        MeshPhongMaterial(color = RGBA(1.0, 1.0, 0, 1.0)))

    # for i = 1:4
    #     setobject!(vis["contact$i"], GeometryBasics.Sphere(Point3f0(0),
    #         convert(Float32, 0.02)),
    #         MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))
    # end

	# force_vis1 = ArrowVisualizer(vis[:force1])
	# setobject!(force_vis1, MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0)))
    #
	# us = u[1] / 2.0
    #
	# settransform!(force_vis1,
	# 			Point(q[1][4] - us[1], 0, q[1][5] - us[1]),
	# 			Vec(us[1], 0, us[2]),
	# 			shaft_radius=0.01,
	# 			max_head_radius=0.025)
    #
	# force_vis2 = ArrowVisualizer(vis[:force2])
	# setobject!(force_vis2, MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0)))
    #
	# settransform!(force_vis2,
	# 			Point(q[1][6] - us[3], 0, q[1][7] - us[4]),
	# 			Vec(us[3], 0, us[4]),
	# 			shaft_radius=0.01,
	# 			max_head_radius=0.025)

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	T = length(q)
    for t = 1:T-1
        MeshCat.atframe(anim, t) do
			if t < T-1

				# # pusher 1
				# if norm(u[t][1:2]) < 1.0e-6
				# 	setvisible!(vis[:force1], false)
				# else
				# 	setvisible!(vis[:force1], true)
                #
				# 	us = u[t] / 10.0
				# 	settransform!(force_vis1,
				# 				Point(q[t+1][4] - us[1], 0, q[t+1][5] - us[2]),
				# 				Vec(us[1], 0, us[2]),
				# 				shaft_radius=0.01,
				# 				max_head_radius=0.025)
				# end

				# # pusher 2
				# if norm(u[t][3:4]) < 1.0e-6
				# 	setvisible!(vis[:force2], false)
				# else
				# 	setvisible!(vis[:force2], true)
                #
				# 	us = u[t] / 2.0
				# 	settransform!(force_vis2,
				# 				Point(q[t+1][6] - us[3], 0.0, q[t+1][7] - us[4]),
				# 				Vec(us[3], 0.0, us[4]),
				# 				shaft_radius=0.01,
				# 				max_head_radius=0.025)
				# end
			end

            settransform!(vis["box"],
				compose(Translation(q[t+1][1], 0.0, q[t+1][2]), LinearMap(RotY(q[t+1][3]))))

            settransform!(vis["pusher1"],
				compose(Translation(q[t+1][4], 0.0, q[t+1][5]), LinearMap(RotY(0.0))))

            settransform!(vis["pusher2"],
                compose(Translation(q[t+1][6], 0.0, q[t+1][7]), LinearMap(RotY(0.0))))
            #
            # for i = 1:4
            #     settransform!(vis["contact$i"],
            #         Translation(([q[t+1][1]; 0.0; q[t+1][2]] + RotY(q[t+1][3]) * [contact_corner_offset[i][1]; 0.0; contact_corner_offset[i][2];])...))
            # end
        end
    end

	settransform!(vis["/Cameras/default"],
		compose(Translation(0.0, -90.0, -1.0),LinearMap(RotZ(pi / 2.0))))
	setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 75)

    MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
function hcat_by_sum(xs::AbstractVector)
    T = eltype(xs[1])
    dim = size(xs[1], 1)
    N = length(xs)

    ns = map(x->size(x,2), xs)
    pushfirst!(ns, 1)
    cumsum!(ns, ns)

    A = similar(xs[1], T, dim, ns[end]-1)
    for i in 1:N
        A[:, ns[i]:(ns[i+1]-1)] .= xs[i]
    end
    return A
end

function ChainRulesCore.rrule(::typeof(hcat_by_sum), xs::AbstractVector)
    N = length(xs)

    ns = map(x->size(x,2), xs)
    pushfirst!(ns, 1)
    cumsum!(ns, ns)

    hcat_by_sum_pullback(Δ) = (NoTangent(), ntuple(i->view(Δ,:,ns[i]:(ns[i+1]-1)), N))
    hcat_by_sum(xs), hcat_by_sum_pullback
end

function ChainRulesCore.rrule(::typeof(parent), A::Base.SubArray)
    parent_pullback(Δ) = (NoTangent(), view(Δ, A.indices...))
    parent(A), parent_pullback
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3d4c79d0-10ec-4ba1-be3e-ce6cd8f83e9e
md"""
### CHEME 5440/7770: Problem Set 2 (PS2) solution
"""

# ╔═╡ 594f9194-91ef-4a41-999b-a553cca2c180
md"""
##### Build the Stoichiometric Matrix
"""

# ╔═╡ d93fad69-7f2f-420c-ba78-0133fc6e1572
md"""
##### Compute the extreme pathway matrix using expa
"""

# ╔═╡ 82ddb8c4-d286-47e4-9aa4-02f95b6792c9
idx_pathway = 13

# ╔═╡ 6389894e-e9b4-4530-b10c-feedb44862f9
md"""
##### Compute binary stoichiometric array
"""

# ╔═╡ 1e1670b6-f87d-463f-bd38-dbb2e6103653
md"""
##### Metabolite Connectivity Array (MCA)
"""

# ╔═╡ 641b3e7d-e0c8-44f8-bbcd-044bff755046
md"""
##### Reaction Connectivity Array (RCA)
"""

# ╔═╡ 5115eed1-546a-4d87-bfb3-d58ecebe653a
md"""
### Is there a relationship between reaction connectivity and usage?
"""

# ╔═╡ d0f77fde-88d5-11ec-3620-fb46e6aebbca
function ingredients(path::String)

	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol("lib")
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ f7e07c69-f8b3-4956-90bc-bca6dbd41ac5
begin
	# import some packages -
	using PlutoUI
	using PrettyTables
	using LinearAlgebra

	# setup paths -
	const _PATH_TO_NOTEBOOK = pwd()
	const _PATH_TO_SRC = joinpath(_PATH_TO_NOTEBOOK,"src")
	const _PATH_TO_CONFIG = joinpath(_PATH_TO_NOTEBOOK,"config")
	const _PATH_TO_FIGS = joinpath(_PATH_TO_NOTEBOOK, "figs")

	# load the PS2 code lib -
	lib = ingredients(joinpath(_PATH_TO_SRC, "Include.jl"));

	# return -
	nothing
end

# ╔═╡ a7cec06e-8a50-4c34-81a3-0f303099cfd5
PlutoUI.LocalResource(joinpath(_PATH_TO_FIGS,"Fig-Urea-cycle.png"))

# ╔═╡ 32f498cf-a379-4178-9ea2-b022c571227d
begin

	# load/parse the network file -
	_PATH_TO_NETWORK_FILE = joinpath(_PATH_TO_CONFIG,"Network_v2.net")
	list_of_reactions = lib.read_reaction_file(_PATH_TO_NETWORK_FILE);

	# build the stoichiometric array -
	(S, mna, rna) = lib.build_stoichiometric_matrix(list_of_reactions; expand=true);

	# get expanded reaction strings -
	expanded_reaction_array = lib.expand_reversible_reactions(list_of_reactions)

	# return -
	nothing
end

# ╔═╡ db048167-7ca6-4845-889c-6960b61a77dd
(ℳ,ℛ) = size(S)

# ╔═╡ 902d9a8b-2f11-425a-91cd-ccb2345503f1
begin

	# compute the extreme pathways -
	T = lib.expa(S)
	P = T[:,1:ℛ]
	C = T[:,(ℛ + 1):end]

	# show -
	nothing
end

# ╔═╡ 1ac2f945-868c-4262-b903-a30db71ce21d
P

# ╔═╡ b95e014c-e486-4608-a645-b62bb5ac2fb6
begin

	# compute the frequency of reaction usage in P -
	(RP,CP) = size(P)
	reaction_usage_array = Array{Int64,1}(undef,CP)
	for j ∈ 1:CP
		usage = 0
		for i ∈ 1:RP
			if (P[i,j] != 0.0)
				usage += 1
			end
		end
		reaction_usage_array[j] = usage
	end

end

# ╔═╡ a5fb1e84-dc81-49b3-aa60-552e2bebc547
let

	# visualize pathways -
	idx_nz = findall(x->x!=0.0,P[idx_pathway,:])

	# initialize storage for the table -
	state_array = Array{Any,2}(undef, length(idx_nz),3)
	for (i,index) ∈ enumerate(idx_nz)
		state_array[i,1] = i
		state_array[i,2] = P[idx_pathway,index]
		state_array[i,3] = expanded_reaction_array[index]
	end

	# setup the header -
	header_data = (["index","coeff","reaction"])

	with_terminal() do
		pretty_table(state_array; header=header_data, alignment=:l)
	end

end

# ╔═╡ 84c09869-1893-450c-8ff7-e5033077a48c
B = lib.binary_stoichiometric_matrix(S)

# ╔═╡ 9786e10f-b1eb-43f4-8e9f-cce5b8f3af1f
# compute the metabolite connectivity array -
MCA = B*transpose(B)

# ╔═╡ e42e400f-0447-40d5-a0bd-7f86f4a1abfa
begin

	# build a ranking table -
	mca_diagonal = diag(MCA)

	# sort -
	idx_sort = sortperm(mca_diagonal; rev=true)

	# build the mca table -
	state_array = Array{Any,2}(undef,ℳ,3)
	for i ∈ 1:ℳ
		state_array[i,1] = i
		state_array[i,2] = mna[idx_sort[i]]
		state_array[i,3] = mca_diagonal[idx_sort[i]]
	end

	# setup header -
	header_data = (["rank","metabolite","connections"])

	with_terminal() do
		pretty_table(state_array; header = header_data)
	end
end

# ╔═╡ ade866b0-b419-45a5-83f7-4a0058adcc55
# compute the reaction connectivity array -
RCA = transpose(B)*B

# ╔═╡ 273d535a-c28e-451f-9a03-5b5397d13798
let
	# build a ranking table -
	rca_diagonal = diag(RCA)

	# sort -
	idx_sort = sortperm(rca_diagonal; rev=true)

	# build the mca table -
	state_array = Array{Any,2}(undef,ℳ,3)
	for i ∈ 1:ℳ
		state_array[i,1] = i
		state_array[i,2] = rna[idx_sort[i]]
		state_array[i,3] = rca_diagonal[idx_sort[i]]
	end

	# setup header -
	header_data = (["rank","reaction","connections"])

	with_terminal() do
		pretty_table(state_array; header = header_data)
	end
end

# ╔═╡ 250b0ed2-e465-4e4c-b2eb-4d5ae80387cc
let

	# build usage table -
	state_table = Array{Any,2}(undef, CP, 5)
	idx_sort = sortperm(reaction_usage_array; rev=true)
	for i ∈ 1:CP
		state_table[i,1] = i
		state_table[i,2] = idx_sort[i]
		state_table[i,3] = reaction_usage_array[idx_sort[i]]
		state_table[i,4] = (1/RP)*reaction_usage_array[idx_sort[i]]
		state_table[i,5] = expanded_reaction_array[idx_sort[i]]
	end

	# setup header -
	header_data = (["rank","index","usage","freq", "reaction"])

	with_terminal() do
		pretty_table(state_table; header=header_data, alignment=:l)
	end
end

# ╔═╡ b0f24d05-af37-477d-a709-b011cba825fc
html"""
<style>
main {
    max-width: 900px;
    width: 75%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}
</style>"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"

[compat]
PlutoUI = "~0.7.34"
PrettyTables = "~1.3.1"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "8979e9802b4ac3d58c503a20f2824ad67f9074dd"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.34"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "bb1064c9a84c52e277f1096cf41434b675cd368b"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.6.1"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─3d4c79d0-10ec-4ba1-be3e-ce6cd8f83e9e
# ╠═a7cec06e-8a50-4c34-81a3-0f303099cfd5
# ╟─594f9194-91ef-4a41-999b-a553cca2c180
# ╠═32f498cf-a379-4178-9ea2-b022c571227d
# ╠═db048167-7ca6-4845-889c-6960b61a77dd
# ╟─d93fad69-7f2f-420c-ba78-0133fc6e1572
# ╠═902d9a8b-2f11-425a-91cd-ccb2345503f1
# ╠═1ac2f945-868c-4262-b903-a30db71ce21d
# ╠═82ddb8c4-d286-47e4-9aa4-02f95b6792c9
# ╠═a5fb1e84-dc81-49b3-aa60-552e2bebc547
# ╟─6389894e-e9b4-4530-b10c-feedb44862f9
# ╠═84c09869-1893-450c-8ff7-e5033077a48c
# ╟─1e1670b6-f87d-463f-bd38-dbb2e6103653
# ╠═9786e10f-b1eb-43f4-8e9f-cce5b8f3af1f
# ╟─e42e400f-0447-40d5-a0bd-7f86f4a1abfa
# ╟─641b3e7d-e0c8-44f8-bbcd-044bff755046
# ╠═ade866b0-b419-45a5-83f7-4a0058adcc55
# ╟─273d535a-c28e-451f-9a03-5b5397d13798
# ╟─5115eed1-546a-4d87-bfb3-d58ecebe653a
# ╠═b95e014c-e486-4608-a645-b62bb5ac2fb6
# ╠═250b0ed2-e465-4e4c-b2eb-4d5ae80387cc
# ╠═f7e07c69-f8b3-4956-90bc-bca6dbd41ac5
# ╠═d0f77fde-88d5-11ec-3620-fb46e6aebbca
# ╠═b0f24d05-af37-477d-a709-b011cba825fc
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using Pkg
Pkg.activate(".")

using HTTP

##########################

init() = HTTP.Router()

function route(resp_fn::Function, app::HTTP.Router, route::String, method::Union{String,Symbol} = :GET)
    HTTP.@register(app,method,route,resp_fn)
end

function run(app::HTTP.Router, host = "127.0.0.1", port = 8081)
    @info "Serving at http://$host:$port ..."
    HTTP.serve(app,host,port)
end

##########################

const APP = init()

function f(r::HTTP.Request)
    return HTTP.Response(200,"{\"Hello\":\"World!\"}")
end

route(f,APP,"/",:GET)

run(APP)


text:
# NOTE: worker processes cannot add more workers, only the client process can.
require("testdefs.jl")

if nprocs() < 2
    remotecall_fetch(1, () -> addprocs(1))
end

id_me = myid()
id_other = filter(x -> x != id_me, procs())[rand(1:(nprocs()-1))]

@test fetch(@spawnat id_other myid()) == id_other
@test @fetchfrom id_other begin myid() end == id_other
@fetch begin myid() end

d = drand((200,200), [id_me, id_other])
s = convert(Matrix{Float64}, d[1:150, 1:150])
a = convert(Matrix{Float64}, d)
@test a[1:150,1:150] == s

@test fetch(@spawnat id_me localpart(d)[1,1]) == d[1,1]
@test fetch(@spawnat id_other localpart(d)[1,1]) == d[1,101]

d=DArray(I->fill(myid(), map(length,I)), (10,10), [id_me, id_other])
d2 = map(x->1, d)
@test reduce(+, d2) == 100

@test reduce(+, d) == ((50*id_me) + (50*id_other))
map!(x->1, d)
@test reduce(+, d) == 100

@unix_only begin

dims = (20,20,20)

@linux_only begin
    S = SharedArray(Int64, dims)
    @test beginswith(S.segname, "/jl")
    @test !ispath("/dev/shm" * S.segname)

    S = SharedArray(Int64, dims; pids=[id_other])
    @test beginswith(S.segname, "/jl")
    @test !ispath("/dev/shm" * S.segname)
end

# TODO : Need a similar test of shmem cleanup for OSX

# SharedArray tests
d = Base.shmem_rand(1:100, dims)
a = convert(Array, d)

partsums = Array(Int, length(procs(d)))
@sync begin
    for (i, p) in enumerate(procs(d))
        @async partsums[i] = remotecall_fetch(p, D->sum(D.loc_subarr_1d), d)
    end
end
@test sum(a) == sum(partsums)

d = Base.shmem_rand(dims)
for p in procs(d)
    idxes_in_p = remotecall_fetch(p, D -> parentindexes(D.loc_subarr_1d)[1], d)
    idxf = first(idxes_in_p)
    idxl = last(idxes_in_p)
    d[idxf] = float64(idxf)
    rv = remotecall_fetch(p, (D,idxf,idxl) -> begin assert(D[idxf] == float64(idxf)); D[idxl] = float64(idxl); D[idxl];  end, d,idxf,idxl)
    @test d[idxl] == rv
end

@test ones(10, 10, 10) == Base.shmem_fill(1.0, (10,10,10))
@test zeros(Int32, 10, 10, 10) == Base.shmem_fill(0, (10,10,10))

d = Base.shmem_rand(dims)
s = Base.shmem_rand(dims)
copy!(s, d)
@test s == d
s = Base.shmem_rand(dims)
copy!(s, sdata(d))
@test s == d

d = SharedArray(Int, dims; init = D->fill!(D.loc_subarr_1d, myid()))
for p in procs(d)
    idxes_in_p = remotecall_fetch(p, D -> parentindexes(D.loc_subarr_1d)[1], d)
    idxf = first(idxes_in_p)
    idxl = last(idxes_in_p)
    @test d[idxf] == p
    @test d[idxl] == p
end

# issue #6362
d = Base.shmem_rand(dims)
s = copy(sdata(d))
ds = deepcopy(d)
@test ds == d
remotecall_fetch(findfirst(id->(id != myid()), procs(ds)), setindex!, ds, 1.0, 1:10)
@test ds != d
@test s == d

# SharedArray as an array
# Since the data in d will depend on the nprocs, just test that these operations work
a = d[1:5]
@test_throws BoundsError d[-1:5]
a = d[1,1,1:3:end]
d[2:4] = 7
d[5,1:2:4,8] = 19

AA = rand(4,2)
A = convert(SharedArray, AA)
B = convert(SharedArray, AA')
@test B*A == AA'*AA

d=SharedArray(Int64, (10,10); init = D->fill!(D.loc_subarr_1d, myid()), pids=[id_me, id_other])
d2 = map(x->1, d)
@test reduce(+, d2) == 100

@test reduce(+, d) == ((50*id_me) + (50*id_other))
map!(x->1, d)
@test reduce(+, d) == 100

# Boundary cases where length(S) <= length(pids)
@test 2.0 == remotecall_fetch(id_other, D->D[2], Base.shmem_fill(2.0, 2; pids=[id_me, id_other]))
@test 3.0 == remotecall_fetch(id_other, D->D[1], Base.shmem_fill(3.0, 1; pids=[id_me, id_other]))

end # @unix_only(SharedArray tests)

# Test @parallel load balancing - all processors should get either M or M+1
# iterations out of the loop range for some M.
if nprocs() < 4
    remotecall_fetch(1, () -> addprocs(4 - nprocs()))
end
workloads = hist(@parallel((a,b)->[a,b], for i=1:7; myid(); end), nprocs())[2]
@test maximum(workloads) - minimum(workloads) <= 1

# @parallel reduction should work even with very short ranges
@test @parallel(+, for i=1:2; i; end) == 3

# Testing timedwait on multiple RemoteRefs
rr1 = RemoteRef()
rr2 = RemoteRef()
rr3 = RemoteRef()

@async begin sleep(0.5); put!(rr1, :ok) end
@async begin sleep(1.0); put!(rr2, :ok) end
@async begin sleep(2.0); put!(rr3, :ok) end

tic()
timedwait(1.0) do
    all(map(isready, [rr1, rr2, rr3]))
end
et=toq()

# assuming that 0.5 seconds is a good enough buffer on a typical modern CPU
try
    @test (et >= 1.0) && (et <= 1.5)
    @test !isready(rr3)
catch
    warn("timedwait tests delayed. et=$et, isready(rr3)=$(isready(rr3))")
end
@test isready(rr1)

# TODO: The below block should be always enabled but the error is printed by the event loop

# Hence in the event of any relevant changes to the parallel codebase,
# please define an ENV variable PTEST_FULL and ensure that the below block is
# executed successfully before committing/merging

if haskey(ENV, "PTEST_FULL")
    println("START of parallel tests that print errors")

    # make sure exceptions propagate when waiting on Tasks
    @test_throws ErrorException (@sync (@async error("oops")))

    # pmap tests
    # needs at least 4 processors (which are being created above for the @parallel tests)
    s = "a"*"bcdefghijklmnopqrstuvwxyz"^100;
    ups = "A"*"BCDEFGHIJKLMNOPQRSTUVWXYZ"^100;
    @test ups == bytestring(Uint8[uint8(c) for c in pmap(x->uppercase(x), s)])
    @test ups == bytestring(Uint8[uint8(c) for c in pmap(x->uppercase(char(x)), s.data)])

    # retry, on error exit
    res = pmap(x->(x=='a') ? error("test error. don't panic.") : uppercase(x), s; err_retry=true, err_stop=true);
    @test length(res) < length(ups)
    @test isa(res[1], Exception)

    # no retry, on error exit
    res = pmap(x->(x=='a') ? error("test error. don't panic.") : uppercase(x), s; err_retry=false, err_stop=true);
    @test length(res) < length(ups)
    @test isa(res[1], Exception)

    # retry, on error continue
    res = pmap(x->iseven(myid()) ? error("test error. don't panic.") : uppercase(x), s; err_retry=true, err_stop=false);
    @test length(res) == length(ups)
    @test ups == bytestring(Uint8[uint8(c) for c in res])

    # no retry, on error continue
    res = pmap(x->(x=='a') ? error("test error. don't panic.") : uppercase(x), s; err_retry=false, err_stop=false);
    @test length(res) == length(ups)
    @test isa(res[1], Exception)

    println("END of parallel tests that print errors")
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
struct CapacityMarket # Assumes a piece-wise linear demand curve

    breakpoints::Vector{Float64}   # demand curve break points
    pricepoints::Vector{Float64}    # $/MW-year

    function CapacityMarket(b, p)
        @assert all(b .>= 0)
        @assert all(p .>= 0)
        @assert length(b) > 1
        @assert length(p) == length(b)
        new(b, p)
    end

end


text:
import Knet.Ops20: conv4, conv4x, conv4w, pool, poolx
using CUDA: CUDA, CUDNN, Mem
using Knet.KnetArrays: DevArray, Cptr

function conv4(w::R,x::R; handle=CUDNN.handle(), alpha=1,
               o...) where {T,R<:DevArray{T}} # padding=0, stride=1, dilation=1, mode=0, group=1
    beta=0 # nonzero beta does not make sense when we create y
    y = similar(x, cdims(w,x;o...))
    (algo,workSpace) = conv4_algo(w, x, y; handle=handle, o...)
    CUDNN.cudnnConvolutionForward(handle,Ref(T(alpha)),TD(x),x,FD(w),w,CD(w,x;o...),algo,workSpace,bytes(workSpace),Ref(T(beta)),TD(y),y)
    return y
end

function conv4x(w::R,x::R,dy::R; handle=CUDNN.handle(), alpha=1,
                   o...) where {T,R<:DevArray{T}} # padding=0, stride=1, dilation=1, mode=0, group=1
    beta = 0
    dx = similar(x)
    (algo,workSpace) = conv4x_algo(w,x,dy,dx; handle=handle, o...)
    CUDNN.cudnnConvolutionBackwardData(handle,Ref(T(alpha)),FD(w),w,TD(dy),dy,CD(w,x;o...),algo,workSpace,bytes(workSpace),Ref(T(beta)),TD(dx),dx)
    return dx
end

function conv4w(w::R,x::R,dy::R; handle=CUDNN.handle(), alpha=1,
                   o...) where {T,R<:DevArray{T}} # padding=0, stride=1, dilation=1, mode=0, group=1
    beta = 0
    dw = similar(w)
    (algo,workSpace) = conv4w_algo(w,x,dy,dw;handle=handle,o...)
    CUDNN.cudnnConvolutionBackwardFilter(handle,Ref(T(alpha)),TD(x),x,TD(dy),dy,CD(w,x;o...),algo,workSpace,bytes(workSpace),Ref(T(beta)),FD(dw),dw)
    return dw
end

function pool(x::R; handle=CUDNN.handle(), alpha=1,
                 o...) where {T,R<:DevArray{T}} # window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0
    y = similar(x, pdims(x; o...))
    beta = 0
    CUDNN.cudnnPoolingForward(handle,PD(x;o...),Ref(T(alpha)),TD(x),x,    Ref(T(beta)),TD(y),y)
    return y
end

function poolx(x::R,y::R,dy::R; handle=CUDNN.handle(), alpha=1, mode=0,
                  o...) where {T,R<:DevArray{T}} # window=2, padding=0, stride=window, maxpoolingNanOpt=0
    dx = similar(x)
    beta = 0
    CUDNN.cudnnPoolingBackward(handle,PD(x;mode=mode,o...),Ref(T(alpha)),TD(y),y,TD(dy),dy,TD(x),x,Ref(T(beta)),TD(dx),dx)
    return dx
end

# cudnn descriptors: need to use mutable structs in order to have finalizers.

mutable struct TD; ptr; end
TD(a::DevArray{T}) where {T} = TD(T,size(a))
TD(T::Type, dims::Integer...) = TD(T, dims)
function TD(T::Type, dims)
    d = Cptr[0]
    CUDNN.cudnnCreateTensorDescriptor(d)
    n = length(dims)
    sz = [Cint(dims[i]) for i=n:-1:1]
    st = similar(sz); st[n] = 1
    for i=(n-1):-1:1; st[i] = st[i+1] * sz[i+1]; end
    dt = CUDNN.cudnnDataType_t(DT(T))
    CUDNN.cudnnSetTensorNdDescriptor(d[1], dt, n, sz, st)
    td = TD(d[1])
    finalizer(x->CUDNN.cudnnDestroyTensorDescriptor(x.ptr), td)
    return td
end

mutable struct FD; ptr; end
FD(a::DevArray{T}) where {T}=FD(T,size(a))
FD(T::Type, dims::Integer...) = FD(T,dims)
function FD(T::Type, dims)
    d = Cptr[0]
    CUDNN.cudnnCreateFilterDescriptor(d)
    n = length(dims)
    sz = [Cint(dims[i]) for i=n:-1:1]
    dt = CUDNN.cudnnDataType_t(DT(T))
    tf = CUDNN.cudnnTensorFormat_t(0)
    CUDNN.cudnnSetFilterNdDescriptor(d[1], dt, tf, n, sz)
    fd = FD(d[1])
    finalizer(x->CUDNN.cudnnDestroyFilterDescriptor(x.ptr), fd)
    return fd
end

mutable struct CD; ptr; end
function CD(w::DevArray,x::DevArray; padding=0, stride=1, dilation=1, mode=0, group=1)
    d = Cptr[0]
    CUDNN.cudnnCreateConvolutionDescriptor(d)
    cd = CD(d[1])
    nd = ndims(x)-2
    dt = CUDNN.cudnnDataType_t(DT(x))
    mode = CUDNN.cudnnConvolutionMode_t(mode)
    CUDNN.cudnnSetConvolutionNdDescriptor(cd,nd,cdsize(padding,nd),cdsize(stride,nd),cdsize(dilation,nd),mode,dt)
    CUDNN.cudnnSetConvolutionGroupCount(cd,group)
    finalizer(x->CUDNN.cudnnDestroyConvolutionDescriptor(x.ptr),cd)
    return cd
end

mutable struct PD; ptr; end
function PD(x::DevArray; window=2, padding=0, stride=window, mode=0, maxpoolingNanOpt=0)
    d = Cptr[0]
    CUDNN.cudnnCreatePoolingDescriptor(d)
    nd = ndims(x)-2
    mode = CUDNN.cudnnPoolingMode_t(mode)
    maxpoolingNanOpt = CUDNN.cudnnNanPropagation_t(maxpoolingNanOpt)
    CUDNN.cudnnSetPoolingNdDescriptor(d[1],mode,maxpoolingNanOpt,nd,cdsize(window,nd),cdsize(padding,nd),cdsize(stride,nd))
    pd = PD(d[1])
    finalizer(x->CUDNN.cudnnDestroyPoolingDescriptor(x.ptr), pd)
    return pd
end

import Base: unsafe_convert
unsafe_convert(::Type{Cptr}, td::TD)=td.ptr
unsafe_convert(::Type{Cptr}, fd::FD)=fd.ptr
unsafe_convert(::Type{Cptr}, cd::CD)=cd.ptr
unsafe_convert(::Type{Cptr}, pd::PD)=pd.ptr

# fill and reverse Cint array with padding etc. for cudnn calls
function cdsize(w, nd)
    if isa(w,Number)
        fill(Cint(w),nd)
    elseif length(w)==nd
        [ Cint(w[nd-i+1]) for i=1:nd ]
    else
        throw(DimensionMismatch("$w $nd"))
    end
end

DT(::DevArray{Float32})=Cint(0)
DT(::DevArray{Float64})=Cint(1)
DT(::DevArray{Float16})=Cint(2)
DT(::Type{Float32}) = Cint(0)
DT(::Type{Float64}) = Cint(1)
DT(::Type{Float16}) = Cint(2)

## Utilities to find a fast algorithm

struct cudnnConvolutionFwdAlgoPerf_t
    algo::Cint
    status::Cint
    time::Cfloat
    memory::Csize_t
    determinism::Cint
    mathType::Cint
    r1::Cint; r2::Cint; r3::Cint
end

const CUDNN_MAX_FIND = 100      # How many times can we call FindAlgorithm
const requestedAlgoCount = 10
const returnedAlgoCount = Cint[0]
bytes(x::DevArray{T}) where {T}=length(x)*sizeof(T)
gpufree() = Mem.info()[1] + (isdefined(CUDA,:pool) ? CUDA.pool[].cached_memory() : CUDA.cached_memory())

# This seems to cover a reasonable subset of the available algorithms
# The user can set this to 0 for a more memory-tight execution
maxWorkspaceSize(w,x,y) = min(gpufree() ÷ 10, bytes(x) * 100)

const conv4_algos = Dict()
function conv4_algo(w::R, x::R, y::R; handle=CUDNN.handle(), o...) where {T,R<:DevArray{T}}
    key = (T,size(w),size(x),o...)
    if haskey(conv4_algos, key)
        p = conv4_algos[key]
    elseif length(conv4_algos) >= CUDNN_MAX_FIND
        p = nothing
    else
        workSpace = similar(w, maxWorkspaceSize(w,x,y) ÷ sizeof(T))
        perfResults = Array{CUDNN.cudnnConvolutionFwdAlgoPerf_t}(undef,requestedAlgoCount)
        wd, xd, yd, cd = FD(w), TD(x), TD(y), CD(w,x;o...)
        CUDNN.cudnnFindConvolutionForwardAlgorithmEx(handle,xd,x,wd,w,cd,yd,y,requestedAlgoCount,returnedAlgoCount,perfResults,workSpace,bytes(workSpace))
        p = perfChoose(perfResults, returnedAlgoCount[1])
        conv4_algos[key] = p
        if p === nothing
            @warn "No good algo found for conv4$o: using default algo=0." maxlog=1
        end
    end
    if p === nothing
        return (CUDNN.cudnnConvolutionFwdAlgo_t(0), cudnnWorkSpace(w))
    else
        return (p.algo, cudnnWorkSpace(w,p.memory))
    end
end

const conv4w_algos = Dict()
function conv4w_algo(w::R,x::R,dy::R,dw::R; handle=CUDNN.handle(), o...) where {T,R<:DevArray{T}}
    key = (T,size(w),size(x),o...)
    if haskey(conv4w_algos, key)
        p = conv4w_algos[key]
    elseif length(conv4w_algos) >= CUDNN_MAX_FIND
        p = nothing
    else
        workSpace = similar(w, maxWorkspaceSize(w,x,dy) ÷ sizeof(T))
        perfResults = Array{CUDNN.cudnnConvolutionBwdFilterAlgoPerf_t}(undef,requestedAlgoCount)
        wd, xd, yd, cd = FD(dw), TD(x), TD(dy), CD(w,x;o...)
        CUDNN.cudnnFindConvolutionBackwardFilterAlgorithmEx(handle,xd,x,yd,dy,cd,wd,dw,requestedAlgoCount,returnedAlgoCount,perfResults,workSpace,bytes(workSpace))
        p = perfChoose(perfResults, returnedAlgoCount[1])
        conv4w_algos[key] = p
        if p === nothing
            @warn "No good algo found for conv4w$o: using default algo=0." maxlog=1
        end
    end
    if p === nothing
        return (CUDNN.cudnnConvolutionBwdFilterAlgo_t(0), cudnnWorkSpace(w))
    else
        return (p.algo, cudnnWorkSpace(w, p.memory))
    end
end

const conv4x_algos = Dict()
function conv4x_algo(w::R,x::R,dy::R,dx::R; handle=CUDNN.handle(), o...) where {T,R<:DevArray{T}}
    key = (T,size(w),size(x),o...)
    if haskey(conv4x_algos, key)
        p = conv4x_algos[key]
    elseif length(conv4x_algos) >= CUDNN_MAX_FIND
        p = nothing
    else
        workSpace = similar(w, maxWorkspaceSize(w,x,dy) ÷ sizeof(T))
        perfResults = Array{CUDNN.cudnnConvolutionBwdDataAlgoPerf_t}(undef,requestedAlgoCount)
        wd, xd, yd, cd = FD(w), TD(dx), TD(dy), CD(w,x;o...)
        CUDNN.cudnnFindConvolutionBackwardDataAlgorithmEx(handle,wd,w,yd,dy,cd,xd,dx,requestedAlgoCount,returnedAlgoCount,perfResults,workSpace,bytes(workSpace))
        p = perfChoose(perfResults, returnedAlgoCount[1])
        conv4x_algos[key] = p
        if p === nothing
            @warn "No good algo found for conv4x$o: using default algo=0." maxlog=1
        end
    end
    if p === nothing
        return (CUDNN.cudnnConvolutionBwdDataAlgo_t(0), cudnnWorkSpace(w))
    else
        return (p.algo, cudnnWorkSpace(w, p.memory))
    end
end

function perfChoose(ps, n)
    if n==length(ps)
        @warn "returnedAlgoCount==requestedAlgoCount" maxlog=1
    end
    (ibest,mbest,tbest) = (0,Inf,Inf)
    for i = 1:n
        # These metrics are written in a sorted fashion where the first element has the lowest compute time.
        if ps[i].status == 0 && ps[i].memory < mbest && ps[i].time < tbest * 1.1
            (ibest,mbest,tbest) = (i,ps[i].memory,ps[i].time)
        end
    end
    if ibest > 0
        return ps[ibest]
    else
        return nothing
    end
end

# Fresh workspace for every op is safer:
cudnnWorkSpace(w)=similar(w, 0)
cudnnWorkSpace(w,len)=similar(w, 1 + (len-1) ÷ sizeof(eltype(w)))

## Dimension helpers:

# outputDim = 1 + ( inputDim + 2*pad - (((filterDim-1)*dilation)+1) )/convolutionStride;

function cdims(w,x; padding=0, stride=1, dilation=1, o...)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            di = (if isa(dilation,Number); dilation; else dilation[i]; end)
            1 + (size(x,i) + 2*pi - (((size(w,i)-1)*di)+1)) ÷ si
        elseif i == N-1
            size(w,N)
        else # i == N
            size(x,N)
        end
    end
end

function pdims(x; window=2, padding=0, stride=window, o...)
    N = ndims(x)
    ntuple(N) do i
        if i < N-1
            wi = (if isa(window,Number); window; else window[i]; end)
            pi = (if isa(padding,Number); padding; else padding[i]; end)
            si = (if isa(stride,Number); stride; else stride[i]; end)
            1 + div(size(x,i) + 2*pi - wi, si)
        else
            size(x,i)
        end
    end
end


text:
function applyhemoics!(system::CVSystem,old=Dict("a"=>0))
    for i = 1:length(system.branches.ID)
        temp = false;
        for j = 1:length(system.hemo.hID)
            if system.branches.ID[i] == system.hemo.hID[j]
                temp = true;
            end
        end
        push!(system.hemo.injured,temp);
    end
    if haskey(old,"hemo")
        hemo = old["hemo"];
        if hemo["Vloss"] == 0
            for j = 1:length(system.hemo.tID)
                push!(system.hemo.bmin,system.branches.beta[system.hemo.tID[j]][end]);
                push!(system.hemo.bmax,2*system.branches.beta[system.hemo.tID[j]][end]);
                push!(system.hemo.Amax,system.branches.A0[system.hemo.tID[j]][end]);
                push!(system.hemo.Amin,0.01*system.branches.A0[system.hemo.tID[j]][end]);
            end
        elseif hemo["Vloss"] > 0
            for j = 1:length(system.hemo.tID)
                push!(system.hemo.bmin,hemo["bmin"][j]);
                push!(system.hemo.bmax,hemo["bmax"][j]);
                push!(system.hemo.Amin,hemo["Amin"][j]);
                push!(system.hemo.Amax,hemo["Amax"][j]);
            end
            system.hemo.Vloss = hemo["Vloss"];
            system.hemo.Vlossinit = system.hemo.Vloss;
        end
    else
        for j = 1:length(system.hemo.tID)
            push!(system.hemo.bmin,system.branches.beta[system.hemo.tID[j]][end]);
            push!(system.hemo.bmax,2*system.branches.beta[system.hemo.tID[j]][end]);
            push!(system.hemo.Amax,system.branches.A0[system.hemo.tID[j]][end]);
            push!(system.hemo.Amin,0.01*system.branches.A0[system.hemo.tID[j]][end]);
        end
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
module ITensorChainRules

using ChainRulesCore
using ITensors

include("zygoterules.jl")

function ChainRulesCore.rrule(::typeof(getindex), x::ITensor, I...)
  y = getindex(x, I...)
  function getindex_pullback(ȳ)
    # TODO: add definition `ITensor(::Tuple{}) = ITensor()`
    # to ITensors.jl so no splatting is needed here.
    x̄ = ITensor(inds(x)...)
    x̄[I...] = ȳ
    Ī = broadcast_notangent(I)
    return (NoTangent(), x̄, Ī...)
  end
  return y, getindex_pullback
end

# Specialized version in order to avoid call to `setindex!`
# within the pullback, should be better for taking higher order
# derivatives in Zygote.
function ChainRulesCore.rrule(::typeof(getindex), x::ITensor)
  y = x[]
  function getindex_pullback(ȳ)
    x̄ = ITensor(ȳ)
    return (NoTangent(), x̄)
  end
  return y, getindex_pullback
end

function setinds_pullback(ȳ, x, a...)
  x̄ = ITensors.setinds(ȳ, inds(x))
  ā = broadcast_notangent(a)
  return (NoTangent(), x̄, ā...)
end

function inv_op(f::Function, args...; kwargs...)
  return error(
    "Trying to differentiate `$f` but the inverse of the operation (`inv_op`) `$f` with arguments $args and keyword arguments $kwargs is not defined.",
  )
end

function inv_op(::typeof(prime), x::ITensor, n::Integer=1; kwargs...)
  return prime(x, -n; kwargs...)
end

function inv_op(::typeof(replaceprime), x::ITensor, n1n2::Pair; kwargs...)
  return replaceprime(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(swapprime), x::ITensor, n1n2::Pair; kwargs...)
  return swapprime(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(addtags), x::ITensor, args...; kwargs...)
  return removetags(x, args...; kwargs...)
end

function inv_op(::typeof(removetags), x::ITensor, args...; kwargs...)
  return addtags(x, args...; kwargs...)
end

function inv_op(::typeof(replacetags), x::ITensor, n1n2::Pair; kwargs...)
  return replacetags(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(swaptags), x::ITensor, n1n2::Pair; kwargs...)
  return swaptags(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(replaceind), x::ITensor, n1n2::Pair; kwargs...)
  return replaceind(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(replaceinds), x::ITensor, n1n2::Pair; kwargs...)
  return replaceinds(x, reverse(n1n2); kwargs...)
end

function inv_op(::typeof(swapind), x::ITensor, args...; kwargs...)
  return swapind(x, reverse(args)...; kwargs...)
end

function inv_op(::typeof(swapinds), x::ITensor, args...; kwargs...)
  return swapinds(x, reverse(args)...; kwargs...)
end

for fname in (
  :prime,
  :setprime,
  :noprime,
  :replaceprime,
  :swapprime,
  :addtags,
  :removetags,
  :replacetags,
  :settags,
  :swaptags,
  :replaceind,
  :replaceinds,
  :swapind,
  :swapinds,
)
  @eval begin
    function ChainRulesCore.rrule(f::typeof($fname), x::ITensor, a...; kwargs...)
      y = f(x, a...; kwargs...)
      function f_pullback(ȳ)
        x̄ = inv_op(f, ȳ, a...; kwargs...)
        if !hassameinds(x, x̄)
          error(
            "Trying to differentiate function `$f` with arguments $a and keyword arguments $kwargs. The forward pass indices $(inds(x)) do not match the reverse pass indices $(inds(x̄)). Likely this is because the priming/tagging operation you tried to perform is not invertible. Please write your code in a way where the index manipulation operation you are performing is invertible. For example, `prime(A::ITensor)` is invertible, with an inverse `prime(A, -1)`. However, `noprime(A)` is in general not invertible since the information about the prime levels of the original tensor are lost. Instead, you might try `prime(A, -1)` or `replaceprime(A, 1 => 0)` which are invertible.",
          )
        end
        ā = broadcast_notangent(a)
        return (NoTangent(), x̄, ā...)
      end
      return y, f_pullback
    end
  end
end

# TODO: This is not being called by Zygote for some reason,
# using a Zygote overload directly instead. Figure out
# why, maybe raise an issue.
#function ChainRulesCore.rrule(::typeof(adjoint), x::ITensor)
#  y = prime(x)
#  function adjoint_pullback(ȳ)
#    return setinds_pullback(ȳ, x)
#  end
#  return y, adjoint_pullback
#end

# Special case for contracting a pair of ITensors
function ChainRulesCore.rrule(::typeof(*), x1::ITensor, x2::ITensor)
  y = x1 * x2
  function contract_pullback(ȳ)
    x̄1 = ȳ * dag(x2)
    x̄2 = dag(x1) * ȳ
    return (NoTangent(), x̄1, x̄2)
  end
  return y, contract_pullback
end

function ChainRulesCore.rrule(::typeof(*), x1::Number, x2::ITensor)
  y = x1 * x2
  function contract_pullback(ȳ)
    x̄1 = ȳ * dag(x2)
    # TODO: define `dag(::Number)` to help with generic code
    x̄2 = dag(ITensor(x1)) * ȳ
    return (NoTangent(), x̄1[], x̄2)
  end
  return y, contract_pullback
end

function ChainRulesCore.rrule(::typeof(*), x1::ITensor, x2::Number)
  y = x1 * x2
  function contract_pullback(ȳ)
    x̄1 = ȳ * dag(ITensor(x2))
    x̄2 = dag(x1) * ȳ
    return (NoTangent(), x̄1, x̄2[])
  end
  return y, contract_pullback
end

function ChainRulesCore.rrule(::typeof(*), x1::ITensor, x2::ITensor, xs::ITensor...)
  y = *(x1, x2, xs...)
  function contract_pullback(ȳ)
    # TODO: use some contraction sequence optimization here
    tn = [x1, x2, xs...]
    N = length(tn)
    env_contracted = Vector{ITensor}(undef, N)
    for n in 1:length(tn)
      tn_left = tn[1:(n - 1)]
      # TODO: define contract([]) = ITensor(1.0)
      env_left = isempty(tn_left) ? ITensor(1.0) : contract(tn_left)
      tn_right = tn[reverse((n + 1):end)]
      env_right = isempty(tn_right) ? ITensor(1.0) : contract(tn_right)
      env_contracted[n] = dag(env_left) * ȳ * dag(env_right)
    end
    return (NoTangent(), env_contracted...)
  end
  return y, contract_pullback
end

function ChainRulesCore.rrule(::typeof(+), x1::ITensor, x2::ITensor)
  y = x1 + x2
  function add_pullback(ȳ)
    return (NoTangent(), ȳ, ȳ)
  end
  return y, add_pullback
end

function ChainRulesCore.rrule(::typeof(itensor), x::Array, a...)
  y = itensor(x, a...)
  function itensor_pullback(ȳ)
    x̄ = array(ȳ)
    ā = broadcast_notangent(a)
    return (NoTangent(), x̄, ā...)
  end
  return y, itensor_pullback
end

function ChainRulesCore.rrule(::typeof(ITensor), x::Array{<:Number}, a::Index...)
  y = ITensor(x, a...)
  function ITensor_pullback(ȳ)
    # TODO: define `Array(::ITensor)` directly
    x̄ = Array(ȳ, inds(ȳ)...)
    ā = broadcast(_ -> NoTangent(), a)
    return (NoTangent(), x̄, ā...)
  end
  return y, ITensor_pullback
end

function ChainRulesCore.rrule(::typeof(ITensor), x::Number)
  y = ITensor(x)
  function ITensor_pullback(ȳ)
    x̄ = ȳ[]
    return (NoTangent(), x̄)
  end
  return y, ITensor_pullback
end

function ChainRulesCore.rrule(::typeof(dag), x::ITensor)
  y = dag(x)
  function dag_pullback(ȳ)
    x̄ = dag(ȳ)
    return (NoTangent(), x̄)
  end
  return y, dag_pullback
end

function ChainRulesCore.rrule(::typeof(permute), x::ITensor, a...)
  y = permute(x, a...)
  function permute_pullback(ȳ)
    x̄ = permute(ȳ, inds(x))
    ā = broadcast_notangent(a)
    return (NoTangent(), x̄, ā...)
  end
  return y, permute_pullback
end

broadcast_notangent(a) = broadcast(_ -> NoTangent(), a)

@non_differentiable broadcast_notangent(::Any)
@non_differentiable Index(::Any...)
@non_differentiable delta(::Any...)
@non_differentiable dag(::Index)
@non_differentiable inds(::Any...)

end


text:
module TheAlgorithms

# Usings/Imports (keep sorted)
using DataFrames
using DifferentialEquations
using GLM
using LinearAlgebra
using Plots
using Random

## Exports
# Please keep the folders/functions sorted

# Exports: math
export collatz_sequence
export euler_method
export line_length
export SIR # TODO: make the name lowercase if possible
export trapazoidal_area
export is_armstrong

# Exports: matrix
export determinant
export lu_decompose

# Exports: project-rosalind
export count_nucleotides
export dna2rna
export reverse_complement

# Exports: searches
export binarysearch
export search

# Exports: statistics
export OLSbeta

## Includes
# Please keep the folders/files sorted (by dependencies then alphabetical order)

# Includes: math
include("math/area_under_curve.jl")
include("math/collatz_sequence.jl")
include("math/euler_method.jl")
include("math/line_length.jl")
include("math/sir_model.jl")
include("math/armstrong_number.jl")

# Includes: matrix
include("matrix/lu_decompose.jl") # used by determinant.jl
include("matrix/determinant.jl")

# Includes: project-rosalind
include("project-rosalind/count_nucleotide.jl")
include("project-rosalind/dna2rna.jl")
include("project-rosalind/reverse_complement.jl")

# Includes: search
include("searches/binary_search.jl")
include("searches/linear_search.jl")

# Includes: statistics
include("statistics/ordinary_least_squares.jl")

end

max_stars_repo_name: madisonestabrook/Julia

max_issues_repo_name: madisonestabrook/Julia

max_forks_repo_name: madisonestabrook/Julia


text:
using BenchmarkTools
H = read_hypergraph("/Users/aaronkelley/.julia/dev/CoexistenceHoles.jl/test/test_data/Lagoon of Venice/Hp.dat")

@benchmark begin
    f = betti_hypergraph_ripscomplex(H; max_dim=4)
end

@benchmark begin
    hypergraph_subdivide(H)
end;

@time begin
    G = hypergraph_subdivide(H; expansion = false);
end;

@time begin
    M = 2.0*ones(length(H), length(H));
    [ M[g[1], g[2]] = 1 for g in G ];
    [ M[g[2], g[1]] = 1 for g in G ];
    M[diagind(M)] .= 0;
end;

@time begin
    flt = Rips(M; threshold=1.5,sparse=true)
end;

@time begin
    barcodes = ripserer(flt; dim_max = 4)
end;

@time begin
    betti = [ count(x -> x[2] .== Inf, barcodes[i]) for i in 1:length(barcodes)]
end;

function betti_hypergraph_ripscomplex(H; max_dim = 3)

    if length(H) == 0
        return Int.(zeros(max_dim + 1))
    elseif length(H) == 1
        return Int.(vcat(1, zeros(max_dim)))
    end

    _max_dim = min(max_dim, length(H) - 2);

    G = hypergraph_subdivide(H; expansion = false)
    # build distance matrix to give as input to Ripserer.
    M = 2.0*ones(length(H), length(H))
    [ M[g[1], g[2]] = 1 for g in G ]
    [ M[g[2], g[1]] = 1 for g in G ]
    M[diagind(M)] .= 0;

    #call Ripser
    flt = Rips(M; threshold=1, sparse=true)
    barcodes = ripserer(flt; dim_max = _max_dim)
    # interpret the outout. Barcodes is an array of tuples of dim = max_dim + 1.
    # All tuples where the second entry is Inf are the persistent ones.
    # Betti numbers just count the number of those persistent ones.
    betti = [ count(x -> x[2] .== Inf, barcodes[i]) for i in 1:length(barcodes)]

    padding = (max_dim+1) - length(betti);
    return vcat(betti, Int.(zeros(padding)))

    ## old function using Ripser library: https://github.com/mtsch/Ripser.jl#master
    # if length(H) == 0
    #     return zeros(max_dim + 1)
    # else
    #
    #     G = hypergraph_subdivide(H; expansion = false)
    #     # build distance matrix to give as input to Ripser.
    #     M = 2.0*ones(length(H), length(H))
    #     [ M[g[1], g[2]] = 1 for g in G ]
    #
    #     #call Ripser
    #     barcodes = ripser(M; dim_max = max_dim, threshold = 1.5)
    #     # interpret the outout. Barcodes is an array of tuples of dim = max_dim + 1.
    #     # All tuples where the second entry is Inf are the persistent ones.
    #     # Betti numbers just count the number of those persistent ones.
    #     betti = [ count(x -> x[2] .== Inf, barcodes[i]) for i in 1:length(barcodes)]
    #
    #     return betti
    # end
end

max_issues_repo_name: mtangulo/coexistenceholes

max_forks_repo_name: mtangulo/coexistenceholes


text:
using Cairo.tex2pango
using Base.Test

fsize = 1.618034
@test tex2pango("ƒ_{Nyquist} [\\mu K]",fsize) == "ƒNyquist [μK]"