text:
using Base.Test

using KyotoCabinet
using KyotoCabinet.c

KyotoCabinet.pack(v::ASCIIString) = convert(Array{Uint8,1}, v)
KyotoCabinet.unpack(T::Type{ASCIIString}, buf::Array{Uint8,1}) = bytestring(buf)

# TODO: add exception type to @test_throws
# See: https://github.com/JuliaLang/julia/commit/6fa50c4183358047c772a508e7a1a44a47c94a95

# TODO
# File creation failures.
# - Invalid file name
# - File name clash with existing dir

# Get/Set failures
# - Read non existing keys
# - Set value when opened in R/O mode

function test_open()
  test_with(empty_db) do db
  end
end

function test_get_set()
  test_with(abc_db) do db
    @assert "1" == get(db, "a")
    @assert "2" == get(db, "b")
    @assert "3" == get(db, "c")
  end
end

function test_length()
  test_with(empty_db) do db
    @assert 0 == length(db)
  end
  test_with(abc_db) do db
    @assert 3 == length(db)
  end
end

# Use cursor to iterate records
function test_iterate_empty()
  test_with(empty_db) do db
    log = string()
    for (k, v) = db
      log = log * " $k:$v"
    end
    @assert "" == log
  end
end

# Test for loop over records
function test_iterate()
  test_with(abc_db) do db
    log = string()
    for (k, v) = db
      log = log * " $k:$v"
    end
    @assert " a:1 b:2 c:3" == log
  end
end

# next should fail on empty db
function test_iterate_nexts_empty()
  test_with(empty_db) do db
    s0 = start(db)
    @assert done(db, s0)
    @test_throws next(db, s0)
  end
end

# Test next->next without done method
# All movements should be in next() method.
function test_iterate_nexts()
  test_with(abc_db) do db
    s0 = start(db)
    @assert !done(db, s0)
    (kv0, s1) = next(db, s0)
    (kv1, s2) = next(db, s1)
    (kv2, s3) = next(db, s2)
    @assert done(db, s3)

    @assert ("a", "1") == kv0
    @assert ("b", "2") == kv1
    @assert ("c", "3") == kv2
  end
end

# Test generator syntax over records
function test_generator()
  test_with(abc_db) do db
    log = join(["$k:$v" for (k, v) = db], " ")
    @assert "a:1 b:2 c:3" == log
  end
end

function test_iterate_next()
  test_with(abc_db) do db
    s0 = start(db)
    @assert !done(db, s0)

    rec, s1 = next(db, s0)
    @assert ("a", "1") == rec
    @assert !done(db, s1)

    rec, s2 = next(db, s1)
    @assert ("b", "2") == rec
    @assert !done(db, s2)

    rec, s3 = next(db, s2)
    @assert ("c", "3") == rec
    @assert done(db, s3)

    @test_throws next(db, s3)
  end
end

# Test keys iterator
function test_keys_empty()
  test_with(empty_db) do db
    log = string()
    for k = keys(db)
      log = log * " $k"
    end
    @assert "" == log
  end
end

function test_keys()
  test_with(abc_db) do db
    log = string()
    for k = keys(db)
      log = log * " $k"
    end
    @assert " a b c" == log
  end
end

# Test values iterator
function test_values_empty()
  test_with(empty_db) do db
    log = string()
    for k = values(db)
      log = log * " $k"
    end
    @assert "" == log
  end
end

function test_values()
  test_with(abc_db) do db
    log = string()
    for k = values(db)
      log = log * " $k"
    end
    @assert " 1 2 3" == log
  end
end

# Test for loop over records
function test_iterate()
  test_with(abc_db) do db
    log = string()
    for (k, v) = db
      log = log * " $k:$v"
    end
    @assert " a:1 b:2 c:3" == log
  end
end

function test_get_set_failures()
  test_with(abc_db) do db
    @test_throws get(db, "z")
  end
end

function test_dict_haskey()
  test_with(abc_db) do db
    @assert !haskey(db, "")

    @assert haskey(db, "a")
    @assert haskey(db, "b")
    @assert !haskey(db, "z")
  end
end

function test_dict_get()
  test_with(abc_db) do db
    @assert "1" == get(db, "a", "0")
    @assert "2" == get(db, "b", "0")
    @assert "0" == get(db, "z", "0")

    f() = "0"
    @assert "1" == get(f, db, "a")
    @assert "2" == get(f, db, "b")
    @assert "0" == get(f, db, "z")
  end
end

function test_dict_get!()
  test_with(abc_db) do db
    @assert "1" == get!(db, "a", "0")
    @assert "2" == get!(db, "b", "0")
    @assert "0" == get!(db, "z", "0")
    @assert "0" == get(db, "z", "z")

    f() = "z"
    @assert "1" == get!(f, db, "a")
    @assert "2" == get!(f, db, "b")
    @assert "z" == get!(f, db, "zz")
    @assert "z" == get(db, "zz", "z")
  end
end

function test_dict_modify()
  test_with(abc_db) do db
    @assert 3 == length(db)
    empty!(db)
    @assert isempty(db)
  end

  test_with(abc_db) do db
    @assert "a" == getkey(db, "a", "0")
    @assert "b" == getkey(db, "b", "0")
    @assert "0" == getkey(db, "z", "0")

    @assert haskey(db, "a")
    delete!(db, "a")
    @assert !haskey(db, "a")

    @assert haskey(db, "b")
    @assert "2" == pop!(db, "b", "0")
    @assert !haskey(db, "b")

    @assert "0" == pop!(db, "z", "0")

    @test_throws pop!(db, "z")
  end
end

function test_associative()
  test_with(abc_db) do db
    @assert "1" == db["a"]
    @assert "2" == db["b"]
    @test_throws db["z"]

    @test "0" == (db["z"] = "0")
    @assert "0" == db["z"]
  end
end

# Generic Associative methods should work with Db as well
function test_associative_merge()
  test_with(empty_db) do db
    @assert !haskey(db, "a")
    @assert !haskey(db, "b")
    merge!(db, ["a"=>"a1", "b"=>"b1"])
    @assert "a1" == db["a"]
    @assert "b1" == db["b"]
  end
end

function test_path()
  file = tempname() * ".kch"
  open(Db{Bytes, Bytes}(), file, KCOWRITER | KCOCREATE) do db
    @assert file == path(db)
  end
end

function test_cas()
  test_with(abc_db) do db
    @assert cas(db, "a", "1", "1a")
    @assert "1a" == get(db, "a")

    @assert !cas(db, "a", "1", "1b")
    @assert "1a" == get(db, "a")

    @assert !cas(db, "z", "0", "0z")

    @assert cas(db, "z", nothing, "0")
    @assert "0" == get(db, "z")
    @assert !cas(db, "z", nothing, "0z")

    @assert cas(db, "z", "0", nothing)
    @assert !haskey(db, "z")
    @assert !cas(db, "z", "0", nothing)
  end
end

function test_bulkset()
  test_with(empty_db) do db
    @assert !haskey(db, "a")
    @assert !haskey(db, "b")
    @assert !haskey(db, "c")

    @assert 3 == bulkset!(db, ["a"=>"a1", "b"=>"b1", "c"=>"c1"], true)

    @assert "a1" == db["a"]
    @assert "b1" == db["b"]
    @assert "c1" == db["c"]

    @assert 2 == bulkset!(db, ["a"=>"a2", "b"=>"b2"], false)

    @assert "a2" == db["a"]
    @assert "b2" == db["b"]
    @assert "c1" == db["c"]
  end
end

function test_bulkdelete()
  test_with(abc_db) do db
    @assert haskey(db, "a")
    @assert haskey(db, "b")
    @assert haskey(db, "c")

    @assert 2 == bulkdelete!(db, ["a", "b"], true)

    @assert !haskey(db, "a")
    @assert !haskey(db, "b")
    @assert haskey(db, "c")

    @assert 1 == bulkdelete!(db, ["a", "c"], true)
    @assert 0 == bulkdelete!(db, ["b", "c"], false)
  end
end

function test_set_get_long_string()
  open(Db{ASCIIString, Bytes}(), tempname() * ".kch", KCOWRITER | KCOCREATE) do db
    bytes = "08 03 22 96 01" * repeat(" 61", 150)
    s = map(s->parseint(Uint8, s, 16), split(bytes, " "))

    db["1"] = s

    @assert s == db["1"]
  end
end

function empty_db(db::Db{ASCIIString,ASCIIString})
end

function abc_db(db::Db{ASCIIString,ASCIIString})
  set!(db, "a", "1")
  set!(db, "b", "2")
  set!(db, "c", "3")
end

function test_with(check::Function, configure::Function)
  file = tempname() * ".kch"
  open(Db{ASCIIString, ASCIIString}(), file, KCOWRITER | KCOCREATE) do db
    configure(db)
    check(db)
  end
end

test_open()
test_length()
test_path()
test_get_set()

test_iterate_empty()
test_iterate()
test_iterate_nexts_empty()
test_iterate_nexts()
test_iterate_next()

test_keys_empty()
test_keys()

test_values_empty()
test_values()

test_generator()
test_get_set_failures()

test_dict_haskey()
test_dict_get()
test_dict_get!()
test_dict_modify()

test_associative()
test_associative_merge()

test_cas()
test_bulkset()
test_bulkdelete()

test_set_get_long_string()


text:
@scope GLOBAL_CONTEXT let
    @sym [Wild{Any}] x y

    @! set x == y         --> y == x
    @! set x == x         --> true
    @! set zero(zero(x))  --> zero(x)
    @! set one(one(x))    --> one(x)
end

@scope GLOBAL_CONTEXT let
    @sym [Wild{Number}] x y [Wild{AbstractMatrix}] X Y

    @! set (+)(x) --> x
    @! set (*)(x) --> x

    @! set x + y --> y + x
    @! set x*y   --> y*x

    @! set commuteswith(+, x, y) --> true
    @! set commuteswith(+, X, y) --> true
    @! set commuteswith(+, x, Y) --> true
    @! set commuteswith(+, X, Y) --> true

    @! set commuteswith(*, x, y) --> true
    @! set commuteswith(*, X, y) --> true
    @! set commuteswith(*, x, Y) --> true
end

@scope GLOBAL_CONTEXT let
    @sym [Wild{Bool}] x y

    @! set (|)(x) --> x
    @! set (&)(x) --> x

    @! set commuteswith(&, x, y) --> true
    @! set commuteswith(|, x, y) --> true
end

ext: jl


text:
"""
    SeisPlotAmplitude(d,fmax,dt ; <keyword arguments>)

Plot amplitude-frequency 2D seismic data `d`

# Arguments
- `d::Array{Real,2}`: 2D data to plot.
- `fmax::Real`: maximum frequency.
- `dt::Real`: time sample interval.

# Keyword arguments

- `title=" "`: title of plot.
- `titlesize=16`: size of title.
- `xlabel=" "`: label on x-axis.
- `xunits=" "`: units of y-axis.
- `ylabel=" "`: label on x-axis.
- `yunits=" "`: units of y-axis.
- `labelsize=14`: size of labels on axis.
- `ox=0`: first point of x-axis.
- `dx=1`: increment of x-axis.
- `oy=0`: first point of y-axis.
- `xticks="NULL"`: ticks on x-axis.
- `yticks="NULL"`: ticks on y-axis.
- `xticklabels="NULL"`: labels on ticks of x-axis.
- `yticklabels="NULL"`: labels on ticks of y-axis.
- `ticksize=11`: size of labels on ticks.
- `fignum="NULL"`: number of figure.
- `wbox=6`: width of figure in inches.
- `hbox=6`: height of figure in inches.
- `dpi=100`: dots-per-inch of figure.
- `name="NULL"`: name of the figure to save (only if `name` is given).

# Example
```julia
julia> d = SeisLinearEvents(); SeisPlotAmplitude(d,100,0.004);
```

Credits: Aaron Stanton, 2015
"""
function SeisPlotAmplitude(d::Array{T,2}, fmax::Real, dt::Real;
                           title=" ", titlesize=16, xlabel=" ", xunits=" ",
                           ylabel=" ", yunits=" ", labelsize=14, ox=0, dx=1,
                           oy=0, xticks="NULL", yticks="NULL",
                           xticklabels="NULL", yticklabels="NULL", ticksize=11,
                           fignum="NULL", wbox=6, hbox=6, dpi=100, name="NULL") where {T<:Real}

pl.ion()
if (fignum == "NULL")
    fig = pl.figure(figsize=(wbox, hbox), dpi=dpi, facecolor="w",
                       edgecolor="k")
else
	fig = pl.figure(num=fignum, figsize=(wbox, hbox), dpi=dpi,
                           facecolor="w", edgecolor="k")
end

xlabel = "Frequency"
xunits = "(Hz)"
ylabel = "Ampltitude"
yunits = ""

nx = size(d, 2)
df = 1/dt/size(d, 1)
FMAX = df*size(d, 1)/2
if fmax > FMAX
    fmax = FMAX
end
nf = convert(Int32, floor((size(d, 1)/2)*fmax/FMAX))
y = fftshift(sum(abs.(fft(d, 1)), dims=2))/nx
y = y[round(Int,end/2):round(Int, end/2)+nf]
norm = maximum(y[:])
if (norm > 0.)
	y = y/norm
end

x = collect(0:df:fmax)
im = pl.plot(x, y)

pl.title(title)
pl.xlabel(join([xlabel " " xunits]))
pl.ylabel(join([ylabel " " yunits]))
pl.axis([0, fmax, 0, 1.1])
pl.title(title, fontsize=titlesize)
pl.xlabel(join([xlabel " " xunits]), fontsize=labelsize)
pl.ylabel(join([ylabel " " yunits]), fontsize=labelsize)
xticks == "NULL" ? nothing : pl.xticks(xticks)
yticks == "NULL" ? nothing : pl.yticks(yticks)
ax = pl.gca()
xticklabels == "NULL" ? nothing : ax.set_xticklabels(xticklabels)
yticklabels == "NULL" ? nothing : ax.set_yticklabels(yticklabels)
pl.setp(ax.get_xticklabels(), fontsize=ticksize)
pl.setp(ax.get_yticklabels(), fontsize=ticksize)

if (name == "NULL")
    pl.show()
else
    pl.savefig(name, dpi=dpi)
    pl.close()
end

return im

end


text:
function work(x)
    x * 2
end

work_single_line(x) = x+5
passed = work_single_line

println(work(4))
println(work_single_line(20))
println(passed(25))

max_stars_repo_name: tarasivashchuk/Cookbooks

max_issues_repo_name: tarasivashchuk/cookbooks

max_forks_repo_name: tarasivashchuk/cookbooks


text:
module PartialFunctions

export $
export <|
include("reversedfunctions.jl")

name = (string ∘ Symbol)

struct PartialFunction{F<:Function, T<:Tuple, N<:NamedTuple} <: Function
    func::F
    args::T
    kwargs::N
    PartialFunction(f, args::T, kwargs::N) where {T<:Tuple, N<:NamedTuple} = let
        new{typeof(f), T, N}(f, args, kwargs)
    end
end

function PartialFunction(f::PartialFunction, newargs::Tuple, newkwargs::NamedTuple)
    allargs = (f.args..., newargs...)
    allkwargs = (;f.kwargs..., newkwargs...)
    PartialFunction(f.func, allargs, allkwargs)
end

(p::PartialFunction)(newargs...; newkwargs...) = p.func(p.args..., newargs...; p.kwargs..., newkwargs...)

"""
    (\$)(f::Function, args...)
Partially apply the given arguments to f. Typically used as infix `f \$ args`

The returned function is of type [`PartialFunctions.PartialFunction{typeof(f), typeof(args)}`](@ref)

# Examples

```jldoctest
julia> using PartialFunctions

julia> simonsays = println \$ "Simon says: "
println("Simon says: ", ...)

julia> simonsays("Partial function application is cool!")
Simon says: Partial function application is cool!
```
"""
($)(f::Function, args::Tuple) = PartialFunction(f, args, (;))
($)(f::Function, arg) = PartialFunction(f, (arg,), (;))
($)(f::Function, kwargs::NamedTuple) = PartialFunction(f, (), kwargs)
($)(f::Function, n::Tuple{<:Tuple, <:NamedTuple}) = let (args, kwargs) = n
    PartialFunction(f, args, kwargs)
end
($)(f::Function, n::Tuple{<:Any, <:NamedTuple}) = let (arg, kwargs) = n
    PartialFunction(f, (arg,), kwargs)
end

($)(f::DataType, args) = ($)(identity∘f, args)
"""
    <|(f, args)

Applies a function to the succeeding argument or tuple of arguments. Acts as the reverse
of [`|>`](@ref), and is especially useful when combined with partial functions for
an alternative, low-parenthese function chaining syntax

# Examples
```@jldoctest
julia> using PartialFunctions

julia> isdigit <| '1'
true

julia> (+) <| (2, 3)...
5

julia> map \$ Int <| [1.0, 2.0, 3.0]
3-element Array{Int64,1}:
 1
 2
 3
```
"""
(<|)(f::Function, args...) = f(args...)
(<|)(f::Function, ::Tuple{}) = f()

function Base.Symbol(pf::PartialFunction)
    func_name = name(pf.func)
    argstring = "(" * join(repr.(pf.args), ", ") *", ..."
    if isempty(pf.kwargs)
        argstring *= ")"
    else
        argstring *= "; " * strip(repr(pf.kwargs), ['(', ')']) * " ...)"
    end
    Base.Symbol("$(func_name)$(argstring)")
end

Base.show(io::IO, pf::PartialFunction) = print(io, name(pf))
Base.show(io::IO, ::MIME"text/plain", pf::PartialFunction) = show(io, pf)

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule FFTW_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("FFTW")
JLLWrappers.@generate_main_file("FFTW", UUID("f5851436-0d7a-5f13-b9de-f02708fd171a"))
end  # module FFTW_jll


text:
module Functions

using ..TopOptProblems, ..FEA, ..CheqFilters
using ..Utilities, ..GPUUtils, CuArrays
using ForwardDiff, LinearAlgebra, GPUArrays
using Parameters: @unpack
import CUDAdrv
using TimerOutputs, CUDAnative, JuAFEM
using StatsFuns, MappedArrays, LazyArrays

export  Objective,
        Constraint,
        VolumeFunction,
        ComplianceFunction,
        AbstractFunction,
        getfevals,
        getmaxfevals,
        maxedfevals,
        getnvars,
        GlobalStress

const to = TimerOutput()
const dev = CUDAdrv.device()
const ctx = CUDAdrv.CuContext(dev)

abstract type AbstractFunction{T} <: Function end

@params struct Objective <: Function
    f
end
@inline function Base.getproperty(o::Objective, s::Symbol)
    s === :f && return getfield(o, :f)
    return getproperty(o.f, s)
end
@inline function Base.setproperty!(o::Objective, s::Symbol, v)
    s === :f && return setfield!(o, :f, v)
    return setproperty!(o.f, s, v)
end

@params struct Constraint <: Function
    f
    s
end
@inline function Base.getproperty(c::Constraint, s::Symbol)
    s === :f && return getfield(c, :f)
    s === :s && return getfield(c, :s)
    return getproperty(c.f, s)
end
@inline function Base.setproperty!(c::Constraint, s::Symbol, v)
    s === :f && return setfield!(c, :f, v)
    s === :s && return setfield!(c, :s, v)
    s === :reuse && return setproperty!(c.f, :reuse, v)
    return setfield!(c.f, s, v)
end

Base.broadcastable(o::Union{Objective, Constraint}) = Ref(o)
getfunction(o::Union{Objective, Constraint}) = o.f
getfunction(f::AbstractFunction) = f
GPUUtils.whichdevice(o::Union{Objective, Constraint}) = o |> getfunction |> whichdevice
Utilities.getsolver(o::Union{Objective, Constraint}) = o |> getfunction |> getsolver
Utilities.getpenalty(o::Union{Objective, Constraint}) = o |> getfunction |> getsolver |> getpenalty
Utilities.setpenalty!(o::Union{Objective, Constraint}, p) = setpenalty!(getsolver(getfunction(o)), p)
Utilities.getprevpenalty(o::Union{Objective, Constraint}) = o |> getfunction |> getsolver |> getprevpenalty

@define_cu(Objective, :f)
(o::Objective)(x, grad) = o.f(x, grad)

@define_cu(Constraint, :f)
(c::Constraint)(x, grad) = c.f(x, grad) - c.s

getfevals(o::Union{Constraint, Objective}) = o |> getfunction |> getfevals
getfevals(f::AbstractFunction) = f.fevals
getmaxfevals(o::Union{Constraint, Objective}) = o |> getfunction |> getmaxfevals
getmaxfevals(f::AbstractFunction) = f.maxfevals
maxedfevals(o::Union{Objective, Constraint}) = maxedfevals(o.f)
maxedfevals(f::AbstractFunction) = getfevals(f) >= getmaxfevals(f)
getnvars(o::Union{Constraint, Objective}) = o |> getfunction |> getnvars
getnvars(f::AbstractFunction) = length(f.solver.vars)

include("compliance.jl")
include("volume.jl")
include("stress.jl")

end


text:
using QuasiMonteCarlo, Distributions, Test
using Test

#1D
lb = 0.0
ub = 5.0
n = 5
d = 1
QuasiMonteCarlo.sample(n,lb,ub,GridSample(0.1))
QuasiMonteCarlo.sample(n,lb,ub,UniformSample())
QuasiMonteCarlo.sample(n,lb,ub,SobolSample())
QuasiMonteCarlo.sample(n,lb,ub,LatinHypercubeSample())
QuasiMonteCarlo.sample(20,lb,ub,LowDiscrepancySample(10))
QuasiMonteCarlo.sample(5,d,Cauchy())
QuasiMonteCarlo.sample(5,d,Normal(0,4))

#ND
lb = [0.1,-0.5]
ub = [1.0,20.0]
n = 5
d = 2

#GridSample{T}
s = QuasiMonteCarlo.sample(n,lb,ub,GridSample([0.1,0.5]))
@test isa(s,Matrix{typeof(s[1][1])}) == true

#UniformSample()
s = QuasiMonteCarlo.sample(n,lb,ub,UniformSample())
@test isa(s,Matrix{typeof(s[1][1])}) == true

#SobolSample()
s = QuasiMonteCarlo.sample(n,lb,ub,SobolSample())
@test isa(s,Matrix{typeof(s[1][1])}) == true

#LHS
s = QuasiMonteCarlo.sample(n,lb,ub,LatinHypercubeSample())
@test isa(s,Matrix{typeof(s[1][1])}) == true

#LDS
s = QuasiMonteCarlo.sample(n,lb,ub,LowDiscrepancySample([10,3]))
@test isa(s,Matrix{typeof(s[1][1])}) == true

#Distribution 1
s = QuasiMonteCarlo.sample(n,d,Cauchy())
@test isa(s,Matrix{typeof(s[1][1])}) == true

#Distribution 2
s = QuasiMonteCarlo.sample(n,d,Normal(3,5))
@test isa(s,Matrix{typeof(s[1][1])}) == true


text:
#tracywidomlis.jl
#Algorithm 2.7 of Random Eigenvalues by Alan Edelman

#Experiment:    Sample random permutations
#Plot:          Histogram of lengths of longest increasing subsequences
#Theory:        The Tracy-Widom law

## Parameters
t = 10000       # number of trials
n = 6^6         # length of permutations
dx = 1/6        # bin size

## Experiment
require("patiencesort.jl")
function tracywidomlis(t, n, dx)
    ## single processor: 1x speed
    #v = [patiencesort(randperm(n)) for i = 1:t]

    ## simple parallelism: Nx speed * 130%
    v = pmap((i)->patiencesort(randperm(n)), 1:t)

    ## maximum parallelism: Nx speed
    #grouped = floor(t/(nprocs()-2))
    #println(grouped)
    #v = vcat(pmap((i)->[patiencesort(randperm(n)) for j = 1:grouped], 1:t/grouped)...)

    w = (v-2sqrt(n))/n^(1/6)
    return hist(w, -5:dx:2)
end
@time grid, count = tracywidomlis(t, n, dx)

## Plot
using Winston
p = FramedPlot()
h = Histogram(count/(t*dx), step(grid))
h.x0 = first(grid)
include("../1/tracywidom.jl")
add(p, h)
if isinteractive()
    Winston.display(p)
else
    file(p, "tracywidomlis.png")
end


text:
#=
# start of assignment 3 code
"""
    compliment(base)

Get the DNA compliment of the provided base:

    A <-> T
    G <-> C

Accepts `String` or `Char`, but always returns `Char`.
If a valid base is not provided, the function throws an error.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> compliment('A')
    'T'

    julia> compliment("G")
    'C'

    julia> compliment("T")
    'A'

    julia> compliment('C')
    'G'
"""
function compliment(base)
    base = uppercase(base)
    if base == "A" || base == 'A'
        return 'T'
    elseif base == "T" || base == 'T'
        return 'A'
    elseif base == "G" || base == 'G'
        return 'C'
    elseif base == "C" || base == 'C'
        return 'G'
    else throw(ErrorException("Invalid base provided."))
    end
end

"""
    ispurine(base)

A boolean function that returns `true` if the base is a purine (A or G)
and `false` if it is not.
The function only supports bases A, C, G, and T (throws an error for other values).
Accepts `String` or `Char`.

Examples
=========

    julia> ispurine('A')
    true

    julia> ispurine("C")
    false

    julia> if ispurine("G")
               println("It's a purine!")
           else
               println("It's a pyrimidine!")
           end
    It's a purine!

    julia> ispurine('B')
    Error: "Base B not supported")
"""
function ispurine(base)
    base = uppercase(base)
    if base == "A" || base == 'A' || base == "G" || base == 'G'
        return true
    elseif base == "T" || base == 'T' || base == "C" || base == 'C'
        return false
    else throw(ErrorException("Base is not supported."))
    end
end

"""
    ispyrimidine(base)

A boolean function that returns `true` if the base is a pyrimidine (C or T)
and `false` if it is not.
The function only supports bases A, C, G, and T (throws an error for other values).
Accepts `String` or `Char`.

Examples
=========

    julia> ispyrimidine('G')
    false

    julia> ispyrimidine("T")
    true

    julia> if ispyrimidine("G")
               println("It's a pyrimidine!")
           else
               println("It's a purine!")
           end
    It's a purine!

    julia> ispyrimidine('X')
    Error: "Base X not supported"
"""
function ispyrimidine(base)
    return !ispurine(base)
end

"""
    base_type(base)

Determines whether a base is a purine (A or G) or pyrimidine (T or C),
and returns a `String`.

Examples
≡≡≡≡≡≡≡≡≡≡

    julia> base_type("G")
    "purine"

    julia> base_type('C')
    "pyrimidine"

    julia> base_type('Z')
    Error: "Base Z not supported"

    julia> x = base_type('A'); println(x)
    purine
"""
function base_type(base)
    base = uppercase(base)
    if ispyrimidine(base)
        return "pyrimidine"
    elseif ispurine(base)
        return "purine"
    else throw(ErrorException("Base is not supported."))
    end
end
# end of assignment 3 code
=#


text:
struct LocalVar
    is_mutable  :: Ref{Bool} # mutability
    is_shared   :: Ref{Bool} # shared between different physical scopes/actual functions.
    sym         :: Symbol
end

GlobalVar = Symbol

readable_var(sym::Symbol) = LocalVar(Ref(false), Ref(false), sym)
global_var(sym::Symbol) = sym


text:
## Eigen Vectors/Values
## ====================
using Sigma

dist(A,B; epsilon = 1.0) = sum([(A[i] - B[i])^2 for i = 1:length(A)]) < epsilon

dist(A,B; epsilon = 1.0) = sum([abs(A[i] - B[i]) for i = 1:length(A)]) < epsilon
dist3(A,B; epsilon = 1.0) = (&)([abs(A[i] - B[i]) < epsilon for i = 1:length(A)]...)
dist2(A,B; epsilon = 1.0) = sum([abs(A[i] - B[i]) for i = 1:length(A)])

function issquare(x)
  size(x)[1] == size(x)[2]
end

function eigen(A, nsamples::Integer; precision=0.001, epsilon = 0.01)
  @assert issquare(A)
  len = size(A)[1]
  v = mvuniform(0.1,1.0, len)
  # V is nonzero
  λ = uniform(0, 1.0)

  sol = rand((A,v,λ), dist3(A*v, λ*v; epsilon=epsilon), nsamples; parallel=true, precision=precision, RandVarType=Sigma.SymbolicRandVar)
  return sol
end

A = mvuniform(0,1,5,5)

sol = eigen(A, 10; epsilon = 0.0001, precision=0.00001)
A = sol[1][1]
v = sol[2][1]
λ = sol[3][1]
dist2(A*v, λ*v)


text:
"""
    config!([logger], level; fmt::AbstractString, levels::Dict{AbstractString, Int}, colorized::Bool) -> Logger

Sets the `Memento._log_levels`, creates a default root logger with a `DefaultHandler`
that prints to stdout.

# Arguments
* 'logger::Union{Logger, AbstractString}`: The logger to configure (optional)
* `level::AbstractString`: the minimum logging level to log message to the root logger (required).

# Keywords
* `fmt::AbstractString`: a format string to pass to the `DefaultFormatter` which describes
    how to log messages (defaults to `Memento.DEFAULT_FMT_STRING`)
* `levels`: the default logging levels to use (defaults to `Memento._log_levels`).
* `colorized`: whether or not the message to stdout should be colorized.
* `recursive`: whether or not to recursive set the level of all child loggers.
* `substitute`: whether or not to substitute the global logger with Memento
  (only supported on julia 0.7).
* `propagate`: whether the logger should also send messages to parent loggers (default: true)

# Returns
* `Logger`: the logger passed in, or the root logger.
"""
config!(level::AbstractString; kwargs...) = config!("root", level; kwargs...)

function config!(logger::AbstractString, level::AbstractString; kwargs...)
    config!(getlogger(logger), level; kwargs...)
end

function config!(
    logger::Logger, level::AbstractString;
    fmt::AbstractString=DEFAULT_FMT_STRING, levels=_log_levels, colorized=true,
    recursive=false, substitute=false, propagate=true
)
    logger.levels = levels
    setlevel!(logger, level; recursive=recursive)
    setpropagating!(logger, propagate)
    handler = DefaultHandler(
        stdout,
        DefaultFormatter(fmt), Dict{Symbol, Any}(:is_colorized => colorized)
    )
    logger.handlers["console"] = handler
    register(logger)

    substitute && substitute!()

    return logger
end

"""
    reset!()

Removes all registered loggers and reinitializes the root logger
without any handlers.
"""
function reset!()
    empty!(_loggers)
    register(Logger("root"))
    nothing
end


text:
using Zygote: @adjoint
using Flux: Optimise

@adjoint function square_metric(x)
	square_metric(x), dy -> (square_metric_back(x, dy),)
end

@adjoint function norm2_metric(x)
	D = norm2_metric(x)
	D, dy -> (norm2_metric_back(x, D, dy),)
end

@adjoint function inner_prod_metric(x)
	inner_prod_metric(x), dy -> (inner_prod_metric_back(x, dy),)
end

@adjoint function positive(x)
	positive(x), dy -> (positive_back(x, dy),)
end

power(x, α) = x.^α
@adjoint function power(x, α)
	y = power(x, α)
	y, function (dy)
		α*(@. dy*x^(α-1)), tr(dy'* (@. y*log(x)))
	end
end

# gradient check for functions
function gradient_check(f, args...; η = 1e-5)
  g = gradient(f, args...)
  dy_expect = η*sum(abs2.(g[1]))
  dy = f(args...)-f([gi == nothing ? arg : arg.-η.*gi for (arg, gi) in zip(args, g)]...)
  @show dy
  @show dy_expect
  isapprox(dy, dy_expect, rtol=1e-3, atol=1e-8)
end

# gradient check for models
function model_gradient_check(f, ps; η = 1e-5)
    g = gradient(f, ps)
    dy_expect = η*sum([gp == nothing ? 0.0 : sum(abs2.(gp)) for gp in values(g.grads)])

    f_original = f()
    for p in ps
			g[p] == nothing ? continue : Optimise.update!(p, -η.*g[p])
    end
    f_new = f()
    dy = f_original - f_new
  @show dy
  @show dy_expect
  isapprox(dy, dy_expect, rtol=1e-2, atol=1e-8)
end

max_issues_repo_name: HamletWantToCode/GPFlux

max_forks_repo_name: HamletWantToCode/GPFlux


text:
using ArgParse
using JSON

include("lib/eval.jl")

function main(args)
    s = ArgParseSettings()
    s.description = "Evaluation script for generations"

    @add_arg_table s begin
        ("--generations"; required=true; help="generations JSON file")
        ("--meta"; action=:store_true; help="print meta data")
    end

    # parse args
    isa(args, AbstractString) && (args=split(args))
    o = parse_args(args, s; as_symbols=true)

    # read JSON file
    results = JSON.parsefile(o[:generations])

    # print META data
    if o[:meta]
        println("META")
        for (k,v) in results["meta"]
            println(k, " => ", v)
        end
        flush(STDOUT)
    end

    hyp, ref = [], []
    for caption in results["captions"]
        push!(hyp, caption["hypothesis"])
        push!(ref, caption["references"])
    end
    results = 0; gc()

    # evaluate
    ti = now()
    @printf("Evaluation started (date=%s)\n", ti)
    scores, bp, hlen, rlen = bleu(hyp, ref)
    @printf("BLEU = %.1f/%.1f/%.1f/%.1f ",
            map(i->i*100,scores)...)
    @printf("(BP=%g, ratio=%g, hyp_len=%d, ref_len=%d)\n",
            bp, hlen/rlen, hlen, rlen)
    tf = now()
    @printf("Time elapsed: %s [%s]\n", tf-ti, tf)
end

!isinteractive() && !isdefined(Core.Main, :load_only) && main(ARGS)

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
using ProgressMeter

include("fig1acegi.jl")
include("supp1acegi.jl")
include("fig1bdfhj.jl")
include("supp1bdfhj.jl")
include("fig2acegi.jl")
include("supp2acegi.jl")
include("fig2bdfhj.jl")
include("supp2bdfhj.jl")
include("fig3acegi.jl")
include("supp3acegi.jl")
include("fig3bdfhj.jl")
include("supp3bdfhj.jl")


text:
module SampleChains

using NestedTuples
using TupleVectors
import Base

include("chain.jl")
include("multichain.jl")
include("utils.jl")

end


text:
import Base: <, <=, ==, !==, >=, >,
             &, |, !, ==,
             +, -, *, /, ^

## a few definitionsn for ValidatedNumerics that don't fit in there:
## Validated numerics doesn't define these, as the order ins't a total order
Base.isless(i::Interval{T}, j::Interval{S}) where {T<:Real, S<:Real} = isless(i.hi, j.lo)
#<=(i::Interval{T}, j::Interval{S}) where {T<:Real, S<:Real} = <=(i.hi, j.lo)

#Base.max(i::Interval, j::Interval) = Interval(max(i.lo,j.lo), max(i.hi,j.hi))
#Base.min(i::Interval, j::Interval) = Interval(min(i.lo,j.lo), min(i.hi,j.hi))

## BInterval represents TRUE (true, true), FALSE (false, false) and MAYBE (false, true)
struct BInterval <: Integer
    lo :: Bool
    hi :: Bool

    function BInterval(a::Bool, b::Bool)
        if a & !b
            a,b=b,a
        end
        new(a, b)
    end
end

const TRUE = BInterval(true, true)
const FALSE = BInterval(false, false)
const MAYBE = BInterval(false, true)

Base.convert(::Type{BInterval}, y::Bool) = y ? TRUE : FALSE
Base.promote_rule(::Type{BInterval}, ::Type{Bool}) = BInterval

(&)(x::BInterval, y::BInterval) = BInterval(x.lo & y.lo, x.hi & y.hi)
(|)(x::BInterval, y::BInterval) = BInterval(x.lo | y.lo, x.hi | y.hi)
(!)(x::BInterval) = BInterval(!x.lo, !x.hi)
==(x::BInterval, y::BInterval) = (x.lo==y.lo)&&(x.hi==y.hi)

## ...
function negate_op(op)
    (op === <)   && return(>=)
    (op === <=)  && return(>)
    (op === ==)  && return(!==)
    (op === !==) && return(==)
    (op === >=)  && return(<)
    (op === >)   && return(<=)
    (op === in)  && return((x,y) -> !in(x,y))
end

## OIinterval includes interval, if defined on interval and if continuous on interval
struct OInterval <: Real
    val::Interval
    def::BInterval
    cont::BInterval
    OInterval(val, def, cont) = new(val, def, cont)
end
(O::OInterval)(o::OInterval) = o
Base.show(io::IO,  o::OInterval)  = print(io, "OInterval: ", o.val, " def=", o.def, " cont=",o.cont)

## some outer constructors...
OInterval(a::T, b::S) where {T <: Real, S <: Real}= OInterval(Interval(a,b), TRUE, TRUE)
OInterval(a::OInterval, b::OInterval) = a === a ? a : error("a is not b?") ## why is this one necessary?
OInterval(a::OInterval) = a
OInterval(a) = OInterval(a,a)   # thin one...
OInterval(i::Interval) = OInterval(i.lo, i.hi)

Base.convert(::Type{OInterval}, i::Interval) = OInterval(i.lo, i.hi)

Base.promote_rule(::Type{OInterval}, ::Type{A}) where {A<:Real} = OInterval
Base.float(i::OInterval) = OInterval(float(i.val), i.def, i.cont)

## A region is two OIntervals.
struct Region
    x::OInterval
    y::OInterval
end

ImplicitEquations.diam(x::OInterval) = diam(x.val)
Base.isfinite(x::OInterval) = isfinite(x.val)

## extend functions for OInterval
## Notice these return BIntervals -- not Bools
##  bypass isless...

<(i::OInterval, j::OInterval)  = i.val <  j.val ? TRUE : (i.val > j.val ? FALSE : MAYBE)
<=(i::OInterval, j::OInterval) = i.val <= j.val ? TRUE : (i.val > j.val ? FALSE : MAYBE)
==(i::OInterval, j::OInterval) =
    if i.val == j.val
        TRUE
    elseif (i < j) == TRUE
        FALSE
    elseif (i > j) == TRUE
        FALSE
    else
        MAYBE
    end
!==(i::OInterval, j::OInterval) = (i < j == TRUE) ? TRUE : ((i > j ==TRUE) ? TRUE : MAYBE)
>=(i::OInterval, j::OInterval) = j < i
>(i::OInterval, j::OInterval) = j <= i

"""

Screen a value using `NaN` values.
Use as  with `f(x,y) = x*y * screen(x > 0)`

Also aliased to I_(x>0)

An expression like `x::OInterval > 0` is not Boolean, but
rather a `BInterval` which allows for a "MAYBE" state. As such, a
simple ternary operator, like `x > 0 ? 1 : NaN` won't work, to screen values.

"""
screen(ex) = (ex == FALSE) ? NaN : 1
const I_ = screen               # indicator function like!

## Functions which are continuous everywhere
## +, -, *, sin, cos, ...
function +(x::OInterval, y::OInterval)
    val = x.val + y.val
    def = x.def & y.def
    cont = x.cont & y.cont
    OInterval(val, def, cont)
end

function -(x::OInterval, y::OInterval)
    val = x.val - y.val
    def = x.def & y.def
    cont = x.cont & y.cont
    OInterval(val, def, cont)
end
-(x::OInterval) = OInterval(-x.val, x.def, x.cont)

function *(x::OInterval, y::OInterval)
    val = x.val * y.val
    def = x.def & y.def
    cont = x.cont & y.cont
    OInterval(val, def, cont)
end

Base.sin(x::OInterval) = OInterval(sin(x.val), x.def, x.cont)
Base.cos(x::OInterval) = OInterval(cos(x.val), x.def, x.cont)
Base.sec(x::OInterval) = 1/cos(x)
Base.csc(x::OInterval) = 1/sin(x)
Base.tan(x::OInterval) = sin(x)/cos(x)
Base.cot(x::OInterval) = 1/tan(x)

Base.sinh(x::OInterval) = (exp(x) - exp(-x))/2
Base.cosh(x::OInterval) = (exp(x) + exp(-x))/2
Base.sech(x::OInterval) = 1/cosh(x)
Base.csch(x::OInterval) = 1/sinh(x)
Base.tanh(x::OInterval) = sinh(x)/cosh(x)
Base.coth(x::OInterval) = 1/tanh(x)

function Base.asin(x::OInterval)
    if x.val.hi < -1.0 ||  x.val.lo > 1.0
        OInterval(x.val, FALSE, FALSE)
    elseif (x.val.lo < -1.0) & (x.val.hi > 1.0)
        OInterval(Interval(-pi/2, pi/2), x.def & MAYBE, x.cont)
    elseif x.val.hi > 1
        OInterval(Interval(asin(x.val.lo), pi/2), x.def & MAYBE, x.cont)
    elseif x.val.lo < -1
        OInterval(Interval(-pi/2, asin(x.val.hi)), x.def & MAYBE, x.cont)
    else
        OInterval(asin(x.val), x.def, x.cont)
    end
end
Base.acsc(x::OInterval) = asin(1/x)

function Base.acos(x::OInterval)
    if x.val.lo > 1.0 || x.val.hi < -1.0
        OInterval(x.val, FALSE, FALSE)
    elseif (x.val.lo < -1.0) & (x.val.hi > 1.0)
        OInterval(Interval(-pi/2, pi/2), x.def & MAYBE, x.cont)
    elseif x.val.hi > 1
        OInterval(Interval(acos(x.val.lo), pi/2), x.def & MAYBE, x.cont)
    elseif x.val.lo < -1
        OInterval(Interval(-pi/2, acos(x.val.hi)), x.def & MAYBE, x.cont)
    else
        OInterval(acos(x.val), x.def, x.cont)
    end
end
Base.asec(x::OInterval) = acos(1/x)

Base.atan(x::OInterval) = OInterval(atan(x.val), x.def, x.cont)
Base.acot(x::OInterval) = OInterval(acot(x.val), x.def, x.cont)

Base.asinh(x::OInterval) = log(x + sqrt(1 - x^2))
Base.acosh(x::OInterval) = log(x + sqrt(x + 1)*sqrt(x - 1))
Base.atanh(x::OInterval) = 1/2*(log(1+x) - log(1-x))
Base.asech(x::OInterval) = acosh(1/x)
Base.acsch(x::OInterval) = asinh(1/x)
Base.acoth(x::OInterval) = atanh(1/x)

Base.exp(x::OInterval) = OInterval(exp(x.val), x.def, x.cont)

function Base.hypot(x::OInterval, y::OInterval)
    val = hypot(x.val, y.val)
    def = x

    OInterval(val, x.def & y.def, x.cont & y.cont)
end

## discontinous functions

## /
## division is slow
function /(x::OInterval, y::OInterval)
    ## 0 is the issue.
    if 0.0 ∈ y.val
        ## maybe defined, maybe continuous
        OInterval(x.val/y.val, x.def & MAYBE, x.cont & MAYBE)
    else
        OInterval(x.val/y.val, x.def & y.def, x.cont & y.cont)
    end
end

## log
function Base.log(x::OInterval)
    if x.val.hi <= 0
        OInterval(x.val, FALSE, FALSE)
    elseif x.val.lo <= 0
        OInterval(log(x.val), BInterval(false, x.def.hi), x.cont)
    else
        OInterval(log(x.val), x.def, x.cont)
    end
end
Base.log(k::Irrational{:e},x::OInterval) = log(x)
Base.log(k::Real, x::OInterval) = log(x)/log(k)

function ^(a::OInterval, x::OInterval)
    OInterval(a.val^x.val, x.def, x.cont)
end

## x^n
function ^(x::OInterval, n::Integer)
    OInterval(x.val^n, x.def, x.cont)
end

## Rational ones can be exact, whereas floating point exponents are not. The main
## example would be `x^(1/3)` and `x^(1//3)`
function ^(x::OInterval, q::Rational)
    q < 0 && return(1/x^(-q))
    ## clean up odd denominator
    if isodd(q.den)
        val = x.val^(q.num)
        val = Interval(sign(val.lo)*abs(val.lo)^(1/q.den), sign(val.hi) * abs(val.hi)^(1/q.den))
    else
        ## nothing to salvage by working over Q.
        return(x^float(q))
    end
    OInterval(val, x.def, x.cont)
end

function ^(x::OInterval, r::Real)
    r < 0 && return(1/x^(-r))
    if x.val.hi < 0
        OInterval(x.val, BInterval(false, false), x.cont)
    elseif x.val.lo < 0
        OInterval(x.val^r, BInterval(false, x.def.hi), x.cont)
    else
        OInterval(x.val^r, x.def, x.cont)
    end
end

Base.sqrt(x::OInterval) = x^(1//2)
Base.cbrt(x::OInterval) = x^(1//3)

## abs, max, min, floor, trunc, ...
Base.abs(x::OInterval) = OInterval(abs(x.val), x.def, x.cont)
function Base.floor(x::OInterval)
    val = floor(x.val)
    def = x.def
    cont = BInterval(x.def.lo & (val.lo == val.hi), x.def.hi)
    OInterval(val, def, cont)
end

function Base.ceil(x::OInterval)
    val = ceil(x.val)
    def = x.def
    cont = BInterval(x.def.lo, (val.lo == val.hi) & x.def.hi)
    OInterval(val, def, cont)
end

Base.max(x::OInterval, y::OInterval) = OInterval(max(x.val, y.val), x.def & y.def, x.cont & y.cont)
Base.min(x::OInterval, y::OInterval) = OInterval(min(x.val, y.val), x.def & y.def, x.cont & y.cont)

## others sign, mod, ...
function Base.sign(x::OInterval)
    OInterval(sign(x.val), x.def, x.cont & ((x.val.lo < 0 < x.val.hi) ? MAYBE : TRUE))
end

## pixels are [0, W) x [0, H) where (0,0) lower left, (W-1, H-1) upper right
## we assume reg is of the form [a,b) x [c,d)
## TODO: Should this depend on the operation? f(I) < a should use wider I, f(I) > a should use narrower?
function xy_region(u, L, R, B, T, W, H)
    px, py = u.x.val,  u.y.val
    a = L + px.lo * (R - L) / W
    b = L + (px.hi) * (R - L) / W
    c = B + py.lo * (T - B) / H
    d = B + (py.hi) * (T - B) / H
    delta = sqrt(eps())

    x, y  = OInterval(a+(u.x.cont==TRUE)*delta,b-delta), OInterval(c+0*delta,d-delta)
    x, y
end

function compute_fxy(p::Pred, u::Region, L, R, B, T, W, H)
    x, y = xy_region(u, L, R, B, T, W, H)
    p.f(x, y)
end

"""

Compute whether predicate holds in a given region. Returns FALSE, TRUE or MAYBE

"""
function compute(p::Pred, u::Region, L, R, B, T, W, H)
    fxy = compute_fxy(p, u, L, R, B, T, W, H)

    (fxy.def == FALSE) && return (FALSE)
    isempty(fxy.val) && return (FALSE & fxy.def)

    if p.op === ==
        return((p.val ∈ fxy.val) ? MAYBE : FALSE)
    elseif negate_op(p.op) === ==
        return((p.val ∈ fxy.val) ? MAYBE : TRUE)
    end

    out = p.op(fxy, p.val)

    ## domain tracking
    out & fxy.def
end

## build up answer
function compute(ps::Preds, u::Region, L, R, B, T, W, H)
    vals = [compute(p, u, L, R, B, T, W, H) for p in ps.ps]
    val = popfirst!(vals)
    for i in 1:length(ps.ops)
        val = ps.ops[i](val, vals[i])
    end
    val
end

"""

Does this function have a zero crossing? Heuristic check.

We return `TRUE` or `MAYBE`. However, that
leaves some functions showing too much red in the case where there is no zero.

"""
function cross_zero(r::Pred, u::Region, L, R, B, T, W, H)
    x, y = xy_region(u, L, R, B, T, W, H)
    dx, dy = diam(x), diam(y)

    n = 20                      # number of random points chosen
    λ1s, λ2s = [0.0; 1.0;rand(n)], [0.0; 1.0; rand(n)]
    β1s, β2s = [1.0; 0.0; rand(n)], [1.0; 0.0; rand(n)]
    for i in 1:(n+2)
        rx, ry = x.val.lo + λ1s[i]*dx, y.val.lo + λ2s[i] * dy
        sx, sy = x.val.lo + β1s[i]*dx, y.val.lo + β2s[i] * dy
        ll = OInterval(rx), OInterval(ry)
        ur = OInterval(sx), OInterval(sy)
        val = (r.f(ll...) - r.val) * (r.f(ur...) - r.val)
        ((val <= 0)==TRUE) && return(TRUE)
    end
    return(MAYBE)
end

"""

Does this function have a value in the pixel satisfying the inequality? Return `TRUE` or `MAYBE`.

"""
function check_inequality(r::Pred, u::Region, L, R, B, T, W, H)
    x, y = xy_region(u, L, R, B, T, W, H)
    dx, dy = diam(x), diam(y)
    # check 10 random points, here fxy.def == TRUE so we can evaluate function
    n = 10
    λ1s, λ2s = [0.0; 1.0; rand(n)], [0.0; 1.0; rand(n)]
    for i in 1:n
        rx, ry = x.val.lo + λ1s[i]*dx, y.val.lo + λ2s[i] * dy
        if r.op(r.f(rx,ry), r.val)
            return(TRUE)
        end
    end
    return(MAYBE)
end


text:
using NURBS
using Base.Test

@testset "Knot Building" begin

    @testset "Open Uniform Knots" begin
        @test typeof(knot(5, 2)) == Array{Float64,1}
        @test knot(5, 2) == [0.0, 0.0, 1.0, 2.0, 3.0, 4.0, 4.0]
        @test knot(6, 3, false, 0.25) == [0.0, 0.0, 0.0, 0.25, 0.5, 0.75, 1.0, 1.0, 1.0]
        @test knot(5, 2, true) == [-2.0, -2.0, -1.0, 0.0, 1.0, 2.0, 2.0]
    end

    @testset "Open NonUniform Proportional Knots" begin
        b = [0. 2 4 6 8; 0 6 3 6 6; 0 0 0 0 0]
        @test typeof(knotc(5, 3, b)) == Array{Float64,1}
        @test isapprox(knotc(5, 3, b), [0.0, 0.0, 0.0, 1.45338, 2.38171, 3.0, 3.0, 3.0]; atol = 1e-6)
    end

    @testset "Periodic Uniform Knots" begin
        @test typeof(knotu(5, 2)) == Array{Float64,1}
        @test knotu(5,2) == [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
        @test knotu(5, 4, 0.25) == [0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0]
    end
end


text:
using GeometricIntegrators.Config
using GeometricIntegrators.Integrators
using GeometricIntegrators.Solvers
using GeometricIntegrators.Tableaus
using GeometricIntegrators.TestProblems.LotkaVolterra2dProblem
using GeometricIntegrators.Utils
using Test

using GeometricIntegrators.TestProblems.LotkaVolterra2dProblem: Δt, nt

set_config(:nls_atol, 8eps())
set_config(:nls_rtol, 2eps())

ode = lotka_volterra_2d_ode()

int  = IntegratorFIRK(ode, getTableauGLRK(8), Δt)
sol  = integrate(int, nt)
refx = sol.q[:,end]

@testset "$(rpad("Special integrators",80))" begin

    pgint = IntegratorPGLRK(ode, getCoefficientsPGLRK(2), Δt)
    pgsol = integrate(pgint, nt)

    # println(rel_err(pgsol.q, refx))
    @test rel_err(pgsol.q, refx) < 4E-12

    pgint = IntegratorPGLRK(ode, getCoefficientsPGLRK(3), Δt)
    pgsol = integrate(pgint, nt)

    # println(rel_err(pgsol.q, refx))
    @test rel_err(pgsol.q, refx) < 2E-16

    pgint = IntegratorPGLRK(ode, getCoefficientsPGLRK(4), Δt)
    pgsol = integrate(pgint, nt)

    # println(rel_err(pgsol.q, refx))
    @test rel_err(pgsol.q, refx) < 2E-16

end


text:
#===============================================================================
# Practice with constructors of composite types or objects.
#
# Author: Carl Cortright
# Date: 12/30/2016
#
===============================================================================#
# This is the type we will construct
type Foo
  bar
  baz
end

# Outer constructors create new instances by calling inner constructors
Foo(x) = Foo(x,x)
Foo() = Foo(0)

# Testing our new constructor
food = Foo()

# Inner constructors allow us to recurse and provide invariants
type OrderedPair
  x::Real
  y::Real

  # Note: the new keyword is locally defined and creates a new instance
  OrderedPair(x, y) = x > y ? error("out of order") : new(x,y)
end

# Seeing if the new constructor works as a longer function
type OrderedPair
  x::Real
  y::Real

  function OrderedPair(x,y)
    if x > y
      error("out of order")
    else
      new(x,y)
    end
  end
end

# Note: once an inner constructor is defined Julia assumes you've defined all the constructors you need
type Foo
  bar
  baz

  Foo(bar, baz) = new(bar, baz)
end

# The new function can be called with less parameters than itself to allow for self-referal
type SelfReferalials
  obj::SelfReferalials

  SelfReferalials() = (x = new(); x.obj = x)
end

x = SelfReferalials()
is(x,x)
is(x,x.obj)
is(x, x.obj.obj.obj)

# Parametric types have new constructors for each type they implement
type Point{T <: Real}
  x::T
  y::T
end

# We can also redefine outer constructors to dispatch on unusual arguments
Point(x::Int64, y::Float64) = Point(convert(x, Float64), y)

max_issues_repo_name: carlcortright/Playground


text:
# This file includes testset for TaskManager

@testset "TaskManager" begin
    @info "Running TaskManagerTestSet ..."

    # Preliminaries.
    gettask(ch) = @async while true
        val = take!(ch)
        val === NaN && break
        val == 0 && error("The task failed.")
        println("Took val" * string(val))
    end

    # TaskManager construction
    struct Mytype1
        x::Int
    end
    comps = [Mytype1(i) for i = 1 : 5]
    chpairs = [Channel(0) for i = 1 : 5]
    comptasks = [gettask(chpair)  for chpair in chpairs]
    ps = Dict(zip(comps, comptasks))
    tm = TaskManager(ps)
    @test checktaskmanager(tm) === nothing  # All tasks are running, nothing is thrown.
    put!(chpairs[1], 0.)     # Fail the trigger task first Mytype1
    put!(chpairs[2], 0.)     # Fail the output task first Mytype1
    @test_throws Exception checktaskmanager(tm)

    @info "Done TaskManagerTestSet."
end # testset


text:
"""
    autotype(X; kw...)

Return a dictionary of suggested scitypes for each column of `X`, a table or
an array based on rules

## Kwargs

* `only_changes=true`:       if true, return only a dictionary of the names for
                              which applying autotype differs from just using
                              the ambient convention. When coercing with
                              autotype, `only_changes` should be true.
* `rules=(:few_to_finite,)`: the set of rules to apply.
"""
autotype(X; kw...) = _autotype(X, Val(ST.trait(X)); kw...)

# For an array object (trait:other)
function _autotype(X::Arr, ::Val{:other};
                   rules::NTuple{N,Symbol} where N=(:few_to_finite,))
    # check that the rules are recognised
    _check_rules(rules)
    # inspect the current element scitype
    sugg_type = elscitype(X)
    # apply rules in sequence
    np = prod(size(X))
    for rule in rules
        if rule === :few_to_finite
            col = vec(X) # needed to compute unique values
            sugg_type = eval(:($rule($sugg_type, $col, $np)))
        else
            col = view(X, :, 1) # only needed for eltype
            sugg_type = eval(:($rule($sugg_type, $col, $np)))
        end
    end
    return sugg_type
end

# For a table object (trait:table)
function _autotype(X, ::Val{:table}; only_changes::Bool=true,
                  rules::NTuple{N,Symbol} where N=(:few_to_finite,))
    # check that the rules are recognised
    _check_rules(rules)
    # recuperate the schema of `X`
    sch = schema(X)
    # dictionary to keep track of the suggested types
    suggested_types = Dict{Symbol,Type{<:Any}}()
    # keep track of the column names for which the suggested type is different
    # than the convention one
    has_changed = Symbol[]
    # go over each column and for each of them apply the rules in order
    # in which they were provided
    columns = Tables.columns(X)
    cols = (Tables.getcolumn(columns, c)
            for c in Tables.columnnames(columns)) |> Tuple
    zipper = zip(sch.names, sch.types, sch.scitypes, cols)
    for (name, machine_type, stype, col) in zipper
        # start with the data type and iterate over the rules to get
        # a suggested type, note that this loop is type unstable but it
        # doesn't really matter, there are few rules and sugg is fast
        sugg_type = stype
        for rule in rules
            sugg_type = eval(:($rule($sugg_type, $col, $sch.nrows)))
        end
        # store the suggested type
        suggested_types[name] = sugg_type
        # check if the suggested type is different than the scitype
        # if so, push it to has_changed to aid an eventual filtering
        if sugg_type != stype
            push!(has_changed, name)
        end
    end
    if only_changes
        filter!(e -> first(e) in has_changed, suggested_types)
    end
    return suggested_types
end

# convenience functions to pass a single rule at the time
autotype(X, rule::Symbol; args...) =
    autotype(X; rules=(rule,), args...)
# convenience function to splat rules
autotype(X, rules::NTuple{N,Symbol} where N; args...) =
    autotype(X; rules=rules, args...)

# -----------------------------------------------------------------
# rules

function _check_rules(rules::NTuple{N,Symbol} where N)
    for rule in rules
        rule in (:few_to_finite,
                 :discrete_to_continuous,
                 :string_to_multiclass) ||
                 throw(ArgumentError("Rule $rule not recognised."))
    end
    return nothing
end

"""
    few_to_finite(type, col, nrows)

For a column `col` with element type `type` and `nrows` rows, check if there
are relatively few values as compared to the number of rows. The heuristic
for "few" is as follows:
1. there are ≤ 3 unique values with more than 5 rows, use `MultiClass{N}` type
2. there are less than 10% unique values out of the number of rows **or**
    there's fewer than 100 unique values (whichever one is smaller):

In both cases:
a. if it's a Real type, return as `OrderedFactor`
b. if it's something else (e.g. a `String`) return as `MultiClass{N}`
"""
function few_to_finite(type::Type, col, nrows::Int)
    nonmissing(type) <: Finite && return type
    unique_vals  = unique(skipmissing(col))
    coltype      = eltype(col)
    nunique_vals = length(unique_vals)
    if nunique_vals ≤ 3 && nrows ≥ 5 ||             # H1
       nunique_vals ≤ max(min(0.1 * nrows, 100), 4) # H2
        # suggest a type
        T = sugg_finite(coltype)
        return T_or_Union_Missing_T(coltype, T)
    end
    return type
end

"""
    discrete_to_continuous(type, _, _)

Return `Continuous` if the current type is either `Count` or `Integer`,
otherwise return the type unchanged.
"""
function discrete_to_continuous(type::Type, _, _)
    if nonmissing(type) <: Union{Count,Integer}
        return T_or_Union_Missing_T(type, Continuous)
    end
    return type
end

"""
    string_to_multiclass(type, col, _)

For a column with element type `<: AbstractString` or `<: AbstractChar` return
Multiclass irrelevant of how many unique values there are. This rule is only
applied on columns which are `Textual` or `Unknown`.
"""
function string_to_multiclass(type::Type, col, _)
    nonmissing(type) in (Textual, Unknown) || return type
    etc = eltype(col)
    if nonmissing(etc) <: Union{AbstractChar,AbstractString}
        return T_or_Union_Missing_T(etc, Multiclass)
    end
    return type
end

# -----------------------------------------------------------------
# Helper functions for the rules

"""
    sugg_finite(type)

Helper function to suggest a finite type corresponding to `T` when there are
few unique values.
"""
function sugg_finite(::Type{<:Union{Missing,T}}) where T
    T <: Real && return OrderedFactor
    return Multiclass
end

"""
    T_or_Union_Missing_T(type, T)

Helper function to return either `T` or `Union{Missing,T}`.
"""
function T_or_Union_Missing_T(type::Type, T::Type)
    return ifelse(type >: Missing, Union{Missing, T}, T)
end


text:
mutable struct TripleBandLinearOp{FM <: FdmMesher}
  direction::Int
  mesher::FM
  i0::Vector{Int}
  i2::Vector{Int}
  reverseIndex::Vector{Int}
  lower::Vector{Float64}
  _diag::Vector{Float64}
  upper::Vector{Float64}
end

function TripleBandLinearOp(direction::Int, mesher::FdmMesher)
  sz = mesher.layout.size
  i0 = Vector{Int}(undef, sz)
  i2 = Vector{Int}(undef, sz)
  reverseIndex = Vector{Int}(undef, sz)
  lower = Vector{Float64}(undef, sz)
  _diag = Vector{Float64}(undef, sz)
  upper = Vector{Float64}(undef, sz)

  newDim = copy(mesher.layout.dim)
  newDim[1], newDim[direction] = newDim[direction], newDim[1] # swap

  newSpacing = FdmLinearOpLayout(newDim).spacing
  newSpacing[1], newSpacing[direction] = newSpacing[direction], newSpacing[1] # swap

  coords = ones(Int, length(mesher.layout.dim))
  @inbounds @simd for i = 1:sz
    ## TripleBandLinearOp part
    i0[i] = neighborhood(mesher.layout, i, coords, direction, -1)
    i2[i] = neighborhood(mesher.layout, i, coords, direction, 1)

    newIndex = dot(coords .- 1, newSpacing) + 1

    reverseIndex[newIndex] = i

    iter_coords!(coords, mesher.layout.dim)
  end

  return TripleBandLinearOp(direction, mesher, i0, i2, reverseIndex, lower, _diag, upper)
end

function FirstDerivativeOp(direction::Int, mesher::FdmMesher)
  sz = mesher.layout.size
  i0 = Vector{Int}(undef, sz)
  i2 = Vector{Int}(undef, sz)
  reverseIndex = Vector{Int}(undef, sz)
  lower = Vector{Float64}(undef, sz)
  _diag = Vector{Float64}(undef, sz)
  upper = Vector{Float64}(undef, sz)

  newDim = copy(mesher.layout.dim)
  newDim[1], newDim[direction] = newDim[direction], newDim[1] # swap

  newSpacing = FdmLinearOpLayout(newDim).spacing
  newSpacing[1], newSpacing[direction] = newSpacing[direction], newSpacing[1] # swap

  coords = ones(Int, length(mesher.layout.dim))
  @inbounds @simd for i = 1:sz
    ## TripleBandLinearOp part
    i0[i] = neighborhood(mesher.layout, i, coords, direction, -1)
    i2[i] = neighborhood(mesher.layout, i, coords, direction, 1)

    newIndex = dot(coords .- 1, newSpacing) + 1

    reverseIndex[newIndex] = i

    ## FirstDerivativeOp part
    hm = get_dminus(mesher, coords, direction)
    hp = get_dplus(mesher, coords, direction)

    zetam1 = hm * (hm + hp)
    zeta0 = hm * hp
    zetap1 = hp * (hm + hp)

    if coords[direction] == 1
      # upwinding scheme
      lower[i] = 0.0
      upper[i] = 1.0 / hp
      _diag[i] = -(upper[i])
    elseif coords[direction] == mesher.layout.dim[direction]
      # downwinding scheme
      _diag[i] = 1.0 / hm
      lower[i] = -_diag[i]
      upper[i] = 0.0
    else
      lower[i] = -hp / zetam1
      _diag[i] = (hp - hm) / zeta0
      upper[i] = hm / zetap1
    end
    iter_coords!(coords, mesher.layout.dim)
  end

  return TripleBandLinearOp(direction, mesher, i0, i2, reverseIndex, lower, _diag, upper)
end

function SecondDerivativeOp(direction::Int, mesher::FdmMesher)
  sz = mesher.layout.size
  i0 = Vector{Int}(undef, sz)
  i2 = Vector{Int}(undef, sz)
  reverseIndex = Vector{Int}(undef, sz)
  lower = Vector{Float64}(undef, sz)
  _diag = Vector{Float64}(undef, sz)
  upper = Vector{Float64}(undef, sz)

  newDim = copy(mesher.layout.dim)
  newDim[1], newDim[direction] = newDim[direction], newDim[1] # swap

  newSpacing = FdmLinearOpLayout(newDim).spacing
  newSpacing[1], newSpacing[direction] = newSpacing[direction], newSpacing[1] # swap

  coords = ones(Int, length(mesher.layout.dim))
  @inbounds @simd for i = 1:sz
    ## TripleBandLinearOp part
    i0[i] = neighborhood(mesher.layout, i, coords, direction, -1)
    i2[i] = neighborhood(mesher.layout, i, coords, direction, 1)

    newIndex = dot(coords .- 1, newSpacing) + 1

    reverseIndex[newIndex] = i

    ## SecondDerivativeOp part
    hm = get_dminus(mesher, coords, direction)
    hp = get_dplus(mesher, coords, direction)

    zetam1 = hm * (hm + hp)
    zeta0 = hm * hp
    zetap1 = hp * (hm + hp)

    co = coords[direction]

    if co == 1 || co == mesher.layout.dim[direction]
      lower[i] = upper[i] = _diag[i] = 0.0
    else
      lower[i] = 2.0 / zetam1
      _diag[i] = -2.0 / zeta0
      upper[i] = 2.0 / zetap1
    end
    iter_coords!(coords, mesher.layout.dim)
  end

  return TripleBandLinearOp(direction, mesher, i0, i2, reverseIndex, lower, _diag, upper)
end

function mult!(trpBandLinOp::TripleBandLinearOp, u::Vector{T}) where {T <: Number}
  # probably should protect for array dims here
  trpBandLinOp.lower = trpBandLinOp.lower .* u
  trpBandLinOp._diag = trpBandLinOp._diag .* u
  trpBandLinOp.upper = trpBandLinOp.upper .* u

  return trpBandLinOp
end

function add!(trpBandLinOp::TripleBandLinearOp, m::TripleBandLinearOp)
  trpBandLinOp.lower += m.lower
  trpBandLinOp._diag += m._diag
  trpBandLinOp.upper += m.upper

  return trpBandLinOp
end

function axpyb!(trpBandLinOp::TripleBandLinearOp, a::Vector{Float64}, x::TripleBandLinearOp, y::TripleBandLinearOp, b::Vector{Float64})
  sz = trpBandLinOp.mesher.layout.size

  if length(a) == 0
    if length(b) == 0
      trpBandLinOp._diag[1:sz] = y._diag
      trpBandLinOp.lower[1:sz] = y.lower
      trpBandLinOp.upper[1:sz] = y.upper
    else
      addB = length(b) > 1 ? b[1:sz] : b[1]
      trpBandLinOp._diag[1:sz] = y._diag + addB
      trpBandLinOp.lower[1:sz] = y.lower
      trpBandLinOp.upper[1:sz] = y.upper
    end
  elseif length(b) == 0
    # this might be improved
    trpBandLinOp._diag[1:sz] = y._diag + (length(a) > 0 ? a .* x._diag : a[1] * x._diag)
    trpBandLinOp.lower[1:sz] = y.lower + (length(a) > 0 ? a .* x.lower : a[1] * x.lower)
    trpBandLinOp.upper[1:sz] = y.upper + (length(a) > 0 ? a .* x.upper : a[1] * x.upper)
  else
    addB = length(b) > 1 ? b[1:sz] : b[1]
    trpBandLinOp._diag[1:sz] = y._diag + (length(a) > 0 ? a .* x._diag : a[1] * x._diag) + addB
    trpBandLinOp.lower[1:sz] = y.lower + (length(a) > 0 ? a .* x.lower : a[1] * x.lower)
    trpBandLinOp.upper[1:sz] = y.upper + (length(a) > 0 ? a .* x.upper : a[1] * x.upper)
  end

  return trpBandLinOp
end

function apply(trpBandLinOp::TripleBandLinearOp, r::Vector{Float64})
  idx = trpBandLinOp.mesher.layout
  length(r) == idx.size || error("inconsistent length of r")

  retArray = Vector{Float64}(undef, length(r))

  for i = 1:idx.size
    retArray[i] = r[trpBandLinOp.i0[i]] * trpBandLinOp.lower[i] + r[i] * trpBandLinOp._diag[i] + r[trpBandLinOp.i2[i]] * trpBandLinOp.upper[i]
  end

  return retArray
end

function solve_splitting(trpBandLinOp::TripleBandLinearOp, r::Vector{Float64}, a::Float64, b::Float64)
  layout = trpBandLinOp.mesher.layout
  length(r) == layout.size || error("inconsistent length of r")

  retArray = Vector{Float64}(undef, length(r))
  tmp = Vector{Float64}(undef, length(r))

  # solving a tridiagonal system (we could use Julia built in f'ns here)
  rim1 = trpBandLinOp.reverseIndex[1]
  bet = 1.0 / (a * trpBandLinOp._diag[rim1] + b)
  bet != 0.0 || error("division by zero")

  retArray[rim1] = r[rim1] * bet

  @inbounds @simd for j = 2:layout.size
    ri = trpBandLinOp.reverseIndex[j]
    tmp[j] = a * trpBandLinOp.upper[rim1] * bet

    bet = b + a * (trpBandLinOp._diag[ri] - tmp[j] * trpBandLinOp.lower[ri])
    bet != 0.0 || error("division by zero")

    bet = 1.0 / bet

    retArray[ri] = (r[ri] - a * trpBandLinOp.lower[ri] * retArray[rim1]) * bet
    rim1 = ri
  end

  @inbounds @simd for j = layout.size - 1:-1:2
    retArray[trpBandLinOp.reverseIndex[j]] -= tmp[j + 1] * retArray[trpBandLinOp.reverseIndex[j + 1]]
  end

  retArray[trpBandLinOp.reverseIndex[1]] -= tmp[2] * retArray[trpBandLinOp.reverseIndex[2]]

  return retArray
end

mutable struct SecondOrderMixedDerivativeOp{FD <: FdmMesher} <: NinePointLinearOp
  d1::Int
  d2::Int
  i00::Vector{Int}
  i10::Vector{Int}
  i20::Vector{Int}
  i01::Vector{Int}
  i21::Vector{Int}
  i02::Vector{Int}
  i12::Vector{Int}
  i22::Vector{Int}
  a00::Vector{Float64}
  a10::Vector{Float64}
  a20::Vector{Float64}
  a01::Vector{Float64}
  a11::Vector{Float64}
  a21::Vector{Float64}
  a02::Vector{Float64}
  a12::Vector{Float64}
  a22::Vector{Float64}
  mesher::FD
end

function SecondOrderMixedDerivativeOp(d1::Int, d2::Int, mesher::FD) where {FD <: FdmMesher}
  sz = mesher.layout.size
  i00 = Vector{Int}(undef, sz)
  i10 = Vector{Int}(undef, sz)
  i20 = Vector{Int}(undef, sz)
  i01 = Vector{Int}(undef, sz)
  i21 = Vector{Int}(undef, sz)
  i02 = Vector{Int}(undef, sz)
  i12 = Vector{Int}(undef, sz)
  i22 = Vector{Int}(undef, sz)
  a00 = Vector{Float64}(undef, sz)
  a10 = Vector{Float64}(undef, sz)
  a20 = Vector{Float64}(undef, sz)
  a01 = Vector{Float64}(undef, sz)
  a11 = Vector{Float64}(undef, sz)
  a21 = Vector{Float64}(undef, sz)
  a02 = Vector{Float64}(undef, sz)
  a12 = Vector{Float64}(undef, sz)
  a22 = Vector{Float64}(undef, sz)

  coords = ones(Int, length(mesher.layout.dim))
  @inbounds @simd for i = 1:sz
    # NinePointLinearOp part
    i10[i] = neighborhood(mesher.layout, i, coords, d2, -1)
    i01[i] = neighborhood(mesher.layout, i, coords, d1, -1)
    i21[i] = neighborhood(mesher.layout, i, coords, d1,  1)
    i12[i] = neighborhood(mesher.layout, i, coords, d2,  1)
    i00[i] = neighborhood(mesher.layout, i, coords, d1, -1, d2, -1)
    i20[i] = neighborhood(mesher.layout, i, coords, d1,  1, d2, -1)
    i02[i] = neighborhood(mesher.layout, i, coords, d1, -1, d2,  1)
    i22[i] = neighborhood(mesher.layout, i, coords, d1,  1, d2,  1)

    # SecondOrderMixedDerivativeOp part
    hm_d1 = get_dminus(mesher, coords, d1)
    hp_d1 = get_dplus(mesher, coords, d1)
    hm_d2 = get_dminus(mesher, coords, d2)
    hp_d2 = get_dplus(mesher, coords, d2)

    zetam1 = hm_d1 * (hm_d1 + hp_d1)
    zeta0 = hm_d1 * hp_d1
    zetap1 = hp_d1 * (hm_d1 + hp_d1)

    phim1 = hm_d2 * (hm_d2 + hp_d2)
    phi0 = hm_d2 * hp_d2
    phip1 = hp_d2 * (hm_d2 + hp_d2)

    c1 = coords[d1]
    c2 = coords[d2]

    if c1 == 1 && c2 == 1
      # lower left corner
      a00[i] = a01[i] = a02[i] = a10[i] = a20[i] = 0.0
      a11[i] = a22[i] = 1.0 / (hp_d1 * hp_d2)
      a21[i] = a12[i] = -a11[i]
    elseif c1 == mesher.layout.dim[d1] && c2 == 1
      # upper left corner
      a22[i] = a21[i] = a20[i] = a10[i] = a00[i] = 0.0
      a01[i] = a12[i] = 1.0 / (hm_d1 * hp_d2)
      a11[i] = a02[i] = -a01[i]
    elseif c1 == 1 && c2 == mesher.layout.dim[d2]
      # lower right corner
      a00[i] = a01[i] = a02[i] = a12[i] = a22[i] = 0.0
      a10[i] = a21[i] = 1.0 / (hp_d1 * hm_d2)
      a20[i] = a11[i] = -a10[i]
    elseif c1 == mesher.layout.dim[d1] && c2 == mesher.layout.dim[d2]
      # upper right corner
      a20[i] = a21[i] = a22[i] = a12[i] = a02[i] = 0.0
      a00[i] = a11[i] = 1.0 / (hm_d1 * hm_d2)
      a10[i] = a01[i] = -a00[i]
    elseif c1 == 1
      # lower side
      a00[i] = a01[i] = a02[i] = 0.0
      a10[i] = hp_d2 / (hp_d1 * phim1)
      a20[i] = -a10[i]
      a21[i] = (hp_d2 - hm_d2) / (hp_d1 * phi0)
      a11[i] = -a21[i]
      a22[i] = hm_d2 / (hp_d1 * phip1)
      a12[i] = -a22[i]
    elseif c1 == mesher.layout.dim[d1]
      # upper side
      a20[i] = a21[i] = a22[i] = 0.0
      a00[i] = hp_d2 / (hm_d1 * phim1)
      a10[i] = -a00[i]
      a11[i] = (hp_d2 - hm_d2) / (hm_d1 * phi0)
      a01[i] = -a11[i]
      a12[i] = hm_d2 / (hm_d1 * phip1)
      a02[i] = -a12[i]
    elseif c2 == 1
      # left side
      a00[i] = a10[i] = a20[i] = 0.0
      a01[i] = hp_d1 / (zetam1 * hp_d2)
      a02[i] = -a01[i]
      a12[i] = (hp_d1 - hm_d1) / (zeta0 * hp_d2)
      a11[i] = -a12[i]
      a22[i] = hm_d1 / (zetap1 * hp_d2)
      a21[i] = -a22[i]
    elseif c2 == mesher.layout.dim[d2]
      # right side
      a22[i] = a12[i] = a02[i] = 0.0
      a00[i] = hp_d1 / (zetam1 * hm_d2)
      a01[i] = -a00[i]
      a11[i] = (hp_d1 - hm_d1) / (zeta0 * hm_d2)
      a10[i] = -a11[i]
      a21[i] = hm_d1 / (zetap1 * hm_d2)
      a20[i] = -a21[i]
    else
      a00[i] = hp_d1 * hp_d2 / (zetam1 * phim1)
      a10[i] = -(hp_d1 - hm_d1) * hp_d2 / (zeta0 * phim1)
      a20[i] = -hm_d1 * hp_d2 / (zetap1 * phim1)
      a01[i] = -hp_d1 * (hp_d2 - hm_d2) / (zetam1 * phi0)
      a11[i] = (hp_d1 - hm_d1) * (hp_d2 - hm_d2) / (zeta0 * phi0)
      a21[i] = hm_d1 * (hp_d2 - hm_d2) / (zetap1 * phi0)
      a02[i] = -hp_d1 * hm_d2 / (zetam1 * phip1)
      a12[i] = hm_d2 * (hp_d1 - hm_d1) / (zeta0 * phip1)
      a22[i] = hm_d1 * hm_d2 / (zetap1 * phip1)
    end
    iter_coords!(coords, mesher.layout.dim)
  end

  return SecondOrderMixedDerivativeOp{FD}(d1, d2, i00, i10, i20, i01, i21, i02, i12, i22, a00, a10, a20, a01, a11, a21, a02, a12, a22, mesher)
end

function mult!(ninePointLin::NinePointLinearOp, u::Vector{T}) where {T <: Number}
  ninePointLin.a11 = ninePointLin.a11 .* u
  ninePointLin.a01 = ninePointLin.a01 .* u
  ninePointLin.a10 = ninePointLin.a10 .* u
  ninePointLin.a21 = ninePointLin.a21 .* u
  ninePointLin.a22 = ninePointLin.a22 .* u
  ninePointLin.a00 = ninePointLin.a00 .* u
  ninePointLin.a02 = ninePointLin.a02 .* u
  ninePointLin.a20 = ninePointLin.a20 .* u
  ninePointLin.a12 = ninePointLin.a12 .* u

  return ninePointLin
end

function apply(ninePointLin::NinePointLinearOp, u::Vector{Float64})
  length(u) == ninePointLin.mesher.layout.size || error("inconsistent length of r")

  retVal = ninePointLin.a00 .* u[ninePointLin.i00] + ninePointLin.a01 .* u[ninePointLin.i01] + ninePointLin.a02 .* u[ninePointLin.i02] + ninePointLin.a10 .* u[ninePointLin.i10] +
          ninePointLin.a11 .* u + ninePointLin.a12 .* u[ninePointLin.i12] + ninePointLin.a20 .* u[ninePointLin.i20] + ninePointLin.a21 .* u[ninePointLin.i21] +
          ninePointLin.a22 .* u[ninePointLin.i22]

  return retVal
end

mutable struct FdmG2Op{FD <: FdmMesher} <: FdmLinearOpComposite
  direction1::Int
  direction2::Int
  x::Vector{Float64}
  y::Vector{Float64}
  dxMap::TripleBandLinearOp{FD}
  dyMap::TripleBandLinearOp{FD}
  corrMap::SecondOrderMixedDerivativeOp{FD}
  mapX::TripleBandLinearOp{FD}
  mapY::TripleBandLinearOp{FD}
  model::G2
end

function FdmG2Op(mesher::FD, model::G2, direction1::Int, direction2::Int) where {FD <: FdmMesher}
  x = get_locations(mesher, direction1)
  y = get_locations(mesher, direction2)

  dxMap = add!(mult!(FirstDerivativeOp(direction1, mesher), (-x * get_a(model))),
              mult!(SecondDerivativeOp(direction1, mesher), (0.5 * get_sigma(model) * get_sigma(model) * ones(mesher.layout.size))))

  dyMap = add!(mult!(FirstDerivativeOp(direction2, mesher), (-y * get_b(model))),
              mult!(SecondDerivativeOp(direction2, mesher), (0.5 * get_eta(model) * get_eta(model) * ones(mesher.layout.size))))

  corrMap = mult!(SecondOrderMixedDerivativeOp(direction1, direction2, mesher), fill(get_rho(model) * get_sigma(model) * get_eta(model), mesher.layout.size))

  mapX = TripleBandLinearOp(direction1, mesher)
  mapY = TripleBandLinearOp(direction2, mesher)

  return FdmG2Op{FD}(direction1, direction2, x, y, dxMap, dyMap, corrMap, mapX, mapY, model)
end

mutable struct FdmHullWhiteOp{FD <: FdmMesher} <: FdmLinearOpComposite
  direction::Int
  x::Vector{Float64}
  dzMap::TripleBandLinearOp{FD}
  mapT::TripleBandLinearOp{FD}
  model::HullWhite
end

function FdmHullWhiteOp(mesher::FD, model::HullWhite, direction::Int) where {FD <: FdmMesher}
  x = get_locations(mesher, direction)

  dzMap = add!(mult!(FirstDerivativeOp(direction, mesher), (-x * get_a(model))),
              mult!(SecondDerivativeOp(direction, mesher), (0.5 * get_sigma(model) * get_sigma(model)) * ones(mesher.layout.size)))

  mapT = TripleBandLinearOp(direction, mesher)

  return FdmHullWhiteOp{FD}(direction, x, dzMap, mapT, model)
end

# OP methods #
function set_time!(op::FdmG2Op, t1::Float64, t2::Float64)
  dynamics = get_dynamics(op.model)

  phi = 0.5 * (short_rate(dynamics, t1, 0.0, 0.0) + short_rate(dynamics, t2, 0.0, 0.0))

  hr = -0.5 * (op.x + op.y .+ phi)
  axpyb!(op.mapX, Vector{Float64}(), op.dxMap, op.dxMap, hr)
  axpyb!(op.mapY, Vector{Float64}(), op.dyMap, op.dyMap, hr)

  return op
end

function set_time!(op::FdmHullWhiteOp, t1::Float64, t2::Float64)
  dynamics = get_dynamics(op.model)

  phi = 0.5 * (short_rate(dynamics, t1, 0.0) + short_rate(dynamics, t2, 0.0))

  axpyb!(op.mapT, Vector{Float64}(), op.dzMap, op.dzMap, -(op.x .+ phi))

  return op
end

function apply_direction(op::FdmG2Op, direction::Int, r::Vector{Float64})
  if direction == op.direction1
    return apply(op.mapX, r)
  elseif direction == op.direction2
    return apply(op.mapY, r)
  else
    return zeros(length(r))
  end
end

function apply_direction(op::FdmHullWhiteOp, direction::Int, r::Vector{Float64})
  if direction == op.direction
    return apply(op.mapT, r)
  else
    return zeros(length(r))
  end
end

function solve_splitting(op::FdmG2Op, direction::Int, r::Vector{Float64}, a::Float64)
  if direction == op.direction1
    return solve_splitting(op.mapX, r, a, 1.0)
  elseif direction == op.direction2
    return solve_splitting(op.mapY, r, a, 1.0)
  else
    return zeros(length(r))
  end
end

function solve_splitting(op::FdmHullWhiteOp, direction::Int, r::Vector{Float64}, a::Float64)
  if direction == op.direction
    return solve_splitting(op.mapT, r, a, 1.0)
  else
    return zeros(length(r))
  end
end

apply_mixed(op::FdmG2Op, r::Vector{Float64}) = apply(op.corrMap, r)
apply_mixed(op::FdmHullWhiteOp, r::Vector{Float64}) = zeros(length(r))

apply(op::FdmG2Op, r::Vector{Float64}) = apply(op.mapX, r) + apply(op.mapY, r) + apply_mixed(op, r)
apply(op::FdmHullWhiteOp, r::Vector{Float64}) = apply(op.mapT, r)

get_size(op::FdmG2Op) = 2
get_size(op::FdmHullWhiteOp) = 1


text:
"""
    wald(obj::EconometricModel, vce::VCE = obj.vce)

Provides the Wald statistic as Wald value, a F-Distribution, and the associated p-value.
"""
function wald(obj::EconometricModel, vce::VCE = obj.vce)
    β = coef(obj)
    V = vcov(obj, vce)
    k = length(β) - 1
    rdf = dof_residual(obj)
    R = hasintercept(obj) ? hcat(zeros(k), I) : I
    Bread = R * β
    Meat = inv(bunchkaufman!(Hermitian(R * V * R')))
    W = (Bread' * Meat * Bread) / size(R, 1)
    if size(R, 1) > 0
        F = FDist(size(R, 1), rdf)
        p = ccdf(F, W)
    else
        F = FDist(1, 1)
        p = NaN
    end
    W, F, p
end


text:
"""
FinEtools (C) 2017-2019, Petr Krysl

Finite Element tools.  Julia implementation  of the finite element method
for continuum mechanics. Package for heat diffusion problems.
"""
module FinEtoolsAcoustics

__precompile__(true)

include("allmodules.jl")

# Exports follow:

###########################################################################
# Acoustics functionality
###########################################################################
using FinEtoolsAcoustics.MatAcoustFluidModule: MatAcoustFluid
# Exported: type of acoustic fluid material
export MatAcoustFluid

using FinEtoolsAcoustics.FEMMAcoustModule: FEMMAcoust, acousticmass, nzebcloadsacousticmass,
acousticstiffness, nzebcloadsacousticstiffness
# Exported: type for linear acoustics  and discretization methods
export FEMMAcoust, acousticmass, nzebcloadsacousticmass,
    acousticstiffness, nzebcloadsacousticstiffness

using FinEtoolsAcoustics.FEMMAcoustSurfModule: FEMMAcoustSurf, acousticABC, pressure2resultantforce, pressure2resultanttorque, acousticcouplingpanels
# Exported: type for acoustic absorbing boundary condition  and  transformation matrices from pressure  to resultants
export FEMMAcoustSurf, acousticABC, pressure2resultantforce, pressure2resultanttorque, acousticcouplingpanels

end # module


text:
"""
    update_and_close_viewers!(wss::Vector{HTTP.WebSockets.WebSocket})

Take a list of viewers, i.e. WebSocket connections from a client,
send a message with data "update" to each of them (to trigger a page reload),
then close the connection. Finally, empty the list since all connections are
closing anyway and clients will re-connect from the re-loaded page.
"""
function update_and_close_viewers!(wss::Vector{HTTP.WebSockets.WebSocket})
    foreach(wss) do wsi
        try
            write(wsi, "update")
            close(wsi)
        catch
        end
    end
    empty!(wss)
    return nothing
end

"""
    file_changed_callback(f_path::AbstractString)

Function reacting to the change of the file at `f_path`. Is set as callback for
the file watcher.
"""
function file_changed_callback(f_path::AbstractString)
    if VERBOSE[]
        println("ℹ [LiveUpdater]: Reacting to change in file '$f_path'...")
    end
    if endswith(f_path, ".html")
        # if html file, update viewers of this file only
        # check the viewer still exists otherwise may error
        # see issue https://github.com/asprionj/LiveServer.jl/issues/95
        if haskey(WS_VIEWERS, f_path)
            update_and_close_viewers!(WS_VIEWERS[f_path])
        end
    else
        # otherwise (e.g. modification to a CSS file), update all viewers
        foreach(update_and_close_viewers!, values(WS_VIEWERS))
    end
    return nothing
end

"""
    get_fs_path(req_path::AbstractString)

Return the filesystem path corresponding to a requested path, or an empty
String if the file was not found.
"""
function get_fs_path(req_path::AbstractString)::String
    # first element after the split is **always** "/"
    r_parts = HTTP.URIs.unescapeuri.(split(HTTP.URI(req_path).path[2:end], "/"))
    fs_path = joinpath(r_parts...)
    if !isempty(CONTENT_DIR[])
        fs_path = joinpath(CONTENT_DIR[], fs_path)
    end
    # if no file is specified, try to append `index.html` and see
    endswith(req_path, "/") && (fs_path = joinpath(fs_path, "index.html"))
    # either the result is a valid file path in which case it's returned otherwise ""
    if isfile(fs_path) || isdir(fs_path)
        return fs_path
    else
        return ""
    end
end

"""
    serve_file(fw, req::HTTP.Request; inject_browser_reload_script::Bool = true)

Handler function for serving files. This takes a file watcher, to which files to be watched can be
added, and a request (e.g. a path entered in a tab of the browser), and converts it to the
appropriate file system path. If the path corresponds to a HTML file, it will inject the reloading
`<script>` (see file `client.html`) at the end of its body, i.e. directly before the `</body>` tag.
All files served are added to the file watcher, which is responsible to check whether they're
already watched or not. Finally the file is served via a 200 (successful) response. If the file
does not exist, a response with status 404 and an according message is sent.
"""
function serve_file(fw, req::HTTP.Request; inject_browser_reload_script::Bool = true)
    ret_code = 200
    fs_path = get_fs_path(req.target)
    # in case the path was not resolved, return a 404
    if isempty(fs_path)
        ret_code = 404
        # Check if /404/ or /404.html exists and serve that as a body
        for f in ("/404/", "/404.html")
            maybe_path = get_fs_path(f)
            if !isempty(maybe_path)
                fs_path = maybe_path
                break
            end
        end
        # If still not found a body, return a generic error message
        if isempty(fs_path)
            return HTTP.Response(404, "404: file not found. The most likely reason " *
                                      "is that the URL you entered has a mistake in it or that " *
                                      "the requested page has been deleted or renamed. Check " *
                                      "also that the server is still running.")
        end
    end

    # Respond with 301 if the path is a directory
    if isdir(fs_path)
        return HTTP.Response(301, ["Location" => req.target * "/"])
    end

    ext = last(splitext(fs_path))[2:end]
    content = read(fs_path, String)

    # build the response with appropriate mime type (this is inspired from Mux
    # https://github.com/JuliaWeb/Mux.jl/blob/master/src/examples/files.jl)
    mime = get(MIME_TYPES, ext, "application/octet-stream")

    # if html, add the browser-sync script to it
    if (inject_browser_reload_script) && (mime in ["text/html", "application/xhtml+xml"])
        end_body_match = match(r"</body>", content)
        if end_body_match === nothing
            # no </body> tag found, trying to add the reload script at the end; this may fail.
            content *= BROWSER_RELOAD_SCRIPT
        else
            end_body = prevind(content, end_body_match.offset)
            # reconstruct the page with the reloading script
            io = IOBuffer()
            write(io, SubString(content, 1:end_body))
            write(io, BROWSER_RELOAD_SCRIPT)
            write(io, SubString(content, nextind(content, end_body):lastindex(content)))
            content = take!(io)
        end
    end

    headers = ["Content-Type" => mime]
    resp = HTTP.Response(ret_code, content)
    isempty(headers) || (resp.headers = HTTP.mkheaders(headers))

    # add the file to the file watcher
    watch_file!(fw, fs_path)

    # return the response
    return resp
end

"""
    ws_upgrade(http::HTTP.Stream)

Upgrade the HTTP request in the stream to a websocket.
"""
function ws_upgrade(http::HTTP.Stream)
    # adapted from HTTP.WebSockets.upgrade; note that here the upgrade will always
    # have  the right format as it always triggered by after a Response
    HTTP.setstatus(http, 101)
    HTTP.setheader(http, "Upgrade" => "websocket")
    HTTP.setheader(http, "Connection" => "Upgrade")
    key = HTTP.header(http, "Sec-WebSocket-Key")
    HTTP.setheader(http, "Sec-WebSocket-Accept" => HTTP.WebSockets.accept_hash(key))
    HTTP.startwrite(http)

    io = http.stream
    return HTTP.WebSockets.WebSocket(io; server=true)
end

"""
    ws_tracker(ws::HTTP.WebSockets.WebSocket, target::AbstractString)

Adds the websocket connection to the viewers in the global dictionary `WS_VIEWERS` to the entry
corresponding to the targeted file.
"""
function ws_tracker(ws::HTTP.WebSockets.WebSocket, target::AbstractString)
    # add to list of html files being "watched"
    # NOTE: this file always exists because the query is generated just after serving it
    fs_path = get_fs_path(target)

    # if the file is already being viewed, add ws to it (e.g. several tabs)
    # otherwise add to dict
    if fs_path ∈ keys(WS_VIEWERS)
        push!(WS_VIEWERS[fs_path], ws)
    else
        WS_VIEWERS[fs_path] = [ws]
    end

    try
        # NOTE: browsers will drop idle websocket connections so this effectively
        # forces the websocket to stay open until it's closed by LiveServer (and
        # not by the browser) upon writing a `update` message on the websocket.
        # See update_and_close_viewers
        while isopen(ws.io)
            sleep(0.1)
        end
    catch err
        # NOTE: there may be several sources of errors caused by the precise moment
        # at which the user presses CTRL+C and after what events. In an ideal world
        # we would check that none of these errors have another source but for now
        # we make the assumption it's always the case (note that it can cause other
        # errors than InterruptException, for instance it can cause errors due to
        # stream not being available etc but these all have the same source).
        # - We therefore do not propagate the error but merely store the information that
        # there was a forcible interruption of the websocket so that the interruption
        # can be guaranteed to be propagated.
        WS_INTERRUPT[] = true
    end
    return nothing
end

"""
    serve(filewatcher; host="127.0.0.1", port=8000, dir="", verbose=false, coreloopfun=(c,fw)->nothing, inject_browser_reload_script::Bool = true)

Main function to start a server at `http://host:port` and render what is in the current
directory. (See also [`example`](@ref) for an example folder).

* `filewatcher` is a file watcher implementing the API described for [`SimpleWatcher`](@ref) (which also is the default) and messaging the viewers (via WebSockets) upon detecting file changes.
* `port` is an integer between 8000 (default) and 9000.
* `dir` specifies where to launch the server if not the current working directory.
* `verbose` is a boolean switch to make the server print information about file changes and connections.
* `coreloopfun` specifies a function which can be run every 0.1 second while the liveserver is going; it takes two arguments: the cycle counter and the filewatcher. By default the coreloop does nothing.

# Example

```julia
LiveServer.example()
serve(host="127.0.0.1", port=8080, dir="example", verbose=true)
```

If you open a browser to `http://127.0.0.1:8080/`, you should see the `index.html` page from the
`example` folder being rendered. If you change the file, the browser will automatically reload the
page and show the changes.
"""
function serve(fw::FileWatcher=SimpleWatcher(file_changed_callback);
               host::String="127.0.0.1", port::Int=8000, dir::AbstractString="", verbose::Bool=false,
               coreloopfun::Function=(c,fw)->nothing,
               inject_browser_reload_script::Bool = true)

    8000 ≤ port ≤ 9000 || throw(ArgumentError("The port must be between 8000 and 9000."))
    setverbose(verbose)

    if !isempty(dir)
        isdir(dir) || throw(ArgumentError("The specified dir '$dir' is not recognised."))
        CONTENT_DIR[] = dir
    end

    start(fw)

    # make request handler
    req_handler = HTTP.RequestHandlerFunction(req -> serve_file(fw, req; inject_browser_reload_script = inject_browser_reload_script))

    server = Sockets.listen(parse(IPAddr, host), port)
    println("✓ LiveServer listening on http://$(host == string(Sockets.localhost) ? "localhost" : host):$port/ ...\n  (use CTRL+C to shut down)")
    @async HTTP.listen(host, port;
                       server=server, readtimeout=0, reuse_limit=0) do http::HTTP.Stream
        # reuse_limit=0 ensures that there won't be an error if killing and restarting the server.
        if HTTP.WebSockets.is_upgrade(http.message)
            # upgrade to websocket
            ws = ws_upgrade(http)
            # add to list of viewers and keep open until written to
            ws_tracker(ws, http.message.target)
        else
            # handle HTTP request
            HTTP.handle(req_handler, http)
        end
    end

    # wait until user interrupts the LiveServer (using CTRL+C).
    try
        counter = 1
        while true
            if WS_INTERRUPT.x || fw.status == :interrupted
                # rethrow the interruption (which may have happened during
                # the websocket handling or during the file watching)
                throw(InterruptException())
            end
            # do the auxiliary function if there is one (by default this does nothing)
            coreloopfun(counter, fw)
            # update the cycle counter and sleep (yields to other threads)
            counter += 1
            sleep(0.1)
        end
    catch err
        if !isa(err, InterruptException)
            throw(err)
        end
    finally
        # cleanup: close everything that might still be alive
        VERBOSE[] && println("\n⋮ shutting down LiveServer")
        # stop the filewatcher
        stop(fw)
        # close any remaining websockets
        for wss ∈ values(WS_VIEWERS), wsi ∈ wss
            close(wsi.io)
        end
        # empty the dictionary of viewers
        empty!(WS_VIEWERS)
        # shut down the server
        close(server)
        VERBOSE[] && println("\n✓ LiveServer shut down.")
        # reset environment variables
        CONTENT_DIR[] = ""
        WS_INTERRUPT[] = false
    end
    return nothing
end


text:
# module Bukdu

function render(T::Type{Val{:Markdown}}, text::String)::Conn
    filtering(render,T,text) do
        Conn(200, Dict("Content-Type"=>"text/html"), chomp(Markdown.html(Markdown.parse(text))))
    end
end


text:
mutable struct NormalMatrixNormLikelihood <: MyXMLElement
    el::XMLOrNothing
    precision_parameter::MyXMLElement
    matrix_parameter::MyXMLElement
    id::String

    function NormalMatrixNormLikelihood(prec::MyXMLElement, mat::MyXMLElement,
                                        id::String)
        return new(nothing, prec, mat, id)
    end
end

function make_xml(nmnl::NormalMatrixNormLikelihood)
    @unpack el, precision_parameter, matrix_parameter, id = nmnl
    el = new_element(bn.NORMAL_MATRIX_NORM)
    set_id!(el, id)

    gp_el = new_child(el, bn.GLOBAL_PRECISION)
    make_xml(precision_parameter)
    add_ref_el(gp_el, precision_parameter.el)

    m_el = new_child(el, bn.MATRIX)
    make_xml(matrix_parameter)
    add_ref_el(m_el, matrix_parameter.el)

    nmnl.el = el
    return el
end

function get_normal_prior(nmnl::NormalMatrixNormLikelihood)
    make_xml(nmnl)
    return nmnl.el
end