text:
function augment_images(im,ll)

    total_augmentations = 6

    im_out = zeros(Float32,size(im,1),size(im,2),size(im,3),size(im,4)*total_augmentations)
    l_out = zeros(Float32,size(ll,1),size(ll,2),size(ll,3),size(ll,4)*total_augmentations)

    count=1

    for i=1:size(im,4)
        im_out[:,:,:,count] = im[:,:,:,i]
        l_out[:,:,:,count] = ll[:,:,:,i]
        count += 1
        #rotations
        for j in [pi/12, pi/6, -pi/12, -pi/6]
            for k=1:size(im,3)
                im_out[:,:,k,count] = imrotate(im[:,:,k,i],j,Reflect())[1:size(im,1),1:size(im,2)]
            end
            for k=1:size(ll,3)
                l_out[:,:,k,count] = imrotate(ll[:,:,k,i],j,Reflect())[1:size(ll,1),1:size(ll,2)]
            end
            count += 1
        end

        #additive noise
        im_out[:,:,:,count] = im[:,:,:,i] .+ 0.1*rand(Float32,size(im_out,1),size(im_out,2),size(im_out,3))
        l_out[:,:,:,count] = ll[:,:,:,i]
        count += 1

        #Set group of pixels to random value x times
    end

    #Shuffle Positions
    myinds=Random.shuffle(collect(1:size(im_out,4)));
    im_out=im_out[:,:,:,myinds]
    l_out=l_out[:,:,:,myinds]

    (im_out,l_out)
end

function ffmpeg_cmd(st,vid_name,num_frames,temp_file)
`$(FFMPEG.ffmpeg) -y -loglevel panic -ss $(st) -i $(vid_name) -frames:v $(num_frames) -f rawvideo -pix_fmt gray $(temp_file)`
end

#=
Convert Discrete points to heatmap for deep learning
=#
function make_heatmap_labels(han::Tracker_Handles,real_w=640,real_h=480,label_img_size=64)
    make_heatmap_labels(han.woi,han.wt.pad_pos,real_w,real_h,label_img_size)
end

function make_heatmap_labels(woi,pad_pos,real_w=640,real_h=480,label_img_size=64)

    d_points=make_discrete_all_whiskers(woi,pad_pos)

    labels=zeros(Float32,label_img_size,label_img_size,size(d_points,1),size(woi,1))

    for i=1:size(labels,4)
        for j=1:size(labels,3)
            this_x = d_points[j,1,i] / real_w * label_img_size
            this_y = d_points[j,2,i] / real_h * label_img_size

            if (this_x !=0.0)&(this_y != 0.0)
                labels[:,:,j,i] = StackedHourglass.gaussian_2d(1.0:1.0:label_img_size,1.0:1.0:label_img_size,this_y,this_x)
                labels[:,:,j,i] = labels[:,:,j,i] ./ maximum(labels[:,:,j,i])
            end
        end
    end

    labels
end

function get_labeled_frames(han::Tracker_Handles,out_hw=256)
    get_labeled_frames(han.wt.vid_name,han.frame_list,out_hw=256)
end

function get_labeled_frames(vid_name,frame_list,out_hw=256)

    (w,h,frame_rate)=get_vid_dims(vid_name)

    imgs=zeros(Float32,out_hw,out_hw,1,length(frame_list))

    temp=zeros(UInt8,w,h)
    for i=1:length(frame_list)
        frame_time = frame_list[i] / frame_rate
        load_single_frame(frame_time,temp,vid_name)
        imgs[:,:,1,i]=Images.imresize(temp',(out_hw,out_hw))
    end
    imgs
end


text:
using Distributions, FastGaussQuadrature

struct Stable <: ContinuousUnivariateDistribution
  α::Float64
  β::Float64
  θ::Float64
  ρ::Float64
  Stable(α::Real,β::Real) = ( β < -1. || β > 1. || α <= 0. || α > 2. ) ?
    error("Parameters' requirements unmet:\n (α,β)∈(0,2]×[-1,1]") :
    α == 2. ? new(2.,0.,0.,.5) :
    new(α,β,β*(α <= 1. ? 1. :
    (α-2.)/α), (1. + β*(α <= 1. ? 1. : (α-2.)/α))/2.)
end

import Distributions.params
function params(d::Stable)
  return (d.α,d.β,d.θ,d.ρ)
end

import Distributions.minimum
function minimum(d::Stable)
  return d.α <= 1. && d.β == 1. ? (d.α == 1. ? 1. : 0.) : -Inf
end

import Distributions.maximum
function maximum(d::Stable)
  return d.α <= 1. && d.β == -1. ? (d.α == 1. ? -1. : 0.) : Inf
end

import Distributions.insupport
function insupport(d::Stable,x::Real)
  return d.α <= 1. && abs(d.β) == 1. ? (d.α == 1. ? x == d.β : sign(x) != -d.β) : isfinite(x)
end

########################
# BASIC SIMULATION TOOLS
########################

import Distributions.rand
function rand(d::Stable) # Chambers. Mellows, and Stuck
  U = rand(Uniform(-pi/2,pi/2))
  return rand(Exponential())^(1-1/d.α) * sin(d.α * (U+(pi/2)*d.θ)) /
    (cos(U) ^ (1/d.α) * cos(U - d.α * (U+(pi/2)*d.θ)) ^ (1-1/d.α))
end

function rand(d::Stable,n::Int) # Chambers. Mellows, and Stuck
  U = rand(Uniform(-pi/2,pi/2),n)
  return rand(Exponential(),n).^(1-1/d.α) .* sin.(d.α .* (U .+ (pi/2)*d.θ)) ./
    (cos.(U) .^ (1/d.α) .* cos.(U .- d.α .* (U .+ (pi/2)*d.θ)) .^ (1-1/d.α))
end
########################
# PDF, CDF
########################

function auxV2(x::AbstractArray,a::Real,th::Real)
  y = (pi/2).*x
  t = (pi/2)*a*th
  return ((sin(a.*y+t).^a)./cos(y)).^(1/(1-a)).*cos((a-1).*y+t)
end

import Distributions.pdf
function pdf(d::Stable,x::Real)
  return x==0 ? gamma(1+1/d.α)*sin(pi*d.ρ)/pi : (x>0 ? pdf(PositiveStable(d.α,d.β),x)*d.ρ : pdf(PositiveStable(d.α,-d.β),-x)*(1-d.ρ))
end

function pdf(d::Stable,X::AbstractArray{Real})
  return pdf(PositiveStable(d.α,d.β),X).*d.ρ + pdf(PositiveStable(d.α,-d.β),-X).*(1-d.ρ)
end

import Distributions.cdf
function cdf(d::Stable,x::Real)
  return d.ρ + cdf(PositiveStable(d.α,d.β),x)*d.ρ - cdf(PositiveStable(d.α,-d.β),-x)*(1-d.ρ)
end

function cdf(d::Stable,X::AbstractArray)
  return d.ρ + cdf(PositiveStable(d.α,d.β),X).*d.ρ - cdf(PositiveStable(d.α,-d.β),-X).*(1-d.ρ)
end

import Distributions.cf
function cf(d::Stable,x::Real)
  if x==0
    return 1
  end
  if d.α==1 && abs(d.β)==1
    return exp(d.β*im*x)
  end
  return exp(-abs(x)^d.α*exp(-pi*d.θ*d.α*sign(x)*im/2))
end

function cf(d::Stable,X::AbstractArray{Real})
  if d.α==1 && abs(d.β)==1
    return exp((d.β*im).*X)
  end
  x = vec(X)
  s1 = exp(-pi*d.θ*d.α*im/2)
  s2 = exp(pi*d.θ*d.α*im/2)
  return reshape(exp(-abs(x).^d.α.*[x[i]>0 ? s1 : s2 for i=1:length(x)]),size(X))
end

import Distributions.mgf
function mgf(d::Stable,x::Real)
  if x==0
    return 1
  end
  if d.α==2
    return exp(x^2)
  elseif d.α==1 && abs(d.β)==1
      return exp(d.β*x)
  elseif x>0
    if d.β==-1 && d.α!=1
        return exp(sign(d.α-1)*x^d.α)
    else
      return Inf
    end
  elseif d.β==1 && d.α!=1
    return exp(sign(d.α-1)*(-x)^d.α)
  end
  return Inf
end

function mgf(d::Stable,X::AbstractArray{Real})
  if d.α==2
    return exp(X.^2)
  elseif d.α==1 && abs(d.β)==1
    return exp(d.β.*X)
  end
  x = vec(X)
  res = exp(sign(d.α-1).*abs(x).^d.α)
  if d.β==-1
    return reshape([x[i]>=0 ? res[i] : Inf for i=1:length(x)],size(X))
  elseif d.β==1
    return reshape([x[i]<=0 ? res[i] : Inf for i=1:length(x)],size(X))
  end
  return reshape([x[i]==0 ? 1 : Inf for i=1:length(x)],size(X))
end

import Distributions.mean
function mean(d::Stable)
  if d.α<=1
    return Inf
  end
  cr=1-d.ρ
  return (sin(pi*d.ρ)-sin(pi*cr))/(d.α*sin(pi/d.α)*gamma(1+1/d.α))
end

import Distributions.var
function var(d::Stable)
  if d.α == 1. && abs(d.β) == 1.
      return 0.
  else
    return Inf
  end
end

function mellin(d::Stable,x::Complex)
  if (real(x) >= d.α && (d.α <= 1. || (d.α < 2. && d.β != -1.))) || real(x) <= -1.
    return Inf
  end
  if (d.α > 1. && d.β == -1.) || d.α == 2.
    return d.ρ * gamma(1. + x) / gamma(1. + x / d.α)
  end
  return d.ρ * (sin(pi * d.ρ * x) * gamma(1. + x)) /
    (d.α * d.ρ * sin(pi * x / d.α) * gamma(1. + x / d.α))
end

function mellin(d::Stable,X::AbstractArray{Complex})
  if (d.α > 1. && d.β == -1.) || d.α == 2.
    return d.ρ .* gamma.(1. .+ X) ./ gamma.(1. .+ X ./ d.α)
  end

  res = Complex{Float64}[]
  for x = X
    push!(res,
      (real(x) >= d.α && (d.α <= 1. || (d.α < 2. && d.β != -1.))) || real(x) <= -1. ?
      Inf : (sin(pi * d.ρ * x) * gamma(1. + x)) / (d.α * d.ρ * sin(pi * x / d.α) * gamma(1. + x / d.α))
    )
  end
  return reshape(d.ρ .* res,size(X))
end

function mellin(d::Stable,x::Real)
  if (real(x) >= d.α && (d.α <= 1. || (d.α < 2. && d.β != -1.))) || real(x) <= -1.
    return Inf
  end
  if (d.α > 1. && d.β == -1.) || d.α == 2.
    return d.ρ * gamma(1. + x) / gamma(1. + x / d.α)
  end
  return d.ρ * (sin(pi * d.ρ * x) * gamma(1. + x)) /
  (d.α * d.ρ * sin(pi * x / d.α) * gamma(1. + x / d.α))
end

function mellin(d::Stable,X::AbstractArray{Real})
  if (d.α > 1. && d.β == -1.) || d.α == 2.
    return d.ρ .* gamma.(1. .+ X) ./ gamma.(1. .+ X ./ d.α)
  end

  res = Float64[]
  for x = X
    push!(res,
      (real(x) >= d.α && (d.α <= 1. || (d.α < 2. && d.β != -1.))) || real(x) <= -1. ?
      Inf : (sin(pi * d.ρ * x) * gamma(1. + x)) / (d.α * d.ρ * sin(pi * x / d.α) * gamma(1. + x / d.α))
    )
  end
  return reshape(d.ρ .* res,size(X))
end

export Stable, rand, minimum, maximum, insupport, pdf, cdf, cf, mgf, mean, mellin, params

ext: jl


text:
using Signals
using StaticArrays
using BenchmarkTools
Signals.async_mode(false)
# write your own tests here
# @test 1

typ = Matrix
typ = SMatrix{4,4,Float64,16}

A = Signal(typ(rand(4,4)))
B = Signal(typ(rand(4,4)))

C = Signal(typ(rand(4,4)))
D = Signal(typ(rand(4,4)))
E = Signal(A,B) do a,b
    a*b
end
F = Signal(C,D) do a,b
    a*b
end
G = Signal(E,F) do e,f
    e*f
end
Z = typ(zeros(4,4))
A(Z)

@benchmark $A($Z)

@benchmark begin begin $A[] = $A[];end; $G(); end

Signals.debug_mode(true)

A("ff")

addprocs(1)
@everywhere using Signals
Signals.async_mode(true)

A = Signal(1)
B = remote_signal(x->x+1,A)

A(10)
A("ff")

@noinline f(a,b) = a+b
a5() = try
    f("ff",1)
catch
    st = catch_stacktrace()
    println(st)
end
b = Task(a5)
A = schedule(b)

using Signals
Signals.async_mode(true)
A_strict = Signal(1; strict_push = true)
B_strict = Signal(A_strict;state = 0) do a,state
    state.x += a
    a + 1
end

A_soft = Signal(1)
B_soft = Signal(A_soft;state = 0) do a,state
    state.x += a
    a + 1
end

A_soft(0);sleep(1) #restart the eventlopp

B_strict.state.x == 1
A_strict(10);A_strict(100);
yield();
B_strict.state.x == 111

A = 5
f() = begin
    global A
    while true
        if A == 6
            break;
        end
        sleep(0.001)
    end
end

B = @async f()
A = 6


text:
using IterativeSolvers
using Base.Test

function randSPD(T, n)
    A = rand(T, n, n) + n * I
    A' * A
end

function approx_eigenvalue_bounds(A)
    λs = eigvals(A)
    mnv, mxv = minimum(real(λs)), maximum(real(λs))
    Δ = (mxv - mnv) / 100
    mnv - Δ, mxv + Δ
end

#Chebyshev
@testset "Chebyshev" begin

n = 10
srand(1234321)

@testset "Matrix{$T}" for T in (Float32, Float64, Complex64, Complex128)
    A = randSPD(T, n)
    b = rand(T, n)
    tol = √(eps(real(T)))

    # Without a preconditioner
    begin
        λ_min, λ_max = approx_eigenvalue_bounds(A)
        x0 = rand(n)
        x, history = chebyshev(A, b, λ_min, λ_max, tol=tol, maxiter=10n, log=true)
        @test isa(history, ConvergenceHistory)
        @test history.isconverged
        @test norm(A * x - b) / norm(b) ≤ tol
    end

    # With an initial guess
    begin
        λ_min, λ_max = approx_eigenvalue_bounds(A)
        x0 = rand(T, n)
        x, history = chebyshev!(x0, A, b, λ_min, λ_max, tol=tol, maxiter=10n, log=true)
        @test isa(history, ConvergenceHistory)
        @test history.isconverged
        @test x == x0
        @test norm(A * x - b) / norm(b) ≤ tol
    end

    # With a preconditioner
    begin
        B = randSPD(T, n)
        B_fact = cholfact!(B)
        λ_min, λ_max = approx_eigenvalue_bounds(B_fact \ A)
        x, history = chebyshev(A, b, λ_min, λ_max, Pl = B_fact, tol=tol, maxiter=10n, log=true)
        @test history.isconverged
        @test norm(A * x - b) / norm(b) ≤ tol
    end
end
end


text:
using IterableTables
using NamedTuples
using Temporal
using DataFrames
using NullableArrays
import TypedTables
using Base.Test

@testset "Temporal" begin

dates  = collect(Date(1999,1,1):Date(1999,1,3))

source_ta1 = TS(collect(1:length(dates)), dates, [:value])

@test isiterable(source_ta1) == true

source_it = getiterator(source_ta1)

@test length(source_it) == 3

df1 = DataFrame(source_ta1)
@test size(df1) == (3,2)
@test df1[:Index] == dates
@test df1[:value] == [1,2,3]

source_ta2 = TS(collect(1:length(dates)), dates, :a)
df2 = DataFrame(source_ta2)
@test size(df2) == (3,2)
@test df2[:Index] == dates
@test df2[:a] == [1,2,3]

source_ta3 = TS(hcat(collect(1:length(dates)), collect(length(dates):-1:1)), dates, [:a, :b])
df3 = DataFrame(source_ta3)
@test size(df3) == (3,3)
@test df3[:Index] == dates
@test df3[:a] == [1,2,3]
@test df3[:b] == [3,2,1]

source_ta4 = TS(hcat(collect(1:length(dates)), collect(length(dates):-1:1)), dates, [:a, :b])
df4 = DataFrame(source_ta4)
@test size(df4) == (3,3)
@test df4[:Index] == dates
@test df4[:a] == [1,2,3]
@test df4[:b] == [3,2,1]

source_tt = TypedTables.@Table(a=[4.,5.], time=[Date(1999,1,1),Date(1999,1,2)], b=[6.,8.], c=[12.,24.])
ta1 = TS(source_tt, index_column=:time)
@test size(ta1) == (2,3)
@test ta1.values == [4. 6. 12.;5. 8. 24]
@test ta1.index == [Date(1999,1,1),Date(1999,1,2)]

source_df = DataFrame(a=[4.,5.], Index=[Date(1999,1,1),Date(1999,1,2)], b=[6.,8.])
ta2 = TS(source_df)
@test size(ta2) == (2,2)
@test ta2.values == [4. 6.;5. 8.]
@test ta2.index == [Date(1999,1,1),Date(1999,1,2)]

end


text:
function HDF5.h5open(filename::AbstractString, mode::AbstractString="r", pv...; swmr=false, fclose_degree = H5F_CLOSE_STRONG, auto_fclose = true)
    checkprops(pv...)
    # pv is interpreted as pairs of arguments
    # the first of a pair is a key of hdf5_prop_get_set
    # the second of a pair is a property value
    fapl = p_create(H5P_FILE_ACCESS, pv...) # file access property list
    # With garbage collection, the other modes don't make sense
    # (Set this first, so that the user-passed properties can overwrite this.)
    fapl["fclose_degree"] = fclose_degree
    fcpl = p_create(H5P_FILE_CREATE, pv...) # file create property list
    modes =
        mode == "r"  ? (true,  false, false, false, false) :
        mode == "r+" ? (true,  true,  false, false, true ) :
        mode == "cw" ? (false, true,  true,  false, true ) :
        mode == "w"  ? (false, true,  true,  true,  false) :
        # mode == "w+" ? (true,  true,  true,  true,  false) :
        # mode == "a"  ? (true,  true,  true,  true,  true ) :
        error("invalid open mode: ", mode)
    if auto_fclose
        try
            disable_dag()
            fid = h5open(filename, modes..., fcpl, fapl; swmr=swmr)
            enable_dag()
            fid
        catch e
            fapl["fclose_degree"] = H5F_CLOSE_DEFAULT
            h5open(filename, modes..., fcpl, fapl; swmr=swmr)
        end
    else
        h5open(filename, modes..., fcpl, fapl; swmr=swmr)
    end
end

function HDF5.h5open(f::Function, args...; ka...)
    fid = h5open(args...; ka...)
    try
        f(fid)
    finally
        close(fid)
    end
end

function HDF5.close(obj::HDF5File)
    if obj.id != -1
        flush(obj)
        h5f_close(obj.id)
        obj.id = -1
    end
    nothing
end

function h5load(src, ::Type{T}, pv...; mode = "r", mmaparrays = true, gc = false, ka...) where T
    gc && @eval GC.gc(true)
    fid = h5open(src, mode, pv...; ka...)
    os = Any[]
    for s in fieldnames(T)
        ft = fieldtype(T, s)
        x = if ft <: AbstractArray
            if string(s) ∈ names(fid)
                !Sys.iswindows() && mmaparrays ? tryreadmmap(fid[string(s)]) : read(fid[string(s)])
            else
                zeros(ft.parameters[1], ntuple(i -> 0, ft.parameters[2]))
            end
        else
            x = read_nonarray(fid, string(s))
            x = try ft(x) catch e x end
        end
        push!(os, x)
    end
    obj = T(os...)
    gc && Threads.nthreads() == 1 && finalizer(x -> HDF5.h5_garbage_collect(), obj)
    return obj
end

function h5load(src, paths = nothing, pv...; mode = "r", mmaparrays = true, gc = false, ka...)
    gc && @eval GC.gc(true)
    fid = h5open(src, mode, pv...; ka...)
    if isnothing(paths)
        obj = !Sys.iswindows() && mmaparrays ? tryreadmmap(fid) : read(fid)
    elseif paths isa AbstractArray
        obj = Dict(path => tryreadmmap(fid[path]) for path in paths)
    end
    gc && Threads.nthreads() == 1 && finalizer(x -> HDF5.h5_garbage_collect(), obj)
    return obj
end

h5load(src, path::AbstractString, pv...; ka...) = h5load(src, [path], pv...; ka...)[path]

function h5save(dst, obj, pv...; exclude = [], ka...)
    @eval GC.gc(true)
    h5open(dst, "w", pv...; ka...) do fid
        if obj isa AbstractDict
            for (k, v) in obj
                k ∈ exclude && continue
                write_batch(fid, k, v)
            end
        else
            for s in propertynames(obj)
                s ∈ exclude && continue
                x = getproperty(obj, s)
                if isa(x, AbstractArray)
                    write_batch(fid, string(s), x)
                elseif isa(x, HDF5Dataset)
                    HDF5.create_external(fid, s, filename(x), name(x))
                else
                    write_nonarray(fid, string(s), x)
                end
            end
        end
    end
    return dst
end

h5loadv(a...; ka...) = h5load(a...; fclose_degree = H5F_CLOSE_DEFAULT, ka...)


text:
export judiModeling, judiDataModeling, judiDataSourceModeling, judiPointSourceModeling, judiJacobian

const adjoint_map = Dict(:forward => :adjoint, :adjoint => :forward, :born => :adjoint_born, :adjoint_born => :born)

# Base abstract types
abstract type judiPropagator{D, O} <: joAbstractLinearOperator{D, D} end
abstract type judiComposedPropagator{D, O} <: judiPropagator{D, O} end
abstract type judiAbstractJacobian{D, O, FT} <: judiComposedPropagator{D, O} end

############################################################################################################################
# Concrete types
struct judiModeling{D, O} <: judiPropagator{D, O}
    m::AbstractSize
    n::AbstractSize
    model::Model
    options::JUDIOptions
end

struct judiPointSourceModeling{D, O} <: judiComposedPropagator{D, O}
    m::AbstractSize
    n::AbstractSize
    F::judiModeling{D, O}
    qInjection::AdjointProjection{D}
    function judiPointSourceModeling{D, O}(F::judiModeling{D, O}, qInjection::AdjointProjection{D}) where {D, O}
        # instantiate un-initialized sizes
        ts = time_space_src(get_nsrc(qInjection), get_nt(qInjection), F.model.n)
        merge!(F.n, ts)
        merge!(F.m, ts)
        update_size(qInjection, F)
        new(F.m, qInjection.n, F, qInjection)
    end
end

struct judiDataModeling{D, O} <: judiComposedPropagator{D, O}
    m::AbstractSize
    n::AbstractSize
    rInterpolation::Projection{D}
    F::judiModeling{D, O}
    function judiDataModeling{D, O}(rInterpolation::Projection{D}, F::judiModeling{D, O}) where {D, O}
        ts = time_space_src(get_nsrc(rInterpolation), get_nt(rInterpolation), F.model.n)
        merge!(F.n, ts)
        merge!(F.m, ts)
        update_size(rInterpolation, F)
        new(rInterpolation.m, F.n, rInterpolation, F)
    end
end

struct judiDataSourceModeling{D, O} <: judiComposedPropagator{D, O}
    m::AbstractSize
    n::AbstractSize
    rInterpolation::Projection{D}
    F::judiModeling{D, O}
    qInjection::AdjointProjection{D}
    function judiDataSourceModeling{D, O}(rInterpolation::Projection{D}, F::judiModeling{D, O}, qInjection::AdjointProjection{D}) where {D, O}
        ts = time_space_src(get_nsrc(rInterpolation), get_nt(rInterpolation), F.model.n)
        merge!(F.n, ts)
        merge!(F.m, ts)
        update_size(rInterpolation, F)
        update_size(qInjection, F)
        new(rInterpolation.m, qInjection.n, rInterpolation, F, qInjection)
    end
end

struct judiJacobian{D, O, FT} <: judiAbstractJacobian{D, O, FT}
    m::AbstractSize
    n::AbstractSize
    F::FT
    q::judiMultiSourceVector
end

struct LazyScal
    s::Number
    P::judiPropagator
end

############################################################################################################################
# Base and JOLI compat
_get_property(J::judiPropagator, ::Val{:fop}) = x -> J*x
_get_property(J::judiPropagator, ::Val{:fop_T}) = x -> J'*x
_get_property(J::judiPropagator, ::Val{:name}) = "$(typeof(J))"
_get_property(J::judiPropagator, ::Val{s}) where {s} = getfield(J, s)

_get_property(C::judiComposedPropagator, ::Val{:model}) = C.F.model
_get_property(C::judiComposedPropagator, ::Val{:options}) = C.F.options

getproperty(J::judiPropagator, s::Symbol) = _get_property(J, Val{s}())

display(P::judiPropagator{D, O}) where {D, O} = println("JUDI $(String(O)){$D} propagator $(repr(P.n)) -> $(repr(P.m))")
show(io::Union{IOBuffer, IOContext}, P::judiPropagator{D, O}) where {D, O} = print(io, "JUDI $(String(O)){$D} propagator $(repr(P.n)) -> $(repr(P.m))")

############################################################################################################################
# Constructors
"""
    judiModeling(model; options=Options())
    judiModeling(model, src_geometry, rec_geometry; options=Options())

Create seismic modeling operator for a velocity model given as a `Model` structure.
The function also takes the source and receiver geometries
as additional input arguments, which creates a combined operator `judiProjection*judiModeling*judiProjection'`.

Example
=======
`Pr` and `Ps` are projection operatos of type `judiProjection` and
`q` is a data vector of type `judiVector`:
    F = judiModeling(model)
    dobs = Pr*F*Ps'*q
    F = judiModeling(model, q.geometry, rec_geometry)
    dobs = F*q
"""
function judiModeling(model::Model; options=Options())
    D = eltype(model.m)
    m = time_space(model.n)
    return judiModeling{D, :forward}(m, m, model, options)
end

judiModeling(model::Model, src_geom::Geometry, rec_geom::Geometry; options=Options()) =
    judiProjection(rec_geom) * judiModeling(model; options=options) * adjoint(judiProjection(src_geom))

"""
    judiJacobian(F,q)

Create a linearized modeling operator from the non-linear modeling operator `F` and
the source `q`. `q` can be either a judiVector (point source Jacobian) or a `judiWeight` for
extended source modeling.
`F` is a full modeling operator including source/receiver projections.

Examples
========
1) `F` is a modeling operator without source/receiver projections:
    J = judiJacobian(Pr*F*Ps',q)
2) `F` is the combined operator `Pr*F*Ps'`:
    J = judiJacobian(F,q)
"""
function judiJacobian(F::judiComposedPropagator{D, O}, q::judiMultiSourceVector; options=nothing) where {D, O}
    update!(F.F.options, options)
    return judiJacobian{D, :born, typeof(F)}(F.m, space(F.model.n), F, q)
end

# Backward compat with giving weights as array. Not recommened
function judiJacobian(F::judiComposedPropagator{D, O}, q::Array{D, N}; options=nothing) where {D, O, N}
    @warn "judiWeights is recommended for judiJacobian(F, weights)"
    nsrc = try length(F.qInjection.data) catch; length(F.rInterpolation.data) end
    update!(F.F.options, options)
    return judiJacobian(F, judiWeights(reshape(q, F.model.n); nsrc=nsrc))
end

############################################################################################################################
# Linear algebra (conj/adj/...)
# Adjoints
conj(F::judiPropagator) = F
transpose(F::judiPropagator) = adjoint(F)

adjoint(s::Symbol) = adjoint_map[s]
adjoint(F::judiModeling{D, O}) where {D, O} = judiModeling{D, adjoint(O)}(F.n, F.m, F.model, F.options)
adjoint(F::judiDataModeling{D, O}) where {D, O} = judiPointSourceModeling{D, adjoint(O)}(adjoint(F.F), adjoint(F.rInterpolation))
adjoint(F::judiPointSourceModeling{D, O}) where {D, O}= judiDataModeling{D, adjoint(O)}(adjoint(F.qInjection), adjoint(F.F))
adjoint(F::judiDataSourceModeling{D, O}) where {D, O} = judiDataSourceModeling{D, adjoint(O)}(adjoint(F.qInjection), adjoint(F.F), adjoint(F.rInterpolation))
adjoint(J::judiJacobian{D, O, FT}) where {D, O, FT} = judiJacobian{D, adjoint(O), FT}(J.n, J.m, J.F, J.q)
adjoint(L::LazyScal) = LazyScal(L.s, adjoint(L.P))

# Composition
*(F::judiModeling{D, O}, P::AdjointProjection{D}) where {D, O} = judiPointSourceModeling{D, O}(F, P)
*(P::Projection{D}, F::judiModeling{D, O}) where {D, O} = judiDataModeling{D, O}(P, F)
*(P::Projection{D}, F::judiPointSourceModeling{D, O}) where {D, O} = judiDataSourceModeling{D, O}(P, F.F, F.qInjection)
*(F::judiDataModeling{D, O}, P::AdjointProjection{D}) where {D, O} = judiDataSourceModeling{D, O}(F.rInterpolation, F.F, P)
*(s::Number, P::judiPropagator) = LazyScal(s, P)
\(P::judiPropagator, s::Number) = LazyScal(1/s, P)
*(L::LazyScal, x::SourceType{T}) where {T<:Number} = L.s * (L.P * x)
*(L::LazyScal, x::Array{T, N}) where {T<:Number, N} = L.s * (L.P * x)

# Propagation via linear algebra `*`
*(F::judiPropagator{T, O}, q::SourceType{T}) where {T<:Number, O} = multi_src_propagate(F, q)
*(F::judiAbstractJacobian{T, O, FT}, q::dmType{T}) where {T<:Number, O, FT} = multi_src_propagate(F, q)

mul!(out::SourceType{T}, F::judiPropagator{T, O}, q::SourceType{T}) where {T<:Number, O} = begin y = F*q; copyto!(out, y) end
mul!(out::SourceType{T}, F::judiPropagator{T, O}, q::Vector{T}) where {T<:Number, O} = begin y = F*q; copyto!(out, y) end
mul!(out::SourceType{T}, F::joLinearFunction{T, T}, q::SourceType{T}) where {T<:Number} = begin y = F*q; copyto!(out, y) end
mul!(out::SourceType{T}, F::judiAbstractJacobian{T, :born, FT}, q::Vector{T}) where {T<:Number, FT} = begin y = F*q[:]; copyto!(out, y) end
mul!(out::SourceType{T}, F::judiAbstractJacobian{T, :born, FT}, q::Array{T, N}) where {T<:Number, FT, N} = begin y = F*q[:]; copyto!(out, y) end
mul!(out::Array{T, 2}, F::judiAbstractJacobian{T, :adjoint_born, FT}, q::SourceType{T}) where {T<:Number, FT} =  begin y = F*q; copyto!(out, y) end
mul!(out::Array{T, 3}, F::judiAbstractJacobian{T, :adjoint_born, FT}, q::SourceType{T}) where {T<:Number, FT} =  begin y = F*q; copyto!(out, y) end

############################################################################################################################
# Propagation input
process_input_data(::judiPropagator, data::judiMultiSourceVector) = data
process_input_data(F::judiModeling, q::Vector) = process_input_data(q, F.model)
process_input_data(F::judiPointSourceModeling, q::Vector) = process_input_data(q, F.qInjection.geometry)
process_input_data(F::judiDataSourceModeling, q::Vector) = process_input_data(q, F.qInjection.data)
process_input_data(F::judiDataModeling, q::Vector) = process_input_data(q, F.model)
process_input_data(J::judiJacobian{D, :adjoint_born, FT}, q::Vector) where {D, FT} =
    process_input_data(q, J.F.qInjection.data)
process_input_data(::judiJacobian{D, :born, FT}, q::dmType{D}) where {D<:Number, FT} = q

make_input(::judiModeling, q::SourceType) = (nothing, make_input(q), nothing, nothing, nothing)
make_input(::judiModeling, rhs::judiRHS) = (make_src(rhs)..., nothing, nothing, nothing)
make_input(F::judiPointSourceModeling, q::SourceType{T}) where {T} = (make_src(q, F.qInjection)..., nothing, nothing, nothing)
make_input(F::judiPointSourceModeling, q::Matrix{T}) where {T} = (F.qInjection.data[1], q, nothing, nothing, nothing)
make_input(F::judiDataModeling, q::SourceType{T}) where {T} = (nothing, make_input(q), F.rInterpolation.geometry, nothing, nothing)
make_input(F::judiDataModeling{T, O}, q::LazyAdd{T}) where {T, O} = (make_src(q)..., F.rInterpolation.geometry, nothing, nothing)
make_input(F::judiDataModeling, rhs::judiRHS) = (make_src(rhs)..., F.rInterpolation.geometry, nothing, nothing)

make_input(F::judiDataSourceModeling, q::SourceType{T}) where {T} = (make_src(q, F.qInjection)..., F.rInterpolation.data[1], nothing, nothing)
make_input(F::judiDataSourceModeling, q::Matrix{T}) where {T} = (F.qInjection.data[1], q, F.rInterpolation.data[1], nothing, nothing)

function make_input(J::judiJacobian{D, :born, FT}, q::dmType{D}) where {D<:Number, FT}
    srcGeom, srcData = make_src(J.q, J.F.qInjection)
    return srcGeom, srcData, J.F.rInterpolation.data[1], nothing, reshape(q, J.model.n)
end

function make_input(J::judiJacobian{D, :adjoint_born, FT}, q::SourceType{D}) where {D, FT}
    srcGeom, srcData = make_src(J.q, J.F.qInjection)
    recGeom, recData = make_src(q, J.F.rInterpolation)
    return srcGeom, srcData, recGeom, recData, nothing
end

############################################################################################################################
# Size update based on linear operator

update_size(w::judiProjection, F::judiPropagator) = set_space_size!(w.n, F.model.n)
update_size(w::jAdjoint{<:judiProjection}, F::judiPropagator) = set_space_size!(w.op.n, F.model.n)

update_size(w::judiWavelet, F::judiPropagator) = set_space_size!(w.m, F.model.n)
update_size(w::jAdjoint{<:judiWavelet}, F::judiPropagator) = set_space_size!(w.op.m, F.model.n)
############################################################################################################################
# indexing
getindex(F::judiModeling{D, O}, i) where {D, O} = judiModeling{D, O}(F.m[i], F.n[i], F.model, F.options[i])
getindex(F::judiDataModeling{D, O}, i) where {D, O} = judiDataModeling{D, O}(F.rInterpolation[i], F.F[i])
getindex(F::judiPointSourceModeling{D, O}, i) where {D, O}= judiPointSourceModeling{D, O}(F.F[i], F.qInjection[i])
getindex(F::judiDataSourceModeling{D, O}, i) where {D, O} = judiDataSourceModeling{D, O}(F.rInterpolation[i], F.F[i], F.qInjection[i])
getindex(J::judiJacobian{D, O, FT}, i) where {D, O, FT} = judiJacobian{D, O, FT}(J.m[i], J.n[i], J.F[i], J.q[i])

############################################################################################################################
# Comparisons
isequal(P1::judiPropagator, P2::judiPropagator) = P1 == P2
==(F1::judiModeling{D, O1}, F2::judiModeling{D, O2}) where {D, O1, O2} = (O1 == O2 && F1.model == F2.model && F1.options == F2.options)
==(F1::judiDataModeling, F2::judiDataModeling) = (F1.F == F2.F && F1.rInterpolation == F2.rInterpolation)
==(F1::judiPointSourceModeling, F2::judiPointSourceModeling) = (F1.F == F2.F && F1.qInjection == F2.qInjection)
==(F1::judiDataSourceModeling, F2::judiDataSourceModeling) = (F1.F == F2.F && F1.qInjection == F2.qInjection && F1.rInterpolation == F2.rInterpolation)
==(F1::judiJacobian{D, O1, FT1}, F2::judiJacobian{D, O2, FT2}) where {D, O1, O2, FT1, FT2} = (O1 == O2 && FT1 == FT2 && F1.F == F2.F && F1.q == F2.q)


text:
using OptimalControl
using Test

@testset "OptimalControl.jl" begin
    # Write your tests here.
end


text:
using Printf
using LinearAlgebra
using LinearAlgebra: Givens, Rotation

# Auxiliary
function adiagmax(A::StridedMatrix)
    adm = zero(typeof(real(A[1])))
    @inbounds begin
        for i = size(A,1)
            adm = max(adm, abs(A[i,i]))
        end
    end
    return adm
end

# Hessenberg Matrix
struct HessenbergMatrix{T,S<:StridedMatrix} <: AbstractMatrix{T}
    data::S
end

Base.copy(H::HessenbergMatrix{T,S}) where {T,S} = HessenbergMatrix{T,S}(copy(H.data))

Base.getindex(H::HessenbergMatrix{T,S}, i::Integer, j::Integer) where {T,S} = i > j + 1 ? zero(T) : H.data[i,j]

Base.size(H::HessenbergMatrix) = size(H.data)
Base.size(H::HessenbergMatrix, i::Integer) = size(H.data, i)

function LinearAlgebra.ldiv!(H::HessenbergMatrix, B::AbstractVecOrMat)
    n = size(H, 1)
    Hd = H.data
    for i = 1:n-1
        G, _ = givens!(Hd, i, i+1, i)
        lmul!(G, view(Hd, 1:n, i+1:n))
        lmul!(G, B)
    end
    ldiv!(Triangular(Hd, :U), B)
end

# Hessenberg factorization
struct HessenbergFactorization{T, S<:StridedMatrix,U} <: Factorization{T}
    data::S
    τ::Vector{U}
end

function _hessenberg!(A::StridedMatrix{T}) where T
    n = LinearAlgebra.checksquare(A)
    τ = Vector{Householder{T}}(undef, n - 1)
    for i = 1:n - 1
        xi = view(A, i + 1:n, i)
        t  = LinearAlgebra.reflector!(xi)
        H  = Householder{T,typeof(xi)}(view(xi, 2:n - i), t)
        τ[i] = H
        lmul!(H', view(A, i + 1:n, i + 1:n))
        rmul!(view(A, :, i + 1:n), H)
    end
    return HessenbergFactorization{T, typeof(A), eltype(τ)}(A, τ)
end
LinearAlgebra.hessenberg!(A::StridedMatrix) = _hessenberg!(A)

Base.size(H::HessenbergFactorization, args...) = size(H.data, args...)

# Schur
struct Schur{T,S<:StridedMatrix} <: Factorization{T}
    data::S
    R::Rotation
end

function Base.getproperty(F::Schur, s::Symbol)
    if s === :T || s === :Schur
        return triu(F.data, -1)
    else
        return getfield(F, s)
    end
end

function wilkinson(Hmm, t, d)
    λ1 = (t + sqrt(t*t - 4d))/2
    λ2 = (t - sqrt(t*t - 4d))/2
    return ifelse(abs(Hmm - λ1) < abs(Hmm - λ2), λ1, λ2)
end

# We currently absorb extra unsupported keywords in kwargs. These could e.g. be scale and permute. Do we want to check that these are false?
function _schur!(H::HessenbergFactorization{T}; tol = eps(real(T)), debug = false, shiftmethod = :Francis, maxiter = 100*size(H, 1), kwargs...) where T
    n = size(H, 1)
    istart = 1
    iend = n
    HH = H.data
    τ = Rotation(Givens{T}[])

    # iteration count
    i = 0

    @inbounds while true
        i += 1
        if i > maxiter
            throw(ArgumentError("iteration limit $maxiter reached"))
        end

        # Determine if the matrix splits. Find lowest positioned subdiagonal "zero"
        for _istart in iend - 1:-1:1
            if abs(HH[_istart + 1, _istart]) <= tol*(abs(HH[_istart, _istart]) + abs(HH[_istart + 1, _istart + 1]))
                    istart = _istart + 1
                debug && @printf("Split! Subdiagonal element is: %10.3e and istart now %6d\n", HH[istart, istart - 1], istart)
                break
            elseif _istart > 1 && abs(HH[_istart, _istart - 1]) <= tol*(abs(HH[_istart - 1, _istart - 1]) + abs(HH[_istart, _istart]))
                    debug && @printf("Split! Next subdiagonal element is: %10.3e and istart now %6d\n", HH[_istart, _istart - 1], _istart)
                istart = _istart
                break
            end
            istart = 1
        end

        # if block size is one we deflate
        if istart >= iend
            debug && @printf("Bottom deflation! Block size is one. New iend is %6d\n", iend - 1)
            iend -= 1

        # and the same for a 2x2 block
        elseif istart + 1 == iend
            debug && @printf("Bottom deflation! Block size is two. New iend is %6d\n", iend - 2)
            iend -= 2

        # run a QR iteration
        # shift method is specified with shiftmethod kw argument
        else
            Hmm = HH[iend, iend]
            Hm1m1 = HH[iend - 1, iend - 1]
            d = Hm1m1*Hmm - HH[iend, iend - 1]*HH[iend - 1, iend]
            t = Hm1m1 + Hmm
            t = iszero(t) ? eps(real(one(t))) : t # introduce a small pertubation for zero shifts
            debug && @printf("block start is: %6d, block end is: %6d, d: %10.3e, t: %10.3e\n", istart, iend, d, t)

            if shiftmethod == :Francis
                # Run a bulge chase
                if iszero(i % 10)
                    # Vary the shift strategy to avoid dead locks
                    # We use a Wilkinson-like shift as suggested in "Sandia technical report 96-0913J: How the QR algorithm fails to converge and how fix it".

                    debug && @printf("Wilkinson-like shift! Subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", HH[iend, iend - 1], HH[iend - 1, iend - 2])
                    _d = t*t - 4d
                    if _d isa Real && _d >= 0
                        # real eigenvalues
                        a = t/2
                        b = sqrt(_d)/2
                        s = a > Hmm ? a - b : a + b
                    else
                        # complex case
                        s = t/2
                    end
                    singleShiftQR!(HH, τ, s, istart, iend)
                else
                    # most of the time use Francis double shifts

                    debug && @printf("Francis double shift! Subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", HH[iend, iend - 1], HH[iend - 1, iend - 2])
                    doubleShiftQR!(HH, τ, t, d, istart, iend)
                end
            elseif shiftmethod == :Rayleigh
                debug && @printf("Single shift with Rayleigh shift! Subdiagonal is: %10.3e\n", HH[iend, iend - 1])

                # Run a bulge chase
                singleShiftQR!(HH, τ, Hmm, istart, iend)
            else
                throw(ArgumentError("only support supported shift methods are :Francis (default) and :Rayleigh. You supplied $shiftmethod"))
            end
        end
        if iend <= 2 break end
    end

    return Schur{T,typeof(HH)}(HH, τ)
end
_schur!(A::StridedMatrix; kwargs...) = _schur!(_hessenberg!(A); kwargs...)
LinearAlgebra.schur!(A::StridedMatrix; kwargs...) = _schur!(A; kwargs...)

function singleShiftQR!(HH::StridedMatrix, τ::Rotation, shift::Number, istart::Integer, iend::Integer)
    m = size(HH, 1)
    H11 = HH[istart, istart]
    H21 = HH[istart + 1, istart]
    if m > istart + 1
        Htmp = HH[istart + 2, istart]
        HH[istart + 2, istart] = 0
    end
    G, _ = givens(H11 - shift, H21, istart, istart + 1)
    lmul!(G, view(HH, :, istart:m))
    rmul!(view(HH, 1:min(istart + 2, iend), :), G')
    lmul!(G, τ)
    for i = istart:iend - 2
        G, _ = givens(HH[i + 1, i], HH[i + 2, i], i + 1, i + 2)
        lmul!(G, view(HH, :, i:m))
        HH[i + 2, i] = Htmp
        if i < iend - 2
            Htmp = HH[i + 3, i + 1]
            HH[i + 3, i + 1] = 0
        end
        rmul!(view(HH, 1:min(i + 3, iend), :), G')
        # mul!(G, τ)
    end
    return HH
end

function doubleShiftQR!(HH::StridedMatrix, τ::Rotation, shiftTrace::Number, shiftDeterminant::Number, istart::Integer, iend::Integer)
    m = size(HH, 1)
    H11 = HH[istart, istart]
    H21 = HH[istart + 1, istart]
    Htmp11 = HH[istart + 2, istart]
    HH[istart + 2, istart] = 0
    if istart + 3 <= m
        Htmp21 = HH[istart + 3, istart]
        HH[istart + 3, istart] = 0
        Htmp22 = HH[istart + 3, istart + 1]
        HH[istart + 3, istart + 1] = 0
    else
        # values doen't matter in this case but variables should be initialized
        Htmp21 = Htmp22 = Htmp11
    end
    G1, r = givens(H11*H11 + HH[istart, istart + 1]*H21 - shiftTrace*H11 + shiftDeterminant, H21*(H11 + HH[istart + 1, istart + 1] - shiftTrace), istart, istart + 1)
    G2, _ = givens(r, H21*HH[istart + 2, istart + 1], istart, istart + 2)
    vHH = view(HH, :, istart:m)
    lmul!(G1, vHH)
    lmul!(G2, vHH)
    vHH = view(HH, 1:min(istart + 3, m), :)
    rmul!(vHH, G1')
    rmul!(vHH, G2')
    lmul!(G1, τ)
    lmul!(G2, τ)
    for i = istart:iend - 2
        for j = 1:2
            if i + j + 1 > iend break end
            # G, _ = givens(H.H,i+1,i+j+1,i)
            G, _ = givens(HH[i + 1, i], HH[i + j + 1, i], i + 1, i + j + 1)
            lmul!(G, view(HH, :, i:m))
            HH[i + j + 1, i] = Htmp11
            Htmp11 = Htmp21
            # if i + j + 2 <= iend
                # Htmp21 = HH[i + j + 2, i + 1]
                # HH[i + j + 2, i + 1] = 0
            # end
            if i + 4 <= iend
                Htmp22 = HH[i + 4, i + j]
                HH[i + 4, i + j] = 0
            end
            rmul!(view(HH, 1:min(i + j + 2, iend), :), G')
            # mul!(G, τ)
        end
    end
    return HH
end

_eigvals!(A::StridedMatrix; kwargs...)           = _eigvals!(_schur!(A; kwargs...))
_eigvals!(H::HessenbergMatrix; kwargs...)        = _eigvals!(_schur!(H; kwargs...))
_eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(_schur!(H; kwargs...))

# Overload methods from LinearAlgebra to make them work generically
if VERSION > v"1.2.0-DEV.0"
    function LinearAlgebra.eigvals!(A::StridedMatrix;
                                    sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
                                    kwargs...)
            if ishermitian(A)
                return LinearAlgebra.sorteig!(eigvals!(Hermitian(A)), sortby)
            end
            LinearAlgebra.sorteig!(_eigvals!(A; kwargs...), sortby)
    end

    LinearAlgebra.eigvals!(H::HessenbergMatrix;
                        sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
                        kwargs...) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby)

    LinearAlgebra.eigvals!(H::HessenbergFactorization;
                        sortby::Union{Function,Nothing}=LinearAlgebra.eigsortby,
                        kwargs...) = LinearAlgebra.sorteig!(_eigvals!(H; kwargs...), sortby)
else

    function LinearAlgebra.eigvals!(A::StridedMatrix; kwargs...)
        if ishermitian(A)
            return eigvals!(Hermitian(A))
        end
        _eigvals!(A; kwargs...)
    end

    LinearAlgebra.eigvals!(H::HessenbergMatrix; kwargs...) = _eigvals!(H; kwargs...)

    LinearAlgebra.eigvals!(H::HessenbergFactorization; kwargs...) = _eigvals!(H; kwargs...)
end

function _eigvals!(S::Schur{T}; tol = eps(real(T))) where T
    HH = S.data
    n = size(HH, 1)
    vals = Vector{complex(T)}(undef, n)
    i = 1
    while i < n
        Hii = HH[i, i]
        Hi1i1 = HH[i + 1, i + 1]
        if abs(HH[i + 1, i]) <= tol*(abs(Hi1i1) + abs(Hii))
            vals[i] = Hii
            i += 1
        else
            d = Hii*Hi1i1 - HH[i, i + 1]*HH[i + 1, i]
            t = Hii + Hi1i1
            x = 0.5*t
            y = sqrt(complex(x*x - d))
            vals[i] = x + y
            vals[i + 1] = x - y
            i += 2
        end
    end
    if i == n
        vals[i] = HH[n, n]
    end
    return vals
end


text:
# identity matrix function
identm(x) = [row == col ? true : false for row in 1:x, col in 1:x]

max_stars_repo_name: jawuku/julia

max_issues_repo_name: jawuku/julia

max_forks_repo_name: jawuku/julia


text:
function UW(M,ndim,nbins)
    results = zeros(size(M))

    for i in 1:ndim

        if length(unique(M[i,:]))<3
            res_temp = zeros(length(M[i,:]))

            for u in 1:length(unique(M[i,:]))
                res_temp[M[i,:] .== unique(M[i,:])[u]] .= u
            end

            results[i,:] = res_temp
        else

            disc = LinearDiscretizer(binedges(DiscretizeUniformWidth(nbins), M[i,:]))
            results[i,:] = [encode(disc,one) for one in M[i,:]]

        end

    end

    return(DataFrame(results))

    return(results)

end


text:
# =========================================================================== #
# General
# =========================================================================== #

function ramp(x_start, x_delta, x_end)
    return collect( x_start:x_delta:x_end )
end

function linspace(x_start, x_end, x_length::Int)
    if x_length <= 1
        return [x_end]
    end
    x_delta = (x_end - x_start) / (x_length - 1)
    x_ = collect( x_start:x_delta:x_end )
    return x_
end

function logspace(a, b, n::Int = 50)
	retval = 10 .^ collect( LinRange(a, b, n) )
end

function mesh_grid(x, y)
    X = [i for i in x, j in 1:length(y)]
    Y = [j for i in 1:length(x), j in y]
    return (X, Y)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
"""

Binning of a vector (of draws).

$(SIGNATURES)

## Positional arguments
```julia
* `x::AbstractVectir` # DataFrame holding the draws
* `n::Int` # Number of bins
* `rng::AbstractRNG=Random.default_rng()` # Random number generator
```

From [Bogumił Kamiński](https://bkamins.github.io/julialang/2020/12/11/binning.html)

Not exported.

"""
function bin_vector(x::AbstractVector, n::Int,
    rng::AbstractRNG=Random.default_rng())

    n > 0 || throw(ArgumentError("number of bins must be positive"))
    l = length(x)

    # find bin sizes
    d, r = divrem(l, n)
    lens = fill(d, n)
    lens[1:r] .+= 1
    # randomly decide which bins should be larger
    shuffle!(rng, lens)

    # ensure that we have data sorted by x, but ties are ordered randomly
    df = DataFrame(id=axes(x, 1), x=x, r=rand(rng, l))
    sort!(df, [:x, :r])

    # assign bin ids to rows
    binids = reduce(vcat, [fill(i, v) for (i, v) in enumerate(lens)])
    df.binids = binids

    # recover original row order
    sort!(df, :id)
    return df.binids
end


text:
#
# inline-issue-minimal.jl -
#
# Demonstrate a reduction of performances with Julia-1.6.0 compared to 1.5.4.
# This version is similar to `inline-issue.jl` but does not use
# `MayOptimize.jl`.
#
module InlineIssue

using BenchmarkTools

# Catmull-Rom spline:
function spline(x::T) where {T<:AbstractFloat}
    abs_x = abs(x)
    if abs_x ≥ 2
        zero(T)
    elseif abs_x ≤ 1 # (((3/2)*x - (5/2))*x^2 + 1)
        ((T(3)/2)*abs_x - (T(5)/2))*abs_x^2 + 1
    else # (((5/2) - (1/2)*x)*x - 4)*x + 2 = (2 - x)^2 (1 - x)/2
        (2 - abs_x)^2*(1 - abs_x)/2
    end
end

# Compute all weights (10 operations).
function compute_weights(::typeof(spline),
                         t::T) where {T<:AbstractFloat}
    u = 1 - t
    v = (T(-1)/2)*t*u
    w1 = v*u
    w4 = v*t
    w = w4 - w1
    w2 = u - w1 + w
    w3 = t - w4 - w
    return (w1, w2, w3, w4)
end

function inbounds_map!(f::Function,
                       dst::Array{T,N},
                       src::Array{T,N}) where {T<:AbstractFloat,N}
    @inbounds for i in eachindex(dst, src)
        dst[i] = f(src[i])
    end
    return dst
end

function simd_map!(f::Function,
                   dst::Array{T,N},
                   src::Array{T,N}) where {T<:AbstractFloat,N}
    @inbounds @simd for i in eachindex(dst, src)
        dst[i] = f(src[i])
    end
    return dst
end

function compute_weights!(f::Function,
                          dst::Array{T,2},
                          src::Array{T,1}) where {T<:AbstractFloat}
    @assert size(dst) == (4,length(src))
    @inbounds @simd for i in eachindex(src)
        w1, w2, w3, w4 = compute_weights(f, src[i])
        dst[1,i] = w1
        dst[2,i] = w2
        dst[3,i] = w3
        dst[4,i] = w4
    end
    return dst
end

@inline function inlined_compute_weights!(f::Function,
                                          dst::Array{T,2},
                                          src::Array{T,1}) where {T<:AbstractFloat}
    @assert size(dst) == (4,length(src))
    @inbounds @simd for i in eachindex(src)
        w1, w2, w3, w4 = compute_weights(f, src[i])
        dst[1,i] = w1
        dst[2,i] = w2
        dst[3,i] = w3
        dst[4,i] = w4
    end
    return dst
end

function runtests(; T::Type{<:AbstractFloat}=Float32, n::Int=1_000)
    x = rand(T, n)
    y = Array{T}(undef, n)
    t = rand(T, n) .- one(T)/2
    z = Array{T}(undef, 4, n)
    print("Tests with Julia-", VERSION, ", T=", T, ", n=", n, "\n")
    print(" ├─ Call spline ($n times):\n");
    print(" │   ├─ map! ───────────────────────"); @btime $(map!)($spline, $y, $x)
    print(" │   ├─ inbounds_map! ──────────────"); @btime $(inbounds_map!)($spline, $y, $x)
    print(" │   └─ simd_map! ──────────────────"); @btime $(simd_map!)($spline, $y, $x)
    print(" └─ Computation of weights with spline ($n times):\n");
    print("     ├─ compute_weights! ───────────"); @btime $(compute_weights!)($spline, $z, $t)
    print("     └─ inlined_compute_weights! ───"); @btime $(inlined_compute_weights!)($spline, $z, $t)

end

runtests(T=Float32)
#println()
#runtests(T=Float64)

end # module


text:
using BarycentricInterpolation
using LinearAlgebra
using Test

# These tests arbitrarily use 5 and 6 degree polynomials (odd and even are used to
# ensure coverage)

equi₅ = Equispaced{5}()
equi₆ = Equispaced{6}()
cheb1₅ = Chebyshev1{5}()
cheb1₆ = Chebyshev1{6}()
cheb2₅ = Chebyshev2{5}()
cheb2₆ = Chebyshev2{6}()
leg₅ = Legendre{5}()
leg₆ = Legendre{6}()

x5 = [0.814723686393179,1.7205156234688,1.8475024397623,2.76087829590132,3.39323754212673,3.49077794712614]
x6 = [0.278498218867048,0.825379738072032,1.78288657350633,2.74777510870561,2.90538819038315,3.87598097214377,4.83314792038672]
arb₅ = ArbitraryPolynomial(x5)
arb₆ = ArbitraryPolynomial(x6)

@testset "Degree" begin
    @test degree(equi₅) == 5
    @test degree(equi₆) == 6
end

@testset "Nodes" begin
    @test nodes(equi₅) ≈ [-1,-0.6,-0.2,0.2,0.6,1]
    @test nodes(equi₆) ≈ [-1,-0.666666666666667,-0.333333333333333,0,0.333333333333333,0.666666666666667,1]

    @test nodes(cheb1₅) ≈ -[0.965925826289068,0.707106781186548,0.258819045102521,-0.258819045102521,-0.707106781186547,-0.965925826289068]
    @test nodes(cheb1₆) ≈ -[0.974927912181824,0.78183148246803,0.433883739117558,6.12323399573677e-17,-0.433883739117558,-0.781831482468029,-0.974927912181824]

    @test nodes(cheb2₅) ≈ -[1,0.809016994374947,0.309016994374947,-0.309016994374947,-0.809016994374947,-1]
    @test nodes(cheb2₆) ≈ -[1,0.866025403784439,0.5,6.12323399573677e-17,-0.5,-0.866025403784439,-1]

    @test nodes(leg₅) ≈ [-0.932469514203152,-0.661209386466264,-0.238619186083197,0.238619186083197,0.661209386466264,0.932469514203152]
    @test nodes(leg₆) ≈ [-0.949107912342758,-0.741531185599394,-0.405845151377397,0,0.405845151377397,0.741531185599394,0.949107912342758]

    @test nodes(arb₅) ≈ x5
    @test nodes(arb₆) ≈ x6

    @test nodes(Equispaced{5, Float64}) ≈ range(-1, stop=1, length=6)
    @test collect(nodes(Equispaced{5, Float32})) isa Vector{Float32}
end

@testset "Weights" begin
    @test weights(equi₅) ≈ [-1,5,-10,10,-5,1]
    @test weights(equi₆) ≈ [1,-6,15,-20,15,-6,1]

    @test weights(cheb1₅) ≈ -[0.258819045102521,-0.707106781186547,0.965925826289068,-0.965925826289068,0.707106781186548,-0.258819045102521]
    @test weights(cheb1₆) ≈ -[0.222520933956314,-0.623489801858733,0.900968867902419,-1,0.900968867902419,-0.623489801858734,0.222520933956314]

    @test weights(cheb2₅) ≈ -[0.5,-1,1,-1,1,-0.5]
    @test weights(cheb2₆) ≈ -[0.5,-1,1,-1,1,-1,0.5]

    leg₅w_arb = weights(ArbitraryPolynomial(nodes(leg₅)))
    leg₅w = weights(leg₅)
    leg₅w .*= leg₅w_arb[3] / leg₅w[3]
    @test leg₅w_arb ≈ leg₅w
    leg₆w_arb = weights(ArbitraryPolynomial(nodes(leg₆)))
    leg₆w = weights(leg₆)
    leg₆w .*= leg₆w_arb[3] / leg₆w[3]
    @test leg₆w_arb ≈ leg₆w

    @test weights(arb₅) ≈ [-0.079601734468871,2.82206348049722,-3.28654127480991,1.17155138701322,-2.4317744194857,1.80430256125404]
    @test weights(arb₆) ≈ [0.0114361913589559,-0.0390633084690722,0.100394094241955,-0.588770483980979,0.552871013562898,-0.0415349717143145,0.00466746500055649]

    @test collect(weights(Chebyshev2{5, Float32})) isa Vector{Float32}
end

xnew = [-0.514624351277159,0.285656117611641,0.427542456238857,0.849303738865132] # randomly generated

@testset "Interpolation matrix" begin
    @test interpolation_matrix(cheb1₅, xnew) ≈ [-0.0289845478236897 0.0959627734159381 -0.207065941756692 0.626076869404405 0.609098232472457 -0.0950873857124182;-0.0105623302174427 0.0465782216908936 0.999201970730136 -0.0492504664603377 0.0197735254445248 -0.00574092118777367;-0.0706583814898062 0.371759729030326 0.841448464571141 -0.206847346420822 0.0915972568975815 -0.0272997225884207;0.362858671520288 0.813049119299893 -0.267459002051222 0.142520710886204 -0.0742818875343108 0.0233123878791484][:, end:-1:1]
    @test interpolation_matrix(cheb1₆, xnew) ≈ [0.0127857533810382 -0.0411607067288632 0.0812980146502686 -0.166310956920905 0.955056229650331 0.199706640020739 -0.0413749740526084;0.041385687119839 -0.161088285741265 0.779202698741526 0.448772538954965 -0.160518366877509 0.0748749109083769 -0.0226291831059332;0.00285521790063521 -0.0123604217994083 0.997916710684666 0.0164279128299457 -0.00734604100263827 0.00362101552128418 -0.00111439413448407;0.184921983997812 0.964706061561875 -0.226419209876787 0.122921499942711 -0.0733011112746609 0.0399052720352278 -0.0127344963861772][:, end:-1:1]
    @test interpolation_matrix(cheb1₅, nodes(cheb1₅)) ≈ I
    @test interpolation_matrix(cheb1₆, nodes(cheb1₆)) ≈ I
end

@testset "Interpolate" begin
    c1 = Chebyshev1{10}()
    @test interpolate(c1, sin.(nodes(c1)), xnew) ≈ sin.(xnew)
    c2 = Chebyshev2{20}(0, 10)
    x = 10*rand(5)
    @test interpolate(c2, sin.(nodes(c2)), x) ≈ sin.(x)
end

@testset "Differentiation matrix" begin
    @test differentiation_matrix(cheb1₅) ≈ -[7.20976852010751 -10.5558337350587 5.27791686752937 -3.04720672422855 1.63299316185545 -0.517638090205042;1.4142135623731 0.707106781186547 -3.04720672422855 1.4142135623731 -0.707106781186548 0.218779599482357;-0.378937381963012 1.63299316185545 0.138700708242031 -1.93185165257814 0.757874763926024 -0.218779599482357;0.218779599482357 -0.757874763926024 1.93185165257814 -0.138700708242031 -1.63299316185545 0.378937381963012;-0.218779599482357 0.707106781186547 -1.4142135623731 3.04720672422855 -0.707106781186545 -1.4142135623731;0.517638090205041 -1.63299316185545 3.04720672422854 -5.27791686752936 10.5558337350587 -7.2097685201075]
    @test differentiation_matrix(cheb1₆) ≈ -[9.84466088119818 -14.5105621059791 7.48352452737972 -4.60952974192497 2.87399478545474 -1.59494677776481 0.512858431636276;1.84827792218219 1.0055981139743 -4.15304279314467 2.05143372654511 -1.18863516903897 0.639524003844967 -0.203155804362921;-0.4564869487803 1.98886610088761 0.267253750914115 -2.55809601537987 1.15238243548124 -0.569229914492255 0.17531059136945;0.22824347439015 -0.797473388882404 2.07652139657234 -1.02695629777827e-15 -2.07652139657234 0.797473388882405 -0.22824347439015;-0.175310591369451 0.569229914492254 -1.15238243548124 2.55809601537987 -0.26725375091411 -1.98886610088761 0.456486948780299;0.203155804362921 -0.639524003844966 1.18863516903897 -2.05143372654511 4.15304279314467 -1.00559811397431 -1.84827792218218;-0.512858431636278 1.59494677776481 -2.87399478545474 4.60952974192498 -7.48352452737972 14.5105621059791 -9.84466088119817]
    c1 = Chebyshev2{20}(0, 10)
    x = 10*rand(5)
    D = differentiation_matrix(c1)
    @test interpolate(c1, D*sin.(nodes(c1)), x) ≈ cos.(x)
end


text:
using JuliaFormatter

kwargs = (;
    indent = 4,
    margin = 80,
    always_for_in = true,
    whitespace_typedefs = true,
    whitespace_ops_in_indices = true,
    remove_extra_newlines = true,
    import_to_using = false,
    pipe_to_function_call = false,
    always_use_return = false, # false plays nicely with KernelAbstractions
    whitespace_in_kwargs = true,
    annotate_untyped_fields_with_any = false,
    short_to_long_function_def = false, # false allows for one-liners
    format_docstrings = false,
    align_assignment = true,
    align_struct_field = true,
    align_conditional = true,
    align_pair_arrow = true,
    align_matrix = true,
    normalize_line_endings = "auto",
    trailing_comma = true,
    join_lines_based_on_source = false, #false forces either single or multiline
    indent_submodule = true,
)

format("."; kwargs...)

max_stars_repo_name: sandreza/GradientLearning

max_issues_repo_name: sandreza/GradientLearning

max_forks_repo_name: sandreza/GradientLearning


text:
############################################################
# judiWeights ##############################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: June 2019

# Updated by Ziyi Yin (ziyi.yin@gatech.edu), Nov 2020

export judiWeights

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

# structure for seismic data as an abstract vector
mutable struct judiWeights{T<:Number} <: judiMultiSourceVector{T}
    nsrc::Integer
    data::Vector{Array{T, N}} where N
end

# Bypass mismatch in naming and fields for backward compat
Base.getproperty(obj::judiWeights, sym::Symbol) = sym == :weights ? getfield(obj, :data) : getfield(obj, sym)

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

## outer constructors

"""
    judiWeights
        nsrc
        weights
Abstract vector for weighting an extended source, which is injected at every grid point, as weighted by this vector.
Constructors
============
Construct vector cell array of weights. The `weights` keyword\\
can also be a single (non-cell) array, in which case the weights are the same for all source positions:
    judiWeights(weights; nsrc=1)
"""
function judiWeights(weights::Array{T, N}; nsrc=1) where {T<:AbstractFloat, N}
    weights = convert(Array{Float32}, weights)
    # length of vector
    weightsCell = [deepcopy(weights) for j=1:nsrc]
    return judiWeights{Float32}(nsrc, weightsCell)
end

# constructor if weights are passed as a cell array
judiWeights(weights::Vector{Array}) =
    judiWeights(convert.(Array{Float32, length(size(weights[1]))}, weights))

function judiWeights(weights::Vector{Array{T, N}}) where {T<:Number, N}
    nsrc = length(weights)
    weights = convert.(Array{Float32, N}, weights)
    return judiWeights{Float32}(nsrc, weights)
end

############################################################
# JOLI conversion
jo_convert(::Type{T}, jw::judiWeights{T}, ::Bool) where {T<:AbstractFloat} = jw
jo_convert(::Type{T}, jw::judiWeights{vT}, B::Bool) where {T<:AbstractFloat, vT} = judiWavefield{T}(jv.nsrc, jo_convert.(T, jw.weights, B))
zero(::Type{T}, v::judiWeights{vT}; nsrc::Integer=v.nsrc) where {T, vT} = judiWeights{T}(nsrc, T(0) .* v.data[1:nsrc])
(w::judiWeights)(x::Vector{<:Array}) = judiWeights(x)

function copy!(jv::judiWeights, jv2::judiWeights)
    jv.data .= jv2.data
    jv
end
copyto!(jv::judiWeights, jv2::judiWeights) = copy!(jv, jv2)

function push!(a::judiWeights{T}, b::judiWeights{T}) where T
	append!(a.weights, b.weights)
	a.nsrc += b.nsrc
end

make_input(w::judiWeights) = w.data[1]

# getindex weights container
"""
    getindex(x,source_numbers)
getindex seismic weights vectors or matrix-free linear operators and extract the entries that correspond\\
to the shot positions defined by `source_numbers`. Works for inputs of type `judiWeights`, `judiModeling`, \\
`judiProjection`, `judiJacobian`, `Geometry`, `judiRHS`, `judiPDE`, `judiPDEfull`.
Examples
========
(1) Extract 2 shots from `judiWeights` vector:
    dsub = getindex(dobs,[1,2])
(2) Extract geometry for shot location 100:
    geometry_sub = getindex(dobs.geometry,100)
(3) Extract Jacobian for shots 10 and 20:
    Jsub = getindex(J,[10,20])
"""
getindex(a::judiWeights{avDT}, srcnum::RangeOrVec) where avDT = judiWeights{avDT}(length(srcnum), a.data[srcnum])


text:
import Colors

function rgba(color::AbstractString, alpha::Float64 = 1.0)
    rgb = parse(Colors.RGB{Float64}, color)
    (Colors.red(rgb), Colors.green(rgb), Colors.blue(rgb), alpha)
end

function rgba(color::Colors.Colorant, alpha::Float64 = 1.0)
    convert(Colors.RGB{Float64}, color)
    (Colors.red(rgb), Colors.green(rgb), Colors.blue(rgb), alpha)
end


text:
export construct_vtx_map

"""
    construct_vtx_map(vtxs,dims)

Returns a matrix `M` such that `MV[i,j] == v`, where v is the index of the
vertex whose coordinates are (i,j)

Arguments:
* vtxs : a list of integer coordinates
* dims : the dimensions of the grid
"""
construct_vtx_map(vtxs, dims) = construct_vtx_grid(vtxs, dims).mtx

export construct_edge_cache

"""
    construct_edge_cache(vtxs,vtx_map)

Returns a `cache` such that `cache[v] = {(0,1),(0,0),...}`, the set of all
valid directions in which a robot may move from vertex `v`.

Arguments:
* vtxs : a list of integer coordinates
* vtx_map : a matrix such that `vtx_map[i,j] = v`, where `vtxs[i,j = v`
"""
function construct_edge_cache(vtxs, vtx_map)
    edge_cache = Vector{Set{Tuple{Int,Int}}}()
    for (v, pt) in enumerate(vtxs)
        edge_list = Set{Tuple{Int,Int}}()
        for d in [(0, 0), (-1, 0), (0, 1), (1, 0), (0, -1)]
            if get(vtx_map, (pt[1] + d[1], pt[2] + d[2]), 0) >= 1
                push!(edge_list, d)
            end
        end
        push!(edge_cache, edge_list)
    end
    edge_cache
end

export validate_edge_cache

"""
    validate_edge_cache(G,vtxs,cache)

Verifies that for a graph `G` embedded in space such that vertex `v` has
coordinates `vtxs[v]`, every edge `v₁ → v₂` is stored in the edge cache. In
other words, for every edge `v₁ → v₂`, the unit direction vector obtained by
`vtxs[v₂] - vtxs[v₁]` is a member of `cache[v₁]`.

Arguments:
* G : a graph on which the functions `vertices` and `outneighbors` can be
    called
* vtxs : a list of integer coordinates
* cache : an edge cache such that e.g. `cache[v] = {(0,1),(0,0),...}`, the
    set of all valid directions in which a robot may move from vertex `v`
"""
function validate_edge_cache(G, vtxs, cache)
    for v in vertices(G)
        vtx = vtxs[v]
        for v2 in outneighbors(G, v)
            vtx2 = vtxs[v2]
            d = [vtx2...] - [vtx...]
            d = tuple(d...)
            try
                @assert(d in cache[v], "$d not in $(edge_cache[v])")
            catch e
                print(e.msg)
                return false
            end
        end
    end
    return true
end

export construct_expanded_zones

"""
    construct_expanded_zones(vtxs,vtx_map,pickup_zones,dropoff_zones;
        shapes=[(1,1),(1,2),(2,1),(2,2)])

A utility for constructing drop-off/pick-up zones where large-footprint
objects need to be delivered. Returns a dictionary mapping vertices to a
dict of shape=>vtxs for expanded size delivery zones. For each starting
vertex, the expanded zone is selected as the appropriately sized region that
overlaps with the original vertex, does not overlap with any obstacles, and
has minimal overlap with all other vertices in the zone list.

Arguments:
* vtxs : a list of integer coordinates
* vtx_map : a matrix such that `vtx_map[i,j] == v`, where `vtxs[v] == (i,j)`
* zones : a list of integer vertex indices identifying the zones to be
    expanded
"""
function construct_expanded_zones(
    vtxs,
    vtx_map,
    zones;
    shapes = [(1, 1), (1, 2), (2, 1), (2, 2)],
    minimize_overlap=true # if false, just go from top left
)
    expanded_zones = Dict{Int,Dict{Tuple{Int,Int},Vector{Int}}}(
        v => Dict{Tuple{Int,Int},Vector{Int}}() for v in zones)
    # populate heatmap with zones
    grid_map = IndicatorGrid(VtxGrid(vtx_map)) #Int.(vtx_map .== 0) # == 1 for all obstacles
    heatmap = zeros(size(grid_map))
    for v in zones
        heatmap[vtxs[v]...] += 1.0
    end
    for s in shapes
        # s = (2,2)
        filtered_grid = IndicatorGrid(Int.(imfilter(grid_map, centered(ones(s))) .> 0))
        for v in zones
            vtx = vtxs[v]
            best_cost = Inf
            best_vtx = (-1, -1)
            for i = 0:s[1]-1
                for j = 0:s[2]-1
                    # start_vtx = clip([vtx...]-[i,j],[1,1],[size(grid_map)...] .- s .+ 1)
                    idx1 = clip(vtx[1] - i, 1, size(grid_map, 1) - s[1] + 1)
                    idx2 = clip(vtx[2] - j, 1, size(grid_map, 2) - s[2] + 1)
                    if filtered_grid[idx1, idx2] == 0
                        if minimize_overlap
                            cost = sum(heatmap[idx1:idx1+s[1]-1, idx2:idx2+s[2]-1])
                            if cost < best_cost
                                best_cost = cost
                                best_vtx = (idx1, idx2)
                            end
                        else
                            if best_cost == Inf
                                best_cost = 0.0
                                best_vtx = (idx1, idx2)
                            end
                        end
                    end
                end
            end
            if any(best_vtx .== -1)
                println("Zone ",v," at ",vtx," unreachable for agents of shape ",s)
            else
                idx1, idx2 = best_vtx
                vtx_list = sort([vtx_map[idx1:idx1+s[1]-1, idx2:idx2+s[2]-1]...])
                expanded_zones[v][s] = vtx_list
            end
        end
    end
    expanded_zones
end
construct_expanded_zones(
    vtxs,
    vtx_map,
    pickup_zones,
    dropoff_zones;
    kwargs...,
) = construct_expanded_zones(
    vtxs,
    vtx_map,
    vcat(pickup_zones, dropoff_zones);
    kwargs...,
)

export validate_expanded_zones

"""
    validate_expanded_zones(vtx_map,expanded_zones)

Verify that all vertices of each expanded zone do not overlap with obstacle
regions, and that each expanded zone contains the original vertex.
"""
function validate_expanded_zones(vtx_map, expanded_zones)
    for (v, shape_dict) in expanded_zones
        for (s, vtx_list) in shape_dict
            try
                @assert(v in vtx_list, "$v not in expanded zone for itself!")
            catch e
                println(e.msg)
                return false
            end
            for vp in vtx_list
                try
                    @assert(
                        vp > 0,
                        "$vp in expanded zone of size $s for $v is out-of-bounds",
                    )
                catch e
                    println(e.msg)
                    return false
                end
            end
        end
    end
    return true
end

export
    SparseDistanceMatrix,
    remap_idx,
    get_val_and_status

"""
    SparseDistanceMatrix{G,M}

Stores a graph and sparse distance matrix. When the distance matrix is queried
for the distance between `start` and `goal`, it first checks that this distance
has been computed. If not, it will first compute the distances from all nodes to
`goal` before returning the queried distance.
"""
struct SparseDistanceMatrix{G,M}
    graph::G
    mtx::M
    SparseDistanceMatrix(g::G,m::S=spzeros(Int,nv(g),nv(g))) where {G,S} = new{G,S}(g,m)
end
Base.size(m::SparseDistanceMatrix) = size(m.mtx)
Base.valtype(m::SparseDistanceMatrix) = valtype(m.mtx)

function get_val_and_status(A::SparseMatrixCSC{T,M},i0,i1) where {T,M}
    if !(1 <= i0 <= size(A, 1) && 1 <= i1 <= size(A, 2)); throw(BoundsError()); end
    r1 = Int(SparseArrays.getcolptr(A)[i1])
    r2 = Int(SparseArrays.getcolptr(A)[i1+1]-1)
    if (r1 > r2); return zero(T), true; end
    r1 = searchsortedfirst(rowvals(A), i0, r1, r2, Base.Sort.Forward)
    empty_flag = ((r1 > r2) || (rowvals(A)[r1] != i0))
    val = empty_flag ? zero(T) : nonzeros(A)[r1]
    return val, empty_flag
end
function get_val_and_status(A::Matrix,i,j)
    if !(1 <= i0 <= size(A, 1) && 1 <= i1 <= size(A, 2)); throw(BoundsError()); end
    return A[i,j], true
end

function (m::SparseDistanceMatrix)(v1::Int,v2::Int)
    val, empty_flag = get_val_and_status(m.mtx,v1,v2)
    if empty_flag
        m.mtx[:,v2] .= gdistances(m.graph,v2;sort_alg=RadixSort)
        val = m.mtx[v1,v2]
    end
    return val
end
Base.getindex(m::SparseDistanceMatrix,v1::Int,v2::Int) = m(v1,v2)
Base.get(m::SparseDistanceMatrix,idxs::Tuple{Int,Int},args...) = m(idxs...)

export recompute_cached_distances!
"""
    recompute_cached_distances!

Recompute all cached distances (important if, e.g., the graph has been modified)
"""
function recompute_cached_distances!(m::SparseDistanceMatrix)
    _, cols, _ = findnz(m.mtx)
    for v in unique(cols)
        m.mtx[:,v] .= gdistances(m.graph,v;sort_alg=RadixSort)
    end
    return m
end

export remove_edges!
function remove_edges!(m::SparseDistanceMatrix,edge_set)
    for e in edge_set
        rem_edge!(m.graph,e)
    end
    recompute_cached_distances!(m)
end

export add_edges!
function add_edges!(m::SparseDistanceMatrix,edge_set)
    for e in edge_set
        add_edge!(m.graph,e)
    end
    recompute_cached_distances!(m)
end

export
    config_index_to_tuple,
    config_tuple_to_index

"""
    config_index_to_tuple(shape::Tuple{Int,Int}, idx::Int)

Convert a config index to a tuple coordinate representation. Assumes the
following layout for config indices given shape `(m,n)`:
          1          2         ... n
         ____________________________
    1   | 1          2         ... n
    2   | n+1        n+2       ... 2n
    ... |
    m   | (m-1)*n+1  (m-1)*n+2 ... m*n
"""
function config_index_to_tuple(shape::Tuple{Int,Int}, idx::Int)
    if !(1 <= idx <= prod(shape)); throw(BoundsError()); end
    (div(idx - 1, shape[2]) + 1, mod(idx - 1, shape[2]) + 1)
end
function config_tuple_to_index(shape::Tuple{Int,Int}, config::Tuple{Int,Int})
    if !all(1 .<= config .<= shape); throw(BoundsError()); end
    shape[2]*(config[1]-1) + config[2]
end

export
    RemappedDistanceMatrix

"""
    RemappedDistanceMatrix

Stores a distance matrix that corresponds to a "shifted" grid graph with a new
configuration space layered over a base graph. The shape `m.s` of
`m::RemappedDistanceMatrix` defines the convolution kernel used to generate the
modified.
"""
struct RemappedDistanceMatrix{M}
    mtx::M
    # relevant for remapping
    base_vtx_map::VtxGrid # VtxGrid for original graph
    base_vtxs::Vector{Tuple{Int,Int}} # VtxList for original graph
    vtx_map::VtxGrid
    vtxs::Vector{Tuple{Int,Int}}
    s::Tuple{Int,Int} # Shape of the large agent
    config::Tuple{Int,Int} # position of sub-agent relative to top left corner

    RemappedDistanceMatrix(m::M,base_map,base_vtxs,vtx_map,vtxs,
        s=(1,1),config=(1,1)) where {M} = new{M}(
        m,base_map,base_vtxs,vtx_map,vtxs,s,config
        )
end
function RemappedDistanceMatrix(m::RemappedDistanceMatrix,config::Tuple{Int,Int})
    RemappedDistanceMatrix(
        m.mtx,
        m.base_vtx_map,
        m.base_vtxs,
        m.vtx_map,
        m.vtxs,
        m.s,
        config
    )
end
RemappedDistanceMatrix(m::RemappedDistanceMatrix,config::Int) = RemappedDistanceMatrix(m,config_index_to_tuple(m.s,config))
function RemappedDistanceMatrix(grid::IndicatorGrid,shape::Tuple{Int,Int},config::Tuple{Int,Int}=(1,1))
    base_vtx_map = VtxGrid(grid)
    base_vtxs = vtx_list_from_vtx_grid(base_vtx_map)
    filtered_grid = convolve_with_occupancy_kernel(grid,shape)
    vtx_map = VtxGrid(filtered_grid)
    vtxs = vtx_list_from_vtx_grid(vtx_map)
    graph = initialize_grid_graph_from_vtx_grid(vtx_map)
    sm = SparseDistanceMatrix(graph)
    m = RemappedDistanceMatrix(sm,base_vtx_map,base_vtxs,vtx_map,vtxs,shape,config)
end

"""
    remap_idx(m::SparseDistanceMatrix,v,config)

Maps an index and configuration from the base grid to the transformed grid
represented by m.graph. `config` represents the position of the query pt
relative to the coordinates that correspond to `v`.
"""
function remap_idx(m::RemappedDistanceMatrix,v::Int,config::Tuple{Int,Int}=m.config)
    (1 <= v <= length(m.base_vtxs)) ? true : return -1
    vtx = m.base_vtxs[v]
    vtx_ = vtx .+ 1 .- config
    # @show config, vtx, vtx_, v
    # v_ = m.vtx_map[vtx_[1],vtx_[2]]
    v_ = get(m.vtx_map, vtx_, -1)
    # @show v_
    return v_
end
function remap_idx(m::RemappedDistanceMatrix,v::Int,config_idx::Int)
    remap_idx(m,v,config_index_to_tuple(m.s,config_idx))
end

function (m::RemappedDistanceMatrix)(v1::Int,v2::Int,config=m.config,
    default_val=valtype(m.mtx)(0)
    )
    i = remap_idx(m,v1,config)
    j = remap_idx(m,v2,config)
    all((1,1) .<= (i,j) .<= size(m.mtx)) ? true : return default_val
    return m.mtx[i,j]
end

function remap_edges(m::RemappedDistanceMatrix,edge_set)
    remapped_edge_set = Set{Edge}()
    for e in edge_set
        ep = Edge(remap_idx(m,e.src),remap_idx(m,e.dst))
        if ep.src != -1 && ep.dst != -1
            push!(remapped_edge_set,ep)
        end
    end
    remapped_edge_set
end
function remove_edges!(m::RemappedDistanceMatrix,edge_set)
    remove_edges!(m.mtx,remap_edges(m,edge_set))
    return m
end
function add_edges!(m::RemappedDistanceMatrix,edge_set)
    add_edges!(m.mtx,remap_edges(m,edge_set))
    return m
end

export
    DistMatrixMap,
    get_distance,
    get_team_config_dist_function

"""
    DistMatrixMap

Maps team size to the effective distance (computed by Djikstra) between
leader (top left) vtxs.
A DistMatrixMap is constructed by starting with a base environment grid
graph, which is represented as a binary occupancy grid. The occupancy grid
is then convolved with kernels of various sizes (which represent
configurations of robots moving as a team). The output of each convolution
represents a new occupancy grid corresponding to the workspace of the robot
team.
It is assumed that the "lead" robot is always in the top left of the
configuration. If a team of robots wishes to query the distance to a
particular target configuration, they pass the leader's current vtx,
the leader's target vtx, and the team configuration (shape) to the
DistMatrixMap, which returns the correct distance.
"""
struct DistMatrixMap{F}
    dist_mtxs::Dict{Tuple{Int,Int},RemappedDistanceMatrix{F}}
end
DistMatrixMap{F}() where {F} = DistMatrixMap{F}(Dict{Tuple{Int,Int},RemappedDistanceMatrix{F}}())
function get_team_config_dist_function(d::DistMatrixMap,shape,config_idx)
    mat = d.dist_mtxs[shape]
    return RemappedDistanceMatrix(mat,config_idx)
end
function add_edges!(m::DistMatrixMap,edge_set)
    for mat in values(m.dist_mtxs)
        add_edges!(mat,edge_set)
    end
end
function remove_edges!(m::DistMatrixMap,edge_set)
    for mat in values(m.dist_mtxs)
        remove_edges!(mat,edge_set)
    end
end

"""
    get_distance(mtx_map::DistMatrixMap,v1::Int,v2::Int,shape::Tuple{Int,Int}=(1,1),config_idx=1)

Returns the length of the minimum distance collision-free path between
vertices `v1` and `v2` for an object of footprint `shape`.

Arguments:
* mtx_map : a `DistMatrixMap`
* v1 : an integer corresponding to the source vertex in the graph
* v2 : an integer corresponding to the destination vertex in the graph
* shape : the footprint of the object that will move between `v1` and `v2`
* config_idx : an integer that identifies the position of the reference
    point within the footprint.
"""
function get_distance(
    mtx_map::DistMatrixMap,
    v1::Int,
    v2::Int,
    shape::Tuple{Int,Int} = (1, 1),
    config_idx = 1,
)
    D = mtx_map.dist_mtxs[shape](v1, v2, config_idx)
end
function get_distance(mtx::Matrix, v1::Int, v2::Int, args...)
    return get(mtx, (v1, v2), 0)
end
function DistMatrixMap(
        base_vtx_map::M,
        base_vtxs::Vector{Tuple{Int,Int}};
        shapes = [(1, 1), (1, 2), (2, 1), (2, 2)],
        ) where {M<:Union{VtxGrid,Matrix{Int}}}
    grid_map = IndicatorGrid(base_vtx_map)
    dist_mtx_map = DistMatrixMap(
        Dict(s=>RemappedDistanceMatrix(grid_map,s) for s in shapes)
    )
end
Base.getindex(d::DistMatrixMap, v1::Int, v2::Int) = get_distance(d, v1, v2, (1, 1), 1)
function (d::DistMatrixMap)(v1::Int,v2::Int,args...)
    get_distance(d,v1,v2,args...)
end

export
    GridFactoryEnvironment,
    # get_graph,
    get_x_dim,
    get_y_dim,
    get_cell_width,
    get_transition_time,
    get_vtxs,
    get_pickup_zones,
    get_dropoff_zones,
    get_obstacles,
    get_pickup_vtxs,
    get_dropoff_vtxs,
    get_obstacle_vtxs,
    get_num_free_vtxs,
    get_free_zones

abstract type AbstractFactoryEnv <:AbstractGraph{Int} end
"""
    GridFactoryEnvironment
"""
@with_kw struct GridFactoryEnvironment{G} <: AbstractFactoryEnv
    graph::G = MetaDiGraph()
    x_dim::Int = 20
    y_dim::Int = 20
    cell_width::Float64 = 0.5
    transition_time::Float64 = 2.0
    vtxs::Vector{Tuple{Int,Int}} = Vector{Tuple{Int,Int}}()
    pickup_zones::Vector{Int}   = collect(1:length(vtxs))
    dropoff_zones::Vector{Int}  = collect(1:length(vtxs))
    free_zones::Vector{Int}     = collect(1:length(vtxs))
    obstacles::Vector{Tuple{Int,Int}} = Vector{Tuple{Int,Int}}()
    vtx_map::VtxGrid = construct_vtx_grid(vtxs, (x_dim, y_dim))
    edge_cache::Vector{Set{Tuple{Int,Int}}} = construct_edge_cache(vtxs,vtx_map)
    expanded_zones::Dict{Int,Dict{Tuple{Int,Int},Vector{Int}}} = construct_expanded_zones(
        vtxs, vtx_map, pickup_zones, dropoff_zones;
        shapes = filter(t->all(t.<=(x_dim,y_dim)),[(1,1),(1,2),(2,1),(2,2)])
        )
    dist_function::DistMatrixMap = DistMatrixMap(vtx_map, vtxs)
end
# get_graph(env::GridFactoryEnvironment) = env.graph
get_x_dim(env::GridFactoryEnvironment) = env.x_dim
get_y_dim(env::GridFactoryEnvironment) = env.y_dim
get_cell_width(env::GridFactoryEnvironment) = env.cell_width
get_transition_time(env::GridFactoryEnvironment) = env.transition_time
get_vtxs(env::GridFactoryEnvironment) = env.vtxs
get_pickup_zones(env::GridFactoryEnvironment) = env.pickup_zones
get_dropoff_zones(env::GridFactoryEnvironment) = env.dropoff_zones
get_obstacles(env::GridFactoryEnvironment) = env.obstacles
get_pickup_vtxs(env::GridFactoryEnvironment) = get_vtxs(env)[get_pickup_zones(env)]
get_dropoff_vtxs(env::GridFactoryEnvironment) =get_vtxs(env)[get_dropoff_zones(env)]
get_obstacle_vtxs(env::GridFactoryEnvironment) = get_obstacles(env)
get_dist_matrix(env::GridFactoryEnvironment) = env.dist_function
function GridFactoryEnvironment(
    env::E,
    graph::G,
) where {E<:GridFactoryEnvironment,G<:AbstractGraph}
    GridFactoryEnvironment(
        graph = graph,
        x_dim = env.x_dim,
        y_dim = env.y_dim,
        cell_width = env.cell_width,
        transition_time = env.transition_time,
        vtxs = env.vtxs,
        pickup_zones = env.pickup_zones,
        dropoff_zones = env.dropoff_zones,
        obstacles = env.obstacles,
        vtx_map = env.vtx_map,
    )
end

Base.zero(env::GridFactoryEnvironment{G}) where {G} =
    GridFactoryEnvironment(env, graph = G())
LightGraphs.edges(env::GridFactoryEnvironment) = edges(env.graph)
for op in [:edgetype,:edges,:has_edge,:has_vertex,:inneighbors,:is_directed,
    :ne,:nv,:outneighbors,:vertices]
    @eval LightGraphs.$op(env::GridFactoryEnvironment, args...) = $op(env.graph,args...)
end
function vtx_to_idx(env::GridFactoryEnvironment,vtx::Tuple{Int,Int})
    env.vtx_map[vtx[1],vtx[2]]
end
function idx_to_vtx(env::GridFactoryEnvironment,idx::Int)
    env.vtxs[idx]
end
function idx_from_offset(env::GridFactoryEnvironment,idx::Int,offset::Tuple{Int,Int}=(0,0))
    vtx_to_idx(env,idx_to_vtx(env,idx) .+ offset)
end
# get_x(env::E,v::Int) where {E<:GridFactoryEnvironment} = get_vtxs(env)[v][1]
# get_y(env::E,v::Int) where {E<:GridFactoryEnvironment} = get_vtxs(env)[v][2]
# get_θ(env::E,v::Int) where {E<:GridFactoryEnvironment} = 0.0
get_num_free_vtxs(env::GridFactoryEnvironment) = length(env.free_zones)
get_free_zones(env::GridFactoryEnvironment) = env.free_zones
get_distance(env::GridFactoryEnvironment,args...) = get_distance(env.dist_function,args...)
get_team_config_dist_function(env::GridFactoryEnvironment,args...) = get_team_config_dist_function(env.dist_function,args...)
function remove_edges!(env::GridFactoryEnvironment,edge_set)
    for e in edge_set
        rem_edge!(env.graph,e)
    end
    remove_edges!(env.dist_function,edge_set)
    return env
end
function add_edges!(env::GridFactoryEnvironment,edge_set)
    for e in edge_set
        add_edge!(env.graph,e)
    end
    add_edges!(env.dist_function,edge_set)
    return env
end

################################################################################
################################ READ AND WRITE ################################
################################################################################
export read_env

function TOML.parse(env::GridFactoryEnvironment)
    toml_dict = Dict()
    toml_dict["title"] = "GridFactoryEnvironment"
    toml_dict["x_dim"] = get_x_dim(env)
    toml_dict["y_dim"] = get_y_dim(env)
    toml_dict["cell_width"] = get_cell_width(env)
    toml_dict["transition_time"] = get_transition_time(env)
    toml_dict["vtxs"] = map(tup -> [tup[1], tup[2]], get_vtxs(env))
    toml_dict["pickup_zones"] = get_pickup_zones(env)
    toml_dict["dropoff_zones"] = get_dropoff_zones(env)
    toml_dict["obstacles"] = map(tup -> [tup[1], tup[2]], get_obstacles(env))
    return toml_dict
end
function TOML.print(io, env::GridFactoryEnvironment)
    TOML.print(io, TOML.parse(env))
end

"""
    read_env(io)

Loads saved environment from a .toml file.
"""
function read_env(io)
    toml_dict = TOML.parsefile(io)
    x_dim = toml_dict["x_dim"]
    y_dim = toml_dict["y_dim"]
    vtxs = map(arr -> (arr[1], arr[2]), toml_dict["vtxs"])
    vtx_grid = construct_vtx_grid(vtxs,(x_dim, y_dim))
    graph = initialize_grid_graph_from_vtx_grid(vtx_grid)
    env = GridFactoryEnvironment(
        graph = graph,
        x_dim = toml_dict["x_dim"],
        y_dim = toml_dict["y_dim"],
        cell_width = toml_dict["cell_width"],
        transition_time = toml_dict["transition_time"],
        vtxs = map(arr -> (arr[1], arr[2]), toml_dict["vtxs"]),
        pickup_zones = toml_dict["pickup_zones"],
        dropoff_zones = toml_dict["dropoff_zones"],
        obstacles = map(arr -> (arr[1], arr[2]), toml_dict["obstacles"]),
    )
    # GridFactoryEnvironment(env, initialize_factory_graph(env))
end

################################################################################
################################ INITIALIZATION ################################
################################################################################

export construct_regular_factory_world

"""
    construct_regular_factory_world()

Returns a `GridFactoryEnvironment` with regularly spaced obstacle regions
surrounded by alternating pick-up and drop-off locations.

Keyword Arguments:
* n_obstacles_x = 2 : number of obstacles in x direction
* n_obstacles_y = 2 : number of obstacles in y direction
* obs_width = [2;2] : obstacle width in both directions
* obs_offset = [1;1] : width of buffer region around each obstacle
* env_pad = [1;1] : env pad
* env_scale = 0.5 : determines the width of each grid cell when the
    coordinates of the environment are transformed to continuous Cartesian
    space.
* transition_time = 2.0 : determines the nominal travel time for a robot
    to move from one grid cell to an adjacent one.
"""
function construct_regular_factory_world(
        ;
        n_obstacles_x = 2,
        n_obstacles_y = 2,
        obs_width = [2; 2],
        obs_offset = [1; 1],
        env_pad = [1; 1],
        env_scale = 0.5,
        transition_time = 2.0,
        exclude_from_free = false,
    )
    PICKUP_FLAG = -1
    DROPOFF_FLAG = -2
    OBSTACLE_FLAG = 0
    NORMAL_VTX_FLAG = 1
    # generate occupancy grid representing the environment
    o = OBSTACLE_FLAG * ones(Int, obs_width[1], obs_width[2]) # obstacle region
    op = pad_matrix(o, (1, 1), NORMAL_VTX_FLAG) # padded obstacles region
    flag = PICKUP_FLAG
    for i = 1:size(o, 1)
        op[i+1, 1] = flag
        flag = flag == PICKUP_FLAG ? DROPOFF_FLAG : PICKUP_FLAG
    end
    for j = 1:size(o, 2)
        op[end, j+1] = flag
        flag = flag == PICKUP_FLAG ? DROPOFF_FLAG : PICKUP_FLAG
    end
    for i in reverse(collect(1:size(o, 1)))
        op[i+1, end] = flag
        flag = flag == PICKUP_FLAG ? DROPOFF_FLAG : PICKUP_FLAG
    end
    for j in reverse(collect(1:size(o, 2)))
        op[1, j+1] = flag
        flag = flag == PICKUP_FLAG ? DROPOFF_FLAG : PICKUP_FLAG
    end
    op = pad_matrix(op, (obs_offset[1] - 1, obs_offset[2] - 1), NORMAL_VTX_FLAG) # padded obstacles region
    A = repeat(op, n_obstacles_x, n_obstacles_y)
    Ap = pad_matrix(A, (env_pad[1], env_pad[2]), NORMAL_VTX_FLAG) # padded occupancy grid
    K = zeros(Int, size(Ap))

    pickup_zones = Vector{Int}()
    dropoff_zones = Vector{Int}()
    obstacles = Vector{Tuple{Int,Int}}()
    vtxs = Vector{Tuple{Int,Int}}()
    k = 0
    for i = 1:size(Ap, 1)
        for j = 1:size(Ap, 2)
            if Ap[i, j] == OBSTACLE_FLAG
                push!(obstacles, (i, j))
            else
                k += 1
                K[i, j] = k
                push!(vtxs, (i, j))
                if Ap[i, j] == PICKUP_FLAG
                    push!(pickup_zones, k)
                elseif Ap[i, j] == DROPOFF_FLAG
                    push!(dropoff_zones, k)
                end
            end
        end
    end
    if exclude_from_free
        free_zones = setdiff(collect(1:length(vtxs)),union(pickup_zones,dropoff_zones))
    else
        free_zones = collect(1:length(vtxs))
    end

    graph = initialize_grid_graph_from_vtx_grid(K)

    env = GridFactoryEnvironment(
        graph = graph,
        x_dim = size(Ap, 1),
        y_dim = size(Ap, 2),
        cell_width = env_scale,
        transition_time = transition_time,
        vtxs = vtxs,
        pickup_zones = pickup_zones,
        dropoff_zones = dropoff_zones,
        free_zones = free_zones,
        obstacles = obstacles,
    )

    return env
end

export construct_factory_env_from_vtx_grid

"""
    construct_factory_env_from_vtx_grid(vtx_grid;kwargs...)

Arguments:
* vtx_grid : a matrix such that `vtx_grid[i,j] > 0` represents free space,
    otherwise an obstacle.

Kewyword Arguments:
* cell_width = 0.5 : determines the width of each grid cell when the
    coordinates of the environment are transformed to continuous Cartesian
    space.
* transition_time = 2.0 : determines the nominal travel time for a robot
    to move from one grid cell to an adjacent one.
* pickup_zones = Int[] : a list of vertices that represent pick-up points
* dropoff_zones = Int[] : a list of vertices that represent drop-off points
"""
function construct_factory_env_from_vtx_grid(
        vtx_grid;
        cell_width = 1.0,
        transition_time = 1.0,
        pickup_zones = collect(1:sum(vtx_grid .> 0)),
        dropoff_zones = collect(1:sum(vtx_grid .> 0)),
    )

    graph = initialize_grid_graph_from_vtx_grid(vtx_grid)
    vtxs = Vector{Tuple{Int,Int}}()
    for i = 1:size(vtx_grid, 1)
        for j = 1:size(vtx_grid, 2)
            if vtx_grid[i, j] > 0
                push!(vtxs, (i, j))
            end
        end
    end

    env = GridFactoryEnvironment(
        graph = graph,
        x_dim = size(vtx_grid, 1),
        y_dim = size(vtx_grid, 2),
        cell_width = cell_width,
        transition_time = transition_time,
        vtxs = vtxs,
        pickup_zones = pickup_zones,
        dropoff_zones = dropoff_zones,
    )
end

export
    construct_factory_env_from_indicator_grid

"""
    construct_factory_env_from_indicator_grid

Args:
    * grid: an indicator grid whose zero entries denote free space (non-zero
    denotes obstacle)
"""
function construct_factory_env_from_indicator_grid(grid;kwargs...)
    construct_factory_env_from_vtx_grid(
        initialize_vtx_grid_from_indicator_grid(grid);kwargs...
    )
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

@testset "Lenses" begin
    test_n = 5000

    """Creates a 3D vector uniformly distributed on the sphere by rejection sampling, i.e., discarding all points with norm > 1.0"""
    function randunit()
        let v = rand(3)
            while (norm(v) > 1.0)
                v = rand(3)
            end
            return normalize(v)
        end
    end

    @testset "Refraction" begin
        Random.seed!(SEED)
        ηᵢ = 1.4
        ηₜ = 1.2
        for i in 1:test_n
            nₛ = randunit()
            rᵢ = randunit()
            rₜ = refractedray(ηᵢ, ηₜ, nₛ, rᵢ)
            if !(rₜ === nothing)
                sinθₜ = norm(cross(-nₛ, rₜ))
                sinθᵢ = norm(cross(nₛ, rᵢ))
                a = isapprox(ηₜ * sinθₜ, ηᵢ * sinθᵢ, rtol = RTOLERANCE, atol = ATOLERANCE) #verify snell's law for the the transmitted and incident ray
                b = isapprox(1.0, norm(rₜ), rtol = RTOLERANCE, atol = ATOLERANCE) #verify transmitted ray is unit
                perp = normalize(cross(nₛ, rᵢ))
                c = isapprox(0.0, rₜ ⋅ perp, rtol = RTOLERANCE, atol = ATOLERANCE) #verify incident and transmitted ray are in the same plane
                @test a && b && c
            end
        end
    end # testset refraction

    # snell
    @testset "Snell" begin
        Random.seed!(SEED)
        for i in 1:test_n
            r = randunit()
            nₛ = randunit()
            n1 = 1.0
            n2 = 1.4
            sθ1, sθ2 = snell(nₛ, r, n1, n2)
            sθ3, sθ4 = snell(nₛ, r, n2, n1)
            @test isapprox(sθ1 * n1, sθ2 * n2, rtol = RTOLERANCE, atol = ATOLERANCE) && isapprox(sθ3 * n2, sθ4 * n1, rtol = RTOLERANCE, atol = ATOLERANCE)
        end
    end # testset snell

    @testset "Reflection" begin
        Random.seed!(SEED)
        # reflection
        for i in 1:test_n
            r = randunit()
            nₛ = randunit()
            if abs(r ⋅ nₛ) < 0.9999 # if r and n are exactly aligned then their sum will be zero, not a vector pointing in the direction of the normal
                reflected = reflectedray(nₛ, r)
                #ensure reflected and refracted rays are in the same plane
                a = isapprox(norm(reflected ⋅ cross(r, nₛ)), 0.0, rtol = RTOLERANCE, atol = ATOLERANCE)
                b = isapprox(norm(reflected), 1.0, rtol = RTOLERANCE, atol = ATOLERANCE)
                c = isapprox(0.0, reflected ⋅ nₛ + r ⋅ nₛ, rtol = RTOLERANCE, atol = ATOLERANCE)
                #ensure sum of reflected and origin ray are in the direction of the normal
                d = isapprox(0.0, norm(nₛ - sign(reflected ⋅ nₛ) * (normalize(reflected - r))), rtol = RTOLERANCE, atol = ATOLERANCE)
                @test a & b & c & d
            end
        end
    end # testset reflection

    @testset "Paraxial" begin
        # check that normally incident rays are focussed across the lens
        l = ParaxialLensEllipse(100.0, 10.0, 10.0, [1.0, 1.0, 0.0], [3.0, 3.0, 3.0])
        r = OpticalRay([2.0, 2.0, 3.0], [1.0, 1.0, 0.0], 1.0, 0.55)
        intsct = halfspaceintersection(surfaceintersection(l, r))
        ref, _, _ = OpticSim.processintersection(OpticSim.interface(intsct), OpticSim.point(intsct), OpticSim.normal(intsct), r, 20.0, 1.0, true, true)
        @test isapprox(ref, normalize([1.0, 1.0, 0.0]), rtol = RTOLERANCE, atol = ATOLERANCE)
        l = ParaxialLensRect(100.0, 10.0, 10.0, [1.0, 1.0, 0.0], [3.0, 3.0, 3.0])
        r = OpticalRay([2.3, 2.4, 3.1], [1.0, 1.0, 0.0], 1.0, 0.55)
        intsct = halfspaceintersection(surfaceintersection(l, r))
        fp = [3.0, 3.0, 3.0] + 100 * normalize([1.0, 1.0, 0.0])
        ref, _, _ = OpticSim.processintersection(OpticSim.interface(intsct), OpticSim.point(intsct), OpticSim.normal(intsct), r, 20.0, 1.0, true, true)
        @test isapprox(ref, normalize(fp - OpticSim.point(intsct)), rtol = RTOLERANCE, atol = ATOLERANCE)
    end # testset Paraxial
end # testset lenses


text:
export initClipWav
export clipWav

"""
    y = clipWav(x::Array, a)
use clipping to create distortion effects
"""
function clipWav(x::Array, a)
    if !(2>=a>=0.5); a=1.0; end
    s = (0.4 + 0.1 * rand());
    c = (s + 1) / 2;
    k = (1 - s) / (s - c)^2;
    x .*= a/(maximum(abs.(x)) + 3e-5);
    L = veclen(x)
    for t = 1:L
        if c < x[t]
            x[t] = 1.0;
        elseif s < x[t] <= c
            x[t] = -k * ( x[t] - c )^2 + 1;
        elseif -c <= x[t] < -s
            x[t] =  k * ( x[t] + c )^2 - 1;
        elseif x[t] < -c
            x[t] = -1.0;
        end
    end
    return x
end

"""
    initClipWav(clipSpan::NTuple{2,Number}) -> clipwav(wav::Array)
init distortion effect function
+ `clipSpan` e.g. (0.5, 2.0)
"""
function initClipWav(clipSpan::NTuple{2,Number})
    clipMin, clipMax = clipSpan
    @assert clipMin <= clipMax
    function clipwav(wav::Array)
        return clipWav(wav, rand()*(clipMax - clipMin) + clipMin)
    end
    return clipwav
end

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
include("registration.jl")

# Classic Control
#=============================================================================#
register("CartPole-v0",
         :CartPoleEnv,
         "/classic_control/CartPole.jl",
	     max_episode_steps=200,
		 reward_threshold=195.0)

register("CartPole-v1",
         :CartPoleEnv,
         "/classic_control/CartPole.jl",
	     max_episode_steps=500,
		 reward_threshold=475.0)

register("Pendulum-v0",
         :PendulumEnv,
         "/classic_control/Pendulum.jl",
		 max_episode_steps=200)

register("MountainCarContinuous-v0",
         :Continuous_MountainCarEnv,
         "/classic_control/Continuous-MountainCar.jl",
		 max_episode_steps=999,
		 reward_threshold=90.0)


text:
using Test
using MotifSequenceGenerator

struct FloatShout
  shout::String
  dur::Float64
  start::Float64
end

using Random

let

N = 5

shouts = [FloatShout(uppercase(randstring(rand(3:5))), rand()+1, rand()) for k in 1:N]
shoutlimits(s::FloatShout) = (s.start, s.start + s.dur);
shouttranslate(s::FloatShout, n) = FloatShout(s.shout, s.dur, s.start + n);
q = 10.0

function shoutlens(R)
    l = 0
    for r in R
        a, b = shoutlimits(r)
        l += b - a
    end
    return l
end

@testset "Float Length δq=$(δq)" for δq in [1.0, 2.0]
    for j in 1:N
        r, s = random_sequence(shouts, q, shoutlimits, shouttranslate, δq)
        ℓ = shoutlens(r)
        @test q - δq ≤ ℓ ≤ q + δq
    end
end

@testset "Float Length, Weights, δq=$(δq)" for δq in [1.0, 2.0]
    weights = rand(1:5, N)
    for j in 1:N
        r, s = random_sequence(shouts, q, shoutlimits, shouttranslate, δq;
        weights = weights)
        ℓ = shoutlens(r)
        @test q - δq ≤ ℓ ≤ q + δq
    end
end

using MotifSequenceGenerator: DeadEndMotifs
@test_throws ArgumentError random_sequence(shouts, q, shoutlimits, shouttranslate, 0.0)
@test_throws DeadEndMotifs random_sequence(shouts, q, shoutlimits, shouttranslate, 0.000001)

end


text:
#          osx4 osx5 ai4  ai5
# linreg   1.03 1.51 0.94 1.08
# housing  0.01 0.01 0.01 0.01
# mnist    4.66 4.96 1.57 1.73
# lenet    -    -    6.55 7.18
# charlm   20.0 22.7 0.33 0.34

using Knet

# Pkg.test("Knet")
load_only = true

for (p,f,o1,o2,o3) =
    (
     (:LinReg, "linreg.jl", "--gcheck 2", "--fast", "--fast"),
     (:Housing, "housing.jl", "--gcheck 2 --atype Array{Float64}", "--fast", "--fast"),
     (:MNIST, "mnist.jl", "--gcheck 2", "--fast", "--fast"),
     (:LeNet, "lenet.jl", "--gcheck 2", "--fast", "--fast"),
     #(:CharLM, "charlm.jl", "--gcheck 2 --winit 0.01", "--fast", "--fast"),
     # TODO: gcheck errors for charlm
     (:CharLM, "charlm.jl", "--gcheck 2", "--fast", "--fast"),
    )
    gpu() < 0 && p == :LeNet && continue
    include(f)
    m = eval(:($p.main))
    m(o1); m(o2); m(o3)
    gc(); gpu() >= 0 && Knet.knetgc()
end


text:
using Statistics
using Logging
using Printf
using DataDeps
using JLD2

using Oceananigans
using Oceananigans.Units
using Oceananigans.Coriolis
using Oceananigans.Models.HydrostaticFreeSurfaceModels
using Oceananigans.TurbulenceClosures

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

Logging.global_logger(OceananigansLogger())

dd = DataDep("cubed_sphere_32_grid",
    "Conformal cubed sphere grid with 32×32 grid points on each face",
    "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/cubed_sphere_grids/cubed_sphere_32_grid.jld2",
    "b1dafe4f9142c59a2166458a2def743cd45b20a4ed3a1ae84ad3a530e1eff538" # sha256sum
)

DataDeps.register(dd)

## Choose a lat-lon grid or cubed sphere face grid

H = 4kilometers

# grid = LatitudeLongitudeGrid(size = (60, 60, 1), longitude = (-40, 40), latitude = (-40, 40), z = (-H, 0))

cs32_filepath = datadep"cubed_sphere_32_grid/cubed_sphere_32_grid.jld2"
grid = ConformalCubedSphereFaceGrid(cs32_filepath, face=1, Nz=1, z=(-H, 0))

## Turbulent diffusivity closure

const νh₀ = 5e3 * (60 / grid.Nx)^2

@inline νh(λ, φ, z, t) = νh₀ * cos(π * φ / 180)

variable_horizontal_diffusivity = HorizontallyCurvilinearAnisotropicDiffusivity(νh=νh)
constant_horizontal_diffusivity = HorizontallyCurvilinearAnisotropicDiffusivity(νh=νh₀)

## Model setup

model = HydrostaticFreeSurfaceModel(
          architecture = CPU(),
                  grid = grid,
    momentum_advection = VectorInvariant(),
          free_surface = ExplicitFreeSurface(gravitational_acceleration=0.1),
        # free_surface = ImplicitFreeSurface(gravitational_acceleration=0.1)
           #  coriolis = nothing,
              coriolis = HydrostaticSphericalCoriolis(scheme = VectorInvariantEnstrophyConserving()),
               closure = nothing,
             # closure = constant_horizontal_diffusivity,
             # closure = variable_horizontal_diffusivity,
               tracers = nothing,
              buoyancy = nothing
)

## Very small sea surface height perturbation so the resulting dynamics are well-described
## by a linear free surface.

A  = 1e-5 * H  # Amplitude of the perturbation
λ₀ = 0   # Central longitude
φ₀ = 20  # Central latitude
Δλ = 10  # Longitudinal width
Δφ = 10  # Latitudinal width

η′(λ, φ, z) = A * exp(- (λ - λ₀)^2 / Δλ^2) * exp(- (φ - φ₀)^2 / Δφ^2)

set!(model, η=η′)

# g = model.free_surface.gravitational_acceleration
# gravity_wave_speed = sqrt(g * H) # hydrostatic (shallow water) gravity wave speed

# # Time-scale for gravity wave propagation across the smallest grid cell
# wave_propagation_time_scale = min(grid.radius * cosd(maximum(abs, grid.ϕᵃᶜᵃ)) * deg2rad(grid.Δλ),
#                                   grid.radius * deg2rad(grid.Δϕ)) / gravity_wave_speed

mutable struct Progress
    interval_start_time :: Float64
end

function (p::Progress)(sim)
    wall_time = (time_ns() - p.interval_start_time) * 1e-9

    @info @sprintf("Time: %s, iteration: %d, max(u⃗): (%.2e, %.2e) m/s, extrema(η): (min=%.2e, max=%.2e), wall time: %s",
                   prettytime(sim.model.clock.time),
                   sim.model.clock.iteration,
                   maximum(abs, sim.model.velocities.u),
                   maximum(abs, sim.model.velocities.v),
                   minimum(abs, sim.model.free_surface.η),
                   maximum(abs, sim.model.free_surface.η),
                   prettytime(wall_time))

    p.interval_start_time = time_ns()

    return nothing
end

simulation = Simulation(model,
                        Δt = 20minutes,
                        stop_time = 30days,
                        iteration_interval = 20,
                        progress = Progress(time_ns()))

output_fields = merge(model.velocities, (η=model.free_surface.η,))

output_prefix = grid isa LatitudeLongitudeGrid ? "lat_lon_waves" : "cubed_sphere_face_waves"

simulation.output_writers[:fields] = JLD2OutputWriter(model, output_fields,
                                                      schedule = TimeInterval(1hour),
                                                      prefix = output_prefix,
                                                      force = true)

run!(simulation)

# #####
# ##### Animation!
# #####

# include("visualize_barotropic_gyre.jl")

# visualize_barotropic_gyre(simulation.output_writers[:fields])