text:
__precompile__()

module BuildUtils

using Compat
using Compat.Sys: iswindows, isapple

export rpath_origin, whole_archive, isexecutable, object_dependencies

"""
    rpath_origin([kernel]; escape=false)

Construct the appropriate linker flags used for setting the runtime search to include
`\$ORIGIN`, the location of the library itself that's looking for a required dependency.
If `escape` is `true`, the leading `\$` in `\$ORIGIN` is doubled. This sometimes useful
since `\$` has a special meaning in the shell.
"""
function rpath_origin(k::Symbol=Sys.KERNEL; escape::Bool=false)
    if isapple(k)
        "-Wl,-rpath,'@loader_path/'"
    elseif iswindows(k)
        ""
    else
        flags = "-Wl,-rpath,'\$ORIGIN' -Wl,-z,origin"
        escape ? replace(flags, '$', "\$\$") : flags
    end
end

"""
    whole_archive([kernel], lib)

Wrap the static library `lib`, specified as a `String`, in "whole archive" and
"no whole archive" flags.
"""
function whole_archive(k::Symbol, lib::String)
    if isapple(k)
        whole = "-Xlinker -all_load"
        nowhole = ""
    else
        whole = "-Wl,--whole-archive"
        nowhole = "-Wl,--no-whole-archive"
    end
    chomp(string(whole, " ", s, " ", nowhole))
end

"""
    isexecutable(program)

Determine whether the given program name or path is executable using the current user's
permissions. This is roughly equivalent to querying `which program` at the command line
and checking that a result is found, but no shelling out occurs.
"""
function isexecutable(k::Symbol, prog::String)
    access = iswindows(k) ? :_access : :access
    X_OK = 1 << 0 # Taken from unistd.h
    # If prog has a slash, we know the user wants to determine whether the given
    # file exists and is executable
    if '/' in prog || (iswindows(k) && '\\' in prog) # Windows can use / too
        isfile(prog) || return false
        return ccall(access, Cint, (Ptr{UInt8}, Cint), prog, X_OK) == 0
    end
    path = get(ENV, "PATH", "")
    # Something is definitely wrong if the user's path is empty...
    @assert !isempty(path)
    sep = iswindows(k) ? ';' : ':'
    for dir in split(path, sep), file in readdir(dir)
        if file == prog || (iswindows(k) && file == prog * ".exe")
            p = joinpath(dir, file)
            @assert isfile(p)
            return ccall(access, Cint, (Ptr{UInt8}, Cint), p, X_OK) == 0
        end
    end
    false
end

"""
    object_dependencies([kernel,] lib)

Determine the shared libraries on which the given library `lib` depends, returning
a `Dict` of library, path pairs.

!!! note
    This function shells out to `otool -L` on macOS, which requires having the Xcode
    command line tools installed. On other systems, it shells out to `ldd`, which
    can be installed with system package managers if it isn't already installed on
    the system.
"""
function object_dependencies(k::Symbol, lib::String)
    isfile(lib) || throw(ArgumentError("library '$lib' does not exist"))
    if isapple(k)
        isexecutable(k, "otool") || error("otool is not available")
        cmd = `otool -L`
    else
        isexecutable(k, "ldd") || error("ldd is not available")
        cmd = `ldd`
    end
    d = Dict{String,String}()
    # The first line is always the input library
    for line in Iterators.drop(eachline(`$cmd $lib`), 1)
        line = lstrip(line)
        # TODO: Move this check out of the loop?
        if isapple(k)
            found = line[1:prevind(line, findfirst(line, ' '))]
            dep = basename(found)
        else
            tokens = split(line)
            length(tokens) == 4 || error("unexpected line format: '$line'")
            dep = tokens[1]
            # Handle 'not found'
            found = tokens[3] == "not" ? join(tokens[3:4], ' ') : tokens[3]
        end
        d[dep] = found
    end
    d
end

for f in [:whole_archive, :isexecutable, :object_dependencies]
    @eval ($f)(s::String) = ($f)(Sys.KERNEL, s)
end

end # module


text:
function branch!(
    x,
    dx,
    output_ode::Vector{T},
    V_r_from,
    V_i_from,
    V_r_to,
    V_i_to,
    current_r_from,
    current_i_from,
    current_r_to,
    current_i_to,
    ix_range::UnitRange{Int},
    ode_range::UnitRange{Int},
    branch::PSY.DynamicBranch,
    inputs::SimulationInputs,
) where {T <: Real}
    #Obtain local device states
    device_states = @view x[ix_range]

    #Obtain references
    sys_f = get_base_frequency(inputs)

    mdl_line_ode!(
        device_states,
        view(output_ode, ode_range),
        V_r_from,
        V_i_from,
        V_r_to,
        V_i_to,
        current_r_from,
        current_i_from,
        current_r_to,
        current_i_to,
        sys_f,
        branch,
    )

    return
end

ext: jl


text:
module InterpolateTest

using ColorPrinting
using ReTest

using ColorPrinting: interpolate, @interpolate

@testset "basic" begin
    x = 1
    @test eval(interpolate("foo")) == "foo"
    @test eval(interpolate("foo$x")) == "foo1"
    @test eval(interpolate("foo$(x + 1)")) == "foo2"
    @test eval(interpolate("foo$(x + 1)bar")) == "foo2bar"
    @test eval(interpolate("foo$(x)bar$x")) == "foo1bar1"
    @test eval(interpolate("foo$(x + (x - 3))")) == "foo-1"
end

z(s) = s*s
macro bar_str(s)
    quote
        x = @interpolate $s
        z(x)
    end
end

@testset "macro interpolate" begin
    x = 3
    let y = 1
        @test bar"a$y" == "a1a1"
        @test bar"a$(x + y)b" == "a4ba4b"
    end
end

end # module


text:
# Very junky script testing performance against direct binary IO
#
# In this test, it's about 5x slower to write to ply than to directly dump the
# bytes, this seems to be due to the many calls to write()

using BufferedStreams

import PlyIO: save_ply, write_header, Ply, Element, ArrayProperty, ListProperty

ply = Ply()

nverts = 100000

vertex = Element("vertex",
                 ArrayProperty("x", randn(nverts)),
                 ArrayProperty("y", randn(nverts)),
                 ArrayProperty("z", randn(nverts)),
                 ArrayProperty("r", rand(nverts)),
                 ArrayProperty("g", rand(nverts)),
                 ArrayProperty("b", rand(nverts)))
push!(ply, vertex)

# Some triangular faces
vertex_index = ListProperty("vertex_index", Int32, Int32)
for i=1:nverts
   push!(vertex_index, rand(0:nverts-1,3))
end
push!(ply, Element("face", vertex_index))

#=
# Some edges
vertex_index = ListProperty("vertex_index", Int32, Int32)
for i=1:nverts
   push!(vertex_index, rand(0:nverts-1,2))
end
push!(ply, Element("edge", vertex_index))
=#

function save_ply_ref(ply, stream::IO)
    write_header(ply, stream, false)
    x = ply["vertex"]["x"].data::Vector{Float64}
    y = ply["vertex"]["y"].data::Vector{Float64}
    z = ply["vertex"]["z"].data::Vector{Float64}
    r = ply["vertex"]["r"].data::Vector{Float64}
    g = ply["vertex"]["g"].data::Vector{Float64}
    b = ply["vertex"]["b"].data::Vector{Float64}

    vertex_index = ply["face"]["vertex_index"]
    viinds = vertex_index.start_inds::Vector{Int32}
    vidata = vertex_index.data::Vector{Int32}

    #=
    for i=1:length(x)
        write(stream, x[i])
        write(stream, y[i])
        write(stream, z[i])
        write(stream, r[i])
        write(stream, g[i])
        write(stream, b[i])
    end

    for i=1:length(viinds)-1
        len = viinds[i+1] - viinds[i]
        write(stream, len)
        write(stream, vidata[viinds[i]:viinds[i+1]-1])
    end
    =#

    # Buffered reordering
    #=
    write(stream, [x y z r g b]')
    write(stream, viinds)
    write(stream, vidata)
    =#

    # Benchmark against direct binary IO (invalid ply!), which should be about
    # as fast as you can hope for.
    write(stream, x)
    write(stream, y)
    write(stream, z)
    write(stream, r)
    write(stream, g)
    write(stream, b)

    write(stream, viinds)
    write(stream, vidata)

end

function save_ply_ref(ply, filename::AbstractString)
    open(filename, "w") do fid
        save_ply_ref(ply, fid)
    end
end

function save_ply_buffered(ply, filename; kwargs...)
    open(filename, "w") do fid
        save_ply(ply, BufferedOutputStream(fid, 2^16); kwargs...)
    end
end

save_ply(ply, "test.ply", ascii=false)
save_ply_buffered(ply, "test.ply", ascii=false)

save_ply_ref(ply, "test2.ply")

Profile.clear_malloc_data()
@time save_ply(ply, "test.ply", ascii=false)
@time save_ply_buffered(ply, "test.ply", ascii=false)

@time save_ply_ref(ply, "test2.ply")

#run(`displaz -script test.ply`)


text:
########################################
## File Name: state_types.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/12
## Description: State Type Definitions for SACBP
########################################

import Base.vec

abstract type PhysState end;

@auto_hash_equals struct PhysPos{T<:Real} <: PhysState
    t::Float64
    pos::Vector{T} # [x1,x2,x3,...]
    dim::Int64
    function PhysPos{T}(t,pos,dim) where {T<:Real}
        if dim != length(pos)
            error(ArgumentError("Invalid position vector length."));
        else
            return new(t,pos,dim);
        end
    end
end
PhysPos(t::Real,pos::Vector{T},dim=length(pos)) where {T<:Real} = PhysPos{T}(t,pos,dim);
PhysPos(t::Real,pos::T,dim=1) where {T<:Real} = PhysPos{T}(t,[pos],dim);
vec(p::PhysPos) = p.pos;

@auto_hash_equals struct PhysManipulate2D{T<:Real} <: PhysState
    t::Float64
    pos::Vector{T} # [x,y]
    θ::T # 0 <= θ < 2pi
    vel::Vector{T} # [vx,vy]
    ω::T
    m::T
    J::T
    r::Vector{T}
    μ::T
    dim::Int64
    function PhysManipulate2D{T}(t,pos,θ,vel,ω,m,J,r,μ,dim) where {T<:Real}
        if dim != 2 || length(pos) != 2 || length(vel) != 2 || length(r) != 2
            error(ArgumentError("Invalid dimension or parameter vector length."));
        else
            return new(t,pos,mod2pi(θ),vel,ω,m,J,r,μ,dim);
        end
    end
end
function PhysManipulate2D(t::Real,pos::Vector{T},θ::T,vel::Vector{T},ω::T,m::T,J::T,r::Vector{T},μ::T,dim::Int64=2) where {T<:Real}
    return PhysManipulate2D{T}(t,pos,θ,vel,ω,m,J,r,μ,dim)
end
function PhysManipulate2D(t::Real,params::Vector{T}) where {T<:Real}
    if length(params) != 11
        error(ArgumentError("Invalid parameter vector length."))
    else
        return PhysManipulate2D(t,params[1:2],params[3],params[4:5],params[6],params[7],params[8],params[9:10],params[11])
    end
end
vec(p::PhysManipulate2D) = vcat(p.pos,p.θ,p.vel,p.ω,p.m,p.J,p.r,p.μ);

@auto_hash_equals struct AugState{P<:PhysState,B<:BelState}
    t::Float64
    p::P
    b::B
    function AugState{P,B}(t,p,b) where {P<:PhysState} where {B<:BelState}
        if isa(b,Belief)
            if t != p.t || t != b.t
                error(ArgumentError("PhysState and BelState have inconsistent time parameters."))
            else
                return new(t,p,b)
            end
        else
            if t != p.t || t != b[1].t
                error(ArgumentError("PhysState and BelState have inconsistent time parameters."))
            else
                return new(t,p,b)
            end
        end
    end
end
AugState(t::Real,p::P,b::B) where {P<:PhysState} where {B<:BelState} = AugState{P,B}(t,p,b);
AugState(p::P,b::B) where {P<:PhysState} where {B<:BelState} = AugState{P,B}(p.t,p,b);


text:
nonsquare(n::Real) = n + floor(typeof(n), 0.5 + sqrt(n))
@show nonsquare.(1:1_000_000) ∩ collect(1:1000) .^ 2

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import UDUnits
using Test

@testset "UDUnits" begin
    system = UDUnits.System()

    unit_m = UDUnits.Unit(system,"m")
    unit_km = UDUnits.Unit(system,"km")
    unit_cm = UDUnits.Unit(system,"cm")
    unit_s = UDUnits.Unit(system,"s")
    unit_h = UDUnits.Unit(system,"h")
    unit_W = UDUnits.Unit(system,"W")
    unit_kg = UDUnits.Unit(system,"kg")
    unit_kg = system["kg"]
    unit_μm = system["μm"]

    @test haskey(system,"kg")
    @test haskey(system,"μm")

    @test_throws ErrorException UDUnits.Unit(system,"badname")

    @test UDUnits.areconvertible(unit_m,unit_cm)
    @test_throws ErrorException UDUnits.Converter(unit_m,unit_kg)

    buf = IOBuffer()
    show(buf,unit_km)
    @test startswith(String(take!(buf)),"<Unit:")

    # converter
    conv = UDUnits.Converter(unit_cm,unit_m)
    @test UDUnits.convert(conv,100.) ≈ 1.
    @test UDUnits.expression(conv) == "0.01*x"

    buf = IOBuffer()
    show(buf,conv)
    @test startswith(String(take!(buf)),"<Converter:")

    unit_m3 = unit_m + 3
    conv = UDUnits.Converter(unit_m,unit_m3)
    @test UDUnits.convert(conv,100.) ≈ 97.

    unit_mm3 = unit_m - 3.
    conv = UDUnits.Converter(unit_m,unit_mm3)
    @test UDUnits.convert(conv,100.) ≈ 103.

    unit_m_per_s = unit_m/unit_s
    unit_km_per_h = unit_km/unit_h
    unit_Hz = 1/unit_s

    @test UDUnits.areconvertible(unit_m_per_s,unit_km_per_h)
    conv = UDUnits.Converter(unit_m_per_s,unit_km_per_h)

    @test UDUnits.convert(conv,100.) ≈ 360.
    @test conv(100.) ≈ 360.

    unit_W2 = unit_kg * unit_m^2/(unit_s^3)
    conv = UDUnits.Converter(unit_W2,unit_W)
    @test UDUnits.convert(conv,100.) ≈ 100.

    @test UDUnits.areconvertible(√(unit_m * unit_m),unit_km)

    unit_logW = log10(unit_W)
    @test UDUnits.areconvertible(unit_W,unit_logW)

    unit_logeW = log(unit_W)
    @test UDUnits.areconvertible(unit_W,unit_logeW)

    unit_s2 = (unit_s * unit_s) ^(1//2)

    @test UDUnits.areconvertible(unit_s2,unit_s)

    @test UDUnits.name(unit_m) == "meter"
    @test UDUnits.symbol(unit_m) == "m"

    # Format strings
    unit_J_per_s = UDUnits.Unit(system,"J/s")

    @test UDUnits.format(unit_J_per_s) == "W"
    @test UDUnits.string(unit_J_per_s) == "W"

    @test UDUnits.format(unit_J_per_s; names = true) == "watt"
    @test UDUnits.format(unit_J_per_s; definition = true) == "m²·kg·s⁻³"
    @test UDUnits.format(unit_J_per_s; definition = true, names = true) ==
        "meter²·kilogram·second⁻³"

    # test clean-up
    unit_m = 0
    system = 0
    conv = 0
    GC.gc()
end


text:
module Metis

using SparseArrays
using LinearAlgebra
import LightGraphs
import TriangleMesh
using METIS_jll: libmetis

# Metis C API
include("metis_h.jl")
const options = fill(Cint(-1), METIS_NOPTIONS)
options[METIS_OPTION_NUMBERING] = 1

# Julia interface
"""
    Metis.Graph

1-based CSR representation of a graph as defined in
section 5.5 "Graph data structure" in the Metis manual.
"""
struct Graph
    nvtxs::idx_t
    xadj::Vector{idx_t}
    adjncy::Vector{idx_t}
    vwgt::Vector{idx_t}
    Graph(nvtxs, xadj, adjncy) = new(nvtxs, xadj, adjncy)
    Graph(nvtxs, xadj, adjncy, vwgt) = new(nvtxs, xadj, adjncy, vwgt)
end

"""
    Metis.Mesh

Representation of a mesh as defined in section 5.6
"Mesh data structure" in the Metis manual.
"""
struct Mesh
    ne::idx_t
    nn::idx_t
    eptr::Vector{idx_t}
    eind::Vector{idx_t}
    Mesh(ne, nn, eptr, eind) = new(ne, nn, eptr, eind)
end

"""
    Mesh(mesh::TriangleMesh.TriMesh)

Constructor of Metis.Mesh using a TriangleMesh.TriMesh.
"""
function Mesh(mesh::TriangleMesh.TriMesh)
    ne = idx_t(mesh.n_cell)
    nn = idx_t(mesh.n_point)
    eptr = Vector{idx_t}(undef, ne+1)
    eind = Vector{idx_t}(undef, length(mesh.cell))
    nn_per_elem, _ = size(mesh.cell)
    eptr[:] = [idx_t(nn_per_elem*i) for i in 0:ne]
    eind[:] = reshape(mesh.cell, length(mesh.cell))
    return Mesh(ne, nn, eptr, eind)
end

"""
    Metis.graph(G::SparseMatrixCSC; check_hermitian=true)

Construct the 1-based CSR representation of the sparse matrix `G`.
If `check_hermitian` is `false` the matrix is not checked for being hermitian
before constructing the graph.
"""
function graph(G::SparseMatrixCSC; check_hermitian=true)
    if check_hermitian
        ishermitian(G) || throw(ArgumentError("matrix must be Hermitian"))
    end
    N = size(G, 1)
    xadj = Vector{idx_t}(undef, N+1)
    xadj[1] = 1
    adjncy = Vector{idx_t}(undef, nnz(G))
    adjncy_i = 0
    @inbounds for j in 1:N
        n_rows = 0
        for k in G.colptr[j] : (G.colptr[j+1] - 1)
            i = G.rowval[k]
            if i != j # don't include diagonal elements
                n_rows += 1
                adjncy_i += 1
                adjncy[adjncy_i] = i
            end
        end
        xadj[j+1] = xadj[j] + n_rows
    end
    resize!(adjncy, adjncy_i)
    return Graph(idx_t(N), xadj, adjncy)
end

"""
    graph(G::LightGraphs.AbstractSimpleGraph)

Construct the 1-based CSR representation of the `LightGraphs` graph `G`.
"""
function graph(G::LightGraphs.AbstractSimpleGraph)
    N = LightGraphs.nv(G)
    xadj = Vector{idx_t}(undef, N+1)
    xadj[1] = 1
    adjncy = Vector{idx_t}(undef, 2*LightGraphs.ne(G))
    adjncy_i = 0
    for j in 1:N
        ne = 0
        for i in LightGraphs.outneighbors(G, j)
            ne += 1
            adjncy_i += 1
            adjncy[adjncy_i] = i
        end
        xadj[j+1] = xadj[j] + ne
    end
    resize!(adjncy, adjncy_i)
    return Graph(idx_t(N), xadj, adjncy)
end

"""
    perm, iperm = Metis.permutation(G)

Compute the fill reducing permutation `perm`
and its inverse `iperm` of `G`.
"""
permutation(G) = permutation(graph(G))

function permutation(G::Graph)
    perm = Vector{idx_t}(undef, G.nvtxs)
    iperm = Vector{idx_t}(undef, G.nvtxs)
    vwgt = isdefined(G, :vwgt) ? G.vwgt : C_NULL
    METIS_NodeND(G.nvtxs, G.xadj, G.adjncy, vwgt, options, perm, iperm)
    return perm, iperm
end

"""
    Metis.partition(G, n; alg = :KWAY)

Partition the graph `G` in `n` parts.
The partition algorithm is defined by the `alg` keyword:
 - :KWAY: multilevel k-way partitioning
 - :RECURSIVE: multilevel recursive bisection
"""
partition(G, nparts; alg = :KWAY) = partition(graph(G), nparts, alg = alg)

function partition(G::Graph, nparts::Integer; alg = :KWAY)
    part = Vector{idx_t}(undef, G.nvtxs)
    vwgt = isdefined(G, :vwgt) ? G.vwgt : C_NULL
    edgecut = fill(idx_t(0), 1)
    if alg === :RECURSIVE
        METIS_PartGraphRecursive(G.nvtxs, idx_t(1), G.xadj, G.adjncy, vwgt, C_NULL, C_NULL,
                                 idx_t(nparts), C_NULL, C_NULL, options, edgecut, part)
    elseif alg === :KWAY
        METIS_PartGraphKway(G.nvtxs, idx_t(1), G.xadj, G.adjncy, vwgt, C_NULL, C_NULL,
                            idx_t(nparts), C_NULL, C_NULL, options, edgecut, part)
    else
        throw(ArgumentError("unknown algorithm $(repr(alg))"))
    end
    return part
end

"""
    Metis.separator(G)

Compute a vertex separator of the graph `G`.
"""
separator(G) = separator(graph(G))

function separator(G::Graph)
    part = Vector{idx_t}(undef, G.nvtxs)
    sepsize = fill(idx_t(0), 1)
    vwgt = isdefined(G, :vwgt) ? G.vwgt : C_NULL
    # METIS_ComputeVertexSeparator segfaults with 1-based indexing
    xadj = G.xadj .- idx_t(1)
    adjncy = G.adjncy .- idx_t(1)
    METIS_ComputeVertexSeparator(G.nvtxs, xadj, adjncy, vwgt, options, sepsize, part)
    part .+= 1
    return part
end

"""
    Metis.mesh_partition(M, n; alg = :DUAL)

Partition the mesh `M` in `n` parts.
The partition algorithm is defined by the `alg` keyword:
 - :DUAL: Partition of the mesh's dual graph
 - :NODAL: Partition of the mesh's nodal graph
"""
function mesh_partition(mesh::TriangleMesh.TriMesh, nparts::Integer; alg = :DUAL)
    M = Mesh(mesh)
    epart = Vector{idx_t}(undef, M.ne)
    npart = Vector{idx_t}(undef, M.nn)
    objval = fill(idx_t(0), 1)
    if alg === :NODAL
        METIS_PartMeshNodal(M.ne, M.nn, M.eptr, M.eind, C_NULL, C_NULL, idx_t(nparts),
                            C_NULL, options, objval, epart, npart)
    elseif alg === :DUAL
        METIS_PartMeshDual(M.ne, M.nn, M.eptr, M.eind, C_NULL, C_NULL, idx_t(2), idx_t(nparts),
                           C_NULL, options, objval, epart, npart)
    else
        throw(ArgumentError("unknown algorithm $(repr(alg))"))
    end
    return epart, npart
end

end # module


text:
using ApproxFun

x=Fun()
u0=0.0x
N=u->[u(-1.),u(1.),0.001u''+6*(1-x^2)*u'+u^2-1.]
u=newton(N,u0)

x=Fun()
u0=0.0x
N=u->[u(-1.),u(1.),0.1u''+sin(u)-1.]
u=newton(N,u0)

max_stars_repo_path: ODEs/Newton iteration for Nonlinear BVPs.jl


text:
# temporary wrapper to pass byte vector through
struct VectorString{T <: AbstractVector{UInt8}} <: AbstractString
    bytes::T
end

Base.codeunits(x::VectorString) = x.bytes

# high-level user API functions
read(io::IO; kw...) = read(Base.read(io, String); kw...)
read(bytes::AbstractVector{UInt8}; kw...) = read(VectorString(bytes); kw...)

function read(str::AbstractString; jsonlines::Bool=false, kw...)
    buf = codeunits(str)
    len = length(buf)
    if len == 0
        error = UnexpectedEOF
        pos = 0
        @goto invalid
    end
    pos = 1
    b = getbyte(buf, pos)
    @wh
    tape = len < 1000 ? Vector{UInt64}(undef, len + 4) :
        Vector{UInt64}(undef, div(len, 10))
    if jsonlines
        pos, tapeidx = jsonlines!(buf, pos, len, b, tape, Int64(1); kw...)
    else
        pos, tapeidx = read!(buf, pos, len, b, tape, Int64(1), Any; kw...)
    end
    @inbounds t = tape[1]
    if isobject(t)
        obj = Object(buf, tape, Dict{Symbol, Int}())
        populateinds!(obj)
        return obj
    elseif isarray(t)
        arr = Array{geteltype(tape[2])}(buf, tape, Int[])
        populateinds!(arr)
        return arr
    else
        return getvalue(Any, buf, tape, 1, t)
    end
@label invalid
    invalid(error, buf, pos, Any)
end

macro check()
    esc(quote
        if (tapeidx + 1) > length(tape)
            newsize = ceil(Int64, ((1 - pos / len) + 1) * tapeidx) + 20
            # println("resizing tape from $(pointer(tape)) $tapeidx to $newsize")
            resize!(tape, newsize)
        end
    end)
end

function read!(buf, pos, len, b, tape, tapeidx, ::Type{Any}, checkint=true; allow_inf::Bool=false)
    if b == UInt8('{')
        return read!(buf, pos, len, b, tape, tapeidx, Object; allow_inf=allow_inf)
    elseif b == UInt8('[')
        return read!(buf, pos, len, b, tape, tapeidx, Array; allow_inf=allow_inf)
    elseif b == UInt8('"')
        return read!(buf, pos, len, b, tape, tapeidx, String)
    elseif b == UInt8('n')
        return read!(buf, pos, len, b, tape, tapeidx, Nothing)
    elseif b == UInt8('t')
        return read!(buf, pos, len, b, tape, tapeidx, True)
    elseif b == UInt8('f')
        return read!(buf, pos, len, b, tape, tapeidx, False)
    elseif (UInt8('0') <= b <= UInt8('9')) || b == UInt8('-') || b == UInt8('+') || (allow_inf && (b == UInt8('N') || b == UInt8('I')))
        float, code, floatpos = Parsers.typeparser(Float64, buf, pos, len, b, Int16(0), Parsers.OPTIONS)
        if code > 0
            !isfinite(float) && !allow_inf && @goto invalid
            @check
            if checkint
                int = unsafe_trunc(Int64, float)
                if int == float
                    @inbounds tape[tapeidx] = INT
                    @inbounds tape[tapeidx + 1] = Core.bitcast(UInt64, int)
                else
                    @inbounds tape[tapeidx] = FLOAT
                    @inbounds tape[tapeidx + 1] = Core.bitcast(UInt64, float)
                end
            else
                @inbounds tape[tapeidx] = FLOAT
                @inbounds tape[tapeidx + 1] = Core.bitcast(UInt64, float)
            end
            return floatpos, tapeidx + 2
        end
        invalid(InvalidNumber, buf, pos, Float64)
    end
@label invalid
    invalid(InvalidChar, buf, pos, Any)
end

function read!(buf, pos, len, b, tape, tapeidx, ::Type{String})
    pos += 1
    @eof
    strpos = pos
    strlen = Int64(0)
    escaped = false
    b = getbyte(buf, pos)
    while b != UInt8('"')
        if b == UInt8('\\')
            escaped = true
            # skip next character
            pos += 2
            strlen += 2
        elseif b < UInt8(' ')
            unescaped_control(b)
        else
            pos += 1
            strlen += 1
        end
        @eof
        b = getbyte(buf, pos)
    end
    @check
    @inbounds tape[tapeidx] = string(strlen)
    @inbounds tape[tapeidx+1] = ifelse(escaped, ESCAPE_BIT | strpos, strpos)
    return pos + 1, tapeidx + 2

@label invalid
    invalid(error, buf, pos, String)
end

struct True end
function read!(buf, pos, len, b, tape, tapeidx, ::Type{True})
    if pos + 3 <= len &&
        b            == UInt8('t') &&
        buf[pos + 1] == UInt8('r') &&
        buf[pos + 2] == UInt8('u') &&
        buf[pos + 3] == UInt8('e')
        @check
        @inbounds tape[tapeidx] = BOOL | UInt64(1)
        return pos + 4, tapeidx + 2
    else
        invalid(InvalidChar, buf, pos, True)
    end
end

struct False end
function read!(buf, pos, len, b, tape, tapeidx, ::Type{False})
    if pos + 4 <= len &&
        b            == UInt8('f') &&
        buf[pos + 1] == UInt8('a') &&
        buf[pos + 2] == UInt8('l') &&
        buf[pos + 3] == UInt8('s') &&
        buf[pos + 4] == UInt8('e')
        @check
        @inbounds tape[tapeidx] = BOOL
        return pos + 5, tapeidx + 2
    else
        invalid(InvalidChar, buf, pos, False)
    end
end

function read!(buf, pos, len, b, tape, tapeidx, ::Type{Nothing})
    if pos + 3 <= len &&
        b            == UInt8('n') &&
        buf[pos + 1] == UInt8('u') &&
        buf[pos + 2] == UInt8('l') &&
        buf[pos + 3] == UInt8('l')
        @check
        @inbounds tape[tapeidx] = NULL
        return pos + 4, tapeidx + 2
    else
        invalid(InvalidChar, buf, pos, Nothing)
    end
end

function read!(buf, pos, len, b, tape, tapeidx, ::Type{Object}; kw...)
    objidx = tapeidx
    eT = EMPTY
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    if b == UInt8('}')
        @check
        @inbounds tape[tapeidx] = object(Int64(2))
        @inbounds tape[tapeidx+1] = eltypelen(eT, Int64(0))
        tapeidx += 2
        pos += 1
        @goto done
    elseif b != UInt8('"')
        error = ExpectedOpeningQuoteChar
        @goto invalid
    end
    pos += 1
    @eof
    tapeidx += 2
    nelem = Int64(0)
    while true
        keypos = pos
        keylen = Int64(0)
        escaped = false
        # read first key character
        b = getbyte(buf, pos)
        # positioned at first character of object key
        while b != UInt8('"')
            if b == UInt8('\\')
                escaped = true
                # skip next character
                pos += 2
                keylen += 2
            else
                pos += 1
                keylen += 1
            end
            @eof
            b = getbyte(buf, pos)
        end
        @check
        @inbounds tape[tapeidx] = string(keylen)
        @inbounds tape[tapeidx+1] = ifelse(escaped, ESCAPE_BIT | keypos, keypos)
        tapeidx += 2
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b != UInt8(':')
            error = ExpectedSemiColon
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        # now positioned at start of value
        prevtapeidx = tapeidx
        pos, tapeidx = read!(buf, pos, len, b, tape, tapeidx, Any; kw...)
        @eof
        b = getbyte(buf, pos)
        @wh
        @inbounds eT = promoteeltype(eT, gettypemask(tape[prevtapeidx]))
        nelem += 1
        if b == UInt8('}')
            @check
            @inbounds tape[objidx] = object(tapeidx - objidx)
            @inbounds tape[objidx+1] = eltypelen(eT, nelem)
            pos += 1
            @goto done
        elseif b != UInt8(',')
            error = ExpectedComma
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
        if b != UInt8('"')
            error = ExpectedOpeningQuoteChar
            @goto invalid
        end
        pos += 1
        @eof
    end

@label done
    return pos, tapeidx
@label invalid
    invalid(error, buf, pos, Object)
end

function read!(buf, pos, len, b, tape, tapeidx, ::Type{Array}; kw...)
    arridx = tapeidx
    eT = EMPTY
    pos += 1
    @eof
    b = getbyte(buf, pos)
    @wh
    if b == UInt8(']')
        @check
        @inbounds tape[tapeidx] = array(Int64(2))
        @inbounds tape[tapeidx+1] = eltypelen(eT, Int64(0))
        tapeidx += 2
        pos += 1
        @goto done
    end
    tapeidx += 2
    nelem = Int64(0)
    while true
        # positioned at start of value
        prevtapeidx = tapeidx
        pos, tapeidx = read!(buf, pos, len, b, tape, tapeidx, Any, eT != FLOAT && eT != (FLOAT | NULL); kw...)
        @eof
        b = getbyte(buf, pos)
        @wh
        @inbounds eT = promoteeltype(eT, gettypemask(tape[prevtapeidx]))
        nelem += 1
        if b == UInt8(']')
            @check
            @inbounds tape[arridx] = array(tapeidx - arridx)
            @inbounds tape[arridx+1] = eltypelen(eT, nelem)
            pos += 1
            @goto done
        elseif b != UInt8(',')
            error = ExpectedComma
            @goto invalid
        end
        pos += 1
        @eof
        b = getbyte(buf, pos)
        @wh
    end

@label done
    return pos, tapeidx
@label invalid
    invalid(error, buf, pos, Array)
end

function jsonlines!(buf, pos, len, b, tape, tapeidx; kw...)
    arridx = tapeidx
    eT = EMPTY
    if pos > len
        @check
        @inbounds tape[tapeidx] = array(Int64(2))
        @inbounds tape[tapeidx+1] = eltypelen(eT, Int64(0))
        tapeidx += 2
        @goto done
    end
    tapeidx += 2
    nelem = Int64(0)
    while true
        @wh
        # positioned at start of value
        prevtapeidx = tapeidx
        pos, tapeidx = read!(buf, pos, len, b, tape, tapeidx, Any, eT != FLOAT && eT != (FLOAT | NULL); kw...)
        @inbounds eT = promoteeltype(eT, gettypemask(tape[prevtapeidx]))
        nelem += 1
        if pos > len
            @check
            @inbounds tape[arridx] = array(tapeidx - arridx)
            @inbounds tape[arridx+1] = eltypelen(eT, nelem)
            @goto done
        end
        b = getbyte(buf, pos)
        if b != UInt8('\n') && b != UInt8('\r')
            error = ExpectedComma
            @goto invalid
        end
        pos += 1
        if pos > len
            @check
            @inbounds tape[arridx] = array(tapeidx - arridx)
            @inbounds tape[arridx+1] = eltypelen(eT, nelem)
            @goto done
        end
        if b == UInt8('\r')
            b = getbyte(buf, pos)
            if b == UInt8('\n')
                pos += 1
                if pos > len
                    @check
                    @inbounds tape[arridx] = array(tapeidx - arridx)
                    @inbounds tape[arridx+1] = eltypelen(eT, nelem)
                    @goto done
                end
            end
        end
        b = getbyte(buf, pos)
    end

@label done
    return pos, tapeidx
@label invalid
    invalid(error, buf, pos, Array)
end


text:
import ..Resampling: resample
import ..UncertainValues: UncertainValue
import ..UncertainDatasets: UncertainIndexDataset

include("add_uncertainindexdataset.jl")
include("subtract_uncertainindexdataset.jl")
include("multiply_uncertainindexdataset.jl")
include("divide_uncertainindexdataset.jl")
include("exponentiate_uncertainindexdataset.jl")


text:
#=
<inpaint> application of the resynthesizer algorithm.

"inpaint" is tradional term.  In GIMP app, AKA "heal selection"

Takes one image and a mask.

Returns new image with the masked region synthesized from its *entire* surroundings.

The one image is first divided into a target and corpus image.
The target image is the input.
The algorithm mutates the target image in the masked region, called the "synth" (a noun) region.

The corpus is inverse of mask applied to input.
The corpus is immutable.
=#

include("../resynthesizer.jl")

function inpaint(
    image::AbstractArray{T,N},
    mask::AbstractArray{Bool,N},
    ) where {N,T}

    # Split one image into two

    targetImage = MaskedImage(image, mask)         # image is mutable

    # The corpus is an immutable copy of the in image, with an inverted mask
    invertedMask = .!mask
    corpusImage = MaskedImage(copy(image), invertedMask)
    @assert isconcretetype(typeof(targetImage.mask))

    resynthesize(targetImage, corpusImage)
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
module HoeffdingD

export hoeffdingd

using StatsBase
using LazyArrays ## for @~ to speed up Qi calculations

"""
    reject_independence(D, n, α)

    Compute the D-test of independence.
    Compute the D-test of independence.
        Given:
        a D (the calculated D measure),
        a n (the number of observations),
        and a α (0 < α < 1) the desired level of significance.
    Return true to indicate rejection of H₀ hypothesis of independence.
    The implementation is based upon chapter 9 in the paper from Wassily Hoeffding (1948).

"""

function reject_independence(D, n, α)
    ρ = √((2 * (n^2 +5n -32)) / (9n * (n-1)*(n-3)*(n-4)*α))
    return D > ρ
end

"""
    hoeffdingd(a::AbstractVector{Union{Missing,T}}, b::AbstractVector{Union{Missing,T}}, pvalue=true, α=0.05) where {T <: Real}

Compute the Hoeffding measure of independence of the two continuous random variables. As defined by Hoeffding [A Non-Parametric Test of Independence](https://projecteuclid.org/journals/annals-of-mathematical-statistics/volume-19/issue-4/A-Non-Parametric-Test-of-Independence/10.1214/aoms/1177730150.full) in chapter 5.

This is a wrapper to remove missing data by pair, and then call the core hoeffdingd to return the adequate value.

"""
function hoeffdingd(a::AbstractVector{Union{Missing,T}}, b::AbstractVector{Union{Missing,T}}) where {T <: Number}

    missa = findall(!ismissing, a)
    missb = findall(!ismissing, b)
    gindexes = findall(in(missa), missb)

    nma = convert(Vector{T}, view(a, view(missa, gindexes)))
    nmb = convert(Vector{T}, view(b, view(missb, gindexes)))

    return hoeffdingd(nma, nmb)
end

"""
    function hoeffdingd(a::AbstractVector{<:T}, b::AbstractVector{<:T}, α) where {T <: Number}

Compute the Hoeffding measure of independence of the two continuous random variables. As defined by Hoeffding [A Non-Parametric Test of Independence](https://projecteuclid.org/journals/annals-of-mathematical-statistics/volume-19/issue-4/A-Non-Parametric-Test-of-Independence/10.1214/aoms/1177730150.full) in chapter 5.

Compute the Hoeffding measure of dependence (D) together with the result of the test of independence when α is presented.
Check that α is between 0 and 1.
Return the results from the core hoeffdingd and the result of the independence test in a tuple.

"""
function hoeffdingd(a::AbstractVector{<:T}, b::AbstractVector{<:T}, α) where {T <: Number}
    if (α <= 0)  || (α >= 1)
        throw(DomainError(α, "α must be: 0 < α < 1"))
    end
    D = hoeffdingd(a,b)
    N = size(a,1)
    return D , reject_independence(D, N, α)
end

"""
    function hoeffdingd(a::AbstractVector{<:T}, b::AbstractVector{<:T}) where {T <: Number}

Compute the Hoeffding measure of independence of the two continuous random variables. As defined by Hoeffding [A Non-Parametric Test of Independence](https://projecteuclid.org/journals/annals-of-mathematical-statistics/volume-19/issue-4/A-Non-Parametric-Test-of-Independence/10.1214/aoms/1177730150.full) in chapter 5.

Returns missing when there is not enough observations (less than 5).

The implementation relies on:
    StatsBase.tiedrank(a) to avoid duplication,
    and LazyArrays (@~) to improve performance.

"""
function hoeffdingd(a::AbstractVector{<:T}, b::AbstractVector{<:T}) where {T <: Number}
    N = size(a, 1)
    if (N != size(b, 1))
        error("Both array should have the same number of (non missing) observations")
    end
    if N < 5
        # error("Not enough observations ($(N) < 5) for Hoeffding's D")
        return missing
    end

    R = StatsBase.tiedrank(a)
    S = StatsBase.tiedrank(b)
    Q = Vector{Float64}(undef, N)

    @inbounds for i in 1:N
        Q[i] = 1 + sum(@~( (R .< R[i]) .& (S .< S[i]) ))
        Q[i] += 0.25 * (sum(@~((R .== R[i]) .& (S .== S[i]))) - 1)
        Q[i] += 0.5 * sum(@~((R .== R[i]) .& (S .< S[i])))
        Q[i] += 0.5 * sum(@~((R .< R[i]) .& (S .== S[i])))
    end

    D1 = sum((Q .- 1) .* (Q .- 2))
    D2 = sum((R .- 1) .* (R .- 2) .* (S .- 1) .* (S .- 2))
    D3 = sum((R .- 2) .* (S .- 2) .* (Q .- 1))
    D = 30 * ((N - 2) * (N - 3) * D1 + D2 - 2 * (N - 2) * D3) / (N * (N - 1) * (N - 2) * (N - 3) * (N - 4))

    return D
end

"""
    hoeffdingd(m::AbstractMatrix{<:T}, α) where {T <: Number}

Compute the Hoeffding measure of independence and the test of independence for each pair of columns in the matrix m.
Returns a matrix with the results.

"""
function hoeffdingd(m::AbstractMatrix{<:T}, α) where {T <: Number}

    cols = size(m)[2]
    results = Array{Tuple{Float64, Bool}}(undef, cols, cols)
    @inbounds for i in 1:cols
        results[i, i] = hoeffdingd(view(m, :, i), view(m, :, i), α)
        @inbounds for j in i + 1:cols
            results[i, j] = hoeffdingd(view(m, :, i), view(m, :, j), α)
            results[j, i] = results[i,j]
        end
    end
    return results

end

"""
    hoeffdingd(m::AbstractMatrix{<:T}) where {T <: Number}

Compute the Hoeffding measure of independence for each pair of columns in the matrix m.
Returns a matrix with the results.

"""
function hoeffdingd(m::AbstractMatrix{<:T}) where {T <: Number}

    cols = size(m)[2]
    results = Array{Float64}(undef, cols, cols)
    @inbounds for i in 1:cols
        results[i, i] = hoeffdingd(view(m, :, i), view(m, :, i))
        @inbounds for j in i + 1:cols
            results[i, j] = hoeffdingd(view(m, :, i), view(m, :, j))
            results[j, i] = results[i,j]
        end
    end
    return results

end

end

max_stars_repo_name: ericqu/HoeffdingD

max_issues_repo_name: ericqu/HoeffdingD

max_forks_repo_name: ericqu/HoeffdingD


text:
include("ncu.jl")
include("grover.jl")
include("qft.jl")
include("rqc.jl")
include("evolution.jl")