text:
using ComplexVisual
@ComplexVisual.import_huge

const CVG = try using ComplexVisualGtk; ComplexVisualGtk; catch nothing; end

const fontface = cv_fontface("sans-serif")

function get_param_n()
    value_min, value_max = 1, 29
    slider_pos = CV_TranslateByOffset(Int)
    slider_pos.value = 3
    set_slider_value_func = z -> begin
        slider_pos.value = round(UInt, max(min(z, value_max), value_min))
        return CV_Response(; redraw_flag=true)
    end
    label_style = cv_black → fontface → cv_fontsize(20)
    rulers=(CV_TickLabelAppearance(; label_style) ↦
                ("%.0f" ⇒ value_min:2:value_max),
            CV_Ruler("" ⇒ value_min+1:2:value_max),
            CV_TickLabelAppearance(; tick_length=0, gap=15, label_style) ↦
                (("n:" ⇒ 0.2),))

    return CV_ParamWithSlider(;
        slider_pos, value_min=value_min-0.1, value_max=value_max+0.1,
        rulers, set_slider_value_func)
end

param_n = get_param_n()

function trafo2(z)
    n = param_n.slider_pos.value
    result = 0.0
    for k in 1:n
      result += (isodd(k) ? +1 : -1)*(z^k)/k
    end
    return result
end

trafo1 = z -> log1p(z)

codomain1 = CV_Math2DCanvas(-2.0 + 2.0im, 2.0 - 2.0im, 100)
codomain2 = CV_Math2DCanvas(-2.0 + 2.0im, 2.0 - 2.0im, 100)

label_style = cv_color(0,0,0) → cv_fontface("sans-serif") → cv_fontsize(20)
rulers=(CV_Ruler(cv_format_ticks("%.0f", -2.0:1.0:2.0...),
    CV_TickLabelAppearance(; label_style)),)

circle = (cv_black → cv_linewidth(3)) ↦ CV_LinePainter(
    cv_arc_lines(0.0, 2π, (1.0,)))

scene = cv_scene_comp_codomains_std((param_n, ), trafo1, trafo2,
    codomain1, codomain2; codomain1_re_rulers=rulers,
    painter1=CV_PortraitPainter(trafo1) → circle,
    painter2=CV_PortraitPainter(trafo2) → circle)
cv_get_redraw_func(scene)()

cv_save_image(cv_get_can_layout(scene), "./Log1pSeries.png")

if CVG !== nothing
    handler = CVG.cvg_visualize(scene)
    CVG.cvg_wait_for_destroy(handler.window)

    CVG.cvg_close(handler); CVG.cv_destroy(scene)
end

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
using FFTW, Plots, SpecialFunctions;
NX=128; NP=64NX; dt = 1/10NX; NT=1024; W=1600; w=W/NP; dx=1/NX;
x=rand(NP); v=2collect(1:NP.>NP/2) .-1.; u()=(x.=mod.(x.+v/2*dt,1))
ik=2π*im*vcat(1,1:NX/2,-NX/2+1:-1); K=zeros(NT, 2); ρ=zeros(NX);
f(g, c)=erf((g-c)/dx)/2; ff(i, c) = f((i+0.5)*dx,c)-f((i-0.5)*dx,c)
d(c)=((mod1(i,NX),ff(i, c)) for i∈(-6:6).+Int(round(c*NX)))
rho()=(fill!(ρ, 0);for j∈d.(x);for k∈j;ρ[k[1]]+=k[2]*w/dx;end;end;ρ)
@gif for t in 1:NT;
  u(); E=real.(ifft((E=fft(rho())./ik;E[1]*=0;E))); u();
  for j∈eachindex(v);v[j]+=sum(k->E[k[1]]*k[2],d(x[j]))*dt;end;
  K[t,:].=(sum(E.^2)*dx,sum(v.^2)*w)./2; scatter(x,v,ylims=(-3,3))
end every 4

ext: jl


text:
using Documenter
import Pkg
using SmartShockFinder

# Define module-wide setups such that the respective modules are available in doctests
DocMeta.setdocmeta!(SmartShockFinder, :DocTestSetup, :(using SmartShockFinder); recursive=true)

# Make documentation
makedocs(
    # Specify modules for which docstrings should be shown
    modules = [SmartShockFinder],
    # Set sitename to SmartShockFinder
    sitename="SmartShockFinder.jl",
    # Provide additional formatting options
    format = Documenter.HTML(
        # Disable pretty URLs during manual testing
        prettyurls = get(ENV, "CI", nothing) == "true",
        # Explicitly add favicon as asset
        # assets = ["assets/favicon.ico"],
        # Set canonical URL to GitHub pages URL
        canonical = "https://trixi-framework.github.io/SmartShockFinder.jl/dev"
    ),
    # Explicitly specify documentation structure
    pages = [
        "Home" => "index.md",
        "Reference" => "reference.md",
        "License" => "license.md"
    ],
    strict = true # to make the GitHub action fail when doctests fail, see https://github.com/neuropsychology/Psycho.jl/issues/34
)

deploydocs(
    repo = "github.com/trixi-framework/SmartShockFinder.jl",
    devbranch = "main",
    push_preview = true
)


text:
##### Beginning of file

import PredictMDFull
import Test

Test.@test(typeof(PredictMDFull._registry_url_list()) <: Vector{String})
Test.@test(length(PredictMDFull._registry_url_list()) > 0)

##### End of file


text:
# This file contains some aliases and generic methods to use neural networks.

export logψ, ∇logψ, ∇logψ!, logψ_and_∇logψ, logψ_and_∇logψ!, grad_cache

"""
    grad_cache(net::N) -> tuple

Creates a tuple holding all derivatives of the network N.
"""
grad_cache(net::NeuralNetwork) = begin
    is_analytic(net) && return RealDerivative(net)
    return WirtingerDerivative(net)
end

# Define various aliases used a bit everywhere in networks and when testing
"""
    logψ(net, σ) -> Number

Evaluates the neural network `net` for the given input `σ`.
This will return the value ⟨σ|ψ⟩ if the network represents a pure state (a ket)
or ⟨σᵣ|ρ̂|σᵪ⟩, where σ=(σᵣ,σᵪ), if the network encodes a mixed state.

This is equivalent to `net(σ)`.

The numerical type of the output is `out_type(net)`.

If `net isa CachedNet` then the computation will be performed efficiently
with minimal allocations.
"""
@inline logψ(net::NeuralNetwork, σ) = net(σ)
@inline log_prob_ψ(net, σ...)       = 2.0*real(net(σ...))
@inline ∇logψ(args...)              = logψ_and_∇logψ(args...)[2]
@inline ∇logψ!(args...)             = logψ_and_∇logψ!(args...)[2]

"""
    logψ_and_∇logψ(net, σ)

Evaluates the neural network `net` for the given input `σ` and returns both the
value and it's gradient w.r.t. to the nework's parameters.

!!  The gradient will be allocated. Don't use this in a hot loop and prefer the
    in-place version `logψ_and_∇logψ!`. The gradient is always allocated by calling
    `grad_cache(net)`.

!!  If `net isa CachedNet` and the gradient has been hand-coded, the hand-coded
    version will be used, otherwise `Zygote` will be used to generate the
    gradient with AD. This case is much slower and incurs a long precompilation
    time.

See also: `logψ`, `logψ_and_∇logψ!`

"""
function logψ_and_∇logψ(net::NeuralNetwork, σ::Vararg{N,V})where {N,V}
    der = grad_cache(net)
    y, der = logψ_and_∇logψ!(der, net, σ)
    return y, der
end

# Autodiff version of the gradient: if you are hitting this, expect
# bad performance. Maybe you should hand-code the gradient of your model?
function logψ_and_∇logψ!(der, net::NeuralNetwork, σ)
    σ = config(σ)
    # Zygote's autodiff: generate the pullback
    y, back = forward(net -> net(σ...), net)

    # This computes the gradient, which is the conjugate of the derivative
    _der = back(Int8(1))[1]

    for key=keys(_der)
        conj!(_der[key])
        copyto!(der[key], _der[key])
    end
    return y, der
end

# Common Operations on networks that should always be defined
# If you hit this section it probably means that you forgot to define one
# of those functions on the type of your neural network.
"""
    input_type(net) -> Type

Returns the numerical `eltype` of the input for efficiently evaluating in a
type-stable way the network.
"""
input_type(net::NeuralNetwork) = error("Not Implemented")

"""
    out_type(net) -> Type

Returns the numerical `eltype` of the output of the network.
"""
out_type(net::NeuralNetwork) = error("Not Implemented")
"""
    input_shape(net)

This returns the shape of the numerical input that should be provided to the
network, inside a tuple. This is basically the shape of the input layer.

If the network encodes a Matrix, then it will be a tuple with two equal-length
states. If the network encodes a Ket it will be a tuple with a single element
of length equal to the number of input units.
"""
input_shape(net::NeuralNetwork) = error("Not Implemented")

"""
    random_input_state(net)

Returns a random input state that can be fed into the network. This is not of
type `State`, but is rather the underlying numerical values.

This method is used for testing and to construct states. The shape of the state
will be the same of `input_shape`.
"""
random_input_state(net::NeuralNetwork) = error("Not Implemented")

"""
    is_analytic(net) -> Bool

Returns true if the network is an analytic function (almost everywhere) or if it
has real weights.
"""
is_analytic(net::NeuralNetwork) = false

"""
    num_params(net) -> Net

Returns the total number of parameters of the neural network.
"""
num_params(net::NeuralNetwork) = sum([length(getfield(net, f)) for f=fieldnames(typeof(net))])

# TODO does this even make sense?!
# the idea was that a shallow-copy of the weights of the net is not
# even a copy....
copy(net::NeuralNetwork) = net


text:
mutable struct BarycentricRefinement{U,D1,T,M} <: AbstractMesh{U,D1,T}
    mesh::M
    parent::M
end

vertextype(br::BarycentricRefinement) = vertextype(br.mesh)
vertices(br::BarycentricRefinement) = vertices(br.mesh)
numvertices(br::BarycentricRefinement) = numvertices(br.mesh)
cells(br::BarycentricRefinement) = cells(br.mesh)
numcells(br::BarycentricRefinement) = numcells(br.mesh)
skeleton(br::BarycentricRefinement, dim::Int) = skeleton(br.mesh, dim)
chart(br::BarycentricRefinement, cell) = chart(br.mesh, cell)

parent(mesh::BarycentricRefinement) = mesh.parent

"""
    barycentric refinement(mesh) -> refined_mesh

Create the mesh obtained by inserting an extra vertex in the barycenters of all cells and
recusively creating fine cells by connecting the barycenter of a k-cell to the already
constructed refined (k-1)-cells on its boundary.
"""
function barycentric_refinement(mesh::Mesh{U,2}) where U

    # Get the points and the faces (segments) we will use in the refinement process
    Edges = skeleton(mesh, 1)

    # Note thier number and use it to find the new number of points and segments(faces)
    # after refienmnts
    NV, NE = numvertices(Edges), numcells(Edges)

    nv = NV + NE
    ne = 2NE

    # define the new vertices array that will hold the old points(coarse) and the
    # new points(finer mesh) of vertices
    verts = similar(mesh.vertices, nv)

    # Now we assign the coarse mesh vetrices (old points) in the upper half excatly
    # as the original coarse mesh
    for i in 1:NV
        verts[i] = mesh.vertices[i]
    end

    # for the second half we store the refienment points and we calculate them using
    # the old points. We insert a point between every pair of points in the original
    # coarse mesh by (point1+ point2)\2
    for (i, Face) in enumerate(cells(mesh))
        verts[NV+i] = cartesian(center(chart(mesh, Face)))
    end

    # Now we create the faces using the new points as well
    # example:
    # oringinal points_index (1,2,3,4) , original faces ((1,2), (2,3), (3,4),(4,5))
    # new points_index(1,2,3,4,5,6,7,8), new_faces((1,5),(5,2),(2,6),(6,3),(3,7),(7,4),(4,8),(8,1))
    #                  a,b      For       E=1 ->   (a,c),(c,b) and so on ..
    edges = similar(Edges.faces, ne)
    for E in 1:NE
        c = NV + E
        a = Edges.faces[E][1]
        b = Edges.faces[E][2]
        edges[2(E-1) + 1] = SVector(a,c)
        edges[2(E-1) + 2] = SVector(c,b)
    end

    # return the new mesh after refienments
    return Mesh(verts, edges)
end

function barycentric_refinement(mesh::Mesh{U,3}) where U

    D1 = 3
    T = coordtype(mesh)

    edges = skeleton(mesh,1)
    faces = skeleton(mesh,2)

    NV, NE, NF = numvertices(mesh), numcells(edges), numcells(faces)

    nv = NV + NE + NF
    verts = Array{vertextype(mesh)}(undef,nv)
    for V in 1 : numvertices(mesh)
        verts[V] = mesh.vertices[V]
    end

    # add a vertex in each edge centroid
    for (E,Edge) in enumerate(cells(edges))
        verts[NV+E] = cartesian(center(chart(edges, Edge)))
    end

    # add a vertex in each face centroid
    for (F,Face) in enumerate(cells(faces))
        verts[NV+NE+F] = cartesian(center(chart(faces, Face)))
    end

    D = copy(transpose(connectivity(edges, faces, identity)))
    rows, vals = rowvals(D), nonzeros(D)

    # add six faces in each coarse face
    nf = 6NF
    fcs = zeros(celltype(mesh), nf)
    for F in 1 : numcells(faces)
        c = NV + NE + F

        for i in nzrange(D, F)
            E = rows[i]
            e = NV + E

            j = abs(vals[i]) # local index of edge E in face F
            a = mesh.faces[F][mod1(j+1,3)]
            b = mesh.faces[F][mod1(j+2,3)]

            fcs[6(F-1)+2(j-1)+1] = index(a,e,c)
            fcs[6(F-1)+2(j-1)+2] = index(b,c,e)
        end
    end

    Nodes = skeleton(mesh, 0)
    node_ctrs = [vertices(Nodes)[node][1] for node in cells(Nodes)]
    Nodes.faces = Nodes.faces[sort_sfc(node_ctrs)]

    fine = Mesh(verts, fcs)
    D = connectivity(Nodes, fine)
    rows, vals = rowvals(D), nonzeros(D)
    sorted_fcs = Vector{celltype(mesh)}()
    for (i,Node) in enumerate(cells(Nodes))
        for k in nzrange(D,i)
            j = rows[k]
            push!(sorted_fcs, fcs[j])
        end
    end

    # sorted_fcs = fcs
    @show length(fcs)
    @show length(sorted_fcs)
    @assert length(fcs) == length(sorted_fcs)

    M = typeof(mesh)
    # BarycentricRefinement{U,3,T,M}(Mesh(verts, sorted_fcs), mesh)
    BarycentricRefinement{U,3,T,M}(Mesh(verts, fcs), mesh)
end

function barycentric_refinement(mesh::Mesh{U,4}) where U

    D1 = dimension(mesh)+1
    T = coordtype(mesh)

    edges = skeleton(mesh,1)
    faces = skeleton(mesh,2)
    tetrs = skeleton(mesh,3)

    NV, NE, NF, NT = numvertices(mesh), numcells(edges), numcells(faces), numcells(tetrs)

    nv = NV + NE + NF + NT
    verts = Array{vertextype(mesh)}(undef,nv)
    for V in 1 : numvertices(mesh)
        verts[V] = mesh.vertices[V]
    end

    # add a vertex in each edge centroid
    for (E,Edge) in enumerate(cells(edges))
        verts[NV+E] = cartesian(center(chart(edges, Edge)))
    end

    # add a vertex in each face centroid
    for (F,Face) in enumerate(cells(faces))
        verts[NV+NE+F] = cartesian(center(chart(faces, Face)))
    end

    # add a vertex in each tet centroid
    for (T,Tetr) in enumerate(cells(tetrs))
        verts[NV+NE+NF+T] = cartesian(center(chart(tetrs, Tetr)))
    end

    D = copy(transpose(connectivity(faces, tetrs, identity)))
    C = copy(transpose(connectivity(edges, faces, identity)))
    rows, vals = rowvals(D), nonzeros(D)

    # add 24 tetrs in each coarse tetr
    nt = 24NT
    idx = 1
    fcs = zeros(celltype(mesh), nt)
    for T in 1 : numcells(tetrs)
        c = NV + NE + NF + T # index of the tetr centroid

        # iterate over all faces of tetr T
        for i in nzrange(D, T)
            F = rows[i]
            f = NV + NE + F # f is the vertex index of the face center

            for k in nzrange(C,F)
                E = rowvals(C)[k]
                e = NV + E # e is the vertex index of the edge center

                a = cells(edges)[E][1]
                b = cells(edges)[E][2]

                if sign(nonzeros(C)[k]) < 0
                    a, b = b, a
                end

                fcs[idx+0] = index(a,e,f,c)
                fcs[idx+1] = index(b,f,e,c)

                idx += 2
            end

            # j = abs(vals[i]) # local index of edge E in face F
            # a = mesh.faces[F][mod1(abs(vals[i])+1,3)]
            # b = mesh.faces[F][mod1(abs(vals[i])+2,3)]
            #
            # fcs[6(F-1)+2(j-1)+1] = index(a,e,c)
            # fcs[6(F-1)+2(j-1)+2] = index(b,c,e)
        end
    end

    # Nodes = skeleton(mesh, 0)
    # node_ctrs = [vertices(Nodes)[node][1] for node in cells(Nodes)]
    # Nodes.faces = Nodes.faces[sort_sfc(node_ctrs)]

    fine = Mesh(verts, fcs)
    for (i,fc) in enumerate(fcs)
        # @show fc
        vs = verts[fc]
        t1 = vs[1]-vs[4]
        t2 = vs[2]-vs[4]
        t3 = vs[3]-vs[4]
        if dot(t1 × t2, t3) < 0
            fcs[i] = @SVector[fc[1],fc[2],fc[4],fc[3]]
        end
    end
    M = typeof(mesh)
    BarycentricRefinement{U,4,T,M}(fine, mesh)
    # D = connectivity(Nodes, fine)
    # rows, vals = rowvals(D), nonzeros(D)
    # sorted_fcs = Vector{celltype(mesh)}()
    # for (i,Node) in enumerate(cells(Nodes))
    #     for k in nzrange(D,i)
    #         j = rows[k]
    #         push!(sorted_fcs, fcs[j])
    #     end
    # end
    #
    # # sorted_fcs = fcs
    # @show length(fcs)
    # @show length(sorted_fcs)
    # @assert length(fcs) == length(sorted_fcs)
    #
    # M = typeof(mesh)
    # BarycentricRefinement{U,3,T,M}(Mesh(verts, sorted_fcs), mesh)
end

children(mesh::Mesh{3,4}, cell) = [24*(cell-1)+1 : 24*cell]
parent(mesh::BarycentricRefinement{3,4}, cell) = div((cell-1),24)+1
parent(mesh::BarycentricRefinement{U,3} where {U}, cell) = div((cell-1),6)+1

"""
    bisecting_refinement(mesh) -> refinement

Construct a refinement of `mesh` by inserting a new vertex on every existing edge.
Every face is subdived in four small faces.

Only defined for 2D meshes.
"""
function bisecting_refinement(mesh::Mesh{U,3}) where U

    D1 = 3

    edges = skeleton(mesh,1)
    faces = skeleton(mesh,2)

    NV, NE, NF = numvertices(mesh), numcells(edges), numcells(faces)

    nv = NV + NE
    verts = Array{vertextype(mesh)}(undef,nv)

    # copy over the existing vertices
    for V in 1 : numvertices(mesh)
        verts[V] = mesh.vertices[V]
    end

    # add a vertex in each edge centroid
    for (E, Edge) in enumerate(cells(edges))
        verts[NV+E] = cartesian(center(chart(edges, Edge)))
    end

    # Build a matrix that given a coarse face gives
    # access to the coarse edges making up its boundary
    D = copy(transpose(connectivity(edges, faces, identity)))
    rows, vals = rowvals(D), nonzeros(D)

    # add four faces in each coarse face
    nf = 4NF
    fcs = zeros(celltype(mesh), nf)

    for F in 1 : numcells(faces)

        # retrieve the indices of the fine vertices in the centroids
        # of the edges in the boundary of coarse face F
        es = zeros(Int,3)
        for i in nzrange(D, F)
            E = rows[i]
            j = abs(vals[i]) # local index of edge E in face F
            @assert 1 <= j <= 3
            e = NV + E
            es[j] = e
        end

        # build the four triangles
        r1 = mesh.faces[F][1]
        r2 = mesh.faces[F][2]
        r3 = mesh.faces[F][3]

        fcs[4(F-1)+1] = index(r1,es[3],es[2])
        fcs[4(F-1)+2] = index(r2,es[1],es[3])
        fcs[4(F-1)+3] = index(r3,es[2],es[1])
        fcs[4(F-1)+4] = index(es[1],es[2],es[3])
    end

    Mesh(verts, fcs)
end

"""
    los_triangle_center(vertices) -> center

Compute the center for a "line-of-sight" refinement.
"""
function los_triangle_center(r)

    n = (r[1]-r[3]) × (r[2]-r[3])
    n /= norm(n)

    l = [
        norm(r[2]-r[3]),
        norm(r[3]-r[1]),
        norm(r[1]-r[2])]

    α = [
        acos((l[2]^2 + l[3]^2 - l[1]^2)/(2*l[2]*l[3])),
        acos((l[3]^2 + l[1]^2 - l[2]^2)/(2*l[3]*l[1])),
        acos((l[1]^2 + l[2]^2 - l[3]^2)/(2*l[1]*l[2]))]

    if α[1] ≥ π/2
        β = 2/3 * sin(α[2]) * sin(α[3]) / cos(α[2] - α[3])
        t = -β * dot(r[2]-r[3], r[3]-r[1]) / dot(n×(r[2]-r[1]), r[3]-r[1])
        r = r[1] + β*(r[2]-r[1]) + t*cross(n, r[2]-r[1])
    elseif α[2] ≥ π/2
        β = 2/3 * sin(α[3]) * sin(α[1]) / cos(α[3] - α[1])
        t = -β * dot(r[3]-r[1], r[1]-r[2]) / dot(n×(r[3]-r[2]), r[1]-r[2])
        r = r[2] + β*(r[3]-r[2]) + t*cross(n, r[3]-r[2])
    elseif α[3] ≥ π/2
        β = 2/3 * sin(α[1]) * sin(α[2]) / cos(α[1] - α[2])
        t = -β * dot(r[1]-r[2], r[2]-r[3]) / dot(n×(r[1]-r[3]), r[2]-r[3])
        r = r[3] + β*(r[1]-r[3]) + t*cross(n, r[1]-r[3])
    else
        r = (r[1]+r[2]+r[3]) / 3
    end
    return r
end

function lineofsight_refinement(mesh::Mesh{U,3}) where U

    D1 = 3

    edges = skeleton(mesh,1)
    faces = skeleton(mesh,2)

    NV, NE, NF = numvertices(mesh), numcells(edges), numcells(faces)

    nv = NV + NE + NF
    verts = Array{vertextype(mesh)}(undef,nv)
    for V in 1 : numvertices(mesh)
        verts[V] = mesh.vertices[V]
    end

    # add a vertex in each face los-centre
    for (F,Face) in enumerate(cells(faces))
        R = vertices(mesh)[Face]
        # @show F
        losc = los_triangle_center(R)
        verts[NV+NE+F] = losc
    end

    # For each edge, find the neighboring triangles, and connect
    # them with a straight line. At a new vertex at the edge crossing
    C = connectivity(edges, faces, identity)
    rows = rowvals(C)
    # vals = nonzeros(EF)
    for (E,Edge) in  enumerate(cells(edges))
        K = nzrange(C,E)
        # Breaks on open meshes
        @assert length(K) == 2
        Fs = rows[K]
        c1 = verts[NV+NE+Fs[1]]
        c2 = verts[NV+NE+Fs[2]]
        e1 = verts[Edge[1]]
        e2 = verts[Edge[2]]
        c1c2 = c2 - c1
        e1e2 = e2 - e1
        c1e1 = e1 - c1
        l12 = norm(c1c2)
        S = c1c2 * c1c2'
        A = one(S) - S / l12^2
        t = -dot(c1e1, A*e1e2) / dot(e1e2, A*e1e2)
        m = e1 + t * e1e2
        verts[NV+E] = m
    end

    # # add a vertex in each edge centroid
    # for (E,Edge) in enumerate(cells(edges))
    #     verts[NV+E] = cartesian(center(chart(edges, Edge)))
    # end

    D = copy(transpose(C))
    rows = rowvals(D)
    vals = nonzeros(D)
    # vals = nonzeros(D)

    # add six faces in each coarse face
    nf = 6NF
    fcs = zeros(celltype(mesh), nf)
    for F in 1 : numcells(faces)
        c = NV + NE + F

        for i in nzrange(D, F)
            E = rows[i]
            e = NV + E

            j = abs(vals[i]) # local index of edge E in face F
            a = mesh.faces[F][mod1(j+1,3)]
            b = mesh.faces[F][mod1(j+2,3)]

            fcs[6(F-1)+2(j-1)+1] = index(a,e,c)
            fcs[6(F-1)+2(j-1)+2] = index(b,c,e)
        end
    end

    D1 = 3
    T = coordtype(mesh)
    M = typeof(mesh)
    BarycentricRefinement{U,D1,T,M}(Mesh(verts, fcs), mesh)
end


text:
export Polygon2,
       PolygonWithHoles2

export area,
       bbox,
       bottom_vertex,
       bounded_side,
       edges,
       has_on_boundary,
       has_on_bounded_side,
       has_on_negative_side,
       has_on_positive_side,
       has_on_unbounded_side,
       holes,
       is_clockwise_oriented,
       is_collinear_oriented,
       is_convex,
       is_counterclockwise_oriented,
       is_simple,
       left_vertex,
       orientation,
       oriented_side,
       outer_boundary,
       right_vertex,
       top_vertex,
       transform,
       vertices

"""
    Polygon2

The type [`Polygon2`](@ref) represents polygons.
"""
Polygon2

"""
    Polygon2()

Creates an empty polygon.
"""
Polygon2()

"""
    Polygon2(polygon::Polygon2)

Copy constructor.
"""
Polygon2(polygon::Polygon2)

"""
    Polygon2(ps::AbstractVector{Point2})
    Polygon2(ps::Point2...)

Creates a polygon with vertices from the sequence defined by the vector `ps`.
"""
Polygon2(ps::AbstractVector{Point2})

Polygon2(ps::AbstractVector) = Polygon2(collect(CxxRef{Point2}, CxxRef.(ps)))
Polygon2(ps::reference_type_union(Point2)...) = Polygon2(collect(CxxRef.(ps)))

"""
    ==(p₁::Polygon2, p₂::Polygon2)

Test for equality: two polygons are equal iff there exists a cyclic permutation
of the vertices of `p₂` such that they are equal to the vertices of `p₁`.
"""
==(p₁::Polygon2, p₂::Polygon2)

"""
    transform(t::AffTransformation2, p::Polygon2)

Returns the image of the polygon `p` under the transformation `t`.
"""
transform(t::AffTransformation2, p::Polygon2)

"""
    push!(p::Polygon2, q::Point2)

Inserts the vertex `q` at the end of polygon `p`.
"""
push!(p::Polygon2, x::Point2)

"""
    empty!(p::Polygon2)

Erases the polygon `p`'s vertices.
"""
empty!(p::Polygon2)

"""
    reverse(p::Polygon2)

Reverses the orientation of the polygon.

The first vertex remains the same.
"""
reverse(p::Polygon2)

"""
    reverse!(p::Polygon2)

In-place version of [`reverse!`](@ref).
"""
reverse!(p::Polygon2)

"""
    vertices(p::Polygon2)

Returns the polygon's vertices.
"""
vertices(p::Polygon2)

"""
    edges(p::Polygon2)

Returns the polygon's edges.
"""
edges(p::Polygon2)

"""
    is_simple(p::Polygon2)

Returns whether `p` is a simple polygon.
"""
is_simple(p::Polygon2)

"""
    is_convex(p::Polygon2)

Returns whether `p` is convex.
"""
is_convex(p::Polygon2)

"""
    orientation(p::Polygon2)

Returns the orientation.

If the number of vertices `size(p) < 3` then [`COLLINEAR`](@ref) is returned.

!!! info "Precondition"

    `is_simple(p)`
"""
orientation(p::Polygon2)

"""
    oriented_side(p::Polygon2, q::Point2)

Returns [`ON_POSITIVE_SIDE`](@ref), [`ON_NEGATIVE_SIDE`](@ref), or
[`ON_ORIENTED_BOUNDARY`](@ref), depending on where point `q` is.

!!! info "Precondition"

    `is_simple(p)`
"""
oriented_side(p::Polygon2, q::Point2)

"""
    bounded_side(p::Polygon2, q::Point2)

Returns [`ON_BOUNDED_SIDE`](@ref), [`ON_BOUNDARY`](@ref), or
[`ON_UNBOUNDED_SIDE`](@ref), depending on where point `q` is.

!!! info "Precondition"

    `is_simple(p)`
"""
bounded_side(p::Polygon2, q::Point2)

"""
    bbox(p::Polygon2)

Returns the smallest bounding box containing this polygon.
"""
bbox(p::Polygon2)

"""
    area(p::Polygon2)

Returns the signed area of the polygon.

This means that the area is positive for counter clockwise polygons and negative
for clockwise polygons.
"""
area(p::Polygon2)

"""
    left_vertex(p::Polygon2)

Returns the leftmost vertex of the polygon with the smallest x-coordinate.
"""
left_vertex(p::Polygon2)

"""
    right_vertex(p::Polygon2)

Returns the rightmost vertex of the polygon with the largest x-coordinate.
"""
right_vertex(p::Polygon2)

"""
    top_vertex(p::Polygon2)

Returns the topmost vertex of the polygon with the largest y-coordinate.
"""
top_vertex(p::Polygon2)

"""
    bottom_vertex(p::Polygon2)

Returns the bottommost vertex of the polygon with the smallest y-coordinate.
"""
bottom_vertex(p::Polygon2)

for name ∈ (:counterclockwise, :clockwise, :collinear)
    fname = Symbol("is_$(name)_oriented")
    NAME = uppercase(string(name))
    @eval begin
        """
            $($fname)(p::Polygon2)

        Returns `orientation(p) == $($NAME)`.
        """
        $fname(p::Polygon2)
    end
end

for (name, es) ∈ zip((:oriented, :bounded),
                     ((:positive_side, :negative_side)
                    , (:boundary, :bounded_side, :unbounded_side)))
    fdel = "$(name)_side"
    for e ∈ es
        CNAME = uppercase("on_$e")
        fname = Symbol("has_on_$e")
        @eval begin
            """
                $($fname)(p::Polygon2, q::Point2)

            Returns `$($fdel)(p, q) == $($CNAME)`.
            """
            $fname(p::Polygon2, q::Point2)
        end
    end
end

"""
    size(p::Polygon2)

Returns the number of vertices of the polygon.
"""
size(p::Polygon2)

"""
    isempty(p::Polygon2)

Returns `size(p) == 0`.
"""
isempty(p::Polygon2)

"""
resize!(p::Polygon2, sz::Integer)

Resizes the polygon.
"""
resize!(p::Polygon2, sz)

"""
    PolygonWithHoles2

The type [`PolygonWithHoles2`](@ref) represents a linear polygon with holes.
"""
PolygonWithHoles2

"""
    PolygonWithHoles2(pgn_boundary::Polygon2)

Constructor from a polygon.
"""
PolygonWithHoles2(pgn_boundary::Polygon2)

"""
    PolygonWithHoles2(pgn_boundary::Polygon2, holes::AbstractVector{Polygon2})
    PolygonWithHoles2(pgn_boundary::Polygon2, holes::Polygon2...)

Constructor from a polygon (outer boundary) and hole polygons.
"""
PolygonWithHoles2(pgn_boundary::Polygon2, holes::AbstractVector{Polygon2})

@cxxdereference PolygonWithHoles2(p::Polygon2, holes::AbstractVector) =
    PolygonWithHoles2(p, collect(CxxRef{Polygon2}, CxxRef.(holes)))
@cxxdereference PolygonWithHoles2(p::Polygon2,
                                  holes::reference_type_union(Polygon2)...) =
    PolygonWithHoles2(p, collect(CxxRef.(holes)))

"""
    is_unbounded(p::PolygonWithHoles2)

Returns `true` if the outer boundary is empty, `false` otherwise.
"""
is_unbounded(p::PolygonWithHoles2)

"""
    outer_boundary(p::PolygonWithHoles2)

Returns the polygon that represents the outer boundary.

Note that this polygon is not necessarily a valid (simple) general polygon
because it may be relatively simple.
"""
outer_boundary(p::PolygonWithHoles2)

"""
    holes(p::PolygonWithHoles2)

Returns the polygon `p`'s hole polygons.
"""
holes(p::PolygonWithHoles2)


text:
"""
    Peters{N,TF,SV,SA} <: AbstractModel

Peter's finite state model with `N` state variables, inputs ``u, \\omega,
\\dot{v}, \\dot{\\omega}`` and parameters ``a, b, a_0, \\alpha_0``
"""
struct Peters{N,TF,TV<:SVector{N,TF},TA<:SMatrix{N,N,TF}} <: AbstractModel
    A::TA
    b::TV
    c::TV
end

# --- Constructors --- #

"""
    Peters{N,TF=Float64}()

Initialize an object of type `Peters` which has `N` aerodynamic degrees of
freedom.
"""
Peters{N}() where N = Peters{N,Float64}()

function Peters{N,TF}() where {N,TF}

    b = zeros(TF, N)
    for n = 1:N-1
        b[n] = (-1)^(n-1)*factorial(big(N + n - 1))/factorial(big(N - n - 1))*
            1/factorial(big(n))^2
    end
    b[N] = (-1)^(N-1)

    c = zeros(TF, N)
    for n = 1:N
        c[n] = 2/n
    end

    d = zeros(TF, N)
    d[1] = 1/2

    D = zeros(TF, N, N)
    for m in 1:N-1
        n = m + 1
        D[n, m] = 1/(2*n)
    end
    for m in 2:N
        n = m - 1
        D[n, m] = -1/(2*n)
    end

    A = D + d*b' + c*d' + 1/2*c*b'

    return Peters(SMatrix{N,N,TF}(A), SVector{N,TF}(b), SVector{N,TF}(c))
end

# --- Traits --- #

number_of_states(::Type{Peters{N,TF,SV,SA}}) where {N,TF,SV,SA} = N
number_of_inputs(::Type{<:Peters}) = 4
number_of_parameters(::Type{<:Peters}) = 6

inplaceness(::Type{<:Peters}) = OutOfPlace()

rate_jacobian_type(::Type{<:Peters}) = Invariant()
state_jacobian_type(::Type{<:Peters}) = Linear()
input_jacobian_type(::Type{<:Peters}) = Nonlinear()
parameter_jacobian_type(::Type{<:Peters}) = Nonlinear()
time_gradient_type(::Type{<:Peters}) = Zeros()

# --- Methods --- #

function get_residual(model::Peters{N,TF,SV,SA}, dx, x, y, p, t) where {N,TF,SV,SA}
    # extract rates
    dλ = SVector{N}(dx)
    # extract states
    λ = SVector{N}(x)
    # extract inputs
    u, ω, vdot, ωdot = y
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    Abar = model.A
    cbar = model.c
    # calculate rates
    return peters_residual(dλ, λ, u, ω, vdot, ωdot, a, b, Abar, cbar)
end

# --- Performance Overloads --- #

function get_rate_jacobian(model::Peters)
    # extract model constants
    Abar = model.A
    # jacobian with respect to aerodynamic states
    return peters_rate_jacobian(Abar)
end

function get_state_jacobian(model::Peters, dx, x, y, p, t)
    # extract inputs
    u, ω, vdot, ωdot = y
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    Abar = model.A
    cbar = model.c
    # jacobian with respect to aerodynamic states
    return peters_state_jacobian(u, b, cbar)
end

function get_input_jacobian(model::Peters{N,TF,SV,SA}, dx, x, y, p, t) where {N,TF,SV,SA}
    # extract states
    λ = SVector{N}(x)
    # extract inputs
    u, ω, vdot, ωdot = y
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    Abar = model.A
    cbar = model.c
    # return jacobian
    return peters_input_jacobian(λ, u, ω, a, b, Abar, cbar)
end

function get_parameter_jacobian(model::Peters{N,TF,SV,SA}, dx, x, y, p, t) where {N,TF,SV,SA}
    # extract states
    λ = SVector{N}(x)
    # extract inputs
    u, ω, vdot, ωdot = y
    # extract parameters
    a, b, a0, α0, cd0, cm0 = p
    # extract model constants
    Abar = model.A
    cbar = model.c
    # return jacobian
    return peters_parameter_jacobian(λ, u, ω, ωdot, a, b, Abar, cbar)
end

# --- Convenience Functions --- #

function set_states!(y, model::Peters; lambda)

    y .= lambda

    return y
end

function set_inputs!(y, model::Peters; u, omega, vdot, omegadot)

    y[1] = u
    y[2] = omega
    y[3] = vdot
    y[4] = omegadot

    return y
end

function set_parameters!(p, model::Peters; a, b, a0, alpha0, cd0, cm0)

    p[1] = a
    p[2] = b
    p[3] = a0
    p[4] = alpha0
    p[5] = cd0
    p[6] = cm0

    return p
end

separate_states(model::Peters, x) = (lambda = x,)

function separate_inputs(model::Peters, y)

    return (u = y[1], omega = y[2], vdot = y[3], omegadot = y[4])
end

function separate_parameters(model::Peters, p)

    return (a = p[1], b = p[2], a0 = p[3], alpha0 = p[4], cd0 = p[5], cm0 = p[6])
end

# --- Internal Methods for Model --- #

function peters_residual(dλ, λ, u, ω, vdot, ωdot, a, b, Abar, cbar)

    return Abar*dλ - cbar*(vdot + u*ω + (b/2-a*b)*ωdot) + u/b*λ
end

peters_rate_jacobian(Abar) = Abar

peters_state_jacobian(u, b, cbar) = u/b*Diagonal(ones(similar_type(cbar)))

function peters_input_jacobian(λ, u, ω, a, b, Abar, cbar)

    dλ_u = -cbar*ω + λ/b
    dλ_ω = -u*cbar
    dλ_vdot = -cbar
    dλ_ωdot = -(b/2-a*b)*cbar

    return hcat(dλ_u, dλ_ω, dλ_vdot, dλ_ωdot)
end

function peters_parameter_jacobian(λ, u, ω, ωdot, a, b, Abar, cbar)

    dλ_a = b*ωdot*cbar
    dλ_b = -cbar*(1/2-a)*ωdot - u/b^2*λ
    dλ_a0 = zero(cbar)
    dλ_α0 = zero(cbar)
    dλ_cd0 = zero(cbar)
    dλ_cm0 = zero(cbar)

    return hcat(dλ_a, dλ_b, dλ_a0, dλ_α0, dλ_cd0, dλ_cm0)
end

# --- Internal Methods for Couplings --- #

function peters_loads(a, b, ρ, a0, α0, cd0, cm0, bbar, u, v, ω, vdot, ωdot, λ)
    # circulatory load factor
    tmp1 = a0*ρ*u*b
    # non-circulatory load factor
    tmp2 = pi*ρ*b^3
    # constant based on geometry
    d = b/2 - a*b
    # induced flow velocity
    λ0 = 1/2 * bbar'*λ
    # normal force at the reference point
    N = tmp1*(v + d*ω - λ0 - u*α0) + tmp2*(vdot/b + u/b*ω - a*ωdot)
    # axial force at the reference point
    A = -a0*ρ*b*(v + d*ω - λ0 - u*α0)^2 + ρ*b*u^2*cd0
    # moment at reference point
    M = -tmp2*(vdot/2 + u*ω + b*(1/8 - a/2)*ωdot) + 2*ρ*b^2*u^2*cm0 + (b/2 + a*b)*N
    # calculate loads
    return SVector(N, A, M)
end

function peters_loads_a(a, b, ρ, a0, α0, bbar, u, v, ω, vdot, ωdot, λ)

    tmp1 = a0*ρ*u*b
    tmp2 = pi*ρ*b^3
    d = b/2 - a*b
    λ0 = 1/2 * bbar'*λ

    N = tmp1*(v + (b/2 - a*b)*ω - λ0 - u*α0) + tmp2*(vdot/b + u/b*ω - a*ωdot)
    N_a = -tmp1*b*ω - tmp2*ωdot
    A_a = 2*a0*ρ*b^2*(v + d*ω - λ0 - u*α0)*ω
    M_a = tmp2*b/2*ωdot + b*N + (b/2 + a*b)*N_a

    return SVector(N_a, A_a, M_a)
end

function peters_loads_b(a, b, ρ, a0, α0, cd0, cm0, bbar, u, v, ω, vdot, ωdot, λ)

    tmp1 = a0*ρ*u*b
    tmp1_b = a0*ρ*u

    tmp2 = pi*ρ*b^3
    tmp2_b = 3*pi*ρ*b^2

    d = b/2 - a*b
    d_b = 1/2 - a

    λ0 = 1/2 * bbar'*λ

    N = tmp1*(v + d*ω - λ0 - u*α0) + tmp2*(vdot/b + u/b*ω - a*ωdot)

    N_b = tmp1_b*(v + d*ω - λ0 - u*α0) + tmp1*d_b*ω +
        tmp2_b*(vdot/b + u/b*ω - a*ωdot) + tmp2*(-vdot/b^2 - u/b^2*ω)
    A_b = -2*a0*ρ*b*(v + d*ω - λ0 - u*α0)*d_b - a0*ρ*(v + d*ω - λ0 - u*α0)^2 +
        ρ*u^2*cd0
    M_b = -tmp2_b*(vdot/2 + u*ω + b*(1/8 - a/2)*ωdot) - tmp2*(1/8 - a/2)*ωdot +
        4*ρ*b*u^2*cm0 + (1/2 + a)*N + (b/2 + a*b)*N_b

    return SVector(N_b, A_b, M_b)
end

function peters_loads_ρ(a, b, ρ, a0, α0, cd0, cm0, bbar, u, v, ω, vdot, ωdot, λ)

    tmp1_ρ = a0*u*b

    tmp2_ρ = pi*b^3

    d = b/2 - a*b

    λ0 = 1/2 * bbar'*λ

    N_ρ = tmp1_ρ*(v + d*ω - λ0 - u*α0) + tmp2_ρ*(vdot/b + u/b*ω - a*ωdot)
    A_ρ = -a0*b*(v + d*ω - λ0 - u*α0)^2 + b*u^2*cd0
    M_ρ = -tmp2_ρ*(vdot/2 + u*ω + (b/8 - a*b/2)*ωdot) + 2*b^2*u^2*cm0 + (b/2 + a*b)*N_ρ

    return SVector(N_ρ, A_ρ, M_ρ)
end

function peters_loads_a0(a, b, ρ, a0, α0, bbar, u, v, ω, λ)

    tmp1_a0 = ρ*u*b

    tmp2 = pi*ρ*b^3

    d = b/2 - a*b

    λ0 = 1/2 * bbar'*λ

    N_a0 = tmp1_a0*(v + d*ω - λ0 - u*α0)
    A_a0 = -ρ*b*(v + d*ω - λ0 - u*α0)^2
    M_a0 = (b/2 + a*b)*N_a0

    return SVector(N_a0, A_a0, M_a0)
end

function peters_loads_α0(a, b, ρ, a0, α0, bbar, u, v, ω, λ)
    tmp1 = a0*ρ*u*b
    d = b/2 - a*b
    λ0 = 1/2 * bbar'*λ
    N_α0 = -tmp1*u
    A_α0 = 2*a0*ρ*b*u*(v + d*ω - λ0 - u*α0)
    M_α0 = (b/2 + a*b)*N_α0
    return SVector(N_α0, A_α0, M_α0)
end

function peters_loads_cd0(b, ρ, u)
    N_cd0 = 0
    A_cd0 = ρ*b*u^2
    M_cd0 = 0
    return SVector(N_cd0, A_cd0, M_cd0)
end

function peters_loads_cm0(b, ρ, u)
    N_cm0 = 0
    A_cm0 = 0
    M_cm0 = 2*ρ*b^2*u^2
    return SVector(N_cm0, A_cm0, M_cm0)
end

function peters_loads_u(a, b, ρ, a0, α0, cd0, cm0, bbar, u, v, ω, λ)
    # circulatory load factor
    tmp1 = a0*ρ*u*b
    tmp1_u = a0*ρ*b
    # non-circulatory load factor
    tmp2 = pi*ρ*b^3
    # constant based on geometry
    d = b/2 - a*b
    # induced flow velocity
    λ0 = 1/2 * bbar'*λ
    # normal force at reference point
    N_u = tmp1_u*(v + d*ω - λ0 - u*α0) - tmp1*α0 + tmp2/b*ω
    # axial force at reference point
    A_u = 2*a0*ρ*b*α0*(v + d*ω - λ0 - u*α0) + 2*ρ*b*u*cd0
    # moment at reference point
    M_u = -tmp2*ω + 4*ρ*b^2*u*cm0 + (b/2 + a*b)*N_u

    return SVector(N_u, A_u, M_u)
end

function peters_loads_v(a, b, ρ, a0, α0, bbar, u, v, ω, λ)
    d = b/2 - a*b
    λ0 = 1/2 * bbar'*λ
    N_v = a0*ρ*u*b
    A_v = -2*a0*ρ*b*(v + d*ω - λ0 - u*α0)
    M_v = (b/2 + a*b)*N_v
    return SVector(N_v, A_v, M_v)
end

function peters_loads_ω(a, b, ρ, a0, α0, bbar, u, v, ω, λ)
    tmp = pi*ρ*b^3
    d = b/2 - a*b
    λ0 = 1/2 * bbar'*λ
    N_ω = a0*ρ*u*b*(b/2 - a*b) + tmp*u/b
    A_ω = 2*a0*ρ*a*b^2*(v + d*ω - λ0 - u*α0)
    M_ω = -tmp*u + (b/2 + a*b)*N_ω
    return SVector(N_ω, A_ω, M_ω)
end

function peters_loads_vdot(a, b, ρ)
    tmp = pi*ρ*b^3
    N_vdot = tmp/b
    A_vdot = 0
    M_vdot = -tmp/2 + (b/2 + a*b)*N_vdot
    return SVector(N_vdot, A_vdot, M_vdot)
end

function peters_loads_ωdot(a, b, ρ)
    tmp = pi*ρ*b^3
    N_ωdot = -tmp*a
    A_ωdot = 0
    M_ωdot = -tmp*(b/8 - a*b/2) + (b/2 + a*b)*N_ωdot
    return SVector(N_ωdot, A_ωdot, M_ωdot)
end

function peters_loads_λ(a, b, ρ, a0, α0, bbar, u, v, ω, λ)
    tmp = a0*ρ*u*b
    d = b/2 - a*b
    λ0 = 1/2 * bbar'*λ
    N_λ = -tmp/2*bbar'
    A_λ = a0*ρ*b*(v + d*ω - λ0 - u*α0)*bbar'
    M_λ = (b/2 + a*b)*N_λ
    return vcat(N_λ, A_λ, M_λ)
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

@testset "bsrsm2" begin
    @testset for elty in [Float32,Float64,Complex64,Complex128]
        @testset "bsrsm2!" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m,n)
            alpha = rand(elty)
            d_X = CudaArray(X)
            d_A = CudaSparseMatrixCSR(sparse(A))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_X = CUSPARSE.bsrsm2!('N','N',alpha,d_A,d_X,'O')
            h_Y = to_host(d_X)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            d_X = CudaArray(rand(elty,n,n))
            @test_throws DimensionMismatch CUSPARSE.bsrsm2!('N','N',alpha,d_A,d_X,'O')
            @test_throws DimensionMismatch CUSPARSE.bsrsm2!('N','T',alpha,d_A,d_X,'O')
            A = sparse(rand(elty,m,n))
            d_A = CudaSparseMatrixCSR(A)
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            @test_throws DimensionMismatch CUSPARSE.bsrsm2!('N','N',alpha,d_A,d_X,'O')
        end

        @testset "bsrsm2" begin
            A = rand(elty,m,m)
            A = triu(A)
            X = rand(elty,m,n)
            alpha = rand(elty)
            d_X = CudaArray(X)
            d_A = CudaSparseMatrixCSR(sparse(A))
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            d_Y = CUSPARSE.bsrsm2('N','N',alpha,d_A,d_X,'O')
            h_Y = to_host(d_Y)
            Y = A\(alpha * X)
            @test Y ≈ h_Y
            A = sparse(rand(elty,m,n))
            d_A = CudaSparseMatrixCSR(A)
            d_A = CUSPARSE.switch2bsr(d_A, convert(Cint,5))
            @test_throws DimensionMismatch CUSPARSE.bsrsm2('N','N',alpha,d_A,d_X,'O')
        end
    end
end


text:
using Primes

function count_prime_cube_pairs(limit)
    cube_limit = Int(ceil((sqrt(12 * limit + 9) - 3)/6))
    prime_list = primes(limit)

    count = 0
    for n in 1:cube_limit - 1
        diff = (n + 1)^3 - n^3

        if diff in prime_list
            count += 1
        end
    end

    println(count)
end

count_prime_cube_pairs(1000000)

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
using MyPkgDemo
using Test

@testset "Tests for add_nums" begin
  @test add_nums(2,3) == 5
  @test add_nums(3,3) == 6
  @test add_nums(2,4) == 6
  @test add_nums(4,4) == 8
  @test add_nums(4,3,2) == 9
  @test add_nums(4,3,4) == 11
end


text:
using Test, KalmanFilter, Random, LinearAlgebra, LazyArrays

Random.seed!(1234)

include("kf.jl")
include("srkf.jl")
include("sigmapoints.jl")
include("augmented_sigmapoints.jl")
include("ukf.jl")
include("srukf.jl")
include("aukf.jl")
include("sraukf.jl")
include("consider.jl")
include("tests.jl")
include("system.jl")


text:
struct Shape{Transform<:Transformation, SDF<:AbstractSignedDistanceField, Mat<:AbstractMaterial}
    # Store the inverse pose since that's what we use to transform input points
    invpose::Transform
    sdf::SDF
    mat::Mat
end

@inline invpose(shape::Shape) = shape.invpose
@inline sdf(shape::Shape) = shape.sdf
@inline mat(shape::Shape) = shape.mat
# Calculate the signed distance within the sdf's coordinate frame by using the inverse pose
@inline function (shape::Shape{Transform, SDF, Mat})(point::Point{3, T}) where
        {T<:Real, Transform<:Transformation, SDF<:AbstractSignedDistanceField, Mat<:AbstractMaterial}
    sdf(shape)(invpose(shape)(point))
end


text:
load("utils/req.jl")
req("utils/utils.jl")

type Unfinished; end
unfinished = Unfinished()

type DictEntry{K,V}
    dict::Dict{K,V}
    key::K
end
DictEntry{K,V}(dict::Dict{K,V}, key) = DictEntry{K,V}(dict, key)

present(entry::DictEntry)   = has(entry.dict, entry.key)
ref(entry::DictEntry)       = entry.dict[entry.key]
assign(entry::DictEntry, x) = (entry.dict[entry.key] = x)

typealias Cache Dict{Function, Dict}

function cacheentry(f::Function, c::Cache, args...)
    DictEntry((@setdefault c[f] Dict()), args)
end

function cached_call(call)
    fullcall = call
    if is_expr(call, doublecolon)
        @expect length(call.args) == 2
        returntype = call.args[2]
        call = call.args[1]
    else
        returntype = Any
    end
    @expect is_expr(call, :call)
    quote
        let entry = cacheentry($call.args...) # args[1] = the function
            local value
            if present(entry)
                value = entry[]
                if is(value, unfinished)
                    error("Reentered evalutation of ( @cached ",
                          ($string(fullcall)), " )\nwith arguments = ",
                          ($expr(:tuple, call.args[2:end])))
                end
                value::($returntype)
            else
                entry[] = unfinished
                entry[] = value = ($fullcall)
                value
            end
        end
    end
end

macro cached(call)
    cached_call(call)
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


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

struct AdminUsergroupsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_adminUsergroupsAddChannels(_api::AdminUsergroupsApi, usergroup_id::String, channel_ids::String; team_id=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.usergroups.addChannels", ["slackAuth"])
    Swagger.set_param(_ctx.form, "usergroup_id", usergroup_id)  # type String
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "channel_ids", channel_ids)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Add one or more default channels to an IDP group.
Param: usergroup_id::String (required)
Param: channel_ids::String (required)
Param: team_id::String
Return: Dict{String, Any}
"""
function adminUsergroupsAddChannels(_api::AdminUsergroupsApi, usergroup_id::String, channel_ids::String; team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsergroupsAddChannels(_api, usergroup_id, channel_ids; team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsergroupsAddChannels(_api::AdminUsergroupsApi, response_stream::Channel, usergroup_id::String, channel_ids::String; team_id=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsergroupsAddChannels(_api, usergroup_id, channel_ids; team_id=team_id, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsergroupsAddTeams(_api::AdminUsergroupsApi, usergroup_id::String, team_ids::String; auto_provision=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.usergroups.addTeams", ["slackAuth"])
    Swagger.set_param(_ctx.form, "usergroup_id", usergroup_id)  # type String
    Swagger.set_param(_ctx.form, "team_ids", team_ids)  # type String
    Swagger.set_param(_ctx.form, "auto_provision", auto_provision)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Associate one or more default workspaces with an organization-wide IDP group.
Param: usergroup_id::String (required)
Param: team_ids::String (required)
Param: auto_provision::Bool
Return: Dict{String, Any}
"""
function adminUsergroupsAddTeams(_api::AdminUsergroupsApi, usergroup_id::String, team_ids::String; auto_provision=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsergroupsAddTeams(_api, usergroup_id, team_ids; auto_provision=auto_provision, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsergroupsAddTeams(_api::AdminUsergroupsApi, response_stream::Channel, usergroup_id::String, team_ids::String; auto_provision=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsergroupsAddTeams(_api, usergroup_id, team_ids; auto_provision=auto_provision, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsergroupsListChannels(_api::AdminUsergroupsApi, usergroup_id::String; team_id=nothing, include_num_members=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.usergroups.listChannels", ["slackAuth"])
    Swagger.set_param(_ctx.query, "usergroup_id", usergroup_id)  # type String
    Swagger.set_param(_ctx.query, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.query, "include_num_members", include_num_members)  # type Bool
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

List the channels linked to an org-level IDP group (user group).
Param: usergroup_id::String (required)
Param: team_id::String
Param: include_num_members::Bool
Return: Dict{String, Any}
"""
function adminUsergroupsListChannels(_api::AdminUsergroupsApi, usergroup_id::String; team_id=nothing, include_num_members=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsergroupsListChannels(_api, usergroup_id; team_id=team_id, include_num_members=include_num_members, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsergroupsListChannels(_api::AdminUsergroupsApi, response_stream::Channel, usergroup_id::String; team_id=nothing, include_num_members=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsergroupsListChannels(_api, usergroup_id; team_id=team_id, include_num_members=include_num_members, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsergroupsRemoveChannels(_api::AdminUsergroupsApi, usergroup_id::String, channel_ids::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.usergroups.removeChannels", ["slackAuth"])
    Swagger.set_param(_ctx.form, "usergroup_id", usergroup_id)  # type String
    Swagger.set_param(_ctx.form, "channel_ids", channel_ids)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Remove one or more default channels from an org-level IDP group (user group).
Param: usergroup_id::String (required)
Param: channel_ids::String (required)
Return: Dict{String, Any}
"""
function adminUsergroupsRemoveChannels(_api::AdminUsergroupsApi, usergroup_id::String, channel_ids::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsergroupsRemoveChannels(_api, usergroup_id, channel_ids; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsergroupsRemoveChannels(_api::AdminUsergroupsApi, response_stream::Channel, usergroup_id::String, channel_ids::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsergroupsRemoveChannels(_api, usergroup_id, channel_ids; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminUsergroupsAddChannels, adminUsergroupsAddTeams, adminUsergroupsListChannels, adminUsergroupsRemoveChannels


text:
using Base.Test
using MPI

# Start workers via `mpiexec` that communicate among themselves via MPI;
# communicate with the workers via TCP
mgr = MPI.MPIManager(np=4)
addprocs(mgr)

refs = []
for w in workers()
    push!(refs, @spawnat w MPI.Comm_rank(MPI.COMM_WORLD))
end
ids = falses(nworkers())
for r in refs
    id = fetch(r)
    @test !ids[id+1]
    ids[id+1] = true
end
for id in ids
    @test id
end

s = @parallel (+) for i in 1:10
    i^2
end
@test s == 385

@mpi_do mgr n = 2
@mpi_do mgr n + 2


text:
## Theory for CAS
using Metatheory
using Metatheory.Library
using Metatheory.EGraphs
using Metatheory.EGraphs.Schedulers
using TermInterface

mult_t = @commutative_monoid (*) 1
plus_t = @commutative_monoid (+) 0

minus_t = @theory begin
    # TODO Jacques Carette's post in zulip chat
    a - a       => 0
    a - b       => a + (-1*b)
    -a          => -1 * a
    a + (-b)    => a + (-1*b)
end

mulplus_t = @theory begin
    # TODO FIXME this rules improves performance and avoids commutative
    # explosion of the egraph
    a + a => 2 * a
    0 * a => 0
    a * 0       => 0
    a * (b + c) == ((a*b) + (a*c))
    a + (b * a) => ((b+1)*a)
end

pow_t = @theory begin
    (y^n) * y   => y^(n+1)
    x^n * x^m   == x^(n+m)
    (x * y)^z   == x^z * y^z
    (x^p)^q     == x^(p*q)
    x^0         => 1
    0^x         => 0
    1^x         => 1
    x^1         => x
    x * x       => x^2
    inv(x)      == x^(-1)
end

div_t = @theory begin
    x / 1 => x
    # x / x => 1 TODO SIGN ANALYSIS
    x / (x / y) => y
    x * (y / x) => y
    x * (y / z) == (x * y) / z
    x^(-1)      == 1 / x
end

trig_t = @theory begin
    sin(θ)^2 + cos(θ)^2     => 1
    sin(θ)^2 - 1            => cos(θ)^2
    cos(θ)^2 - 1            => sin(θ)^2
    tan(θ)^2 - sec(θ)^2     => 1
    tan(θ)^2 + 1            => sec(θ)^2
    sec(θ)^2 - 1            => tan(θ)^2

    cot(θ)^2 - csc(θ)^2     => 1
    cot(θ)^2 + 1            => csc(θ)^2
    csc(θ)^2 - 1            => cot(θ)^2
end

# Dynamic rules
fold_t = @theory begin
    -(a::Number)            |> -a
    a::Number + b::Number   |> a + b
    a::Number * b::Number   |> a * b
    a::Number ^ b::Number   |> begin b < 0 && a isa Int && (a = float(a)) ; a^b end
    a::Number / b::Number   |> a/b
end

using Calculus: differentiate
diff_t = @theory begin
    ∂(y, x::Symbol) |> begin
        z = extract!(_egraph, simplcost; root=y.id)
        @show z
        zd = differentiate(z, x)
        @show zd
        zd
    end
end

cas = fold_t ∪ mult_t ∪ plus_t ∪ minus_t ∪
    mulplus_t ∪ pow_t ∪ div_t ∪ trig_t ∪ diff_t


text:
using Test
import IJulia: Comm, comm_target

target = :notarget
comm_id = "6BA197D8A67A455196279A59EB2FE844"
comm = Comm(target, comm_id, false)
@test :notarget == comm_target(comm)
@test !comm.primary

# comm_info_request in comm_manager.jl
const comms = Dict{String, Comm}(
    "id" => Comm(Symbol("jupyter.widget"), "id", false)
)
msg_content = Dict("target_name" => "jupyter.widget")
reply = if haskey(msg_content, "target_name")
    let t, cb
        t = Symbol(msg_content["target_name"])
        cb(k, v) = comm_target(v) == t # For 0.6
        cb(kv) = comm_target(kv[2]) == t # For 0.7
        filter(cb, comms)
    end
else
    comms
end
@test Dict{String,Comm} == typeof(reply)
_comms = Dict{String,Dict{Symbol,Symbol}}()
for (comm_id,comm) in reply
    local comm_id, comm
    _comms[comm_id] = Dict(:target_name => comm_target(comm))
end
@test Dict("id"=>Dict(:target_name=>Symbol("jupyter.widget"))) == _comms
content = Dict(:comms => _comms)
@test Dict(:comms=>Dict("id"=>Dict(:target_name=>Symbol("jupyter.widget")))) == content


text:
using Pkg
Pkg.add(
    ["Parameters", "CSV", "DataFrames", "DataFramesMeta", "Query", "Printf",
    "Statistics", "Random", "Distributions", "Combinatorics", "Flux", "Tracker"]
    )


text:
module TestIndex
using Test, DataFrames
using DataFrames: Index

i = Index()
push!(i, :A)
push!(i, :B)

inds = Any[1,
           big(1),
           :A,
           [true, false],
           [1],
           [big(1)],
           big(1):big(1),
           [:A],
           Union{Int, Missing}[1],
           Union{BigInt, Missing}[big(1)],
           Union{Symbol, Missing}[:A],
           Any[1],
           Any[:A]]

for ind in inds
    if ind == :A || ndims(ind) == 0
        @test i[ind] == 1
    else
        @test (i[ind] == [1])
    end
end

@test_throws MethodError i[1.0]
@test_throws ArgumentError i[true]
@test_throws ArgumentError i[false]
@test_throws ArgumentError i[Union{Bool, Missing}[true, false]]
@test_throws ArgumentError i[Any[1, missing]]
@test_throws ArgumentError i[[1, missing]]
@test_throws ArgumentError i[[true, missing]]
@test_throws ArgumentError i[Any[true, missing]]
@test_throws MethodError i[[:A, missing]]
@test_throws MethodError i[Any[:A, missing]]
@test_throws ArgumentError i[1.0:1.0]
@test_throws ArgumentError i[[1.0]]
@test_throws ArgumentError i[Any[1.0]]

@test i[1:1] == 1:1

@test_throws BoundsError i[[true]]
@test_throws BoundsError i[true:true]
@test_throws BoundsError i[[true, false, true]]

@test_throws ArgumentError i[["a"]]
@test_throws ArgumentError i[Any["a"]]

@test i[[]] == Int[]
@test i[Int[]] == Int[]
@test i[Symbol[]] == Int[]

@test names(i) == [:A,:B]
@test names!(i, [:a,:a], makeunique=true) == Index([:a,:a_1])
@test_throws ArgumentError names!(i, [:a,:a])
@test names!(i, [:a,:b]) == Index([:a,:b])
@test rename(i, Dict(:a=>:A, :b=>:B)) == Index([:A,:B])
@test rename(i, :a => :A) == Index([:A,:b])
@test rename(i, :a => :a) == Index([:a,:b])
@test rename(i, [:a => :A]) == Index([:A,:b])
@test rename(i, [:a => :a]) == Index([:a,:b])
@test rename(x->Symbol(uppercase(string(x))), i) == Index([:A,:B])
@test rename(x->Symbol(lowercase(string(x))), i) == Index([:a,:b])

@test delete!(i, :a) == Index([:b])
push!(i, :C)
@test delete!(i, 1) == Index([:C])

i = Index([:A, :B, :C, :D, :E])
i2 = copy(i)
names!(i2, reverse(names(i2)))
names!(i2, reverse(names(i2)))
@test names(i2) == names(i)
for name in names(i)
  i2[name] # Issue #715
end

end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
using .Pipeline

# =========================================================================================
"""
    tess_pipeline_lstm_box(
        pipe::TessPipeline,
        filename::AbstractString
    )::Bool

Generate a LSTM BOX file from the pipeline and save it to the specified file.  Returns
`false` if there is a problem adding the LSTM BOX generator to the output.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | pipe     |         | The pipline to collect the LSTM BOX from.
| R | filename |         | The file to write the LSTM BOX to.

__Details:__

If the file exists it will be overwritten.

__Examples:__

```jldoctest
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

tess_pipeline_lstm_box(pipeline, "My Book.box")

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
    true
end

for line in readlines("My Book.box")[1:10]
    println(line)
end

# output

N 11 577 422 591 0
o 11 577 422 591 0
  11 577 422 591 0
o 11 577 422 591 0
n 11 577 422 591 0
e 11 577 422 591 0
  11 577 422 591 0
w 11 577 422 591 0
o 11 577 422 591 0
u 11 577 422 591 0
```

See also: [`tess_run_pipeline`](@ref), [`tess_pipeline_word_box`](@ref),
          [`tess_pipeline_unlv`](@ref)
"""
function tess_pipeline_lstm_box(
            pipe::TessPipeline,
            filename::AbstractString
        )::Bool
    local dir         = dirname(filename)
    local tmpfile, io = mktemp(dir;cleanup=false); close(io)
    local tmpname     = basename(tmpfile)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_lstm_box(tmpfile)

    rm(tmpfile)

    if ptr == C_NULL
        return false
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local task = pipeline_start_task(
        joinpath(dir, "$(tmpname).box"),
        filename
    )

    pipeline_add_renderer(pipe, ptr, task, RENDERER_LSTM_BOX)

    return true
end

# =========================================================================================
"""
    tess_pipeline_lstm_box(
        pipe::TessPipeline
    )::Union{TessOutput{String}, Nothing}

Generate a LSTM BOX file from the pipeline and save it to a string.  Returns `nothing` if
there is a problem adding the LSTM BOX generator to the output.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | pipe  |         | The pipline to collect the LSTM BOX data from.

__Examples:__

```jldoctest
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

book = tess_pipeline_lstm_box(pipeline)

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
end

for line in split(book[], "\\n")[1:10]
    println(line)
end

# output

N 11 577 422 591 0
o 11 577 422 591 0
  11 577 422 591 0
o 11 577 422 591 0
n 11 577 422 591 0
e 11 577 422 591 0
  11 577 422 591 0
w 11 577 422 591 0
o 11 577 422 591 0
u 11 577 422 591 0
```

See also: [`tess_run_pipeline`](@ref), [`tess_pipeline_word_box`](@ref),
          [`tess_pipeline_unlv`](@ref)
"""
function tess_pipeline_lstm_box(
            pipe::TessPipeline
        )::Union{TessOutput{String}, Nothing}
    local path = mktempdir(;cleanup=false)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_lstm_box(
        joinpath(path, "output")
    )

    if ptr == C_NULL
        rm(path;recursive = true)
        return nothing
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local out = TessOutput(String)
    local task = pipeline_start_task(path, "output.box", out)

    pipeline_add_renderer(pipe, ptr, task, RENDERER_LSTM_BOX)

    return out
end

# =========================================================================================
"""
    tess_pipeline_lstm_box(
        func::Function,
        pipe::TessPipeline
    )::Bool

Generate a LSTM BOX file from the pipeline and pass it back to the client via a callback.
Returns `false` if there is a problem adding the LSTM BOX generator to the output.

__Arguments:__

| T | Name  | Default | Description
|:--| :---- | :------ | :----------
| R | func  |         | The function to call with the lines of text.
| R | pipe  |         | The pipline to collect the LSTM BOX data from.

__Details:__

The text will be passed to the caller one line at a time. The "\\n" line terminator will be
included with the text.

__Examples:__

```jldoctest
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

count = 0
tess_pipeline_lstm_box(pipeline) do line
    global count
    if count < 10
        print(line)
    end
    count += 1
end

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
end

# output

N 11 577 422 591 0
o 11 577 422 591 0
  11 577 422 591 0
o 11 577 422 591 0
n 11 577 422 591 0
e 11 577 422 591 0
  11 577 422 591 0
w 11 577 422 591 0
o 11 577 422 591 0
u 11 577 422 591 0
true
```

See also: [`tess_run_pipeline`](@ref), [`tess_pipeline_word_box`](@ref),
          [`tess_pipeline_unlv`](@ref)
"""
function tess_pipeline_lstm_box(
            func::Function,
            pipe::TessPipeline
        )::Bool
    local path = mktempdir(;cleanup=false)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_lstm_box(
        joinpath(path, "output")
    )

    if ptr == C_NULL
        rm(path;recursive = true)
        return false
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local task = pipeline_start_task(path, "output.box", func, "")

    pipeline_add_renderer(pipe, ptr, task, RENDERER_LSTM_BOX)

    return true
end


text:
x = 3
y = 2

println(2x)

# same as y / x
println(x \ y)

# pipe
println([1,2,3] |> x -> x * 10 |> sum)

# like lisp
println(*(10, 2, 20))

# like natural formula
println(2(10 * 10)x)

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
include("common.jl")
@testset "Unconstrained time step" begin
    # Check that numerical solutions approximate analytical solutions,
    # independent of problem structure
    @testset "standard history" begin
        # standard algorithm
        alg = MethodOfSteps(BS3(); constrained=false)

        # different tolerances
        alg1 = MethodOfSteps(Tsit5(), constrained=false, max_fixedpoint_iters=100,
                             fixedpoint_abstol=1e-12, fixedpoint_reltol=1e-12)
        alg2 = MethodOfSteps(DP8(), constrained=false, max_fixedpoint_iters=10,
                             fixedpoint_abstol=1e-8, fixedpoint_reltol=1e-10)
        alg3 = MethodOfSteps(Tsit5(), constrained=true, max_fixedpoint_iters=10,
                             fixedpoint_abstol=1e-8, fixedpoint_reltol=1e-10)
        alg4 = MethodOfSteps(DP5(), constrained=false, max_fixedpoint_iters=100,
                             fixedpoint_abstol=1e-12, fixedpoint_reltol=1e-12)

        ## Single constant delay
        @testset "single constant delay" begin
            @testset "short time span" begin
                ### Not in-place function with scalar history function
                prob = prob_dde_1delay_scalar_notinplace
                sol = solve(prob, alg)

                @test sol.errors[:l∞] < 3.0e-5
                @test sol.errors[:final] < 2.1e-5
                @test sol.errors[:l2] < 1.2e-5

                ### Not in-place function with vectorized history function
                prob = prob_dde_1delay_notinplace
                sol2 = solve(prob, alg)

                @test sol.t ≈ sol2.t && sol.u ≈ sol2[1, :]

                ### In-place function
                prob = prob_dde_1delay
                sol2 = solve(prob, alg)

                @test sol.t ≈ sol2.t && sol.u ≈ sol2[1, :]
            end

            @testset "long time span" begin
                prob = prob_dde_1delay_long_scalar_notinplace

                sol1 = solve(prob, alg1)
                sol2 = solve(prob, alg2)
                sol3 = solve(prob, alg3)
                sol4 = solve(prob, alg4)

                @test abs(sol1[end] - sol2[end]) < 5.3e-4
                @test abs(sol1[end] - sol3[end]) < 1.1e-8
                @test abs(sol1[end] - sol4[end]) < 2.9e-4
            end
        end

        ## Two constant delays
        @testset "two constant delays" begin
            @testset "short time span" begin
                ### Not in-place function with scalar history function
                prob = prob_dde_2delays_scalar_notinplace
                sol = solve(prob, alg)

                @test sol.errors[:l∞] < 1.9e-6
                @test sol.errors[:final] < 1.2e-6
                @test sol.errors[:l2] < 1.0e-6

                ### Not in-place function with vectorized history function
                prob = prob_dde_2delays_notinplace
                sol2 = solve(prob, alg)

                @test sol.t ≈ sol2.t && sol.u ≈ sol2[1, :]

                ### In-place function
                prob = prob_dde_2delays
                sol2 = solve(prob, alg)

                @test sol.t ≈ sol2.t && sol.u ≈ sol2[1, :]
            end

            @testset "long time span" begin
                prob = prob_dde_2delays_long_scalar_notinplace

                sol1 = solve(prob, alg1)
                # sol2 = solve(prob, alg) # aborted because of dt <= dtmin
                sol3 = solve(prob, alg3)
                sol4 = solve(prob, alg4)

                # relaxed tests to prevent floating point issues
                @test abs(sol1[end] - sol3[end]) < 7.9e-13 # 7.9e-15
                @test abs(sol1[end] - sol4[end]) < 1.6e-12 # 1.6e-14
            end
        end
    end

    ## Non-standard history functions
    @testset "non-standard history" begin
        alg = MethodOfSteps(Tsit5(), constrained=false, max_fixedpoint_iters=100,
                            fixedpoint_abstol=1e-12, fixedpoint_reltol=1e-12)

        @testset "idxs" begin
            function f(du,u,h,p,t)
                du[1] = -h(p, t-0.2;idxs=1) + u[1]
            end
            h(p, t; idxs=nothing) = typeof(idxs) <: Number ? 0.0 : [0.0]

            prob = DDEProblem(f, [1.0], h, (0.0, 100.), constant_lags = [0.2])
            solve(prob, alg)
        end

        @testset "in-place" begin
            function f(du,u,h,p,t)
                h(du, p, t-0.2)
                du[1] = -du[1] + u[1]
            end
            h(val, p, t) = (val .= 0.0)

            prob = DDEProblem(f, [1.0], h, (0.0, 100.0), constant_lags = [0.2])
            solve(prob, alg)
        end
    end
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using CUDA
using StaticArrays

# %%
θ = range(0, 2π; length=361)[1:end-1]
θ_cu = cu(collect(θ))

# %%
r_smf32(t) = @SMatrix(Float32[
    cos(t) -sin(t) 0
    sin(t)  cos(t) 0
    0      0       1
])

f(t) = r_smf32(t) * @SVector Float32[1, 0, 0]

# %%
f.(θ_cu)

# %%
using CUDA, NNlib, NNlibCUDA

# %%
?batched_vec

# %%
R(t) = [
    cos(t) -sin(t) 0
    sin(t)  cos(t) 0
    0      0       1
]

θ = range(0, 2π; length=361)[1:end-1]
A = mapslices(t -> R(t[1, 1]), reshape(θ, 1, 1, :); dims=(1, 2))
A_cu = cu(A)

# %%
batched_vec(A_cu, cu([1.0, 0.0, 0.0]))

# %%
using CUDA, NNlib, NNlibCUDA

R(t) = [
    cos(t) -sin(t) 0
    sin(t)  cos(t) 0
    0      0       1
]

θ = range(0, 2π; length=361)[1:end-1]
A = mapslices(t -> R(t[1, 1]), reshape(θ, 1, 1, :); dims=(1, 2))
batched_vec(cu(A), cu([1.0, 0.0, 0.0]))

# %% [markdown]
# See also https://discourse.julialang.org/t/how-to-broadcast-or-batch-multiply-a-batch-of-matrices-with-another-matrix-on-the-gpu/67259

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
using Test, Statistics

pdir = joinpath(@__DIR__, "..", "examples", "m10.4s")
cd(pdir) do

  isdir("tmp") && rm("tmp", recursive=true);

  include(joinpath(pdir, "m10.4s.jl"))

  @test 0.8 < mean(m10_4s[:bp].value) < 0.9

  isdir("tmp") && rm("tmp", recursive=true);

end


text:
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

# rep   = parse(Int,ARGS[1])
# ip    = parse(Int,ARGS[2])
# dir = ARGS[3]

alpha = 0.1:0.1:1
pArr = [50, 100, 150]
n = 300
NUM_REP = 1000

bnm = Binomial(1, 0.05)

for ip=1:3
  for ialpha=1:10

    p = pArr[ip]
    Random.seed!(54298)

    # generate X
    Ωx = Matrix{Float64}(I, p, p)
    for i = 1:(p-1)
      for j = (i+1):p
        Ωx[i, j] = rand(bnm) * 0.8
        Ωx[j, i] = Ωx[i, j]
      end
    end
    d = Vector{Float64}(undef, p)
    rand!(Uniform(0.5, 2.5), d)
    d .= sqrt.(d)
    D = Diagonal(d)
    Ωx = D * Ωx * D

    level = maximum( diag(Ωx) ) * sqrt(2. * log(p) / n)

    bc = Bernoulli()
    Ωy = copy(Ωx)
    rd = sample(1:p-1, 8; replace=false)
    sort!(rd)
    Ωy[rd[1], rd[5]] += (2 * rand(bc) - 1) * level
    Ωy[rd[2], rd[6]] += (2 * rand(bc) - 1) * level
    Ωy[rd[3], rd[7]] += (2 * rand(bc) - 1) * level
    Ωy[rd[4], rd[8]] += (2 * rand(bc) - 1) * level

    Ωy[rd[5], rd[1]] = Ωy[rd[1], rd[5]]
    Ωy[rd[6], rd[2]] = Ωy[rd[2], rd[6]]
    Ωy[rd[7], rd[3]] = Ωy[rd[3], rd[7]]
    Ωy[rd[8], rd[4]] = Ωy[rd[4], rd[8]]

    de = max( abs(minimum(eigvals(Ωx))), abs(minimum(eigvals(Ωy))) ) + 0.05
    Ωx = (Ωx + de*I)/(1+de)
    Ωy = (Ωy + de*I)/(1+de)

    Ωy = (1. - alpha[ialpha]) * Ωx + alpha[ialpha] * Ωy
    Δ = Ωx - Ωy

    Σx = inv(Symmetric(Ωx))
    Σy = inv(Symmetric(Ωy))
    dist_X = MvNormal(convert(Matrix, Σx))
    dist_Y = MvNormal(convert(Matrix, Σy))

    rp = div(p * (p+1), 2)
    res = Matrix{DiffPrecResultNormal}(undef, rp, NUM_REP)

    @show p, ialpha
    @time for rep=1:NUM_REP

      # generate data
      Random.seed!(12346 + rep*100 + ialpha)
      X = rand(dist_X, n)'
      Y = rand(dist_Y, n)'

      Sx = Symmetric( X'X / n )
      Sy = Symmetric( Y'Y / n )

      it = 0
      for col=1:p
        for row=col:p
          it = it + 1
          res[it, rep] = DiffPrecTest.estimate(SeparateNormal(), X, Y, row, col; Sx=Sx, Sy=Sy)
        end
      end

    end

    @save "YinXia_power_res_$(ip)_$(ialpha).jld" res
  end
end


text:
# Simple Quadratic Regression
# Originally written by Keegan Go for lsqpy
# Translated into LinearLeastSquares.jl by Karanveer Mohan and David Zeng

using LinearLeastSquares
using Gadfly
# Set the random seed to get consistent data
srand(1)

# Number of examples to use
n = 100

# Specify the true value of the variable
true_coeffs = [2; -2; 0.5]

# Generate data
x_data = rand(n, 1) * 5
x_data_expanded = hcat([x_data .^ i for i in 1 : 3]...)
y_data = x_data_expanded * true_coeffs + 0.5 * rand(n, 1)

quadratic_coeff = Variable()
slope = Variable()
offset = Variable()
quadratic = offset + x_data * slope + quadratic_coeff * x_data .^ 2
residuals = quadratic - y_data
fit_error = sum_squares(residuals)
optval = minimize!(fit_error)

# Create some evenly spaced points for plotting, again replicate powers
t = reshape([0 : 0.1 : 5], length([0 : 0.1 : 5]), 1)
t_squared = t .^ 2

# Plot our regressed function
p = plot(
  layer(x=x_data, y=y_data, Geom.point),
  layer(x=t, y=evaluate(offset) + t * evaluate(slope) + t_squared * evaluate(quadratic_coeff), Geom.line),
  Theme(panel_fill=color("white"))
)
# draw(PNG("quadratic_regression.png", 16cm, 12cm), p)


text:
module Day13

export day13

using CombinedParsers
import CombinedParsers.Regexp: word, whitespace_maybe

function parseinput(filename)
    input = Vector{String}()
    for line in eachline(filename)
        push!(input, line)
    end
    earliest = parse(Int, input[1])
    buses = join(Repeat(Either('x', Numeric(Int))), ',')
    (earliest, input[2] |> buses)
end

# Bus with timestamp t arrives at dt = l*t - earliest.
# If earliest = k*t + r, then dt = (k+1)*t - earliest
# = t - r = t - mod(earliest, t)
function day13(::Val{:a}, filename)
    (earliest, buses::Vector{Union{Char,Int}}) = parseinput(filename)
    buses = filter(buses) do x
        isa(x, Int)
    end
    later = map(buses) do t
        r = mod(earliest, t)
        r == 0 ? 0 : t - r
    end
    (dt,ind) = findmin(later)
    dt * buses[ind]
end

# The time t is such that t = -i+1 mod n_i, where
# n_i is the timestamp of the i-th bus in the list
# including the 'x'. We solve using the CRT.
function day13(::Val{:b}, filename)
    (earliest, buses::Vector{Union{Char,Int}}) = parseinput(filename)
    N = prod(filter(c -> isa(c, Int), buses))
    suma = 0
    for (i,m) in enumerate(buses)
        if isa(m, Int)
            n = div(N,m)
            suma = mod(suma + (n * (-i+1) * gcdx(n,m)[2]), N)
        end
    end
    suma
end

end


text:
using Bijectors
using Random

rng = MersenneTwister(123456)

@testset "parameter wrapper" begin
    @testset "bijector" begin
        # scalar
        l = rand(rng)
        # vector
        N = 5
        v = rand(rng, N)
        # matrix
        M, M′ = 4, 7
        W = rand(rng, M, M′)

        # identity bijector
        pl = Parameter(l)
        @test l == value(pl)
        @test length(pl) == 1 == num_of_param(pl)
        @test size(pl) == ()
        pv = Parameter(v)
        @test v == value(pv)
        @test length(pv) == N == num_of_param(pv)
        @test size(pv) == (N,)
        pW = Parameter(W)
        @test W == value(pW)
        @test length(pW) == M*M′ == num_of_param(pW)
        @test size(pW) == (M, M′)

        # Exp/Log bijector
        pl_pos = Parameter(l, Val(:pos))
        @test first(pl_pos.x) == Bijectors.Log{0}()(l)
        @test value(pl_pos) == Bijectors.Exp{0}()(Bijectors.Log{0}()(l))
        pv_pos = Parameter(v, Val(:pos))
        @test pv_pos.x == Bijectors.Log{1}()(v)
        @test value(pv_pos) == Bijectors.Exp{1}()(Bijectors.Log{1}()(v))
        pW_pos = Parameter(W, Val(:pos))
        @test pW_pos.x == Bijectors.Log{2}()(W)
        @test value(pW_pos) == Bijectors.Exp{2}()(Bijectors.Log{2}()(W))

        # check input domain
        @test_throws DomainError Parameter(-3.0, Val(:pos))
        @test_throws DomainError Parameter(9.0, inv(Bijectors.Logit(1.3, 6.0)))
    end
end

@testset "constant" begin
    n, m = 3, 5
    V = randn(rng, n)
    M = randn(rng, n, m)
    cV = Constant(V)
    cM = Constant(M)
    @test length(cV) == n
    @test num_of_param(cV) == 0
    @test size(cV) == (n,)
    @test length(cM) == n*m
    @test num_of_param(cM) == 0
    @test size(cM) == (n, m)
    @test value(cV) == V
    @test value(cM) == M
end

@testset "multiple parameters" begin
    π_num = float(π)
    r = rand(rng)
    θ = rand(rng)*π_num
    ϕ = rand(rng)*2*π_num

    sph_coord = Parameter([r, θ, ϕ], Bijectors.Exp{0}(),
                          inv(Bijectors.Logit(0.0, π_num)),
                          inv(Bijectors.Logit(0.0, 2*π_num)))
    @test length(sph_coord) == 3
    @test num_of_param(sph_coord) == 3
    @test size(sph_coord) == (3,)
    @test value(sph_coord) ≈ [r, θ, ϕ]
end


text:
using Revise
using DataFrames, CSV, JDF
using WeakRefStrings

a[!, :stringarr] = StringArray(rand(["a", "a", "b"], size(a,1)))
a[!, :cate] = categorical(a[!, :stringarr])

@time a = CSV.read("c:/data/feature_matrix_cleaned.csv");

@time savejdf("c:/data/feature_matrix_cleaned.csv.jdf", a)

a = nothing
@time a = loadjdf("c:/data/feature_matrix_cleaned.csv.jdf")

type_compress!(a, compress_float=true)
@time savejdf("c:/data/feature_matrix_cleaned.csv.compressed.jdf", a)

using BenchmarkTools
@benchmark a = loadjdf("c:/data/feature_matrix_cleaned.csv.jdf")

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
using Plots, Test

@testset "Subplot sclicing" begin
    pl = @test_nowarn plot(
        rand(4, 8),
        layout = 4,
        yscale = [:identity :identity :log10 :log10],
    )
    @test pl[1][:yaxis][:scale] == :identity
    @test pl[2][:yaxis][:scale] == :identity
    @test pl[3][:yaxis][:scale] == :log10
    @test pl[4][:yaxis][:scale] == :log10
end

@testset "Plot title" begin
    pl = plot(rand(4, 8), layout = 4, plot_title = "My title")
    @test pl[:plot_title] == "My title"
    @test pl[:plot_titleindex] == 5

    plot!(pl)
    @test pl[:plot_title] == "My title"
    @test pl[:plot_titleindex] == 5

    plot!(pl, plot_title = "My new title")
    @test pl[:plot_title] == "My new title"
    @test pl[:plot_titleindex] == 5
end

@testset "Plots.jl/issues/4083" begin
    p = plot(plot(1:2), plot(1:2); border = :grid, plot_title = "abc")
    @test p[1][:framestyle] === :grid
    @test p[2][:framestyle] === :grid
    @test p[3][:framestyle] === :none
end

@testset "Coverage" begin
    p = plot((plot(i) for i in 1:4)..., layout = (2, 2))

    sp = p[end]
    @test sp isa Plots.Subplot
    @test size(sp) == (1, 1)
    @test length(sp) == 1
    @test sp[1, 1] == sp
    @test Plots.get_subplot(p, UInt32(4)) == sp
    @test Plots.series_list(sp) |> first |> Plots.get_subplot isa Plots.Subplot
    @test Plots.get_subplot(p, keys(p.spmap) |> first) isa Plots.Subplot

    gl = p[2, 2]
    @test gl isa Plots.GridLayout
    @test length(gl) == 1
    @test size(gl) == (1, 1)
    @test Plots.layout_args(gl) == (gl, 1)

    @test size(p, 1) == 2
    @test size(p, 2) == 2
    @test size(p) === (2, 2)
    @test ndims(p) == 2

    @test p[1][end] isa Plots.Series
    show(devnull, p[1])

    @test Plots.getplot(p) == p
    @test Plots.getattr(p) == p.attr
    @test Plots.backend_object(p) == p.o
    @test occursin("Plot", string(p))
    print(devnull, p)

    @test Plots.to_pixels(1Plots.mm) isa AbstractFloat
    @test Plots.ispositive(1Plots.mm)
    @test size(Plots.defaultbox) == (0Plots.mm, 0Plots.mm)
    show(devnull, Plots.defaultbox)
    show(devnull, p.layout)

    @test Plots.make_measure_hor(1Plots.mm) == 1Plots.mm
    @test Plots.make_measure_vert(1Plots.mm) == 1Plots.mm

    @test Plots.parent(p.layout) isa Plots.RootLayout
    show(devnull, Plots.parent_bbox(p.layout))

    rl = Plots.RootLayout()
    show(devnull, rl)
    @test parent(rl) === nothing
    @test Plots.parent_bbox(rl) == Plots.defaultbox
    @test Plots.bbox(rl) == Plots.defaultbox

    el = Plots.EmptyLayout()
    @test Plots.update_position!(el) === nothing
    @test size(el) == (0, 0)
    @test length(el) == 0
    @test el[1, 1] === nothing

    @test Plots.left(el) == 0Plots.mm
    @test Plots.top(el) == 0Plots.mm
    @test Plots.right(el) == 0Plots.mm
    @test Plots.bottom(el) == 0Plots.mm

    @test_throws ErrorException Plots.layout_args(nothing)
end


text:
export SPSystem

"""
    SPSystem(polynomials, vars) <: AbstractSystem

Create a system using the [`StaticPolynomials`](https://github.com/JuliaAlgebra/StaticPolynomials.jl) package.
Note that `StaticPolynomials` leverages Julias metaprogramming capabilities to automatically generate
functions to evaluate the system and its Jacobian. These generated functions are *very fast* but
at the cost of possibly large compile times. The compile time depends on the size of the support of the polynomial system.
If you intend to solve a large system or you need to solve a system
with the *same support* but different coefficients even large compile times can be worthwile.
As a general rule of thumb this usually is twice as fast as solving the same system using [`FPSystem`](@ref).

## Example
You can use `SPSystem` as follows with solve
```julia
@polyvar x y
F = [x^2+3y^4-2, 2y^2+3x*y+4]
solve(F, system=SPSystem)
```
"""
struct SPSystem{S<:SP.PolynomialSystem} <: AbstractSystem
    system::S
end

SPSystem(polys::Vector{<:MP.AbstractPolynomial}, vars) = SPSystem(SP.PolynomialSystem(polys, variables=vars))
SPSystem(polys::Vector{<:MP.AbstractPolynomial}; kwargs...) = SPSystem(SP.PolynomialSystem(polys; kwargs...))
function SPSystem(exponents::Vector{<:Matrix}, coeffs::Vector{<:Vector})
    n = length(exponents)
    @assert n == length(coeffs)
    SPSystem(SP.PolynomialSystem((SP.Polynomial(coeffs[i], exponents[i]) for i in 1:n)...))
end

Base.size(F::SPSystem) = (SP.npolynomials(F.system), SP.nvariables(F.system))

cache(F::SPSystem, x, p=nothing) = SystemNullCache()
Base.@propagate_inbounds evaluate!(u, F::SPSystem, x, ::SystemNullCache) = SP.evaluate!(u, F.system, x)
Base.@propagate_inbounds evaluate!(u, F::SPSystem, x, p, ::SystemNullCache) = SP.evaluate!(u, F.system, x, p)
evaluate(F::SPSystem, x, ::SystemNullCache) = SP.evaluate(F.system, x)
evaluate(F::SPSystem, x, p, ::SystemNullCache) = SP.evaluate(F.system, x, p)
Base.@propagate_inbounds jacobian!(U, F::SPSystem, x, ::SystemNullCache) = SP.jacobian!(U, F.system, x)
Base.@propagate_inbounds jacobian!(U, F::SPSystem, x, p, ::SystemNullCache) = SP.jacobian!(U, F.system, x, p)
jacobian(F::SPSystem, x, ::SystemNullCache) = SP.jacobian(F.system, x)
jacobian(F::SPSystem, x, p, ::SystemNullCache) = SP.jacobian(F.system, x, p)
Base.@propagate_inbounds evaluate_and_jacobian!(u, U, F::SPSystem, x, ::SystemNullCache) = SP.evaluate_and_jacobian!(u, U, F.system, x)
Base.@propagate_inbounds evaluate_and_jacobian!(u, U, F::SPSystem, x, p, ::SystemNullCache) = SP.evaluate_and_jacobian!(u, U, F.system, x, p)
evaluate_and_jacobian(F::SPSystem, x, ::SystemNullCache) = SP.evaluate_and_jacobian(F.system, x)
evaluate_and_jacobian(F::SPSystem, x, p, ::SystemNullCache) = SP.evaluate_and_jacobian(F.system, x, p)
Base.@propagate_inbounds differentiate_parameters!(U, F::SPSystem, x, p, ::SystemNullCache) = SP.differentiate_parameters!(U, F.system, x, p)
differentiate_parameters(F::SPSystem, x, p, ::SystemNullCache) = SP.differentiate_parameters(F.system, x, p)

(F::SPSystem)(x...) = evaluate(F, x...)

"""
    set_coefficients!(F::SPSystem, coeffs::Vector{<:Vector})

Set the coefficients of `F` to the provided coefficients
"""
@generated function set_coefficients!(G::SP.PolynomialSystem{N}, coefficients::AbstractVector{<:AbstractVector}) where {N}
    quote
        $((:(copyto!(G.polys[$i].coefficients, coefficients[$i])) for i=1:N)...)
        G
    end
end
function set_coefficients!(G::SPSystem, coefficients::AbstractVector{<:AbstractVector})
    set_coefficients!(G.system, coefficients)
    G
end


text:
"""
    load_bert_pretrain(path, sym = :all)

Loading bert data/model from pretrain bson data. use `sym` to determine which data to load.
possible value: `:all`(default), `bert_model`, `wordpiece`, `tokenizer`.
"""
function load_bert_pretrain(path::AbstractString, sym = :all)
  @info "loading pretrain bert model: $(basename(path)) $(sym == :all ? "" : sym)"
  @assert sym ∈ (:all, :bert_model, :wordpiece, :tokenizer) "sym only support :all, :bert_model, :wordpiece, :tokenizer"
  if istfbson(path)

    bson = BSON.parse(path)

    config = BSON.raise_recursive(bson[:config], Main)
    tokenizer = named2tokenizer(config["filename"])

    sym == :tokenizer && return tokenizer

    wordpiece = WordPiece(convert(Vector{String}, BSON.raise_recursive(bson[:vocab], Main)))

    sym == :wordpiece && return wordpiece

    weights = BSON.raise_recursive(bson[:weights], Main)
    bert_model = load_bert_from_tfbson(config, weights)

    sym == :bert_model && return bert_model

    return bert_model, wordpiece, tokenizer
  elseif isbson(path)
    bson = BSON.parse(path)
    tokenizer = BSON.raise_recursive(bson[:tokenizer], Main)
    sym == :tokenizer && return tokenizer
    wordpiece = BSON.raise_recursive(bson[:wordpiece], Main)
    sym == :wordpiece && return wordpiece
    bert_model = BSON.raise_recursive(bson[:bert_model], Main)
    sym == :bert_model && return bert_model
    return bert_model, wordpiece, tokenizer
  else
    error("""
          Unknown data type: Is it a bert model?
          If it is the weight file get from google-bert, then run `tfckpt2bson` on that file beforehand.
    """)
  end
end

# function load_pretrain_file(path::AbstractString)
#   if istfbson(path)
#     @info "loading raw data: $path"
#     bson = BSON.load(path)
#     bson[:config], bson[:weights], bson[:vocab]
#   elseif isbson(path)
#     @info "loading model"
#     bson = BSON.load(path)
#     bson[:bert_model], bson[:wordpiece], bson[:tokenizer]
#   elseif iszip(path)
#     data = ZipFile.Reader(path)
#     config, weights, vocab = readckptfolder(data)
#   elseif isdir(path)
#     config, weights, vocab = readckptfolder(path)
#   else
#     error("""Unknown file format""")
#   end
# end


text:
const YEAR  = 365.242
const GNEWT = 39.4845/YEAR^2  # Units of MSUN*AU^3/YEAR^2
const NDIM  = 3
const third = 1.0/3.0

include("../src/init_nbody.jl")

#@testset "init_nbody" begin

elements = readdlm("elements.txt",',')

#elements[2:3,1] /= 100.0
#n_body = 4
n_body = 3
t0 = 7257.93115525-7300.0
elements[:,3] .-= 7300.0

jac_init     = zeros(Float64,7*n_body,7*n_body)
jac_init_big = zeros(BigFloat,7*n_body,7*n_body)
jac_init_num = zeros(BigFloat,7*n_body,7*n_body)
x,v = init_nbody(elements,t0,n_body,jac_init)
elements_big = big.(elements)
t0big = big(t0)
xbig,vbig = init_nbody(elements_big,t0big,n_body,jac_init_big)
elements0 = copy(elements)
#dq = big.([1e-10,1e-5,1e-6,1e-6,1e-6,1e-5,1e-5])
#dq = big.([1e-10,1e-8,1e-8,1e-8,1e-8,1e-8,1e-8])
dq = big.([1e-12,1e-10,1e-10,1e-10,1e-10,1e-10,1e-10])

#dq = big.([1e-15,1e-15,1e-15,1e-15,1e-15,1e-15,1e-15])
t0big = big(t0)
# Now, compute derivatives numerically:
for j=1:n_body
  for k=1:7
    elementsbig = big.(elements0)
    dq0 = dq[k]; if j==1 && k==1 ; dq0 = big(1e-15); end
    elementsbig[j,k] += dq0
    xp,vp = init_nbody(elementsbig,t0big,n_body)
    elementsbig[j,k] -= 2dq0
    xm,vm = init_nbody(elementsbig,t0big,n_body)
    for l=1:n_body, p=1:3
      i1 = (l-1)*7+p
      if k == 1
        j1 = j*7
      else
        j1 = (j-1)*7+k-1
      end
      jac_init_num[i1,  j1] = (xp[p,l]-xm[p,l])/dq0*.5
      jac1 = jac_init[i1,j1]; jac2 = jac_init_num[i1,j1]
      if abs(jac1-jac2) > 1e-4*abs(jac1+jac2) && abs(jac1+jac2) > 1e-14
        println(l," ",p," ",j," ",k," ",jac_init_num[i1,j1]," ",jac_init[i1,j1]," ",jac_init_num[i1,j1]/jac_init[i1,j1])
      end
      jac_init_num[i1+3,j1] = (vp[p,l]-vm[p,l])/dq0*.5
      jac1 = jac_init[i1+3,j1]; jac2 = jac_init_num[i1+3,j1]
      if abs(jac1-jac2) > 1e-4*abs(jac1+jac2) && abs(jac1+jac2) > 1e-14
        println(l," ",p+3," ",j," ",k," ",jac_init_num[i1+3,j1]," ",jac_init[i1+3,j1]," ",jac_init_num[i1+3,j1]/jac_init[i1+3,j1])
      end
    end
  end
  jac_init_num[j*7,j*7]=1.0
end

#println("Maximum jac_init-jac_init_num: ",maximum(abs.(jac_init-jac_init_num)))
println("Maximum jac_init-jac_init_num: ",maximum(abs.(asinh.(jac_init)-asinh.(jac_init_num))))
println("Maximum jac_init-jac_init_big: ",maximum(abs.(asinh.(jac_init)-asinh.(jac_init_big))))
println("Maximum jac_init/jac_init_big-1: ",maximum(abs.(jac_init[jac_init_big .!= 0.0]./jac_init_big[jac_init_big .!= 0.0].-1)))
println("Maximum jac_init_big-jac_init_num: ",maximum(abs.(asinh.(jac_init_num)-asinh.(jac_init_big))))
println(convert(Array{Float64,2},jac_init_big-jac_init_num))
#@test isapprox(jac_init_num,jac_init)
@test isapprox(jac_init_num,jac_init;norm=maxabs)
@test isapprox(jac_init,convert(Array{Float64,2},jac_init_big);norm=maxabs)
#end

for i=1:21, j=1:21
  if jac_init_big[i,j] != 0.0
    dij = abs(jac_init[i,j]/jac_init_big[i,j]-1)
    if dij > 1e-12
      jac_max = dij
      println(i," ",j," ",convert(Float64,jac_max)," ",jac_init[i,j]," ",convert(Float64,jac_init_big[i,j])," ",convert(Float64,jac_init[i,j]-jac_init_big[i,j]))
    end
  end
end


text:
module TestPlots

using Plasmo
using Plots, PlasmoPlots
gr()
using Test

function test_plot_1()
    graph = OptiGraph()

    @optinode(graph,n1)
    @variable(n1, y >= 2)
    @variable(n1,x >= 0)
    @constraint(n1,x + y >= 3)
    @objective(n1, Min, y)

    @optinode(graph,n2)
    @variable(n2, y)
    @variable(n2,x >= 0)
    @constraint(n2,x + y >= 3)
    @objective(n2, Min, y)

    @optinode(graph,n3)
    @variable(n3, y )
    @variable(n3,x >= 0)
    @constraint(n3,x + y >= 3)
    @objective(n3, Min, y)

    @linkconstraint(graph, n1[:x] + n2[:x] + n3[:x] == 3)

    plt_graph1 = layout_plot(graph,node_labels = true,markersize = 60,labelsize = 30,
    linewidth = 4,layout_options = Dict(:tol => 0.01,:iterations => 2));

    plt_matrix1 = matrix_plot(graph,node_labels = true,markersize = 30);

    return true
end

function test_plot_2()
    graph1 = OptiGraph()

    @optinode(graph1,n1)
    @variable(n1, x >= 0)
    @variable(n1, y >= 2)
    @constraint(n1,x + y >= 3)
    @objective(n1, Min, y)

    @optinode(graph1,n2)
    @variable(n2,x >= 0)
    @variable(n2, y >= 2)
    @constraint(n2,x + y >= 3)
    @objective(n2, Min, y)

    @optinode(graph1,n3)
    @variable(n3,x >= 0)
    @variable(n3, y >= 2)
    @constraint(n3,x + y >= 3)
    @objective(n3, Min, y)

    @linkconstraint(graph1, n1[:x] + n2[:x] + n3[:x] == 3)

    graph2 = OptiGraph()

    @optinode(graph2,n4)
    @variable(n4, x >= 0)
    @variable(n4, y >= 2)
    @constraint(n4,x + y >= 5)
    @objective(n4, Min, y)

    @optinode(graph2,n5)
    @variable(n5, x >= 0)
    @variable(n5, y >= 2)
    @constraint(n5,x + y >= 5)
    @objective(n5, Min, y)

    @optinode(graph2,n6)
    @variable(n6, x >= 0)
    @variable(n6, y >= 2 )
    @constraint(n6,x + y >= 5)
    @objective(n6, Min, y)

    @linkconstraint(graph2, n4[:x] + n5[:x] + n6[:x] == 5)

    graph3 = OptiGraph()

    @optinode(graph3,n7)
    @variable(n7, x >= 0)
    @variable(n7, y >= 2)
    @constraint(n7,x + y >= 7)
    @objective(n7, Min, y)

    @optinode(graph3,n8)
    @variable(n8, x >= 0)
    @variable(n8, y >= 2)
    @constraint(n8,x + y >= 7)
    @objective(n8, Min, y)

    @optinode(graph3,n9)
    @variable(n9,x >= 0)
    @variable(n9, y >= 2)
    @constraint(n9,x + y >= 7)
    @objective(n9, Min, y)

    @linkconstraint(graph3, n7[:x] + n8[:x] + n9[:x] == 7)

    graph0 = OptiGraph()
    @optinode(graph0,n0)
    @variable(n0,x)
    @constraint(n0,x >= 0)

    add_subgraph!(graph0,graph1)
    add_subgraph!(graph0,graph2)
    add_subgraph!(graph0,graph3)
    @linkconstraint(graph0,n3[:x] + n5[:x] + n7[:x] == 10)

    @linkconstraint(graph0,n0[:x] + n3[:x] == 3)
    @linkconstraint(graph0,n0[:x] + n5[:x] == 5)
    @linkconstraint(graph0,n0[:x] + n7[:x] == 7)

    for (i,node) in enumerate(all_nodes(graph0))
        node.label = "n$i"
    end

    plt_graph2 = layout_plot(graph0,node_labels = true,markersize = 60,labelsize = 30,linewidth = 4,subgraph_colors = true,
    layout_options = Dict(:tol => 0.001,:C => 2, :K => 4, :iterations => 5));

    plt_matrix2 = matrix_plot(graph0,node_labels = true,subgraph_colors = true,markersize = 16);

    return true
end

function run_tests()
    for name in names(@__MODULE__; all = true)
        if !startswith("$(name)", "test_")
            continue
        end
        @testset "$(name)" begin
            getfield(@__MODULE__, name)()
        end
    end
end

end

TestPlots.run_tests()


text:
using SurrogateModelOptim
using PlotlyJS
using Statistics
using Random

# Optimize the test function Rosenbrock
function rosenbrock_2D(x)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end
search_range=[(-5.0,5.0),(-5.0,5.0)]
fun_max = 90036.0
fun_min = 0.0
noiselevel = (fun_max-fun_min)*0.1 # 10% noise

noisy_rosenbrock_2D = function (x)
    rosenbrock_2D(x)-noiselevel/2+noiselevel*rand()
end

result = smoptimize(noisy_rosenbrock_2D, search_range;
                    options=SurrogateModelOptim.Options(
                    smooth=:single,
                    iterations=20,
                    num_interpolants=10, #Preferably even number of added processes
                    num_start_samples=4,
                    create_final_surrogate=true, #Use the results from last iteration to
                                                 #re-create the surrogate before using it for plotting
                        ))
show(result)

function plot_fun_2D(fun,sr,title)
    N = 51
    x = range(sr[1][1], stop = sr[1][2], length = N)
    y = range(sr[2][1], stop = sr[2][2], length = N)
    grid(x,y) = [x,y]
    z = @. fun(grid(x,y'))

    trace = surface(x=x,y=y,z=z, colorscale="Viridis")
    layout = Layout(title=title)
    p = plot(trace,layout)
end

display(plot_fun_2D(noisy_rosenbrock_2D,search_range,"Original function with added noise"))
display(plot_fun_2D(rosenbrock_2D,search_range,"Original function without noise"))
display(plot_fun_2D(x->median(result.sm_interpolant(x)),search_range,"Estimated function"))


text:
module IterableTables

using Requires, IteratorInterfaceExtensions, TableTraits, TableTraitsUtils

# include("integrations/datastreams.jl")
# include("integrations/gadfly.jl")
# include("integrations/timeseries.jl")
# include("integrations/juliadb.jl")
include("integrations/generators.jl")
# include("integrations/temporal.jl")

function __init__()
    @require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" if !isdefined(DataFrames, :Tables)
        include("integrations/dataframes.jl")
    end
    @require StatsModels="3eaba693-59b7-5ba5-a881-562e759f1c8d" include("integrations/statsmodels.jl")
end

end # module