text:
# Test CG complete
@testset "CG small problem" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = rand(T,(n,n))
        A = sqrt(A*A')
        b = rand(T,n)
        alg = CG(maxiter = 2n, tol = 10n*eps(real(T))*norm(b)) # because of loss of orthogonality, we choose maxiter = 2n
        x, info = @inferred linsolve(wrapop(A), wrapvec(b), wrapvec(zero(b)), alg)
        @test b ≈ A*unwrapvec(x)
        @test info.converged > 0

        A = rand(T,(n,n))
        A = sqrt(A*A')
        α₀ = rand(real(T)) + 1
        α₁ = rand(real(T))
        x, info = @inferred linsolve(wrapop(A), wrapvec(b), wrapvec(zero(b)), alg, α₀, α₁)
        @test b ≈ (α₀*I+α₁*A)*unwrapvec(x)
        @test info.converged > 0
    end
end

# Test CG complete
@testset "CG large problem" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = rand(T,(N,N))
        A = sqrt(sqrt(A*A'))/N
        b = rand(T,N)
        x, info = @inferred linsolve(wrapop(A), wrapvec(b);
                    isposdef = true, maxiter = 1, krylovdim = N, rtol = 10*N*eps(real(T)))
        @test b ≈ A*unwrapvec(x) + unwrapvec(info.residual)

        α₀ = rand(real(T)) + 1
        α₁ = rand(real(T))
        x, info = @inferred linsolve(wrapop(A), wrapvec(b), α₀, α₁;
                    isposdef = true, maxiter = 1, krylovdim = N, rtol = 10*N*eps(real(T)))
        @test b ≈ (α₀*I+α₁*A)*unwrapvec(x) + unwrapvec(info.residual)
    end
end

# Test GMRES complete
@testset "GMRES full factorization" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = rand(T,(n,n)).-one(T)/2
        b = rand(T,n)
        alg = GMRES(krylovdim = n, maxiter = 2, tol = 2*n*eps(real(T))*norm(b))
        x, info = @inferred linsolve(wrapop(A), wrapvec(b), wrapvec(zero(b)), alg)
        @test info.converged > 0
        @test b ≈ A*unwrapvec(x)

        A = rand(T,(n,n))
        α₀ = rand(T)
        α₁ = -rand(T)
        x, info = @inferred(linsolve(wrapop(A), wrapvec(b), wrapvec(zero(b)), alg, α₀, α₁))
        @test b ≈ (α₀*I+α₁*A)*unwrapvec(x)
        @test info.converged > 0
    end
end

# Test GMRES with restart
@testset "GMRES with restarts" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = rand(T,(N,N)).-one(T)/2
        A = I-T(9/10)*A/maximum(abs, eigvals(A))
        b = rand(T,N)
        x, info = @inferred linsolve(wrapop(A), wrapvec(b); krylovdim = 3*n, maxiter = 50, rtol = 10*N*eps(real(T)))
        @test b ≈ A*unwrapvec(x) + unwrapvec(info.residual)

        A = rand(T,(N,N)).-one(T)/2
        α₀ = maximum(abs, eigvals(A))
        α₁ = -rand(T)
        α₁ *= T(9)/T(10)/abs(α₁)
        x, info = @inferred linsolve(wrapop(A), wrapvec(b), α₀, α₁; krylovdim = 3*n, maxiter = 50, rtol = 10*N*eps(real(T)))
        @test b ≈ (α₀*I+α₁*A)*unwrapvec(x) + unwrapvec(info.residual)
    end
end

# Test BICGStab
@testset "BiCGStab" begin
    @testset for T in (Float32, Float64, ComplexF32, ComplexF64)
        A = rand(T,(N,N)).-one(T)/2
        A = I-T(9/10)*A/maximum(abs, eigvals(A))
        b = rand(T,N)
        alg = BiCGStab(; maxiter=2N, tol=10N*eps(real(T))*norm(b))
        x, info = @inferred linsolve(wrapop(A), wrapvec(b), wrapvec(zero(b)), alg)
        @test b ≈ A*unwrapvec(x)

        A = rand(T,(N,N)).-one(T)/2
        α₀ = maximum(abs, eigvals(A))
        α₁ = -rand(T)
        α₁ *= T(9)/T(10)/abs(α₁)
        x, info = @inferred linsolve(wrapop(A), wrapvec(b), wrapvec(zero(b)), alg, α₀, α₁)
        @test b ≈ (α₀*I+α₁*A)*unwrapvec(x) # + unwrapvec(info.residual)
    end
end


text:
include(joinpath(Sys.BINDIR, "..", "share", "julia", "stdlib", "v$(VERSION.major).$(VERSION.minor)", "REPL", "src", "emoji_symbols.jl"));

📁 = open("emoji_symbols.py", "w")

println(📁, "from __future__ import unicode_literals\n\n")
println(📁, "emoji_symbols = [")
for (📞, (🔑, 🍺)) in enumerate(emoji_symbols)
    print(📁, "  (\"", escape_string(🔑), "\", u\"",  🍺, "\")")
    📞 < length(emoji_symbols) && print(📁, ",")
    println(📁, "")
end
println(📁, "]")

close(📁)

ext: jl

max_stars_repo_name: albertas-jn/Julia-sublime

max_issues_repo_name: albertas-jn/Julia-sublime

max_forks_repo_name: albertas-jn/Julia-sublime


text:
export sigmoid

doc"""
    sigmoid(x)

Sigmoid logistic function.

```math
f(x) = (1 + \exp(-x))^{-1}
```
"""
sigmoid(x::Var) = Var(sigmoid(x.data), ∇sigmoid!, (x,))
sigmoid(x::Array) = sigmoid.(x)
sigmoid(x::T) where T<:AbstractFloat = T(1 / (1 + exp(-x)))
sigmoid(x::CuArray) = CUDNN.sigmoid(x)
sigmoid(x::Node) = Node(sigmoid, (x,))

function ∇sigmoid!(y::Var, x::Var)
    isnothing(x.grad) && return
    ∇sigmoid!(y.data, y.grad, x.data, x.grad)
end

function ∇sigmoid!(y::Array{T}, gy::Array{T}, x::Array{T}, gx::Array{T}) where T
    @inbounds for i = 1:length(gx)
        gx[i] += gy[i] * y[i] * (T(1) - y[i])
    end
end

∇sigmoid!(y::CuArray, gy::CuArray, x::CuArray, gx::CuArray) = CUDNN.∇sigmoid!(y, gy, x, gx)


text:
# To implement
function lenscap(lensmap, f, args...)

end


text:
#===================================================================================================
This is a container type for bringing together CPU and Memory with a common clock

There's a good chance this implementation will change, so try not to make functions written
for this thing as generic as possible!
===================================================================================================#

# this will serve as the default. Note this was the NES mater clock / 16
# this works out to about 559ns
const NES_CLOCK_PERIOD = 1.0/1.662607 * 10.0^(-6)

mutable struct Chipset
    cpu::CPU
    ram::Memory

    δt::Float64

    clock::Int64

    Chipset(c::CPU, m::Memory, δt::AbstractFloat) = new(c, m, δt, 0)
end

Chipset(c::CPU, m::Memory) = Chipset(c, m, NES_CLOCK_PERIOD)
Chipset() = Chipset(CPU(), Memory())

ν(cs::Chipset) = 1.0/cs.δt
ω(cs::Chipset) = 2π/cs.δt

# a tick where the CPU is forced to sleep
idletick!(cs::Chipset) = (cs.clock += 1)


text:
using PkgBenchmark
results = benchmarkpkg("ReverseMcCormick")
show(results)

#=
# specify tag and uncommit to benchmark versus prior tagged version
tag =
results = judge("ReverseMcCormick", tag)
show(results)
=#

export_markdown("results.md", results)


text:
__precompile__(false)

module ShaderManager

using ModernGL

using RessourceManager
using LogManager
using FileManager
using GraphicsManager
using GLLists

const GPU = GraphicsManager

mutable struct ShaderProperty
  value::Union{Nothing,Any}
  default::Union{Nothing,Any}
  typ::Union{Nothing,GLenum}
  layout::Union{Nothing,Dict{Symbol,Any}}
  uniform::Bool

  ShaderProperty(;value=nothing,default=nothing,typ=nothing,layout=nothing,uniform=false) = new(value,default,typ,layout,uniform)
end
export ShaderProperty

Uniform(;value=nothing,default=nothing,typ=nothing,layout=nothing) = ShaderProperty(value=value,default=default,typ=typ,layout=layout,uniform=true)
export Uniform

mutable struct ShaderProgram
  id::Integer
  locations::Dict{String, Integer}

  ShaderProgram(id::Integer) = new(id,Dict())
end
export ShaderProgram

abstract type IShaderManager end

mutable struct JShaderManager <: IShaderManager
  programs::Dict{Symbol, ShaderProgram}
  active_program::Union{Nothing, Symbol}

  JShaderManager() = new(Dict(),nothing)
end
export JShaderManager

""" TODO """
getActiveProgram(this::IShaderManager) = this.active_program
export getActiveProgram

function clean(this::IShaderManager)
  delPrograms(this)
end
export clean

""" TODO """
function doProgram(this::IShaderManager, f::Function)
  for (k,p) in this.programs f((k,p)) end
end
export doProgram

""" TODO """
function delPrograms(this::IShaderManager)
  for (k,p) in this.programs
    if this.active_program == k useProgram(this, nothing) end
    glDeleteProgram(p.id)
  end
  this.programs = typeof(this.programs)()
end
export delPrograms

""" TODO """
function useProgram(this::IShaderManager, program::Union{Nothing, Symbol})
  this.active_program = program
  p = program != nothing ? getProgram(this, program) : nothing
  @GLCHECK glUseProgram(program != nothing ? p.id : 0)
  p
end
export useProgram

""" TODO """
hasProgram(this::IShaderManager, program::Symbol) = haskey(this.programs, program)
export hasProgram

""" TODO """
getProgram(this::IShaderManager, program::Symbol) = hasProgram(this, program) ? this.programs[program] : nothing
export getProgram

""" TODO """
setProgram(this::IShaderManager, program::Symbol, id::Integer) = this.programs[program] = ShaderProgram(id)
export setProgram

nested_eltype(x) = nested_eltype(typeof(x))
nested_eltype(::Type{T}) where T<:AbstractArray = nested_eltype(eltype(T))
nested_eltype(::Type{T}) where T = T

#@deprecated
#flatten(a::AbstractArray) = begin try; while true a = vec(a...) end; catch ex; end; a; end

# faster version
function flatten(arr::AbstractArray)
    rst = nested_eltype(arr)[]
    grep(v) = for x in v
        if isa(x, Array) grep(x) else push!(rst, x) end
    end
    grep(arr)
    rst
end

""" TODO """
function delProgram(this::IShaderManager, program::Symbol)
  p = getProgram(this, program)
  if p == nothing return end
  if this.active_program == program useProgram(this, nothing) end
  this.programs = filter((x)->x[1] != program, this.programs)
  @GLCHECK glDeleteProgram(p.id)
  p
end
export delProgram

""" TODO """
function reloadProgram(this::IShaderManager, program::Symbol, shaders::AbstractArray)
  result=false

  #try
    id = createShaderProgram(shaders; name=program)
    if id <= 0 return false end

    if id > 0 && id != getProgram(this, program)
      delProgram(this, program)
      setProgram(this, program, id)
    end
  #catch ex
  #  error("Could not load $program: $ex.")
  #end
  result
end
export reloadProgram

""" TODO """
function setProperty(this::IShaderManager, name::String, value::Any; program::Union{Nothing,Symbol}=this.active_program)
  if program == nothing end
  p = getProgram(this, program)
  if p == nothing return end

  l = glGetUniformLocation(p.id, name)

  if l>-1
    elems=length(value)
    typ=typeof(value)
    dims=0

    if isa(value,AbstractArray)
      sz = size(value)
      typ = nested_eltype(value)
      dims = length(sz) > 1 ? sz[2] : 1
      value = flatten(value) #typ[value...]
      #dims = isa(value, AbstractArray{typ,1}) ? 1 : 2
    end

    t = nothing
    val = typ(0)

    if dims == 0
      t = (l, value)
      elems = 0
    elseif dims == 1 && isa(val, AbstractFloat)
      t = (l, 1, value)
    elseif dims == 2 && isa(val, AbstractFloat)
      t = (l, 1, false, value)
    end

    list = GLLists.UNIFORMS
    f = haskey(list, dims) && haskey(list[dims], elems) && haskey(list[dims][elems], typ) ? list[dims][elems][typ] : nothing

    if t != nothing && f != nothing
      (f)(t...)
    else
      warn("Shader Program $program: Shader Property $name (typ:$(string(typ)), dims:$dims, elems:$elems) is not defined!")
    end
  end
end
export setMode

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

function getShaderTypeString(value::Any)
  len=length(value)
	isArray = isa(value,AbstractArray)
  typ = isArray ? nested_eltype(value) : typeof(value)
	base_typ = GLLists.BASE_TYPES[typ].name
	type_str=base_typ #float
	if isArray
		if len>=2 && len <= 4
			type_str = (base_typ[1] != 'f' ? base_typ[1] : "")*"vec"*string(len)
			if isa(value[1],AbstractArray) type_str*="["*string(len)*"]" end #vec2[3]
		else type_str*="["*string(len)*"]"  #float[3]
		end
	end
	type_str
end

"""
load all content from shaders located in shaders folder in root dir
"""
function loadShaders(shader_vars=Dict{Symbol,Any}();auto_inputs=false,dir=RessourceManager.getPath(:SHADERS)) #joinpath(@__DIR__,"../shaders/"))
  types=["VSH","FSH","GSH","CSH"]
  uniforms = haskey(shader_vars, :UNIFORMS) ? shader_vars[:UNIFORMS] : Dict()
  globals = haskey(shader_vars, :GLOBALS) ? shader_vars[:GLOBALS] : Dict()

  #shader_files = filter(x->isfile(dir*x) && uppercase(replace(splitext(x)[end],"."=>"")) == "GLSL",readdir(dir))

  list, list_paths = findFiles(dir, "GLSL")

  # replace vars and set types
  for (key,file) in list
    shader=list[key]
    content=file[:content]

    keystr=string(key)
    typ = :NOTHING
    for t in types
      if occursin(t, keystr)
        typ = Symbol(t)
        break
      end
    end

    inputs = ""

    if auto_inputs
      if typ == :FSH
        inputs *="layout(location = 0) out vec4 fragColor;\n"
      else typ == :CSH
        inputs *= "layout(local_size_x = \$iDispatchX, local_size_y = \$iDispatchY, local_size_z = \$iDispatchZ) in;\n"
      end
      if length(inputs)>0 inputs*="\n" end
    end

    if length(uniforms)>0
      for (uname, shp) in uniforms
        ulayout=""
        if shp.layout != nothing
          ulayout = "layout("
          c=0
          for (k,v) in shp.layout
            c+=1
            ulayout*=(c>1 ? ", " : "")*string(k)*(v == nothing ? "" : "="*string(v))
          end
          ulayout *= ") "
        end

        utype = ""
        try
          utype = shp.typ != nothing ? GLLists.LIST_TYPE_STRING[shp.typ] : getShaderTypeString(shp.value)
        catch ex
          warn("Shader Property $uname::$(typeof(shp.value)) is not valid!")
          continue
        end

        uvalue = "" #!haskey(uform, :default) || uform.default == nothing ? "" : *" = "*string(uform.default)
        inputs*=ulayout*"uniform "*utype*" "*string(uname)*uvalue*";\n"
      end
      inputs*="\n"
    end

    vars=(x->x[1]).(collect(eachmatch(r"\$(\w+)",inputs*content)))

    for var in vars
      entry = nothing
      svar=Symbol(var)
      if haskey(globals, svar) entry=string(globals[svar])
      elseif haskey(list, var) entry="\""*list[var][:content]*"\""
      end
      if entry != nothing
        inputs=replace(inputs,"\$"*var=>entry)
        content=replace(content,"\$"*var=>entry)
      end
    end

    shader[:type] = typ
    shader[:inputs] = inputs
    shader[:content] = content
  end

  # replace imports
  found=true
  while found
    found=false
    for (key,file) in list
      content=file[:content]

      if occursin("#import ", content)
        r=collect(eachmatch(r"(\#import\s+\"(\w+(\.\w+)*)\")",content))
        paths=[[x[1],x[2]] for x in r]

        for path in paths
          f=path[2]
          if dirname(f) == "" f=dir*f end
          kkey=Symbol(abspath(f))
          line="// FILE \""*f*"\" NOT FOUND!"
          if haskey(list_paths, kkey) line=list_paths[kkey][:content] end
          content = replace(content,path[1]=>line)
        end

        list[key][:content]=content

        found=true
      end
    end
  end

  # add glsl version
  for (key,file) in list
    file[:content] = get_glsl_version_string()*"\n"*file[:inputs]*file[:content]
  end

  list
end
export loadShaders

#stat("nodes.txt").mtime

end #ShaderManager

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
function tridiag(Nvec, gamma, a, b, c, r, x)

    gamma .= 0.0

    beta = b[1]
    x[1] = r[1] / beta

    for n = 2:Nvec
        gamma[n] = c[n - 1] / beta
        beta = b[n] - a[n] * gamma[n]
        x[n] = (r[n] - a[n] * x[n - 1]) / beta
    end

    for n = (Nvec - 1):-1:1
        x[n] = x[n] - gamma[n + 1] * x[n + 1]
    end

    return nothing

end


text:
"""
The Hbv type contains the state variables (sm, suz, slz, st_uh), the inputs
(epot, prec) for one time step, the parameters (param) and the time step
length (tstep) for the HBV model.
"""
type Hbv <: Hydro

  sm::Float64
  suz::Float64
  slz::Float64
  st_uh::Array{Float64,1}
  hbv_ord::Array{Float64,1}
  epot::Float64
  prec::Float64
  q_sim::Float64
  param::Array{Float64,1}
  tstep::Float64
  time::DateTime

end

"""
    Hbv(tstep, time)

Constructor for HBV with predefined state variables, parameters and inputs.
The time step (tstep) is given as a fraction of one day. Thus, for hourly input
data tstep should be set to 1/24.
"""
function Hbv(tstep::Float64, time::DateTime)

  sm      = 0.0
  suz     = 0.0
  slz     = 0.0
  st_uh   = zeros(Float64, 20)

  epot   = 0.0
  prec   = 0.0

  q_sim = 0.0

  # fc, lp, k0, k1, k2, beta, perc, ulz, maxbas

  param = [100., 0.8, 0.05, 0.05, 0.01, 1., 2., 30., 2.5]

  hbv_ord = compute_hbv_ord(param[9])

  Hbv(sm, suz, slz, st_uh, hbv_ord, epot, prec, q_sim, param, tstep, time)

end

"""
    Hbv(tstep, time, param)

Constructor for HBV with predefined state variables and inputs. The parameter
values are given as input. The time step (tstep) is given as a fraction of one
day. Thus, for hourly input data tstep should be set to 1/24.
"""
function Hbv(tstep::Float64, time::DateTime, param::Array{Float64,1})

  sm      = 0.0
  suz     = 0.0
  slz     = 0.0
  st_uh   = zeros(Float64, 20)

  epot   = 0.0
  prec   = 0.0

  q_sim = 0.0

  hbv_ord = compute_hbv_ord(param[9])

  Hbv(sm, suz, slz, st_uh, hbv_ord, epot, prec, q_sim, param, tstep, time)

end

"""
    init_states(mdata::Hbv)

Initilize the state variables of the model.
"""
function init_states(mdata::Hbv)

  mdata.sm      = 0.0
  mdata.suz     = 0.0
  mdata.slz     = 0.0

  for i in eachindex(mdata.st_uh)
    mdata.st_uh[i] = 0.0
  end

end

"""
    get_param_range(mdata::Hbv)

Get allowed parameter ranges for the calibration of the model.
"""
function get_param_range(mdata::Hbv)

  param_range_hydro = [(1., 1000.),      # fc
                       (0.5, 0.99),      # lp
                       (0.001, 0.999),   # k0
                       (0.001, 0.999),   # k1
                       (0.001, 0.999),   # k2
                       (1., 5.),         # beta
                       (0.1, 1000.),     # perc
                       (1., 1000.),      # ulz
                       (1., 20.)]       # maxbas

end

"""
    get_states(mdata::Hbv)

Get state variables for computing penelty during calibration.

"""
function get_states(mdata::Hbv)

  return [mdata.suz; mdata.slz]

end

"""
    assign_param(mdata::Hbv, param::Array{Float64,1})

Assign parameter values to the Hbv type.
"""
function assign_param(mdata::Hbv, param::Array{Float64,1})

  for i in eachindex(mdata.param)
    mdata.param[i] = param[i]
  end

  mdata.hbv_ord = compute_hbv_ord(param[9])

end

"""
    enkf_hydro(mdata::Array{Hbv,1}, obs_ens, q_sim)

Implementation of the ensemble Kalman filter for the HBV model.
"""
function enkf_hydro(mdata::Array{Hbv,1}, obs_ens, q_sim)

  nens = length(mdata)

  # Allocate arrays

  sm    = zeros(Float64, length(mdata[1].sm), nens)
  suz   = zeros(Float64, length(mdata[1].suz), nens)
  slz   = zeros(Float64, length(mdata[1].slz), nens)
  st_uh = zeros(Float64, length(mdata[1].st_uh), nens)

  # Add states to arrays

  for iens = 1:nens

    sm[:, iens]    = mdata[iens].sm
    suz[:, iens]   = mdata[iens].suz
    slz[:, iens]   = mdata[iens].slz
    st_uh[:, iens] = mdata[iens].st_uh

  end

  # Run ensemble kalman filter

  sm    = enkf(sm, obs_ens, q_sim)
  suz   = enkf(suz, obs_ens, q_sim)
  slz   = enkf(slz, obs_ens, q_sim)
  st_uh = enkf(st_uh, obs_ens, q_sim)

  # Check limits of states

  sm[sm .< 0] = 0.
  sm[sm .> mdata[1].param[1]] = mdata[1].param[1]

  suz[suz .< 0] = 0.
  slz[slz .< 0] = 0.
  st_uh[st_uh .< 0] = 0.

  # Add arrays to states

  for iens = 1:nens

    mdata[iens].sm = sm[:, iens][1]
    mdata[iens].suz = suz[:, iens][1]
    mdata[iens].slz = slz[:, iens][1]
    mdata[iens].st_uh = st_uh[:, iens]

  end

end

"""
    run_timestep(mdata::Hbv)

Propagate the model one time step and return simulated dischage.
"""
function run_timestep(mdata::Hbv)

  # mdata and parameters

  sm = mdata.sm
  suz = mdata.suz
  slz = mdata.slz
  st_uh = mdata.st_uh
  hbv_ord = mdata.hbv_ord
  epot = mdata.epot
  prec = mdata.prec

  fc     = mdata.param[1]
  lp     = mdata.param[2]
  k0     = mdata.param[3]
  k1     = mdata.param[4]
  k2     = mdata.param[5]
  beta   = mdata.param[6]
  perc   = mdata.param[7]
  ulz    = mdata.param[8]
  maxbas = mdata.param[9]

  # Input for current time step

  prec_now = prec
  epot_now = epot

  # Soil moisture zone (assume no evaporation during rainfall)

  if prec_now > 0.

    # Beta function

    f_recharge = (sm / fc) ^ beta

    # Groundwater recharge

    recharge = f_recharge * prec_now

    # Update soil moisture zone

    sm = sm + prec_now - recharge

    # Add excess soil moisture to groundwater recharge

    if sm > fc
      recharge += sm - fc
      sm = fc
    end

    # No evapotranspiration

    eact = 0.

  else

    # Compute actual evapotranspiration

    eact = epot_now * min(sm/(fc*lp), 1.)

    # Update soil moisture zone

    sm = sm - eact

    # Check limits for soil moisture zone

    if sm < 0.
      eact = max(eact + sm, 0.)
      sm = 0.
    end

    # No groundwater recharge

    recharge = 0.

  end

  # Add recharge to upper groundwater box

  suz = suz + recharge

  # Remove percolation from upper groundwater box

  perc_now = min(perc, suz)

  suz = suz - perc_now

  # Compute runoff from upper groundwater box and update storage

  q_suz = k1 * suz + k0 * max(suz-ulz, 0.)

  suz = suz - q_suz

  if suz < 0.
    q_suz = max(q_suz + suz, 0.)
    suz = 0.
  end

  # Add precolation to lower groundwater box

  slz = slz + perc_now

  # Compute runoff from lower groundwater box and update storage

  q_slz = k2 * slz

  slz = slz - q_slz

  # Convolution of unit hydrograph

  q_tmp = q_suz + q_slz

  nh = length(hbv_ord)

  for k = 1:nh-1
    st_uh[k] = st_uh[k+1] + hbv_ord[k]*q_tmp
  end

  st_uh[nh] = hbv_ord[nh] * q_tmp

  # Compute total runoff

  q_tot = st_uh[1]

  # Assign states

  mdata.time  += Dates.Hour(mdata.tstep)
  mdata.sm = sm
  mdata.suz = suz
  mdata.slz = slz
  mdata.st_uh = st_uh
  mdata.q_sim = q_tot

  return nothing

end

# Function for computing ordinates of unit hydrograph

function compute_hbv_ord(maxbas)

  triang = Distributions.TriangularDist(0, maxbas)

  triang_cdf = Distributions.cdf(triang, 0:20)

  hbv_ord = diff(triang_cdf)

  return(hbv_ord)

end

max_stars_repo_name: jmgnve/Vann

max_forks_repo_name: jmgnve/Vann


text:
function DataFrames.DataFrame(c::PerformanceCurve, args... ; kwargs...)
    df = DataFrames.DataFrame(x=xvalues(c), y=yvalues(c), args...; kwargs...)
    DataFrames.rename!(df, :x => xlabel(c), :y => ylabel(c))
end


text:
export GroupActions, apply_actions

################
# Group actions
################
"""
    GroupActions(actions::Function...)

Store a bunch of group actions `(f1, f2, f3, ...)`.
Each action has to return a tuple.
The actions are applied in the following sense
1) f1 is applied on the original solution `s`
2) f2 is applied on `s` and the results of 1
3) f3 is applied on `s` and the results of 1) and 2)
and so on

## Example
```julia-repl
julia> f1(s) = (s * s,);

julia> f2(s) = (2s, -s, 5s);

julia> f3(s) = (s + 1,);

julia> GroupActions(f1)(3)
(3, 9)

julia> GroupActions(f1, f2)(3)
(3, 9, 6, -3, 15, 18, -9, 45)

julia> GroupActions(f1,f2, f3)(3)
(3, 9, 6, -3, 15, 18, -9, 45, 4, 10, 7, -2, 16, 19, -8, 46)
```
"""
struct GroupActions{T<:Tuple}
    actions::T
end
GroupActions(::Nothing) = GroupActions(())
GroupActions(actions::GroupActions) = actions
GroupActions(actions::Function...) = GroupActions(actions)
GroupActions(actions) = GroupActions(actions...)

function (actions::GroupActions)(s)
    S = [s]
    T = typeof(s)
    apply_actions(actions, s) do sᵢ
        sⱼ = convert(T,sᵢ)
    	if sⱼ != s
    		push!(S,sⱼ)
    	end
    	false
    end
    S
end

apply_actions(cb, action::GroupActions, s) = apply_actions(action.actions, s, cb)
function apply_actions(actions::Tuple, x, cb)
    f, rest = first(actions), Base.tail(actions)
    y = f(x)
    if isa(x, AbstractVector{<:Number}) && isa(y, AbstractVector{<:Number})
        cb(y) && return true
        if apply_actions(rest, y, cb)
            return true
        end
    else
        for yᵢ in f(x)
            cb(yᵢ) && return true
            if apply_actions(rest, yᵢ, cb)
                return true
            end
        end
    end
    apply_actions(rest, x, cb)
end
apply_actions(::Tuple{}, s, cb) = false


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module MiscTest
using ParallelAccelerator
using Compat

#ParallelAccelerator.DomainIR.set_debug_level(4)
#ParallelAccelerator.ParallelIR.set_debug_level(3)
#ParallelAccelerator.CGen.set_debug_level(4)
#ParallelAccelerator.set_debug_level(4)
#using CompilerTools
#CompilerTools.LivenessAnalysis.set_debug_level(3)
#CompilerTools.LambdaHandling.set_debug_level(3)

@acc function for_ret()
    for f = 1:10
    end
end

function test1()
    return for_ret()
end

@acc function for_ret1(n)
    for f = 1:n
    end
end

function test2()
    return for_ret1(10)
end

@acc function opt_At_mul_B!(X, W)
     X' * W
end

function test_At_mul_B(m::Int, k::Int, n::Int)
    W = Array{Float64}(m, k)
    X = Array{Float64}(m, n)
    fill!(W, 3)
    fill!(X, 5)
    opt_At_mul_B!(X, W)
end

function test3()
  all(Bool[ x == 150.0 for x in MiscTest.test_At_mul_B(10,10,10) ])
end

@acc function f()
    W = zeros(Int, 5, 5)
    s = Int[sum(W[:,j]) for j in 1:5]
end

function test4()
    f()
end

@acc function const_array_init()
    Int[1,2,3] .+ Int[4,5,6]
end

function test5()
    const_array_init() == [5,7,9]
end

@acc function mod_rem_test(x, y)
    Int[ mod(x, y), rem(x, y) ]
end

@acc function end_test(a)
    a[3:end] .+ 1
end

function test6()
    end_test([1;2;3;4]) == [4;5]
end

@acc function end_test_numeric(a)
    a[3:4] .+ 1
end

function test7()
    end_test_numeric([1;2;3;4]) == [4;5]
end

@acc function accum1()
    a = zeros(Int, 3)
    b = a
    for i = 1:3
       a += ones(Int, 3)
    end
    a .+ b
end

@acc function accum2()
    a = zeros(Int, 3)
    b = a
    for i = 1:3
       a .+= ones(Int, 3)
    end
    a .+ b
end

@acc function test8_acc1(a :: Array{Int64, 3})
    a[:, 1, 1] + 2 * 3
end
@acc function test8_acc2(a :: Array{Int64, 3})
    a[1, :, 1] + 2 * 3
end
@acc function test8_acc3(a :: Array{Int64, 3})
    a[1, 1, :] + 2 * 3
end
@acc function test8_acc4(a :: Array{Int64, 3})
    a[:, :, 1] + 2 * 3
end
@acc function test8_acc5(a :: Array{Int64, 3})
    a[:, 1, :] + 2 * 3
end
@acc function test8_acc6(a :: Array{Int64, 3})
    a[1, :, :] + 2 * 3
end

function test8()
    a = ones(Int64, 5,6,7)
    sub_tests = [test8_acc1;
                 test8_acc2;
                 test8_acc3;
                 test8_acc4;
                 test8_acc5;
                 test8_acc6]
    res = true
    for subtest in sub_tests
        b = subtest(a)
        c = @noacc subtest(a)
        this_res = (ndims(b) == ndims(c) && sum(b) == sum(c))
        if !this_res
            println("sub-test ", subtest, " failed. correct = ", c, " accelerated = ", b)
        end
        res = res && this_res
    end
    res
end

end

using Base.Test
println("Testing miscellaneous features...")
#@test MiscTest.test1() == nothing
#@test MiscTest.test2() == nothing
@test MiscTest.test3()
@test MiscTest.test4() == [0.0, 0.0, 0.0, 0.0, 0.0]
@test MiscTest.test5()
#@test MiscTest.test6()
@test MiscTest.test7()
@test MiscTest.test8()
@test MiscTest.mod_rem_test(7,3) == [1, 1]
@test MiscTest.mod_rem_test(7,-3) == [-2, 1]
@test MiscTest.mod_rem_test(-7,3) == [2, -1]
@test MiscTest.mod_rem_test(-7,-3) == [-1, -1]
@test MiscTest.accum1() == [3,3,3]
if VERSION > v"0.5.0-" && !Compat.is_windows()
@test MiscTest.accum2() == [6,6,6]
else
@test MiscTest.accum2() == [3,3,3]
end
println("Done testing miscellaneous features...")


text:
module tst_cov
    using Test
    using CloudCovErr

    @testset "PerStar" begin
        ttt_residimIn = ones(51,51)
        ttt_maskim = zeros(Bool,51,51)
        ttt_maskim[26,26] = true
        ttt_starim = ones(51,51)
        data_in, stars_in, kmasked2d = stamp_cutter(26,26,ttt_residimIn,ttt_starim,ttt_maskim)
        @test kmasked2d[17,17]

        function psfmodel_test(x,y,Np)
            psftest = [
             0.02  0.02  0.03  0.03  0.04  0.04  0.04  0.03  0.03  0.02  0.02;
             0.02  0.03  0.04  0.05  0.06  0.06  0.06  0.05  0.04  0.03  0.02;
             0.03  0.04  0.05  0.07  0.09  0.1   0.09  0.07  0.05  0.04  0.03;
             0.03  0.05  0.07  0.11  0.17  0.2   0.17  0.11  0.07  0.05  0.03;
             0.04  0.06  0.09  0.17  0.33  0.5   0.33  0.17  0.09  0.06  0.04;
             0.04  0.06  0.1   0.2   0.5   1.0   0.5   0.2   0.1   0.06  0.04;
             0.04  0.06  0.09  0.17  0.33  0.5   0.33  0.17  0.09  0.06  0.04;
             0.03  0.05  0.07  0.11  0.17  0.2   0.17  0.11  0.07  0.05  0.03;
             0.03  0.04  0.05  0.07  0.09  0.1   0.09  0.07  0.05  0.04  0.03;
             0.02  0.03  0.04  0.05  0.06  0.06  0.06  0.05  0.04  0.03  0.02;
             0.02  0.02  0.03  0.03  0.04  0.04  0.04  0.03  0.03  0.02  0.02;
            ]
            return psftest
        end
        kmasked2d = zeros(Bool,11,11)
        kmasked2d[6,end] = true
        kmasked2d[6,6] = true
        circmask = kstar_circle_mask(11;rlim=25)
        x_star = 10
        y_star = 10
        flux_star = 100
        psft, kstar, kpsf2d, kcond0, kcond, kpred, dnt = gen_pix_mask(kmasked2d,psfmodel_test,circmask,x_star,y_star,flux_star,Np=11,thr=2000)
        @test reshape(kstar,11,11)[6,end]
        @test !reshape(kstar,11,11)[6,1]
        @test (kcond0, kcond, kpred, dnt) == (71, 71, 9, 0)
        psft, kstar, kpsf2d, kcond0, kcond, kpred, dnt = gen_pix_mask(kmasked2d,psfmodel_test,circmask,x_star,y_star,flux_star,Np=11,thr=200)
        @test !reshape(kstar,11,11)[1,1]
        @test !reshape(kpsf2d,11,11)[1,1]
        function psfmodel_test(x,y,Np)
            psftest = [
             0.02  0.02  0.03  0.03  0.04  0.04  0.04  0.03  0.03  0.02  0.02;
             0.02  0.03  0.04  0.05  0.06  0.06  0.06  0.05  0.04  0.03  0.02;
             0.03  0.04  0.05  0.07  0.09  0.1   0.09  0.07  0.05  0.04  0.03;
             0.03  0.05  0.07  0.11  0.17  0.2   0.17  0.11  0.07  0.05  0.03;
             0.04  0.06  0.09  0.17  0.33  0.5   0.33  0.17  0.09  0.06  0.04;
             0.04  0.06  0.1   0.2   0.5   1.0   0.5   0.2   0.1   0.06  0.04;
             0.04  0.06  0.09  0.17  0.33  0.5   0.33  0.17  0.09  0.06  0.04;
             0.03  0.05  0.07  0.11  0.17  0.2   0.17  0.11  0.07  0.05  0.03;
             0.03  0.04  0.05  0.07  0.09  0.1   0.09  0.07  0.05  0.04  0.03;
             0.02  0.03  0.04  0.05  0.06  0.06  0.06  0.05  0.04  0.03  0.02;
             -0.2  0.02  0.03  0.03  0.04  0.04  0.04  0.03  0.03  0.02  0.02;
            ]
            return psftest
        end
        flux_star = 1e5
        psft, kstar, kpsf2d, kcond0, kcond, kpred, dnt = gen_pix_mask(kmasked2d,psfmodel_test,circmask,x_star,y_star,flux_star,Np=11,thr=200)
        @test (dnt .& 2^5) != 0

        Np = 5
        widx = 3
        widy = 3
        px0 = 1
        py0 = 1
        sx0 = 8
        sy0 = 8
        tilex = 2
        tiley = 2
        T = Float32
        Δx = (widx-1)÷2
        Δy = (widy-1)÷2
        padx = Np+Δx+px0
        pady = Np+Δy+py0
        stepx = (sx0+2) ÷ tilex
        stepy = (sy0+2) ÷ tiley
        in_subimage = zeros(T,stepx+2*padx,stepy+2*pady)
        in_subimage = convert(Array{Float32,2},reshape(1:size(in_subimage)[1]*size(in_subimage)[2],size(in_subimage)))
        ism = zeros(T,stepx+2*padx,stepy+2*pady)
        bimage = zeros(T,stepx+2*padx-2*Δx,stepy+2*pady-2*Δy)
        bism = zeros(T,stepx+2*padx-2*Δx,stepy+2*pady-2*Δy,2*Np-1, Np);
        cov1 = zeros(T,Np*Np,Np*Np)
        μ1 = zeros(T,Np*Np)
        cov_avg!(bimage, ism, bism, in_subimage, Np=Np, widx=widx, widy=widy)
        build_cov!(cov1,μ1,6,6,bimage,bism,Np,widx,widy)
        @test abs(cov1[1,1]-271.500183105)<1e-7
        @test abs(cov1[1,end]-271.499633789)<1e-7
        @test abs(cov1[end,end]-271.500732421875)<1e-7

        kpsf2d=zeros(Bool,3,3)
        kpsf2d[2,2] = true
        kpsf2d[2,3] = true
        km = kpsf2d[:]
        km[1] = true
        km[end] = true
        data_in = [
          0.516654    0.128205   1.94308;
         -0.0805772  -0.920908   0.212453;
         -0.774471    0.165229  -0.363535;
        ]
        stars_in = 2*ones(3,3)
        cov_loc = [
        47.2245   1.33946   0.628186   0.841306   0.288469  -0.437706   0.754434  -0.245601  -0.150857;
         0.0     47.2602    1.39354    0.570245   0.871201   0.243023  -0.457453   0.806491  -0.209147;
         0.0      0.0      47.166      1.40991    0.471935   0.767516   0.251835  -0.398556   0.772153;
         0.0      0.0       0.0       47.1832     1.44052    0.461676   0.757649   0.282371  -0.314066;
         0.0      0.0       0.0        0.0       47.2298     1.34296    0.471258   0.766785   0.218381;
         0.0      0.0       0.0        0.0        0.0       47.2845     1.36993    0.492573   0.715498;
         0.0      0.0       0.0        0.0        0.0        0.0       47.1459     1.3015     0.425096;
         0.0      0.0       0.0        0.0        0.0        0.0        0.0       47.1451     1.30432;
         0.0      0.0       0.0        0.0        0.0        0.0        0.0        0.0       47.2068;
        ]
        μ = [0.005037597380578518
             0.003629490500316024
             0.0026571564376354218
             0.0032599247060716152
             0.002772299339994788
             0.0034678715746849775
             0.0051109688356518745
             0.0037797277327626944
             0.003628455102443695]
        psft = [
        0.0274769  0.0360258  0.032605;
        0.0330318  0.0403823  0.0339776;
        0.0301219  0.033388   0.0255433;
        ]
        out = condCovEst_wdiag(cov_loc,μ,km,kpsf2d,data_in,stars_in,psft;Np=3,export_mean=false,n_draw=0,diag_on=true)
        @test length(out) == 1
        @test length(out[1]) == 6
        @test abs.(out[1][1] .- 133.81487945671913) .< 1e-7
        out = condCovEst_wdiag(cov_loc,μ,km,kpsf2d,data_in,stars_in,psft;Np=3,export_mean=true,n_draw=0,diag_on=true)
        @test length(out) == 2
        @test size(out[2]) == (3,3)
        out = condCovEst_wdiag(cov_loc,μ,km,kpsf2d,data_in,stars_in,psft;Np=3,export_mean=true,n_draw=2,diag_on=true)
        @test length(out) == 3
        @test size(out[3]) == (9,2)
    end
end


text:
using Test, LinearMaps, LinearAlgebra

@testset "transpose/adjoint" begin
    A = 2 * rand(ComplexF64, (20, 10)) .- 1
    v = rand(ComplexF64, 10)
    w = rand(ComplexF64, 20)
    V = rand(ComplexF64, 10, 3)
    W = rand(ComplexF64, 20, 3)
    Av = A * v
    AV = A * V
    M = @inferred LinearMap(A)
    N = @inferred LinearMap(M)

    @test @inferred M' * w == A' * w
    @test @inferred mul!(copy(V), adjoint(M), W) ≈ A' * W
    @test @inferred transpose(M) * w == transpose(A) * w
    @test @inferred transpose(LinearMap(transpose(M))) * v == A * v
    @test @inferred LinearMap(M')' * v == A * v
    @test @inferred transpose(transpose(M)) == M
    @test (M')' == M
    Mherm = @inferred LinearMap(A'A)
    @test @inferred ishermitian(Mherm)
    @test @inferred !issymmetric(Mherm)
    @test @inferred !issymmetric(transpose(Mherm))
    @test @inferred ishermitian(transpose(Mherm))
    @test @inferred ishermitian(Mherm')
    @test @inferred isposdef(Mherm)
    @test @inferred isposdef(transpose(Mherm))
    @test @inferred isposdef(adjoint(Mherm))
    @test @inferred !(transpose(M) == adjoint(M))
    @test @inferred !(adjoint(M) == transpose(M))
    @test @inferred transpose(M') * v ≈ transpose(A') * v
    @test @inferred transpose(LinearMap(M')) * v ≈ transpose(A') * v
    @test @inferred LinearMap(transpose(M))' * v ≈ transpose(A)' * v
    @test @inferred transpose(LinearMap(transpose(M))) * v ≈ Av
    @test @inferred adjoint(LinearMap(adjoint(M))) * v ≈ Av

    @test @inferred mul!(copy(w), transpose(LinearMap(M')), v) ≈ transpose(A') * v
    @test @inferred mul!(copy(w), LinearMap(transpose(M))', v) ≈ transpose(A)' * v
    @test @inferred mul!(copy(w), transpose(LinearMap(transpose(M))), v) ≈ Av
    @test @inferred mul!(copy(w), adjoint(LinearMap(adjoint(M))), v) ≈ Av
    @test @inferred mul!(copy(V), transpose(M), W) ≈ transpose(A) * W
    @test @inferred mul!(copy(V), adjoint(M), W) ≈ A' * W

    B = @inferred LinearMap(Symmetric(rand(10, 10)))
    @test transpose(B) == B
    @test B == transpose(B)

    B = @inferred LinearMap(Hermitian(rand(ComplexF64, 10, 10)))
    @test adjoint(B) == B
    @test B == B'
end


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
SUITE["Build"] = BenchmarkGroup()

# ==============================================================================
# Decomposition-based approach
# ==============================================================================
include("building_BFFPSV18.jl")

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=false)

# return the sample with the smallest time value in each test
println("minimum: ", minimum(results))

# return the median for each test
println("median: ", median(results))

# ==============================================================================
# Create plots
# ==============================================================================

function plot_dashed_line!(v, time_horizon=20.)
    plot!(x->x, x->v, 0.0, time_horizon, line=2, color="red", linestyle=:dash,
          legend=nothing)
end

𝑂_BLDF01[:mode] = "reach"
𝑂_BLDF01[:plot_vars] = [0, 25]
𝑂_BLDF01[:project_reachset] = true
res = solve(build_TV, 𝑂_BLDF01, op=BFFPSV18(𝑂_dense_BLDF01))

plot(res,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t\raisebox{-0.5mm}{\textcolor{white}{.}}",
     ylab=L"x_{25}\raisebox{2mm}{\textcolor{white}{.}}",
     xtick=[0., 0.5, 1.], ytick=[-6e-3, -3e-3, 0., 3e-3, 6e-3],
     xlims=(0., 1.), ylims=(-Inf, 6e-3),
     bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
     size=(1000, 1000))
plot_dashed_line!(4e-3)
plot_dashed_line!(5.1e-3)
savefig(@relpath "BLDF01_time_horizon_1.png")

plot(res,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t\raisebox{-0.5mm}{\textcolor{white}{.}}",
     ylab=L"x_{25}\raisebox{2mm}{\textcolor{white}{.}}",
     xtick=[0., 10., 20.], ytick=[-6e-3, -3e-3, 0., 3e-3, 6e-3],
     bottom_margin=6mm, left_margin=2mm, top_margin=3mm,
     ylims=(-Inf, 6e-3),
     size=(1000, 1000))
plot_dashed_line!(4e-3)
plot_dashed_line!(5.1e-3)
savefig(@relpath "BLDF01_time_horizon_20.png")

𝑂_BLDC01[:mode] = "reach"
𝑂_BLDC01[:plot_vars] = [0, 25]
𝑂_BLDC01[:project_reachset] = true
res = solve(build_CONST, 𝑂_BLDC01, op=BFFPSV18(𝑂_dense_BLDC01))

plot(res,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t\raisebox{-0.5mm}{\textcolor{white}{.}}",
     ylab=L"x_{25}\raisebox{2mm}{\textcolor{white}{.}}",
     xtick=[0., 0.5, 1.], ytick=[-6e-3, -3e-3, 0., 3e-3, 6e-3],
     xlims=(0., 1.), ylims=(-Inf, 6e-3),
     bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
     size=(1000, 1000))
plot_dashed_line!(4e-3)
plot_dashed_line!(5.1e-3)
savefig(@relpath "BLDC01_time_horizon_1.png")

plot(res,
     tickfont=font(30, "Times"), guidefontsize=45,
     xlab=L"t\raisebox{-0.5mm}{\textcolor{white}{.}}",
     ylab=L"x_{25}\raisebox{2mm}{\textcolor{white}{.}}",
     xtick=[0., 10., 20.], ytick=[-6e-3, -3e-3, 0., 3e-3, 6e-3],
     bottom_margin=6mm, left_margin=2mm, top_margin=3mm,
     ylims=(-Inf, 6e-3),
     size=(1000, 1000))
plot_dashed_line!(4e-3)
plot_dashed_line!(5.1e-3)
savefig(@relpath "BLDC01_time_horizon_20.png")

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
# using Languages
using Test
using SimilaritySearch
using TextSearch
using LinearAlgebra
fit = TextSearch.fit

const text0 = "@user;) #jello.world"
const text1 = "hello world!! @user;) #jello.world :)"
const text2 = "a b c d e f g h i j k l m n o p q"
const corpus = ["hello world :)", "@user;) excellent!!", "#jello world."]

@testset "Character q-grams" begin
    config = TextConfig()
    config.group_usr = false
    config.nlist = []
    config.qlist = [3]
    config.slist = []
    a = [(p.first, p.second) for p in compute_bow(config, text0)[1]]
    b = [(:rld, 1), (Symbol("@us"), 1), (:orl, 1), (:wor, 1), (Symbol("lo."), 1), (:use, 1), (:ser, 1), (Symbol("o.w"), 1), (:llo, 1), (Symbol("r;)"), 1), (:jel, 1), (Symbol(" #j"), 1), (Symbol("er;"), 1), (Symbol("ld "), 1), (Symbol(";) "), 1), (Symbol(" @u"), 1), (:ell, 1), (Symbol(") #"), 1), (Symbol(".wo"), 1), (Symbol("#je"), 1)]
    @test sort(a) == sort(b)
end

@testset "Word n-grams" begin
    config = TextConfig()
    config.group_usr = false
    config.nlist = [1, 2]
    config.qlist = []
    config.slist = []
    a = [(p.first, p.second) for p in compute_bow(config, text0)[1]]
    b = [(Symbol("#jello"), 1), (Symbol("#jello ."), 1), (:., 1), (Symbol(". world"), 1), (Symbol(";)"), 1), (Symbol(";) #jello"), 1), (Symbol("@user"), 1), (Symbol("@user ;)"), 1), (:world, 1)]
    @test sort(a) == sort(b)
 end

@testset "Skip-grams" begin
    config = TextConfig()
    config.nlist = []
    config.qlist = []
    config.del_punc = true
    config.slist = [(2,1), (2, 2), (3, 1), (3, 2)]
    #L = collect(compute_bow(text2, config))
    #sort!(L)
    a = [(p.first, p.second) for p in compute_bow(config, text2)[1]]
    b = [(Symbol("a c"), 1), (Symbol("a c e"), 1), (Symbol("a d"), 1), (Symbol("a d g"), 1), (Symbol("b d"), 1), (Symbol("b d f"), 1), (Symbol("b e"), 1), (Symbol("b e h"), 1), (Symbol("c e"), 1), (Symbol("c e g"), 1), (Symbol("c f"), 1), (Symbol("c f i"), 1), (Symbol("d f"), 1), (Symbol("d f h"), 1), (Symbol("d g"), 1), (Symbol("d g j"), 1), (Symbol("e g"), 1), (Symbol("e g i"), 1), (Symbol("e h"), 1), (Symbol("e h k"), 1), (Symbol("f h"), 1), (Symbol("f h j"), 1), (Symbol("f i"), 1), (Symbol("f i l"), 1), (Symbol("g i"), 1), (Symbol("g i k"), 1), (Symbol("g j"), 1), (Symbol("g j m"), 1), (Symbol("h j"), 1), (Symbol("h j l"), 1), (Symbol("h k"), 1), (Symbol("h k n"), 1), (Symbol("i k"), 1), (Symbol("i k m"), 1), (Symbol("i l"), 1), (Symbol("i l o"), 1), (Symbol("j l"), 1), (Symbol("j l n"), 1), (Symbol("j m"), 1), (Symbol("j m p"), 1), (Symbol("k m"), 1), (Symbol("k m o"), 1), (Symbol("k n"), 1), (Symbol("k n q"), 1), (Symbol("l n"), 1), (Symbol("l n p"), 1), (Symbol("l o"), 1), (Symbol("m o"), 1), (Symbol("m o q"), 1), (Symbol("m p"), 1), (Symbol("n p"), 1), (Symbol("n q"), 1), (Symbol("o q"), 1)]
    @test sort(a) == sort(b)
end

@testset "Tokenizer, BOW, and vectorize" begin # test_vmodel
    config = TextConfig()
    config.nlist = [1]
    config.qlist = []
    config.slist = []
    config.group_usr = false

    @test tokenize(config, text1) == [Symbol(h) for h in ["hello", "world", "!!",  "@user", ";)", "#jello", ".", "world", ":)"]]
    model = fit(VectorModel, config, corpus)
    @test length(vectorize(model, TfModel, text1)) == 8
    @test length(vectorize(model, TfModel, text2)) == 0
end

const labeled_corpus = [("me gusta", 1), ("me encanta", 1), ("lo odio", 2), ("odio esto", 2), ("me encanta esto LOL!", 1)]
const sentiment_text = "lol, esto me encanta"

@testset "DistModel tests" begin
    config = TextConfig()
    config.nlist = [1]
    dmodel = fit(DistModel, config, [x[1] for x in labeled_corpus], [x[2] for x in labeled_corpus])
    @show sentiment_text
    @show dmodel
    #a = vectorize(dmodel, TfIdf, sentiment_text)
    #b = [(:me1,1.0),(:me2,0.0),(:encanta1,1.0),(:encanta2,0.0),(:esto1,0.4),(:esto2,0.6),(:lol1,1.0),(:lol2,0.0)]
    #@test string(a) == string(b)

end

@testset "DistModel-normalize! tests" begin
    config = TextConfig()
    config.nlist = [1]
    X = [x[1] for x in labeled_corpus]
    y = [x[2] for x in labeled_corpus]
    dmodel = fit(DistModel, config, X, y, norm_by=minimum)
    #dmap = id2token(dmodel)
    #@show sentiment_text
    #@show dmodel
    #d1 = [(dmap[t.id], t.weight) for t in vectorize(dmodel, sentiment_text).tokens]
    #d2 = [(:me1, 1.0), (:me2, 0.0), (:encanta1, 1.0), (:encanta2, 0.0), (:esto1, 0.4), (:esto2, 0.6), (:lol1, 1.0), (:lol2, 0.0)]
    #@test string(d1) == string(d2)
end

@testset "EntModel tests" begin
    config = TextConfig()
    config.nlist = [1]
    X = [x[1] for x in labeled_corpus]
    y = [x[2] for x in labeled_corpus]
    dmodel = fit(DistModel, config, X, y)
    emodel = fit(EntModel, dmodel)
    emodel_ = fit(EntModel, config, X, y)
    a = vectorize(emodel, X)
    b = vectorize(emodel_, X)
    @test 0.999 < dot(a, b)

    # a = [(emap[t.id], t.weight) for t in .tokens]
    # b = [(:esto,0.0290494),(:encanta,1.0),(:me,1.0),(:lol,1.0)]
    # sort!(a)
    # sort!(b)
    # @test string(a) == string(b)

    # @show [(maptoken[term.id], term.id, term.weight) for term in vectorize(sentiment_text, emodel).terms]
    # @show vectorize(text4, vmodel)
end
#@test
# @test TextConfig()

@testset "distances" begin
    u = Dict(:el => 0.9, :hola => 0.1, :mundo => 0.2) |> normalize!
    v = Dict(:el => 0.4, :hola => 0.2, :mundo => 0.4) |> normalize!
    w = Dict(:xel => 0.4, :xhola => 0.2, :xmundo => 0.4) |> normalize!

    dist = angle_distance
    @test dist(u, v) ≈ 0.5975474808029686
    @test dist(u, u) <= eps(Float32)
    @test dist(w, u) ≈ 1.5707963267948966
end

@testset "operations" begin
    u = Dict(:el => 0.1, :hola => 0.2, :mundo => 0.4)
    v = Dict(:el => 0.2, :hola => 0.4, :mundo => 0.8)
    w = Dict(:el => 0.1^2, :hola => 0.2^2, :mundo => 0.4^2)
    y = Dict(:el => 0.1/9, :hola => 0.2/9, :mundo => 0.4/9)
    @test u == u
    @test u != v
    @test u + u == v
    @test u * u == w
    @test u * (1/9) == y
    @test (1/9) * u == y
    @test dot(normalize!(u + v - v), normalize!(u)) > 0.99
end

@testset "io" begin
    buff = IOBuffer("""{"key1": "value1a", "key2c": "value2a"}
{"key1": "value1b", "key2c": "value2b"}
{"key1": "value1c", "key2b": "value2c"}
{"key1": "value1d", "key2a": "value2d"}""")
    itertweets(buff) do x
        @info x
    end
end

_corpus = [
    "la casa roja",
    "la casa verde",
    "la casa azul",
    "la manzana roja",
    "la pera verde esta rica",
    "la manzana verde esta rica",
    "la hoja verde",
]

@testset "transpose bow" begin
    config = TextConfig()
    config.nlist = [1]
    config.qlist = []
    config.slist = []
    model = fit(VectorModel, config, _corpus)
    @show _corpus
    X = [vectorize(model, FreqModel, x) for x in _corpus]
    dX = transpose(X)
end

@testset "invindex" begin
    config = TextConfig()
    config.qlist = [3, 4]
    config.nlist = [1, 2, 3]

    model = fit(VectorModel, config, _corpus)
    invindex = InvIndex()
    for c in _corpus
        push!(invindex, invindex.n + 1, vectorize(model, TfidfModel, c))
    end

    q = vectorize(model, TfidfModel, "la casa roja")
    res = search(invindex, cosine_distance, q, KnnResult(4))
    ires = [r.objID for r in res]
    @show res, _corpus[ires]
    @test ires == [1, 2, 4, 3]
    shortindex = prune(invindex, 3)
    res = search(shortindex, cosine_distance, q, KnnResult(4))
    ires = [r.objID for r in res]
    @show res, _corpus[ires]
    @test sort(ires) == [1, 2, 3, 4]
end

@testset "centroid computing" begin
    config = TextConfig()
    config.nlist = [1]
    config.qlist = []
    config.slist = []
    model = fit(VectorModel, config, _corpus)
    @show _corpus
    X = [vectorize(model, FreqModel, x) for x in _corpus]
    x = sum(X) |> normalize!
    @test 0.999 < dot(x, Dict(:la=>0.736665,:verde=>0.39922,:azul=>0.112482,:pera=>0.087128,:esta=>0.174256,:roja=>0.224964,:hoja=>0.112482,:casa=>0.337445,:rica=>0.174256,:manzana=>0.19961))
end

@testset "rocchio" begin
    config = TextConfig()
    config.nlist = [1,2]
    config.qlist = [3]
    config.slist = []
    corpus = [x[1] for x in labeled_corpus]

    model = fit(VectorModel, config, corpus)
    X = [vectorize(model, TfidfModel, x) for x in corpus]
    y = [x[2] for x in labeled_corpus]
    rocchio = fit(Rocchio, X, y)
    @show rocchio.protos rocchio.pops
    @test sum(predict.(rocchio, X) .== y)/length(y) == 1.0

    for p in transform.(rocchio, X)
        println(p)
    end
end

@testset "rocchio bagging" begin
    config = TextConfig()
    config.nlist = [1, 2]
    config.qlist = [3]
    config.slist = []
    corpus = [x[1] for x in labeled_corpus]
    model = fit(VectorModel, config, corpus)
    X = [vectorize(model, TfidfModel, x) for x in corpus]
    y = [x[2] for x in labeled_corpus]
    rocchio = fit(RocchioBagging, X, y)
    for p in transform.(rocchio, X)
        println("===>", p)
    end

    #@show rocchio.protos rocchio.pops
    @show predict.(rocchio, X)
    @show y
end

@testset "neardup" begin
    config = TextConfig()
    config.nlist = [1]
    config.qlist = []
    config.slist = []

    function create_corpus()
        alphabet = Char.(97:100)
        corpus = String[]
        L = Char[]
        for i in 1:10000
            resize!(L, 0)
            for j in 1:7
                append!(L, rand(alphabet, 2))
                push!(L, ' ')
            end

            push!(corpus, join(L))
        end

        corpus
    end

    corpus = create_corpus()
    model = fit(VectorModel, config, corpus)
    @show corpus[1:10]
    X = [vectorize(model, TfModel, x) for x in corpus]
    L, D  = neardup(X, 0.2)
    @test length(X) > length(unique(L))
end


text:
module test_jive_runtests_target2

using Test

@test true

end # module test_jive_runtests_target2


text:
"""
    load_file_line_by_line(filename::AbstractString) -> Array{SubString{String},1}
Load and return the content of the text file `filename` in `Array{SubString{String},1}`, i.e., one element for each line.
# Arguments
- `filename::AbstractString`: the filename of the file to load,
# Keywords
# Returns
- `Array{SubString{String},1}`: the array of the file content,
# Throws
- `Error`: in the case of file does not exist.
"""
function load_file_line_by_line(filename::AbstractString)
    if !isfile(filename)
        error("File $(filename) does not exist.")
    end
    f=open(filename) do file
        read(file,String)
    end
    return split(f,"\n")
end


text:
module Autoencoder

export trainandsave, load

using Base.Iterators: partition
using BSON
using BSON: @save
using CUDA
using Images
using NNlib
using Parameters
using Flux
using Flux.Data.MNIST
using Flux: @epochs, mse, throttle
using Statistics

@with_kw mutable struct Args
    α::Float64 = 1e-3      # learning rate
    epochs::Int = 20       # number of epochs
    N::Int = 64            # size of the encoding (i.e. hidden layer)
    batchsize::Int = 100   # batch size for training
    num::Int = 20          # number of random digits in the sample image (UNUSED)
    throttle::Int = 1      # throttle timeout (called once every X seconds)
end

global X = MNIST.images(:test) # train using the same validation dataset

function get_processed_data(args)
    # load images and convert image of type RBG to Float
    imgs = channelview.(X)
    # partition into batches of size `batchsize`
    batches = partition(imgs, args.batchsize)
    traindata = [float(hcat(vec.(imgs)...)) for imgs in batches]
    return gpu.(traindata)
end

function train(; kwargs...)
    args = Args(; kwargs...)

    traindata = get_processed_data(args)

    @info "Constructing model..."

    # You can try to make the encoder/decoder network larger
    # Also, the output of encder is a coding of the given input.
    # In this case, the input dimension is 28^2 and the output dimension of the
    # encoder is 32. This implies that the coding is a compressed representation.
    # We can make lossy compression via this `encoder`.
    intermediate = 28^2 ÷ 2
    encoder = Chain(
        Dense(28^2, intermediate, leakyrelu) |> gpu,
        Dense(intermediate, args.N, leakyrelu) |> gpu)
    decoder = Chain(
        Dense(args.N, intermediate, leakyrelu) |> gpu,
        Dense(intermediate, 28^2, leakyrelu) |> gpu)

    # define main model as a Chain of encoder and decoder models
    m = Chain(encoder, decoder)

    @info "Training model..."
    loss(x) = mse(m(x), x)

    # callback, optimizer, and training
    callback = throttle(() -> @show(mean(loss.(traindata))), args.throttle)
    opt = ADAM(args.α)

    @epochs args.epochs Flux.train!(loss, params(m), zip(traindata), opt, cb=callback)

    return m, encoder, decoder
end

img(x::Vector) = Gray.(reshape(clamp.(x, 0, 1), 28, 28))

"""
Sample a compressed image from the model.

```julia
sample(m, 20)
```
"""
function sample(m, num=1)
    # convert image of type RGB to Float
    imgs = channelview.(X)
    # number of random digits (truth)
    before = [imgs[i] for i in rand(1:length(imgs), num)]
    # after applying autoencoder to `before` input image
    after = img.(map(x->cpu(m)(float(vec(x))), before))
    # stack `before` and `after` images them all together
    Gray.(hcat(vcat.(before, after)...))
end

function trainandsave(; kwargs...)
    m, encoder, decoder = train(; kwargs...)
    m = cpu(m)
    @save "models/autoencoder.bson" m
    encoder = cpu(encoder)
    @save "models/encoder.bson" encoder
    decoder = cpu(decoder)
    @save "models/decoder.bson" decoder
end

function load(model_filename_bson, sym=:m)
    model = BSON.load(model_filename_bson)[sym]
    return model
end

end # module


text:
# ---------- #
# State type #
# ---------- #

struct State
    board::Board
    player::Int
    round::Int
end
State() = State(Board(0, 0), 1, 0)

function Base.show(io::IO, mime::MIME"text/plain", s::State)
    if s.player == 1 || s.player == 2
        msg = "It is Player $(s.player)'s ($(MARKERS[s.player])) turn"
    else
        msg = "It is player $(s.player)'s turn"
    end

    println(io, msg, " on round $(s.round)")
    println(io, "The score is $(score(s)) and the board looks like:\n")
    show(io, mime, s.board)
end
Base.show(io::IO, s::State) = show(io, MIME"text/plain"(), s)

# --------------------------- #
# bitboard api for State type #
# --------------------------- #

function legal_actions_bits(s::State)
    if s.player == 1
        return legal_actions_bits(s.board.p1_placed, s.board.p2_placed)
    else
        return legal_actions_bits(s.board.p2_placed, s.board.p1_placed)
    end
end

score(s::State) = score(s.board.p1_placed), score(s.board.p2_placed)
legal_actions(s::State) = decode(legal_actions_bits(s))
function game_over(s::State)
    if s.round  < 4
        return false, 0, 0
    end
    game_over(s.board.p1_placed, s.board.p2_placed)
end

function do_action(s::State, action::UInt64)::State
    p1 = s.board.p1_placed
    p2 = s.board.p2_placed
    if s.player == 1
        new_p1, new_p2 = next_state_bits(action, p1, p2)
        new_player = 2
    else
        new_p2, new_p1 = next_state_bits(action, p2, p1)
        new_player = 1
    end
    State(Board(new_p1, new_p2), new_player, s.round+1)
end

do_action(s::State, a::Tuple{Int,Int}) = do_action(s, encode(a))
do_pass(s::State) = State(s.board, (s.player == 1) + 1, s.round)

opponent(i::Integer) = (i == 1) + 1
opponent(s::State) = opponent(state.player)


text:
# TODO: define VertexDataGraph, a graph with only data on the
# vertices, and EdgeDataGraph, a graph with only data on the edges.
struct DataGraph{VD,ED,V,E,G<:AbstractGraph} <: AbstractDataGraph{VD,ED,V,E}
  underlying_graph::G
  vertex_data::Dictionary{V,VD}
  edge_data::Dictionary{E,ED}
end
underlying_graph(graph::DataGraph) = graph.underlying_graph
vertex_data(graph::DataGraph) = graph.vertex_data
edge_data(graph::DataGraph) = graph.edge_data

function copy(graph::DataGraph)
  return DataGraph(
    copy(underlying_graph(graph)), copy(vertex_data(graph)), copy(edge_data(graph))
  )
end

function is_directed(::Type{<:DataGraph{VD,ED,V,E,G}}) where {VD,ED,V,E,G}
  return is_directed(G)
end

function DataGraph{VD,ED}(
  underlying_graph::AbstractGraph,
  vertex_data::Dictionary{V,VD},
  edge_data::Dictionary{E,ED},
) where {VD,ED,V,E}
  G = typeof(underlying_graph)
  return DataGraph{VD,ED,V,E,G}(underlying_graph, vertex_data, edge_data)
end

function DataGraph{VD,ED}(underlying_graph::AbstractGraph) where {VD,ED}
  V = eltype(underlying_graph)
  E = edgetype(underlying_graph)
  vertex_data = Dictionary{V,VD}()
  edge_data = Dictionary{E,ED}()
  return DataGraph{VD,ED}(underlying_graph, vertex_data, edge_data)
end

DataGraph(underlying_graph::AbstractGraph) = DataGraph{Any,Any}(underlying_graph)
function DataGraph{VD}(underlying_graph::AbstractGraph) where {VD}
  return DataGraph{VD,Any}(underlying_graph)
end


text:
function pluralize end

function pluralize(word::String,lang_code::String)
    lang_code = lowercase(lang_code)
    if lang_code == "eng"
        return Pluralize.ENG.pluralize(word)
    elseif lang_code == "fra"
        return Pluralize.FRA.pluralize(word)
    else
        error("Unsupported language[$lang_code]")
    end
end


text:
using JuMP, Mosek, MosekTools, LinearAlgebra, Suppressor, StatsBase#, SCS#, Ipopt

# Solve the perspective relaxation (naive but cheaper formulation)
function solvepersprelax(X, y, m, n, p, lambda, gamma)
    mod=Model(Mosek.Optimizer)

    @variable(mod, beta[1:p, 1:n])
    @variable(mod, W[1:n, 1:n], Symmetric)
    @variable(mod, theta[1:p, 1:p], Symmetric)
    @variable(mod, t) # Epigraph variable for rotated SOC constraint

    @constraint(mod, Symmetric(Matrix(1.0I, n, n)-W) in PSDCone());

    @constraint(mod, Symmetric([theta beta; beta' W]) in PSDCone());

    @constraint(mod, [t; 1.0; vec(y-X*beta)] in RotatedSecondOrderCone());

    @objective(mod, Min, (0.5/m)*t+0.5*1.0/(2.0*gamma)*sum(theta[i,i] for i=1:p)+lambda*sum(W[i,i] for i=1:n))

    @suppress optimize!(mod)

    ofv_dual=objective_value(mod)
    u,v,=svd(value.(W))
    # Use a cutoff of 10^-4 for singular values, to avoid e.g. Julia counting singular values of 10^-12 as part of the rank
    ofv_primal=(0.5/m)*LinearAlgebra.dot(y.-X*value.(beta), y.-X*value.(beta))+1.0/(2.0*gamma)*norm(value.(beta))^2+lambda*sum(abs.(v).>1e-4)
    # Todo: round greedily and resolve second stage problem??

    return value.(beta), ofv_dual, ofv_primal, value.(W)
end

# Solve the perspective relaxation (expensive DCL version)
function solvepersprelax_DCL(X, y, m, n, p, lambda, gamma)
    mod=Model(Mosek.Optimizer)

    @variable(mod, beta[1:p, 1:n])
    @variable(mod, B[1:p, 1:p], Symmetric)
    @variable(mod, W[1:n, 1:n], Symmetric)

    @constraint(mod, Symmetric(Matrix(1.0I, n, n)-W) in PSDCone());

    @constraint(mod, Symmetric([B beta; beta' W]) in PSDCone());

    @objective(mod, Min, (0.5/m)*sum(y[i,j]^2 for i=1:m for j=1:n)-(1.0/m)*LinearAlgebra.dot(y, X*beta)+0.5*LinearAlgebra.dot(B, 1.0/gamma*Matrix(1.0I, p, p)+(1.0/m)X'*X)+lambda*sum(W[i,i] for i=1:n))

    @suppress optimize!(mod)

    @show ofv_dual=objective_value(mod)

    u,v,=svd(value.(W))
    # Use a cutoff of 10^-4 for singular values, to avoid e.g. Julia counting singular values of 10^-12 as part of the rank
    ofv_primal=(0.5/m)*LinearAlgebra.dot(y.-X*value.(beta), y.-X*value.(beta))+1.0/(2.0*gamma)*norm(value.(beta))^2+lambda*sum(abs.(v).>1e-4)

    return value.(beta), ofv_dual, ofv_primal, value.(W)
end

# Solve the nuclear norm "relaxation" (not actually a relaxation without big-M constraints and projection matrices)
function solveNNrelax(X, y, m, n, p, lambda, gamma)
    mod=Model(Mosek.Optimizer)

    @variable(mod, beta[1:p, 1:n])
    @variable(mod, t) # Epigraph variable for rotated SOC constraint
    @variable(mod, U[1:p, 1:p])
    @variable(mod, V[1:n, 1:n])

    @constraint(mod, [t; 1.0; vec(y-X*beta)] in RotatedSecondOrderCone());

    @constraint(mod, Symmetric([U beta; beta' V]) in PSDCone());

    @objective(mod, Min, (0.5/m)*t+1.0/(2.0*gamma)*sum(beta[i,j]^2 for i=1:p for j=1:n)+0.5*lambda*sum(U[i,i] for i=1:p)+0.5*lambda*sum(V[i,i] for i=1:n))

    @suppress optimize!(mod)

    # Todo: round greedily and resolve second stage problem??

    return value.(beta), objective_value(mod)
end

max_stars_repo_path: Matrix Regression/convexpenalties.jl

max_issues_repo_path: Matrix Regression/convexpenalties.jl

max_forks_repo_path: Matrix Regression/convexpenalties.jl


text:
using SafeTestsets

const GROUP = get(ENV, "GROUP", "All")
const is_APPVEYOR = Sys.iswindows() && haskey(ENV,"APPVEYOR")
const is_TRAVIS = haskey(ENV,"TRAVIS")

#Start Test Script

@time begin

if GROUP == "All" || GROUP == "OperatorInterface"
    @time @safetestset "Poisson example" begin include("DerivativeOperators/poisson.jl") end
    @time @safetestset "Heat equation example" begin include("DerivativeOperators/heat_equation.jl") end
    @time @safetestset "Robin Boundary Condition Operators" begin include("DerivativeOperators/robin.jl") end
    @time @safetestset "Composite Operators Interface" begin include("DerivativeOperators/composite_operators_interface.jl") end
    @time @safetestset "BC and Coefficient Compositions" begin include("DerivativeOperators/bc_coeff_compositions.jl") end
    @time @safetestset "Derivative Operators Interface" begin include("DerivativeOperators/derivative_operators_interface.jl") end
    @time @safetestset "Ghost Derivative Operators Interface" begin include("DerivativeOperators/ghost_derivative_operators_interface.jl") end
    @time @safetestset "Validate Regular Derivative Operators" begin include("DerivativeOperators/regular_operator_validation.jl") end
    @time @safetestset "Validate and Compare Generic Operators" begin include("DerivativeOperators/generic_operator_validation.jl") end
    @time @safetestset "Validate Boundary Padded Array Concretization" begin include("DerivativeOperators/boundary_padded_array.jl") end
    #@time @safetestset "Validate Higher Dimensional Boundary Extension" begin include("DerivativeOperators/multi_dim_bc_test.jl") end
    @time @safetestset "2nd order check" begin include("DerivativeOperators/2nd_order_check.jl") end
    @time @safetestset "Non-linear Diffusion" begin include("DerivativeOperators/Fast_Diffusion.jl") end
    @time @safetestset "KdV" begin include("DerivativeOperators/KdV.jl") end # 2-Soliton case needs implementation
    @time @safetestset "Heat Equation" begin include("DerivativeOperators/heat_eqn.jl") end
    @time @safetestset "Matrix-Free Operators" begin include("DerivativeOperators/matrixfree.jl") end
    @time @safetestset "Convolutions" begin include("DerivativeOperators/convolutions.jl") end
    @time @safetestset "Differentiation Dimension" begin include("DerivativeOperators/differentiation_dimension.jl") end
    @time @safetestset "Higher Dimensional Concretization" begin include("DerivativeOperators/concretization.jl") end
    @time @safetestset "Coefficient Functions" begin include("DerivativeOperators/coefficient_functions.jl") end
    @time @safetestset "Upwind Operator Interface" begin include("DerivativeOperators/upwind_operators_interface.jl") end
    @time @safetestset "Basic SDO Examples" begin include("DerivativeOperators/BasicSDOExamples.jl") end
    @time @safetestset "3D laplacian Test" begin include("DerivativeOperators/3D_laplacian.jl") end
    # @time @safetestset "Linear Complementarity Problem Examples" begin include("DerivativeOperators/lcp.jl"); include("DerivativeOperators/lcp_split.jl") end
end

if GROUP == "All" || GROUP == "MOLFiniteDifference"
    @time @safetestset "MOLFiniteDifference Interface" begin include("MOL/MOLtest.jl") end
    #@time @safetestset "MOLFiniteDifference Interface: Linear Convection" begin include("MOL/MOL_1D_Linear_Convection.jl") end
    @time @safetestset "MOLFiniteDifference Interface: 1D Linear Diffusion" begin include("MOL/MOL_1D_Linear_Diffusion.jl") end
    @time @safetestset "MOLFiniteDifference Interface: 1D Non-Linear Diffusion" begin include("MOL/MOL_1D_NonLinear_Diffusion.jl") end
    @time @safetestset "MOLFiniteDifference Interface: 2D Diffusion" begin include("MOL/MOL_2D_Diffusion.jl") end
    @time @safetestset "MOLFiniteDifference Interface: 1D HigherOrder" begin include("MOL/MOL_1D_HigherOrder.jl") end
    @time @safetestset "MOLFiniteDifference Interface: 1D Partial DAE" begin include("MOL/MOL_1D_PDAE.jl") end
end

if GROUP == "All" || GROUP == "Misc"
    @time @safetestset "Utilities Tests" begin include("Misc/utils.jl") end
    @time @safetestset "JacVec Operators Interface" begin include("Misc/jacvec_operators.jl") end
    @time @safetestset "JacVec Operator Integration Test" begin include("Misc/jacvec_integration_test.jl") end
end

if !is_APPVEYOR && (GROUP == "All" || GROUP == "Multithreading")
    @time @safetestset "2D and 3D fast multiplication" begin include("DerivativeOperators/2D_3D_fast_multiplication.jl") end
end

if GROUP == "GPU"

end
end


text:
"""
    Client

    Client(url::URIs.URI)
    Client(url::AbstractString)
    Client(host, port, scheme)

Druid queryable host connection information store.

Do not provide the endpoint, it will be set based on the type of query. For
example, "druid/v2/sql" for SQL queries and "druid/v2" for all other queries.

# Examples
```julia-repl
julia> client = Client("http://localhost:8888/")
Client(URI("http://localhost:8888/"))

julia> client = Client("localhost", 8082, "http")
Client(URI("http://localhost:8082"))
```
"""
struct Client
    url::URI
end
Client(host, port, scheme) =
    Client(URI(;scheme=scheme, host=host, port=port))
Client(url::AbstractString) = Client(URI(url))


text:
struct RadiationPattern
    gains_db::Matrix{Float64}   # gains in decibels (I think)
    gains::Matrix{Float64}      # gains in real units
    thetas::Vector{Float64}     # azimuth angles (I think)
    phis::Vector{Float64}       # elevation angles (I think)
end

# the following is just so we can pass in reals or ranges
RealOrRange = Union{Real,Range}
istart(r::Real) = r
istart(r::Range) = start(r)
Base.step(r::Real) = 1       # doesn't matter, could be any number

function RadiationPattern(a::AbstractAntenna, fmhz::Real;
                          thetas::RealOrRange=0:180,
                          phis::RealOrRange=0:360
                         )
    vs = VoltageSource(a)
    return RadiationPattern(a.wires, vs, fmhz; thetas=thetas, phis=phis)
end

# returns:
#  max_db = maximum gain, in db
#  theta = inclination in degrees (0 along +ve z-axis, to zenith)
#  phi = azimuth in degrees (0 along +ve x-axis, 90 along +ve y-axis)
function Base.max(rp::RadiationPattern)
    max_db = maximum(rp.gains_db)
    max_idx = indmax(rp.gains_db)
    theta_idx, phi_idx = ind2sub( size(rp.gains_db), max_idx )

    theta = rp.thetas[theta_idx]
    phi = mod(rp.phis[phi_idx], 360.0)      # ensure 360.0 treated as zero

    return max_db, theta, phi
end

# TODO: find the horizontal plane, don't just assume it is at index 91
function hpat(rp::RadiationPattern)
    return rp.thetas, vec(rp.gains_db[91,:])
end

# TODO: dont just assume it is at index 91
# TODO: this is a train wreck
function f2b(rp::RadiationPattern)
    thetas, gains_db = hpat(rp)
    return gains_db[1] , gains_db[181]
end


text:
using Compat

data_dir = joinpath(Pkg.dir("Word2Vec"), "test", "data")
model_dir = joinpath(Pkg.dir("Word2Vec"), "test", "models")

test_filename = isempty(ARGS) ? "text8" : ARGS[1]
test_file = joinpath(data_dir, test_filename)
model_file = joinpath(model_dir, test_filename * ".model")

np = (length(ARGS) > 1) ? @compat(parse(Int, ARGS[2])) : 4
addprocs(np)

@everywhere using Word2Vec
@everywhere using Base.Test
@everywhere using Blocks
@everywhere using Base.FS
@everywhere using Compat

function parallel_word_embedding(filename, savemodel="")
    b = Block(File(filename), nworkers())
    embed = WordEmbedding(100, Word2Vec.random_inited, Word2Vec.huffman_tree, subsampling = 1e-4)
    @time train(embed, b, "/tmp/emb")
    isempty(savemodel) || save(embed, savemodel)
    embed
end

function test_nearest_words(embed)
    println("nearest words:")
    inp = ["king", "queen", "prince"]
    for w in inp
        nwords = find_nearest_words(embed, w)
        println("\t$w - $(join(nwords,','))")
    end
    nwords = find_nearest_words(embed, join(inp, ' '))
    println("\t$(join(inp,',')) - $(join(nwords,','))")
end

function interactive_nearest_words(embed)
    cont = true
    while cont
        println("\nEnter words to find nearest word:")
        l = strip(readline(STDIN))
        nwords = find_nearest_words(embed, strip(l))
        println("\t$(l) - $(join(nwords,','))")
    end
end

function restore_or_train()
    if isfile(model_file)
        restore(model_file)
    else
        parallel_word_embedding(test_file, model_file)
    end
end

embed = restore_or_train()
show(embed)
test_nearest_words(embed)
interactive_nearest_words(embed)


text:
using Rasters, Test, ArchGDAL, Dates, Statistics, GeoInterface
using Rasters.LookupArrays, Rasters.Dimensions
using GeoInterface: Point, Polygon, LineString

include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))

A = [missing 7.0f0; 2.0f0 missing]
B = [1.0 0.4; 2.0 missing]
ga = Raster(A, (X, Y); missingval=missing)
ga99 = replace_missing(ga, -9999)
gaNaN = replace_missing(ga, NaN32)
gaMi = replace_missing(ga)
st = RasterStack((a=A, b=B), (X, Y); missingval=(a=missing,b=missing))

polygon = [[-20.0, 30.0],
           [-20.0, 10.0],
           [0.0, 10.0],
           [0.0, 30.0],
           [-20.0, 30.0]]
vals = [1, 2, 3, 4, 5]

@testset "replace_missing" begin
    @test all(isequal.(ga99, [-9999.0f0 7.0f0; 2.0f0 -9999.0f0]))
    @test missingval(ga99) === -9999.0f0
    @test all(isequal.(gaNaN, [NaN32 7.0f0; 2.0f0 NaN32]))
    @test missingval(gaNaN) === NaN32
    @test all(isequal.(gaMi, ga))
    @test missingval(gaMi) === missing
    # The second layer NaN32 will be promoted to NaN to match the array type
    @test missingval(replace_missing(st, NaN32)) === (a=NaN32, b=NaN)
    @test all(map(values(replace_missing(st, NaN32)), (a=[NaN32 7.0f0; 2.0f0 NaN32], b=[1.0 0.4; 2.0 NaN])) do x, y
        all(x .=== y)
    end)
    dNaN = replace_missing(ga, NaN32; filename="test.tif")
    @test all(isequal.(dNaN, [NaN32 7.0f0; 2.0f0 NaN32]))
    @test Rasters.isdisk(dNaN)
    rm("test.tif")
    stNaN = replace_missing(st, NaN32; filename="teststack.tif")
    @test all(map(stNaN[Band(1)], (a=[NaN32 7.0f0; 2.0f0 NaN32], b=[1.0 0.4; 2.0 NaN])) do x, y
        all(x .=== y)
    end)
    rm("teststack_a.tif")
    rm("teststack_b.tif")
end

@testset "boolmask" begin
    @test boolmask(ga) == [false true; true false]
    @test boolmask(ga99) == [false true; true false]
    @test boolmask(gaNaN) == [false true; true false]
    @test dims(boolmask(ga)) == (X(NoLookup(Base.OneTo(2))), Y(NoLookup(Base.OneTo(2))))
end

@testset "missingmask" begin
    @test all(missingmask(ga) .=== [missing true; true missing])
    @test all(missingmask(ga99) .=== [missing true; true missing])
    @test all(missingmask(gaNaN) .=== [missing true; true missing])
    @test dims(missingmask(ga)) == (X(NoLookup(Base.OneTo(2))), Y(NoLookup(Base.OneTo(2))))
end

@testset "mask" begin
    A1 = [missing 1; 2 3]
    A2 = view([0 missing 1; 0 2 3], :, 2:3)
    ga1 = Raster(A1, (X, Y); missingval=missing)
    ga2 = Raster(A2, (X, Y); missingval=missing)
    @test all(mask(ga1; with=ga) .=== mask(ga2; with=ga) .=== [missing 1; 2 missing])
    ga2 = replace_missing(ga1 .* 1.0; missingval=NaN)
    @test all(mask(ga2; with=ga) .=== [NaN 1.0; 2.0 NaN])
    ga3 = replace_missing(ga1; missingval=-9999)
    mask!(ga3; with=ga)
    @test all(ga3 .=== [-9999 1; 2 -9999])
    dmask = mask(ga3; with=ga, filename="mask.tif")
    @test Rasters.isdisk(dmask)
    rm("mask.tif")
    stmask = mask(replace_missing(st, NaN); with=ga, filename="mask.tif")
    @test Rasters.isdisk(stmask)
    rm("mask_a.tif")
    rm("mask_b.tif")
    @testset "to polygon" begin
        a1 = Raster(ones(X(-20:5), Y(0:30)))
        st1 = RasterStack(a1, a1)
        ser1 = RasterSeries([a1, a1], Ti(1:2))
        @test all(mask(a1; with=polygon) .===
            mask(st1; with=polygon)[:layer1] .===
            mask(ser1; with=polygon)[1]
        )
        # TODO: investigate this more for Points/Intervals
        # Exactly how do we define when boundary values are inside/outside a polygon
        @test sum(skipmissing(mask(a1; with=polygon, boundary=:inside))) == 19 * 19
        @test sum(skipmissing(mask(a1; with=polygon, boundary=:center))) == 20 * 20
        @test sum(skipmissing(mask(a1; with=polygon, boundary=:touches))) == 21 * 21
    end
end

@testset "classify" begin
    A1 = [missing 1; 2 3]
    ga1 = Raster(A1, (X, Y); missingval=missing)
    @test all(classify(ga1, 1=>99, 2=>88, 3=>77) .=== [missing 99; 88 77])
    @test all(classify(ga1, 1=>99, 2=>88, 3=>77; others=0) .=== [missing 99; 88 77])
    @test all(classify(ga1, 1=>99, 2=>88; others=0) .=== [missing 99; 88 0])
    A2 = [1.0 2.5; 3.0 4.0]
    ga2 = Raster(A2, (X, Y); missingval=missing)
    @test classify(ga2, (2, 3)=>:x, >(3)=>:y) == [1.0 :x; 3.0 :y]
    @test classify(ga2, (>=(1), <(2))=>:x, >=(3)=>:y) == [:x 2.5; :y :y]
    classify!(ga2, (1, 2.5)=>0.0, >=(3)=>-1.0; lower=(>), upper=(<=))
    @test ga2 == [1.0 0.0; -1.0 -1.0]
    classify!(ga2, [1 2.5 0.0; 2.5 4.0 -1.0]; lower=(>), upper=(<=))
    @test ga2 == [1.0 0.0; -1.0 -1.0]
    @test all(classify(ga1, [1 99; 2 88; 3 77]) .=== [missing 99; 88 77])
    @test_throws ArgumentError classify(ga1, [1, 2, 3])
end

@testset "points" begin
    ga = Raster(A, (X(9.0:1.0:10.0), Y(0.1:0.1:0.2)); missingval=missing)
    @test all(
              collect(
                      points(ga; order=(Y, X)))
              .=== [missing (0.2, 9.0); (0.1, 10.0) missing])
    @test all(collect(points(ga; order=(X, Y))) .=== [missing (9.0, 0.2); (10.0, 0.1) missing])
    @test all(points(ga; order=(X, Y), ignore_missing=true) .===
              [(9.0, 0.1) (9.0, 0.2); (10.0, 0.1) (10.0, 0.2)])
end

@testset "extract" begin
    A1 = [1 2; 3 4]
    A2 = [5 6; 7 8]
    dimz = (X(9.0:1.0:10.0), Y(0.1:0.1:0.2))
    ga = Raster(A1, dimz; name=:test, missingval=missing)
    ga2 = Raster(A2, dimz; name=:test2, missingval=missing)
    st = RasterStack(ga, ga2)
    @testset "from Raster" begin
        @test all(extract(ga, [missing, [9.0, 0.1], [10.0, 0.2], [10.0, 0.3]]) .===
                  [missing, (X=9.0, Y=0.1, test=1), (X=10.0, Y=0.2, test=4), (X=10.0, Y=0.3, test=missing)])
        @test all(extract(ga, ([9.0, 10.0, 10.0], [0.1, 0.2, 0.3])) .===
                  [(X=9.0, Y=0.1, test=1), (X=10.0, Y=0.2, test=4), (X=10.0, Y=0.3, test=missing)])
        @test all(extract(ga, Polygon([[9.0, 0.1], [10.0, 0.2], [10.0, 0.3]])) .===
                  [(X=9.0, Y=0.1, test=1), (X=10.0, Y=0.2, test=4), (X=10.0, Y=0.3, test=missing)])
        @test all(extract(ga, [(0.1, 9.0), (0.2, 10.0), (0.3, 10.0), missing]; order=(Y, X)) .===
                  [(Y=0.1, X=9.0, test=1), (Y=0.2, X=10.0, test=4), (Y=0.3, X=10.0, test=missing), missing])
        @test all(extract(ga, ([0.1, 0.2, 0.3], [9.0, 10.0, 10.0]); order=(Y, X)) .===
                  [(Y=0.1, X=9.0, test=1), (Y=0.2, X=10.0, test=4), (Y=0.3, X=10.0, test=missing)])
        @test all(extract(ga, Polygon([[0.1, 9.0], [0.2, 10.0], [0.3, 10.0]]); order=(Y, X)) .===
                  [(Y=0.1, X=9.0, test=1), (Y=0.2, X=10.0, test=4), (Y=0.3, X=10.0, test=missing)])
        @test all(extract(ga, [(missing, 9.0), (0.2, 10.0), (0.3, 10.0)]; order=(Y, X)) .===
                  [(Y=missing, X=missing, test=missing), (Y=0.2, X=10.0, test=4), (Y=0.3, X=10.0, test=missing)])
        @test all(extract(ga, [[missing, 9.0], [0.2, 10.0], [0.3, 10.0]]; order=(Y, X)) .===
                  [(Y=missing, X=missing, test=missing), (Y=0.2, X=10.0, test=4), (Y=0.3, X=10.0, test=missing)])
        @test_throws ArgumentError extract(ga, [missing, [9.0, 0.1], [10.0, 0.2], [10.0, 0.3]]; order=())
        @test_throws ArgumentError extract(ga, [missing, [9.0, 0.1], [10.0, 0.2], [10.0, 0.3]]; order=(Z, X))
    end
    @testset "from stack" begin
        @test all(extract(st, [missing, [9.0, 0.1], [10.0, 0.2], [10.0, 0.3]]) .===
              [missing, (X=9.0, Y=0.1, test=1, test2=5), (X=10.0, Y=0.2, test=4, test2=8), (X=10.0, Y=0.3, test=missing, test2=missing)])
    end
    @testset "Tables.jl compatible" begin
        @test all(extract(ga, [(X=9.0, Y=0.1), (X=10.0, Y=0.2), (X=10.0, Y=0.3)]) .===
              [(X=9.0, Y=0.1, test=1), (X=10.0, Y=0.2, test=4), (X=10.0, Y=0.3, test=missing)])
        @test all(extract(ga, [(X=9.0, Y=0.1), (X=10.0, Y=0.2), (X=10.0, Y=0.3)]; order=(Y=>:Y, X=>:X)) .===
              [(Y=0.1, X=9.0, test=1), (Y=0.2, X=10.0, test=4), (Y=0.3, X=10.0, test=missing)])
    end
end

@testset "trim, crop, extend" begin
    A = [missing missing missing
         missing 2.0     0.5
         missing 1.0     missing]

    r_fwd = Raster(A, (X(1.0:1.0:3.0), Y(1.0:1.0:3.0)); missingval=missing)
    r_revX = reverse(r_fwd; dims=X)
    r_revY = reverse(r_fwd; dims=Y)

    ga = r_revY
    for ga in (r_fwd, r_revX, r_revY)
        # Test with missing on all sides
        ga_r = rot180(ga)
        trimmed = trim(ga)
        trimmed_r = trim(ga_r)
        trimmed_d = trim(ga_r)
        if ga === r_fwd
            @test all(trimmed .=== [2.0 0.5; 1.0 missing])
            @test all(trimmed_r .=== [missing 1.0; 0.5 2.0])
        end
        cropped = crop(ga; to=trimmed)
        _, cropped2 = crop(trimmed, ga)
        cropped_r = crop(ga_r; to=trimmed_r)
        @test all(cropped .=== trimmed)
        @test all(cropped_r .=== trimmed_r)
        extended = extend(cropped, ga)[1]
        extended_r = extend(cropped_r; to=ga_r)
        extended_d = extend(cropped; to=ga, filename="extended.tif")
        @test all(extended .=== ga)
        @test all(extended_r .=== ga_r)

        @testset "to polygons" begin
            A1 = Raster(zeros(X(-20:-5; sampling=Points()), Y(0:30; sampling=Points())))
            A2 = Raster(ones(X(-20:-5; sampling=Points()), Y(0:30; sampling=Points())))
            A1crop1 = crop(A1; to=polygon)
            A1crop2, A2crop = crop(A1, A2; to=polygon)
            size(A1crop1)
            @test size(A1crop1) == size(A1crop2) == size(A2crop) == (16, 21)
            @test bounds(A1crop1) == bounds(A1crop2) == bounds(A2crop) == ((-20, -5), (10, 30))
            A1extend1 = extend(A1; to=polygon)
            A1extend2, A2extend = extend(A1, A2; to=polygon)
            @test all(A1extend1 .=== A1extend2)
            @test size(A1extend1) == size(A1extend2) == size(A2extend) == (21, 31)
            @test bounds(A1extend1) == bounds(A1extend2) == bounds(A2extend) == ((-20.0, 0.0), (0, 30))
        end
        @testset "to table" begin
            A1 = Raster(zeros(X(-20:-5; sampling=Points()), Y(0:30; sampling=Points())))
            table = map(polygon, vals) do p, v
                (x=p[1], y=p[2], val1=v, val2=2.0f0v)
            end
            A1crop1 = crop(A1; to=table, order=(X=>:x, Y=>:y))
        end
    end

end

@testset "mosaic" begin
    reg1 = Raster([0.1 0.2; 0.3 0.4], (X(2.0:1.0:3.0), Y(5.0:1.0:6.0)))
    reg2 = Raster([1.1 1.2; 1.3 1.4], (X(3.0:1.0:4.0), Y(6.0:1.0:7.0)))
    irreg1 = Raster([0.1 0.2; 0.3 0.4], (X([2.0, 3.0]), Y([5.0, 6.0])))
    irreg2 = Raster([1.1 1.2; 1.3 1.4], (X([3.0, 4.0]), Y([6.0, 7.0])))

    span_x1 = Explicit(vcat((1.5:1.0:2.5)', (2.5:1.0:3.5)'))
    span_x2 = Explicit(vcat((2.5:1.0:3.5)', (3.5:1.0:4.5)'))
    exp1 = Raster([0.1 0.2; 0.3 0.4], (X(Sampled([2.0, 3.0]; span=span_x1)), Y([5.0, 6.0])))
    exp2 = Raster([1.1 1.2; 1.3 1.4], (X(Sampled([3.0, 4.0]; span=span_x2)), Y([6.0, 7.0])))
    @test val(span(mosaic(first, exp1, exp2), X)) == [1.5 2.5 3.5; 2.5 3.5 4.5]
    @test all(mosaic(first, [reg1, reg2]) .===
              mosaic(first, irreg1, irreg2) .===
              mosaic(first, (irreg1, irreg2)) .===
              [0.1 0.2 missing;
               0.3 0.4 1.2;
               missing 1.3 1.4])
    @test all(mosaic(last, reg1, reg2) .===
              mosaic(last, irreg1, irreg2) .===
              mosaic(last, exp1, exp2) .=== [0.1 0.2 missing;
                                             0.3 1.1 1.2;
                                             missing 1.3 1.4])

    # 3 dimensions
    A1 = Raster(ones(2, 2, 2), (X(2.0:-1.0:1.0), Y(5.0:1.0:6.0), Ti(DateTime(2001):Year(1):DateTime(2002))))
    A2 = Raster(zeros(2, 2, 2), (X(3.0:-1.0:2.0), Y(4.0:1.0:5.0), Ti(DateTime(2002):Year(1):DateTime(2003))))
    @test all(mosaic(mean, A1, A2) |> parent .===
              first(mosaic(mean, RasterStack(A1), RasterStack(A2))) .===
              cat([missing missing missing
                   missing 1.0     1.0
                   missing 1.0     1.0    ],
                   [0.0     0.0 missing
                   0.0     0.5     1.0
                   missing 1.0     1.0    ],
                   [0.0     0.0     missing
                   0.0     0.0     missing
                   missing missing missing], dims=3))
end

@testset "inpolygon" begin
    gi_polygon = Polygon([polygon])
    poly = gi_polygon
    for poly in (polygon, gi_polygon)
        @test inpolygon((-10.0, 20.0), poly) == true
        @test inpolygon((-19.0, 29.0), poly) == true
        @test inpolygon((-30.0, 20.0), poly) == false
        @test inpolygon([(-10.0, 20.0), (-30.0, 40.0)], poly) == [true, false]
        @test inpolygon(Point([-20.0, 50.0]), poly) == false
        @test inpolygon(LineString([[-10.0, 20.0], [-30.0, 40.0]]), poly) == [true, false]
        @test inpolygon(Polygon([[[-10.0, 20.0], [-30.0, 40.0]]]), poly) == [true, false]
    end
end

@testset "rasterize" begin
    gi_polygon = Polygon([polygon])
    rev_polygon = reverse.(polygon)
    A1 = Raster(zeros(X(-20:5; sampling=Intervals()), Y(0:30; sampling=Intervals())))
    A2 = Raster(zeros(Y(0:30; sampling=Intervals()), X(-20:5; sampling=Intervals())))
    st = RasterStack((A1, A1))

    A = A1
    poly = polygon
    ord = (X, Y)
    for A in (A1, A2), (ord, poly) in (((X, Y), polygon), ((X, Y), gi_polygon), ((Y, X), rev_polygon))
        A .= 0
        rasterize!(A, poly, vals; order=ord)
        @test sum(A) == 14 # The last value overwrites the first
        A .= 0
        rasterize!(A, poly; shape=:point, fill=1, order=ord)
        @test sum(A) == 4
        A .= 0
        rasterize!(A, poly; shape=:line, fill=1, order=ord)
        @test sum(A) == 20 + 20 + 20 + 20
        A .= 0
        rasterize!(A, poly; shape=:polygon, fill=1, boundary=:center, order=ord)
        @test sum(A) == 20 * 20
        A .= 0
        rasterize!(A, poly; shape=:polygon, fill=1, boundary=:touches, order=ord)
        @test sum(A) == 21 * 21
        A .= 0
        rasterize!(A, poly; shape=:polygon, fill=1, boundary=:inside, order=ord)
        @test sum(A) == 19 * 19
        A = Raster(zeros(X(-20:5; sampling=Intervals()), Y(0:30; sampling=Intervals())))
        R = rasterize(poly, vals; to=A, order=ord)
        @test sum(skipmissing(R)) == 14 # The last value overwrites the first
        R = rasterize(poly; fill=1, to=A, order=ord)
        @test sum(skipmissing(R)) == 20 * 20
        @test_throws ArgumentError rasterize(poly, vals; to=A, order=())
        @test_throws ArgumentError rasterize!(A, poly; shape=:notashape, fill=1, order=ord)
        @test_throws ArgumentError rasterize!(A, poly; shape=:polygon, fill=1, boundary=:notaboundary, order=ord)

        st[:layer1] .= st[:layer2] .= 0
        rasterize!(st, poly; shape=:point, fill=1, order=ord)
        @test sum(st[:layer1]) == sum(st[:layer2]) == 4
        st[:layer1] .= st[:layer2] .= 0
        ra1 = rasterize(poly, vals; to=st, order=ord)
        @test name(ra1) == :layer1
        @test sum(skipmissing(ra1)) == 14 # The last value overwrites the first
        stvals = [(1, 2), (3, 4), (5, 6), (7, 8), (9, 10)]
        st2 = rasterize(poly, stvals; to=st, order=ord)
        @test sum(skipmissing(st2[:layer1])) == 24 # The last value overwrites the first
        @test sum(skipmissing(st2[:layer2])) == 28
    end

    @testset "table" begin
        table = map(polygon, vals) do p, v
            (x=p[1], y=p[2], val1=v, val2=2.0f0v)
        end
        @test sum(skipmissing(rasterize(table; to=A, order=(X=>:x, Y=>:y), name=:val1))) == 14
        rst = rasterize(table; to=A, order=(X=>:x, Y=>:y), name=(:val1, :val2))
        @test map(sum ∘ skipmissing, rst) === (val1=14, val2=28.0f0)
        R = rasterize(table; to=A, order=(X=>:x, Y=>:y))
        @test sum(skipmissing(R[:val1])) === 14
        @test sum(skipmissing(R[:val2])) === 28.0f0
        @test keys(R) == (:val1, :val2)
        @test dims(R) == dims(A)
        @test_throws ArgumentError rasterize(table; to=A, order=(), name=:val1)
        @test_throws ArgumentError rasterize(table; to=A, order=(Z,), name=:val1)
        @test_throws ArgumentError rasterize(table; to=A, name=:val1)
    end

end

@testset "resample" begin
    raster_path = maybedownload("https://download.osgeo.org/geotiff/samples/gdal_eg/cea.tif")

    output_res = 0.0027
    output_crs = EPSG(4326)
    resample_method = "near"

    ## Resample cea.tif manually with ArchGDAL
    wkt = convert(String, convert(WellKnownText, output_crs))
    AG_output = ArchGDAL.read(raster_path) do dataset
        ArchGDAL.gdalwarp([dataset], ["-t_srs", "$(wkt)",
                                "-tr", "$(output_res)", "$(output_res)",
                                "-r", "$(resample_method)"]) do warped
            ArchGDAL.read(ArchGDAL.getband(warped, 1))
        end
    end

    ## Resample cea.tif using resample
    cea = read(geoarray(raster_path))
    raster_output = resample(cea, output_res; crs=output_crs, method=resample_method)
    disk_output = resample(cea, output_res; crs=output_crs, method=resample_method, filename="resample.tif")

    cea_permuted = permutedims(read(geoarray(raster_path)), (Y, X, Band))
    permuted_output = resample(cea_permuted, output_res; crs=output_crs, method=resample_method)

    # Compare ArchGDAL, resample and permuted resample
    @test AG_output ==
        raster_output[Band(1)] ==
        disk_output[Band(1)] ==
        permutedims(permuted_output, (X, Y, Band))[Band(1)]
    @test abs(step(dims(raster_output, Y))) ≈
        abs(step(dims(raster_output, X))) ≈
        abs(step(dims(disk_output, X))) ≈
        abs(step(dims(permuted_output, X))) ≈ output_res

    rm("resample.tif")

    @testset "snapped size and dim index match" begin
        snaptarget = raster_output
        snapped = resample(cea; to=snaptarget)
        disk_snapped = resample(cea; to=snaptarget, filename="raster.tif")
        @test size(snapped) == size(disk_snapped) == size(snaptarget)
        @test isapprox(index(snaptarget, Y), index(snapped, Y))
        @test isapprox(index(snaptarget, X), index(snapped, X))
        @test isapprox(index(snaptarget, Y), index(disk_snapped, Y))
        @test isapprox(index(snaptarget, X), index(disk_snapped, X))
    end
end


text:
@testset "basic" begin
    @testset "GNNChain" begin
        n, din, d, dout = 10, 3, 4, 2

        g = GNNGraph(random_regular_graph(n, 4),
                    graph_type=GRAPH_T,
                    ndata= randn(Float32, din, n))

        gnn = GNNChain(GCNConv(din => d),
                       BatchNorm(d),
                       x -> tanh.(x),
                       GraphConv(d => d, tanh),
                       Dropout(0.5),
                       Dense(d, dout))

        testmode!(gnn)

        test_layer(gnn, g, rtol=1e-5, exclude_grad_fields=[:μ, :σ²])

        @testset "Parallel" begin
            AddResidual(l) = Parallel(+, identity, l)

            gnn = GNNChain(ResGatedGraphConv(din => d, tanh),
                           BatchNorm(d),
                           AddResidual(ResGatedGraphConv(d => d, tanh)),
                           BatchNorm(d),
                           Dense(d, dout))

            testmode!(gnn)

            test_layer(gnn, g, rtol=1e-5, exclude_grad_fields=[:μ, :σ²])
        end

        @testset "Only graph input" begin
            nin, nout = 2, 4
            ndata = rand(nin, 3)
            edata = rand(nin, 3)
            g = GNNGraph([1,1,2], [2, 3, 3], ndata=ndata, edata=edata)
            m = NNConv(nin => nout, Dense(2, nin*nout, tanh))
            chain = GNNChain(m)
            y = m(g, g.ndata.x, g.edata.e)
            @test m(g).ndata.x == y
            @test chain(g).ndata.x == y
        end
    end

    @testset "WithGraph" begin
        x = rand(Float32, 2, 3)
        g = GNNGraph([1,2,3], [2,3,1], ndata=x)
        model = SAGEConv(2 => 3)
        wg = WithGraph(model, g)
        # No need to feed the graph to `wg`
        @test wg(x) == model(g, x)
        @test Flux.params(wg) == Flux.params(model)
        g2 = GNNGraph([1,1,2,3], [2,4,1,1])
        x2 = rand(Float32, 2, 4)
        # WithGraph will ignore the internal graph if fed with a new one.
        @test wg(g2, x2) == model(g2, x2)

        wg = WithGraph(model, g, traingraph=false)
        @test length(Flux.params(wg)) == length(Flux.params(model))

        wg = WithGraph(model, g, traingraph=true)
        @test length(Flux.params(wg)) == length(Flux.params(model)) + length(Flux.params(g))
    end
end


text:
"A wrapper around taxonkit to get matching data dir"
function taxid_to_species(taxids::Array{String,1}, db_dir::String)
    db = db_dir
    ids = taxids
    outdata=Array{String,1}()

    bindir = Conda.bin_dir(:vsearch)

    tmpdir = mktempdir("/dev/shm"; prefix="taxid_", cleanup=true)
    in_taxid_lis = "$tmpdir/taxidlist.txt"
    lineage_out = "$tmpdir/lineage.txt"
    reformat_out = "$tmpdir/reformat.txt"
    open(in_taxid_lis,"w") do idsf
        for i in ids
            println(idsf, i)
        end
    end

    lineage = run(`$bindir/taxonkit lineage --data-dir $db -o $lineage_out  $in_taxid_lis`)
    reformat = read(`$bindir/taxonkit reformat --miss-rank-repl "NA" --data-dir $db $lineage_out`, String)
    for l in split(reformat, "\n")
        if l != ""
            species = split(l, ";")[end]
            push!(outdata, species)
        end
    end
    return(outdata)
end

"A wrapper around taxonkit to get matching data dir"
function taxid_to_genus(taxids::Array{String,1}, db_dir::String)
    db = db_dir
    ids = taxids
    outdata=Array{String,1}()

    bindir = Conda.bin_dir(:vsearch)

    tmpdir = mktempdir("/dev/shm"; prefix="taxid_", cleanup=true)
    in_taxid_lis = "$tmpdir/taxidlist.txt"
    lineage_out = "$tmpdir/lineage.txt"
    reformat_out = "$tmpdir/reformat.txt"
    open(in_taxid_lis,"w") do idsf
        for i in ids
            println(idsf,i )
        end
    end

    lineage = run(`$bindir/taxonkit lineage --data-dir $db -o $lineage_out  $in_taxid_lis`)
    reformat = read(`$bindir/taxonkit reformat --miss-rank-repl "NA" --data-dir $db $lineage_out`, String)
    for l in split(reformat, "\n")
        if l != ""
            genus = split(l, ";")[end-1]
            push!(outdata, genus)
        end
    end
    return(outdata)
end


text:
using InfiniteLinearAlgebra, LinearAlgebra, BandedMatrices, InfiniteArrays, MatrixFactorizations, LazyArrays, FillArrays, SpecialFunctions, Test
import LazyArrays: colsupport, rowsupport, MemoryLayout, DenseColumnMajor, TriangularLayout, resizedata!
import LazyBandedMatrices: BroadcastBandedLayout
import BandedMatrices: _BandedMatrix, _banded_qr!, BandedColumns
import InfiniteLinearAlgebra: partialqr!, AdaptiveQRData, AdaptiveLayout

@testset "Adaptive QR" begin
    @testset "test partialqr!" begin
        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1)
        l,u = bandwidths(A)
        F = AdaptiveQRData(A);
        @test bandwidths(F.data) == (1,2)
        n = 3
        partialqr!(F,n);
        F̃ = qrunblocked(Matrix(A[1:n+l,1:n]))
        @test F̃.factors ≈ F.data[1:n+l,1:n]
        @test F̃.τ  ≈ F.τ
        @test triu!(F.data[1:n+1,1:n+2]) ≈ F̃.Q'*A[1:n+1,1:n+2] # test extra columns have been modified
        n  = 6
        partialqr!(F,n);
        F̃ = qrunblocked(Matrix(A[1:n+l,1:n]))
        @test F̃.factors ≈ F.data[1:n+l,1:n]
        @test F̃.τ  ≈ F.τ
        @test triu!(F.data[1:n+1,1:n+2]) ≈ F̃.Q'*A[1:n+1,1:n+2] # test extra columns have been modified
    end

    @testset "AdaptiveQRFactors" begin
        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1)
        F = qr(A);
        @test F.factors[1,1] ≈ -sqrt(2)
        @test F.factors[100,100] ≈ qrunblocked(A[1:101,1:100]).factors[100,100]
        @test F.τ[1] ≈ 1+sqrt(2)/2
        @test F.τ[100] ≈ qrunblocked(A[1:101,1:100]).τ[100]
        Q,R = F;
        @test bandwidths(R) == (0,2)
        @test size(R) == (∞,∞)
        @test R[1,1] == -sqrt(2)
    end

    @testset "col/rowsupport" begin
        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1)
        F = qr(A);
        @test MemoryLayout(typeof(F.factors)) isa AdaptiveLayout{BandedColumns{DenseColumnMajor}}
        @test bandwidths(F.factors) == (1,2)
        @test colsupport(F.factors,1) ==  1:2
        @test colsupport(F.factors,5) ==  3:6
        @test rowsupport(F.factors,1) ==  1:3
        @test rowsupport(F.factors,5) ==  4:7
        Q,R = F;
        @test MemoryLayout(typeof(R)) isa TriangularLayout
        @test colsupport(R,1) == 1:1
        @test colsupport(R,5) == 3:5
        @test rowsupport(R,1) == 1:3
        @test rowsupport(R,5) == 5:7
        @test colsupport(Q,1) == 1:2
        @test colsupport(Q,5) == 3:6
        @test rowsupport(Q,1) == 1:3
        @test rowsupport(Q,5) == 4:7
    end

    @testset "Qmul" begin
        A = _BandedMatrix(Vcat(Ones(1,∞), (1:∞)', Ones(1,∞)), ∞, 1, 1)
        Q,R = qr(A);
        b = Vcat([1.,2,3],Zeros(∞))
        @test lmul!(Q, Base.copymutable(b)).datasize[1] == 4
        @test lmul!(Q, Base.copymutable(b)).data[1:4] ≈ qr(A[1:4,1:3]).Q*[1,2,3]

        @test Q[1,1] ≈ -1/sqrt(2)
        @test Q[200_000,200_000] ≈ -1.0
        @test Q[1:101,1:100] ≈ qr(A[1:101,1:100]).Q[:,1:100]

        @test lmul!(Q, Base.copymutable(b)) == Q*b

        r = lmul!(Q', Base.copymutable(b))
        nr = length(r.data)
        @test qr(A[1:nr+1,1:nr]).Q'b[1:nr+1] ≈ r[1:nr+1]

        @test Q'*b == r

        @test factorize(A) isa typeof(qr(A))
        @test qr(A)\b == A\b
        @test (A*(A\b))[1:100] ≈ [1:3; Zeros(97)]
    end

    @testset "Bessel J" begin
        z = 1000; # the bigger z the longer before we see convergence
        A = BandedMatrix(0 => -2*(0:∞)/z, 1 => Ones(∞), -1 => Ones(∞))
        b = Vcat([besselj(1,z)], Zeros(∞))
        F = qr(A)
        @test qr(A[1:3000,1:3000]).Q'b[1:3000] ≈ (F.Q'b)[1:3000]
        @time J = A \ Vcat([besselj(1,z)], Zeros(∞))
        @test J[1:2000] ≈ [besselj(k,z) for k=0:1999]

        z = 10_000; # the bigger z the longer before we see convergence
        A = BandedMatrix(0 => -2*(0:∞)/z, 1 => Ones(∞), -1 => Ones(∞))
        @time J = A \ Vcat([besselj(1,z)], Zeros(∞))
        @test J[1:20_000] ≈ [besselj(k,z) for k=0:20_000-1]
    end

    @testset "5-band" begin
        A = BandedMatrix(-2 => Ones(∞), -1 => Vcat(1, Zeros(∞)), 0 => Vcat([1,2,3],Zeros(∞)).+3, 1 => Vcat(1, Zeros(∞)), 2 => Ones(∞))
        b = Vcat([3,4,5],Zeros(∞))
        x = qr(A) \ b
        @test x[1:2000] ≈ (A[1:2000,1:2000]\b[1:2000])
    end

    @testset "broadcast" begin
        A = BandedMatrix(0 => -2*(0:∞), 1 => Ones(∞), -1 => Ones(∞))
        B = BandedMatrix(-2 => Ones(∞), -1 => Vcat(1, Zeros(∞)), 0 => Vcat([1,2,3],Zeros(∞)).+3, 1 => Vcat(1, Zeros(∞)), 2 => Ones(∞))

        AB = BroadcastArray(+,A,B)
        @test MemoryLayout(typeof(AB)) isa BroadcastBandedLayout{typeof(+)}

        C = cache(AB);
        resizedata!(C,103,100); resizedata!(C,203,200);
        @test C[103,104] ≈ 1.0
        C = qr(AB).factors.data.data;
        resizedata!(C,102,104); resizedata!(C,202,204);
        @test C[202,204] == AB[202,204]
        F = qr(AB);
        partialqr!(F.factors.data, 100); partialqr!(F.factors.data, 200);
        @test norm(F.factors.data.data.data[Base.OneTo.(F.factors.data.data.datasize)...]) ≤ 4000
        b = Vcat([3,4,5],Zeros(∞))
        @time x = qr(AB) \ b;
        @test x[1:300] ≈ AB[1:300,1:300] \ b[1:300]
    end
end