text:
using Zygote, NNlib, Test, Random, LinearAlgebra
using Zygote: gradient
using NNlib: conv
import Random

function ngradient(f, xs::AbstractArray...)
  grads = zero.(xs)
  for (x, Δ) in zip(xs, grads), i in 1:length(x)
    δ = sqrt(eps())
    tmp = x[i]
    x[i] = tmp - δ/2
    y1 = f(xs...)
    x[i] = tmp + δ/2
    y2 = f(xs...)
    x[i] = tmp
    Δ[i] = (y2-y1)/δ
  end
  return grads
end

gradcheck(f, xs...) =
  all(isapprox.(ngradient(f, xs...),
                gradient(f, xs...), rtol = 1e-5, atol = 1e-5))

gradtest(f, xs::AbstractArray...) = gradcheck((xs...) -> sum(sin.(f(xs...))), xs...)
gradtest(f, dims...) = gradtest(f, rand.(Float64, dims)...)

Random.seed!(0)

@test gradient(//, 2, 3) === (1//3, -2//9)

@test gradtest((x, W, b) -> σ.(W*x .+ b), 5, (2,5), 2)
@test gradtest((x, W, b) -> σ.(W*x .+ b), (5,3), (2,5), 2)
@test gradtest((x, W, b) -> logσ.(W*x .+ b), 5, (2,5), 2)
@test gradtest((x, W, b) -> logσ.(W*x .+ b), (5,3), (2,5), 2)

@test gradtest((w, x) -> w'*x, randn(10, 2), randn(10))
@test gradtest((w, x) -> transpose(w)*x, randn(5,5), randn(5,5))

@test gradtest(x -> sum(x, dims = (2, 3)), (3,4,5))
@test gradtest(x -> prod(x, dims = (2, 3)), (3,4,5))
@test gradtest(x -> prod(x), (3,4,5))

@test gradtest(x -> softmax(x).*(1:3), 3)
@test gradtest(x -> softmax(x).*(1:3), (3,5))
@test gradtest(x -> logsoftmax(x).*(1:3), 3)
@test gradtest(x -> logsoftmax(x).*(1:3), (3,5))

@test gradtest(conv, rand(10, 3, 2), randn(Float64,2, 3, 2))
@test gradtest(conv, rand(10, 10, 3, 2), randn(Float64,2, 2, 3, 2))
@test gradtest(conv, rand(10, 10, 10, 3, 2), randn(Float64,2, 2, 2, 3, 2))

@test gradtest(x -> maxpool(x, (2,2)), rand(10, 10, 3, 2))
@test gradtest(x -> maxpool(x, (2,2,2)), rand(10, 10, 10, 3, 2))

@test gradtest(x -> meanpool(x, (2,2)), rand(10, 10, 3, 2))
@test gradtest(x -> meanpool(x, (2,2,2)), rand(5, 5, 5, 3, 2))

@test gradtest(x -> permutedims(x, [3,1,2]), rand(4,5,6))

@test gradtest(x -> repeat(x; inner=2), rand(5))
@test gradtest(x -> repeat(x; inner=2, outer=3), rand(5))
@test gradtest(x -> repeat(x; inner=(2,2,1), outer=(1,1,3)), rand(5,4,3))

@test gradtest(kron, rand(5), rand(3))
@test gradtest(kron, rand(5), rand(3), rand(8))
@test gradtest(kron, rand(5,1), rand(3,1))
@test gradtest(kron, rand(5,1), rand(3,1), rand(8,1))
@test gradtest(kron, rand(5,2), rand(3,2), rand(8,2))

@testset "maximum" begin
  @test gradtest(maximum, rand(2, 3))

  @test gradtest(x -> maximum(x, dims=1), rand(2, 3))
  @test gradtest(x -> maximum(x, dims=2), rand(2, 3))
  @test gradtest(x -> maximum(x, dims=3), rand(2, 3, 4))

  @test gradtest(x -> maximum(x, dims=[1, 2]), rand(2, 3, 4))
end

@testset "minimum" begin
  @test gradtest(minimum, rand(2, 3))

  @test gradtest(x -> minimum(x, dims=1), rand(2, 3))
  @test gradtest(x -> minimum(x, dims=2), rand(2, 3))
  @test gradtest(x -> minimum(x, dims=3), rand(2, 3, 4))

  @test gradtest(x -> minimum(x, dims=[1, 2]), rand(2, 3, 4))
end

@testset "backsolve" begin
  rng, P, Q = MersenneTwister(123456), 10, 9
  X, Y, y = randn(rng, P, P), randn(rng, P, Q), randn(rng, P)

  # \
  @test gradtest(X -> X \ Y, X)
  @test gradtest(Y -> X \ Y, Y)
  @test gradtest(X -> X \ y, X)
  @test gradtest(y -> X \ y, y)

  # /
  @test gradtest(X -> Y' / X, X)
  @test gradtest(Y -> Y' / X, Y)
  @test gradtest(X -> y' / X, X)
  @test gradtest(y -> y' / X, y)
end

@testset "Symmetric" begin
  rng, P = MersenneTwister(123456), 7
  A = randn(rng, P, P)
  @test gradtest(Symmetric, A)
end

@testset "diag" begin
  rng, P = MersenneTwister(123456), 10
  A = randn(rng, P, P)
  @test gradtest(diag, A)
end

@testset "dense + UniformScaling" begin
  rng = MersenneTwister(123456)
  A, λ = randn(rng, 10, 10), randn(rng)
  @test gradtest(A->A + 5I, A)
  @test gradtest(λ->A + λ[1] * I, [λ])
end

@testset "cholesky" begin
  rng, N = MersenneTwister(123456), 5
  A = randn(rng, N, N)
  @test gradtest(A->logdet(cholesky(A' * A + 1e-6I)), A)
  @testset "cholesky - scalar" begin
    y, back = Zygote.forward(cholesky, 5.0 * ones(1, 1))
    y′, back′ = Zygote.forward(cholesky, 5.0)
    C̄ = randn(rng, 1, 1)
    @test back′((factors=C̄,))[1] isa Real
    @test back′((factors=C̄,))[1] ≈ back((factors=C̄,))[1][1, 1]
  end
end

using Distances

Zygote.refresh()

@testset "distances" begin
  rng, P, Q, D = MersenneTwister(123456), 10, 9, 8

   # Check sqeuclidean.
  let
    x, y = randn(rng, D), randn(rng, D)
    @test gradtest(x->sqeuclidean(x, y), x)
    @test gradtest(y->sqeuclidean(x, y), y)
  end

   # Check binary colwise.
  let
    X, Y = randn(rng, D, P), randn(rng, D, P)
    @test gradtest(X->colwise(SqEuclidean(), X, Y), X)
    @test gradtest(Y->colwise(SqEuclidean(), X, Y), Y)
  end

   # Check binary pairwise.
  let
    X, Y = randn(rng, D, P), randn(rng, D, Q)
    @test gradtest(X->pairwise(SqEuclidean(), X, Y), X)
    @test gradtest(Y->pairwise(SqEuclidean(), X, Y), Y)
  end

   # Check unary pairwise.
  @test gradtest(X->pairwise(SqEuclidean(), X), randn(rng, D, P))
end

function cat_test(f, A::Union{AbstractVector, AbstractMatrix}...)
  @test gradtest(f, A...)
  Z, back = Zygote.forward(f, A...)
  Ā = back(randn(size(Z)))
  @test all(map((a, ā)->ā isa typeof(a), A, Ā))
end

@testset "vcat" begin

  # Vector-only.
  cat_test(vcat, randn(1))
  cat_test(vcat, randn(10))

  cat_test(vcat, randn(1), randn(1))
  cat_test(vcat, randn(5), randn(1))
  cat_test(vcat, randn(1), randn(6))
  cat_test(vcat, randn(5), randn(6))

  # Matrix-only.
  for c in [1, 2], r1 in [1, 2], r2 in [1, 2]
    cat_test(vcat, randn(r1, c))
    cat_test(vcat, randn(r1, c), randn(r2, c))
  end

  # Matrix-Vector / Vector-Matrix / Vector-Matrix-Vector.
  for r in [1, 2]
    cat_test(vcat, randn(r, 1), randn(r))
    cat_test(vcat, randn(r), randn(r, 1))
    cat_test(vcat, randn(r), randn(r, 1), randn(r))
  end
end

@testset "hcat" begin

  # Vector-only.
  for r in [1, 2]
    cat_test(hcat, randn(r))
    cat_test(hcat, randn(r), randn(r))
  end

  # Matrix-only.
  for r in [1, 2], c1 in [1, 2], c2 in [1, 2]
    cat_test(hcat, randn(r, c1), randn(r, c2))
  end

  # Matrix-Vector / Vector-Matrix / Vector-Matrix-Vector.
  for r in [1, 2], c in [1, 2]
    cat_test(hcat, randn(r, c), randn(r))
    cat_test(hcat, randn(r), randn(r, c))
    cat_test(hcat, randn(r), randn(r, c), randn(r))
  end
end

@testset "one(s) and zero(s)" begin
  @test Zygote.gradient(x->sum(ones(size(x))), randn(5))[1] isa Nothing
  @test Zygote.gradient(x->sum(one(x)), randn(3, 3))[1] isa Nothing
  @test Zygote.gradient(x->sum(zeros(size(x))), randn(7))[1] isa Nothing
  @test Zygote.gradient(x->sum(zero(x)), randn(3))[1] isa Nothing
end

@testset "* sizing" begin
  @test size(Zygote.gradient((x, y)->sum(x * y), randn(1, 1), randn(1, 10))[1]) == (1, 1)
  @test size(Zygote.gradient((x, y)->sum(x * y), randn(1, 1), randn(1, 10))[2]) == (1, 10)
end


text:
## Testing hand-coded trait implementations
###########################################

# impl will do:
#  - check supertraits are implemented
#  - import bindings
#  - make methods

# @traitimpl D2{T1, T2} begin
# +(x::T1, y::T2) = x.t1 + y.t2
# -(x::T1, y::T2) = x.t1 - y.t2
# end

# becomes:

# helper functions
function get_fnname(ex::Expr)
    out = ex.args[1].args[1]
    if typeof(out)==Symbol
        return out
    elseif out.head==:curly
        out = out.args[1]
    else
        error("Could not parse function definition: $ex")
    end
end
get_fnsig(ex::Expr) = ex.args[1].args[2:end]
function prefix_module!(ex::Expr, modname::Symbol)
    # prefix the module of the function:
    # :(+(x::T1, y::T2) = x.t1 + y.t2) -> :(Base.+(x::T1, y::T2) = x.t1 + y.t2)
    ex.head== :(=) || error("Not a function definition: $ex")
    ex.args[1].head==:call || error("Not a function definition: $ex")

    fnname = get_fnname(ex)
    if typeof(fnname)==Symbol
        ex.args[1].args[1] = :($modname.$fnname)
    elseif fnname.head==:(.)
        # has a module already prefixed
    else
        error("something went wrong, $fnname is not a Symbol")
    end
end

type T1
t1::Float64
end
type T2
t2::Float64
end

# functions in trait definition
implfns = [:(+(x::T1, y::T2) = x.t1 + y.t2),
           :(-(x::T1, y::T2) = x.t1 - y.t2)]
#fn3 = :(sin(x::T1) = sin(x.t1))

# check supertraits are satisfied
Base.sin(x::T1) = sin(x.t1)
Base.sin(x::T2) = sin(x.t2)
Base.cos(x::T1) = cos(x.t1)
Base.cos(x::T2) = cos(x.t2)
@assert istrait(traitgetsuper(D2{T1, T2}))

# make dict of implfns
tmp = implfns
implfns = Dict()
for implfn in tmp
    implfns[eval(get_fnname(implfn))] = implfn
end
# check that there are no double defs
length(tmp)==length(implfns) || error("Duplicate method definition(s)")

# check right number of defs
length(D2{T1,T2}().methods)==length(implfns) || error("Not right number of method definitions")
# check that the signature of fns agrees with D2{T1,T2}().methods
for (fn,_fn) in D2{T1,T2}().methods
    # for now just check length
    if length(_fn.env.defs.sig)-1!=length(get_fnsig(implfns[fn])) # -1 because of return type
        error("""Method definition:
                 $fn  $sig
                 does not match implementation:
                 $(implfns[fn])""")
    end
end

# make method definitions
for (fn, fndef) in implfns
    modname = module_name(Base.function_module(fn, Tuple{Vararg{Any}}))
    # I think meta-programming is needed to make the new function def
    prefix_module!(fndef, modname)
    eval(fndef)
end

# checks
@assert istrait(D2{T1, T2})
@assert T1(5) + T2(4)==5+4

ext: jl


text: # CUDA related types


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

using Markdown
using InteractiveUtils

# ╔═╡ ce74e976-4be2-11eb-29e9-4f2e622164a7
using HePPCAT

# ╔═╡ 21e07958-4be4-11eb-3f15-791d9e54ed12
using LinearAlgebra

# ╔═╡ 16f41744-4be3-11eb-1b34-4d680b206e76
md"""
# Quick start notebook for [`HePPCAT.jl`](https://github.com/dahong67/HePPCAT.jl)
"""

# ╔═╡ 261c7946-4be3-11eb-22a1-2fcb36cd9de6
md"""
**Step 1:** load the package.
"""

# ╔═╡ 13a9c3ae-4be3-11eb-0400-d94006c590e7
md"""
**Step 2:** form data into a `Vector` with `Matrix` elements -
each matrix is a group of samples (columns) having equal noise variance.
"""

# ╔═╡ 494732f8-4be3-11eb-22e3-9b1ac3913f41
begin # Example: generating synthetic data with heteroscedastic noise
   d = 100        # number of features (ambient dimension)
   k = 1          # number of factors (latent dimension)

   n = [800,200]  # number of samples in each group/block
   v = [10,0.01]  # noise variance for each group/block
   L = length(n)  # number of groups/blocks

   F = randn(d,k)/sqrt(d) # synthetic factor matrix
   X = [F*randn(k,n[l]) + sqrt(v[l])*randn(d,n[l]) for l in 1:L]  # synthetic data
end

# ╔═╡ 4cc04d34-4be3-11eb-31d5-0fb11143809b
md"""
This code generates an example synthetic dataset with heteroscedastic noise:
the `n[1]` = $(n[1]) samples in `X[1]` have noise variance `v[1]` = $(v[1]),
and the `n[2]` = $(n[2]) samples in `X[2]` have noise variance `v[2]` = $(v[2]).
"""

# ╔═╡ 707ac7cc-4be3-11eb-3091-4b5397319891
md"""
**Step 3:** run the `heppcat` method.
"""

# ╔═╡ a223fd5c-4be3-11eb-097b-4db63a3a9aa3
model = heppcat(X,k,100)  # run 100 iterations

# ╔═╡ ad0256d8-4be3-11eb-084b-2ffb4fece64d
md"""
For more info, use `?heppcat` to access the [docstring](https://docs.julialang.org/en/v1/manual/documentation/#Accessing-Documentation).
"""

# ╔═╡ c7db5bf8-4be3-11eb-3460-15b039b47f07
md"""
**Step 4:** extract factor matrix and noise variance estimates.
"""

# ╔═╡ e6a22602-4be3-11eb-0532-276b6a63a00b
model.F

# ╔═╡ e73c22c0-4be3-11eb-276a-afb37f0146ea
model.v

# ╔═╡ e9a19a04-4be3-11eb-0f4e-a72c20cf75e4
md"""
For more info about the returned datatype use `?HePPCATModel`
for the [docstring](https://docs.julialang.org/en/v1/manual/documentation/#Accessing-Documentation).
"""

# ╔═╡ 1fcfcefc-4be4-11eb-0041-cf96d0098e53
md"""
**Evaluating recovery:**
check recovery of the factor covariance and noise variances.
"""

# ╔═╡ 2841c4aa-4be4-11eb-206e-df2fd2696ac6
norm(model.F*model.F' - F*F')/norm(F*F')  # normalized estimation error

# ╔═╡ 2d2f42bc-4be4-11eb-043d-2723dad42b10
[model.v v]

# ╔═╡ 393e4b20-4be4-11eb-3d69-6fa84bae8d49
md"""
Compared with homoscedastic PPCA:
"""

# ╔═╡ 5b273d08-4be4-11eb-1d7b-7373113b2ca6
ppca = heppcat(X,k,0)  # initialization is homoscedastic PPCA

# ╔═╡ 5f50f2c0-4be4-11eb-1939-518d61b5fa73
norm(ppca.F*ppca.F' - F*F')/norm(F*F')  # normalized estimation error

# ╔═╡ 6281428a-4be4-11eb-0b2e-e7f9423afdb4
[ppca.v v]

# ╔═╡ 65b6c7c4-4be4-11eb-0950-b71e94e85221
md"""
`HePPCAT` accounts for heterogeneous quality among the samples
and is generally more robust.
"""

# ╔═╡ Cell order:
# ╟─16f41744-4be3-11eb-1b34-4d680b206e76
# ╟─261c7946-4be3-11eb-22a1-2fcb36cd9de6
# ╠═ce74e976-4be2-11eb-29e9-4f2e622164a7
# ╟─13a9c3ae-4be3-11eb-0400-d94006c590e7
# ╠═494732f8-4be3-11eb-22e3-9b1ac3913f41
# ╟─4cc04d34-4be3-11eb-31d5-0fb11143809b
# ╟─707ac7cc-4be3-11eb-3091-4b5397319891
# ╠═a223fd5c-4be3-11eb-097b-4db63a3a9aa3
# ╟─ad0256d8-4be3-11eb-084b-2ffb4fece64d
# ╟─c7db5bf8-4be3-11eb-3460-15b039b47f07
# ╠═e6a22602-4be3-11eb-0532-276b6a63a00b
# ╠═e73c22c0-4be3-11eb-276a-afb37f0146ea
# ╟─e9a19a04-4be3-11eb-0f4e-a72c20cf75e4
# ╟─1fcfcefc-4be4-11eb-0041-cf96d0098e53
# ╠═21e07958-4be4-11eb-3f15-791d9e54ed12
# ╠═2841c4aa-4be4-11eb-206e-df2fd2696ac6
# ╠═2d2f42bc-4be4-11eb-043d-2723dad42b10
# ╟─393e4b20-4be4-11eb-3d69-6fa84bae8d49
# ╠═5b273d08-4be4-11eb-1d7b-7373113b2ca6
# ╠═5f50f2c0-4be4-11eb-1939-518d61b5fa73
# ╠═6281428a-4be4-11eb-0b2e-e7f9423afdb4
# ╟─65b6c7c4-4be4-11eb-0950-b71e94e85221


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random, TimeSeries, MarketData, Dates
CairoMakie.activate!() #HIDE
let
    tempo = string.(timestamp(ohlc))
    lentime = length(tempo)
    slice_dates = range(1,lentime, step = lentime÷20)

    colors = values(ohlc.Close) .> values(ohlc.Open)
    lowV = values(ohlc.Low)
    highV = values(ohlc.High)
    linesegs = []
    for i in 1:lentime
        push!(linesegs, Point2f0(i, lowV[i]))
        push!(linesegs, Point2f0(i, highV[i]))
    end
    linesegs = Point2f0.(linesegs)
    cmap = ("#f64325", "#78f518")

    fig = Figure(resolution = (1200, 700), font = "sans", fontsize = 20)
    ax = Axis(fig, backgroundcolor = :black, ygridcolor = "#65866b",
        xgridcolor = "#65866b", xgridstyle=:dash, ygridstyle=:dash)
    barplot!(ax, 1:lentime, values(ohlc.Open), fillto = values(ohlc.Close),
        color = colors,strokewidth = 0.5, strokecolor = colors, colormap = cmap)
    linesegments!(ax, linesegs, color = colors, colormap = cmap)
    xlims!(ax, 200, 300)
    ylims!(ax, 13,24)
    ax.xticks = (slice_dates, tempo[slice_dates])
    ax.xticklabelrotation = π/4
    fig[1,1] = ax
    fig
    save(joinpath(@__DIR__, "output", "candlestick.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random","TimeSeries", "MarketData", "Dates"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Tomoko
using Random
using Test
using Distributions

p = PopRates()
x = bitrand(p.loci)
df = runsim(p,1:1:1000)
sample_betabin(α,β,n,N) = mapreduce(x-> [x,n-x],hcat,rand(BetaBinomial(n,α,β),N))
αtest = 0.05
dat = [sample_betabin(αtest,αtest,100,200) for ii in 1:20]

@testset "Tomoko.jl" begin
    0.0 < p.f(x)
    (mean(df.D) - Tomoko.pop_D(p))/Tomoko.pop_D(p) < 0.05 # check that D is near the equilibrium value
    abs((mean(df.pop_size) - p.κ)/p.κ) < 0.05 # check that the pop_size approaches the eq value
    abs((median(Tomoko.neutral_D.(dat)) - 2 αtest)/(2*αtest)) < 0.1 # check that minka algorithm is correct
end

max_stars_repo_name: StatPhysBio/Tomoko

max_issues_repo_name: StatPhysBio/Tomoko

max_forks_repo_name: StatPhysBio/Tomoko


text:
struct SuperOpConnection{A,B,C} <: AbsOpConnection
    op_conn_l_id::A
    op_conn_r_id::B

    op_conn_l_r::C
end

SuperOpConnection{A,B,C}() where {A,B,C} = SuperOpConnection(A(), B(), C())

@inline length_l_id(c::SuperOpConnection) = length(c.op_conn_l_id)
@inline length_r_id(c::SuperOpConnection) = length(c.op_conn_r_id)
@inline length_l_r(c::SuperOpConnection)  = length(c.op_conn_l_r)
@inline Base.length(c::SuperOpConnection) = length_l_id(c) + length_r_id(c) + length_l_r(c)
@inline Base.size(c::SuperOpConnection) = (length(c), )
function Base.eltype(::Type{SuperOpConnection{L,R,LR}}) where {L,R,LR}
    T1, SC1 = eltype(L).parameters
    T2, SC2 = eltype(R).parameters
    T3, SC3 = eltype(R).parameters
    T = promote_type(T1, T2, T3)

    return Tuple{T, SC1}
end

Base.:(==)(a::SuperOpConnection, b::SuperOpConnection) = (a.op_conn_l_id == b.to_change &&
                                                          a.op_conn_r_id == b.op_conn_r_id &&
                                                          a.op_conn_l_r == b.op_conn_l_r )
# ad
#Base.@propagate_inbounds Base.getindex(c::OpConnection, i) = (c.mel[i], c.changes[i])

function Base.iterate(iter::SuperOpConnection, state=(1,(1,1)))

    el, inner_state = state

    if el == 1
        res = iterate(iter.op_conn_l_id, inner_state)

        if !isnothing(res)
            conn, inner = res
            return conn, (el, inner)
        else
            el = 2
            inner_state = (1,1)
        end
    end

    if el == 2
        res = iterate(iter.op_conn_r_id, inner_state)

        if !isnothing(res)
            conn, inner = res
            return conn, (el, inner)
        else
            el = 3
            inner_state = (1,1)
        end
    end

    if el == 3
        res = iterate(iter.op_conn_l_r, inner_state)

        if !isnothing(res)
            conn, inner = res
            return conn, (el, inner)
        else
            el = 4
            inner_state = (1,1)
        end
    end

    return nothing
end

@inline Base.getindex(c::SuperOpConnection, i) = begin
    if i <= length_l_id(c)
        return getindex(c.op_conn_l_id, i)
    elseif i <= length_l_id(c) + length_r_id(c)
        return getindex(c.op_conn_r_id, i - length_l_id(c))
    else
        return getindex(c.op_conn_l_r, i - length_l_id(c) - length_r_id(c))
    end
end

# showing
function Base.show(io::IO, c::SuperOpConnection)
    print(io, "$(length(c))-elements - SuperOpConnection{...}:\n")
    if length(c) > 0
        Base.print_matrix(IOContext(io, :compact=>true), collect(c))
    end
    return io
end


text:
include("makeMaze.jl")

@time g = makeMaze(200, 200, [1, 1], [200, 200], 1000)

writecsv("E:\\GitHub\\makeMaze\\maze_out.csv", getMap(g))
writecsv("E:\\GitHub\\makeMaze\\maze_patha_out.csv", getPath(g))

p = getPath(g);

g.path_b

g.path_b[10]

g.pa
g.pb

quit()

max_stars_repo_name: RandomMaze/makeMaze

max_issues_repo_name: RandomMaze/makeMaze

max_forks_repo_name: RandomMaze/makeMaze


text:
# Tseng, "On Accelerated Proximal Gradient Methods for Convex-Concave
# Optimization" (2008).
#
# Beck, Teboulle, "A Fast Iterative Shrinkage-Thresholding Algorithm
# for Linear Inverse Problems", SIAM Journal on Imaging Sciences, vol. 2,
# no. 1, pp. 183-202 (2009).

using Base.Iterators
using ProximalAlgorithms.IterationTools
using ProximalOperators: Zero
using LinearAlgebra
using Printf

"""
    ForwardBackwardIteration(; <keyword-arguments>)

Instantiate the forward-backward splitting algorithm (see [1, 2]) for solving
optimization problems of the form

    minimize f(Ax) + g(x),

where `f` is smooth and `A` is a linear mapping (for example, a matrix).

# Arguments
- `x0`: initial point.
- `f=Zero()`: smooth objective term.
- `A=I`: linear operator (e.g. a matrix).
- `g=Zero()`: proximable objective term.
- `Lf=nothing`: Lipschitz constant of the gradient of x ↦ f(Ax).
- `gamma=nothing`: stepsize to use, defaults to `1/Lf` if not set (but `Lf` is).
- `adaptive=false`: forces the method stepsize to be adaptively adjusted.
- `minimum_gamma=1e-7`: lower bound to `gamma` in case `adaptive == true`.
- `fast=false`: enables Nesterov acceleration.

# References
- [1] Tseng, "On Accelerated Proximal Gradient Methods for Convex-Concave
Optimization" (2008).
- [2] Beck, Teboulle, "A Fast Iterative Shrinkage-Thresholding Algorithm
for Linear Inverse Problems", SIAM Journal on Imaging Sciences, vol. 2, no. 1,
pp. 183-202 (2009).
"""

@Base.kwdef struct ForwardBackwardIteration{R,C<:Union{R,Complex{R}},Tx<:AbstractArray{C},Tf,TA,Tg}
    f::Tf = Zero()
    A::TA = I
    g::Tg = Zero()
    x0::Tx
    Lf::Maybe{R} = nothing
    gamma::Maybe{R} = Lf === nothing ? nothing : (1 / Lf)
    adaptive::Bool = false
    minimum_gamma::R = real(eltype(x0))(1e-7)
    fast::Bool = false
end

Base.IteratorSize(::Type{<:ForwardBackwardIteration}) = Base.IsInfinite()

Base.@kwdef mutable struct ForwardBackwardState{R,Tx,TAx}
    x::Tx             # iterate
    Ax::TAx           # A times x
    f_Ax::R           # value of smooth term
    grad_f_Ax::TAx    # gradient of f at Ax
    At_grad_f_Ax::Tx  # gradient of smooth term
    gamma::R          # stepsize parameter of forward and backward steps
    y::Tx             # forward point
    z::Tx             # forward-backward point
    g_z::R            # value of nonsmooth term (at z)
    res::Tx           # fixed-point residual at iterate (= z - x)
    theta::R = one(real(eltype(x)))
    z_prev::Tx = copy(x)
end

f_model(state::ForwardBackwardState) = f_model(state.f_Ax, state.At_grad_f_Ax, state.res, state.gamma)

function Base.iterate(iter::ForwardBackwardIteration{R}) where {R}
    x = copy(iter.x0)
    Ax = iter.A * x
    grad_f_Ax, f_Ax = gradient(iter.f, Ax)

    gamma = iter.gamma

    if gamma === nothing
        # compute lower bound to Lipschitz constant of the gradient of x ↦ f(Ax)
        xeps = x .+ R(1)
        grad_f_Axeps, f_Axeps = gradient(iter.f, iter.A * xeps)
        L = norm(iter.A' * (grad_f_Axeps - grad_f_Ax)) / R(sqrt(length(x)))
        gamma = R(1) / L
    end

    # compute initial forward-backward step
    At_grad_f_Ax = iter.A' * grad_f_Ax
    y = x - gamma .* At_grad_f_Ax
    z, g_z = prox(iter.g, y, gamma)

    # compute initial fixed-point residual
    res = x - z

    state = ForwardBackwardState(
        x=x, Ax=Ax, f_Ax=f_Ax, grad_f_Ax=grad_f_Ax, At_grad_f_Ax=At_grad_f_Ax,
        gamma=gamma, y=y, z=z, g_z=g_z, res=res,
    )

    return state, state
end

function Base.iterate(iter::ForwardBackwardIteration{R}, state::ForwardBackwardState{R,Tx,TAx}) where {R,Tx,TAx}
    Az, f_Az, grad_f_Az, At_grad_f_Az = nothing, nothing, nothing, nothing
    a, b, c = nothing, nothing, nothing

    # backtrack gamma (warn and halt if gamma gets too small)
    while iter.gamma === nothing || iter.adaptive == true
        if state.gamma < iter.minimum_gamma
            @warn "parameter `gamma` became too small ($(state.gamma)), stopping the iterations"
            return nothing
        end
        f_Az_upp = f_model(state)
        Az = iter.A * state.z
        grad_f_Az, f_Az = gradient(iter.f, Az)
        tol = 10 * eps(R) * (1 + abs(f_Az))
        if f_Az <= f_Az_upp + tol
            break
        end
        state.gamma *= 0.5
        state.y .= state.x .- state.gamma .* state.At_grad_f_Ax
        state.g_z = prox!(state.z, iter.g, state.y, state.gamma)
        state.res .= state.x .- state.z
    end

    if iter.fast == true
        theta1 = (R(1) + sqrt(R(1) + 4 * state.theta^2)) / R(2)
        extr = (state.theta - R(1)) / theta1
        state.theta = theta1
        state.x .= state.z .+ extr .* (state.z .- state.z_prev)
        state.z_prev, state.z = state.z, state.z_prev
    else
        state.x, state.z = state.z, state.x
    end

    # TODO: if iter.fast == true, in the adaptive case we should still be able
    # to save some computation by extrapolating Ax and (if f is quadratic)
    # f_Ax, grad_f_Ax, At_grad_f_Ax.
    if iter.fast == false && (iter.gamma === nothing || iter.adaptive == true)
        state.Ax = Az
        state.f_Ax = f_Az
        state.grad_f_Ax = grad_f_Az
    else
        mul!(state.Ax, iter.A, state.x)
        state.f_Ax = gradient!(state.grad_f_Ax, iter.f, state.Ax)
    end

    mul!(state.At_grad_f_Ax, adjoint(iter.A), state.grad_f_Ax)
    state.y .= state.x .- state.gamma .* state.At_grad_f_Ax
    state.g_z = prox!(state.z, iter.g, state.y, state.gamma)

    state.res .= state.x .- state.z

    return state, state
end

# Solver

struct ForwardBackward{R, K}
    maxit::Int
    tol::R
    verbose::Bool
    freq::Int
    kwargs::K
end

function (solver::ForwardBackward)(x0; kwargs...)
    stop(state::ForwardBackwardState) = norm(state.res, Inf) / state.gamma <= solver.tol
    disp((it, state)) =
        @printf("%5d | %.3e | %.3e\n", it, state.gamma, norm(state.res, Inf) / state.gamma)
    iter = ForwardBackwardIteration(; x0=x0, solver.kwargs..., kwargs...)
    iter = take(halt(iter, stop), solver.maxit)
    iter = enumerate(iter)
    if solver.verbose
        iter = tee(sample(iter, solver.freq), disp)
    end
    num_iters, state_final = loop(iter)
    return state_final.z, num_iters
end

ForwardBackward(; maxit=10_000, tol=1e-8, verbose=false, freq=100, kwargs...) =
    ForwardBackward(maxit, tol, verbose, freq, kwargs)


text:
@testset "ResNet" begin
  for model in [ResNet18, ResNet34, ResNet50, ResNet101, ResNet152]
    m = model()
    @test size(m(rand(Float32, 256, 256, 3, 50))) == (1000, 50)
  end
end


text:
include("header.jl")

#=
# Basic Tools for Quasi-1D Steady Compressible Flow
This notebook demonstrates the basic syntax for some tools for computing
quasi-1d steady compressible flow.
=#

# ### Set up the module
using MAE103

#=
### Setting basic properties and states
We can set thermodynamic properties and states in a straightforward manner.
However, it is important to remember that we have to explicitly define the
type of property or state we are setting. Examples below will show how this works.
=#

#=
Let's say we wish to set the pressure to 10000 Pa. Pascals are the default units
of pressure, as can be verified by using the `default_unit` function:
=#
default_unit(Pressure)

#=
So if we do not specify the unit, it is **automatically set to the default unit**:
=#
Pressure(10000)

#=
We can set a quantity with another unit using the syntax u"[unit]". For example,
if we set pressure to 1 atm, it will still convert it to the default
unit.
=#
p = Pressure(1u"atm")

#=
However, we can always report the quantity in some desired units with the `value`
function:
=#
value(p,u"atm")
#-
value(p,u"psi")
#-
value(p,u"kPa")

#=
#### Other thermodynamic quantities
We can set most any other thermodynamic quantity in similar fashion:
=#
T = Temperature(20u"°C")
#-
T0 = StagnationTemperature(20)
#-
MachNumber(2.0)
#-
Enthalpy(50)
#-
Entropy(10)
#-
Area(50u"cm^2")
#-
Length(5)
#=
and others...
=#

#=
#### Gas properties
We can set the properties of the gas that we are analyzing. (Note: It is
assumed that the gas is perfect.)
=#
SpecificHeatRatio(1.3)
#-
GasConstant(320)

#=
and we can define a gas with these values:
=#
gas = PerfectGas(γ=SpecificHeatRatio(1.3),R=GasConstant(320))

#=
We have **pre-defined gases** (at standard conditions), as well, for convenience:
=#
Air
#-
He
#-
O2
#-
CO2
#-
H2
#-
N2
#-
Ar
#=
#### Equations of state
We can apply the equation of state for a perfect gas to determine other quantities.
For example, suppose we have carbon dioxide at 1.2 kg/m^3 and 80 kPa. What is the temperature?
=#
T = Temperature(Density(1.2),Pressure(80u"kPa"),gas=CO2)

#=
You can switch the order of the arguments and it will still work:
=#
T = Temperature(Pressure(80u"kPa"),Density(1.2),gas=CO2)

#=
Then we can calculate the enthalpy, for example:
=#
Enthalpy(T,gas=CO2)

#=
What is the speed of sound of air at 20 degrees Celsius? Let's find out:
=#
SoundSpeed(Temperature(20u"°C"),gas=Air)

#=
How about oxygen?
=#
SoundSpeed(Temperature(20u"°C"),gas=O2)

#=
**Note: the default gas is air. So if you do not put the `gas=` argument in,
it will assume air at standard conditions.**
=#


text:
import Base.show
show(io::IO,eq::Equation)=print(io,ps(eq))
show(io::IO,ex::Factor)=print(io,ps(ex))

ps(eq::Equation)="$(ps(eq.lhs)) ≖ $(ps(eq.rhs))"
function ps(ex::Expression)
	s=""
	for t in ex
		s*=ps(t)*"+"
	end
	s[1:end-1]
end
function ps(t::Term)
	s=""
	for f in t
		s*=ps(f)*"*"
	end
	s[1:end-1]
end
function ps{T}(ar::Array{T,1})
	if isempty(ar)
		return "$T[]"
	end
	s="$T["
	for a in ar
		s*=ps(a)*","
	end
	s[1:end-1]*"]"
end
#function ps{T}(a::T)
#	"$T"
#end
function ps(a)
	io=IOBuffer()
	show(io,a)
	str=takebuf_string(io)
end
function ps(c::Component)
	ar=getargs(c)
	s="$(typeof(c))("
	for a in ar
		s*="$(ps(a)),"
	end
	s[1:end-1]*")"
end
ps(s::Symbol)=":$s"
ps(n::Number)="$n"

function print(io::IO,u::U)
	print(io,u.x)
	print(io," [")
	print(io,u.units)
	print(io,"]")
end
ps(a::Alt)="Alt($(a.indices))"
function show(io::IO,a::Alt)
	print(io,ps(a))
end
function print(io::IO,a::Alt)
	print(io,"ϵ$(a.indices)")
end


text:
# wrappers of the low-level CUBLAS functionality

import CUDAdrv: CuStream

# Note: we don't implement padded storage dimensions
function _mkplan(xtype, xdims, region)
    nrank = length(region)
    sz = [xdims[i] for i in region]
    csz = copy(sz)
    csz[1] = div(sz[1],2) + 1
    batch = prod(xdims) ÷ prod(sz)

    pp = Ref{cufftHandle_t}()
    if (nrank == 1) && (batch == 1)
        cufftPlan1d(pp, sz[1], xtype, 1)
    elseif (nrank == 2) && (batch == 1)
        cufftPlan2d(pp, sz[2], sz[1], xtype)
    elseif (nrank == 3) && (batch == 1)
        cufftPlan3d(pp, sz[3], sz[2], sz[1], xtype)
    else
        rsz = (length(sz) > 1) ? rsz = reverse(sz) : sz
        if ((region...,) == ((1:nrank)...,))
            # handle simple case ... simply! (for robustness)
           cufftPlanMany(pp, nrank, Cint[rsz...], C_NULL, 1, 1, C_NULL, 1, 1,
                         xtype, batch)
        else
            if nrank==1 || all(diff(collect(region)) .== 1)
                # _stride: successive elements in innermost dimension
                # _dist: distance between first elements of batches
                if region[1] == 1
                    istride = 1
                    idist = prod(sz)
                    cdist = prod(csz)
                else
                    if region[end] != length(xdims)
                        throw(ArgumentError("batching dims must be sequential"))
                    end
                    istride = prod(xdims[1:region[1]-1])
                    idist = 1
                    cdist = 1
                end
                inembed = Cint[rsz...]
                cnembed = (length(csz) > 1) ? Cint[reverse(csz)...] : Cint[csz[1]]
                ostride = istride
                if xtype == CUFFT_R2C || xtype == CUFFT_D2Z
                    odist = cdist
                    onembed = cnembed
                else
                    odist = idist
                    onembed = inembed
                end
                if xtype == CUFFT_C2R || xtype == CUFFT_Z2D
                    idist = cdist
                    inembed = cnembed
                end
            else
                if any(diff(collect(region)) .< 1)
                    throw(ArgumentError("region must be an increasing sequence"))
                end
                cdims = collect(xdims)
                cdims[region[1]] = div(cdims[region[1]],2)+1

                if region[1] == 1
                    istride = 1
                    ii=1
                    while (ii < nrank) && (region[ii] == region[ii+1]-1)
                        ii += 1
                    end
                    idist = prod(xdims[1:ii])
                    cdist = prod(cdims[1:ii])
                    ngaps = 0
                else
                    istride = prod(xdims[1:region[1]-1])
                    idist = 1
                    cdist = 1
                    ngaps = 1
                end
                nem = ones(Int,nrank)
                cem = ones(Int,nrank)
                id = 1
                for ii=1:nrank-1
                    if region[ii+1] > region[ii]+1
                        ngaps += 1
                    end
                    while id < region[ii+1]
                        nem[ii] *= xdims[id]
                        cem[ii] *= cdims[id]
                        id += 1
                    end
                    @assert nem[ii] >= sz[ii]
                end
                if region[end] < length(xdims)
                    ngaps += 1
                end
                # CUFFT represents batches by a single stride (_dist)
                # so we must verify that region is consistent with this:
                if ngaps > 1
                    throw(ArgumentError("batch regions must be sequential"))
                end

                inembed = Cint[reverse(nem)...]
                cnembed = Cint[reverse(cem)...]
                ostride = istride
                if xtype == CUFFT_R2C || xtype == CUFFT_D2Z
                    odist = cdist
                    onembed = cnembed
                else
                    odist = idist
                    onembed = inembed
                end
                if xtype == CUFFT_C2R || xtype == CUFFT_Z2D
                    idist = cdist
                    inembed = cnembed
                end
            end
            cufftPlanMany(pp, nrank, Cint[rsz...],
                          inembed, istride, idist, onembed, ostride, odist,
                          xtype, batch)
        end
    end
    pp[]
end

# this is used implicitly in the unsafe_execute methods below:
unsafe_convert(::Type{cufftHandle_t}, p::CuFFTPlan) = p.plan

convert(::Type{cufftHandle_t}, p::CuFFTPlan) = p.plan

destroy_plan(plan::CuFFTPlan) = cufftDestroy(plan)

set_stream(plan::CuFFTPlan, stream::CuStream) = cufftSetStream(plan, stream)

function assert_applicable(p::CuFFTPlan{T,K}, X::CuArray{T}) where {T,K}
    (size(X) == p.sz) ||
        throw(ArgumentError("CuFFT plan applied to wrong-size input"))
end

function assert_applicable(p::CuFFTPlan{T,K}, X::CuArray{T}, Y::CuArray{Ty}) where {T,K,Ty}
    assert_applicable(p, X)
    (size(Y) == p.osz) ||
        throw(ArgumentError("CuFFT plan applied to wrong-size output"))
    # type errors should be impossible by dispatch, but just in case:
    if p.xtype ∈ [CUFFT_C2R, CUFFT_Z2D]
        (Ty == real(T)) ||
            throw(ArgumentError("Type mismatch for argument Y"))
    elseif p.xtype ∈ [CUFFT_R2C, CUFFT_D2Z]
        (Ty == complex(T)) ||
            throw(ArgumentError("Type mismatch for argument Y"))
    else
        (Ty == T) ||
            throw(ArgumentError("Type mismatch for argument Y"))
    end
end

function unsafe_execute!(plan::cCuFFTPlan{cufftComplex,K,true,N},
                         x::CuArray{cufftComplex,N}) where {K,N}
    @assert plan.xtype == CUFFT_C2C
    cufftExecC2C(plan, x, x, K)
end
function unsafe_execute!(plan::rCuFFTPlan{cufftComplex,K,true,N},
                         x::CuArray{cufftComplex,N}) where {K,N}
    @assert plan.xtype == CUFFT_C2R
    cufftExecC2R(plan, x, x)
end

function unsafe_execute!(plan::cCuFFTPlan{cufftComplex,K,false,N},
                         x::CuArray{cufftComplex,N}, y::CuArray{cufftComplex}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_C2C
    cufftExecC2C(plan, x, y, K)
end
function unsafe_execute!(plan::rCuFFTPlan{cufftComplex,K,false,N},
                         x::CuArray{cufftComplex,N}, y::CuArray{cufftReal}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_C2R
    cufftExecC2R(plan, x, y)
end

function unsafe_execute!(plan::rCuFFTPlan{cufftReal,K,false,N},
                         x::CuArray{cufftReal,N}, y::CuArray{cufftComplex,N}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_R2C
    cufftExecR2C(plan, x, y)
end

# double prec.
function unsafe_execute!(plan::cCuFFTPlan{cufftDoubleComplex,K,true,N},
                         x::CuArray{cufftDoubleComplex,N}) where {K,N}
    @assert plan.xtype == CUFFT_Z2Z
    cufftExecZ2Z(plan, x, x, K)
end
function unsafe_execute!(plan::rCuFFTPlan{cufftDoubleComplex,K,true,N},
                         x::CuArray{cufftDoubleComplex,N}) where {K,N}
    @assert plan.xtype == CUFFT_Z2D
    cufftExecZ2D(plan, x, x)
end

function unsafe_execute!(plan::cCuFFTPlan{cufftDoubleComplex,K,false,N},
                         x::CuArray{cufftDoubleComplex,N}, y::CuArray{cufftDoubleComplex}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_Z2Z
    cufftExecZ2Z(plan, x, y, K)
end
function unsafe_execute!(plan::rCuFFTPlan{cufftDoubleComplex,K,false,N},
                         x::CuArray{cufftDoubleComplex,N}, y::CuArray{cufftDoubleReal}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_Z2D
    cufftExecZ2D(plan, x, y)
end

function unsafe_execute!(plan::rCuFFTPlan{cufftDoubleReal,K,false,N},
                         x::CuArray{cufftDoubleReal,N}, y::CuArray{cufftDoubleComplex,N}
                         ) where {K,N}
    @assert plan.xtype == CUFFT_D2Z
    cufftExecD2Z(plan, x, y)
end


text:
module VTIModuli

using NewtonMethods: quasinewton
using LinearAlgebra: I, Diagonal, Symmetric
using Roots: fzero

export vphaseP, vphaseSV, vphaseSH, moduli2thomsen, velocities0, thomsen2moduli, vgroup, groupangle, dvPdthetav, dvSVdthetav, vgroupSH, groupangleSH, phaseangleP, phaseangleSV, phaseangleSH, vpvsYoungsPoisson, stiffmatrixTI, cracks2moduli
export findthomsen, findmoduli, findcracks
export VMeasure

include("base.jl")
include("types.jl")
include("invert.jl")
include("cracks.jl")

end


text:
module LssTests

using Revise
using Subsystems
using BlockArrays
using Test
using Infiltrator

@testset "LSS" begin
    @testset "has correct" begin
        N = 5;
        ne = rand(1:5, N)
        me = rand(0:2, N)
        pe = rand(1:5, N)

        As = [rand(ne[i],ne[i]) for i in 1:N]
        Bs = [rand(ne[i],me[i]) for i in 1:N]
        Cs = [rand(pe[i],ne[i]) for i in 1:N]
        K = [[rand(ne[i],ne[j]) for j in 1:N] for i in 1:N]
        Inc = rand(Bool, N, N);

        subs = [LdtiSubsystem(i, As[i], Bs[i], Cs[i]) for i in 1:N]

        Abe = BlockArray{Real}(zeros(sum(ne), sum(ne)), ne, ne)
        Bbe = BlockArray{Real}(zeros(sum(ne), sum(me)), ne, me)
        Cbe = BlockArray{Real}(zeros(sum(pe), sum(ne)), pe, ne)

        for i in 1:N
            Abe[Block(i,i)] = As[i]
            Bbe[Block(i,i)] = Bs[i]
            Cbe[Block(i,i)] = Cs[i]

            for j in vcat(1:i-1, i+1:N)
                if Inc[i,j]
                    addNeighbour(subs[i], subs[j], K[i][j], false)
                    view(Abe, Block(i,j)) .= K[i][j]
                end
            end
        end

        @testset "dimensions" begin
            (ng, mg, pg) = _getdimarrays(subs)
            @test ng == ne
            @test mg == me
            @test pg == pe
        end

        @testset "Block matrices" begin
            ms = tomonolith(subs)
            @test Abe == ms.A
            @test Bbe == ms.B
            @test Cbe == ms.C
        end
    end
end

end


text:
#! /usr/bin/julia

# Rosetta Code, Brownian tree

using Color, Images, FixedPointNumbers

const W = 512
const W0 = W>>1
const H = 512
const H0 = H>>1
const N = iceil(1.0*W*H)
const SIDESTICK = false

function motecolor(x::Int, y::Int)
    h = clamp(180*(atan2(y-H0, x-W0)/pi + 1.0), 0.0, 360.0)
    return HSV(h, 0.5, 0.5)
end

img = Image(zeros(RGB{Ufixed8}, H, W))
img["x", W0, "y", H0] = RGB(1, 1, 1)
isfree = trues(W, H)
isfree[W0, H0] = false
for i in 1:N
    x = rand(1:W)
    y = rand(1:H)
    isfree[x, y] || continue
    mc = motecolor(x, y)
    for j in 1:10^3
        xp = x + rand(-1:1)
        yp = y + rand(-1:1)
        iscontained = 0 < xp <= W && 0 < yp <= H
        if  iscontained && isfree[xp, yp]
            x = xp
            y = yp
            continue
        else
            if SIDESTICK || (iscontained && !isfree[xp, yp])
                isfree[x, y] = false
                img["x", x, "y", y] = mc
            end
            break
        end
    end
end

imwrite(img, "brownian_tree.png")

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
struct AggregatedResults <: Results
    variables::Dict{Symbol, DataFrames.DataFrame}
    total_cost::Dict
    optimizer_log::Dict
    time_stamp::DataFrames.DataFrame
    check_sum::Array
    duals::Dict{Symbol, Any}
end

get_res_variables(result::AggregatedResults) = result.variables
get_cost(result::AggregatedResults) = result.total_cost
get_log(result::AggregatedResults) = result.optimizer_log
get_time_stamp(result::AggregatedResults) = result.time_stamp
get_duals(result::AggregatedResults) = result.duals
get_variables(result::AggregatedResults) = result.dual_variables
"""This function creates the correct results struct for the context"""
function _make_results(variables::Dict,
                      total_cost::Dict,
                      optimizer_log::Dict,
                      time_stamp::DataFrames.DataFrame,
                      check_sum::Array,
                      duals::Dict)
    return AggregatedResults(variables, total_cost, optimizer_log, time_stamp, check_sum, duals)
end

# internal function to parse through the reference dictionary and grab the file paths
function _reading_references(results::Dict, duals::Array, stage::String, step::Array,
                             references::Dict, extra_time::Int64)

    for name in (dual)
        date_df = references[stage][name]
        step_df = DataFrames.DataFrame(Date = Dates.DateTime[], Step = String[], File_Path = String[])
        for n in 1:length(step)
            step_df = vcat(step_df, date_df[date_df.Step .== step[n], :])
        end
        results[name] = DataFrames.DataFrame()
        for (ix,time) in enumerate(step_df.Date)
            file_path = step_df[ix, :File_Path]
            var = Feather.read("$file_path")
            correct_var_length = size(1:(size(var, 1) - extra_time), 1)
            results[name] = vcat(results[name],var[1:correct_var_length,:])
        end
    end
    return results
end
# internal function to parse through the reference dictionary and grab the file paths
function _reading_references(results::Dict, dual::Array, stage::String,
                             references::Dict, extra_time::Int64)
    for name in dual
        date_df = references[stage][name]
        results[name] = DataFrames.DataFrame()
        for (ix,time) in enumerate(date_df.Date)
            file_path = date_df[ix, :File_Path]
            var = Feather.read(file_path)
            correct_var_length = size(1:(size(var, 1) - extra_time_length), 1)
            results[name] = vcat(results[name],var[1:correct_var_length,:])
        end
    end
    return results
end
# internal function to remove the overlapping results and only use the most recent
function _reading_time(file_path::String, time_length::Number)
    time_file_path = joinpath(dirname(file_path), "time_stamp.feather")
    temp_time_stamp = Feather.read("$time_file_path")
    time_stamp = temp_time_stamp[(1:time_length),:]
    return time_stamp
end

""" This sums all of the rows in a result dataframe """
function rowsum(variable::DataFrames.DataFrame, name::String)
    variable = DataFrames.DataFrame(Symbol(name) => sum.(eachcol(variable)))
    return variable
end
""" This sums each column in a result dataframe """
function columnsum(variable::DataFrames.DataFrame)
    shortvar = DataFrames.DataFrame()
    varnames = collect(names(variable))
    eachsum = (sum.(eachrow(variable)))
    for i in 1: size(variable, 1)
        df = DataFrames.DataFrame(Symbol(varnames[i]) => eachsum[i])
        shortvar = hcat(shortvar, df)
    end
    return shortvar
end
# internal function to check for duals
function _find_duals(variables::Array)
    duals = []
    for i in 1:length(variables)
        if occursin("dual", String.(variables[i]))
            duals = vcat(duals, variables[i])
        end
    end
    return duals
end
# internal function for differentiating variables from duals in file names
function _concat_string(duals::Vector{Symbol})
    duals = (String.(duals)).*"_dual"
    return duals
end


text:
function _pass_abstract_jump(optimizer::Union{Nothing, JuMP.OptimizerFactory},
                              parameters::Bool,
                              JuMPmodel::Union{JuMP.AbstractModel,Nothing})

    if isa(optimizer, Nothing)
        @info("The optimization model has no optimizer attached")
    end

    if !isnothing(JuMPmodel) && parameters

        if !haskey(JuMPmodel.ext, :params)
            @info("Model doesn't have Parameters enabled. Parameters will be enabled")
        end

        PJ.enable_parameters(JuMPmodel)

        return JuMPmodel

    end

    JuMPmodel = JuMP.Model(optimizer)

    if parameters
        PJ.enable_parameters(JuMPmodel)
    end

    return JuMPmodel

end

function _make_container_array(V::DataType, parameters::Bool, ax...)

    if parameters
        return JuMP.Containers.DenseAxisArray{PGAE{V}}(undef, ax...)
    else
        return JuMP.Containers.DenseAxisArray{GAE{V}}(undef, ax...)
    end

    return

end

function _make_expressions_dict(transmission::Type{S},
                                V::DataType,
                                bus_numbers::Vector{Int64},
                                time_steps::UnitRange{Int64},
                                parameters::Bool) where {S<:PM.AbstractPowerModel}

    return DSDA(:nodal_balance_active =>  _make_container_array(V,
                                                                parameters,
                                                                bus_numbers,
                                                                time_steps),
                :nodal_balance_reactive => _make_container_array(V,
                                                                 parameters,
                                                                 bus_numbers,
                                                                 time_steps))

end

function _make_expressions_dict(transmission::Type{S},
                                V::DataType,
                                bus_numbers::Vector{Int64},
                                time_steps::UnitRange{Int64},
                                parameters::Bool) where {S<:PM.AbstractActivePowerModel}

    return DSDA(:nodal_balance_active =>  _make_container_array(V,
                                                                parameters,
                                                                bus_numbers,
                                                                time_steps))
end

function _canonical_init(bus_numbers::Vector{Int64},
                        jump_model::JuMP.AbstractModel,
                        optimizer::Union{Nothing,JuMP.OptimizerFactory},
                        transmission::Type{S},
                        time_steps::UnitRange{Int64},
                        resolution::Dates.Period,
                        initial_time::Dates.DateTime,
                        parameters::Bool,
                        sequential_runs::Bool,
                        ini_con::DICKDA) where {S<:PM.AbstractPowerModel}

    V = JuMP.variable_type(jump_model)

    canonical = CanonicalModel(jump_model,
                              optimizer,
                              parameters,
                              sequential_runs,
                              time_steps,
                              resolution,
                              initial_time,
                              DSDA(),
                              DSDA(),
                              zero(JuMP.GenericAffExpr{Float64, V}),
                              _make_expressions_dict(transmission,
                                                     V,
                                                     bus_numbers,
                                                     time_steps,
                                                     parameters),
                              parameters ? DRDA() : nothing,
                              ini_con,
                              nothing);

    return canonical

end

mutable struct CanonicalModel
    JuMPmodel::JuMP.AbstractModel
    optimizer_factory::Union{Nothing, JuMP.OptimizerFactory}
    parametrized::Bool
    sequential_runs::Bool
    time_steps::UnitRange{Int64}
    resolution::Dates.Period
    initial_time::Dates.DateTime
    variables::Dict{Symbol, JuMP.Containers.DenseAxisArray}
    constraints::Dict{Symbol, JuMP.Containers.DenseAxisArray}
    cost_function::JuMP.AbstractJuMPScalar
    expressions::Dict{Symbol, JuMP.Containers.DenseAxisArray}
    parameters::Union{Nothing, Dict{UpdateRef, JuMP.Containers.DenseAxisArray}}
    initial_conditions::DICKDA
    pm_model::Union{Nothing, PM.AbstractPowerModel}

    function CanonicalModel(JuMPmodel::JuMP.AbstractModel,
                            optimizer_factory::Union{Nothing, JuMP.OptimizerFactory},
                            parametrized::Bool,
                            sequential_runs::Bool,
                            time_steps::UnitRange{Int64},
                            resolution::Dates.Period,
                            initial_time::Dates.DateTime,
                            variables::Dict{Symbol, JuMP.Containers.DenseAxisArray},
                            constraints::Dict{Symbol, JuMP.Containers.DenseAxisArray},
                            cost_function::JuMP.AbstractJuMPScalar,
                            expressions::Dict{Symbol, JuMP.Containers.DenseAxisArray},
                            parameters::Union{Nothing, Dict{UpdateRef, JuMP.Containers.DenseAxisArray}},
                            initial_conditions::DICKDA,
                            pm_model::Union{Nothing, PM.AbstractPowerModel})

        #prevents having empty parameters and parametrized canonical model
        @assert isnothing(parameters) == !parametrized

        if (sequential_runs && sequential_runs != parametrized)
            throw(ArgumentError("Sequential simulations can't run when the specified OperationModel
                                 is not parametrized"))
        end

        new(JuMPmodel,
            optimizer_factory,
            parametrized,
            sequential_runs,
            time_steps,
            resolution,
            initial_time,
            variables,
            constraints,
            cost_function,
            expressions,
            parameters,
            initial_conditions,
            pm_model)

    end

end

function CanonicalModel(::Type{T},
                         sys::PSY.System,
                         optimizer::Union{Nothing,JuMP.OptimizerFactory};
                         kwargs...) where {T<:PM.AbstractPowerModel}

    sequential_runs = get(kwargs, :sequential_runs, false)
    user_defined_model = get(kwargs, :JuMPmodel, nothing)
    ini_con = get(kwargs, :initial_conditions, DICKDA())
    parameters = get(kwargs, :parameters, false)
    forecast = get(kwargs, :forecast, true)
    jump_model = _pass_abstract_jump(optimizer, parameters, user_defined_model)
    initial_time = get(kwargs,
                       :initial_time,
                       PSY.get_forecasts_initial_time(sys))

    if forecast
        horizon = PSY.get_forecasts_horizon(sys)
        time_steps = 1:horizon
        resolution = PSY.get_forecasts_resolution(sys)
    else
        resolution = PSY.get_forecasts_resolution(sys)
        time_steps = 1:1
    end

    bus_numbers = sort([PSY.get_number(b) for b in PSY.get_components(PSY.Bus, sys)])

    return _canonical_init(bus_numbers,
                           jump_model,
                           optimizer,
                           T,
                           time_steps,
                           resolution,
                           initial_time,
                           parameters,
                           sequential_runs,
                           ini_con)

end

function InitialCondition(canonical::CanonicalModel,
                          device::T,
                          access_ref::Symbol,
                          value::Float64,
                          cache::Union{Nothing, Type{<:AbstractCache}}=nothing) where T <: PSY.Device

    if model_has_parameters(canonical)
        return InitialCondition(device,
                                UpdateRef{JuMP.VariableRef}(access_ref),
                                PJ.add_parameter(canonical.JuMPmodel, value),
                                cache)
    else
        !hasfield(T, access_ref) && error("Device of of type $(T) doesn't contain
                                            the field $(access_ref)")
        return InitialCondition(device,
                                UpdateRef{T}(access_ref),
                                value,
                                cache)
    end

end

function get_ini_cond(canonical_model::CanonicalModel, key::ICKey)
    return get(canonical_model.initial_conditions, key, Vector{InitialCondition}())
end

# Var_ref
function get_value(canonical::CanonicalModel, ref::UpdateRef{JuMP.VariableRef})
    return var(canonical, ref.access_ref)
end

# param_ref
function get_value(canonical::CanonicalModel, ref::UpdateRef{PJ.ParameterRef})
    for (k, v) in canonical.parameters
        if k.access_ref == ref.access_ref
            return v
        end
    end
    return
end

_variable_type(cm::CanonicalModel) = JuMP.variable_type(cm.JuMPmodel)
model_time_steps(canonical_model::CanonicalModel) = canonical_model.time_steps
model_resolution(canonical_model::CanonicalModel) = canonical_model.resolution
model_has_parameters(canonical_model::CanonicalModel) = canonical_model.parametrized
model_runs_sequentially(canonical_model::CanonicalModel) = canonical_model.sequential_runs
model_initial_time(canonical_model::CanonicalModel) = canonical_model.initial_time
#Internal Variables, Constraints and Parameters accessors
vars(canonical_model::CanonicalModel) = canonical_model.variables
cons(canonical_model::CanonicalModel) = canonical_model.constraints
var(canonical_model::CanonicalModel, name::Symbol) = canonical_model.variables[name]
con(canonical_model::CanonicalModel, name::Symbol) = canonical_model.constraints[name]
par(canonical_model::CanonicalModel, param_reference::UpdateRef) = canonical_model.parameters[param_reference]
exp(canonical_model::CanonicalModel, name::Symbol) = canonical_model.expressions[name]
get_initial_conditions(canonical_model::CanonicalModel) = canonical_model.initial_conditions


text:
function add_solver_xml!(::NEOSSolver{:CPLEX, :NL}, m::NLModel)
	# Add user options
	param_string = ""
	for key in keys(m.solver.params)
		param_string *= "set $(key) $(m.solver.params[key])\n"
	end
	m.xmlmodel = replace(m.xmlmodel, r"<options>.*</options>"s, "<options><![CDATA[$(param_string)]]></options>")
end


text:
using DashBootstrapComponents

badges = html_span([
    dbc_badge("Primary", color = "primary", className = "me-1"),
    dbc_badge("Secondary", color = "secondary", className = "me-1"),
    dbc_badge("Success", color = "success", className = "me-1"),
    dbc_badge("Warning", color = "warning", className = "me-1"),
    dbc_badge("Danger", color = "danger", className = "me-1"),
    dbc_badge("Info", color = "info", className = "me-1"),
    dbc_badge("Light", color = "light", text_color = "dark", className = "me-1"),
    dbc_badge("Dark", color = "dark"),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
module Slider

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Slider (Clockwise)" => Ahorn.EntityPlacement(
        Maple.Slider,
        "point",
        Dict{String, Any}(
            "clockwise" => true
        )
    )
)

Ahorn.editingOptions(entity::Maple.Slider) =  Dict{String, Any}(
    "surface" => Maple.slider_surfaces
)

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Slider, room::Maple.Room)
    clockwise = get(entity.data, "clockwise", false)
    dir = clockwise ? 1 : -1

    x, y = Int(entity.data["x"]), Int(entity.data["y"])
    radius = 12

    Ahorn.drawArrow(ctx, x + radius, y, x + radius, y + 0.001 * dir, Ahorn.colors.selection_selected_fc, headLength=6)
    Ahorn.drawArrow(ctx, x - radius, y, x - radius, y + 0.001 * -dir, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.selection(entity::Maple.Slider)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.Slider, room::Maple.Room) = Ahorn.drawCircle(ctx, 0, 0, 12, (1.0, 0.0, 0.0, 1.0))

end


text:
"""
    AbstractContext

An imputation context records summary information about missing data for an imputation algorithm.
All `AbstractContext`s are callable with a function, which allows us to write code like:

```julia
context() do c
    # My imputation code using a clean context
end
```

This do-block will pass a fresh context to your code and apply the `on_complete` function on
the resulting data and context state. By default, `on_complete` will throw an
[ImputeError](@ref) if we have too many missing values.
"""
abstract type AbstractContext end

# We implement a version of copy for all contexts which reconstructs the context from the
# raw fields.
Base.copy(ctx::T) where {T <: AbstractContext} = T(fieldvalues(ctx)...)

"""
    ismissing(ctx::AbstractContext, x) -> Bool

Uses `ctx.is_missing` to determine if x is missing. If x is a `NamedTuple` or an `AbstractArray`
then `ismissing` will return true if `ctx.is_missing` returns true for any element.
The ctx.count is increased whenever whenever we return true and if `ctx.count / ctx.num`
exceeds our `ctx.limit` we throw an `ImputeError`

# Arguments
* `ctx::Context`: the contextual information about missing information.
* `x`: the value to check (may be an single values, abstract array or row)
"""
function Base.ismissing(ctx::AbstractContext, x)
    was_missing = if isa(x, NamedTuple)
        any(ctx.is_missing, Tuple(x))
    elseif isa(x, AbstractArray)
        any(ctx.is_missing, x)
    else
        ctx.is_missing(x)
    end

    missing_update!(ctx, was_missing)

    return was_missing
end

"""
    findfirst(ctx::AbstractContext, data::AbstractVector) -> Int

Returns the first non-missing index in `data`.

# Arguments
* `ctx::AbstractContext`: the context to pass into `ismissing`
* `data::AbstractVector`: the data array to search

# Returns
* `Int`: the first index in `data` that isn't missing
"""
function Base.findfirst(ctx::AbstractContext, data::AbstractVector)
    return findfirst(x -> !ismissing(ctx, x), data)
end

"""
    findlast(ctx::AbstractContext, data::AbstractVector) -> Int

Returns the last non-missing index in `data`.

# Arguments
* `ctx::AbstractContext`: the context to pass into `ismissing`
* `data::AbstractVector`: the data array to search

# Returns
* `Int`: the last index in `data` that isn't missing
"""
function Base.findlast(ctx::AbstractContext, data::AbstractVector)
    return findlast(x -> !ismissing(ctx, x), data)
end

"""
    findnext(ctx::AbstractContext, data::AbstractVector) -> Int

Returns the next non-missing index in `data`.

# Arguments
* `ctx::AbstractContext`: the context to pass into `ismissing`
* `data::AbstractVector`: the data array to search

# Returns
* `Int`: the next index in `data` that isn't missing
"""
function Base.findnext(ctx::AbstractContext, data::AbstractVector, idx::Int)
    return findnext(x -> !ismissing(ctx, x), data, idx)
end

@auto_hash_equals mutable struct Context <: AbstractContext
    num::Int
    count::Int
    limit::Float64
    is_missing::Function
    on_complete::Function
end

"""
    Context

Records base information about the missing data and assume all observations are equally
weighted.

# Keyword Arguments
* `n::Int`: number of observations
* `count::Int`: number of missing values found
* `limit::Float64`: portion of total values allowed to be imputed (should be between 0.0 and 1.0).
* `is_missing::Function`: must return a Bool indicating if the value counts as missing
* `on_complete::Function`: a function to run when imputation is complete
"""
function Context(;
    limit::Float64=0.1,
    is_missing::Function=ismissing,
    on_complete::Function=complete
)
    return Context(0, 0, limit, is_missing, on_complete)
end

function Base.empty(ctx::Context)
    _ctx = copy(ctx)
    _ctx.num = 0
    _ctx.count = 0

    return _ctx
end

function missing_update!(ctx::Context, was_missing)
    ctx.num += 1

    if was_missing
        ctx.count += 1
    end
end

function complete(ctx::Context, data)
    missing_ratio = ctx.count / ctx.num

    if missing_ratio > ctx.limit
        throw(ImputeError(
            "More than $(ctx.limit * 100)% of values were missing ($missing_ratio)."
        ))
    end

    return data
end

@auto_hash_equals mutable struct WeightedContext <: AbstractContext
    num::Int
    s::Float64
    limit::Float64
    is_missing::Function
    on_complete::Function
    wv::AbstractWeights
end

"""
    WeightedContext(wv; limit=1.0, is_missing=ismissing, on_complete=complete)

Records information about the missing data relative to a set of weights.
This context type can be useful if some missing observation are more important than others
(e.g., more recent observations in time series datasets)

# Arguments
* `wv::AbstractWeights`: a set of statistical weights to use when evaluating the importance
  of each observation. Will be accumulated during imputation.

# Keyword Arguments
* `num::Int`: number of observations
* `s::Float64`: sum of the weights of missing values
* `limit::Float64`: portion of total values allowed to be imputed (should be between 0.0 and 1.0).
* `is_missing::Function`: returns a Bool if the value counts as missing
* `on_complete::Function`: a function to run when imputation is complete
"""
function WeightedContext(
    wv::AbstractWeights;
    limit::Real=1.0,
    is_missing::Function=ismissing,
    on_complete::Function=complete
)
    return WeightedContext(0, 0.0, limit, is_missing, on_complete, wv)
end

function Base.empty(ctx::WeightedContext)
    _ctx = copy(ctx)
    _ctx.num = 0
    _ctx.s = 0.0

    return _ctx
end

function missing_update!(ctx::WeightedContext, was_missing)
    ctx.num += 1

    if was_missing
        ctx.s += ctx.wv[ctx.num]
    end
end

function complete(ctx::WeightedContext, data)
    missing_ratio = ctx.s / sum(ctx.wv)

    if missing_ratio > ctx.limit
        throw(ImputeError(
            "More than $(ctx.limit * 100)% of weighted values were missing ($missing_ratio)."
        ))
    end

    return data
end

#=
Define our callable methods for each context. Once we drop 1.0 we should be able to just
define this on the `AbstractContext`.
=#
for T in (Context, WeightedContext)
    @eval begin
        function (ctx::$T)(f::Function)
            _ctx = empty(ctx)
            return ctx.on_complete(_ctx, f(_ctx))
        end
    end
end


text:
@testset "Test composition order" begin
  Qx, x = FlintQQ["x"]
  K, a = NumberField(x^6 + 108, "a")
  A = automorphisms(K)
  for f in A
    for g in A
      @test g(f(a)) == compose(f, g)(a)
      @test g(f(a)) == (f * g)(a)
    end
  end
end


text:
using Test, DiscreteExteriorCalculus, DiscretePDEs
using UniqueVectors: UniqueVector
using LinearAlgebra: eigen

file_name = joinpath(@__DIR__, "sphere.geo")
geo_write!(file_name, characteristic_length_factor=.3,
    footer="""
    Sphere(1) = {0, 0, 0, 1, -Pi/2, Pi/2, 2*Pi};
    """)
initialize!()
gmsh_open!(file_name)
N, K = 3, 3
mesh!(K)

node_tags, points, tcomp = get_triangulated_complex(N, K)
@test typeof(tcomp) <: TriangulatedComplex{N, K}
comp = tcomp.complex
orient!(comp)
m = Metric(N)
mesh = Mesh(tcomp, circumcenter(m))

laplacian = differential_operator(m, mesh, "Δ", 1, true)

vals, vects = eigen(collect(laplacian))
inds = sortperm(vals)
vals, vects = vals[inds], vects[:, inds]

# solutions are spherical harmonics (or combinations of them since they
# are degenerate)
correct_modes = vcat([[l*(l+1) for _ in -l:l] for l in 0:3]...)

@testset "laplacian on a sphere" begin
    @test isapprox(vals[1:length(correct_modes)], correct_modes, rtol=1e-2)
end

if false # set to true to display plots with Gmsh
    comp_points = UniqueVector([c.points[1] for c in comp.cells[1]])
    ordering = [findfirst(isequal(p), comp_points) for p in points]
    for i in 1:length(correct_modes)
        v = vects[:,i]
        add_field!("Eigenvector", node_tags, v[ordering])
    end
    gui!()
end
rm(file_name)


text:
using ROOTFramework, ROOT, Cxx
import JSON

# Init GUI, open output file and create canvas

rootgui()
set_default_palette(:viridis)

tfile = open(TFile, "out.root", "recreate")

canvas = @cxxnew TCanvas(pointer("histplots"), pointer("Histogram Plots"))
@cxx canvas->Divide(2, 2)

# TH1

hist1 = THxx(-20:0.2:20, "hist1", "Hist")

@time for i in 1:10000
    push!(hist1, 5.0 * randn())
end
@cxx canvas->cd(1)
draw(hist1)

ROOTFramework.fit(hist1, "gaus")
draw(hist1)

write(tfile, hist1)
#hist_dict1 = JSON.parse(rootjson(hist1))

# TH2

set_default_palette(:viridis)

hist2 = THxx((-10:0.5:10, -10:0.5:10), "hist2", "Hist2")

@time for i in 1:100000
    push!(hist2, (3.0 * randn(), 3.0 * randn()))
end
@cxx canvas->cd(2)
draw(hist2, "lego2")
@cxx canvas->cd(3)
draw(hist2, "colz1")

#write(tfile, hist2)
#hist_dict2 = JSON.parse(rootjson(hist2))

# TH3

hist3 = THxx((-10:0.5:10, -10:0.5:10, -10:0.5:10), "hist3", "Hist")

@time for i in 1:100000
    push!(hist3, (3.0 * randn(), 3.0 * randn(), 3.0 * randn()))
end
@cxx canvas->cd(4)
draw(hist3)

#write(tfile, hist3)
#hist3_dict = JSON.parse(rootjson(hist3))

# Done

#write(tfile, canvas)

#close(tfile)

#=
new_TBrowser();
=#


text:
# An N-dimensional point
struct Point{N,T} <: StaticVector{N,T}
    data::NTuple{N,T}

    function Point{N,T}(x::NTuple{N,T}) where {N,T}
        return new{N,T}(x)
    end

    function Point{N,T}(x::NTuple{N,Any}) where {N,T}
        return new{N,T}(StaticArrays.convert_ntuple(T, x))
    end
end

const Point_2D = Point{2}
const Point_3D = Point{3}

# NOTE: This is pretty much an exact copy of the Point in GeometryBasics.
# (https://github.com/JuliaGeometry/GeometryBasics.jl/blob/master/src/fixed_arrays.jl)

# Constructors & Conversions
# -------------------------------------------------------------------------------------------------
# Construct from vector with known size
function Point{N}(x::AbstractVector{T}) where {N,T}
    @assert N <= length(x)
    return Point{N,T}(ntuple(i -> x[i], Val(N)))
end

# Construct from vector with known size, and convert to type T₁
function Point{N,T₁}(x::AbstractVector{T₂}) where {N,T₁,T₂}
    @assert N <= length(x)
    return Point{N,T₁}(ntuple(i -> convert(T₁, x[i]), Val(N)))
end

# Construct from a Tuple
function Point(x::T) where {N,T <: Tuple{Vararg{Any,N}}}
    return Point{N,StaticArrays.promote_tuple_eltype(T)}(x)
end

# Construct from a Tuple with known size
function Point{N}(x::T) where {N,T <: Tuple}
    return Point{N,StaticArrays.promote_tuple_eltype(T)}(x)
end

# Construct from a Point
@generated function (::Type{Point{N,T}})(p::Point) where {N,T}
    idx = [:(p[$i]) for i in 1:N]
    return quote
        $(Point){N,T}($(idx...))
    end
end

# Convert from a Point
@generated function Base.convert(::Type{Point{N,T}}, p::Point) where {N,T}
    idx = [:(p[$i]) for i in 1:N]
    return quote
        $(Point){N,T}($(idx...))
    end
end

# Convert Tuple to Point
function Base.convert(::Type{Point{N,T}}, x::NTuple{N,T}) where {N,T}
    return Point{N,T}(x)
end
function Base.convert(::Type{Point{N,T}}, x::Tuple) where {N,T}
    return Point{N,T}(convert(NTuple{N,T}, x))
end

# Base (and similar_type)
# -------------------------------------------------------------------------------------------------
@generated function StaticArrays.similar_type(::Type{SV}, ::Type{T},
                                              s::Size{N}) where {SV <: Point,T,N}
    return if length(N) === 1
        Point{N[1],T}
    else
        StaticArrays.default_similar_type(T, s(), Val{length(N)})
    end
end
function Base.broadcasted(f, a::AbstractArray{T}, b::Point) where {T <: Point}
    return Base.broadcasted(f, a, (b,))
end
Base.@propagate_inbounds function Base.getindex(p::Point{N,T}, i::Int) where {N,T}
    return p.data[i]
end
Base.Tuple(p::Point) = p.data

# Methods
# -------------------------------------------------------------------------------------------------
@inline distance(𝐩₁::Point, 𝐩₂::Point) = norm(𝐩₁ - 𝐩₂)
@inline distance²(𝐩₁::Point, 𝐩₂::Point) = norm²(𝐩₁ - 𝐩₂)
@inline Base.isapprox(𝐩₁::Point, 𝐩₂::Point) = distance²(𝐩₁, 𝐩₂) < (1e-5)^2 # 100 nm
@inline midpoint(𝐩₁::Point, 𝐩₂::Point) = (𝐩₁ + 𝐩₂)/2
@inline norm²(𝐩::Point) = 𝐩 ⋅ 𝐩

# Sort points based on their distance from a given point
sortpoints(p::Point, points::Vector{<:Point}) = points[sortperm(distance².(Ref(p), points))]
function sortpoints!(p::Point_2D, points::Vector{<:Point_2D})
    permute!(points, sortperm(distance².(Ref(p), points)))
    return nothing
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
@testset "simple" begin
    @testset "simple.gallery.`Ast Pattern`" begin
        isnothing(x) = x !== nothing

        rmlines(a::Expr) =
            Expr(a.head, filter(isnothing, map(rmlines, a.args))...)

        rmlines(::LineNumberNode) = nothing

        rmlines(a::Any) = a

        ex = quote
            function f(x, y=5)
                x + 1
                x + 2
                x + "3"
                x + "4"
                y + x
            end
        end

        @test @match rmlines(ex) begin
            quote
                function f(x, y=$default)
                    x + $a
                    x + $b
                    x + $c
                    x + $d
                    y + x
                end
            end => (a, b, c, d, default) == (1, 2, "3", "4", 5)
        end

        ex = quote
            struct A{T} end
            struct B end
            struct C end
            struct D{G} end
        end

        @match rmlines(ex) begin
            Do(names=[]) &&
                quote
                    $(
                        Many(
                            (
                                :(struct $name end)
                                || :(struct $name{$(_...)} end)
                            ) &&
                            Do(push!(names, name))
                        )...
                    )
                end => names == [:A, :B, :C, :D]
        end
    end
    @testset "gallery.simple.`Array & Tuple`" begin

        @test @match [1, 2, 3] begin
            [1, a..., b] => (a, b) == ([2], 3)
        end
        # Attention:
        # `a` is not a copy of a part of original array.
        # In Python, `_, *a, b = [1, 2, 3]` means `a` is a copy,
        #   while in Julia, `a` here is a `view`[1].
        # [1] view: https://docs.julialang.org/en/latest/devdocs/subarrays/index.html

        @test @match (1, 2, 3) begin
            (1, 2, a) => a == 3
        end
        # Note that, you cannot pack a tuple,
        # and this is invalid
        # ```
        # @match (1, 2, 3) begin
        #   (1, a..., b) => ...
        # end
        # ```

        # nested match
        @test @match [1, [2, 3], (4, 5)] begin
            [1, [2, a], (b, 5), tail...] =>
                a == 3 && b == 4 && isempty(tail)
        end
    end
    @testset "gallery.simple.`Or Pattern`" begin
        # If I want number 1 or 2:
        function is_it(x)
            @match x begin
                1 || 2 => true
                _     => false
            end
        end
        @test is_it(1)
        @test is_it(2)
        @test !is_it(3)

        # Patterns can be nested
        function nested_or(x)
            @match x begin
                1 || [1 || 2] =>  true
                _             =>  false
            end
        end
        @test nested_or(1)
        @test nested_or([1])
        @test nested_or([2])
        @test !nested_or([1, 1])
    end

    @testset "gallery.simple.`And Pattern`" begin
        @test @match 1 begin
            a && 1 : 10 => a == 1
            _           => false
        end
    end
end


text:
using Distributions
using Roots

# Some global variables that will stay constant
alpha = 0.013
alpha_q = (1-(1-alpha)^3)
b = 0.0124
d = 0.00822
beta = 0.98
gamma = 1.0
sigma = 2.0

# The default wage distribution: a discretized log normal
log_wage_mean, wage_grid_size, max_wage = 20, 200, 170
w_vec = linspace(1e-3, max_wage, wage_grid_size + 1)
logw_dist = Normal(log(log_wage_mean), 1)
cdf_logw = cdf(logw_dist, log(w_vec))
pdf_logw = cdf_logw[2:end] - cdf_logw[1:end-1]
p_vec = pdf_logw ./ sum(pdf_logw)
w_vec = (w_vec[1:end-1] + w_vec[2:end]) / 2

"""
Compute the reservation wage, job finding rate and value functions of the
workers given c and tau.

"""
function compute_optimal_quantities(c::Float64, tau::Float64)
  mcm = McCallModel(alpha_q,
                    beta,
                    gamma,
                    c-tau,           # post-tax compensation
                    sigma,
                    collect(w_vec-tau),   # post-tax wages
                    p_vec)

  w_bar, V, U = compute_reservation_wage(mcm, return_values=true)
  lmda = gamma * sum(p_vec[w_vec-tau .> w_bar])

  return w_bar, lmda, V, U
end

"""
Compute the steady state unemployment rate given c and tau using optimal
quantities from the McCall model and computing corresponding steady state
quantities

"""
function compute_steady_state_quantities(c::Float64, tau::Float64)
  w_bar, lmda, V, U = compute_optimal_quantities(c, tau)

  # Compute steady state employment and unemployment rates
  lm = LakeModel(lambda=lmda, alpha=alpha_q, b=b, d=d)
  x = rate_steady_state(lm)
  e_rate, u_rate = x

  # Compute steady state welfare
  w = sum(V .* p_vec .* (w_vec - tau .> w_bar)) / sum(p_vec .* (w_vec - tau .> w_bar))
  welfare = e_rate .* w + u_rate .* U

  return e_rate, u_rate, welfare
end

"""
Find tax level that will induce a balanced budget.

"""
function find_balanced_budget_tax(c::Float64)
  function steady_state_budget(t::Float64)
    e_rate, u_rate, w = compute_steady_state_quantities(c, t)
    return t - u_rate * c
  end

  tau = fzero(steady_state_budget, 0.0, 0.9 * c)

  return tau
end

# Levels of unemployment insurance we wish to study
Nc = 60
c_vec = linspace(5.0, 140.0, Nc)

tax_vec = Array(Float64, Nc,)
unempl_vec = Array(Float64, Nc,)
empl_vec = Array(Float64, Nc,)
welfare_vec = Array(Float64, Nc,)

for i = 1:Nc
  t = find_balanced_budget_tax(c_vec[i])
  e_rate, u_rate, welfare = compute_steady_state_quantities(c_vec[i], t)
  tax_vec[i] = t
  unempl_vec[i] = u_rate
  empl_vec[i] = e_rate
  welfare_vec[i] = welfare
end

titles = ["Unemployment" "Employment" "Tax" "Welfare"]
plot(c_vec, [unempl_vec empl_vec tax_vec welfare_vec],
    color=:blue, lw=2, alpha=0.7, title=titles, legend=:none, layout=(2,2))


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 9285d6f0-9e71-11eb-0411-1928431042b7
begin
	import Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ f455bcc9-02bb-4d20-961d-4da33624a575
begin
	Pkg.add("Plots")
	Pkg.add("PlutoUI")
	Pkg.add("StaticArrays")
	Pkg.add("BenchmarkTools")
	using Revise
	using Plots
	using PlutoUI
	using WormlikeChain
	using StaticArrays
	using BenchmarkTools
end

# ╔═╡ f4577ff0-f627-4f54-87d3-27f4518e51db
md"""
## Defining a bead
A bead is defined by a potential energy function on it's neighbors and its dimension
"""

# ╔═╡ 362e0a98-ea1a-49b9-a3b4-a50dcdf6d10a
perbead_param_keys= (:k,:kθ)

# ╔═╡ cb945668-43cf-48aa-893c-2b968af8549e
function chain_pe(p,b,n,perbead_params,global_params)
	rp= (p .- b)
	rn= (n .- b)
	dp= sqrt(sum(rp.^2))
	dn= sqrt(sum(rn.^2))
	cosθ= sum(rp .* rn)/dp/dn
	perbead_params.k*1//2*((dn-1)^2) + perbead_params.kθ*(cosθ+1)
end

# ╔═╡ 4e5f3d5e-0cab-4fea-905d-4be4acef6526
const Ndims=2

# ╔═╡ cca6e1b6-84a2-4a1e-ab1c-a31d4b8ecb74
beaddef= BeadDefinition(chain_pe,  perbead_param_keys, (), Val(Ndims))

# ╔═╡ 3a27b649-d3f5-47b7-860d-fb67f1ac9f1f
md"""
BeadDefinition automatically differentiates `chain_pe` to get a `chain_force`
`chain_force` can then be compiled into optimized assembly

"""

# ╔═╡ 4242953f-0f86-47a1-beab-bd2583a3dbd5
md"""
## Defining a chain
A chain is just a circle of beads.
"""

# ╔═╡ 09e78f9a-20f4-4ba6-8cb5-018897e96be1
const Nbeads=50

# ╔═╡ 8a6b422a-f288-4cbf-b8f0-ca9d41d8282a
begin
	r= (Nbeads/2π)
	θs= range(0.0, step= 2π/Nbeads, length=Nbeads)
	xs= (Nbeads/2π) .* cos.(θs)
	ys= (Nbeads/2π) .* sin.(θs)
	pos= [xs ys]
	# pos= zeros(Nbeads,Ndims)
	# pos[:,2] .= range(0.0, step= 1.0, length=Nbeads)
	vel= zero(pos)# .+ [0.0 -0.1]
	scatter(pos[:,1],pos[:,2])
	#beadparams= fill((k=0.0, kθ=0.0),Nbeads)
	beadparams= fill((k=10.0, kθ=4.4),Nbeads)
	# beadparams[Nbeads]= (k=0.0, kθ= 0.0)# turn off interation to split ends
	# beadparams[1]= (k=3.5, kθ= 0.0)
end

# ╔═╡ e12c77f3-870c-40ea-aec8-ae01750e26ee
chain1= Chain(beaddef,beadparams, (), pos,vel)

# ╔═╡ 440db016-35da-4e7e-9255-af9c42995a13
md"""
## Defining an additional potential energy term
"""

# ╔═╡ 78bfb5e5-9ec3-4c9f-9a43-329bf2087fb8
wallparams= (k= 1.0, L= 3*Nbeads/π)

# ╔═╡ 21a861fc-06d6-45bd-9a18-5f3573582693
function harmonicwallpe(pos,time,params)
	pe= 0
	for x in (pos[1])
		pe+= 1//2*params.k*(max(x-params.L/2,0))^2
		pe+= 1//2*params.k*(min(x+params.L/2,0))^2
	end
	pe
end

# ╔═╡ 06a89ca3-2768-4551-85bd-b24160b9b0e8
interactions = [((1,i),) for i in 1:Nbeads]

# ╔═╡ 3736c303-1e9c-44f2-b92a-09e632c19bfd
walldef= SpecificForce(harmonicwallpe, wallparams, interactions, Val(Ndims))

# ╔═╡ ed818b30-5e04-472b-b54f-6d2d5e77280f
md"""
## Putting it all together in a ChainSystem
First create an empty system, then add in chains and extra forces
"""

# ╔═╡ 7b971a08-e603-47e9-bf73-1c53b2f15d38
empty_system= ChainSystem(0.0,beaddef,())

# ╔═╡ c84c7c2b-1791-44ed-861d-b7076f3311b2
chain1system= append(empty_system,chain1)

# ╔═╡ fa9b09f7-6ab7-44a1-9e88-9ec1e108510f
system= append(chain1system,walldef)

# ╔═╡ b2f82994-caed-42dd-8fb8-3ad965446e35
sum(force_pe(system, system.init_pos,0.0)[1][:,2])

# ╔═╡ 4cd004e3-beeb-47a9-94a2-a4738c637ecc
randn(10,2)

# ╔═╡ eb54332d-3f8e-4884-92b1-949267412239
md"""
## Compiling the kernel and running a sim
Kernel get recompiled and speciallized for the new types created
"""

# ╔═╡ bc1f2f08-707d-4899-8366-21ddd07dfda3

# ╔═╡ 7379a001-450a-420f-a012-e2a7b92d1e1c
begin
	invmass= 1.0
    γ= 0.0
    invβ= 1.0
    VNbeads= Val(size(system.init_pos)[1])
    VNdims= Val(size(system.init_pos)[2])
    step0= UInt64(0)
    time0= 0.0
    Vrngkey= Val(UInt64(123))
    Δt= 0.01
    VChainbounds= Val(system.chainbounds)
    VNsteps= Val(1000)
end

# ╔═╡ 5b49a595-f83e-4084-9277-8a72e6d1e9be
VNbeads

# ╔═╡ bb8e8fcb-ebc2-495b-a938-31a303a30278
VNdims

# ╔═╡ a88ade2a-66a6-443d-91a3-d48526905c86
totalframes= 100

# ╔═╡ eb2dcad9-3ec6-4d43-9588-6f0c1d2239ed
begin
	trajectory= zeros(Nbeads,Ndims,totalframes)
	pes= zeros(totalframes)
 	kes= zeros(totalframes)
	simpos= copy(system.init_pos)
    simvel= copy(system.init_vel)
	simvel[1,:] = [0.1 0.1]
 	simvel .+= [0 1]

	for step in 1:totalframes
		trajectory[:,:,step] .= simpos
		f, pe= force_pe(system,simpos,0.0)
		kes[step] = 1//2*sum((simvel .+ (Δt/2).*f) .^ 2)
		pes[step] = pe
    	refcpukernel!(simpos, simvel,
                system.sp_forces, system.sp_interactions, system.sp_params,
                system.beaddef.chain_force, system.perbead_params,
			    system.global_params,
                invmass, γ, invβ, VNbeads, step0, time0,
                Vrngkey, Δt, VNdims, VChainbounds, VNsteps)
 	end

end

# ╔═╡ 587def7f-a48b-4140-8571-bcf9151e0409
@bind seeframe Slider(1:totalframes)

# ╔═╡ 036a7fd3-11a0-4aae-a5d5-b819c1bf114b
scatter(trajectory[:,1,seeframe],trajectory[:,2,seeframe]; xlims = (-wallparams.L/2, wallparams.L/2), ylims = (-wallparams.L/2,wallparams.L/2))

# ╔═╡ 39fe110d-20b0-4362-b4d7-1f51a9fbf128
seeframe

# ╔═╡ e3db67ad-11f6-42f1-aaad-d81f3a772397
begin
	plot(kes)
	plot!(pes)
	plot!(kes+pes)
end

# ╔═╡ c3c2b1b4-4e73-4b68-8f2d-39fe1db6f550
plotly()

# ╔═╡ Cell order:
# ╟─f4577ff0-f627-4f54-87d3-27f4518e51db
# ╠═362e0a98-ea1a-49b9-a3b4-a50dcdf6d10a
# ╠═cb945668-43cf-48aa-893c-2b968af8549e
# ╠═4e5f3d5e-0cab-4fea-905d-4be4acef6526
# ╠═cca6e1b6-84a2-4a1e-ab1c-a31d4b8ecb74
# ╟─3a27b649-d3f5-47b7-860d-fb67f1ac9f1f
# ╟─4242953f-0f86-47a1-beab-bd2583a3dbd5
# ╠═09e78f9a-20f4-4ba6-8cb5-018897e96be1
# ╠═8a6b422a-f288-4cbf-b8f0-ca9d41d8282a
# ╠═e12c77f3-870c-40ea-aec8-ae01750e26ee
# ╟─440db016-35da-4e7e-9255-af9c42995a13
# ╠═78bfb5e5-9ec3-4c9f-9a43-329bf2087fb8
# ╠═21a861fc-06d6-45bd-9a18-5f3573582693
# ╠═06a89ca3-2768-4551-85bd-b24160b9b0e8
# ╠═3736c303-1e9c-44f2-b92a-09e632c19bfd
# ╟─ed818b30-5e04-472b-b54f-6d2d5e77280f
# ╠═7b971a08-e603-47e9-bf73-1c53b2f15d38
# ╠═c84c7c2b-1791-44ed-861d-b7076f3311b2
# ╠═fa9b09f7-6ab7-44a1-9e88-9ec1e108510f
# ╠═b2f82994-caed-42dd-8fb8-3ad965446e35
# ╠═4cd004e3-beeb-47a9-94a2-a4738c637ecc
# ╟─eb54332d-3f8e-4884-92b1-949267412239
# ╠═bc1f2f08-707d-4899-8366-21ddd07dfda3
# ╠═7379a001-450a-420f-a012-e2a7b92d1e1c
# ╠═5b49a595-f83e-4084-9277-8a72e6d1e9be
# ╠═bb8e8fcb-ebc2-495b-a938-31a303a30278
# ╠═eb2dcad9-3ec6-4d43-9588-6f0c1d2239ed
# ╠═a88ade2a-66a6-443d-91a3-d48526905c86
# ╠═036a7fd3-11a0-4aae-a5d5-b819c1bf114b
# ╠═39fe110d-20b0-4362-b4d7-1f51a9fbf128
# ╠═587def7f-a48b-4140-8571-bcf9151e0409
# ╠═e3db67ad-11f6-42f1-aaad-d81f3a772397
# ╠═9285d6f0-9e71-11eb-0411-1928431042b7
# ╠═f455bcc9-02bb-4d20-961d-4da33624a575
# ╠═c3c2b1b4-4e73-4b68-8f2d-39fe1db6f550


text:
import .Flux
Flux.trainable(bn::InvertibleBatchNorm) = (bn.b, bn.logs)
Flux.@functor InvertibleBatchNorm


text:
# ---
# title: 1530. Number of Good Leaf Nodes Pairs
# id: problem1530
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/number-of-good-leaf-nodes-pairs/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree and an integer `distance`. A pair of two
# different **leaf** nodes of a binary tree is said to be good if the length of
# **the shortest path** between them is less than or equal to `distance`.
#
# Return _the number of good leaf node pairs_ in the tree.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/09/e1.jpg)
#
#
#
#     Input: root = [1,2,3,null,4], distance = 3
#     Output: 1
#     Explanation: The leaf nodes of the tree are 3 and 4 and the length of the shortest path between them is 3. This is the only good pair.
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/07/09/e2.jpg)
#
#
#
#     Input: root = [1,2,3,4,5,6,7], distance = 3
#     Output: 2
#     Explanation: The good pairs are [4,5] and [6,7] with shortest path = 2. The pair [4,6] is not good because the length of ther shortest path between them is 4.
#
#
# **Example 3:**
#
#
#
#     Input: root = [7,1,4,6,null,5,3,null,null,null,null,null,2], distance = 3
#     Output: 1
#     Explanation: The only good pair is [2,5].
#
#
# **Example 4:**
#
#
#
#     Input: root = [100], distance = 1
#     Output: 0
#
#
# **Example 5:**
#
#
#
#     Input: root = [1,1,1], distance = 2
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the `tree` is in the range `[1, 2^10].`
#   * Each node's value is between `[1, 100]`.
#   * `1 <= distance <= 10`
#
#
## @lc code=start
using LeetCode

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


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Plan <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    publisher::Any # spec type: Union{ Nothing, String } # spec name: publisher
    product::Any # spec type: Union{ Nothing, String } # spec name: product
    promotionCode::Any # spec type: Union{ Nothing, String } # spec name: promotionCode

    function Plan(;name=nothing, publisher=nothing, product=nothing, promotionCode=nothing)
        o = new()
        validate_property(Plan, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(Plan, Symbol("publisher"), publisher)
        setfield!(o, Symbol("publisher"), publisher)
        validate_property(Plan, Symbol("product"), product)
        setfield!(o, Symbol("product"), product)
        validate_property(Plan, Symbol("promotionCode"), promotionCode)
        setfield!(o, Symbol("promotionCode"), promotionCode)
        o
    end
end # type Plan

const _property_map_Plan = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("publisher")=>Symbol("publisher"), Symbol("product")=>Symbol("product"), Symbol("promotionCode")=>Symbol("promotionCode"))
const _property_types_Plan = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("publisher")=>"String", Symbol("product")=>"String", Symbol("promotionCode")=>"String")
Base.propertynames(::Type{ Plan }) = collect(keys(_property_map_Plan))
Swagger.property_type(::Type{ Plan }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Plan[name]))}
Swagger.field_name(::Type{ Plan }, property_name::Symbol) =  _property_map_Plan[property_name]

function check_required(o::Plan)
    true
end

function validate_property(::Type{ Plan }, name::Symbol, val)
end


text:
function templateTestTmapfoldl1(arr::AbstractArray)::Bool
	init = randInit();
	return mapfoldl(x->2x, +, arr; init = init) == tmapfoldl(x->2x, +, arr; init = init);
end

function testTmapfoldl1Stable()::Bool
	return templateTestTmapfoldl1(randArray());
end

function testTmapfoldl1Unstable()::Bool
	arr = randStrArray();
	str = lipsum();
	init = randstring(rand(1:1000));
	return mapfoldl(x->str[1:x], *, arr, init = init) == tmapfoldl(x->str[1:x], *, arr; init = init);
end

function testTmapfoldl1Offset()::Bool
	return templateTestTmapfoldl1(offsetArray());
end

function testTmapfoldl1Empty()::Bool
	return templateTestTmapfoldl1([]);
end

#####

function templateTestTmapfoldl2(arr::AbstractArray)::Bool
	return mapfoldl(x->2x, +, arr) == tmapfoldl(x->2x, +, arr);
end

function testTmapfoldl2Stable()::Bool
	return templateTestTmapfoldl2(randArray());
end

function testTmapfoldl2Unstable()::Bool
	arr = randStrArray();
	str = lipsum();
	return mapfoldl(x->str[1:x], *, arr) == tmapfoldl(x->str[1:x], *, arr);
end

function testTmapfoldl2Offset()::Bool
	return templateTestTmapfoldl2(offsetArray());
end

function testTmapfoldl2Empty()
	tmapfoldl(x->2x, (x, y)->x, []);
end


text:
### Exports ###

export ConvNConfig, ConvNState, ConvNeurons

### Types ###

ConvNConfig = NConfig{Tuple{Int,Int,Int}}
const ConvNState = NState{Array{Float32, 3}, Array{Bool, 3}}
const ConvNeurons = Neurons{ConvNConfig, ConvNState}
ConvNeurons(size::Tuple{Int,Int,Int}; a=0.02, b=0.2, c=-65, d=8, thresh=30, θRate=0.5, traceRate=0.75, boostRate=0.01) =
  ConvNeurons(
    ConvNConfig(
      size,
      a,
      b,
      c,
      d,
      thresh,
      θRate,
      traceRate,
      boostRate),
    ConvNState(
      fill(c, size...),
      zeros(Float32, size...),
      zeros(Float32, size...),
      zeros(Bool, size...),
      zeros(Float32, size...),
      zeros(Float32, size...),
      zeros(Float32, size...)))

### Methods ###

function Base.show(io::IO, neurons::ConvNeurons)
  print(io, "ConvNeurons $(neurons.conf.size)")
end
# TODO: Make this actually random?
function randnode(::Type{ConvNeurons}, size::Tuple{Int,Int,Int}=(rand(1:20,3)...,))
  ConvNeurons(size)
end


text:
include("./genericOrdering.jl")

@enum RRBSetting sop
include("./sopOrdering.jl")
include("./sopRRB.jl")


text:
# Tests of rankings

using Stats
using Base.Test

a = [1.0, 2.0, 2.0, 3.0, 4.0, 4.0, 4.0, 5.0]
x = [3.0, 1.0, 2.0, 4.0, 4.0, 2.0, 5.0, 4.0]  # x is a permutated version of a

@test ordinalrank(a) == [1, 2, 3, 4, 5, 6, 7, 8]
@test ordinalrank(x) == [4, 1, 2, 5, 6, 3, 8, 7]

@test competerank(a) == [1, 2, 2, 4, 5, 5, 5, 8]
@test competerank(x) == [4, 1, 2, 5, 5, 2, 8, 5]

@test denserank(a) == [1, 2, 2, 3, 4, 4, 4, 5]
@test denserank(x) == [3, 1, 2, 4, 4, 2, 5, 4]

@test tiedrank(a) == [1.0, 2.5, 2.5, 4.0, 6.0, 6.0, 6.0, 8.0]
@test tiedrank(x) == [4.0, 1.0, 2.5, 6.0, 6.0, 2.5, 8.0, 6.0]


text:
@testset "Symmetric/Hermitian rules" begin
    @testset "$(SymHerm)(::AbstractMatrix{$T}, :$(uplo))" for
        SymHerm in (Symmetric, Hermitian),
        T in (Float64, ComplexF64),
        uplo in (:U, :L)

        N = 3
        @testset "frule" begin
            x = randn(T, N, N)
            Δx = randn(T, N, N)
            # can't use frule_test here because it doesn't yet ignore nothing tangents
            Ω = SymHerm(x, uplo)
            Ω_ad, ∂Ω_ad = frule((Zero(), Δx, Zero()), SymHerm, x, uplo)
            @test Ω_ad == Ω
            ∂Ω_fd = jvp(_fdm, z -> SymHerm(z, uplo), (x, Δx))
            @test ∂Ω_ad ≈ ∂Ω_fd
        end
        @testset "rrule" begin
            x = randn(T, N, N)
            ∂x = randn(T, N, N)
            ΔΩ = randn(T, N, N)
            @testset "back(::$MT)" for MT in (Matrix, LowerTriangular, UpperTriangular)
                rrule_test(SymHerm, MT(ΔΩ), (x, ∂x), (uplo, nothing))
            end
            @testset "back(::Diagonal)" begin
                rrule_test(SymHerm, Diagonal(ΔΩ), (x, Diagonal(∂x)), (uplo, nothing))
            end
        end
    end
    @testset "$(f)(::$(SymHerm){$T}) with uplo=:$uplo" for f in (Matrix, Array),
        SymHerm in (Symmetric, Hermitian),
        T in (Float64, ComplexF64),
        uplo in (:U, :L)

        N = 3
        x = SymHerm(randn(T, N, N), uplo)
        Δx = randn(T, N, N)
        ∂x = SymHerm(randn(T, N, N), uplo)
        ΔΩ = f(SymHerm(randn(T, N, N), uplo))
        frule_test(f, (x, Δx))
        frule_test(f, (x, SymHerm(Δx, uplo)))
        rrule_test(f, ΔΩ, (x, ∂x))
    end
end


text:
# This is for testing display in IJulia. A huge image should
# be automatically restricted before sending to IJulia for display.

using Images

url = "http://upload.wikimedia.org/wikipedia/commons/9/97/The_Earth_seen_from_Apollo_17.jpg"

const savedir = joinpath(tempdir(), "Images")
if !isdir(savedir)
    mkdir(savedir)
end

const savename = joinpath(savedir, "The_Earth_seen_from_Apollo_17.jpg")
if !isfile(savename)
    download(url, savename)
end

img = imread(savename)