text:
# ---
# title: 1401. Circle and Rectangle Overlapping
# id: problem1401
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Geometry
# link: <https://leetcode.com/problems/circle-and-rectangle-overlapping/description/>
# hidden: true
# ---
#
# Given a circle represented as (`radius`, `x_center`, `y_center`) and an axis-
# aligned rectangle represented as (`x1`, `y1`, `x2`, `y2`), where (`x1`, `y1`)
# are the coordinates of the bottom-left corner, and (`x2`, `y2`) are the
# coordinates of the top-right corner of the rectangle.
#
# Return True if the circle and rectangle are overlapped otherwise return False.
#
# In other words, check if there are **any** point (xi, yi) such that belongs to
# the circle and the rectangle at the same time.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/20/sample_4_1728.png)
#
#
#
#     Input: radius = 1, x_center = 0, y_center = 0, x1 = 1, y1 = -1, x2 = 3, y2 = 1
#     Output: true
#     Explanation: Circle and rectangle share the point (1,0)
#
#
# **Example 2:**
#
# **![](https://assets.leetcode.com/uploads/2020/02/20/sample_2_1728.png)**
#
#
#
#     Input: radius = 1, x_center = 0, y_center = 0, x1 = -1, y1 = 0, x2 = 0, y2 = 1
#     Output: true
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2020/03/03/sample_6_1728.png)**
#
#
#
#     Input: radius = 1, x_center = 1, y_center = 1, x1 = -3, y1 = -3, x2 = 3, y2 = 3
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: radius = 1, x_center = 1, y_center = 1, x1 = 1, y1 = -3, x2 = 2, y2 = -1
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= radius <= 2000`
#   * `-10^4 <= x_center, y_center, x1, y1, x2, y2 <= 10^4`
#   * `x1 < x2`
#   * `y1 < y2`
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
""" Data structure for views of attributed C-sets
"""
module ACSetViews
export ACSetView, backing, @compute_prop, @select_where,
  ACSetViewMorphism

using ...CSetDataStructures
using ..CSets
using ..FinSets
using ...Theories: CatDesc, AttrDesc, dom, codom_num

using MLStyle: @match
using PrettyTables: pretty_table
using Tables

# Data types
############

""" One useful way of thinking of ACSetViews is that they are
dataframes backed by ACSets.

They are also useful for computing limits

P is a symbol, the object in the schema that the view focuses onto
"""
struct ACSetView{A <: ACSet, P, Attrs <: NamedTuple} <: AbstractArray{Attrs,1}
  backing :: A
  parts :: AbstractVector{Int}
  function ACSetView(backing::A, P::Symbol, parts::AbstractVector{Int}) where
      {CD <: CatDesc,AD <: AttrDesc{CD},Ts <: Tuple, A <: ACSet{CD,AD,Ts}}
    attr_names = filter(a -> dom(AD,a) == P, AD.attr)
    attr_types = map(a -> Ts.parameters[codom_num(AD,a)], attr_names)
    row_type = NamedTuple{attr_names, Tuple{attr_types...}}
    new{A,P,row_type}(backing,parts)
  end
  function ACSetView(backing::A, P::Symbol) where
      {CD <: CatDesc,AD <: AttrDesc{CD},Ts <: Tuple, A <: ACSet{CD,AD,Ts}}
    ACSetView(backing,P,1:nparts(backing,P))
  end
end

# We're going to overload getproperty, so we use these instead

backing(av::ACSetView) = getfield(av,:backing)

CSetDataStructures.parts(av::ACSetView) = getfield(av,:parts)

# Overloaded accessors
######################

Base.size(av::ACSetView) = size(parts(av))

@generated function Base.getindex(av::ACSetView{<:ACSet,P,Attrs}, i::Int) where
    {P, Attrs <: NamedTuple}
  tuple_args = map(Attrs.parameters[1]) do attr
    :($attr = av[i,$(Expr(:quote, attr))])
  end
  Expr(:tuple,tuple_args...)
end

function Base.getindex(av::ACSetView{<:ACSet,P,<:NamedTuple},i::Int,name::Union{Symbol,Vector{Symbol}}) where {P}
  subpart(backing(av), parts(av)[i], name)
end

function Base.getindex(av::ACSetView{<:ACSet,P,<:NamedTuple},i,name::Union{Symbol,Vector{Symbol}}) where {P}
  subpart(backing(av), view(parts(av),i), name)
end

function Base.getindex(av::ACSetView{<:ACSet,P,<:NamedTuple},name::Union{Symbol,Vector{Symbol}}) where {P}
  subpart(backing(av), parts(av), name)
end

function Base.getindex(av::ACSetView{<:ACSet,P,<:NamedTuple},indicator::BitArray{1}) where {P}
  ACSetView(backing(av),P,parts(av)[indicator])
end

Base.getproperty(av::ACSetView, name::Symbol) =
  subpart(backing(av), parts(av), name)

function Base.show(io::IO, ::MIME"text/plain", av::T) where {ACS <: AbstractACSet, P, T<:ACSetView{ACS,P}}
  print(io, "ACSetView of object $P")
  println(io)
  table = tables(backing(av))[P][parts(av)]
  if !(eltype(table) <: NamedTuple{(),Tuple{}} || isempty(table))
    pretty_table(io, table, nosubheader=true,
                 show_row_number=true, row_number_column_title=string(P))
  end
end

# Tables.jl interface
#####################

Tables.istable(av::ACSetView) = true

Tables.columnaccess(av::ACSetView) = true
Tables.columns(av::ACSetView) = av

function Tables.getcolumn(av::ACSetView{<:ACSet,P,Attrs}, i::Int) where {P, Attrs <: NamedTuple}
  av[:,Attrs.parameters[1][i]]
end

Tables.getcolumn(av::ACSetView, nm::Symbol) = av[:,nm]

Tables.getcolumn(av::ACSetView, i::Int, nm::Symbol) = av[i,nm]

Tables.columnnames(av::ACSetView{<:ACSet,P,Attrs}) where {P,Attrs<:NamedTuple} = Attrs.parameters[1]

# Macros for Querying/filtering
###############################

""" This macro takes an expression using attributes of the object selected by the ACSetView,
and returns an array of the result of that expression applied to the attribute of each part
in the ACSetView.

Example: @compute_prop(edges_table, 2 * src.weight)

Eventually, this should be superceded by Query.jl or AlgebraicRelations.jl functionality, so I have
purposely made the name a bit obtuse and awkward.

However, currently these have functionality not replicated by Query.jl: the ability to use compound attributes.
TODO: Add this functionality to Query.jl
"""
macro compute_prop(av, ex)
  :(@. $(replace_symbols(esc(av), ex)))
end

replace_symbols(av, ex::Symbol) = :($av[$(Expr(:quote, ex))])

function replace_symbols(av, ex::Expr)
  @match ex begin
    Expr(:call, f, args...) => Expr(:call, esc(f), map(arg -> replace_symbols(av,arg), args)...)
    Expr(:$, arg) => esc(arg)
    Expr(:., args...) => :($av[$(flatten_dot_expr(ex))])
    _ => error("unsupported syntax in @query")
  end
end

replace_symbols(av, ex) = ex

function flatten_dot_expr(ex::Expr)
  @match ex begin
    Expr(:., x, y) => [flatten_dot_expr(x)..., y.value]
    _ => error("unsupported syntax in @query")
  end
end

function flatten_dot_expr(ex::Symbol)
  [ex]
end

""" This macro takes in a similar expression to @compute_prop, except that the expression
should return a boolean, and then it returns a new ACSetView containing only the parts
that return true when the expression is applied to them.

Again, eventually this should be superceded by Query.jl or AlgebaricRelations.jl, but again
there is current functionality that is not available in either.
"""
macro select_where(av, ex)
  :($(esc(av))[@. $(replace_symbols(esc(av), ex))])
end

# Morphisms
###########

struct ACSetViewMorphism{AV <: ACSetView}
  f::FinFunction{Int,Int}
  dom::AV
  codom::AV
  function ACSetViewMorphism{AV}(f::FinFunction, dom::AV, codom::AV) where {AV <: ACSetView}
    @assert length(dom(f)) == size(dom)
    @assert length(codom(f)) == size(codom)
    new{AV}(f,dom,codom)
  end
  function ACSetViewMorphism(α::ACSetTransformation,P::Symbol)
    dom = ACSetView(α.dom,P)
    codom = ACSetView(α.codom,P)
    f = α.components[P]
    new{typeof(dom)}(f,dom,codom)
  end
end

function CSets.is_natural(α::ACSetViewMorphism)
  all(α.dom[i] == α.codom[α.f(i)] for i in eachindex(α.dom))
end

end


text:
using Oceananigans.Utils: cell_advection_timescale
using Oceananigans.TurbulenceClosures: cell_diffusion_timescale

"""
    CFL{D, S}

An object for computing the Courant-Freidrichs-Lewy (CFL) number.
"""
struct CFL{D, S}
           Δt :: D
    timescale :: S
end

"""
    CFL(Δt [, timescale=Oceananigans.cell_advection_timescale])

Returns an object for computing the Courant-Freidrichs-Lewy (CFL) number
associated with time step or `TimeStepWizard` `Δt` and `timescale`.

See also `AdvectiveCFL` and `DiffusiveCFL`.
"""
CFL(Δt) = CFL(Δt, cell_advection_timescale)

(c::CFL)(model) = c.Δt / c.timescale(model)

"""
    AdvectiveCFL(Δt)

Returns an object for computing the Courant-Freidrichs-Lewy (CFL) number
associated with time step or `TimeStepWizard` `Δt` and the time scale
for advection across a cell.

Example
=======
```julia
julia> model = NonhydrostaticModel(grid=RegularRectilinearGrid(size=(16, 16, 16), length=(8, 8, 8)));

julia> cfl = AdvectiveCFL(1.0);

julia> data(model.velocities.u) .= π;

julia> cfl(model)
6.283185307179586
```
"""
AdvectiveCFL(Δt) = CFL(Δt, cell_advection_timescale)

"""
    DiffusiveCFL(Δt)

Returns an object for computing the diffusive Courant-Freidrichs-Lewy (CFL) number
associated with time step or `TimeStepWizard` `Δt` and the time scale for diffusion
across a cell associated with `model.closure`.

The maximum diffusive CFL number among viscosity and all tracer diffusivities is
returned.

Example
=======
```julia
julia> model = NonhydrostaticModel(grid=RegularRectilinearGrid(size=(16, 16, 16), length=(1, 1, 1)));

julia> dcfl = DiffusiveCFL(0.1);

julia> dcfl(model)
2.688e-5
```
"""
DiffusiveCFL(Δt) = CFL(Δt, cell_diffusion_timescale)

#####
##### Accurate CFL via reduction
#####

using CUDA, CUDAKernels, KernelAbstractions, Tullio

using Oceananigans.Models
using Oceananigans.Grids: halo_size
using Oceananigans.Operators: Δxᶠᶜᵃ, Δyᶜᶠᵃ, Δzᵃᵃᶠ

accurate_cell_advection_timescale(model) = accurate_cell_advection_timescale(model.grid, model.velocities)

function accurate_cell_advection_timescale(grid, velocities)
    Nx, Ny, Nz = size(grid)
    Hx, Hy, Hz = halo_size(grid)

    is = 1+Hx:Nx+Hx
    js = 1+Hy:Ny+Hy
    ks = 1+Hz:Nz+Hz

    u = view(velocities.u.data.parent, is, js, ks)
    v = view(velocities.v.data.parent, is, js, ks)
    w = view(velocities.w.data.parent, is, js, ks)

    min_timescale = minimum(
        @tullio (min) timescale[k] := 1 / (  abs(u[i, j, k]) / Δxᶠᶜᵃ(i, j, k, grid)
                                           + abs(v[i, j, k]) / Δyᶜᶠᵃ(i, j, k, grid)
                                           + abs(w[i, j, k]) / Δzᵃᵃᶠ(i, j, k, grid))
    )

    return min_timescale
end


text:
# generate test plots of scores
#    matching https://en.wikipedia.org/wiki/Scoring_rule
#
#
using RatPack
using PyPlot

x = collect( 0 : 0.01 : 1.0 )
norm = false
norm = true

rule_log = ScoreLogarithmic(; normalise=norm )
y_log = zeros(size(x))
for i=1:length(x)
    y_log[i] = scoring_function( rule_log, [x[i], 1-x[i]], 1 )
end

rule_brier = ScoreBrier(; normalise=norm )
y_brier = zeros(size(x))
for i=1:length(x)
    y_brier[i] = scoring_function( rule_brier, [x[i], 1-x[i]], 1 )
end

rule_quad = ScoreQuadratic(; normalise=norm )
y_quad = zeros(size(x))
for i=1:length(x)
    y_quad[i] = scoring_function( rule_quad, [x[i], 1-x[i]], 1 )
end

rule_sphere = ScoreSpherical(; normalise=norm )
y_sphere = zeros(size(x))
for i=1:length(x)
    y_sphere[i] = scoring_function( rule_sphere, [x[i], 1-x[i]], 1 )
end

figure(1)
clf()
plot(x, y_log)
plot(x, y_brier)
plot(x, y_quad)
plot(x, y_sphere)
ylim([-3, 1])

# test scoring on cases


text:
using DashBootstrapComponents, DashHtmlComponents

badges = html_div([
    html_h1(["Example heading", dbc_badge("New", className = "ms-1")]),
    html_h2(["Example heading", dbc_badge("New", className = "ms-1")]),
    html_h3(["Example heading", dbc_badge("New", className = "ms-1")]),
    html_h4(["Example heading", dbc_badge("New", className = "ms-1")]),
    html_h5(["Example heading", dbc_badge("New", className = "ms-1")]),
    html_h6(["Example heading", dbc_badge("New", className = "ms-1")]),
]);

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
#=
SubString tests

Copyright 2018 Gandalf Software, Inc., Scott P. Jones, and other contributors to the Julia language
Licensed under MIT License, see LICENSE.md
Based initially on julia/test/strings/types.jl
=#

function testuni(T)
    str = T(u8str)
    str2 = str^2
    str2plain = u8str^2

    @test length(str2) == 2 * length(str)
    @test sizeof(str2) == 2 * sizeof(str)

    for i1 = 1:length(str2)
        if is_valid(str2, i1)
            for i2 = i1:length(str2)
                if is_valid(str2, i2)
                    @eval @test length($str2[$i1:$i2]) == length($str2plain[$i1:$i2])
                    @eval @test $str2[$i1:$i2] == $str2plain[$i1:$i2]
                end
            end
        end
    end

    # tests that SubString of a single multibyte `Char` string, like "∀" which takes 3 bytes
    # gives the same result as `getindex` (except that it is a view not a copy)

    s1 = T("∀")
    s2 = T("∀∀")

    for idx in 0:1
        @test SubString(s1, 1, idx) == s1[1:idx]
    end

    # Substring provided with invalid end index throws BoundsError
    @test_throws StringIndexError  SubString(s1, 1, 2)
    @test_throws StringIndexError  SubString(s1, 1, 3)
    @test_throws BoundsError SubString(s1, 1, 4)

    # Substring provided with invalid start index throws BoundsError
    @test SubString(s2, 1, 1) == s1
    @test SubString(s2, 1, 4) == s2
    @test SubString(s2, 4, 4) == s1

    @test_throws StringIndexError  SubString(s2, 1, 2)
    @test_throws StringIndexError  SubString(s2, 1, 5)
    @test_throws StringIndexError  SubString(s2, 2, 4)
    @test_throws BoundsError SubString(s2, 0, 1)
    @test_throws BoundsError SubString(s2, 0, 4)
    @test_throws BoundsError SubString(s2, 1, 7)
    @test_throws BoundsError SubString(s2, 4, 7)

    # tests for SubString of more than one multibyte `Char` string
    # we are consistent with `getindex` for `String`
    for idx in [0, 1, 4]
        @test SubString(s2, 1, idx) == s2[1:idx]
        @test SubString(s2, 4, idx) == s2[4:idx]
    end

    # index beyond lastindex(s2)
    for idx in [2:3; 5:6]
        @test_throws StringIndexError SubString(s2, 1, idx)
    end
    for idx in 7:8
        @test_throws BoundsError SubString(s2, 1, idx)
    end

    str = T("aa\u2200\u2222bb")
    u = SubString(str, 3, 6)
    @test length(u) == 2
    b = IOBuffer()
    write(b, u)
    @test String(take!(b)) == "\u2200\u2222"

    @test_throws StringIndexError SubString(str, 4, 5)
    @test_throws BoundsError iterate(u, 0)
    @test iterate(u, 7) === nothing
    @test_throws BoundsError getindex(u, 0)
    @test_throws BoundsError getindex(u, 7)
    @test_throws BoundsError getindex(u, 0:1)
    @test_throws BoundsError getindex(u, 7:7)
    @test reverseind(u, 1) == 4
    if T === String
        @test typeof(Base.cconvert(Ptr{Int8}, u)) == SubString{String}
    else
        @test typeof(Base.cconvert(Ptr{Int8}, u)) == String
    end
    @test Base.cconvert(Ptr{Int8}, u) == u

    str = T("føøbar")
    @test_throws BoundsError SubString(str, 10, 10)
    u = SubString(str, 4, 3)
    @test length(u) == 0
    b = IOBuffer()
    write(b, u)
    @test String(take!(b)) == ""

    # sizeof
    @test sizeof(SubString(T("abc\u2222def"), 4, 4)) == 3

    # proper nextind/prevind/thisind for SubString{String}
    let rng = MersenneTwister(1),
        strs = ["∀∃∀" * String(rand(rng, Char, 40)) * "∀∃∀", String(rand(rng, Char, 50))]
        for sa in strs
            try
                s = T(sa)
            catch ex
                println("$T(\"$sa\") = $(collect(codeunits(sa)))")
                s = sa
            end
            a = 0
            while !str_done(s, a)
                a = nextind(s, a)
                b = a - 1
                while !str_done(s, b)
                    ss = SubString(s, a, b)
                    s2 = s[a:b]
                    @test ncodeunits(ss) == ncodeunits(s2)
                    for i in 1:ncodeunits(ss)
                        @test thisind(ss, i) == thisind(s2, i)
                    end
                    for i in 0:ncodeunits(ss)
                        @test nextind(ss, i) == nextind(s2, i)
                        for j in 0:ncodeunits(ss)+5
                            if j > 0 || is_valid(ss, i)
                                @test nextind(ss, i, j) == nextind(s2, i, j)
                            end
                        end
                    end
                    for i in 1:ncodeunits(ss)+1
                        @test prevind(ss, i) == prevind(s2, i)
                        for j in 0:ncodeunits(ss)+5
                            if j > 0 || is_valid(ss, i)
                                @test prevind(ss, i, j) == prevind(s2, i, j)
                            end
                        end
                    end
                    b = nextind(s, b)
                end
            end
        end
    end

    # for is_valid(SubString{String})
    s = T("Σx + βz - 2")
    for i in -1:ncodeunits(s)+2
        if checkbounds(Bool, s, i)
            if is_valid(s, i)
                ss = SubString(s, 1, i)
                for j = 1:ncodeunits(ss)
                    @test is_valid(ss, j) == is_valid(s, j)
                end
            else
                @test_throws StringIndexError SubString(s, 1, i)
            end
        elseif i > 0
            @test_throws BoundsError SubString(s, 1, i)
        else
            @test SubString(s, 1, i) == ""
        end
    end

    # length(SubString{String}) performance specialization
    s = T("|η(α)-ϕ(κ)| < ε")
    @test length(SubString(s, 1, 0))  == length(s[1:0])
    @test length(SubString(s, 4, 4))  == length(s[4:4])
    @test length(SubString(s, 1, 7))  == length(s[1:7])
    @test length(SubString(s, 4, 11)) == length(s[4:11])

    if false
    @testset "reverseind" for S in (T, SubString, GenericString)
        for prefix in ("", "abcd", "\U0001d6a4\U0001d4c1", "\U0001d6a4\U0001d4c1c",
                       " \U0001d6a4\U0001d4c1"),
            suffix in ("", "abcde", "\U0001d4c1β\U0001d6a4", "\U0001d4c1β\U0001d6a4c",
                       " \U0001d4c1β\U0001d6a4"),
            c in ('X', 'δ', '\U0001d6a5')

            s = convert(S, T(string(prefix, c, suffix)))
            r = reverse(s)
            ri = fnd(First, ==(c), r)
            @test c == s[reverseind(s, ri)] == r[ri]
            s = convert(S, T(string(prefix, prefix, c, suffix, suffix)))
            pre = convert(S, T(prefix))
            sb = SubString(s, nextind(pre, lastindex(pre)),
                           lastindex(convert(S, T(string(prefix, prefix, c, suffix)))))
            r = reverse(sb)
            ri = fnd(First, ==(c), r)
            @test c == sb[reverseind(sb, ri)] == r[ri]
        end
    end
    end
end

function teststr(T)
    str = T("tempus fugit")              #length(str)==12

    ss = SubString(str, 1, lastindex(str)) #match source string
    @test length(ss) == length(str)

    ss = SubString(str, 1:lastindex(str))
    @test length(ss) == length(str)

    ss = SubString(str,1,0)    #empty SubString
    @test length(ss)==0

    ss = SubString(str,1:0)
    @test length(ss)==0

    @test_throws BoundsError SubString(str, 14, 20)  #start indexing beyond source string length
    @test_throws BoundsError SubString(str, 10, 16)  #end indexing beyond source string length

    @test_throws BoundsError SubString(T(""), 1, 4)  #empty source string
    @test_throws BoundsError SubString(T(""), 1, 1)  #empty source string, identical start and end index
    @test_throws BoundsError SubString(T(""), 10, 12)
    @test SubString(T(""), 12, 10) == ""

    @test SubString(T("foobar"), big(1), big(3)) == "foo"

    # search and SubString (issue #5679)
    let str = T("Hello, world!"),
        u = SubString(str, 1, 5)
        @test fnd(Last, T("World"), u) == 0:-1
        @test fnd(Last, ==('z'), u) == 0
        @test fnd(Last, T("ll"), u) == 3:4
    end

    # SubString created from SubString
    let str = T("Hello, world!"),
        u = SubString(str, 2, 5)
        for idx in 1:4
            @test SubString(u, 2, idx) == u[2:idx]
            @test SubString(u, 2:idx) == u[2:idx]
        end
        @test_throws BoundsError SubString(u, 1, 10)
        @test_throws BoundsError SubString(u, 1:10)
        @test_throws BoundsError SubString(u, 20:30)
        @test_throws BoundsError SubString(u, -1:10)
        @test SubString(u, 20:15) == ""
        @test SubString(u, -1, -10) == ""
        @test SubString(SubString(T("123"), 1, 2), -10, -20) == ""
    end

    # sizeof
    @test sizeof(SubString("abc\u2222def",4,4)) == 3

    # issue #3710
    @test prevind(SubString(T("{var}"),2,4),4) == 3

    # issue #4183
    @test split(SubString(T("x"), 2, 0), "y") == [""]

    # issue #6772
    if T != UTF8Str # Todo: Broken for UTF8Str?
        @test parse(Float64, SubString(T("10"),1,1)) === 1.0
        @test parse(Float64, SubString(T("1 0"),1,1)) === 1.0
        @test parse(Float32, SubString(T("10"),1,1)) === 1.0f0
    end

    # issue #5870
    @test !occurs_in(Regex(T("aa")), SubString(T(""),1,0))
    @test occurs_in(Regex(T("")), SubString(T(""),1,0))

    # is_valid, length, prevind, nextind for SubString{String}
    let s = T("lorem ipsum"),
        sdict = Dict(SubString(s, 1, 11)  => T("lorem ipsum"),
                     SubString(s, 1, 6)   => T("lorem "),
                     SubString(s, 1, 0)   => T(""),
                     SubString(s, 2, 4)   => T("ore"),
                     SubString(s, 2, 11)  => T("orem ipsum"),
                     SubString(s, 15, 14) => T(""),
                     )
        for (ss, s) in sdict
            @test ncodeunits(ss) == ncodeunits(s)
            for i in -2:13
                @test is_valid(ss, i) == is_valid(s, i)
            end
            for i in 1:ncodeunits(ss), j = i-1:ncodeunits(ss)
                @test length(ss, i, j) == length(s, i, j)
            end
        end
        for (ss, s) in sdict
            @test length(ss) == length(s)
            for i in 0:ncodeunits(ss), j = 0:length(ss)+1
                @test prevind(ss, i+1, j) == prevind(s, i+1, j)
                @test nextind(ss, i, j) == nextind(s, i, j)
            end
            @test_throws BoundsError prevind(s, 0)
            @test_throws BoundsError prevind(ss, 0)
            @test_throws BoundsError nextind(s, ncodeunits(ss)+1)
            @test_throws BoundsError nextind(ss, ncodeunits(ss)+1)
        end
    end

    ss = SubString(T("hello"), 1, 5)
    @test length(ss, 1, 0) == 0
    @test_throws BoundsError length(ss, 1, -1)
    @test_throws BoundsError length(ss, 1, 6)
    @test_throws BoundsError length(ss, 1, 10)
    #println("$(typeof(ss))(\"$ss\")")
    @test_throws BoundsError prevind(ss, 0, 1)
    @test prevind(ss, 1, 1) == 0
    @test prevind(ss, 6, 1) == 5
    @test_throws BoundsError prevind(ss, 7, 1)
    @test_throws BoundsError nextind(ss, -1, 1)
    @test nextind(ss, 0, 1) == 1
    @test nextind(ss, 5, 1) == 6
    @test_throws BoundsError nextind(ss, 6, 1)

    @testset "reverseind of empty strings" begin
        for s in ("",
                  SubString(T(""), 1, 0),
                  SubString(T("ab"), 1, 0),
                  SubString(T("ab"), 2, 1),
                  SubString(T("ab"), 3, 2),
                  GenericString(T("")))
            @test reverseind(s, 0) == 1
            @test reverseind(s, 1) == 0
        end
    end

    ## Cstring tests ##

    # issue #13974: comparison against pointers
    str = T("foobar")
    ptr = pointer(str)
    cstring = Cstring(ptr)
    @test ptr == cstring
    @test cstring == ptr

    # convenient NULL string creation from Ptr{Cvoid}
    nullstr = Cstring(C_NULL)

    # Comparisons against NULL strings
    @test ptr != nullstr
    @test nullstr != ptr

    # Short-hand comparison against C_NULL
    @test nullstr == C_NULL
    @test C_NULL == nullstr
    @test cstring != C_NULL
    @test C_NULL != cstring

end

## SubString tests ##
@testset "SubString tests" begin
    for T in (UTF8Str, ) #UCS2Str
        @testset "Test Unicode substrings $T" begin
            testuni(T)
        end
    end
    for T in (UTF8Str, ASCIIStr, LatinStr) #UCS2Str
        @testset "Test ASCII substrings $T" begin
            teststr(T)
        end
    end
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(γ::Variogram; nlags=100, maxlag=3.)
  # start at 1e-6 instead of 0 to avoid nugget artifact in plot
  h = range(1e-6, stop=maxlag, length=nlags)

  seriestype --> :path
  xlim --> (0, maxlag)
  ylim --> (0, Inf)
  xguide --> "Lag h"
  yguide --> "Variogram(h)"
  label --> "variogram"

  h, γ.(h)
end


text:
@testset ExtendedTestSet "add.jl" begin
    @testset "add.jl" begin
        @test add(4, 4) == 8
    end

    @testset "add.jl" begin
        @test add(1, 2) == 3
    end

    @testset "add.jl" begin
        @test add(-1, -1) == -2
    end
end


text:
"""
# module JuliaToFortran

# Examples

```jldoctest
julia>
```
"""
module JuliaToFortran

using Fortran90Namelists.FortranToJulia

export to_fortran

function to_fortran(v::Int)
    FortranData(string(v))
end
function to_fortran(v::Float32; scientific::Bool = false)
    str = string(v)
    scientific && return FortranData(replace(str, r"f"i => "e"))
    return FortranData(str)
end
function to_fortran(v::Float64; scientific::Bool = false)
    str = string(v)
    scientific && return FortranData(replace(str, r"e"i => "d"))
    return FortranData(string(v))
end
function to_fortran(v::Bool)
    v ? FortranData(".true.") : FortranData(".false.")
end
function to_fortran(v::AbstractString)
    return FortranData("'$v'")
end

function Base.string(s::FortranData)
    return string(s.data)
end

end


text:
using MathLink
import MathLink: put, get, endpacket, nextpacket, getnextraw

link = MathLink._defaultlink()

X = rand(80,120)

put(link, X)
endpacket(link)
pkt = nextpacket(link)
t = getnextraw(link)

import MathLink: loopbackopen, CLink, mlib, getfunction

rmeterp = Ref(C_NULL)
rleaftoken = Ref(C_NULL)
rdepth = Ref{Clong}(0)

loopback = loopbackopen()

ccall((:MLGetArrayTypeWithDepthAndLeafType, mlib), Cint,
      (CLink, CLink, Ptr{Ptr{Cvoid}}, Ptr{Clong}, Ptr{Ptr{Cvoid}}),
      link, loopback, rmeterp, rdepth, rleaftoken)

sz = ntuple(i -> getfunction(loopback)[2], rdepth[])

tk = MathLink.getrawtype(link)
Y = Array{Float64}(undef, sz)

ccall((:MLGetBinaryNumberArrayData, mlib), Cint,
      (CLink, Ptr{Cvoid}, Ptr{Cvoid}, Clong, Clong),
      link, rmeterp[], Y, sizeof(Y), tk)

@test X == Y

w = WExpr(WSymbol("Sqrt"), Any[2])
MathLink.meval(w)

u = MathLink.WReal("2.000000000000000000")
w = WExpr(WSymbol("Sqrt"), Any[u])

MathLink.put(link, WExpr(WSymbol("EvaluatePacket"), Any[w]))
MathLink.EndPacket(link)
packet, _ = MathLink.getfunction(link)
t = MathLink.GetRawType(link)
MathLink.get(link, Float64)

# Integers 228 (Int32) / 230 (Int64)
# Real 246 (Float64)
MathLink.meval(w)

using MathLink
import MathLink: WExpr, WSymbol

link = MathLink._defaultlink()

MathLink.put(link, WExpr(WSymbol("EvaluatePacket"), Any[WExpr(WSymbol("ToExpression"), Any["Function[x,x*2][3]", WSymbol("InputForm"), WSymbol("Hold")])]))
MathLink.EndPacket(link)
packet, _ = MathLink.getfunction(link)

t = MathLink.GetNextRaw(link)
MathLink.getargcount(link)
t = MathLink.GetNextRaw(link)
MathLink.get(link, WSymbol)

t = MathLink.GetNext(link)
MathLink.getargcount(link)
MathLink.get(link, WSymbol)

MathLink.getfunction(link)

t = MathLink.GetRawType(link)
MathLink.getargcount(link)


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSVCProjPO

using Test, RLinearAlgebra, LinearAlgebra, Random

@testset "LSVC Projection Partial Orthogonal -- Procedural" begin

    # Supertype
    @test supertype(LinSysVecColProjPO) == LinSysVecColProjection

    # Testing context
    Random.seed!(1010)

    A = rand(5, 10)
    x = rand(10)
    b = A * x

    z = rand(10)

    rsub = LinSysVecColProjPO()

    # Verify that Z (storage) is nothing
    @test isnothing(rsub.Z)

    # Verify that Z (storage) is appropriately initialized
    @test let
        e = zeros(10)
        e[1] = 1.0
        res = dot(A[:,1], A * z - b)
        RLinearAlgebra.rsubsolve!(rsub, z, (e, A, res), 1)
        (length(rsub.Z) == rsub.m) & reduce(&, [length(z)==10 for z in rsub.Z])
    end

    # Verify that solution has zero residual against previous directions
    @test let
        for i = 1:6
            e = zeros(10)
            e[i] = 1.0
            res = dot(A[:,i], A * z - b)
            RLinearAlgebra.rsubsolve!(rsub, z, (e, A, res), i)
        end

        norm( A[:,2:6]'*(A * z - b)) < 1e-13 #Ajdust for condition number
    end

    @test let
        for i = 1:8
            e = zeros(10)
            e[i] = 1.0
            res = dot(A[:,i], A * z - b)
            RLinearAlgebra.rsubsolve!(rsub, z, (e, A, res), i)
        end

        norm( A[:,4:8]'*(A * z - b)) < 1e-13 #Adjust for condition number
    end
end

end # End module


text:
using KernelAbstractions
using Statistics
using Oceananigans.Architectures: architecture, device_event
using Oceananigans.Fields: location, ZReducedField, Field

instantiate(X) = X()

mask_immersed_field!(field::Field, value=zero(eltype(field.grid))) =
    mask_immersed_field!(field, field.grid, location(field), value)

mask_immersed_field!(field, grid, loc, value) = NoneEvent()

function mask_immersed_field!(field::Field, grid::ImmersedBoundaryGrid, loc, value)
    arch = architecture(field)
    loc = instantiate.(loc)
    return launch!(arch, grid, :xyz, _mask_immersed_field!, field, loc, grid, value; dependencies = device_event(arch))
end

@kernel function _mask_immersed_field!(field, loc, grid, value)
    i, j, k = @index(Global, NTuple)
    @inbounds field[i, j, k] = scalar_mask(i, j, k, grid, grid.immersed_boundary, loc..., value, field)
end

mask_immersed_reduced_field_xy!(field::ZReducedField, value=zero(eltype(field.grid)); k) =
    mask_immersed_reduced_field_xy!(field, field.grid, location(field), value; k)

mask_immersed_reduced_field_xy!(::Nothing, args...; kwargs...) = NoneEvent()

mask_immersed_reduced_field_xy!(field, grid, loc, value; k) = NoneEvent()

function mask_immersed_reduced_field_xy!(field::ZReducedField, grid::ImmersedBoundaryGrid, loc, value; k)
    arch = architecture(field)
    loc = instantiate.(loc)
    return launch!(arch, grid, :xy,
                   _mask_immersed_reduced_field_xy!, field, loc, grid, value, k;
                   dependencies = device_event(arch))
end

@kernel function _mask_immersed_reduced_field_xy!(field, loc, grid, value, k)
    i, j = @index(Global, NTuple)
    @inbounds field[i, j, 1] = scalar_mask(i, j, k, grid, grid.immersed_boundary, loc..., value, field)
end

#####
##### mask_immersed_velocities for NonhydrostaticModel
#####

mask_immersed_velocities!(U, arch, grid) = tuple(NoneEvent())

#####
##### Masking for GridFittedBoundary
#####

@inline function scalar_mask(i, j, k, grid, ::AbstractGridFittedBoundary, LX, LY, LZ, value, field)
    return @inbounds ifelse(solid_interface(LX, LY, LZ, i, j, k, grid),
                            value,
                            field[i, j, k])
end

mask_immersed_velocities!(U, arch, grid::ImmersedBoundaryGrid) = Tuple(mask_immersed_field!(q) for q in U)


text:
# Check conversion and output
for i = 1:10000
    a = {rand(Int8), rand(Int16), rand(Int32), rand(Int64), rand(Int128),
         rand(Uint8), rand(Uint16), rand(Uint32), rand(Uint64), rand(Uint128),
         rand(Float32), rand(Float64)}
    b = map(string, a)
    for j in a
        if j != 0
            @test "[$(BigFloat(j)), $(BigFloat(j))]" == string(Interval(j))
        else
            @test "[$(BigFloat(j)), -$(BigFloat(j))]" == string(Interval(j))
        end
    end
end

# Check precision
a = Interval(1)
@test precision(a) == get_bigfloat_precision()

# Check squaring
a = Interval(-1, 1)
@test (left(square(a)), right(square(a))) == (0, 1)
@test (left(a * a), right(a * a)) == (-1, 1)

# Check -
a = Interval(0, 1)
@test (left(-a), right(-a)) == (-1, 0)
@test signbit(right(-a)) == 1

# Check rand
a = Interval(-1, 1)
@test -1 <= rand(a) <= 1
r = rand(a, 1, 3)
@test size(r) == (1, 3)
@test typeof(r) == Matrix{BigFloat}


text:
struct ColorizedArray{C<:Colorant,N,A<:AbstractArray,L<:AbstractArray} <: AbstractArray{C,N}
    intensity::A
    label::L
end

"""
    ColorizedArray(intensity, label::IndirectArray) -> A

Create an array, combining a `label` array (where each pixel is
assigned one of a list of discrete colors) and an `intensity` array
(where each pixel has a scalar value). `A` satisfies

    A[i,j,...] = intensity[i,j,...] * label[i,j,...]

The label array "tinges" the grayscale intensity with the color
associated with that point's label.

This computation is performed lazily, as to be suitable even for large arrays.
"""
function ColorizedArray(intensity, label::IndirectArray{C,N}) where {C<:Colorant,N}
    depwarn("`ColorizedArray(intensity, label)` is deprecated, use `mappedarray(*, intensity, label)` from `MappedArrays` instead", :ColorizedArray)

    axes(intensity) == axes(label) || throw(DimensionMismatch("intensity and label must have the same axes, got $(axes(intensity)) and $(axes(label))"))
    CI = typeof(zero(C)*zero(eltype(intensity)))
    ColorizedArray{CI,N,typeof(intensity),typeof(label)}(intensity, label)
end

# TODO: an implementation involving AxisArray that matches the shared
# axes, and therefore allows `label` to be of lower dimensionality
# than `intensity`.

intensitytype(::Type{ColorizedArray{C,N,A,L}}) where {C<:Colorant,N,A<:AbstractArray,L<:AbstractArray} = A
labeltype(::Type{ColorizedArray{C,N,A,L}}) where {C<:Colorant,N,A<:AbstractArray,L<:AbstractArray} = L

Base.size(A::ColorizedArray) = size(A.intensity)
Base.axes(A::ColorizedArray) = axes(A.intensity)
Base.IndexStyle(::Type{CA}) where {CA<:ColorizedArray} = IndexStyle(IndexStyle(intensitytype(CA)), IndexStyle(labeltype(CA)))

@inline function Base.getindex(A::ColorizedArray, i::Integer)
    @boundscheck checkbounds(A, i)
    @inbounds ret = A.intensity[i]*A.label[i]
    ret
end
@inline function Base.getindex(A::ColorizedArray{C,N}, I::Vararg{Int,N}) where {C,N}
    @boundscheck checkbounds(A, I...)
    @inbounds ret = A.intensity[I...]*A.label[I...]
    ret
end


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for the Combinatorial API

############################################################################################################################

@testset "Combinatorial " begin
   g = completegraph(10)

   @test nv(g) == 10
   @test ne(g) == 90
   @test size(g) == (10,90)

   @test vertices(g) == 1 : 10

   @test hasvertex(g, rand(1:10))
   @test all(hasvertex(g, rand(1:10, 4)))

   @test hasedge(g, 1=>2)
   @test all(hasedge(g, edges(g)))

   @test fadj(g, 1) == collect(2:10)
   @test fadj!(g, 10, zeros(Int, 9)) == collect(1:9)

   @test outdegree(g, rand(1:10)) == 9
   @test outdegree(g, [1,2,3]) == [9,9,9]
   @test outdegree(g) == fill(9, 10)

   @test indegree(g, rand(1:10)) == 9
   @test indegree(g, [1,2,3]) == [9,9,9]
   @test indegree(g) == fill(9, 10)

   v = rand(vertices(g))
   @test g[v] == fadj(g, v)
end


text:
include("setup.jl")

@testset "BasisFunctions" begin
    include("basisfunctions.jl")
end

@testset "Fitting" begin
    include("fit.jl")
end

@testset "Designmatrix" begin
    include("designmatrix.jl")
end

@testset "Splines" begin
    include("splines.jl")
end

@testset "Predict" begin
    include("predict.jl")
end

@testset "Effects" begin
    include("effects.jl")
end

@testset "Statistics" begin
    include("statistics.jl")
end

@testset "Utilities" begin
    include("utilities.jl")
end
#@testset "ClusterPermutation" begin
#include("clusterpermutation.jl")
#end

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
function _transform(x::Array{nf_elem}, y::fmpz_mat)
  n = length(x)
  @assert n == nrows(y)
  m = ncols(y)
  z = Array{nf_elem}(m)
  for i in 1:m
    z[i] = x[1]^y[1, i]
    for j in 2:n
      z[i] = z[i]*x[j]^y[j, i]
    end
  end
  return z
end

function _make_row_primitive(x::fmpz_mat, j::Int)
  y = x[j, 1]
  for i in 1:ncols(x)
    y = gcd(y, x[j, i])
  end
  if y > 1
    for i in 1:ncols(x)
      x[j, i] = div(x[j, i], y)
    end
  end
end

function _make_row_primitive!(x::Array{fmpz, 1})
  y = x[1]
  for i in 2:length(x)
    y = gcd(y, x[i])
    if y == 1
      return x
    end
  end
  if y > 1
    for i in 1:ncols(x)
      x[i] = div(x[i], y)
    end
  end
end

function _primitive_element(F::Union{FqNmodFiniteField, FqFiniteField})
  #println("Computing primitive element of $F")
  #println("Have to factor $(order(F) - 1)")
  oF = order(F) - 1
  fac = factor(oF)
  exps = Vector{fmpz}()
  for (l, ll) in fac
    push!(exps, div(oF, l))
  end

  while true
    a = rand(F)
    if iszero(a)
      continue
    end
    is_primitive = true
    for e in exps
      if isone(a^e)
        is_primitive = false
        break
      end
    end
    if is_primitive
      return a
    end
  end
end


text:
using LSMH: filter_collisions!

minhash = LSMH.minhash

@testset "hashing.jl" begin
    @testset "minhash" begin
        s = "Hello, world!"

        @test minhash(s, 5) == 0x12da77c8
        @test minhash(s, 5, 0) == 0x12da77c8
        @test minhash(s, 5, 1) == 0x2011211e

        @test_throws DomainError minhash(s, 1)
        @test_throws DomainError minhash(s, 2, -1)
    end
    @testset "signature" begin
        s = "Hello, world!"

        @test signature(s, 5, 3) == UInt32[0x2011211e,0x1134f986,0x291d3c45]
        @test signature(s, 2, 2) == UInt32[0x09cf7d74,0x0a6f416e]

        @test_throws DomainError signature(s, 1, 3)
        @test_throws DomainError signature(s, 2, 0)
    end
    @testset "AbstractSignature" begin
        sig = Signature("ABC", UInt32[1,2,3])

        @test identifier(sig) == "ABC"
        @test signature(sig) == UInt32[1,2,3]

        @testset "Signature" begin
            @test isdefined(sig, :id)
            @test isdefined(sig, :signature)
        end
    end
    @testset "hashband" begin
        sig = Signature("ABC", Array{UInt32}(1:100))
        @test hashband(sig, 1, 10) == 0x83af1d94
        @test hashband(sig, 11, 20) == 0x531e9145
    end
    @testset "HashTable" begin
        ht = HashTable(UInt32, Int, 1)
        @test isdefined(ht, :band)
        @test isdefined(ht, :hashtable)
        @test band(ht) == 1
        ht[UInt32(1)] = 100
        @test hashtable(ht) == Dict{UInt32, Set{Int}}(1 => Set([100]))
        ht[UInt32(1)] = 100
        @test hashtable(ht) == Dict{UInt32, Set{Int}}(1 => Set([100]))
        ht[UInt32(1)] = 101
        @test hashtable(ht) == Dict{UInt32, Set{Int}}(1 => Set([100,101]))
        @test length(ht) == 1
        ht[UInt32(2)] = 101
        @test length(ht) == 2
        @test haskey(ht, UInt32(1))
    end
    @testset "lsh" begin
        A = Signature("A", UInt32[1,2,3,4])
        B = Signature("B", UInt32[1,2,5,6])
        @test lsh(A, 2) == UInt32[0xc3642e86, 0xe68b97c2]
        @test lsh(B, 2) == UInt32[0xc3642e86, 0x584cfdb6]

        C = Signature("C", UInt32[2,3,5,6])
        signatures = [A, B, C]
        candiate_duplicates = lsh(signatures, 2)

        candiate_duplicates[1].band
        @test collect(values(candiate_duplicates[1]))[1] == Set(["A", "B"])
        @test collect(values(candiate_duplicates[2]))[1] == Set(["B", "C"])
    end
    @testset "filter_collisions!" begin
        d = Dict("A" => Set([1,2]), "B" => Set([3]))
        ht = HashTable(1, d)
        filter_collisions!(ht)
        @test hashtable(ht) == Dict("A" => Set([1,2]))
        filter_collisions!(d)
        @test d == Dict("A" => Set([1,2]))
    end
    @testset "jaccard" begin
        @test jaccard(1:100, 2:100) == .99
        @test jaccard([], []) == 1.
    end
    @testset "candidates" begin
        d = Dict(nothing => Set([1,2,3]))
        @test candidates(d) == Vector{Set{Int}}([Set([2,3]),Set([2,1]),Set([3,1])])
    end
end


text:
module_path = splitdir(@__FILE__)[1]

# example 1 - single turbine still image

# This file will create VTK files that can be uploaded into Paraview to visualize
# a single wind turbine with no flow field and no animations.

using JLD
using LinearAlgebra
# Load WFVisual module
import WFVisual
using Statistics
wfv=WFVisual

# Load GeometricTools: https://github.com/byuflowlab/GeometricTools.jl
import GeometricTools
gt=GeometricTools

data_path = joinpath(module_path, "../datav07/")

# Create save path where to store vtks
save_path = "example_1_still/"         # Save path of this example
gt.create_path(save_path, true)

# --------------------- WIND FARM LAYOUT ---------------------------------------

turbine_x = [0.0]
turbine_y = [0.0]
turbine_z = [0.0]

nturbs = length(turbine_x)
hub_height = zeros(nturbs) .+ 100.0
rotor_diameter = zeros(nturbs) .+ 100.0
nBlades = zeros(Int8,nturbs) .+ 3

wind_direction = 228.0

yaw = [0.0] .- wind_direction

# --------------------- PERIMETER AND FLUID DOMAIN -----------------------------
NDIVSx = 100              # Cells in the parametric x-direction
NDIVSy = 100            # Cells in the parametric y-direction
NDIVSz = 25               # Cells in the geometric z-direction

# Dummy perimeter
Rper = 1000.0
perimeter_points = Rper.*[ [cos(a), sin(a), 0] for a in range(0, 2*pi, 179)]

# --------------------- GENERATE FLOW FIELD ------------------------------------

_zmin = 0
_zmax = maximum(hub_height) + 1.25*maximum(rotor_diameter)/2
z_off = 0.0
spl_s=0.01
fdom = wfv.generate_perimetergrid(perimeter_points,
                                  NDIVSx, NDIVSy, NDIVSz;
                                  z_min=_zmin+z_off, z_max=_zmax+z_off,
                                  verify_spline=false,
                                  spl_s=spl_s,
                                  save_path=nothing
                                )

# Guassian wake function
magVinf = 12.0                          # (m/s) wind speed
Vinf_Oaxis = gt.rotation_matrix2(0, 0, wind_direction) # Free-stream coordinate system
invVinf_Oaxis = Vinf_Oaxis'             # Inverse transformation
refH = mean(hub_height)                 # (m) reference height for shear layer

function wake(X)

  k = 0.0325 # we can let this vary
  CT = 8.0/9.0
  loss = 0
  alpha = 0.12

  # Calculate loss in freestream direction
  for i = 1:length(turbine_x)     # Iterate over turbines

      # Distance from this turbine in global coordinates
      DX = X - [turbine_x[i], turbine_y[i], turbine_z[i]+hub_height[i]]

      # Distance in free-stream coordinates
      dX = invVinf_Oaxis*DX
      dx = dX[1]

      if dx > 0.
          dy = dX[2]
          dz = dX[3]
          sigma = k*dx + rotor_diameter[i]/sqrt(8.0)
          this_loss = (1.0-sqrt(1.0-CT/(8*sigma^2/(rotor_diameter[i]^2))))*
                          exp(-1.0/2.0*(dy/sigma)^2)*exp(-1.0/2.0*(dz/sigma)^2)
          loss += this_loss^2
      end
  end
  loss = sqrt(loss)

  windspeed= magVinf
  # Uncomment this line adds the shear layer
  windspeed *= (X[3]/refH)^alpha

  # Return velocity vector in global coordinate system
  return Vinf_Oaxis*[ windspeed*(1-loss), 0, 0 ]
end

gt.calculate_field(fdom, wake, "wake", "vector", "node")

# --------------------- GENERATE WIND FARM -------------------------------------
wfv.generate_windfarm(rotor_diameter, hub_height, nBlades,
                              turbine_x, turbine_y, turbine_z,
                              yaw,
                              perimeter_points, fdom;
                              NDIVSx=NDIVSx, NDIVSy=NDIVSy, NDIVSz=NDIVSz,
                              save_path=save_path, spl_s=0.01,
                              data_path=data_path, paraview=false);

nothing


text:
"""
    refmesh_to_mesh!(node)
    refmesh_to_mesh(node)

Compute a node mesh based on the reference mesh, the transformation matrix and the tapering.
The mutating version adds the new mesh to the `mesh` field of the geometry attribute of the
node.

# Examples

```julia
using PlantGeom
file = joinpath(dirname(dirname(pathof(PlantGeom))),"test","files","simple_plant.opf")
opf = read_opf(file)

node = opf[1][1][1]

new_mesh = refmesh_to_mesh(node)

using MeshViz, GLMakie
viz(new_mesh)
```
"""
refmesh_to_mesh!, refmesh_to_mesh

function refmesh_to_mesh(node)
    if node[:geometry] !== nothing

        ref_mesh = node[:geometry].ref_mesh.mesh

        # Get the reference mesh and taper it in z and y (the principal axis is following x already):
        if node[:geometry].ref_mesh.taper # Taper only if enableScale="true" in the OPF: taper == true
            ref_mesh = taper(ref_mesh, node[:geometry].dUp, node[:geometry].dDwn)
        end

        scaled_mesh = Array{Meshes.Point3}(undef, Meshes.nvertices(ref_mesh))
        for (i, p) in enumerate(ref_mesh.points)
            scaled_mesh[i] = Meshes.Point3(node[:geometry].transformation(p.coords))
        end

        return Meshes.SimpleMesh(scaled_mesh, Meshes.topology(ref_mesh))
    else
        return nothing
    end
end

function refmesh_to_mesh!(node)
    if node[:geometry] !== nothing

        ref_mesh = node[:geometry].ref_mesh.mesh

        # Get the reference mesh and taper it in z and y (the principal axis is following x already):
        if node[:geometry].ref_mesh.taper # Taper only if enableScale="true" in the OPF: taper == true
            ref_mesh = taper(ref_mesh, node[:geometry].dUp, node[:geometry].dDwn)
        end

        scaled_mesh = Array{Meshes.Point3}(undef, Meshes.nvertices(ref_mesh))
        for (i, p) in enumerate(ref_mesh.points)
            scaled_mesh[i] = Meshes.Point3(node[:geometry].transformation(p.coords))
        end

        node[:geometry].mesh = Meshes.SimpleMesh(scaled_mesh, Meshes.topology(ref_mesh))
        return node[:geometry].mesh
    else
        return nothing
    end
end


text:
using Random
using LinearAlgebra
using Distances

"""
Gaussian/RBF/Squared Exponential correlation function
"""
function rbf(x, y, θ=1.0)
    1/sqrt(2*pi)*exp.(-θ*0.5*(norm(x .- y))^2)
end

"""
Gaussian/RBF/Squared Exponential correlation function
"""
function rbf_single(dist, θ=1.0)
    1/sqrt(2*pi)*exp.(-θ*0.5*(dist)^2)
end

"""
Derivative of Gaussian/RBF/Squared Exponential correlation function
with respect to hyperparameter theta
"""
function rbf_prime(x, y, θ=1.0)
    - 1/sqrt(2*pi) * 0.5 * norm(x.-y)^2 * exp.(-θ*0.5*(norm(x .- y))^2)
end

"""
Derivative of Gaussian/RBF/Squared Exponential correlation function
with respect to hyperparameter theta
"""
function rbf_prime_single(dist, θ=1.0)
    - 1/sqrt(2*pi) * 0.5 * dist^2 * exp.(-θ*0.5*(dist)^2)
end

"""
Second derivative of Gaussian/RBF/Squared Exponential correlation function
with respect to hyperparameter theta
"""
function rbf_prime2(x, y, θ=1.0)
    1/sqrt(2*pi) * 0.25 * norm(x.-y)^4 * exp.(-θ*0.5*(norm(x .- y))^2)
end

"""
Second derivative of Gaussian/RBF/Squared Exponential correlation function
with respect to hyperparameter theta
"""
function rbf_prime2_single(dist, θ=1.0)
    1/sqrt(2*pi) * 0.25 * dist^4 * exp.(-θ*0.5*(dist)^2)
end

"""
Builds cross-covariance matrix using corr kernel function.
Ex. RBF, exponential, Matern

Inputs are assumed to be 1D arrays or 2D column vectors.
Returns matrix of size len(s1) times len(s2)
"""
function K(s1, s2, θ, corr=rbf)
    K = zeros(size(s1, 1), size(s2, 1))
    if s1 != s2
        for i = 1:size(s1, 1)
            for j = 1:size(s2, 1)
                cur = corr(s1[i, :], s2[j, :], θ)
                K[i, j] = isa(cur, Array) ? cur[1] : cur
            end
        end
    else #s1==s2 can compute half the matrix and fill in other half
        for i = 1:size(s1, 1)
            for j = 1:i
                cur = corr(s1[i, :], s2[j, :], θ)
                K[i, j] = K[j, i] = isa(cur, Array) ? cur[1] : cur
            end
        end
    end
    return size(s1, 1)==size(s2, 1) ? K+UniformScaling(1e-8) : K
end

#rows of s1 and s2 contain the coordinates of each location
function fastK(s1, s2, θ, corr=rbf_single)
    corrtheta = dist -> corr(dist, θ)
    res = corrtheta.(pairwise(Euclidean(), s1, s2, dims=1))
    return size(s1, 1)==size(s2, 1) ? res + UniformScaling(1e-8) : res
end

"""
Samples from a Normal distribution with mean 0 and covariance matrix K
Scales result by 40. Works in 1D.
"""
function sample(K)
    c = cholesky(K)
    v = 40*randn(size(K, 1), 1)
    return c.L*(v)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module TestStats

using Bootstrap
using FactCheck

x = collect(100:-1:0) ## avoid sorting

## reference values, computed with boot:::norm_inter (R)
qr = (
      (1.0,      Inf),
      (0.0,     -Inf),
      (0.5,     50.0),
      (0.4713,  47.07265),
      (0.1,      9.205641),
      (0.052,    4.319606)
      )

facts("Normal-interpolated quantiles") do

    context("Reference values: alpha scalar") do
        for (alpha, ref) in qr
            @fact Bootstrap.iquantile(x, alpha) --> roughly(ref, 1e-5)
        end
    end

    context("Reference values: alpha vector") do
        a = Float64[x[1] for x in qr]
        r = Float64[x[2] for x in qr]
        y = Bootstrap.iquantile(x, a)
        ## workaround, since roughly cannot handle 'Inf' here
        y = round.(y, 4)
        r = round.(r, 4)
        @fact y --> r
    end

    context("Monotone quantiles") do
        x = collect(0:100)
        ## check that quantiles are monotone
        y0 = -Inf
        for alpha in 0:0.001:1
            y1 = Bootstrap.iquantile(x, alpha)
            @fact y1 --> greater_than_or_equal(y0)
            y0 = y1
        end
    end

end

end


text:
function ant(width, height)
    y, x = fld(height, 2), fld(width, 2)
    M = falses(height, width)

    dir = im
    for i in 0:100000
        x in 1:width && y in 1:height || break
        dir *= M[y, x] ? im : -im
        M[y, x] = !M[y, x]
        x, y = reim(x + im * y + dir)
    end

    for row in 1:size(M,1)
        println(mapreduce(x -> x ? 'x' : '.', *, M[row,:]))
    end
end

ant(100, 100)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# This file is part of the ValidatedNumerics.jl package; MIT licensed

type IntervalParameters

    precision_type::Type
    precision::Int
    rounding::Symbol
    pi::Interval{BigFloat}

    IntervalParameters() = new(BigFloat, 256, :narrow)  # leave out pi
end

const parameters = IntervalParameters()

## Precision:

doc"`big53` creates an equivalent `BigFloat` interval to a given `Float64` interval."
function big53(a::Interval{Float64})
    x = setprecision(Interval, 53) do  # precision of Float64
        convert(Interval{BigFloat}, a)
    end
end

setprecision(::Type{Interval}, ::Type{Float64}) = parameters.precision_type = Float64
# does not change the BigFloat precision

function setprecision{T<:AbstractFloat}(::Type{Interval}, ::Type{T}, prec::Integer)
    #println("SETTING BIGFLOAT PRECISION TO $precision")
    setprecision(BigFloat, prec)

    parameters.precision_type = T
    parameters.precision = prec
    parameters.pi = convert(Interval{BigFloat}, pi)

    prec
end

setprecision{T<:AbstractFloat}(::Type{Interval{T}}, prec) = setprecision(Interval, T, prec)

setprecision(::Type{Interval}, prec::Integer) = setprecision(Interval, BigFloat, prec)

function setprecision(f::Function, ::Type{Interval}, prec::Integer)

    old_precision = precision(Interval)
    setprecision(Interval, prec)

    try
        return f()
    finally
        setprecision(Interval, old_precision)
    end
end

# setprecision(::Type{Interval}, precision) = setprecision(Interval, precision)
setprecision(::Type{Interval}, t::Tuple) = setprecision(Interval, t...)

precision(::Type{Interval}) = (parameters.precision_type, parameters.precision)

const float_interval_pi = convert(Interval{Float64}, pi)  # does not change

pi_interval(::Type{BigFloat}) = parameters.pi
pi_interval(::Type{Float64})  = float_interval_pi

# Rounding for rational intervals, e.g for sqrt of rational interval:
# Find the corresponding AbstractFloat type for a given rational type

Base.float{T}(::Type{Rational{T}}) = typeof(float(one(Rational{T})))

# better to just do the following ?
# Base.float(::Type{Rational{Int64}}) = Float64
# Base.float(::Type{Rational{BigInt}}) = BigFloat

# Use that type for rounding with rationals, e.g. for sqrt:

if VERSION < v"0.5.0-dev+1182"

    function Base.with_rounding{T}(f::Function, ::Type{Rational{T}},
        rounding_mode::RoundingMode)
        setrounding(f, float(Rational{T}), rounding_mode)
    end

else
    function Base.setrounding{T}(f::Function, ::Type{Rational{T}},
        rounding_mode::RoundingMode)
        setrounding(f, float(Rational{T}), rounding_mode)
    end
end

float(x::Interval) =
    # @round(BigFloat, convert(Float64, x.lo), convert(Float64, x.hi))
    convert(Interval{Float64}, x)

## Change type of interval rounding:

doc"""`rounding(Interval)` returns the current interval rounding mode.
There are two possible rounding modes:

- :narrow  -- changes the floating-point rounding mode to `RoundUp` and `RoundDown`.
This gives the narrowest possible interval.

- :wide -- Leaves the floating-point rounding mode in `RoundNearest` and uses
`prevfloat` and `nextfloat` to achieve directed rounding. This creates an interval of width 2`eps`.
"""

rounding(::Type{Interval}) = parameters.rounding

function setrounding(::Type{Interval}, mode)
    if mode ∉ [:wide, :narrow]
        throw(ArgumentError("Only possible interval rounding modes are `:wide` and `:narrow`"))
    end

    parameters.rounding = mode  # a symbol
end

big{T}(x::Interval{T}) = convert(Interval{BigFloat}, x)


text:
using Pathfinder
using PSIS
using Random
using Test

@testset "resample" begin
    x = randn(100)

    rng = Random.seed!(Random.default_rng(), 42)
    xsub = Pathfinder.resample(rng, x, 500)
    @test all(in.(xsub, Ref(x)))

    # weight the first 20% much higher than the remaining 80%
    log_weights = randn(100)
    log_weights[1:20] .+= 1_000
    xsub, result = Pathfinder.resample(rng, x, log_weights, 500)
    @test all(in.(xsub, Ref(x[1:20])))
    @test result isa PSIS.PSISResult
end


text:
using SubsetSelection, DataFrames, LinearAlgebra

function split_data(n, holdout)
    batch = ceil(Int, n*holdout)
    perm = shuffle(collect(1:n))
    train = perm[1:(n-batch)]
    val = perm[(n-batch+1): n]
    return train, val
end

function error(ℓ::Regression, Y, X, indices, w)
  Z = Y .- X[:,indices]*w
  return dot(Z,Z)/size(Z,1)
end

function error(ℓ::Classification, Y, X, indices, w)
    df = DataFrame(Output = (Y.>0), Probability = 1. ./ (1 .+ exp.(X[:,indices]*w)))
    sort!(df, :Probability)
    # sort!(df, cols = :Probability)
    # df = sort(df, cols = :Probability)
    df[!,:Above]=(1:length(Y)) .- cumsum(df[:,:Output])

    nP = length(findall(o->o>0, Y)); nN = length(findall(o->o<0, Y));

    AUC = 1-sum(df[:,:Above].*df[:,:Output])/nP/nN
    return 1-AUC
end

function error2(ℓ::Classification, Y, X, indices, w)
    Y_truth = Y .> 0
    Y_pred = X[:,indices]*w .> 0

    MR = sum(1 .- Y_truth.*Y_pred .- (1 .- Y_truth).*(1 .- Y_pred)) / size(Y,1)
    return MR
end

function grad_primal(ℓ::LossFunction, Y, X, w, γ)
  n,k = size(X)
  g = zeros(k)
  for j in 1:k
    g[j] += w[j] / γ
  end
  for i in 1:n
    g += grad_loss(ℓ, Y[i], X[i,:], w)
  end
  return g
end

function grad_loss(ℓ::OLS, y, x, w)
  return -(y-dot(x,w))*x
end
function grad_loss(ℓ::L1SVM, y, x, w)
  return -y*(1-y*dot(x,w) > 0)*x
end
function grad_loss(ℓ::LogReg, y, x, w)
  z = exp.(y*dot(x,w))
  return -y/(1+z)*x
end

function value_primal(ℓ::LossFunction, Y, X, w, γ)
  n,k = size(X)
  v = dot(w,w)/2/γ
  for i in 1:n
    v += value_loss(ℓ, Y[i], X[i,:], w)
  end
  return v
end

function value_loss(ℓ::OLS, y, x, w)
  return .5*(y-dot(x,w))^2
end
function value_loss(ℓ::L1SVM, y, x, w)
  return max(0, 1-y*dot(x,w))
end
function value_loss(ℓ::LogReg, y, x, w)
  z = exp.(-y*dot(x,w))
  return log(1+z)
end


text:
using myJuliaUtils
using Test
using LinearAlgebra
using Random

# testing some helper function

# one matrix case (qq is a vector)
n = round(Int, 2 + rand() * 100)
qq = genPDmat(n)
qq = vech(qq)

A, sigma = myJuliaUtils.decVCV(qq)

@test vec2sym(qq) ≈ vec2ltriW1(A) * Diagonal(sigma).^2 * vec2ltriW1(A)'

# many matrix case (qq is a matrix)

n = 5
T = 100
c = round(Int, n * (n + 1) / 2)
qq = Array{Float64}(undef, c, T)

for i = 1:T
    q = genPDmat(n)
    qq[:,i] = vech(q)
end

A, sigma = myJuliaUtils.decVCV(qq)

for i = 1:T
    @test vec2sym(qq[:,i]) ≈ vec2ltriW1(A[:,i]) * Diagonal(sigma[:,i]).^2 * vec2ltriW1(A[:,i])'
end

A, sigma = myJuliaUtils.decVCV(qq,true)

for i = 1:T
    @test vec2sym(qq[:,i]) != vec2ltriW1(A[:,i]) * Diagonal(sigma[:,i]).^2 * vec2ltriW1(A[:,i])'
end

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/optimize.jl
# Load subproblem type optimization routines and define general optimize!
# throw_optimize_hook!.
#############################################################################

include("optimize_subsolve.jl")
include("nonconvex/presolve.jl")
include("nonconvex/branch/branch.jl")
include("nonconvex/optimize_local_nlp.jl")
include("nonconvex/optimize_relax_lp.jl")
include("nonconvex/preprocess.jl")
include("nonconvex/postprocess.jl")
include("nonconvex/optimize_nonconvex.jl")

throw_optimize_hook!(m::Optimizer) = optimize_hook!(_ext_type(m), m)

function MOI.optimize!(m::Optimizer)

    m._start_time = time()
    m._time_left = m._time_limit

    if !_enable_optimize_hook(m)

        _parse_classify_problem!(m)

        new_time = time() - m._start_time
        m._parse_time = new_time
        m._run_time = new_time

        optimize!(Val{m._problem_type}(), m)
    else
        throw_optimize_hook!(m)
    end

    return nothing
end


text:
using Test
using TightlyBound
using Suppressor
#=
include("../Utility.jl")
include("../BandTools.jl")
include("../Atomic.jl")
include("../Atomdata.jl")
include("../Crystal.jl")
include("../DFToutMod.jl")
include("../TB.jl")
include("../CalcTB.jl")
include("../FitTB.jl")
=#

#basic loading

#include("../includes_laguerre.jl")

function test1()

    @testset "testing laguerre fitting fake example" begin

        c1 = makecrys([7.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0], ["Li"], units="Bohr");
        c2 = makecrys([8.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0], ["Li"], units="Bohr");
        c3 = makecrys([8.5 0 0; 0 8.5 0; 0 0 7.0], [0 0 0], ["Li"], units="Bohr");
        c4 = makecrys([6.0 0 0; 0 6.0 0; 0 0 6.0], [0 0 0], ["Li"], units="Bohr");

        c5 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.49 0 0 ], ["Li", "Li"], units="Bohr");
        c6 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.48 0 0 ], ["Li", "Li"], units="Bohr");
        c7 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.47 0 0 ], ["Li", "Li"], units="Bohr");
        c8 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.46 0 0 ], ["Li", "Li"], units="Bohr");
        c9 = makecrys([14.0 0 0; 0 7.0 0; 0 0 7.0], [0 0 0; 0.45 0 0 ], ["Li", "Li"], units="Bohr");

        @suppress begin
#        if true

            database = Dict()
            database[(:Li, :Li)] = TightlyBound.CalcTB.make_coefs(Set(["Li", "Li"]), 2)

            tbc1 = TightlyBound.CalcTB.calc_tb_fast(c1, database, use_threebody=false);
            tbc2 = TightlyBound.CalcTB.calc_tb_fast(c2, database, use_threebody=false);
            tbc3 = TightlyBound.CalcTB.calc_tb_fast(c3, database, use_threebody=false);
            tbc4 = TightlyBound.CalcTB.calc_tb_fast(c4, database, use_threebody=false);
            tbc5 = TightlyBound.CalcTB.calc_tb_fast(c5, database, use_threebody=false);
            tbc6 = TightlyBound.CalcTB.calc_tb_fast(c6, database, use_threebody=false);
            tbc7 = TightlyBound.CalcTB.calc_tb_fast(c7, database, use_threebody=false);
            tbc8 = TightlyBound.CalcTB.calc_tb_fast(c8, database, use_threebody=false);
            tbc9 = TightlyBound.CalcTB.calc_tb_fast(c9, database, use_threebody=false);

            tbc_list = [tbc1 tbc2 tbc3 tbc4 tbc5 tbc6 tbc7 tbc8 tbc9]
            newdatabase = TightlyBound.FitTB.do_fitting(tbc_list, fit_threebody=false, do_plot=false)
            #        println("newdartabase")
            #        println(newdatabase[("Li", "Li")])
            @test sum(abs.(newdatabase[(:Li, :Li)].datH .- database[(:Li, :Li)].datH)) ≤ 1e-5
        end

    end
end

test1()


text:
#####################
# Compilers to LogicCircuits data structures starting from already parsed line objects
#####################

"""
Compile lines into a unstructured logical circuit
"""
compile_logical(lines::CircuitFormatLines)::UnstLogicalΔ =
    compile_logical_m(lines)[1]

"""
Compile lines into a unstructured logical circuit,
while keeping track of id-to-node mappings
"""
function compile_logical_m(lines::CircuitFormatLines)

    # linearized circuit nodes
    circuit = Vector{UnstLogicalΔNode}()

    # mapping from circuit node ids to node objects
    id2node = Dict{ID,UnstLogicalΔNode}()

    # literal cache is responsible for making leaf literal nodes unique and adding them to `circuit`
    lit_cache = Dict{Lit,LogicalLeafNode}()
    literal_node(l::Lit) = get!(lit_cache, l) do
        leaf = LiteralNode(l)
        push!(circuit,leaf) # also add new leaf to linearized circuit before caller
        leaf
    end

    seen_true = false
    seen_false = false

    true_node = TrueNode()
    false_node = FalseNode()

    function compile(ln::CircuitFormatLine)
        error("Compilation of line $ln is not supported")
    end
    function compile(::Union{CircuitHeaderLine,CircuitCommentLine})
         # do nothing
    end
    function compile(ln::WeightedLiteralLine)
        # Weighted Literal are implicitly an OR node
        # Here making that explicit in the Circuit
        lit_node = literal_node(literal(ln))
        or_node = ⋁Node([lit_node])
        push!(circuit, lit_node)
        push!(circuit, or_node)
        id2node[ln.node_id] = or_node
    end
    function compile(ln::LiteralLine)
        id2node[ln.node_id] = literal_node(literal(ln))
    end
    function compile(ln::ConstantLine)
        if constant(ln) == true
            n = true_node
            seen_true || (push!(circuit,n); seen_true = true)
        else
            n = false_node
            seen_false || (push!(circuit,n); seen_false = true)
        end
        id2node[ln.node_id] = n
    end
    function compile_elements(e::Element)
        n = ⋀Node([id2node[e.prime_id],id2node[e.sub_id]])
        push!(circuit,n)
        n
    end
    function compile(ln::DecisionLine)
        n = ⋁Node(map(compile_elements, ln.elements))
        push!(circuit,n)
        id2node[ln.node_id] = n
    end
    function compile(ln::BiasLine)
        n = ⋁Node([circuit[end]])
        push!(circuit,n)
        id2node[ln.node_id] = n
    end

    for ln in lines
        compile(ln)
    end

    return circuit, id2node
end

#TODO add compile_struct_logical

"""
Compile lines into a smooth unstructured logical circuit
"""
compile_smooth_logical(lines::CircuitFormatLines)::UnstLogicalΔ =
    compile_smooth_logical_m(lines)[1]

"""
Compile lines into a smooth unstructured logical circuit,
while keeping track of id-to-node mappings
"""
function compile_smooth_logical_m(lines::CircuitFormatLines)

    # linearized circuit nodes
    circuit = Vector{UnstLogicalΔNode}()

    # mapping from circuit node ids to node objects
    id2node = Dict{ID,UnstLogicalΔNode}()

    # literal cache is responsible for making leaf literal nodes unique and adding them to `circuit`
    lit_cache = Dict{Lit,LogicalLeafNode}()
    literal_node(l::Lit) = get!(lit_cache, l) do
        leaf = LiteralNode(l)
        push!(circuit,leaf) # also add new leaf to linearized circuit before caller
        leaf
    end

    smoothing_warning = "Cannot compile a smooth logical circuit from lines that are not normalized: there is no way to smooth without knowing the variable scope. Instead compile a non-smooth logical circuit and smooth it afterwards."

    function compile(ln::CircuitFormatLine)
        error("Compilation of line $ln is not supported")
    end
    function compile(::Union{CircuitHeaderLine,CircuitCommentLine})
         # do nothing
    end
    function compile(ln::WeightedLiteralLine)
        # Weighted Literal are implicitly an OR node
        # Here making that explicit in the Circuit
        @assert is_normalized(ln) " $smoothing_warning"
        lit_node = literal_node(literal(ln))
        or_node = ⋁Node([lit_node])
        push!(circuit, lit_node)
        push!(circuit, or_node)
        id2node[ln.node_id] = or_node
    end
    function compile(ln::LiteralLine)
        @assert is_normalized(ln) " $smoothing_warning"
        id2node[ln.node_id] = literal_node(literal(ln))
    end

    function compile(ln::WeightedNamedConstantLine)
        @assert constant(ln) == true
        # because we promise to compile a smooth circuit, here we need to add a "smoothing or gate"
        n = ⋁Node([literal_node(var2lit(variable(ln))),
                    literal_node(-var2lit(variable(ln)))])
        push!(circuit,n)
        id2node[ln.node_id] = n
    end
    function compile(ln::AnonymousConstantLine)
        error(smoothing_warning)
    end
    function compile_elements(e::TrimmedElement)
        error(smoothing_warning)
    end
    function compile_elements(e::NormalizedElement)
        n = ⋀Node([id2node[e.prime_id],id2node[e.sub_id]])
        push!(circuit,n)
        n
    end
    function compile(ln::DecisionLine)
        n = ⋁Node(map(compile_elements, ln.elements))
        push!(circuit,n)
        id2node[ln.node_id] = n
    end
    function compile(ln::BiasLine)
        n = ⋁Node([circuit[end]])
        push!(circuit,n)
        id2node[ln.node_id] = n
    end

    for ln in lines
        compile(ln)
    end

    return circuit, id2node
end

"""
Compile circuit and vtree lines into a structured logical circuit with its vtree
"""
function compile_smooth_struct_logical(circuit_lines::CircuitFormatLines,
                                vtree_lines::VtreeFormatLines)
    compile_smooth_struct_logical_m(circuit_lines,vtree_lines)[1:2]
end

"""
Compile circuit and vtree lines into a structured logical circuit with its vtree,
while keeping track of id-to-node mappings
"""
function compile_smooth_struct_logical_m(circuit_lines::CircuitFormatLines,
                                  vtree_lines::VtreeFormatLines)
    vtree, id2vtree = compile_vtree_format_lines_m(vtree_lines)
    circuit, id2circuit = compile_smooth_struct_logical_m(circuit_lines, id2vtree)
    return circuit, vtree, id2vtree, id2circuit
end

"""
Compile circuit lines and vtree node mapping into a structured logical circuit,
while keeping track of id-to-node mappings
"""
function compile_smooth_struct_logical_m(lines::CircuitFormatLines,
                                         id2vtree::Dict{ID, PlainVtreeNode})

    # linearized circuit nodes
    circuit = Vector{StructLogicalΔNode{PlainVtreeNode}}()

    # mapping from node ids to node objects
    id2node = Dict{ID,StructLogicalΔNode{PlainVtreeNode}}()

    # literal cache is responsible for making leaf literal nodes unique and adding them to `circuit`
    lit_cache = Dict{Lit,StructLogicalLeafNode{PlainVtreeNode}}()
    literal_node(l::Lit, v::PlainVtreeNode) = get!(lit_cache, l) do
        leaf = StructLiteralNode{PlainVtreeNode}(l,v)
        push!(circuit,leaf) # also add new leaf to linearized circuit before caller
        leaf
    end

    smoothing_warning = "Cannot compile a smooth logical circuit from lines that are not normalized: functionality to determine variable scope and perform smoothing not implemented in the line compiler.  Instead compile a non-smooth logical circuit and smooth it afterwards."

    function compile(ln::CircuitFormatLine)
        error("Compilation of line $ln is not supported")
    end
    function compile(::Union{CircuitHeaderLine,CircuitCommentLine})
         # do nothing
    end
    function compile(ln::WeightedLiteralLine)
        # Weighted Literal are implicitly an OR node
        # Here making that explicit in the Circuit
        @assert is_normalized(ln) smoothing_warning
        lit_node = literal_node(ln.literal, id2vtree[ln.vtree_id])
        or_node = Struct⋁Node{PlainVtreeNode}([lit_node], id2vtree[ln.vtree_id])

        push!(circuit, lit_node)
        push!(circuit, or_node)
        id2node[ln.node_id] = or_node
    end
    function compile(ln::LiteralLine)
        @assert is_normalized(ln) smoothing_warning
        id2node[ln.node_id] = literal_node(ln.literal, id2vtree[ln.vtree_id])
    end

    function compile(ln::ConstantLine)
        vtree = id2vtree[ln.vtree_id]
        if is_normalized(ln)
            variable = (vtree::PlainVtreeLeafNode).var
            @assert !(ln isa WeightedNamedConstantLine) || variable == ln.variable "PlainVtree mapping must agree with variable field of circuit line"
        else
            error(smoothing_warning)
        end
        if constant(ln) == true
            # because we promise to compile a smooth circuit, here we need to add an or gate
            n = Struct⋁Node{PlainVtreeNode}([literal_node(var2lit(variable), vtree), literal_node(-var2lit(variable), vtree)], vtree)
        else
            error("False leaf logical circuit nodes not yet implemented")
        end
        push!(circuit,n)
        id2node[ln.node_id] = n
    end
    function compile_elements(e::TrimmedElement, ::PlainVtreeNode)
        error(smoothing_warning)
    end
    function compile_elements(e::NormalizedElement, v::PlainVtreeNode)
        n = Struct⋀Node{PlainVtreeNode}([id2node[e.prime_id], id2node[e.sub_id]], v)
        push!(circuit,n)
        n
    end
    function compile(ln::DecisionLine)
        vtree = id2vtree[ln.vtree_id]
        n = Struct⋁Node{PlainVtreeNode}(map(e -> compile_elements(e, vtree), ln.elements), vtree)
        push!(circuit,n)
        id2node[ln.node_id] = n
    end
    function compile(ln::BiasLine)
        n = Struct⋁Node{PlainVtreeNode}([circuit[end]], circuit[end].vtree)
        push!(circuit,n)
        id2node[ln.node_id] = n
    end

    for ln in lines
        compile(ln)
    end

    return circuit, id2node
end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

"""
    ObjectiveContainer{T}

A helper struct to simplify the handling of objective functions in
Utilities.Model.
"""
mutable struct ObjectiveContainer{T} <: MOI.ModelLike
    is_sense_set::Bool
    sense::MOI.OptimizationSense
    is_function_set::Bool
    single_variable::Union{Nothing,MOI.VariableIndex}
    scalar_affine::Union{Nothing,MOI.ScalarAffineFunction{T}}
    scalar_quadratic::Union{Nothing,MOI.ScalarQuadraticFunction{T}}
    function ObjectiveContainer{T}() where {T}
        o = new{T}()
        MOI.empty!(o)
        return o
    end
end

function MOI.empty!(o::ObjectiveContainer{T}) where {T}
    o.is_sense_set = false
    o.sense = MOI.FEASIBILITY_SENSE
    o.is_function_set = false
    o.single_variable = nothing
    o.scalar_affine =
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], zero(T))
    o.scalar_quadratic = nothing
    return
end

function MOI.is_empty(o::ObjectiveContainer)
    return !o.is_sense_set && !o.is_function_set
end

###
### ObjectiveSense
###

MOI.supports(::ObjectiveContainer, ::MOI.ObjectiveSense) = true

MOI.get(o::ObjectiveContainer, ::MOI.ObjectiveSense) = o.sense

function MOI.set(o::ObjectiveContainer, ::MOI.ObjectiveSense, value)
    if value == MOI.FEASIBILITY_SENSE
        MOI.empty!(o)
    end
    o.is_sense_set = true
    o.sense = value
    return
end

###
### ObjectiveFunctionType
###

function MOI.get(
    o::ObjectiveContainer{T},
    ::MOI.ObjectiveFunctionType,
) where {T}
    if o.single_variable !== nothing
        return MOI.VariableIndex
    elseif o.scalar_quadratic !== nothing
        return MOI.ScalarQuadraticFunction{T}
    end
    @assert o.scalar_affine !== nothing
    return MOI.ScalarAffineFunction{T}
end

###
### ObjectiveFunction
###

function MOI.supports(
    ::ObjectiveContainer{T},
    ::MOI.ObjectiveFunction{
        <:Union{
            MOI.VariableIndex,
            MOI.ScalarAffineFunction{T},
            MOI.ScalarQuadraticFunction{T},
        },
    },
) where {T}
    return true
end

function MOI.get(
    o::ObjectiveContainer{T},
    ::MOI.ObjectiveFunction{F},
) where {T,F}
    if o.single_variable !== nothing
        return convert(F, o.single_variable)
    elseif o.scalar_quadratic !== nothing
        return convert(F, o.scalar_quadratic)
    end
    @assert o.scalar_affine !== nothing
    return convert(F, o.scalar_affine)
end

function MOI.set(
    o::ObjectiveContainer,
    ::MOI.ObjectiveFunction{MOI.VariableIndex},
    f::MOI.VariableIndex,
)
    o.is_function_set = true
    o.single_variable = copy(f)
    o.scalar_affine = nothing
    o.scalar_quadratic = nothing
    return
end

function MOI.set(
    o::ObjectiveContainer{T},
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}},
    f::MOI.ScalarAffineFunction{T},
) where {T}
    o.is_function_set = true
    o.single_variable = nothing
    o.scalar_affine = copy(f)
    o.scalar_quadratic = nothing
    return
end

function MOI.set(
    o::ObjectiveContainer{T},
    ::MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{T}},
    f::MOI.ScalarQuadraticFunction{T},
) where {T}
    o.is_function_set = true
    o.single_variable = nothing
    o.scalar_affine = nothing
    o.scalar_quadratic = copy(f)
    return
end

###
### MOI.ListOfModelAttributesSet
###

function MOI.get(o::ObjectiveContainer, ::MOI.ListOfModelAttributesSet)
    ret = MOI.AbstractModelAttribute[]
    if o.is_sense_set
        push!(ret, MOI.ObjectiveSense())
    end
    if o.is_function_set
        F = MOI.get(o, MOI.ObjectiveFunctionType())
        push!(ret, MOI.ObjectiveFunction{F}())
    end
    return ret
end

###
### MOI.modify
###

function MOI.modify(
    o::ObjectiveContainer,
    ::MOI.ObjectiveFunction,
    change::MOI.AbstractFunctionModification,
)
    if o.single_variable !== nothing
        o.single_variable = modify_function(o.single_variable, change)
    elseif o.scalar_quadratic !== nothing
        o.scalar_quadratic = modify_function(o.scalar_quadratic, change)
    else
        @assert o.scalar_affine !== nothing
        o.is_function_set = true
        o.scalar_affine = modify_function(o.scalar_affine, change)
    end
    return
end

###
### MOI.delete
###

function MOI.delete(o::ObjectiveContainer, x::MOI.VariableIndex)
    if o.single_variable !== nothing
        if x == o.single_variable
            sense = o.sense
            MOI.empty!(o)
            if o.is_sense_set
                MOI.set(o, MOI.ObjectiveSense(), sense)
            end
        end
    elseif o.scalar_quadratic !== nothing
        o.scalar_quadratic = remove_variable(o.scalar_quadratic, x)
    else
        @assert o.scalar_affine !== nothing
        o.scalar_affine = remove_variable(o.scalar_affine, x)
    end
    return
end

function MOI.delete(o::ObjectiveContainer, x::Vector{MOI.VariableIndex})
    keep = v -> !(v in x)
    if o.single_variable !== nothing
        if o.single_variable in x
            sense = o.sense
            MOI.empty!(o)
            if o.is_sense_set
                MOI.set(o, MOI.ObjectiveSense(), sense)
            end
        end
    elseif o.scalar_quadratic !== nothing
        o.scalar_quadratic = filter_variables(keep, o.scalar_quadratic)
    else
        @assert o.scalar_affine !== nothing
        o.scalar_affine = filter_variables(keep, o.scalar_affine)
    end
    return
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains types and methods to display vertex and edge dataframes.

###
# TODO: Avoid reliance on DataFrames internal methods for display. Try to avoid label computation for all vertices
###
################################################# IMPORT/EXPORT ############################################################

export VertexDescriptor, EdgeDescriptor

function Base.show(io::IO, g::Graph)
   write(io, "Graph($(nv(g)) vertices, $(ne(g)) edges, $(listvprops(g)) vertex properties, $(listeprops(g)) edge properties)")
end

################################################# VERTEXDESCRIPTOR ########################################################

immutable VertexDescriptor
   g::Graph
end

function Base.show(io::IO, V::VertexDescriptor)
   g = V.g
   N = nrow(vdata(g))
   if N == 0
      println("No vertex properties to show")
   else
      Vdata = hcat(DataFrame(Labels=encode(g)), vdata(g))
      show(io, Vdata, true, :VertexID, false)
   end
end

function Base.getindex(V::VertexDescriptor, x)
   g = V.g
   v = decode(g, x)
   println("Vertex $x")
   for prop in listvprops(g)
      println("$prop => $(getvprop(g, v, prop))")
   end
end

Base.getindex(V::VertexDescriptor, x::AbstractArray) = error("Vector getindex isn't supported, try a filter instead")

################################################# EDGEDESCRIPTOR ##########################################################

immutable EdgeDescriptor
   g::Graph
end

function Base.show(io::IO, E::EdgeDescriptor)
   g = E.g
   N = nrow(edata(g))

   if N == 0
      println("No edge properties to show")
   else
      eit = edges(g)
      uls = encode(g, eit.us)
      vls = encode(g, eit.vs)
      Edata = hcat(DataFrame(Source=uls, Target=vls), edata(g))
      show(io, Edata, true, :Index, false)
   end
end

function Base.getindex(E::EdgeDescriptor, x, y)
   g = E.g
   e = decode(g, x=>y)
   if hasedge(g, e)
      println("Edge $x => $y")
      for prop in listeprops(g)
         println("$prop => $(geteprop(g, e, prop))")
      end
   else
      error("Invalid edge $x => $y")
   end
end

Base.getindex(E::EdgeDescriptor, x::AbstractArray) = error("Vector getindex isn't supported, try a filter instead")
Base.getindex(E::EdgeDescriptor, x::AbstractArray, y::AbstractArray) = error("Vector getindex isn't supported, try a filter instead")
################################################# SPLICED ASSIGNMENT ######################################################

# Make the graph type iterable
Base.start(g::Graph) = 1
Base.done(g::Graph, i) = i == 3

function Base.next(g::Graph, i)
   i == 1 && return VertexDescriptor(g), 2
   i == 2 && return EdgeDescriptor(g), 3
end


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AbstractPilotResponseImpl

export AbstractPilotResponse

abstract AbstractPilotResponse

end

module AbstractPilotResponseInterfaces

export updatePilotResponse

function updatePilotResponse()
end

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
#!/usr/bin/julia
function value(list)
    result = 0
    for v in list
        result = result * 2
        result = result + v
    end
    return result
end

function find_most(candidates, index)
    count = sum([v[index] for v in candidates])
    return count >= length(candidates) / 2 ? 1 : 0
end

function find_least(candidates, index)
    count = sum([v[index] for v in candidates])
    return count < length(candidates) / 2 ? 1 : 0
end

input_file = ARGS[1]
open(input_file) do io
    data = map.(x -> parse(Int64,x), map(collect, readlines(io)))
    candidates = data
    index = 1
    while length(candidates) > 1
        most = find_most(candidates, index)
        candidates = filter(v -> v[index] == most, candidates)
        index = index + 1
    end
    oxy = value(candidates[1])
    candidates = data
    index = 1
    while length(candidates) > 1
        least = find_least(candidates, index)
        candidates = filter(v -> v[index] == least, candidates)
        index = index + 1
    end
    co2 = value(candidates[1])
    println(oxy * co2)
end


text:
"""
    EmptyShape{T,Dim}

An empty shape with no points.
"""
struct EmptyShape{T,Dim} <: Shape{T,Dim} end

EmptyShape(s::Shape{T,Dim}) where {T,Dim} = EmptyShape{T,Dim}()

name(shape::EmptyShape) = "EmptyShape"

outer_radius(c::EmptyShape{T,Dim}) where {T,Dim} = zero(T)
volume(shape::EmptyShape{T,Dim}) where {T,Dim} = zero(T)

import Base.issubset
issubset(c1::Shape, c2::EmptyShape) = false
issubset(c1::EmptyShape, c2::Shape) = false

import Base.in
in(x::AbstractVector, c::EmptyShape) = false


text:
# All functionality corresponding to the correlation coefficient of the
# convergent cross mapping (CCM) as Loss function Γ.

## Constructors:

"""
    CCM_ρ <: AbstractLoss

    Constructor for the CCM_ρ loss function (correlation coefficient of the
    convergent cross mapping) based on Sugihara et al.[^Sugihara2012], see also
    [`ccm`](@ref). In this case MCDTS tries to maximize the correlation coefficient
    of the convergent cross mapping from the input `data` and `Y_CCM`, the time
    series CCM should cross map to from `data` (see [`mcdts_embedding`](@ref)).

    ## Fieldnames
    * `timeseries`: The time series CCM should cross map to.
    * `threshold:: = 1: A threshold for the sufficient correlation of the
      cross-mapped values and the true values from `Y_CMM` for the current embedding.
      When the correlation coefficient exeeds this threshold in an embedding cycle
      the embedding stops.
    * `samplesize::Real = 1.`: the fraction of all phase space points
     to be considered in the computation of CCM_ρ.

    ## Defaults
    * When calling `CCM_ρ(timeseries)`, a CCM_ρ-object is created, storing
      `timeseries`, which is considered to be causally depended and the
      `threshold=1` is used, i.e. no threshold, since the correlation coefficient
      can not exceed 1.

    [^Sugihara2012]: Sugihara et al., [Detecting Causality in Complex Ecosystems. Science 338, 6106, 496-500](https://doi.org/10.1126/science.1227079)
"""
struct CCM_ρ <: AbstractLoss
    timeseries::Vector
    threshold::AbstractFloat
    samplesize::Real

    # Constraints and Defaults
    CCM_ρ(x,y=1.0,z=1) = begin
        @assert 0 < y ≤ 1 "Threshold must be a ∈ (0 1]"
        @assert 0 < z ≤ 1. "The samplesize must be in the interval (0 1]."
        typeof(y) <: Int ? new(x,convert(AbstractFloat, -y),z) : new(x,-y,z)
    end
end

## Functions:

"""
    Return the loss based on the negative correlation coefficient for CCM.
"""
function compute_loss(Γ::CCM_ρ, Λ::AbstractDelayPreselection, dps::Vector{P}, Y_act::Dataset{D, T}, Ys, τs, w::Int, ts::Int, τ_vals, ts_vals; metric=Euclidean(), kwargs...) where {P, D, T}

    Y_other = DelayEmbeddings.standardize(Γ.timeseries)
    samplesize = Γ.samplesize

    max_idx = get_max_idx(Λ, dps, τ_vals, ts_vals, ts) # get the candidate delays
    isempty(max_idx) && return Float64[], Int64[], []

    ρ_CCM = zeros(Float64, length(max_idx))

    for (i,τ_idx) in enumerate(max_idx)
        # create candidate phase space vector for this peak/τ-value
        tau_trials = (τ_vals...,τs[τ_idx-1],)
        ts_trials = (ts_vals...,ts,)
        Y_trial = genembed(Ys, tau_trials.*(-1), ts_trials)
        # account for value-shift due to negative lags
        Ys_other = Y_other[1+maximum(tau_trials.*(-1)):length(Y_trial)+maximum(tau_trials.*(-1))]
        # compute ρ_CCM for Y_trial and Y_other
        ρ_CCM[i], _, _ = MCDTS.ccm(Y_trial, Ys_other; metric = metric, w = w, samplesize = samplesize)
    end
    return -ρ_CCM, max_idx, [nothing for i in max_idx]
end

"""
    ccm(X, y; kwargs...) → ρ, y_hat, y_idx

    Compute the convergent crossmapping (CCM) [^Sugihara2012](Sugihara et al. 2012) of a
    vector time series `X` (an embedded time series `x`) on the time series `y`
    NOTE: 'X' and 'y' must have the same length and you have to make sure that
    'y' starts at the same time index as 'X' does. - When using [`genembed`](@ref)
    with negative delays to construct `X` from `x`, which is mandatory here, then
    'y' needs to be shifted by the largest negative delay value, which has been
    used to construct `X`.

    Returns the correlation coefficient of `y` and its predicted values for `y_hat`,
    based on the nearest neighbour structure of `X`. `y_idx` are the corresponding
    indices, which have been used for computing `y_hat`.
    It is said that 'y' causes 'x', if ρ increases with increasing time series
    length AND ρ is "quite high".

    Keyword arguments:
    *`metric = Euclidean()`: The metric for vector distance computation.
    *`w::Int = 1`: The Theiler window in sampling units.
    *`lags::Array = [0]`: The lag for the cross mapping, in order to detect time lagged
                          causal relationships. The output ρ is an array of size
                          `length(lags)`, the output Y_hat is the one corresponding
                          to a lag of zero.
    * `samplesize::Real = 0.1`: fraction of all phase space points (=`length(X)`)
                                to be considered (fiducial points v)

"""
function ccm(X::Dataset{D,T},Y::Vector{T}; metric = Euclidean(), w::Int = 1,
    lags::AbstractArray = [0], samplesize::Real = 1.) where {D,T<:Real}

    K = D+1
    @assert length(X)==length(Y)
    @assert 0 < samplesize ≤ 1. "`samplesize` for computing ρ-CCM must be ∈ (0,1]"
    XX = Matrix(X)
    N = length(X)

    # consider subset of state space points
    NN = floor(Int, samplesize*N)
    if samplesize == 1
        ns = 1:N # the fiducial point indices
    else
        ns = sample(1:N, NN; replace=false) # the fiducial point indices
    end

    vxs = X[ns] # the fiducial points in the data set
    vtree = KDTree(X, metric)
    allNNidxs, allNNdist = DelayEmbeddings.all_neighbors(vtree, vxs, ns, K, w)
    Y_hat = zeros(T, NN) # preallocation
    # loop over each fiducial point
    for (i,v) in enumerate(vxs)
        NNidxs = allNNidxs[i] # indices of k nearest neighbors to v
        NNdist = allNNdist[i] # indices of k nearest neighbors to v
        # determine weights
        u = zeros(K)
        ws = zeros(K)
        @inbounds for (k, j) in enumerate(NNdist)
            u[k] = exp(-(j/NNdist[1]))
        end
        ws = u ./ sum(u)
        # compute Y_hat as a wheighted mean
        Y_hat[i] = sum(ws .* Y[NNidxs])
    end
    ρ = Statistics.cor(Y_hat, Y[ns])
    return ρ, Dataset(Y_hat), ns
end

max_stars_repo_path: src/Cost functions/CCM.jl

max_issues_repo_path: src/Cost functions/CCM.jl

max_forks_repo_path: src/Cost functions/CCM.jl


text:
################################################################################
# Initialization
#

using CUDA, Base.Test

include("perfutil.jl")

#
# Set-up
#

dev = CuDevice(0)
ctx = CuContext(dev)

cgctx = CuCodegenContext(ctx, dev)
ENV["ARCH"] = cgctx.arch
include("kernels/load.jl")

################################################################################
# Smoke testing
#

#
# ptx loading
#

launch(reference_dummy(), 1, 1, ())

#
# @cuda macro
#

@cuda (1, 1) kernel_dummy()

# kernel dims
@target ptx kernel_empty() = return nothing
@test_throws ErrorException @eval begin
    @cuda (0, 0) kernel_empty()
end
@eval begin
    @cuda (1, 1) kernel_empty()
end

# kernel name
@test_throws ErrorException @eval begin
    @cuda (1, 1) Module.kernel_foo()
end
@test_throws ErrorException @eval begin
    @cuda (1, 1) InvalidPrefixedKernel()
end

# external kernel
module KernelModule
    export kernel_empty2
    @target ptx kernel_empty2() = return nothing
end
@eval begin
    using KernelModule
    @cuda (1, 1) kernel_empty2()
end

# unlowered call
i = 0
@timeit begin # setup
        i += 1
        fname = symbol("kernel_dummy_$i")
        @eval @target ptx $fname() = return nothing
    end begin # benchmark
        @eval @cuda (1, 1) $fname()
    end begin # verification
    end begin # teardown
    end "macro_full" "unlowered call to @cuda (staged function execution, kernel compilation, runtime API interactions and asynchronous kernel launch)"
synchronize(ctx)

# pre-lowered call
@target ptx kernel_dummy() = return nothing
@timeit begin # setup
    end begin # benchmark
        @cuda (1, 1) kernel_dummy()
    end begin # verification
    end begin # teardown
    end "macro_lowered" "lowered call to @cuda (runtime API interactions and asynchronous kernel launch)"
synchronize(ctx)

#
# Argument passing
#

dims = (512, 512)
len = prod(dims)

# manual allocation
@timeit begin # setup
        input = round(rand(Float32, dims) * 100)
    end begin # benchmark
        input_dev = CuArray(input)
        output_dev = CuArray(Float32, dims)

        launch(reference_copy(), len, 1, (input_dev.ptr, output_dev.ptr))
        output = to_host(output_dev)

        free(input_dev)
        free(output_dev)
    end begin # verification
        @test_approx_eq input output
    end begin # teardown
    end "copy_reference" "vector copy reference execution"
@timeit begin # setup
        input = round(rand(Float32, dims) * 100)
    end begin # benchmark
        input_dev = CuArray(input)
        output_dev = CuArray(Float32, dims)

        @cuda (len, 1) kernel_copy(input_dev, output_dev)
        output = to_host(output_dev)

        free(input_dev)
        free(output_dev)
    end begin # verification
        @test_approx_eq input output
    end begin # teardown
    end "copy_manual" "vector copy on manually allocated GPU arrays"

# auto-managed host data
@timeit begin # setup
        input = round(rand(Float32, dims) * 100)
        output = Array(Float32, dims)
    end begin # benchmark
        @cuda (len, 1) kernel_copy(CuIn(input), CuOut(output))
    end begin # verification
        @test_approx_eq input output
    end begin # teardown
    end "copy_managed" "vector copy on managed GPU arrays"

# auto-managed host data, without specifying type
let
    input = round(rand(Float32, dims) * 100)
    output = Array(Float32, dims)

    @cuda (len, 1) kernel_copy(input, output)
    @test_approx_eq input output
end

# auto-managed host data, without specifying type, not using containers
let
    input = rand(Float32, dims)
    output = Array(Float32, dims)

    @cuda (len, 1) kernel_copy(round(input*100), output)
    @test_approx_eq round(input*100) output
end

# scalar through single-value array
let
    arr = round(rand(Float32, dims) * 100)
    val = Float32[0]

    @cuda (len, 1) kernel_lastvalue(CuIn(arr), CuOut(val))
    @test_approx_eq arr[dims...] val[1]
end

destroy(cgctx)


text:
export graph_degopt, get_topo_order_degopt

"""
    (graph,crefs)=graph_degopt(k;T=ComplexF64,input=:A)
    (graph,crefs)=graph_degopt(x,z;input=:A)
    (graph,crefs)=graph_degopt(d::Degopt;input=:A)

Corresponds to the (for a fixed number of multiplications) degree-optimal
polynomial

    B1=A
    B2=(x *I+x *A)(x *I+x *A)
    B3=(x *I+x *A+x*B2)(x *I+x *A+x *B2)
     ..

and

    Out=z*I+z*A1+z*B2+z*B3+z*B4...

The `x`-values are given in the argument `x`, which is a
`Vector{Tuple{Vector,Vector}}`, containing the elements of each sum. The
`z`-vector contains the elements to form the output, and `input` determines the
name of the matrix A above. If the parameter `k` is supplied instead of the
coefficients, all coeffs will be set to one. The general recursion is given in
(9) in the paper referenced below.

**Reference**

1. P. Bader, S. Blanes, and F. Casas, "Computing the matrix exponential
   with an optimized Taylor polynomial approximation", Mathematics, 7(12), 2019.
   DOI: [10.3390/math7121174](https://doi.org/10.3390/math7121174)
"""
function graph_degopt(x, z; input = :A)
    T = promote_type(eltype(eltype(eltype(x))), eltype(z))
    (graph, crefs) = graph_degopt_B(x, T, input = input)

    k = size(x, 1)

    # Add the polynomial in the Bk coeffs

    z_nodes = [:I; input]
    for s = 2:k+1
        push!(z_nodes, Symbol("B$s"))
    end
    key = Symbol("T2k$(k+3)")
    crefs_new = add_sum!(graph, key, z, z_nodes, Symbol("T2k"))
    append!(crefs, crefs_new)

    # Set the output
    empty!(graph.outputs)
    push!(graph.outputs, key)

    return (graph, crefs)
end
function graph_degopt(k; T = ComplexF64, input = :A)
    x = Vector{Tuple{Vector{T},Vector{T}}}(undef, k)
    for i = 1:k
        x[i] = (ones(T, i + 1), ones(T, i + 1))
    end

    z = ones(T, k + 2)

    return graph_degopt(x, z; input = :A)
end
function graph_degopt(degopt::Degopt; input = :A)
    return graph_degopt(degopt.x, degopt.y, input = input)
end

# Normally x::Vector{Tuple{Vector{Number},Vector{Number}}}
# Containing the coefficients in the recursion
function graph_degopt_B(x, T; input = :A)
    k = size(x, 1)

    graph = Compgraph(T)
    useful_syms = [:I, input]

    cref = []
    for s = 2:k+1
        base0 = "B$s"
        base_a = "Ba$s"
        base_b = "Bb$s"

        # First poly
        c_a = x[s-1][1]
        crefs_a = add_sum!(
            graph,
            Symbol(base_a),
            c_a,
            useful_syms,
            Symbol("$(base_a)_"),
        )

        append!(cref, crefs_a)

        # Second poly
        c_b = x[s-1][2]
        crefs_b = add_sum!(
            graph,
            Symbol(base_b),
            c_b,
            useful_syms,
            Symbol("$(base_b)_"),
        )
        append!(cref, crefs_b)

        # Multiply them together
        add_mult!(graph, Symbol(base0), Symbol(base_a), Symbol(base_b))
        push!(useful_syms, Symbol(base0))
        empty!(graph.outputs)
        push!(graph.outputs, Symbol(base0))
    end

    return (graph, cref)
end

"""
    order=get_topo_order_degopt(k)

A special implementation of [`get_topo_order`](@ref) for degree-optimal polynomials
generated with [`graph_degopt`](@ref). The natural order of computation is to compute
row by row.

See also [`get_degopt_crefs`](@ref).
"""
function get_topo_order_degopt(k)
    (x, z) = get_degopt_crefs(k)
    computation_order = Vector{Symbol}(undef, 0)
    for i = 1:k
        for n = 1:2
            for j = 2:(i+1)
                push!(computation_order, x[i][n][j][1])
            end
        end
        push!(computation_order, Symbol("B$(i+1)"))
    end
    for i = 2:(2+k)
        push!(computation_order, z[i][1])
    end
    return computation_order
end


text:
using BenchmarkTools
using DataFrames
using PartiallyTypedTables: partiallytyped
using Tables

# Monkey patch!
Tables.rows(df::DataFrames.AbstractDataFrame) = eachrow(df)

function filtered_sum(rows)
    s = 0.0
    for r in rows
        if r.a > 3
            s += r.b
        end
    end
    return s
end

SUITE = BenchmarkGroup()

n = 1_000_000
df = DataFrame(a = randn(n), b = randn(n))

SUITE["none"] = @benchmarkable filtered_sum($(eachrow(df)))
SUITE["partial"] = @benchmarkable filtered_sum($(partiallytyped(df, :a)))
SUITE["full"] = @benchmarkable filtered_sum($(Tables.rows(DataFrames.columntable(df))))


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Leptonica_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Leptonica")
JLLWrappers.@generate_main_file("Leptonica", UUID("6a1430e4-294a-53a5-a485-ec66ef6b843c"))
end  # module Leptonica_jll