text:
abstract type TrainingMethods end

struct Batch <: TrainingMethods end
struct BatchL2 <: TrainingMethods
    lambda::AbstractFloat
end

struct BatchCM <: TrainingMethods end
struct BatchCML2 <: TrainingMethods
    lambda::AbstractFloat
end

function train_core(W::AbstractMatrix{T}, b::AbstractMatrix{T}, X::AbstractMatrix{T}, Y::AbstractMatrix{T},
                    mf::MappingFunctions, tm::Batch) where {T<:BaseType}
    H = mapping(W, b, X, mf)
    return Y*pinv(H);
end

function train_core(W::AbstractMatrix{T}, b::AbstractMatrix{T}, X::AbstractMatrix{T}, Y::AbstractMatrix{T},
                    mf::MappingFunctions, tm::BatchL2) where {T<:BaseType}
    H = mapping(W, b, X, mf)
    L = size(H,1)
    N = size(H,2)
    tH = ctranspose(H)
    if N > L
        return Y*tH*pinv(H*tH + tm.lambda*eye(T,L))
    else
        return Y*(pinv(tH*H + tm.lambda*eye(T,N))*tH)
    end
end

function train_core(W::AbstractMatrix{T}, b::AbstractMatrix{T}, X::AbstractMatrix{T}, Y::AbstractMatrix{T},
                    mf::MappingFunctions, tm::BatchCM) where {T<:BaseType}
    A, B = cormat(W, b, X, Y, mf)
    return B*pinv(A)'
end

function train_core(W::AbstractMatrix{T}, b::AbstractMatrix{T}, X::AbstractMatrix{T}, Y::AbstractMatrix{T},
                    mf::MappingFunctions, tm::BatchCML2) where {T<:BaseType}
    A, B = cormat(W, b, X, Y, mf)
    return B*pinv(A+tm.lambda*eye(T,size(A,1)))'
end

function cormat(W::AbstractMatrix{T}, b::AbstractMatrix{T}, X::AbstractMatrix{T}, Y::AbstractMatrix{T},
                mf::MappingFunctions) where {T<:BaseType}
    P = size(X,2)
    Nh = size(W, 1)
    No = size(Y, 1)
    Hj = zeros(T, Nh, 1)
    A = zeros(T, Nh, Nh)
    B = zeros(T, No, Nh)
    for i=1:P
        for j=1:Nh
            Hj[j] = mapping(W[j,:], b[j], X[:,i], mf);
        end
        tHj = ctranspose(Hj)
        A += Hj * tHj
        B += Y[:,i] * tHj
    end

    return A, B
end


text:
"""
    MeshImportModule

Module for mesh import  functions.
"""
module MeshImportModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ..FENodeSetModule: FENodeSet
import ..FESetModule: setlabel!, nodesperelem
import ..FESetModule: AbstractFESet, FESetL2, FESetL3, FESetP1
import ..FESetModule: FESetQ4, FESetQ8, FESetQ9
import ..FESetModule: FESetT3, FESetT6, FESetT4, FESetT10
import ..FESetModule: FESetH8, FESetH20, FESetH27
import ..MeshModificationModule: renumberconn!
import ..MeshExportModule: H5MESH
import ..MeshExportModule: VTK
import Unicode: uppercase, isdigit
import LinearAlgebra: norm
using DelimitedFiles

"""
!!! note

The arrays are reallocated as the data is read from files. The size of the
arrays is increased by this much. If the number of entities to be read  is
large, the chunk should be probably increased so that only a few reallocations
are needed.
"""
const chunk = 1000

# Fix up an old style floating-point number without the exponent letter.
function _fixupdecimal(s)
    os = ""
    for i = length(s):-1:1
        os = s[i] * os
        if (s[i] == '-') && (i > 1) && (uppercase(s[i-1]) != "E") && isdigit(s[i-1])
            os = "E" * os
        end
        if (s[i] == '+') && (i > 1) && (uppercase(s[i-1]) != "E") && isdigit(s[i-1])
            os = "E" * os
        end
    end
    return os
end

"""
    import_NASTRAN(filename)

Import tetrahedral (4- and 10-node) NASTRAN mesh (.nas file).

Limitations:
1. only the GRID and CTETRA  sections are read.
2. Only 4-node and 10-node tetrahedra  are handled.
3. The file should be free-form (data separated by commas).
Some fixed-format files can also be processed (large-field, but not small-field).

# Output
Data dictionary, with keys
- "`fens`" (finite element nodes),
- "`fesets`" (array of finite element sets).
"""
function import_NASTRAN(filename; allocationchunk=chunk, expectfixedformat = false)
    lines = readlines(filename)

    nnode = 0
    node = zeros(allocationchunk, 4)
    maxnodel = 10
    nelem = 0
    elem = zeros(FInt, allocationchunk, maxnodel + 3)

    current_line = 1
    while true
        if current_line >= length(lines)
            break
        end
        temp  = lines[current_line]
        current_line = current_line + 1
        if (length(temp) >= 4) && (uppercase(temp[1:4]) == "GRID")
            fixedformat = expectfixedformat
            largefield = false
            if (length(temp) >= 5) && (uppercase(temp[1:5]) == "GRID*")
                largefield = true; fixedformat = true
            end
            @assert (!fixedformat) || (fixedformat && largefield) "Can handle either free format or large-field fixed format"
            nnode = nnode + 1
            if size(node, 1) < nnode
                node = vcat(node, zeros(allocationchunk, 4))
            end
            if fixedformat
                # $------1-------2-------3-------4-------5-------6-------7-------8-------9-------0
                # GRID*                  5               0-1.66618812195+1-3.85740337853+0
                # *       1.546691269367+1               0
                node[nnode, 1] = parse(Float64, _fixupdecimal(temp[9:24]))
                node[nnode, 2] = parse(Float64, _fixupdecimal(temp[41:56]))
                node[nnode, 3] = parse(Float64, _fixupdecimal(temp[57:72]))
                temp  = lines[current_line]
                current_line = current_line + 1
                node[nnode, 4] = parse(Float64, _fixupdecimal(temp[9:24]))
            else
                # Template:
                #   GRID,1,,-1.32846E-017,3.25378E-033,0.216954
                A = split(replace(temp, "," => " "))
                for  six = 1:4
                    node[nnode, six] = parse(Float64, A[six+1])
                end
            end
        end # GRID
        if (length(temp) >= 6) && (uppercase(temp[1:6]) == "CTETRA")
            # Template:
            #   CTETRA,1,3,15,14,12,16,8971,4853,8972,4850,8973,4848
            nelem = nelem + 1
            if size(elem, 1) < nelem
                elem = vcat(elem, zeros(FInt, allocationchunk, maxnodel + 3))
            end
            continuation = ""
            fixedformat = (length(temp) == 72) || expectfixedformat # Is this fixed-format record? This is a guess based on the length of the line.
            if fixedformat
                continuation = temp[min(length(temp), 72):length(temp)]
                temp = temp[1:min(length(temp), 72)]
            end
            A = split(replace(temp, "," => " "))
            elem[nelem, 1] = parse(FInt, A[2])
            elem[nelem, 2] = parse(FInt, A[3])
            if length(A) == 7  #  nodes per element  equals  4
                nperel = 4
            else
                nperel = 10
                if length(A) < 13 # the line is continued: read the next line
                    temp  = lines[current_line]
                    current_line = current_line + 1
                    if fixedformat
                        temp = temp[length(continuation)+1:min(length(temp), 72)]
                    end
                    temp = strip(temp)
                    Acont = split(replace(temp, "," => " "))
                    A = vcat(A, Acont)
                end
            end
            elem[nelem, 3] = nperel
            for  six = 1:nperel
                elem[nelem, six+3] = parse(FInt, A[six+3])
            end
        end # CTETRA

    end # while

    # Trim the arrays
    node = node[1:nnode,:]
    elem = elem[1:nelem,:]

    # The nodes need to be in serial order:  if they are not,  the element
    # connectivities  will not point at the right nodes
    @assert  norm(collect(1:nnode)-node[:,1]) == 0 "Nodes are not in serial order"

    # Process output arguments
    # Extract coordinates
    xyz=node[:,2:4]
    # Cleanup element connectivities
    ennod = unique(elem[:,3])
    @assert length(ennod) == 1 "Cannot handle mixture of element types"

    @assert ((ennod[1] == 4) || (ennod[1] == 10)) "Unknown element type"
    conn = elem[:,4:3+convert(FInt, ennod[1])]
    pids = elem[:,2] # property identifier: different identifiers, different materials
    upids = unique(pids)

    # Create output arguments. First the nodes
    fens = FENodeSet(xyz)

    # Now for each distinct property identifier, create one finite element set
    fesets = AbstractFESet[]
    for pid in upids
        ix = [i for i in 1:length(pids) if pids[i] == pid]
        C = conn[ix, :]
        if ennod[1] == 4
            fes = FESetT4(C)
        else
            fes = FESetT10(C)
        end
        push!(fesets, fes)
    end

    output = FDataDict("fens"=>fens, "fesets"=>fesets, "property_ids"=>upids)
    return output
end

mutable struct _AbaqusElementSection
    ElementLine::AbstractString
    nelem::FInt
    elem::Array{FInt,2}
end

mutable struct _AbaqusNSetSection
    NSetLine::AbstractString
    nsetname::AbstractString
    generate::Bool
    nodes::Vector{FInt}
end

"""
    import_ABAQUS(filename)

Import Abaqus mesh (.inp file).

Limitations:
1. Only the `*NODE` and `*ELEMENT`  sections are read
2. Only 4-node and 10-node tetrahedra, 8-node or 20-node  hexahedra, 4-node
   quadrilaterals, 3-node triangles are handled.

# Output
Data dictionary, with keys
- "`fens`" = finite element nodes.
- "`fesets`" = array of finite element sets.
- "`nsets`" = dictionary of "node sets", the keys are the names of the sets.
"""
function import_ABAQUS(filename; allocationchunk=chunk)
    lines = readlines(filename)

    maxelnodes = 20

    warnings = String[]
    nnode = 0
    node = zeros(allocationchunk, 4)
    Reading_nodes = false
    next_line = 1
    while true
        if next_line > length(lines)
            break
        end
        temp = uppercase(strip(lines[next_line]))
        next_line = next_line + 1
        if (length(temp) >= 5) && (temp[1:5] == "*NODE")
            Reading_nodes = true
            nnode = 0
            node = zeros(allocationchunk, 4)
            temp = uppercase(strip(lines[next_line]))
            next_line = next_line + 1
        end
        if Reading_nodes
            if temp[1:1] == "*" # another section started
                Reading_nodes = false
                break
            end
            nnode = nnode + 1
            if size(node, 1) < nnode # if needed, allocate more space
                node = vcat(node, zeros(allocationchunk, 4))
            end
            A = split(replace(temp, "," => " "))
            for  six = 1:length(A)
                node[nnode, six] = parse(Float64, A[six])
            end
        end
    end # while

    nelemset = 0
    elemset = _AbaqusElementSection[]
    Reading_elements = false
    next_line = 1
    while true
        if next_line > length(lines)
            break
        end
        temp = uppercase(strip(lines[next_line]))
        next_line = next_line + 1
        if (length(temp) >= 8) && (temp[1:8] == "*ELEMENT")
            Reading_elements = true
            nelemset = nelemset + 1
            nelem = 0
            a = _AbaqusElementSection(temp, nelem, zeros(FInt, allocationchunk, maxelnodes+1))
            push!(elemset, a)
            temp = uppercase(strip(lines[next_line]))
            next_line = next_line + 1
        end
        if Reading_elements
            if temp[1:1] == "*" # another section started
                Reading_elements = false
                break
            end
            elemset[nelemset].nelem = elemset[nelemset].nelem + 1
            if size(elemset[nelemset].elem, 1) < elemset[nelemset].nelem
                elemset[nelemset].elem = vcat(elemset[nelemset].elem,
                zeros(FInt, allocationchunk, maxelnodes+1))
            end
            A = split(temp, ",")
            if (A[end] == "") # the present line is continued on the next one
                temp = uppercase(strip(lines[next_line]))
                next_line = next_line + 1
                Acont = split(temp, ",")
                A = vcat(A[1:end-1], Acont)
            end
            for ixxxx = 1:length(A)
                elemset[nelemset].elem[elemset[nelemset].nelem, ixxxx] = parse(FInt, A[ixxxx])
            end
        end
    end # while

    node = node[1:nnode, :] # truncate the array to just the lines read
    # The nodes need to be in serial order:  if they are not,  the element
    # connectivities  will not point at the right nodes. So,  if that's the case we
    # will figure out the sequential numbering of the nodes  and then we will
    # renumbered the connectivvities of the elements.
    newnumbering = collect(1:nnode)
    if norm(collect(1:nnode)-node[:,1]) != 0
        newnumbering = zeros(FInt, convert(FInt, maximum(node[:,1])))
        jn = 1
        for ixxxx = 1:size(node, 1)
            if node[ixxxx,1] != 0
                on = convert(FInt, node[ixxxx,1])
                newnumbering[on] = jn
                jn = jn + 1
            end
        end
    end

    nsetsarr = _AbaqusNSetSection[]
    Reading_nsetsarr = false
    next_line = 1
    while true
        if next_line > length(lines)
            break
        end
        temp = uppercase(strip(lines[next_line]))
        next_line = next_line + 1
        if (length(temp) >= 5) && (temp[1:5] == "*NSET")
            Reading_nsetsarr = true
            # *NSET, NSET=THE-REGION, GENERATE
            l = split(temp, ",")
            sn = ""
            generate = false
            for t in l
                t = rstrip(lstrip(t))
                ts = split(t, "=")
                if ts[1] == "NSET"
                    sn = ts[2]
                end
                if ts[1] == "GENERATE"
                    generate = true
                end
            end
            a = _AbaqusNSetSection(temp, sn, generate, FInt[])
            push!(nsetsarr, a)
            temp = uppercase(strip(lines[next_line]))
            next_line = next_line + 1
        end
        if Reading_nsetsarr
            if temp[1:1] == "*" # another section started
                Reading_nsetsarr = false
            else
                A = split(temp, ",")
                if nsetsarr[end].generate
                    fn = parse(FInt, A[1])
                    ln = parse(FInt, A[2])
                    st = parse(FInt, A[3])
                    nsetsarr[end].nodes = collect(fn:st:ln)
                else
                    for ixxxx = 1:length(A)
                        if !isempty(A[ixxxx])
                            nn = parse(FInt, A[ixxxx])
                            push!(nsetsarr[end].nodes,  nn)
                        end
                    end
                end
            end
        end
    end # while

    # Process output arguments
    # Nodes
    xyz = node[:,2:4]
    fens = FENodeSet(xyz)

    # Element sets
    fesets = AbstractFESet[]

    function feset_construct(elemset1)
        temp = uppercase(strip(elemset1.ElementLine))
        b  = split(temp, ",")
        for ixxx = 1:length(b)
            c = split(b[ixxx], "=")
            if (uppercase(strip(c[1])) == "TYPE") && (length(c) > 1)
                TYPE = uppercase(strip(c[2]))
                if (length(TYPE) >= 4) && (TYPE[1:4] == "C3D8")
                    return FESetH8(elemset1.elem[:, 2:9])
                elseif (length(TYPE) >= 5) && (TYPE[1:5] == "C3D20")
                    return FESetH20(elemset1.elem[:, 2:21])
                elseif ((length(TYPE) >= 4) && (TYPE[1:4] == "C3D4")) || ((length(TYPE) >= 5) && (TYPE[1:5] == "DC3D4"))
                    return FESetT4(elemset1.elem[:, 2:5])
                elseif (length(TYPE) >= 5) && (TYPE[1:5] == "C3D10")
                    return FESetT10(elemset1.elem[:, 2:11])
                elseif (length(TYPE) >= 5) && (TYPE[1:5] == "DC2D3")
                    return FESetT3(elemset1.elem[:, 2:4])
                elseif (length(TYPE) >= 2) && (TYPE[1:2] == "S3")
                    return FESetT3(elemset1.elem[:, 2:4])
                elseif (length(TYPE) >= 5) && ((TYPE[1:5] == "CPS4R") || (TYPE[1:4] == "CPS4"))
                    return FESetQ4(elemset1.elem[:, 2:5])
                elseif (length(TYPE) >= 4) && (TYPE[1:4] == "CPS3")
                    return FESetT3(elemset1.elem[:, 2:4])
                elseif (length(TYPE) >= 4) && (TYPE[1:4] == "CPS6")
                    return FESetT6(elemset1.elem[:, 2:7])
                else
                    return nothing
                end
            end
        end
    end

    for ixxxx = 1:length(elemset)
        elemset[ixxxx].elem = elemset[ixxxx].elem[1:elemset[ixxxx].nelem, :]
        fes = feset_construct(elemset[ixxxx])
        if (fes == nothing)
            push!(warnings, "Don't know how to handle " * elemset[ixxxx].ElementLine)
        else
            fes = renumberconn!(fes, newnumbering)
            push!(fesets, fes)
        end
    end

    nsets = Dict()
    for ixxxx = 1:length(nsetsarr)
        nsets[nsetsarr[ixxxx].nsetname] = nsetsarr[ixxxx].nodes
    end

    output = FDataDict("fens"=>fens, "fesets"=>fesets, "nsets"=>nsets, "warnings"=>warnings)
    return output
end

"""
    import_MESH(filename)

Import mesh in the MESH format (.mesh, .xyz, .conn triplet of files).

# Output
Data dictionary, with keys
- "`fens`" = finite element nodes.
- "`fesets`" = array of finite element sets.
"""
function import_MESH(filename)
    warnings = String[]
    meshfilebase, ext = splitext(filename)
    if ext == ""
        ext = ".mesh"
    end
    meshfile = meshfilebase * ext
    meshfiledir = dirname(meshfile)
    # Mesh file in the format of the FinEtools .mesh file
    datinfo = open(meshfile, "r") do file
        readdlm(file)
    end
    # The first line is the name of the file with the coordinates
    Xfile = isfile(datinfo[1]) ? datinfo[1] : joinpath(meshfiledir, datinfo[1])
    X = open(Xfile, "r") do file
        readdlm(file, ' ', Float64)
    end
    # The second line is the name of the element name tag.
    etype = datinfo[2]
    Cfile = isfile(datinfo[3]) ? datinfo[3] : joinpath(meshfiledir, datinfo[3])
    C = open(Cfile, "r") do file
        readdlm(file, ' ', Int64)
    end

    fens = FENodeSet(X)

    if     etype == "H8"
        fes = FESetH8(C)
    elseif etype == "H20"
        fes = FESetH20(C)
    elseif etype == "T4"
        fes = FESetT4(C)
    elseif etype == "T10"
        fes = FESetT10(C)
    elseif etype == "Q4"
        fes = FESetQ4(C)
    elseif etype == "Q8"
        fes = FESetQ8(C)
    elseif etype == "T3"
        fes = FESetT3(C)
    elseif etype == "T6"
        fes = FESetT6(C)
    else
        push!(warnings, "Don't know how to handle " * etype)
        fes = nothing
    end

    # Is an optional label file present?
    if length(datinfo) >= 4
        labelfile = isfile(datinfo[4]) ? datinfo[4] : joinpath(meshfiledir, datinfo[4])
        lab = open(labelfile, "r") do file
            readdlm(file, ' ', Int64)
        end
        setlabel!(fes, vec(lab))
    end

    output = FDataDict("fens"=>fens, "fesets"=>[fes], "warnings"=>warnings)
    return output
end

"""
    import_H5MESH(filename)

Import mesh in the H5MESH format (.h5mesh file).

# Output
Data dictionary, with keys
- "`fens`" = finite element nodes.
- "`fesets`" = array of finite element sets.
"""
function import_H5MESH(filename)
    warnings = String[]
    meshfilebase, ext = splitext(filename)
    if ext == ""
        ext = ".h5mesh"
    end
    fname = meshfilebase

    X = H5MESH.retrieve_matrix(fname, ext, "xyz")
    fens = FENodeSet(X)

    etype = H5MESH.retrieve_number(fname, ext, "etype")

    C = H5MESH.retrieve_matrix(fname, ext, "conn")
    if eltype(C) != FInt
        C = FInt.(C)
    end
    if     etype == VTK.H8
        fes = FESetH8(C)
    elseif etype == VTK.H20
        fes = FESetH20(C)
    elseif etype == VTK.T4
        fes = FESetT4(C)
    elseif etype == VTK.T10
        fes = FESetT10(C)
    elseif etype == VTK.Q4
        fes = FESetQ4(C)
    elseif etype == VTK.Q8
        fes = FESetQ8(C)
    elseif etype == VTK.T3
        fes = FESetT3(C)
    elseif etype == VTK.T6
        fes = FESetT6(C)
    else
        push!(warnings, "Don't know how to handle " * etype)
        fes = nothing
    end

    lab = H5MESH.retrieve_matrix(fname, ext, "label")
    setlabel!(fes, vec(lab))

    output = FDataDict("fens"=>fens, "fesets"=>[fes], "warnings"=>warnings)
    return output
end

end #Module

ext: jl


text:
"""
$(SIGNATURES)

Return whether `vertices` are ordered according to vertex order type `O` (a subtype of [`VertexOrder`](@ref)),
and as a result *strongly* convex (see e.g. [CGAL documentation](https://doc.cgal.org/latest/Convex_hull_2/index.html)
for a definition of strong convexity).
"""
function is_ordered_and_strongly_convex(vertices::AbstractVector{<:PointLike}, order::Type{O}) where {O<:VertexOrder}
    op = orientation_comparator(O)
    n = length(vertices)
    n <= 2 && return true
    @inbounds begin
        δprev = vertices[n] - vertices[n - 1]
        δnext = vertices[1] - vertices[n]
        for i in Base.OneTo(n - 1)
            op(cross2(δprev, δnext), 0) || return false
            δprev = δnext
            δnext = vertices[i + 1] - vertices[i]
        end
        return op(cross2(δprev, δnext), 0)
    end
end


text:
@testset "config_namespace" begin
    @testset "user namespace" begin
        # Validate that we can process whatever the current system's namespace is
        result = @test K8sClusterManagers.config_namespace() isa Union{String,Nothing}

        if !(result isa Test.Pass)
            @info "kubectl config view:\n" * read(`$(kubectl()) config view`, String)
        end
    end

    @testset "multiple contexts" begin
        config_path = tempname()
        write(config_path,
              """
              apiVersion: v1
              kind: Config
              contexts:
              - context:
                  namespace: extra-namespace
                name: extra
              - context:
                  namespace: test-namespace
                name: test
              current-context: test
              """)

        withenv("KUBECONFIG" => config_path) do
            @test K8sClusterManagers.config_namespace() == "test-namespace"
        end
    end

    # Note: Mirrors the config returned when running in the CI environment
    @testset "no current context" begin
        config_path = tempname()
        write(config_path,
              """
              apiVersion: v1
              kind: Config
              current-context: ""
              """)

        withenv("KUBECONFIG" => config_path) do
            @test K8sClusterManagers.config_namespace() === nothing
        end
    end

    @testset "no namespace for context" begin
        config_path = tempname()
        write(config_path,
              """
              apiVersion: v1
              kind: Config
              contexts:
              - context:
                  cluster: ""
                  user: ""
                name: test
              current-context: test
              """)

        withenv("KUBECONFIG" => config_path) do
            @test K8sClusterManagers.config_namespace() == nothing
        end
    end

    # Note: Can occur when the current context is deleted
    @testset "bad active context" begin
        config_path = tempname()
        write(config_path,
              """
              apiVersion: v1
              kind: Config
              contexts:
              - context:
                  namespace: test-namespace
                name: test
              current-context: foo
              """)

        withenv("KUBECONFIG" => config_path) do
            @test K8sClusterManagers.config_namespace() == nothing
        end
    end
end

@testset "pod_namespace" begin
    @testset "inside of pod" begin
        patches = [@patch isfile(f) = f == NAMESPACE_FILE
                   @patch read(f, ::Type{String}) = "pod-namespace"]

        apply(patches) do
            @test K8sClusterManagers.pod_namespace() == "pod-namespace"
        end
    end

    @testset "outside of pod" begin
        patches = [@patch isfile(f) = false]

        apply(patches) do
            @test K8sClusterManagers.pod_namespace() === nothing
        end
    end
end

@testset "current_namespace" begin
    @testset "fallback namespace" begin
        config_path = touch(tempname())
        patches = [@patch isfile(f) = false]

        withenv("KUBECONFIG" => config_path) do
            apply(patches) do
                @test K8sClusterManagers.current_namespace() === DEFAULT_NAMESPACE
            end
        end
    end
end


text:
"""
    ProductArray(f, a1, a2)

Takes the Cartesian product of two arrays, `a1` and `a2`, by applying `f` to each pairwise
combination of elements of `a1` and `a2`. The dimensionality of the output is the sum of
the dimensionality of the inputs.

# Examples

julia> ProductArray(tuple, 1:3, 1:4)
3×4 ProductArray{Tuple{Int64,Int64},2,typeof(tuple),UnitRange{Int64},UnitRange{Int64}}:
 (1, 1)  (1, 2)  (1, 3)  (1, 4)
 (2, 1)  (2, 2)  (2, 3)  (2, 4)
 (3, 1)  (3, 2)  (3, 3)  (3, 4)
"""
struct ProductArray{T, N, F, A1, A2} <: AbstractArray{T, N}
    f::F
    a1::A1
    a2::A2
end

function ProductArray(f::F, a1::AbstractArray{<:Any, n}, a2::AbstractArray{<:Any, m}) where {F, n, m}
    T = Core.Compiler.return_type(f, Tuple{eltype(a1), eltype(a2)})
    return ProductArray{T, _add(n, m), F, typeof(a1), typeof(a2)}(f, a1, a2)
end

size(a::ProductArray) = (size(a.a1)..., size(a.a2)...)
axes(a::ProductArray) = (axes(a.a1)..., axes(a.a2)...)

@propagate_inbounds function getindex(a::ProductArray{T, n}, i::Vararg{Int, n}) where {T, n}
    n1 = ndims(a.a1)
    n2 = ndims(a.a2)

    (i1, i2) = split_inds(Val(n1), Val(n2), i)

    x1 = a.a1[i1...]
    x2 = a.a2[i2...]

    return a.f(x1, x2)::T
end

@generated function split_inds(::Val{n1}, ::Val{n2}, i::Tuple{Vararg{Int, n3}}) where {n1, n2, n3}
    @assert n3 == n1 + n2

    exprs1 = Any[ :(i[$j]) for j = 1:n1 ]
    exprs2 = Any[ :(i[$j]) for j = (1+n1):n3 ]

    return quote
        @_inline_meta
        return ( tuple($(exprs1...)), tuple($(exprs2...)) )
    end
end


text:
"""
onesite infinite dmrg
"""
@with_kw struct Idmrg1{} <: Algorithm
    tol_galerkin::Float64 = Defaults.tol
    tol_gauge::Float64 = Defaults.tolgauge
    maxiter::Int = Defaults.maxiter
    verbose::Bool = Defaults.verbose
end

function find_groundstate(st::InfiniteMPS, ham::Hamiltonian,alg::Idmrg1,opars=params(st,ham))
    pars = SimpleEnv(st,opars);

    curu = [st.AR[i] for i in 1:length(st)];
    prevc = st.CR[0];

    err = 0.0;

    for topit in 1:alg.maxiter
        curc = copy(prevc);

        for i in 1:length(st)
            @tensor curu[i][-1 -2;-3] := curc[-1,1]*curu[i][1,-2,-3]

            (eigvals,vecs) = let st=st,pars=pars
                eigsolve(curu[i],1,:SR,Lanczos()) do x
                    ac_prime(x,i,st,pars)
                end
            end

            (curu[i],curc)=leftorth!(vecs[1])

            #partially update pars
            setleftenv!(pars,i+1,st,transfer_left(leftenv(pars,i,st),ham,i,curu[i]))
        end

        for i in length(st):-1:1

            @tensor curu[i][-1 -2;-3] := curu[i][-1,-2,1]*curc[1,-3]

            (eigvals,vecs) = let st=st,pars=pars
                eigsolve(curu[i],1,:SR,Lanczos()) do x
                    ac_prime(x,i,st,pars)
                end
            end

            (curc,temp)=rightorth(vecs[1],(1,),(2,3,))
            curu[i] = permute(temp,(1,2),(3,))

            #partially update pars
            setrightenv!(pars,i-1,st,transfer_right(rightenv(pars,i,st),ham,i,curu[i]))
        end

        err = norm(curc-prevc)
        prevc = curc;
        err<alg.tol_galerkin && break;

        alg.verbose && println("idmrg iter $(topit) err $(err)")
    end

    nst = InfiniteMPS(curu,tol=alg.tol_gauge);
    return nst,params(nst,ham,tol=opars.tol,maxiter=opars.maxiter),err;
end


text:
abstract type Boundary end
abstract type BoundaryMedium end

abstract type BoundaryFlatSubtypes end

struct PeriodicBoundaryCondition<:BoundaryMedium end

defaultFlatBoundary(x,t) = 0
defaultFlatBoundary(x,y,t) = 0
defaultFlatBoundary(x,y,z,t) = 0

struct DirichletBoundaryCondition<:BoundaryMedium
    fMin::Function
    fMax::Function
end
function DirichletBoundaryCondition(f::Function)
    return DirichletBoundaryCondition(f,f)
end
function DirichletBoundaryCondition()
    return DirichletBoundaryCondition(defaultFlatBoundary,defaultFlatBoundary)
end
struct DirichletBoundaryCondition_NewmannBoundaryCondition<:BoundaryMedium
    fMin::Function
    fMax::Function
end
function DirichletBoundaryCondition_NewmannBoundaryCondition()
    return DirichletBoundaryCondition_NewmannBoundaryCondition(defaultFlatBoundary,defaultFlatBoundary)
end

struct NewmannBoundaryCondition<:BoundaryMedium
    fMin::Function
    fMax::Function
end
function NewmannBoundaryCondition(f::Function)
    return NewmannBoundaryCondition(f,f)
end
function NewmannBoundaryCondition()
    return NewmannBoundaryCondition(defaultFlatBoundary,defaultFlatBoundary)
end
struct NewmannBoundaryCondition_DirichletBoundaryCondition<:BoundaryMedium
    fMin::Function
    fMax::Function
end
function NewmannBoundaryCondition_DirichletBoundaryCondition()
    return NewmannBoundaryCondition_DirichletBoundaryCondition(defaultFlatBoundary,defaultFlatBoundary)
end


text:
@doc raw"""
    ExponentiatedKernel([α=1])

The exponentiated kernel is a Mercer kernel given by:

```
    κ(x,y) = exp(α⋅xᵀy)   α > 0
```
where `α` is a positive scaling parameter.

# Examples

```jldoctest; setup = :(using MLKernels)
julia> ExponentiatedKernel()
ExponentiatedKernel{Float64}(1.0)

julia> ExponentiatedKernel(2)
ExponentiatedKernel{Float64}(2.0)

julia> ExponentiatedKernel(2.0f0)
ExponentiatedKernel{Float32}(2.0)
```
"""
struct ExponentiatedKernel{T<:AbstractFloat} <: MercerKernel{T}
    α::T
    function ExponentiatedKernel{T}(α::Real=T(1)) where {T<:AbstractFloat}
        @check_args(ExponentiatedKernel, α, α > zero(T), "α > 0")
        return new{T}(α)
    end
end
ExponentiatedKernel(α::T=1.0) where {T<:Real} = ExponentiatedKernel{promote_float(T)}(α)

@inline basefunction(::ExponentiatedKernel) = ScalarProduct()

@inline kappa(κ::ExponentiatedKernel{T}, xᵀy::T) where {T} = exp(κ.α*xᵀy)

function convert(
        ::Type{K},
        κ::ExponentiatedKernel
    ) where {K>:ExponentiatedKernel{T}} where T
    return ExponentiatedKernel{T}(κ.α)
end


text:
# Autogenerated wrapper script for MPIR_jll for x86_64-apple-darwin
export libmpir

JLLWrappers.@generate_wrapper_header("MPIR")
JLLWrappers.@declare_library_product(libmpir, "@rpath/libmpir.23.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libmpir,
        "lib/libmpir.23.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MagicDicts
using Test

@testset "MagicDicts.jl" begin

    let
        D = MagicDict()
        D['A', 'A'] = rand()
        D['A', 'B'] = rand()
        D['B', 'A'] = rand()
        @test all(D['A', 'A':'B'] .== [D['A', 'A'],D['A', 'B']])
        @test all(D['A', ['A','B']] .== [D['A', 'A'],D['A', 'B']])
        @test all(D['A':'B', 'A'] .== [D['A', 'A'],D['B', 'A']])
        @test all(D[['A','B'], 'A'] .== [D['A', 'A'],D['B', 'A']])
        D['C', 'A':'B'] = rand()
        @test all(D['C', 'A':'B'] .== D['C', 'A'])
        D[collect(1:10), 'A':'B'] = rand()
        @test all(D[collect(1:10), 'A':'B'] .== D[1, 'A'])
        D = MagicDict()
        D[1,2,3,4] = rand()
        @test haskey(D, 1)
        @test haskey(D, 1, 2)
        @test haskey(D, 1, 2, 3)
        @test haskey(D, 1, 2, 3, 4)
        @test !haskey(D, 1, 3, 4)
    end
end


text:
# This file was generated, do not modify it. # hide
train, test = partition(eachindex(y), 0.70, shuffle=true, rng=52)
@load XGBoostClassifier
xgb_model = XGBoostClassifier()


text:
s = """
import Base: show
import Base.Math
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
importall Lint
"""
msgs = lintstr(s)
@test isempty(msgs)

# Avoid warning users about dynamic includes.
s = """
script = \"test.jl\"; include(script)
"""
msgs = lintstr(s)
@test msgs[1].code == :I372
@test isempty(filter(x -> !isinfo(x), msgs))


text:
using Test
using StaticArrays

model = Dynamics.DubinsCar()
x,u = rand(model)
n,m = size(model)
N = 101
dt = 0.1
tf = 10.0

X = [@SVector rand(n) for k = 1:N]
U = [@SVector rand(m) for k = 1:N-1]
Z = Traj(x,u,dt,N)
@test Z[1].dt == dt
@test Z[1].t == 0
@test Z[end].dt == 0
@test Z[end].t == tf

Z = Traj(X,U,fill(dt,N))
@test Z[1].dt == dt
@test Z[1].t == 0
@test Z[end].dt == 0
@test Z[end].t ≈ tf


text:
#############################################################
#
# Tweets Functions
#
#############################################################

function get_retweets_id(id::AbstractString; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/statuses/retweets/$(id).json", options)

    #Return array of type TWEETS
    #return to_TWEETS(r)
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function get_single_tweet_id(id::AbstractString; options=Dict{AbstractString, AbstractString}())

    #Add required parameter(s) to options dict
    options["id"] = id

    r = get_oauth("https://api.twitter.com/1.1/statuses/show.json", options)

    #Return array of type TWEETS
    #return to_TWEETS(r)
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function post_destroy_single_tweet_id(id::AbstractString; options=Dict{AbstractString, AbstractString}())

    r = post_oauth("https://api.twitter.com/1.1/statuses/destroy/$(id).json", options)

    #Return array of type TWEETS
    #return to_TWEETS(r)
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function post_status_update(status::AbstractString; options=Dict{AbstractString, AbstractString}())

    #Add status into options Dict
    options["status"] = status

    r = post_oauth("https://api.twitter.com/1.1/statuses/update.json", options)

    #Return array of type TWEETS
    #return to_TWEETS(r)
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function post_status_retweet_id(id::AbstractString; options=Dict{AbstractString, AbstractString}())

    r = post_oauth("https://api.twitter.com/1.1/statuses/retweet/$(id).json", options)

    #Return array of type TWEETS
    #return to_TWEETS(r)
    return r.status == 200 ? to_TWEETS(Requests.json(r)) : r

end

function post_status_update_media()
    error("Twitter API not fully implemented")
end

function get_retweeters_id(id::AbstractString; options=Dict{AbstractString, AbstractString}())

    #Add required parameter(s) to options dict
    options["id"] = id

    #Return a Dict for now, only useful data value is array of ids
    r = get_oauth("https://api.twitter.com/1.1/statuses/retweeters/ids.json", options)

    return r.status == 200 ? Requests.json(r) : r

end


text:
s = string(
"73167176531330624919225119674426574742355349194934",
"96983520312774506326239578318016984801869478851843",
"85861560789112949495459501737958331952853208805511",
"12540698747158523863050715693290963295227443043557",
"66896648950445244523161731856403098711121722383113",
"62229893423380308135336276614282806444486645238749",
"30358907296290491560440772390713810515859307960866",
"70172427121883998797908792274921901699720888093776",
"65727333001053367881220235421809751254540594752243",
"52584907711670556013604839586446706324415722155397",
"53697817977846174064955149290862569321978468622482",
"83972241375657056057490261407972968652414535100474",
"82166370484403199890008895243450658541227588666881",
"16427171479924442928230863465674813919123162824586",
"17866458359124566529476545682848912883142607690042",
"24219022671055626321111109370544217506941658960408",
"07198403850962455444362981230987879927244284909188",
"84580156166097919133875499200524063689912560717606",
"05886116467109405077541002256983155200055935729725",
"71636269561882670428252483600823257530420752963450")

a = [Int32(s[i]-48) for i in 1:length(s)]
println(maximum([reduce(*, a[i:i+12]) for i in 1:length(a)-13]))

max_stars_repo_path: 008 - Largest Product in a Series/008 - Largest_Product_in_a_Series.jl

max_stars_repo_name: LanguageAdventures/Project-Euler-Julia

max_issues_repo_name: LanguageAdventures/Project-Euler-Julia

max_forks_repo_name: LanguageAdventures/Project-Euler-Julia


text:
function clean(phone_number)
    cleaned = replace(phone_number, r"[^0-9]" => "")
    try
        match(r"^1?([2-9][0-9]{2}[2-9][0-9]{6})$", cleaned).captures[1]
    catch
        nothing
    end
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
Random.seed!(1)
const d = 2
const Point = SArray{Tuple{d},Float32,1,d}
using Makie
T = (0:0.001f0:2)[2:end]
n = 200
σ = 0.01f0
μ = 0.2f0
a = 4*0.8f0
b = 4*0.15f0
B = sparse((-μ - a + b)*I, n, n)
B1 = sparse(-1f0*I, n, n)
B2 = sparse(1f0*I, n, n)

for i in 1:n
    k1 = rand(1:n-1)
    k1 += (i <= k1)
    k2 = rand(1:n-2)
    k2 += (min(i, k1) <= k2)
    k2 += (max(i, k1) <= k2)
    @assert i ≠ k1
    @assert i ≠ k2
    @assert k1 ≠ k2

    B[i, k1] = a
    B[i, k2] = -b

    B1[i,k1] = 1
    B2[i, k2] = -1

end
#x0 = 0.1*randn(Point, n)
x = copy(x0)
xs = [x]
ts = [0.0]
t = 0.0
xp = Node(x)
r = .5
p1 = scatter(xp, markersize=3, color=1:n)
xlims!( -r, r)
ylims!( -r, r)

display(p1)
nlz(x) = x/norm(x)
for i in eachindex(T)
    global t, x, xs
    local dt = T[i] - t
    @assert dt > 0
    x = x + (B*x)*dt + σ*sqrt(dt)*randn(Point, n)
     xp[] = x
    t = T[i]
    yield()
    sleep(0.0001)
    push!(ts, t)
    push!(xs, copy(x))
end
error("stop")
p1 = scatter(xp, markersize=3, color=1:n)
xlims!( -r, r)
ylims!( -r, r)
using Base.Iterators: take
framerate=30
record(p1, "boids_animation.mp4", 1:200; framerate = framerate) do i
    xp[] = xs[i*10]
end


text:
@testset "preview" begin
    root = "assets"
    abs_root = abspath(pwd(), root)

    @testset "page structure" begin
        page_dir = @suppress_err preview_demos(joinpath(abs_root, "page", "template"), theme="grid", require_html=false)
        md_index = joinpath(page_dir, "index.md")
        @test_reference joinpath("references", "preview", "index_page.md") read(md_index, String) by=ignore_CR

        @test readdir(page_dir) == ["covers", "index.md", "subsection_1", "subsection_2"]
        @test readdir(joinpath(page_dir, "subsection_2")) == ["card_3.md", "card_4.md"]

        for demo_file in filter(DemoCards.is_democard, readdir(joinpath(page_dir, "subsection_2")))
            demo_file = joinpath(joinpath(page_dir, "subsection_2"), demo_file)
            @test_reference joinpath("references", "preview", basename(demo_file)) read(demo_file, String) by=ignore_CR
        end
        @test "gridtheme.css" in readdir(joinpath(dirname(page_dir), "democards"))
    end

    @testset "section with non page structure" begin
        page_dir = @suppress_err preview_demos(joinpath(abs_root, "preview", "scripts"), theme="grid", require_html=false)
        md_index = joinpath(page_dir, "index.md")
        @test_reference joinpath("references", "preview", "index_sec_nopage.md") read(md_index, String) by=ignore_CR

        @test readdir(page_dir) == ["assets", "covers", "index.md", "scripts"]
        @test readdir(joinpath(page_dir, "covers")) == ["demo2.svg", "democards_logo.svg"]
        @test readdir(joinpath(page_dir, "scripts")) == ["assets", "demo1.ipynb", "demo1.jl", "demo1.md", "demo2.md"]
    end

    @testset "section with page structure" begin
        page_dir = @suppress_err preview_demos(joinpath(abs_root, "page", "template", "subsection_2"), theme="grid", require_html=false)
        md_index = joinpath(page_dir, "index.md")
        @test_reference joinpath("references", "preview", "index_sec_page.md") read(md_index, String) by=ignore_CR

        @test readdir(page_dir) == ["covers", "index.md", "subsection_2"]
        @test readdir(joinpath(page_dir, "subsection_2")) == ["card_3.md", "card_4.md"]
    end

    @testset "section with invalid page structure" begin
        @test_throws ArgumentError @suppress_err preview_demos(joinpath(abs_root, "page"))
    end

    @testset "file with non page structure" begin
        page_dir = @suppress_err preview_demos(joinpath(abs_root, "preview", "scripts", "demo2.md"), theme="grid", require_html=false)
        md_index = joinpath(page_dir, "index.md")
        @test_reference joinpath("references", "preview", "index_file.md") read(md_index, String) by=ignore_CR

        @test readdir(page_dir) == ["assets", "covers", "index.md", "scripts"]
        @test readdir(joinpath(page_dir, "covers")) == ["demo2.svg"]
        @test readdir(joinpath(page_dir, "scripts")) == ["assets", "demo2.md"]
        @test readdir(joinpath(page_dir, "scripts", "assets")) == ["logo.svg"]

        page_dir = @suppress_err preview_demos(joinpath(abs_root, "preview_demo1.jl"), theme="grid", require_html=false)
        md_index = joinpath(page_dir, "index.md")
        @test_reference joinpath("references", "preview", "index_file2.md") read(md_index, String) by=ignore_CR

        @test readdir(page_dir) == ["covers", "index.md", "preview_section"]
        @test readdir(joinpath(page_dir, "covers")) == ["democards_logo.svg"]
        @test readdir(joinpath(page_dir, "preview_section")) == ["preview_demo1.ipynb", "preview_demo1.jl", "preview_demo1.md"]
    end

    @testset "file with page structure" begin
        page_dir = @suppress_err preview_demos(joinpath(abs_root, "page", "template", "subsection_2", "card_3.md"), theme="grid", require_html=false)
        md_index = joinpath(page_dir, "index.md")
        @test_reference joinpath("references", "preview", "index_file_page.md") read(md_index, String) by=ignore_CR

        @test readdir(page_dir) == ["covers", "index.md", "subsection_2"]
        @test readdir(joinpath(page_dir, "subsection_2")) == ["card_3.md"]

        page_dir = @suppress_err preview_demos(joinpath(abs_root, "page", "template", "subsection_2", "card_3.md"), theme=nothing, require_html=false)
        @test !isfile(joinpath(page_dir, "index.md"))
        @test readdir(page_dir) == ["subsection_2"]
        @test readdir(joinpath(page_dir, "subsection_2")) == ["card_3.md"]
    end

    index_page = @suppress_err preview_demos(joinpath(abs_root, "preview"), theme="grid", require_html=true)
    @test isfile(index_page)
    index_page = @suppress_err preview_demos(joinpath(abs_root, "preview"), require_html=true)
    @test isfile(index_page)

    @testset "don't copy unrelated assets" begin
        src_page_dir = joinpath(abs_root, "preview")
        @test readdir(src_page_dir) == ["assets", "other_scripts", "scripts"]
        page_dir = @suppress_err preview_demos(joinpath(src_page_dir, "other_scripts", "nested", "demo3.jl"), theme="grid", require_html=false)

        # `assets` is copied, but `scripts` is not
        @test readdir(page_dir) == ["assets", "covers", "index.md", "other_scripts"]
        @test readdir(joinpath(page_dir, "other_scripts")) == ["assets", "nested"]
        @test readdir(joinpath(page_dir, "other_scripts", "nested")) == ["assets", "demo3.ipynb", "demo3.jl", "demo3.md"]
    end
end


text:
extend(a::Array, ::Type{Val{:count}}) = length(a)
extend(a::Array, ::Type{Val{:append}}) = x->push!(a, x)

function Base.getproperty(a::Array, name::Symbol)
    extend(a, Val{name})
end


text:
# julia -i -q -L replcxx.jl

const __initcxx = """
#include <iostream>
#include <eigen3/Eigen/Eigen>
#include <eigen3/Eigen/Geometry>

using namespace std;
using namespace Eigen;
"""

__replcxx_init() = begin
    while !isdefined(Base,:active_repl)
        sleep(0.1)
    end

    @eval using Cxx
    @eval @cxx_str($__initcxx)
    write(STDIN.buffer,'<');
end

println("/// Cxx.jl REPL ///", "\n\n", __initcxx)
@async __replcxx_init()


text:
#### Chart option functions

#1: colorscheme!
names = ["sales", "administration", "information technology", "customer support", "development",
"marketing"]
max = [6500, 16000, 30000, 38000, 52000, 25000]
values = [4300, 10000, 28000, 35000, 50000, 19000]
values2 = 0.7 .* values
r = radar(names, hcat(values, values2), max, fill = [false, true])
colorscheme!(r, ("acw", "JapaneseGarden"))
@test typeof(r) == EChart

v = [100,60,80,70,50]
n = ["A", "B", "C", "D", "E"]
fn = funnel(n, v)
colorscheme!(fn,  ("acw", "VitaminC"), reversePalette = true)
colorscheme!(fn, "purple")
colorscheme!(fn, ["purple", "yellow"])
@test typeof(fn) == EChart

#2: flip!
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
bm = bar(x, hcat(0.95 .* y, 1.25 .* y, y), color = ["red", "gray", "blue"], stack = true)
flip!(bm)
@test typeof(bm) == EChart
flip!(bm, rotatedims = true)
@test typeof(bm) == EChart

#3: lineargradient
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
ar = area(x, y, color = lineargradient("purple", "cyan"))
smooth!(ar)
xaxis!(ar, name = "Day of Week")
yaxis!(ar, name = "Daily High Temperature °C")
@test typeof(ar) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
ar = area(x, y, color = lineargradient("purple", "cyan", x0 = 0.4, y0 = 0.9, x2 = .5, y2 = 1))
smooth!(ar)
@test typeof(ar) == EChart

#4: radialgradient
pp = pie(["all"], [1], color = radialgradient("red", "purple"))
@test typeof(pp) == EChart

pp2 = pie(["all"], [1], color = radialgradient("red", "purple", x = 0.25, y = 0.3, r = 0.4))
@test typeof(pp2) == EChart

#5: seriesnames!
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
bm = bar(x, hcat(0.95 .* y, 1.25 .* y,y))
seriesnames!(bm, ["Pink Series", "Blue Series", "Tan Series"])
@test typeof(bm) == EChart

#6: slider!
dt = ["2016/08/26", "2016/08/29", "2016/08/30", "2016/08/31", "2016/09/01", "2016/09/02",
"2016/09/06", "2016/09/07", "2016/09/08", "2016/09/09", "2016/09/12", "2016/09/13",
"2016/09/14", "2016/09/15", "2016/09/16", "2016/09/19", "2016/09/20", "2016/09/21",
"2016/09/22", "2016/09/23", "2016/09/26"]

close_ = [12.3800, 12.4700, 12.5500, 12.6000, 12.4400, 12.5000, 12.6700, 12.7000, 12.7300,
12.3800, 12.7000, 12.3800, 12.1400, 12.1100, 12.1100, 12.1100, 12.0000, 12.0900, 12.1800,
12.1700, 12.0100]

open_ = [12.4700, 12.3800, 12.4600, 12.4800, 12.6600, 12.5300, 12.4900, 12.6200, 12.5000,
12.6100, 12.3200, 12.5300, 12.1800, 12.1400, 12.0500, 12.1200, 12.0900, 12.0800, 12.1700,
12.1200, 12.1200]

high_ = [12.5500, 12.5000, 12.5600, 12.6100, 12.7200, 12.5700, 12.6700, 12.7500, 12.7500,
12.6850, 12.7700, 12.6800, 12.3100, 12.1800, 12.1300, 12.3300, 12.1900, 12.1100, 12.2900,
12.2200, 12.1200]

low_ = [12.3400, 12.3800, 12.4300, 12.4800, 12.3500, 12.4600, 12.4300, 12.6200, 12.5000,
12.3800, 12.2800, 12.3300,12.1100, 12.0600, 12.0100, 12.0586, 11.9600, 12.0100, 12.1600,
12.1200, 12.0000]

c = candlestick(dt, open_, close_, low_, high_)
slider!(c)
@test typeof(c) == EChart

#7: smooth!
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
ar = area(x, y, color = lineargradient("purple", "cyan"))
smooth!(ar)
smooth!(ar, 1)
smooth!(ar, [1])
xaxis!(ar, name = "Day of Week")
yaxis!(ar, name = "Daily High Temperature °C")
@test typeof(ar) == EChart

#8: text!
sgrp = scatter(mtcars, :MPG, :HP, :Cyl)
title!(sgrp, text = "Horsepower vs. Miles Per Gallon")
text!(sgrp, subtext = "Source: mtcars dataset from R", left = "0%", bottom = "0%")
@test typeof(sgrp) == EChart

text!(sgrp, 2, sublink = "https://stat.ethz.ch/R-manual/R-devel/library/datasets/html/mtcars.html")
@test typeof(sgrp) == EChart

#9: title!
x = 1:20
y = 3 .+ 5x
a = bar(x, y)
title!(a, text = "Bar Plot Title", subtext = "Secondary Title")
@test typeof(a) == EChart

#10: toolbox!
x = 1:20
y = 3 .+ 5x
a = bar(x, y)
toolbox!(a, chartTypes = ["bar", "line"])
@test typeof(a) == EChart

#11: xarea!/yarea!
sgrp = scatter(mtcars, :MPG, :HP, :Cyl)
xarea!(sgrp, 28, 33)
xarea!(sgrp, 10, 14)
@test typeof(sgrp) == EChart

yarea!(sgrp, 100, 150, series = 2)
@test typeof(sgrp) == EChart

#12: xaxis!/yaxis!
x = 1:10
y = [2^x for x in x]
l = line(x, y)
yaxis!(l, _type = "log", scale = true)
@test typeof(l) == EChart

x = 1:10
y = [2^x for x in x]
l = line(x, y)
yaxis!(l, _type = "log", scale = true, formatter = "{value} kg")
@test typeof(l) == EChart

#13: xgridlines!/ygridlines!
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
l = line(x, y)
ygridlines!(l, show = false)
@test typeof(l) == EChart

#14: xline!/yline!
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
y2 = 3.7 .* y
as = area(x, hcat(y, y2))
yline!(as, "average")
yline!(as, 40.362)
xline!(as, "Tuesday")
@test typeof(as) == EChart

#15: radial!
x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
b = bar(x, y)
radial!(b)
@test typeof(b) == EChart

x = ["Monday","Tuesday","Wednesday","Thursday","Friday","Saturday","Sunday"]
y = [11, 11, 15, 13, 12, 13, 10]
b = bar(x, hcat(0.95 .* y, 1.25 .* y,y), stack = [1,1,2])
radial!(b)
@test typeof(b) == EChart

#16: jitter!
x = rand(1:10, 500)
y = rand(0:10, 500) + x
a = scatter(x, y)
jitter!(a)
@test typeof(a) == EChart

jitter!(a, pctxrange = 0, pctyrange = 0.05)
@test typeof(a) == EChart

jitter!(a, pctxrange = 0.1, pctyrange = 0.05)
@test typeof(a) == EChart

jitter!(a, 1, pctxrange = 0.1, pctyrange = 0.05)
@test typeof(a) == EChart

sm = scatter(1:30, hcat(rand(30), rand(30)))
jitter!(sm)
@test typeof(sm) == EChart

#17: theme!
@test typeof(theme!(sm, ECharts.vintage)) == EChart

@test typeof(Theme(JSON.parsefile(joinpath(dirname(@__FILE__), "..", "docs/js/themes/chalk.json")))) == Theme


text:
"""Convert a `Colorant` to an `(r,g,b,a)` tuple ∈ [0,1]⁴, as accepted by Matplotlib."""
toRGBAtuple(c) = toRGBAtuple(RGBA(c))
toRGBAtuple(c::RGBA) = (c.r, c.g, c.b, c.alpha)

mplcolors = C0, C1, C2, C3, C4, C5, C6, C7, C9, C10 = parse.(RGB,
    ["#1f77b4", "#ff7f0e", "#2ca02c", "#d62728", "#9467bd",
     "#8c564b", "#e377c2", "#7f7f7f", "#bcbd22", "#17becf" ]
)

black = Gray(0)
white = Gray(1)
lightgrey = lightgray = Gray(0.77)

"""
Mix a color with white. `original` specifies how much is left of the original color.
`0`: output is pure white. `1`: output is the original color.
Visually equivalent to setting the color's alpha to `original`, on a white background.
"""
lighten(c::C, original = 0.8) where {C<:Color} = C(mix(RGB(c), RGB(white), 1 - original))
#   White goes last, because colors can't have negative channels (`b - a` in `mix`), even
#   with ColorVectorSpace.

"""
Mix a color with black. `original` specifies how much is left of the original color.
`0`: output is pure black. `1`: output is the original color.
"""
darken(c::C, original = 0.8) where {C<:Color} = C(mix(RGB(black), RGB(c), original))

"""Linearly interpolate ("lerp") between `a` (`t = 0`) and `b` (`t = 1`)."""
mix(a, b, t=0.5) = a + t * (b - a)

"""
De-emphasise part of an Axes by colouring it light grey.
`part` is one of {:xlabel, :ylabel, :xaxis, :yaxis}.
"""
function deemph(part::Symbol, ax; color = lightgrey)
    color = toRGBAtuple(color)
    if part == :xlabel
        ax.xaxis.get_label().set_color(color)
    elseif part == :ylabel
        ax.yaxis.get_label().set_color(color)
        hasproperty(ax, :hylabel) && ax.hylabel.set_color(color)
    elseif part == :xaxis
        foreach(loc -> ax.spines[loc].set_color(color), ["top", "bottom"])
        ax.tick_params(; axis = "x", which = "both", color, labelcolor = color)
    elseif part == :yaxis
        foreach(loc -> ax.spines[loc].set_color(color), ["left", "right"])
        ax.tick_params(; axis = "y", which = "both", color, labelcolor = color)
    end
end


text:
# The license for this code is available at https://github.com/avik-pal/FastStyleTransfer.jl/blob/master/LICENSE.md

#---------------------------Residual Block-----------------------------------

struct ResidualBlock
    conv_layers
    norm_layers
end

ResidualBlock(chs::Int) =
   ResidualBlock((Conv((3,3), chs=>chs, pad = (1,1)), Conv((3,3), chs=>chs, pad = (1,1))), (BatchNorm(chs), BatchNorm(chs)))

function (r::ResidualBlock)(x)
    value = relu.(r.norm_layers[1](r.conv_layers[1](x)))
    r.norm_layers[2](r.conv_layers[2](value)) + x
end

#--------------------------Reflection Pad-------------------------------------

# Paper suggests using Reflection Padding. However normal padding is being used until this layer is implemented

struct ReflectionPad
    dim::Int
end

#----------------------Convolution Block--------------------------------------

ConvBlock(chs::Pair{<:Int,<:Int}, kernel::Tuple{Int,Int}, stride::Tuple{Int,Int} = (1,1), pad::Tuple{Int,Int} = (0,0)) =
    Chain(Conv(kernel, chs, stride = stride, pad = pad), ReflectionPad(kernel[1]÷2))

#-------------------------Upsample--------------------------------------------

Upsample(x) = repeat(x, inner = (2,2,1,1))

#----------------------Upsampling BLock---------------------------------------

# TODO: Use reflection padding instead of zero padding once its implemented

UpsamplingBlock(chs::Pair{<:Int,<:Int}, kernel::Tuple{Int,Int}, stride::Tuple{Int,Int}, upsample::Int, pad::Tuple{Int,Int} = (0,0)) =
    Chain(Conv(kernel, chs, stride = stride, pad = (kernel[1]÷2, kernel[2]÷2)), x -> Upsample(x))


text:
function possemideftri1(opt)
    TOL = 1e-4
    m = JuMP.Model(opt)

    JuMP.@variable(m, x, Int)
    JuMP.@constraint(m, x >= 0)
    JuMP.@variable(m, y >= 0)
    JuMP.@variable(m, Z[1:2, 1:2], Symmetric)
    JuMP.@constraint(m, svec(1.0 * Z) in Hypatia.PosSemidefTriCone{Float64, Float64}(3))
    JuMP.@objective(m, Max, 3x + y - Z[1, 1])
    JuMP.@constraint(m, 3x + 2y <= 10)
    JuMP.@constraint(m, svec([2 x; x 2]) in Hypatia.PosSemidefTriCone{Float64, Float64}(3))
    c1 = JuMP.@constraint(m, Z[1, 2] >= 1)
    c2 = JuMP.@constraint(m, y >= Z[2, 2])
    JuMP.optimize!(m)

    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    @test isapprox(JuMP.objective_value(m), 7.5, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), 7.5, atol = TOL)
    @test isapprox(JuMP.value(x), 2, atol = TOL)
    @test isapprox(JuMP.value(y), 2, atol = TOL)
    @test isapprox(JuMP.value.(vec(Z)), [0.5, 1, 1, 2], atol = TOL)

    JuMP.delete(m, c1)
    JuMP.delete(m, c2)
    JuMP.@constraint(m, x >= 2)
    JuMP.set_lower_bound(Z[1, 2], 2)
    c2 = JuMP.@constraint(m, y >= Z[2, 2] + Z[1, 1])
    JuMP.optimize!(m)
    @test JuMP.termination_status(m) == MOI.INFEASIBLE
    @test JuMP.primal_status(m) == MOI.NO_SOLUTION
    return
end

function possemideftri2(opt)
    TOL = 1e-4
    d = 3
    for is_complex in (false, true)
        R = (is_complex ? ComplexF64 : Float64)
        dim = svec_length(R, d)
        mat = Hermitian(smat(R, Float64.(1:dim)), :U)
        λ₁ = eigmax(mat)
        m = JuMP.Model(opt)

        x = JuMP.@variable(m, [1:dim])
        JuMP.@constraint(m, x in Hypatia.PosSemidefTriCone{Float64, R}(dim))
        JuMP.@objective(m, Max, dot(svec(mat), x))
        x_diag = [x[svec_idx(R, i, i)] for i in 1:d]
        JuMP.@constraint(m, sum(x_diag) == 1)
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.objective_value(m), λ₁, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), λ₁, atol = TOL)
        x_val = JuMP.value.(x)
        x_mat = Hermitian(smat(R, x_val), :U)
        @test isapprox(tr(x_mat), 1, atol = TOL)
        @test eigmin(x_mat) >= -TOL

        JuMP.set_binary.(x_diag)
        JuMP.optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        x_val = JuMP.value.(x)
        x_mat = Hermitian(smat(R, x_val), :U)
        @test isapprox(tr(x_mat), 1, atol = TOL)
        @test eigmin(x_mat) >= -TOL
        @test isapprox(x_mat[d, d], 1, atol = TOL)
    end
    return
end

function possemideftrisparse1(opt)
    # arrow PSD formulation of second-order cone constraint
    TOL = 1e-4
    row_idxs = [1, 2, 3, 4, 2, 3, 4]
    col_idxs = [1, 1, 1, 1, 2, 3, 4]

    for is_complex in (false, true)
        R = (is_complex ? ComplexF64 : Float64)
        KT = Hypatia.PosSemidefTriSparseCone{Hypatia.Cones.PSDSparseDense, Float64, R}
        K = KT(4, row_idxs, col_idxs, false)
        x_dim = vec_length(R, 3)
        x0 = vec_copyto!(zeros(x_dim), -Real.(1:x_dim) .+ 0.5)
        opt_val = norm(Real.(1:x_dim), 2)
        m = JuMP.Model(opt)

        JuMP.@variable(m, y)
        JuMP.@objective(m, Min, y)
        JuMP.@variable(m, x[1:x_dim], Int)
        JuMP.@constraint(m, x .<= x0)
        aff = vcat(y, rt2 * x, y, y, y)
        JuMP.@constraint(m, aff in K)
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(JuMP.value.(x), -Float64.(1:x_dim), atol = TOL)

        JuMP.@constraint(m, y <= opt_val - 0.2)
        JuMP.optimize!(m)
        @test JuMP.termination_status(m) == MOI.INFEASIBLE
        @test JuMP.primal_status(m) == MOI.NO_SOLUTION
    end
    return
end

function possemideftrisparse2(opt)
    TOL = 1e-4
    row_idxs = [1, 2, 3, 2, 3, 3, 4, 4, 5, 5]
    col_idxs = [1, 1, 1, 2, 2, 3, 3, 4, 4, 5]
    KT = Hypatia.PosSemidefTriSparseCone{Hypatia.Cones.PSDSparseDense, Float64, Float64}
    K = KT(5, row_idxs, col_idxs, true)
    m = JuMP.Model(opt)

    JuMP.@variable(m, y[1:3], Int)
    JuMP.@constraint(m, y .<= 20)
    JuMP.@objective(m, Min, sum(y))
    aff = [8, -2, -4, y[1], 7, y[2], 5, y[3], -6, 5]
    JuMP.@constraint(m, y .>= 3)
    JuMP.@constraint(m, aff in K)

    JuMP.optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    @test isapprox(JuMP.objective_value(m), 15, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), 15, atol = TOL)

    JuMP.@constraint(m, y[3] == 3)
    JuMP.optimize!(m)
    @test JuMP.termination_status(m) == MOI.INFEASIBLE
    @test JuMP.primal_status(m) == MOI.NO_SOLUTION
    return
end

function epinormeucl1(opt)
    TOL = 1e-4
    m = JuMP.Model(opt)

    JuMP.@variable(m, x)
    JuMP.@variable(m, y)
    JuMP.@variable(m, z, Int)
    JuMP.@constraint(m, z <= 2.5)
    JuMP.@objective(m, Min, x + 2y)
    JuMP.@constraint(m, [z, x, y] in Hypatia.EpiNormEuclCone{Float64}(3))
    JuMP.set_integer(x)
    JuMP.optimize!(m)

    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    opt_obj = -1 - 2 * sqrt(3)
    @test isapprox(JuMP.objective_value(m), opt_obj, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), opt_obj, atol = TOL)
    @test isapprox(JuMP.value(x), -1, atol = TOL)
    @test isapprox(JuMP.value(y), -sqrt(3), atol = TOL)
    @test isapprox(JuMP.value(z), 2, atol = TOL)

    JuMP.unset_integer(x)
    JuMP.optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    opt_obj = -2 * sqrt(5)
    @test isapprox(JuMP.objective_value(m), opt_obj, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), opt_obj, atol = TOL)
    @test isapprox(abs2(JuMP.value(x)) + abs2(JuMP.value(y)), 4, atol = TOL)
    @test isapprox(JuMP.value(z), 2, atol = TOL)
    return
end

function epinormeucl2(opt)
    TOL = 1e-4
    m = JuMP.Model(opt)

    JuMP.@variable(m, x[1:3], Bin)
    K = Hypatia.EpiNormEuclCone{Float64}(4)
    c1 = JuMP.@constraint(m, vcat(inv(sqrt(2)), x .- 0.5) in K)
    JuMP.optimize!(m)
    @test JuMP.termination_status(m) == MOI.INFEASIBLE
    @test JuMP.primal_status(m) == MOI.NO_SOLUTION

    JuMP.@objective(m, Min, 1 + sum(x))
    JuMP.delete(m, c1)
    JuMP.@constraint(m, vcat(1.5, x .- 1) in Hypatia.EpiNormEuclCone{Float64}(4))
    JuMP.set_start_value.(x, [1, 0, 0])
    JuMP.optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    @test isapprox(JuMP.objective_value(m), 2, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), 2, atol = TOL)
    return
end

function epipersquare1(opt)
    m = JuMP.Model(opt)
    JuMP.@variable(m, x[1:3], Int)
    K = Hypatia.EpiPerSquareCone{Float64}(5)
    JuMP.@constraint(m, vcat(0.25, 1, x .- 0.5) in K)
    JuMP.optimize!(m)
    @test JuMP.termination_status(m) == MOI.INFEASIBLE
    @test JuMP.primal_status(m) == MOI.NO_SOLUTION
    return
end

function epipersquare2(opt)
    TOL = 1e-4
    m = JuMP.Model(opt)

    JuMP.@variable(m, v, Bin)
    JuMP.@variable(m, w[1:2], Int)
    JuMP.@constraint(m, vcat(1, v, w .- 1.5) in Hypatia.EpiPerSquareCone{Float64}(4))
    JuMP.@objective(m, Min, v)
    JuMP.optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    @test isapprox(JuMP.objective_value(m), 1, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), 1, atol = TOL)
    @test isapprox(JuMP.value(v), 1, atol = TOL)

    JuMP.unset_integer.(w)
    JuMP.@objective(m, Min, v + sum(w) - 1)
    JuMP.optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    @test isapprox(JuMP.objective_value(m), 1, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), 1, atol = TOL)
    @test isapprox(JuMP.value.(w), [0.5, 0.5], atol = 10TOL)
    return
end

function epinorminf1(opt)
    TOL = 1e-4
    for use_dual in (false, true), is_complex in (false, true)
        R = (is_complex ? ComplexF64 : Float64)
        u_fix = (is_complex ? rt2 : 1.0) * (use_dual ? 0.5 : 0.4)
        w_fix = (is_complex ? [0.3, 0.3, 0.7, 0.7, 1.2, 1.2] : [0.3, 0.7, 1.2])
        w_opt = round.(w_fix)
        dim = vec_length(R, 3)
        m = JuMP.Model(opt)

        JuMP.@variable(m, w[1:dim], Bin)
        K = Hypatia.EpiNormInfCone{Float64, R}(1 + dim, use_dual)
        c1 = JuMP.@constraint(m, vcat(u_fix, w .- 0.5) in K)
        JuMP.optimize!(m)
        @test JuMP.termination_status(m) == MOI.INFEASIBLE
        @test JuMP.primal_status(m) == MOI.NO_SOLUTION

        JuMP.@objective(m, Min, 1 + sum(w))
        JuMP.delete(m, c1)
        K = Hypatia.EpiNormInfCone{Float64, R}(1 + dim, use_dual)
        if use_dual
            u_fix *= 2
        end
        JuMP.@constraint(m, vcat(u_fix, w - w_fix) in K)
        JuMP.set_start_value.(w, w_opt)
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.objective_value(m), 1 + sum(w_opt), atol = TOL)
        @test isapprox(JuMP.objective_bound(m), 1 + sum(w_opt), atol = TOL)
        @test isapprox(JuMP.value.(w), w_opt, atol = TOL)
    end
    return
end

function epinorminf2(opt)
    TOL = 1e-4
    for use_dual in (false, true), is_complex in (false, true)
        R = (is_complex ? ComplexF64 : Float64)
        dim = vec_length(R, 2)
        m = JuMP.Model(opt)

        JuMP.@variable(m, w[1:dim], Int)
        K = Hypatia.EpiNormInfCone{Float64, R}(1 + dim, use_dual)
        aff = [i * w[i] for i in eachindex(w)]
        JuMP.@constraint(m, vcat(10.5, aff) in K)
        JuMP.@objective(m, Min, sum(w))
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        w_opt = JuMP.value.(w)
        @test isapprox(JuMP.objective_value(m), sum(w_opt), atol = TOL)
        @test isapprox(JuMP.objective_bound(m), sum(w_opt), atol = TOL)
        aff_opt = vec_copyto!(zeros(R, 2), JuMP.value.(aff))
        @test norm(aff_opt, (use_dual ? 1 : Inf)) <= 10.5 + TOL
    end
    return
end

function epinormspectraltri1(opt)
    TOL = 1e-4
    d = 3
    for use_dual in (false, true), is_complex in (false, true)
        R = (is_complex ? ComplexF64 : Float64)
        dim = svec_length(R, d)
        vec = collect(1:dim)
        mat = smat(R, fill(0.1, dim))
        σ = eigvals(Hermitian(mat, :U))
        opt_val = (use_dual ? sum(abs, σ) : maximum(abs, σ))
        m = JuMP.Model(opt)

        JuMP.@variable(m, y)
        x = JuMP.@variable(m, [1:dim], Int)
        JuMP.@constraint(m, [i in 1:dim], i - 1 <= x[i] <= i + 1)
        K = Hypatia.EpiNormSpectralTriCone{Float64, R}(1 + dim, use_dual)
        x_vec = scale_svec(R, vec - x, rt2)
        JuMP.@constraint(m, vcat(y, x_vec .+ 0.1) in K)
        JuMP.@objective(m, Min, y)
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(JuMP.value.(x), vec, atol = TOL)
    end
    return
end

function epinormspectraltri2(opt)
    TOL = 1e-4
    d = 3
    for use_dual in (false, true), is_complex in (false, true)
        R = (is_complex ? ComplexF64 : Float64)
        dim = svec_length(R, d)
        opt_val = (use_dual ? 1 : d)
        m = JuMP.Model(opt)

        x = JuMP.@variable(m, [1:dim], Bin)
        K = Hypatia.EpiNormSpectralTriCone{Float64, R}(1 + dim, use_dual)
        x_vec = scale_svec(R, 1.0 * x, rt2)
        JuMP.@constraint(m, vcat(1, x_vec) in K)
        JuMP.@objective(m, Max, sum(x))
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(sum(JuMP.value.(x)), opt_val, atol = TOL)
    end
    return
end

function epinormspectral1(opt)
    TOL = 1e-4
    d1 = 2
    d2 = 3
    for use_dual in (false, true), is_complex in (false, true)
        R = (is_complex ? ComplexF64 : Float64)
        dim = vec_length(R, d1 * d2)
        vec = collect(1:dim)
        mat = vec_copyto!(zeros(R, d1, d2), fill(0.1, dim))
        σ = svdvals(mat)
        opt_val = (use_dual ? sum(σ) : σ[1])
        m = JuMP.Model(opt)

        JuMP.@variable(m, y)
        x = JuMP.@variable(m, [1:dim], Int)
        JuMP.@constraint(m, [i in 1:dim], i - 1 <= x[i] <= i + 1)
        K = Hypatia.EpiNormSpectralCone{Float64, R}(d1, d2, use_dual)
        JuMP.@constraint(m, vcat(y, vec - x .+ 0.1) in K)
        JuMP.@objective(m, Min, y)
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(JuMP.value.(x), vec, atol = TOL)
    end
    return
end

function epinormspectral2(opt)
    TOL = 1e-4
    d1 = d2 = 2
    for use_dual in (false, true), is_complex in (false, true)
        R = (is_complex ? ComplexF64 : Float64)
        dim = vec_length(R, d1 * d2)
        opt_val = (use_dual ? 1 : d1)
        m = JuMP.Model(opt)

        x = JuMP.@variable(m, [1:dim], Bin)
        K = Hypatia.EpiNormSpectralCone{Float64, R}(d1, d2, use_dual)
        JuMP.@constraint(m, vcat(1, x) in K)
        JuMP.@objective(m, Max, sum(x))
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(sum(JuMP.value.(x)), opt_val, atol = TOL)
    end
    return
end

function hypogeomean1(opt)
    TOL = 1e-4
    m = JuMP.Model(opt)

    JuMP.@variable(m, u)
    JuMP.@variable(m, w[1:3], Int)
    JuMP.@constraint(m, vcat(u, w) in Hypatia.HypoGeoMeanCone{Float64}(4))
    JuMP.@objective(m, Max, u)
    w_max = [1.1, 2.3, 3.5]
    JuMP.@constraint(m, w .<= w_max)
    JuMP.optimize!(m)

    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    opt_obj = PajaritoExtras.geomean(Float64.(1:3))
    @test isapprox(JuMP.objective_value(m), opt_obj, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), opt_obj, atol = TOL)
    @test isapprox(JuMP.value(u), opt_obj, atol = TOL)
    @test isapprox(JuMP.value.(w), 1:3, atol = TOL)

    JuMP.unset_integer.(w)
    JuMP.optimize!(m)
    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    opt_obj = PajaritoExtras.geomean(w_max)
    @test isapprox(JuMP.objective_value(m), opt_obj, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), opt_obj, atol = TOL)
    @test isapprox(JuMP.value.(w), w_max, atol = TOL)
    return
end

function hypogeomean2(opt)
    TOL = 1e-4
    m = JuMP.Model(opt)
    JuMP.@variable(m, x, Bin)
    JuMP.@variable(m, y[1:3], Int)
    JuMP.@constraint(m, vcat(0.5 * x, y) in Hypatia.HypoGeoMeanCone{Float64}(4))
    JuMP.@objective(m, Max, 5x - sum(y))
    JuMP.optimize!(m)

    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    @test isapprox(JuMP.objective_value(m), 2, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), 2, atol = TOL)
    @test isapprox(JuMP.value(x), 1, atol = TOL)
    @test isapprox(JuMP.value.(y), ones(3), atol = TOL)
    return
end

function hyporootdettri1(opt)
    TOL = 1e-4
    for d in (2, 3), is_complex in (false, true)
        R = (is_complex ? ComplexF64 : Float64)
        w_dim = svec_length(R, d)
        m = JuMP.Model(opt)

        w = JuMP.@variable(m, [1:w_dim], Bin)
        K = Hypatia.HypoRootdetTriCone{Float64, R}(1 + w_dim)
        JuMP.@constraint(m, vcat(1, w) in K)
        JuMP.@objective(m, Max, sum(w))
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.objective_value(m), d, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), d, atol = TOL)
        w_val = JuMP.value.(w)
        w_mat = Hermitian(smat(R, w_val), :U)
        @test isapprox(w_mat, I, atol = TOL)

        JuMP.@constraint(m, sum(w) <= d - 1)
        JuMP.optimize!(m)
        @test JuMP.termination_status(m) == MOI.INFEASIBLE
        @test JuMP.primal_status(m) == MOI.NO_SOLUTION
    end
    return
end

function hyporootdettri2(opt)
    TOL = 1e-4
    (m, x, Q, opt_x, opt_Q) = _setup_expdesign(opt)
    opt_val = sqrt(det(opt_Q))
    K = Hypatia.HypoRootdetTriCone{Float64, Float64}(4)

    JuMP.@variable(m, y)
    JuMP.@objective(m, Max, y)
    JuMP.@constraint(m, vcat(y, svec(Q)) in K)
    JuMP.optimize!(m)

    @test JuMP.termination_status(m) == MOI.OPTIMAL
    @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
    @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
    @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
    @test isapprox(JuMP.value.(x), opt_x, atol = TOL)
    return
end

function vector_epipersepspectral1(opt)
    TOL = 1e-4
    sep_spectral_funs = [
        Hypatia.Cones.NegLogSSF(),
        Hypatia.Cones.NegEntropySSF(),
        Hypatia.Cones.NegSqrtSSF(),
        Hypatia.Cones.NegPower01SSF(0.3),
        Hypatia.Cones.Power12SSF(3 // 2),
    ]

    for h_fun in sep_spectral_funs, use_dual in (false, true)
        D = (use_dual ? Dual : Prim)
        Q = Hypatia.Cones.VectorCSqr{Float64}
        K = Hypatia.EpiPerSepSpectralCone{Float64}(h_fun, Q, 3, use_dual)
        m = JuMP.Model(opt)

        JuMP.@variable(m, u)
        JuMP.@variable(m, w[1:3], Int)
        JuMP.@constraint(m, w .<= [3.2, 5.5, 6.9])
        JuMP.@constraint(m, sum(w) == 8)
        JuMP.@objective(m, Min, u)
        epiper = PajaritoExtras.swap_epiper(D, [u, 1]...)
        JuMP.@constraint(m, vcat(epiper..., w) in K)
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        w_sol = JuMP.value.(w)
        opt_val = PajaritoExtras.val_or_conj(D)(w_sol, h_fun)
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(JuMP.value(u), opt_val, atol = TOL)
        @test isapprox(w_sol, round.(w_sol), atol = TOL)
    end
    return
end

function vector_epipersepspectral2(opt)
    TOL = 1e-4
    # only functions for which the perspective is zero when the perspective variable is zero
    sep_spectral_funs = [
        Hypatia.Cones.NegLogSSF(),
        Hypatia.Cones.NegSqrtSSF(),
        Hypatia.Cones.NegPower01SSF(3 // 10),
    ]

    for h_fun in sep_spectral_funs, use_dual in (false, true)
        D = (use_dual ? Dual : Prim)
        Q = Hypatia.Cones.VectorCSqr{Float64}
        w = [0.5, 2.3]
        m = JuMP.Model(opt)

        JuMP.@variable(m, x[1:2], Bin)
        JuMP.@constraint(m, sum(x) == 1)
        JuMP.@variable(m, u[1:2])
        JuMP.@objective(m, Min, sum(u))
        for i in 1:2
            K = Hypatia.EpiPerSepSpectralCone{Float64}(h_fun, Q, 2, use_dual)
            epiper = PajaritoExtras.swap_epiper(D, u[i], x[i])
            JuMP.@constraint(m, vcat(epiper..., w .+ i) in K)
        end
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        x_sol = JuMP.value.(x)
        round_x_sol = round.(x_sol)
        @test isapprox(x_sol, round_x_sol, atol = TOL)
        @test round_x_sol in ([0.0, 1.0], [1.0, 0.0])
        i = findfirst(Bool.(round_x_sol))
        u_sol = JuMP.value.(u)
        opt_val = PajaritoExtras.val_or_conj(D)(w .+ i, h_fun)
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(u_sol[i], opt_val, atol = TOL)
        @test isapprox(sum(u_sol), opt_val, atol = TOL)
    end
    return
end

function matrix_epipersepspectral1(opt)
    TOL = 1e-4
    sep_spectral_funs = [
        Hypatia.Cones.NegLogSSF(),
        Hypatia.Cones.NegEntropySSF(),
        Hypatia.Cones.Power12SSF(1.2),
    ]

    for h_fun in sep_spectral_funs, use_dual in (false, true), is_complex in (false, true)
        R = (is_complex ? ComplexF64 : Float64)
        D = (use_dual ? Dual : Prim)
        Q = Hypatia.Cones.MatrixCSqr{Float64, R}
        K = Hypatia.EpiPerSepSpectralCone{Float64}(h_fun, Q, 2, use_dual)
        w_dim = svec_length(R, 2)
        vec_diag = svec(Matrix{R}(Diagonal([1.2, 2.5])))
        m = JuMP.Model(opt)

        JuMP.@variable(m, u)
        JuMP.@variable(m, w[1:w_dim])
        W_diag = [w[svec_idx(R, i, i)] for i in 1:2]
        JuMP.set_integer.(W_diag)
        JuMP.@constraint(m, sum(W_diag) == 6)
        JuMP.@constraint(m, sum(w) >= 7)
        JuMP.@objective(m, Min, u)
        K_vec = w - vec_diag
        epiper = PajaritoExtras.swap_epiper(D, [u, 1]...)
        JuMP.@constraint(m, vcat(epiper..., K_vec) in K)
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        diag_sol = JuMP.value.(W_diag)
        @test isapprox(diag_sol, round.(diag_sol), atol = TOL)
        ω = eigvals(Hermitian(smat(R, JuMP.value.(K_vec)), :U))
        if PajaritoExtras.dom_pos(D, h_fun)
            @test minimum(ω) > -TOL
            ω = max.(ω, 1e-9)
        end
        opt_val = PajaritoExtras.val_or_conj(D)(ω, h_fun)
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(JuMP.value(u), opt_val, atol = TOL)
    end
    return
end

function matrix_epipersepspectral2(opt)
    TOL = 1e-4
    # only functions that are decreasing
    dual_and_h = [
        (false, Hypatia.Cones.NegSqrtSSF()),
        (false, Hypatia.Cones.NegPower01SSF(3 // 10)),
        (true, Hypatia.Cones.NegEntropySSF()),
        (true, Hypatia.Cones.NegSqrtSSF()),
        (true, Hypatia.Cones.NegPower01SSF(0.5)),
    ]

    for (use_dual, h_fun) in dual_and_h
        D = (use_dual ? Dual : Prim)
        Q = Hypatia.Cones.MatrixCSqr{Float64, Float64}
        K = Hypatia.EpiPerSepSpectralCone{Float64}(h_fun, Q, 2, use_dual)
        (m, x, Q, opt_x, opt_Q) = _setup_expdesign(opt)
        opt_val = PajaritoExtras.val_or_conj(D)(eigvals(Hermitian(opt_Q, :U)), h_fun)

        JuMP.@variable(m, y)
        JuMP.@objective(m, Min, y)
        epiper = PajaritoExtras.swap_epiper(D, [y, 1]...)
        JuMP.@constraint(m, vcat(epiper..., svec(Q)) in K)
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(JuMP.value.(x), opt_x, atol = TOL)
    end
    return
end

function _setup_expdesign(opt)
    m = JuMP.Model(opt)
    JuMP.@variable(m, x[1:4], Int)
    JuMP.@constraint(m, x[1:2] .>= 0.5) # avoids ill-posedness
    JuMP.@constraint(m, x[3:4] .>= 0)
    JuMP.@constraint(m, sum(x) <= 8.5)

    V = [1 1 -0.2 -0.5; 1 -1 0.5 -0.2]
    Q = V * diagm(x) * V'
    opt_x = [4, 4, 0, 0]
    opt_Q = Hermitian(V * Diagonal(opt_x) * V')
    return (m, x, Q, opt_x, opt_Q)
end

function wsosinterpnonnegative1(opt)
    TOL = 1e-4
    # real case
    dom = PolyUtils.BoxDomain{Float64}(-ones(1), ones(1))
    fs = Function[
        x -> (x - 0.5)^2 + 0.3, # min 0.3
        x -> 2x + 1.2, # min -0.8
        x -> -1 - x, # min -2
        x -> 3 - 2x^2, # min 1
    ]
    (U, pts, Ps) = PolyUtils.interpolate(dom, 2)
    pts = vec(pts)
    f_pts = [f_i.(pts) for f_i in fs]

    for use_dual in (false, true)
        opt_val = (use_dual ? -0.8 : 0.3)
        K = Hypatia.WSOSInterpNonnegativeCone{Float64, Float64}(U, Ps, use_dual)
        (m, x, y) = _setup_polymin(opt, use_dual, K, f_pts, 3.0)
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.value.(x), [1, 0, 0, 1], atol = TOL)
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(JuMP.value(y), opt_val, atol = TOL)
    end
    return
end

function wsosinterpnonnegative2(opt)
    TOL = 1e-4
    # complex case
    dom = [x -> 1 - abs2(x[1])]
    fs = Function[
        x -> abs2(x - 0.5) - 0.4, # min -0.4
        x -> 2 * real(x) - imag(x) + 1, # min -1.2361
        x -> 2 * abs2(x) + imag(x), # min -0.125
        x -> real(x) + 0.1, # min -0.9
    ]
    (pts, Ps) = PolyUtils.interpolate(ComplexF64, 1, 1, dom, [1])
    pts = first.(pts)
    U = length(pts)
    f_pts = [f_i.(pts) for f_i in fs]

    for use_dual in (false, true)
        opt_val = (use_dual ? -0.9 : -0.4)
        K = Hypatia.WSOSInterpNonnegativeCone{Float64, ComplexF64}(U, Ps, use_dual)
        (m, x, y) = _setup_polymin(opt, use_dual, K, f_pts, 3.0)
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        @test isapprox(JuMP.value.(x), [1, 0, 1, 0], atol = TOL)
        @test isapprox(JuMP.objective_value(m), opt_val, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_val, atol = TOL)
        @test isapprox(JuMP.value(y), opt_val, atol = TOL)
    end
    return
end

function _setup_polymin(opt, use_dual, K, f_pts, bound)
    # pick 2/4 polynomials to either
    # maximize the minimum value (primal) or minimize the maximum value (dual)
    m = JuMP.Model(opt)
    JuMP.@variable(m, x[1:4], Bin)
    JuMP.@constraint(m, sum(x) == 2)
    JuMP.@variable(m, y)

    if use_dual
        JuMP.@objective(m, Min, y)
        U = MOI.dimension(K)
        for i in 1:4
            z_i = JuMP.@variable(m, [1:U])
            JuMP.@constraint(m, z_i in K)
            JuMP.@constraint(m, sum(z_i) == 1)
            JuMP.@constraint(m, y >= JuMP.dot(f_pts[i], z_i) - bound * x[i])
        end
    else
        JuMP.@objective(m, Max, y)
        for i in 1:4
            JuMP.@constraint(m, f_pts[i] .+ (bound * (1 - x[i]) - y) in K)
        end
    end
    return (m, x, y)
end

function specialorderedset(opt)
    TOL = 1e-4
    for (S, sol) in [(MOI.SOS1, [1, 0, 0]), (MOI.SOS2, [1, 1, 0])]
        if !MOI.supports_constraint(MOI.instantiate(opt), MOI.VectorOfVariables, S{Float64})
            # only test if MOIPajarito supports SOS1/2
            continue
        end
        m = JuMP.Model(opt)

        JuMP.@variable(m, u)
        JuMP.@variable(m, w[1:3])
        w0 = w + [1.1, 2.3, 3.5]
        JuMP.@constraint(m, vcat(u, w0) in Hypatia.HypoGeoMeanCone{Float64}(4))
        JuMP.@objective(m, Max, u)
        JuMP.@constraint(m, w .<= 1.5)
        JuMP.@constraint(m, w in S([1.0, 2, 3]))
        JuMP.optimize!(m)

        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        opt_obj = PajaritoExtras.geomean(JuMP.value.(w0))
        @test isapprox(JuMP.objective_value(m), opt_obj, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_obj, atol = TOL)
        @test isapprox(JuMP.value(u), opt_obj, atol = TOL)
        @test isapprox(JuMP.value.(w), 1.5 * sol, atol = TOL)

        JuMP.set_integer.(w)
        JuMP.optimize!(m)
        @test JuMP.termination_status(m) == MOI.OPTIMAL
        @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        opt_obj = PajaritoExtras.geomean(JuMP.value.(w0))
        @test isapprox(JuMP.objective_value(m), opt_obj, atol = TOL)
        @test isapprox(JuMP.objective_bound(m), opt_obj, atol = TOL)
        @test isapprox(JuMP.value(u), opt_obj, atol = TOL)
        @test isapprox(JuMP.value.(w), sol, atol = TOL)
    end
    return
end


text:
using SimplePosets
using SimpleGraphs: SimpleGraph, AbstractSimpleGraph, SimpleDigraph, vertex_type, NV, elist, in_neighbors, add_edges!

import Graphs
import Graphs.simple_graph, Graphs.add_edge!, Graphs.topological_sort_by_dfs
export convert_simple

# From https://github.com/scheinerman/SimpleGraphs.jl/blob/1396758729f95912d7f245f9c70957f4993be417/src/simple_converters.jl
function convert_simple(G::AbstractSimpleGraph)
    T = vertex_type(G)
    n = NV(G)
    has_dir = isa(G,SimpleDigraph)

    d = SimpleGraphs.vertex2idx(G)
    dinv = Dict{Int,T}()
    for k in keys(d)
        v = d[k]
        dinv[v] = k
    end

    H = simple_graph(n,is_directed=has_dir)

    EE = elist(G)
    for e in EE
        u = d[e[1]]
        v = d[e[2]]
        add_edge!(H,u,v)
    end
    return (H,d,dinv)
end

export toposort
function toposort(m::Model)
    (g, _, names) = poset(m).D |> convert_simple
    setdiff(map(v -> names[v], Graphs.topological_sort_by_dfs(g)), arguments(m))
end


text:
# # Double Gyre

# ## Model

module doublegyre   #jl

using ReachabilityAnalysis

@taylorize function doublegyre!(dx, x, params, t)
    local A = 0.1
    dx[1] = -π*A*sin(π*x[1])*cos(π*x[2])
    dx[2] = π*A*cos(π*x[1])*sin(π*x[2])
    return dx
end

function model(X0)
    S = @system(x' = doublegyre!(x), dim:2)
    return IVP(S, X0)
end

end #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
import VoronoiCells: LEFT, RIGHT, LOWER, UPPER
using VoronoiDelaunay

@testset "Area of polygons" begin
	# Area of square
	S = [Point2D(LEFT, LOWER) ; Point2D(RIGHT, LOWER) ; Point2D(RIGHT, UPPER) ; Point2D(LEFT, UPPER)]
	@test polyarea(S) ≈ 1.0

	# Area of triangle
	T1 = [Point2D(LEFT, LOWER) ; Point2D(RIGHT, LOWER) ; 0.5*(Point2D(RIGHT, UPPER) + Point2D(LEFT, UPPER))]
	@test polyarea(T1) ≈ 0.5

	T2 = [Point2D(LEFT, LOWER) ; Point2D(LEFT+0.5, LOWER) ;
	Point2D(LEFT, LOWER+0.5)]
	@test polyarea(T2) ≈ 0.5^3
end


text:
# 2D FWI on a small version of the Overthrust model with SPG from the minConf library
# Author: pwitte.slim@gmail.com
# Date: December 2018
#

using Statistics, Random, LinearAlgebra
using JUDI, SlimOptim, HDF5, SegyIO, PyPlot

# Load starting model
n,d,o,m0 = read(h5open("$(JUDI.JUDI_DATA)/small_overthrust_model.h5","r"), "n", "d", "o", "m0")
model0 = Model((n[1],n[2]), (d[1],d[2]), (o[1],o[2]), m0)

# Bound constraints
v0 = sqrt.(1f0 ./ model0.m)
vmin = ones(Float32,model0.n) .* 1.3f0
vmax = ones(Float32,model0.n) .* 6.5f0
vmin[:,1:21] .= v0[:,1:21]   # keep water column fixed
vmax[:,1:21] .= v0[:,1:21]

# Slowness squared [s^2/km^2]
mmin = vec((1f0 ./ vmax).^2)
mmax = vec((1f0 ./ vmin).^2)

# Load data
block = segy_read("$(JUDI.JUDI_DATA)/small_overthrust_shot_records.segy")
d_obs = judiVector(block)

# Set up wavelet
src_geometry = Geometry(block; key="source")
wavelet = ricker_wavelet(src_geometry.t[1],src_geometry.dt[1],0.008f0)  # 8 Hz wavelet
q = judiVector(src_geometry,wavelet)

############################### FWI ###########################################

# Optimization parameters
fevals = 10
batchsize = 8
fvals = []

# Objective function for minConf library
function objective_function(x)
    model0.m .= x

    # fwi function value and gradient
    i = randperm(d_obs.nsrc)[1:batchsize]
    fval, grad = fwi_objective(model0, q[i], d_obs[i])
    grad = .125f0*grad/norm(grad, Inf)  # scale for line search

    global fvals; fvals = [fvals; fval]
    return fval, vec(grad.data)
end

# Bound projection
ProjBound(x) = median([mmin x mmax]; dims=2)[1:end]

# FWI with SPG
options = spg_options(verbose=3, maxIter=fevals, memory=3)
x, fsave, funEvals= spg(objective_function, vec(m0), ProjBound, options)

# Save results
h5open("result_2D_small_overthrust_spg.h5", "w") do file
    write(file, "x", sqrt.(1f0 ./ reshape(x, model0.n)), "fsave", fsave, "fhistory", convert(Array{Float32, 1}, fvals))
end

# Plot convergence and final result
figure(); plot(fvals/norm(fvals, Inf));
xlabel("Iteration no."); ylabel("Normalized residual"); title("Convergence of FWI w/ SPG")
figure(); imshow(sqrt.(1f0 ./ adjoint(reshape(x, model0.n))), extent=(0, 20.0, 5.15, 0)); title("FWI with SPG");
xlabel("Lateral position [km]"); ylabel("Depth [km]")


text:
function background_θ(z, p::Params)
    return p.θ0 + (p.π0 - z)
end

function background_ϕ(z, p::Params)
    θ = background_θ(z, p)
    return 0.5*(p.θ0 + θ)*(p.π0 - z)
end

function background_ψ(z, p::Params)
    return background_ϕ(z, p)
end


text:
var = 7
if var > 10
    println("var has value $var and is bigger than 10.")
elseif var < 10    # This elseif clause is optional.
    println("var has value $var and is smaller than 10.")
else                    # The else clause is optional too.
    println("var has value $var and is 10.")
end
#> prints "var has value 7 and is smaller than 10."

# if 0 print("ok") end
# ERROR: LoadError: TypeError: non-boolean (Int64) used in boolean context
a = 10
b = 15
z = if a > b a
    else  b
    end
println(z) #> 15
z = a > b ? a : b # with ternary operator

var = 7
varout = "var has value $var"
cond = var > 10 ? "and is bigger than 10." : var < 10 ? "and is smaller than 10" : "and is 10."
println("$varout $cond") # var has value 7 and is smaller than 10

function sqroot(n::Int)
    n >= 0 || error("n must be non-negative")
    n == 0 && return 0
    sqrt(n)
end
sqroot(4) #> 2.0
sqroot(0) #> 0.0
# sqroot(-6) #> ERROR: LoadError: n must be non-negative


text:
@testset "function defs" begin
    @test CSTParser.defines_function(CSTParser.parse("function f end"))
    @test CSTParser.defines_function(CSTParser.parse("function f() end"))
    @test CSTParser.defines_function(CSTParser.parse("function f()::T end"))
    @test CSTParser.defines_function(CSTParser.parse("function f(x::T) where T end"))
    @test CSTParser.defines_function(CSTParser.parse("function f{T}() end"))
    @test CSTParser.defines_function(CSTParser.parse("f(x) = x"))
    @test CSTParser.defines_function(CSTParser.parse("f(x)::T = x"))
    @test CSTParser.defines_function(CSTParser.parse("f{T}(x)::T = x"))
    @test CSTParser.defines_function(CSTParser.parse("f{T}(x)::T = x"))
    @test CSTParser.defines_function(CSTParser.parse("*(x,y) = x"))
    @test CSTParser.defines_function(CSTParser.parse("*(x,y)::T = x"))
    @test CSTParser.defines_function(CSTParser.parse("!(x::T)::T = x"))
    @test CSTParser.defines_function(CSTParser.parse("a + b = a"))
    @test CSTParser.defines_function(CSTParser.parse("a/b = x"))
    @test !CSTParser.defines_function(CSTParser.parse("a.b = x"))
end

@testset "datatype defs" begin
    @test CSTParser.defines_struct(CSTParser.parse("struct T end"))
    @test CSTParser.defines_struct(CSTParser.parse("mutable struct T end"))
    @test CSTParser.defines_mutable(CSTParser.parse("mutable struct T end"))
    @test CSTParser.defines_abstract(CSTParser.parse("abstract type T end"))
    # @test CSTParser.defines_abstract(CSTParser.parse("abstract T"))
    @test CSTParser.defines_primitive(CSTParser.parse("primitive type a b end"))
end

@testset "get_name" begin
    @test valof(CSTParser.get_name(CSTParser.parse("struct T end"))) == "T"
    @test valof(CSTParser.get_name(CSTParser.parse("struct T{T} end"))) == "T"
    @test valof(CSTParser.get_name(CSTParser.parse("struct T <: T end"))) == "T"
    @test valof(CSTParser.get_name(CSTParser.parse("struct T{T} <: T end"))) == "T"

    @test valof(CSTParser.get_name(CSTParser.parse("mutable struct T end"))) == "T"
    @test valof(CSTParser.get_name(CSTParser.parse("mutable struct T{T} end"))) == "T"
    @test valof(CSTParser.get_name(CSTParser.parse("mutable struct T <: T end"))) == "T"
    @test valof(CSTParser.get_name(CSTParser.parse("mutable struct T{T} <: T end"))) == "T"

    @test valof(CSTParser.get_name(CSTParser.parse("abstract type T end"))) == "T"
    @test valof(CSTParser.get_name(CSTParser.parse("abstract type T{T} end"))) == "T"
    @test valof(CSTParser.get_name(CSTParser.parse("abstract type T <: T end"))) == "T"
    @test valof(CSTParser.get_name(CSTParser.parse("abstract type T{T} <: T end"))) == "T"
    # NEEDS FIX: v0.6 dep
    # @test CSTParser.get_name(CSTParser.parse("abstract T")).val == "T"
    # @test CSTParser.get_name(CSTParser.parse("abstract T{T}")).val == "T"
    # @test CSTParser.get_name(CSTParser.parse("abstract T <: T")).val == "T"
    # @test CSTParser.get_name(CSTParser.parse("abstract T{T} <: T")).val == "T"

    @test valof(CSTParser.get_name(CSTParser.parse("function f end"))) == "f"
    @test valof(CSTParser.get_name(CSTParser.parse("function f() end"))) == "f"
    @test valof(CSTParser.get_name(CSTParser.parse("function f()::T end"))) == "f"
    @test valof(CSTParser.get_name(CSTParser.parse("function f(x::T) where T end"))) == "f"
    @test valof(CSTParser.get_name(CSTParser.parse("function f{T}() end"))) == "f"

    # Operators
    @test CSTParser.str_value(CSTParser.get_name(CSTParser.parse("function +() end"))) == "+"
    @test CSTParser.str_value(CSTParser.get_name(CSTParser.parse("function (+)() end"))) == "+"
    @test CSTParser.str_value(CSTParser.get_name(CSTParser.parse("+(x,y) = x"))) == "+"
    @test CSTParser.str_value(CSTParser.get_name(CSTParser.parse("+(x,y)::T = x"))) == "+"
    @test CSTParser.str_value(CSTParser.get_name(CSTParser.parse("!(x)::T = x"))) == "!"
    @test CSTParser.str_value(CSTParser.get_name(CSTParser.parse("!(x) = x"))) == "!"
end

@testset "get_sig_params" begin
    @test CSTParser.get_where_params(CSTParser.parse("a where T")) == ["T"]
    @test CSTParser.get_where_params(CSTParser.parse("a where {T}")) == ["T"]
    @test CSTParser.get_where_params(CSTParser.parse("a where {T,S}")) == ["T", "S"]
    @test CSTParser.get_where_params(CSTParser.parse("a where T <: S")) == ["T"]
    @test CSTParser.get_where_params(CSTParser.parse("a where T{S}")) == ["T"]
    @test CSTParser.get_where_params(CSTParser.parse("a where T{S} <: R")) == ["T"]
    @test CSTParser.get_sig_params(CSTParser.parse("f() where T")) == ["T"]
    @test CSTParser.get_sig_params(CSTParser.parse("f{T}()")) == ["T"]
    @test CSTParser.get_sig_params(CSTParser.parse("f{T}() where S")) == ["S", "T"]
    @test CSTParser.get_sig_params(CSTParser.parse("f{T}() where S{R} where J")) == ["J", "S", "T"]
end

@testset "get_sig_params" begin
    f = x -> CSTParser.str_value.(CSTParser.get_args(CSTParser.parse(x)))
    @test f("function f(a) end") == ["a"]
    @test f("function f(a::T) end") == ["a"]
    @test f("function f(a,b) end") == ["a", "b"]
    @test f("function f(a::T,b::T) end") == ["a", "b"]
    @test f("function f(a::T,b::T) where T end") == ["a", "b"]
    @test f("function f{T}(a::T,b::T) where T end") == ["a", "b"]
    @test f("function f{T}(a::T,b::T;c = 1) where T end") == ["a", "b", "c"]

    @test f("a -> a") == ["a"]
    @test f("a::T -> a") == ["a"]
    @test f("(a::T) -> a") == ["a"]
    @test f("(a,b) -> a") == ["a", "b"]

    @test f("map(1:10) do a
        a
    end") == ["a"]
    @test f("map(1:10) do a,b
        a
    end") == ["a", "b"]
end


text:
# Note that the output files are not true JSON. Each file must be manually edited to be a JSON object.

# Small values:
stream = open( "bits_1e-200_1e-308.json", "w" );

for i in linspace( 1e-200, 1e-308, 1007 )
    write( stream, string( "\"", bits( i ), "\":\"", i, "\",\n" ) );
end

close( stream );

# Medium values:
stream = open( "bits_-1e3_1e3.json", "w" );

for i in linspace( -1e3, 1e3, 1007 )
    write( stream, string( "\"", bits( i ), "\":\"", i, "\",\n" ) );
end

close( stream );

# Large values:
stream = open( "bits_1e200_1e308.json", "w" );

for i in linspace( 1e200, 1e308, 1007 )
    write( stream, string( "\"", bits( i ), "\":\"", i, "\",\n" ) );
end

close( stream );

# Subnormal values:
stream = open( "bits_1e-310_5e-324.json", "w" );

for i in linspace( 1e-310, 5e-324, 1007 )
    write( stream, string( "\"", bits( i ), "\":\"", i, "\",\n" ) );
end

close( stream );


text:
# Minimal CxxWrap dummy module for the tests to run
module CxxWrap

# Base type for wrapped C++ types
abstract type CppAny end
abstract type CppBits <: CppAny end
abstract type CppDisplay <: AbstractDisplay end
abstract type CppArray{T,N} <: AbstractArray{T,N} end
abstract type CppAssociative{K,V} <: AbstractDict{K,V} end
abstract type CppEnum end
Base.convert(::Type{Int32}, x::CppEnum) = reinterpret(Int32, x)
import Base: +, |
+(a::T, b::T) where {T <: CppEnum} = reinterpret(T, Int32(a) + Int32(b))
|(a::T, b::T) where {T <: CppEnum} = reinterpret(T, Int32(a) | Int32(b))
cxxdowncast(x) = error("No downcast for type $(supertype(typeof(x))). Did you specialize SuperType to enable automatic downcasting?")
abstract type SmartPointer{T} <: CppAny end

mutable struct CppFunctionInfo
  name::Any
  argument_types::Array{Type,1}
  reference_argument_types::Array{Type,1}
  return_type::Type
  function_pointer::Int
  thunk_pointer::Int
end

mutable struct ConstructorFname
  _type::DataType
end

end

max_stars_repo_name: RelationalAI-oss/libcxxwrap-julia

max_issues_repo_name: RelationalAI-oss/libcxxwrap-julia

max_forks_repo_name: RelationalAI-oss/libcxxwrap-julia


text:
#=
    Cheating model.

    The model is a port of the Cheating model (the first version)
    https://github.com/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Chapter2_MorePyMC/Ch2_MorePyMC_PyMC3.ipynb
    """
    Social data has an additional layer of interest as people are not always honest with responses, which
    adds a further complication into inference. For example, simply asking individuals "Have you ever cheated
    on a test?" will surely contain some rate of dishonesty. What you can say for certain is that the true
    rate is less than your observed rate (assuming individuals lie only about not cheating; I cannot
    imagine one who would admit "Yes" to cheating when in fact they hadn't cheated).

    To present an elegant solution to circumventing this dishonesty problem, and to demonstrate Bayesian
    modeling, we first need to introduce the binomial distribution.

    ...

    We will use the binomial distribution to determine the frequency of students cheating during an exam.
    If we let $N$ be the total number of students who took the exam, and assuming each student is interviewed
    post-exam (answering without consequence), we will receive integer $X$ "Yes I did cheat" answers. We
    then find the posterior distribution of $p$, given $N$, some specified prior on $p$, and observed data $X$.

    This is a completely absurd model. No student, even with a free-pass against punishment, would admit to
    cheating. What we need is a better algorithm to ask students if they had cheated. Ideally the algorithm
    should encourage individuals to be honest while preserving privacy. The following proposed algorithm
    is a solution I greatly admire for its ingenuity and effectiveness:

        In the interview process for each student, the student flips a coin, hidden from the interviewer.
        The student agrees to answer honestly if the coin comes up heads. Otherwise, if the coin comes
        up tails, the student (secretly) flips the coin again, and answers "Yes, I did cheat" if the
        coin flip lands heads, and "No, I did not cheat", if the coin flip lands tails. This way, the
        interviewer does not know if a "Yes" was the result of a guilty plea, or a Heads on a second
        coin toss. Thus privacy is preserved and the researchers receive honest answers.
    """

    This model:

    Summary Statistics
    parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

            p    0.2273    0.1249     0.0040    0.0087   290.4227    1.0058      121.5666
     obs_prop    0.3559    0.0460     0.0015    0.0023   403.9539    0.9991      169.0891

    Quantiles
    parameters      2.5%     25.0%     50.0%     75.0%     97.5%
        Symbol   Float64   Float64   Float64   Float64   Float64

            p    0.0255    0.1364    0.2166    0.3005    0.4916
     obs_prop    0.2700    0.3200    0.3500    0.3800    0.4500

   Cf cheating2.jl for a port of a variant of the PyMC3 model (also from the "Hacker" book).
=#

using Turing, StatsPlots, Distributions, StatsBase
using CSV
include("jl_utils.jl")

@model function cheating_model(N,X)
    p ~ Uniform(0,1) # freq cheating

    truths       ~ filldist(Bernoulli(p),N)
    first_flips  ~ filldist(Bernoulli(0.5), N)
    second_flips ~ filldist(Bernoulli(0.5), N)

    val = tzeros(N) # We must defined val to be used as LHS in @.
    @. val = first_flips*truths + (1 - first_flips)*second_flips
    obs_prop ~ Dirac(sum(val) / N)
    X ~ Binomial(N, obs_prop )

    p_prior ~ Uniform(0,1)

end

N = 100
yes_response = 35
model = cheating_model(N, yes_response)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(5),  1_000)
# chns = sample(model, PG(5),  MCMCThreads(), 10_000, 4)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns[[:p,:obs_prop,:p_prior]])
# display(plot(chns[[:p,:obs_prop,:p_prior]]))

# Another way:
println("mean(truths): $(mean(group(chns, :truths).value.data,dims=1) |> mean)")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank