text:
struct OrbitGrid{T}
    energy::AbstractVector{T}
    pitch::AbstractVector{T}
    r::AbstractVector{T}
    counts::Vector{Int}
    orbit_index::Array{Int,3}
    class::Array{Symbol,3}
    tau_p::Array{T,3}
    tau_t::Array{T,3}
end

function Base.show(io::IO, og::OrbitGrid)
    print(io, "OrbitGrid: $(length(og.energy))×$(length(og.pitch))×$(length(og.r)):$(length(og.counts))")
end

function orbit_grid(M::AxisymmetricEquilibrium, eo::AbstractVector, po::AbstractVector, ro::AbstractVector;
                    q = 1, amu = H2_amu, kwargs...)

    nenergy = length(eo)
    npitch = length(po)
    nr = length(ro)

    orbit_index = zeros(Int,nenergy,npitch,nr)
    class = fill(:incomplete,(nenergy,npitch,nr))
    tau_t = zeros(Float64,nenergy,npitch,nr)
    tau_p = zeros(Float64,nenergy,npitch,nr)

    norbs = nenergy*npitch*nr
    subs = CartesianIndices((nenergy,npitch,nr))

    p = Progress(norbs)
    channel = RemoteChannel(()->Channel{Bool}(norbs), 1)
    orbs = fetch(@sync begin
        @async while take!(channel)
            ProgressMeter.next!(p)
        end
        @async begin
            orbs = @distributed (vcat) for i=1:norbs
                ie,ip,ir = Tuple(subs[i])
                c = EPRCoordinate(M,eo[ie],po[ip],ro[ir],q=q,amu=amu)
                try
                    o = get_orbit(M, c; kwargs...)
                catch
                    o = Orbit(EPRCoordinate(;q=q,amu=amu),:incomplete)
                end

                if o.class in (:incomplete,:invalid,:lost)
                    o = Orbit(o.coordinate,:incomplete)
                end
                put!(channel, true)
                o
            end
            put!(channel, false)
            orbs
        end
    end)

    for i=1:norbs
        class[subs[i]] = orbs[i].class
        tau_p[subs[i]] = orbs[i].tau_p
        tau_t[subs[i]] = orbs[i].tau_t
    end

    grid_index = filter(i -> orbs[i].class != :incomplete, 1:norbs)
    orbs = filter(x -> x.class != :incomplete, orbs)
    norbs = length(orbs)
    orbit_index[grid_index] = 1:norbs

    return orbs, OrbitGrid(eo,po,ro,fill(1,norbs),orbit_index,class,tau_p,tau_t)

end

function segment_orbit_grid(M::AxisymmetricEquilibrium, orbit_grid::OrbitGrid, orbits::Vector;
                        norbits=1000, combine=(length(orbits[1].path) != 0),
                        q = 1, amu = H2_amu, kwargs...)

    eo = orbit_grid.energy
    po = orbit_grid.pitch
    ro = orbit_grid.r

    nenergy = length(eo)
    npitch = length(po)
    nr = length(ro)
    norbs = length(orbits)

    e_range = extrema(eo)
    p_range = extrema(po)
    r_range = extrema(ro)
    orbs_index = zeros(Int,norbs)
    for i = 1:length(orbit_grid.orbit_index)
        ii = orbit_grid.orbit_index[i]
        ii == 0 && continue
        orbs_index[ii] != 0 && continue
        orbs_index[ii] = i
    end

    orbit_index = zeros(Int,nenergy,npitch,nr)

    norm = abs.([-(e_range...), -(p_range...), -(r_range...)])
    mins = [e_range[1],p_range[1],r_range[1]]
    oclasses = [:potato, :stagnation, :trapped, :co_passing, :ctr_passing]

    orbit_counts = Dict{Symbol,Int}(o=>count(x -> x.class == o, orbits)
                                    for o in oclasses)

    nclusters = 0
    orbs = eltype(orbits)[]
    orbit_num = 0
    for oc in oclasses
        nk = max(Int(ceil(norbits*orbit_counts[oc]/norbs)),1)
        if (nclusters + nk) > norbits
            nk = norbits - nclusters
        end
        nk == 0 && continue

        inds_oc = findall([o.class == oc for o in orbits])
        coords = hcat((([o.coordinate.energy, o.coordinate.pitch, o.coordinate.r] .- mins)./norm
                       for o in orbits if o.class == oc)...)

        if nk == 1
            if !combine
                c = coords .* norm .+ mins
                cc = EPRCoordinate(M,mean(c,dims=2)...;q = q, amu=amu)
                try
                    o = get_orbit(M, GCParticle(cc.energy,cc.pitch,cc.r,cc.z,cc.m,cc.q); kwargs...)
                    push!(orbs,o)
                    orbit_num = orbit_num + 1
                catch
                    o = Orbit(cc)
                    push!(orbs,o)
                    orbit_num = orbit_num + 1
                end
            else
                o = combine_orbits(orbits[inds_oc])
                push!(orbs,o)
                orbit_num = orbit_num + 1
            end
            orbit_index[orbs_index[inds_oc]] .= orbit_num
            nclusters = nclusters + nk
            continue
        end

        k = kmeans(coords,nk)
        if !combine
            coords = k.centers.*norm .+ mins
            for i=1:size(coords,2)
                w = k.assignments .== i
                sum(w) == 0 && continue
                cc = EPRCoordinate(M,coords[1,i],coords[2,i],coords[3,i], q=q, amu=amu)
                try
                    o = get_orbit(M, GCParticle(cc.energy,cc.pitch,cc.r,cc.z,cc.m,cc.q); kwargs...)
                    push!(orbs,o)
                    orbit_num = orbit_num + 1
                catch
                    o = Orbit(cc)
                    push!(orbs,o)
                    orbit_num = orbit_num + 1
                end
                orbit_index[orbs_index[inds_oc[w]]] .= orbit_num
            end
        else
            for i=1:nk
                w = k.assignments .== i
                sum(w) == 0 && continue
                o = combine_orbits(orbits[inds_oc[w]])
                push!(orbs,o)
                orbit_num = orbit_num + 1
                orbit_index[orbs_index[inds_oc[w]]] .= orbit_num
            end
        end
        nclusters = nclusters + nk
    end

    counts = [count(x -> x == i, orbit_index) for i=1:length(orbs)]
    return orbs, OrbitGrid(eo,po,ro,counts,orbit_index,orbit_grid.class,orbit_grid.tau_p,orbit_grid.tau_t)

end

function orbit_matrix(M::AxisymmetricEquilibrium, grid::OrbitGrid, energy, pitch, r, z; kwargs...)
    nenergy = length(energy)
    npitch = length(pitch)
    nr = length(r)
    nz = length(z)
    norbits = length(grid.counts)
    subs = CartesianIndices((nenergy,npitch,nr,nz))
    nsubs = length(subs)

    p = Progress(nsubs)
    channel = RemoteChannel(()->Channel{Bool}(nsubs), 1)
    R = fetch(@sync begin
        @async while take!(channel)
            ProgressMeter.next!(p)
        end
        @async begin
            R = @distributed (hcat) for i=1:nsubs
                ie,ip,ir,iz = Tuple(subs[i])
                gcp = GCParticle(energy[ie],pitch[ip],r[ir],z[iz])
                o = get_orbit(M,gcp;store_path=false,kwargs...)
                Rcol = spzeros(norbits)
                if !(o.class in (:lost,:incomplete)) && o.coordinate.r > M.axis[1]
                    oi = orbit_index(grid,o.coordinate)
                    (oi > 0) && (Rcol[oi] = 1.0)
                end
                put!(channel,true)
                Rcol
            end
            put!(channel,false)
            R
        end
    end)
    return R
end

function write_orbit_grid(grid::OrbitGrid;filename="orbit_grid.h5")
    h5open(filename,"w") do file
        file["energy"] = collect(grid.energy)
        file["pitch"] = collect(grid.pitch)
        file["r"] = collect(grid.r)
        file["counts"] = grid.counts
        file["orbit_index"] = grid.orbit_index
        file["class"] = String.(grid.class)
        file["tau_p"] = grid.tau_p
        file["tau_t"] = grid.tau_p
    end
    nothing
end

function read_orbit_grid(filename)
    isfile(filename) || error("File does not exist")

    f = h5open(filename)
    energy = read(f["energy"])
    pitch = read(f["pitch"])
    r = read(f["r"])
    counts = read(f["counts"])
    orbit_index = read(f["orbit_index"])
    class = Symbol.(read(f["class"]))
    tau_p = read(f["tau_p"])
    tau_t = read(f["tau_t"])
    close(f)

    return OrbitGrid(energy,pitch,r,counts,orbit_index,class,tau_p,tau_t)
end

function map_orbits(grid::OrbitGrid, f::Vector)
    if length(grid.counts) != length(f)
        throw(ArgumentError("Incompatible sizes"))
    end
    dorb = abs((grid.r[2]-grid.r[1])*(grid.energy[2]-grid.energy[1])*(grid.pitch[2]-grid.pitch[1]))
    return [i == 0 ? zero(f[1]) : f[i]/(grid.counts[i]*dorb) for i in grid.orbit_index]
end

function orbit_index(grid::OrbitGrid, o::EPRCoordinate; nearest=false)

    if !nearest
        if (grid.energy[1] <= o.energy <= grid.energy[end]) &&
           (grid.pitch[1] <= o.pitch <= grid.pitch[end]) &&
           (grid.r[1] <= o.r <= grid.r[end])

            i = argmin(abs.(o.energy .- grid.energy))
            j = argmin(abs.(o.pitch .- grid.pitch))
            k = argmin(abs.(o.r .- grid.r))
            ind = grid.orbit_index[i,j,k]
        else
            ind = 0
        end
    else
        inds = filter(x->grid.orbit_index[x] != 0, CartesianIndices(grid.orbit_index))
        data = hcat( ([grid.energy[I[1]], grid.pitch[I[2]], grid.r[I[3]]] for I in inds)...)
        tree = KDTree(data)
        idxs, dists = knn(tree,[o.energy,o.pitch,o.r],1,false)
        ind = grid.orbit_index[inds[idxs[1]]]
    end

    return ind
end

function bin_orbits(grid::OrbitGrid, orbits; weights::Union{Nothing,Vector},nearest=false)

    if weights != nothing
        length(weights) == length(orbits) || error("Incompatible weight vector size")
        w = weights
    else
        w = fill(1.0,length(orbits))
    end

    f = zeros(length(grid.counts))
    if !nearest
        for (io,o) in enumerate(orbits)
            i = argmin(abs.(o.energy .- grid.energy))
            j = argmin(abs.(o.pitch .- grid.pitch))
            k = argmin(abs.(o.r .- grid.r))
            ind = grid.orbit_index[i,j,k]
            if ind != 0
                f[ind] += w[io]
            end
        end
    else
        inds = filter(x->grid.orbit_index[x] != 0, CartesianIndices(grid.orbit_index))
        data = hcat( ([grid.energy[I[1]], grid.pitch[I[2]], grid.r[I[3]]] for I in inds)...)
        tree = KDTree(data)
        for (io, o) in enumerate(orbits)
            idxs, dists = knn(tree,[o.energy,o.pitch,o.r],1,false)
            f[grid.orbit_index[inds[idxs[1]]]] += w[io]
        end
    end

    return f
end

function combine_orbits(orbits)
    norbits = length(orbits)
    norbits == 1 && return orbits[1]

    o = orbits[1]
    r = o.path.r
    z = o.path.z
    phi = o.path.phi
    pitch = o.path.pitch
    energy = o.path.energy
    dt = o.path.dt
    #dl = o.path.dl

    c = o.coordinate
    isa(c, EPRCoordinate) || error("Wrong orbit coordinate. Expected EPRCoordinate")
    ec = c.energy
    pc = c.pitch
    rc = c.r
    zc = c.z
    tau_p = o.tau_p
    tau_t = o.tau_t

    for i=2:norbits
        oo = orbits[i]
        ec = ec + oo.coordinate.energy
        pc = pc + oo.coordinate.pitch
        rc = rc + oo.coordinate.r
        zc = zc + oo.coordinate.z
        tau_t = tau_t + oo.tau_t
        tau_p = tau_p + oo.tau_p
        append!(r, oo.path.r)
        append!(z, oo.path.z)
        append!(phi, oo.path.phi)
        append!(pitch, oo.path.pitch)
        append!(energy, oo.path.energy)
        append!(dt, oo.path.dt)
        #append!(dl, oo.path.dl)
    end
    ec = ec/norbits
    pc = pc/norbits
    rc = rc/norbits
    zc = zc/norbits
    tau_p = tau_p/norbits
    tau_t = tau_p/norbits

    cc = EPRCoordinate(ec,pc,rc,zc,zero(zc),c.m,c.q)
    path = OrbitPath(o.path.vacuum, o.path.drift, energy,pitch,r,z,phi,dt)#,dl)

    if all(x -> x.class == orbits[1].class, orbits)
        class = orbits[1].class
    else
        class = :meta
    end

    return Orbit(cc, class, tau_p, tau_t, path)
end

function mc2orbit(M::AxisymmetricEquilibrium, d::FIDASIMGuidingCenterParticles, GCP::T;  kwargs...) where T <: Function
    p = Progress(d.npart)
    channel = RemoteChannel(()->Channel{Bool}(d.npart),1)

    t = @sync begin
        @async while take!(channel)
            ProgressMeter.next!(p)
        end

        @async begin
            orbs = @distributed (vcat) for i=1:d.npart
                o = get_orbit(M,GCP(d.energy[i],M.sigma*d.pitch[i],d.r[i]/100,d.z[i]/100); kwargs...,store_path=false)
                put!(channel,true)
                o.coordinate
            end
            put!(channel,false)
            orbs
        end
    end
    return fetch(t)
end

function fbm2orbit(M::AxisymmetricEquilibrium,d::FIDASIMGuidingCenterFunction; GCP=GCDeuteron, n=1_000_000, kwargs...)
    dmc = fbm2mc(d,n=n)
    return mc2orbit(M, dmc, GCP; kwargs...)
end

struct OrbitSpline{T<:Function}
    n::Int
    itp::T
end

@inline (os::OrbitSpline)(x) = os.itp(x)
@inline Base.length(os::OrbitSpline) = os.n

function OrbitSpline(p::OrbitPath, t)
    length(p) == 0 && return OrbitSpline(0, x -> S4(zeros(4)))
    eprz = hcat(p.energy,p.pitch,p.r,p.z)
    oi = scale(interpolate(eprz, (BSpline(Cubic(Periodic(OnGrid()))), NoInterp())), t, 1:4)
    return OrbitSpline(length(t), x -> S4(oi.(x,1:4)))
end

function OrbitSpline(o::Orbit, t)
    return OrbitSpline(o.path, t)
end

function OrbitSpline(o::Orbit)
    t = range(0.0, 1.0, length=length(o))
    return OrbitSpline(o.path, t)
end

function OrbitSpline(p::OrbitPath)
    t = range(0.0, 1.0, length=length(p))
    return OrbitSpline(p, t)
end

OrbitSpline(o::OrbitSpline) = o
OrbitSpline(o::OrbitSpline, t) = o


text:
using FractionalDiffEq, Plots
h = 0.5; t = collect(0:h:1); u0 = 1
fun(t)=0
prob = DODEProblem([1, 0.1], [x->6*x*(1-x), 0], (0, 1), fun, u0, t)
sol = solve(prob, h, DOMatrixDiscrete())
plot(sol)

ext: jl


text:
module KhepriBase
using ColorTypes
using LinearAlgebra
using StaticArrays
using Dierckx
using Dates
using Sockets
using Base.Iterators
using Base: @kwdef

#=
Dierckx does not yet support derivatives. Until they update it,
we extend it here to have that support.
=#
export evaluate_derivative
function evaluate_derivative(spline::Spline2D, x::AbstractVector, y::AbstractVector, dx::Int, dy::Int)
    m = length(x)
    @assert length(y) == m
    xin = convert(Vector{Float64}, x)
    yin = convert(Vector{Float64}, y)
    zout = Vector{Float64}(undef, m)
    nx = length(spline.tx)
    ny = length(spline.ty)
    lwrk = (nx*ny)+(spline.kx+1)*m+(spline.ky+1)*m
    wrk = Vector{Float64}(undef, lwrk)
    ier = Ref{Int32}()
    kwrk = 2*m
    iwrk = Vector{Int32}(undef, kwrk)
    ccall((:pardeu_, Dierckx.libddierckx), Nothing,
          (Ref{Float64}, Ref{Int32},  # ty, ny
           Ref{Float64}, Ref{Int32},  # tx, nx
           Ref{Float64},              # c
           Ref{Int32}, Ref{Int32},    # ky, kx
           Ref{Int32}, Ref{Int32},    # nuy, nux
           Ref{Float64},              # y
           Ref{Float64},              # x
           Ref{Float64},              # z
           Ref{Int32},                # m
           Ref{Float64}, Ref{Int32},  # wrk, lwrk
           Ref{Int32}, Ref{Int32},    # iwrk, kwrk
           Ref{Int32}),               # ier
          spline.ty, ny,
          spline.tx, nx,
          spline.c,
          spline.ky, spline.kx,
          dy, dx,
          yin, xin, zout, m,
          wrk, lwrk,
          iwrk, kwrk,
          ier)
    ier[] == 0 || error(Dierckx._eval2d_message)
    return zout
end
function evaluate_derivative(spline::Spline2D, x::Number, y::Number, dx::Int, dy::Int)
  if dx >= spline.kx
    ϵ = 1e-9
    (evaluate(spline, x+ϵ, y) - evaluate(spline, x-ϵ, y))/2/ϵ
  elseif dy >= spline.kx
    ϵ = 1e-9
    (evaluate(spline, x, y+ϵ) - evaluate(spline, x, y-ϵ))/2/ϵ
  else
    evaluate_derivative(spline, [x], [y], dx, dy)[1]
  end
end

import Base:
    +, -, *, /, length,
    show, showerror,
    zero, iterate, convert,
    getindex, firstindex, lastindex, broadcastable,
    union, fill
import LinearAlgebra:
    cross, dot, norm

include("Parameters.jl")
include("Utils.jl")
include("Coords.jl")
include("Regions.jl")
include("Paths.jl")
include("Geometry.jl")
include("Backends.jl")
include("Frontend.jl")
include("Shapes.jl")
include("Primitives.jl")
include("BIM.jl")
include("Camera.jl")

export khepribase_interface_file
khepribase_interface_file() = joinpath(@__DIR__, "Interface.jl")

# From ColorTypes
export RGB, RGBA, rgb, rgba

export and_mark_deleted,
       @remote_functions,
       parse_signature,
       encode,
       decode,
       decode_or_error,
       @encode_decode_as,
       SocketBackend,
       create_backend_connection,
       current_backends,
       @remote,
       @get_remote,
       connection,
       reset_backend,
       IOBufferBackend,
       backend_frame_at,
       acad_revolution,
       backend_bounding_box,
       b_delete_shapes,
       backend_fill_curves,
       backend_loft_curve_point,
       backend_loft_curves,
       backend_loft_surface_point,
       backend_loft_surfaces,
       backend_map_division,
       backend_node_displacement_function,
       backend_panel,
       backend_revolve_curve,
       backend_revolve_point,
       backend_revolve_surface,
       backend_surface_boundary,
       backend_surface_domain,
       captured_shape,
       captured_shapes,
       create_block,
       dimension,
       generate_captured_shape,
       generate_captured_shapes,
       pre_selected_shapes_from_se,
       ref,
       save_file,
       shape_from_ref,
       show_truss_deformation,
       slice_ref

# Ports for socket-based backends are defined here to avoid conflicts
export autocad_port, revit_port, rhino_port, unity_port, unreal_port, blender_port
const autocad_port = 11000
const revit_port = 11001
const rhino_port = 12000
const unity_port = 11002
const unreal_port = 11010
const blender_port = 11003

end


text:
using GeoPhyInv
using Statistics
using Test

zgrid=range(-1000.0,stop=1000.0,length=201)
xgrid=range(-1000.0,stop=1000.0,length=201)
mgrid = [zgrid, xgrid];

vpb = [1500., 3500.] # bounds for vp
rhob = [1500., 3500.] # density bounds
medium=Medium(mgrid);
update!(medium, [:vp,:rho], [vpb,rhob]);
fill!(medium);
update!(medium, [:vp,:rho], randn_perc=5); # add some random noise

medium_new=similar(medium)
update!(medium_new, [:vp,:rho], randn_perc=5.) # add some noise
update!(medium_new, [:vp], constant_pert=500.) # perturb velocity

ageom=AGeom(medium.mgrid,:surf, SSrcs(3), Recs(30)); # surface seismic

tgrid = range(0.0,stop=1.0,length=1000) # generate a temporal grid
wav = ricker(10.0, tgrid, tpeak=0.25,); # Choose a source wavelet
srcwav = SrcWav(tgrid, ageom, [:P]) # initialize
update!(srcwav, [:P], wav) # distribute to all supersources

pa=SeisForwExpt(Fdtd(),medium=medium, ageom=ageom, srcwav=srcwav, tgrid=tgrid, verbose=true);

@time update!(pa);
d1=copy(pa[:data][:P])

update!(pa, medium_new)

@time update!(pa);
d2=copy(pa[:data][:P])

@test d1≠d2


text:
using ControlSystems
using Interpolations
using SpecialFunctions

"""
    snr, amp, k0, k1, sigma_e2 = predictSNR(ntf, R=64, amp=..., f0=0)

Predict the SNR curve of a binary delta-sigma modulator by using the describing
function method of Ardalan and Paulos.

The modulator is specified by a noise transfer function `ntf`. The band of
interest is defined by the oversampling ratio `R` and the center frequency `f0`.

The input signal is characterized by the `amp` vector. `amp` defaults to
[-120, -110, ... -20, -15, -10, -9, -8, ... 0] dB, where 0 dB means a full-scale
(peak value = 1) sine wave.

The algorithm assumes that the `amp` vector is sorted in increasing order;
once instability is detected, the remaining SNR values are set to -Inf.

Output:
- `snr`: a vector of SNR values (in dB)
- `amp`: a vector of amplitudes (in dB)
- `k0`: the quantizer signal gain
- `k1`: the quantizer noise gain
- `sigma_e2`: the power of the quantizer noise (not in dB)

The describing function method of A&P assumes that the quantizer processes
signal and noise components separately. The quantizer is modelled as two
(not necessarily equal) linear gains, `k0` and `k1`, and an additive white
gaussian noise source of power `sigma_e2`. `k0`, `k1` and `sigma_e2` are
calculated as functions of the input.

The modulator's loop filter is assumed to have nearly infinite gain at
the test frequency.
"""
function predictSNR(ntf, R=64, amp=[-120:10:-20; -15; -10:0], f0=0)
    Nb = 100
    if f0 == 0
        band_of_interest = range(0, pi/R; length=Nb)
    else
        band_of_interest = range(2*pi*(f0-0.25/R), 2*pi*(f0+0.25/R); length=Nb)
        XTAB = range(-2, 0; length=21)
        YTAB = [0.46575960516930   0.67366999387741;
                0.47904652357101   0.68426650762558;
                0.49316295981407   0.69527947902679;
                0.50817364454269   0.70673173666000;
                0.52414894104004   0.71864765882492;
                0.54116523265839   0.73105299472809;
                0.55930554866791   0.74397552013397;
                0.57866013050079   0.75744456052780;
                0.59932720661163   0.77149158716202;
                0.62141352891922   0.78615015745163;
                0.64503526687622   0.80145609378815;
                0.67031890153885   0.81744754314423;
                0.69740217924118   0.83416539430618;
                0.72643494606018   0.85165339708328;
                0.75758063793182   0.86995816230774;
                0.79101717472076   0.88912981748581;
                0.82693856954575   0.90922164916992;
                0.86555624008179   0.93029111623764;
                0.90710091590881   0.95239937305450;
                0.95182400941849   0.97561222314835;
                1.00000000000000   1.00000000000000]
        itp1 = scale(interpolate(YTAB[:,1], BSpline(Cubic(Line(OnGrid())))), XTAB)
        itp2 = scale(interpolate(YTAB[:,2], BSpline(Cubic(Line(OnGrid())))), XTAB)
    end

    ntf = tf(ntf)
    num = ntf.matrix[1].num[end:-1:0]
    den = ntf.matrix[1].den[end:-1:0]
    num1 = num - den

    N = length(amp)
    snr = zeros(N) .- Inf
    k0 = zeros(N)
    k1 = zeros(N)
    sigma_e2 = zeros(N)

    u = 10 .^ (amp/20)

    Nimp = 100
    unstable = false
    rho = 0
    fprime = 0
    m0 = 0
    for n = 1:N
        # Calculate sigma_e2
        if f0 == 0
            erfinvu = erfinv(u[n])
            sigma_e2[n] = 1 - u[n]^2 - 2/pi*exp(-2*erfinvu^2)
        else # Sinusoidal input.
            # Solve sqrt(pi)*u/2 = rho * hypergeo(0.5,2,-rho^2);
            # Formulate as solve f(rho)=0, f = rho*M(0.5,2,-rho^2)-K
            # and use the secant method.
            K = 0.5 * sqrt(pi) * u[n]
            if n == 1
                rho = u[n]^2    # Initial guess; otherwise use previous value.
                fprime = 1
            end
            drho = 1
            f_prev = 0
            for itn = 1:20
                m0 = itp2(-rho^2)
                f = rho*m0 - K
                if itn > 1
                    fprime = max((f-f_prev)/drho, 0.5)  # Secant approx.
                end
                if abs(f) < 1e-8
                    break   # Converged
                end
                drho = -f/fprime
                if abs(drho) > 0.2
                    drho = sign(drho)*0.2
                end
                if abs(drho) < 1e-6
                    break
                end
                rho = rho + drho
                f_prev = f
            end
            m1 = itp1(-rho^2)
            sigma_e2[n] = 1 - u[n]^2/2 - 2/pi*m1^2
        end

        # Iterate to solve for k1 and sigma_1.
        if n > 1
            k1[n] = k1[n-1] # Use the previous value of k1 as the initial guess.
        else
            k1[n] = 1.2
        end
        k1_prev = 0
        itn = 0
        if f0 == 0
            k1sigma1 = sqrt(2/pi)*exp(-erfinvu^2)
        else
            k1sigma1 = sqrt(2/pi)*m1
        end
        sigma_1 = 0
        while abs(k1[n]-k1_prev) > 1e-6*(1+k1[n]) && itn < 100
            # Create the function: H_hat = L1/(1-k1*L1)=(H-1)/(H*(1-k1)+k1).
            den1 = (1-k1[n])*num + den*k1[n]
            # Calculate pGain, the square of the 2-norm of H_hat.
            pGain, Nimp = powerGain(num1, den1)
            if isinf(pGain)
                unstable = true
                break
            end

            sigma_1 = sqrt(pGain*sigma_e2[n])
            k1_prev = k1[n]
            k1[n] = k1sigma1/sigma_1
            itn = itn+1
        end
        if unstable
            break
        end

        if f0 == 0
            y0 = sqrt(2)*erfinvu*sigma_1
            k0[n] = u[n]/y0
        else
            k0[n] = sqrt(2/pi)*m0/sigma_1
        end

        h = freqz(PolynomialRatio(num, (1-k1[n])*num + k1[n]*den), band_of_interest)
        # For both DC and sine wave inputs, use u^2/2 as the signal
        # power since true DC measurements are usually impossible.
        snr[n] = dbp(0.5*u[n]^2 / (sum(h.^2) / (R*Nb) * sigma_e2[n]))
    end

    return snr, amp, k0, k1, sigma_e2
end

"""
    pGain, Nimp = powerGain(num, den, Nimp0=100)

Calculate the power gain of a TF given in coefficient form. `Nimp` is the
recommended number of impulse response samples for use in future calls and
`Nimp0` is the suggested number to use.
"""
function powerGain(num, den, Nimp0=100)
    Nimp = Nimp0
    unstable = false

    sys = tf(num, den, 1)
    imp, = impulse(sys, Nimp)
    if sum(abs.(imp[Nimp-10:Nimp])) < 1e-8 && Nimp > 80  # Use fewer samples
        Nimp = round(Nimp/1.3)
    else
        while sum(abs.(imp[Nimp-10:Nimp])) > 1e-6
            Nimp = Nimp*2
            imp, = impulse(sys, Nimp)
            if sum(abs.(imp[Nimp-10:Nimp])) >= 50 || Nimp >= 1e4
                # H is close to being unstable
                unstable = true
                break
            end
        end
    end
    if !unstable
        pGain = sum(imp.^2)
    else
        pGain = Inf
    end

    return pGain, Nimp
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
include("unit.jl")
include("functionality.jl")
include("verification.jl")


text:
using DocOpt

doc="""WORLD-based F0 extraction for audio signals.
Usage:
    f0.jl [options] <src_dir> <dst_dir>
    f0.jl --version
    f0.jl -h | --help

Options:
    -h --help        show this message
    --period=PERIOD  frame period in msec [default: 5.0]
    --max=MAX        Maximum number that will be processed [default: 200]
"""

# TODO(ryuichi) allow to specify DioOption as command line options?

using VoiceConversion.Tools
using WAV

using Logging
@Logging.configure(level=DEBUG, output=STDOUT)

searchdir(path, key) = filter(x -> contains(x, key), readdir(path))

let
    args = docopt(doc, version=v"0.0.2")

    srcdir = args["<src_dir>"]
    dstdir = args["<dst_dir>"]
    mkdir_if_not_exist(dstdir)

    period = parse(Float64, args["--period"])
    nmax = parse(Int, args["--max"])

    files = searchdir(srcdir, ".wav")
    @info("$(length(files)) data found.")

    # perform feature extraction for each file
    for (count, filename) in enumerate(files)
        path = joinpath(srcdir, filename)
        savepath = joinpath(dstdir, string(splitext(basename(path))[1],
                                          "_f0.jld"))

        @info("Start processing $(path)")
        elapsed = @elapsed begin
            x, fs = wavread(path)
            size(x, 2) != 1 && error("The input data must be monoral.")
            x = map(Float64, vec(x))

            f0, timeaxis = dio(x, fs, DioOption(period=period))
            f0 = stonemask(x, fs, timeaxis, f0)

            save(savepath,
                 "description", "WORLD-based F0",
                 "type", "f0",
                 "fs", fs,
                 "period", period,
                 "feature_vector", f0
                 )
        end
        @info("Elapsed time in F0 estimation is $(elapsed) sec.")
        @info("Dumped to $(savepath)")

        count += 1
        count > nmax && break
    end

    println("Finished")
end


text:
module AnyAscii

export anyascii

const blocks = Dict{UInt16, Vector{String}}()

function anyascii(string::AbstractString)::String
	buf = IOBuffer(sizehint=div(ncodeunits(string), 2))
	for c in string
		if isascii(c)
			write(buf, UInt8(c))
		else
			print(buf, anyascii(c))
		end
	end
	String(take!(buf))
end

function anyascii(char::AbstractChar)::String
	cp = codepoint(char)
	blocknum = cp >> 8
	block = get!(blocks, blocknum) do
		file = joinpath(@__DIR__, "data", string(blocknum, base=16, pad=3))
		isfile(file) ? splitblock(read(file)) : String[]
	end
	lo = (cp & 0xff) + 1
	lo <= lastindex(block) ? block[lo] : ""
end

function splitblock(bytes::Vector{UInt8})::Vector{String}
	block = String[]
	d = 1
	while true
		i = d
		d = findnext(b -> b == 0xff, bytes, d)
		push!(block, String(bytes[i:d - 1]))
		if d == lastindex(bytes)
			break
		end
		d += 1
	end
	block
end

end

max_stars_repo_name: RastaKolobar/anyascii

max_issues_repo_name: RastaKolobar/anyascii

max_forks_repo_name: RastaKolobar/anyascii


text:
using SimpleWebsockets
using Test, Suppressor
import Sockets: getaddrinfo
include("testservers.jl")
include("testclients.jl")

@testset "SimpleWebsockets" begin
    @info "Testing SimpleWebsockets"
    @test_nowarn WebsocketServer()
    @test_nowarn WebsocketClient()

    server1 = WebsocketServer()
    client = WebsocketClient()
    @suppress begin
        @test_nowarn listen(server1, :invalidlistener) do noop end
        @test_nowarn listen(client, :invalidlistener) do noop end
    end

    @testset "Unit Tests" begin
        @info "Unit Tests"
        include("unittests.jl")
    end

    @testset "Server listens and closes" begin
        @info "Server listens and closes"
        for config in [(; ssl = false), (; ssl = true)]
            server2 = WebsocketServer(; config...)
            details = @test_nowarn servercanlisten(server2, 8080)
            @test details isa NamedTuple
            @test haskey(details, :port)
            @test haskey(details, :host)
            @test details.port === 8080
            @test details.host === getaddrinfo("localhost")
            @test !isopen(server2)
        end
    end

    @testset "Client connects and disconnects on ws and wss" begin
        @info "Client connects and disconnects on ws and wss"
        for config in [
            (; server = (; ssl = false,), client = (; url = "ws://localhost",)),
            (; server = (; ssl = true,), client = (; url = "wss://localhost",))
        ]
            server3 = WebsocketServer(; config.server...)
            @test_nowarn echoserver(server3, 8080)
            client = WebsocketClient()
            @suppress begin
                closed = @test_nowarn clientconnects(client, 8080, config.client.url)
                @test !isopen(client)
                @test closed isa NamedTuple
                @test haskey(closed, :code)
                @test haskey(closed, :description)
                @test closed.code === 1000
                @test closed.description === "Normal connection closure"
                @test length(server3) === 0
                @test_nowarn close(server3)
            end
        end
    end
    @testset "Client passes connection errors to callback handler" begin
        @info "Client passes connection errors to callback handler" wait = "wait for socket to timeout..."
        client = WebsocketClient()
        err = clientconnects(client, 8080, "ws://badurl.bad")
        @test err isa ConnectError
        @test err.msg === "Sockets.DNSError"
    end
    @testset "Client sends and receives messages up to max payload" begin
        @info "Client sends and receives messages up to max payload" (;
            description = "This incrementally scales data payload to the echo server up to the maximum allowed",
            purpose = "Scale through [8, 16, 32]bit payloads and fragmentation in text / binary combo's",
            time = "Please wait, moving a lot of test data..."
        )...
        for binary in [true, false]
            @info "Testing $(binary ? "binary" : "text") server."
            local server4 = WebsocketServer(; binary = binary)
            echoserver(server4, 8080)
            count = 0
            @sync for clientbinary in [true, false]
                count += 1
                @info "Opening $(clientbinary ? "binary" : "text") client $count"
                local client = WebsocketClient(; binary = clientbinary)
                @async begin
                    closed = echoclient(client, 8080; server4.config...)
                    @test !isopen(client)
                    @test closed.code === 1000
                    @test closed.description === "Normal connection closure"
                end
            end
            close(server4)
        end
        @info "...Done"
    end

    @testset "Server client feedback" begin
        server5 = WebsocketServer()
        echoserver(server5, 8080)

        @testset "ping pong" begin
            @info "ping pong"
            client = WebsocketClient()
            closed = pingclient(client, 8080)
            @test closed === "testping"
            @test !isopen(client)
            @test isopen(server5)
        end

        @testset "Server rejects clients with bad payloads" begin
            @info "Server rejects clients with bad payloads"
            client = WebsocketClient()
            closed = badclient(client, 8080)
            @test closed.code === 1009
            @test closed.description === "Maximum message size of 1048576 Bytes exceeded"
            @test !isopen(client)
            @test isopen(server5)
            client = WebsocketClient(; fragmentOutgoingMessages = false)
            closed = badclient(client, 8080)
            @test closed.code === 1008
            @test closed.description === "frame size exceeds maximum of 65536 Bytes."
            @test !isopen(client)
            @test isopen(server5)
        end

        @testset "Client times out" begin
            @info "Client times out"
            client = WebsocketClient()
            closed = timeoutclient(client, 8080)
            @test closed.code === 1006
            @test closed.description === "could not ping the server."
            @test !isopen(client)
            @test isopen(server5)
            close(server5)
        end

        server6 = WebsocketServer()
        echoserver(server6, 8080)

        @testset "server gracefully closes clients and exits" begin
            @info "server gracefully closes clients and exits"
            client = WebsocketClient()
            client2 = WebsocketClient()
            local reason, reason2
            @async reason = wsclient(client)
            @async reason2 = wsclient(client2)
            while length(server6) < 2
                sleep(0.1)
            end
            @test length(server6) === 2
            close(server6)

            @test !isopen(client)
            @test !isopen(client2)
            @test !isopen(server6)
            for r in [reason, reason2]
                @test r.code === 1001
                @test r.description === "Remote peer is going away"
            end
        end
    end
end

max_stars_repo_name: citkane/SimpleWebsockets

max_issues_repo_name: citkane/SimpleWebsockets

max_forks_repo_name: citkane/SimpleWebsockets


text:
    # MPT type
    mutable struct MPT{N}
        cores::Vector{Array{Float64,N}}
        normcore::Int64
        function MPT(cores::Vector{Array{Float64,N}},normcore::Int64) where N
            new{ndims(cores[1])}(cores,normcore)
        end
    end
    MPT(cores) = MPT(cores,0);

    # aliases for MPS and MPO
    const MPS = MPT{3};
    const MPO = MPT{4};

    # middle sizes of each core
    function Base.size(mpt::MPT)
        sizes = [size(core) for core in mpt.cores];
        [sizes[i][2:end-1] for i in 1:length(sizes)]
    end

    # number of cores
    function order(mpt::MPT)
        collect(size(mpt.cores))[1]
    end

    # number of elements in each core and overall
    function Base.length(mpt::MPT)
        lengths = [length(core) for core in mpt.cores];
        [lengths,[sum(lengths)]]
    end

    # ranks of each core
    function LinearAlgebra.rank(mpt::MPT)
        sizes = [size(core) for core in mpt.cores];
        [[sizes[i][1], sizes[i][end]] for i in 1:length(sizes)]
    end

    function LinearAlgebra.norm(mps::MPS)
        N   = order(mps);
        mptnorm = contractcores(mps[1],mps[1]);
        for i = 2:N
            mptnorm = mptnorm*contractcores(mps[i],mps[i]);
        end
        return sqrt(mptnorm[1])
    end

    function LinearAlgebra.norm(mpo::MPO)
        mps = mpo2mps(mpo);
        norm(mps)
    end

    # indexing to get and set a core in an MPT
    Base.IndexStyle(::Type{<:MPT}) = IndexLinear()
    Base.getindex(mpt::MPT, i::Int) = mpt.cores[i] # tt[1] gives the first core
    Base.getindex(mpt::MPT, range::UnitRange{Int64}) = [mpt.cores[i] for i in range]
    Base.setindex!(mpt::MPT,v,i::Int) = setindex!(mpt.cores, v, i) # tt[1] = rand(1,5,3) sets the first core

    # indexing to get an element or row/column of the reconstructed vector or matrix

    function getelement(tt::MPT{3},row::Int)
        sz = size(tt);
        d  = ();
        el = 1;
        for i = 1:length(sz)
            d = (d...,1:sz[i][1])
        end
        linear = LinearIndices(d);
        index = findall(x->x==row,linear)[1];
        for i = 1:length(sz)
            el = el*tt[i][:,index[i],:];
        end
        return el[1]
    end

    function getelement(ttm::MPT{4},row::Int,col::Int)
        sz = size(ttm);
        d1  = ();
        d2  = ();
        el = 1;
        for i = 1:length(sz)
            d1 = (d1...,1:sz[i][1])
            d2 = (d2...,1:sz[i][2])
        end
        linear1 = LinearIndices(d1);
        index1 = findall(x->x==row,linear1)[1];
        linear2 = LinearIndices(d2);
        index2 = findall(x->x==col,linear2)[1];
        for i = 1:length(sz)
            el = el*ttm[i][:,index1[i],index2[i],:];
        end
        return el[1]
    end

    function getcol(ttm::MPT{4},col::Int)
        sz = size(ttm);
        d  = ();
        for i = 1:length(sz)
            d = (d...,1:sz[i][1])
        end
        linear = LinearIndices(d);
        index = findall(x->x==col,linear)[1];

        cores = Vector{Array{Float64,3}}(undef,order(ttm));
        for i = 1:length(sz)
            cores[i] = ttm[i][:,:,index[i],:];
        end
        return MPT(cores)
    end

    function getrow(ttm::MPT{4},row::Int)
        sz = size(ttm);
        d  = ();
        for i = 1:length(sz)
            d = (d...,1:sz[i][2])
        end
        linear = LinearIndices(d);
        index = findall(x->x==row,linear)[1];

        cores = Vector{Array{Float64,3}}(undef,order(ttm));
        for i = 1:length(sz)
            cores[i] = ttm[i][:,index[i],:,:];
        end
        return MPT(cores)
    end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
module ValueGraphFromDataFrame

using SimpleValueGraphs
using DataFrames

export ValDiGraph

function vattr_types_from_df(vertex_df)
    vatypes = eltype.(eachcol(vertex_df))
    vanames = names(vertex_df)
    (; (Symbol.(vanames) .=> vatypes)...)
end

function make_edgelist_nodeidx(edgelist, vertex_df, vertex_df_col=1)
    vind = DataFrame(name=vertex_df[:,vertex_df_col], idx=1:size(vertex_df, 1))
    elnames = Symbol.(names(edgelist, [1,2]))
    res = leftjoin(edgelist, vind, on=elnames[1] => :name)
    rename!(res, :idx => :idx1)
    res = leftjoin(res, vind, on=elnames[2] => :name)
    select!(res, :idx1, :idx, Not(elnames))
    rename!(res, :idx1 => elnames[1], :idx => elnames[2])
    return res
end

function add_vertices_from_df!(vg, vertex_df)
    for r in eachrow(vertex_df)
       add_vertex!(vg, NamedTuple(r))
    end
end

function add_edges_from_df!(vg, edgelist)
    if size(edgelist, 2) > 2
        for r in eachrow(edgelist)
            add_edge!(vg, r[1], r[2], NamedTuple(r[Not([1,2])]))
        end
    else
        for r in eachrow(edgelist)
            add_edge!(vg, r[1], r[2])
        end
    end
end

function SimpleValueGraphs.ValDiGraph(edgelist::DataFrame, vertex_df::DataFrame; vertex_df_col=1)
    vg = ValDiGraph(0, vertexval_types=vattr_types_from_df(vertex_df), vertexval_init=undef)
    el = make_edgelist_nodeidx(edgelist, vertex_df, vertex_df_col)
    add_vertices_from_df!(vg, vertex_df)
    add_edges_from_df!(vg, el)
    return vg
end

end


text:
module CUTENSOR

using ..CuArrays
using ..CuArrays: active_context, @argout, @workspace

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

using CEnum
const cudaDataType_t = cudaDataType

const libcutensor = Ref("libcutensor")

# core library
include("libcutensor_common.jl")
include("error.jl")
include("libcutensor.jl")

# low-level wrappers
include("tensor.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

const _handles = Dict{CuContext,Ref{cutensorHandle_t}}()
const _handle = Ref{Union{Ref{cutensorHandle_t},Nothing}}(nothing)

function handle()
    if _handle[] == nothing
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _handle[] = get!(_handles, active_context[]) do
            context = active_context[]
            handle = Ref{cutensorHandle_t}()
            cutensorInit(handle)
            handle
        end
    end
    return _handle[]
end

end


text:
using Test
import GDAL
using DiskArrays: eachchunk, haschunks, Chunked, GridChunks, readblock!
import ArchGDAL; const AG = ArchGDAL

@testset "RasterDataset Type" begin
    AG.readraster("ospy/data4/aster.img") do ds
        @testset "Test forwarded methods" begin
            @test ds isa AG.RasterDataset{UInt8}
            @test AG.getgeotransform(ds) == [419976.5, 15.0, 0.0, 4.6624225e6, 0.0, -15.0]
            @test AG.nraster(ds) == 3
            @test AG.getband(ds,1) isa AG.AbstractRasterBand
            @test startswith(AG.getproj(ds),"PROJCS")
            @test AG.width(ds)==5665
            @test AG.height(ds)==5033
            @test AG.getdriver(ds) isa AG.Driver
            @test splitpath(AG.filelist(ds)[1]) == ["ospy", "data4", "aster.img"]
            @test splitpath(AG.filelist(ds)[2]) == ["ospy", "data4", "aster.rrd"]
            @test AG.listcapability(ds) isa Dict
            @test AG.ngcp(ds)==0
            @test AG.write(ds,tempname()) == C_NULL
            @test AG.testcapability(ds,"ODsCCreateLayer") == false
        end
        @testset "DiskArray chunk interface" begin
            b = AG.getband(ds,1)
            @test eachchunk(ds) == GridChunks(size(ds),(64,64,1))
            @test eachchunk(b) == GridChunks(size(b),(64,64))
            @test haschunks(ds) == Chunked()
            @test haschunks(b) == Chunked()
        end
        @testset "Reading into non-arrays" begin
            data1 = view(zeros(3,3,3), 1:3, 1:3, 1:3)
            readblock!(ds, data1, 1:3, 1:3, 1:3)
            @test data1 == ds[1:3,1:3,1:3]
        end
    end
end

@testset "Test Array getindex" begin
    AG.readraster("ospy/data4/aster.img") do ds
        @testset "Dataset indexing" begin
            @testset "dims dropped correctly" begin
                @test typeof(ds[:, :, :]) <: Array{UInt8,3}
                @test typeof(ds[:, :, 1]) <: Array{UInt8,2}
                @test typeof(ds[:, 1, 1]) <: Array{UInt8,1}
                @test typeof(ds[1, 1, 1]) <: UInt8
            end
            @testset "range indexing" begin
                buffer = ds[1:AG.width(ds), 1:AG.height(ds), 1:1]
                @test typeof(buffer) <: Array{UInt8,3}
                total = sum(buffer)
                count = sum(buffer .> 0)
                @test total / count ≈ 76.33891347095299
                @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
            end
            @testset "colon indexing" begin
                buffer = ds[:, :, 1]
                total = sum(buffer)
                count = sum(buffer .> 0)
                @test total / count ≈ 76.33891347095299
                @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
            end
            @testset "int indexing" begin
                @test ds[755, 2107, 1] == 0xff
            end
        end
        @testset "RasterBand indexing" begin
            band = AG.getband(ds, 1)
            @testset "dims dropped correctly" begin
                @test typeof(band[:, :]) <: Array{UInt8,2}
                @test typeof(band[:, 1]) <: Array{UInt8,1}
                @test typeof(band[1, 1]) <: UInt8
            end
            @testset "range indexing" begin
                buffer = band[1:AG.width(band), 1:AG.height(band)]
                @test typeof(buffer) <: Array{UInt8,2}
                total = sum(buffer)
                count = sum(buffer .> 0)
                @test total / count ≈ 76.33891347095299
                @test total / (AG.height(band) * AG.width(band)) ≈ 47.55674749653172
            end
            @testset "colon indexing" begin
                buffer = band[:, :]
                total = sum(buffer)
                count = sum(buffer .> 0)
                @test total / count ≈ 76.33891347095299
                @test total / (AG.height(band) * AG.width(band)) ≈ 47.55674749653172
            end
            @testset "int indexing" begin
                @test band[755, 2107] == 0xff
            end
        end
    end
end

cp("ospy/data4/aster.img", "ospy/data4/aster_write.img"; force=true)

@testset "Test Array setindex" begin
    AG.readraster("ospy/data4/aster_write.img"; flags=AG.OF_Update) do ds
        @testset "Dataset setindex" begin
            @test ds[755, 2107, 1] == 0xff
            ds[755, 2107, 1] = 0x00
            @test ds[755, 2107, 1] == 0x00
            ds[755:755, 2107:2107, 1:1] = reshape([0x01], 1, 1, 1)
            @test ds[755, 2107, 1] == 0x01
            ds[755:755, 2107:2107, 1] = reshape([0x02], 1, 1)
            @test ds[755, 2107, 1] == 0x02
            ds[755:755, 2107, 1] = [0x03]
            @test ds[755, 2107, 1] == 0x03
        end
        @testset "RasterBand setindex" begin
            band = AG.getband(ds, 1)
            band[755:755, 2107] = [0x00]
            @test band[755, 2107] == 0x00
            band[755, 2107] = 0x01
            @test band[755, 2107] == 0x01
            band[755:755, 2107:2107] = reshape([0xff], 1, 1)
            @test band[755, 2107] == 0xff
            buffer = band[:, :]
            band[:, :] = buffer .* 0x00
            @test sum(band[:, :]) == 0x00
            band[:, 1:500] = buffer[:, 1:500]
            band[:, 501:end] = buffer[:, 501:end]
            @test sum(buffer) / sum(buffer .> 0) ≈ 76.33891347095299
            @test_throws DimensionMismatch band[:, 501:end] = [1, 2, 3]
        end
    end
end

@testset "Test Array constructor" begin
    AG.readraster("ospy/data4/aster_write.img"; flags=AG.OF_Update) do ds
        buffer = Array(ds)
        typeof(buffer) <: Array{UInt8,3}
        total = sum(buffer[:, :, 1:1])
        count = sum(buffer[:, :, 1:1] .> 0)
        @test total / count ≈ 76.33891347095299
        @test total / (AG.height(ds) * AG.width(ds)) ≈ 47.55674749653172
    end
end


text:
function img2colConvolve(cLayer::CL, Ai::AbstractArray{T,N}) where {T,N, CL <: ConvLayer}

    Aip = padding(cLayer, Ai)
    axBT = axes(cLayer.B)
    axB = axBT[1:end-2]
    axBend = axBT[end]
    Z = col2img(cLayer.K * img2col(Aip), cLayer.outputS) .+
            view(cLayer.B, axB..., 1, axBend)

    actFun = cLayer.actFun
    Ao = eval(:($actFun($Z)))

    return Dict(:Z => Z, :A => Ao)
end #function img2colConvolve(cLayer::CL

export img2colConvolve

### dimg2colConvolve!

export dimg2colConvolve!

function dimg2colConvolve!(
    cLayer::Conv1D,
    Ai::AbstractArray{T1,3},
    dAi::AbstractArray{T2,3},
    dZ::AbstractArray{T3,3},
) where {T1,T2,T3}

    padS = paddingSize(cLayer, Ai)
    Aip = padding(cLayer, Ai)
    Aipv = img2col(Aip)
    dZv = img2col(dZ)

    n_Hi, c, m = size(Aip)
    s_H = cLayer.s
    f_H = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1

    dAi .= col2img(cLayer.K' * dZv, (n_Hi,c,m))[1+padS[1]:end-padS[2], :, :]
    cLayer.dK = dZv * Aipv'
    cLayer.dB = sum(permutedims(dZ, [1, 3, 2]), dims = 1:2)

    return nothing
end # function dimg2colConvolve!(
    #     cLayer::Conv1D,
    #     Ai::AbstractArray{T,3},
    #     dA::AbstractArray{T,3},
    #     dZ::AbstractArray{T,3},
    # ) where {T}

function dimg2colConvolve!(
    cLayer::Conv2D,
    Ai::AbstractArray{T1,4},
    dAi::AbstractArray{T2,4},
    dZ::AbstractArray{T3,4},
) where {T1,T2,T3}

    padS = paddingSize(cLayer, Ai)
    Aip = padding(cLayer, Ai)
    Aipv = img2col(Aip)
    dZv = img2col(dZ)
    n_Hi, n_Wi, c, m = size(Aip)
    s_H, s_W = cLayer.s
    f_H, f_W = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1

    dAi .= col2img(cLayer.K' * dZv, (n_Hi,n_Wi,c,m))[1+padS[1]:end-padS[2], 1+padS[3]:end-padS[4], :, :]
    cLayer.dK = dZv * Aipv'
    cLayer.dB = sum(permutedims(dZ, [1, 2, 4, 3]), dims = 1:3)

    return nothing
end # function dimg2colConvolve!(
    #     cLayer::Conv2D,
    #     Ai::AbstractArray{T,4},
    #     dA::AbstractArray{T,4},
    #     dZ::AbstractArray{T,4},
    # ) where {T}

function dimg2colConvolve!(
    cLayer::Conv3D,
    Ai::AbstractArray{T1,5},
    dAi::AbstractArray{T2,5},
    dZ::AbstractArray{T3,5},
) where {T1,T2,T3}

    padS = paddingSize(cLayer, Ai)
    Aip = padding(cLayer, Ai)
    Aipv = img2col(Aip)
    dZv = img2col(dZ)
    n_Hi, n_Wi, n_Di, c, m = size(Aip)
    s_H, s_W, s_D = cLayer.s
    f_H, f_W, f_D = cLayer.f
    n_H = (n_Hi - f_H) ÷ s_H + 1
    n_W = (n_Wi - f_W) ÷ s_W + 1
    n_D = (n_Di - f_D) ÷ s_D + 1

    dAi .= col2img(cLayer.K' * dZv, (n_Hi, n_Wi, n_Di, c, m))[1+padS[1]:end-padS[2],
        1+padS[3]:end-padS[4],
        1+padS[5]:end-padS[6],
        :,
        :,
    ]
    cLayer.dK = dZv * Aipv'
    cLayer.dB = sum(permutedims(dZ, [1, 2, 3, 5, 4]), dims = 1:4)

    return nothing
end # function dimg2colConvolve!(
    #     cLayer::Conv3D,
    #     Ai::AbstractArray{T,5},
    #     dA::AbstractArray{T,5},
    #     dZ::AbstractArray{T,5},
    # ) where {T}


text:
## ExoplanetsSysSims/src/abc_distance.jl
## (c) 2015 Eric B. Ford

dist_L1_fractional(x::Real, y::Real) =  (x==y==0.0) ? 0.0 : abs( (x-y)*2/(x+y) )
dist_L1_abs(x::Real, y::Real) =  abs(x-y)
dist_L2_fractional(x::Real, y::Real) =  (x==y==0.0) ? 0.0 : abs( (x-y)^2*2/(x^2+y^2) )
dist_L2_abs(x::Real, y::Real) =  (x-y)^2

# Library of functions to convert distance vector into a scalar distance
calc_scalar_distance_sum(d::Vector{Float64}) = sum(d)
calc_scalar_distance_rms(d::Vector{Float64}) = sqrt(sumsq(d))
calc_scalar_distance_max(d::Vector{Float64}) = maximum(d)

# How to combine two distances based on subsets of the full distance vector (should always be greater than either for two pass algorithm to work)
combine_scalar_distances_sum(d1::Float64, d2::Float64) = d1+d2
combine_scalar_distances_rms(d1::Float64, d2::Float64) = sqrt(d1*d1+d2*d2)
combine_scalar_distances_max(d1::Float64, d2::Float64) = max(d1,d2)

# Pick which one to be used
calc_scalar_distance(d::Vector{Float64}) = calc_scalar_distance_sum(d)
combine_scalar_distances(d1::Float64, d2::Float64) = combine_scalar_distances_sum(d1,d2)

# compute supremum of differences between empirical cdfs.
# Borrowed from JuliaStats/HypothesisTests.jl
function ksstats(x::AbstractVector{T}, y::AbstractVector{S}) where {T<:Real, S<:Real}
    n_x, n_y = length(x), length(y)
    sort_idx = sortperm([x; y], alg=MergeSort)
    pdf_diffs = [ones(n_x)/n_x; -ones(n_y)/n_y][sort_idx]
    cdf_diffs = cumsum(pdf_diffs)
    deltap = maximum(cdf_diffs)
    deltan = -minimum(cdf_diffs)
    delta = max(deltap,deltan)
    (n_x, n_y, deltap, deltan, delta)
end
# weighted version   # WARNING:  Function is untested
function ksstats(x::AbstractVector{T}, y::AbstractVector{S}, wx::AbstractVector{T}, wy::AbstractVector{T}) where {T<:Real, S<:Real}
    n_x, n_y = length(x), length(y)
    wx .*= 1.0/sum(wx)
    wy .*= 1.0/sum(wy)
    sort_idx = sortperm([x; y], alg=MergeSort)
    pdf_diffs = [wx; -wy][sort_idx]
    cdf_diffs = cumsum(pdf_diffs)
    deltap = maximum(cdf_diffs)
    deltan = -minimum(cdf_diffs)
    delta = max(deltap,deltan)
    (n_x, n_y, deltap, deltan, delta)  # should the first two values returned here be sum(wx) and sum(wy) before normalizing?  For now, we ignore all but the final delta anyway.
end

dist_KS(x::AbstractVector{T}, y::AbstractVector{S}) where {T<:Real, S<:Real}  = ksstats(x,y)[5]
dist_KS(x::AbstractVector{T}, y::AbstractVector{S}, wx::AbstractVector{T}, wy::AbstractVector{T}) where {T<:Real, S<:Real} = ksstats(x,y,wx,wy)[5]

##### Function to compute the KS distance for discrete arrays (i.e. of the planet multiplicities in the systems):
function ksstats_ints(x::AbstractVector{T}, y::AbstractVector{S}) where {T <: Int, S <: Int}
    n_x, n_y = length(x), length(y)
    xy_max = max(maximum(x), maximum(y)) #maximum of both x and y
    x_counts = [sum(x .== n) for n in 1:xy_max]
    y_counts = [sum(y .== n) for n in 1:xy_max]
    pdf_diffs = x_counts./n_x .- y_counts./n_y
    cdf_diffs = cumsum(pdf_diffs)
    deltap = maximum(cdf_diffs)
    deltan = -minimum(cdf_diffs)
    delta = max(deltap,deltan)
    (n_x, n_y, deltap, deltan, delta)
end

##### Function to compute the 2-sample Anderson-Darling (AD) distance between two empirical CDFs (continuous distributions):
function ADstats(x::AbstractVector{T}, y::AbstractVector{S}) where {T <: Real, S <: Real}
    #This function computes the AD distance according to A. N. Pettitt (1976) Eq. (1.2)
    n, m = length(x), length(y)
    N = n + m
    sort_idx = sortperm([x; y], alg=MergeSort) #array of indices that would sort the combined array
    M_i_diffs = [ones(n); zeros(m)][sort_idx]
    M_i_array = cumsum(M_i_diffs)[1:end-1] #array of M_i except for last element, i.e. from i=1 to i=N-1
    i_array = 1:(N-1) #array of i from i=1 to i=N-1
    AD_dist = (1/(n*m))*sum(((M_i_array*N .- n*i_array).^2)./(i_array.*(N .- i_array))) #AD distance
    return AD_dist
end

function ADstats_mod(x::AbstractVector{T}, y::AbstractVector{S}) where {T <: Real, S <: Real}
    #This function is the same as 'ADstats' except without the factor of 'nm/N' before the integral
    n, m = length(x), length(y)
    N = n + m
    sort_idx = sortperm([x; y], alg=MergeSort) #array of indices that would sort the combined array
    M_i_diffs = [ones(n); zeros(m)][sort_idx]
    M_i_array = cumsum(M_i_diffs)[1:end-1] #array of M_i except for last element, i.e. from i=1 to i=N-1
    i_array = 1:(N-1) #array of i from i=1 to i=N-1
    AD_dist = (N/((n*m)^2))*sum(((M_i_array*N .- n*i_array) .^2) ./(i_array .*(N .- i_array))) #AD distance
    return AD_dist
end

##### Function to compute the Cressie-Read Power Divergence (CRPD) statistic for the distributions of observed planet multiciplities:
function CRPDstats(En::AbstractVector{T}, On::AbstractVector{S}) where {T <: Int, S <: Int}
    #En and On must be arrays of the total numbers of systems with 1,2,3,... observed planets, in the simulated (i.e. expected) and the actual (i.e. observed Kepler) data, respectively
    n_max = max(length(En), length(On))
    En = append!(En, zeros(n_max - length(En)))
    On = append!(On, zeros(n_max - length(On)))

    E_array = En./sum(En) #normalized numbers (fractions) of simulated systems with 1,2,3,... observed planets
    O_array = On./sum(On) #normalized numbers (fractions) of actual Kepler systems with 1,2,3,... observed planets
    rho = 0.
    for (i,E_i) in enumerate(E_array)
        if En[i] != 0
            rho += O_array[i]*((O_array[i]/E_array[i])^(2/3) - 1)
        end
    end
    rho = (9/5)*rho
    return rho
end

# lambda:  rate for Poisson process, i.e., expected value for number of events
# k:  number of events observed
function calc_num_events_to_maximize_poisson_pdf(lambda::Real)
  kstar = floor(Int64,lambda)
  if lambda > kstar+1
     kstar +=1
  end
  return kstar
end

function log_prob_poisson_num_events_given_rate(lambda::Real,k::Real)
   k*log(lambda)-lambda-lgamma(k+1)
end

function delta_log_prob_poisson_num_events_given_rate(lambda::Real, k::Integer, kstar::Integer = calc_num_events_to_maximize_poisson_pdf(lambda) )
  delta = (k-kstar)*log(lambda)
  if kstar>=k+1
	delta += sum(log.((k+1):kstar))
  elseif k>=kstar+1
     delta -= sum(log.((kstar+1):k))
  end
  return delta
end

function distance_poisson_likelihood(lambda::Real, k::Integer)
  kstar = calc_num_events_to_maximize_poisson_pdf(lambda)
  khi = round(Int64,lambda+sqrt(lambda))
  delta_logprob_one_sigma_hi = delta_log_prob_poisson_num_events_given_rate(lambda,khi,kstar)
  delta_logp_k = delta_log_prob_poisson_num_events_given_rate(lambda,k,kstar)
  delta_logp_k/delta_logprob_one_sigma_hi
end

function distance_poisson_draw(lambda::Real, k::Integer)
  d =  Distributions.Poisson(lambda)
  simulated_number_of_detections = rand(d)
  abs( simulated_number_of_detections -k)
end

"""
    distance_sum_of_bernoulli_draws(num_pl_obs, num_targets_obs, prob_detect_list, num_targets_sim, bin_match_list)

Calculate L2 distance and the number of simulated planets detected for a binned occurrence rate where each planet's detection is drawn from a Bernoulli distribution using the planet's detection probability.

# Arguments:
- `num_pl_obs::Integer`: Number of planets in (true) observed catalog
- `num_targets_obs::Integer`: Number of target stars in (true) observed catalog
- `prob_detect_list::Vector{Float}`: Detection probabilities for planets in simulated catalog
- `num_targets_sim::Integer`: Number of target stars in simulated catalog
- `bin_match_list::Vector{Integer}`: List of simulated catalog planet indices associated with current occurrence rate bin

# Returns:
- L2 distance between provided (true) observed catalog and simulated catalog summary statistics
- Number of detected simulated planets
"""
function distance_sum_of_bernoulli_draws(num_pl_obs::Integer, num_targets_obs::Integer, prob_detect_list::Vector{TReal}, num_targets_sim::Integer, bin_match_list::Vector{TInt}) where {TReal<:Real, TInt<:Integer}
   @assert(0<=num_pl_obs<=num_targets_obs)
   num_pl_match = length(bin_match_list)
   @assert(0<=length(bin_match_list))

   num_detect_sim = 0
   if num_pl_match >= 1
      num_draws_all = min(max(1,floor(Int64, num_targets_obs/num_targets_sim)),1000)
      @assert(1<=num_draws_all<=1000)
      for i in 1:num_pl_match
         pl_id = bin_match_list[i]
         @assert(1<=pl_id<=length(prob_detect_list))
         prob_detect = min(prob_detect_list[pl_id],1.0)
         num_detect_sim += sum(rand(Bernoulli(prob_detect),num_draws_all))
      end
      # If number of targets observed is not a multiple of number of targets simulated, then pick a random set to make total number of draws equal (as long as there are some planets to choose from)
      for i in (num_pl_match*num_draws_all+1):(floor(Int64,num_pl_match*num_targets_obs/num_targets_sim))
          pl_id = bin_match_list[rand(1:num_pl_match)]
          prob_detect = min(prob_detect_list[pl_id],1.0)
          num_detect_sim += rand(Bernoulli(prob_detect))
      end
   end
   distance = dist_L2_abs(num_pl_obs/num_targets_obs, num_detect_sim/num_targets_obs)
   return distance, num_detect_sim
end

# compute Canberra distance.
function distance_canberra_modified(x::AbstractVector{T}, y::AbstractVector{S}) where {T<:Real, S<:Real}
    @assert length(x) == length(y)
    dist_sum = 0.0
    for i in 1:length(x)
        numer = abs(x[i]-y[i])
        denom = sqrt(x[i] + y[i])
        if denom == 0.0
            continue
        else
            dist_sum += numer/denom
        end
    end
    return dist_sum
end

# compute original Canberra distance.
function distance_canberra_orig(x::AbstractVector{T}, y::AbstractVector{S}) where {T<:Real, S<:Real}
    @assert length(x) == length(y)
    dist_sum = 0.0
    for i in 1:length(x)
        numer = abs(x[i]-y[i])
        denom = abs(x[i]) + abs(y[i])
        if denom == 0.0
            continue
        else
            dist_sum += numer/denom
        end
    end
    return dist_sum
end

distance_canberra(x::AbstractVector{T}, y::AbstractVector{S}) where {T<:Real, S<:Real} = distance_canberra_modified(x,y)  # WARNING: Defaults to modified Canberra distance.  At some point, update Danley's code to call distance_canberra_modified, so can remove this

# compute Cosine distance.
function distance_cosine(x::AbstractVector{T}, y::AbstractVector{S}) where {T<:Real, S<:Real}
    @assert length(x) == length(y)
    numer = 0.0
    denom_1 = 0.0
    denom_2 = 0.0
    for i in 1:length(x)
        numer += x[i]*y[i]
        denom_1 += x[i]^2
        denom_2 += y[i]^2
    end
    return numer/(sqrt(denom_1)*sqrt(denom_2))
end

function wasserstein_distance(u_values::AbstractVector{T1}, v_values::AbstractVector{T2}, p::Integer) where { T1<:Real, T2<:Real}
  # code adapted from SciPy.stats._cdf_distance
  u_sorter = issorted(u_values) ? (1:length(u_values)) : sortperm(u_values)
  v_sorter = issorted(v_values) ? (1:length(v_values)) : sortperm(v_values)
  all_values = sort(vcat(u_values,v_values),alg=MergeSort)
  deltas = all_values[2:end] .- all_values[1:end-1]
  u_cdf = map(x->searchsortedlast(u_values[u_sorter],x)/ length(u_values),all_values[1:end-1])
  v_cdf = map(x->searchsortedlast(v_values[v_sorter],x)/ length(v_values),all_values[1:end-1])
  if p==1
     return sum(abs.(u_cdf.-v_cdf).*deltas)
    elseif p==2
      return sqrt(sum((u_cdf.-v_cdf).^2 .* deltas))
    else
      return pow(sum(abs.(u_cdf.-v_cdf).^p .* deltas),1/p)
    end
end

earth_mover_distance(u_values::AbstractVector{T1}, v_values::AbstractVector{T2}) where { T1<:Real, T2<:Real} = wasserstein_distance(u_values,v_values,1)

# TODO USER SCI: IMPORTANT:  Replace the distance function with something well thought out for your particular scientific application.  See examples
function calc_distance_vector_demo(summary1::CatalogSummaryStatistics, summary2::CatalogSummaryStatistics, pass::Int64, sim_param::SimParam ; verbose::Bool = false)
  d = Array{Float64}(undef,0)
  if pass == 1
    if verbose
      println("# Summary 1, pass 1: ",summary1)
      println("# Summary 2, pass 1: ",summary2)
    end
    d = zeros(3)
    # Since observed and simulated catalogs can have different summary statistics for the number of planets, prefer detections if avaliable (e.g., after pass2), otherwise use expected (e.g., from pass 1)
    np1 = haskey(summary1.stat,"planets detected") ? summary1.stat["planets detected"] : summary1.stat["expected planets detected"]
    np2 = haskey(summary2.stat,"planets detected") ? summary2.stat["planets detected"] : summary2.stat["expected planets detected"]
    d[1] = dist_L1_abs(np1/summary1.stat["num targets"],np2/summary2.stat["num targets"])    #  Normalize so different statistics weighted appropriately and not dominated by this one
    #println("np1 = ",np1,", np2 = ",np2)
    #println("np1 (normalized) = ",np1/summary1.stat["num targets"],", np2 (normalized) = ",np2/summary2.stat["num targets"],", d[1] = ",d[1])
    #d[2] = dist_L1_abs(summary1.stat["mean log10 P"],summary2.stat["mean log10 P"])
    #d[3] = dist_L1_abs(summary1.stat["mean log10 depth"],summary2.stat["mean log10 depth"])
    #d[4] = dist_L1_abs(summary1.stat["std log10 P"],summary2.stat["std log10 P"])
    #d[5] = dist_L1_abs(summary1.stat["std log10 depth"],summary2.stat["std log10 depth"])
    #d[2] = dist_KS(summary1.stat["P list"], summary2.stat["P list"])
    #d[3] = dist_KS(summary1.stat["depth list"], summary2.stat["depth list"])
    # EDITS for christiansen-single-bin
    #d[2] = dist_KS(summary1.stat["P list"], summary2.stat["P list"],summary1.stat["weight list"],summary2.stat["weight list"])
    #d[3] = dist_KS(summary1.stat["depth list"], summary2.stat["depth list"],summary1.stat["weight list"],summary2.stat["weight list"])
    # END EDITS
  elseif pass == 2
    max_tranets_in_sys = get_int(sim_param,"max_tranets_in_sys")
    d = zeros(max_tranets_in_sys)
    for n in 1:max_tranets_in_sys
      d[n] = n*dist_L1_abs(summary1.stat["num_sys_tranets"][n]/summary1.stat["num targets"],summary2.stat["num_sys_tranets"][n]/summary2.stat["num targets"])
    end
  else
    println("# calc_distance_vector_demo doesn't know what to do for pass= ", pass)
  end
  return d
end

function test_abc_distance(cat_obs::KeplerObsCatalog, cat_phys::KeplerPhysicalCatalog, sim_param::SimParam)
  ss_pass1 = calc_summary_stats_sim_pass_one_demo(cat_obs,cat_phys,sim_param)
  ss_pass2 = calc_summary_stats_sim_pass_two_demo(cat_obs,cat_phys,ss_pass1,sim_param)
  d1 = calc_distance_vector_demo(ss_pass1,ss_pass1, 1, sim_param)
  d2 = calc_distance_vector_demo(ss_pass2,ss_pass2, 2, sim_param)
end


text:
#########
# UTILS #
#########
listed_deps(project::Project) =
    append!(collect(keys(project.deps)), collect(keys(project.extras)))

###########
# READING #
###########
read_project_uuid(::Nothing) = nothing
function read_project_uuid(uuid::String)
    try uuid = UUID(uuid)
    catch err
        err isa ArgumentError || rethrow()
        pkgerror("Could not parse project UUID as a UUID")
    end
    return uuid
end
read_project_uuid(uuid) = pkgerror("Expected project UUID to be a string")

read_project_version(::Nothing) = nothing
function read_project_version(version::String)
    try version = VersionNumber(version)
    catch err
        err isa ArgumentError || rethrow()
        pkgerror("Could not parse project version as a version")
    end
end
read_project_version(version) = pkgerror("Expected project version to be a string")

read_project_deps(::Nothing, section::String) = Dict{String,UUID}()
function read_project_deps(raw::Dict{String,Any}, section_name::String)
    deps = Dict{String,UUID}()
    for (name, uuid) in raw
        try
            uuid = UUID(uuid)
        catch err
            err isa ArgumentError || rethrow()
            pkgerror("Malfomed value for `$name` in `$(section_name)` section.")
        end
        deps[name] = uuid
    end
    return deps
end
function read_project_deps(raw, section_name::String)
    pkgerror("Expected `$(section_name)` section to be a key-value list")
end

read_project_targets(::Nothing, project::Project) = Dict{String,Vector{String}}()
function read_project_targets(raw::Dict{String,Any}, project::Project)
    for (target, deps) in raw
        deps isa Vector{String} || pkgerror("""
            Expected value for target `$target` to be a list of dependency names.
        """)
    end
    return raw
end
read_project_targets(raw, project::Project) =
    pkgerror("Expected `targets` section to be a key-value list")

read_project_compat(::Nothing, project::Project) = Dict{String,String}()
function read_project_compat(raw::Dict{String,Any}, project::Project)
    for (name, version) in raw
        try VersionSpec(semver_spec(version))
        catch err
            pkgerror("Could not parse compatibility version for dependency `$name`")
        end
    end
    return raw
end
read_project_compat(raw, project::Project) =
    pkgerror("Expected `compat` section to be a key-value list")

function validate(project::Project)
    # deps
    dep_uuids = collect(values(project.deps))
    if length(dep_uuids) != length(unique(dep_uuids))
        pkgerror("Two different dependencies can not have the same uuid")
    end
    # extras
    extra_uuids = collect(values(project.extras))
    if length(extra_uuids) != length(unique(extra_uuids))
        pkgerror("Two different `extra` dependencies can not have the same uuid")
    end
    dep_names = keys(project.deps)
    # TODO decide what to do in when `add`ing a dep that is already in `extras`
    #   also, reintroduce test files for this
    #=
    for (name, uuid) in project.extras
        name in dep_names && pkgerror("name `$name` is listed in both `deps` and `extras`")
        uuid in dep_uuids && pkgerror("uuid `$uuid` is listed in both `deps` and `extras`")
    end
    =#
    # targets
    listed = listed_deps(project)
    for (target, deps) in project.targets, dep in deps
        if length(deps) != length(unique(deps))
            pkgerror("A dependency was named twice in target `$target`")
        end
        dep in listed || pkgerror("""
            Dependency `$dep` in target `$target` not listed in `deps` or `extras` section.
            """)
    end
    # compat
    for (name, version) in project.compat
        name == "julia" && continue
        name in listed ||
            pkgerror("Compat `$name` not listed in `deps` or `extras` section.")
    end
end

function Project(raw::Dict)
    project = Project()
    project.other    = raw
    project.name     = get(raw, "name", nothing)
    project.manifest = get(raw, "manifest", nothing)
    project.uuid     = read_project_uuid(get(raw, "uuid", nothing))
    project.version  = read_project_version(get(raw, "version", nothing))
    project.deps     = read_project_deps(get(raw, "deps", nothing), "deps")
    project.extras   = read_project_deps(get(raw, "extras", nothing), "extras")
    project.compat   = read_project_compat(get(raw, "compat", nothing), project)
    project.targets  = read_project_targets(get(raw, "targets", nothing), project)
    validate(project)
    return project
end

function read_project(io::IO; path=nothing)
    raw = nothing
    try
        raw = TOML.parse(io)
    catch err
        if err isa TOML.ParserError
            pkgerror("Could not parse project $(something(path,"")): $(err.msg)")
        elseif all(x -> x isa TOML.ParserError, err)
            pkgerror("Could not parse project $(something(path,"")): $err")
        else
            rethrow()
        end
    end
    return Project(raw)
end

read_project(path::String) =
    isfile(path) ? open(io->read_project(io;path=path), path) : Project()

###########
# WRITING #
###########
function destructure(project::Project)::Dict
    raw = deepcopy(project.other)

    should_delete(x::Dict) = isempty(x)
    should_delete(x)       = x === nothing
    entry!(key::String, src) = should_delete(src) ? delete!(raw, key) : (raw[key] = src)

    entry!("name",     project.name)
    entry!("uuid",     project.uuid)
    entry!("version",  project.version)
    entry!("manifest", project.manifest)
    entry!("deps",     project.deps)
    entry!("extras",   project.extras)
    entry!("compat",   project.compat)
    entry!("targets",  project.targets)
    return raw
end

_project_key_order = ["name", "uuid", "keywords", "license", "desc", "deps", "compat"]
project_key_order(key::String) =
    something(findfirst(x -> x == key, _project_key_order), length(_project_key_order) + 1)

write_project(project::Project, project_file::AbstractString) =
    write_project(destructure(project), project_file)
write_project(project::Project, io::IO) =
    write_project(destructure(project), io)
write_project(project::Dict, project_file::AbstractString) =
    open(io -> write_project(project, io), project_file; truncate=true)
write_project(project::Dict, io::IO) =
    TOML.print(io, project, sorted=true, by=key -> (project_key_order(key), key))
function write_project(project::Project, env, old_env, ctx::Context; display_diff=true)
    project = destructure(ctx.env.project)
    if !isempty(project) || ispath(env.project_file)
        if display_diff && !(ctx.currently_running_target)
            printpkgstyle(ctx, :Updating, pathrepr(env.project_file))
            Pkg.Display.print_project_diff(ctx, old_env, env)
        end
        if !ctx.preview
            mkpath(dirname(env.project_file))
            write_project(project, env.project_file)
        end
    end
end


text:
x = [1,2,3]
function add!(x...)
    r = x[1]
    s = x[3]
    a = sum(r, s)
    return a
end
add(x...)

function get_big_matrix_item(big_m)

    item = big_m[3,3]
    return item

end

function get_little_matrix_item(lil_m)

    item = lil_m[2,1]
    return item
end

function add!(a, b, c, x...)

    a[1] = b + get_big_matrix_item(x[1])
    a[2] = c + get_little_matrix_item(x[2])

end

big_m = rand(Int8, 5,5)
lil_m = rand(Int8, 2,2)

a = [2,3]
y = [big_m, lil_m]
add!
a

# note:
# haven't really solved problem - it runs properly, but I don't think its
# actually pulling the objects properly, I need to explicitly pass the arguments
# to the function, but I'm not sure if the ODE problem will let me do that
# that is, I need to put something in add!(...) but I'm not sure if the ODEproblem
# will accept that?

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
# Inspired by
# * https://github.com/fonsp/pluto-on-binder/blob/compiler/postBuild
# * https://github.com/bat/batjl-environment/blob/master/build_sysimage.jl

import Pkg

srcdir = @__DIR__
Pkg.activate(srcdir)

@info("Building custom system image in project \"$srcdir\".")

@info("Instantiating project.")
Pkg.instantiate()

@info("Precompiling packages.")
Pkg.precompile()

proj = Pkg.API.project();
pkgs = Symbol.(keys(proj.dependencies))
proj.ispackage && push!(pkgs, Symbol(proj.name))

@info("Package to include in system image: $(join(pkgs, " "))")

@info("Building system image.")

import PackageCompiler, Libdl

sysimage_path = joinpath(srcdir, "JuliaSysimage." * Libdl.dlext)

PackageCompiler.create_sysimage(
    Symbol.(pkgs),
    sysimage_path = sysimage_path,
    precompile_execution_file = joinpath(srcdir, "pkgcompile_exec.jl"),
    cpu_target = PackageCompiler.default_app_cpu_target(),
    replace_default = false
)

@info("Created custom Julia system image \"$sysimage_path\".")

max_stars_repo_name: greimel/pluto-binder

max_issues_repo_name: greimel/pluto-binder

max_forks_repo_name: greimel/pluto-binder


text:
using ApproxFun, SingularIntegralEquations, Test

@testset "GreensFun" begin
    cts = 100
    Γ,ψ = ones(cts),ones(cts)
    [(Γ[k],ψ[k]) = (k*Γ[k-1],ψ[k-1]+1/k) for k=2:cts]

    function FK1(z)
        if abs(z) <=.25
            upd = 1.
            ret = upd
            l=1
            zt2 = (.5z)^2
            while norm(upd) >= eps()
                upd = zt2^l/(Γ[l]*Γ[l+1])
                ret += upd
                l+=1
            end
        else
            ret = 2besseli(1,abs(z))/abs(z)
        end
        .5ret
    end
    function GK1(z)
        if abs(z) <= .25
            upd = 1.
            ret = upd
            l=1
            zt2 = (.5z)^2
            while norm(upd) >= eps()
                upd = zt2^l/(Γ[l]*Γ[l+1])*(ψ[l]+ψ[l+1])
                ret += upd
                l+=1
            end
        else
            abz = abs(z)
            ret = 4/abz^2 - 4(besselk(1,abz)-besseli(1,abz)*(log(abz/2)+γ))/abz
        end
        .5ret
    end

    d = Segment(-1.,1.)
    sp = Space(d)
    wsp = JacobiWeight(-0.5,-0.5,sp)
    ⨍ = DefiniteLineIntegral(wsp)
    x = Fun(identity,d)

    β = -2.0
    G = exp(β*x)

    Kf = (x,y)->(1/π*exp(β/2*(y-x))*β^2/2^2*(FK1(β*(y-x)/2)*(log(abs(β)/4) + γ) - GK1(β*(y-x)/2)/2))/sqrt(1-y^2)
    K = GreensFun(Kf,sp⊗wsp)
    @test K(0.1,0.2) ≈ Kf(0.1,0.2)

    K0f = (x,y)-> x == y ? abs(β)/(4sqrt(1-y^2)) : exp(β/2*(y-x))*abs(β)/2*besseli(1,abs(β*(y-x)/2))/abs(y-x)/sqrt(1-y^2)
    K0 = GreensFun(K0f,CauchyWeight(sp⊗wsp,0))
    @test K0(0.1,0.2) ≈ K0f(0.1,0.2)*logabs(0.1-0.2)/π

    K2f = (x,y)->exp(β/2*(y-x))/sqrt(1-y^2)
    K2 = GreensFun(K2f,CauchyWeight(sp⊗wsp,2))
    @test K2(0.1,0.2) ≈ K2f(0.1,0.2)/(0.1-0.2)^2/π
end


text:
@testset "Partitioning" begin
  setify(lists) = Set(Set.(lists))

  d = CartesianGrid{T}(10,10)
  p = partition(d, UniformPartition(100))
  @test sprint(show, p) == "100 Partition{2,$T}"
  @test sprint(show, MIME"text/plain"(), p) == "100 Partition{2,$T}\n  └─1 View{10×10 CartesianGrid{2,$T}}\n  └─1 View{10×10 CartesianGrid{2,$T}}\n  └─1 View{10×10 CartesianGrid{2,$T}}\n  └─1 View{10×10 CartesianGrid{2,$T}}\n  └─1 View{10×10 CartesianGrid{2,$T}}\n  ⋮\n  └─1 View{10×10 CartesianGrid{2,$T}}\n  └─1 View{10×10 CartesianGrid{2,$T}}\n  └─1 View{10×10 CartesianGrid{2,$T}}\n  └─1 View{10×10 CartesianGrid{2,$T}}\n  └─1 View{10×10 CartesianGrid{2,$T}}"

  @testset "UniformPartition" begin
    Random.seed!(123)

    grid = CartesianGrid{T}(3,3)
    p = partition(grid, UniformPartition(3, false))
    @test setify(indices(p)) == setify([[1,2,3], [4,5,6], [7,8,9]])
    p = partition(grid, UniformPartition(3))
    if VERSION < v"1.7"
      @test setify(indices(p)) == setify([[8,6,9], [4,1,7], [2,3,5]])
    else
      @test setify(indices(p)) == setify([[4,6,1], [9,8,3], [5,7,2]])
    end

    grid = CartesianGrid{T}(2,3)
    p = partition(grid, UniformPartition(3, false))
    @test setify(indices(p)) == setify([[1,2], [3,4], [5,6]])
  end

  @testset "DirectionPartition" begin
    grid = CartesianGrid{T}(3,3)

    # basic checks on small regular grid data
    p = partition(grid, DirectionPartition(T.((1,0))))
    @test setify(indices(p)) == setify([[1,2,3], [4,5,6], [7,8,9]])

    p = partition(grid, DirectionPartition(T.((0,1))))
    @test setify(indices(p)) == setify([[1,4,7], [2,5,8], [3,6,9]])

    p = partition(grid, DirectionPartition(T.((1,1))))
    @test setify(indices(p)) == setify([[1,5,9], [2,6], [3], [4,8], [7]])

    p = partition(grid, DirectionPartition(T.((1,-1))))
    @test setify(indices(p)) == setify([[1], [2,4], [3,5,7], [6,8], [9]])

    # opposite directions produce same partition
    dir1 = (rand(T), rand(T)); dir2 = .-dir1
    p1 = partition(grid, DirectionPartition(dir1))
    p2 = partition(grid, DirectionPartition(dir2))
    @test setify(indices(p1)) == setify(indices(p2))

    # partition of arbitrarily large regular grid always
    # returns the "lines" and "columns" of the grid
    for n in [10,100,200]
      grid = CartesianGrid{T}(n,n)

      p = partition(grid, DirectionPartition(T.((1,0))))
      @test setify(indices(p)) == setify([collect((i-1)*n+1:i*n) for i in 1:n])
      ns = [nelements(d) for d in p]
      @test all(ns .== n)

      p = partition(grid, DirectionPartition(T.((0,1))))
      @test setify(indices(p)) == setify([collect(i:n:n*n) for i in 1:n])
      ns = [nelements(d) for d in p]
      @test all(ns .== n)
    end
  end

  @testset "FractionPartition" begin
    grid = CartesianGrid{T}(10,10)

    p = partition(grid, FractionPartition(T(0.5)))
    @test nelements(p[1]) == nelements(p[2]) == 50
    @test length(p) == 2

    p = partition(grid, FractionPartition(T(0.7)))
    @test nelements(p[1]) == 70
    @test nelements(p[2]) == 30

    p = partition(grid, FractionPartition(T(0.3)))
    @test nelements(p[1]) == 30
    @test nelements(p[2]) == 70
  end

  @testset "BlockPartition" begin
    grid = CartesianGrid{T}(10,10)

    p = partition(grid, BlockPartition(T(5),T(5)))
    @test length(p) == 4
    @test all(nelements.(p) .== 25)

    p = partition(grid, BlockPartition(T(5),T(2)))
    @test length(p) == 12
    @test Set(nelements.(p)) == Set([5,10])

    grid = CartesianGrid{T}(50, 50, 50)

    p = partition(grid, BlockPartition(T(1.), T(1.), T(1.), neighbors = false))
    @test length(p) == 125000
    @test Set(nelements.(p)) == Set(1)
    @test metadata(p) == Dict{Any,Any}()

    p = partition(grid, BlockPartition(T(5.), T(5.), T(5.), neighbors = true))
    @test length(p) == 1000
    @test Set(nelements.(p)) == Set(125)
    n = metadata(p)[:neighbors]
    @test length(n) == length(p)
    @test all(0 .< length.(n) .< 27)
  end

  @testset "BisectPointPartition" begin
    grid = CartesianGrid((10,10), T.((-0.5,-0.5)), T.((1.0, 1.0)))

    p = partition(grid, BisectPointPartition(T.((0.0,1.0)), T.((5.0,5.1))))
    p1, p2 = p[1], p[2]
    @test nelements(p1) == 60
    @test nelements(p2) == 40

    # all points in p1 are below those in p2
    pts1 = [centroid(p1, i) for i in 1:nelements(p1)]
    pts2 = [centroid(p2, i) for i in 1:nelements(p2)]
    X1 = reduce(hcat, coordinates.(pts1))
    X2 = reduce(hcat, coordinates.(pts2))
    M1 = maximum(X1, dims=2)
    m2 = minimum(X2, dims=2)
    @test all(X1[2,j] < m2[2] for j in 1:size(X1,2))
    @test all(X2[2,j] > M1[2] for j in 1:size(X2,2))

    # flipping normal direction is equivalent to swapping subsets
    p₁ = partition(grid, BisectPointPartition(T.(( 1.0,0.0)), T.((5.1,5.0))))
    p₂ = partition(grid, BisectPointPartition(T.((-1.0,0.0)), T.((5.1,5.0))))
    @test nelements(p₁[1]) == nelements(p₂[2]) == 60
    @test nelements(p₁[2]) == nelements(p₂[1]) == 40
  end

  @testset "BisectFractionPartition" begin
    grid = CartesianGrid((10,10), T.((-0.5,-0.5)), T.((1.0,1.0)))

    p = partition(grid, BisectFractionPartition(T.((1.0,0.0)), T(0.2)))
    p1, p2 = p[1], p[2]
    @test nelements(p1) == 20
    @test nelements(p2) == 80

    # all points in p1 are to the left of p2
    pts1 = [centroid(p1, i) for i in 1:nelements(p1)]
    pts2 = [centroid(p2, i) for i in 1:nelements(p2)]
    X1 = reduce(hcat, coordinates.(pts1))
    X2 = reduce(hcat, coordinates.(pts2))
    M1 = maximum(X1, dims=2)
    m2 = minimum(X2, dims=2)
    @test all(X1[1,j] < m2[1] for j in 1:size(X1,2))
    @test all(X2[1,j] > M1[1] for j in 1:size(X2,2))

    # flipping normal direction is equivalent to swapping subsets
    p₁ = partition(grid, BisectFractionPartition(T.(( 1.0,0.0)), T(0.2)))
    p₂ = partition(grid, BisectFractionPartition(T.((-1.0,0.0)), T(0.8)))
    @test nelements(p₁[1]) == nelements(p₂[2]) == 20
    @test nelements(p₁[2]) == nelements(p₂[1]) == 80
  end

  @testset "BallPartition" begin
    pset = PointSet(T[
      0 1 1 0 0.2
      0 0 1 1 0.2
    ])

    # 3 balls with 1 point, and 1 ball with 2 points
    p = partition(pset, BallPartition(T(0.5)))
    n = nelements.(p)
    @test length(p) == 4
    @test count(i->i==1, n) == 3
    @test count(i->i==2, n) == 1
    @test setify(indices(p)) == setify([[1,5],[2],[3],[4]])

    # 5 balls with 1 point each
    p = partition(pset, BallPartition(T(0.2)))
    @test length(p) == 5
    @test all(nelements.(p) .== 1)
    @test setify(indices(p)) == setify([[1],[2],[3],[4],[5]])
  end

  @testset "PlanePartition" begin
    grid = CartesianGrid((3,3), T.((-0.5,-0.5)), T.((1.0,1.0)))
    p = partition(grid, PlanePartition(T.((0,1))))
    @test setify(indices(p)) == setify([[1,2,3],[4,5,6],[7,8,9]])

    grid = CartesianGrid((4,4), T.((-0.5,-0.5)), T.((1.0,1.0)))
    p = partition(grid, PlanePartition(T.((0,1))))
    @test setify(indices(p)) == setify([1:4,5:8,9:12,13:16])
  end

  @testset "PredicatePartition" begin
    grid = CartesianGrid((3,3), T.((-0.5,-0.5)), T.((1.0,1.0)))

    # partition even from odd locations
    pred(i,j) = iseven(i+j)
    p = partition(grid, PredicatePartition(pred))
    @test setify(indices(p)) == setify([1:2:9,2:2:8])
  end

  @testset "SpatialPredicatePartition" begin
    g = CartesianGrid((10,10), T.((-0.5,-0.5)), T.((1.0,1.0)))

    # check if there are 100 partitions, each one having only 1 point
    sp = SpatialPredicatePartition((x,y) -> norm(x-y) < T(1))
    s = indices(partition(g, sp))
    @test length(s) == 100
    nelms = [nelements(d) for d in partition(g, sp)]
    @test all(nelms .== 1)

    # defining a predicate to check if points x and y belong to the square [0.,5.]x[0.,5.]
    pred(x, y) = all(T[0,0] .<= x .<= T[5,5]) && all(T[0,0] .<= y .<= T[5,5])
    sp = SpatialPredicatePartition(pred)
    p = partition(g, sp)
    s = indices(p)
    n = nelements.(p)

    # There will be 65 subsets:
    # 1 subset with 36 points (inside square [0.,5.]x[0.,5.])
    # 64 subsets with only 1 point inside each of them
    @test length(s) == 65
    @test maximum(length.(s)) == 36
    @test count(i->i==1, n) == 64
    @test count(i->i==36, n) == 1
  end

  @testset "ProductPartition" begin
    g = CartesianGrid((100,100), T.((-0.5,-0.5)), T.((1.0,1.0)))
    bm  = BlockPartition(T(10), T(10))
    bn  = BlockPartition(T(5), T(5))
    bmn = ProductPartition(bm, bn)

    # Bm x Bn = Bn with m > n
    s1 = indices(partition(g, bmn))
    s2 = indices(partition(g, bn))
    @test setify(s1) == setify(s2)

    # pXp=p (for deterministic p)
    for p in [BlockPartition(T(10), T(10)),
              BisectFractionPartition(T.((0.1,0.1)))]
      pp = ProductPartition(p, p)
      s1 = indices(partition(g, pp))
      s2 = indices(partition(g, p))
      @test setify(s1) == setify(s2)
    end
  end

  @testset "HierarchicalPartition" begin
    g = CartesianGrid((100,100), T.((-0.5,-0.5)), T.((1.0,1.0)))
    bm  = BlockPartition(T(10), T(10))
    bn  = BlockPartition(T(5), T(5))
    bmn = HierarchicalPartition(bm, bn)

    # Bn -> Bm = Bm with m > n
    s1 = indices(partition(g, bmn))
    s2 = indices(partition(g, bn))
    @test setify(s1) == setify(s2)
  end

  @testset "Mixed Tests" begin
    g = CartesianGrid((100,100), T.((-0.5,-0.5)), T.((1.0,1.0)))
    bm  = BlockPartition(T(10), T(10))
    bn  = BlockPartition(T(5), T(5))
    bmn = ProductPartition(bm, bn)
    hmn = HierarchicalPartition(bm, bn)

    # Bm*Bn = Bm->Bn
    s1 = indices(partition(g, bmn))
    s2 = indices(partition(g, hmn))
    @test setify(s1) == setify(s2)
  end

  @testset "Utilities" begin
    d = CartesianGrid{T}(10,10)
    l, r = split(d, T(0.5))
    @test nelements(l) == 50
    @test nelements(r) == 50
    l, r = split(d, T(0.5), T.((1,0)))
    @test nelements(l) == 50
    @test nelements(r) == 50
    lpts = [centroid(l, i) for i in 1:nelements(l)]
    rpts = [centroid(r, i) for i in 1:nelements(r)]
    cl = mean(coordinates.(lpts))
    cr = mean(coordinates.(rpts))
    @test cl[1] < cr[1]
    @test cl[2] == cr[2]
    l, r = split(d, T(0.5), T.((0,1)))
    @test nelements(l) == 50
    @test nelements(r) == 50
    lpts = [centroid(l, i) for i in 1:nelements(l)]
    rpts = [centroid(r, i) for i in 1:nelements(r)]
    cl = mean(coordinates.(lpts))
    cr = mean(coordinates.(rpts))
    @test cl[1] == cr[1]
    @test cl[2] < cr[2]

    d = CartesianGrid{T}(10,10)
    s = slice(d, T(0.5):T(11), T(0.5):T(11))
    @test s == CartesianGrid(P2(1,1), P2(10,10), dims=(9,9))
  end
end


text:
function jacobian_batched(f, x::AbstractMatrix{T}, array_mover::Function)::Tuple where {T <: AbstractFloat}
    y, back = Zygote.pullback(f, x)
    z = zeros(eltype(x), size(x)) |> array_mover
    res = zeros(size(x, 1), size(x, 1), size(x, 2)) |> array_mover
    for i in 1:size(y, 1)
        z[i, :] .= one(eltype(x))
        res[i, :, :] .= only(back(z))
        z[i, :] .= zero(eltype(x))
    end
    y, res
end

function make_mover(acceleration::AbstractResource, data_type::DataType)
    (data_type <: AbstractFloat) || error("data_type must be a float type")

    if acceleration isa CUDALibs
        arr_t = CuArray
    else
        arr_t = Array
    end

    function mover(x)
        convert(arr_t{data_type}, x)
    end
    mover
end


text:
module Light

using CthulhuVision.Math

export Ray, RGB, origin, direction, pointat

struct RGB
    r::Float32
    g::Float32
    b::Float32
end

@inline Base.:*(t::Float32, rgb::RGB) :: RGB = RGB(t*rgb.r, t*rgb.g, t*rgb.b)
@inline Base.:*(a::RGB, b::RGB) :: RGB = RGB(a.r * b.r, a.g * b.g, a.b * b.b)
@inline Base.:/(rgb::RGB, t::Float32) :: RGB = RGB(rgb.r/t, rgb.g/t, rgb.b/t)
@inline Base.:+(a::RGB, b::RGB) :: RGB = RGB(a.r + b.r, a.g + b.g, a.b + b.b)

struct Ray
    a::Vector3
    b::Vector3
end

@inline origin(r::Ray) = r.a
@inline direction(r::Ray) = r.b
@inline pointat(r::Ray, t::Float32) = r.a + t*r.b

end