text:
"""
    kinship_pruning(grm::Matrix; method=:bottom_up, cutoff=0.125, min_samples=0)

Return a BitVector indicating which samples remains after kinship pruning.
Valid values for `method` are `:bottom_up`, `:top_down`, `:gcta`, and `:plink`.
When `min_samples` is defined, pruning stops if that number of samples remain.

- `:bottom_up` runs bottom-up greedy method, iteratively eliminating neighbors of sample with the lowest degree.
- `:top_down` runs top-down greedy method, iteratively eliminating sample with the highest degree.
- `:gcta` runs the method implemented in GCTA.
- `:plink` runs the method implemented in PLINK.

In general, `:bottom_up` greedy method performs the best (keeping the hightest number of subjects) and has a guaranteed approximamtion ratio.
With `min_samples` defined, `:top_down` greedy method keeps the least number of off-diagonal values.
"""
function kinship_pruning(grm::AbstractMatrix; method=:bottom_up, cutoff=0.125, min_samples=0)
    if method == :bottom_up
        r = _kinship_pruning_bottomup(grm; cutoff=cutoff, min_samples=min_samples)
    elseif method == :top_down
        r = _kinship_pruning_topdown(grm; cutoff=cutoff, min_samples=min_samples)
    elseif method == :gcta
        r = _kinship_pruning_gcta(grm; cutoff=cutoff)
    elseif method == :plink
        r = _kinship_pruning_plink(grm; cutoff=cutoff, min_samples=min_samples)
    else
        @error "Invalid value of argument method. Valid values are: :bottom_up, :top_down, :gcta, and :plink."
    end
    r
end

function _find_neighbors!(neighbors::Vector{<:Integer}, i::Integer, grm::AbstractMatrix, r::BitVector, degree::Integer; cutoff=0.125)
    m = size(grm, 1)
    n_neighbors = 0
    @inbounds for j in 1:m
        if r[j] && i != j && grm[i, j] > cutoff
            n_neighbors += 1
            neighbors[n_neighbors] = j
        end
    end
    @assert n_neighbors == degree
    neighbors
end

function _count_degrees(grm::AbstractMatrix; cutoff=0.125)
    m = size(grm, 1)
    degrees = zeros(Int, m)
    @inbounds for i in 1:m
        for j in 1:m
            if j == i; continue; end
            if grm[i, j] > cutoff
                degrees[i] += 1
            end
        end
    end
    degrees
end

function _eliminate_node!(i::Integer, grm::AbstractMatrix, r::BitVector, degrees::Vector{<:Integer}; cutoff=0.125)
    m = size(grm, 1)
    r[i] = false
    degrees[i] = 0
    @inbounds for j in 1:m
        if r[j] && grm[i, j] > cutoff
            degrees[j] -= 1
        end
    end
    nothing
end

function _kinship_pruning_gcta(grm::AbstractMatrix; cutoff=0.125)
    m = size(grm, 1)
    degrees = _count_degrees(grm; cutoff=cutoff)

    r = trues(m)
    @inbounds for i in 1:m, j in (i+1):m
        if j == i; continue; end
        if grm[i, j] > cutoff
            if degrees[i] > degrees[j]
                r[i] = false
            else
                r[j] = false
            end
        end
    end
    r
end

function _kinship_pruning_plink(grm::AbstractMatrix; cutoff=0.125, min_samples=0)
    m = size(grm, 1)
    degrees = _count_degrees(grm; cutoff=cutoff)
    r = trues(m)
    nbs = [0]
    while sum(degrees) > 0 && count(r) > min_samples
        if any(degrees .== 1)
            to_add = findfirst(x -> x == 1, degrees)
            to_eliminate = _find_neighbors!(nbs, to_add, grm, r, degrees[to_add]; cutoff=cutoff)
            _eliminate_node!(to_eliminate[1], grm, r, degrees; cutoff=cutoff)
        else
           to_eliminate = findfirst(x -> x == maximum(degrees), degrees)
            _eliminate_node!(to_eliminate, grm, r, degrees; cutoff=cutoff)
        end
    end
    if sum(degrees) != 0
        @warn "not all edges were removed: $(sum(degrees) ÷ 2) edges remain"
    end
    r
end

function _kinship_pruning_topdown(grm::AbstractMatrix; cutoff=0.125, min_samples=0)
    m = size(grm, 1)
    degrees = _count_degrees(grm; cutoff=cutoff)
    r = trues(m)
    while sum(degrees) > 0 && count(r) > min_samples
        to_eliminate = findfirst(x -> x == maximum(degrees), degrees)
        _eliminate_node!(to_eliminate, grm, r, degrees; cutoff=cutoff)
    end
    if sum(degrees) != 0
        @warn "not all edges were removed: $(sum(degrees) ÷ 2) edges remain"
    end
    r
end

function _kinship_pruning_bottomup(grm::AbstractMatrix; cutoff=0.125, min_samples=0)
    m = size(grm, 1)
    degrees = _count_degrees(grm; cutoff=cutoff)
    nbs = Array{Int}(undef, m)
    r = trues(m)
    while sum(degrees) > 0 && count(r) > min_samples
        minpositive = minimum(degrees[degrees .> 0])
        to_add = findfirst(x -> x == minpositive, degrees)
        to_eliminate = _find_neighbors!(nbs, to_add, grm, r, degrees[to_add]; cutoff=cutoff)
        @inbounds for i in 1:degrees[to_add]
            _eliminate_node!(to_eliminate[i], grm, r, degrees; cutoff=cutoff)
        end
    end
    if sum(degrees) != 0
        @warn "not all edges were removed: $(sum(degrees) ÷ 2) edges remain"
    end
    r
end


text:
"""
    get_id(x)

Get the IDENTIFIER name of a variable, possibly in the presence of
type declaration operators.
"""

function get_id(x::EXPR{BinarySyntaxOpCall})
    if x.args[2] isa EXPR{OPERATOR{ComparisonOp,Tokens.ISSUBTYPE,false}} || x.args[2] isa EXPR{OPERATOR{DeclarationOp,Tokens.DECLARATION,false}} || x.args[2] isa EXPR{OPERATOR{WhereOp,Tokens.WHERE,false}}
        return get_id(x.args[1])
    else
        return x
    end
end

get_id(x::EXPR{Curly}) = get_id(x.args[1])
get_id(x) = x

"""
    get_t(x)

Basic inference in the presence of type declarations.
"""
get_t(x) = :Any
function get_t(x::EXPR{BinarySyntaxOpCall})
    if x.args[2] isa EXPR{OPERATOR{DeclarationOp,Tokens.DECLARATION,false}}
        return Expr(x.args[3])
    else
        return :Any
    end
end

# NEEDS FIX
# function func_sig(x::EXPR)
#     name = x.args[1]
#     args = x.args[2:end]
#     if name isa EXPR && name.head == CURLY
#         params = name.args[2]
#         name = name.args[1]
#     end
#     if name isa EXPR && name.head isa OPERATOR{DotOp,Tokens.DOT}
#         mod = name.args[1]
#         name = name.args[2]
#     end
#     if name isa EXPR{Quotenode}
#         name = name.val
#     end
# end

"""
    _track_assignment(ps, x, val, defs = [])

When applied to the lhs of an assignment returns a vector of the
newly defined variables.
"""
function _track_assignment(ps::ParseState, x, val, defs = [])
    return defs
end
function _track_assignment(ps::ParseState, x::EXPR{IDENTIFIER}, val, defs = [])
    t = infer_t(val)
    push!(defs, Variable(Expr(x), t, val))
    return defs
end

function _track_assignment(ps::ParseState, x::EXPR{BinarySyntaxOpCall}, val, defs = [])
    if x.args[2] isa EXPR{OPERATOR{DeclarationOp,Tokens.DECLARATION,false}}
        t = Expr(x.args[3])
        push!(defs, Variable(Expr(get_id(x.args[1])), t, val))
    end
    return defs
end

function _track_assignment(ps::ParseState, x::EXPR{Curly}, val, defs = [])
    t = infer_t(val)
    push!(defs, Variable(Expr(get_id(x)), t, val))
    return defs
end

function _track_assignment(ps::ParseState, x::EXPR{TupleH}, val, defs = [])
    for a in x.args
        _track_assignment(ps, a, val, defs)
    end
    return defs
end

infer_t(x) = :Any
infer_t(x::EXPR{LITERAL{Tokens.INTEGER}}) = :Int
infer_t(x::EXPR{LITERAL{Tokens.FLOAT}}) = :Float64
infer_t(x::EXPR{LITERAL{Tokens.STRING}}) = :String
infer_t(x::EXPR{LITERAL{Tokens.TRIPLE_STRING}}) = :String
infer_t(x::EXPR{LITERAL{Tokens.CHAR}}) = :Char
infer_t(x::EXPR{LITERAL{Tokens.TRUE}}) = :Bool
infer_t(x::EXPR{LITERAL{Tokens.FALSE}}) = :Bool
infer_t(x::EXPR{LITERAL{Tokens.CMD}}) = :Cmd

infer_t(x::EXPR{Vect}) = :(Array{Any,1})
infer_t(x::EXPR{Vcat}) = :(Array{Any,N})
infer_t(x::EXPR{TypedVcat}) = :(Array{$(Expr(x.args[1])),N})
infer_t(x::EXPR{Hcat}) = :(Array{Any,2})
infer_t(x::EXPR{TypedHcat}) = :(Array{$(Expr(x.args[1])),2})
infer_t(x::EXPR{Quote}) = :Expr
infer_t(x::EXPR{StringH}) = :String
infer_t(x::EXPR{Quotenode}) = :QuoteNode

function get_symbols(x, offset = 0, symbols = []) end

function get_symbols(x::EXPR, offset = 0, symbols = [])
    for a in x.args
        for v in a.defs
            push!(symbols, (v, offset + (1:a.span)))
        end
        if contributes_scope(a)
            get_symbols(a, offset, symbols)
        end
        if a isa EXPR{ModuleH} || a isa EXPR{BareModule}
            m_scope = get_symbols(a.args[3])
            offset2 = offset + a.args[1].span + a.args[2].span
            for mv in m_scope
                push!(symbols, (Variable(Expr(:(.), a.defs[1].id, QuoteNode(mv[1].id)), mv[1].t, mv[1].val), mv[2] + offset2))
            end
        end
        offset += a.span
    end
    return symbols
end

"""
    contributes_scope(x)
Checks whether the body of `x` is included in the toplevel namespace.
"""
contributes_scope(x) = false
contributes_scope(x::EXPR{FileH}) = true
contributes_scope(x::EXPR{Begin}) = true
contributes_scope(x::EXPR{Block}) = true
contributes_scope(x::EXPR{Const}) = true
contributes_scope(x::EXPR{Global}) = true
contributes_scope(x::EXPR{Local}) = true
contributes_scope(x::EXPR{If}) = true
contributes_scope(x::EXPR{MacroCall}) = true
contributes_scope(x::EXPR{TopLevel}) = true

ext: jl


text:
module Periods

import Base: *, /, get, isapprox, show

export TimeUnit, Second, Minute, Hour, Day, Year, Century,
    J2000, J1950, MJD,
    seconds, minutes, hours, days, years, centuries,
    Period, *, /, get

include("constants.jl")

"""
All time units are subtypes of the abstract type `TimeUnit`.
The following time units are defined:

* `Second`
* `Minute`
* `Hour`
* `Day`
* `Year`
* `Century`
"""
abstract type TimeUnit end

struct Second <: TimeUnit end
struct Minute <: TimeUnit end
struct Hour <: TimeUnit end
struct Day <: TimeUnit end
struct Year <: TimeUnit end
struct Century <: TimeUnit end

const seconds = Second()
const minutes = Minute()
const hours = Hour()
const days = Day()
const years = Year()
const centuries = Century()

struct Period{U<:TimeUnit,T<:Number}
    unit::Type{U}
    Δt::T
    Period{U}(Δt::T) where {U<:TimeUnit,T<:Number} = new{U,T}(U, Δt)
end

get(p::Period) = p.Δt

show(io::IO, p::Period{Second}) = print(io, "$(get(p)) seconds")
show(io::IO, p::Period{Minute}) = print(io, "$(get(p)) minutes")
show(io::IO, p::Period{Hour}) = print(io, "$(get(p)) hours")
show(io::IO, p::Period{Day}) = print(io, "$(get(p)) days")
show(io::IO, p::Period{Year}) = print(io, "$(get(p)) years")
show(io::IO, p::Period{Century}) = print(io, "$(get(p)) centuries")

(*)(Δt::T, ::U) where {T<:Number, U<:TimeUnit} = Period{U}(Δt)

(*)(x::T, p::Period) where {T<:Number} = Period{p.unit}(p.Δt * x)
(*)(p::Period, x::T) where {T<:Number} = Period{p.unit}(p.Δt * x)
(/)(x::T, p::Period) where {T<:Number} = Period{p.unit}(x / p.Δ)
(/)(p::Period, x::T) where {T<:Number} = Period{p.unit}(p.Δt / x)

isapprox(p1::Period{U}, p2::Period{U}) where {U<:TimeUnit} = get(p1) ≈ get(p2)

(::Second)(p::Period{Second})  = p
(::Second)(p::Period{Minute})  = Period{Second}(p.Δt * SECONDS_PER_MINUTE)
(::Second)(p::Period{Hour})    = Period{Second}(p.Δt * SECONDS_PER_HOUR)
(::Second)(p::Period{Day})     = Period{Second}(p.Δt * SECONDS_PER_DAY)
(::Second)(p::Period{Year})    = Period{Second}(p.Δt * SECONDS_PER_YEAR)
(::Second)(p::Period{Century}) = Period{Second}(p.Δt * SECONDS_PER_CENTURY)

(::Minute)(p::Period{Second})  = Period{Minute}(p.Δt / SECONDS_PER_MINUTE)
(::Minute)(p::Period{Minute})  = p
(::Minute)(p::Period{Hour})    = Period{Minute}(p.Δt * MINUTES_PER_HOUR)
(::Minute)(p::Period{Day})     = Period{Minute}(p.Δt * MINUTES_PER_DAY)
(::Minute)(p::Period{Year})    = Period{Minute}(p.Δt * MINUTES_PER_YEAR)
(::Minute)(p::Period{Century}) = Period{Minute}(p.Δt * MINUTES_PER_CENTURY)

(::Hour)(p::Period{Second})  = Period{Hour}(p.Δt / SECONDS_PER_HOUR)
(::Hour)(p::Period{Minute})  = Period{Hour}(p.Δt / MINUTES_PER_HOUR)
(::Hour)(p::Period{Hour})    = p
(::Hour)(p::Period{Day})     = Period{Hour}(p.Δt * HOURS_PER_DAY)
(::Hour)(p::Period{Year})    = Period{Hour}(p.Δt * HOURS_PER_YEAR)
(::Hour)(p::Period{Century}) = Period{Hour}(p.Δt * HOURS_PER_CENTURY)

(::Day)(p::Period{Second})  = Period{Day}(p.Δt / SECONDS_PER_DAY)
(::Day)(p::Period{Minute})  = Period{Day}(p.Δt / MINUTES_PER_DAY)
(::Day)(p::Period{Hour})    = Period{Day}(p.Δt / HOURS_PER_DAY)
(::Day)(p::Period{Day})     = p
(::Day)(p::Period{Year})    = Period{Day}(p.Δt * DAYS_PER_YEAR)
(::Day)(p::Period{Century}) = Period{Day}(p.Δt * DAYS_PER_CENTURY)

(::Year)(p::Period{Second})  = Period{Year}(p.Δt / SECONDS_PER_YEAR)
(::Year)(p::Period{Minute})  = Period{Year}(p.Δt / MINUTES_PER_YEAR)
(::Year)(p::Period{Hour})    = Period{Year}(p.Δt / HOURS_PER_YEAR)
(::Year)(p::Period{Day})     = Period{Year}(p.Δt / DAYS_PER_YEAR)
(::Year)(p::Period{Year})    = p
(::Year)(p::Period{Century}) = Period{Year}(p.Δt * YEARS_PER_CENTURY)

(::Century)(p::Period{Second})  = Period{Century}(p.Δt / SECONDS_PER_CENTURY)
(::Century)(p::Period{Minute})  = Period{Century}(p.Δt / MINUTES_PER_CENTURY)
(::Century)(p::Period{Hour})    = Period{Century}(p.Δt / HOURS_PER_CENTURY)
(::Century)(p::Period{Day})     = Period{Century}(p.Δt / DAYS_PER_CENTURY)
(::Century)(p::Period{Year})    = Period{Century}(p.Δt / YEARS_PER_CENTURY)
(::Century)(p::Period{Century}) = p

end


text:
using Turing, ContinuousBenchmarks
using Mamba: describe

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "school8-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "school8.model.jl"))

data = deepcopy(schools8data[1])
delete!(data, "tau")

# chn = sample(school8(data=data), HMC(2000, 0.75, 5))

_data = map(ind->data[ind], ["J", "y", "sigma"])
@tbenchmark(HMC(20, 0.75, 5), school8, (_data...))
bench_res = @tbenchmark(HMC(2000, 0.75, 5), school8, (_data...))
# bench_res[4].names = ["phi[1]", "phi[2]", "phi[3]", "phi[4]"]
logd = build_logd("School 8", bench_res...)

# describe(chn)

#include(joinpath(ContinuousBenchmarks.BENCH_DIR, "school8-stan.run.jl"))

#logd["stan"] = stan_d
#logd["time_stan"] = get_stan_time("schools8")

print_log(logd)
send_log(logd)


text:
# Multi dimensional array broadcast
# this kernel can handle all arrays with different dimensions and broadcasting multiple dimensions
# x and y expected to be an array
#in stride_x and stride_y, values corresponding to broadcasting dim should be zero
#in the one that is being broadcasted
# N_z is the size of the z array
# dimlen_z is the number of the dimensions in z

#performance limitations
#stride values and coordinates kept as arrays, causing them to replaced in heap,shared memory
# (slow: http://stackoverflow.com/a/13485322)
# TODO-enis understand heap memory and dynamic allocation for possible bug cause with arrays of >100 dims
# http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#heap-memory-allocation

fp = open("cuda17.cu","w")
#using Knet: broadcast_ops

function cuda17src(f, j=f, ex="$f(xi,yi)")
  sprint() do s

    for (T,F) in [("float","$(f)_32"),("double","$(f)_64")]
        print(s,

"""
__global__ void _$(F)_17($T *x,$T *y, $T *z,
                                int *stride_x, int *stride_y,
                                int *stride_z,int N_z,int dimlen_z) {

    int index_z = threadIdx.x + (blockIdx.x * blockDim.x);
    int index_x,index_y,coords;
    int temp_index;

    while (index_z < N_z) {
        temp_index = index_z;

        // TODO replace (i/n) == (i>>log2(n)) also %
        index_x =0;
        index_y = 0;
        for (int i=dimlen_z-1; i>0; i--)
        {
            coords = temp_index / stride_z[i];
            index_x+= stride_x[i]*coords;
            index_y+= stride_y[i]*coords;
            temp_index = temp_index % stride_z[i];
        }
        index_x+= temp_index;
        index_y+= temp_index;

        $T xi = x[index_x];
        $T yi = y[index_y];
        z[index_z]=$ex;;
        index_z+=(blockDim.x * gridDim.x);
    }
}

extern "C" {
  $DLLEXPORT void $(F)_17($T *x,$T *y,$T *z, int *stride_x, int *stride_y,int *stride_z,int N_z,int dimlen_z) {

    _$(F)_17<<<256,256>>>(x,y,z,stride_x,stride_y,stride_z,N_z,dimlen_z);

  }
}
""")
    end
  end
end

for a in broadcast_ops
    if !isa(a,Tuple); a=(a,); end
    print(fp,cuda17src(a...))
end

close(fp)


text:
# Forward- and Reverse-Mode Differentiation
# =========================================

include("utils.jl");

# Differentiation tools are frequently described as implementing "forward mode"
# or "reverse mode" AD. This distinction was briefly covered in the intro
# notebook, but here we'll go into more detail. We'll start with an intuition
# for what the distinction *means* in terms of the differentiation process; then
# we'll discuss why it's an important consideration in practice.

# Consider a simple mathematical expression:

y = :(sin(x^2) * 5)

# Written as a Wengert list:

Wengert(y)

# The ability to take derivatives mechanically relies on two things: Firstly, we
# know derivatives for each basic function in our program (e.g.
# $\frac{dy_2}{dy_1}=cos(y_1)$). Secondly, we have a rule of composition called the
# *chain rule* which lets us compose these basic derivatives together.
#
# $$
# \frac{dy}{dx} = \frac{dy_1}{dx} \times \frac{dy_2}{dy_1} \times \frac{dy}{dy_2}
# $$
#
# More specifically:
#
# $$
# \begin{align}
# \frac{dy}{dx} &= 2x \times cos(y_1) \times 5 \\
#               &= 2x \times cos(x^2) \times 5
# \end{align}
# $$
#
# The forward/reverse distinction basically amounts to: given that we do
# multiplications one at a time, do we evaluate $\frac{dy_1}{dx} \times
# \frac{dy_2}{dy_1}$ first, or $\frac{dy_2}{dy_1} \times \frac{dy}{dy_2}$? (This seems
# like a pointless question right now, given that either gets us the same
# results, but bear with me.)
#
# It's easier to see the distinction if we think algorithmically. Given some
# enormous Wengert list with $n$ instructions, we have two ways to differentiate
# it:
#
# **(1)**: Start with the known quantity $\frac{dy_0}{dx} = \frac{dx}{dx} = 1$
# at the beginning of the list. Look up the derivative for the next instruction
# $\frac{dy_{i+1}}{dy_i}$ and multiply out the top, getting $\frac{dy_1}{dx}$,
# $\frac{dy_2}{dx}$, ... $\frac{dy_{n-1}}{dx}$, $\frac{dy}{dx}$. Because we
# walked forward over the Wengert list, this is called *forward mode*. Each
# intermediate derivative $\frac{dy_i}{dx}$ is known as a *perturbation*.
#
# **(2)**: Start with the known quantity $\frac{dy}{dy_n} = \frac{dy}{dy} = 1$
# at the end of the list. Look up the derivative for the previous instruction
# $\frac{dy_i}{dy_{i-1}}$ and multiply out the bottom, getting
# $\frac{dy}{dy_n}$, $\frac{dy}{dy_{n-1}}$, ... $\frac{dy}{dy_1}$,
# $\frac{dy}{dx}$. Because we walked in reverse over the Wengert list, this is
# called *reverse mode*. Each intermediate derivative $\frac{dy}{dy_i}$ is known
# as a *sensitivity*.
#
# This all seems very academic, so we need to explain why it might make a
# difference to performance. I'll give two related explanations: dealing with
# mulitple variables, and working with vectors rather than scalars.

# Explanation 1: Multiple Variables
# ---------------------------------
#
# So far we have dealt only with simple functions that take a number, and return
# a number. But more generally we'll deal with functions that take, or produce,
# multiple numbers of interest.
#
# For example what if we have a function that returns *two* numbers, and we want
# derivatives for both? Do we have to do the whole thing twice over?

y = quote
  y2 = sin(x^2)
  y3 = y2 * 5
end

# Let's say we want both of the derivatives $\frac{dy_2}{dx}$ and
# $\frac{dy_3}{dx}$. You can probably see where this is going now; the Wengert
# list representation of this expression has not actually changed!

Wengert(y)

# Now, we discussed that when doing forward mode differentiation, we actually
# calculate *every* intermediate derivative $\frac{dy_i}{dx}$, which means we get
# $\frac{dy_2}{dx}$ for free. This property goes all the way back to our
# original, recursive formulation of differentiation, which calculated the
# derivatives of a complex expression by combining the derivatives of simpler
# ones.

derive(Wengert(y), :x)

# In our output, $y_7 = \frac{dy_2}{dx}$ and $y_8 = \frac{dy_3}{dx}$.
#
# Let's consider the opposite situation, a function of two variables $a$ and
# $b$, where we'd like to get $\frac{dy}{da}$ and $\frac{dy}{db}$.

y = :(sin(a) * b)
#-
Wengert(y)

# This one is a bit tougher. We can start the forward-mode differentiation
# process with $\frac{da}{da} = 1$ or with $\frac{db}{db} = 1$, but if we want
# both we'll have to go through the entire multiplying-out process twice.
#
# But both variables ultimately end up at the same place, $y$, and we know that
# $\frac{dy}{dy} = 1$. Aha, so perhaps we can use reverse mode for this
# instead!
#
# Exactly opposite to forward mode, reverse mode gives us every intermediate
# gradient $\frac{dy_i}{dy}$ for free, ultimately leading back in the inputs
# $\frac{da}{dy}$ and $\frac{db}{dy}$.
#
# It's easy to see, then, why reverse-mode differentiation – or backpropagation
# – is so effective for machine learning. In general we have a large computation
# with millions of parameters, yet only a single scalar loss to optimise. We can
# get gradients even for these millions of inputs in a single pass, enabling ML
# to scale to complex tasks like image and voice recognition.

# Explanation 2: Vector Calculus
# ------------------------------
#
# So far we have dealt only with simple functions that take a number, and return
# a number. But more generally we'll deal with functions that take, or produce,
# *vectors* containing multiple numbers of interest.
#
# It's useful to consider how our idea of differentiation works when we have
# vectors. For example, a function that takes a vector of length $2$ to another
# vector of length $2$:

f(x) = [x[1] * x[2], cos(x[1])]

x = [2, 3]
y = f(x)

# We now need to talk about what we mean by $\frac{d}{dx}f(x)$, given that we
# can't apply the usual limit rule. What we *can* do is take the derivative of
# any scalar *element* of $y$ with respect to any element of $x$. (We'll use
# subscripts $x_n$ to refer to the $n^{th}$ index of $x$.) For example:
#
# $$
# \begin{align}
# \frac{dy_1}{dx_1} &= \frac{d}{dx_1} x_1 \times x_2 = x_2 \\
# \frac{dy_1}{dx_2} &= \frac{d}{dx_2} x_1 \times x_2 = x_1 \\
# \frac{dy_2}{dx_1} &= \frac{d}{dx_1} \cos(x_1) = \sin(x_1) \\
# \frac{dy_2}{dx_2} &= \frac{d}{dx_2} \cos(x_1) = 0 \\
# \end{align}
# $$
#
# It's a little easier if we organise all of these derivatives into a matrix.
#
# $$
# J_{ij} = \frac{dy_i}{dx_j}
# $$
#
# This $2\times2$ matrix is called the *Jacobian*, and in general it's what we mean by
# $\frac{dy}{dx}$. (The Jacobian for a scalar function like $y = \sin(x)$ only
# has one element, so it's consistent with our current idea of the derivative
# $\frac{dy}{dx}$.) The key point here is that the Jacobian is a potentially
# large object: it has a size `length(y) * length(x)`. Now, we discussed that
# the distinction between forward and reverse mode is whether we propagate
# $\frac{dy_i}{dx}$ or $\frac{dy}{dy_i}$, which can have a size of either
# `length(y_i) * length(x)` or `length(y) * length(y_i)`.
#
# It should be clear, then, what mode is better if we have a gazillion inputs
# and one output. In forward mode we need to carry around a gazillion
# "perturbations" *for each* element of $y_i$, whereas in reverse we only need a
# gradient of the same size of $x$. And vice versa.


text:
using Latexify
using Test

@test latexify(((1.0, 2), (3, 4)); env=:align) == replace(
raw"\begin{align}
1.0 =& 3 \\
2 =& 4
\end{align}
", "\r\n"=>"\n")

@test latexify(((1.0, 2), (3, 4)); separator = [" =& ", " ∈& "], env = :align) == replace(
raw"\begin{align}
1.0 =& 3 \\
2 ∈& 4
\end{align}
", "\r\n"=>"\n")

lhs = [:a, :b]
rhs = [1, 2]

@test latexify(lhs, rhs; env = :aligned) == replace(
raw"\[
\begin{aligned}
a =& 1 \\
b =& 2
\end{aligned}
\]
", "\r\n"=>"\n")

# @test_throws MethodError latexify(rn; bad_kwarg="should error")

# Latexify.@generate_test latexify(["a=1"], env=:align)
@test latexify(["a=1"], env = :align) == replace(
raw"\begin{align}
a =& 1
\end{align}
", "\r\n"=>"\n")


text:
using Test
using PooledArrays

let a = rand(10), b = rand(10,10), c = rand(1:10, 1000)
    @test PooledArray(a) == a
    @test PooledArray(b) == b
    pc = PooledArray(c)
    @test pc == c
    @test convert(Array, pc) == c
    #@test issorted(pc.pool)

    @test copy(pc) == pc

    @test isa(similar(pc), typeof(pc))

    @test reverse(pc) == reverse(c)

    @test pc[34:250] == c[34:250]

    @test sort(pc) == sort(c)
    @test sortperm(pc) == sortperm(c)

    push!(pc, -10)
    push!(c, -10)
    @test pc == c
    @test maximum(values(pc.pool)) == 11
    @test maximum(keys(pc.revpool)) == 11

    #@test issorted(pc.pool)

    @test map(identity, pc) == pc
    @test map(x->2x, pc) == map(x->2x, c)

    # case where the outputs are one-to-many
    pa = PooledArray([1,2,3,4,5,6])
    @test map(isodd, pa) == [true,false,true,false,true,false]

    px = PooledArray(rand(128))
    py = PooledArray(rand(200))

    @test isa(vcat(px, py), PooledArray{Float64, UInt16})
    @test vcat(px, Array(py)) == vcat(px, py)

    px2 = copy(px)
    resize!(px2, 100)
    @test px2 == px[1:100]

    @test findall(PooledArray([true,false,true])) == [1,3]
end


text:
module TestExamples2DAcousticPerturbation

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "Acoustic Perturbation" begin
  @trixi_testset "elixir_acoustics_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_convergence.jl"),
      l2 = [0.0019921138796370834, 0.002090394698052287, 0.0006091925854593805, 0.0, 0.0, 0.0, 0.0],
      linf = [0.00769282588065634, 0.008276649669227254, 0.004196479023954813, 0.0, 0.0, 0.0, 0.0])
  end

  @trixi_testset "elixir_acoustics_gauss.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_gauss.jl"),
      l2 = [0.08005276517890283, 0.08005276517890268, 0.4187202920734123, 0.0, 0.0, 0.0, 0.0],
      linf = [0.17261097190220992, 0.17261097190220973, 1.13601894068238, 0.0, 0.0, 0.0, 0.0])
  end

  @trixi_testset "elixir_acoustics_gaussian_source.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_gaussian_source.jl"),
      l2 = [0.004296394903650806, 0.004241280404758938, 0.006269684906035964, 0.0, 0.0, 0.0, 0.0],
      linf = [0.03970270697049378, 0.04151096349298151, 0.0640019829058819, 0.0, 0.0, 0.0, 0.0])
  end

  @trixi_testset "elixir_acoustics_gauss_wall.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_gauss_wall.jl"),
      l2 = [0.019419398248465843, 0.019510701017551826, 0.04818246051887614,
            7.382060834820337e-17, 0.0, 1.4764121669640674e-16, 1.4764121669640674e-16],
      linf = [0.18193631937316496, 0.1877464607867628, 1.0355388011792845,
              2.220446049250313e-16, 0.0, 4.440892098500626e-16, 4.440892098500626e-16])
  end

  @trixi_testset "elixir_acoustics_monopole.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_acoustics_monopole.jl"),
      l2 = [0.006816790293009947, 0.0065068948357351625, 0.008724512056168938,
            0.0009894398191644543, 0.0, 7.144325530679576e-17, 7.144325530679576e-17],
      linf = [1.000633375007386, 0.5599788929862504, 0.5738432957070382,
              0.015590137026938428, 0.0, 2.220446049250313e-16, 2.220446049250313e-16],
      maxiters=50)
  end
end

end # module


text:
abstract type BoundaryConditions
end

struct PeriodicBoundaryConditions{cType <: Real} <: BoundaryConditions
    boundary::SVector{6,cType}
end

PeriodicBoundaryConditions(L::Real) = PeriodicBoundaryConditions(SVector(0, L, 0, L, 0, L))

function Base.iterate(pbc::PeriodicBoundaryConditions,state=1)
  state > length(pbc.boundary) && return nothing
  pbc.boundary[state], state + 1
end

function Base.getindex(pbc::PeriodicBoundaryConditions, i::Integer)
    1 <= i <= length(pbc.boundary) || throw(BoundsError(pbc, i))
    pbc.boundary[i]
end

struct InfiniteBox{cType <: Real} <: BoundaryConditions
    boundary::SVector{6,<:cType}
end

InfiniteBox() = InfiniteBox(SVector(-Inf, Inf, -Inf, Inf, -Inf, Inf))

struct CubicPeriodicBoundaryConditions{cType <: Real} <: BoundaryConditions
    L::cType
end

function get_interparticle_distance(ri, rj, pbc::PeriodicBoundaryConditions)
    rij = ri - rj
    x, y, z = rij
    while x  < pbc[1]   x += pbc[2]-pbc[1] end
    while x >= pbc[2]   x -= pbc[2]-pbc[1] end
    while y  < pbc[3]   y += pbc[4]-pbc[3] end
    while y >= pbc[4]   y -= pbc[4]-pbc[3] end
    while z  < pbc[5]   z += pbc[6]-pbc[5] end
    while z >= pbc[6]   z -= pbc[6]-pbc[5] end
    rij = @SVector [x, y, z]
    r2 = rij[1]^2 + rij[2]^2 + rij[3]^2
    r = sqrt(r2)
    return (rij, r, r2)
end

function get_interparticle_distance(ri, rj, bc::CubicPeriodicBoundaryConditions)
    rij = ri - rj
    x, y, z = rij
    size = bc.L
    radius = 0.5 * size
    while x >= radius    x -= size end
    while x < -radius    x += size end
    while y >= radius    y -= size end
    while y < -radius    y += size end
    while z >= radius    z -= size end
    while z < -radius    z += size end
    rij = @SVector [x, y, z]
    r2 = rij[1]^2 + rij[2]^2 + rij[3]^2
    r = sqrt(r2)
    return (rij, r, r2)
end

function get_interparticle_distance(ri, rj, ::BoundaryConditions)
    rij = ri - rj
    r2 = rij[1]^2 + rij[2]^2 + rij[3]^2
    r = sqrt(r2)
    (rij, r, r2)
end


text:
abstract type AbstractObject end

struct ExampleObject1 <: AbstractObject end
const EXAMPLE_OBJECT_1 = ExampleObject1()

struct ExampleObject2 <: AbstractObject end
const EXAMPLE_OBJECT_2 = ExampleObject2()

struct ExampleObject3 <: AbstractObject end
const EXAMPLE_OBJECT_3 = ExampleObject3()


text:
struct QuickWrapper{E<:AbstractEnv, D<:NamedTuple} <: AbstractWrapper
    env::E
    data::D
end

"""
    QuickWrapper(env; kwargs...)

Create a wrapper to override specific behavior of the environment with keywords.

Each keyword argument corresponds to a CommonRLInterface function to be overridden. The keyword arguments can either be static objects or functions. If a keyword argument is a function, the arguments will be the wrapped environment and any other arguments. `provided` is automatically handled.

# Examples

Override the action space statically:
```julia
w = QuickWrapper(env; actions=[-1, 1])
observe(w) # returns the observation from env
actions(w) # returns [-1, 1]
```

Override the `act!` function to return the reward squared:
```julia
w = QuickWrapper(env; act! = (env, a) -> act!(env, a).^2)
act!(w, a) # returns the squared reward for taking action a in env
```
"""
QuickWrapper(e; kwargs...) = QuickWrapper(e, values(kwargs))

wrapped_env(w::QuickWrapper) = w.env

macro quick_forward(f)
    quote
        function $f(w::QuickWrapper, args...)
            if haskey(w.data, nameof($f))
                _call(w.data[nameof($f)], w.env, args...)
            else
                $f(w.env, args...)
            end
        end

        function CommonRLInterface.provided(::typeof($f), w::QuickWrapper, args...)
            if haskey(w.data, nameof($f))
                return true
            else
                return provided($f, w.env, args...)
            end
        end

        function CommonRLInterface.provided(::typeof($f), TT::Type{<:Tuple{QuickWrapper{E,D}, Vararg}}) where {E,D}
            if hasfield(D, nameof($f))
                return true
            else
                return provided($f, Tuple{E, TT.parameters[2:end]...})
            end
        end
    end
end

_call(f::Function, args...) = f(args...)
_call(other, args...) = other

@quick_forward CommonRLInterface.reset!
@quick_forward CommonRLInterface.actions
@quick_forward CommonRLInterface.observe
@quick_forward CommonRLInterface.act!
@quick_forward CommonRLInterface.terminated

@quick_forward CommonRLInterface.render
@quick_forward CommonRLInterface.state
@quick_forward CommonRLInterface.setstate!
@quick_forward CommonRLInterface.valid_actions
@quick_forward CommonRLInterface.valid_action_mask
@quick_forward CommonRLInterface.observations
@quick_forward CommonRLInterface.player

function CommonRLInterface.clone(w::QuickWrapper, args...)
    if haskey(w.data, :clone)
        QuickWrapper(_call(w.data[:clone], w.env), w.data)
    else
        QuickWrapper(clone(w.env), w.data)
    end
end

function CommonRLInterface.provided(f::typeof(CommonRLInterface.clone), w::QuickWrapper, args...)
    if haskey(w.data, :clone)
        return true
    else
        return provided(f, w.env, args...)
    end
end

function CommonRLInterface.provided(f::typeof(CommonRLInterface.clone), TT::Type{<:Tuple{QuickWrapper{E,D}, Vararg}}) where {E,D}
    if hasfield(D, :clone)
        return true
    else
        return provided(f, Tuple{E, TT.parameters[2:end]...})
    end
end


text:
"""
    PoissonProcess{T<:Real} <: AbstractMarkovChain

##### Fields
- `rate::T` : the rate; i.e., the parameter "lambda".
"""
struct PoissonProcess{T<:Real} <: AbstractMarkovChain
    rate::T
end

"""
    Base.show(io::IO, pp::PoissonProcess) -> nothing

Show a `PoissonProcess` in a human-friendly manner.
"""
function Base.show(io::IO, pp::PoissonProcess)
    println(io, typeof(pp))
    println(io, " Rate (i.e. lambda):")
    println(io, pp.rate)
    return nothing
end

"""
    Base.eltype(pp::PoissonProcess) -> DataType

Return the type associated with `pp`'s `rate` parameter.
"""
function Base.eltype(pp::PoissonProcess)
    return typeof(pp.rate)
end

"""
    nexpectedarrivals(pp::PoissonProcess, intervallength::Real) -> Int

Return the expected number of arrivals `E[N(intervallength)]` over the time interval
`[0, intervallength]` for `pp`.

##### Complexity
-Time complexity: `O(1)`.
"""
function nexpectedarrivals(pp::PoissonProcess, intervallength::Real)
    return pp.rate*intervallength
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
@inline get_transformation(R, t) = SMatrix{3, 4, Float64, 12}(R..., t...)

"""
For polynom, with coefficients as matrices, compute traces of the coefficients.
"""
function trace(p::SMatrix{3, 3, T, 9})::T where T
    mapreduce(mi -> tr(coefficient.(p, mi)) * mi, +, monomials(p[1]))
end

function null_space!(F, p1, p2)
    @inbounds @simd for i in 1:length(p1)
        p1i, p2i = p1[i], p2[i]
        F[i, 1] = p2i[1] * p1i[1] # x2 * x1
        F[i, 2] = p2i[1] * p1i[2] # x2 * y1
        F[i, 3] = p2i[1]          # x2

        F[i, 4] = p2i[2] * p1i[1] # y2 * x1
        F[i, 5] = p2i[2] * p1i[2] # y2 * y1
        F[i, 6] = p2i[2]          # y2

        F[i, 7] = p1i[1]          # x1
        F[i, 8] = p1i[2]          # y1
        F[i, 9] = 1.0
    end
    svd(F; full=true).V
end

@inbounds function to_polynoms(v1, v2, z)
    p1 = -v2[1] * z^3 + (v1[1] - v2[2]) * z^2 + (v1[2] - v2[3]) * z + v1[3]
    p2 = -v2[4] * z^3 + (v1[4] - v2[5]) * z^2 + (v1[5] - v2[6]) * z + v1[6]
    p3 = -v2[7] * z^4 + (v1[7] - v2[8]) * z^3 + (v1[8] - v2[9]) * z^2 + (v1[9] - v2[10]) * z + v1[10]
    p1, p2, p3
end

function compute_rref!(M, E1, E2, E3, E4)
    @polyvar x y z

    # One equation from rank constraint.
    @inbounds c11 = E1[1, 1] * x + E2[1, 1] * y + E3[1, 1] * z + E4[1, 1]
    @inbounds c12 = E1[1, 2] * x + E2[1, 2] * y + E3[1, 2] * z + E4[1, 2]
    @inbounds c13 = E1[1, 3] * x + E2[1, 3] * y + E3[1, 3] * z + E4[1, 3]

    @inbounds c21 = E1[2, 1] * x + E2[2, 1] * y + E3[2, 1] * z + E4[2, 1]
    @inbounds c22 = E1[2, 2] * x + E2[2, 2] * y + E3[2, 2] * z + E4[2, 2]
    @inbounds c23 = E1[2, 3] * x + E2[2, 3] * y + E3[2, 3] * z + E4[2, 3]

    @inbounds c31 = E1[3, 1] * x + E2[3, 1] * y + E3[3, 1] * z + E4[3, 1]
    @inbounds c32 = E1[3, 2] * x + E2[3, 2] * y + E3[3, 2] * z + E4[3, 2]
    @inbounds c33 = E1[3, 3] * x + E2[3, 3] * y + E3[3, 3] * z + E4[3, 3]

    # row1 - vector containing 20 coefficients for the first equation.
    row1 =
        c11 * c22 * c33 + c12 * c23 * c31 +
        c13 * c21 * c32 - c13 * c22 * c31 -
        c12 * c21 * c33 - c11 * c23 * c32

    # 9 equations from trace constraint.
    e1 = E1 .* x + E2 .* y + E3 .* z + E4
    e2 = E1' .* x + E2' .* y + E3' .* z + E4'

    mat_part = (e1 * e2) * e1
    trace_part = trace(e1 * e2) * e1
    row33 = mat_part .- 0.5 .* trace_part

    @inbounds monoms = monomials(row33[1])
    @inbounds @simd for i in 1:length(monoms)
        M[1:9, i] .= reshape(coefficient.(row33, monoms[i]), 9)
    end
    @inbounds M[10, :] .= coefficients(row1)

    order = UInt8[1,7,2,4,3,11,8,14,5,12,6,13,17,9,15,18,10,16,19,20]
    Base.permutecols!!(M, order)
    rref!(M)
end


text:
using Faker
using Base.Test

@test  ismatch(r"[0-9]", Faker.random_digit())
@test  ismatch(r"([0-9])|([0-9]{2})|([0-9]{3})|([0-9]{4})", Faker.random_int())
@test  ismatch(r"[0-9]", Faker.random_digit() )
@test  ismatch(r"[0-9]", Faker.random_digit_not_null() )
#print(Faker.random_digit_or_empty() )
#print(Faker.random_digit_not_null_or_empty() )
@test  ismatch(r"[0-9]", Faker.random_number() )
@test  ismatch(r"[0-9]", Faker.randomize_nb_elements() )
@test  ismatch(r"[a-zA-Z]", Faker.random_letter() )
@test  ismatch(r"[a-zA-Z]", Faker.random_element() )
@test  ismatch(r"[0-9]{2}[a-zA-Z]{2}", Faker.bothify("##??") )

#Faker.state_abbr()
#Faker.latitude()
#Faker.street_name()
#Faker.address()
#Faker.street_address()
#Faker.postcode()
#Faker.longitude()
#Faker.country()
#Faker.geo_coordinate("None", 0.001)
#Faker.secondary_address()
#Faker.city_prefix()
#Faker.city_suffix()
#Faker.building_number()
#Faker.city()
#Faker.state()

#Faker.company()
#Faker.company_suffix()
#Faker.company_prefix()
#Faker.catch_phrase()
#Faker.bs()

#Faker.date_time_ad()
#Faker.months()
#Faker.am_pm()
#Faker.iso8601()
#Faker.date_time_this_century()
#Faker.date_time()
#Faker.month_name()
#Faker.date_time_this_month()
#Faker.date_time_this_decade()
#Faker.day_of_week()
#Faker.day_of_month()
#Faker.time("%H:%M:%S")
#Faker.date_time_between("-30y", "now")
#Faker.unix_time()
#Faker.date_time_this_year()
#Faker.timezone()
#Faker.century()
#Faker.date("%Y-%m-%d")
#Faker.year()

#Faker.email()
#Faker.free_email()
#Faker.safe_email()
#Faker.user_name()
#Faker.domain_name()
#Faker.domain_word()
#Faker.domain_suffix()
#Faker.ipv4()
#Faker.ipv6()
#Faker.uri()
#Faker.url()

#Faker.text()
#Faker.sentence()
#Faker.word()
#Faker.paragraphs()
#Faker.words()
#Faker.paragraph()
#Faker.sentences()

#Faker.phone_number())

#Faker.password()
#Faker.locale()
#Faker.md5()
#Faker.sha1()
#Faker.null_boolean()
#Faker.sha256()
#Faker.country_code()
#Faker.language_code()
#Faker.boolean()

#colores
#Faker.color_name()
#Faker.safe_color_name()
#Faker.hex_color()
#Faker.safe_hex_color()
#Faker.rgb_color_list()
#Faker.rgb_color()
#Faker.rgb_css_color()

#CREDITCARD
#Faker.credit_card_provider()
#Faker.credit_card_number()
#Faker.credit_card_full()
#Faker.credit_card_security_code()
#Faker.credit_card_type()
#Faker.generate_number("52",16)

#file
#Faker.mime_type()
#Faker.file_extension()
#Faker.file_name()

#person
#Faker.first_name()
#Faker.last_name()
#Faker.prefixe()

#profile
#Faker.simple_profile()
#Faker.profile()

#job
#Faker.job()

#useragent
#Faker.mac_processor()
#Faker.linux_processor()
#Faker.user_agent()
#Faker.windows_platform_token()
#Faker.linux_platform_token()
#Faker.mac_platform_token()
#Faker.chrome()
#Faker.firefox()
#Faker.safari()
#Faker.opera()
#Faker.internet_explorer()

# write your own tests here
#@test 1 == 1


text:
using AQFED, Test
using StatsBase
import AQFED.American:
    AndersenLakeRepresentation,
    AndersenLakeNRepresentation,
    priceAmerican,
    americanBoundaryPutQDP
import AQFED.TermStructure: ConstantBlackModel

@testset "ALNegative" begin
    spot = 80.0
    strike = 100.0
    r = -0.005
    q = -0.05
    tte = 1.0
    vol = 0.2

    model = ConstantBlackModel(vol, r, q)
    pricer = AndersenLakeNRepresentation(model, tte, 1e-8, 5, 16, 21, 201)
    price = priceAmerican(pricer, strike, spot)
    println(price)
    @test isapprox(20.030946627004035, price, atol = 1e-4)
end

@testset "ALQ0" begin
    strike = 100.0
    spot = 80.0
    vol = 0.3
    tte = 0.75
    q = 0.0
    r = 0.04
    model = ConstantBlackModel(vol, r, q)
    pricer = AndersenLakeRepresentation(model, tte, 1e-8, 10, 32, 101, 101)
    refPrice = priceAmerican(pricer, strike, spot)
    pricer = AndersenLakeRepresentation(model, tte, 1e-8, 7, 8, 15, 31)
    price = priceAmerican(pricer, strike, spot)
    println(price, " ", refPrice, " ", price - refPrice)
    @test isapprox(refPrice, price, atol = 1e-5)
    pricer = AndersenLakeRepresentation(model, tte, 1e-8, 16, 16, 31, 63)
    price = priceAmerican(pricer, strike, spot)
    println(price, " ", refPrice, " ", price - refPrice)
    @test isapprox(refPrice, price, atol = 1e-5)

    spot = 100.0
    vol = 0.1
    tte = 0.25
    q = 0.0
    r = 0.02
    model = ConstantBlackModel(vol, r, q)
    pricer = AndersenLakeRepresentation(model, tte, 1e-8, 10, 32, 101, 101)
    refPrice = priceAmerican(pricer, strike, spot)
    pricer = AndersenLakeRepresentation(model, tte, 1e-8, 7, 8, 15, 31)
    price = priceAmerican(pricer, strike, spot)
    println(price, " ", refPrice, " ", price - refPrice)
    @test isapprox(refPrice, price, atol = 1e-5)

end

@testset "ALInstability" begin
    spot = 100.0
    strike = 100.0
    vol = 0.1
    r = 0.1
    q = 0.001
    tte = 3.0
    isCall = false
    model = ConstantBlackModel(vol, r, q)
    pricer = AndersenLakeRepresentation(model, tte, 1e-8, 10, 32, 101, 101)
    price = priceAmerican(pricer, strike, spot)
    println(price)
    @test isapprox(1.7910185909001428, price, atol = 1e-12)
end

skipRInfQ(rInfq::Int, r::Float64, q::Float64)::Bool =
    (rInfq == 1 && r >= q) || (rInfq == -1 && r <= q)

function batchAL(threshold, rInfq, isBatch, m, n, l, p)
    rs = [0.02, 0.04, 0.06, 0.08, 0.10]
    qs = [0.0, 0.04, 0.08, 0.12]
    spots = [25.0, 50.0, 80.0, 90.0, 100.0, 110.0, 120.0, 150.0, 175.0, 200.0]
    ttes = [1.0 / 12, 0.25, 0.5, 0.75, 1.0]
    sigmas = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6]
    strike = 100.0

    refPrices = zeros((length(rs), length(qs), length(ttes), length(sigmas), length(spots)))
    @elapsed for (ir, r) in enumerate(rs), (iq, q) in enumerate(qs)
        if !skipRInfQ(rInfq, r, q)
            for (it, tte) in enumerate(ttes), (iv, vol) in enumerate(sigmas)
                model = ConstantBlackModel(vol, r, q)
                pricer = AndersenLakeRepresentation(model, tte, 1e-8, 16, 16, 31, 63)
                for (is, spot) in enumerate(spots)
                    refPrices[ir, iq, it, iv, is] = priceAmerican(pricer, strike, spot)
                end

            end
        end

    end
    prices = zeros((length(rs), length(qs), length(ttes), length(sigmas), length(spots)))
    elap = @elapsed for (ir, r) in enumerate(rs), (iq, q) in enumerate(qs)
        if !skipRInfQ(rInfq, r, q)
            for (it, tte) in enumerate(ttes), (iv, vol) in enumerate(sigmas)
                model = ConstantBlackModel(vol, r, q)
                local pricer
                if isBatch
                    pricer = AndersenLakeRepresentation(model, tte, 1e-8, m, n, l, p)
                end
                for (is, spot) in enumerate(spots)
                    if refPrices[ir, iq, it, iv, is] > threshold
                        if !isBatch
                            pricer =
                                AndersenLakeRepresentation(model, tte, 1e-8, m, n, l, p)
                        end
                        prices[ir, iq, it, iv, is] = priceAmerican(pricer, strike, spot)
                        if abs(prices[ir, iq, it, iv, is] - refPrices[ir, iq, it, iv, is]) > 0.1
                            println(
                                ir,
                                " ",
                                iq,
                                " ",
                                it,
                                " ",
                                iv,
                                " ",
                                is,
                                " ",
                                refPrices[ir, iq, it, iv, is],
                                " ",
                                prices[ir, iq, it, iv, is],
                                " ",
                                spot,
                                " ",
                                vol,
                                " ",
                                tte,
                                " ",
                                q,
                                " ",
                                r,
                            )
                        end
                    end
                end

            end
        end

    end
    return prices, refPrices, elap
end
@testset "ALPositiveSet" begin
    threshold = 0.5
    rInfq = 0 #1, -1 or 0
    isBatch = false
    prices, refPrices, elap = batchAL(threshold, rInfq, isBatch, 7, 8, 15, 31)
    thrIndices = findall(z -> z > threshold, refPrices)
    mae = maxad(prices[thrIndices], refPrices[thrIndices])
    mre = maxad(prices[thrIndices] ./ refPrices[thrIndices], ones(length(thrIndices)))
    rmse = rmsd(prices[thrIndices], refPrices[thrIndices])
    println(
        "AL ",
        7,
        " ",
        rmse,
        " ",
        mae,
        " ",
        mre,
        " ",
        elap,
        " ",
        length(thrIndices) / elap,
    )

end

@testset "QDHalley" begin
    r = 0.02
    q = 0.04
    vol = 0.4
    T = 0.015
    K = 100.0
    Szero = 90.0
    #Szero = K * min(r / q, 1)
    T = 5.0
    model = ConstantBlackModel(vol, r, q)
    for i = 1:100
        tte = T * i / 100
        # americanBoundaryPutQDP(false, Szero, K, r, q, tte, vol, 1e-6, Halley)
        # americanBoundaryPutQDP(false, Szero, K, r, q, tte, vol, 1e-6, InverseQuadratic)
        Szero = americanBoundaryPutQDP(false, model, Szero, K, tte, 1e-6)
        # Szero = americanBoundaryPutQDP(false, Szero, K, r, q, tte, vol, 1e-6, SuperHalley)
        #        println("Szero = ", Szero)
        @test isless(Szero, 50.0)
    end
    @test isapprox(Szero, 22.665, atol = 1e-3)

end


text:
using POMDPs
using POMDPModels
using POMDPModelTools
using StaticArrays
using JLD2, FileIO
using Random
using Logging
using Distributions
using Statistics
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
using CMSSPs

rng = MersenneTwister(5678)

inhor_fn = ARGS[1]
outhor_fn = ARGS[2]
beta = parse(Float64, ARGS[3])

arguments = ["./paramsets/scale-1.toml","./paramsets/simtime-1.toml",
             "./paramsets/cost-1.toml", "30"]

scale_file = arguments[1]
simtime_file = arguments[2]
cost_file = arguments[3]
num_trials = parse(Int64, arguments[4])

# First, parse parameter files to get filenames and construct params object
params = parse_params(scale_file=scale_file, simtime_file=simtime_file, cost_file=cost_file)

# Get MDP and policy - DIFFERENT SEEDS
cf_mdp = get_cf_mdp(MultiRotorUAVState, MultiRotorUAVAction, params, beta)
hopon_policy = load_modal_horizon_policy_localapprox(inhor_fn, outhor_fn)

costs = Vector{Float64}(undef, 0)
successes = Vector{Int64}(undef, 0)

sim = HopOnOffSingleCarSimulator(params, rng)

for i = 1:num_trials

    tot_cost = 0.0

    # Go to 0.0, 0.0
    curr_state = generate_start_state(MultiRotorUAVState, params, rng)
    goal_state = get_state_at_rest(MultiRotorUAVState, Point(0.0, 0.0))

    reset_sim(sim)

    @show curr_state, mean(sim.time_to_finish)

    while true

        # @show curr_state, mean(sim.time_to_finish)

        tp_dist = get_arrival_time_distribution(0, mean(sim.time_to_finish), params, rng)

        a = get_best_intramodal_action(hopon_policy, 0, tp_dist, curr_state, goal_state)
        # @show a
        # readline()

        if a == nothing
            println("ABORT!")
            break
        end

        (new_state, reward) = generate_sr(cf_mdp, curr_state, a.action, rng)
        curr_state = new_state
        tot_cost += -1.0*reward

        is_done = step_sim(sim)

        if is_done
            if isterminal(cf_mdp, get_relative_state(cf_mdp, curr_state, goal_state))
                @info "Success!"
                push!(successes, 1)
            else
                @info "Failure!"
                push!(successes, 0)
            end
            break
        end
    end

    push!(costs, tot_cost)
end

@show costs
@show successes


text:
# This file contains the implementation of explicit Runkge-Kutta
# solver from (Hairer & Wanner 1992 p.134, p.165-169).

include("tableaus.jl")

# intermediate level interface

"""

A general Runge-Kutta stepper (it cen represent either, a fixed step
or an adaptive step algorithm).

"""
immutable RKStepper{Step,T} <: AbstractStepper
    # m3: this is an abstract type.  Is that ok?

    # pwl: I think this is fine, otherwise we would have to add even
    # more type parameters to RKStepper
    tableau::Tableau
    function RKStepper(tab)
        if Step == :fixed && isadaptive(tab)
            error("Cannot construct a fixed step method from an adaptive step tableau")
        elseif Step == :adaptive && !isadaptive(tab)
            error("Cannot construct an adaptive step method from an fixed step tableau")
        end
        new(convert(T,tab))
    end
end

typealias RKStepperFixed{T}    RKStepper{:fixed,   T}
typealias RKStepperAdaptive{T} RKStepper{:adaptive,T}

order(stepper::RKStepper) = minimum(order(stepper.tableau))

# TODO: possibly handle the initial stepsize and the tableau conversion here?
solve{S,T}(ode::ExplicitODE, stepper::RKStepper{S,T}, options::Options{T}) =
    Solver{RKStepper{S,T}}(ode,stepper,options)

# lower level interface

# explicit RK stepper

"""

Pre allocated arrays to store temporary data.  Used only by
Runge-Kutta stepper.

"""
type RKWorkArrays{T}
    y   ::T
    ynew::T
    yerr::T
    ks  ::Vector{T}
end

"""
State for the Runge-Kutta stepper.
"""
type RKState{T,S} <: AbstractState
    step    ::Step{T,S}
    dt      ::T
    work    ::RKWorkArrays{S}
    timeout ::Int
    # This is not currently incremented with each step
    iters   ::Int
end

function show(io::IO, state::RKState)
    show(state.step)
    println("dt      = $(state.dt)")
    println("timeout = $(state.timeout)")
    println("work    = $(state.work)")
end

function start{S,T}(s::Solver{RKStepper{S,T}})
    t0, dt0, y0 = s.ode.t0, s.options.initstep, s.ode.y0

    # TODO: we should do the Butcher table conversion somewhere
    lk = lengthks(s.stepper.tableau)
    work = RKWorkArrays(zero(y0), # y
                        zero(y0), # ynew
                        zero(y0), # yerr
                        Array(typeof(y0), lk)) # ks

    # we have to allocate each component separately
    for i = 1:lk
        work.ks[i]=zero(y0)
    end

    # pre-initialize work.ks[1]
    s.ode.F!(t0,y0,work.ks[1])

    step = Step(t0,deepcopy(y0),deepcopy(work.ks[1]))

    timeout = 0 # for step control
    return RKState(step,dt0,work,timeout,0)
end

#####################
# Fixed step method #
#####################

function next{T}(s::Solver{RKStepperFixed{T}}, state::RKState)
    step = state.step
    work = state.work

    dof  = length(step.y)
    b    = s.stepper.tableau.b
    dt   = min(state.dt,s.options.tstop-step.t)

    # m3: why is it necessary to copy here and then copy back below?

    # pwl: to my understanding calc_next_k! needs the starting value
    # (step.y) but work.ynew is changing in the inner loop, so we need
    # two distinct arrays, both starting as step.y.

    copy!(work.ynew,step.y)

    for k=1:length(b)
        calc_next_k!(work, k, s.ode, s.stepper.tableau, step, dt)
        for d=1:dof
            work.ynew[d] += dt * b[k]*work.ks[k][d]
        end
    end
    step.t += dt
    copy!(step.y,work.ynew)
    return ((step.t,step.y), state)
end

########################
# Adaptive step method #
########################

function next{T}(sol::Solver{RKStepperAdaptive{T}}, state::RKState)

    const timeout_const = 5

    # the initial values
    dt      = state.dt          # dt is the previous stepisze, it is
    # modified inside the loop
    timeout = state.timeout
    work    = state.work
    step    = state.step
    tableau = sol.stepper.tableau

    # The while loop continues until we either find a stepsize which
    # leads to a small enough error or the stepsize reaches
    # prob.minstep

    # trim the inital stepsize to avoid overshooting
    dt      = min(dt, sol.options.tstop-state.step.t)

    while true

        # Do one step (assumes ks[1]==f0).  After calling work.ynew
        # holds the new step.
        # TODO: return ynew instead of passing it as work.ynew?

        # work.y and work.yerr and work.ks are updated after this step
        rk_embedded_step!(work, sol.ode, tableau, step, dt)

        # changes work.yerr
        err, newdt, timeout = stepsize_hw92!(work, step, tableau, dt, timeout, sol.options)

        # trim again in case newdt > dt
        newdt = min(newdt, sol.options.tstop-state.step.t)

        if abs(newdt) < sol.options.minstep  # minimum step size reached, break
            # passing the newdt to state will result in done()
            state.dt = newdt
            break
        end

        if err > 1 # error is too large, repeat the step with smaller dt
            # redo step with smaller dt and reset the timeout
            dt      = newdt
            timeout = timeout_const
        else
            # step is accepted

            # preload ks[1] for the next step
            if isFSAL(sol.stepper.tableau)
                copy!(work.ks[1],work.ks[end])
            else
                sol.ode.F!(step.t+dt, work.ynew, work.ks[1])
            end

            # Swap bindings of y and ytrial, avoids one copy
            step.y, work.ynew = work.ynew, step.y

            # Update state with the data from the step we have just
            # made:
            step.t += dt
            state.dt = newdt
            state.timeout = timeout
            break
        end
    end
    return ((step.t,step.y),state)
end

##########################
# Lower level algorithms #
##########################

function rk_embedded_step!(work      ::RKWorkArrays,
                           ode       ::ExplicitODE,
                           tableau   ::Tableau,
                           last_step ::Step,
                           dt)
    # Does one embedded R-K step updating work.ynew, work.yerr and work.ks.
    # Assumes that work.ks[:,1] is already calculated!
    # Modifies work.y, work.ynew and work.yerr only

    y      = last_step.y
    dof    = length(y)
    b      = tableau.b

    fill!(work.ynew, zero(eltype(y)))
    fill!(work.yerr, zero(eltype(y)))

    for s=1:lengthks(tableau)
        # we skip the first step beacause we assume that work.ks[1] is
        # already computed
        if s > 1
            calc_next_k!(work, s, ode, tableau, last_step, dt)
        end
        for d=1:dof
            work.ynew[d] += b[1,s]*work.ks[s][d]
            work.yerr[d] += b[2,s]*work.ks[s][d]
        end
    end

    for d=1:dof
        work.yerr[d] = dt*(work.ynew[d]-work.yerr[d])
        work.ynew[d] = y[d] + dt*work.ynew[d]
    end

end

function stepsize_hw92!{T}(work,
                           last_step ::Step,
                           tableau   ::TableauRKExplicit,
                           dt        ::T,
                           timeout,
                           options   ::Options)
    # Estimates the error and a new step size following Hairer &
    # Wanner 1992, p167 (with some modifications)
    #
    # If timeout>0 no step size increase is allowed, timeout is
    # decremented in here.
    #
    # Returns the error, newdt and the number of timeout-steps
    #
    # TODO:
    # - allow component-wise reltol and abstol?
    # - allow other norms

    ord = minimum(order(tableau))
    timout_after_nan = 5
    fac = [T(8//10), T(9//10), T(1//4)^(1//(ord+1)), T(38//100)^(1//(ord+1))][1]
    facmax = T(5) # maximal step size increase. 1.5-5
    facmin = 1./facmax  # maximal step size decrease. ?
    dof = length(last_step.y)

    # in-place calculate yerr./tol
    for d=1:dof

        # if outside of domain (usually NaN) then make step size smaller by maximum
        if options.isoutofdomain(work.y[d])
            return T(10), dt*facmin, timout_after_nan
        end

        y0 = last_step.y[d] # TODO: is this supposed to be the last successful step?
        y1 = work.ynew[d]    # the approximation to the next step
        sci = (options.abstol + options.reltol*max(options.norm(y0),options.norm(y1)))
        work.yerr[d] = work.yerr[d]/sci # Eq 4.10
    end

    # TOOD: should we use options.norm here as well?
    err   = norm(work.yerr) # Eq. 4.11
    newdt = min(options.maxstep, dt*max(facmin, fac*(1/err)^(1//(ord+1)))) # Eq 4.13 modified

    if timeout > 0
        newdt = min(newdt, dt)
        timeout -= 1
    end

    return err, newdt, timeout
end

# For clarity we pass the RKWorkArrays part of the state separately,
# this is the only part of state that can be changed here
function calc_next_k!(work      ::RKWorkArrays,
                      i         ::Int,
                      ode       ::ExplicitODE,
                      tableau   ::Tableau,
                      last_step ::Step,
                      dt)
    dof = length(last_step.y)
    t, a, c = last_step.t, tableau.a, tableau.c

    copy!(work.y,last_step.y)
    for j=1:i-1
        for d=1:dof
            work.y[d] += dt * work.ks[j][d] * a[i,j]
        end
    end
    ode.F!(t + c[i]*dt, work.y, work.ks[i])
end


text:
"""
Julia package for creating new packages quickly. See [`PkgSkeleton.generate`](@ref).
"""
module PkgSkeleton

using ArgCheck: @argcheck
import Dates
using DocStringExtensions: SIGNATURES
import LibGit2
import UUIDs

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

####
#### Template values
####

"""
$(SIGNATURES)

Populate a vector of replacement pairs, either from arguments or by querying global settings
and state.
"""
function get_replacement_values(; pkg_name)
    c = LibGit2.GitConfig()     # global configuration
    #_getgitopt(opt, type = AbstractString) = LibGit2.get(type, c, opt)
    ["{UUID}" => UUIDs.uuid4(),
     "{PKGNAME}" => pkg_name,
     "{GHUSER}" => "tbeason", # _getgitopt("github.user"),
     "{USERNAME}" => "Tyler Beason", #_getgitopt("user.name"),
     "{USEREMAIL}" => "tbeas12@gmail.com", #_getgitopt("user.email"),
     "{YEAR}" => Dates.year(Dates.now())]
end

####
#### Template application
####

"""
$(SIGNATURES)

Replace multiple pairs in `str`, using `replace` iteratively.
"""
replace_multiple(str, replacements) = foldl(replace, replacements; init = str)

"""
$(SIGNATURES)

Copy from `src_dir` to `dest_dir` recursively, making the substitutions of

1. file contents and

2. filenames

using `replacements`.

Existing files are not overwritten when `skip_existing_files = true` (the default).

Directory names are *not* replaced.

Return a list of `source => dest` path pairs, with `source ≡ nothing` when `dest` was not
overwritten.
"""
function copy_and_substitute(src_dir, dest_dir, replacements;
                             skip_existing_files::Bool = true)
    results = Vector{Pair{Union{String,Nothing},String}}()
    for (root, dirs, files) in walkdir(src_dir)
        sub_dir = relpath(root, src_dir)
        mkpath(normpath(dest_dir, sub_dir))
        for file in files
            srcfile = joinpath(root, file)
            destfile = normpath(joinpath(dest_dir, sub_dir,
                                         replace_multiple(file, replacements)))
            if isfile(destfile) && skip_existing_files
                push!(results, nothing => destfile)
            else
                push!(results, srcfile => destfile)
                srcstring = read(srcfile, String)
                deststring = replace_multiple(srcstring, replacements)
                write(destfile, deststring)
            end
        end
    end
    results
end

"""
$(SIGNATURES)

Return the template directory for `name`. Symbols refer to built-in templates, while strings
are considered paths. Directories are always verified to exist.
"""
function resolve_template_dir(name::Symbol)
    dir = abspath(joinpath(@__DIR__, "..", "templates", String(name)))
    @argcheck isdir(dir) "Could not find built-in template $(name)."
    dir
end

function resolve_template_dir(dir::AbstractString)
    @argcheck isdir(dir) "Could not find directory $(dir)."
    dir
end

"""
$(SIGNATURES)

Extract a package name using the last component of a path.

The following all result in `"Foo"`:

```julia
pkg_name_from_path("/tmp/Foo")
pkg_name_from_path("/tmp/Foo/")
pkg_name_from_path("/tmp/Foo.jl")
pkg_name_from_path("/tmp/Foo.jl/")
```
"""
function pkg_name_from_path(path::AbstractString)
    base, ext = splitext(basename(isdirpath(path) ? dirname(path) : path))
    @argcheck isempty(ext) || ext == ".jl" "Invalid extension $(ext), specify package name manually."
    base
end

####
#### exposed API
####

"""
$(SIGNATURES)

Generate the skeleton for a Julia package in `dest_dir`.

The directory is transformed with `expanduser`, replacing `~` in paths.

# Arguments

`template` specifies the template to use. Symbols (eg `:default`, which is the default)
refer to *built-in* templates delivered with this package. Strings are considered paths.

`skip_existing_dir = true` (the default) aborts package generation for existing directories.

`skip_existing_files = true` (the default) prevents existing files from being overwritten.

`pkg_name` can be used to specify a package name. Note that it is derived from `dest_dir`:
the package name is `"Foo"` for all of

1. `"/tmp/Foo"`,
2. `"/tmp/Foo/"`,
3. `"/tmp/Foo.jl"`,
4. `"/tmp/Foo.jl/"`,

Use a different name only when you know what you are doing.

`git_init = true` (the default) ensures that an *empty* repository is generated in
`dest_dir`. You still have to commit files yourself.

`docs_manifest` completes the `Manifest.toml` in the `docs` subdirectory. You usually want
this.
"""
function generate(dest_dir; template = :default,
                  skip_existing_dir::Bool = true,
                  skip_existing_files::Bool = true,
                  pkg_name = pkg_name_from_path(dest_dir),
                  git_init::Bool = true, docs_manifest::Bool = true)
    dest_dir = expanduser(dest_dir)
    # preliminary checks
    @argcheck !isfile(dest_dir) "destination $(dest_dir) is a file."
    if skip_existing_dir && isdir(dest_dir)
        @warn "destination $(dest_dir) exists, skipping package generation.\nConsider `skip_existing_dir = false`."
        return false
    end

    # copy and substitute
    @info "getting template values"
    replacements = get_replacement_values(; pkg_name = pkg_name)
    @info "copy and substitute"
    results = copy_and_substitute(resolve_template_dir(template), dest_dir, replacements;
                                  skip_existing_files = skip_existing_files)
    for (src, dest) in results
        src ≡ nothing && println(stderr, "not overwriting $(dest)")
    end

    # git initialization
    if git_init
        @info "initializing git repository"
        LibGit2.init(dest_dir)
    end

    # docs manifest
    if docs_manifest
        @info "adding documenter (completing the Manifest.toml for docs)"
        docs = joinpath(dest_dir, "docs")
        cd(docs) do
            run(`$(Base.julia_cmd()) --project=$(docs) -e 'import Pkg; Pkg.add("Documenter"); Pkg.develop(Pkg.PackageSpec(; path = ".."))'`)
        end
    end

    # done
    @info "successfully generated $(pkg_name)" dest_dir
    true
end

end # module


text:
# matlab like function definitions

module odesample

# harmonic ode with damping
# Lotka–Volterra equations
# duffin eq
# hopf
# van der pol
# heat eq
# logistic
# lorenz
# euler rigid body

# =========================================================================== #
# General
# =========================================================================== #

function std_zero(x)
    return 0.0
end

# =========================================================================== #
# Harmonic Oscillator
# =========================================================================== #

struct HarmonicOscillatorOptions{T<:Function}
    m::Float64
    d::Float64
    c::Float64
    f::T # function
    function HarmonicOscillatorOptions{T}(m, d, c, f::T) where {T <: Function}
        return new(m, d, c, f)
    end
end
function HarmonicOscillatorOptions(m, d, c, f)
    HarmonicOscillatorOptions{typeof(f)}(m, d, c, f)
end
function HarmonicOscillatorOptions(m, d, c)
    HarmonicOscillatorOptions{typeof(std_zero)}(m, d, c, std_zero)
end

function HarmonicOscillatorOptions()
    return HarmonicOscillatorOptions(1.0, 0.3, 2.0)
end

function HarmonicOscillatorOptions(;
    m = 1.0,
    d = 0.3,
    c = 2.0,
    f = std_zero)
    return HarmonicOscillatorOptions(m, d, c, f)
end

function harmonic_oscillator!(der_q, t, q, options::HarmonicOscillatorOptions)

    x = q[1]
    v = q[2]

    dxdt = v
    dvdt = 1/options.m * (options.f(t) - options.c * x - options.d * v )

    der_q[1] = dxdt
    der_q[2] = dvdt

    return der_q
end

function harmonic_oscillator_controller!(der_q, t, q, options::HarmonicOscillatorOptions)

    x = q[1]
    v = q[2]

    dxdt = v
    dvdt = 1/options.m * (options.f(t, q) - options.c * x - options.d * v )

    der_q[1] = dxdt
    der_q[2] = dvdt

    return der_q
end

# =========================================================================== #
# Brusselator equations
# =========================================================================== #

struct BrusselatorOptions
    a::Float64
    b::Float64
end
function BrusselatorOptions()
    return BrusselatorOptions(1.0, 3.0)
end

function brusselator!(der_q, t, q, options::BrusselatorOptions)
    x = q[1]
    y = q[2]

    dxdt = options.a + x^2 * y - options.b * x - x
    dydt = options.b * x - x^2 * y

    der_q[1] = dxdt
    der_q[2] = dydt

    return der_q
end

# =========================================================================== #
# Lotka–Volterra equations
# =========================================================================== #

struct LotkaVolterraOptions
    a::Float64
    b::Float64
    c::Float64
    d::Float64
end
function LotkaVolterraOptions()
    return LotkaVolterraOptions(1.1, 0.4, 0.4, 0.1)
end

function lotkavolterra!(der_q, t, q, options::LotkaVolterraOptions)
    x = q[1]
    y = q[2]

    dxdt = options.a * x - options.b * x * y
    dydt = options.d * x * y - options.c * y

    der_q[1] = dxdt
    der_q[2] = dydt

    return der_q
end

# =========================================================================== #
# Duffin equations
# =========================================================================== #

struct DuffinOptions
    d::Float64
    c1::Float64
    c3::Float64
    F::Float64
    omega::Float64
end
function DuffinOptions()
    return DuffinOptions(0.02, 1.0, 5.0, 8.0, 0.5)
end

function duffin!(der_q, t, q, options)
    x = q[1]
    v = q[2]

    dxdt = v
    dvdt = omega.F * cos(options.omega * t) - options.d * v - options.c1 * x - options.c3 * x^3

    der_q[1] = dxdt
    der_q[2] = dvdt

    return der_q
end

# =========================================================================== #
# Hopf equations
# =========================================================================== #

struct HopfOptions
    a::Float64
    b::Float64
end
function HopfOptions()
    return HopfOptions(0.1, 0.75)
end

function hopf!(der_q, t, q, options::HopfOptions)
    x = q[1]
    y = q[2]

    dxdt = - x + options.a * y + x^2 * y
    dydt = options.b - options.a *y - x^2 * y

    der_q[1] = dxdt
    der_q[2] = dydt

    return der_q
end

# =========================================================================== #
# Van der Pol equations
# =========================================================================== #

struct VanDerPolOptions
    d::Float64
end
function VanDerPolOptions()
    return VanDerPolOptions(1.0)
end

function van_der_pol!(der_q, t, q, options::VanDerPolOptions)

    x = q[1]
    v = q[2]

    dxdt = v
    dvdt = options.d * (1 - x^2) * v - x

    der_q[1] = dxdt
    der_q[2] = dvdt

    return der_q
end

# =========================================================================== #
# Heat equations
# =========================================================================== #

struct HeatOptions
    tau::Float64
end
function HeatOptions()
    return HeatOptions(2.0)
end

function heat_eq!(der_q, t, q, options::HeatOptions)
    der_q[1] = -options.tau * q
    return der_q
end

# =========================================================================== #
# Logistic equations
# =========================================================================== #

struct LogisticOptions
    tau::Float64
end
function LogisticOptions()
    return LogisticOptions(1.0)
end

function logistic!(der_q, t, q, options::LogisticOptions)
    der_q[1] = options.tau * q * (1 - q)
    return der_q
end

# =========================================================================== #
# Lorenz equations
# =========================================================================== #

struct LorenzOptions
    sigma::Float64
    rho::Float64
    beta::Float64
end
function LorenzOptions()
    return LorenzOptions(10.0, 28.0, 8/3)
end

function lorenz!(der_q, t, q, options::LorenzOptions)

    x = q[1]
    y = q[2]
    z = q[3]

    dxdt = options.sigma * (y - x)
    dydt = x * (options.rho - z) - y
    dzdt = x * y - options.beta * z

    der_q[1] = dxdt
    der_q[2] = dydt
    der_q[3] = dzdt

    return der_q
end

# =========================================================================== #
# Euler Rigid Body equations
# =========================================================================== #
# not implemented yet

# =========================================================================== #
# Stiff test 1
# =========================================================================== #

function stiff_test_v1(der_q, t, q)

    der_q[1] = -0.013 * q[1] - 1000 * q[1] * q[3]
    der_q[2] = -2500 * q[2] * q[3]
    der_q[3] = -0.013 * q[1] - 1000 * q[1] * q[3] -2500 * q[2] * q[3]

end

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module ClassicalOrthogonalPolynomials
using LazyBandedMatrices: LazyBandedLayout
using InfiniteArrays: parentindices
using IntervalSets: UnitRange
using ContinuumArrays, QuasiArrays, LazyArrays, FillArrays, BandedMatrices, BlockArrays,
    IntervalSets, DomainSets, ArrayLayouts, SpecialFunctions,
    InfiniteLinearAlgebra, InfiniteArrays, LinearAlgebra, FastGaussQuadrature, FastTransforms, FFTW,
    LazyBandedMatrices, HypergeometricFunctions

import Base: @_inline_meta, axes, getindex, unsafe_getindex, convert, prod, *, /, \, +, -,
                IndexStyle, IndexLinear, ==, OneTo, tail, similar, copyto!, copy,
                first, last, Slice, size, length, axes, IdentityUnitRange, sum, _sum, cumsum,
                to_indices, _maybetail, tail, getproperty, inv, show, isapprox, summary
import Base.Broadcast: materialize, BroadcastStyle, broadcasted, Broadcasted
import LazyArrays: MemoryLayout, Applied, ApplyStyle, flatten, _flatten, colsupport, adjointlayout,
                sub_materialize, arguments, sub_paddeddata, paddeddata, PaddedLayout, resizedata!, LazyVector, ApplyLayout, call,
                _mul_arguments, CachedVector, CachedMatrix, LazyVector, LazyMatrix, axpy!, AbstractLazyLayout, BroadcastLayout,
                AbstractCachedVector, AbstractCachedMatrix, paddeddata, cache_filldata!
import ArrayLayouts: MatMulVecAdd, materialize!, _fill_lmul!, sublayout, sub_materialize, lmul!, ldiv!, ldiv, transposelayout, triangulardata,
                        subdiagonaldata, diagonaldata, supdiagonaldata
import LazyBandedMatrices: SymTridiagonal, Bidiagonal, Tridiagonal, unitblocks, BlockRange1, AbstractLazyBandedLayout
import LinearAlgebra: pinv, factorize, qr, adjoint, transpose, dot
import BandedMatrices: AbstractBandedLayout, AbstractBandedMatrix, _BandedMatrix, bandeddata
import FillArrays: AbstractFill, getindex_value, SquareEye

import QuasiArrays: cardinality, checkindex, QuasiAdjoint, QuasiTranspose, Inclusion, SubQuasiArray,
                    QuasiDiagonal, MulQuasiArray, MulQuasiMatrix, MulQuasiVector, QuasiMatMulMat,
                    ApplyQuasiArray, ApplyQuasiMatrix, LazyQuasiArrayApplyStyle, AbstractQuasiArrayApplyStyle,
                    LazyQuasiArray, LazyQuasiVector, LazyQuasiMatrix, LazyLayout, LazyQuasiArrayStyle,
                    _getindex, layout_getindex, _factorize, AbstractQuasiArray, AbstractQuasiMatrix, AbstractQuasiVector

import InfiniteArrays: OneToInf, InfAxes, Infinity, AbstractInfUnitRange, InfiniteCardinal, InfRanges
import ContinuumArrays: Basis, Weight, basis, @simplify, Identity, AbstractAffineQuasiVector, ProjectionFactorization,
    inbounds_getindex, grid, plotgrid, transform, transform_ldiv, TransformFactorization, QInfAxes, broadcastbasis, Expansion,
    AffineQuasiVector, AffineMap, WeightLayout, WeightedBasisLayout, WeightedBasisLayouts, demap, AbstractBasisLayout, BasisLayout,
    checkpoints, weight, unweightedbasis, MappedBasisLayouts, __sum, invmap
import FastTransforms: Λ, forwardrecurrence, forwardrecurrence!, _forwardrecurrence!, clenshaw, clenshaw!,
                        _forwardrecurrence_next, _clenshaw_next, check_clenshaw_recurrences, ChebyshevGrid, chebyshevpoints

import FastGaussQuadrature: jacobimoment

import BlockArrays: blockedrange, _BlockedUnitRange, unblock, _BlockArray, block, blockindex, BlockSlice
import BandedMatrices: bandwidths

export OrthogonalPolynomial, Normalized, orthonormalpolynomial, LanczosPolynomial,
            Hermite, Jacobi, Legendre, Chebyshev, ChebyshevT, ChebyshevU, ChebyshevInterval, Ultraspherical, Fourier, Laguerre,
            HermiteWeight, JacobiWeight, ChebyshevWeight, ChebyshevGrid, ChebyshevTWeight, ChebyshevUWeight, UltrasphericalWeight, LegendreWeight, LaguerreWeight,
            WeightedUltraspherical, WeightedChebyshev, WeightedChebyshevT, WeightedChebyshevU, WeightedJacobi,
            ∞, Derivative, .., Inclusion,
            chebyshevt, chebyshevu, legendre, jacobi, ultraspherical,
            legendrep, jacobip, ultrasphericalc, laguerrel,hermiteh, normalizedjacobip,
            jacobimatrix, jacobiweight, legendreweight, chebyshevtweight, chebyshevuweight, Weighted, PiecewiseInterlace

import Base: oneto

include("interlace.jl")

cardinality(::FullSpace{<:AbstractFloat}) = ℵ₁
cardinality(::EuclideanDomain) = ℵ₁

transform_ldiv(A::AbstractQuasiArray{T}, f::AbstractQuasiArray{V}, ::Tuple{<:Any,InfiniteCardinal{0}}) where {T,V}  =
    adaptivetransform_ldiv(convert(AbstractQuasiArray{promote_type(T,V)}, A), f)

function chop!(c::AbstractVector, tol::Real)
    @assert tol >= 0

    for k=length(c):-1:1
        if abs(c[k]) > tol
            resize!(c,k)
            return c
        end
    end
    resize!(c,0)
    c
end

setaxis(c, ::OneToInf) = c
setaxis(c, ax::BlockedUnitRange) = PseudoBlockVector(c, (ax,))

function adaptivetransform_ldiv(A::AbstractQuasiArray{U}, f::AbstractQuasiArray{V}) where {U,V}
    T = promote_type(U,V)

    r = checkpoints(A)
    fr = f[r]
    maxabsfr = norm(fr,Inf)

    tol = 20eps(real(T))

    for n = 2 .^ (4:∞)
        An = A[:,oneto(n)]
        cfs = An \ f
        maxabsc = maximum(abs, cfs)
        if maxabsc == 0 && maxabsfr == 0
            return zeros(T,∞)
        end

        un = A * [cfs; Zeros{T}(∞)]
        # we allow for transformed coefficients being a different size
        ##TODO: how to do scaling for unnormalized bases like Jacobi?
        if maximum(abs,@views(cfs[n-2:end])) < 10tol*maxabsc &&
                all(norm.(un[r] - fr, 1) .< tol * n * maxabsfr*1000)
            return setaxis([chop!(cfs, tol); zeros(T,∞)], axes(A,2))
        end
    end
    error("Have not converged")
end

abstract type OrthogonalPolynomial{T} <: Basis{T} end

# OPs are immutable
copy(a::OrthogonalPolynomial) = a
copy(a::SubQuasiArray{<:Any,N,<:OrthogonalPolynomial}) where N = a

"""
    jacobimatrix(P)

returns the Jacobi matrix `X` associated to a quasi-matrix of orthogonal polynomials
satisfying
```julia
x = axes(P,1)
x*P == P*X
```
Note that `X` is the transpose of the usual definition of the Jacobi matrix.
"""
jacobimatrix(P) = error("Override for $(typeof(P))")

"""
    recurrencecoefficients(P)

returns a `(A,B,C)` associated with the Orthogonal Polynomials P,
satisfying for `x in axes(P,1)`
```julia
P[x,2] == (A[1]*x + B[1])*P[x,1]
P[x,n+1] == (A[n]*x + B[n])*P[x,n] - C[n]*P[x,n-1]
```
Note that `C[1]`` is unused.

The relationship with the Jacobi matrix is:
```julia
1/A[n] == X[n+1,n]
-B[n]/A[n] == X[n,n]
C[n+1]/A[n+1] == X[n,n+1]
```
"""
function recurrencecoefficients(Q::AbstractQuasiMatrix{T}) where T
    X = jacobimatrix(Q)
    c,a,b = subdiagonaldata(X), diagonaldata(X), supdiagonaldata(X)
    inv.(c), -(a ./ c), Vcat(zero(T), b) ./ c
end

const WeightedOrthogonalPolynomial{T, A<:AbstractQuasiVector, B<:OrthogonalPolynomial} = WeightedBasis{T, A, B}

function isorthogonalityweighted(wS::WeightedOrthogonalPolynomial)
    w,S = wS.args
    w == orthogonalityweight(S)
end

"""
    singularities(f)

gives the singularity structure of an expansion, e.g.,
`JacobiWeight`.
"""
singularities(::WeightLayout, w) = w
singularities(lay::BroadcastLayout, a) = singularitiesbroadcast(call(a), map(singularities, arguments(lay, a))...)
singularities(::WeightedBasisLayouts, a) = singularities(BroadcastLayout{typeof(*)}(), a)
singularities(w) = singularities(MemoryLayout(w), w)
singularities(f::Expansion) = singularities(basis(f))
singularities(S::WeightedOrthogonalPolynomial) = singularities(S.args[1])

singularities(S::SubQuasiArray) = singularities(parent(S))[parentindices(S)[1]]

struct NoSingularities end

singularities(::Number) = NoSingularities()
singularities(r::Base.RefValue) = r[] # pass through

orthogonalityweight(P::SubQuasiArray{<:Any,2,<:Any,<:Tuple{AbstractAffineQuasiVector,Slice}}) =
    orthogonalityweight(parent(P))[parentindices(P)[1]]

_weighted(w, P) = w .* P
weighted(P::AbstractQuasiMatrix) = _weighted(orthogonalityweight(P), P)

OrthogonalPolynomial(w::Weight) =error("Override for $(typeof(w))")

@simplify *(B::Identity, C::OrthogonalPolynomial) = C*jacobimatrix(C)

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::OrthogonalPolynomial)
    x == axes(C,1) || throw(DimensionMismatch())
    C*jacobimatrix(C)
end

# function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::BroadcastQuasiVector, C::OrthogonalPolynomial)
#     axes(a,1) == axes(C,1) || throw(DimensionMismatch())
#     # re-expand in OP basis
#     broadcast(*, C * (C \ a), C)
# end

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), a::AbstractAffineQuasiVector, C::OrthogonalPolynomial)
    x = axes(C,1)
    axes(a,1) == x || throw(DimensionMismatch())
    broadcast(*, C * (C \ a), C)
end

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::WeightedOrthogonalPolynomial)
    x == axes(C,1) || throw(DimensionMismatch())
    w,P = C.args
    P2, J = (x .* P).args
    @assert P == P2
    (w .* P) * J
end

##
# Multiplication for mapped and subviews x .* view(P,...)
##

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
    T = promote_type(eltype(x), eltype(C))
    x == axes(C,1) || throw(DimensionMismatch())
    P = parent(C)
    kr,jr = parentindices(C)
    y = axes(P,1)
    Y = P \ (y .* P)
    X = kr.A \ (Y     - kr.b * Eye{T}(∞))
    P[kr, :] * view(X,:,jr)
end

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), x::Inclusion, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
    T = promote_type(eltype(x), eltype(C))
    x == axes(C,1) || throw(DimensionMismatch())
    P = parent(C)
    kr,_ = parentindices(C)
    y = axes(P,1)
    Y = P \ (y .* P)
    X = kr.A \ (Y     - kr.b * Eye{T}(∞))
    P[kr, :] * X
end

# function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::AbstractQuasiVector, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Slice}})
#     T = promote_type(eltype(f), eltype(C))
#     axes(f,1) == axes(C,1) || throw(DimensionMismatch())
#     P = parent(C)
#     kr,jr = parentindices(C)
#     (f[invmap(kr)] .* P)[kr,jr]
# end

# function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::AbstractQuasiVector, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}})
#     T = promote_type(eltype(f), eltype(C))
#     axes(f,1) == axes(C,1) || throw(DimensionMismatch())
#     P = parent(C)
#     kr,jr = parentindices(C)
#     (f[invmap(kr)] .* P)[kr,jr]
# end

broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), f::Broadcasted, C::SubQuasiArray{<:Any,2,<:Any,<:Tuple{<:AbstractAffineQuasiVector,<:Any}}) =
    broadcast(*, materialize(f), C)

function jacobimatrix(C::SubQuasiArray{T,2,<:Any,<:Tuple{AbstractAffineQuasiVector,Slice}}) where T
    P = parent(C)
    kr,jr = parentindices(C)
    Y = jacobimatrix(P)
    kr.A \ (Y - kr.b * Eye{T}(size(Y,1)))
end

function recurrencecoefficients(C::SubQuasiArray{T,2,<:Any,<:Tuple{AbstractAffineQuasiVector,Slice}}) where T
    P = parent(C)
    kr,jr = parentindices(C)
    A,B,C = recurrencecoefficients(P)
    A * kr.A, A*kr.b + B, C
end

_vec(a) = vec(a)
_vec(a::InfiniteArrays.ReshapedArray) = _vec(parent(a))
_vec(a::Adjoint{<:Any,<:AbstractVector}) = a'

include("clenshaw.jl")
include("ratios.jl")
include("normalized.jl")
include("lanczos.jl")

function _tritrunc(_, X, n)
    c,a,b = subdiagonaldata(X),diagonaldata(X),supdiagonaldata(X)
    Tridiagonal(c[OneTo(n-1)],a[OneTo(n)],b[OneTo(n-1)])
end

function _tritrunc(::SymTridiagonalLayout, X, n)
    a,b = diagonaldata(X),supdiagonaldata(X)
    SymTridiagonal(a[OneTo(n)],b[OneTo(n-1)])
end

_tritrunc(X, n) = _tritrunc(MemoryLayout(X), X, n)

jacobimatrix(V::SubQuasiArray{<:Any,2,<:Any,<:Tuple{Inclusion,OneTo}}) =
    _tritrunc(jacobimatrix(parent(V)), maximum(parentindices(V)[2]))

grid(P::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,Any}}) =
    eigvals(symtridiagonalize(jacobimatrix(P)))

function golubwelsch(X)
    D, V = eigen(symtridiagonalize(X))  # Eigenvalue decomposition
    D, V[1,:].^2
end

function golubwelsch(V::SubQuasiArray)
    x,w = golubwelsch(jacobimatrix(V))
    w .*= sum(orthogonalityweight(parent(V)))
    x,w
end

function factorize(L::SubQuasiArray{T,2,<:Normalized,<:Tuple{Inclusion,OneTo}}) where T
    x,w = golubwelsch(L)
    TransformFactorization(x, L[x,:]'*Diagonal(w))
end

function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{Inclusion,OneTo}}) where T
    Q = Normalized(parent(L))[parentindices(L)...]
    D = L \ Q
    F = factorize(Q)
    TransformFactorization(F.grid, D*F.plan)
end

function factorize(L::SubQuasiArray{T,2,<:OrthogonalPolynomial,<:Tuple{<:Inclusion,<:AbstractUnitRange}}) where T
    _,jr = parentindices(L)
    ProjectionFactorization(factorize(parent(L)[:,oneto(maximum(jr))]), jr)
end

function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial}, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial})
    axes(A,1) == axes(B,1) || throw(DimensionMismatch())
    _,jA = parentindices(A)
    _,jB = parentindices(B)
    (parent(A) \ parent(B))[jA, jB]
end

function \(A::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Any,Slice}}, B::SubQuasiArray{<:Any,2,<:OrthogonalPolynomial,<:Tuple{Any,Slice}})
    axes(A,1) == axes(B,1) || throw(DimensionMismatch())
    parent(A) \ parent(B)
end

# assume we can expand w_B in wA to reduce to polynomial multiplication
function \(wA::WeightedOrthogonalPolynomial, wB::WeightedOrthogonalPolynomial)
    _,A = arguments(wA)
    w_B,B = arguments(wB)
    A \ ((A * (wA \ w_B)) .* B)
end

include("classical/hermite.jl")
include("classical/jacobi.jl")
include("classical/chebyshev.jl")
include("classical/ultraspherical.jl")
include("classical/laguerre.jl")
include("classical/fourier.jl")
include("stieltjes.jl")
include("decompOPs.jl")

export cholesky_jacobimatrix

end # module


text:
##--- Identity ---##
type IdentityTransformation <: Transformation end
transform!(::IdentityTransformation, x) = x
invert!(::IdentityTransformation, x) = x
transform(::IdentityTransformation, x) = copy(x)
invert(::IdentityTransformation, x) = copy(x)
isinvertible(::IdentityTransformation) = true
Base.inv(::IdentityTransformation) = IdentityTransformation()

##--- Invertible Transforms ---##
for (t1,f1,t2,f2) in ( (:ExpTransformation, :exp, :LogTransformation, :log),
                       (:LogisticTransformation, :logistic, :LogitTransformation, :logit) )
    @eval begin
        type $(t1) <: Transformation end
        type $(t2) <: Transformation end
        transform!(::$(t1), x) = map!($(f1),x)
        transform!(::$(t2), x) = map!($(f2),x)
        transform(::$(t1), x) = $(f1)(x)
        transform(::$(t2), x) = $(f2)(x)
        invert!(::$(t1), x) = map!($(f2),x)
        invert!(::$(t2), x) = map!($(f1),x)
        invert(::$(t1), x) = $(f2)(x)
        invert(::$(t2), x) = $(f1)(x)
        Base.inv(::$(t1)) = $(t2)()
        Base.inv(::$(t2)) = $(t1)()
        isinvertible(::$(t1)) = true
        isinvertible(::$(t2)) = true
    end
end

## --- Scaling --- ##
immutable ScaleTransformation{T<:Number} <: Transformation
    val::T
end
ScaleTransformation() = ScaleTransformation(1.0)
isinvertible(::ScaleTransformation) = true
function invert!{T}(tfm::ScaleTransformation{T}, x::AbstractArray{T})
    scale!(one(T)/tfm.val, x)
end
function transform!{T}(tfm::ScaleTransformation, x::AbstractArray{T})
    scale!(tfm.val, x)
end
Base.inv{T<:Number}(tfm::ScaleTransformation{T}) = ScaleTransformation(one(T)/tfm.val)

## --- Shifting --- ##
immutable ShiftTransformation{T<:Number} <: Transformation
    val::T
end
ShiftTransformation() = ShiftTransformation(0.0)
isinvertible(::ShiftTransformation) = true
function transform!{T}(tfm::ShiftTransformation{T}, x::AbstractArray{T})
    for i in eachindex(x); x[i] += tfm.val; end
    x
end
function invert!{T}(tfm::ShiftTransformation{T}, x::AbstractArray{T})
    for i in eachindex(x); x[i] -= tfm.val; end
    x
end
Base.inv(tfm::ShiftTransformation) = ShiftTransformation(-tfm.val)


text:
using PyCall
using .PyPlot
import .PyPlot: loglog, plot, semilogx, semilogy

### plotting Cℓs

for plot in (:plot, :loglog, :semilogx, :semilogy)
    @eval function ($plot)(ic::InterpolatedCℓs, args...; kwargs...)
		($plot)(ic.ℓ, ic.Cℓ, args...; kwargs...)
	end
	@eval function ($plot)(ic::InterpolatedCℓs{<:AbstractExtrapolation{<:Measurement}}, args...; kwargs...)
		errorbar(ic.ℓ, Measurements.value.(ic.Cℓ), Measurements.uncertainty.(ic.Cℓ), args...; marker=".", ls="", capsize=2, kwargs...)
		($plot) in [:loglog,:semilogx] && xscale("log")
		($plot) in [:loglog,:semilogy] && yscale("log")
	end
end

### plotting FlatFields

plotsize₀ = 4

pretty_name(s::Symbol) = pretty_name(Val.(Symbol.(split(string(s),"")))...)
pretty_name(::Val{s},::Val{:x}) where {s} = "$s Map"
pretty_name(::Val{s},::Val{:l}) where {s} = "$s Fourier"
pretty_name(::Val{:T},::Val{:x}) where {s} = "Map"
pretty_name(::Val{:T},::Val{:l}) where {s} = "Fourier"

# generic plotting some components of a FlatField
function _plot(f::FlatField{P}, ax, k, title, vlim; units=:deg, ticklabels=true, axeslabels=false, kwargs...) where {N,θ,P<:Flat{N,θ}}
    if string(k)[2] == 'x'
        x = θ*N/Dict(:deg=>60,:arcmin=>1)[units]/2
    elseif string(k)[2] == 'l'
        x = fieldinfo(f).nyq
    else
        throw(ArgumentError("Invalid `which`: $k"))
    end
    extent = [-x,x,-x,x]
    (title == nothing) && (title="$(pretty_name(k)) ($(N)x$(N) @ $(θ)')")
    (vlim == nothing) && (vlim=:sym)
    _plot(Array(f[k]); ax=ax, extent=extent, title=title, vlim=vlim, kwargs...)
    if ticklabels
        if string(k)[2] == 'x'
            @pydef mutable struct MyFmt <: pyimport(:matplotlib).ticker.ScalarFormatter
                __call__(self,v,p=nothing) = py"super"(MyFmt,self).__call__(v,p)*Dict(:deg=>"°",:arcmin=>"′")[units]
            end
            ax.xaxis.set_major_formatter(MyFmt())
            ax.yaxis.set_major_formatter(MyFmt())
            if axeslabels
                ax.set_xlabel("RA")
                ax.set_ylabel("Dec")
            end
        else
            ax.set_xlabel(raw"$\ell_x$")
            ax.set_ylabel(raw"$\ell_y$")
            ax.tick_params(axis="x", rotation=45)
        end
        ax.tick_params(labeltop=false, labelbottom=true)
    else
        ax.tick_params(labeltop=false, labelleft=false)
    end
end

# plotting a map
function _plot(m::AbstractMatrix{<:Real}; ax=gca(), title=nothing, vlim=:sym, cmap="RdBu_r", vscale=:linear, cbar=true, kwargs...)

    # some logic to automatically get upper/lower limits
    if vlim==:sym
        vmax = quantile(abs.(m[@. !isnan(m)][:]),0.999)
        vmin = -vmax
    elseif vlim==:asym
        vmin, vmax = (quantile(m[@. !isnan(m)][:],q) for q=(0.001,0.999))
    elseif isa(vlim,Tuple)
        vmin, vmax = vlim
    else
        vmax = vlim
        vmin = -vmax
    end

    m[isinf.(m)] .= NaN

    cax = ax.matshow(clamp.(m,vmin,vmax); vmin=vmin, vmax=vmax, cmap=cmap, rasterized=true, kwargs...)
    cbar && gcf().colorbar(cax,ax=ax)
    title!=nothing && ax.set_title(title, y=1)
    ax
end

# plotting fourier coefficients
function _plot(m::AbstractMatrix{<:Complex}; vscale=:log, kwargs...)
    dat = ifftshift(unfold(m))
    if vscale==:log
        dat .= log10.(abs.(dat))
    end
    _plot(real.(dat); vlim=(nothing,nothing), cmap=nothing, kwargs...)
end

@doc doc"""
    plot(f::Field; kwargs...)
    plot(fs::VecOrMat{\<:Field}; kwarg...)

Plotting fields.
"""
plot(f::Field; kwargs...) = plot([f]; kwargs...)
function plot(fs::AbstractVecOrMat{F}; plotsize=plotsize₀, which=default_which(fs), title=nothing, vlim=nothing, return_all=false, kwargs...) where {F<:Field}
    (m,n) = size(tuple.(fs, which)[:,:])
    fig,axs = subplots(m, n; figsize=plotsize.*[1.4*n,m], squeeze=false)
    axs = getindex.(Ref(axs), 1:m, (1:n)') # see https://github.com/JuliaPy/PyCall.jl/pull/487#issuecomment-456998345
    _plot.(fs,axs,which,title,vlim; kwargs...)
    tight_layout(w_pad=-10)
    return_all ? (fig,axs,which) : isjuno ? fig : nothing
end
default_which(::AbstractVecOrMat{<:FlatS0}) = [:Ix]
default_which(::AbstractVecOrMat{<:FlatS2}) = [:Ex :Bx]
default_which(::AbstractVecOrMat{<:FlatS02}) = [:Ix :Ex :Bx]
function default_which(fs::AbstractVecOrMat{<:Field})
    try
        ensuresame((default_which([f]) for f in fs)...)
    catch x
        x isa AssertionError ? throw(ArgumentError("Must specify `which` argument by hand for this combination of fields to plot.")) : rethrow()
    end
end

@doc doc"""
    animate(fields::Vector{\<:Vector{\<:Field}}; interval=50, motionblur=false, kwargs...)

"""
animate(f::AbstractVecOrMat{<:Field}; kwargs...) = animate([f]; kwargs...)
animate(annonate::Function, args...; kwargs...) = animate(args...; annonate=annonate, kwargs...)
function animate(fields::AbstractVecOrMat{<:AbstractVecOrMat{<:Field}}; fps=25, motionblur=false, annonate=nothing, filename=nothing, kwargs...)
    fig, axs, which = plot(first.(fields); return_all=true, kwargs...)
    motionblur = (motionblur == true) ? [0.1, 0.5, 1, 0.5, 0.1] : (motionblur == false) ? [1] : motionblur

    if (annonate!=nothing); annonate(fig,axs,which); end

    ani = pyimport("matplotlib.animation").FuncAnimation(fig,
        i->begin
            for (f,ax,k) in tuple.(fields,axs,which)
                if length(f)>1
                    img = ax.images[1]
                    img.set_data(sum(x*getindex(f[mod1(i-j+1,length(f))],k) for (j,x) in enumerate(motionblur)) / sum(motionblur))
                end
            end
            first.(getproperty.(axs,:images))[:]
        end,
        1:maximum(length.(fields)[:]),
        interval=1000/fps, blit=true
        )
    close()
    if filename!=nothing
        ani.save(filename,writer="imagemagick",savefig_kwargs=Dict(:facecolor=>fig.get_facecolor()))
        if endswith(filename,".gif")
            run(`convert -layers Optimize $filename $filename`)
        end
    end
    HTML(ani.to_html5_video())
end


text:
module TestWDISearchWDI

using Base.Test
using WorldBankDataTd
using Dates

country_data = { { "total"=>4,"per_page"=>"25000","pages"=>1,"page"=>1 },
                 {
                    ["latitude"=>"42.5075","id"=>"AND","iso2Code"=>"AD","incomeLevel"=>["id"=>"NOC","value"=>"High income: nonOECD"],"adminregion"=>["id"=>"","value"=>""],"lendingType"=>["id"=>"LNX","value"=>"Not classified"],"region"=>["id"=>"ECS","value"=>"Europe & Central Asia (all income levels)"],"capitalCity"=>"Andorra la Vella","name"=>"Andorra","longitude"=>"1.5218"],
                    ["latitude"=>"","id"=>"ARB","iso2Code"=>"1A","incomeLevel"=>["id"=>"NA","value"=>"Aggregates"],"adminregion"=>["id"=>"","value"=>""],"lendingType"=>["id"=>"","value"=>"Aggregates"],"region"=>["id"=>"NA","value"=>"Aggregates"],"capitalCity"=>"","name"=>"Arab World","longitude"=>""],
                    ["latitude"=>"24.4764","id"=>"ARE","iso2Code"=>"AE","incomeLevel"=>["id"=>"NOC","value"=>"High income: nonOECD"],"adminregion"=>["id"=>"","value"=>""],"lendingType"=>["id"=>"LNX","value"=>"Not classified"],"region"=>["id"=>"MEA","value"=>"Middle East & North Africa (all income levels)"],"capitalCity"=>"Abu Dhabi","name"=>"United Arab Emirates","longitude"=>"54.3705"],
                    ["latitude"=>"-34.6118","id"=>"ARG","iso2Code"=>"AR","incomeLevel"=>["id"=>"UMC","value"=>"Upper middle income"],"adminregion"=>["id"=>"LAC","value"=>"Latin America & Caribbean (developing only)"],"lendingType"=>["id"=>"IBD","value"=>"IBRD"],"region"=>["id"=>"LCN","value"=>"Latin America & Caribbean (all income levels)"],"capitalCity"=>"Buenos Aires","name"=>"Argentina","longitude"=>"-58.4173"]
                 }
               }
df_country = WorldBankDataTd.parse_country(country_data)

indicator_data = { { "total"=>5,"per_page"=>"25000","pages"=>1,"page"=>1 },
                   { ["sourceOrganization"=>"World Bank and International Energy Agency (IEA Statistics © OECD/IEA, http://www.iea.org/stats/index.asp).  ","id"=>"12.1_TD.LOSSES","topics"=>{},"sourceNote"=>"Transmission and distribution losses (%): Transmission and distribution (T&D) losses measure power lost in the transmission of (high-voltage) electricity from power generators to distributors and in the distribution of (medium- and low-voltage) electricity from distributors to end-users. T&D losses are represented as a percentage of gross electricity production. They include both technical and nontechnical (or commercial) losses. Included in the latter are unmetered, unbilled, and unpaid electricity, including theft, which could be significant in developing countries. Aggregate T&D system indicators may be dominated by factors other than losses. The location of primary energy resources (such as hydro lakes and coal seams) and large loads (cities and industries) may be more significant factors in T&D efficiency indicators than the losses or efficiency of the transmission system itself. Properly separating true losses (and hence the efficiency potential of transmission systems) from exogenous location and scale factors and nontechnical losses would require detailed studies of system-dynamic interactions and real operating requirements that are not practical for global tracking purposes.","name"=>"Transmission and distribution losses (%)","source"=>["id"=>"35","value"=>"Sustainable Energy for All"]],
                     ["sourceOrganization"=>"World Bank and International Energy Agency (IEA Statistics © OECD/IEA, http://www.iea.org/stats/index.asp).  ","id"=>"13.1_INDUSTRY.ENERGY.INTENSITY","topics"=>{},"sourceNote"=>"Energy intensity of industrial sector (MJ/\$2005):  A ratio between energy consumption in industry (including energy industry own use) and industry sector value added measured at purchasing power parity. Energy intensity is an indication of how much energy is used to produce one unit of economic output. Lower ratio indicates that less energy is used to produce one unit of output. ","name"=>"Energy intensity of industrial sector (MJ/\$2005)","source"=>["id"=>"35","value"=>"Sustainable Energy for All"]],
                     ["sourceOrganization"=>"World Bank and International Energy Agency (IEA Statistics © OECD/IEA, http://www.iea.org/stats/index.asp).  ","id"=>"14.1_AGR.ENERGY.INTENSITY","topics"=>{},"sourceNote"=>"Energy intensity of agricultural sector (MJ/\$2005):  A ratio between energy consumption in agricultural sector (including forestry and fishing) and agricultural sector value added measured at purchasing power parity. Energy intensity is an indication of how much energy is used to produce one unit of economic output. Lower ratio indicates that less energy is used to produce one unit of output. ","name"=>"Energy intensity of agricultural sector (MJ/\$2005)","source"=>["id"=>"35","value"=>"Sustainable Energy for All"]],
                     ["sourceOrganization"=>"World Bank and International Energy Agency (IEA Statistics © OECD/IEA, http://www.iea.org/stats/index.asp).  ","id"=>"15.1_OTHER.SECT.ENER.INTENS","topics"=>{},"sourceNote"=>"Energy intensity of other sectors (MJ/\$2005): A ratio between energy consumption in other sectors (including services, residential and transport) and services sector value added measured at purchasing power parity. Energy intensity is an indication of how much energy is used to produce one unit of economic output. Lower ratio indicates that less energy is used to produce one unit of output. ","name"=>"Energy intensity of other sectors (MJ/\$2005)","source"=>["id"=>"35","value"=>"Sustainable Energy for All"]],
                     ["sourceOrganization"=>"World Bank and International Energy Agency (IEA Statistics © OECD/IEA, http://www.iea.org/stats/index.asp).  ","id"=>"16.1_DECOMP.EFFICIENCY.IND","topics"=>{},"sourceNote"=>"The Index is created from the series of year to year Energy Intensity component of the Decomposition Analysis with the DIVISIA LMDI I method. The Significant Period is defined as the period for which the available data is consistent in the energy and value added databases such that all the sectors have energy and value added for those years. The CAGR is calculated over that Significant Period.","name"=>"Divisia Decomposition Analysis - Energy Intensity component Index","source"=>["id"=>"35","value"=>"Sustainable Energy for All"]]
                   }
                 }

df_indicator = WorldBankDataTd.parse_indicator(indicator_data)

WorldBankDataTd.set_country_cache(df_country)

@test search_wdi("countries",:name,r"(Andorr|Arg)"i)[:capital] ==  UTF8String["Andorra la Vella", "Buenos Aires"]
@test search_wdi("countries",:capital,r"(Andorr|Abu)"i)[:iso2c] ==  ASCIIString["AD", "AE"]

WorldBankDataTd.set_indicator_cache(df_indicator)

@test search_wdi("indicators",:description,r"energy intensity of"i)[:name] == UTF8String["Energy intensity of industrial sector (MJ/\$2005)", "Energy intensity of agricultural sector (MJ/\$2005)", "Energy intensity of other sectors (MJ/\$2005)"]
@test search_wdi("indicators",:description,r"transmiss"i)[:indicator] == ASCIIString["12.1_TD.LOSSES"]

end