text:
import MathematicalSets.AbstractSet

# ===========================================================================
# This is a collection of concrete subtypes of `AbstractSet` for the purposes
# of testing and code coverage.
#
# For use cases of `MathematicalSets` in external libraries, see the section
# `Usage` in the documentation.
# ===========================================================================

"""
    DummySet <: MathematicalSets.AbstractSet

The dummy set is a set with no fields. It is used to cover the `MathematicalSets`
functions.
"""
struct DummySet <: MathematicalSets.AbstractSet end

"""
    Singleton{T<:AbstractVector} <: AbstractSet

A singleton is a set with one element.

### Fields

- `element` -- element in the singleton
"""
struct Singleton{T<:AbstractVector} <: MathematicalSets.AbstractSet
    element::T
end

MathematicalSets.space_dimension(s::Singleton) = length(s.element)

# the dimension of a point is conventionally 0
MathematicalSets.dimension(s::Singleton) = 0

# the following metric properties for a singleton are usually zero, although
# depending on the applicaton one could be interested to have e.g. a Dirac δ
# if commented, defaults to returning `nothing` in `MathematicalSets`:
MathematicalSets.volume(s::Singleton) = 0.
MathematicalSets.surface_area(s::Singleton) = 0.

"""
    Ball{T} <: AbstractSet

Euclidean ``n``-dimensional ball in ``\\mathbb{R}^n``.

### Fields

- `center` -- vector that represents the center of the ball
- `radius` -- scalar that represents the radius of the ball
"""
struct Ball{T, VT <: AbstractVector{T}} <: MathematicalSets.AbstractSet
    center::VT
    radius::T
end
radius(S::Ball) = S.radius
center(S::Ball) = S.center

MathematicalSets.space_dimension(S::Ball) = length(S.center)
MathematicalSets.dimension(S::Ball) = space_dimension(S)

# convenience constructor: if only the radius and dimension are supplied,
# the ball is centered in the origin
Ball(R::T, n::Int) where T = Ball(zeros(T, n), R)

using SpecialFunctions # needed for the gamma function

function MathematicalSets.volume(S::Ball)
    n, R = dimension(S), radius(S)
    return π^(n/2) * R^n / gamma(n/2 + 1)
end

function MathematicalSets.surface_area(S::Ball)
    n, R = dimension(S), radius(S)
    return π^(n/2) * n * R^(n-1) / gamma(n/2 + 1)
end


text:
module TestArrays

include("preamble.jl")

@testset begin
    A = reshape(UnionVector(Any[1, 2.0]), 1, :)
    @test A isa UnionArray
    @test A[1] === 1
    @test A[2] === 2.0
    @test A[1, 1] === 1
    @test A[1, 2] === 2.0
    @test (A[1, 1] = 3.0) isa Any
    @test A[1, 1] === 3.0
    @test (A[1] = 3.0) isa Any
    @test A[1] === 3.0
end

end  # module

ext: jl


text:
#
# Unit tests for UnitfulAstrodynamics.jl
#

# include("R2BP.jl")
# include("CR3BP.jl")
# include("Propagation.jl")


text:
function Base.show(io::IO, design::PlackettBurman)
    print(io, "$(typeof(design))\n",
          "Dimension: $(size(design.matrix))\n",
          "Factors: $(design.factors)\n",
          "Dummy Factors: $(design.dummy_factors)\n",
          "Formula: $(design.formula)\n",
          "Design Matrix:\n")
    show(io, design.matrix)
end

function Base.show(io::IO, design::BoxBehnken)
    print(io, "$(typeof(design))\n",
          "Dimension: $(size(design.matrix))\n",
          "Factors: $(design.factors)\n",
          "Formula: $(design.formula)\n",
          "Design Matrix:\n")
    show(io, design.matrix)
end

function Base.show(io::IO, design::CentralComposite)
    print(io, "$(typeof(design))\n",
          "Dimension: $(size(design.matrix))\n",
          "Factors: $(design.factors)\n",
          "Formula: $(design.formula)\n",
          "Alpha: $(design.alpha)\n",
          "Face: $(design.face)\n",
          "Design Matrix:\n")
    show(io, design.matrix)
end

function Base.show(io::IO, design::FullFactorial)
    if typeof(design.matrix) == Missing
        print(io, "$(typeof(design))\n",
              "Factors: $(design.factors)\n",
              "Formula: $(design.formula)\n",
              "Design Matrix: On demand")
    else
        print(io, "$(typeof(design))\n",
              "Dimension: $(size(design.matrix))\n",
              "Factors: $(design.factors)\n",
              "Formula: $(design.formula)\n",
              "Design Matrix:\n")
        show(io, design.matrix)
    end
end

function Base.show(io::IO, design::FractionalFactorial2Level)
    if typeof(design.matrix) == Missing
        print(io, "$(typeof(design))\n",
              "Factors: $(design.factors)\n",
              "Formula: $(design.formula)\n",
              "Design Matrix: On demand")
    else
        print(io, "$(typeof(design))\n",
              "Dimension: $(size(design.matrix))\n",
              "Factors: $(design.factors)\n",
              "Formula: $(design.formula)\n",
              "Design Matrix:\n")
        show(io, design.matrix)
    end
end

function Base.show(io::IO, design::DesignDistribution)
    print(io, "$(typeof(design))\n",
          "Formula: $(design.formula)\n",
          "Factor Distributions:")

    for factor in keys(design.factors)
        print(io, "\n$(factor): $(design.factors[factor])")
    end
end

function Base.show(io::IO, design::RandomDesign)
    print(io, "$(typeof(design))\n",
          "Dimension: $(size(design.matrix))\n",
          "Factors: $(design.factors)\n",
          "Formula: $(design.formula)\n",
          "Design Matrix:\n")
    show(io, design.matrix)
end

function Base.show(io::IO, design::RandomLHCDesign)
    print(io, "$(typeof(design))\n",
          "Dimension: $(size(design.matrix))\n",
          "Factors: $(design.factors)\n",
          "Design Matrix:\n")
    show(io, design.matrix)
end

function Base.show(io::IO, design::OptimLHCDesign)
    print(io, "$(typeof(design))\n",
          "Dimension: $(size(design.matrix))\n",
          "Factors: $(design.factors)\n",
          "Fitness: $(design.fitness)\n",
          "Design Matrix:\n")
    show(io, design.matrix)
end

function Base.show(io::IO, design::OptimalDesign)
    print(io, "$(typeof(design))\n",
          "Dimension: $(size(design.matrix))\n",
          "Factors: $(design.factors)\n",
          "Formula: $(design.formula)\n",
          "Selected Candidate Rows: $(design.selected_experiments)\n",
          "Optimality Criteria: $(design.criteria)\n",
          "Design Matrix:\n")
    show(io, design.matrix)
end


text:
#----------------------------------------------------------------
#----------------------------------------------------------------
# Block Matrix multiplication

@inline function block_mat_mat_mul!(C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    @_sync begin
        @_spawn begin
            #_mul!(    C11, A11, B11, sz)
            gemm_kernel!(C11, A11, B11)
            _mul_add!(C11, A12, B21, sz)
        end
        @_spawn begin
            _mul!(    C12, A11, B12, sz)
            _mul_add!(C12, A12, B22, sz)
        end
        @_spawn begin
            _mul!(    C21, A21, B11, sz)
            _mul_add!(C21, A22, B21, sz)
        end
        _mul!(    C22, A21, B12, sz)
        _mul_add!(C22, A22, B22, sz)
    end
end

function block_mat_vec_mul!(C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:sz,     1:end];
        C21 = C[sz+1:end, 1:end];

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:end];
        B21 = B[sz+1:end, 1:end];
    end
    @_sync begin
        @_spawn begin
            #_mul!(    C11, A11, B11, sz)
            gemm_kernel!(C11, A11, B11)
            _mul_add!(C11, A12, B21, sz)
        end
        _mul!(    C21, A21, B11, sz)
        _mul_add!(C21, A22, B21, sz)
    end
end

function block_covec_mat_mul!(C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:end,    1:sz]; C12 = C[1:end,    sz+1:end]

        A11 = A[1:end,    1:sz]; A12 = A[1:end,    sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    @_sync begin
        @_spawn begin
            #_mul!(    C11, A11, B11, sz)
            gemm_kernel!(C11, A11, B11)
            _mul_add!(C11, A12, B21, sz)
        end
        _mul!(    C12, A11, B12, sz)
        _mul_add!(C12, A12, B22, sz)
    end
end

function block_vec_covec_mul!(C, A, B, sz)
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:end];
        A21 = A[sz+1:end, 1:end];

        B11 = B[1:end,     1:sz]; B12 = B[1:end,     sz+1:end]
    end
    @_sync begin
        @_spawn begin
            #_mul!(    C11, A11, B11, sz)
            gemm_kernel!(C11, A11, B11)
        end
        @_spawn begin
            _mul!(C12, A11, B12, sz)
        end
        @_spawn begin
            _mul!(C21, A21, B11, sz)
        end
        _mul!(C22, A21, B12, sz)
    end
end

function block_covec_vec_mul!(C, A, B, sz)
    @inbounds @views begin
        A11 = A[1:end,    1:sz]; A12 = A[1:sz,     sz+1:end]

        B11 = B[1:sz,     1:end];
        B21 = B[sz+1:end, 1:end];
    end
    gemm_kernel!(C, A11, B11)
    #_mul!(    C, A11, B11, sz)
    _mul_add!(C, A12, B21, sz)
end

#----------------------------------------------------------------
#----------------------------------------------------------------
# Block Matrix addition-multiplication

@inline function block_mat_mat_mul_add!(C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    @_sync begin
        @_spawn begin
            add_gemm_kernel!(C11, A11, B11, Val(factor))
            _mul_add!(C11, A12, B21, sz, Val(factor))
        end
        @_spawn begin
            _mul_add!(C12, A11, B12, sz, Val(factor))
            _mul_add!(C12, A12, B22, sz, Val(factor))
        end
        @_spawn begin
            _mul_add!(C21, A21, B11, sz, Val(factor))
            _mul_add!(C21, A22, B21, sz, Val(factor))
        end
        _mul_add!(C22, A21, B12, sz, Val(factor))
        _mul_add!(C22, A22, B22, sz, Val(factor))
    end
end

function block_mat_vec_mul_add!(C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz,     1:end];
        C21 = C[sz+1:end, 1:end];

        A11 = A[1:sz,     1:sz]; A12 = A[1:sz,     sz+1:end]
        A21 = A[sz+1:end, 1:sz]; A22 = A[sz+1:end, sz+1:end]

        B11 = B[1:sz,     1:end];
        B21 = B[sz+1:end, 1:end];
    end
    @_sync begin
        @_spawn begin
            add_gemm_kernel!(C11, A11, B11, Val(factor))
            _mul_add!(C11, A12, B21, sz, Val(factor))
        end
        _mul_add!(C21, A21, B11, sz, Val(factor))
        _mul_add!(C21, A22, B21, sz, Val(factor))
    end
end

function block_covec_mat_mul_add!(C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:end,    1:sz]; C12 = C[1:end,     sz+1:end]

        A11 = A[1:end,    1:sz]; A12 = A[1:end,    sz+1:end]

        B11 = B[1:sz,     1:sz]; B12 = B[1:sz,     sz+1:end]
        B21 = B[sz+1:end, 1:sz]; B22 = B[sz+1:end, sz+1:end]
    end
    @_sync begin
        @_spawn begin
            add_gemm_kernel!(C11, A11, B11, Val(factor))
            _mul_add!(C11, A12, B21, sz, Val(factor))
        end
        _mul_add!(C12, A11, B12, sz, Val(factor))
        _mul_add!(C12, A12, B22, sz, Val(factor))
    end
end

function block_vec_covec_mul_add!(C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        C11 = C[1:sz,     1:sz]; C12 = C[1:sz,     sz+1:end]
        C21 = C[sz+1:end, 1:sz]; C22 = C[sz+1:end, sz+1:end]

        A11 = A[1:sz,     1:end];
        A21 = A[sz+1:end, 1:end];

        B11 = B[1:end,    1:sz]; B12 = B[1:end,    sz+1:end]
    end
    @_sync begin
        @_spawn begin
            add_gemm_kernel!(C11, A11, B11, Val(factor))
        end
        @_spawn begin
            _mul_add!(C12, A11, B12, sz, Val(factor))
        end
        @_spawn begin
            _mul_add!(C21, A21, B11, sz, Val(factor))
        end
        _mul_add!(C22, A21, B12, sz, Val(factor))
    end
end

function block_covec_vec_mul_add!(C, A, B, sz, ::Val{factor} = Val(1)) where {factor}
    @inbounds @views begin
        A11 = A[1:end,    1:sz]; A12 = A[1:end,    sz+1:end]

        B11 = B[1:sz,     1:end];
        B21 = B[sz+1:end, 1:end];
    end
    add_gemm_kernel!(C, A11, B11, Val(factor))
    _mul_add!(C, A12, B21, sz, Val(factor))
end


text:
semantic_token_test() = LanguageServer.textDocument_semanticTokens_full_request(LanguageServer.SemanticTokensParams(LanguageServer.TextDocumentIdentifier("testdoc"), missing, missing), server, server.jr_endpoint)

@testset "function calls" begin
        settestdoc("""
        function hello()
            println("hello world")
        end
        """)
        @test semantic_token_test() == SemanticToken(0, 9, 5, SemanticTokenKinds.Function)
end


text:
#!/usr/bin/env julia

using Gtk
using Cairo
using GtkBuilderAid

# Test using canvas features of DrawingArea
mutable struct CanvasData
  surface::CairoSurface
end

# Clear the surface to white
function clear_surface(surf::CairoSurface)
  ctx = CairoContext(surf)
  set_source_rgb(ctx, 1, 1, 1)
  paint(ctx)
end

function draw_brush(
    canvas,
    x,
    y,
    surf)
  px = Int(round(x - 3))
  py = Int(round(y - 3))
  ctx = CairoContext(surf)
  rectangle(ctx, px, py, 6, 6)
  fill(ctx)
  reveal_area(canvas, px, py, 6, 6)
end

ft = @GtkFunctionTable begin

@guarded Cint(0) function configure_event_cb(
    canvas,
    configure_event,
    qsdata)
  userdata = qsdata.userdata
  userdata.surface = create_similar_surface(canvas, Gtk.GEnum(Cairo.CONTENT_COLOR_ALPHA))
  clear_surface(userdata.surface)
  return Cint(1)
end

@guarded Cint(1) function draw_cb(
    canvas,
    ctx_ptr,
    qsdata)
  userdata = qsdata.userdata

  ctx = CairoContext(ctx_ptr)

  set_source_surface(ctx, userdata.surface, 0, 0)
  paint(ctx)
  return Cint(0)
end

@guarded Cint(0) function button_press_event_cb(
    canvas,
    event_ptr,
    qsdata)
  userdata = qsdata.userdata

  event = Gtk.GdkEvent(event_ptr)

  # Make sure the surface is valid
  if userdata.surface.ptr == C_NULL
    return Cint(0)
  end

  if event.button == 1
    draw_brush(canvas, event.x, event.y, userdata.surface)
  elseif event.button == 3
    clear_surface(userdata.surface)
    reveal(GObject(canvas)::Gtk.GtkWidget, false)
  end

  return Cint(1)
end

@guarded Cint(0) function motion_notify_event_cb(
    canvas,
    event_ptr,
    qsdata)
  userdata = qsdata.userdata

  event = Gtk.GdkEvent(event_ptr)

  # Make sure the surface is valid
  if userdata.surface.ptr == C_NULL
    return Cint(0)
  end

  if (event.state & Gtk.GdkModifierType.BUTTON1) != 0
    draw_brush(canvas, event.x, event.y, userdata.surface)
  end

  return Cint(1)

end

end

start_application(ft, "io.github.matt5sean3.GtkBuilderAid.second", "drawing_window", "resources/main.ui", CanvasData(CairoSurface(C_NULL, -1, -1)))


text:
#=
/* Switches for the triangulator.                                            */

-p Triangulates a Planar Straight Line Graph (.poly file).
-r Refines a previously generated mesh.
-q Quality mesh generation with no angles smaller than 20 degrees. An alternate minimum angle may be specified after the `q'.
-a Imposes a maximum triangle area constraint. A fixed area constraint (that applies to every triangle) may be specified after the `a', or varying area constraints may be read from a .poly file or .area file.
-u Imposes a user-defined constraint on triangle size.
-A Assigns a regional attribute to each triangle that identifies what segment-bounded region it belongs to.
-c Encloses the convex hull with segments.
-D Conforming Delaunay: use this switch if you want all triangles in the mesh to be Delaunay, and not just constrained Delaunay; or if you want to ensure that all Voronoi vertices lie within the triangulation.
-j Jettisons vertices that are not part of the final triangulation from the output .node file (including duplicate input vertices and vertices ``eaten'' by holes).
-e Outputs (to an .edge file) a list of edges of the triangulation.
-v Outputs the Voronoi diagram associated with the triangulation. Does not attempt to detect degeneracies, so some Voronoi vertices may be duplicated.
-n Outputs (to a .neigh file) a list of triangles neighboring each triangle.
-g Outputs the mesh to an Object File Format (.off) file, suitable for viewing with the Geometry Center's Geomview package.
-B Suppresses boundary markers in the output .node, .poly, and .edge output files.
-P Suppresses the output .poly file. Saves disk space, but you lose the ability to maintain constraining segments on later refinements of the mesh.
-N Suppresses the output .node file.
-E Suppresses the output .ele file.
-I Suppresses mesh iteration numbers.
-O Suppresses holes: ignores the holes in the .poly file.
-X Suppresses exact arithmetic.
-z Numbers all items starting from zero (rather than one). Note that this switch is normally overrided by the value used to number the first vertex of the input .node or .poly file. However, this switch is useful when calling Triangle from another program.
-o2 Generates second-order subparametric elements with six nodes each.
-Y Prohibits the insertion of Steiner points on the mesh boundary. If specified twice (-YY), it prohibits the insertion of Steiner points on any segment, including internal segments.
-S Specifies the maximum number of added Steiner points.
-i Uses the incremental algorithm for Delaunay triangulation, rather than the divide-and-conquer algorithm.
-F Uses Steven Fortune's sweepline algorithm for Delaunay triangulation, rather than the divide-and-conquer algorithm.
-l Uses only vertical cuts in the divide-and-conquer algorithm. By default, Triangle uses alternating vertical and horizontal cuts, which usually improve the speed except with vertex sets that are small or short and wide. This switch is primarily of theoretical interest.
-s Specifies that segments should be forced into the triangulation by recursively splitting them at their midpoints, rather than by generating a constrained Delaunay triangulation. Segment splitting is true to Ruppert's original algorithm, but can create needlessly small triangles. This switch is primarily of theoretical interest.
-C Check the consistency of the final mesh. Uses exact arithmetic for checking, even if the -X switch is used. Useful if you suspect Triangle is buggy.
-Q Quiet: Suppresses all explanation of what Triangle is doing, unless an error occurs.
-V Verbose: Gives detailed information about what Triangle is doing. Add more `V's for increasing amount of detail. `-V' gives information on algorithmic progress and detailed statistics.
-h Help: Displays complete instructions.

/*  If the size of the object file is important to you, you may wish to      */
/*  generate a reduced version of triangle.o.  The REDUCED symbol gets rid   */
/*  of all features that are primarily of research interest.  Specifically,  */
/*  the -DREDUCED switch eliminates Triangle's -i, -F, -s, and -C switches.  */
/*  The CDT_ONLY symbol gets rid of all meshing algorithms above and beyond  */
/*  constrained Delaunay triangulation.  Specifically, the -DCDT_ONLY switch */
/*  eliminates Triangle's -r, -q, -a, -u, -D, -Y, -S, and -s switches.       */
=#

mutable struct TriangulateOptions
    pslg::Bool #p
    regionattrib::Bool # A
    convex::Bool # c
    jettison::Bool # j
    firstnumberiszero::Bool # z
    edgesout::Bool # e
    voronoi::Bool # v
    neighbors::Bool # n
    nobound::Bool # B
    nopolywritten::Bool # P
    nonodewritten::Bool # N
    noelewritten::Bool # E
    noiterationnum::Bool # I
    noholes::Bool # O
    noexactaritmetic::Bool # X
    order::Bool # o
    orderHow::Int64 # 1...2...3
    dwyer::Bool # l
    quiet::Bool # Q
    verbose::Bool # V
    TriangulateOptions() = new(false,
    false, false, false, false, false, false, false,
    # No_xyz_ selector(s)
    true, true, false, false, true, false, false,
    # order
    false, 0,
    false,
    # Quiet Verbose
    true, false)
end

function getTriangulateStringOptions(self::TriangulateOptions)
    output_stri = ""

    if self.pslg
        output_stri = output_stri * "p"
    end

    if self.regionattrib
        output_stri = output_stri * "A"
    end

    if self.convex
        output_stri = output_stri * "c"
    end

    if self.jettison
        output_stri = output_stri * "j"
    end

    if self.firstnumberiszero
        output_stri = output_stri * "z"
    end

    if self.edgesout
        output_stri = output_stri * "e"
    end

    if self.voronoi
        output_stri = output_stri * "v"
    end

    if self.neighbors
        output_stri = output_stri * "n"
    end

    if self.nobound
        output_stri = output_stri * "B"
    end

    if self.nopolywritten
        output_stri = output_stri * "P"
    end

    if self.nonodewritten
        output_stri = output_stri * "N"
    end

    if self.noelewritten
        output_stri = output_stri * "E"
    end

    if self.noiterationnum
        output_stri = output_stri * "I"
    end

    if self.noholes
        output_stri = output_stri * "O"
    end

    if self.noexactaritmetic
        output_stri = output_stri * "X"
    end

    if self.order && self.orderHow > 0
        output_stri = output_stri * "o" * string(self.orderHow)
    end

    if self.dwyer
        output_stri = output_stri * "l"
    end

    if self.quiet
        output_stri = output_stri * "Q"
    end

    if self.verbose
        output_stri = output_stri * "V"
    end

    return output_stri
end


text:
module TestConditional
using Test
using LinearAlgebra
using GaussianDistributions
using GaussianDistributions: ConditionalMean
using ForwardDiff: derivative

const nsigma = 4

@testset "conditional" begin
    @testset "multivariate" begin
        d = 3
        μ = randn(d)
        Σ = randn(d, d)
        Σ = Σ'Σ
        G = Gaussian(μ, Σ)

        x = randn()
        tol = 1e-13

        C = G | (1, x)
        @test C == conditional(G, 1, x)
        @test mean(C[1]) ≈ x
        @test cov(C[1]) < tol # collapse of uncertainty
        @test all(vc < vg for (vc, vg) in zip(var(C), var(G))) # uniform reduction in uncertainty

        x = randn(2)
        C = G | (1:2, x)
        @test mean(C[1:2]) ≈ x
        @test sum(std(C[1:2])) < tol # collapse of uncertainty

        # play around with posterior
        x = randn()
        σ = .1
        C = G | (1, x, σ^2)
        @test abs(mean(C[1])-x) < abs(mean(G[1])-x) # posterior mean is closer to data than prior mean
        @test cov(C[1]) < cov(G[1])  # posterior uncertainty is smaller than prior
        @test sum(var(C)) < sum(var(G))

        # test conversion of vector of noise variances to diagonal covariance
        n = 2
        x = randn(n)
        σ² = σ^2 * exp.(randn(n))
        C = G | (1:n, x, σ²)
        E = G | (1:n, x, Diagonal(σ²))
        @test C isa Gaussian
        @test mean(C) ≈ mean(E)
        @test cov(C) ≈ cov(E)
    end

    @testset "bayesian linear regression" begin
        # i.e. conditioning on linear operator data
        n, m = 4, 8
        A = randn(n, m)
        σx = 2
        x = σx * randn(m)
        b = A*x

        x = Gaussian(σx^2*I(m))

        # noiseless conditioning
        P = x | (A, b)
        PC = Matrix(cov(P))
        λ = sort!(eigvals(PC))
        @test all(<(1e-12), λ[1:n])
        @test all(>(1e-12), λ[n+1:end])
        @test A*mean(P) ≈ b

        # test iterative conditioning
        Q = x | (A[1:n-1, :], b[1:n-1])
        Q = Q | (A[n, :], b[n])
        @test mean(Q) ≈ mean(P)
        @test cov(Q) ≈ cov(P)

        # with noise
        σ = 1e-6
        y = Gaussian(b, σ^2*I(n))
        P = x | (A, y)

        # n observations only uniquely determine n components
        PC = Hermitian(Matrix(cov(P)))
        λ = sort!(eigvals(PC))
        @test all(<(σ), λ[1:n])
        @test all(>(σ), λ[n+1:end])

        # if we observe a complete set, we can reconstruct the solution
        @test norm(A*mean(P)-mean(y)) < 3σ

        # test iterative conditioning
        Q = x | (A[1:n-1, :], y[1:n-1])
        Q = Q | (A[n, :], y[n])
        @test mean(Q) ≈ mean(P)
        @test cov(Q) ≈ cov(P)
    end

    @testset "process mean" begin
        # create dataset
        G = 4
        x = sort(randn(G))
        f(x) = sin(2π*x)

        l = .1
        k(x, y) = exp(-(x-y)^2/2l^2)
        G = Gaussian(k)
        K = k.(x, x')
        cm = ConditionalMean(G, x, f.(x), K)

        @test marginal(cm, x) ≈ f.(x)

        # test differentibility
        using ForwardDiff: derivative
        g(x) = derivative(cm, x)
    end

    @testset "process" begin
        # create dataset
        n = 16
        x = sort(randn(n))
        f(x) = sin(2π*x)
        y = f.(x)

        # setup Gaussian process
        l = .1
        k(x, y) = exp(-sum(abs2, x-y)/2l^2)
        G = Gaussian(k)
        K = k.(x, x')
        C = G | (x, y) # noiseless conditioning
        @test mean(C, x) ≈ f.(x)
        tol = 10eps(Float64)

        v = var(C, x)
        @test all(v .< tol) # collapse of uncertainty

        # noisy conditioning
        n = 128
        x = sort(randn(n))
        y = f.(x)
        σ = .01
        @. y += σ*randn() # adding noise
        C = G | (x, y, σ^2)
        d = mean(C, x)-f.(x)
        @test all(abs.(d) .< nsigma*σ)
        v = var(C, x)
        @test all(v .< σ^2) # collapse of uncertainty

        # test differentibility
        g(x) = derivative(mean(C), x)
        @test isapprox(g(0), 2π, rtol = .1)

        # TODO: iterative conditioning
        # xs = randn()
        # ys = f(xs)
        # C = C | (xs, ys, σ^2)

        # d = mean(C)(xs) - f(xs)
        # @test abs(d) < nsigma*σ
        # tol = σ
        # v = std(C)(xs)
        # @test v < tol

        @testset "single datapoint conditioning" begin
            x = randn()
            y = f(x)
            C = G | (x, y)
            @test mean(C)(x) ≈ y
            @test std(C)(x) < 1e-12

            # with noise
            y += σ*randn()
            C = G | (x, y, σ^2)
            @test isapprox(mean(C)(x), y, atol = nsigma*σ)
            @test std(C)(x) < nsigma*σ
            @test std(C)(x) > σ/2 # uncertainty does not collapse completely

            # single point conditioning with 2d input
            x = randn(2)
            y = randn()
            C = G | (x, y)
            @test mean(C)(x) ≈ y
            @test std(C)(x) < 1e-12

            # with noise
            y += σ*randn()
            C = G | (x, y, σ^2)
            @test isapprox(mean(C)(x), y, atol = nsigma*σ)
            @test std(C)(x) < nsigma*σ
        end
    end
end

end


text:
using Replay

instructions = ["using UnicodeFun", "\"\\\\pi\" |> to_latex"]

replay(instructions, julia_project=@__DIR__)


text:
EIN = (
    ("asmUnit.Psi1", [0.0, -1.0] ),
    ("asmUnit.Psi2", [0.3, -0.9] ),
    ("transformer.X.ix", [0.0, 0.0] ),
    ("J.w", 1.0),
    (:LOAD, 0.5),

)

AUS = (
    ("der(asmUnit.Psi1)", [0.0, 0.0]),
    ("der(asmUnit.Psi2)", [0.0, 0.0]),
#     ("der(transformer.X.ix)", [0.0, 0.0]),  vector init coming soon
    ("der(J.w)", 0.0),
    ("transformer.X.derINI1", 0.0),
    ("transformer.X.derINI2", 0.0),
    ("meterGRID.powerRe", -0.8),

)

max_stars_repo_name: johhell/EPSphasor

max_issues_repo_name: johhell/EPSphasor

max_forks_repo_name: johhell/EPSphasor


text:
t1 = read_tree("$(dirname(pathof(TreeKnit)))/..//test/splitgraph/3solutions/t1.nwk")
t2 = read_tree("$(dirname(pathof(TreeKnit)))/..//test/splitgraph/3solutions/t2.nwk")

treelist = deepcopy(Any[t1, t2])

mcc = naive_mccs(treelist)
mcc_names = TreeKnit.name_mcc_clades!(treelist, mcc)
for (i,t) in enumerate(treelist)
	treelist[i] = TreeKnit.reduce_to_mcc(t, mcc)
end

g = trees2graph(treelist)

oa = OptArgs(resolve = false, likelihood_sort = false, γ = 2)
Trange = reverse(oa.Tmin:oa.dT:oa.Tmax)

@testset "No likelihood" begin
	local out = SplitGraph.sa_opt(g; Trange=Trange, γ=2, M = 100, resolve = false)[1]
	sort!(out)
	@test out[1] == Bool[0,1,1]
	@test out[2] == Bool[1,0,1]
	@test out[3] == Bool[1,1,0]
end

@testset "With likelihood" begin
	local out = []
	for rep in 1:20
		tmp = SplitGraph.opttrees(
			t1, t2;
			γ=2, M=10, Trange=Trange, likelihood_sort = true
		)
		if !in(tmp[1], out)
			push!(out, tmp[1])
		end
	end
	@test out[1] == [["C1","C2"]]
end

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/RecombTools

max_forks_repo_name: PierreBarrat/RecombTools


text:
struct RNASeq <: GeneticSeq
  data::BitArray{2}

  function RNASeq(x; checkinput::Bool=true)
    if checkinput
      if _checkinput(x)
        return new(x)
      else
        throw(ErrorException("Invalid input"))
      end
    else
      return new(x)
    end
  end
end

function onehot(::Type{RNASeq}, x::Char)
  if x     == 'A'
    return BitArray{1}([1, 0, 0, 0])
  elseif x == 'C'
    return BitArray{1}([0, 1, 0, 0])
  elseif x == 'G'
    return BitArray{1}([0, 0, 1, 0])
  elseif x == 'U'
    return BitArray{1}([0, 0, 0, 1])
  else
    throw(ErrorException("Invalid `Char`"))
  end
end

function onehotinv(::Type{RNASeq}, x::BitArray{1})
  return ['A', 'C', 'G', 'U'][x][1]
end

function RNASeq(seq::Vector{Char})
  x = BitArray{2}(fill(0, (4, length(seq))))
  for i in 1:length(seq)
    if seq[i]     == 'A'
      x[1, i] = 1
    elseif seq[i] == 'C'
      x[2, i] = 1
    elseif seq[i] == 'G'
      x[3, i] = 1
    elseif seq[i] == 'U'
      x[4, i] = 1
    else
      throw(ErrorException("Invalid `Char` at index $i"))
    end
  end
  return RNASeq(x, checkinput=false)
end

function RNASeq(x::String)
  return RNASeq(Vector{Char}(x))
end

function RNASeq(x::BitArray{1})
  return RNASeq(reshape(x, (4, 1)))
end


text:
"""
This is Tutorial_ETOPO1.jl. It contains all the necessary commands to import the netCDF file from ETOPO1,
to convert it to the GeoData format and to export it to a Paraview format. The following steps are performed:
1. Read data from this file.
2. Select a longitude and latitude range.
3. Put the data in a GeoData format (this is the format that is used internally in the GMG).
4. Export that data to a format readable by Paraview.

ETOPO1 data files used in this example can be downloaded here:
[https://ngdc.noaa.gov/mgg/global/global.html](https://ngdc.noaa.gov/mgg/global/global.html)

doi:10.7289/V5C8276M
"""

using GMT, GeoParams, GeophysicalModelGenerator, Shapefile, Luxor

# 1. define where the file is located on your computer
filename_topo = "./ETOPO1/ETOPO1_Ice_g_gmt4.grd" # this only works if you are in the directory where the data is located
filename_geo  = "./tectonic_maps_4dmb_2020_09_17/shape_files/tect_units_alcapadi.shp"
filename_faults = "./tectonic_maps_4dmb_2020_09_17/shape_files/faults_alcapadi.shp"
filename_gmt = "./tectonic_maps_4dmb_2020_09_17/GMT_example/alcapadi_polygons.gmt"
plot("alcapadi_polygons.gmt", show=true)
# 3. Select latitude and longitude range
lat_min = 37.0
lat_max = 49.0
lon_min = 4.0
lon_max = 20.0
region = [lon_min lon_max lat_min lat_max]

# 2. load desired topography data using GMT.jl and adapt the data to make it compatible with GeoData
G = gmtread(filename_topo, limits=[lon_min,lon_max,lat_min,lat_max]);
Lon,Lat,Depth    =   LonLatDepthGrid(G.x[1:end],G.y[1:end],0);
Depth[:,:,1]     =   1e-3*G.z';

# 3. load data for the geological map
TectUnitData = zeros(size(Lon));

# assign "tectonic units" as phases to surface grid using shapefile import
table = Shapefile.Table(filename_geo)
tectunits = unique(table.tect_unit); # get tectonic units

    for ishape = 1:length(table)
        phase = findall(==(0),cmp.(table.tect_unit[ishape], tectunits)); # assign a certain phase to the respective tectonic unit

        #xv = zeros(length(table.geometry[ishape].points),1); # preallocate the vector to hold the different points
        #yv = zeros(length(table.geometry[ishape].points),1); # preallocate the vector to hold the different points
        xv = Vector{Float64}(undef,length(table.geometry[ishape].points));
        yv = Vector{Float64}(undef,length(table.geometry[ishape].points));

        for ipoint = 1:length(table.geometry[ishape].points)
            xv[ipoint] = table.geometry[ishape].points[ipoint].x;
            yv[ipoint] = table.geometry[ishape].points[ipoint].y;
        end

        xa = vec(Lon);
        ya = vec(Lat);

        polygon = Point.(xv,yv);
        points  =  Point.(xa,ya);

        inside = [isinside(p, polygon; allowonedge=true) for p in points]

        inside = inside .& (vec(TectUnitData).==0);

        TectUnitData[inside] .= phase;
    end

# Set up the Data structure
Data_set        =   GeoData(Lon, Lat, Depth, (Topography=Depth*km,tect_unit=TectUnitData))

# Export the data structure to Paraview format
Write_Paraview(Data_set, "test_GeoMap")

# Set up the Data structure
data_Topo        =   GeoData(Lon, Lat, Depth, (Topography=Depth*km,))

# Export the data structure to Paraview format
Write_Paraview(Data_set, "test_netcdf_ETOPO1")


text:
@time include("0_setup.jl")
using IterableTables, FastGroupBy

@everywhere import JuliaDB.extractarray

if !isdefined(:data)
    @time data = load("lc_juliadb")
end

# fields = fieldnames(eltype(data))
fields = names(data)
ftypes = eltype(data).parameters

xx = JuliaDB.extractarray(data, x -> x.loan_status)

numeric_cols = [fields[i] for i = 1:length(fields) if (ftypes[i] <: Number ||
            ftypes[i] <: Nullable{<:Number}) && !(fields[i] in [:id, :member_id, :dti_joint])];

const bad_statuses = ("Late (16-30 days)","Late (31-120 days)","Default","Charged Off")
good_loans = filter(data) do row
    !(row.loan_status in bad_statuses)
end
bad_loans = filter(x->x.loan_status in bad_statuses, data)

using Gadfly
import NullableArrays: dropnull

# Density plot for bad and good loans
plots = Gadfly.Plot[]

good_numbers = good_loans[:loan_amnt]
@which columns(good_loans,:loan_amnt)

for (name, g, b) in zip(numeric_cols, columns(good_loans), columns(bad_loans))
    g = NullableArrays.dropnull(g)
    b = NullableArrays.dropnull(b)
    p = plot(layer(x=g, Geom.density, Theme(default_color=colorant"green")),
             layer(x=b, Geom.density, Theme(default_color=colorant"red")),
              Guide.title(string(name)), Guide.ylabel("density"))
    push!(plots, p)
end

perm = randperm(length(data))
train_till = round(Int, length(data) * 3/4)

training_subidx = sort!(perm[1:train_till])
testing_subidx = sort!(perm[train_till+1:end]);

training_subset = data[training_subidx]
testing_subset = data[testing_subidx]

features_train = [revol_util_train int_rate_train mths_since_last_record_train annual_inc_joint_train total_rec_prncp_train all_util_train]

labels = collect(map(x->x in bad_statuses, data))

import NullableArrays: dropnull
findnonnulls(xs::Columns) = find(x->!any(map(isnull, x)), xs)

function input_matrix(table, fields)
    float_features = collect(values(table, fields))
    tmp = findnonnulls(float_features) # indices of the rows where all fields are non-null
    nzidxs = collect(keys(table,1)[tmp]) # corresponding indices in the table
    reduce(hcat, map(dropnull, columns(float_features[tmp]))), values(labels[nzidxs]) # matrix, label vector
end

training_matrix, train_labels = input_matrix(training_subset, numeric_features)
model = build_forest(train_labels, training_matrix, 3, 10, 0.8, 6)

#     0.512169 seconds (7.67 k allocations: 440.766 KiB)
#   Out[98]:
#   Ensemble of Decision Trees
#   Trees:      10
#   Avg Leaves: 41.8
#   Avg Depth:  8.0

f = open("  loanmodel.jls", "w")
serialize(f, model)
close(f)

features_test, test_labels = input_matrix(testing_subset, numeric_features)

@time predictions = mapslices(features_test, 2) do fs
    DecisionTree.apply_forest(model, fs)
end;

# Receiver Operating Characteristics curve
using ROC
curve = roc(convert(Vector{Float64}, predictions[:]),convert(BitArray{1}, test_labels))

# An ROC plot in Gadfly with data calculuated using ROC.jl

plot(layer(x=curve.FPR, y=curve.TPR, Geom.line),
     layer(x = linspace(0.0,1.0,101), y = linspace(0.0,1.0,101),
       Geom.line, Theme(default_color=colorant"red")), Guide.title("ROC"),
       Guide.xlabel("False Positive Rate"),Guide.ylabel("True Positive Rate"))

# Area Under Curve
AUC(curve)

# => 0.7291976225854385

max_stars_repo_path: src/julia/lending club/3_load_data.jl

max_issues_repo_path: src/julia/lending club/3_load_data.jl

max_forks_repo_path: src/julia/lending club/3_load_data.jl


text:
mutable struct Network <: NeuralNet

    function Network()
    """
    TODO:    prompts user for all input

    #read data
    #save data
    #show data (tabulated form?)
    """

    end

function Network(id)
    """
        TODO:    prompts user for input

        # Satisfied with only being given an id
            Used for Debugging
        """

end
    #recognize, track, memorize, and perform many other useful functions

    function recognize(id)
        """ gets model from an id
        compares a 'model' with already existing info (in a JSON database?)
        """
    end

    function find(id)
        """
        retrieve a neural Network object, from a given id
        """
    end
    function memorize(id)
        """
        adds to the memory database (JSON file) - connfiured for NoSQL Data Entry
        """
    end

    """
    further NoSQL commands, here, below
    """

end

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text:
push!(LOAD_PATH,"../src/")
using Documenter, ImageContrastAdjustment, ColorVectorSpace, ColorTypes
makedocs(sitename="Documentation",
            format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"))
deploydocs(repo = "github.com/zygmuntszpak/ImageContrastAdjustment.jl.git")


text:
import Base.:+
import Base.:-
import Base.:*

export ℝ³
export cross

"""
    Represents a point in ℝ³.

field: a.
"""
struct ℝ³ <: VectorSpace
    a::Array{Float64}
    ℝ³(x₁::Real, x₂::Real, x₃::Real) = new(float.([x₁; x₂; x₃]))
    ℝ³(a::Array) = begin
        @assert(length(a) == 3, "The input vector must contain exactly three elements.")
        ℝ³(a...)
    end
end

## Unary Operators ##

+(r::ℝ³) = r
-(r::ℝ³) = ℝ³(-vec(r))

## Binary Operators ##

+(r1::ℝ³, r2::ℝ³) = ℝ³(vec(r1) + vec(r2))
-(r1::ℝ³, r2::ℝ³) = ℝ³(vec(r1) - vec(r2))
*(r::ℝ³, λ::Real) = ℝ³(λ .* vec(r))
*(λ::Real, r::ℝ³) = r * λ

"""
    cross(r1, r2)

Perform a cross product with the given vectors `r1` and `r2`.
"""
cross(r1::ℝ³, r2::ℝ³) = begin
    M = transpose(reshape([vec(r1); vec(r1); vec(r2)], :, length(r1)))
    M = convert(Array{Float64}, M)
    ℝ³(map(x -> cofactor(M, 1, x), 1:length(r1)))
end

Base.isapprox(a::Array{ℝ³},
              b::Array{ℝ³};
              atol::Float64 = TOLERANCE) = begin
    for (elementa, elementb) in zip(a, b)
        if isapprox(elementa, elementb, atol = atol) == false
            return false
        end
    end
    return true
end

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
# FUNCTIONS

function sum_of_squares(vec)
    return sum(n^2 for n in vec)
end

function square_of_sum(vec)
    return sum(vec)^2
end

# SCRIPT

function main()
    vec=1:100
    println(square_of_sum(vec) - sum_of_squares(vec))
end

main()

max_stars_repo_name: shanedrabing/project-euler

max_issues_repo_name: shanedrabing/project-euler

max_forks_repo_name: shanedrabing/project-euler


text:
"""
    HDF5Logger

Creates an object for logging frames of data to an HDF5 file. The frames are
expected to have the same size, and the total number of frames to log is
expected to be known in advance. It currently works with scalars, vectors, and
matrices of basic types that the HDF5 package supports.

# Example 1: Logging a Vector over Time

```julia
log = Log("my_file.h5") # Create a logger.
num_frames = 100        # Set the total number of frames to log.
example_data = [1., 2., 3.] # Create an example of a frame of data.
add!(log, "/a_vector", example_data, num_frames) # Create the stream.
log!(log, "/a_vector", [4., 5., 6.]) # Log a single frame of data.
log!(log, "/a_vector", [7., 8., 9.]) # Log the next frame.
close!(log) # Always clean up when done. Use try-catch to make sure.

# Read it back out using the regular HDF5 library.
using HDF5
x = HDF5.h5open("my_file.h5", "r") do logs
    read(logs, "/a_vector")
end
```

# Example 2: Logging a Scalar over Time

The underlying HDF5 library (HDF5.jl) logs things as arrays, so passing in
scalars will result in 1-by-n arrays in the HDF5 file, instead of an n-length
Vector. One can retrieve the n-length Vector by using `squeeze`. E.g.:

```julia
log = Log("my_file.h5") # Create a logger.
add!(log, "/a_scalar", 0., 1000) # Create stream with space for 1000 samples.
log!(log, "/a_scalar", 0.1) # Log a single frame of data.
log!(log, "/a_scalar", 0.2) # Log the next frame.
close!(log) # Always clean up when done. Use try-catch to make sure.

# Read it back out using the regular HDF5 library.
using HDF5
t = HDF5.h5open("my_file.h5", "r") do logs
    read(logs, "/a_scalar")
end # t is now a 1-by-1000 Array.
t = squeeze(t, 1) # Make an n-element Vector.
```

# Notes

One often needs to `log!` immediately after an `add!`, so the `add!` function
can also log its example data as the first frame for convenience. Just use:

```julia
add!(log, "/group/name", data, num_frames, true) # Add stream; log data.
```

"""
module HDF5Logger

using HDF5

export Log, add!, log!, close!

mutable struct Stream
    count::Int64
    length::Int64
    rank::Int64
    dataset::HDF5Dataset
end

struct Log
    streams::Dict{String,Stream}
    file_id::HDF5File
    function Log(file_name::String)
        new(Dict{String,Stream}(), h5open(file_name, "w"))
    end
end

function prepare_group!(log::Log, slug::String)
    groups = filter(x->!isempty(x), split(slug, '/')) # Explode to group names
    group_id = g_open(log.file_id, "/") # Start at top group
    for k = 1:length(groups)-1 # For each group up to dataset
        if exists(group_id, String(groups[k]))
            group_id = g_open(group_id, String(groups[k]))
        else
            group_id = g_create(group_id, String(groups[k]))
        end
    end
    return (group_id, String(groups[end])) # Group ID and name of dataset
end

function add!(log::Log, slug::String, data, num_samples::Int64, keep::Bool = false)
    dims                    = isbits(data) ? 1 : size(data)
    group_id, group_name    = prepare_group!(log, slug)
    dataset_id              = d_create(group_id, group_name,
                                       datatype(eltype(data)),
                                       dataspace(dims..., num_samples))
    log.streams[slug]       = Stream(0, num_samples, length(dims), dataset_id)
    if keep
        log!(log, slug, data)
    end
end

function log!(log::Log, slug::String, data)
    @assert(haskey(log.streams, slug),
            "The logger doesn't have that key. Perhaps you need to `add` it?")
    @assert(log.streams[slug].count < log.streams[slug].length,
            "We've already used up all of the allocated space for this stream!")
    log.streams[slug].count += 1
    colons = (Colon() for i in 1:log.streams[slug].rank)
    log.streams[slug].dataset[colons..., log.streams[slug].count] = data
end

function close!(log::Log)
    close(log.file_id)
end

end # module


text:
"""Boundary Bools"""

#=
   Boundary Bools
                                            lo  hi       lo     hi
        false  ⇋  a  Closed bound            ClOp   ⇋  (false, true)
        true   ⇋  an Open   bound            OpOp   ⇋  (true , true)
=#

const PairedBoundaryTypes = [ ClCl, OpCl, ClOp, OpOp ]

""" obtains a pair of boundary bools from a paired boundary type """
@inline boundaries(::Type{ClCl}) = (false, false)
@inline boundaries(::Type{ClOp}) = (false, true )
@inline boundaries(::Type{OpCl}) = (true,  false)
@inline boundaries(::Type{OpOp}) = (true,  true )

""" obtains a paired boundary type from a pair of boundary bools """
boundaries{B<:Bool}(lo::B,hi::B) =
   PairedBoundaryTypes[ one(Int8)+reinterpret(Int8,lo)+(reinterpret(Int8,hi)<<1) ]

""" obtains a paired boundary type from an Rvl [type] """
@inline boundaries{R<:Real}(::Type{Rvl{ClCl,R}}) = ClCl
@inline boundaries{R<:Real}(::Type{Rvl{ClOp,R}}) = ClOp
@inline boundaries{R<:Real}(::Type{Rvl{OpCl,R}}) = OpCl
@inline boundaries{R<:Real}(::Type{Rvl{OpOp,R}}) = OpOp
boundaries{R<:Real}(x::Rvl{ClCl,R}) = boundaries(Rvl{ClCl,R})
boundaries{R<:Real}(x::Rvl{ClOp,R}) = boundaries(Rvl{ClCl,R})
boundaries{R<:Real}(x::Rvl{OpCl,R}) = boundaries(Rvl{ClCl,R})
boundaries{R<:Real}(x::Rvl{OpOp,R}) = boundaries(Rvl{ClCl,R})

""" obtains a paired boundary type, the unary negation of a paired boundary type"""
@inline negate(::Type{ClCl}) = ClCl
@inline negate(::Type{ClOp}) = OpCl
@inline negate(::Type{OpCl}) = ClOp
@inline negate(::Type{OpOp}) = OpOp

""" obtains a Real type from an Rvl [type] """
@inline realtype{G<:Grasp,R<:Real}(::Type{Rvl{G,R}}) = R
realtype{G<:Grasp,R<:Real}(x::Rvl{G,R}) = realtype(Rvl{G,R})


text:
module BarycentricLagrange
# Reference:
# Barycentric Lagrange Interpolation.
# Berrut, J.-P., & Trefethen, L. N. (2004).
# SIAM Review, 46(3), 501–517. doi:10.1137/S36144502417716

export bclag_interp_eval
export bclag_interp_matrix
export bclag_interp_weights

using Compat

function bclag_interp_weight(x, j, n)
    p = 1.0
    for k=1:j-1
        p *= x[j] - x[k]
    end
    for k=j+1:n
        p *= x[j] - x[k]
    end
    return 1/p
end

"""
    w = bclag_interp_weights(x)

Compute weights for Barycentric Lagrange interpolation.
"""
function bclag_interp_weights(x::Vector)
    n = length(x)
    w = zeros(n)
    for j=1:n
        w[j] = bclag_interp_weight(x, j, n)
    end
    return w
end

function bclag_interp_eval(x::Vector, f::Vector{T}, xx::Number,
                           w::Vector = bclag_interp_weights(x)) where T
    n = length(x)
    @assert length(f) == length(w) == n
    numer = zero(T)
    denom = 0.0
    for j=1:n
        xdiff = xx-x[j]
        if xdiff==0
            return f[j]
        end
        tmp = w[j] / xdiff
        numer += f[j]*tmp
        denom += tmp
    end
    ff = numer/denom
    return ff
end

"""
    ff = bclag_interp_eval(x, f, xx [, w])

Barycentric Lagrange interpolation of `f` from points `x` to points `xx`

`w` is vector of weights, `w = bclag_interp_weights(x)`
"""
function bclag_interp_eval(x::Vector, f::Vector{T}, xx::Vector,
                           w::Vector = bclag_interp_weights(x)) where T
    n = length(xx)
    ff = Array{T}(undef, n)
    for j=1:n
        ff[j] = bclag_interp_eval(x, f, xx[j], w)
    end
    return ff
end

"""
    M = bclag_interp_matrix(x, xx [, w])

Compute the matrix `M` such that `ff = M*f` gives the barycentric Lagrange interpolation of `f` from points `x` to points `xx`.

`w` is vector of weights, `w = bclag_interp_weights(x)`
"""
function bclag_interp_matrix(x::Vector, xx::Vector,
                             w::Vector = bclag_interp_weights(x))
    n = length(x)
    N = length(xx)
    @assert length(w)==n
    B, denom, exact = matrix_alloc(n, N)
    compute_matrix!(B, denom, exact, x, xx, w)
    return B
end

function matrix_alloc(n, N)
    B = Array{Float64}(undef, N, n)
    denom = Array{Float64}(undef, N)
    exact = Array{Int64}(undef, N)
    return B, denom, exact
end

function compute_matrix!(B, denom, exact, x, xx, w)
    # No allocations or size checks done here
    N, n = size(B)
    fill!(denom, 0.0)
    fill!(exact, 0)
    for j=1:n
        for k=1:N
            xdiff = xx[k]-x[j]
            if xdiff != 0
                temp = w[j]/xdiff
                B[k,j] = temp
                denom[k] += temp
            else
                exact[k] = j
            end
        end
    end
    B ./= denom
    for jj=1:N
        if exact[jj] != 0
            @. B[jj,:] = 0.0
            B[jj + N*(exact[jj]-1)] = 1.0
        end
    end
end

end # module


text:
# ---
# title: 1071. Greatest Common Divisor of Strings
# id: problem1071
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/greatest-common-divisor-of-strings/description/>
# hidden: true
# ---
#
# For two strings `s` and `t`, we say "`t` divides `s`" if and only if `s = t +
# ... + t`  (`t` concatenated with itself 1 or more times)
#
# Given two strings str1 and str2, return the largest string `x` such that `x`
# divides both `str1` and `str2`.
#
#
#
# **Example 1:**
#
#
#
#     Input: str1 = "ABCABC", str2 = "ABC"
#     Output: "ABC"
#
#
# **Example 2:**
#
#
#
#     Input: str1 = "ABABAB", str2 = "ABAB"
#     Output: "AB"
#
#
# **Example 3:**
#
#
#
#     Input: str1 = "LEET", str2 = "CODE"
#     Output: ""
#
#
# **Example 4:**
#
#
#
#     Input: str1 = "ABCDEF", str2 = "ABC"
#     Output: ""
#
#
#
#
# **Constraints:**
#
#   * `1 <= str1.length <= 1000`
#   * `1 <= str2.length <= 1000`
#   * `str1` and `str2` consist of English uppercase letters.
#
#
## @lc code=start
using LeetCode

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


text:
# Julia wrapper for header: /Volumes/BOOTCAMP/programs/compa_libs/iup/include/wd.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@static Sys.iswindows() ? (const libwd_ = "cd") : (const libwd_ = "libcd")  # Name of IUP shared lib.

# These used to be in wd_h.jl (meanwhile removed)
const wdVectorTextTransform = cdVectorTextTransform
const wdActivate = cdActivate
# Skipping MacroDefinition: wdClip(mode)cdClip(mode)
# Skipping MacroDefinition: wdBegin(mode)cdBegin(mode)
# Skipping MacroDefinition: wdEnd()cdEnd();
# Skipping MacroDefinition: wdMM2Pixel(mm_dx,mm_dy,dx,dy)cdMM2Pixel(mm_dx,mm_dy,dx,dy)
# Skipping MacroDefinition: wdPixel2MM(dx,dy,mm_dx,mm_dy)cdPixel2MM(dx,dy,mm_dx,mm_dy)

function wdCanvasWindow(canvas::Ptr{cdCanvas}, xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdCanvasWindow, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xmin, xmax, ymin, ymax)
end

function wdCanvasGetWindow(canvas::Ptr{cdCanvas}, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdCanvasGetWindow, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, xmin, xmax, ymin, ymax)
end

function wdCanvasViewport(canvas::Ptr{cdCanvas}, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdCanvasViewport, libwd_), Cvoid, (Ptr{cdCanvas}, Cint, Cint, Cint, Cint), canvas, xmin, xmax, ymin, ymax)
end

function wdCanvasGetViewport(canvas::Ptr{cdCanvas}, xmin::Ptr{Cint}, xmax::Ptr{Cint}, ymin::Ptr{Cint}, ymax::Ptr{Cint})
  ccall( (:wdCanvasGetViewport, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), canvas, xmin, xmax, ymin, ymax)
end

function wdCanvasWorld2Canvas(canvas::Ptr{cdCanvas}, xw::Cdouble, yw::Cdouble, xv::Ptr{Cint}, yv::Ptr{Cint})
  ccall( (:wdCanvasWorld2Canvas, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{Cint}, Ptr{Cint}), canvas, xw, yw, xv, yv)
end

function wdCanvasWorld2CanvasSize(canvas::Ptr{cdCanvas}, hw::Cdouble, vw::Cdouble, hv::Ptr{Cint}, vv::Ptr{Cint})
  ccall( (:wdCanvasWorld2CanvasSize, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{Cint}, Ptr{Cint}), canvas, hw, vw, hv, vv)
end

function wdCanvasCanvas2World(canvas::Ptr{cdCanvas}, xv::Integer, yv::Integer, xw::Ptr{Cdouble}, yw::Ptr{Cdouble})
  ccall( (:wdCanvasCanvas2World, libwd_), Cvoid, (Ptr{cdCanvas}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}), canvas, xv, yv, xw, yw)
end

function wdCanvasSetTransform(canvas::Ptr{cdCanvas}, sx::Cdouble, sy::Cdouble, tx::Cdouble, ty::Cdouble)
  ccall( (:wdCanvasSetTransform, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, sx, sy, tx, ty)
end

function wdCanvasGetTransform(canvas::Ptr{cdCanvas}, sx::Ptr{Cdouble}, sy::Ptr{Cdouble}, tx::Ptr{Cdouble}, ty::Ptr{Cdouble})
  ccall( (:wdCanvasGetTransform, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, sx, sy, tx, ty)
end

function wdCanvasTranslate(canvas::Ptr{cdCanvas}, dtx::Cdouble, dty::Cdouble)
  ccall( (:wdCanvasTranslate, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, dtx, dty)
end

function wdCanvasScale(canvas::Ptr{cdCanvas}, dsx::Cdouble, dsy::Cdouble)
  ccall( (:wdCanvasScale, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, dsx, dsy)
end

function wdCanvasClipArea(canvas::Ptr{cdCanvas}, xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdCanvasClipArea, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xmin, xmax, ymin, ymax)
end

function wdCanvasGetClipArea(canvas::Ptr{cdCanvas}, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdCanvasGetClipArea, libwd_), Cint, (Ptr{cdCanvas}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, xmin, xmax, ymin, ymax)
end

function wdCanvasIsPointInRegion(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble)
  ccall( (:wdCanvasIsPointInRegion, libwd_), Cint, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, x, y)
end

function wdCanvasOffsetRegion(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble)
  ccall( (:wdCanvasOffsetRegion, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, x, y)
end

function wdCanvasGetRegionBox(canvas::Ptr{cdCanvas}, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdCanvasGetRegionBox, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, xmin, xmax, ymin, ymax)
end

function wdCanvasHardcopy(canvas::Ptr{cdCanvas}, ctx::Ptr{Cint}, data::Ptr{Cvoid}, draw_func::Ptr{Cvoid})
  ccall( (:wdCanvasHardcopy, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{Cint}, Ptr{Cvoid}, Ptr{Cvoid}), canvas, ctx, data, draw_func)
end

function wdCanvasPixel(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, color::Clong)
  ccall( (:wdCanvasPixel, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Clong), canvas, x, y, color)
end

function wdCanvasMark(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble)
  ccall( (:wdCanvasMark, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, x, y)
end

function wdCanvasLine(canvas::Ptr{cdCanvas}, x1::Cdouble, y1::Cdouble, x2::Cdouble, y2::Cdouble)
  ccall( (:wdCanvasLine, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, x1, y1, x2, y2)
end

function wdCanvasVertex(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble)
  ccall( (:wdCanvasVertex, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble), canvas, x, y)
end

function wdCanvasRect(canvas::Ptr{cdCanvas}, xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdCanvasRect, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xmin, xmax, ymin, ymax)
end

function wdCanvasBox(canvas::Ptr{cdCanvas}, xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdCanvasBox, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xmin, xmax, ymin, ymax)
end

function wdCanvasArc(canvas::Ptr{cdCanvas}, xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdCanvasArc, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xc, yc, w, h, angle1, angle2)
end

function wdCanvasSector(canvas::Ptr{cdCanvas}, xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdCanvasSector, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xc, yc, w, h, angle1, angle2)
end

function wdCanvasChord(canvas::Ptr{cdCanvas}, xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdCanvasChord, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), canvas, xc, yc, w, h, angle1, angle2)
end

function wdCanvasText(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdCanvasText, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}), canvas, x, y, s)
end

function wdCanvasPutImageRect(canvas::Ptr{cdCanvas}, image::Ptr{Cint}, x::Cdouble, y::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdCanvasPutImageRect, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{Cint}, Cdouble, Cdouble, Cint, Cint, Cint, Cint), canvas, image, x, y, xmin, xmax, ymin, ymax)
end

function wdCanvasPutImageRectRGB(canvas::Ptr{cdCanvas}, iw::Integer, ih::Integer, r::Ptr{Cuchar}, g::Ptr{Cuchar}, b::Ptr{Cuchar}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdCanvasPutImageRectRGB, libwd_), Cvoid, (Ptr{cdCanvas}, Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), canvas, iw, ih, r, g, b, x, y, w, h, xmin, xmax, ymin, ymax)
end

function wdCanvasPutImageRectRGBA(canvas::Ptr{cdCanvas}, iw::Integer, ih::Integer, r::Ptr{Cuchar}, g::Ptr{Cuchar}, b::Ptr{Cuchar}, a::Ptr{Cuchar}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdCanvasPutImageRectRGBA, libwd_), Cvoid, (Ptr{cdCanvas}, Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), canvas, iw, ih, r, g, b, a, x, y, w, h, xmin, xmax, ymin, ymax)
end

function wdCanvasPutImageRectMap(canvas::Ptr{cdCanvas}, iw::Integer, ih::Integer, index::Ptr{Cuchar}, colors::Ptr{Clong}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdCanvasPutImageRectMap, libwd_), Cvoid, (Ptr{cdCanvas}, Cint, Cint, Ptr{Cuchar}, Ptr{Clong}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), canvas, iw, ih, index, colors, x, y, w, h, xmin, xmax, ymin, ymax)
end

function wdCanvasPutBitmap(canvas::Ptr{cdCanvas}, bitmap::Ptr{Cint}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble)
  ccall( (:wdCanvasPutBitmap, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{Cint}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, bitmap, x, y, w, h)
end

function wdCanvasLineWidth(canvas::Ptr{cdCanvas}, width::Cdouble)
  ccall( (:wdCanvasLineWidth, libwd_), Cdouble, (Ptr{cdCanvas}, Cdouble), canvas, width)
end

function wdCanvasFont(canvas::Ptr{cdCanvas}, type_face::String, style::Integer, size::Cdouble)
  ccall( (:wdCanvasFont, libwd_), Cint, (Ptr{cdCanvas}, Ptr{UInt8}, Cint, Cdouble), canvas, type_face, style, size)
end

function wdCanvasGetFont(canvas::Ptr{cdCanvas}, type_face::String, style::Ptr{Cint}, size::Ptr{Cdouble})
  ccall( (:wdCanvasGetFont, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{UInt8}, Ptr{Cint}, Ptr{Cdouble}), canvas, type_face, style, size)
end

function wdCanvasMarkSize(canvas::Ptr{cdCanvas}, size::Cdouble)
  ccall( (:wdCanvasMarkSize, libwd_), Cdouble, (Ptr{cdCanvas}, Cdouble), canvas, size)
end

function wdCanvasGetFontDim(canvas::Ptr{cdCanvas}, max_width::Ptr{Cdouble}, height::Ptr{Cdouble}, ascent::Ptr{Cdouble}, descent::Ptr{Cdouble})
  ccall( (:wdCanvasGetFontDim, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, max_width, height, ascent, descent)
end

function wdCanvasGetTextSize(canvas::Ptr{cdCanvas}, s::String, width::Ptr{Cdouble}, height::Ptr{Cdouble})
  ccall( (:wdCanvasGetTextSize, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, s, width, height)
end

function wdCanvasGetTextBox(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdCanvasGetTextBox, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, x, y, s, xmin, xmax, ymin, ymax)
end

function wdCanvasGetTextBounds(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String, rect::Ptr{Cdouble})
  ccall( (:wdCanvasGetTextBounds, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}, Ptr{Cdouble}), canvas, x, y, s, rect)
end

function wdCanvasStipple(canvas::Ptr{cdCanvas}, w::Integer, h::Integer, fgbg::Ptr{Cuchar}, w_mm::Cdouble, h_mm::Cdouble)
  ccall( (:wdCanvasStipple, libwd_), Cvoid, (Ptr{cdCanvas}, Cint, Cint, Ptr{Cuchar}, Cdouble, Cdouble), canvas, w, h, fgbg, w_mm, h_mm)
end

function wdCanvasPattern(canvas::Ptr{cdCanvas}, w::Integer, h::Integer, color::Ptr{Clong}, w_mm::Cdouble, h_mm::Cdouble)
  ccall( (:wdCanvasPattern, libwd_), Cvoid, (Ptr{cdCanvas}, Cint, Cint, Ptr{Clong}, Cdouble, Cdouble), canvas, w, h, color, w_mm, h_mm)
end

function wdCanvasVectorTextDirection(canvas::Ptr{cdCanvas}, x1::Cdouble, y1::Cdouble, x2::Cdouble, y2::Cdouble)
  ccall( (:wdCanvasVectorTextDirection, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Cdouble, Cdouble), canvas, x1, y1, x2, y2)
end

function wdCanvasVectorTextSize(canvas::Ptr{cdCanvas}, size_x::Cdouble, size_y::Cdouble, s::String)
  ccall( (:wdCanvasVectorTextSize, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}), canvas, size_x, size_y, s)
end

function wdCanvasGetVectorTextSize(canvas::Ptr{cdCanvas}, s::String, x::Ptr{Cdouble}, y::Ptr{Cdouble})
  ccall( (:wdCanvasGetVectorTextSize, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, s, x, y)
end

function wdCanvasVectorCharSize(canvas::Ptr{cdCanvas}, size::Cdouble)
  ccall( (:wdCanvasVectorCharSize, libwd_), Cdouble, (Ptr{cdCanvas}, Cdouble), canvas, size)
end

function wdCanvasVectorText(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdCanvasVectorText, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}), canvas, x, y, s)
end

function wdCanvasMultiLineVectorText(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdCanvasMultiLineVectorText, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}), canvas, x, y, s)
end

function wdCanvasGetVectorTextBounds(canvas::Ptr{cdCanvas}, s::String, x::Cdouble, y::Cdouble, rect::Ptr{Cdouble})
  ccall( (:wdCanvasGetVectorTextBounds, libwd_), Cvoid, (Ptr{cdCanvas}, Ptr{UInt8}, Cdouble, Cdouble, Ptr{Cdouble}), canvas, s, x, y, rect)
end

function wdCanvasGetVectorTextBox(canvas::Ptr{cdCanvas}, x::Cdouble, y::Cdouble, s::String, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdCanvasGetVectorTextBox, libwd_), Cvoid, (Ptr{cdCanvas}, Cdouble, Cdouble, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), canvas, x, y, s, xmin, xmax, ymin, ymax)
end

function wdWindow(xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdWindow, libwd_), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble), xmin, xmax, ymin, ymax)
end

function wdGetWindow(xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdGetWindow, libwd_), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), xmin, xmax, ymin, ymax)
end

function wdViewport(xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdViewport, libwd_), Cvoid, (Cint, Cint, Cint, Cint), xmin, xmax, ymin, ymax)
end

function wdGetViewport(xmin::Ptr{Cint}, xmax::Ptr{Cint}, ymin::Ptr{Cint}, ymax::Ptr{Cint})
  ccall( (:wdGetViewport, libwd_), Cvoid, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), xmin, xmax, ymin, ymax)
end

function wdWorld2Canvas(xw::Cdouble, yw::Cdouble, xv::Ptr{Cint}, yv::Ptr{Cint})
  ccall( (:wdWorld2Canvas, libwd_), Cvoid, (Cdouble, Cdouble, Ptr{Cint}, Ptr{Cint}), xw, yw, xv, yv)
end

function wdWorld2CanvasSize(hw::Cdouble, vw::Cdouble, hv::Ptr{Cint}, vv::Ptr{Cint})
  ccall( (:wdWorld2CanvasSize, libwd_), Cvoid, (Cdouble, Cdouble, Ptr{Cint}, Ptr{Cint}), hw, vw, hv, vv)
end

function wdCanvas2World(xv::Integer, yv::Integer, xw::Ptr{Cdouble}, yw::Ptr{Cdouble})
  ccall( (:wdCanvas2World, libwd_), Cvoid, (Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}), xv, yv, xw, yw)
end

function wdClipArea(xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdClipArea, libwd_), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble), xmin, xmax, ymin, ymax)
end

function wdGetClipArea(xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdGetClipArea, libwd_), Cint, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), xmin, xmax, ymin, ymax)
end

function wdGetClipPoly(n::Ptr{Cint})
  ccall( (:wdGetClipPoly, libwd_), Ptr{Cdouble}, (Ptr{Cint},), n)
end

function wdPointInRegion(x::Cdouble, y::Cdouble)
  ccall( (:wdPointInRegion, libwd_), Cint, (Cdouble, Cdouble), x, y)
end

function wdOffsetRegion(x::Cdouble, y::Cdouble)
  ccall( (:wdOffsetRegion, libwd_), Cvoid, (Cdouble, Cdouble), x, y)
end

function wdRegionBox(xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdRegionBox, libwd_), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), xmin, xmax, ymin, ymax)
end

function wdHardcopy(ctx::Ptr{Cint}, data::Ptr{Cvoid}, cnv::Ptr{Cint}, draw_func::Ptr{Cvoid})
  ccall( (:wdHardcopy, libwd_), Cvoid, (Ptr{Cint}, Ptr{Cvoid}, Ptr{Cint}, Ptr{Cvoid}), ctx, data, cnv, draw_func)
end

function wdPixel(x::Cdouble, y::Cdouble, color::Clong)
  ccall( (:wdPixel, libwd_), Cvoid, (Cdouble, Cdouble, Clong), x, y, color)
end

function wdMark(x::Cdouble, y::Cdouble)
  ccall( (:wdMark, libwd_), Cvoid, (Cdouble, Cdouble), x, y)
end

function wdLine(x1::Cdouble, y1::Cdouble, x2::Cdouble, y2::Cdouble)
  ccall( (:wdLine, libwd_), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble), x1, y1, x2, y2)
end

function wdVertex(x::Cdouble, y::Cdouble)
  ccall( (:wdVertex, libwd_), Cvoid, (Cdouble, Cdouble), x, y)
end

function wdRect(xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdRect, libwd_), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble), xmin, xmax, ymin, ymax)
end

function wdBox(xmin::Cdouble, xmax::Cdouble, ymin::Cdouble, ymax::Cdouble)
  ccall( (:wdBox, libwd_), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble), xmin, xmax, ymin, ymax)
end

function wdArc(xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdArc, libwd_), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), xc, yc, w, h, angle1, angle2)
end

function wdSector(xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdSector, libwd_), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), xc, yc, w, h, angle1, angle2)
end

function wdChord(xc::Cdouble, yc::Cdouble, w::Cdouble, h::Cdouble, angle1::Cdouble, angle2::Cdouble)
  ccall( (:wdChord, libwd_), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble), xc, yc, w, h, angle1, angle2)
end

function wdText(x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdText, libwd_), Cvoid, (Cdouble, Cdouble, Ptr{UInt8}), x, y, s)
end

function wdPutImageRect(image::Ptr{Cint}, x::Cdouble, y::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdPutImageRect, libwd_), Cvoid, (Ptr{Cint}, Cdouble, Cdouble, Cint, Cint, Cint, Cint), image, x, y, xmin, xmax, ymin, ymax)
end

function wdPutImageRectRGB(iw::Integer, ih::Integer, r::Ptr{Cuchar}, g::Ptr{Cuchar}, b::Ptr{Cuchar}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdPutImageRectRGB, libwd_), Cvoid, (Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), iw, ih, r, g, b, x, y, w, h, xmin, xmax, ymin, ymax)
end

function wdPutImageRectRGBA(iw::Integer, ih::Integer, r::Ptr{Cuchar}, g::Ptr{Cuchar}, b::Ptr{Cuchar}, a::Ptr{Cuchar}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdPutImageRectRGBA, libwd_), Cvoid, (Cint, Cint, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), iw, ih, r, g, b, a, x, y, w, h, xmin, xmax, ymin, ymax)
end

function wdPutImageRectMap(iw::Integer, ih::Integer, index::Ptr{Cuchar}, colors::Ptr{Clong}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble, xmin::Integer, xmax::Integer, ymin::Integer, ymax::Integer)
  ccall( (:wdPutImageRectMap, libwd_), Cvoid, (Cint, Cint, Ptr{Cuchar}, Ptr{Clong}, Cdouble, Cdouble, Cdouble, Cdouble, Cint, Cint, Cint, Cint), iw, ih, index, colors, x, y, w, h, xmin, xmax, ymin, ymax)
end

function wdPutBitmap(bitmap::Ptr{Cint}, x::Cdouble, y::Cdouble, w::Cdouble, h::Cdouble)
  ccall( (:wdPutBitmap, libwd_), Cvoid, (Ptr{Cint}, Cdouble, Cdouble, Cdouble, Cdouble), bitmap, x, y, w, h)
end

function wdLineWidth(width::Cdouble)
  ccall( (:wdLineWidth, libwd_), Cdouble, (Cdouble,), width)
end

function wdFont(type_face::Integer, style::Integer, size::Cdouble)
  ccall( (:wdFont, libwd_), Cvoid, (Cint, Cint, Cdouble), type_face, style, size)
end

function wdGetFont(type_face::Ptr{Cint}, style::Ptr{Cint}, size::Ptr{Cdouble})
  ccall( (:wdGetFont, libwd_), Cvoid, (Ptr{Cint}, Ptr{Cint}, Ptr{Cdouble}), type_face, style, size)
end

function wdMarkSize(size::Cdouble)
  ccall( (:wdMarkSize, libwd_), Cdouble, (Cdouble,), size)
end

function wdFontDim(max_width::Ptr{Cdouble}, height::Ptr{Cdouble}, ascent::Ptr{Cdouble}, descent::Ptr{Cdouble})
  ccall( (:wdFontDim, libwd_), Cvoid, (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), max_width, height, ascent, descent)
end

function wdTextSize(s::String, width::Ptr{Cdouble}, height::Ptr{Cdouble})
  ccall( (:wdTextSize, libwd_), Cvoid, (Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}), s, width, height)
end

function wdTextBox(x::Cdouble, y::Cdouble, s::String, xmin::Ptr{Cdouble}, xmax::Ptr{Cdouble}, ymin::Ptr{Cdouble}, ymax::Ptr{Cdouble})
  ccall( (:wdTextBox, libwd_), Cvoid, (Cdouble, Cdouble, Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}), x, y, s, xmin, xmax, ymin, ymax)
end

function wdTextBounds(x::Cdouble, y::Cdouble, s::String, rect::Ptr{Cdouble})
  ccall( (:wdTextBounds, libwd_), Cvoid, (Cdouble, Cdouble, Ptr{UInt8}, Ptr{Cdouble}), x, y, s, rect)
end

function wdStipple(w::Integer, h::Integer, stipple::Ptr{Cuchar}, w_mm::Cdouble, h_mm::Cdouble)
  ccall( (:wdStipple, libwd_), Cvoid, (Cint, Cint, Ptr{Cuchar}, Cdouble, Cdouble), w, h, stipple, w_mm, h_mm)
end

function wdPattern(w::Integer, h::Integer, pattern::Ptr{Clong}, w_mm::Cdouble, h_mm::Cdouble)
  ccall( (:wdPattern, libwd_), Cvoid, (Cint, Cint, Ptr{Clong}, Cdouble, Cdouble), w, h, pattern, w_mm, h_mm)
end

function wdVectorTextDirection(x1::Cdouble, y1::Cdouble, x2::Cdouble, y2::Cdouble)
  ccall( (:wdVectorTextDirection, libwd_), Cvoid, (Cdouble, Cdouble, Cdouble, Cdouble), x1, y1, x2, y2)
end

function wdVectorTextSize(size_x::Cdouble, size_y::Cdouble, s::String)
  ccall( (:wdVectorTextSize, libwd_), Cvoid, (Cdouble, Cdouble, Ptr{UInt8}), size_x, size_y, s)
end

function wdGetVectorTextSize(s::String, x::Ptr{Cdouble}, y::Ptr{Cdouble})
  ccall( (:wdGetVectorTextSize, libwd_), Cvoid, (Ptr{UInt8}, Ptr{Cdouble}, Ptr{Cdouble}), s, x, y)
end

function wdVectorCharSize(size::Cdouble)
  ccall( (:wdVectorCharSize, libwd_), Cdouble, (Cdouble,), size)
end

function wdVectorText(x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdVectorText, libwd_), Cvoid, (Cdouble, Cdouble, Ptr{UInt8}), x, y, s)
end

function wdMultiLineVectorText(x::Cdouble, y::Cdouble, s::String)
  ccall( (:wdMultiLineVectorText, libwd_), Cvoid, (Cdouble, Cdouble, Ptr{UInt8}), x, y, s)
end

function wdGetVectorTextBounds(s::String, x::Cdouble, y::Cdouble, rect::Ptr{Cdouble})
  ccall( (:wdGetVectorTextBounds, libwd_), Cvoid, (Ptr{UInt8}, Cdouble, Cdouble, Ptr{Cdouble}), s, x, y, rect)
end


text:
function attach_constraints_rows(cellvec,cellconstr,cellmask=Fill(true,length(cellconstr)))
  lazy_map(ConstrainRowsMap(),cellvec,cellconstr,cellmask)
end

function attach_constraints_cols(cellmat,cellconstr,cellmask=Fill(true,length(cellconstr)))
  cellconstr_t = lazy_map(transpose,cellconstr)
  lazy_map(ConstrainColsMap(),cellmat,cellconstr_t,cellmask)
end

struct ConstrainRowsMap <: Map end

function Arrays.return_cache(k::ConstrainRowsMap,array,constr,mask)
  return_cache(*,constr,array)
end

@inline function Arrays.evaluate!(cache,k::ConstrainRowsMap,array,constr,mask)
  if mask
    evaluate!(cache,*,constr,array)
  else
    array
  end
end

function Arrays.return_cache(k::ConstrainRowsMap,matvec::Tuple,constr,mask)
  mat, vec = matvec
  cmat = return_cache(k,mat,constr,mask)
  cvec = return_cache(k,vec,constr,mask)
  (cmat,cvec)
end

@inline function Arrays.evaluate!(cache,k::ConstrainRowsMap,matvec::Tuple,constr,mask)
  if mask
    cmat, cvec = cache
    mat, vec = matvec
    _mat = evaluate!(cmat,k,mat,constr,mask)
    _vec = evaluate!(cvec,k,vec,constr,mask)
    (_mat,_vec)
  else
    matvec
  end
end

struct ConstrainColsMap <: Map end

function Arrays.return_cache(k::ConstrainColsMap,array,constr_t,mask)
  return_cache(*,array,constr_t)
end

@inline function Arrays.evaluate!(cache,k::ConstrainColsMap,array,constr_t,mask)
  if mask
    evaluate!(cache,*,array,constr_t)
  else
    array
  end
end

function Arrays.return_cache(k::ConstrainColsMap,matvec::Tuple,constr_t,mask)
  mat, vec = matvec
  return_cache(k,mat,constr_t,mask)
end

@inline function Arrays.evaluate!(cache,k::ConstrainColsMap,matvec::Tuple,constr_t,mask)
  if mask
    mat, vec = matvec
    _mat = evaluate!(cache,k,mat,constr_t,mask)
    (_mat,vec)
  else
    matvec
  end
end

function identity_constraints(cell_axes)
  lazy_map(IdentityConstraintMap(),cell_axes)
end

struct IdentityConstraintMap <: Map end

function Arrays.return_cache(k::IdentityConstraintMap,axs)
  n = length(axs[1])
  a = zeros(n,n)
  CachedArray(a)
end

function Arrays.evaluate!(cache,k::IdentityConstraintMap,axs)
  n = length(axs[1])
  setsize!(cache,(n,n))
  a = cache.array
  fill!(a,zero(eltype(a)))
  o = one(eltype(a))
  @inbounds for i in 1:size(a,1)
    a[i,i] = o
  end
  a
end


text:
using MathOptFormat
using MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

DIR_PATH = dirname(@__FILE__)

push!(Base.LOAD_PATH, joinpath(DIR_PATH, "..",".."))

using ModelCompare

cd(DIR_PATH)

model = MathOptFormat.MOF.Model()
x = MOI.add_variable(model)
MOI.set(model, MOI.VariableName(), x, "x")
c = MOI.add_constraint(model,
    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}.([1],[x]), 0.0),
    MOI.GreaterThan(1.0)
)
MOI.set(model, MOI.ConstraintName(), c, "c")
c2 = MOI.add_constraint(model,
    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0),
    MOI.Interval(-1.0, 1.0)
)
MOI.set(model, MOI.ConstraintName(), c2, "c2")
MOI.write_to_file(model, "orig_A.mps")

remove_intervals(model)
MOI.write_to_file(model, "orig_A_no_interval.mps")

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

model_ = MathOptFormat.MOF.Model()
z = MOI.add_variable(model_)
MOI.set(model_, MOI.VariableName(), z, "z")
x = MOI.add_variable(model_)
MOI.set(model_, MOI.VariableName(), x, "x")
c = MOI.add_constraint(model_,
    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}.([1,2],[z,x]), 0.0),
    MOI.GreaterThan(1.0)
)
MOI.set(model_, MOI.ConstraintName(), c, "c")
c2 = MOI.add_constraint(model_,
    MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0),
    MOI.Interval(-1.0, 1.0)
)
MOI.set(model_, MOI.ConstraintName(), c2, "c2")
MOI.write_to_file(model_, "orig_B.mps")
remove_intervals(model_)

compare_models(model, model_)


text:
using MechanicalUnits
using Test
@info "Starting test"
@testset "MechanicalUnits.jl" begin
    @testset "pretty-print.jl"  begin
        include("pretty-print.jl")
    end
    @testset "conversion_promotion.jl"  begin
        include("conversion_promotion.jl")
    end
    @testset "prefixes.jl"  begin
        include("prefixes.jl")
    end
end
nothing

max_stars_repo_name: hustf/MechanicalUnits

max_issues_repo_name: hustf/MechanicalUnits


text:
module BenchCartesian

using BenchmarkTools
using Transducers

function copyto_manual!(ys::AbstractMatrix, xs::AbstractMatrix)
    @assert axes(ys) == axes(xs)
    for j in 1:size(xs, 2), i in 1:size(xs, 1)
        @inbounds ys[i, j] = xs[i, j]
    end
    return ys
end

function copyto_iter!(ys, xs)
    @assert axes(ys) == axes(xs)
    for I in CartesianIndices(xs)
        @inbounds ys[I] = xs[I]
    end
    return ys
end

function copyto_xf!(ys, xs)
    foreach(Map(identity), CartesianIndices(xs)) do I
        @inbounds ys[I] = xs[I]
        nothing
    end
    return ys
end

const SUITE = BenchmarkGroup()

let xs = randn(3, 10^3)
    ys = zero(xs)
    s1 = SUITE["copyto!"] = BenchmarkGroup()
    s1["man"] = @benchmarkable(copyto_manual!($ys, $xs))
    s1["iter"] = @benchmarkable(copyto_iter!($ys, $xs))
    s1["xf"] = @benchmarkable(copyto_xf!($ys, $xs))
end

end  # module
BenchCartesian.SUITE


text:
"""
    tracenorm(rho)

Trace norm of `rho`.

It is defined as

```math
T(ρ) = Tr\\{\\sqrt{ρ^† ρ}\\}.
```

Depending if `rho` is hermitian either [`tracenorm_h`](@ref) or
[`tracenorm_nh`](@ref) is called.
"""
function tracenorm(rho::DenseOperator{B,B}) where B<:Basis
    ishermitian(rho) ? tracenorm_h(rho) : tracenorm_nh(rho)
end
function tracenorm(rho::T) where T<:AbstractOperator
    throw(ArgumentError("tracenorm not implemented for $(typeof(rho)). Use dense operators instead."))
end

"""
    tracenorm_h(rho)

Trace norm of `rho`.

It uses the identity

```math
T(ρ) = Tr\\{\\sqrt{ρ^† ρ}\\} = \\sum_i |λ_i|
```

where ``λ_i`` are the eigenvalues of `rho`.
"""
function tracenorm_h(rho::DenseOperator{B,B}) where B<:Basis
    s = eigvals(Hermitian(rho.data))
    sum(abs.(s))
end
function tracenorm_h(rho::T) where T<:AbstractOperator
    throw(ArgumentError("tracenorm_h not implemented for $(typeof(rho)). Use dense operators instead."))
end

"""
    tracenorm_nh(rho)

Trace norm of `rho`.

Note that in this case `rho` doesn't have to be represented by a square
matrix (i.e. it can have different left-hand and right-hand bases).

It uses the identity

```math
    T(ρ) = Tr\\{\\sqrt{ρ^† ρ}\\} = \\sum_i σ_i
```

where ``σ_i`` are the singular values of `rho`.
"""
tracenorm_nh(rho::DenseOperator) = sum(svdvals(rho.data))
function tracenorm_nh(rho::T) where T<:AbstractOperator
    throw(ArgumentError("tracenorm_nh not implemented for $(typeof(rho)). Use dense operators instead."))
end

"""
    tracedistance(rho, sigma)

Trace distance between `rho` and `sigma`.

It is defined as

```math
T(ρ,σ) = \\frac{1}{2} Tr\\{\\sqrt{(ρ - σ)^† (ρ - σ)}\\}.
```

It calls [`tracenorm`](@ref) which in turn either uses [`tracenorm_h`](@ref)
or [`tracenorm_nh`](@ref) depending if ``ρ-σ`` is hermitian or not.
"""
tracedistance(rho::T, sigma::T) where T<:DenseOperator = 0.5*tracenorm(rho - sigma)
function tracedistance(rho::AbstractOperator, sigma::AbstractOperator)
    throw(ArgumentError("tracedistance not implemented for $(typeof(rho)) and $(typeof(sigma)). Use dense operators instead."))
end

"""
    tracedistance_h(rho, sigma)

Trace distance between `rho` and `sigma`.

It uses the identity

```math
T(ρ,σ) = \\frac{1}{2} Tr\\{\\sqrt{(ρ - σ)^† (ρ - σ)}\\} = \\frac{1}{2} \\sum_i |λ_i|
```

where ``λ_i`` are the eigenvalues of `rho` - `sigma`.
"""
tracedistance_h(rho::T, sigma::T) where {B<:Basis,T<:DenseOperator{B,B}}= 0.5*tracenorm_h(rho - sigma)
function tracedistance_h(rho::AbstractOperator, sigma::AbstractOperator)
    throw(ArgumentError("tracedistance_h not implemented for $(typeof(rho)) and $(typeof(sigma)). Use dense operators instead."))
end

"""
    tracedistance_nh(rho, sigma)

Trace distance between `rho` and `sigma`.

Note that in this case `rho` and `sigma` don't have to be represented by square
matrices (i.e. they can have different left-hand and right-hand bases).

It uses the identity

```math
    T(ρ,σ) = \\frac{1}{2} Tr\\{\\sqrt{(ρ - σ)^† (ρ - σ)}\\}
         = \\frac{1}{2} \\sum_i σ_i
```

where ``σ_i`` are the singular values of `rho` - `sigma`.
"""
tracedistance_nh(rho::T, sigma::T) where {T<:DenseOperator} = 0.5*tracenorm_nh(rho - sigma)
function tracedistance_nh(rho::AbstractOperator, sigma::AbstractOperator)
    throw(ArgumentError("tracedistance_nh not implemented for $(typeof(rho)) and $(typeof(sigma)). Use dense operators instead."))
end

"""
    entropy_vn(rho)

Von Neumann entropy of a density matrix.

The Von Neumann entropy of a density operator is defined as

```math
S(ρ) = -Tr(ρ \\log(ρ)) = -\\sum_n λ_n\\log(λ_n)
```

where ``λ_n`` are the eigenvalues of the density matrix ``ρ``, ``\\log`` is the
natural logarithm and ``0\\log(0) ≡ 0``.

# Arguments
* `rho`: Density operator of which to calculate Von Neumann entropy.
* `tol=1e-15`: Tolerance for rounding errors in the computed eigenvalues.
"""
function entropy_vn(rho::DenseOperator{B,B}; tol::Float64=1e-15) where B<:Basis
    evals::Vector{ComplexF64} = eigvals(rho.data)
    evals[abs.(evals) .< tol] .= 0.0im
    sum([d == 0.0im ? 0.0im : -d*log(d) for d=evals])
end
entropy_vn(psi::StateVector; kwargs...) = entropy_vn(dm(psi); kwargs...)

"""
    fidelity(rho, sigma)

Fidelity of two density operators.

The fidelity of two density operators ``ρ`` and ``σ`` is defined by

```math
F(ρ, σ) = Tr\\left(\\sqrt{\\sqrt{ρ}σ\\sqrt{ρ}}\\right),
```

where ``\\sqrt{ρ}=\\sum_n\\sqrt{λ_n}|ψ⟩⟨ψ|``.
"""
fidelity(rho::T, sigma::T) where {B<:Basis,T<:DenseOperator{B,B}} = tr(sqrt(sqrt(rho.data)*sigma.data*sqrt(rho.data)))

"""
    ptranspose(rho, index)

Partial transpose of rho with respect to subspace specified by index.
"""
function ptranspose(rho::DenseOperator{B,B}, index::Int=1) where B<:CompositeBasis

    # Define permutation
    N = length(rho.basis_l.bases)
    perm = [1:N;]
    perm[index] = N
    perm[N] = index

    # Permute indexed subsystem to last position
    rho_perm = permutesystems(rho, perm)

    # Transpose corresponding blocks
    m = Int(prod(rho_perm.basis_l.shape[1:N-1]))
    n = rho_perm.basis_l.shape[N]
    for i=1:n, j=1:n
        rho_perm.data[m*(i-1)+1:m*i, m*(j-1)+1:m*j] = permutedims(rho_perm.data[m*(i-1)+1:m*i, m*(j-1)+1:m*j])
    end

    return permutesystems(rho_perm, perm)
end

"""
    PPT(rho, index)

Peres-Horodecki criterion of partial transpose.
"""
PPT(rho::DenseOperator{B,B}, index::Int) where B<:CompositeBasis = all(real.(eigvals(ptranspose(rho, index).data)) .>= 0.0)

"""
    negativity(rho, index)

Negativity of rho with respect to subsystem index.

The negativity of a density matrix ρ is defined as

```math
N(ρ) = \\|ρᵀ\\|,
```
where `ρᵀ` is the partial transpose.
"""
negativity(rho::DenseOperator{B,B}, index::Int) where B<:CompositeBasis = 0.5*(tracenorm(ptranspose(rho, index)) - 1.0)

"""
    logarithmic_negativity(rho, index)

The logarithmic negativity of a density matrix ρ is defined as

```math
N(ρ) = \\log₂\\|ρᵀ\\|,
```
where `ρᵀ` is the partial transpose.
"""
logarithmic_negativity(rho::DenseOperator{B,B}, index::Int) where B<:CompositeBasis = log(2, tracenorm(ptranspose(rho, index)))

"""
    avg_gate_fidelity(x, y)

The average gate fidelity between two superoperators x and y.
"""
function avg_gate_fidelity(x::T, y::T) where T <: Union{PauliTransferMatrix{B, B} where B, SuperOperator{B, B} where B, ChiMatrix{B, B} where B}
    dim = 2 ^ length(x.basis_l)
    return (tr(transpose(x.data) * y.data) + dim) / (dim^2 + dim)
end

"""
    entanglement_entropy(state, partition, [entropy_fun=entropy_vn])

Computes the entanglement entropy of `state` between the list of sites `partition`
and the rest of the system. The state must be defined in a composite basis.

If `state isa AbstractOperator` the operator-space entanglement entropy is
computed, which has the property
```julia
entanglement_entropy(dm(ket)) = 2 * entanglement_entropy(ket)
```

By default the computed entropy is the Von-Neumann entropy, but a different
function can be provided (for example to compute the entanglement-renyi entropy).
"""
function entanglement_entropy(psi::Ket{B}, partition::Vector{Int}, entropy_fun=entropy_vn) where B<:CompositeBasis
    # check that sites are within the range
    @assert all(partition .<= length(psi.basis.bases))

    rho = ptrace(psi, partition)
    return entropy_fun(rho)
end

function entanglement_entropy(rho::DenseOperator{B,B}, partition::Vector{Int}, args...) where {B<:CompositeBasis}
    # check that sites is within the range
    hilb = rho.basis_l
    all(partition .<= length(hilb.bases)) || throw(ArgumentError("Indices in partition must be within the bounds of the composite basis."))
    length(partition) <= length(hilb.bases) || throw(ArgumentError("Partition cannot include the whole system."))

    # build the doubled hilbert space for the vectorised dm, normalized like a Ket.
    b_doubled = hilb^2
    rho_vec = normalize!(Ket(b_doubled, vec(rho.data)))

    return entanglement_entropy(rho_vec,
                                vcat(partition, partition.+length(hilb.bases)),
                                args...)
end

entanglement_entropy(state, partition::Number, args...) =
    entanglement_entropy(state, [partition], args...)
entanglement_entropy(state, partition, args...) =
    entanglement_entropy(state, collect(partition), args...)


text:
using Test
using EmergencyBrakingSystem
using AutomotiveDrivingModels
using AutoViz
using AutomotiveSensors
using AutomotivePOMDPs
using Parameters
using StaticArrays
using Reel
using Random

@testset "EmergencyBrakingSystem.jl" begin

Random.seed!(1)

ped_v = 5.0 / 3.6
ped_theta = π/2
ped_right_side = true
hitpoint = 0
ped_t_collision = 1.0
ego_y = 0.0
ped_x = 100.0
ped_y = 0.0

obstacle_right_offset = -VehicleDef().width/2 - 1.0
obstacle_1 = ConvexPolygon([VecE2(44, obstacle_right_offset), VecE2(44, obstacle_right_offset-2), VecE2(48.5, obstacle_right_offset-2), VecE2(48.5, obstacle_right_offset)],4)
obstacles = [obstacle_1]

# unavoidable situation
ego_v = 50.0 / 3.6
ped_v = 0.0
ego_x = ped_x - ego_v * ped_t_collision - VehicleDef().length/2;
(scenes, timestep, env, ego_vehicle, sensor, sensor_observations, risk, ttc, collision_rate, brake_request, prediction, collision) = EmergencyBrakingSystem.evaluate_scenario(ego_x, ego_y, ego_v, ped_x, ped_y, ped_v, ped_theta, obstacles);
@test collision[end] == true

# avoidable situation
ego_v = 30.0 / 3.6
ego_x = ped_x - ego_v * ped_t_collision - VehicleDef().length/2;
(scenes, timestep, env, ego_vehicle, sensor, sensor_observations, risk, ttc, collision_rate, brake_request, prediction, collision) = EmergencyBrakingSystem.evaluate_scenario(ego_x, ego_y, ego_v, ped_x, ped_y, ped_v, ped_theta, obstacles);
@test collision[end] == false

end


text:
"""
    obstacles
"""

struct ObstacleConstraints <: Constraints
    n
    ineq
    n_stage
end

function constraints!(c, Z, con::ObstacleConstraints, model, idx, h, T)
    for t = 1:T
        x = view(Z, idx.x[t])
        obstacles!(view(c, (t - 1) * con.n_stage .+ (1:con.n_stage)), x)
    end
    nothing
end

function constraints_jacobian!(∇c, Z, con::ObstacleConstraints, model, idx, h, T)
    n = model.n

    c_tmp = zeros(con.n_stage)
    shift = 0

    for t = 1:T
        r_idx = (t - 1) * con.n_stage .+ (1:con.n_stage)
        c_idx = idx.x[t]

        len = length(r_idx) * length(c_idx)
        x = view(Z, idx.x[t])

        ForwardDiff.jacobian!(reshape(view(∇c, shift .+ (1:len)),
            con.n_stage, n),
            obstacles!, c_tmp, x)

        shift += len
    end

    return nothing
end

function constraints_sparsity(con::ObstacleConstraints, model, idx, T;
    shift_row = 0, shift_col = 0)
    n = model.n

    row = []
    col = []

    for t = 1:T
        r_idx = shift_row + (t - 1) * con.n_stage .+ (1:con.n_stage)
        c_idx = shift_col .+ idx.x[t]

        row_col!(row, col, r_idx, c_idx)
    end

    return collect(zip(row, col))
end

function circle_obs(x, y, xc, yc, r)
    (x - xc)^2.0 + (y - yc)^2.0 - r^2.0
end

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
"""
    randSeisSrc()

Generate a SeisSrc structure filled with random values.
"""
function randSeisSrc(; mw::Float32=0.0f0)
  R = SeisSrc()

  m0 = 1.0e-7*(10^(1.5*((mw == zero(Float32) ? 10^(rand(Float64)-1.0) : Float64(mw))+10.7)))

  setfield!(R, :id, string(rand(1:2^12)))                           # :id
  setfield!(R, :m0, m0)                                             # :m0
  setfield!(R, :gap, abs(rand_lon()))                               # :gap
  setfield!(R, :misc, rand_misc(rand(4:24)))                        # :misc
  setfield!(R, :mt, m0.*(rand(6).-0.5))                             # :mt
  setfield!(R, :dm, m0.*(rand(6).-0.5))                             # :dm
  setfield!(R, :npol, rand(Int64(6):Int64(120)))                    # :npol
  [note!(R, randstring(rand(16:256))) for i = 1:rand(1:6)]          # :notes
  setfield!(R, :pax, vcat(rand_lon(2,3),
      (m0*rand() < 0.5 ? -1 : 1) .* [rand() -1.0*rand() rand()]))   # :pax
  setfield!(R, :planes, rand_lon(3,2))                              # :planes
  setfield!(R, :src, R.id * ",randSeisSrc")                         # :src
  note!(R, "+origin ¦ " * R.src)
  setfield!(R, :st, SourceTime(randstring(2^rand(6:8)), rand(), rand(), rand()))
  return R
end

"""
    randSeisHdr()

Generate a SeisHdr structure filled with random values.
"""
function randSeisHdr()
  H = SeisHdr()

  # moment magnitude and m0 in SI units
  mw = Float32(log10(10.0^rand(0.0:0.1:7.0)))
  m0 = 1.0e-7*(10^(1.5*(Float64(mw)+10.7)))

  # a random earthquake location
  dmin = 10.0^(rand()+0.5)
  dmax = dmin + exp(mw)*10.0^rand()

  loc = EQLoc(
           (rand(0.0:1.0:89.0) + rand())*-1.0^rand(1:2),
           (rand(0.0:1.0:179.0)+rand())*-1.0^rand(1:2),
           min(10.0*randexp(Float64), 660.0),
           rand()*4.0,
           rand()*4.0,
           rand()*8.0,
           rand(),
           rand(),
           rand(),
           10.0*10.0^rand(),
           dmin,
           dmax,
           4 + round(Int64, rand()*exp(mw)),
           rand(0x00:0x10:0xf0),
           rand_datum(),
           rand(loc_types),
           "",
           rand(loc_methods))

  # setfield!(loc, :sig, "1σ")
  # event type
  rtyp = rand()
  typ = rtyp > 0.3 ? "earthquake" : rand(evtypes)

  # Generate random header
  setfield!(H, :id, string(rand(1:2^62)))                           # :id
  setfield!(H, :int, (mw < 0.0f0 ? 0x00 : floor(UInt8, mw),
    randstring(rand(2:4))))                                         # :int
  setfield!(H, :loc, loc)                                           # :loc
  setfield!(H, :mag, EQMag( val   = mw,          # :mag
                            scale = (rand() < 0.5 ? "m" : "M") * randstring(2),
                            nst   = rand(3:100),
                            gap   = abs(rand_lon()),
                            src   = randstring(24)
                            )
            )
  setfield!(H, :misc, rand_misc(rand(4:24)))                        # :misc
  [note!(H, randstring(rand(16:256))) for i = 1:rand(1:6)]          # :notes
  setfield!(H, :ot, now())                                          # :ot
  setfield!(H, :typ, typ)                                           # :typ

  # header :src
  H.src = "randSeisHdr:" * H.id
  note!(H, "+origin ¦ " * H.src)
  return H
end

"""
    randPhaseCat()

Generate a random seismic phase catalog suitable for testing EventChannel,
EventTraceData, and SeisEvent objects.
"""
function randPhaseCat(n::Int64)
  npha = (n <= 0) ? rand(3:18) : n
  phase = Array{String, 1}(undef, npha)
  for j = 1:length(phase)
    phase[j] = rand(phase_list)
  end
  unique!(phase)

  P = PhaseCat()
  for j in phase
    P[j] = SeisPha(rand(), rand(), rand(), rand(), rand(), rand(), rand(), rand(), rand(pol_list), Char(rand(0x30:0x39)))
  end
  return P
end
randPhaseCat() = randPhaseCat(0)

"""
    randSeisEvent([, c=0.2, s=0.6])

Generate a SeisEvent structure filled with random header and channel data.
* 100*c is the percentage of :data channels _after the first_ with irregularly-sampled data (fs = 0.0)
* 100*s is the percentage of :data channels _after the first_ with guaranteed seismic data.

See also: `randSeisChannel`, `randSeisData`, `randSeisHdr`, `randSeisSrc`
"""
function randSeisEvent(N::Int64; c::Float64=0.0, s::Float64=1.0, nx::Int64=0)
  V = SeisEvent(hdr=randSeisHdr(),
                source=randSeisSrc(),
                data=convert(EventTraceData, randSeisData(N, c=c, s=s, nx=nx)))
  V.source.eid = V.hdr.id
  V.hdr.loc.src = V.hdr.id * "," * V.hdr.loc.src

  for i = 1:V.data.n
    setindex!(getfield(getfield(V, :data), :pha),  randPhaseCat(),  i)
    setindex!(getfield(getfield(V, :data), :az),   rand_lat(),      i)
    setindex!(getfield(getfield(V, :data), :baz),  rand_lat(),      i)
    setindex!(getfield(getfield(V, :data), :dist), abs(rand_lon()), i)
  end
  return V
end
randSeisEvent(; c::Float64=0.0, s::Float64=1.0, nx::Int64=0) = randSeisEvent(rand(8:24), c=c, s=s, nx=nx)


text:
export LennardJonesModel,
       LennardJonesAtom

import CUDA

struct LennardJonesModel
    rc²::Float32
    rs²::Float32
    δ⁻²::Float32
    LennardJonesModel(cutoff, switch) = new(cutoff^2, switch^2, 1/(cutoff^2 - switch^2))
end

LennardJonesAtom(ε, σ) = LJAtom(0.5σ, 2*sqrt(ε))

struct LJAtom
    half_σ::Float32
    twice_sqrt_ε::Float32
end

@inline shfl_sync(mask, val::LJAtom, src) = LJAtom(
    CUDA.shfl_sync(mask, val.half_σ, src),
    CUDA.shfl_sync(mask, val.twice_sqrt_ε, src)
)

@inline function interaction(r²::Float32,
                             model::LennardJonesModel,
                             atom_i::LJAtom,
                             atom_j::LJAtom)::Tuple{Float32,Float32}
    σ = atom_i.half_σ + atom_j.half_σ
    ε4 = atom_i.twice_sqrt_ε*atom_j.twice_sqrt_ε
    s⁻² = σ*σ/r²
    s⁻⁶ = s⁻²*s⁻²*s⁻²
    ε4s⁻⁶ = atom_i.twice_sqrt_ε*atom_j.twice_sqrt_ε*s⁻⁶
    E = ε4s⁻⁶*(s⁻⁶ - 1)
    minus_E′r = 6ε4s⁻⁶*(2s⁻⁶ - 1)
    x = (r² - model.rs²)*model.δ⁻²
    x *= 0.5f0(sign(x) - sign(x-1))
    x² = x*x
    g = 1 + x*x²*(15x - 6x² - 10)
    minus_g′r = 60x²*(1 - 2x + x²)*model.δ⁻²*r²
    return E*g, minus_E′r*g + E*minus_g′r
end


text:
using HTTP, JSON

function yangsimoptapi(req::HTTP.Request)
    d=String(req.body)
    println(d)
    d = JSON.parse(d)
    #timeout = pop!(d, "timeout_seconds")
    #tol = pop!(d, "tolerance")

    @info "Starting YangSimOpt API ..."
	error_response = Dict()
	results = Dict()
    results = d
    println(d)

	if isempty(error_response)
    	@info "Hahaha."
    	return HTTP.Response(200, JSON.json(Dict("Jab"=>"dfa")))
	else
		@info "An error occured in the Julia code."
		return HTTP.Response(500, JSON.json(error_response))
	end
end

function health(req::HTTP.Request)
    return HTTP.Response(200, JSON.json(Dict("Julia-api"=>"healthy!")))
end

# define REST endpoints to dispatch to "service" functions
const ROUTER = HTTP.Router()

HTTP.@register(ROUTER, "POST", "/yangsimoptapi", yangsimoptapi)
HTTP.@register(ROUTER, "GET", "/health", health)
HTTP.serve(ROUTER, "0.0.0.0", 9527, reuseaddr=true)

max_stars_repo_name: mingtaoli/YangSimOptOnlineServer

max_issues_repo_name: mingtaoli/YangSimOptOnlineServer

max_forks_repo_name: mingtaoli/YangSimOptOnlineServer


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

@doc raw"""Adapts a secret into a projected volume.  The contents of the target Secret's Data field will be presented in a projected volume as files using the keys in the Data field as the file names. Note that this is identical to a secret volume source without the default mode.

    IoK8sApiCoreV1SecretProjection(;
        items=nothing,
        name=nothing,
        optional=nothing,
    )

    - items::Vector{IoK8sApiCoreV1KeyToPath} : If unspecified, each key-value pair in the Data field of the referenced Secret will be projected into the volume as a file whose name is the key and content is the value. If specified, the listed keys will be projected into the specified paths, and unlisted keys will not be present. If a key is specified which is not present in the Secret, the volume setup will error unless it is marked optional. Paths must be relative and may not contain the '..' path or start with '..'.
    - name::String : Name of the referent. More info: https://kubernetes.io/docs/concepts/overview/working-with-objects/names/#names
    - optional::Bool : Specify whether the Secret or its key must be defined
"""
mutable struct IoK8sApiCoreV1SecretProjection <: SwaggerModel
    items::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1KeyToPath} } # spec name: items
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    optional::Any # spec type: Union{ Nothing, Bool } # spec name: optional

    function IoK8sApiCoreV1SecretProjection(;items=nothing, name=nothing, optional=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1SecretProjection, Symbol("items"), items)
        setfield!(o, Symbol("items"), items)
        validate_property(IoK8sApiCoreV1SecretProjection, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiCoreV1SecretProjection, Symbol("optional"), optional)
        setfield!(o, Symbol("optional"), optional)
        o
    end
end # type IoK8sApiCoreV1SecretProjection

const _property_map_IoK8sApiCoreV1SecretProjection = Dict{Symbol,Symbol}(Symbol("items")=>Symbol("items"), Symbol("name")=>Symbol("name"), Symbol("optional")=>Symbol("optional"))
const _property_types_IoK8sApiCoreV1SecretProjection = Dict{Symbol,String}(Symbol("items")=>"Vector{IoK8sApiCoreV1KeyToPath}", Symbol("name")=>"String", Symbol("optional")=>"Bool")
Base.propertynames(::Type{ IoK8sApiCoreV1SecretProjection }) = collect(keys(_property_map_IoK8sApiCoreV1SecretProjection))
Swagger.property_type(::Type{ IoK8sApiCoreV1SecretProjection }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1SecretProjection[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1SecretProjection }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1SecretProjection[property_name]

function check_required(o::IoK8sApiCoreV1SecretProjection)
    true
end

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