text:
module ArrayTimings

using Printf
using BenchmarkTools
using BenchmarkTools: prettymemory
using ArrayTools
using LinearAlgebra
using LazyAlgebra

#
# The benchmarks show that `all_indices` is as fast as `eachindex`.
#

flatten(A::AbstractArray{T}) where {T} =
    copyto!(Vector{T}(undef, length(A)), A)

function prettytime(ns)
    a = abs(ns)
    if a < 1e3
        scale, units = 1e0, "ns"
    elseif a < 1e6
        scale, units = 1e-3, "μs"
    elseif a < 1e9
        scale, units = 1e-6, "ms"
    else
        scale, units = 1e-9, "s "
    end
    return string(@sprintf("%8.3f", scale*ns), " ", units)
end

function prettyflops(Gflops)
    a = abs(Gflops)
    if true || a >= 1
        scale, units = 1e0, "Gflops"
    else
        scale, units = 1e3, "Mflops"
    end
    return string(@sprintf("%8.3f", scale*Gflops), " ", units)
end

function prt(str::AbstractString,
             nops::Integer,
             trial::BenchmarkTools.Trial; kwds...)
    prt(stdout, str, nops, trial; kwds...)
end

function prt(io::IO,
             str::AbstractString,
             nops::Integer,
             trial::BenchmarkTools.Trial;
             width::Integer = 0, mark::Bool=false)
    best = minimum(trial)
    print(io, str, repeat(" ", max(1, width - length(str))),
          prettytime(time(best)), "  ",
          prettyflops(nops/time(best)),
          " (", allocs(best), " allocations: ",
          prettymemory(memory(best)), ")")
    mark && print(io,  " [*]")
    print(io,  "\n")
end

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

sum1(A::AbstractArray) = begin
    s = zero(eltype(A))
    for v in A
        s += v
    end
    return s
end

sum2(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds for v in A
        s += v
    end
    return s
end

sum3(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds @simd for v in A
        s += v
    end
    return s
end

sum4(A::AbstractArray) = begin
    s = zero(eltype(A))
    for i in eachindex(A)
        s += A[i]
    end
    return s
end

sum5(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds for i in eachindex(A)
        s += A[i]
    end
    return s
end

sum6(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds @simd for i in eachindex(A)
        s += A[i]
    end
    return s
end

sum7(A::AbstractArray) = begin
    s = zero(eltype(A))
    for i in all_indices(A)
        s += A[i]
    end
    return s
end

sum8(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds for i in all_indices(A)
        s += A[i]
    end
    return s
end

sum9(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds @simd for i in all_indices(A)
        s += A[i]
    end
    return s
end

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

function dot4(x::AbstractArray{T,N},
              y::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    s = zero(T)
    for i in eachindex(x, y)
        s += x[i]*y[i]
    end
    return s
end

function dot5(x::AbstractArray{T,N},
              y::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    s = zero(T)
    @inbounds for i in eachindex(x, y)
        s += x[i]*y[i]
    end
    return s
end

function dot6(x::AbstractArray{T,N},
              y::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    s = zero(T)
    @inbounds @simd for i in eachindex(x, y)
        s += x[i]*y[i]
    end
    return s
end

function dot7(x::AbstractArray{T,N},
              y::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    s = zero(T)
    for i in all_indices(x, y)
        s += x[i]*y[i]
    end
    return s
end

function dot8(x::AbstractArray{T,N},
              y::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    s = zero(T)
    @inbounds for i in all_indices(x, y)
        s += x[i]*y[i]
    end
    return s
end

function dot9(x::AbstractArray{T,N},
              y::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    s = zero(T)
    @inbounds @simd for i in all_indices(x, y)
        s += x[i]*y[i]
    end
    return s
end

function unsafe_dot(x::AbstractArray{T,N},
                    y::AbstractArray{T,N}) where {T<:AbstractFloat,N}
    s = zero(T)
    @inbounds @simd for i in eachindex(x)
        s += x[i]*y[i]
    end
    return s
end

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

norm2_1(A::AbstractArray) = begin
    s = zero(eltype(A))
    for v in A
        s += abs2(v)
    end
    return sqrt(s)
end

norm2_2(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds for v in A
        s += abs2(v)
    end
    return sqrt(s)
end

norm2_3(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds @simd for v in A
        s += abs2(v)
    end
    return sqrt(s)
end

norm2_4(A::AbstractArray) = begin
    s = zero(eltype(A))
    for i in eachindex(A)
        s += abs2(A[i])
    end
    return sqrt(s)
end

norm2_5(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds for i in eachindex(A)
        s += abs2(A[i])
    end
    return sqrt(s)
end

norm2_6(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds @simd for i in eachindex(A)
        s += abs2(A[i])
    end
    return sqrt(s)
end

norm2_7(A::AbstractArray) = begin
    s = zero(eltype(A))
    for i in all_indices(A)
        s +=abs2(A[i])
    end
    return sqrt(s)
end

norm2_8(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds for i in all_indices(A)
        s += abs2(A[i])
    end
    return sqrt(s)
end

norm2_9(A::AbstractArray) = begin
    s = zero(eltype(A))
    @inbounds @simd for i in all_indices(A)
        s += abs2(A[i])
    end
    return sqrt(s)
end

norm2_9(::Type{T}, A::AbstractArray) where {T<:AbstractFloat} = begin
    local s::T = zero(T)
    @inbounds @simd for i in all_indices(A)
        s += abs2(A[i])
    end
    return sqrt(s)
end
#------------------------------------------------------------------------------

function test_sum()
    dims = (100, 100)
    inds = (2:99, 2:99)
    for T in (Float32, Float64)
        A = rand(T, dims)
        V = view(A, inds...)
        X = rand(T, dims)
        Y = rand(T, dims)
        W = view(X, inds...)
        Z = view(Y, inds...)
        a = flatten(A)
        x = flatten(X)
        y = flatten(Y)

        println("\nTest SUM on a regular array (T=$T):")
        n = length(A)  # number of operations
        prt("  sum(a) ", n, @benchmark(sum($a)))
        prt("  sum(A) ", n, @benchmark(sum($A)))
        prt("  sum1(A)", n, @benchmark(sum1($A)))
        prt("  sum2(A)", n, @benchmark(sum2($A)))
        prt("  sum3(A)", n, @benchmark(sum3($A)))
        prt("  sum4(A)", n, @benchmark(sum4($A)))
        prt("  sum5(A)", n, @benchmark(sum5($A)))
        prt("  sum6(A)", n, @benchmark(sum6($A)))
        prt("  sum7(A)", n, @benchmark(sum7($A)))
        prt("  sum8(A)", n, @benchmark(sum8($A)))
        prt("  sum9(a)", n, @benchmark(sum9($a)); mark=true)
        prt("  sum9(A)", n, @benchmark(sum9($A)); mark=true)

        println("\nTest SUM on a sub-array (T=$T):")
        n = length(V)  # number of operations
        prt("  sum ", n, @benchmark(sum($V)))
        prt("  sum1", n, @benchmark(sum1($V)))
        prt("  sum2", n, @benchmark(sum2($V)))
        prt("  sum3", n, @benchmark(sum3($V)))
        prt("  sum4", n, @benchmark(sum4($V)))
        prt("  sum5", n, @benchmark(sum5($V)))
        prt("  sum6", n, @benchmark(sum6($V)))
        prt("  sum7", n, @benchmark(sum7($V)))
        prt("  sum8", n, @benchmark(sum8($V)))
        prt("  sum9", n, @benchmark(sum9($V)); mark=true)
    end
end

function test_dot()
    dims = (100, 100)
    inds = (2:99, 2:99)
    for T in (Float32, Float64)
        A = rand(T, dims)
        V = view(A, inds...)
        X = rand(T, dims)
        Y = rand(T, dims)
        W = view(X, inds...)
        Z = view(Y, inds...)
        a = flatten(A)
        x = flatten(X)
        y = flatten(Y)

        println("\nTest DOT on regular arrays (T=$T):")
        n = 2*length(X)  # number of operations
        prt("  x'*y      ", n, @benchmark($x'*$y))
        prt("  dot(x,y)  ", n, @benchmark(dot($x,$y)))
        prt("  dot4      ", n, @benchmark(dot4($X,$Y)))
        prt("  dot5      ", n, @benchmark(dot5($X,$Y)))
        prt("  dot6      ", n, @benchmark(dot6($X,$Y)))
        prt("  dot7      ", n, @benchmark(dot7($X,$Y)))
        prt("  dot8      ", n, @benchmark(dot8($X,$Y)))
        prt("  dot9      ", n, @benchmark(dot9($X,$Y)); mark=true)
        prt("  unsafe_dot", n, @benchmark(unsafe_dot($X,$Y)))
        prt("  vdot      ", n, @benchmark(vdot($X,$Y)))

        println("\nTest DOT on sub-arrays (T=$T):")
        n = 2*length(W)  # number of operations
        prt("  dot4      ", n, @benchmark(dot4($W,$Z)))
        prt("  dot5      ", n, @benchmark(dot5($W,$Z)))
        prt("  dot6      ", n, @benchmark(dot6($W,$Z)))
        prt("  dot7      ", n, @benchmark(dot7($W,$Z)))
        prt("  dot8      ", n, @benchmark(dot8($W,$Z)))
        prt("  dot9      ", n, @benchmark(dot9($W,$Z)); mark=true)
        prt("  unsafe_dot", n, @benchmark(unsafe_dot($W,$Z)))
        prt("  vdot      ", n, @benchmark(vdot($W,$Z)))
    end
end

function test_norm()
    dims = (100, 100)
    inds = (2:99, 2:99)
    for T in (Float32, Float64)
        A = rand(T, dims)
        V = view(A, inds...)
        X = rand(T, dims)
        Y = rand(T, dims)
        W = view(X, inds...)
        Z = view(Y, inds...)
        a = flatten(A)
        x = flatten(X)
        y = flatten(Y)

        println("\nTest NORM2 on a regular array (T=$T):")
        n = 2*length(A)  # number of operations
        prt("  norm(a,2)         ", n, @benchmark(norm($a,2)))
        prt("  norm(A,2)         ", n, @benchmark(norm($A,2)))
        prt("  norm2_1(A)        ", n, @benchmark(norm2_1($A)))
        prt("  norm2_2(A)        ", n, @benchmark(norm2_2($A)))
        prt("  norm2_3(A)        ", n, @benchmark(norm2_3($A)))
        prt("  norm2_4(A)        ", n, @benchmark(norm2_4($A)))
        prt("  norm2_5(A)        ", n, @benchmark(norm2_5($A)))
        prt("  norm2_6(A)        ", n, @benchmark(norm2_6($A)))
        prt("  norm2_7(A)        ", n, @benchmark(norm2_7($A)))
        prt("  norm2_8(A)        ", n, @benchmark(norm2_8($A)))
        prt("  norm2_9(a)        ", n, @benchmark(norm2_9($a)); mark=true)
        prt("  norm2_9(A)        ", n, @benchmark(norm2_9($A)); mark=true)
        prt("  vnorm2(A)         ", n, @benchmark(vnorm2($A)))
        prt("  norm2_9(Float32,A)", n, @benchmark(norm2_9(Float32,$A)); mark=true)
        prt("  vnorm2(Float32,A) ", n, @benchmark(vnorm2(Float32,$A)))
        prt("  norm2_9(Float64,A)", n, @benchmark(norm2_9(Float64,$A)); mark=true)
        prt("  vnorm2(Float64,A) ", n, @benchmark(vnorm2(Float64,$A)))

        println("\nTest NORM2 on a sub-array (T=$T):")
        n = 2*length(V)  # number of operations
        prt("  norm(V,2)         ", n, @benchmark(norm($V,2)))
        prt("  norm2_1(V)        ", n, @benchmark(norm2_1($V)))
        prt("  norm2_2(V)        ", n, @benchmark(norm2_2($V)))
        prt("  norm2_3(V)        ", n, @benchmark(norm2_3($V)))
        prt("  norm2_4(V)        ", n, @benchmark(norm2_4($V)))
        prt("  norm2_5(V)        ", n, @benchmark(norm2_5($V)))
        prt("  norm2_6(V)        ", n, @benchmark(norm2_6($V)))
        prt("  norm2_7(V)        ", n, @benchmark(norm2_7($V)))
        prt("  norm2_8(V)        ", n, @benchmark(norm2_8($V)))
        prt("  norm2_9(V)        ", n, @benchmark(norm2_9($V)); mark=true)
        prt("  vnorm2(V)         ", n, @benchmark(vnorm2($V)))
        prt("  norm2_9(Float32,V)", n, @benchmark(norm2_9(Float32,$V)); mark=true)
        prt("  vnorm2(Float32,V) ", n, @benchmark(vnorm2(Float32,$V)))
        prt("  norm2_9(Float64,V)", n, @benchmark(norm2_9(Float64,$V)); mark=true)
        prt("  vnorm2(Float64,V) ", n, @benchmark(vnorm2(Float64,$V)))
    end
end

function test_update()
    dims = (100, 100)
    inds = (2:99, 2:99)
    for T in (Float32, Float64)
        X = rand(T, dims)
        Y = similar(X)
        U = view(X, inds...)
        V = view(Y, inds...)
        α = T(0.01)

        println("\nTest UPDATE on regular arrays (T=$T):")
        n = 2*length(X)  # number of operations
        vfill!(Y,0)
        prt("  Y .= Y .+ α.*X   ", n, @benchmark($Y .+= $α.*$X))
        vfill!(Y,0)
        prt("  vupdate!(Y,α,X)  ", n, @benchmark(vupdate!($Y,$α,$X)))

        println("\nTest UPDATE on sub-arrays (T=$T):")
        n = 2*length(U)  # number of operations
        vfill!(V,0)
        prt("  V .= V .+ α.*V   ", n, @benchmark($V .+= $α.*$U))
        vfill!(V,0)
        prt("  vupdate!(V,α,U)  ", n, @benchmark(vupdate!($V,$α,$U)))
    end
end

function test_combine()
    dims = (100, 100)
    inds = (2:99, 2:99)
    for T in (Float32, Float64)
        X = rand(T, dims)
        Y = rand(T, dims)
        Z = similar(X)
        U = view(X, inds...)
        V = view(Y, inds...)
        W = view(Z, inds...)
        α = T(0.3)
        β = T(π)

        println("\nTest COMBINE on regular arrays (T=$T):")
        n = 3*length(X)  # number of operations
        prt("  α*X + β*Y             ", n, @benchmark($α*$X + $β*$Y))
        prt("  vcombine(α,X,β,Y)     ", n, @benchmark(vcombine($α,$X,$β,$Y)))
        prt("  vcombine!(Z,α,X,β,Y)  ", n, @benchmark(vcombine!($Z,$α,$X,$β,$Y)); mark=true)

        println("\nTest COMBINE on sub-arrays (T=$T):")
        n = 3*length(U)  # number of operations
        prt("  α*U + β*V             ", n, @benchmark($α*$U + $β*$V))
        prt("  vcombine(α,U,β,V)     ", n, @benchmark(vcombine($α,$U,$β,$V)))
        prt("  vcombine!(W,α,U,β,V)  ", n, @benchmark(vcombine!($W,$α,$U,$β,$V)); mark=true)
    end
end

end # module

if !isinteractive()
    println("\nIn these tests, the asterisk [*] marks the version ",
            "which should be among the fastest ones and which should ",
            "require no additional memory.")
    ArrayTimings.test_sum()
    ArrayTimings.test_dot()
    ArrayTimings.test_norm()
    ArrayTimings.test_update()
    ArrayTimings.test_combine()
end

nothing


text:
module TestMissing

using Test
using RLEVectors

standard_rle = RLEVector(["a","d","c","b"],[2,4,6,8])
missing_ok_rle = RLEVector(Vector{Union{Missing, String}}(["a","d","c","b"]), [2,4,6,8])
with_missing_rle = RLEVector(Vector{Union{Missing, String}}(["a",missing,"c","b"]), [2,4,6,8])

@testset "Staring with simple vector" begin
    x = standard_rle
    @test disallowmissing(x) == standard_rle
    @test allowmissing(x) == missing_ok_rle
end

@testset "Staring with missings possible but not present" begin
    x = missing_ok_rle
    @test allowmissing(x) == x
    @test disallowmissing(x) == standard_rle
end

@testset "Staring with missing possible" begin
    x = with_missing_rle
    @test_throws MethodError disallowmissing(x)
    x = with_missing_rle
    @test typeof(allowmissing(x)) == typeof(x)
end

end # module

ext: jl


text:
#!/usr/local/bin/julia

function linear_search(array::Array{Int64,1}, key::Int64)
    for i in array
        if i == key
            return true
    return false
        end
    end
end

# Benchmarks
# julia> @timev linear_search([1:100000000...], 10000000)
#  11.617622 seconds (100.00 M allocations: 3.033 GiB, 54.06% gc time)
# elapsed time (ns): 11617622167
# gc time (ns):      6280118240
# bytes allocated:   3256270960
# pool allocs:       99999500
# non-pool GC allocs:9
# malloc() calls:    2
# realloc() calls:   15
# GC pauses:         6
# full collections:  4
# true


text:
using Jchemo

n = 6 ; p = 7 ; q = 2 ; m = 3
Xtrain = rand(n, p) ; Ytrain = rand(n, q)
Xtest = rand(m, p) ; Ytest = rand(m, q)

nlv = 5 ;
fm = plskern(Xtrain, Ytrain; nlv = nlv) ;
pnames(fm)
fm.T
fm.P

Jchemo.summary(fm, Xtrain).explvar

fm.T
Jchemo.transform(fm, Xtrain)
Jchemo.transform(fm, Xtrain; nlv = 1)
Jchemo.transform(fm, Xtest)

Jchemo.coef(fm).B
Jchemo.coef(fm).int
Jchemo.coef(fm; nlv = 2).B
Jchemo.coef(fm; nlv = 2).int
Jchemo.coef(fm; nlv = 0).B
Jchemo.coef(fm; nlv = 0).int

Jchemo.predict(fm, Xtest).pred
Jchemo.predict(fm, Xtest; nlv = nlv).pred

Jchemo.predict(fm, Xtest; nlv = 0:3).pred
Jchemo.predict(fm, Xtest; nlv = 0).pred

pred = Jchemo.predict(fm, Xtest).pred ;
msep(pred, Ytest)

### Weighted PLS

w = collect(1:n)
fm = plskern(Xtrain, Ytrain, w; nlv = nlv) ;


text:
# Creating the images that are used in the
# Readme/documentation.

# ImageMagick, together with FileIO is used to save the Images
using ComplexPortraits
using ImageMagick
using FileIO

using Colors

@ComplexPortraits.import_normal

function save_image(name, img)
  name = joinpath("images", name * ".png")
  save(File(format"PNG", name), img)
end

docfct(z) = (z-1)/(z*z + z +1)
docid(z) = z

function doc_portrait(name; kwargs...)
  img = portrait(
    -2.0 + 2.0im, 2.0 - 2.0im, docfct; kwargs..., no_pixels=(800,800))
  save_image(name, img)
  img = portrait(
    -2.0 + 2.0im, 2.0 - 2.0im, docid; kwargs..., no_pixels=(800,800))
  save_image(name * "_id", img)
  return nothing
end

function cs_fade_out(base_color_func)
  return function(z, fz)
    color = HSV(base_color_func(z,fz))
    absz = abs(z)
    return HSV(color.h,
      absz ≤ 1.5 ? color.s : color.s*max(1-2*(absz-1.5), 0.0),
      absz ≤ 1.75 ? color.v : min(1.0, color.v + 2*(absz-1.75)))
  end
end

if true

cs_intro = cs_fade_out(cs_j(colormap=nist_colors()))
save_image("intro01",
  portrait(-2.0 + 2.0im, 2.0 - 2.0im, z -> z;
  no_pixels=300, point_color=cs_intro))
save_image("intro02",
  portrait(-2.0 + 2.0im, 2.0 - 2.0im, z -> z^2;
  no_pixels=300, point_color=cs_intro))

save_image("wheel01",
  portrait(-2.0 + 2.0im, 2.0 - 2.0im, z -> z;
  no_pixels=300, point_color=cs_fade_out(cs_p())))

save_image("wheel02",
  portrait(-2.0 + 2.0im, 2.0 - 2.0im, z -> z;
  no_pixels=300, point_color=cs_fade_out(cs_p(colormap=nist_colors()))))

doc_portrait("cs_c01", point_color=cs_c())
doc_portrait("cs_c02", point_color=cs_c(phaseres=10))

doc_portrait("cs_d01", point_color=cs_d())

doc_portrait("cs_e01", point_color=cs_e())
doc_portrait("cs_e02", point_color=cs_e(logabsres=10))

doc_portrait("cs_j01", point_color=cs_j())
doc_portrait("cs_j02", point_color=cs_j(jumps=[60/180*pi, 240/180*pi]))

doc_portrait("cs_m01", point_color=cs_m())
doc_portrait("cs_m02", point_color=cs_m(logabsres=10))

doc_portrait("cs_p01", point_color=cs_p())

doc_portrait("cs_q01", point_color=cs_q())

doc_portrait("cs_gridReIm01", point_color=cs_gridReIm())
doc_portrait("cs_gridReIm02", point_color=cs_gridReIm(reim_b=0.8))

doc_portrait("cs_gridReIm_logabs01",
  point_color=cs_gridReIm_logabs())
doc_portrait("cs_gridReIm_logabs02",
  point_color=cs_gridReIm_logabs(abs_b=0.9, reim_b=0.8))

doc_portrait("cs_gridReIm_abs01", point_color=cs_gridReIm_abs())

test_img = load("./images/Portraits.png")

doc_portrait("cs_useImage01",
  point_color=cs_useImage(
    test_img;
    img_upperleft=0.0 + 0.20im,
    img_lowerright=0.80 + 0.0im,
    is_hv_periodic=(false, true)
    ))

doc_portrait("cs_useImage02",
  point_color=cs_useImage(
    test_img;
    img_upperleft=0.0 + 0.20im,
    img_lowerright=0.80 + 0.0im,
    is_hv_periodic=(true, true)
    ))

end

doc_portrait("cs_stripes01",
  point_color=cs_stripes(
   [5.0 + 0.0im],
   [HSV(0.0, 0.0, 1.0), HSV(0.0, 0.0, 0.0)]))

doc_portrait("cs_stripes02",
  point_color=cs_stripes(
   [5.0 + 0.0im, 0.0 + 5.0im],
   [HSV(0.0, 0.0, 1.0), HSV(0.0, 0.0, 0.0),
    HSV(0.0, 1.0, 1.0), HSV(240.0, 1.0, 1.0)]))

doc_portrait("cs_stripes03",
  point_color=cs_stripes(
    exp(-1.0im*pi/4).*[2.0 + 0.0im, 0.0 + 2.0im, 2.0 + 2.0im],
    [HSV(0.0, 0.0, 1.0), HSV(0.0, 0.0, 0.0),
     HSV(0.0, 1.0, 1.0), HSV(60.0, 1.0, 1.0),
     HSV(120.0, 1.0, 1.0), HSV(180.0, 1.0, 1.0),
     HSV(240.0, 1.0, 1.0), HSV(300.0, 1.0, 1.0)]))

doc_portrait("cs_stripes04",
  point_color=cs_angle_abs_stripes(
   [5.0 + 0.0im],
   [HSV(0.0, 0.0, 1.0), HSV(0.0, 0.0, 0.0)]))

doc_portrait("cs_stripes05",
  point_color=cs_angle_abs_stripes(
   [5.0 + 0.0im, 0.0 + 5.0im],
   [HSV(0.0, 0.0, 1.0), HSV(0.0, 0.0, 0.0),
    HSV(0.0, 1.0, 1.0), HSV(240.0, 1.0, 1.0)]))

doc_portrait("cs_stripes06",
  point_color=cs_angle_abs_stripes(
    exp(-1.0im*pi/4).*[2.0 + 0.0im, 0.0 + 2.0im, 2.0 + 2.0im],
    [HSV(0.0, 0.0, 1.0), HSV(0.0, 0.0, 0.0),
     HSV(0.0, 1.0, 1.0), HSV(60.0, 1.0, 1.0),
     HSV(120.0, 1.0, 1.0), HSV(180.0, 1.0, 1.0),
     HSV(240.0, 1.0, 1.0), HSV(300.0, 1.0, 1.0)]))

# vim:syn=julia:cc=79:fdm=indent:sw=2:ts=2:


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LAMMPS_jll
using Base
using Base: UUID
include(path) = Base.include(@__MODULE__, path)
include(joinpath("..", ".pkg", "platform_augmentation.jl"))
import JLLWrappers

JLLWrappers.@generate_main_file_header("LAMMPS")
JLLWrappers.@generate_main_file("LAMMPS", UUID("5b3ab26d-9607-527c-88ea-8fe5ba57cafe"))
end  # module LAMMPS_jll


text:
const MIPHyperParamConfig = NamedTuple{(:M, :λ),Tuple{Float64,Float64}}

struct MIP <: Algorithm
    optimizer_factory
    hypograph::Bool
    regularizer::Union{Nothing,ObjectivePenalty}
    cp_warm_start::Bool
    hyperparameter_set::Array{MIPHyperParamConfig}
end
function MIP(
    optimizer_factory;
    hypograph::Bool=false,
    regularizer=nothing,
    cp_warm_start=true)
    return MIP(
        optimizer_factory,
        hypograph,
        regularizer,
        cp_warm_start,
        [(M = 2.0^M_exp_val, λ = λ_val) for λ_val in (0.0,), M_exp_val in -1:9]
    )
end

has_regularizer(algo::MIP) = algo.regularizer !== nothing
get_regularizer(algo::MIP) = algo.regularizer

function train(learner::Learner{MIP})
    @assert learner.validation_data !== nothing

    formulation_fn = learner.algorithm.hypograph ? formulate_data_point_hypograph_big_m : formulate_data_point_graph_big_m
    if learner.algorithm.regularizer !== nothing
        _add_objective_penalty!(learner, learner.algorithm.regularizer)
    end

    hyperparam_results = Dict{Any,TrainingResults}()
    for (M, λ) in learner.algorithm.hyperparameter_set
        learner.param_lb .= -M
        learner.param_ub .= M
        learner.offset_lb = -M
        learner.offset_ub = M
        if λ > 0
            # TODO: Refactor so that this is less hacky
            push!(learner.objective_penalties, L2Regularizer(λ))
        end
        opt_model = formulate_big_m(learner, formulation_fn)
        if learner.algorithm.cp_warm_start
            contextless_train_data = DataSet([DataPoint(Float64[], dp.bid_price, dp.pay_price) for dp in learner.train_data.samples])
            cp_learner = Learner(contextless_train_data, nothing, nothing, MIP(learner.algorithm.optimizer_factory), Float64[], Float64[], learner.offset_lb, learner.offset_ub, ObjectivePenalty[])
            cp_opt_model = formulate_big_m(cp_learner)
            cp_jump_model = cp_opt_model.jump_model
            cp_solved_opt_model = optimize!(cp_learner, cp_opt_model)
            @assert cp_solved_opt_model.jump_model === cp_jump_model
            if JuMP.primal_status(cp_jump_model) == MOI.FEASIBLE_POINT
                for i in 1:feature_count(learner.train_data)
                    JuMP.set_start_value(opt_model.linear_model.weights[i], 0.0)
                end
                JuMP.set_start_value(opt_model.linear_model.offset, cp_solved_opt_model.linear_model.offset)
                for i in 1:sample_count(learner.train_data)
                    JuMP.set_start_value(opt_model.vs[i], cp_solved_opt_model.vs[i])
                    JuMP.set_start_value(opt_model.ys[i], cp_solved_opt_model.ys[i])
                    for j in 1:size(opt_model.zs, 2)
                        JuMP.set_start_value(opt_model.zs[i,j], cp_solved_opt_model.zs[i,j])
                    end
                end
            else
                @warn "CP warmstart did not produce a feasible solution"
            end
        end
        solved_opt_model = optimize!(learner, opt_model)
        _verify_solution(learner.train_data, solved_opt_model)

        hyperparam_results[(M, λ)] = TrainingResults(learner, solved_opt_model.linear_model, reward_bound=JuMP.objective_bound(solved_opt_model.jump_model))
    end
    best_validation_reward = -Inf
    best_validation_key = nothing

    iter_results = IterationResult[]
    for (key, result) in hyperparam_results
        reward = result.final_validation_eval.reward
        if reward > best_validation_reward
            best_validation_reward = reward
            best_validation_key = key
        end
        push!(iter_results, IterationResult(result.final_is_eval, result.final_os_eval, nothing, Dict("key" => key, "validation_eval" => result.final_validation_eval)))
    end
    best_result = hyperparam_results[best_validation_key]
    best_result.iteration_evals = iter_results
    return best_result
end

max_stars_repo_name: joehuchette/reserve-price-optimization

max_issues_repo_name: joehuchette/reserve-price-optimization

max_forks_repo_name: joehuchette/reserve-price-optimization


text:
function compareSolversOut(n, nruns; maxtime=Inf)

    t1 = time()

    # first, force a compile
    tab = compareSolvers(1000,2)

    tab = compareSolvers(n, nruns, maxtime=maxtime)
    fn = string("compSolvers", n , "x", nruns, ".csv")
    writecsv(fn,tab)

    t2 = time() - t1
    println("Took ", t2, " seconds.")
end

function compareSolvers(n, nruns; maxtime=Inf)

    akpwBuild = Array(Float64,nruns)
    akpwSolve = Array(Float64,nruns)
    primBuild = Array(Float64,nruns)
    primSolve = Array(Float64,nruns)
    augBuild = Array(Float64,nruns)
    augSolve = Array(Float64,nruns)

    nedges = Array(Int64,nruns)
    nlist = Array(Int64,nruns)
    ilist = Array(Int64,nruns)

    for i in 1:nruns
        a = wtedChimera(n,i)
        la = lap(a)

        b = randn(n)
        b = b - mean(b)

        nlist[i] = n
        ilist[i] = i
        nedges[i] = nnz(a)

        tic()
        f = treeSolver(a,akpw)
        t = toq()
        akpwBuild[i] = t

        tic()
        x = f(b, maxtime=maxtime)
        t = toq()
        akpwSolve[i] = t

        tic()
        f = treeSolver(a,prim)
        t = toq()
        primBuild[i] = t

        tic()
        x = f(b, maxtime=maxtime)
        t = toq()
        primSolve[i] = t

        tic()
        f = augTreeLapSolver(la)
        t = toq()
        augBuild[i] = t

        tic()
        x = f(b, maxtime=maxtime)
        t = toq()
        augSolve[i] = t

        print(".")
    end

    tab = ["nlist"; nlist]
    tab = [tab ["ilist"; ilist]]
    tab = [tab ["nedges"; nedges]]
    tab = [tab ["akpwBuild"; akpwBuild]]
    tab = [tab ["akpwSolve"; akpwSolve]]
    tab = [tab ["primBuild"; primBuild]]
    tab = [tab ["primSolve"; primSolve]]
    tab = [tab ["augBuild"; augBuild]]
    tab = [tab ["augSolve"; augSolve]]

    #=
    tab = makeColumn(:nlist)
    tab = [tab makeColumn(:ilist)]
    tab = [tab makeColumn(:nedges)]
    tab = [tab makeColumn(:akpwBuild)]
    tab = [tab makeColumn(:akpwSolve)]
    tab = [tab makeColumn(:primBuild)]
    tab = [tab makeColumn(:primSolve)]
    tab = [tab makeColumn(:augBuild)]
    tab = [tab makeColumn(:augSolve)]
    =#

    return tab
end

function treeSolver(a,treeAlg;maxtime=Inf)
    t = treeAlg(a)
    la = lap(a)
    lt = lap(t)
    f = pcgLapSolver(la,lt;maxtime=maxtime)
end


text:
#################### Chains ####################

## Constructors ##

# Constructor to handle a vector of vectors.
Chains(val::AbstractVector{<:AbstractVector{<:Union{Missing, Real}}}, args...; kwargs...) =
	Chains(copy(reduce(hcat, val)'), args...; kwargs...)

# Constructor to handle a 1D array.
Chains(val::AbstractVector{<:Union{Missing, Real}}, args...; kwargs...) =
	Chains(reshape(val, :, 1, 1), args...; kwargs...)

# Constructor to handle a 2D array
Chains(val::AbstractMatrix{<:Union{Missing, Real}}, args...; kwargs...) =
    Chains(reshape(val, size(val, 1), size(val, 2), 1), args...; kwargs...)

# Constructor to handle parameter names that are not Symbols.
function Chains(
    val::AbstractArray{<:Union{Missing,Real},3},
    parameter_names::AbstractVector,
    args...;
    kwargs...
)
    return Chains(val, Symbol.(parameter_names), args...; kwargs...)
end

# Generic chain constructor.
function Chains(
    val::AbstractArray{<:Union{Missing, Real},3},
    parameter_names::AbstractVector{Symbol} = Symbol.(:param_, 1:size(val, 2)),
    name_map = (parameters = parameter_names,);
    start::Int = 1,
    thin::Int = 1,
    evidence = missing,
    info::NamedTuple = NamedTuple()
)
    # Make sure that we have a `:parameters` index and # Copying can avoid state mutation.
    _name_map = initnamemap(name_map)

    # Preclean the name_map of names that aren't in the
    # parameter_names vector.
    for names in _name_map
        filter!(x -> x ∈ parameter_names, names)
    end

    # Store unassigned variables.
    unassigned = Set(Symbol[])

    # Check that all parameters are assigned.
    for param in parameter_names
        if all(param ∉ names for names in _name_map)
            push!(unassigned, param)
        end
    end

    # Assign all unassigned parameter names.
    append!(_name_map[:parameters], unassigned)

    # Construct the AxisArray.
    arr = AxisArray(val;
                    iter = range(start, step=thin, length=size(val, 1)),
                    var = parameter_names,
                    chain = 1:size(val, 3))

    # Create the new chain.
    return Chains(arr, evidence, _name_map, info)
end

"""
    Chains(c::Chains, section::Union{Symbol,String})
    Chains(c::Chains, sections)

Return a new chain with only a specific `section` or multiple `sections` pulled out.

# Examples
```jldoctest
julia> chn = Chains(rand(100, 2, 1), [:a, :b], Dict(:internals => [:a]));

julia> names(chn)
2-element Vector{Symbol}:
 :a
 :b

julia> chn2 = Chains(chn, :internals);

julia> names(chn2)
1-element Vector{Symbol}:
 :a
```
"""
Chains(c::Chains, section::Union{Symbol,String}) = Chains(c, (section,))
function Chains(chn::Chains, sections)
    # Make sure the sections exist first.
	  all(haskey(chn.name_map, Symbol(x)) for x in sections) ||
		    error("some sections are not present in the chain")

	  # Create the new section map.
    name_map = (; (Symbol(k) => chn.name_map[Symbol(k)] for k in sections)...)

    # Extract wanted values.
    value = chn.value[:, reduce(vcat, name_map), :]

    # Create the new chain.
    return Chains(value, chn.logevidence, name_map, chn.info)
end
Chains(chain::Chains, ::Nothing) = chain

# Groups of parameters

"""
    namesingroup(chains::Chains, sym::Union{String,Symbol})

Return the names of all parameters in a chain that belong to the group `sym`.

This is based on the MCMCChains convention that parameters with names of the form `:sym[index]`
belong to one group of parameters called `:sym`.

If the chain contains a parameter of name `:sym` it will be returned as well.

# Example
```jldoctest
julia> chn = Chains(rand(100, 2, 2), ["A[1]", "A[2]"]);

julia> namesingroup(chn, :A)
2-element Vector{Symbol}:
 Symbol("A[1]")
 Symbol("A[2]")
```
"""
namesingroup(chains::Chains, sym::String) = namesingroup(chains, Symbol(sym))
function namesingroup(chains::Chains, sym::Symbol)
    # Start by looking up the symbols in the list of parameter names.
    names_of_params = names(chains)
    regex = Regex("^$sym\$|^$sym\\[")
    indices = findall(x -> match(regex, string(x)) !== nothing, names(chains))
    return names_of_params[indices]
end

"""
    group(chains::Chains, name::Union{String,Symbol})

Return a subset of the chain chain with all parameters in the group `Symbol(name)`.
"""
function group(chains::Chains, name::Union{String,Symbol})
    return chains[:, namesingroup(chains, name), :]
end

#################### Indexing ####################

Base.getindex(c::Chains, i::Integer) = c[i, :, :]
Base.getindex(c::Chains, i::AbstractVector{<:Integer}) = c[i, :, :]

Base.getindex(c::Chains, v::String) = c[:, Symbol(v), :]
Base.getindex(c::Chains, v::AbstractVector{String}) = c[:, Symbol.(v), :]

Base.getindex(c::Chains, v::Symbol) = c[:, v, :]
Base.getindex(c::Chains, v::AbstractVector{Symbol}) = c[:, v, :]

Base.getindex(chn::Chains, i, j, k) = _getindex(chn, chn.value[_toindex(i, j, k)...])

_getindex(::Chains, data) = data
function _getindex(chains::Chains, data::AxisArray{<:Any,3})
    names = data.axes[2].val
    namemap = namemap_intersect(chains.name_map, names)
    return Chains(data, chains.logevidence, namemap, chains.info)
end

# convert strings to symbols but try to keep all dimensions for multiple parameters
_toindex(i, j, k) = (i, string2symbol(j), k)
_toindex(i::Integer, j, k) = (i:i, string2symbol(j), k)
_toindex(i, j, k::Integer) = (i, string2symbol(j), k:k)
_toindex(i::Integer, j, k::Integer) = (i:i, string2symbol(j), k:k)

# return an array or a number if a single parameter is specified
const SingleIndex = Union{Symbol,String,Integer}
_toindex(i, j::SingleIndex, k) = (i, string2symbol(j), k)
_toindex(i::Integer, j::SingleIndex, k) = (i, string2symbol(j), k)
_toindex(i, j::SingleIndex, k::Integer) = (i, string2symbol(j), k)
_toindex(i::Integer, j::SingleIndex, k::Integer) = (i, string2symbol(j), k)

Base.setindex!(c::Chains, v, i...) = setindex!(c.value, v, i...)
Base.lastindex(c::Chains) = lastindex(c.value, 1)
Base.lastindex(c::Chains, d::Integer) = lastindex(c.value, d)

"""
    Base.get(c::Chains, v::Symbol; flatten=false)
    Base.get(c::Chains, vs::Vector{Symbol}; flatten=false)

Return a `NamedTuple` with `v` as the key, and matching paramter
names as the values.

Passing `flatten=true` will return a `NamedTuple` with keys ungrouped.

# Example

```jldoctest
julia> chn = Chains([1:2 3:4]);

julia> get(chn, :param_1)
(param_1 = [1; 2],)

julia> get(chn, [:param_2])
(param_2 = [3; 4],)

julia> get(chn, :param_1; flatten=true)
(param_1 = 1,)
```
"""
Base.get(c::Chains, v::Symbol; flatten = false) = get(c, [v], flatten=flatten)
function Base.get(c::Chains, vs::Vector{Symbol}; flatten = false)
    pairs = Dict()
    for v in vs
        syms = namesingroup(c, v)
        len = length(syms)
        val = ()
        if len > 1
            val = ntuple(i -> c.value[:,syms[i],:], length(syms))
        elseif len == 1
            val = c.value[:,syms[1],:]
        else
            continue
        end

        if flatten
            for i in eachindex(syms)
                pairs[syms[i]] = val[i]
            end
        else
            pairs[v] = val
        end
    end
    return _dict2namedtuple(pairs)
end

"""
    get(c::Chains; section::Union{Vector{Symbol}, Symbol; flatten=false}

Return all parameters in a given section(s) as a `NamedTuple`.

Passing `flatten=true` will return a `NamedTuple` with keys ungrouped.

# Example

```jldoctest
julia> chn = Chains([1:2 3:4], [:a, :b], Dict(:internals => [:a]));

julia> get(chn; section=:parameters)
(b = [3; 4],)

julia> get(chn; section=[:internals])
(a = [1; 2],)
```
"""
function Base.get(
    c::Chains;
    section::Union{Symbol,AbstractVector{Symbol}},
    flatten = false
)
    names = Set(Symbol[])
    regex = r"[^\[]*"
    _section = section isa Symbol ? (section,) : section
    for v in _section
        v in keys(c.name_map) || error("section $v does not exist")

        # If the name contains a bracket,
        # split it so get can group them correctly.
        if flatten
            append!(names, c.name_map[v])
        else
            for name in c.name_map[v]
                m = match(regex, string(name))
                push!(names, Symbol(m.match))
            end
        end
    end

    return get(c, collect(names); flatten = flatten)
end

"""
    get_params(c::Chains; flatten=false)

Return all parameters packaged as a `NamedTuple`. Variables with a bracket
in their name (as in "P[1]") will be grouped into the returned value as P.

Passing `flatten=true` will return a `NamedTuple` with keys ungrouped.

# Example

```jldoctest
julia> chn = Chains(1:5);

julia> x = get_params(chn);

julia> x.param_1
2-dimensional AxisArray{Int64,2,...} with axes:
    :iter, 1:1:5
    :chain, 1:1
And data, a 5×1 Matrix{Int64}:
 1
 2
 3
 4
 5
```
"""
get_params(c::Chains; flatten = false) = get(c, section = sections(c), flatten=flatten)

#################### Base Methods ####################

function Base.show(io::IO, chains::Chains)
    print(io, "MCMC chain (", summary(chains.value.data), ")")
end

function Base.show(io::IO, mime::MIME"text/plain", chains::Chains)
    print(io, "Chains ", chains, ":\n\n", header(chains))

    # Show summary stats.
    summaries = describe(chains)
    for summary in summaries
        println(io)
        show(io, mime, summary)
    end
end

Base.keys(c::Chains) = names(c)
Base.size(c::Chains) = size(c.value)
Base.size(c::Chains, ind) = size(c)[ind]
Base.length(c::Chains) = length(range(c))
Base.first(c::Chains) = first(c.value[Axis{:iter}].val)
Base.step(c::Chains) = step(c.value[Axis{:iter}].val)
Base.last(c::Chains) = last(c.value[Axis{:iter}].val)

Base.convert(::Type{Array}, chn::Chains) = convert(Array, chn.value)

#################### Auxilliary Functions ####################

"""
    range(chains::Chains)

Return the range of iteration indices of the `chains`.
"""
Base.range(chains::Chains) = chains.value[Axis{:iter}].val

"""
    setrange(chains::Chains, range)

Generate a new chain from `chains` with iterations indexed by `range`.

The new chain and `chains` share the same data in memory.
"""
function setrange(chains::Chains, range::AbstractRange{<:Integer})
    if length(chains) != length(range)
        error("length of `range` (", length(range),
              ") is not equal to the number of iterations (", length(chains), ")")
    end

    value = AxisArray(chains.value.data;
                      iter = range, var = names(chains), chain = MCMCChains.chains(chains))

    return Chains(value, chains.logevidence, chains.name_map, chains.info)
end

"""
    resetrange(chains::Chains)

Generate a new chain from `chains` with iterations indexed by `1:n`, where `n` is the number
of samples per chain.

The new chain and `chains` share the same data in memory.
"""
resetrange(chains::Chains) = setrange(chains, 1:size(chains, 1))

"""
    chains(c::Chains)

Return the names or symbols of each chain in a `Chains` object.
"""
function chains(c::Chains)
    return c.value[Axis{:chain}].val
end

"""
    names(chains::Chains)

Return the parameter names in the `chains`.
"""
Base.names(chains::Chains) = chains.value[Axis{:var}].val

"""
    names(chains::Chains, section::Symbol)

Return the parameter names of a `section` in the `chains`.
"""
Base.names(chains::Chains, section::Symbol) = convert(Vector{Symbol}, chains.name_map[section])

"""
    names(chains::Chains, sections)

Return the parameter names of the `sections` in the `chains`.
"""
function Base.names(c::Chains, sections)
    names = Symbol[]
    for section in sections
        append!(names, c.name_map[section])
    end
    return names
end

"""
    get_sections(chains[, sections])

Return multiple `Chains` objects, each containing only a single section.
"""
function get_sections(chains::Chains, sections = keys(chains.name_map))
    return [Chains(chains, section) for section in sections]
end
get_sections(chains::Chains, section::Union{Symbol, String}) = Chains(chains, section)

"""
    sections(c::Chains)

Retrieve a list of the sections in a chain.
"""
sections(c::Chains) = collect(keys(c.name_map))

"""
    header(c::Chains; section=missing)

Return a string containing summary information for a `Chains` object.
If the `section` keyword is used, this function prints only the relevant section
header.

# Example
```julia
# Printing the whole header.
header(chn)

# Print only one section's header.
header(chn, section = :parameter)
```
"""
function header(c::Chains; section=missing)
    rng = range(c)

    # Function to make section strings.
    section_str(sec, arr) = string(
        "$sec",
        repeat(" ", 18 - length(string(sec))),
        "= $(join(map(string, arr), ", "))\n"
    )

    # Set up string array.
    section_strings = String[]

    # Get section lines.
    if section isa Missing
        for (sec, nms) in pairs(c.name_map)
            section_string = section_str(sec, nms)
            push!(section_strings, section_string)
        end
    else
        section in keys(c.name_map) ||
            throw(ArgumentError("$section not found in name map."))
        section_string = section_str(section, c.name_map[section])
        push!(section_strings, section_string)
    end

    #

    # Return header.
    return string(
        ismissing(c.logevidence) ? "" : "Log evidence      = $(c.logevidence)\n",
        "Iterations        = $(first(c)):$(last(c))\n",
        "Thinning interval = $(step(c))\n",
        "Chains            = $(join(map(string, chains(c)), ", "))\n",
        "Samples per chain = $(length(range(c)))\n",
        section_strings...
    )
end

function indiscretesupport(c::Chains,
                           bounds::Tuple{Real, Real}=(0, Inf))
  nrows, nvars, nchains = size(c.value)
  result = Array{Bool}(undef, nvars * (nrows > 0))
  for i in 1:nvars
    result[i] = true
    for j in 1:nrows, k in 1:nchains
      x = c.value[j, i, k]
      if !isinteger(x) || x < bounds[1] || x > bounds[2]
        result[i] = false
        break
      end
    end
  end
  result
end

function link(c::Chains)
  cc = copy(c.value.data)
  for j in axes(cc, 2)
    x = cc[:, j, :]
    if minimum(x) > 0.0
      cc[:, j, :] = maximum(x) < 1.0 ? logit.(x) : log.(x)
    end
  end
  cc
end

### Chains specific functions ###

"""
    sort(c::Chains[; lt=NaturalSort.natural])

Return a new column-sorted version of `c`.

By default the columns are sorted in natural sort order.
"""
function Base.sort(c::Chains; lt = NaturalSort.natural)
    v = c.value
    x, y, z = size(v)
    unsorted = collect(zip(1:y, v.axes[2].val))
    sorted = sort(unsorted, by = x -> string(x[2]), lt=lt)
    new_axes = (v.axes[1], Axis{:var}([n for (_, n) in sorted]), v.axes[3])
    new_v = copy(v.data)
    for i in eachindex(sorted)
        new_v[:, i, :] = v[:, sorted[i][1], :]
    end

    aa = AxisArray(new_v, new_axes...)

    # Sort the name map too:
    namemap = deepcopy(c.name_map)
    for names in namemap
        sort!(names, by=string, lt=lt)
    end

    return Chains(aa, c.logevidence, namemap, c.info)
end

"""
    setinfo(c::Chains, n::NamedTuple)

Return a new `Chains` object with a `NamedTuple` type `n` placed in the `info` field.

# Example
```julia
new_chn = setinfo(chn, NamedTuple{(:a, :b)}((1, 2)))
```
"""
function setinfo(c::Chains, n::NamedTuple)
    return Chains(c.value, c.logevidence, c.name_map, n)
end

"""
    set_section(chains::Chains, namemap)

Create a new `Chains` object from `chains` with the provided `namemap` mapping of parameter
names.

Both chains share the same underlying data. Any parameters in the chain that are unassigned
will be placed into the `:parameters` section.
"""
function set_section(chains::Chains, namemap)
    # Initialize the name map.
    _namemap = initnamemap(namemap)

    # Make sure all the names are in the new name map.
    newnames = Set(Symbol[])
    names_of_params = names(chains)
    for names in _namemap
        filter!(x -> x ∈ names_of_params, names)
        for name in names
            push!(newnames, name)
        end
    end
    missingnames = setdiff(names_of_params, newnames)

    # Assign everything that is missing to :parameters.
    if !isempty(missingnames)
        @warn "Section mapping does not contain all parameter names, " *
            "$missingnames assigned to :parameters."
        for name in missingnames
            push!(_namemap.parameters, name)
        end
    end

    return Chains(chains.value, chains.logevidence, _namemap, chains.info)
end

_default_sections(c::Chains) = haskey(c.name_map, :parameters) ? :parameters : nothing

function _clean_sections(chains::Chains, sections)
    return filter(sections) do section
        haskey(chains.name_map, Symbol(section))
    end
end
function _clean_sections(chains::Chains, section::Union{String,Symbol})
    return haskey(chains.name_map, Symbol(section)) ? section : ()
end
_clean_sections(::Chains, ::Nothing) = nothing

#################### Concatenation ####################

Base.cat(c::Chains, cs::Chains...; dims = Val(1)) = _cat(dims, c, cs...)
Base.cat(c::T, cs::T...; dims = Val(1)) where T<:Chains = _cat(dims, c, cs...)

Base.vcat(c::Chains, cs::Chains...) = _cat(Val(1), c, cs...)
Base.vcat(c::T, cs::T...) where T<:Chains = _cat(Val(1), c, cs...)

Base.hcat(c::Chains, cs::Chains...) = _cat(Val(2), c, cs...)
Base.hcat(c::T, cs::T...) where T<:Chains = _cat(Val(2), c, cs...)

AbstractMCMC.chainscat(c::Chains, cs::Chains...) = _cat(Val(3), c, cs...)

_cat(dim::Int, cs::Chains...) = _cat(Val(dim), cs...)

function _cat(::Val{1}, c1::Chains, args::Chains...)
    # check inputs
    thin = step(c1)
    all(c -> step(c) == thin, args) || throw(ArgumentError("chain thinning differs"))
    nms = names(c1)
    all(c -> names(c) == nms, args) || throw(ArgumentError("chain names differ"))
    chns = chains(c1)
    all(c -> chains(c) == chns, args) || throw(ArgumentError("sets of chains differ"))

    # concatenate all chains
    data = mapreduce(c -> c.value.data, vcat, args; init = c1.value.data)
    value = AxisArray(data;
                      iter = range(first(c1); length = size(data, 1), step = thin),
                      var = nms,
                      chain = chns)

    return Chains(value, missing, c1.name_map, c1.info)
end

function _cat(::Val{2}, c1::Chains, args::Chains...)
    # check inputs
    rng = range(c1)
    all(c -> range(c) == rng, args) || throw(ArgumentError("chain ranges differ"))
    chns = chains(c1)
    all(c -> chains(c) == chns, args) || throw(ArgumentError("sets of chains differ"))

    # combine names and sections of parameters
    nms = names(c1)
    n = length(nms)
    for c in args
        nms = union(nms, names(c))
        n += length(names(c))
        n == length(nms) || throw(ArgumentError("non-unique parameter names"))
    end

    name_map = mapreduce(c -> c.name_map, merge_union, args; init = c1.name_map)

    # concatenate all chains
    data = mapreduce(c -> c.value.data, hcat, args; init = c1.value.data)
    value = AxisArray(data; iter = rng, var = nms, chain = chns)

    return Chains(value, missing, name_map, c1.info)
end

function _cat(::Val{3}, c1::Chains, args::Chains...)
    # check inputs
    rng = range(c1)
    all(c -> range(c) == rng, args) || throw(ArgumentError("chain ranges differ"))
    nms = names(c1)
    all(c -> names(c) == nms, args) || throw(ArgumentError("chain names differ"))

    # concatenate all chains
    data = mapreduce(c -> c.value.data, (x, y) -> cat(x, y; dims = 3), args;
                     init = c1.value.data)
    value = AxisArray(data; iter = rng, var = nms, chain = 1:size(data, 3))

    return Chains(value, missing, c1.name_map, c1.info)
end

function pool_chain(c::Chains)
    data = c.value.data
    pool_data = reshape(permutedims(data, [1, 3, 2]), :, size(data, 2), 1)
    return Chains(pool_data, names(c), c.name_map; info=c.info)
end

"""
    replacenames(chains::Chains, dict::AbstractDict)

Replace parameter names by creating a new `Chains` object that shares the same underlying data.

# Examples
```jldoctest
julia> chn = Chains(rand(100, 2, 2), ["one", "two"]);

julia> chn2 = replacenames(chn, "one" => "A");

julia> names(chn2)
2-element Vector{Symbol}:
 :A
 :two

julia> chn3 = replacenames(chn2, Dict("A" => "one", "two" => "B"));

julia> names(chn3)
2-element Vector{Symbol}:
 :one
 :B
```
"""
replacenames(chains::Chains, dict::AbstractDict) = replacenames(chains, pairs(dict)...)
function replacenames(chains::Chains, old_new::Pair...)
    isempty(old_new) && error("you have to specify at least one replacement")

    # Set new parameter names and a new name map.
    names_of_params = copy(names(chains))
    namemap = deepcopy(chains.name_map)
    for (old, new) in old_new
        symold_symnew = Symbol(old) => Symbol(new)

        replace!(names_of_params, symold_symnew)
        for names in namemap
            replace!(names, symold_symnew)
        end
    end

    value = AxisArray(
        chains.value.data;
        iter = range(chains), var = names_of_params, chain = 1:size(chains, 3)
    )

    return Chains(value, chains.logevidence, namemap, chains.info)
end


text:
"""
sigmoid activation function
"""
sigmoid(z::Number) = 1 / (1 + e^(-z))
sigmoid{T<:Number}(arr::Array{T}) = map(sigmoid, arr)

partial_sigmoid(z::Number) = z * (1 - z)

"""
tanh activation function
"""
tanh = Base.tanh

partial_tanh(z::Number) = 2 / (e^z + e^-z)

"""
ReLU activation function
"""
relu(x::Real) = max(0, x)
relu{T<:Real}(arr::Array{T}) = map(relu, arr)

max_stars_repo_name: zxteloiv/julia-ex

max_issues_repo_name: zxteloiv/julia-ex

max_forks_repo_name: zxteloiv/julia-ex


text:
"""
```julia
simplify(x; expand=false,
            threaded=false,
            thread_subtree_cutoff=100,
            nonzero_denominators=true,
            rewriter=nothing)
```

Simplify an expression (`x`) by applying `rewriter` until there are no changes.
`expand=true` applies [`expand`](/api/#expand) in the beginning of each fixpoint iteration.

By default, simplify will assume denominators are not zero and allow cancellation in fractions.
Pass `simplify_fractions=false` to prevent this.
"""
function simplify(x;
                  expand=false,
                  polynorm=nothing,
                  threaded=false,
                  simplify_fractions=true,
                  thread_subtree_cutoff=100,
                  rewriter=nothing)
    if polynorm !== nothing
        Base.depwarn("simplify(..; polynorm=$polynorm) is deprecated, use simplify(..; expand=$polynorm) instead",
                        :simplify)
    end

    f = if rewriter === nothing
        if threaded
            threaded_simplifier(thread_subtree_cutoff)
        elseif expand
            serial_expand_simplifier
        else
            serial_simplifier
        end
    else
        Fixpoint(rewriter)
    end

    x = PassThrough(f)(x)
    simplify_fractions && has_operation(x, /) ?
        SymbolicUtils.simplify_fractions(x) : x
end

has_operation(x, op) = (istree(x) && (operation(x) == op ||
                                      any(a->has_operation(a, op),
                                          unsorted_arguments(x))))

Base.@deprecate simplify(x, ctx; kwargs...)  simplify(x; rewriter=ctx, kwargs...)


text:
module CxxWrap

module CxxWrapCore

import Libdl
import MacroTools

if ccall(:jl_is_debugbuild, Cint, ()) != 0 && get(ENV, "JLCXX_DIR", "") == ""
  @warn "Julia debug-build may be incompatible with the libcxxwrap-julia binaries, see https://github.com/JuliaInterop/libcxxwrap-julia#building-libcxxwrap-julia for instructions on building compatible binaries."
end

export @wrapmodule, @readmodule, @wraptypes, @wrapfunctions, @safe_cfunction, @initcxx,
ConstCxxPtr, ConstCxxRef, CxxRef, CxxPtr,
CppEnum, ConstArray, CxxBool, CxxLong, CxxULong,
ptrunion, gcprotect, gcunprotect, isnull

const libcxxwrap_version_range = (v"0.6.6",  v"0.7")

# Convert path if it contains lib prefix on windows
function lib_path(so_path::AbstractString)
  path_copy = so_path
  @static if Sys.iswindows()
    basedir, libname = splitdir(so_path)
    libdir_suffix = Sys.WORD_SIZE == 32 ? "32" : ""
    if startswith(libname, "lib") && !isfile(so_path)
      path_copy = joinpath(basedir*libdir_suffix, libname[4:end])
    end
  end
  return path_copy
end

const depsfile = joinpath(dirname(dirname(@__FILE__)), "deps", "deps.jl")
if !isfile(depsfile)
  error("$depsfile not found, CxxWrap did not build properly")
end
include(depsfile)
const jlcxx_path = libcxxwrap_julia

prefix_path() = dirname(dirname(jlcxx_path))

# Internal function names
const cxxwrap_cfnames = [
  :cxxwrap_version_string,
  :has_cxx_module,
  :initialize_cxxwrap,
  :register_julia_module,
  :get_module_functions,
  :bind_module_constants,
  :get_box_types
]

for fname in cxxwrap_cfnames
  @eval const $(Symbol(fname,:_p)) = Ref(C_NULL)
end

function load_cxxwrap_lib()
  loaded_lib = Libdl.dllist()[findall(x -> occursin("libcxxwrap_julia.",x), Libdl.dllist())]
  if length(loaded_lib) == 1
    return Libdl.dlopen(loaded_lib[1], Libdl.RTLD_GLOBAL)
  end

  return Libdl.dlopen(jlcxx_path, Libdl.RTLD_GLOBAL)
end

function load_cxxwrap_symbols()
  if cxxwrap_version_string_p[] != C_NULL
    return
  end

  libcxxwrap = load_cxxwrap_lib()
  for fname in cxxwrap_cfnames
    getproperty(@__MODULE__, Symbol(fname,:_p))[] = Libdl.dlsym(libcxxwrap, string(fname))
  end
end

function checkversion()
  jlcxxversion = VersionNumber(unsafe_string(ccall(cxxwrap_version_string_p[], Cstring, ())))
  if !(libcxxwrap_version_range[1] <= jlcxxversion < libcxxwrap_version_range[2])
    error("This version of CxxWrap requires a libcxxwrap-julia in the range $(libcxxwrap_version_range), but version $jlcxxversion was found")
  end
end

# Must also be called during precompile
load_cxxwrap_symbols()
checkversion()

# Welcome to the C/C++ integer type mess

abstract type CxxSigned <: Signed end
abstract type CxxUnsigned <: Unsigned end

# long is a special case, because depending on the platform it overlaps with int or long long. See https://en.cppreference.com/w/cpp/language/types
primitive type CxxLong <: CxxSigned 8*sizeof(Clong) end
primitive type CxxULong <: CxxUnsigned 8*sizeof(Culong) end
primitive type CxxBool <: CxxUnsigned 8*sizeof(Cuchar) end
const CharSigning = supertype(Cchar) == Signed ? CxxSigned : CxxUnsigned
primitive type CxxChar <: CharSigning 8*sizeof(Cchar) end
primitive type CxxUChar <: CxxUnsigned 8*sizeof(Cuchar) end
const WCharSigning = supertype(Cwchar_t) == Signed ? CxxSigned : CxxUnsigned
primitive type CxxWchar <: WCharSigning 8*sizeof(Cwchar_t) end

# This macro adds the fixed integer types described on https://en.cppreference.com/w/cpp/types/integer
# Names are e.g. CxxInt32 or CxxUInt64
macro add_int_types()
  result = quote end
  for signed in (Symbol(), :U)
    super = signed == :U ? CxxUnsigned : CxxSigned
    for nbits in (8,16,32,64)
      push!(result.args, esc(:(primitive type $(Symbol(:Cxx, signed, :Int, nbits)) <: $super $nbits end)))
    end
  end
  return result
end
@add_int_types

# Get the equivalen Julia type for a Cxx integer type
@generated julia_int_type(::Type{T}) where {T<:CxxSigned} = Symbol(:Int, 8*sizeof(T))
@generated julia_int_type(::Type{T}) where {T<:CxxUnsigned} = Symbol(:UInt, 8*sizeof(T))

to_julia_int(x::Union{CxxSigned,CxxUnsigned}) = reinterpret(julia_int_type(typeof(x)),x)

# Conversion to and from the equivalent Julia type
Base.convert(::Type{T}, x::Number) where {T<:Union{CxxSigned,CxxUnsigned}} = reinterpret(T, convert(julia_int_type(T), x))
Base.convert(::Type{JT}, x::CT) where {JT<:Number,CT<:Union{CxxSigned,CxxUnsigned}} = convert(JT,reinterpret(julia_int_type(CT), x))
Base.convert(::Type{T}, x::CT) where {T <: Union{CxxWrapCore.CxxSigned, CxxWrapCore.CxxUnsigned}, CT <: Union{CxxWrapCore.CxxSigned, CxxWrapCore.CxxUnsigned}}  = convert(T,reinterpret(julia_int_type(CT), x))

# Convenience constructors
(::Type{T})(x) where {T<:Union{CxxWrapCore.CxxSigned,CxxWrapCore.CxxUnsigned}} = convert(T,x)

Base.flipsign(x::CxxLong, y::CxxLong) = reinterpret(CxxLong, flipsign(to_julia_int(x), to_julia_int(y)))

# Trait type to indicate a type is a C++-wrapped type
struct IsCxxType end
struct IsNormalType end

@inline cpp_trait_type(::Type) = IsNormalType

# Enum type interface
abstract type CppEnum end
Base.convert(::Type{Int32}, x::CppEnum) = reinterpret(Int32, x)
import Base: +, |
+(a::T, b::T) where {T <: CppEnum} = reinterpret(T, convert(Int32,a) + convert(Int32,b))
|(a::T, b::T) where {T <: CppEnum} = reinterpret(T, convert(Int32,a) | convert(Int32,b))

"""
Base class for smart pointers
"""
abstract type SmartPointer{T} end
@inline cpp_trait_type(::Type{SmartPointer{T}}) where {T} = IsCxxType

Base.show(io::IO, x::SmartPointer) = print(io, "C++ smart pointer of type ", typeof(x))

allocated_type(t::Type) = Any
dereferenced_type(t::Type) = Any

__cxxwrap_smartptr_dereference(p::SmartPointer{T}) where {T} = __cxxwrap_smartptr_dereference(CxxRef(p))
__cxxwrap_smartptr_construct_from_other(t::Type{<:SmartPointer{T}}, p::SmartPointer{T}) where {T} = __cxxwrap_smartptr_construct_from_other(t,CxxRef(p))
__cxxwrap_smartptr_cast_to_base(p::SmartPointer{T}) where {T} = __cxxwrap_smartptr_cast_to_base(CxxRef(p))

function Base.getindex(p::SmartPointer{T}) where {T}
  return __cxxwrap_smartptr_dereference(p)
end

# No conversion if source and target type are identical
Base.convert(::Type{T}, p::T) where {PT,T <: SmartPointer{PT}} = p

# Construct from a related pointer, e.g. a std::weak_ptr from std::shared_ptr
function Base.convert(::Type{T1}, p::T2) where {T, T1 <: SmartPointer{T}, T2 <: SmartPointer{T}}
  return __cxxwrap_smartptr_construct_from_other(T1, p)
end

# upcast to base class
function Base.convert(::Type{T1}, p::T2) where {BaseT,DerivedT, T1 <: SmartPointer{BaseT}, T2 <: SmartPointer{DerivedT}}
  if !(DerivedT <: BaseT)
    error("$DerivedT does not inherit from $BaseT in smart pointer convert")
  end
  # First convert to base type
  base_p = __cxxwrap_smartptr_cast_to_base(p)
  return convert(T1, base_p)
end

struct StrictlyTypedNumber{NumberT}
  value::NumberT
end
function Base.convert(t::Type{<:StrictlyTypedNumber}, n::Number)
  @assert t == StrictlyTypedNumber{typeof(n)}
  return StrictlyTypedNumber{typeof(n)}(n)
end

abstract type CxxBaseRef{T} <: Ref{T} end

struct CxxPtr{T} <: CxxBaseRef{T}
  cpp_object::Ptr{T}
  CxxPtr{T}(x::Ptr) where {T} = new{T}(x)
  CxxPtr{T}(x::CxxBaseRef) where {T} = new{T}(x.cpp_object)
end

struct ConstCxxPtr{T} <: CxxBaseRef{T}
  cpp_object::Ptr{T}
  ConstCxxPtr{T}(x::Ptr) where {T} = new{T}(x)
  ConstCxxPtr{T}(x::CxxBaseRef) where {T} = new{T}(x.cpp_object)
end

struct CxxRef{T} <: CxxBaseRef{T}
  cpp_object::Ptr{T}
  CxxRef{T}(x::Ptr) where {T} = new{T}(x)
  CxxRef{T}(x::CxxBaseRef) where {T} = new{T}(x.cpp_object)
end

struct ConstCxxRef{T} <: CxxBaseRef{T}
  cpp_object::Ptr{T}
  ConstCxxRef{T}(x::Ptr) where {T} = new{T}(x)
  ConstCxxRef{T}(x::CxxBaseRef) where {T} = new{T}(x.cpp_object)
end

_ref_type(::Type{RefT}, ::Type{<:CxxBaseRef{T}}) where {RefT,T} = _ref_type(RefT,T)
_ref_type(::Type{RefT}, ::Type{T}) where {RefT,T} = _ref_type(RefT, T, cpp_trait_type(T))
_ref_type(::Type{RefT}, ::Type{T}, ::Type) where {RefT,T} = error("Manual Cxx Reference creation is only for C++ types")
function _ref_type(::Type{RefT}, ::Type{T}, ::Type{IsCxxType}) where {RefT,T}
  if isconcretetype(T)
    return RefT{supertype(T)}
  end
  return RefT{T}
end
_make_ref(::Type{RefT}, x::T) where {RefT,T} = _ref_type(RefT,T)(x.cpp_object)

CxxPtr(x) = _make_ref(CxxPtr,x)
ConstCxxPtr(x) = _make_ref(ConstCxxPtr,x)
CxxRef(x) = _make_ref(CxxRef,x)
ConstCxxRef(x) = _make_ref(ConstCxxRef,x)

Base.convert(t::Type{<:CxxBaseRef{T}}, x::Ptr{NT}) where {T <: Union{CxxUnsigned, CxxSigned}, NT <: Integer} = t(reinterpret(Ptr{T}, x))

Base.:(==)(a::Union{CxxPtr,ConstCxxPtr}, b::Union{CxxPtr,ConstCxxPtr}) = (a.cpp_object == b.cpp_object)
Base.:(==)(a::CxxBaseRef, b::Ptr) = (a.cpp_object == b)
Base.:(==)(a::Ptr, b::CxxBaseRef) = (b == a)
Base.:(==)(a::Union{CxxRef,ConstCxxRef}, b) = (a[] == b)
Base.:(==)(a, b::Union{CxxRef,ConstCxxRef}) = (b == a)
Base.:(==)(a::Union{CxxRef,ConstCxxRef}, b::Union{CxxRef,ConstCxxRef}) = (a[] == b[])

_deref(p::CxxBaseRef, ::Type) = unsafe_load(p.cpp_object)
_deref(p::CxxBaseRef{T}, ::Type{IsCxxType}) where {T} = dereferenced_type(T)(p.cpp_object)

_store_to_cxxptr(::Any,::Any,::Type) = error("Resetting the value to a C++ pointer or reference is only supported for non-wrapped types")
_store_to_cxxptr(r::Union{ConstCxxPtr{T},ConstCxxRef{T}}, x::T, ::Type{IsNormalType}) where {T} = error("Setting the value of a const reference or pointer is not allowed")
_store_to_cxxptr(r::Union{CxxPtr{T},CxxRef{T}}, x::T, ::Type{IsNormalType}) where {T} = unsafe_store!(r.cpp_object, x)

Base.unsafe_load(p::CxxBaseRef{T}) where {T} = _deref(p, cpp_trait_type(T))
Base.unsafe_string(p::CxxBaseRef) = unsafe_string(p.cpp_object)
Base.getindex(r::CxxBaseRef) = unsafe_load(r)
Base.setindex!(r::CxxBaseRef{T}, x::T) where {T}  = _store_to_cxxptr(r, x, cpp_trait_type(T))

# Delegate iteration to the contained type
Base.iterate(x::CxxBaseRef) = Base.iterate(x[])
Base.iterate(x::CxxBaseRef, state) = Base.iterate(x[], state)
Base.length(x::CxxBaseRef) = Base.length(x[])
Base.size(x::CxxBaseRef, d) = Base.size(x[], d)

# Delegate broadcast operations to the contained type
Base.BroadcastStyle(::Type{<:CxxBaseRef{T}}) where {T} = Base.BroadcastStyle(T)
Base.axes(x::CxxBaseRef) = Base.axes(x[])
Base.broadcastable(x::CxxBaseRef) = Base.broadcastable(x[])

# Delegate indexing to the wrapped type
Base.getindex(x::CxxBaseRef, i::Int) = Base.getindex(x[], i)
Base.setindex!(x::CxxBaseRef, val, i::Int) = Base.setindex!(x[], val, i)

Base.convert(::Type{RT}, p::SmartPointer{T}) where {T, RT <: CxxBaseRef{T}} = p[]
Base.cconvert(::Type{RT}, p::SmartPointer{T}) where {T, RT <: CxxBaseRef{T}} = p[]
function Base.convert(::Type{T1}, p::SmartPointer{DerivedT}) where {BaseT,T1 <: BaseT, DerivedT <: BaseT}
  return cxxupcast(T1, p[])[]
end
Base.convert(to_type::Type{Any}, x::CxxWrapCore.SmartPointer{DerivedT}) where {DerivedT} = x
function Base.convert(to_type::Type{<:Ref{T1}}, p::T2) where {BaseT,DerivedT, T1 <: BaseT, T2 <: SmartPointer{DerivedT}}
  return to_type(convert(T1,p))
end

Base.unsafe_convert(to_type::Type{<:CxxBaseRef}, x) = to_type(x.cpp_object)

# This is defined on the C++ side for each wrapped type
cxxupcast(x) = cxxupcast(CxxRef(x))
cxxupcast(x::CxxRef) = error("No upcast for type $(supertype(typeof(x))). Did you specialize SuperType to enable automatic upcasting?")
function cxxupcast(::Type{T}, x) where {T}
  cxxupcast(T, cxxupcast(x))
end
cxxupcast(::Type{T}, x::CxxBaseRef{T}) where {T} = x

struct ConstArray{T,N} <: AbstractArray{T,N}
  ptr::ConstCxxPtr{T}
  size::NTuple{N,Int}
end

ConstArray(ptr::ConstCxxPtr{T}, args::Vararg{Int,N}) where {T,N} = ConstArray{T,N}(ptr, (args...,))

Base.IndexStyle(::ConstArray) = IndexLinear()
Base.size(arr::ConstArray) = arr.size
Base.getindex(arr::ConstArray, i::Integer) = unsafe_load(arr.ptr.cpp_object, i)

function __delete end
function delete(x)
  __delete(CxxPtr(x))
  x.cpp_object = C_NULL
end

# Encapsulate information about a function
mutable struct CppFunctionInfo
  name::Any
  argument_types::Array{Type,1}
  return_type::Type
  julia_return_type::Type
  function_pointer::Int
  thunk_pointer::Int
  override_module::Union{Nothing,Module}
end

const _gc_protected = Dict{UInt64,Tuple{Any, Int}}()

function protect_from_gc(x)
  id = objectid(x)
  (_,n) = get(_gc_protected, id, (x,0))
  _gc_protected[id] = (x,n+1)
  return
end

function unprotect_from_gc(x)
  id = objectid(x)
  (_,n) = get(_gc_protected, id, (x,0))
  if n == 0
    println("warning: attempt to unprotect non-protected object $x")
  end
  if n == 1
    delete!(_gc_protected, id)
  else
    _gc_protected[id] = (x,n-1)
  end
  return
end

function initialize_cxx_lib()
  _c_protect_from_gc = @cfunction protect_from_gc Nothing (Any,)
  _c_unprotect_from_gc = @cfunction unprotect_from_gc Nothing (Any,)
  ccall(initialize_cxxwrap_p[], Cvoid, (Any, Any, Ptr{Cvoid}, Ptr{Cvoid}), @__MODULE__, CppFunctionInfo, _c_protect_from_gc, _c_unprotect_from_gc)
end

# Must also be called during precompile
initialize_cxx_lib()

function __init__()
  load_cxxwrap_symbols()
  checkversion()
  initialize_cxx_lib()
end

function has_cxx_module(mod::Module)
  r = ccall(has_cxx_module_p[], Cuchar, (Any,), mod)
  return r != 0
end

function register_julia_module(mod::Module, fptr::Ptr{Cvoid})
  ccall(register_julia_module_p[], Cvoid, (Any,Ptr{Cvoid}), mod, fptr)
end

function register_julia_module(mod::Module)
  fptr = Libdl.dlsym(Libdl.dlopen(mod.__cxxwrap_sopath, mod.__cxxwrap_flags), mod.__cxxwrap_wrapfunc)
  if !has_cxx_module(mod)
    empty!(mod.__cxxwrap_pointers)
    register_julia_module(mod, fptr)
  end
  if length(mod.__cxxwrap_pointers) != mod.__cxxwrap_nbpointers
    error("Binary part of module was changed since last precompilation, please rebuild.")
  end
end

function get_module_functions(mod::Module)
  ccall(get_module_functions_p[], Any, (Any,), mod)
end

function bind_constants(m::Module, symbols::Array, values::Array)
  ccall(bind_module_constants_p[], Cvoid, (Any,Any,Any), m, symbols, values)
end

function box_types(mod::Module)
  ccall(get_box_types_p[], Any, (Any,), mod)
end

"""
Protect a variable from garbage collection by adding it to the global array kept by CxxWrap
"""
function gcprotect(x)
  protect_from_gc(x)
end

"""
Unprotect a variable from garbage collection by removing it from the global array kept by CxxWrap
"""
function gcunprotect(x)
  unprotect_from_gc(x)
end

# Interpreted as a constructor for Julia  > 0.5
mutable struct ConstructorFname
  _type::DataType
end

# Interpreted as an operator call overload
mutable struct CallOpOverload
  _type::DataType
end

process_fname(fn::Symbol) = fn
process_fname(fn::Tuple{Symbol,Module}) = :($(fn[2]).$(fn[1]))
process_fname(fn::ConstructorFname) = :(::$(Type{fn._type}))
function process_fname(fn::CallOpOverload)
  return :(arg1::$(fn._type))
end

make_func_declaration(fn, argmap) = :($(process_fname(fn))($(argmap...)))
function make_func_declaration(fn::CallOpOverload, argmap)
  return :($(process_fname(fn))($((argmap[2:end])...)))
end

# By default, no argument overloading happens
argument_overloads(t::Type) = Type[]
@static if Int != Cint
  argument_overloads(t::Type{Cint}) = [Int]
end
@static if UInt != Cuint
  argument_overloads(t::Type{Cuint}) = [UInt, Int]
else
  argument_overloads(t::Type{Cuint}) = [Int]
end
argument_overloads(t::Type{Float64}) = [Int, Irrational]
function argument_overloads(t::Type{Array{AbstractString,1}})
  return [Array{String,1}]
end
argument_overloads(t::Type{Ptr{T}}) where {T <: Number} = [Array{T,1}]

"""
Create a Union containing the type and a smart pointer to any type derived from it
"""
function ptrunion(::Type{T}) where {T}
  result{T2 <: T} = Union{T2, SmartPointer{T2}}
  return result
end

smart_pointer_type(t::Type) = smart_pointer_type(cpp_trait_type(t), t)
smart_pointer_type(::Type{IsNormalType}, t::Type) = t
smart_pointer_type(::Type{IsCxxType}, x::Type{T}) where {T} = ptrunion(x)

function smart_pointer_type(::Type{<:SmartPointer{T}}) where {T}
  result{T2 <: T} = SmartPointer{T2}
  return result
end

map_julia_arg_type(t::Type) = Union{Base.invokelatest(smart_pointer_type,t),argument_overloads(t)...}
map_julia_arg_type(a::Type{StrictlyTypedNumber{T}}) where {T} = T
map_julia_arg_type(x::Type{T}) where {T<:Integer} = map_julia_arg_type(x, Base.invokelatest(cpp_trait_type, T))
map_julia_arg_type(x::Type{<:Integer}, ::Type{IsNormalType}) = Integer
map_julia_arg_type(x::Type{<:Integer}, ::Type{IsCxxType}) = x

const PtrTypes{T} = Union{CxxPtr{T}, Array{T}, CxxRef{T}, Base.RefValue{T}, Ptr{T},T}
const ConstPtrTypes{T} = Union{Ref{T}, Array{T}}

map_julia_arg_type(t::Type{<:CxxBaseRef{T}}) where {T} = map_julia_arg_type(t, Base.invokelatest(cpp_trait_type, T))

map_julia_arg_type(t::Type{ConstCxxRef{T}}, ::Type{IsNormalType}) where {T} = Union{Ref{T},map_julia_arg_type(T)}
map_julia_arg_type(t::Type{ConstCxxPtr{T}}, ::Type{IsNormalType}) where {T} = Union{ConstPtrTypes{T},Ptr{Cvoid}}
map_julia_arg_type(t::Type{CxxRef{T}}, ::Type{IsNormalType}) where {T} = PtrTypes{T}
map_julia_arg_type(t::Type{CxxPtr{T}}, ::Type{IsNormalType}) where {T} = Union{PtrTypes{T},Ptr{Cvoid}}

map_julia_arg_type(t::Type{ConstCxxRef{T}}, ::Type{IsNormalType}) where {T<:Union{CxxSigned,CxxUnsigned}} = Union{ConstPtrTypes{julia_int_type(T)},map_julia_arg_type(T)}
map_julia_arg_type(t::Type{ConstCxxPtr{T}}, ::Type{IsNormalType}) where {T<:Union{CxxSigned,CxxUnsigned}} = Union{ConstPtrTypes{julia_int_type(T)},Ptr{Cvoid}}
map_julia_arg_type(t::Type{CxxRef{T}}, ::Type{IsNormalType}) where {T<:Union{CxxSigned,CxxUnsigned}} = PtrTypes{julia_int_type(T)}
map_julia_arg_type(t::Type{CxxPtr{T}}, ::Type{IsNormalType}) where {T<:Union{CxxSigned,CxxUnsigned}} = Union{PtrTypes{julia_int_type(T)},Ptr{Cvoid}}

map_julia_arg_type(t::Type{ConstCxxRef{T}}, ::Type{IsCxxType}) where {T} = Union{map_julia_arg_type(T),ConstCxxRef{<:T},CxxRef{<:T}}
map_julia_arg_type(t::Type{ConstCxxPtr{T}}, ::Type{IsCxxType}) where {T} = Union{CxxPtr{<:T},ConstCxxPtr{<:T}, Ptr{Cvoid}}
map_julia_arg_type(t::Type{CxxRef{T}}, ::Type{IsCxxType}) where {T} = Union{map_julia_arg_type(T),CxxRef{<:T}}
map_julia_arg_type(t::Type{CxxPtr{T}}, ::Type{IsCxxType}) where {T} = Union{CxxPtr{<:T},Ptr{Cvoid}}

map_julia_arg_type(t::Type{CxxPtr{CxxChar}}) = Union{PtrTypes{Cchar}, String}
map_julia_arg_type(t::Type{<:Array{T}}) where {T <: Union{CxxSigned,CxxUnsigned}} = Union{t, Array{julia_int_type(T)}}
map_julia_arg_type(t::Type{<:Array{Ptr{T}}}) where {T <: Union{CxxSigned,CxxUnsigned}} = Union{t, Array{Ptr{julia_int_type(T)}}}
map_julia_arg_type(t::Type{ConstCxxPtr{CxxChar}}) = Union{ConstPtrTypes{Cchar}, String}

# names excluded from julia type mapping
const __excluded_names = Set([
      :cxxupcast,
      :__cxxwrap_smartptr_dereference,
      :__cxxwrap_smartptr_construct_from_other,
      :__cxxwrap_smartptr_cast_to_base
])

function Base.cconvert(to_type::Type{<:CxxBaseRef{T}}, x) where {T}
  return cxxconvert(to_type, x, cpp_trait_type(T))
end
function Base.cconvert(to_type::Type{<:CxxBaseRef{T}}, x::Ptr{PT}) where {T<:Integer,PT<:Integer}
  @assert T == PT || julia_int_type(T) == PT
  return to_type(x)
end
@inline Base.cconvert(::Type{T}, v::T) where {T2,T <: CxxBaseRef{T2}} = v
Base.cconvert(to_type::Type{<:CxxBaseRef{T}}, v::CxxBaseRef{T}) where {T} = to_type(v.cpp_object)
Base.unsafe_convert(to_type::Type{<:CxxBaseRef{T}}, v::CxxBaseRef) where {T} = to_type(v.cpp_object)
Base.unsafe_convert(to_type::Type{<:CxxBaseRef}, v::Base.RefValue) = to_type(pointer_from_objref(v))

cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x, ::Type{IsNormalType}) where {T} = Ref{T}(convert(T,x))
cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x::Base.RefValue{T}, ::Type{IsNormalType}) where {T} = x
cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x::Ptr{T}, ::Type{IsNormalType}) where {T} = to_type(x)
cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x::Ptr{Cvoid}, ::Type{IsNormalType}) where {T} = to_type(x)
cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x::Union{Array,String}, ::Type{IsNormalType}) where {T} = to_type(pointer(x))
function cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x, ::Type{IsCxxType}) where {T}
  return convert(T,x)
end

function cxxconvert(to_type::Type{<:CxxBaseRef{T}}, x::CxxBaseRef, ::Type{IsCxxType}) where {T}
  return to_type(convert(T,x[]))
end

# Build the expression to wrap the given function
function build_function_expression(func::CppFunctionInfo, mod=nothing)
  # Arguments and types
  argtypes = func.argument_types
  argsymbols = map((i) -> Symbol(:arg,i[1]), enumerate(argtypes))

  # These are actually indices into the module-global function pointer table
  fpointer = func.function_pointer
  thunk = func.thunk_pointer

  map_c_arg_type(t::Type) = t
  map_c_arg_type(::Type{Array{T,1}}) where {T <: AbstractString} = Any
  map_c_arg_type(::Type{Type}) = Any
  map_c_arg_type(::Type{T}) where {T <: Tuple} = Any
  map_c_arg_type(::Type{ConstArray{T,N}}) where {T,N} = Any
  map_c_arg_type(::Type{T}) where {T<:Union{CxxSigned,CxxUnsigned}} = julia_int_type(T)

  # Builds the return type passed to ccall
  map_c_return_type(t) = t
  map_c_return_type(::Type{T}) where {T <: Tuple} = Any
  map_c_return_type(::Type{ConstArray{T,N}}) where {T,N} = Any
  map_c_return_type(::Type{T}) where {T<:Union{CxxSigned,CxxUnsigned}} = map_c_arg_type(T)

  # Build the types for the ccall argument list
  c_arg_types = map_c_arg_type.(func.argument_types)
  c_return_type = map_c_return_type(func.return_type)

  # Builds the return-type annotation for the Julia function
  map_julia_return_type(t) = t
  map_julia_return_type(::Type{T}) where {T<:Union{CxxSigned,CxxUnsigned}} = map_c_arg_type(T)
  map_julia_return_type(::Type{CxxBool}) = Bool

  # Build the final call expression
  call_exp = quote end
  if thunk == 0
    push!(call_exp.args, :(ccall(__cxxwrap_pointers[$fpointer], $c_return_type, ($(c_arg_types...),), $(argsymbols...)))) # Direct pointer call
  else
    push!(call_exp.args, :(ccall(__cxxwrap_pointers[$fpointer], $c_return_type, (Ptr{Cvoid}, $(c_arg_types...)), __cxxwrap_pointers[$thunk], $(argsymbols...)))) # use thunk (= std::function)
  end

  function map_julia_arg_type_named(fname, t)
    if fname ∈ __excluded_names
      return t
    end
    return map_julia_arg_type(t)
  end

  # Build an array of arg1::Type1... expressions
  function argmap(signature)
    result = Expr[]
    for (t, s) in zip(signature, argsymbols)
      push!(result, :($s::$(map_julia_arg_type_named(func.name, t))))
    end
    return result
  end

  fname = mod === nothing ? func.name : (func.name,mod)
  function_expression = :($(make_func_declaration(fname, argmap(argtypes)))::$(map_julia_return_type(func.julia_return_type)) = $call_exp)
  return function_expression
end

function makereftype(::Type{T}, mod) where {T}
  basename = T.name.name
  refname = Symbol(basename,"Dereferenced")
  tmod = T.name.module
  params = (T.parameters...,)
  if isempty(params)
    # If the type is non-parametric, this should be hit only once
    @assert !isdefined(mod, refname)
    @assert mod == tmod
    return Core.eval(mod, :(struct $refname <: $T cpp_object::Ptr{Cvoid} end; $refname))
  end
  if !isdefined(tmod, refname)
    #@assert mod == tmod
    parameternames = (Symbol(:T,i) for i in 1:length(params))
    expr = :(struct $refname{$(parameternames...)} <: $basename{$(parameternames...)} cpp_object::Ptr{Cvoid} end)
    Core.eval(mod, expr)
  end
  expr = :($tmod.$refname{$(params...)})
  return Core.eval(mod, expr)
end

function wrap_reference_converters(julia_mod)
  boxtypes = box_types(julia_mod)
  for bt in boxtypes
    st = supertype(bt)
    Core.eval(julia_mod, :(@inline $(@__MODULE__).cpp_trait_type(::Type{<:$st}) = $(@__MODULE__).IsCxxType))
    Core.eval(julia_mod, :(Base.convert(t::Type{$st}, x::T) where {T <: $st} = $(cxxupcast)($st,x)))
    Core.eval(julia_mod, :($(@__MODULE__).allocated_type(::Type{$st}) = $bt))
    reftype = makereftype(st, julia_mod)
    Core.eval(julia_mod, :($(@__MODULE__).dereferenced_type(::Type{$st}) = $reftype))
    Core.eval(julia_mod, :(Base.convert(::Type{$st}, x::$bt) = x))
    Core.eval(julia_mod, :(Base.convert(::Type{$st}, x::$reftype) = x))
  end
end

# Wrap functions from the cpp module to the passed julia module
function wrap_functions(functions, julia_mod)
  basenames = Set([
    :getindex,
    :setindex!,
    :convert,
    :deepcopy_internal,
    :+,
    :*,
    :(==)
  ])

  for func in functions
    if func.override_module != nothing
      Core.eval(julia_mod, build_function_expression(func, func.override_module))
    elseif func.name ∈ basenames
      @warn "Adding method to $(func.name) in the Base module by default. This default will change in the next release of CxxWrap, please use the C++ function `set_override_module` to explicitly set the method module."
      Core.eval(julia_mod, build_function_expression(func, Base))
    else
      Core.eval(julia_mod, build_function_expression(func))
    end
  end
end

# Place the types for the module with the name corresponding to the current module name in the current module
function wraptypes(jlmod)
  symbols = Any[]
  values = Any[]
  bind_constants(jlmod, symbols, values)
  for (sym,val) in zip(symbols, values)
    Core.eval(jlmod, :(const $sym = $val))
  end
end

function wrapfunctions(jlmod)
  module_functions = get_module_functions(jlmod)
  wrap_reference_converters(jlmod)
  wrap_functions(module_functions, jlmod)
end

function readmodule(so_path::AbstractString, funcname, m::Module, flags)
  if flags === nothing
    flags = Libdl.RTLD_LAZY | Libdl.RTLD_DEEPBIND
  end
  Core.eval(m, :(const __cxxwrap_pointers = Ptr{Cvoid}[]))
  Core.eval(m, :(const __cxxwrap_sopath = $so_path))
  Core.eval(m, :(const __cxxwrap_wrapfunc = $(QuoteNode(funcname))))
  Core.eval(m, :(const __cxxwrap_flags = $flags))
  fptr = Libdl.dlsym(Libdl.dlopen(so_path, flags), funcname)
  register_julia_module(m, fptr)
  nb_pointers = length(m.__cxxwrap_pointers)
  Core.eval(m, :(const __cxxwrap_nbpointers = $nb_pointers))
end

function wrapmodule(so_path::AbstractString, funcname, m::Module, flags)
  readmodule(so_path, funcname, m, flags)
  wraptypes(m)
  wrapfunctions(m)
end

"""
  @wrapmodule libraryfile [functionname]

Place the functions and types from the C++ lib into the module enclosing this macro call
Calls an entry point named `define_julia_module`, unless another name is specified as
the second argument.
"""
macro wrapmodule(libraryfile, register_func=:(:define_julia_module), flags=:(nothing))
  return :(wrapmodule($(esc(libraryfile)), $(esc(register_func)), $__module__, $(esc(flags))))
end

"""
  @readmodule libraryfile [functionname]

Read a C++ module and associate it with the Julia module enclosing the macro call.
"""
macro readmodule(libraryfile, register_func=:(:define_julia_module), flags=:(nothing))
  return :(readmodule($(esc(libraryfile)), $(esc(register_func)), $__module__, $(esc(flags))))
end

"""
  @wraptypes

Wrap the types defined in the C++ side of the enclosing module. Requires that
`@readmodule` was called first.
"""
macro wraptypes()
  return :(wraptypes($__module__))
end

"""
  @wrapfunctions

Wrap the functions defined in the C++ side of the enclosing module. Requires that
`@readmodule` and `@wraptypes` was called first.
"""
macro wrapfunctions()
  return :(wrapfunctions($__module__))
end

"""
  @initcxx

Initialize the C++ pointer tables in a precompiled module using CxxWrap. Must be called from within
`__init__` in the wrapped module
"""
macro initcxx()
  return :(register_julia_module($__module__))
end

struct SafeCFunction
  fptr::Ptr{Cvoid}
  return_type::Type
  argtypes::Array{Type,1}
end

macro safe_cfunction(f, rt, args)
  return esc(:($(@__MODULE__).SafeCFunction(@cfunction($f, $rt, $args), $rt, [$(args.args...)])))
end

isnull(x::CxxBaseRef) = (x.cpp_object == C_NULL)

reference_type_union(::Type{T}) where {T} = reference_type_union(T, Base.invokelatest(cpp_trait_type, T))
reference_type_union(::Type{T}, ::Type{IsNormalType}) where {T} = T
reference_type_union(::Type{T}, ::Type{IsCxxType}) where {T} = Union{T, CxxBaseRef{T}, SmartPointer{T}}
reference_type_union(t::TypeVar) = t

dereference_argument(x) = x
dereference_argument(x::CxxBaseRef) = x[]
dereference_argument(x::SmartPointer) = x[]

macro cxxdereference(f)
  fdict = MacroTools.splitdef(f)

  function maparg(a)
    (argname, argtype, slurp, default) = MacroTools.splitarg(a)
    return MacroTools.combinearg(argname, :($(@__MODULE__).reference_type_union($(argtype))), slurp, default)
  end

  # Adapt the signature
  fdict[:args] .= maparg.(fdict[:args])
  fdict[:kwargs] .= maparg.(fdict[:kwargs])

  # Dereference the arguments
  deref_expr = quote end
  for arg in vcat(fdict[:args], fdict[:kwargs])
    (argname, _, slurp, _) = MacroTools.splitarg(arg)
    if argname === nothing
      continue
    end
    if !slurp
      push!(deref_expr.args, :($argname = $(@__MODULE__).dereference_argument($argname)))
    else
      push!(deref_expr.args, :($argname = (($(@__MODULE__).dereference_argument.($argname))...,)))
    end
  end
  insert!(fdict[:body].args, 1, deref_expr)
  fdict[:body] = MacroTools.flatten(fdict[:body])

  # Reassemble the function
  return esc(MacroTools.combinedef(fdict))
end

end

include("StdLib.jl")

using .CxxWrapCore
using .CxxWrapCore: CxxBaseRef, jlcxx_path, argument_overloads, SafeCFunction, reference_type_union, dereference_argument, @cxxdereference

export @wrapmodule, @readmodule, @wraptypes, @wrapfunctions, @safe_cfunction, @initcxx,
ConstCxxPtr, ConstCxxRef, CxxRef, CxxPtr,
CppEnum, ConstArray, CxxBool, CxxLong, CxxULong,
ptrunion, gcprotect, gcunprotect, isnull

using .StdLib: StdVector, StdString, StdWString, StdValArray

export StdVector, StdString, StdWString, StdValArray, @cxxdereference

end # module


text:
module LangrethsRules

using LinearAlgebra
import KadanoffBaym: calculate_weights

export TimeOrderedGF, TimeOrderedConvolution, ⋆
export greater, lesser, advanced, retarded

function skew_hermitify!(x)
    for i in 1:size(x, 1)
        for j in 1:(i - 1)
            x[j,i] = -conj(x[i,j])
        end

        x[i,i] = eltype(x) <: Real ? 0.0 : im * imag(x[i,i])
    end
    x
end

include("common.jl")
include("rules.jl")

end # module


text:
using MathTeXEngine
using CairoMakie
using Colors
using LaTeXStrings

import MathTeXEngine: TeXChar, VLine, HLine, leftinkbound, descender

draw_texelement!(args... ; size=64) = nothing

function draw_texelement!(ax, texchar::TeXChar, position, scale ; size=64)
    x = (position[1] + leftinkbound(texchar)) * size
    y = (position[2] + descender(texchar.font) * scale) * size
    text!(ax, string(texchar.char), font=texchar.font,
        position=Point2f0(x, y),
        textsize=size*scale,
        space=:data)
end

function draw_texelement!(ax, line::VLine, position, scale ; size=64)
    lw = line.thickness * scale / 2
    xmid, y0 = position
    x0 = xmid - lw
    x1 = xmid + lw
    y1 = y0 + line.height
    points = Point2f0[(x0, y0), (x0, y1), (x1, y1), (x1, y0)]
    poly!(ax, points .* size, color=:black, shading=false, linewidth=0)
end

function draw_texelement!(ax, line::HLine, position, scale ; size=64)
    lw = line.thickness * scale / 2
    x0, ymid = position
    x1 = x0 + line.width
    y0 = ymid - lw
    y1 = ymid + lw
    points = Point2f0[(x0, y0), (x0, y1), (x1, y1), (x1, y0)]
    mesh!(ax, points .* size, color=:black, shading=false)
end

draw_texelement_helpers!(args... ; size=64) = nothing

# TODO Move helper rect to the module
function draw_texelement_helpers!(ax, texchar::TeXChar, position, scale ; size=64)
    x = position[1] * size
    # Characters are drawn from the bottom left of the font bounding box but
    # their position is relative to the baseline, so we need to offset them
    y = position[2] * size
    w = inkwidth(texchar) * size * scale
    h = topinkbound(texchar) * size * scale
    a = advance(texchar) * size * scale
    d = bottominkbound(texchar) * size * scale
    left = leftinkbound(texchar) * size * scale
    right = rightinkbound(texchar) * size * scale

    # The space between th origin and the left ink bound
    mesh!(ax,
        Point2f0[
            (x, y + d),
            (x + left, y + d),
            (x + left, y + h),
            (x, y + h)
        ],
        color=RGBA(1, 1, 0, 0.6),
        shading=false
    )

    # The advance after the right inkbound
    mesh!(ax,
        Point2f0[
            (x + right, y + d),
            (x + a, y + d),
            (x + a, y + h),
            (x + right, y + h)
        ],
        color=RGBA(0, 1, 0, 0.3),
        shading=false
    )

    # The descender
    mesh!(ax,
        Point2f0[
            (x + left, y),
            (x + right, y),
            (x + right, y + d),
            (x + left, y + d)
        ],
        color=RGBA(0, 0, 1, 0.3),
        shading=false
    )

    # The inkbound above the baseline
    mesh!(ax,
        Point2f0[
            (x + left, y),
            (x + right, y),
            (x + right, y + h),
            (x + left, y + h)
        ],
        color=RGBA(1, 0, 0, 0.5),
        shading=false
    )
end

function makie_tex!(ax, latex::LaTeXString ; debug=false, size=64)
    for (elem, pos, scale) in generate_tex_elements(latex)
        draw_texelement!(ax, elem, pos, scale ; size=size)
        if debug
            draw_texelement_helpers!(ax, elem, pos, scale ; size=size)
        end
    end
end

begin  # Quick test
    fig = Figure()
    fig[1, 1] = Label(fig, "LaTeX in Makie.jl", tellwidth=false, textsize=64)
    ax = Axis(fig[2, 1])
    ax.aspect = DataAspect()
    tex = L"\lim_{α →\infty} A^j v_{(a + b)_k}^i \sqrt{2} x!= \sqrt{\frac{1+2}{4+a+x}}\int_{0}^{2π} \sin(x) dx"
    tex = L"Momentum $p_x$"
    makie_tex!(ax, tex, debug=false, size=300)
    fig
end


text:
function parse_commandline()
    @eval begin
        using ArgParse;
    end
    @eval begin
        s = ArgParseSettings(;suppress_warnings=true);
        @add_arg_table s begin
            "--outdir"
                help = "specifies directory to store results in; defaults to '~/MORBIT_BENCHMARKS'"
                arg_type = String
                default = joinpath( ENV["HOME"], "MORBIT_BENCHMARKS" )
            "--filename"
                help = """
                name of results file;
                    possible extensions are 'jld2' and 'csv',
                    if no extension is given, results are stored in both 'jld2' and 'csv'
                """
                arg_type = String
                default = ""
                required = false
            "--runs-per-setting"
                help = "number of random starts per setting."
                default = Threads.nthreads()
                arg_type = Int
            "--resume-from"
                help = "old result file to resume from."
                default = ""
                arg_type = String
            "--save-every"
                help = "save results every X iterations."
                default = 100
                arg_type = Int

            # optinal positional arguments if VSCode does something dumb
            "pos1"
                required = false
            "pos2"
                required = false
            "pos3"
                required = false
            "pos4"
                required = false
            "pos5"
                required = false
        end
        return parse_args(s)
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_glyphrepresentation

"""
    vtk_glyphrepresentation(;kwargs...)
    vtk_glyphrepresentation(children::Any;kwargs...)
    vtk_glyphrepresentation(children_maker::Function;kwargs...)

A GlyphRepresentation component.
GlyphRepresentation using a source on port=1 as Glyph and the points of the source on port=0 to position the given glyphs
It takes the following set of properties:
- actor: Properties to assign to the vtkActor
- mapper: Properties to assign to the vtkGlyph3DMapper
- property: Properties to assign to the vtkProperty (actor.getProperty())
- colorMapPreset: Name of the preset to use for controlling the color mapping
- colorDataRange: Range to use for the color scale
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
- `id` (String; optional): The ID used to identify this component.
- `actor` (Dict; optional): Properties to set to the actor
- `mapper` (Dict; optional): Properties to set to the vtkGlyph3DMapper
- `property` (Dict; optional): Properties to set to the actor.property
- `colorMapPreset` (String; optional): Preset name for the lookup table color map
- `colorDataRange` (Array of Reals; optional): Data range use for the colorMap
"""
function vtk_glyphrepresentation(; kwargs...)
        available_props = Symbol[:children, :id, :actor, :mapper, :property, :colorMapPreset, :colorDataRange]
        wild_props = Symbol[]
        return Component("vtk_glyphrepresentation", "GlyphRepresentation", "dash_vtk", available_props, wild_props; kwargs...)
end

vtk_glyphrepresentation(children::Any; kwargs...) = vtk_glyphrepresentation(;kwargs..., children = children)
vtk_glyphrepresentation(children_maker::Function; kwargs...) = vtk_glyphrepresentation(children_maker(); kwargs...)

max_stars_repo_name: ericbdaniels/dash-vtk

max_issues_repo_name: ericbdaniels/dash-vtk

max_forks_repo_name: ericbdaniels/dash-vtk


text:
using ModelicaScriptingTools
using Test
using OMJulia: sendExpression

refdir = "subprojects/shm-ref"
# only simulate full models for two seconds in order to reduce execution time
# of tests and size of reference data files
twosec = Dict("stopTime" => 2, "numberOfIntervals" => 2000)

withOMC("test-output", ".") do omc
    sendExpression(omc, "setCommandLineOptions(\"-d=newInst,nfAPI\")")
    installAndLoad(omc, "Modelica"; version="3.2.3")
    @testset "Simulations" begin
        @testset "Kotani2005" begin
            testmodel(omc, "SHM.Kotani2005.Examples.FullModel.KotaniFullExample"; refdir=refdir, override=twosec)
        end
        @testset "SeidelThesis" begin
            testmodel(omc, "SHM.SeidelThesis.Examples.FullModel.SeidelThesisFullExample"; refdir=refdir, override=twosec)
        end
        @testset "SchoelzelThesis" begin
            @testset "AVDelayExample" begin
                installAndLoad(omc, "SHMConduction") # SHMConduction is required for AVDelayExample
                testmodel(omc, "SHM.SchoelzelThesis.Examples.AVDelayExample"; refdir=refdir)
            end
            @testset "MultiDelayExample" begin
                testmodel(omc, "SHM.SchoelzelThesis.Examples.MultiDelayExample"; refdir=refdir)
            end
            @testset "PacemakerExample" begin
                testmodel(omc, "SHM.SchoelzelThesis.Examples.PacemakerExample"; refdir=refdir)
            end
            @testset "RefractoryDelayExample" begin
                testmodel(omc, "SHM.SchoelzelThesis.Examples.RefractoryDelayExample"; refdir=refdir)
            end
            @testset "RefractoryExample" begin
                testmodel(omc, "SHM.SchoelzelThesis.Examples.RefractoryExample"; refdir=refdir)
            end
            @testset "UnidirectionalModularExample" begin
                testmodel(omc, "SHM.SchoelzelThesis.Examples.UnidirectionalModularExample"; refdir=refdir)
            end
            @testset "BidirectionalContractionExample" begin
                testmodel(omc, "SHM.SchoelzelThesis.Examples.BidirectionalContractionExample"; refdir=refdir)
            end
            @testset "BidirectionalContractionMixinExample" begin
                testmodel(omc, "SHM.SchoelzelThesis.Examples.BidirectionalContractionMixinExample"; refdir=refdir)
            end
        end
        @testset "shm-conduction" begin
            testmodel(omc, "SHMConduction.Examples.ModularExample"; refdir=refdir)
        end
    end
end

max_stars_repo_name: CSchoel/shm

max_issues_repo_name: CSchoel/shm

max_forks_repo_name: CSchoel/shm


text:
using PtFEM

n = 5           # No of nodes
els = n - 1     # No of elements

data = Dict(
  # Frame(nels, nn, ndim, nst, nip, finite_element(nod, nodof))
  :struc_el => Frame(els, n, 3, 1, 1, Line(2, 3)),
  :properties => [
    4.0e6 1.0e6 0.3e6 0.3e6],
  :x_coords => Float64[i for i in 0:n-1],
  :y_coords => zeros(n),
  :z_coords => zeros(n),
  :gamma => zeros(n),     # No rotation of element about local axis (dip)
  :g_num => hcat(1:n-1,2:n)',
  :support => [
    (1, [0 0 0 0 0 0])
  ],
  :loaded_nodes => [
    (n, [0.0 0.0 100.0 0.0 0.0 0.0])
  ]
)

data |> display
println()

@time fem, dis_df, fm_df = p44(data)
println()

dis_df |> display
println()

println("Actions:")
fm_df |> display
println()


text:
using MCMCChains
using Test

## CHAIN TESTS
# Define the experiment
n_iter = 4000
n_name = 3
n_chain = 2

# some sample experiment results
val = randn(n_iter, n_name, n_chain) .+ [1, 2, 3]'
val = hcat(val, rand(1:2, n_iter, 1, n_chain))

# construct a Chains object
chn = Chains(val, start = 1, thin = 2)

# Chains object for discretediag
val_disc = rand(Int16, 200, n_name, n_chain)
chn_disc = Chains(val_disc, start = 1, thin = 2)

@testset "basic chains functions" begin
    @test first(chn) == 1
    @test step(chn) == 2
    @test last(chn) == 7999
    @test size(chn) == (4000, 4, 2)
    @test size(chn[1:1000, :, :]) == (1000, 4, 2)
    @test keys(chn) == ["Param1", "Param2", "Param3", "Param4"]
    @test isa(chn[:,1,:], MCMCChains.AbstractChains)
    @test isa(chn[200:300,"Param1",:], MCMCChains.AbstractChains)
    @test isa(chn[200:300,["Param1", "Param3"],:], MCMCChains.AbstractChains)
    @test isa(chn[200:300,"Param1",1], MCMCChains.AbstractChains)
    @test length(vec(chn[:,1,:].value)) == n_chain * n_iter
    @test all(collect(skipmissing(chn[:,1,1].value)) .== val[:,1,1])
    @test all(chn[:,1,2].value .== val[:,1,2])
    @test all(MCMCChains.indiscretesupport(chn) .== [false, false, false, true])
    @test setinfo(chn, NamedTuple{(:A, :B)}((1,2))).info == NamedTuple{(:A, :B)}((1,2))
    @test isa(set_section(chn, Dict(:internals => ["Param1"])), MCMCChains.AbstractChains)
end

@testset "function tests" begin
    # the following tests only check if the function calls work!
    @test MCMCChains.diag_all(rand(50, 2), :weiss, 1, 1, 1) != nothing
    @test MCMCChains.diag_all(rand(50, 2), :hangartner, 1, 1, 1) != nothing
    @test MCMCChains.diag_all(rand(50, 2), :billingsley, 1, 1, 1) != nothing

    @test isa(discretediag(chn_disc[:,2,:]), Vector{MCMCChains.ChainDataFrame})
    @test isa(gelmandiag(chn[:,1,:]), MCMCChains.ChainDataFrame)
    @test isa(gewekediag(chn[:,1,:]), Vector{MCMCChains.ChainDataFrame})
    @test isa(heideldiag(chn[:,1,:]), Vector{MCMCChains.ChainDataFrame})
    @test isa(rafterydiag(chn[:,1,:]), Vector{MCMCChains.ChainDataFrame})
end

@testset "concatenation tests" begin
    v1 = rand(500, 5, 1)
    v2 = rand(500, 5, 1)
    v3 = rand(500, 4, 1)
    v4 = rand(500, 5, 4)

    c1 = Chains(v1, ["a1", "a2", "a3", "a4", "a5"])
    c2 = Chains(v2, ["a1", "a2", "a3", "a4", "a5"], start = 501)
    c3 = Chains(v3, ["z1", "z2", "z3", "z4"])
    c4 = Chains(v4, ["a1", "a2", "a3", "a4", "a5"])

    # Test dim 1
    c1_2 = cat(c1, c2, dims = 1)
    @test c1_2.value.data == cat(v1, v2, dims=1)
    @test range(c1_2) == 1:1:1000
    @test names(c1_2) == names(c1) == names(c2)
    @test chains(c1_2) == chains(c1) == chains(c2)
    @test c1_2.value == vcat(c1, c2).value

    # Test dim 2
    c1_3 = cat(c1, c3, dims = 2)
    @test c1_3.value.data == cat(v1, v3, dims=2)
    @test range(c1_3) == 1:1:500
    @test names(c1_3) == cat(names(c1), names(c3), dims=1)
    @test chains(c1_3) == chains(c1) == chains(c3)
    @test c1_3.value == hcat(c1, c3).value

    # Test dim 3
    c1_4 = cat(c1, c4, dims = 3)
    @test c1_4.value.data == cat(v1, v4, dims=3)
    @test range(c1_4) == 1:1:500
    @test names(c1_4) == names(c1) == names(c4)
    @test length(chains(c1_4)) == length(chains(c1)) + length(chains(c4))
    @test c1_4.value == chainscat(c1, c4).value
end


text:
#=

  Place number puzzle in Julia.

  http://ai.uwaterloo.ca/~vanbeek/Courses/Slides/introduction.pdf
  """
  Place numbers 1 through 8 on nodes
  - each number appears exactly once
  - no connected nodes have consecutive numbers
       2 - 5
     / | X | \
   1 - 3 - 6 - 8
     \ | X | /
       4 - 7
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function place_number_puzzle(problem,print_solutions=true,all_solutions=true)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>6,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))

    graph = problem[:graph]
    rows,cols = size(graph)
    n     = problem[:n]

    @variable(model, 1 <= x[1:n] <= n, Int)

    @constraint(model, x in CS.AllDifferentSet())

    #=
    for i in 1:rows
        # Here we constraint the domain of t (the difference)
        # to >= 1 (i.e. 2..n)
        t = @variable(model, [1:1], CS.Integers(2:n))
        my_abs(model,x[graph[i,1]],x[graph[i,2]],t[1])
    end
    =#
    # Alternative:
    for (i,j) in eachrow(graph)
        t = @variable(model, [1:1], CS.Integers(2:n))
        my_abs(model,x[i],x[j],t[1])
    end
    # symmetry breaking
    @constraint(model,x[1] <= x[n])

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                # println("solution #$sol")
                x_val = convert.(Integer,JuMP.value.(x; result=sol))
                println("x:$x_val")

            end
        end
    else
        println("status:$status")
    end

    return status
end

problem = Dict(
   :graph => resize_matrix([[1,3], [1,4], [1,5],
                            [2,4], [2,5], [2,6],
                            [3,1], [3,4], [3,7],
                            [4,1], [4,2], [4,3], [4,5], [4,7], [4,8],
                            [5,1], [5,2], [5,4], [5,6], [5,7], [5,8],
                            [6,2], [6,5], [6,8],
                            [7,3], [7,4], [7,5],
                            [8,4], [8,5], [8,6]]),
    :n => 8
)

@time place_number_puzzle(problem)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
__precompile__(true)
module Packing
using GeometryTypes

include("rectangle.jl")
export RectanglePacker

end # module