text:
using Test, Interpolations, DualNumbers, Unitful

@testset "ExtrapTypeStability" begin

    # Test type-stability of 1-dimensional extrapolation
    f(x) = sin((x-3)*2pi/9 - 1)
    xmax = 10
    A = Float64[f(x) for x in 1:xmax]
    itpg = interpolate(A, BSpline(Linear()))

    schemes = (
        Flat,
        Line,
        Reflect,
        Periodic
    )

    for etp in map(E -> @inferred(extrapolate(itpg, E())), schemes),
        x in (
            # In-bounds evaluation
            3.4, 3, dual(3.1),
            # Out-of-bounds evaluation
            -3.4, -3, dual(-3,1),
            13.4, 13, dual(13,1)
        )
        @inferred(etp(x))
    end

    # Test type-stability of 2-dimensional extrapolation with homogeneous scheme
    g(y) = (y/100)^3
    ymax = 4
    A = Float64[f(x)*g(y) for x in 1:xmax, y in 1:ymax]
    itp2 = interpolate(A, BSpline(Linear()))

    for (etp2,E) in map(E -> (extrapolate(itp2, E()), E), schemes),
        x in (
            # In-bounds evaluation
            3.4, 3, dual(3.1),
            # Out-of-bounds evaluation
            -3.4, -3, dual(-3,1),
            13.4, 13, dual(13,1)
        ),
        y in (
            # In-bounds evaluation
            2.1, 2, dual(2.3, 1),
            # Out-of-bounds evaluation
            -2.1, -2, dual(-2.3, 1),
            12.1, 12, dual(12.1, 1)
        )
        @inferred(etp2(x, y))
    end

    A = [1 2; 3 4]
    Af = Float64.(A)
    for B in (A, Af)
        itpg2 = interpolate(B, BSpline(Linear()))
        etp = extrapolate(itpg2, NaN)
        @test typeof(@inferred(etp(dual(1.5,1), dual(1.5,1)))) ==
              typeof(@inferred(etp(dual(6.5,1), dual(3.5,1))))
    end

    # issue #258
    t = [2017, 2018]
    u = u"V"
    y = [1.f0, 2.f0] .* u
    how = Interpolations.Gridded(Interpolations.Linear())
    itp = interpolate((t,), y, Gridded(Linear()))
    etp = extrapolate(itp, 1.f0 * u)
    @test @inferred(itp(2018)) === 2.0f0 * u
    @test @inferred(etp(2018)) === 2.0f0 * u
    @test @inferred(etp(2019)) === 1.0f0 * u

    # Unitful vector knots
    t = [0.0, 1.0] * u"m"
    itp = LinearInterpolation(t, y)
    @test itp(0.5u"m") == 1.5f0 * u
    @test itp([0.5, 0.75]u"m") == [1.5f0, 1.75f0] * u

    # issues #370, #424
    itp = LinearInterpolation([1.0, 2.0, 3.0], [1.0, 2.0, 3.0], extrapolation_bc=missing)
    @test ismissing(itp(4.0))
    @test @inferred(eltype(itp),itp(3.0)) === 3.0
    @test @inferred(eltype(itp),itp(4.0)) === missing
end

ext: jl


text:
struct topic_modeling_dist <: distibution_sample
    α::AbstractArray{Float64,1}
end

#topic_modeling_dist
function log_likelihood!(r::AbstractArray,x::AbstractArray, distibution_sample::topic_modeling_dist , group::Int64 = -1)
    if length(distibution_sample.α) == 0
        return
    end
    @inbounds for i in eachindex(r)
        r[i] = distibution_sample.α[Int64.(x[i])]
    end
end


text:
using MPI, MPIArrays

MPI.Init()
const comm = MPI.COMM_WORLD
const rank = MPI.Comm_rank(comm)
const nb_procs = MPI.Comm_size(comm)
@assert nb_procs == 4

N = 4 # size of the matrix

# Sequentially display the partitons
function display_partitions(a::MPIArray, message)
  if rank != 0
    return
  end

  println("---- $message ----")

  for r in 0:nb_procs-1
    println("local matrix for rank $r:")
    display(getblock(a[localindices(a,r)...]))
    println()
  end
end

# Initial tiled distribution:
A = MPIArray{Int}(comm, (2, 2), N, N)
forlocalpart!(lp -> fill!(lp, rank), A)
sync(A)
display_partitions(A, "Initial distribution")

# Rows only distro
redistribute!(A, 4, 1)
display_partitions(A, "Rows")

#  non-uniform
redistribute!(A, [N], [1,0,0,3])
display_partitions(A, "Non-uniform")

redistribute!(A)
display_partitions(A, "Restored default")

if rank == 0
  println("-------Complete array -------")
  display(A)
  println()
end

# Clean up
free(A)
MPI.Finalize()


text:
"""
Function that replaces a `for` loop by a corresponding `while` loop saving explicitely the *iterator* and its *state*.
"""
function transform_for(expr, ui8::BoxedUInt8)
  @capture(expr, for element_ in iterator_ body__ end) || return expr
  ui8.n += one(UInt8)
  next = Symbol("_iteratornext_", ui8.n)
  state = Symbol("_iterstate_", ui8.n)
  iterator_value = Symbol("_iterator_", ui8.n)
  quote
    $iterator_value = $iterator
    $next = iterate($iterator_value)
    while $next !== nothing
      ($element, $state) = $next
      $(body...)
      $next = iterate($iterator_value, $state)
    end
  end
end

"""
Function that replaces a variable `x` in an expression by `_fsmi.x` where `x` is a known slot.
"""
function transform_slots(expr, symbols::Base.KeySet{Symbol, Dict{Symbol,Any}})
  expr isa Expr || return expr
  expr.head == :let && return transform_slots_let(expr, symbols)
  for i in 1:length(expr.args)
    expr.head == :kw && i == 1 && continue
    expr.head == Symbol("quote") && continue
    expr.args[i] = expr.args[i] isa Symbol && expr.args[i] in symbols ? :(_fsmi.$(expr.args[i])) : expr.args[i]
  end
  expr
end

"""
Function that handles `let` block
"""
function transform_slots_let(expr::Expr, symbols::Base.KeySet{Symbol, Dict{Symbol,Any}})
  @capture(expr, let vars__; body__ end)
  locals = Set{Symbol}()
  for var in vars
    sym = var.args[1].args[2].value
    push!(locals, sym)
    var.args[1] = sym
  end
  body = postwalk(x->transform_let(x, locals), :(begin $(body...) end))
  :(let $((:($var) for var in vars)...); $body end)
end

"""
Function that replaces a variable `_fsmi.x` in an expression by `x` where `x` is a variable declared in a `let` block.
"""
function transform_let(expr, symbols::Set{Symbol})
  expr isa Expr || return expr
  expr.head == :. || return expr
  expr = expr.args[2].value in symbols ? :($(expr.args[2].value)) : expr
end

"""
Function that replaces a `arg = @yield ret` statement by
```julia
  @yield ret;
  arg = arg_
```
where `arg_` is the argument of the function containing the expression.
"""
function transform_arg(expr)
  @capture(expr, arg_ = ex_) || return expr
  _is_yield(ex) || return expr
  ret = length(ex.args) > 2 ? ex.args[3:end] : [nothing]
  quote
    @yield $(ret...)
    $arg = _arg
  end
end

"""
Function that replaces a `@yield ret` or `@yield` statement by
```julia
  @yield ret
  _arg isa Exception && throw(_arg)
```
to allow that an `Exception` can be thrown into a `@resumable function`.
"""
function transform_exc(expr)
  _is_yield(expr) || return expr
  ret = length(expr.args) > 2 ? expr.args[3:end] : [nothing]
  quote
    @yield $(ret...)
    _arg isa Exception && throw(_arg)
  end
end

"""
Function that replaces a `try`-`catch`-`finally`-`end` expression having a top level `@yield` statement in the `try` part
```julia
  try
    before_statements...
    @yield ret
    after_statements...
  catch exc
    catch_statements...
  finally
    finally_statements...
  end
```
with a sequence of `try`-`catch`-`end` expressions:
```julia
  try
    before_statements...
  catch
    catch_statements...
    @goto _TRY_n
  end
  @yield ret
  try
    after_statements...
  catch
    catch_statements...
  end
  @label _TRY_n
  finally_statements...
```
"""
function transform_try(expr, ui8::BoxedUInt8)
  @capture(expr, (try body__ catch exc_; handling__ end) | (try body__ catch exc_; handling__ finally always__ end)) || return expr
  ui8.n += one(UInt8)
  new_body = []
  segment = []
  handling = handling == nothing ? [nothing] : handling
  for ex in body
    if _is_yield(ex)
      ret = length(ex.args) > 2 ? ex.args[3:end] : [nothing]
      exc == nothing ? push!(new_body, :(try $(segment...) catch; $(handling...); @goto $(Symbol("_TRY_", :($(ui8.n)))) end)) : push!(new_body, :(try $(segment...) catch $exc; $(handling...) ; @goto $(Symbol("_TRY_", :($(ui8.n)))) end))
      push!(new_body, quote @yield $(ret...) end)
      segment = []
    else
      push!(segment, ex)
    end
  end
  if segment != []
    exc == nothing ? push!(new_body, :(try $(segment...) catch; $(handling...) end)) : push!(new_body, :(try $(segment...) catch $exc; $(handling...) end))
  end
  push!(new_body, :(@label $(Symbol("_TRY_", :($(ui8.n))))))
  always == nothing || push!(new_body, quote $(always...) end)
  quote $(new_body...) end
end

"""
Function that replaces a `@yield ret` or `@yield` statement with
```julia
  _fsmi._state = n
  return ret
  @label _STATE_n
  _fsmi._state = 0xff
```
"""
function transform_yield(expr, ui8::BoxedUInt8)
  _is_yield(expr) || return expr
  ret = length(expr.args) > 2 ? expr.args[3:end] : [nothing]
  ui8.n += one(UInt8)
  quote
    _fsmi._state = $(ui8.n)
    return $(ret...)
    @label $(Symbol("_STATE_", :($(ui8.n))))
    _fsmi._state = 0xff
  end
end

"""
Function that replaces a `@yield ret` or `@yield` statement with
```julia
  return ret
```
"""
function transform_yield(expr)
  _is_yield(expr) || return expr
  ret = length(expr.args) > 2 ? expr.args[3:end] : [nothing]
  quote
    $(ret...)
  end
end

"""
Function returning whether an expression is a `@yield` macro
"""
_is_yield(ex) = false

function _is_yield(ex::Expr)
  ex.head == :macrocall && ex.args[1] == Symbol("@yield")
end


text:
using Test
using FiniteVolumes

@testset "Interface with other packages" begin
    @testset "Measurements" begin
        using Measurements
        mesh = CartesianMesh(10)
        w0 = map(x -> sin(2π*x), cell_centers(mesh))
        i_face = (FiniteVolumes.CartesianMeshes.Half(11),)

        flux = LinearAdvectionFlux(1.0 ± 1.0)
        @test (Upwind())(flux, mesh, w0, i_face) |> Measurements.value == w0[5]

        dt = 0.01
        w = w0 .- dt*FiniteVolumes.div(flux, mesh, w0, Upwind())
    end

    @testset "Unitful" begin
        using Unitful: kg, m, s
        mesh = CartesianMesh(0.0m, 1.0m, 10)
        i_face = (FiniteVolumes.CartesianMeshes.Half(11),)
        dt = 0.01s

        flux = LinearAdvectionFlux(1m/s)
        w1 = map(x -> sin(2π*x/1m)*1kg/m, cell_centers(mesh))
        @test (Upwind())(flux, mesh, w1, i_face) |> typeof == typeof(w1[1]* 1m/s)
        w2 = w1 .- dt*FiniteVolumes.div(flux, mesh, w1, Upwind())
        @test eltype(w2) == eltype(w1)

    #     flux = Wave1DFlux(1m/s)
    #     w1 = map(x -> (sin(2π*x/1m)*kg/m, 0kg/m/s), cell_centers(mesh))
    #     Δw = map(x -> (0kg/m/s, 0kg/m/s^2), cell_centers(mesh))
    #     w2 = w1 .- dt.*FiniteVolumes.div(flux, mesh, w1, Upwind())
    #     @test eltype(w2) == eltype(w1)
    end

    @testset "ForwardDiff" begin
        using ForwardDiff
        function run(v)
            f = LinearAdvectionFlux(v)
            mesh = CartesianMesh(0.0, 1.0, 100)
            w = [sin(2π*x) for x in cell_centers(mesh)]
            dt = 0.01
            w2 = w .- dt*FiniteVolumes.div(f, mesh, w, Upwind())
        end
        ForwardDiff.derivative(run, 1.0)
    end

    @testset "OrdinaryDiffEq" begin
        using OrdinaryDiffEq
        mesh = CartesianMesh(10)
        flux = LinearAdvectionFlux(1.0)

        w₀ = map(x -> x < 0.5 ? 1.0 : 0.0, cell_centers(mesh))
        dwdt_upwind(w, p, t) = -FiniteVolumes.div(flux, mesh, w, Upwind())
        prob = ODEProblem(dwdt_upwind, w₀, 0.4)
        sol1 = solve(prob, Euler(), dt=0.005, saveat=0.1)
        sol2 = solve(prob, RK4(), dt=0.005, adaptive=false, saveat=0.1)
        sol3 = solve(prob, ImplicitEuler(), dt=0.1, saveat=0.1)

        @test all(0.0 .<= sol1.u[end] .<= 1.0)
        @test all(0.0 .<= sol2.u[end] .<= 1.0)
        @test all(0.0 .<= sol2.u[end] .<= 1.0)
    end

    @testset "Implicit heat equation" begin
        using OrdinaryDiffEq
        mesh = PeriodicCartesianMesh(10)
        w₀ = map(x -> exp(-1000*(x - 0.5)^2), cell_centers(mesh))
        dwdt(w, p, t) = -FiniteVolumes.div(FiniteVolumes.∇, mesh, w, FiniteVolumes.Centered())
        prob = ODEProblem(dwdt, w₀, (0, 0.1))
        sol = solve(prob, ImplicitEuler(), dt=0.02, saveat=0.02)
    end
end


text:
@testset "Auxiliary Functions Test" begin
    @testset "check constant columns" begin
		@test_throws Exception PartialLeastSquaresRegressor.check_constant_cols([1.0 1;1 2;1 3])
		@test_throws Exception PartialLeastSquaresRegressor.check_constant_cols([1.0;1;1][:,:])
		@test_throws Exception PartialLeastSquaresRegressor.check_constant_cols([1.0 2 3])
		@test_throws Exception PartialLeastSquaresRegressor.check_constant_cols([1.0; 1; 1][:,:])

		@test PartialLeastSquaresRegressor.check_constant_cols([1.0 1;2 2;3 3])
		@test PartialLeastSquaresRegressor.check_constant_cols([1.0;2;3][:,:])
	end

	@testset "checkparams" begin
	     #@test_logs PLSRegressor.check_params(2,1,"linear")
		 @test_throws Exception PartialLeastSquaresRegressor.check_params(-1,2,"linear")
		 @test_throws Exception PartialLeastSquaresRegressor.check_params(1,2,"x")

		 @test PartialLeastSquaresRegressor.check_params(1,2,"linear")
	end

	@testset "checkdata" begin
		 @test_throws Exception PartialLeastSquaresRegressor.check_data(zeros(0,0), 0)
		 @test_throws Exception PartialLeastSquaresRegressor.check_data(zeros(1,1), 10)
		 @test PartialLeastSquaresRegressor.check_data(zeros(1,1), 1)
	end
end


text:
using PlasmoAlgorithms
using Base.Test

include("benders_example.jl")

r = bendersolve(g)

@test isapprox(r.objval, 5.28, atol=0.01)


text:
@recipe(Poly) do scene
    Theme(;
        color = theme(scene, :color),
        visible = theme(scene, :visible),
        strokecolor = RGBAf0(0,0,0,0),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        strokewidth = 0.0,
        linestyle = nothing,
    )
end
AbstractPlotting.convert_arguments(::Type{<: Poly}, v::AbstractVector{<: VecTypes}) = convert_arguments(Scatter, v)
AbstractPlotting.convert_arguments(::Type{<: Poly}, v::AbstractVector{<: Union{Circle, Rectangle}}) = (v,)
AbstractPlotting.convert_arguments(::Type{<: Poly}, args...) = convert_arguments(Scatter, args...)
AbstractPlotting.convert_arguments(::Type{<: Poly}, vertices::AbstractArray, indices::AbstractArray) = convert_arguments(Mesh, vertices, indices)

function plot!(plot::Poly{<: Tuple{Union{AbstractMesh, GeometryPrimitive}}})
    mesh!(
        plot, plot[1],
        color = plot[:color], colormap = plot[:colormap], colorrange = plot[:colorrange],
        shading = false, visible = plot[:visible]
    )
    wireframe!(
        plot, plot[1],
        color = plot[:strokecolor], linestyle = plot[:linestyle],
        linewidth = plot[:strokewidth], visible = plot[:visible]
    )
end

function plot!(plot::Poly{<: Tuple{<: AbstractVector{P}}}) where P
    positions = plot[1]
    bigmesh = lift(positions) do p
        polys = GeometryTypes.split_intersections(p)
        merge(GLPlainMesh.(polys))
    end
    mesh!(plot, bigmesh, color = plot[:color], visible = plot[:visible])
    outline = lift(positions) do p
        push!(copy(p), p[1]) # close path
    end
    lines!(
        plot, outline, visible = plot[:visible],
        color = plot[:strokecolor], linestyle = plot[:linestyle],
        linewidth = plot[:strokewidth],
    )
end

function plot!(plot::Poly{<: Tuple{<: AbstractVector{T}}}) where T <: Union{Circle, Rectangle, Rect}
    positions = plot[1]
    position = lift(positions) do rects
        Point.(minimum.(rects))
    end
    markersize = lift(positions, name = "markersize") do rects
        widths.(rects)
    end
    scatter!(
        plot, position,
        marker = T, markersize = markersize, transform_marker = true,
        marker_offset = Vec2f0(0),
        color = plot[:color],
        strokecolor = plot[:strokecolor],
        colormap = plot[:colormap],
        colorrange = plot[:colorrange],
        strokewidth = plot[:strokewidth], visible = plot[:visible]
    )
end
function data_limits(p::Poly{<: Tuple{<: AbstractVector{T}}}) where T <: Union{Circle, Rectangle, Rect}
    xyz = p.plots[1][1][]
    msize = p.plots[1][:markersize][]
    xybb = FRect3D(xyz)
    mwidth = FRect3D(xyz .+ msize)
    union(mwidth, xybb)
end

@recipe(Arrows, points, directions) do scene
    theme = Theme(
        arrowhead = automatic,
        arrowtail = nothing,
        linecolor = :black,
        linewidth = 1,
        arrowsize = 0.3,
        linestyle = nothing,
        scale = Vec3f0(1),
        normalize = false,
        lengthscale = 1.0f0
    )
    # connect arrow + linecolor by default
    get!(theme, :arrowcolor, theme[:linecolor])
    theme
end

# For the matlab/matplotlib users
const quiver = arrows
const quiver! = arrows!
export quiver, quiver!

arrow_head(N, marker) = marker
arrow_head(N, marker::Automatic) = N == 2 ? '▲' : Pyramid(Point3f0(0, 0, -0.5), 1f0, 1f0)

scatterfun(N) = N == 2 ? scatter! : meshscatter!

convert_arguments(::Type{<: Arrows}, x, y, u, v) = (Point2f0.(x, y), Vec2f0.(u, v))
function convert_arguments(::Type{<: Arrows}, x::AbstractVector, y::AbstractVector, u::AbstractMatrix, v::AbstractMatrix)
    (vec(Point2f0.(x, y')), vec(Vec2f0.(u, v)))
end
convert_arguments(::Type{<: Arrows}, x, y, z, u, v, w) = (Point3f0.(x, y, z), Vec3f0.(u, v, w))

function plot!(arrowplot::Arrows{<: Tuple{AbstractVector{<: Point{N, T}}, V}}) where {N, T, V}
    @extract arrowplot (points, directions, lengthscale, arrowhead, arrowsize, arrowcolor)
    headstart = lift(points, directions, lengthscale) do points, directions, s
        map(points, directions) do p1, dir
            dir = arrowplot[:normalize][] ? StaticArrays.normalize(dir) : dir
            Point{N, Float32}(p1) => Point{N, Float32}(p1 .+ (dir .* Float32(s)))
        end
    end
    linesegments!(
        arrowplot, headstart,
        color = arrowplot[:linecolor], linewidth = arrowplot[:linewidth],
        linestyle = arrowplot[:linestyle],
    )
    scatterfun(N)(
        arrowplot,
        map(x-> last.(x), headstart),
        marker = lift(x-> arrow_head(N, x), arrowhead), markersize = arrowsize,
        color = arrowcolor, rotations = directions
    )
end

@recipe(Wireframe) do scene
    default_theme(scene, LineSegments)
end

function argument_conversion(::Type{Wireframe}, x::AbstractVector, y::AbstractVector, z::AbstractMatrix)
    (ngrid(x, y)..., z)
end

xvector(x::AbstractVector, len) = x
xvector(x::ClosedInterval, len) = range(minimum(x), stop=maximum(x), length=len)
xvector(x::AbstractMatrix, len) = x

yvector(x, len) = xvector(x, len)'
yvector(x::AbstractMatrix, len) = x

"""
    `wireframe(x, y, z)`, `wireframe(positions)`, or `wireframe(mesh)`

Draws a wireframe, either interpreted as a surface or as a mesh.
"""
function plot!(plot::Wireframe{<: Tuple{<: Any, <: Any, <: AbstractMatrix}})
    points_faces = lift(plot[1:3]...) do x, y, z
        T = eltype(z); M, N = size(z)
        points = vec(Point3f0.(xvector(x, M), yvector(y, N), z))
        # Connect the vetices with faces, as one would use for a 2D Rectangle
        # grid with M,N grid points
        faces = decompose(Face{2, GLIndex}, SimpleRectangle(0, 0, 1, 1), (M, N))
        view(points, faces)
    end
    linesegments!(plot, Theme(plot), points_faces)
end

Theme(x::AbstractPlot) = x.attributes
# function wireframe(
#         geometry, data::Dict
#     )
#     points = const_lift(geometry) do g
#         decompose(Point3f0, g) # get the point representation of the geometry
#     end
#     # Get line index representation
#     indices = decompose(Face{2, GLIndex}, value(geometry))
#     data[:indices] = reinterpret(GLuint, indices)
#     _default(points, style"linesegment"(), data)
# end
function plot!(plot::Wireframe{Tuple{T}}) where T
    points = lift(plot[1]) do g
        # get the point representation of the geometry
        indices = decompose(Face{2, GLIndex}, g)
        points = decompose(Point3f0, g)
        view(points, indices)
    end
    linesegments!(plot, Theme(plot), points)
end

function sphere_streamline(linebuffer, ∇ˢf, pt, h, n)
    push!(linebuffer, pt)
    df = normalize(∇ˢf(pt[1], pt[2], pt[3]))
    push!(linebuffer, normalize(pt .+ h*df))
    for k=2:n
        cur_pt = last(linebuffer)
        push!(linebuffer, cur_pt)
        df = normalize(∇ˢf(cur_pt...))
        push!(linebuffer, normalize(cur_pt .+ h*df))
    end
    return
end

@recipe(StreamLines, points, directions) do scene
    Theme(
        h = 0.01f0,
        n = 5,
        color = :black,
        linewidth = 1
    )
end

function plot!(plot::StreamLines{<: AbstractVector{T}}) where T
    @extract plot (points, directions)
    linebuffer = T[]
    lines = map(directions, points, plot[:h], plot[:n]) do ∇ˢf, origins, h, n
        empty!(linebuffer)
        for point in origins
            sphere_streamline(linebuffer, ∇ˢf, point, h, n)
        end
        linebuffer
    end
    linesegments!(plot, Theme(plot), lines)
end

@recipe(Series, series) do scene
    Theme(
        seriescolors = :Set1,
        seriestype = :lines
    )
end
convert_arguments(::Type{<: Series}, A::AbstractMatrix{<: Number}) = (A,)
function plot!(sub::Series)
    A = sub[1]
    colors = map_once(sub[:seriescolors], A) do colors, A
        cmap = to_colormap(colors)
        if size(A, 2) > length(cmap)
            @info("Colormap doesn't have enough distinctive values. Please consider using another value for seriescolors")
            cmap = interpolated_getindex.((cmap,), range(0, stop=1, length=M))
        end
        cmap
    end
    map_once(A, sub[:seriestype]) do A, stype
        empty!(sub.plots)
        N, M = size(A)
        for i = 1:M
            c = map(getindex, colors, Node(i))
            attributes = Theme(color = c)
            a_view = view(A, :, i)
            if stype in (:lines, :scatter_lines)
                lines!(sub, attributes, a_view)
            end
            # if stype in (:scatter, :scatter_lines)
            #     scatter!(subsub, attributes, a_view)
            # end
            # subsub
        end
    end
    labels = get(sub, :labels) do
        map(i-> "y $i", 1:size(A[], 2))
    end
    legend!(sub, copy(sub.plots), labels)
    sub
end

"""
    `annotations(strings::Vector{String}, positions::Vector{Point})`

Plots an array of texts at each position in `positions`.
"""
@recipe(Annotations, text, position) do scene
    default_theme(scene, Text)
end

function plot!(plot::Annotations)
    position = plot[2]
    sargs = (
        plot[:model], plot[:font],
        plot[1], position,
        getindex.(Ref(plot), (:color, :textsize, :align, :rotation))...,
    )
    N = value(position) |> eltype |> length
    tp = map(sargs...) do model, font, args...
        if length(args[1]) != length(args[2])
            error("For each text annotation, there needs to be one position. Found: $(length(t)) strings and $(length(p)) positions")
        end
        atlas = get_texture_atlas()
        io = IOBuffer(); combinedpos = Point{N, Float32}[]; colors = RGBAf0[]
        scales = Vec2f0[]; fonts = NativeFont[]; rotations = Quaternionf0[]; alignments = Vec2f0[]
        broadcast_foreach(1:length(args[1]), args...) do idx, text, startpos, color, tsize, alignment, rotation
            # the fact, that Font == Vector{FT_FreeType.Font} is pretty annoying for broadcasting.
            # TODO have a better Font type!
            f = to_font(font)
            f = isa(f, NativeFont) ? f : f[idx]
            c = to_color(color)
            rot = to_rotation(rotation)
            ali = to_align(alignment)
            pos, s = layout_text(text, startpos, tsize, f, alignment, rot, model)
            print(io, text)
            n = length(pos)
            append!(combinedpos, pos)
            append!(scales, s)
            append!(colors, repeated(c, n))
            append!(fonts,  repeated(f, n))
            append!(rotations, repeated(rot, n))
            append!(alignments, repeated(ali, n))
        end
        (String(take!(io)), combinedpos, colors, scales, fonts, rotations)
    end
    t_attributes = Attributes()
    t_attributes[:position] = map(x-> x[2], tp)
    t_attributes[:color] = map(x-> x[3], tp)
    t_attributes[:textsize] = map(x-> x[4], tp)
    t_attributes[:font] = map(x-> x[5], tp)
    t_attributes[:rotation] = map(x-> x[6], tp)
    t_attributes[:align] = Vec2f0(0)
    t_attributes[:model] = Mat4f0(I)
    text!(plot, t_attributes, map(x-> x[1], tp))
    plot
end

is2d(scene::SceneLike) = widths(limits(scene)[])[3] == 0.0

function plot!(scene::SceneLike, subscene::AbstractPlot, attributes::Attributes)
    plot_attributes, rest = merged_get!(:plot, scene, attributes) do
        Theme(
            show_axis = true,
            show_legend = false,
            scale_plot = true,
            center = true,
            axis = Attributes(),
            legend = Attributes(),
            camera = automatic,
            limits = automatic,
            padding = Vec3f0(0.1),
            raw = false
        )
    end
    if plot_attributes[:raw][] == false
        s_limits = limits(scene)
        map_once(plot_attributes[:limits], plot_attributes[:padding]) do limit, padd
            if limit == automatic
                @info("calculating limits")
                @log_performance "calculating limits" begin
                    x = data_limits(scene)
                    # for when scene is empty
                    dlimits = if x == FRect3D(Vec3f0(0), Vec3f0(0))
                        data_limits(subscene)
                    else
                        union(x, data_limits(subscene))
                    end
                    lim_w = widths(dlimits)
                    padd_abs = lim_w .* Vec3f0(padd)
                    s_limits[] = FRect3D(minimum(dlimits) .- padd_abs, lim_w .+  2padd_abs)
                end
            else
                s_limits[] = FRect3D(limit)
            end
        end
        area_widths = RefValue(widths(pixelarea(scene)[]))
        map_once(pixelarea(scene), s_limits, plot_attributes[:scale_plot]) do area, limits, scaleit
            # not really sure how to scale 3D scenes in a reasonable way
            if scaleit && is2d(scene) # && area_widths[] != widths(area)
                area_widths[] = widths(area)
                mini, maxi = minimum(limits), maximum(limits)
                l = ((mini[1], maxi[1]), (mini[2], maxi[2]))
                xyzfit = fit_ratio(area, l)
                s = to_ndim(Vec3f0, xyzfit, 1f0)
                @info("calculated scaling: ", Tuple(s))
                scale!(scene, s)
            end
            return
        end
        if plot_attributes[:show_axis][] && !(any(isaxis, plots(scene)))
            axis_attributes = plot_attributes[:axis]
            if is2d(scene)
                @info("Creating axis 2D")
                axis2d!(scene, axis_attributes, s_limits)
            else
                limits3d = map(s_limits) do l
                    mini, maxi = minimum(l), maximum(l)
                    tuple.(Tuple.((mini, maxi))...)
                end
                @info("Creating axis 3D")
                axis3d!(scene, axis_attributes, limits3d)
            end
        end
        # if plot_attributes[:show_legend][] && haskey(p.attributes, :colormap)
        #     legend_attributes = plot_attributes[:legend][]
        #     colorlegend(scene, p.attributes[:colormap], p.attributes[:colorrange], legend_attributes)
        # end
        if plot_attributes[:camera][] == automatic
            cam = cameracontrols(scene)
            if cam == EmptyCamera()
                if is2d(scene)
                    @info("setting camera to 2D")
                    cam2d!(scene)
                else
                    @info("setting camera to 3D")
                    cam3d!(scene)
                end
            end
        end
    end
    push!(scene.plots, subscene)
    value(plot_attributes[:center]) && center!(scene)
    scene
end

@recipe(Arc, origin, radius, start_angle, stop_angle) do scene
    Theme(;
        default_theme(scene, Lines)...,
        resolution = 361,
    )
end

function plot!(p::Arc)
    args = getindex.(p, (:origin, :radius, :start_angle, :stop_angle, :resolution))
    positions = lift(args...) do origin, radius, start_angle, stop_angle, resolution
        map(range(start_angle, stop=stop_angle, length=resolution)) do angle
            origin .+ (Point2f0(sin(angle), cos(angle)) .* radius)
        end
    end
    lines!(p, Theme(p), positions)
end

function AbstractPlotting.plot!(plot::Plot(AbstractVector{<: Complex}))
    plot[:axis, :labels] = ("Re(x)", "Im(x)")
    lines!(plot, lift(im-> Point2f0.(real.(im), imag.(im)), x[1]))
end

@recipe(BarPlot, x, y) do scene
    Theme(;
        fillto = 0.0,
        color = theme(scene, :color),
        colormap = theme(scene, :colormap),
        colorrange = automatic,
        marker = Rect,
        width = nothing
    )
end

function data_limits(p::BarPlot)
    xy = p.plots[1][1][]
    msize = p.plots[1][:markersize][]
    xybb = FRect3D(xy)
    y = last.(msize) .+ last.(xy)
    bb = AbstractPlotting._boundingbox(first.(xy), y)
    union(bb, xybb)
end

convert_arguments(::Type{<: BarPlot}, x::AbstractVector{<: Number}, y::AbstractVector{<: Number}) = (x, y)
convert_arguments(::Type{<: BarPlot}, y::AbstractVector{<: Number}) = (1:length(y), y)

function AbstractPlotting.plot!(p::BarPlot)
    pos_scale = lift(p[1], p[2], p[:fillto], p[:width]) do x, y, fillto, hw
        nx, ny = length(x), length(y)
        cv = x
        x = if nx == ny
            cv
        elseif nx == ny + 1
            0.5diff(cv) + cv[1:end-1]
        else
            error("bar recipe: x must be same length as y (centers), or one more than y (edges).\n\t\tlength(x)=$(length(x)), length(y)=$(length(y))")
        end
        # compute half-width of bars
        if hw == nothing
            hw = mean(diff(x)) # TODO ignore nan?
        end
        # make fillto a vector... default fills to 0
        positions = Point2f0.(cv, Float32.(fillto))
        scales = Vec2f0.(abs.(hw), y)
        offset = Vec2f0.(hw ./ -2f0, 0)
        positions, scales, offset
    end
    scatter!(
        p, lift(first, pos_scale),
        marker = p[:marker], marker_offset = lift(last, pos_scale),
        markersize = lift(getindex, pos_scale, Node(2)),
        color = p[:color], colormap = p[:colormap], colorrange = p[:colorrange],
        transform_marker = true
    )
end

convert_arguments(P::Type{<:AbstractPlot}, f::Function) = convert_arguments(P, f, -5, 5)
convert_arguments(P::Type{<:AbstractPlot}, f::Function, r) = convert_arguments(P, r, f.(r))
convert_arguments(P::Type{<:AbstractPlot}, f::Function, min, max) =
    convert_arguments(P, f, PlotUtils.adapted_grid(f, (min, max)))

@recipe(ScatterLines) do scene
    Theme()
end

function plot!(scene::SceneLike, ::Type{<:ScatterLines}, attributes::Attributes, p...)
    plot!(scene, Lines, attributes, p...)
    plot!(scene, Scatter, attributes, p...)
    scene
end


text:
insts = OrderedDict()
insts["minimal"] = [
    ((3, true, false),),
    ((3, false, true),),
    ((3, false, true), :SOCExpPSD),
    ]
insts["fast"] = [
    ((10, true, false),),
    ((10, false, true),),
    ((10, false, true), :SOCExpPSD),
    ((50, true, false),),
    ((50, false, true),),
    ((50, false, true), :SOCExpPSD),
    ((400, true, false),),
    ((400, false, true),),
    ((400, true, false),),
    ((400, false, true),),
    ((400, false, true), :SOCExpPSD),
    ]
insts["various"] = [
    ((50, false, true), :SOCExpPSD),
    ((1000, true, false),),
    ((1000, false, true),),
    ((2000, true, false),),
    ((2000, false, true),),
    ((4000, true, false),),
    ((4000, false, true),),
    ((8000, true, false),),
    ((8000, false, true),),
    ]
return (PortfolioJuMP, insts)


text:
"""
    interpolation_matrix(x, ξ)

Create interpolation matrix from grid `ξ` to grid `x`.
"""
function interpolation_matrix(x, ξ)
    A = spzeros(length(x), length(ξ))
    ξᵀ = ξ'
    n = ξᵀ[[1], 1:end-1] .≤ x .< ξᵀ[[1], 2:end]
    for m = 1:length(x)-1
        nₘ = findfirst(n[m, :])
        if !isnothing(nₘ)
            A[m, nₘ] = (ξᵀ[nₘ+1] - x[m]) / (ξᵀ[nₘ+1] - ξᵀ[nₘ])
            A[m, nₘ+1] = (x[m] - ξᵀ[nₘ]) / (ξᵀ[nₘ+1] - ξᵀ[nₘ])
        end
    end
    A[x.<ξ[1], 1] .= 1
    A[ξ[end].≤x, end] .= 1
end


text:
generic2arg = Dict{Symbol,Any}(
# :/ => (:(1./x2),:(-x1./abs2(x2))), # (*,N) (V,V) (M,M)
# :\ => (:(-x2./abs2(x1)),:(1./x1)), # (N,*) (V,V) (M,M) (M,V) (V,M)
)

# TODO:

# eval
# scale
# generic_scale!: Not exported
# scale!
# cross
# triu
# tril
# triu!
# tril!
# diff
# gradient
# diag
# generic_vecnormMinusInf: Not exported
# generic_vecnormInf: Not exported
# generic_vecnorm1: Not exported
# norm_sqr: Not exported
# generic_vecnorm2: Not exported
# generic_vecnormp: Not exported
# vecnormMinusInf: Not exported
# vecnormInf: Not exported
# vecnorm1: Not exported
# vecnorm2: Not exported
# vecnormp: Not exported
# vecnorm
# norm
# norm1: Not exported
# norm2: Not exported
# normInf: Not exported
# vecdot
# dot
# rank
# trace
# inv
# \
# /
# cond
# condskeel
# issym
# ishermitian
# istriu
# istril
# isdiag
# linreg
# peakflops
# axpy!: Not exported
# reflector!: Not exported
# reflectorApply!: Not exported
# det
# logdet
# logabsdet
# isapprox


text:
"""
This file will be called when running `make test` from the repo top-level. It
reads all the ".jl" files in the test directory into an array, and runs
`include()` on each of them.
"""

dir = readdir(dirname(@__FILE__))
tests = filter(x -> (endswith(x, "jl") && x != "runtests.jl"), dir)
include.(tests)

max_stars_repo_name: ryapric/boilerplate

max_issues_repo_name: ryapric/boilerplate

max_forks_repo_name: ryapric/boilerplate


text:
@testset "Betweenness" begin
    # self loops
    s2 = DiGraph(3)
    add_edge!(s2, 1, 2); add_edge!(s2, 2, 3); add_edge!(s2, 3, 3)
    s1 = Graph(s2)
    g3 = PathGraph(5)

    gint = loadgraph(joinpath(testdir, "testdata", "graph-50-500.jgz"), "graph-50-500")

    c = vec(readcsv(joinpath(testdir, "testdata", "graph-50-500-bc.txt")))
    for g in testdigraphs(gint)
        z  = @inferred(betweenness_centrality(g))
        zp = @inferred(parallel_betweenness_centrality(g))
        @test all(isapprox.(z, zp))
        @test map(Float32, z)  == map(Float32, c)

        y  = @inferred(betweenness_centrality(g, endpoints=true, normalize=false))
        yp = parallel_betweenness_centrality(g, endpoints=true, normalize=false)
        @test all(isapprox(y, yp))
        @test round.(y[1:3], 4) ==
            round.([122.10760591498584, 159.0072453120582, 176.39547945994505], 4)

        x  = @inferred(betweenness_centrality(g, 3))
        xp = parallel_betweenness_centrality(g, 3)

        x2  = @inferred(betweenness_centrality(g, collect(1:20)))
        xp2 = parallel_betweenness_centrality(g, collect(1:20))

        @test length(x) == 50
        @test length(xp) == 50
        @test length(x2) == 50
        @test length(xp2) == 50
    end

    @test @inferred(betweenness_centrality(s1)) == [0, 1, 0]
    @test parallel_betweenness_centrality(s1)   == [0, 1, 0]

    @test @inferred(betweenness_centrality(s2)) == [0, 0.5, 0]
    @test parallel_betweenness_centrality(s2)   == [0, 0.5, 0]

    g = Graph(2)
    add_edge!(g, 1, 2)
    z  = @inferred(betweenness_centrality(g; normalize=true))
    zp = parallel_betweenness_centrality(g; normalize=true)
    @test all(isapprox(z, zp))
    @test z[1] == z[2] == 0.0
    z2  = @inferred(betweenness_centrality(g, vertices(g)))
    zp2 = parallel_betweenness_centrality(g, vertices(g))
    @test all(isapprox(z2, zp2))
    z3  = @inferred(betweenness_centrality(g, [vertices(g);]))
    zp3 = parallel_betweenness_centrality(g, [vertices(g);])
    @test all(isapprox(z3, zp3))

    @test z == z2 == z3

    z  = @inferred(betweenness_centrality(g3; normalize=false))
    zp = parallel_betweenness_centrality(g3; normalize=false)
    @test all(isapprox(z, zp))
    @test z[1] == z[5] == 0.0
end


text:
@testset "SimpleANN" begin
    @testset "Construction" begin
        let net = SimpleANN(2,1)
            @test net.insize == 2
            @test net.outsize == 1
            @test net.weights == zeros(1,2)
            @test net.thresholds == zeros(1)
        end

        let net = SimpleANN(zeros(3,2), zeros(3))
            @test net.insize == 2
            @test net.outsize == 3
            @test net.weights == zeros(3,2)
            @test net.thresholds == zeros(3)
        end

        @test_throws ArgumentError SimpleANN(0,2)
        @test_throws ArgumentError SimpleANN(2,0)

        @test_throws DimensionMismatch SimpleANN(zeros(3,2), zeros(2))
    end

    @testset "Fire" begin
        let net = SimpleANN([1. -1.], [0.])
            xs = [0,0]
            @test fire(net, xs) == [0]
            @test xs == [0,0]

            xs = [1,0]
            @test fire(net, xs) == [1]
            @test xs == [1,0]

            xs = [0,1]
            @test fire(net, xs) == [0]
            @test xs == [0,1]

            xs = [1,1]
            @test fire(net, xs) == [0]
            @test xs == [1,1]
        end
    end

    @testset "Fire!" begin
        let net = SimpleANN([1. -1.], [0.])
            xs, ys, expect = [0,0], [0], [0]
            @test fire!(net, xs, ys) == expect
            @test ys == expect
            @test xs == [0,0]

            xs, ys, expect = [1,0], [0], [1]
            @test fire!(net, xs, ys) == expect
            @test ys == expect
            @test xs == [1,0]

            xs, ys, expect = [0,1], [0], [0]
            @test fire!(net, xs, ys) == [0]
            @test xs == [0,1]

            xs, ys, expect = [1,1], [0], [0]
            @test fire!(net, xs, ys) == [0]
            @test xs == [1,1]
        end
    end
end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_sastre_m7_opt_rho6_0(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sastre_m7_opt_rho6_0!(A_copy)
end

@inline function exp_sastre_m7_opt_rho6_0!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=10
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    memslots10 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 Ba8 Bb8 B8 T2k10
    # Computing Bb2 = x*I+x*A
    coeff1=0.0015638303136970531
    coeff2=0.5131399756353182
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.0015638303136970531
    coeff2=0.5131399756353182
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.0015678308682144713
    coeff2=-0.001644406122563135
    coeff3=1.0000648973491737
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=-0.009556418198624978
    coeff2=0.0008384641570218439
    coeff3=0.00012764063512709018
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.0025650924538428376
    coeff2=0.0751072955748711
    coeff3=0.013533154589413957
    coeff4=1.1136335440836176
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.029899374961962634
    coeff2=0.17661587215939123
    coeff3=0.05543256972765084
    coeff4=1.0216298885266226
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.0021189804037908145
    coeff2=0.001037076938904245
    coeff3=1.0216723193069674
    coeff4=0.022323328476928236
    coeff5=-0.052946917852281214
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.016491731477218847
    coeff2=0.0008905549732892586
    coeff3=9.855583452895887e-5
    coeff4=-0.0008621015249925553
    coeff5=-2.1080367849877218e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.012857678730185047
    coeff2=0.20416569487077543
    coeff3=0.04205470017452324
    coeff4=0.01518879866885416
    coeff5=0.044060889132843686
    coeff6=1.0071150699762386
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.017396334722190973
    coeff2=0.04276963803433957
    coeff3=0.03219045127574837
    coeff4=-0.07421310609041179
    coeff5=0.09089674435798574
    coeff6=0.9929884164279147
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=3.1087497397974926
    coeff2=0.3694214555352366
    coeff3=0.15935169296586307
    coeff4=7.9059413858030805
    coeff5=1.1026236327353052
    coeff6=-0.048941742425579186
    coeff7=0.24374173790191617
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.009133614148579657
    coeff2=0.14669027898511353
    coeff3=0.0775106388395176
    coeff4=0.020068924005133555
    coeff5=0.2544089375564568
    coeff6=1.9371031112856034
    coeff7=0.9819279402861139
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing Ba8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9717659119346435
    coeff2=0.03708072735488672
    coeff3=-0.0032737750163317585
    coeff4=0.021723450159739942
    coeff5=-0.13300038220488442
    coeff6=-0.005257171591711241
    coeff7=0.02345129553527665
    coeff8=0.9810542134903524
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb8 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9717659119346435
    coeff2=0.03708072735488672
    coeff3=-0.0032737750163317585
    coeff4=0.021723450159739942
    coeff5=-0.13300038220488442
    coeff6=-0.005257171591711241
    coeff7=0.02345129553527665
    coeff8=0.9810542134903524
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B8 with operation: mult
    mul!(memslots10,memslots2,memslots3)
    # Deallocating Ba8 in slot 2
    # Deallocating Bb8 in slot 3
    # Computing T2k10 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7+x*B8
    coeff1=-0.004180474678744686
    coeff2=-0.002570469257093357
    coeff3=-0.0009454274909569563
    coeff4=-0.003167848437529552
    coeff5=-0.003321896755863846
    coeff6=-0.015009539610819524
    coeff7=-0.02988218764945254
    coeff8=0.004141999495986352
    coeff9=1.002949767210698
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9 .+ coeff9.*memslots10
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    # Deallocating B8 in slot 10
    return memslots1 # Returning T2k10
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
module ReactionNetworkImporter

using DataStructures, DiffEqBiological

abstract type NetworkFileFormat end

# exported data types
struct RSSANetwork <: NetworkFileFormat end
struct BNGNetwork <: NetworkFileFormat end
struct ParsedReactionNetwork
    "DiffEqBiological Network"
    rn

    "Initial Conditions"
    u₀

    "Parameters"
    p

    "Parameter Expressions"
    paramexprs

    "Dict from short sym in rn.syms to full sym for species name"
    symstonames

    "Dict from lumped species name (symbol) to group of species ids"
    groupstoids

    "string representation of the generated reaction_network macro"
    rnstr
end
ParsedReactionNetwork(rn, u₀; p=nothing,
                              paramexprs=nothing,
                              symstonames=nothing,
                              groupstoids=nothing,
                              rnstr=nothing) =
                              ParsedReactionNetwork(rn, u₀, p,
                                                    paramexprs,
                                                    symstonames,
                                                    groupstoids,
                                                    rnstr)

export RSSANetwork, BNGNetwork, ParsedReactionNetwork

# parsers
include("parsing_routines_rssafiles.jl")
include("parsing_routines_bngnetworkfiles.jl")

export loadrxnetwork

end # module

max_issues_repo_name: isaacsas/ParseRxNetwork


text:
using PSLQ, LinearAlgebra
using Test

@testset "PSLQ" begin include("pslq.jl") end


text:
# Activation functions

# Some of activation functions need a wrapper for GPU support
# https://github.com/JuliaGPU/CuArrays.jl/issues/614

# @cufunc softplus(x::Real) = ifelse(x > 0, x + log1p(exp(-x)), log1p(exp(x)))

# @cufunc logσ(x::Real) = -softplus(-x)

# @cufunc function gelu(x::Real)
#     p = oftype(x / 1, π)
#     λ = oftype(x / 1, √(2 / p))
#     α = oftype(x / 1, 0.044715)
#     h = oftype(x / 1, 0.5)
#     h * x * (one(x) + tanh(λ * (x + α * x^3)))
# end

# @cufunc lisht(x::Real) = x * tanh(x)

# @cufunc logcosh(x::Real) = x + softplus(-2x) - log(oftype(x, 2))

# @cufunc mish(x::Real) = x * tanh(softplus(x))

# @cufunc tanhshrink(x::Real) = x - tanh(x)


text:
#
# Visualizations related to orbital energy
#

# These are not fully implemented!

"""
Plot the zero velocity curves for the Synodic, normalized CR3BP system.
"""
function zerovelocityplot(orbit::CR3BPOrbit;
                          nondimensional_range = range(-2; stop=2, length=1000),
                          kwargs...)

    curves = zerovelocity_curves(orbit; nondimensional_range = nondimensional_range)

    default = (; title  = "Zero Velocity Curves" * (name(system(orbit)) == (:Primary, :Secondary) ? "" : " for the $(name(system(orbit))[1])-$(name(system(orbit))[2]) System"),
                 xlabel = "X ($(lengthunit(orbit)))",
                 ylabel = "Y ($(lengthunit(orbit)))",
                 labels = :none,
                 formatter = :plain,
                 grid      = :on,
                 linewidth = 2,
                 dpi       = 150)

    options = merge(default, kwargs)

    fig = plot(; options...)
    for curve ∈ curves
        plot!(fig, curve[:,1], curve[:,2]; kwargs...)
    end

    return fig
end

"""
Plot the zero velocity curves for the Synodic, normalized CR3BP
system to the last figure.
"""
function zerovelocityplot!(fig, orbit::CR3BPOrbit;
                          nondimensional_range = range(-2; stop=2, length=1000),
                          kwargs...)

    curves = zerovelocity_curves(orbit; nondimensional_range = nondimensional_range)

    default = (; )
    options = merge(default, kwargs)

    fig = plot!(; options...)
    for (i,curve) ∈ zip(1:length(curves), curves)
        plot!(fig, curve[:,1], curve[:,2]; kwargs...)
    end

    return fig
end


text:
module Synchronizers

using Sockets

function stack(io::IO,msg::Vector{UInt8})
    frontbytes = reinterpret(UInt8,Int16[length(msg)])
    item = UInt8[frontbytes...,msg...]
    write(io,item)
end

function unstack(io::IO)
    sizebytes = [read(io,UInt8),read(io,UInt8)]
    size = reinterpret(Int16,sizebytes)[1]

    msg = UInt8[]
    for i in 1:size
        push!(msg,read(io,UInt8))
    end
    return msg
end

stack(io::IO,n::Int) = stack(io,UInt8[reinterpret(UInt8,Int[n])...])
unstack(io::IO,x::Type{Int}) = reinterpret(Int,unstack(io))[1]

struct Record
    fname::AbstractString
    data::Vector{UInt8}
end

import Base.==
==(a::Record,b::Record) = a.fname==b.fname && a.data==b.data

function stack(socket::IO,record::Record)
    binary = UInt8[record.fname...,UInt8('\n'),record.data...]
    stack(socket,binary)
end

function unstack(socket::IO,x::Type{Record})
    binary = unstack(socket)
    n = findfirst(x->x==UInt8('\n'),binary)
    fname = String(binary[1:n-1])
    data = binary[n+1:end]
    return Record(fname,data)
end

### I could test whether the stuff works for the buffer.

mutable struct Ledger
    dir::AbstractString
    records::Vector{Record}
end

==(a::Ledger,b::Ledger) = a.records==b.records

import Base.length
length(ledger::Ledger) = length(ledger.records)

function readfile(fname::AbstractString)
    file = open(fname,"r")
    data = UInt8[]
    while !eof(file)
        push!(data,read(file,UInt8))
    end
    close(file)
    return data
end

function Ledger(rootdir::AbstractString)
    mkpath(rootdir)
    records = Record[]
    time = []

    ### I actually need to sort them by date
    for (dir, _, files) in walkdir(rootdir)
        for fname in files
            fullpath = "$dir/$fname"
            data = readfile(fullpath)
            record = Record(relpath(fullpath,rootdir),data)
            push!(records,record)
            push!(time,ctime(fullpath))
        end
    end

    if length(records)>1
        sp = sortperm(time)
        records = records[sp]
    end

    return Ledger(rootdir,records)
end

### One can latter extend this method to specify what to do with a specific type
import Base.push!
function push!(l::Ledger,r::Record)
    push!(l.records,r)
    mkpath(dirname(l.dir * r.fname))
    write(l.dir * r.fname,r.data)
end

# Now the question is how to send and receive the record

function serve(port,ledger::Ledger)
    server = listen(port)
    while true
        socket = accept(server)
        @async begin
            n = length(ledger)::Int
            stack(socket,n)

            while true
                m = unstack(socket,Int)
                if m==0
                    n = length(ledger)::Int
                    stack(socket,n)
                else
                    stack(socket,ledger.records[m])
                end
            end
        end
    end
end

# If I make a daemon, how should I name this?
struct Synchronizer
    socket
    n ### The global state
    ledger::Ledger
end

function Synchronizer(port,ledger::Ledger)
    socket = connect(port)
    n = unstack(socket,Int)
    return Synchronizer(socket,n,ledger)
end

import Base.isready
function isready(s::Synchronizer)
    if length(s.ledger)<s.n
        return true
    else
        stack(s.socket,0)
        n = unstack(s.socket,Int) ### Perhaps I could have unstack with a Datatype!!!
        if n>s.n
            s.n = n
            return true
        else
            return false
        end
    end
end

function sync(s::Synchronizer)
    while isready(s)
        m  = length(s.ledger) + 1
        stack(s.socket,m)
        record = unstack(s.socket,Record)
        push!(s.ledger,record)
    end
end

end # module


text:
data=vec(readdlm("../kepler62ef_planets.txt",',',Float64) )
include("../test_ttv.jl")
@time ttv1,ttv2=test_ttv(5,40,20,data); # This call includes time to compile
Profile.clear_malloc_data()
@time ttv1,ttv2=test_ttv(5,40,20,data,WriteOutput=false,num_evals=100000);

println("# Ignore this: ",hash(ttv1)+hash(ttv2)) # This just makes sure optimizer doesn't optimize away important calculations.

# Test that results haven't changed much
data_inner_out=readdlm("inner_ttv.txt",Float64)
data_outer_out=readdlm("outer_ttv.txt",Float64)
data_inner_ref=readdlm("inner_ttv.txt.ref",Float64)
data_outer_ref=readdlm("outer_ttv.txt.ref",Float64)
max_abs_diff = max(maximum(abs(data_inner_out.-data_inner_ref)),maximum(abs(data_outer_out.-data_outer_ref)))
println("# Maximum difference: ",max_abs_diff)

max_stars_repo_name: ericagol/TTVFaster

max_issues_repo_name: ericagol/TTVFaster

max_forks_repo_name: ericagol/TTVFaster


text:
module FieldsOperations

using Gridap
using Gridap.Helpers

export change_basis
import Gridap: evaluate!
import Gridap: return_size
import Gridap: gradient
import Gridap.MapApply: MapFromKernel

function change_basis(basis::Basis,changeofbasis::Matrix)
  BasisFromChangeOfBasis(basis,changeofbasis)
end

mutable struct FieldLikeAndGradient{D,T,N,V,G} <: FieldLike{D,T,N}
  val::V
  grad::G
end

function FieldLikeAndGradient(
  val::FieldLike{D,T,N},
  grad::FieldLike{D,TG,N}) where {D,T,TG,N}
  V = typeof(val)
  G = typeof(grad)
  FieldLikeAndGradient{D,T,N,V,G}(val,grad)
end

function evaluate!(
  this::FieldLikeAndGradient{D,T,N},
  points::AbstractVector{<:Point{D}},
  v::AbstractArray{T,N}) where {D,T,N}
  evaluate!(this.val,points,v)
end

function return_size(
  f::FieldLikeAndGradient,s::NTuple{N,Int} where N)
  return_size(f.val,s)
end

gradient(f::FieldLikeAndGradient) = f.grad

struct BasisFromChangeOfBasisGrad{D,T,C,B<:Basis} <: Basis{D,T}
  basis::B
  changeofbasis::Matrix{C}
  v::Vector{T}
end

function BasisFromChangeOfBasisGrad(
  basis::Basis{D,T},changeofbasis::Matrix{C}) where {D,T,C}
  B = typeof(basis)
  S = Base._return_type(*,(C,T))
  v = zeros(S,length(basis))
  BasisFromChangeOfBasisGrad{D,S,C,B}(basis,changeofbasis,v)
end

struct BasisFromChangeOfBasis{D,T,C,B<:Basis,G} <: Basis{D,T}
  basis::B
  changeofbasis::Matrix{C}
  g::G
  v::Vector{T}
end

function BasisFromChangeOfBasis(
  basis::Basis{D,T},changeofbasis::Matrix{C}) where {D,T,C}
  B = typeof(basis)
  g = BasisFromChangeOfBasisGrad(gradient(basis),changeofbasis)
  G = typeof(g)
  S = Base._return_type(*,(C,T))
  v = zeros(S,length(basis))
  BasisFromChangeOfBasis{D,S,C,B,G}(basis,changeofbasis,g,v)
end

const _BasisFromChangeOfBasis{D,T} = Union{
  BasisFromChangeOfBasis{D,T}, BasisFromChangeOfBasisGrad{D,T}}

function evaluate!(
  this::_BasisFromChangeOfBasis{D,T},
  points::AbstractVector{<:Point{D}},
  v::AbstractArray{T,N}) where {D,T,N}
  evaluate!(this.basis,points,v)
  ndofs, npoints = size(v)
  for j in 1:npoints
    for i in 1:ndofs
      a = zero(T)
      for k in 1:ndofs
        @inbounds a += this.changeofbasis[i,k]*v[k,j]
      end
      this.v[i] = a
    end
    for i in 1:ndofs
      @inbounds v[i,j] = this.v[i]
    end
  end
end

function return_size(
  f::_BasisFromChangeOfBasis,s::NTuple{N,Int} where N)
  return_size(f.basis,s)
end

gradient(f::BasisFromChangeOfBasis) = f.g

gradient(f::BasisFromChangeOfBasisGrad) = @notimplemented

function gradient(m::MapFromKernel)
  _gradient(m.kernel,m.inputs...)
end

_gradient(k,ms...) = @notimplemented

for op in (:+, :-)
  @eval begin
    function _gradient(k::ArrayKernelFromBroadcastedFunction{typeof($op)},ms...)
      g = [ gradient(m) for m in ms  ]
      apply(k,g...)
    end
  end
end

end # module


text:
# Copyright (C) 2021 Ilya Baran.  This program is distributed under the terms of the MIT license.

# AbstractReaction type
abstract type AbstractReaction end

function apply!(reaction::AbstractReaction, derivatives::Variables, t::Number, variables::Variables, parameters::Variables, extraVariables::Dict{String, Any})::Nothing
    error("Function apply! must be defined for concrete reactions");
end

function reactionName(reaction::AbstractReaction)::String reaction.name end

# SimpleReaction is of the form aA + bB + ... --> cC + dD + ...
struct SimpleReaction <: AbstractReaction
    name::String
    reactants::Array{String}
    products::Array{String}
    rateParameter::String
end

function apply!(reaction::SimpleReaction, derivatives::Variables,  t::Number, variables::Variables, parameters::Variables, extraVariables::Dict{String, Any})::Nothing
    rate = parameters[reaction.rateParameter];
    for reactant = reaction.reactants
        rate *= variables[reactant];
    end
    #println("Reaction $(reaction.name) has rate $rate and parameter $(parameters[reaction.rateParameter])");
    for reactant = reaction.reactants
        derivatives[reactant] -= rate;
    end
    for product = reaction.products
        derivatives[product] += rate;
    end
end

# GeneralRateReaction has a function that computes the rate, then subtracts it from reactants and adds it to products
struct GeneralRateReaction <: AbstractReaction
    name::String
    reactants::Array{String}
    products::Array{String}
    rate::Function # Maps time, variables, parameters, and extraVariables to a rate.  time and extraVariables are optional, but if extraVariables are provided, time must be also
    usesT::Bool
    usesTAndExtraVariables::Bool
    function GeneralRateReaction(name, reactants, products, rate)
        variables = Variables(Array{Variable, 1}());
        usesT = !applicable(rate, variables, variables);
        usesTAndExtraVariables = usesT && !applicable(rate, 0, variables, variables);
        new(name, reactants, products, rate, usesT, usesTAndExtraVariables);
    end
end

function apply!(reaction::GeneralRateReaction, derivatives::Variables, t::Number, variables::Variables, parameters::Variables, extraVariables::Dict{String, Any})::Nothing
    if reaction.usesTAndExtraVariables
        rate = reaction.rate(t, variables, parameters, extraVariables);
    elseif reaction.usesT
        rate = reaction.rate(t, variables, parameters);
    else
        rate = reaction.rate(variables, parameters);
    end
    #println("Reaction $(reaction.name) has rate $rate");
    for reactant = reaction.reactants
        derivatives[reactant] -= rate;
    end
    for product = reaction.products
        derivatives[product] += rate;
    end
end

# GeneralReaction has a function that modifies the derivatives
struct GeneralReaction <: AbstractReaction
    name::String
    apply!::Function # Modifies derivatives (first argument).  Also takes time, variables, parameters, and extraVariables (which it can modify).
end

function apply!(reaction::GeneralReaction, derivatives::Variables, t::Number, variables::Variables, parameters::Variables, extraVariables::Dict{String, Any})::Nothing
    reaction.apply!(derivatives, t, variables, parameters, extraVariables);
    nothing;
end

max_stars_repo_name: ilya-baran-personal/modelBlocks

max_issues_repo_name: ilya-baran-personal/modelBlocks

max_forks_repo_name: ilya-baran-personal/modelBlocks


text:
using BlockArrays
using ControlSystems
using LinearAlgebra
using PredictiveControl

const MPC = PredictiveControl

# Some utilities for testing
include( "../../src/utilities.jl" )
include( "../testUtils.jl" )

# Create a sample system
A = [1.0 0.0 0.0 0.0;
     0.0 2.0 0.0 0.0;
     0.0 0.0 3.0 0.0;
     0.0 0.0 0.0 4.0];
B = [1.0 0.0;
     0.0 0.0;
     0.0 1.0;
     0.0 0.0];
C = to_matrix( Float64, I, 4 )

sys = StateSpace( A, B, C, 0, 0.1 )

N  = 5
nx = ControlSystems.nstates( sys )
nu = ControlSystems.ninputs( sys )

###################################################################
# Test the propagation of the input through the system
###################################################################

# Horizon must be greater than 1
@test_throws DomainError MPC.prediction( sys, -1 )
@test_throws DomainError MPC.prediction( sys,  0 )

# Form a prediction matrix to test
mat = MPC.prediction( sys, N )

@test blocksize( mat ) == (N, N)
@test size( mat ) == (nx*N, nu*N)

# Make sure the blocks have the proper values for the system
for i = 1:N, j = i:N
    b = view( mat, Block( j, i ) )

    @test size( b ) == (nx, nu)
    @test b == A^(j-i) * B
end

###################################################################
# Test the propagation of the initial state through the controlled
# system to form the new input sequence
###################################################################
K = [1.0 0.0 1.0 0.0;
     0.0 1.0 0.0 1.0]

csys = StateSpace( A - B*K, B, C, 0, 0.1 )

# Horizon must be greater than 1
@test_throws DomainError MPC.inputprediction( csys, K, -1 )
@test_throws DomainError MPC.inputprediction( csys, K,  0 )

# Form a prediction matrix to test
mat = MPC.inputprediction( csys, K, N )

@test blocksize( mat ) == (N, N)
@test size( mat ) == (nu*N, nu*N)

# Make sure the blocks have the proper values for the system
for i = 1:N, j = i:N
    b = view( mat, Block( j, i ) )

    @test size( b ) == (nu, nu)

    if i == j
        @test b == 1.0*I(nu)
    else
        @test b == -K * csys.A^((j-1)-i) * csys.B
    end
end


text:
using Test
using PolynomialBasis
# using Revise
using CutCellDG
include("../useful_routines.jl")

x0 = [1.0, 5.0]
L, W = 4.0, 1.0
nelmtsx = 2
nelmtsy = 1
numghostlayers = 2
nodesperelement = 9

mesh = CutCellDG.CGMesh(x0, [L, W], [nelmtsx, nelmtsy], nodesperelement)
paddedmesh = CutCellDG.BoundaryPaddedMesh(mesh, numghostlayers)

dx = 1.0
dy = 0.5

testbottomcoords = [
    1.0 1.0 2.0 2.0 3.0 3.0 4.0 4.0 5.0 5.0
    4.0 4.5 4.0 4.5 4.0 4.5 4.0 4.5 4.0 4.5
]
testrightcoords = [
    6.0 6.0 6.0 7.0 7.0 7.0
    5.0 5.5 6.0 5.0 5.5 6.0
]
testtopcoords = [
    1.0 1.0 2.0 2.0 3.0 3.0 4.0 4.0 5.0 5.0
    6.5 7.0 6.5 7.0 6.5 7.0 6.5 7.0 6.5 7.0
]
testleftcoords = [
    -1.0 -1.0 -1.0 0.0 0.0 0.0
    5.0 5.5 6.0 5.0 5.5 6.0
]

@test allapprox(
    CutCellDG.bottom_ghost_coordinates(paddedmesh),
    testbottomcoords,
)
@test allapprox(CutCellDG.right_ghost_coordinates(paddedmesh), testrightcoords)
@test allapprox(CutCellDG.top_ghost_coordinates(paddedmesh), testtopcoords)
@test allapprox(CutCellDG.left_ghost_coordinates(paddedmesh), testleftcoords)

x0 = [0.0, 0.0]
L, W = 2.0, 1.0
nelmtsx, nelmtsy = 2, 1
numghostlayers = 1
polyorder = 2

xI = [0.0, 0.5]
normal = [0.0, 1.0]
tol = 1e-8
boundingradius = 6.0

levelsetbasis = LagrangeTensorProductBasis(2, polyorder)
dim, nf = size(interpolation_points(levelsetbasis))
mesh = CutCellDG.CGMesh(x0, [L, W], [nelmtsx, nelmtsy], nf)

levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, normal, xI),
    mesh,
    levelsetbasis,
)

paddedmesh = CutCellDG.BoundaryPaddedMesh(
    CutCellDG.background_mesh(levelset),
    numghostlayers,
)
cutmesh = CutCellDG.CutMesh(mesh, levelset)

refseedpoints, seedcellids = CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
spatialseedpoints =
    CutCellDG.map_to_spatial(refseedpoints, seedcellids, cutmesh)

paddedlevelset = CutCellDG.BoundaryPaddedLevelSet(
    paddedmesh,
    spatialseedpoints,
    seedcellids,
    levelset,
    tol,
    boundingradius,
)

Dmx = CutCellDG.first_order_horizontal_backward_difference(paddedlevelset)
Dpx = CutCellDG.first_order_horizontal_forward_difference(paddedlevelset)
Dmy = CutCellDG.first_order_vertical_backward_difference(paddedlevelset)
Dpy = CutCellDG.first_order_vertical_forward_difference(paddedlevelset)

numnodes = CutCellDG.number_of_nodes(mesh)
@test allapprox(Dmy, ones(length(numnodes)))
@test allapprox(Dpy, ones(length(numnodes)))
@test allapprox(Dmx, zeros(length(numnodes)))
@test allapprox(Dpx, zeros(length(numnodes)))

x0 = [0.0, 0.0]
L, W = 2.0, 1.0
nelmtsx, nelmtsy = 2, 1
numghostlayers = 1
polyorder = 2

xI = [0.5, 0.0]
normal = [1.0, 0.0]
tol = 1e-8
boundingradius = 6.0

levelsetbasis = LagrangeTensorProductBasis(2, polyorder)
dim, nf = size(interpolation_points(levelsetbasis))
mesh = CutCellDG.CGMesh(x0, [L, W], [nelmtsx, nelmtsy], nf)
numnodes = CutCellDG.number_of_nodes(mesh)

levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, normal, xI),
    mesh,
    levelsetbasis,
)
paddedmesh = CutCellDG.BoundaryPaddedMesh(
    CutCellDG.background_mesh(levelset),
    numghostlayers,
)
cutmesh = CutCellDG.CutMesh(mesh, levelset)

refseedpoints, seedcellids = CutCellDG.seed_zero_levelset(2, levelset, cutmesh)
spatialseedpoints =
    CutCellDG.map_to_spatial(refseedpoints, seedcellids, cutmesh)

paddedlevelset = CutCellDG.BoundaryPaddedLevelSet(
    paddedmesh,
    spatialseedpoints,
    seedcellids,
    levelset,
    tol,
    boundingradius,
)

Dmx = CutCellDG.first_order_horizontal_backward_difference(paddedlevelset)
Dpx = CutCellDG.first_order_horizontal_forward_difference(paddedlevelset)
Dmy = CutCellDG.first_order_vertical_backward_difference(paddedlevelset)
Dpy = CutCellDG.first_order_vertical_forward_difference(paddedlevelset)

@test allapprox(Dmx, ones(length(numnodes)))
@test allapprox(Dpx, ones(length(numnodes)))
@test allapprox(Dmy, zeros(length(numnodes)))
@test allapprox(Dpy, zeros(length(numnodes)))


text:
bootstrapsample(xs::AbstractVector) = xs[rand(1:length(xs), length(xs))]

function twotailedquantile(quantile::AbstractFloat)
    0. ≤ quantile ≤ 1. || throw(DomainError(quantile))
    lq = (1-quantile)/2
    uq = quantile+lq
    lq, uq
end

"""
    AbstractConfidenceInterval{T<:Real}

A type representing a confidence interval.
"""
abstract type AbstractConfidenceInterval{T<:Real} end

"""
    lower(ci::AbstractConfidenceInterval)

Returns the lower bound of a confidence interval.
"""
lower(ci::AbstractConfidenceInterval) = ci.lower

"""
    upper(ci::AbstractConfidenceInterval)

Returns the upper bound of a confidence interval.
"""
upper(ci::AbstractConfidenceInterval) = ci.upper

"""
    confint(ci::AbstractConfidenceInterval)

Returns the lower and upper bounds of a confidence interval.
"""
HypothesisTests.confint(ci::AbstractConfidenceInterval) = lower(ci), upper(ci)

"""
    quantile(ci::AbstractConfidenceInterval)

Returns the quantile of a confidence interval.
"""
Distributions.quantile(ci::AbstractConfidenceInterval) = ci.quantile

function Base.show(io::IO, ci::AbstractConfidenceInterval)
    print(io, quantile(ci), "CI (", round(lower(ci), digits=PRECISION), ", ",
          round(upper(ci), digits=PRECISION), ")")
end

"""
    ConfidenceInterval(lower, upper, quantile)

A type representing the `lower` lower and `upper` upper bounds of an effect size confidence
interval at the `quantile` quantile.

    ConfidenceInterval(xs, ys, es; quantile)

Calculates a confidence interval for the effect size `es` between two vectors `xs` and `ys`
at the `quantile` quantile.
"""
struct ConfidenceInterval{T<:Real} <: AbstractConfidenceInterval{T}
    lower::T
    upper::T
    quantile::Float64

    function ConfidenceInterval(l::T,u::T,q::Float64) where T<:Real
        l ≤ u || throw(DomainError((l,u), "l > u"))
        0. ≤ q ≤ 1. || throw(DomainError(q))
        new{T}(l,u,q)
    end
end

function ConfidenceInterval(xs::AbstractVector{T}, ys::AbstractVector{T}, es::T;
                            quantile::AbstractFloat) where T<:Real
    0. ≤ quantile ≤ 1. || throw(DomainError(quantile))
    nx = length(xs)
    ny = length(ys)
    σ² = √(((nx+ny)/(nx*ny)) + (es^2 / 2(nx+ny)))
    _, uq = twotailedquantile(quantile)
    z = Distributions.quantile(Normal(), uq)
    ci = z*√σ²
    ConfidenceInterval(es-ci, es+ci, quantile)
end

"""
    BootstrapConfidenceInterval(lower, upper, quantile, bootstrap)

A type representing the `lower` lower and `upper` upper bounds of an effect size confidence
interval at the `quantile` quantile with `bootstrap` boostrap resamples.

    BootstrapConfidenceInterval(xs, ys, f, bootstrap; quantile)

Calculates the effect size confidence interval between two vectors `xs` and `ys` at the
`quantile` quantile by applying `f` to `bootstrap` bootstrap resamples of `xs` and `ys`.
"""
struct BootstrapConfidenceInterval{T<:Real} <: AbstractConfidenceInterval{T}
    lower::T
    upper::T
    quantile::Float64
    bootstrap::Int64

    function BootstrapConfidenceInterval(l::T,u::T,q::Float64,b::Int64) where T<:Real
        l ≤ u || throw(DomainError((l,u), "l > u"))
        0. ≤ q ≤ 1. || throw(DomainError(q))
        b > 1 || throw(DomainError(b))
        new{T}(l,u,q,b)
    end
end

function BootstrapConfidenceInterval(xs::AbstractVector{T}, ys::AbstractVector{T},
                                     f::Function, bootstrap::Integer=1000;
                                     quantile::AbstractFloat) where T<:Real
    0. ≤ quantile ≤ 1. || throw(DomainError(quantile))
    bootstrap > 1 || throw(DomainError(bootstrap))
    es = map(_->f(bootstrapsample(xs), bootstrapsample(ys)), 1:bootstrap)
    # es = @distributed vcat for _ = 1:bootstrap
    #     f(bootstrapsample(xs), bootstrapsample(ys))
    # end
    lq, uq = twotailedquantile(quantile)
    BootstrapConfidenceInterval(Distributions.quantile(es, lq),
        Distributions.quantile(es, uq), quantile, bootstrap)
end


text:
using egads

width = 5.0
rad   = 1.0

(context, status) = EG_open( )
if (status < 0) error("Can't open, failure code: %i", status) end

# build cylinder
datacyl = [0.0, 0.0,  width/2.0,
           0.0, 0.0, -width/2.0,
           rad]
(ebody1, status) = EG_makeSolidBody(context, CYLINDER, datacyl)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# build end points of sphere
datasph = [0.0, 0.0, -width/2.0, rad]
(ebody2, status) = EG_makeSolidBody(context, SPHERE, datasph)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# boolean operations
(emodel, status) = EG_solidBoolean(ebody1, ebody2, FUSION)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# save model
status = EG_saveModel(emodel, "test.egads")
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_saveModel(emodel, "test.step")
if (status != EGADS_SUCCESS) cleanup(status, context) end

# delete
status = EG_deleteObject(ebody1)
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_deleteObject(ebody2)
if (status != EGADS_SUCCESS) cleanup(status, context) end
status = EG_deleteObject(emodel)
if (status != EGADS_SUCCESS) cleanup(status, context) end

# close everything
cleanup( context )


text: struct MarsType <: abstractPlanet end


text:
module BoundingVolumeTrees

using StaticArrays
import StaticArrays.insert
const Point{T} = SVector{3,T}
const Vector3 = Point

# some WIP tree structures intended to replace the older DynamicAABBTree
include("BoundingVolumeTrees/BinaryTrees.jl")
include("BoundingVolumeTrees/IndexedBinaryTrees.jl")
include("BoundingVolumeTrees/BalancedTrees.jl")
include("BoundingVolumeTrees/AABBs.jl")
include("BoundingVolumeTrees/AABBTrees.jl")

#!me ideally this isn't needed.  attempt to have a pool of leaves at indices 1:nleaves reserved
function index_entities(aabbtree::B) where {T,S,B<:Union{IndexedBinaryTree{AVLData{AABBNodeData{T,S}}},
                                                         IndexedBinaryTree{AABBNodeData{T,S}}}}
  [x[2] for x in sort(map(n->(entity(value(data(n))),index(n)),Leaves(bt)), by=first)]
end

#!me this one is optimized, but lacking features like balancing
#!me also the design is a bit grim
# however, it works, so it's the one we use in HalfEdges for now
include("BoundingVolumeTrees/DynamicAABBTrees.jl")

end


text:
# This file was generated, do not modify it. # hide
X0, y0 = @load_reduced_ames;

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
# Broadcasting of AbstractData objects
# https://docs.julialang.org/en/v1/manual/interfaces/#Broadcast-Styles

abstract type DataStyle <: Base.BroadcastStyle end

abstract type DataColumnStyle <: DataStyle end
struct VFStyle{A} <: DataColumnStyle end
DataStyle(::Type{VF{S, A}}) where {S, A} = VFStyle{parent_array_type(A)}()
DataColumnStyle(::Type{VFStyle{A}}) where {A} = VFStyle{A}

abstract type Data1DStyle{Ni} <: DataStyle end
struct IFHStyle{Ni, A} <: Data1DStyle{Ni} end
DataStyle(::Type{IFH{S, Ni, A}}) where {S, Ni, A} =
    IFHStyle{Ni, parent_array_type(A)}()

abstract type DataSlab1DStyle{Ni} <: DataStyle end
DataSlab1DStyle(::Type{IFHStyle{Ni, A}}) where {Ni, A} = IFStyle{Ni, A}

struct IFStyle{Ni, A} <: DataSlab1DStyle{Ni} end
DataStyle(::Type{IF{S, Ni, A}}) where {S, Ni, A} =
    IFStyle{Ni, parent_array_type(A)}()

abstract type DataSlab2DStyle{Nij} <: DataStyle end
struct IJFStyle{Nij, A} <: DataSlab2DStyle{Nij} end
DataStyle(::Type{IJF{S, Nij, A}}) where {S, Nij, A} =
    IJFStyle{Nij, parent_array_type(A)}()

abstract type Data2DStyle{Nij} <: DataStyle end
struct IJFHStyle{Nij, A} <: Data2DStyle{Nij} end
DataStyle(::Type{IJFH{S, Nij, A}}) where {S, Nij, A} =
    IJFHStyle{Nij, parent_array_type(A)}()
DataSlab2DStyle(::Type{IJFHStyle{Nij, A}}) where {Nij, A} = IJFStyle{Nij, A}

abstract type Data1DXStyle{Ni} <: DataStyle end
struct VIFHStyle{Ni, A} <: Data1DXStyle{Ni} end
DataStyle(::Type{VIFH{S, Ni, A}}) where {S, Ni, A} =
    VIFHStyle{Ni, parent_array_type(A)}()
Data1DXStyle(::Type{VIFHStyle{Ni, A}}) where {Ni, A} = VIFHStyle{Ni, A}
DataColumnStyle(::Type{VIFHStyle{Ni, A}}) where {Ni, A} = VFStyle{A}
DataSlab1DStyle(::Type{VIFHStyle{Ni, A}}) where {Ni, A} = IFStyle{Ni, A}

abstract type Data2DXStyle{Nij} <: DataStyle end
struct VIJFHStyle{Nij, A} <: Data2DXStyle{Nij} end
DataStyle(::Type{VIJFH{S, Nij, A}}) where {S, Nij, A} =
    VIJFHStyle{Nij, parent_array_type(A)}()
Data2DXStyle(::Type{VIJFHStyle{Nij, A}}) where {Nij, A} = VIJFHStyle{Nij, A}
DataColumnStyle(::Type{VIJFHStyle{Nij, A}}) where {Nij, A} = VFStyle{A}
DataSlab2DStyle(::Type{VIJFHStyle{Nij, A}}) where {Nij, A} = IJFStyle{Nij, A}

abstract type Data3DStyle <: DataStyle end

Base.Broadcast.BroadcastStyle(::Type{D}) where {D <: AbstractData} =
    DataStyle(D)

# precedence rules

# scalars are broadcast over the data object
Base.Broadcast.BroadcastStyle(
    ::Base.Broadcast.AbstractArrayStyle{0},
    ds::DataStyle,
) = ds

Base.Broadcast.BroadcastStyle(
    ::VFStyle{A1},
    ::IFHStyle{Ni, A2},
) where {Ni, A1, A2} = VIFHStyle{Ni, promote_parent_array_type(A1, A2)}()

Base.Broadcast.BroadcastStyle(
    ::VFStyle{A1},
    ::IJFHStyle{Nij, A2},
) where {Nij, A1, A2} = VIJFHStyle{Nij, promote_parent_array_type(A1, A2)}()

Base.Broadcast.BroadcastStyle(
    ::VFStyle{A1},
    ::VIFHStyle{Ni, A2},
) where {Ni, A1, A2} = VIFHStyle{Ni, promote_parent_array_type(A1, A2)}()

Base.Broadcast.BroadcastStyle(
    ::VFStyle{A1},
    ::VIJFHStyle{Nij, A2},
) where {Nij, A1, A2} = VIJFHStyle{Nij, promote_parent_array_type(A1, A2)}()

Base.Broadcast.BroadcastStyle(
    ::IFHStyle{Ni, A1},
    ::VIFHStyle{Ni, A2},
) where {Ni, A1, A2} = VIFHStyle{Ni, promote_parent_array_type(A1, A2)}()

Base.Broadcast.BroadcastStyle(
    ::IJFHStyle{Nij, A1},
    ::VIJFHStyle{Nij, A2},
) where {Nij, A1, A2} = VIJFHStyle{Nij, promote_parent_array_type(A1, A2)}()

Base.Broadcast.broadcastable(data::AbstractData) = data

@inline function slab(
    bc::Base.Broadcast.Broadcasted{DS},
    inds...,
) where {Ni, DS <: Union{Data1DStyle{Ni}, Data1DXStyle{Ni}}}
    _args = slab_args(bc.args, inds...)
    _axes = (SOneTo(Ni),)
    Base.Broadcast.Broadcasted{DataSlab1DStyle(DS)}(bc.f, _args, _axes)
end

@inline function slab(
    bc::Base.Broadcast.Broadcasted{DS},
    inds...,
) where {Nij, DS <: Union{Data2DStyle{Nij}, Data2DXStyle{Nij}}}
    _args = slab_args(bc.args, inds...)
    _axes = (SOneTo(Nij), SOneTo(Nij))
    Base.Broadcast.Broadcasted{DataSlab2DStyle(DS)}(bc.f, _args, _axes)
end

@inline function column(
    bc::Base.Broadcast.Broadcasted{DS},
    inds...,
) where {N, DS <: Union{Data1DXStyle{N}, Data2DXStyle{N}}}
    _args = column_args(bc.args, inds...)
    _axes = nothing
    Base.Broadcast.Broadcasted{DataColumnStyle(DS)}(bc.f, _args, _axes)
end

@inline function column(
    bc::Base.Broadcast.Broadcasted{DS},
    inds...,
) where {DS <: DataColumnStyle}
    bc
end

@propagate_inbounds function column(
    bc::Union{Data1D, Base.Broadcast.Broadcasted{<:Data1D}},
    i,
    h,
)
    Ref(slab(bc, h)[i])
end

@propagate_inbounds function column(
    bc::Union{Data2D, Base.Broadcast.Broadcasted{<:Data2D}},
    i,
    j,
    h,
)
    Ref(slab(bc, h)[i, j])
end

function Base.similar(
    bc::Union{IJFH{<:Any, Nij, A}, Broadcast.Broadcasted{IJFHStyle{Nij, A}}},
    ::Type{Eltype},
) where {Nij, A, Eltype}
    _, _, _, _, Nh = size(bc)
    PA = parent_array_type(A)
    array = similar(PA, (Nij, Nij, typesize(eltype(A), Eltype), Nh))
    return IJFH{Eltype, Nij}(array)
end

function Base.similar(
    bc::Union{IFH{<:Any, Ni, A}, Broadcast.Broadcasted{IFHStyle{Ni, A}}},
    ::Type{Eltype},
) where {Ni, A, Eltype}
    _, _, _, _, Nh = size(bc)
    PA = parent_array_type(A)
    array = similar(PA, (Ni, typesize(eltype(A), Eltype), Nh))
    return IFH{Eltype, Ni}(array)
end

function Base.similar(
    ::Union{IJF{<:Any, Nij, A}, Broadcast.Broadcasted{IJFStyle{Nij, A}}},
    ::Type{Eltype},
) where {Nij, A, Eltype}
    Nf = typesize(eltype(A), Eltype)
    array = MArray{Tuple{Nij, Nij, Nf}, eltype(A), 3, Nij * Nij * Nf}(undef)
    return IJF{Eltype, Nij}(array)
end

function Base.similar(
    ::Union{IF{<:Any, Ni, A}, Broadcast.Broadcasted{IFStyle{Ni, A}}},
    ::Type{Eltype},
) where {S, Ni, A, Eltype}
    Nf = typesize(eltype(A), Eltype)
    array = MArray{Tuple{Ni, Nf}, eltype(A), 2, Ni * Nf}(undef)
    return IF{Eltype, Ni}(array)
end

function Base.similar(
    bc::Union{VF{<:Any, A}, Broadcast.Broadcasted{VFStyle{A}}},
    ::Type{Eltype},
) where {A, Eltype}
    _, _, _, Nv, _ = size(bc)
    PA = parent_array_type(A)
    array = similar(PA, (Nv, typesize(eltype(A), Eltype)))
    return VF{Eltype}(array)
end

function Base.similar(
    bc::Union{VIFH{<:Any, Ni, A}, Broadcast.Broadcasted{VIFHStyle{Ni, A}}},
    ::Type{Eltype},
) where {Ni, A, Eltype}
    _, _, _, Nv, Nh = size(bc)
    PA = parent_array_type(A)
    array = similar(PA, (Nv, Ni, typesize(eltype(A), Eltype), Nh))
    return VIFH{Eltype, Ni}(array)
end

function Base.similar(
    bc::Union{VIJFH{<:Any, Nij, A}, Broadcast.Broadcasted{VIJFHStyle{Nij, A}}},
    ::Type{Eltype},
) where {Nij, A, Eltype}
    _, _, _, Nv, Nh = size(bc)
    PA = parent_array_type(A)
    array = similar(PA, (Nv, Nij, Nij, typesize(eltype(A), Eltype), Nh))
    return VIJFH{Eltype, Nij}(array)
end

function Base.mapreduce(
    fn::F,
    op::Op,
    bc::Union{
        IJFH{<:Any, Nij, A},
        Base.Broadcast.Broadcasted{IJFHStyle{Nij, A}},
    },
) where {F, Op, Nij, A}
    # mapreduce across DataSlab2D
    _, _, _, _, Nh = size(bc)
    mapreduce(op, 1:Nh) do h
        Base.@_inline_meta
        slabview = @inbounds slab(bc, h)
        mapreduce(fn, op, slabview)
    end
end

function Base.mapreduce(
    fn::F,
    op::Op,
    bc::Union{IFH{<:Any, Ni, A}, Base.Broadcast.Broadcasted{IFHStyle{Ni, A}}},
) where {F, Op, Ni, A}
    # mapreduce across DataSlab1D
    _, _, _, _, Nh = size(bc)
    mapreduce(op, 1:Nh) do h
        Base.@_inline_meta
        slabview = @inbounds slab(bc, h)
        mapreduce(fn, op, slabview)
    end
end

function Base.mapreduce(fn::F, op::Op, bc::IJF{S, Nij}) where {F, Op, S, Nij}
    # mapreduce across DataSlab2D nodes
    mapreduce(op, Iterators.product(1:Nij, 1:Nij)) do (i, j)
        Base.@_inline_meta
        idx = CartesianIndex(i, j, 1, 1, 1)
        node = @inbounds bc[idx]
        fn(node)
    end
end

function Base.mapreduce(fn::F, op::Op, bc::IF{S, Ni}) where {F, Op, S, Ni}
    # mapreduce across DataSlab1D nodes
    mapreduce(op, 1:Ni) do i
        Base.@_inline_meta
        idx = CartesianIndex(i, 1, 1, 1, 1)
        node = @inbounds bc[idx]
        fn(node)
    end
end

function Base.mapreduce(
    fn::F,
    op::Op,
    bc::Union{VF{<:Any, A}, Base.Broadcast.Broadcasted{VFStyle{A}}},
) where {F, Op, A}
    # mapreduce across DataColumn levels
    _, _, _, Nv, _ = size(bc)
    mapreduce(op, 1:Nv) do v
        Base.@_inline_meta
        idx = CartesianIndex(1, 1, 1, v, 1)
        level = @inbounds bc[idx]
        fn(level)
    end
end

function Base.mapreduce(
    fn::F,
    op::Op,
    bc::Union{VIFH{<:Any, Ni, A}, Base.Broadcast.Broadcasted{VIFHStyle{Ni, A}}},
) where {F, Op, Ni, A}
    # mapreduce across columns
    _, _, _, _, Nh = size(bc)
    mapreduce(op, Iterators.product(1:Ni, 1:Nh)) do (i, h)
        Base.@_inline_meta
        columnview = @inbounds column(bc, i, h)
        mapreduce(fn, op, columnview)
    end
end

function Base.mapreduce(
    fn::F,
    op::Op,
    bc::Union{
        VIJFH{<:Any, Nij, A},
        Base.Broadcast.Broadcasted{VIJFHStyle{Nij, A}},
    },
) where {F, Op, Nij, A}
    # mapreduce across columns
    _, _, _, _, Nh = size(bc)
    mapreduce(op, Iterators.product(1:Nij, 1:Nij, 1:Nh)) do (i, j, h)
        Base.@_inline_meta
        columnview = @inbounds column(bc, i, j, h)
        mapreduce(fn, op, columnview)
    end
end

# broadcasting scalar assignment
# Performance optimization for the common identity scalar case: dest .= val
@inline function Base.copyto!(
    dest::AbstractData,
    bc::Base.Broadcast.Broadcasted{<:Base.Broadcast.AbstractArrayStyle{0}},
)
    # TODO: we can write an optimized fill! method here that directly computes offsets for set_struct!
    # for now fallback to the default implementation
    DS = typeof(DataStyle(typeof(dest)))
    return copyto!(
        dest,
        Base.Broadcast.instantiate(
            Base.Broadcast.Broadcasted{DS}(bc.f, bc.args, axes(dest)),
        ),
    )
end

function Base.copyto!(
    dest::IJFH{S, Nij},
    bc::Union{IJFH{S, Nij}, Base.Broadcast.Broadcasted{IJFHStyle{Nij, A}}},
) where {S, Nij, A}
    _, _, _, _, Nh = size(bc)
    @inbounds for h in 1:Nh
        slab_dest = slab(dest, h)
        slab_bc = slab(bc, h)
        copyto!(slab_dest, slab_bc)
    end
    return dest
end

function Base.copyto!(
    dest::IFH{S, Ni},
    bc::Union{IFH{S, Ni}, Base.Broadcast.Broadcasted{IFHStyle{Ni, A}}},
) where {S, Ni, A}
    _, _, _, _, Nh = size(bc)
    @inbounds for h in 1:Nh
        slab_dest = slab(dest, h)
        slab_bc = slab(bc, h)
        copyto!(slab_dest, slab_bc)
    end
    return dest
end

# inline inner slab(::DataSlab2D) copy
@inline function Base.copyto!(
    dest::IJF{S, Nij},
    bc::Union{IJF{S, Nij, A}, Base.Broadcast.Broadcasted{IJFStyle{Nij, A}}},
) where {S, Nij, A}
    @inbounds for j in 1:Nij, i in 1:Nij
        idx = CartesianIndex(i, j, 1, 1, 1)
        dest[idx] = convert(S, bc[idx])
    end
    return dest
end

# inline inner slab(::DataSlab1D) copy
@inline function Base.copyto!(
    dest::IF{S, Ni},
    bc::Base.Broadcast.Broadcasted{IFStyle{Ni, A}},
) where {S, Ni, A}
    @inbounds for i in 1:Ni
        idx = CartesianIndex(i, 1, 1, 1, 1)
        dest[idx] = convert(S, bc[idx])
    end
    return dest
end

# inline inner column(::DataColumn) copy
@inline function Base.copyto!(
    dest::VF{S},
    bc::Union{VF{S, A}, Base.Broadcast.Broadcasted{VFStyle{A}}},
) where {S, A}
    _, _, _, Nv, _ = size(dest)
    @inbounds for v in 1:Nv
        idx = CartesianIndex(1, 1, 1, v, 1)
        dest[idx] = convert(S, bc[idx])
    end
    return dest
end

function Base.copyto!(
    dest::VIFH{S, Ni},
    bc::Union{VIFH{S, Ni, A}, Base.Broadcast.Broadcasted{VIFHStyle{Ni, A}}},
) where {S, Ni, A}
    # copy contiguous columns
    _, _, _, _, Nh = size(dest)
    @inbounds for h in 1:Nh, i in 1:Ni
        col_dest = column(dest, i, h)
        col_bc = column(bc, i, h)
        copyto!(col_dest, col_bc)
    end
    return dest
end

function Base.copyto!(
    dest::VIJFH{S, Nij},
    bc::Union{VIJFH{S, Nij, A}, Base.Broadcast.Broadcasted{VIJFHStyle{Nij, A}}},
) where {S, Nij, A}
    # copy contiguous columns
    _, _, _, _, Nh = size(dest)
    @inbounds for h in 1:Nh, j in 1:Nij, i in 1:Nij
        col_dest = column(dest, i, j, h)
        col_bc = column(bc, i, j, h)
        copyto!(col_dest, col_bc)
    end
    return dest
end


text:
function assemble_cut_cell_displacement_bilinear_form!(
    sysmatrix,
    basis,
    cellquads,
    stiffness,
    mesh,
    cellsign,
    cellid,
    vectosymmconverter,
)

    dim = dimension(mesh)
    jac = jacobian(mesh)
    quad = cellquads[cellsign, cellid]
    cellmatrix = vec(
        displacement_bilinear_form(
            basis,
            quad,
            stiffness[cellsign],
            jac,
            dim,
            vectosymmconverter,
        ),
    )
    nodeids = nodal_connectivity(mesh, cellsign, cellid)
    assemble_cell_matrix!(sysmatrix, nodeids, dim, cellmatrix)
end

function assemble_displacement_bilinear_forms!(
    sysmatrix,
    basis,
    cellquads,
    stiffness,
    mesh,
)

    ncells = number_of_cells(mesh)
    jac = jacobian(mesh)
    dim = dimension(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()

    uniformquad = uniform_cell_quadrature(cellquads)
    uniformbf1 = vec(
        displacement_bilinear_form(
            basis,
            uniformquad,
            stiffness[+1],
            jac,
            dim,
            vectosymmconverter,
        ),
    )
    uniformbf2 = vec(
        displacement_bilinear_form(
            basis,
            uniformquad,
            stiffness[-1],
            jac,
            dim,
            vectosymmconverter,
        ),
    )

    uniformcellmatrices = [uniformbf1, uniformbf2]

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        check_cellsign(cellsign)
        if cellsign == +1
            nodeids = nodal_connectivity(mesh, +1, cellid)
            assemble_cell_matrix!(sysmatrix, nodeids, dim, uniformbf1)
        elseif cellsign == -1
            nodeids = nodal_connectivity(mesh, -1, cellid)
            assemble_cell_matrix!(sysmatrix, nodeids, dim, uniformbf2)
        else
            assemble_cut_cell_displacement_bilinear_form!(
                sysmatrix,
                basis,
                cellquads,
                stiffness,
                mesh,
                +1,
                cellid,
                vectosymmconverter,
            )
            assemble_cut_cell_displacement_bilinear_form!(
                sysmatrix,
                basis,
                cellquads,
                stiffness,
                mesh,
                -1,
                cellid,
                vectosymmconverter,
            )
        end
    end
end


text:
module MyPkg
	export testing
# package code goes here

	function testing(x, y)
		x + y
	end
end # module


text:
"""
    is_leap_year(year)

Return `true` if `year` is a leap year in the gregorian calendar.

"""
function is_leap_year(year)
    if year % 4 == 0
        if year % 100 == 0
            if year % 400 == 0
                return true
            else
                return false
            end
        else
            return true
        end
    else
        return false
    end
end

max_stars_repo_name: drkndl/Julia-Exercism

max_issues_repo_name: drkndl/Julia-Exercism

max_forks_repo_name: drkndl/Julia-Exercism


text:
struct Handler3D
    c_ptr::Ptr{Cvoid}
end

struct OpenGlRenderState
    c_ptr::Ptr{Cvoid}
end

struct View
    c_ptr::Ptr{Cvoid}
end

# void pangolin_create_window(char* window_name, int width, int height, int panel) {
function create_window(name::String, width::Int, height::Int; panel::Int=0)::Nothing
    ccall((:pangolin_create_window, :libpangolin_c_api),
          Cvoid,
          (Cstring, Cint, Cint, Cint),
          name, Cint(width), Cint(height), Cint(panel))
    return
end

function OpenGlRenderState(
        projection_mat::Matrix{Float64},
        model_view_mat::Matrix{Float64})
    # void* pangolin_create_render_state(
    #         float* proj_mat_ptr,
    #         float* model_view_ptr
    #         )
    proj_ptr = Base.unsafe_convert(Ptr{Cdouble}, projection_mat)
    model_view_ptr = Base.unsafe_convert(Ptr{Cdouble}, model_view_mat)
    s_cam = ccall((:pangolin_create_render_state, :libpangolin_c_api),
                  Ptr{Cvoid},
                  (Ptr{Cdouble}, Ptr{Cdouble}),
                  proj_ptr, model_view_ptr)
    OpenGlRenderState(s_cam)
end

# void* pangolin_create_render_state_no_params()
function OpenGlRenderState()
    s_cam = ccall((:pangolin_create_render_state_no_params, :libpangolin_c_api),
                  Ptr{Cvoid},
                  ())
    OpenGlRenderState(s_cam)
end

function get_projection_view_mat(s_cam::OpenGlRenderState)
    out = Matrix{Cdouble}(undef, 4, 4)
    out_ptr = Base.unsafe_convert(Ptr{Cdouble}, out)
    # pangolin_get_projection_model_view_matrix

        ccall((:pangolin_get_projection_model_view_matrix, :libpangolin_c_api),
              Cvoid,
              (Ptr{Cvoid}, Ptr{Cdouble}),
              s_cam.c_ptr, out_ptr)
    return out
end

# void pangolin_set_model_view_matrix(void* s_cam_void_ptr, float* in){
function set_model_view_mat!(s_cam::OpenGlRenderState,
                             model_view::Matrix{Float64};
                             conversion::Bool=true)
    in_ptr = Base.unsafe_convert(Ptr{Cdouble}, model_view)
    # pangolin_get_projection_model_view_matrix
    if conversion
        ccall((:pangolin_set_model_view_matrix_with_conversion, :libpangolin_c_api),
              Cvoid,
              (Ptr{Cvoid}, Ptr{Cdouble}),
              s_cam.c_ptr, in_ptr)
    else
        ccall((:pangolin_set_model_view_matrix, :libpangolin_c_api),
              Cvoid,
              (Ptr{Cvoid}, Ptr{Cdouble}),
              s_cam.c_ptr, in_ptr)
    end
end

function Handler3D(rs::OpenGlRenderState)
    handler = ccall((:pangolin_create_handler3d, :libpangolin_c_api),
                    Ptr{Cvoid},
                    (Ptr{Cvoid},),
                    rs.c_ptr)

    Handler3D(handler)
end

# void* pangolin_create_display(void* renderstate_ptr)
function View(handler::Handler3D)
    d_cam = ccall((:pangolin_create_display, :libpangolin_c_api),
                  Ptr{Cvoid},
                  (Ptr{Cvoid},),
                  handler.c_ptr)
    View(d_cam)
end

function show(v::View; show_window::Bool=true)
    # pangolin_view_show
    ccall((:pangolin_view_show, :libpangolin_c_api),
          Cvoid,
          (Ptr{Cvoid}, Cint),
          v.c_ptr, Cint(show_window))
end

function should_quit()
    status = ccall((:pangolin_should_quit, :libpangolin_c_api),Cint, ())
    return status != 0
end

function activate(v::View, s_cam::OpenGlRenderState)
    ccall((:pangolin_view_active_render_sate, :libpangolin_c_api),
          Cvoid,
          (Ptr{Cvoid}, Ptr{Cvoid}),
          v.c_ptr, s_cam.c_ptr)
end

function glDrawColouredCube()
    ccall((:pangolin_glDrawColouredCube, :libpangolin_c_api),
          Cvoid, ())
end

# void pangolin_finish_frame()
function finish_frame()
    ccall((:pangolin_finish_frame, :libpangolin_c_api), Cvoid, ())
end

# # pangolin_destroy_window()
function destroy_window(name::String)
    ccall((:pangolin_destroy_window, :libpangolin_c_api), Cvoid, (Cstring,), name)
end

"""
int make_window_context_current(char* window_name)
return true if successful
"""
function make_window_context_current(name::String)
    res = ccall(
        (:make_window_context_current, :libpangolin_c_api),
        Cint,
        (Cstring,),
        name
    )
    return (res != 0)
end


text:
# Simple utils for incrementally reading basic SGF files

type SGF
    parts::Vector{AbstractString}
    cindex::Int
    finished::Bool
    error::Bool  # If we encountered an error in parsing
    cplayer::Color  # Track here so we can read without updating board state
end

# Example reading in SGF and updating board:
function playout_sgf(filecontents::AbstractString)
    sgf = go.SGF(filecontents)
    board = go.Board()
    while true
        move = get_next_move(sgf)
        if move == nothing
            break
        end
        play_move(board, move)
    end
    board
end

function load_sgf(filename::AbstractString)
    f = open(filename, "r")
    contents = readall(f)
    SGF(contents)
end

# Returns nothing if the game contains unsupported features
function SGF(contents::AbstractString; debug=false)
    # Check if there is a handicap - we want to ignore these
    handicap_regex = r"A[W|B]\[(|..)\]"
    if match(handicap_regex, contents) != nothing
        debug && println(STDERR, "Contains handicaps")
        return nothing
    end

    moves = split(contents, ";")
    sgf = SGF(moves, 1, false, false, BLACK)
    sgf
end

"find sgf files in a folder"
function find_sgf(folder::AbstractString)
    files = Vector{AbstractString}()
    folders = Vector{AbstractString}()
    push!(folders, folder)
    while length(folders) > 0
        cpath = pop!(folders)
        contents = readdir(cpath)
        for entry in contents
            path = joinpath(cpath, entry)
            if isdir(path)
                push!(folders, path)
            elseif endswith(path, ".sgf")
                push!(files, path)
            end
        end
    end
    files
end

function get_next_move(sgf::SGF)
    # Assumes the move is always right after semicolon
    regex = r"^(W|B)\[(|..)\]"  # W[xx] / B[xx] / W[] / B[]
    move = nothing
    while move == nothing
        if sgf.cindex > length(sgf.parts)
            sgf.finished = true
            return nothing
        end
        part = sgf.parts[sgf.cindex]
        m = match(regex, part)
        sgf.cindex += 1
        if m == nothing
            # Try the next one if this move has no matches
            continue
        else
            color = m.captures[1] == "B" ? BLACK : WHITE
            coords = m.captures[2]
            if color != sgf.cplayer
                println(STDERR, "Incorrect color move in SGF: $(part)")
                sgf.error = true
                return
            end
            if coords == ""
                move = PASS_MOVE
            else
                move = Point(N - (coords[2] - 'a'), coords[1] - 'a' + 1)
            end
            sgf.cplayer = -sgf.cplayer
        end
    end
    move
end

"""
Utility function
Use to filter a list of sgf files to games where both players are above `minrating`
dan: Whether the ratings are in dan
"""
function players_over_rating(filename::AbstractString, minrating::Int; dan=false)
    open(filename) do f
        lines = readall(f)
        if dan
            wrating_regex = r"WR\[([3-9])d\]"
            brating_regex = r"BR\[([3-9])d\]"
        else
            wrating_regex = r"WR\[([0-9]{4})\]"
            brating_regex = r"BR\[([0-9]{4})\]"
        end
        mw = match(wrating_regex, lines)
        mb = match(brating_regex, lines)
        if mw == nothing || mb == nothing
            return false
        end
        if parse(Int, mw.captures[1]) < minrating || parse(Int, mb.captures[1]) < minrating
            return false
        end
        return true
    end
end