text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "LibCURL"
version = v"7.64.1"

# Collection of sources required to build LibCURL
sources = [
    "https://curl.haxx.se/download/curl-7.64.1.tar.gz" =>
    "432d3f466644b9416bc5b649d344116a753aeaa520c8beaf024a90cba9d3d35d",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/curl-7.64.1

# Configure and build
./configure \
    --prefix=$prefix \
    --host=$target \
    --with-mbedtls \
    --without-ssl \
    --disable-manual

if [[ $target == *-w64-mingw32 ]]; then
    LDFLAGS="$LDFLAGS -L$prefix/bin"
elif [[ $target == x86_64-apple-darwin14 ]]; then
    LDFLAGS="$LDFLAGS -L$prefix/lib -Wl,-rpath,$prefix/lib"
else
    LDFLAGS="$LDFLAGS -L$prefix/lib -Wl,-rpath-link,$prefix/lib"
fi
make -j${nproc} LDFLAGS="$LDFLAGS" CPPFLAGS="$CPPFLAGS -I$prefix/include"
make install-exec
"""

# Build for ALL THE PLATFORMS!
platforms = supported_platforms()

# The products that we will ensure are always built
# Note: Avoid including the `curl` executable as it is unneeded and requires additional
# configuration to work.
products(prefix) = [
    LibraryProduct(prefix, "libcurl", :libcurl),
    # ExecutableProduct(prefix, "curl", :curl)
]

# Dependencies that must be installed before this package can be built
gh = "https://github.com"
dependencies = [
    "$gh/JuliaWeb/MbedTLSBuilder/releases/download/v0.20.0/build_MbedTLS.v2.6.1.jl",
    "$gh/bicycle1885/ZlibBuilder/releases/download/v1.0.4/build_Zlib.v1.2.11.jl",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
@testset "Q approximators" begin
    include("neural_network.jl")
    include("tabular.jl")
end

ext: jl


text:
# AbstractVec
#   - VecSeq: wrap
#   - VecMPI (TODO)
#   - VecGhost (TODO)
# for the MPI variants we won't be able to attach finalizers, as destroy needs to be called collectively.

const CVec = Ptr{Cvoid}

abstract type AbstractVec{T} <: AbstractVector{T} end
scalartype(::AbstractVec{T}) where {T} = T

# allows us to pass XXVec objects directly into CVec ccall signatures
Base.cconvert(::Type{CVec}, obj::AbstractVec) = obj.ptr
# allows us to pass XXVec objects directly into Ptr{CVec} ccall signatures
Base.unsafe_convert(::Type{Ptr{CVec}}, obj::AbstractVec) =
    convert(Ptr{CVec}, pointer_from_objref(obj))

"""
    VecSeq(v::Vector)

A standard, sequentially-stored serial PETSc vector, wrapping the Julia vector `v`.

This reuses the array `v` as storage, and so `v` should not be `resize!`-ed or otherwise have its length modified while the PETSc object exists.

This should only be need to be called for more advanced uses, for most simple usecases, users should be able to pass `Vector`s directly and have the wrapping performed automatically

# External Links
$(_doc_external("Vec/VecCreateSeqWithArray"))
"""
mutable struct VecSeq{T} <: AbstractVec{T}
    ptr::CVec
    array::Vector{T}
end

"""
    Vec(v::CVec)

Container for an abstract PETSc vector

# External Links
$(_doc_external("Vec/Vec"))
"""
mutable struct Vec{T} <: AbstractVec{T}
    ptr::CVec
end

Base.eltype(::Type{V}) where {V<:AbstractVec{T}} where T = T
Base.eltype(v::AbstractVec{T}) where {T} = T
Base.size(v::AbstractVec) = (length(v),)
Base.parent(v::AbstractVec) = v.array

# this allows setting V[1:2] = 3:4 on a PetscVec (more convenient)
function Base.setindex!(v::AbstractVec, val, I)
    v.array[I]=val
end
Base.getindex(v::AbstractVec, I) = v.array[I]

@for_libpetsc begin
    function VecSeq(comm::MPI.Comm, X::Vector{$PetscScalar}; blocksize=1)
        @assert initialized($petsclib)
        v = VecSeq(C_NULL, X)
        @chk ccall((:VecCreateSeqWithArray, $libpetsc), PetscErrorCode,
                (MPI.MPI_Comm, $PetscInt, $PetscInt, Ptr{$PetscScalar}, Ptr{CVec}),
                comm, blocksize, length(X), X, v)
        finalizer(destroy, v)
        return v
    end
    function destroy(v::AbstractVec{$PetscScalar})
        finalized($petsclib) ||
        @chk ccall((:VecDestroy, $libpetsc), PetscErrorCode, (Ptr{CVec},), v)
        return nothing
    end
    function Base.length(v::AbstractVec{$PetscScalar})
        r_sz = Ref{$PetscInt}()
        @chk ccall((:VecGetSize, $libpetsc), PetscErrorCode,
          (CVec, Ptr{$PetscInt}), v, r_sz)
        return r_sz[]
    end
    function LinearAlgebra.norm(v::AbstractVec{$PetscScalar}, normtype::NormType=NORM_2)
        r_val = Ref{$PetscReal}()
        @chk ccall((:VecNorm, $libpetsc), PetscErrorCode,
                   (CVec, NormType, Ptr{$PetscReal}),
                   v, normtype,r_val)
        return r_val[]
    end

    function assemblybegin(V::AbstractVec{$PetscScalar})
        @chk ccall((:VecAssemblyBegin, $libpetsc), PetscErrorCode, (CVec,), V)
        return nothing
    end
    function assemblyend(V::AbstractVec{$PetscScalar})
        @chk ccall((:VecAssemblyEnd, $libpetsc), PetscErrorCode, (CVec,), V)
        return nothing
    end
    function assemble(V::AbstractVec{$PetscScalar})
        assemblybegin(V)
        assemblyend(V)
    end

    function ownershiprange(vec::AbstractVec{$PetscScalar})
        r_lo = Ref{$PetscInt}()
        r_hi = Ref{$PetscInt}()
        @chk ccall((:VecGetOwnershipRange, $libpetsc), PetscErrorCode,
          (CVec, Ptr{$PetscInt}, Ptr{$PetscInt}), vec, r_lo, r_hi)
        r_lo[]:(r_hi[]-$PetscInt(1))
    end

    function view(vec::AbstractVec{$PetscScalar}, viewer::AbstractViewer{$PetscLib}=ViewerStdout($petsclib, getcomm(vec)))
        @chk ccall((:VecView, $libpetsc), PetscErrorCode,
                    (CVec, CPetscViewer),
                vec.ptr, viewer);
        return nothing
    end

    function localsize(cv::CVec)
        r_sz = Ref{$PetscInt}()
        @chk ccall((:VecGetLocalSize, $libpetsc), PetscErrorCode,
            (CVec, Ptr{$PetscInt}), cv, r_sz)
        return r_sz[]
    end

    function setvalues!(vec::AbstractVec{$PetscScalar},idxs,vals, insertmode::InsertMode)
        idxs = Vector(Int32,idxs);
        vals = Vector(vals);

        @chk ccall((:VecSetValues, $libpetsc), PetscErrorCode,
                 (CVec, $PetscInt, Ptr{$PetscInt}, Ptr{$PetscScalar},InsertMode), vec, length(idxs), idxs, vals, insertmode)
        return nothing

    end

    function unsafe_localarray(::Type{$PetscScalar}, cv::CVec; read::Bool=true, write::Bool=true)
        r_pv = Ref{Ptr{$PetscScalar}}()
        if write
            if read
                @chk ccall((:VecGetArray, $libpetsc), PetscErrorCode,
                    (CVec, Ptr{Ptr{$PetscScalar}}), cv, r_pv)
            else
                @chk ccall((:VecGetArrayWrite, $libpetsc), PetscErrorCode,
                    (CVec, Ptr{Ptr{$PetscScalar}}), cv, r_pv)
            end
        else
            @chk ccall((:VecGetArrayRead, $libpetsc), PetscErrorCode,
                (CVec, Ptr{Ptr{$PetscScalar}}), cv, r_pv)
        end
        r_sz = Ref{$PetscInt}()

        @chk ccall((:VecGetLocalSize, $libpetsc), PetscErrorCode,
            (CVec, Ptr{$PetscInt}), cv, r_sz)
        v = unsafe_wrap(Array, r_pv[], r_sz[]; own = false)

        if write
            if read
                finalizer(v) do v
                    @chk ccall((:VecRestoreArray, $libpetsc), PetscErrorCode, (CVec, Ptr{Ptr{$PetscScalar}}), cv, Ref(pointer(v)))
                    return nothing
                end
            else
                finalizer(v) do v
                    @chk ccall((:VecRestoreArrayWrite, $libpetsc), PetscErrorCode, (CVec, Ptr{Ptr{$PetscScalar}}), cv, Ref(pointer(v)))
                    return nothing
                end
            end
        else
            finalizer(v) do v
                @chk ccall((:VecRestoreArrayRead, $libpetsc), PetscErrorCode, (CVec, Ptr{Ptr{$PetscScalar}}), cv, Ref(pointer(v)))
                return nothing
            end
        end
        return v
    end

    function Base.fill!(v::AbstractVec{$PetscScalar}, x)
        @chk ccall((:VecSet, $libpetsc),
                   PetscErrorCode,
                   (CVec, $PetscScalar),
                   v, $PetscScalar(x))
        return v
    end

    function Base.setindex!(
        v::AbstractVec{$PetscScalar},
        val,
        i::Integer,
    )
        @chk ccall(
            (:VecSetValues, $libpetsc),
            PetscErrorCode,
            (CVec, $PetscInt, Ptr{$PetscInt}, Ptr{$PetscScalar}, InsertMode),
            v,
            1,
            Ref{$PetscInt}(i - 1),
            Ref{$PetscScalar}(val),
            INSERT_VALUES,
        )

        return val
    end

    function Base.getindex(
        v::AbstractVec{$PetscScalar},
        i::Integer,
    )
        vals = [$PetscScalar(0)]
        @chk ccall(
            (:VecGetValues, $libpetsc),
            PetscErrorCode,
            (CVec, $PetscInt, Ptr{$PetscInt}, Ptr{$PetscScalar}),
            v,
            1,
            Ref{$PetscInt}(i - 1),
            vals,
        )

        return vals[1]
    end
end

"""
    unsafe_localarray(PetscScalar, ptr:CVec; read=true, write=true)
    unsafe_localarray(ptr:AbstractVec; read=true, write=true)

Return an `Array{PetscScalar}` containing local portion of the PETSc data.

Use `read=false` if the array is write-only; `write=false` if read-only.

!!! note
    `Base.finalize` should be called on the `Array` before the data can be used.
"""
unsafe_localarray

unsafe_localarray(v::AbstractVec{T}; kwargs...) where {T} =
    unsafe_localarray(T, v.ptr; kwargs...)

"""
    map_unsafe_localarray!(f!, x::AbstractVec{T}; read=true, write=true)

Convert `x` to an `Array{T}` and apply the function `f!`.

Use `read=false` if the array is write-only; `write=false` if read-only.

# Examples
```julia-repl
julia> map_unsafe_localarray(x; write=true) do x
   @. x .*= 2
end

!!! note
    `Base.finalize` should is automatically called on the array.
"""
function map_unsafe_localarray!(f!, v::AbstractVec{T}; kwargs...) where {T}
    array = unsafe_localarray(T, v.ptr; kwargs...)
    f!(array)
    Base.finalize(array)
end

function Base.show(io::IO, ::MIME"text/plain", vec::AbstractVec)
    _show(io, vec)
end

VecSeq(X::Vector{T}; kwargs...) where {T} = VecSeq(MPI.COMM_SELF, X; kwargs...)
AbstractVec(X::AbstractVector) = VecSeq(X)

"""
    ownershiprange(vec::AbstractVec)

The range of indices owned by this processor, assuming that the vectors are laid out with the first n1 elements on the first processor, next n2 elements on the second, etc. For certain parallel layouts this range may not be well defined.

Note: unlike the C function, the range returned is inclusive (`idx_first:idx_last`)

# External Links
$(_doc_external("Vec/VecGetOwnershipRange"))
"""
ownershiprange

"""
    setvalues!(
        vector::AbstractVec{PetscScalar},
        indices::Vector{PetscInt},
        vals::Vector{PetscScalar},
        mode::InsertMode;
        num_vals = length(ind)
    )

Insert a set of values into the `vector`. Equivalent to one of the following
depending on the `mode`
```julia
vector[indices[1:num_vals]] .= vals[1:num_vals]
vector[indices[1:num_vals]] .+= vals[1:num_vals]
```

!!! warning
    `indices` should use 0-based indexing!

# External Links
$(_doc_external("Vec/VecSetValues"))
"""
function setvalues!(::AbstractVec) end

@for_libpetsc function setvalues!(
    vec::AbstractVec{$PetscScalar},
    inds::Vector{$PetscInt},
    vals::Vector{$PetscScalar},
    mode::InsertMode;
    num_vals = length(inds)
)
    @chk ccall(
         (:VecSetValues, $libpetsc),
         PetscErrorCode,
         (CVec, $PetscInt, Ptr{$PetscInt}, Ptr{$PetscScalar}, InsertMode),
         vec,
         num_vals,
         inds,
         vals,
         mode,
    )
    return vals
end

"""
    getvalues!(
        vector::AbstractVec{PetscScalar},
        indices::Vector{PetscInt},
        vals::Vector{PetscScalar};
        num_vals = length(inds)
    )

Get a set of values from the `vector`. Equivalent to one of the following
```julia
vals[1:num_vals] .= vector[indices[1:num_vals]]
```

!!! warning
    `indices` should use 0-based indexing!

# External Links
$(_doc_external("Vec/VecGetValues"))
"""
function getvalues!(::AbstractVec) end

@for_libpetsc function getvalues!(
    vec::AbstractVec{$PetscScalar},
    inds::Vector{$PetscInt},
    vals::Vector{$PetscScalar};
    num_vals = length(inds)
)
    @chk ccall(
         (:VecGetValues, $libpetsc),
         PetscErrorCode,
         (CVec, $PetscInt, Ptr{$PetscInt}, Ptr{$PetscScalar}),
         vec,
         num_vals,
         inds,
         vals,
    )
    return vals
end


text:
@enum(SocketStatus,
SOCKET_DONE,
SOCKET_NOT_READY,
SOCKET_PARTIAL,
SOCKET_DISCONNECTED,
SOCKET_ERROR)


text:
using Plots
plotlyjs()
include("FoggyMountainObj.jl")

n = 200
x = range(-20., stop=20., length=n)
y = x

##
p1 = surface(x,y,(x,y)->-FoggyMountainObj(x,y))
savefig("FoggySurface.png")
gui()

##
p2 = contour(x,y,(x,y)->-FoggyMountainObj(x,y))
gui()
savefig("FoggyContour.png")

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
module LSH

using Primes
using SparseArrays
using Statistics
using Base.Threads

import Base: match, show

export LSHModel
export match, signature

const DEFAULT_THRESHOLD = 0.75

struct LSHModel{T <: Integer}
    A::Vector{T}
    B::Vector{T}
    c::T
    termhashes::Array{T,2}
    dochashes::Array{T,2}
    bandhash::Dict{Float64,Any}
end

function LSHModel(nhashes::Integer, dtm::SparseMatrixCSC{T,T}) where {T <: Integer}
    ndocs,nterms = size(dtm)
    c = T(nextprime(nterms))
    A = rand(T(1):T(c), nhashes)
    B = rand(T(1):T(c), nhashes)

    termhashes = Array{T,2}(undef,nterms,nhashes)
    dochashes = Array{T,2}(undef,ndocs,nhashes)

    termhash!(termhashes, A, B, c)
    apply_term_weights!(termhashes, dtm)
    dochash!(dochashes, termhashes, dtm)

    LSHModel{T}(A, B, c, termhashes, dochashes, Dict{Float64,Any}())
end

function show(io::IO, model::LSHModel{T}) where T
    nterms,nhashes = size(model.termhashes)
    ndocs,_ = size(model.dochashes)
    print(io, "LSHModel{$T}(nhashes=$nhashes, ndocs=$ndocs, nterms=$nterms)")
end

# Weighted Minhash
function termhash!(termhashes::Array{T,2}, A::Vector{T}, B::Vector{T}, c::T) where {T <: Integer}
    @info("generating term hashes...")
    nterms,nhashes = size(termhashes)
    terms = 1:nterms
    for h in 1:nhashes
        a = A[h]
        b = B[h]
        termhashes[:,h] = (x->T((a*x+b)%c)).(terms)
    end
    nothing
end

function apply_term_weights!(termhashes::Array{T,2}, dtm::SparseMatrixCSC{T,T}) where {T <: Integer}
    @info("applying term weights...")
    ndocs,nterms = size(dtm)
    twf = [log(1+ndocs/sum(dtm[:,t])) for t in 1:nterms]
    mw = minimum(twf)
    termweights = [ceil(T, w/mw) for w in twf]
    termhashes .*= termweights
    nothing
end

function dochash!(dochashes::Array{T,2}, termhashes::Array{T,2}, dtm::SparseMatrixCSC{T,T}) where {T <: Integer}
    @info("generating signatures...")
    ndocs,nhashes = size(dochashes)

    Threads.@threads for h in 1:nhashes
        th = termhashes[:,h]
        for d in 1:ndocs
            tm = dtm[d,:]
            nzinds = SparseArrays.nonzeroinds(tm)
            dochashes[d,h] = minimum(th[nzinds] .* tm[nzinds])
        end
    end
    nothing
end

### LSH
similarity(sig1::Vector{T}, sig2::Vector{T}) where {T <: Integer} = mean(.==(sig1, sig2))
threshold(nbands::T, nrows::T) where {T <: Integer} = (1/nbands) ^ (1/nrows)

function factors(n::T) where {T <: Integer}
    allfactors = T[]

    for i in T(1):n
        if n % i == 0
            push!(allfactors, i)
        end
    end

    allfactors
end

# nbands and nrows are across hash algorithms
function best_partition(nhashes::T, threslim::Float64) where {T <: Integer}
    allfactors = factors(nhashes)
    mindiff = Inf
    bestpart = (T(0),T(0))

    for nbands in allfactors
        nrows = div(nhashes, nbands)
        t = threshold(nbands, nrows)
        tdiff = abs(t - threslim)
        if tdiff < mindiff
            mindiff = tdiff
            bestpart = (nbands, nrows)
        end
    end

    bestpart
end

function bandhash(model::LSHModel{T}, threshold::Float64) where {T <: Integer}
    get!(model.bandhash, threshold) do
        dochashes = model.dochashes
        ndocs,nhashes = size(dochashes)
        nbands,nrows = best_partition(T(nhashes), threshold)

        ht = Dict{Vector{T},Vector{T}}()

        for d in T(1):ndocs
            dh = dochashes[d,:]
            for b in 1:nbands
                brange = (1:nrows) .+ ((b-1)*nrows)
                bsig = dh[brange]
                bsigdocs = get!(ht, bsig) do
                    T[]
                end
                append!(bsigdocs, d)
            end
        end
        return ht
    end
end

match(model::LSHModel{T}, sig::Array{T,2}, threshold::Float64=DEFAULT_THRESHOLD) where {T <: Integer} = [match(model, sig[d,:], threshold) for d in 1:size(sig,1)]

function match(model::LSHModel{T}, sig::Vector{T}, threshold::Float64=DEFAULT_THRESHOLD) where {T <: Integer}
    nhashes = size(model.dochashes, 2)
    @assert(length(sig) == nhashes)

    nbands,nrows = best_partition(nhashes, threshold)
    bh = bandhash(model, threshold)
    similars = Set{Int}()

    for h in 1:nhashes
        for b in 1:nbands
            brange = (1:nrows) .+ ((b-1)*nrows)
            bsig = sig[brange]
            (x->push!(similars,x)).(get(bh, bsig, Int[]))
        end
    end
    result = []
    for d in similars
        sim = similarity(sig, model.dochashes[d,:])
        (sim < threshold) && continue
        push!(result, (similarity=sim, id=d))
    end
    sort!(result, by=x->x.similarity; rev=true)
end

function signature(model::LSHModel{T}, tm::SparseVector{T,T}) where {T <: Integer}
    termhashes = model.termhashes
    nterms,nhashes = size(termhashes)
    nzinds = SparseArrays.nonzeroinds(tm)
    T[minimum(termhashes[nzinds,h] .* tm[nzinds]) for h in 1:nhashes]
end

function signature(model::LSHModel{T}, tm::SparseMatrixCSC{T,T}) where {T <: Integer}
    nterms1,nhashes = size(model.termhashes)
    ndocs,nterms2 = size(tm)
    @assert(nterms1 == nterms2)
    dochashes = Array{T,2}(undef, ndocs, nhashes)

    for d in 1:ndocs
        dochashes[d,:] = signature(model, tm[d,:])
    end

    dochashes
end

function selftest(model::LSHModel{T}) where {T <: Integer}
    @info("running selftest...")
    ndocs,nhashes = size(model.dochashes)
    nfailures = 0
    for d in 1:ndocs
        sig = model.dochashes[d,:]
        matches = match(model, sig)
        (d in [x.id for x in matches]) || (nfailures += 1; @warn("docid did not match", docid=d))

        if (length(matches) > 1) && (d != first(matches).id)
            # check if there are multiple matches with same score
            topscore = first(matches).similarity
            docscore = matches[findfirst(x->x.id == d, matches)].similarity
            (docscore == topscore) || (nfailures += 1; @warn("doc not best batch for its own signature", docid=d, matches=matches))
        end
    end
    @info("selftest done", nfailures)
    nfailures
end

end # module LSH


text:
using Distributed

addprocs(4)

using JuliaDBMeta
using JuliaDB, Dagger
using IndexedTables: transform
using Test

iris1 = loadtable(joinpath(@__DIR__, "tables", "iris.csv"), distributed = false)
iris2 = table(iris1, chunks = 5)

@testset "utils" begin
    @test JuliaDBMeta.isquotenode(QuoteNode(:test))
    @test !JuliaDBMeta.isquotenode(Expr(:call, exp, 3))
    @test !JuliaDBMeta.isquotenode(3)
    @test !JuliaDBMeta.ispair(2)
    @test JuliaDBMeta.ispair(:(1 => 3))
    @test !JuliaDBMeta.istuple(2)
    @test JuliaDBMeta.istuple(:((1,2,3)))
end

@testset "with" begin
    t = table([1,2,3], [4,5,6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    @test (@with t :x .+ :y) == column(t, :x) .+ column(t, :y)
    @test @with(:x .+ :y)(t) == column(t, :x) .+ column(t, :y)

    f(t) = @with(t, :x .+ :z)
    s = @inferred f(t)
    @test s == [1.1, 2.2, 3.3]
    @with t :x .= :x + :y .+ 1
    @test column(t, :x) == [6,8,10]
    @test @with(t, length(_))== length(t)
    @test @with t Base.string.(:y) == ["4", "5", "6"]
    v = @with t begin
        l = :x
        {l}
    end
    @test v == (l = [6,8,10],)
    c = :x
    @with(t, cols(c)) == [6,8,10]
end

@testset "byrow" begin
    t = table([1, 2, 3], [4, 5, 6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    @test (@map t :x + :y) == column(t, :x) .+ column(t, :y)
    @test @map(:x + :y)(t) == column(t, :x) .+ column(t, :y)
    s = @map(:x + :z)
    @test s(t) == [1.1, 2.2, 3.3]
    @test @map(t, {:z}) == select(t, (:z,))
    @test @map(t, :x + :x) == [2, 4, 6]
    @test @map(t, _.y) == @map(t, :y)
    @test @map(t, :x + ^(:s isa Symbol ? 1 : 0)) == [2, 3, 4]
    c = :x
    @test @map(t, cols(c) + ^(:s isa Symbol ? 1 : 0)) == [2, 3, 4]

    @byrow! t :x = :x + :y + 1
    @test t == table([6, 8, 10], [4, 5, 6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    @byrow!(:x = 1)(t)
    @test column(t, :x) == [1,1,1]
    f! = @byrow!(:x = :x + :y + 1)
    @inferred f!(t)
    t = table([1, 2, 3], [4, 5, 6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    c = :x
    @byrow! t cols(c) = _.x + _.y + 1
    @test t == table([6, 8, 10], [4, 5, 6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    t1 = table([6, 8, 10], [4, 5, 6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    t2 = table(t1, chunks = 2)
    @test collect(@byrow! t2 :x = 1) == @byrow! t1 :x = 1
end

@testset "filter" begin
    t = table([1, 2, 3], [4, 5, 6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    @test (@filter t 4*:x == :y) == t[1:1]
    @test (@filter :z < 0.25)(t) == t[1:2]
    @test (@filter :z < 0.25)(t) == @filter t _.z < 0.25
end

@testset "transform" begin
    t = table([1,2,3], [4,5,6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    @test transform(t, :x => [2,3,4], :y => [4,5,6]) ==
        JuliaDBMeta._transform(t, :x => [2,3,4], :y => [4,5,6]) ==
        JuliaDBMeta._transform(t, Columns(x = [2,3,4], y = [4,5,6]))

    @test (@transform_vec t (a = :x .+ :y,)) == transform(t, :a => [1,2,3] .+ [4,5,6])
    @test @transform_vec((a = :x .+ :y,))(t) == transform(t, :a => [1,2,3] .+ [4,5,6])
    @test (@transform_vec t (z = :x .+ :y,)) == transform(t, :z => [1,2,3] .+ [4,5,6])

    @test (@transform t (a = :x .+ :y,))  == transform(t, :a => [1,2,3] .+ [4,5,6])
    @test @transform((a = :x .+ :y,))(t)  == transform(t, :a => [1,2,3] .+ [4,5,6])
    @test (@transform t (z = :x + :y,))  == transform(t, :z => [1,2,3] .+ [4,5,6])
    @test collect(@transform table(t, chunks = 2) (z = :x + :y,)) == (@transform t (z = :x + :y,))
end

@testset "where" begin
    t = table([1,2,3], [4,5,6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    @test (@where_vec t (:x .< 3) .& (:z .== 0.2)) == view(t, [2])
    @test @where_vec(t, 1:2) == view(t, 1:2)
    @test @where_vec(rows(t), 1:2) == view(t, 1:2)
    @test @where_vec((:x .< 3) .& (:z .== 0.2))(t) == view(t, [2])
    @test (@where t (:x < 3) .& (:z == 0.2)) == view(t, [2])
    @test @where((:x < 3) .& (:z == 0.2))(t) == view(t, [2])

    t = table([1,1,3], [4,5,6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    grp = groupby(@map((z = :z,))∘@where(:y != 5), t, :x, flatten = true)
    @test grp == table([1, 3], [0.1, 0.3], names = [:x, :z], pkey = :x)
    @test collect(@where iris2 :SepalLength > 4) == @where iris1 :SepalLength > 4
end

@testset "pair" begin
    t = table([1,2,3], [4,5,6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    @test select(t, (:x, :y) => i -> i.x+i.y) == select(t, @=>(:x+:y)) == @select(t, :x+:y) == @select(:x+:y)(t)
    @test select(t, (:x,) => i -> i.x+i.x) == select(t, @=>(:x+:x)) == @select(t, :x+:x)
    @test select(t, @=>(:x, :x+:y)) == select(t, @=>((:x, :x+:y))) == table([1,2,3], [5,7,9], names = [:x, Symbol("x + y")]) == @select(t, (:x, :x+:y))
    @test select(t, @=>(:x, :a => :x+:y)) == table([1,2,3], [5,7,9], names = [:x, :a]) == @select(t, (:x, :a => :x+:y))
end

@testset "apply" begin
    t = table([1,2,3], [4,5,6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    s = t |> @where(:x >= 2) |> @transform((s = :x + :y,))
    expected = table([2, 3], [5, 6], [0.2, 0.3], [7, 9], names = [:x, :y, :z, :s])
    @test s == expected
    s2 = @apply t begin
        @where :x >= 2
        @transform {s = :x+:y}
    end
    @test s2 == expected
    s3 = @apply t begin
        @where :x >= 2
        @transform {s = :x+:y}
        map(i -> i.s, _)
    end
    @test s3 == [7, 9]
    s4 = @apply rows(t) begin
        @where :x >= 2
        @transform {s = :x+:y}
        map(i -> i.s, _)
    end
    @test s4 == [7, 9]
    @test @apply(sort(_, :y))(t) == sort(t, :y)
    @test @apply(t, sort(_, :y))  == sort(t, :y)
    @test @apply(t, sort) == sort(t)

    t = table([1,2,2], [4,5,6], [0.1, 0.2, 0.3], names = [:x, :y, :z])
    t1 = @apply t :x flatten=true begin
        partialsort(_, 1:1, by = i -> i.y)
        @map {:y}
    end
    @test t1 ==  table([1,2], [4,5], names = [:x, :y], pkey = :x)
end

@testset "applychunked" begin
    t1 = @apply iris1 begin
        @where :SepalLength != 5.1
        @transform {Ratio = :SepalLength / :SepalWidth}
    end
    t2 = @applychunked iris2 begin
        @where :SepalLength != 5.1
        @transform {Ratio = :SepalLength / :SepalWidth}
    end
    @test t1 == collect(t2)
end

@testset "groupby" begin
    t = table([1,2,1,2], [4,5,6,7], [0.1, 0.2, 0.3,0.4], names = [:x, :y, :z])
    t1 = @groupby t :x {maximum(:y - :z)}
    @test t1 == table([1,2], [5.7,6.6], names = [:x, Symbol("maximum(y - z)")], pkey = :x)
    outcome = table([1,2], [5.7, 3.3], names = [:x, :m], pkey = :x)
    @test @groupby(t, :x, {m = maximum(:y - :z) / _.key.x}) == outcome
    c = :z
    @test @groupby(t, :x, {m = maximum(cols(:y) - cols(c)) / _.key.x}) == outcome
    @test @groupby(:x, {m = maximum(:y - :z) / _.key.x})(t) == outcome
    @test @groupby(reindex(t, :x), {m = maximum(:y - :z) / _.key.x}) == outcome
    @test @groupby({m = maximum(:y - :z) / _.key.x})(reindex(t, :x)) == outcome
    @test @groupby(t, :x, {l = length(_)}) == table([1,2], [2,2], names = [:x, :l], pkey = :l)
    @test @groupby(t, :x, {l = length(_)}) == t |> @groupby(:x, {l = length(_)})
    g = @groupby t :x {
                       l = length(_)

                      }
    @test g == @groupby(t, :x, {l = length(_)})
    @test @groupby(t, :x, flatten = true, _) == reindex(t, :x)
    @test @groupby(t, :x, {identity = _}) == groupby(identity, t, :x)

end


text:
import Base: length, step, rand, iterate, eltype, getindex, firstindex, lastindex
import StatsBase: var, mean
import Statistics: cov
import Distributions: invcov # remove this when you remove gaussianvectors.jl file.
using LinearAlgebra, SparseArrays, StatsBase, Distributions

# model.jl implements AbstractModel and SDEModel.
include("model.jl")
# samplepath.jl implements AbstractSamplePath and SamplePath and
# extends Base.step, Base.length, to the appropriate parametric
# types.
include("samplepath.jl")
# basisfunctions.jl implements the Fourier and Faber-Schauder basis.
include("basisfunctions.jl")
# SDE.jl implements AbstractSDE, SDE and extends
#Distribution.rand to sample an SDE sample path.
include("sde.jl")
# GaussianVectors.jl implements Gaussian vectors. Is depricated.
include("gaussianvectors.jl")
# gaussianprocess.jl implements the abstract type AbstractGaussianProcess with
# subtypes GaussianProcess and FaberSchauderExpansionWithGaussianCoefficients.
# Extends Base.length and Distributions.rand to AbstractGaussianProcess subtypes.
include("gaussianprocess.jl")

# Check of matrices die positief definiet moeten zijn ook daadwerkelijk positief
# definiet zijn.

# calculateposterior.jl implements calculateposterior methods to calculate the
# posterior.
include("calculateposterior.jl")


text:
import DataStructures: OrderedDict

import Iterators: chain
import Dates
import URIParser: URI, unescape, escape
import HttpCommon: parsequerystring

export headers,
    header,
    parseqs,
    parseqsr,
    params,
    fileinfo,
    isiter,
    RFC1123_datetime

function emptyfn(args...; kwargs...)
end

const str = string

RFC1123_datetime(t::Dates.DateTime) = Dates.format(t, Dates.RFC1123Format) * " GMT"
RFC1123_datetime() = RFC1123_datetime(Dates.now(Dates.UTC))

typealias Headers{K<:STR_TYPE, V<:Union(STR_TYPE, Array{STR_TYPE})} Associative{K, V}

headers(args...) =  OrderedDict{STR_TYPE,Union(STR_TYPE,Array{STR_TYPE})}(args...)
params(args...) = OrderedDict{STR_TYPE,Any}(args...)

function header(h::Headers, key::STR_TYPE, value::STR_TYPE)
    if haskey(h, key)
        if isa(h[key], STR_TYPE)
            h[key] = String[h[key], value]
        else
            push!(h[key], value)
        end
    else
        h[key] = value
    end
end

isiter(obj) = method_exists(Base.next, (typeof(obj), Integer))

function fileinfo(filename::STR_TYPE)
    if isfile(filename)
        (_, ext) = splitext(filename)

        mime = "application/octet-stream"

        if length(ext) > 1
            ext = ext[2:end]
            if haskey(mimetypes,ext)
                mime = mimetypes[ext]
            end
        end

        return ext, mime, filesize(filename)
    end
    return N
end

function parseqs(query::STR_TYPE)
    if !('=' in query)
        return throw("Not a valid query string: $query, must contain at least one key=value pair.")
    end

    q = Dict{STR_TYPE,Any}()

    for set in split(query, "&")
        key, val = split(set, "=")
        val = unescape(val)
        key = unescape(key)

        if haskey(q, key)
            push!(q[key], val)
        else
            q[key] = Any[val]
        end
    end

    q
end

function parseqsr{T<:Any}(allq::Associative{STR_TYPE, T})
    akeys = reverse(sort(collect(keys(allq))))

    q = OrderedDict{STR_TYPE,Any}()

    for key in akeys
        skeys = split(key, ".")
        d = q
        pkeys = skeys[1:end-1]

        if !isempty(pkeys)
            for (i, k) in enumerate(pkeys)
                if !haskey(d, k)
                    d[k] = OrderedDict{STR_TYPE,Any}()
                end
                d = d[k]
            end
        end

        k = skeys[end]
        d[k] = allq[key]
    end
    q
end

parseqsr(query::STR_TYPE) = parseqsr(parseqs(query))

import Base.mktemp

if !method_exists(mktemp, (STR_TYPE,))
    # Create and return the name of a temporary file along with an IOStream
    function mktemp(parent::STR_TYPE)
        b = joinpath(parent, "tmpXXXXXX")
        p = ccall(:mkstemp, Int32, (Ptr{UInt8},), b) # modifies b
        systemerror(:mktemp, p == -1)
        return (b, fdio(p, true))
    end
end

if !method_exists(mktemp, (Function,))
    function mktemp(fn::Function, parent=tempdir())
        (tmp_path, tmp_io) = mktemp(parent)
        try
            fn(tmp_path, tmp_io)
        finally
            close(tmp_io)
            rm(tmp_path)
        end
    end
end

export Event, add, pop, target, targets, listener, trigger, Events

type Event
    name::STR_TYPE
    target
    args
    kwargs
    stop::Bool
    data
end

Event(name::STR_TYPE; args=(), kwargs...) = Event(name, N, args, kwargs, false, Dict())

add(e::Event, target) = (e.target = (e.target, target); e)
pop(e::Event) = (e.target = e.target[1]; e)
target(e::Event) = e.target[2]

function targets(e::Event)
    local targets = {}
    t = e.target
    while t != N
        push!(t)
        t = t[1]
    end
    targets
end

immutable type Events
    listeners
end

Events(args::(Union(STR_TYPE,Symbol), Union(Array{Function}, Function))...; kwargs...) = begin
    evs = Events(Dict{STR_TYPE, Array{Function}}())

    for (key, cb) in chain(args, kwargs)
        if isa(cb, Function)
            listener(evs, key, cb)
        else
            for c in cb
                listener(evs, key, c)
            end
        end
    end
    evs
end

function listener(events::Events, key::Union(Symbol,STR_TYPE), cb::Function)
    key = string(key)
    l = get(events.listeners, key, N)
    if l == N
        l = Function[]
        events.listeners[key] = l
    end
    push!(l, cb)
    events
end

function trigger(events::Events, event::Event, target)
    add(event, target)
    try
        if haskey(events.listeners, event.name)
            for cb in events.listeners[event.name]
                cb(event, event.args...; event.kwargs...)
                if event.stop
                    break
                end
            end
        end
        event.data
    catch e
        event.data[:targets] = targets(event)
        rethrow()
    finally
        pop(event)
    end
    event.data
end


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

mutable struct PolicyDefinitionGroup <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    displayName::Any # spec type: Union{ Nothing, String } # spec name: displayName
    category::Any # spec type: Union{ Nothing, String } # spec name: category
    description::Any # spec type: Union{ Nothing, String } # spec name: description
    additionalMetadataId::Any # spec type: Union{ Nothing, String } # spec name: additionalMetadataId

    function PolicyDefinitionGroup(;name=nothing, displayName=nothing, category=nothing, description=nothing, additionalMetadataId=nothing)
        o = new()
        validate_property(PolicyDefinitionGroup, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(PolicyDefinitionGroup, Symbol("displayName"), displayName)
        setfield!(o, Symbol("displayName"), displayName)
        validate_property(PolicyDefinitionGroup, Symbol("category"), category)
        setfield!(o, Symbol("category"), category)
        validate_property(PolicyDefinitionGroup, Symbol("description"), description)
        setfield!(o, Symbol("description"), description)
        validate_property(PolicyDefinitionGroup, Symbol("additionalMetadataId"), additionalMetadataId)
        setfield!(o, Symbol("additionalMetadataId"), additionalMetadataId)
        o
    end
end # type PolicyDefinitionGroup

const _property_map_PolicyDefinitionGroup = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("displayName")=>Symbol("displayName"), Symbol("category")=>Symbol("category"), Symbol("description")=>Symbol("description"), Symbol("additionalMetadataId")=>Symbol("additionalMetadataId"))
const _property_types_PolicyDefinitionGroup = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("displayName")=>"String", Symbol("category")=>"String", Symbol("description")=>"String", Symbol("additionalMetadataId")=>"String")
Base.propertynames(::Type{ PolicyDefinitionGroup }) = collect(keys(_property_map_PolicyDefinitionGroup))
Swagger.property_type(::Type{ PolicyDefinitionGroup }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PolicyDefinitionGroup[name]))}
Swagger.field_name(::Type{ PolicyDefinitionGroup }, property_name::Symbol) =  _property_map_PolicyDefinitionGroup[property_name]

function check_required(o::PolicyDefinitionGroup)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

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


text:
using BenchmarkTools
using Logging
using TSXPlaygroundBenchmarks: BenchRandomIncrements
using TerminalLoggers

global_logger(TerminalLogger())

let path = get(ARGS, 1, joinpath(@__DIR__, "build", "results.json"))
    mkpath(dirname(path))
    open(versioninfo, joinpath(dirname(path), "versioninfo.txt"), write = true)
    @info "Start benchmarks. Result will be saved to: $path"
    @info "Loading benchmarks..."
    ntasks_list = 2 .^ (0:floor(Int, log2(Threads.nthreads())))
    if ntasks_list[end] !== Threads.nthreads()
        push!(ntasks_list, Threads.nthreads())
    end
    suite = BenchRandomIncrements.setup(
        # ntasks_list = 1:Threads.nthreads(),
        ntasks_list = ntasks_list,
        # include_tsx_static = false,
        nbatches_list = [100],
        ncounters_list = 2 .^ (15:2:25),
    )
    @info "Running benchmarks..."
    results = run(suite; verbose = true)
    @info "Saving benchmarks to $path..."
    BenchmarkTools.save(path, results)
    @info "Saving benchmarks to $path...DONE"
end


text:
using MacroTools
using Test

@testset "MacroTools" begin

include("match.jl")
include("split.jl")
include("destruct.jl")

end


text:
function computus(year)

    # Tracking lunar orbiting
    # Year's position on the 19 year lunar phase cycle
    # This tracks our location on the metonic cycle, which is a 19 year cycle
    # until lunar phases are on the exact same day as the start of the cycle
    a = mod(year, 19)

    # Century index
    k = floor(year/100)

    # shift of metonic cycle, add a day offset every 300 years
    # 8/25 is roughly 1/3rd, so every 300 years, that will add a day.
    # 13 is an offset.
    # 8 days off every 2500 years.
    # This might be why the alg. breaks down at 4200
    p = floor((13 + 8 * k) / 25)

    # Every 400 years, there are 97 leap days. k is already a century index,
    # so every 400 years, it will be 4. q will then subtract 1 from that and
    # give us 3, so that 100-3=97
    # (k-q) should be thought of as its own variable.
    q = floor(k / 4)

    # Correction to starting point of calculation each century, mod 30 because
    # the lunar calendar is about 30. This is a combination of previous terms.
    M = mod(15 - p + k - q, 30)

    # Finding the next Sunday
    # The difference in the number of leap days between the Gregorian and
    # Julian calendars. 4 = Number of leap days per year. 4 is an offset.
    # k-q is the number of additional leap days, and mod 7 because week
    N = mod(4 + k - q, 7)

    # Number of days from March 21st until the full moon
    # moon year = 12 lunar months = 354 days = 365 - 11 and mod(-11, 30) = 19
    # 0.25 year offset is handled by leap days later
    d = mod(19 * a + M, 30)

    # 52x7 = 364, but there are 365 days a year, so we are offsetting for this.
    # 2*b + 4*c  will drop by 1 every year, but 2 on leap years
    b = mod(year, 4)
    c = mod(year, 7)

    # Ignoring the 6*d term, this is the number of days until the next sunday
    # from March 22. The 6*d adds the appropriate offset to the next Sunday.
    e = mod(2 * b + 4 * c + 6 * d + N, 7)

    # Historical corrections for April 26 and 25 if we are on the last 8 years
    # of the metonic cycle. There is some ambiguity here I could not understand
    # related to the length of each lunar month being ~29.5 days, which could
    # mean that this is actually a mathematical offset.
    # Could be an equinox correction where the second half of the metonic cycle
    # has an earlier equinox?
    if (d == 29 && e == 6) || (d == 28 && e == 6 && a > 10)
        e = -1
    end

    # uncomment to recreate Servois's table
    # return mod(21 + d, 31)
    if(22+ d + e > 31)
        return string(d + e - 9)*" April"
    else
        return string(22+d+e)*" March"
    end
end

easter = computus(1800)


text:
struct InconsistentVRep{N, T, AT} <: VRepresentation{N, T}
    points::Polyhedra.PointsHull{N, T, AT}
    rays::Polyhedra.RaysHull{N, T, AT}
    function InconsistentVRep{N, T, AT}(points, lines, rays) where {N, T, AT}
        new{N, T, AT}(Polyhedra.PointsHull{N, T, AT}(points), Polyhedra.RaysHull(lines, rays))
    end
end
Polyhedra.dualtype(::Type{InconsistentVRep{N,T,AT}}, ::Type{AT}) where {N, T, AT} = Polyhedra.Intersection{N, T, AT}
Polyhedra.hvectortype(::Type{InconsistentVRep{N, T, AT}}) where {N, T, AT} = AT
Polyhedra.vvectortype(::Type{InconsistentVRep{N, T, AT}}) where {N, T, AT} = AT
Polyhedra.similar_type(PT::Type{<:InconsistentVRep}, d::FullDim{N}, ::Type{T}) where {N, T} = InconsistentVRep{N, T, Polyhedra.similar_type(Polyhedra.hvectortype(PT), d, T)}
Polyhedra.fulltype(::Type{InconsistentVRep{N, T, AT}}) where {N, T, AT} = InconsistentVRep{N, T, AT}
#Polyhedra.@subrepelem InconsistentVRep SymPoint points
Polyhedra.@subrepelem InconsistentVRep Point points
Polyhedra.@subrepelem InconsistentVRep Line rays
Polyhedra.@subrepelem InconsistentVRep Ray rays
@testset "Representation tests" begin
    @testset "MixMatRep with bad arguments" begin
        A = [1 1; -1 0; 0 -1]
        b = [1, 0, 0]
        linset = BitSet([1])

        @test_throws ErrorException hrep(A, [0, 0], linset)
        @test_throws ErrorException hrep(A, b, BitSet([4]))
        @test_throws ErrorException MixedMatHRep{3, Int}(A, b, BitSet())
        ine = hrep(A, b, linset)
        @test fulldim(ine) == 2
        @test (@inferred FullDim(ine)) == FullDim{2}()
        @test coefficienttype(ine) == Int
        @test translate(ine, [1, 0]).b == [2, -1, 0]

        V = [0 1; 1 0]
        @test_throws ErrorException MixedMatVRep{3, Int}(V, [1 0], BitSet())
        @test_throws ErrorException vrep(zeros(0, 2), [1 0]) # V-consistency
        @test_throws ErrorException vrep(V, [1 0 0], BitSet())
        @test_throws ErrorException vrep(V, [1 1], BitSet([2]))
        ext = vrep(V)
        @test fulldim(ext) == 2
        @test (@inferred FullDim(ine)) == FullDim{2}()
        @test MP.coefficienttype(ext) == Int
        @test translate(ext, [1, 0]).V == [1 1; 2 0]
    end

    @testset "LPHRepresentation with bad arguments" begin
        @test_throws DimensionMismatch LPHRepresentation(ones(2, 2), [1], [1], [1, 2], [1, 2])
        @test_throws DimensionMismatch LPHRepresentation(ones(2, 2), [1, 2], [1, 2], [1], [1])
        @test_throws DimensionMismatch LPHRepresentation{2, Int, Matrix{Int}}(ones(Int, 1, 1), [1], [1], [1], [1])
    end

    @testset "Lifted Representation with bad arguments" begin
        A = [1 -1 -1; 0 1 0; 0 0 1]
        ls = BitSet([1])

        @test_throws ErrorException LiftedHRepresentation(A, BitSet([4]))
        @test_throws ErrorException LiftedHRepresentation{3, Int}(A, ls)
        ine = copy(LiftedHRepresentation(A, ls))
        @test ine.A == A
        @test ine.A !== A
        #@test linset(ine) == ls
        @test ine.linset !== ls
        @test Polyhedra.similar_type(LiftedHRepresentation{2, Int, Matrix{Int}}, Float64) == LiftedHRepresentation{2, Float64, Matrix{Float64}}
        @test Polyhedra.similar_type(LiftedHRepresentation{2, Int, SparseMatrixCSC{Int, Int}}, FullDim{3}(), Float64) == LiftedHRepresentation{3, Float64, SparseMatrixCSC{Float64, Int}}

        A2 = [1 1; -1 0; 0 -1]
        b2 = [1, 0, 0]
        linset2 = BitSet([1])
        ine2 = hrep(A2, b2, linset2)

        ine = LiftedHRepresentation(ine2)
        @test ine.A == A
        @test ine.linset == ls

        V = [1 0 1; 1 1 0]
        Vlinset = BitSet(2)
        @test_throws ErrorException LiftedVRepresentation{3, Int}(V, Vlinset)
        @test_throws ErrorException LiftedVRepresentation(V, BitSet([4]))
        ext = copy(LiftedVRepresentation(V, Vlinset))
        @test ext.R == V
        @test ext.R !== V
        #@test linset(ext) == Vlinset
        @test ext.linset !== Vlinset

        @test Polyhedra.similar_type(LiftedVRepresentation{2, Int, SparseMatrixCSC{Int, Int}}, Float64) == LiftedVRepresentation{2, Float64, SparseMatrixCSC{Float64, Int}}
        @test Polyhedra.similar_type(LiftedVRepresentation{2, Int, Matrix{Int}}, FullDim{3}(), Float64) == LiftedVRepresentation{3, Float64, Matrix{Float64}}
    end

    @testset "eltype for some iterators is incorrect #7" begin
        function collecttest(it, exp_type)
            @test MP.coefficienttype(it) == exp_type
            a = collect(it)
            @test typeof(a) = Vector{exp_type}
        end
        hps = [HyperPlane([1, 2, 3], 7.)]
        shps = [@inferred HyperPlane((@SVector [1, 2, 3]), 7.)]
        @test eltype(shps) == HyperPlane{3, Float64, SVector{3, Float64}}
        hss = [HalfSpace([4, 5, 6.], 8)]
        shss = [@inferred HalfSpace((@SVector [4., 5., 6.]), 8)]
        @test eltype(shss) == HalfSpace{3, Float64, SVector{3, Float64}}
        for (hr, AT) in (((@inferred hrep(hps)), Vector{Float64}),
                         ((@inferred hrep(shps)), SVector{3, Float64}),
                         ((@inferred hrep(hss)), Vector{Float64}),
                         ((@inferred hrep(shss)), SVector{3, Float64}),
                         ((@inferred hrep(hps, hss)), Vector{Float64}),
                         ((@inferred hrep(shps, shss)), SVector{3, Float64}),
                         (hrep([1 2 3; 4 5 6], [7., 8], BitSet([1])), Vector{Float64}),
                         (hrep(spzeros(2, 3), [7., 8], BitSet([1])), SparseVector{Float64, Int}),
                         (SimpleHRepresentation([1 2 3; 4 5 6], [7., 8], BitSet([1])), Vector{Float64}),
                         (SimpleHRepresentation([1 2 3; 4 5 6], [7., 8]), Vector{Float64}))
            @test (@inferred coefficienttype(hr)) == Float64
            @test                                               (@inferred eltype(allhalfspaces(hr)))  == HalfSpace{3, Float64, AT}
            @test                                               (@inferred collect(allhalfspaces(hr))) isa Vector{HalfSpace{3, Float64, AT}}
            @test isempty(allhalfspaces(hr)) == iszero(nallhalfspaces(hr))
            @test (@inferred Polyhedra.halfspacetype(hr))    == (@inferred eltype(halfspaces(hr)))     == HalfSpace{3, Float64, AT}
            @test                                               (@inferred collect(halfspaces(hr)))    isa Vector{HalfSpace{3, Float64, AT}}
            @test isempty(halfspaces(hr)) == iszero(nhalfspaces(hr))
            @test (@inferred Polyhedra.hyperplanetype(hr))   == (@inferred eltype(hyperplanes(hr)))    == HyperPlane{3, Float64, AT}
            @test                                               (@inferred collect(hyperplanes(hr)))   isa Vector{HyperPlane{3, Float64, AT}}
            @test isempty(hyperplanes(hr)) == iszero(nhyperplanes(hr))
            @test_throws DimensionMismatch hr * ones(2, 3) # TODO replace by ones(2, 3) \ hr
        end
        #symps = [SymPoint([0, 1])]
        #ssymps = [SymPoint(@SVector [0, 1])]
        ps = [[1, 2], [3, 4]]
        sps = [(@SVector [1, 2]), (@SVector [3, 4])]
        rs = [Ray([0, 1])]
        srs = [Ray(@SVector [0, 1])]
        ls = [Line([0, 1])]
        sls = [Line(@SVector [0, 1])]
        @test eltype(sps) == SVector{2, Int}
        for (vr, AT) in (#((@inferred vrep(symps)), Vector{Int}),
                         #((@inferred vrep(ssymps)), SVector{2, Int}),
                         (vrep(ps), Vector{Int}),
                         ((@inferred vrep(sps)), SVector{2, Int}),
                         #((@inferred vrep(symps, ps)), Vector{Int}),
                         #(convexhull(symps..., ps...), Vector{Int}),
                         #(convexhull(ps..., symps...), Vector{Int}),
                         (convexhull(ps...), Vector{Int}),
                         #((@inferred vrep(ssymps, sps)), SVector{2, Int}),
                         #((@inferred convexhull(ssymps..., sps...)), SVector{2, Int}),
                         #((@inferred convexhull(sps..., ssymps...)), SVector{2, Int}),
                         ((@inferred convexhull(sps...)), SVector{2, Int}),
                         ((@inferred vrep(ls)), Vector{Int}),
                         ((@inferred vrep(sls)), SVector{2, Int}),
                         ((@inferred vrep(rs)), Vector{Int}),
                         ((@inferred vrep(srs)), SVector{2, Int}),
                         ((@inferred vrep(ls, rs)), Vector{Int}),
                         ((@inferred vrep(sls, srs)), SVector{2, Int}),
                         ((@inferred vrep(ps, ls, rs)), Vector{Int}),
                         ((@inferred vrep(sps, sls, srs)), SVector{2, Int}),
                         (vrep([1 2; 3 4]), Vector{Int}),
                         (vrep(spzeros(Int, 2, 2)), SparseVector{Int, Int}),
                         (SimpleVRepresentation([1 2; 3 4], zeros(Int, 0, 0), BitSet()), Vector{Int}),
                         (SimpleVRepresentation([1 2; 3 4], zeros(Int, 0, 0)), Vector{Int}),
                         (SimpleVRepresentation([1 2; 3 4]), Vector{Int}))
            @test (@inferred coefficienttype(vr)) == Int
#            @test (@inferred Polyhedra.sympointtype(vr)) == (@inferred eltype(sympoints(vr)))  == SymPoint{2, Int, AT}
#            @test                                           (@inferred collect(sympoints(vr))) isa Vector{SymPoint{2, Int, AT}}
#            @test isempty(sympoints(vr)) == iszero(nsympoints(vr))
            @test (@inferred Polyhedra.pointtype(vr))    == (@inferred eltype(points(vr)))     == AT
            @test                                           (@inferred collect(points(vr)))    isa Vector{AT}
            @test isempty(points(vr)) == iszero(npoints(vr))
            @test                                           (@inferred eltype(allrays(vr)))    == Ray{2, Int, AT}
            @test                                           (@inferred collect(allrays(vr)))   isa Vector{Ray{2, Int, AT}}
            @test isempty(allrays(vr)) == iszero(nallrays(vr))
            @test (@inferred Polyhedra.linetype(vr))     == (@inferred eltype(lines(vr)))      == Line{2, Int, AT}
            @test                                           (@inferred collect(lines(vr)))     isa Vector{Line{2, Int, AT}}
            @test isempty(lines(vr)) == iszero(nlines(vr))
            @test (@inferred Polyhedra.raytype(vr))      == (@inferred eltype(rays(vr)))       == Ray{2, Int, AT}
            @test                                           (@inferred collect(rays(vr)))      isa Vector{Ray{2, Int, AT}}
            @test isempty(rays(vr)) == iszero(nrays(vr))
            @test_throws DimensionMismatch ones(2, 1) * vr
        end
    end

    @testset "Iterating over halfspaces of a MixedMatHRep broken #9" begin
        A = [1 2; 3 4; 5 6]
        b = [1, 2, 3]
        halfspace = [1, 3]
        hyperplane = [2]
        linset = BitSet(2)
        hr = hrep(A, b, linset)
        Aall = [3 4; -3 -4; 1 2; 5 6]
        ball = [2, -2, 1, 3]
        for (i, h) in enumerate(allhalfspaces(hr))
            @test h.a == Aall[i, :]
            @test h.β == ball[i]
            @test isa(h, HalfSpace{2, Int})
        end
        for (i, h) in enumerate(halfspaces(hr))
            @test h.a == A[halfspace[i], :]
            @test h.β == b[halfspace[i]]
            @test isa(h, HalfSpace{2, Int})
        end
        for (i, h) in enumerate(hyperplanes(hr))
            @test h.a == A[hyperplane[i], :]
            @test h.β == b[hyperplane[i]]
            @test isa(h, HyperPlane{2, Int})
        end
    end

    @testset "Change FullDim" begin
        N = 5
        M = 10
        T = Int64
        reps = [MixedMatHRep{N, T, Matrix{T}}, MixedMatVRep{N, T, Matrix{T}}, LiftedHRepresentation{N, T, Matrix{T}}, LiftedVRepresentation{N, T, Matrix{T}}]
        for rep in reps
            changedrep = Polyhedra.similar_type(rep, FullDim{M}())
            @test fulldim(changedrep) == M
            @test (@inferred FullDim(changedrep)) == FullDim{M}()
            @test MP.coefficienttype(changedrep) == T
        end
    end

    @testset "Cartesian product" begin
        A = [1 2; 3 4; 5 6]
        a = [7, 8, 9]
        B = [10 11 12; 13 14 15]
        b = [16, 17]
        p1 = hrep(A, a, BitSet([2]))
        p2 = hrep(B, b, BitSet([1]))
        p = p1 * p2
        @test p.A == [A[2,:]' zeros(1, 3)
                      zeros(1, 2) B[1, :]'
                      A[[1,3],:] zeros(2, 3)
                      zeros(1, 2) B[2, :]']
        @test p.b == [a[2]; b[1]; a[[1,3]]; b[2]]
        @test p.linset == BitSet([1, 2])
    end

    @testset "isempty not working correctly for iterators #17" begin
        function vtest(vr, nr, np)
            hasr = nr > 0
            hasp = np > 0
            @test npoints(vr) == length(points(vr)) == np
            @test nallrays(vr) == length(allrays(vr)) == nr
            @test nlines(vr) == length(lines(vr)) == 0
            @test nrays(vr) == length(rays(vr)) == nr
            @test haspoints(vr) == !isempty(points(vr)) == hasp
            @test hasallrays(vr) == !isempty(allrays(vr)) == hasr
            @test haslines(vr) == !isempty(lines(vr)) == false
            @test hasrays(vr) == !isempty(rays(vr)) == hasr
        end
        vtest(vrep(zeros(0, 3)), 0, 0)
        vtest(vrep(zeros(1, 2)), 0, 1)
        vtest(vrep(zeros(1, 4), ones(2, 4)), 2, 1)
        vtest(vrep(zeros(2, 1), ones(1, 1)), 1, 2)
        vtest(LiftedVRepresentation(zeros(0, 2)), 0, 0)
        vtest(LiftedVRepresentation([1 0; 1 1]), 0, 2)
        vtest(LiftedVRepresentation([0 1; 0 2]), 2, 0)
        vtest(LiftedVRepresentation([0 0; 1 1]), 1, 1)
        function htest(hr, ne, ni)
            hase = ne > 0
            hasi = ni > 0
            @test nallhalfspaces(hr) == length(allhalfspaces(hr)) == 2ne + ni
            @test nhyperplanes(hr) == length(hyperplanes(hr)) == ne
            @test nhalfspaces(hr) == length(halfspaces(hr)) == ni
            @test hasallhalfspaces(hr) == !isempty(allhalfspaces(hr)) == hase || hasi
            @test hashyperplanes(hr) == !isempty(hyperplanes(hr)) == hase
            @test hashalfspaces(hr) == !isempty(halfspaces(hr)) == hasi
        end
        htest(hrep(ones(0, 3), zeros(0)), 0, 0)
        htest(hrep(ones(1, 2), zeros(1)), 0, 1)
        htest(hrep(ones(2, 4), zeros(2), BitSet(1:2)), 2, 0)
        htest(hrep(ones(3, 1), zeros(3), BitSet([2])), 1, 2)
        htest(LiftedHRepresentation(ones(0, 2)), 0, 0)
        htest(LiftedHRepresentation([0 1; 0 2]), 0, 2)
        htest(LiftedHRepresentation([0 1; 0 2], BitSet(1:2)), 2, 0)
        htest(LiftedHRepresentation([0 1; 0 2], BitSet([2])), 1, 1)
    end

    @testset "Building rep with different type" begin
        @test coefficienttype(MixedMatHRep{2, Float64}([1 2; 3 4], [1, 2], BitSet())) == Float64
        @test coefficienttype(MixedMatVRep{2, Float64}([1 2; 3 4], [1 2; 3 4], BitSet())) == Float64
        @test coefficienttype(LiftedHRepresentation{1, Float64}([1 2; 3 4])) == Float64
        @test coefficienttype(LiftedVRepresentation{1, Float64}([1 2; 3 4])) == Float64
    end

    @testset "Chebyshev center" begin
        p = hrep(eye(2), zeros(2))
        @test_throws ErrorException chebyshevcenter(p, lp_solver) # unbounded

        p = hrep([1 1; -1 -1], [0, -1])
        @test_throws ErrorException chebyshevcenter(p, lp_solver) # empty

        # examples/chebyshevcenter.ipynb
        A = [ 2  1
              2 -1
             -1  2
             -1 -2]
        b = ones(4)
        p = hrep(A, b)
        c, r = chebyshevcenter(p, lp_solver)
        @test c ≈ [0, 0] atol=1e-6
        @test r ≈ 0.4472135955 atol=1e-6

        p = convexhull([0, 0], [0, 1], [1, 0])
        @test_throws ErrorException chebyshevcenter(p) # Not yet implemented
    end

    @testset "V-consistency with iterator constructor" begin
        T = Int
        AT = Vector{Int}
        for VRepType in (Polyhedra.LiftedVRepresentation{2, T, Matrix{T}},
                         Polyhedra.MixedMatVRep{2, T, Matrix{T}},
                         Polyhedra.Hull{2, T, AT})
            @test_throws ErrorException VRepType(AT[], [Line([1, 2])])
            @test_throws ErrorException VRepType(AT[], Line{2, T, AT}[], [Ray([1, 2])])
            @test_throws ErrorException VRepType(AT[], [Line([1, 2])], [Ray([1, 2])])
            @test isempty(VRepType(AT[], Line{2, T, AT}[], Ray{2, T, AT}[]))
            @test isempty(VRepType(Line{2, T, AT}[], Ray{2, T, AT}[]))
            v = VRepType([Line([1, 2])])
            @test collect(points(v)) == [[0, 0]]
            @test collect(lines(v)) == [Line([1, 2])]
            @test !hasrays(v)
        end
        for vinc in (InconsistentVRep{2, T, AT}(AT[], Line{2, T, AT}[], [Ray([1, 2])]),
                     InconsistentVRep{2, T, AT}(AT[], [Line([1, 2])], Ray{2, T, AT}[]),
                     InconsistentVRep{2, T, AT}(AT[], [Line([1, 2])], [Ray([1, 2])]))
            @test_throws ErrorException Polyhedra.checkvconsistency(vinc)
            pinc = polyhedron(vinc)
            @test_throws ErrorException Polyhedra.checkvconsistency(pinc)
        end
    end
    @testset "Combination of different coefficient type" begin
        @testset "V-representation" begin
            generator_fulltest(convexhull([1, 0], Line([0, 1.])), convexhull(Line([0, 1]), [1, 0.]))
            @test conichull(convexhull([1, 0.]), conichull([0, 1])) isa Polyhedra.RaysHull{2, Float64}
            @test convexhull(convexhull([1, 0]), Line([0, 1.])) isa Polyhedra.Hull{2, Float64}
            @test convexhull(Line([0, 1.]), convexhull([1, 0])) isa Polyhedra.Hull{2, Float64}
            @test convexhull(convexhull(Line([0, 1.])), [1, 0]) isa Polyhedra.Hull{2, Float64}
            @test convexhull(convexhull(Line([1, 0])), Line([0, 1.])) isa Polyhedra.LinesHull{2, Float64}
            @test convexhull(conichull([1, 0.]), Line([0, 1])) isa Polyhedra.RaysHull{2, Float64}
            @test convexhull(conichull([1, 0.]), [0, 1]) isa Polyhedra.Hull{2, Float64}
        end
        @testset "H-representation" begin
            #FIXME inference fails @test (@inferred (HyperPlane([1, 1], 0) ∩ HyperPlane([1, 0], 1)) ∩ (HyperPlane([1, 1], 0) ∩ HyperPlane([1., 0.], 1))) isa Polyhedra.HyperPlanesIntersection{2, Float64}
            @test ((HyperPlane([1, 1], 0) ∩ HyperPlane([1, 0], 1)) ∩ (HyperPlane([1, 1], 0) ∩ HyperPlane([1., 0.], 1))) isa Polyhedra.HyperPlanesIntersection{2, Float64}
            @test HyperPlane([1, 1], 0) ∩ HalfSpace([1., 0.], 1.) isa Polyhedra.Intersection{2, Float64}
            #FIXME inference fails @test (@inferred (HalfSpace([1., 0.], 1.) ∩ (HyperPlane([1, 1], 0) ∩ HyperPlane([1, 0], 1)))) isa Polyhedra.Intersection{2, Float64}
            @test ((HalfSpace([1., 0.], 1.) ∩ (HyperPlane([1, 1], 0) ∩ HyperPlane([1, 0], 1)))) isa Polyhedra.Intersection{2, Float64}
        end
    end
    @testset "Conversion with different array type" begin
        @testset "V-representation" begin
            vv = convexhull(@SVector [0, 1]) + conichull((@SVector [1, 1]), Line(@SVector [1, 0]))
            mv = vrep([0 1], [1 1; 1 0], BitSet([2]))
            generator_fulltest(vv, mv)
            mvv = @inferred typeof(mv)(vv)
            vmv = @inferred typeof(vv)(mv)
            generator_fulltest(vv, mvv)
            generator_fulltest(vmv, vv)
        end
        @testset "H-representation" begin
            vh = HalfSpace((@SVector [1, 0]), 0) ∩ HyperPlane((@SVector [0, 1]), 1)
            mh = hrep([0 1; 1 0], [1, 0], BitSet(1))
            inequality_fulltest(vh, mh)
            mvh = @inferred typeof(mh)(vh)
            vmh = @inferred typeof(vh)(mh)
            inequality_fulltest(vh, mvh)
            inequality_fulltest(vmh, vh)
        end
    end
    @testset "Copy test" begin
        @testset "V-representation" begin
            vr = convexhull(@SVector [0, 1])
            vc = copy(vr)
            @test vc !== vr
            @test vc.points !== vr.points
            generator_fulltest(vc, vr)
            vr = conichull(Line(@SVector [1, 1]), Line(@SVector [1, 0]))
            vc = copy(vr)
            @test vc !== vr
            @test vc.lines !== vr.lines
            generator_fulltest(vc, vr)
            vr = conichull(Line(@SVector [1, 1]), Line(@SVector [1, 0]), @SVector [1, 1])
            vc = copy(vr)
            @test vc !== vr
            @test vc.rays !== vr.rays
            @test vc.lines !== vr.lines
            @test vc.lines.lines !== vr.lines.lines
            generator_fulltest(vc, vr)
            vr = convexhull(@SVector [0, 1]) + conichull(Line(@SVector [1, 1]), (@SVector [1, 1]), Line(@SVector [1, 0]))
            vc = copy(vr)
            @test vc !== vr
            @test vc.points !== vr.points
            @test vc.points.points !== vr.points.points
            @test vc.rays !== vr.rays
            @test vc.rays.lines !== vr.rays.lines
            @test vc.rays.lines.lines !== vr.rays.lines.lines
            @test vc.rays.rays !== vr.rays.rays
            generator_fulltest(vc, vr)
        end
        @testset "H-representation" begin
            hr = HyperPlane([1, 0], 1) ∩ HyperPlane([0, 1], -1)
            hc = copy(hr)
            @test hc !== hr
            @test hc.hyperplanes !== hr.hyperplanes
            inequality_fulltest(hc, hr)
            hr = HyperPlane([1, 0], 1) ∩ HyperPlane([1, 1], 0) ∩ HalfSpace([0, 1], -1)
            hc = copy(hr)
            @test hc !== hr
            @test hc.hyperplanes !== hr.hyperplanes
            @test hc.hyperplanes.hyperplanes !== hr.hyperplanes.hyperplanes
            @test hc.halfspaces !== hr.halfspaces
            inequality_fulltest(hc, hr)
        end
    end
    @testset "Preserving sparsity" begin
        for h in (HalfSpace(sparsevec([1], [1], 2), 1) ∩ HyperPlane(sparsevec([2], [-1]), 3), hrep(sparse([1, 2], [1, 2], [1, -1]), [1, 3]))
            @test Polyhedra.Intersection(h) isa Polyhedra.Intersection{2, Int, SparseVector{Int, Int}}
            @test LPHRepresentation(h) isa LPHRepresentation{2, Int, SparseMatrixCSC{Int, Int}}
            @test MixedMatHRep(h) isa MixedMatHRep{2, Int, SparseMatrixCSC{Int, Int}}
            @test LiftedHRepresentation(h) isa LiftedHRepresentation{2, Int, SparseMatrixCSC{Int, Int}}
            @test Polyhedra.Intersection{2, Float64}(h) isa Polyhedra.Intersection{2, Float64, SparseVector{Float64, Int}}
            @test LPHRepresentation{2, Float64}(h) isa LPHRepresentation{2, Float64, SparseMatrixCSC{Float64, Int}}
            @test MixedMatHRep{2, Float64}(h) isa MixedMatHRep{2, Float64, SparseMatrixCSC{Float64, Int}}
            @test LiftedHRepresentation{2, Float64}(h) isa LiftedHRepresentation{2, Float64, SparseMatrixCSC{Float64, Int}}
        end
        for v in (convexhull(sparsevec([1], [1], 2), sparsevec([2], [-1])), vrep(sparse([1, 2], [1, 2], [1, -1])))
            @test Polyhedra.Hull(v) isa Polyhedra.Hull{2, Int, SparseVector{Int, Int}}
            @test MixedMatVRep(v) isa MixedMatVRep{2, Int, SparseMatrixCSC{Int, Int}}
            @test LiftedVRepresentation(v) isa LiftedVRepresentation{2, Int, SparseMatrixCSC{Int, Int}}
            @test Polyhedra.Hull{2, Float64}(v) isa Polyhedra.Hull{2, Float64, SparseVector{Float64, Int}}
            @test MixedMatVRep{2, Float64}(v) isa MixedMatVRep{2, Float64, SparseMatrixCSC{Float64, Int}}
            @test LiftedVRepresentation{2, Float64}(v) isa LiftedVRepresentation{2, Float64, SparseMatrixCSC{Float64, Int}}
        end
    end
end


text:
using Dash

index_string = "
<!DOCTYPE html>

        {%metas%}
        <title>{%title%}</title>
        {%favicon%}
        {%css%}

        {%app_entry%}

            {%config%}
            {%scripts%}
            {%renderer%}

"

app = dash(assets_folder = "./dash_docs/assets", meta_tags = [Dict(["name"=>"description", "content" => "Dash for Julia User Guide and Documentation. Dash is a framework for building analytical web apps in Julia, R, and Python."])], index_string = index_string, suppress_callback_exceptions = true)

app.title = "Dash for Julia User Guide and Documentation | Plotly"

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
const TRANSITION = [1.0 1.0; 0.0 1.0]
const POS_NOISE = 1.5
const NOISE = [POS_NOISE 0.0; 0.0 sqrt(2)*POS_NOISE]
const PREDICT_NOISE = [0.1 0.0; 0.0 0.3]

struct State
    p::Float64
    v::Float64
    Σ::Array{Float64, 2}
end
State(p, v) = State(p, v, NOISE)

function predict(s::State, noise=PREDICT_NOISE)
    p, v = TRANSITION * [s.p ; s.v]
    Σ = TRANSITION * s.Σ * TRANSITION' + noise
    State(p, v, Σ)
end

function update(predicted::State, observed::State)
    K = predicted.Σ / (predicted.Σ + observed.Σ)
    # x = [predicted.p ; predicted.v]
    # z = [observed.p; observed.v]
    # p, v = x + K * (z - x)
    p = observed.p
    v = observed.v
    Σ = predicted.Σ - K * predicted.Σ
    State(p, v, Σ)
end

function smoothed_velocity_estimate(old::State, new::State, λ)
    diff = new.p - old.p
    scale = λ
    #scale = 1.0 - 0.5 * (1 - λ) * (1 + 1.0 / sqrt(2π * old.Σ[2, 2]^2))
    return scale * diff + (1 - scale) * old.v
end

function pdf(s::State, p::Float64)
    σ = s.Σ[1, 1]
    exp(-((s.p - p)^2) / (2σ^2)) / (sqrt(2π * σ^2))
end


text:
import FractionalCalculus.FracDiffAlg

"""
Hadamard sense fractional derivative algorithms.
"""
abstract type Hadamard <: FracDiffAlg end

"""
@article{article,
author = {Yin, Chuntao and Li, Changpin and Bi, Qinsheng},
year = {2018},
month = {01},
pages = {},
title = {Approximating Hadamard Derivative and Fractional Differential Equation via the Finite Part Integral},
journal = {SSRN Electronic Journal},
doi = {10.2139/ssrn.3281675}
}
"""

"""
    fracdiff(f, α, a, x, h, Hadamard_LRect())

Using finite part integral of **left rectangular formula** to compute the fractional hadamard derivative.
"""
struct Hadamard_LRect <: Hadamard end

"""
    fracdiff(f, α, a, x, h, Hadamard_Trap())

Using finite part integral of **right rectangular formula** to compute the fractional hadamard derivative.
"""
struct Hadamard_RRect <: Hadamard end

"""
    fracdiff(f, α, a, x, h, Hadamard_Trap())

Using finite part integral of **trapezoidal formula** to compute the fractional hadamard derivative.
"""
struct Hadamard_Trap <: Hadamard end

################################################################
###                    Type defination done                  ###
################################################################

#=
Hadamard Left Rectangular computing algorithms
=#
function fracdiff(f, α, x₀, x, h, ::Hadamard_LRect)
    N = Int64((x-x₀)/h)
    result = zero(Float64)

    for i ∈ 0:N-1
        result += LRectCoeff(i, N, h, α, x₀)*f(x₀+i*h)
    end

    return result
end
function LRectCoeff(i, n, h, α, x₀)
    if 0 ≤ i ≤ n-2
        return 1/gamma(1-α)*((log((x₀+n*h)/(x₀+i*h)))^(-α) - (log((x₀+n*h)/(x₀+(i+1)*h)))^(-α))
    elseif i == n-1
        return 1/gamma(1-α)*(log((x₀+n*h)/(x₀+i*h)))^(-α)
    end
end

#=
Hadamard Right Rectangular computing algorithm
=#
function fracdiff(f, α, x₀, x, h, ::Hadamard_RRect)
    N = Int64((x-x₀)/h)
    result = zero(Float64)

    for i ∈ 0:N-1
        result += LRectCoeff(i, N, h, α, x₀)*f(x₀+(i+1)*h)
    end

    return result
end
function RRectCoeff(i, n, h, α, x₀)
    if i == 0
        return 1/2*LRectCoeff(i, n, h, α, x₀)
    elseif 1 ≤ i ≤ n-1
        return 1/2*(LRectCoeff(i-1, n, h, α, x₀) + LRectCoeff(i, n, h, α, x₀))
    elseif i == n
        return 1/2*LRectCoeff(i-1, n, h, α, x₀)
    end
end

#=
Hadamard trapezoidal computing algorithm
=#
function fracdiff(f, α, x₀, x, h, ::Hadamard_Trap)
    N=Int64((x-x₀)/h)

    result = zero(Float64)

    for i ∈ 0:N
        result += RRectCoeff(i, N, h, α, x₀)*f(x₀+i*h)
    end

    return result
end


text:
module Celeste

# submodules
include("Log.jl")

include("SensitiveFloats.jl")

include("Model.jl")
include("Transform.jl")
include("PSF.jl")
include("SDSSIO.jl")

include("DeterministicVI.jl")
include("StochasticVI.jl")
include("MCMC.jl")

include("Infer.jl")
include("ParallelRun.jl")

include("Stripe82Score.jl")

end # module


text:
#=
Part of this work is derived from astropy/photutils. The relevant derivations
are considered under a BSD 3-clause license. =#

"""
    CircularAperture(x, y, r)
    CircularAperture(position, r)

A circular aperture.

A circular aperture with radius `r`. `r` must be greater than or equal to 0.

# Examples
```jldoctest
julia> ap = CircularAperture(0, 0, 10)
CircularAperture(0, 0, r=10)
```
"""
struct CircularAperture{T <: Number} <: AbstractAperture
    x::T
    y::T
    r::T

    function CircularAperture(x::T, y::T, r::T) where T <: Number
        r ≥ 0 || error("Invalid radius r=$r. r must be greater than or equal to 0.")
        new{T}(x, y, r)
    end
end

CircularAperture(center::AbstractVector, r) = CircularAperture(center..., r)
CircularAperture(x, y, r) = CircularAperture(promote(x, y, r)...)

function Base.show(io::IO, c::CircularAperture)
    print(io, "CircularAperture($(c.x), $(c.y), r=$(c.r))")
end

function bbox(c::CircularAperture{T}) where T
    xmin = ceil(Int, c.x - c.r - 0.5)
    ymin = ceil(Int, c.y - c.r - 0.5)
    xmax = ceil(Int, c.x + c.r - 0.5)
    ymax = ceil(Int, c.y + c.r - 0.5)
    return (xmin, xmax, ymin, ymax)
end

function mask(c::CircularAperture; method = :exact)
    bounds = edges(c)
    ny, nx = size(c)
    return circular_overlap(bounds..., nx, ny, c.r, method = method)
end

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

"""
    CircularAnnulus(x, y, r_in, r_out)
    CircularAnnulus(position, r_in, r_out)

A circular annulus with inner radius `r_in` and outer radius `r_out`. 0 ≤ `r_in` ≤ `r_out`.

# Examples
```jldoctest
julia> ap = CircularAnnulus(0, 0, 5, 10)
CircularAnnulus(0, 0, r_in=5, r_out=10)
```
"""
struct CircularAnnulus{T <: Number} <: AbstractAperture
    x::T
    y::T
    r_in::T
    r_out::T

    function CircularAnnulus(x::T, y::T, r_in::T, r_out::T) where T <: Number
        0 ≤ r_in ≤ r_out || error("Invalid radii ($r_in, $r_out). r_out must be greater than r_in which must be greater than or equal to 0.")
        new{T}(x, y, r_in, r_out)
end
end

CircularAnnulus(center::AbstractVector, r_in, r_out) = CircularAnnulus(center..., r_in, r_out)
CircularAnnulus(x, y, r_in, r_out) = CircularAnnulus(promote(x, y, r_in, r_out)...)

function Base.show(io::IO, c::CircularAnnulus)
    print(io, "CircularAnnulus($(c.x), $(c.y), r_in=$(c.r_in), r_out=$(c.r_out))")
end

function bbox(c::CircularAnnulus)
    xmin = ceil(Int, c.x - c.r_out - 0.5)
    ymin = ceil(Int, c.y - c.r_out - 0.5)
    xmax = ceil(Int, c.x + c.r_out - 0.5)
    ymax = ceil(Int, c.y + c.r_out - 0.5)
    return (xmin, xmax, ymin, ymax)
end

function mask(c::CircularAnnulus; method = :exact)
    bounds = edges(c)
    ny, nx = size(c)
    out = circular_overlap(bounds..., nx, ny, c.r_out, method = method)
    out .-= circular_overlap(bounds..., nx, ny, c.r_in, method = method)
end


text:
module Videre

# reserved

end # module

max_stars_repo_name: Gnimuc/Videre

max_issues_repo_name: Gnimuc/Videre

max_forks_repo_name: Gnimuc/Videre


text:
using CompositionalNetworks
using ConstraintDomains
using Documenter

makedocs(;
    modules=[CompositionalNetworks, ConstraintDomains],
    authors="Jean-François Baffier",
    repo="https://github.com/JuliaConstraints/CompositionalNetworks.jl/blob/{commit}{path}#L{line}",
    sitename="CompositionalNetworks.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", nothing) == "true",
        canonical="https://JuliaConstraints.github.io/CompositionalNetworks.jl",
        assets = ["assets/favicon.ico"; "assets/github_buttons.js"; "assets/custom.css"],
    ),
    pages=[
        "Home" => "index.md",
        "ICNs" => "icn.md",
        "Layers" => [
            "Transformation" => "transformation.md",
            "Arithmetic" => "arithmetic.md",
            "Aggregation" => "aggregation.md",
            "Comparison" => "comparison.md",
        ],
        "Dependencies" => [
            "ConstraintDomains.jl" => "domain.md",
        ],
        "Library" => [
            "public.md",
            "internal.md",
        ],
    ],
)

deploydocs(;
    repo="github.com/JuliaConstraints/CompositionalNetworks.jl.git",
    devbranch="main",
)


text:
# CodeGenOptions
function getArgv0(x::CodeGenOptions)
    @check_ptrs x
    return unsafe_string(clang_CodeGenOptions_getArgv0(x))
end

function getCommandLineArgs(x::CodeGenOptions)
    @check_ptrs x
    n = clang_CodeGenOptions_getCommandLineArgsNum(x)
    files = Vector{Ptr{Cuchar}}(undef, n)
    clang_CodeGenOptions_getCommandLineArgs(x, files, n)
    return unsafe_string.(files)
end

function PrintStats(x::CodeGenOptions)
    @check_ptrs x
    return clang_CodeGenOptions_PrintStats(x)
end


text: abstract type AbstractCumulant end


text:
const BZ = 1.38065e-23 # Boltzmann
const PL = 6.62607e-34 # Planck
const PM = 8.8541878128e-12 # permittivity

const M_He = 6.65e-27
const M_Ne = 3.35e-26
const M_Ar = 6.63e-26
const M_Kr = 1.39e-25
const M_Xe = 2.18e-25

const D_He = 2.19e-10
const D_Ne = 2.60e-10
const D_Ar = 3.66e-10
const D_Kr = 4.20e-10
const D_Xe = 4.94e-10

const ω_He = 0.66
const ω_Ne = 0.66
const ω_Ar = 0.81
const ω_Kr = 0.8
const ω_Xe = 0.85


text:
function NNCostGradient(X,y,LayerStruct;λ=0)
    m=size(X,1)                  #Number of training samples
    n=size(X,2)                  #Number of Features
    k=size(LayerStruct,1)-1           #Number of additional layers ( It does not include the first layer)

    ΘNparams=ΘNparams=Array{Int64,1}(undef,k)

    for i=1:k
       ΘNparams[i]=(LayerStruct[i]+1)*LayerStruct[i+1]
    end

    CostFunction=function (θ)
                    Θ=Array{Array{Float64,2},1}(undef, k)   #Initialize the shape of the Θ Array of Arrays
                    a=Array{Array{Float64,2},1}(undef, k+1)
                    z=Array{Array{Float64,2},1}(undef, k)
                    a[1]=[ones(m) X]
                    r=zeros(k)
                    c=0
                    for i=1:k
                        Θ[i]=reshape(θ[c+1:c+ΘNparams[i]],(LayerStruct[i+1],LayerStruct[i]+1))
                        c=ΘNparams[i]
                        z[i]=a[i]*Θ[i]'
                        a[i+1]=Sigmoid(z[i])
                        a[i+1]=[ones(m) a[i+1]]
                        r[i]=sum(Θ[i].^2)
                    end

                    labels=l=Vector(1:LayerStruct[end])'
                    Y=(y.==labels).*1.0
        J=(-1/m)*sum(Y.*log.(a[end][:,2:end]).+(1 .-Y).*log.(1 .-a[end][:,2:end]))+((0.5*λ/m).*sum(r))
                    return J
    end

        Gradient = function (θ)
                            Θ=Array{Array{Float64,2},1}(undef, k)   #Initialize the shape of the Θ Array of Arrays
                    a=Array{Array{Float64,2},1}(undef, k+1)
                    z=Array{Array{Float64,2},1}(undef, k)
                    a[1]=[ones(m) X]
                    r=zeros(k)
                    c=0
                    for i=1:k
                        Θ[i]=reshape(θ[c+1:c+ΘNparams[i]],(LayerStruct[i+1],LayerStruct[i]+1))
                        c=ΘNparams[i]
                        z[i]=a[i]*Θ[i]'
                        a[i+1]=Sigmoid(z[i])
                        a[i+1]=[ones(m) a[i+1]]

                    end
                    δ=Array{Array{Float64,2},1}(undef, k)
                    ∇=Array{Array{Float64,2},1}(undef, k)
                    labels=l=Vector(1:LayerStruct[end])'
                    Y=(y.==labels).*1.0
                    δ[end]=a[end][:,2:end].-Y
                    ∇[end]=((δ[end]'*a[end-1])./m).+(λ./m).*[zeros(size(Θ[end],1)) Θ[end][:,2:end]]
                    for i=k-1:-1:1
                    δ[i]=(δ[i+1]*Θ[i+1]).*SigmoidGradient([ones(m) z[i]])
                    ∇[i]=((δ[i][:,2:end]'*a[i])./m).+(λ./m).*[zeros(size(Θ[i],1)) Θ[i][:,2:end]]
                     end

                    ∇f=∇[1][:]
                    for i=2:k

                        ∇f=vcat(∇f,∇[i][:])

                    end

        return ∇f

                end
    return CostFunction, Gradient,ΘNparams
end


text:
#!/usr/bin/env julia

const patterns = split("""
    *.1
    *.c
    *.cpp
    *.h
    *.inc
    *.jl
    *.lsp
    *.make
    *.md
    *.mk
    *.rst
    *.scm
    *.sh
    *.yml
    *Makefile
""")

const errors = Set{Tuple{String,Int,String}}()

for path in eachline(`git ls-files -- $patterns`)
    lineno = 0
    non_blank = 0

    file_err(msg) = push!(errors, (path, 0, msg))
    line_err(msg) = push!(errors, (path, lineno, msg))

    for line in eachline(path, keep=true)
        lineno += 1
        contains(line, '\r')   && file_err("non-UNIX line endings")
        contains(line, '\ua0') && line_err("non-breaking space")
        endswith(line, '\n')   || line_err("no trailing newline")
        line = chomp(line)
        endswith(line, r"\s")  && line_err("trailing whitespace")
        contains(line, r"\S")  && (non_blank = lineno)
    end
    non_blank < lineno         && line_err("trailing blank lines")
end

if isempty(errors)
    println(stderr, "Whitespace check found no issues.")
    exit(0)
else
    println(stderr, "Whitespace check found $(length(errors)) issues:")
    for (path, lineno, msg) in sort!(collect(errors))
        if lineno == 0
            println(stderr, "$path -- $msg")
        else
            println(stderr, "$path:$lineno -- $msg")
        end
    end
    exit(1)
end

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
module FEMfunctions

  using LinearAlgebra
  using SparseArrays

  using WriteVTK

  export onedBilinear
  export onedLinForm
  export onedMesh
  export onedQuadratureRule
  export onedShape

  export saveFEMasVTK

  export TriMesh_ElementAdjacency
  export TriMesh_Interpolate
  export TriMesh_ProjectDerivatives
  export TriMesh_PromoteL2Q
  export TriMesh_Search
  export TriMesh_SquareMesh

  export twodBilinear
  export twodLinForm
  export twodMassMatrix
  export twodQuadratureRule
  export twodShape

  export TetMesh_QuadratureRule

  export threedBilinear

  include("onedBilinear.jl")
  include("onedLinForm.jl")
  include("onedMesh.jl")
  include("onedQuadratureRule.jl")
  include("onedShape.jl")

  include("saveFEMasVTK.jl")

  include("TriMesh_ElementAdjacency.jl")
  include("TriMesh_Interpolate.jl")
  include("TriMesh_ProjectDerivatives.jl")
  include("TriMesh_PromoteL2Q.jl")
  include("TriMesh_Search.jl")
  include("TriMesh_SquareMesh.jl")

  include("twodBilinear.jl")
  include("twodLinForm.jl")
  include("twodMassMatrix.jl")
  include("twodQuadratureRule.jl")
  include("twodShape.jl")

  include("TetMesh_QuadratureRule.jl")

  include("threedBilinear.jl")

end


text:
include("header.jl")

#=
These Jupyter notebooks and associated code serve as a companion to the course.

These are not static notebooks. **You are encouraged to play with them: run with different parameter values, create other plots, etc.!** The whole point is to enable an understanding of fluid mechanics through interactiveness.

*The code is released under the [MIT license](https://opensource.org/licenses/MIT).*
<hr />
=#

#=
## Table of Notebooks

- [Introduction](0.1-WelcomeToJulia.ipynb)

#### Chapter 1 - Plotting, fields, and flow lines
=#

#src - [Flow lines](1.1-FlowLines1.ipynb)
#src - [Plotting flow fields](1.2-PlottingFlowFields.ipynb)
#src - [Vorticity and streamfunction](1.3-VorticityAndStreamfunction.ipynb)

#=
#### Chapter 2 - Viscous flows at Re << 1

-

#### Chapter 3 - Wall-bounded flows at Re >> 1
=#
# - [Rayleigh problem](3.0-RayleighProblem.ipynb)
#src - [Flat plate flow](3.1-FlatPlateFlow.ipynb)
# - [Blasius and Falkner-Skan](3.2-BlasiusAndFalknerSkan.ipynb)
#src - [Computational solution](3.3-ComputationalSolution.ipynb)


text:
using Pkg

Pkg.activate(".") # activate the Counter environment from Julia

import Counter: Domain, API, Infrastructure

using .Domain, .API, .Infrastructure

const FILE_COUNTER = "./test_counter.txt" # Julia gives a warning when the value is changed

item = create("ar", 1000)

add_to_file(FILE_COUNTER, [item])

result = read_from_file(FILE_COUNTER)

println(result)


text:
"""
	AcousticPropagation

Calculates the sound pressure field for a cylindrical slice of the ocean.

TODO:
2. Define mature plotting function.
	a. Change output type to own struct `Ray`.
	b. Extend base plotting function to plot with `Ray` input?
3. Test ray tracing on scenarios.
	a. Change computation condition to range instead of arc length.
	b. Create mature plotting function.
	c. Convert to structs? Maybe later, for 2.0.
4. Implement amplitude calculation.
5. Produce pressure field output.
6. Combine with sonar equations.
7. Test pressure and TL field on scenarios.
8. Create as formal Julia package 1.0 (or 0.1?) as proof of concept.

Questions:
* How to rename library
* How to do tests for varying parameters?
"""
module AcousticPropagation
using ForwardDiff
using DifferentialEquations
using LinearAlgebra
using Interpolations

export Ray

"""

"""
function InterpolatingFunction(rng, val)
	Itp = LinearInterpolation(rng, val, extrapolation_bc = Flat())
	return ItpFcn(r) = Itp(r)
end

"""

"""
function InterpolatingFunction(rng, dpt, val)
	Itp = LinearInterpolation((dpt, rng), val, extrapolation_bc = Flat())
	return ItpFcn(r, z) = Itp(z, r)
end

"""
	t_rfl::Vector = boundary_reflection(t_inc::Vector, t_bnd::Vector)

`t_inc`    tangent vector of incident ray
`t_bnd`    tangent vector of boundary
`t_rfl`    tangent vector of reflected ray

TODO:
* Reconcile the two versions via error checking.
"""
function boundary_reflection(t_inc::Vector, t_bnd::Vector)
	# works for parabolic boundary
	MyAngle(tng) = atan(tng[2]/tng[1])
	θ_inc = MyAngle(t_inc)
	θ_bnd = MyAngle(t_bnd)

	c = cos(θ_inc)/t_inc[1]

	θ_inc_flat = θ_inc - θ_bnd
	θ_rfl_flat = -θ_inc_flat
	θ_rfl = θ_rfl_flat + θ_bnd

	return [cos(θ_rfl), sin(θ_rfl)]/c
end
# function boundary_reflection(t_inc::Vector, t_bnd::Vector)
# 	# works for generic boundary
# 	n_bnd = [-t_bnd[2], t_bnd[1]]
# # 	t_rfl = t_inc - 2(t_inc ⋅ n_bnd)*n_bnd
# 	t_rfl = t_inc - 2LinearAlgebra.dot(t_inc, n_bnd)*n_bnd

# 	MyAngle(tng) = atand(tng[2]/tng[1])
# 	θ_inc = MyAngle(t_inc)
# 	θ_bnd = MyAngle(t_bnd)
# 	θ_rfl = MyAngle(t_rfl)

# 	# if !isapprox(θ_inc - θ_bnd, θ_bnd - θ_rfl; atol = 1e-6)
# 	# 	println(t_inc)
# 	# 	println(θ_inc)
# 	# 	println(t_bnd)
# 	# 	println(θ_bnd)
# 	# 	println(t_rfl)
# 	# 	println(θ_rfl)
# 	# 	println(θ_inc - θ_bnd)
# 	# 	println(θ_bnd - θ_rfl)
# 	# 	error("Angle calculated incorrectly.")
# 	# end

# 	return t_rfl
# end

"""
	Boundary(z)

`z(r)`    depth of boundary in terms of range `r`
"""
struct Boundary
	z::Function
	dz_dr::Function
	condition::Function
	affect!::Function
	function Boundary(z::Function)
		dz_dr(r) = ForwardDiff.derivative(z, r)
		condition(u, t, ray) = z(u[1]) - u[2]
		affect!(ray) = ray.u[3], ray.u[4] = boundary_reflection([ray.u[3], ray.u[4]], [1, dz_dr(ray.u[1])])
		return new(z, dz_dr, condition, affect!)
	end
end
function Boundary(rz::AbstractArray)
	r_ = [rng for rng ∈ rz[:, 1]]
	z_ = [dpt for dpt ∈ rz[:, 2]]
	zFcn = InterpolatingFunction(r_, z_)
	return Boundary(zFcn)
end
function Boundary(z::Real)
	zFcn(r) = z
	return Boundary(zFcn)
end

"""
	Medium(c, R)

`c(r, z)`    sound speed at range `r` and depth `z`
`R`          maximum range of slice
"""
struct Medium
	c::Function
	∂c_∂r::Function
	∂c_∂z::Function
	∂²c_∂r²::Function
	∂²c_∂r∂z::Function
	∂²c_∂z²::Function
	R::Real
	function Medium(c::Function, R::Real)
		c_(x) = c(x[1], x[2])
		∇c_(x) = ForwardDiff.gradient(c_, x)
		∇c(r, z) = ∇c_([r, z])
		∂c_∂r(r, z) = ∇c(r, z)[1]
		∂c_∂z(r, z) = ∇c(r, z)[2]

		∂c_∂r_(x) = ∂c_∂r(x[1], x[2])
		∇∂c_∂r_(x) = ForwardDiff.gradient(∂c_∂r_, x)
		∇∂c_∂r(r, z) = ∇∂c_∂r_([r, z])

		∂c_∂z_(x) = ∂c_∂z(x[1], x[2])
		∇∂c_∂z_(x) = ForwardDiff.gradient(∂c_∂r_, x)
		∇∂c_∂z(r, z) = ∇∂c_∂z_([r, z])

		∂²c_∂r²(r, z) = ∇∂c_∂r(r, z)[1]
		∂²c_∂r∂z(r, z) = ∇∂c_∂r(r, z)[2]
		∂²c_∂z²(r, z) = ∇∂c_∂z(r, z)[2]

		return new(c, ∂c_∂r, ∂c_∂z, ∂²c_∂r², ∂²c_∂r∂z, ∂²c_∂z², R)
	end
end
function Medium(c::AbstractArray, R::Real = c[end, 1])
	r_ = [rc for rc ∈ c[1, 2:end]]
	z_ = [zc for zc ∈ c[2:end, 1]]
	c_ = c[2:end, 2:end]

	cFcn = InterpolatingFunction(r_, z_, c_)
	return Medium(cFcn, R)
end
function Medium(z::AbstractVector, c::AbstractVector, R = z[end])
	cMat = vcat([0 0 R], hcat(z, c, c))
	return Medium(cMat, R)
end
function Medium(c::Real, R::Real)
	cFcn(r, z) = c
	return Medium(cFcn, R)
end

"""
	Position(r, z, θ)

`r`    range of entity
`z`    depth of entity
"""
struct Position
	r::Real
	z::Real
end

struct Signal
	f::Real
end

struct Source
	Pos::Position
	Sig::Signal
end

function acoustic_propagation_problem(
	θ₀::Real,
	Src::Source,
	Ocn::Medium,
	Bty::Boundary,
	Ati::Boundary)

	function eikonal!(du, u, p, s)
		r = u[1]
		z = u[2]
		ξ = u[3]
		ζ = u[4]
		τ = u[5]
		pʳ = u[6]
		pⁱ = u[7]
		qʳ = u[8]
		qⁱ = u[9]

		∂²c_∂n²(r, z) = Ocn.c(r, z)^2*(
			Ocn.∂²c_∂r²(r, z)*ζ^2
			- 2Ocn.∂²c_∂r∂z(r, z)*ξ*ζ
			+ Ocn.∂²c_∂z²(r, z)*ξ^2
		)

		du[1] = dr_ds = Ocn.c(r, z)*ξ
		du[2] = dz_ds = Ocn.c(r, z)*ζ
		du[3] = dξ_ds = -Ocn.∂c_∂r(r, z)/Ocn.c(r, z)^2
		du[4] = dζ_ds = -Ocn.∂c_∂z(r, z)/Ocn.c(r, z)^2
		du[5] = dτ_ds = 1/Ocn.c(r, z)
		du[6] = dpʳ_ds = ∂²c_∂n²(r, z)/Ocn.c(r, z)^2*qʳ
		du[7] = dpⁱ_ds = ∂²c_∂n²(r, z)/Ocn.c(r, z)^2*qⁱ
		du[8] = dqʳ_ds = Ocn.c(r, z)*pʳ
		du[9] = dqⁱ_ds = Ocn.c(r, z)*pⁱ
	end

	rng_condition(u, t, ray) = Ocn.R/2 - abs(u[1] - Ocn.R/2)
	rng_affect!(ray) = terminate!(ray)
	CbRng = ContinuousCallback(rng_condition, rng_affect!)
	CbBty = ContinuousCallback(Bty.condition, Bty.affect!)
	CbAti = ContinuousCallback(Ati.condition, Ati.affect!)
	CbBnd = CallbackSet(CbRng, CbBty, CbAti)

	r₀ = Src.Pos.r
	z₀ = Src.Pos.z
	ξ₀ = cos(θ₀)/Ocn.c(r₀, z₀)
	ζ₀ = sin(θ₀)/Ocn.c(r₀, z₀)
	τ₀ = 0.0

	λ₀ = Ocn.c(r₀, z₀)/Src.Sig.f
	ω = Src.Sig.f
	p₀ʳ = 1.0
	p₀ⁱ = 0.0
	W₀ = 100λ₀ # 10..50
	q₀ʳ = 0.0
	q₀ⁱ = ω*W₀^2/2

	u₀ = [r₀, z₀, ξ₀, ζ₀, τ₀, p₀ʳ, p₀ⁱ, q₀ʳ, q₀ⁱ]

	TLmax = 100
	S = 10^(TLmax/10)
	sSpan = (0., S)

	prob_eikonal = ODEProblem(eikonal!, u₀, sSpan)

	return prob_eikonal, CbBnd
end

function solve_acoustic_propagation(prob_eikonal, CbBnd)
	@time RaySol = solve(prob_eikonal, callback = CbBnd, reltol=1e-8, abstol=1e-8)
	return RaySol
end

struct Ray
	θ₀
	Sol
	S
	r
	z
	ξ
	ζ
	τ
	p
	q
	θ
	c
	function Ray(θ₀::Real, Src::Source, Ocn::Medium, Bty::Boundary; Ati::Boundary = Boundary(0))
		Prob, CbBnd = acoustic_propagation_problem(θ₀, Src, Ocn, Bty, Ati)
		Sol = solve_acoustic_propagation(Prob, CbBnd)

		S = Sol.t[end]
		r(s) = Sol(s, idxs = 1)
		z(s) = Sol(s, idxs = 2)
		ξ(s) = Sol(s, idxs = 3)
		ζ(s) = Sol(s, idxs = 4)
		τ(s) = Sol(s, idxs = 5)
		p(s) = Sol(s, idxs = 6) + im*Sol(s, idxs = 7)
		q(s) = Sol(s, idxs = 8) + im*Sol(s, idxs = 9)
		θ(s) = atan(ζ(s)/ξ(s))
		c(s) = cos(θ(s))/ξ(s)

		return new(θ₀, Sol, S, r, z, ξ, ζ, τ, p, q, θ, c)
	end
end

struct Beam
	θ₀::Real
	ray
	b::Function
	S::Real
	function Beam(θ₀::Real, Src::Source, Ocn::Medium, Bty::Boundary, Ati::Boundary = Boundary(0))

		ray = Ray(θ₀, Src, Ocn, Bty, Ati)

		r(s) = ray.r(s)
		z(s) = ray.z(s)
		τ(s) = ray.τ(s)
		p(s) = ray.p(s)
		q(s) = ray.q(s)
		c(s) = ray.c(s)
		W(s) = sqrt(-2/ω/imag(p(s)/q(s)))

		c₀ = c(0)
		ω = 2π*Src.Sig.f
		λ₀ = c₀/Src.Sig.f
		W₀ = W(0)
		q₀ = q(0)

		A = 1/c₀ * exp(im*π/4)*sqrt(q₀*ω*cos(θ₀)/2π)
		b_(s, n) = A * sqrt(c(s)/r(s)/q(s)) * exp(-im*ω * (τ(s) + p(s)/q(s)*n^2/2))
		b(s, n) = max(100, b_(s, n))

		return new(θ₀, ray, b, ray.S)
	end
end

# struct Field
# 	θ₀::Union{Real,Vector}
# 	p::Func
# 	function Field(θ₀, Rcv, Src, Ocn, Bty, Before, After)
# 		return new(θ₀, p, TL)
# 	end
# end

# struct FieldCoherent
# 	θ₀::Real
# 	function FieldCoherent(θ₀, Src::Position, Ocn::Medium, Bty::Boundary; Ati::Boundary = Boundary(0))
# 		BeforeSum(p) = p
# 		AfterSum(p) = p
# 		return Field(θ₀, Src, Ocn, Bty, BeforeSum, AfterSum)
# 	end
# end

Base.broadcastable(m::Position) = Ref(m)
Base.broadcastable(m::Medium) = Ref(m)
Base.broadcastable(m::Boundary) = Ref(m)
Base.broadcastable(m::Source) = Ref(m)

end


text:
@testset "GeometricPrimitives" begin

S = GRAYTR.Sphere()
mat = GRAYTR.Lambert(NoLight())

G = GeometricPrimitive(S, mat, nothing, 1)
@test G.id == 1
@test G.shape == S
@test G.material == mat

shapes = [Sphere(), Disk()]
prims = apply_material(shapes, mat)
@test typeof(prims) == Array{GeometricPrimitive, 1}
@test length(prims) == 2
@test typeof(prims[1].shape) == Sphere
@test typeof(prims[2].shape) == Disk

end # main testset

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
@testset "Banana" begin
    using MLUtilities
    dm = BananaModel()

    g = DataModel.getgrad(dm)
    l = DataModel.getllik(dm)
    for i in 1:10
        @test checkgrad(randn(2),l,g, eps = 1e-8) < 1e-4
    end
end


text:
using CSV
using DataFrames
using DelimitedFiles
using YAML
using Query

using Base

function share!(d::Dict, set::Dict; save = true, overwrite = false)

    # If there is already sharing data, read it and return.
    d_read = read_build("share"; save = save, overwrite = overwrite);
    if !isempty(d_read)
        [d[k] = v for (k,v) in d_read]
        set[:notrd] = setdiff(set[:s], d[:utd][:,:s])
        return d
    end

    # READ SHARING DATA.
    d_read = read_build("share_i"; save = save, overwrite = overwrite);
    if isempty(d_read)
        y = read_from(joinpath("src","readfiles","build","shareinp.yml"))
        d_read = Dict(k => edit_with(v) for (k,v) in y)
        d_read = Dict(k => sort(filter_with(df, set; extrapolate = true)) for (k, df) in d_read)

        write_build("share_i", d_read; save = save)
    end
    [d[k] = v for (k,v) in d_read]

    # Do the sharing.
    share_pce!(d)
    share_sgf!(d)
    share_utd!(d, set)
    share_region!(d, set)
    share_labor!(d, set)
    share_rpc!(d, set)

    d_save = Dict(k => d[k] for k in [:gsp,:labor,:pce,:rpc,:region,:sgf,:utd])
    write_build("share", d_save; save = save)

    return d
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
module JinjaTemplates

export LazyTemplateLoader, get_template_name, render, get_searchpaths, add_searchpath!

using PyCall
@pyimport jinja2

""" Object to store in memory PyObjects representing Jinja templates """
type LazyTemplateLoader
    jinjaenv
    templates

    function LazyTemplateLoader(paths::Array; kwargs...)
        env = jinja2.Environment(;
            loader=jinja2.FileSystemLoader(paths),
            kwargs...
        )
        return new(env, Dict())
    end
    LazyTemplateLoader(path::AbstractString; kwargs...) = LazyTemplateLoader([path]; kwargs...)
end

""" Returns a jinja template. Onload loads the template once requested """
function Base.getindex(loader::LazyTemplateLoader, key)
    try
        return loader.templates[key]
    catch x
        if isa(x, KeyError)
            loader.templates[key] = loader.jinjaenv[:get_template](key)
            return loader.templates[key]
        else
            throw(x)
        end
    end
end

"""
Returns an array of the curret template sources
"""
get_searchpaths(loader::LazyTemplateLoader) = loader.jinjaenv[:loader][:searchpath]

"""
Insert location into search path

`at` can be :start, :end, :index
index is a 1 based index representing the position
"""
function add_searchpath!(loader::LazyTemplateLoader, path::AbstractString; at=:start, index=1)
    searchpath = loader.jinjaenv[:loader]["searchpath"]
    if is(at, :start)
        searchpath[:insert](0, path)
    elseif is(at, :end)
        searchpath[:append](path)
    elseif is(at, :index)
        searchpath[:insert](index-1, path)
    end

    # Clear cache
    loader.templates=Dict()
    return loader
end

"""
Takes a variable and converts its type into a string fit for the filesystem
"""
function get_template_name(item)
    t = string(typeof(item))
    return replace(t, r"[^\w\.]", s"_")
end

""" Renders given kwargs into the specified template """
function render(loader::LazyTemplateLoader, template::AbstractString; kwargs...)
    template = loader[template]
    return template[:render](; kwargs...)
end

render(loader::LazyTemplateLoader, obj, extension; kwargs...) = render(loader, "$get_template_name(obj).$extension"; kwargs...)

end # module


text:
# dnabin.jl
# =========
#
# A nucleic acid type compatible with R's APE package.
#
# This file is a part of the BioJulia ecosystem.
# License is MIT: https://github.com/BioJulia/BioBridgeR.jl/blob/master/LICENSE.md

primitive type DNAbin <: NucleicAcid 8 end

const char_to_dnabin = [0x02 for _ in 0x00:0xff]
const dnabin_to_char = ['?' for _ in 0x01:0xed]

# derived from "The DDBJ/ENA/GenBank Feature Table Definition"
# §7.4.1 Nucleotide base code (IUPAC)
# http://www.insdc.org/documents/feature_table.html#7.4.1

for (char, doc, bits) in [
        ('-', "DNAbin Gap",                                   0x04),
        ('A', "DNAbin Adenine",                               0x88),
        ('C', "DNAbin Cytosine",                              0x28),
        ('G', "DNAbin Guanine",                               0x48),
        ('T', "DNAbin Thymine",                               0x18),
        ('M', "DNAbin Adenine or Cytosine",                   0xA0),
        ('R', "DNAbin Adenine or Guanine",                    0xC0),
        ('W', "DNAbin Adenine or Thymine",                    0x90),
        ('S', "DNAbin Cytosine or Guanine",                   0x60),
        ('Y', "DNAbin Cytosine or Thymine",                   0x30),
        ('K', "DNAbin Guanine or Thymine",                    0x50),
        ('V', "DNAbin Adenine, Cytosine or Guanine",          0xE0),
        ('H', "DNAbin Adenine, Cytosine or Thymine",          0xB0),
        ('D', "DNAbin Adenine, Guanine or Thymine",           0xD0),
        ('B', "DNAbin Cytosine, Guanine or Thymine",          0x70),
        ('N', "DNAbin Adenine, Cytosine, Guanine or Thymine", 0xF0)]
    var = Symbol("DNAbin_", char != '-' ? char : "Gap")
    @eval begin
        @doc $(doc) const $(var) = reinterpret(DNAbin, $(bits))
        char_to_dnabin[$(Int(char) + 1)] = char_to_dnabin[$(Int(lowercase(char)) + 1)] = $(bits)
        dnabin_to_char[$(Int(bits) - 3)] = $(char)
    end
end

# Conversion
# ----------

# Conversion between DNAbin and DNA.

@inline function _DNA_to_raw(nt::DNA)
    bits = UInt8(nt)
    db = swapbits(bits, 0, 3) << 4
    db |= ifelse(iscertain(nt), 0x08, 0x00)
    db |= ifelse(isgap(nt), 0x04, 0x00)
    return db
end

@inline function Base.convert(::Type{DNAbin}, nt::DNA)
    return DNAbin(_DNA_to_raw(nt))
end

@inline function _raw_to_DNA(x::UInt8)
    return swapbits(x, 4, 7) >> 4
end

@inline function Base.convert(::Type{DNA}, nt::DNAbin)
    return BioSymbols.DNA(_raw_to_DNA(UInt8(nt)))
end

# Conversion between DNAbin and characters.

@inline function Base.convert(::Type{DNAbin}, c::Char)
    if c > '\xff'
        throw(InexactError())
    end
    @inbounds nt = char_to_dnabin[Int(c) + 1]
    if (nt & 0x02) > 0
        throw(InexactError())
    end
    return DNAbin(nt)
end

@inline function Base.convert(::Type{Char}, nt::DNAbin)
    @inbounds return dnabin_to_char[Int(nt) - 3]
end

# Alphabet
# --------

@eval function BioSymbols.alphabet(::Type{DNAbin})
    return $(tuple([DNAbin(BioSymbols.DNA(x)) for x in 0b0000:0b1111]...))
end

# Show / print
# ------------

BioSymbols.prefix(::Type{DNAbin}) = "DNAbin"

# Basic operators
# ---------------

function Base.:~(x::DNAbin)
    return reinterpret(DNAbin, ~reinterpret(UInt8, x))
end

function Base.:+(x::DNAbin, y::Integer)
    return reinterpret(DNAbin, (UInt8(x) + y % UInt8))
end

@inline function BioSymbols.gap(::Type{DNAbin})
    return reinterpret(DNAbin,0x04)
end

@inline function BioSymbols.isGC(nt::DNAbin)
    return UInt8(nt) & 0x96 == 0x00
end

@inline function BioSymbols.ispurine(nt::DNAbin)
    return (UInt8(nt) & 0x36) == 0x00
end

@inline function BioSymbols.ispyrimidine(nt::DNAbin)
    return (UInt8(nt) & 0xC6) == 0x00
end

@inline function BioSymbols.isambiguous(nt::DNAbin)
    return (UInt8(nt) & 0x0C) == 0x00
end

@inline function BioSymbols.iscertain(nt::DNAbin)
    return (UInt8(nt) & 0x08) == 0x08
end

@inline function BioSymbols.isgap(nt::DNAbin)
    return nt === DNAbin_Gap
end

@inline function BioSymbols.iscompatible(x::DNAbin, y::DNAbin)
    return (UInt8(x) >> 4) & (UInt8(y) >> 4) != 0
end

@inline function BioSymbols.complement(x::DNAbin)
    bits = UInt8(x)
    return DNAbin((bits & 0x10) << 3 | (bits & 0x80) >> 3 |
        (bits & 0x20) << 1 | (bits & 0x40) >> 1 | (bits & 0x0F))
end

@inline function BioSymbols.isvalid(x::DNAbin)
    return (UInt8(x) & 0x03) == 0x00
end

# Tranfer of data to and from R
# -----------------------------

@inline function check_class(rs::Ptr{RawSxp})
    if !isObject(rs) || isS4(rs) || rcopy(rcall_p(:class, rs)) != "DNAbin"
        error("Object is not of S3 class DNAbin, aborting")
    end
end

function rcopy(::Type{Array{DNAbin, 2}}, rs::Ptr{RawSxp})
    # Protect the RawSxp from GC.
    protect(rs)
    try
        # Check class of RawSxp.
        check_class(rs)
        return convert(Array{DNAbin, 2}, unsafe_array(rs))
    finally
        unprotect(1)
    end
end

function rcopy(::Type{Array{DNA, 2}}, rs::Ptr{RawSxp})
    # Protect the RawSxp from GC.
    protect(rs)
    try
        # Check class of RawSxp.
        check_class(rs)
        ra = unsafe_array(rs)
        ja = similar(ra, DNA)
        ja .= _raw_to_DNA.(ra)
        return ja
    finally
        unprotect(1)
    end
end

function rcopy(::Type{Vector{DNASequence}}, rs::Ptr{RawSxp})
    # Protect the RawSxp from GC.
    protect(rs)
    try
        nucs = rcopy(Array{DNA, 2}, rs)
        nSeq, seqLen = size(nucs)

        # Initialize the output vector of biological sequences.
        seqs = Vector{DNASequence}(nSeq)
        @inbounds for i in 1:nSeq
            seqs[i] = DNASequence(nucs[i,:])
        end
        return seqs
    finally
        unprotect(1)
    end
end

@inline function rcopytype(::Type{RClass{:DNAbin}}, s::Ptr{RawSxp})
    return Vector{DNASequence}
end

function sexp(::Type{RawSxp}, seq::Array{DNAbin,2})
    ra = protect(allocArray(RawSxp, size(seq)...))
    @inbounds for i ∈ 1:endof(seq)
        ra[i] = UInt8(seq[i])
    end
    setclass!(ra, sexp("DNAbin"))
    unprotect(1)
    return ra
end

function sexp(::Type{RawSxp}, seq::Vector{BioSequence{A}}) where {A<:DNAAlphabet}
    ncols = length(seq[1])
    @inbounds for i ∈ 2:endof(seq)
        if length(seq[i]) != ncols
            error("Sequences must all be the same length")
        end
    end
    ra = protect(allocArray(RawSxp, length(seq), ncols))
    z = 1
    @inbounds for j ∈ 1:ncols, i ∈ 1:endof(seq)
        ra[z] = _DNA_to_raw(seq[i][j])
        z += 1
    end
    setclass!(ra, sexp("DNAbin"))
    unprotect(1)
    return ra
end

sexp(seq::Vector{BioSequence{A}}) where {A<:DNAAlphabet} = sexp(RawSxp, seq)