text:
module Constructions

export Construction
export dependency_order
export @place
export @modify
export @construct
export @replace
export @remove

abstract type AbstractGeometricElement end

struct PlacedElement <: AbstractGeometricElement
    name::String
    representation::Any
    required_by::Set{String}
end

struct ConstructedElement <: AbstractGeometricElement
    name::String
    representation::Any
    required_by::Set{String}
    requires::Set{String}
    construct::Function
end

name(g::AbstractGeometricElement) = g.name
representation(g::AbstractGeometricElement) = g.representation
required_by(g::AbstractGeometricElement) = g.required_by

requires(g::ConstructedElement) = g.requires
requires(::PlacedElement) = Set{String}()

construct(g::ConstructedElement) = g.construct
construct(g::PlacedElement) = C -> representation(g)

replace_representation(g::PlacedElement, new_representation) = PlacedElement(name(g), new_representation, required_by(g))

replace_representation(g::ConstructedElement, new_representation) = ConstructedElement(name(g), new_representation, required_by(g), requires(g), construct(g))

replace_required_by(g::PlacedElement, new_required_by) = PlacedElement(name(g), representation(g), new_required_by)

replace_required_by(g::ConstructedElement, new_required_by) = ConstructedElement(name(g), representation(g), new_required_by, requires(g), construct(g))

struct Construction
    elements::Dict{String,AbstractGeometricElement}
end

Construction() = Construction(Dict{String,AbstractGeometricElement}())

function dependency_order(C::Construction)
    result = String[]
    elements = Set(keys(C.elements))
    while !isempty(elements)
        free_elements = filter(e -> isdisjoint(requires(C.elements[e]), elements), elements)
        append!(result, free_elements)
        setdiff!(elements, free_elements)
    end
    result
end

import Base.show

function show(io::IO, C::Construction)
    for en in dependency_order(C)
        e = C.elements[en]
        if is_placed(e)
            println(io, name(e), ": ", representation(e), "; ")
        elseif is_constructed(e)
            r = requires(e)
            print(io, "{")
            for s in r
                print(io, s, ", ")
            end
            println(io, "} => ", name(e), ": ", representation(e), "; ")
        end
    end
end

import Base.getindex

function getindex(C::Construction, ename::String)
    if !haskey(C.elements, ename)
        throw(ArgumentError("Element $ename not found."))
    else
        representation(C.elements[ename])
    end
end

function is_constructed(element)
    element isa ConstructedElement
end

function is_placed(element)
    element isa PlacedElement
end

function representation(ename::String)
    if !haskey(C.elements, ename)
        throw(ArgumentError("Element $ename not found."))
    else
        representation(C.elements[ename])
    end
end

function place_element!(C::Construction, ename::String, representation )
    if haskey(C.elements, ename)
        throw(ArgumentError("Name $ename is already used in the construction."))
    else
        C.elements[ename] = PlacedElement(ename, representation, Set{String}())
    end
    nothing
end

function construct_element!(C::Construction, ename::String, construct::Function, depends_on::Set{String})
    if haskey(C.elements, ename)
        throw(ArgumentError("Name $ename is already used in the construction."))
    else
        unknown_elements = filter(e -> !haskey(C.elements,e), depends_on)
        if !isempty(unknown_elements)
            throw(ArgumentError("Unknown dependencies: $unknown_elements"))
        end
        try
            representation = construct(C)
            C.elements[ename] = ConstructedElement(ename, representation, Set{String}(), depends_on, construct)
        catch err
            throw(ArgumentError("Construction rule failed with error $err."))
        end
        foreach( e -> push!(C.elements[e].required_by, ename) , depends_on)
    end
    nothing
end

function collect_dependencies(C::Construction, ename::String, result::Set{String})
    if !haskey(C.elements, ename)
        throw(ArgumentError("Element $ename not found."))
    else
        followingelements = required_by(C.elements[ename])
        union!(result, followingelements)
        for e in followingelements
            push!(result, e)
            collect_dependencies(C, e, result)
        end
    end
    result
end

function update_dependencies!(C::Construction, ename::String, requires_elements::Set{String}, required_by_elements::Set{String})
    if haskey(C.elements, ename)
        # element has been updated but the dependencies have not changed
        affectedelements = collect_dependencies(C, ename, Set{String}() )
        while !isempty(affectedelements)
            for ae in affectedelements
                element = C.elements[ae]
                if isdisjoint(requires(element), affectedelements)
                    if element isa ConstructedElement
                        try
                            C.elements[ae] = replace_representation(element, element.construct(C))
                        catch err
                            throw(ErrorException("Construction rule failed with error $err."))
                        end
                    end
                    delete!(affectedelements, ae)
                end
            end
        end
    else
        # element has been deleted
        for e in requires_elements
            if haskey(C.elements, e)
                delete!(C.elements[e].required_by, ename)
            end
        end
        for e in required_by_elements
            if haskey(C.elements, e)
                be = C.elements[e]
                delete!(C.elements, e)
                update_dependencies!(C, name(be), requires(be), required_by(be))
            end
        end
    end
    nothing
end

function modify_placed_element!(C::Construction, ename::String, newrepresentation)
    if !haskey(C.elements, ename)
        throw(ArgumentError("Element $ename not found."))
    else
        element = C.elements[ename]
        if element isa ConstructedElement
            throw(ArgumentError("Element $ename is a constructed element."))
        else
            C.elements[ename] = replace_representation(element, newrepresentation)
            update_dependencies!(C, ename, requires(element), required_by(element))
        end
    end
    nothing
end

function remove_element!(C::Construction, ename::String)
    if !haskey(C.elements, ename)
        throw(ArgumentError("Element $ename not found."))
    else
        element = C.elements[ename]
        delete!(C.elements, ename)
        update_dependencies!(C, ename, requires(element), required_by(element))
    end
    nothing
end

function replace_element!(C::Construction, ename::String, representation)
    if !haskey(C.elements, ename)
        throw(ArgumentError("Element $ename not found."))
    else
        element = C.elements[ename]
        re = requires(element)
        for e in re
            delete!(C.elements[e].required_by, ename)
        end
        delete!(C.elements, ename)
        place_element!(C, ename, representation)
        element = replace_required_by(C.elements[ename], required_by(element))
        C.elements[ename] = element
        update_dependencies!(C, ename, requires(element), required_by(element))
    end
    nothing
end

function replace_element!(C::Construction, ename::String, construct::Function, depends_on::Set{String})
    if !haskey(C.elements, ename)
        throw(ArgumentError("Element $ename not found."))
    else
        element = C.elements[ename]
        re = requires(element)
        for e in re
            delete!(C.elements[e].required_by, ename)
        end
        delete!(C.elements, ename)
        construct_element!(C, ename, construct, depends_on)
        element = replace_required_by(C.elements[ename], required_by(element))
        C.elements[ename] = element
        update_dependencies!(C, ename, requires(element), required_by(element))
    end
    nothing
end

macro place(construction, name, representation)
    quote
        place_element!($(esc(construction)), $(esc(name)), $(esc(representation)))
        $(esc(construction))[$(esc(name))]
    end
end

macro modify(construction, name, new_representation)
    quote
        modify_placed_element!($(esc(construction)), $(esc(name)), $(esc(new_representation)))
        $(esc(construction))[$(esc(name))]
    end
end

function make_construct(rule, dependencies...)
    arguments = map( d ->:(C[$(esc(d))]) , dependencies )
    quote
        C -> $(esc(rule))($(arguments...))
    end
end

macro construct(construction, name, rule, dependencies...)
    construct = make_construct(rule, dependencies...)
    quote
        construct_element!($(esc(construction)), $(esc(name)), $construct, Set{String}([$(map(e->esc(e), dependencies)...)]))
        $(esc(construction))[$(esc(name))]
    end
end

macro remove(construction, name)
    quote
        remove_element!($(esc(construction)), $(esc(name)))
    end
end

macro replace(construction, name, representation)
    quote
        replace_element!($(esc(construction)), $(esc(name)), $(esc(representation)))
        $(esc(construction))[$(esc(name))]
    end
end

macro replace(construction, name, rule, dependencies...)
    construct = make_construct(rule, dependencies...)
    quote
        replace_element!($(esc(construction)), $(esc(name)), $construct, Set{String}([$(map(e->esc(e), dependencies)...)]))
        $(esc(construction))[$(esc(name))]
    end
end

using Plots

@recipe function construction_plot_recipe(C::Construction)
    plotorder = dependency_order(C)
    for en in plotorder
        label := en
        @series C[en]
    end
end

end


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

JLLWrappers.@generate_main_file_header("MUMPS_seq")
JLLWrappers.@generate_main_file("MUMPS_seq", UUID("d7ed1dd3-d0ae-5e8e-bfb4-87a502085b8d"))
end  # module MUMPS_seq_jll

ext: jl


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

import Base.show

##### Observatory
"""
Type holding information about an observing site.  Its fields are:

* `name`: the name of the site
* `latitude`: North-ward latitude of the site in degrees
* `longitude`: East-ward longitude of the site in degrees
* `altitude`: altitude of the site in meters
* `tz`: the number of hours of offset from UTC
"""
immutable Observatory
    name::AbstractString
    latitude::Real
    longitude::Real
    altitude::Real
    tz::Real # There are non-integer time zones
    # Define constructor that automatically converts longitude and latitude with
    # "ten", for convenience.
    Observatory(name, lat, long, alt, tz) =
        new(name, ten(lat), ten(long), float(alt), ten(tz))
end

# New type representation
function show(io::IO, obs::Observatory)
    println(io, "Observatory: ", obs.name)
    println(io, "latitude:    ", obs.latitude, "°N")
    println(io, "longitude:   ", obs.longitude, "°E")
    println(io, "altitude:    ", obs.altitude, " m")
    tzdec, tzint = modf(obs.tz)
    print(io,   "time zone:   ", @sprintf("UTC%+d", tzint),
          tzdec == 0 ? "" : @sprintf(":%d", abs(tzdec*60)))
end

##### Planet
"""
Type holding information about a planet.  Its fields are:

Designation:

* `name`: the name

Physical characteristics:

* `radius`: mean radius in meters
* `eqradius`: equatorial radius in meters
* `polradius`: polar radius in meters
* `mass`: mass in kilogram

Orbital characteristics (epoch J2000):

* `ecc`: eccentricity of the orbit
* `axis`: semi-major axis of the orbit in meters
* `period`: sidereal orbital period in seconds
"""
immutable Planet
    name::AbstractString
    radius::Real
    eqradius::Real
    polradius::Real
    mass::Real
    ecc::Real
    axis::Real
    period::Real
    Planet(name, radius, eqradius, polradius, mass, ecc, axis, period) =
        new(name, float(radius), float(eqradius),
            float(polradius), float(mass), float(ecc),
            float(axis), float(period))
end

# New type representation
function show(io::IO, pl::Planet)
    println(io, "Planet:            ", ucfirst(pl.name))
    println(io, "mean radius:       ", pl.radius, " m")
    println(io, "equatorial radius: ", pl.eqradius, " m")
    println(io, "polar radius:      ", pl.polradius, " m")
    println(io, "mass:              ", pl.mass, " kg")
    println(io, "eccentricity:      ", pl.ecc)
    println(io, "semi-major axis:   ", pl.axis, " m")
    print(io,   "period:            ", pl.axis, " s")
end

export Observatory, Planet


text:
# Reference: Cormen, p. 18
# Notes: I ran into the problem of Julia's parametric typing system, because I
# wanted the algorithm to work on any array of real numbers. This meant I wanted
# concrete subtypes of Int and AbstractFloat to work correctly. Naively, I wrote
# the function signature:
#    insertionsort(A::Array{Real, 1})
# but this fails because, for example,
#    Array{Int64, 1} <: Array{Real, 1} => false
# I think that this is due to compiler optimization, if I remember the lecture I
# watched correctly. Regardless, for parametric type abstraction for methods,
# the correct syntax is given below:
function insertionsort{T<:Real}(A::Array{T,1})
    for j = 2:length(A)
        key = A[j]
        # Insert A[j] into the sorted sequence A[1..j-1].
        i = j - 1
        while i > 0 && A[i] > key
            A[i+1] = A[i]
            i -= 1
        end
        A[i+1] = key
    end
end

# Reference: Cormen, p. 34
function mergesort{T<:Real}(A::Array{T,1})
    # Reference: Cormen, p. 31
    # This was more straightforward to implement, but I developed a bit of a
    # "hacked" solution for introducing the sentinel '∞' that is used in the
    # textbook. I essentially set ∞ to the maximum allowable value for the
    # specified type, so it acts like infinity for every practical purpose.
    function merge{T<:Real}(A::Array{T,1}, p::Int, q::Int, r::Int)
        ∞ = typemax(T) - 1
        n₁, n₂ = q - p + 1, r - q
        L, R = zeros(T,n₁+1), zeros(T, n₂+1)
        L[1:end-1] = A[p:p+n₁-1]
        R[1:end-1] = A[q+1:q+n₂]
        L[n₁+1], R[n₂+1] = ∞, ∞
        i, j = 1, 1
        for k = p:r
            if L[i] <= R[j]
                A[k] = L[i]
                i += 1
            else
                A[k] = R[j]
                j += 1
            end
        end
    end

    function kernel{T<:Real}(A::Array{T,1}, p::Int, r::Int)
        if p < r
            q = Int(floor((p + r) / 2))
            kernel(A, p, q)
            kernel(A, q + 1, r)
            merge(A, p, q, r)
        end
    end

    kernel(A, 1, length(A))
end

# Reference: Cormen, p. 161
function heapsort{T<:Real}(A::Array{T,1})
    # Reference: Cormen, p. 152
    parent(i::Int) = Int(floor(i/2))
    leftchild(i::Int) = 2i
    rightchild(i::Int) = 2i + 1

    # Reference: Cormen, p. 154
    # The modification from the book is that maxheapify also takes an additional
    # parameter, 'size', which gives the size of the heap.
    function maxheapify{T<:Real}(A::Array{T,1}, i::Int, size::Int)
        # get the children of i
        l, r = leftchild(i), rightchild(i)

        # check if the left child exists and is greater than i
        if l <= size && A[l] > A[i]
            largest = l
        else
            largest = i
        end

        # check if the right child exists and is greater than i
        if r <= size && A[r] > A[largest]
            largest = r
        end

        # if one of the children is greater than the node A[i], swap them and
        # recurse
        if largest != i
            A[i], A[largest] = A[largest], A[i]
            maxheapify(A, largest, size)
        end
    end

    # Reference: Cormen, p. 157
    # The modification made here is that the heap size is calculated and
    # returned instead of being made an attribute of A. This is to avoid making
    # new types.
    function buildmaxheap{T<:Real}(A::Array{T, 1})
        size = length(A)
        for i = Int(floor(length(A)/2)):-1:1
            maxheapify(A, i, size) end
        size
    end

    # The actual heap sort is down here:
    size = buildmaxheap(A)
    for i = length(A):-1:2
        A[1], A[i] = A[i], A[1]
        size -= 1
        maxheapify(A, 1, size)
    end
end

function quicksort{T<:Real}(A::Array{T,1})

    # Reference: Cormen, p. 171
    function partition{T<:Real}(A::Array{T,1}, p::Int, r::Int)
        x = A[r]
        i = p - 1
        for j = p:-1:(r-1)
            if A[j] <= x
                i += 1
                A[i], A[j] = A[j], A[i]
            end
        end
        A[i+1], A[r] = A[r], A[i+1]

        i + 1
    end

    # Reference: Cormen, p. 179
    function randomizedpartition{T<:Real}(A::Array{T,1}, p::Int, r::Int)
        i = rand(p:r)
        A[r], A[i] = A[i], A[r]
        partition(A, p, r)
    end

    # Reference: Cormen, p. 171
    function kernel{T<:Real}(A::Array{T,1}, p::Int, r::Int)
        if p < r
            q = randomizedpartition(A, p, r)
            kernel(A, p, q - 1)
            kernel(A, q + 1, r)
        end
    end

    kernel(A, 1, length(A))
end

function countingsort{T<:Real}(A::Array{T,1}, B::Array{T,1}, k::Int)
    C = zeros(k+1)
    for j = 1:length(A)
        C[A[j]] += 1
    end
    for i = 2:(k+1)
        C[i] = C[i] + C[i-1]
    end
    for j = length(A):-1:1
        B[C[A[j]]] = A[j]
        C[A[j]] = C[A[j]] - 1
    end
end


text:
# extract all sites
const LIQUID_STR = (:liquid,:LIQUID,:L,:l)

"""
    is_liquid(x::Union{Symbol,String})

Returns `true` if the symbol is in `(:liquid,:LIQUID,:L,:l)`.

If a string is passed, it is converted to symbol.
"""
is_liquid(sym::Symbol) = sym in LIQUID_STR
is_liquid(str::String) = is_liquid(Symbol(str))

const VAPOUR_STR = (:vapor,:VAPOR,:VAPOUR,:vapour,:g,:G,:v,:V,:gas,:GAS)

"""
    is_vapour(x::Union{Symbol,String})

Returns `true` if the symbol is in `(:vapor,:VAPOR,:VAPOUR,:vapour,:g,:G,:v,:V,:gas,:GAS)`.

If a string is passed, it is converted to symbol.
"""
is_vapour(sym::Symbol) = sym in VAPOUR_STR
is_vapour(str::String) = is_vapour(Symbol(str))

const SUPERCRITICAL_STR = (:sc,:SC,:supercritical,:SUPERCRITICAL)

"""
    is_supercritical(x::Union{Symbol,String})

Returns `true` if the symbol is in `(:sc,:SC,:supercritical,:SUPERCRITICAL)`.

If a string is passed, it is converted to symbol.
"""
is_supercritical(sym::Symbol) = sym in SUPERCRITICAL_STR
is_supercritical(str::String) = is_vapour(Symbol(str))

"""
    ∑(iterator)

equivalent to `sum(iterator,init=0.0)`.

"""
function ∑(iterator) #not collecting is faster
    #_0 = ifelse(length(iterator)>0,zero(first(iterator)),0.0)
    _0 = 0.0
    return reduce(Base.add_sum,iterator,init=_0)
end

"""
    ∑(fn,iterator)

equivalent to `sum(fniterator,init=0.0)`.

"""
function ∑(fn,iterator) #not collecting is faster
    #_0 = ifelse(length(iterator)>0,zero(first(iterator)),0.0)
    _0 = 0.0
    return mapreduce(fn,Base.add_sum,iterator,init=_0)
end

"""
    xlogx(x::Number)
Return `x * log(x)` for `x ≥ 0`, handling ``x = 0`` by taking the downward limit.

copied from LogExpFunctions.jl
"""
function xlogx(x::Number)
    result = x * log(x)
    ifelse(iszero(x), zero(result), result)
end

struct FractionVector{T,V} <: AbstractVector{T}
    vec::V
    val::T
end
#=

Fraction Vector
useful when expressing fractions in n-1 components.
the last component is calculated at build time.
it allocates less than creating a new vector or appending.
=#
##
function FractionVector(v::AbstractVector)
    a = one(eltype(v))
    any(x->x<0,v) && throw(DomainError(v,"all elements of a fraction vector should be positive."))
    a -=sum(v)
    a < 0 && throw(DomainError(a,"the values of the input vector add to more than one"))
    return FractionVector(v,a)
end

function FractionVector(v::Real)
    a = one(v)
    (v < zero(v)) && throw(DomainError(v,"all elements of a fraction vector should be positive."))
    a -= v
    a < 0 && throw(DomainError(a,"the values of the input vector add to more than one"))
    return FractionVector(v,a)
end

@inline Base.eltype(v::FractionVector{T}) where T = T
@inline Base.length(v::FractionVector)::Int = Base.length(v.vec) + 1

@inline function Base.length(v::FractionVector{T,<:Real})::Int where T
    return 2
end

@inline Base.size(v::FractionVector) = (length(v),)
@inline function Base.getindex(v::FractionVector,i::Int)
    @boundscheck checkbounds(v, i)
    return ifelse(length(v)==i,v.val,v.vec[min(length(v.vec),i)])
end

@inline function Base.getindex(v::FractionVector{T,<:Real},i::Int) where T
    @boundscheck checkbounds(v, i)
    return ifelse(i==1,v.vec,v.val)
end

Base.IndexStyle(::Type{<:FractionVector}) = IndexLinear()

unicode_subscript(a::Int) = join('₀'+d for d in reverse(digits(a)))


text:
# ------------------- #
# BasisMatrix Type #
# ------------------- #

abstract type AbstractBasisMatrixRep end
const ABSR = AbstractBasisMatrixRep

struct Tensor <: ABSR end
struct Direct <: ABSR end
struct Expanded <: ABSR end

mutable struct BasisMatrix{BST<:ABSR, TM<:AbstractMatrix}
    order::Matrix{Int}
    vals::Matrix{TM}
end

Base.eltype(bm::BasisMatrix{BST,TM}) where {BST, TM} = eltype(TM)

Base.show(io::IO, b::BasisMatrix{BST}) where {BST} =
    print(io, "BasisMatrix{$BST} of order $(b.order)")

Base.ndims(bs::BasisMatrix) = size(bs.order, 2)

# not the same if either type parameter is different
function ==(::BasisMatrix{BST1}, ::BasisMatrix{BST2}) where {BST1<:ABSR,BST2<:ABSR}
    false
end

function ==(::BasisMatrix{BST,TM1},
            ::BasisMatrix{BST,TM2}) where {BST<:ABSR,TM1<:AbstractMatrix,TM2<:AbstractMatrix}
    false
end

# if type parameters are the same, then it is the same if all fields are the
# same
function ==(b1::BasisMatrix{BST,TM},
            b2::BasisMatrix{BST,TM}) where {BST<:ABSR,TM<:AbstractMatrix}
    b1.order == b2.order && b1.vals == b2.vals
end

# -------------- #
# Internal Tools #
# -------------- #

@inline function _checkx(N, x::AbstractMatrix)
    size(x, 2) != N && error("Basis is $N dimensional, x must have $N columns")
    x
end

@inline function _checkx(N, x::AbstractVector{T}) where T
    # if we have a 1d basis, we can evaluate at each point
    if N == 1
        return x
    end

    # If Basis is > 1d, one evaluation point and reshape to (1,N) if possible...
    if length(x) == N
        return reshape(x, 1, N)
    end

    # ... or throw an error
    error("Basis is $N dimensional, x must have $N elements")
end

@inline function _checkx(N, x::TensorX)
    # for BasisMatrix{Tensor} family. Need one vector per dimension
    if length(x) == N
        return x
    end

    # otherwise throw an error
    error("Basis is $N dimensional, need one Vector per dimension")
end

# common checks for all convert methods
function check_convert(bs::BasisMatrix, order::Matrix)
    d = ndims(bs)
    numbas, d1 = size(order)

    d1 != d && error("ORDER incompatible with basis functions")  # 35-37

    # 39-41
    if any(minimum(order, dims=1) .< bs.order)
        error("Order of derivative operators exceeds basis")
    end
    return d, numbas, d1
end

"""
Do common transformations to all constructor of `BasisMatrix`

##### Arguments

- `N::Int`: The number of dimensions in the corresponding `Basis`
- `x::AbstractArray`: The points for which the `BasisMatrix` should be
constructed
- `order::Array{Int}`: The order of evaluation for each dimension of the basis

##### Returns

- `m::Int`: the total number of derivative order basis functions to compute.
This will be the number of rows in the matrix form of `order`
- `order::Matrix{Int}`: A `m × N` matrix that, for each of the `m` desired
specifications, gives the derivative order along all `N` dimensions
- `minorder::Matrix{Int}`: A `1 × N` matrix specifying the minimum desired
derivative order along each dimension
- `numbases::Matrix{Int}`: A `1 × N` matrix specifying the total number of
distinct derivative orders along each dimension
- `x::AbstractArray`: The properly transformed points at which to evaluate
the basis

"""
function check_basis_structure(N::Int, x, order)
    order = _check_order(N, order)

    # initialize basis structure (66-74)
    m = size(order, 1)  # by this time order is a matrix
    if m > 1
        minorder = minimum(order, dims=1)
        numbases = (maximum(order, dims=1) .- minorder) .+ 1
    else
        minorder = order + zeros(Int, 1, N)
        numbases = fill(1, 1, N)
    end

    x = _checkx(N, x)

    return m, order, minorder, numbases, x
end

# --------------- #
# convert methods #
# --------------- #

# no-op. Don't worry about the order argument.
Base.convert(::Type{T}, bs::BasisMatrix{T}, order=bs.order) where {T<:ABSR} = bs

# funbconv from direct to expanded
function Base.convert(::Type{Expanded}, bs::BasisMatrix{Direct,TM},
                  order=fill(0, 1, size(bs.order, 2))) where TM
    d, numbas, d1 = check_convert(bs, order)

    vals = Array{TM}(undef, numbas, 1)

    for i=1:numbas
        vals[i] = bs.vals[order[i, d] - bs.order[d]+1, d]  # 63
        for j=d-1:-1:1
            vals[i] = row_kron(vals[i],
                               bs.vals[order[i, j] - bs.order[j]+1, j])  #65
        end
    end
    BasisMatrix{Expanded,TM}(order, vals)
end

# funbconv from tensor to expanded
function Base.convert(::Type{Expanded}, bs::BasisMatrix{Tensor,TM},
                  order=fill(0, 1, size(bs.order, 2))) where TM
    d, numbas, d1 = check_convert(bs, order)

    vals = Array{TM}(undef, numbas, 1)

    for i=1:numbas  # 54
        vals[i] = bs.vals[order[i, d] - bs.order[d]+1, d]  # 55
        for j=d-1:-1:1  # 56
            # 57
            vals[i] = kron(vals[i], bs.vals[order[i, j] - bs.order[j]+1, j])
        end
    end

    BasisMatrix{Expanded,TM}(order, vals)
end

# funbconv from tensor to direct
# HACK: there is probably a more efficient way to do this, but since I don't
#       plan on doing it much, this will do for now. The basic point is that
#       we need to expand the rows of each element of `vals` so that all of
#       them have prod([size(v, 1) for v in bs.vals])) rows.
function Base.convert(::Type{Direct}, bs::BasisMatrix{Tensor,TM},
                      order=fill(0, 1, size(bs.order, 2))) where TM
    d, numbas, d1 = check_convert(bs, order)
    vals = Array{TM}(undef, numbas, d)
    raw_ind = Array{Vector{Int}}(undef, d)

    for j=1:d
        for i=1:size(bs.vals, 1)
            if !(isempty(bs.vals[i, j]))  # 84
                raw_ind[j] = collect(1:size(bs.vals[i, j], 1))  # 85
                break
            end
        end
    end

    ind = gridmake(raw_ind...)  # 90

    for j=1:d, i=1:numbas
        if !(isempty(bs.vals[i, j]))  # 93
            vals[i, j] = bs.vals[i, j][ind[:, j], :]  # 94
        end
    end

    BasisMatrix{Direct,TM}(order, vals)
end

# ------------ #
# Constructors #
# ------------ #

# method to construct BasisMatrix in direct or expanded form based on
# a matrix of `x` values  -- funbasex
function BasisMatrix(::Type{T2}, basis::Basis{N,BF}, ::Direct,
                     x::AbstractArray=nodes(basis)[1], order=0) where {N,BF,T2}
    m, order, minorder, numbases, x = check_basis_structure(N, x, order)
    # 76-77
    out_order = minorder
    out_format = Direct()
    val_type = bmat_type(T2, basis, x)
    vals = Array{val_type}(undef, maximum(numbases), N)

    # now do direct form, will convert to expanded later if needed
    for j=1:N
        # 126-130
        if (m > 1)
            orderj = unique(order[:, j])
        else
            orderj = order[1, j]
        end

        #131-135
        if length(orderj) == 1
            vals[1, j] = evalbase(T2, basis.params[j], x[:, j], orderj[1])
        else
            vals[orderj-minorder[j]+1, j] =
                evalbase(T2, basis.params[j], x[:, j], orderj)
        end
    end

    # construct Direct Format
    BasisMatrix{Direct,val_type}(out_order, vals)
end

function BasisMatrix(::Type{T2}, basis::Basis, ::Expanded,
                     x::AbstractArray=nodes(basis)[1], order=0) where T2  # funbasex
    # create direct form, then convert to expanded
    bsd = BasisMatrix(T2, basis, Direct(), x, order)
    convert(Expanded, bsd, bsd.order)
end

function BasisMatrix(::Type{T2}, basis::Basis{N,BT}, ::Tensor,
                     x::TensorX=nodes(basis)[2], order=0) where {N,BT,T2}

    m, order, minorder, numbases, x = check_basis_structure(N, x, order)
    out_order = minorder
    out_format = Tensor()
    val_type = bmat_type(T2, basis, x[1])  # TODO: reduce(promot_type, eltype.(x))
    vals = Array{val_type}(undef, maximum(numbases), N)

    # construct tensor base
    for j in 1:N
        # 113-117
        if (m > 1)
            orderj = unique(order[:, j])
        else
            orderj = order[1, j]
        end

        #118-122
        if length(orderj) == 1
            vals[1, j] = evalbase(T2, basis.params[j], x[j], orderj[1])
        else
            vals[orderj-minorder[j]+1, j] = evalbase(T2, basis.params[j], x[j], orderj)
        end
    end

    BasisMatrix{Tensor,val_type}(out_order, vals)
end

# When the user doesn't supply a ABSR, we pick one for them.
# for x::AbstractMatrix we pick direct
# for x::TensorX we pick Tensor
function BasisMatrix(::Type{T2}, basis::Basis, x::AbstractArray, order=0) where T2
    BasisMatrix(T2, basis, Direct(), x, order)
end

function BasisMatrix(::Type{T2}, basis::Basis, x::TensorX, order=0) where T2
    BasisMatrix(T2, basis, Tensor(), x, order)
end

# method to allow passing types instead of instances of ABSR
function BasisMatrix(::Type{T2}, basis, ::Type{BST},
                     x::Union{AbstractArray,TensorX}, order=0) where {BST<:ABSR,T2}
    BasisMatrix(T2, basis, BST(), x, order)
end

function BasisMatrix(basis, ::Type{BST},
                     x::Union{AbstractArray,TensorX}, order=0) where BST<:ABSR
    BasisMatrix(basis, BST(), x, order)
end

# method without vals eltypes
function BasisMatrix(basis::Basis, tbm::TBM,
                     x::Union{AbstractArray,TensorX}, order=0) where TBM<:ABSR
    BasisMatrix(Nothing, basis, tbm, x, order)
end

function BasisMatrix(basis::Basis, x::Union{AbstractArray,TensorX}, order=0)
    BasisMatrix(Nothing, basis, x, order)
end

# method without x
function BasisMatrix(basis::Basis, tbm::Union{Type{TBM},TBM}) where TBM<:ABSR
    BasisMatrix(Nothing, basis, tbm)
end


text:
function ≒(a,b)
    if a ≈ b
        return true
    elseif abs(a-b) < 1e-10
        return true
    else
        return false
    end
end

@testset "BSplineBasis" begin

    Δt = 1.0e-8
    ε = 1.0e-8
    ∞ = Inf

    @testset "0th degree basis" begin
        Random.seed!(42)

        p = 0
        v = rand(10)
        k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
        P = BSplineSpace(p, k)
        fP = FastBSplineSpace(p, k)
        n = dim(P)

        @test degree(P) == degree(fP) == 0

        # Check the values of B-spline basis funcitons
        for i in 1:n, t in rand(5)
            for t in rand(5)
                @test (bsplinebasis₊₀(P, i, t)
                     ≈ bsplinebasis₊₀(fP, i, t)
                     ≈ bsplinebasis₋₀(P, i, t)
                     ≈ bsplinebasis₋₀(fP, i, t)
                     ≈ bsplinebasis(P, i, t)
                     ≈ bsplinebasis(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
            end
        end

        # Check the values of the derivative of B-spline basis funcitons
        for i in 1:n
            for t in rand(5)
                @test (bsplinebasis′₊₀(P, i, t)
                     ≈ bsplinebasis′₊₀(fP, i, t)
                     ≈ bsplinebasis′₋₀(P, i, t)
                     ≈ bsplinebasis′₋₀(fP, i, t)
                     ≈ bsplinebasis′(P, i, t)
                     ≈ bsplinebasis′(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
            end
        end

        # Check the derivative
        @testset "derivative" begin
            Random.seed!(57)
            for i in 1:n, t in rand(5)
                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
                b′   = bsplinebasis′(fP, i, t)
                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
                @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
            end
        end
    end

    @testset "1st degree basis" begin
        Random.seed!(42)

        p = 1
        v = rand(10)
        k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
        P = BSplineSpace(p, k)
        fP = FastBSplineSpace(p, k)
        n = dim(P)

        @test degree(P) == degree(fP) == 1

        # Check the values of B-spline basis funcitons
        for i in 1:n, t in rand(5)
            for t in rand(5)
                @test (bsplinebasis₊₀(P, i, t)
                     ≈ bsplinebasis₊₀(fP, i, t)
                     ≈ bsplinebasis₋₀(P, i, t)
                     ≈ bsplinebasis₋₀(fP, i, t)
                     ≈ bsplinebasis(P, i, t)
                     ≈ bsplinebasis(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
            end
        end

        # Check the values of the derivative of B-spline basis funcitons
        for i in 1:n
            for t in rand(5)
                @test (bsplinebasis′₊₀(P, i, t)
                     ≈ bsplinebasis′₊₀(fP, i, t)
                     ≈ bsplinebasis′₋₀(P, i, t)
                     ≈ bsplinebasis′₋₀(fP, i, t)
                     ≈ bsplinebasis′(P, i, t)
                     ≈ bsplinebasis′(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
            end
        end

        # Check the derivative
        @testset "derivative" begin
            Random.seed!(57)
            for i in 1:n, t in rand(5)
                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
                b′   = bsplinebasis′(fP, i, t)
                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
                @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
            end
        end
    end

    @testset "2nd degree basis" begin
        Random.seed!(42)

        p = 2
        v = rand(10)
        k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
        P = BSplineSpace(p, k)
        fP = FastBSplineSpace(p, k)
        n = dim(P)

        @test degree(P) == degree(fP) == 2

        # Check the values of B-spline basis funcitons
        for i in 1:n, t in rand(5)
            for t in rand(5)
                @test (bsplinebasis₊₀(P, i, t)
                     ≈ bsplinebasis₊₀(fP, i, t)
                     ≈ bsplinebasis₋₀(P, i, t)
                     ≈ bsplinebasis₋₀(fP, i, t)
                     ≈ bsplinebasis(P, i, t)
                     ≈ bsplinebasis(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
            end
        end

        # Check the values of the derivative of B-spline basis funcitons
        for i in 1:n
            for t in rand(5)
                @test (bsplinebasis′₊₀(P, i, t)
                     ≈ bsplinebasis′₊₀(fP, i, t)
                     ≈ bsplinebasis′₋₀(P, i, t)
                     ≈ bsplinebasis′₋₀(fP, i, t)
                     ≈ bsplinebasis′(P, i, t)
                     ≈ bsplinebasis′(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
            end
        end

        # Check the derivative
        @testset "derivative" begin
            Random.seed!(57)
            for i in 1:n, t in rand(5)
                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
                b′   = bsplinebasis′(fP, i, t)
                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
                @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
            end
        end
    end

    @testset "3rd degree basis" begin
        Random.seed!(42)

        p = 3
        v = rand(10)
        k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
        P = BSplineSpace(p, k)
        fP = FastBSplineSpace(p, k)
        n = dim(P)

        @test degree(P) == degree(fP) == 3

        # Check the values of B-spline basis funcitons
        for i in 1:n, t in rand(5)
            for t in rand(5)
                @test (bsplinebasis₊₀(P, i, t)
                     ≈ bsplinebasis₊₀(fP, i, t)
                     ≈ bsplinebasis₋₀(P, i, t)
                     ≈ bsplinebasis₋₀(fP, i, t)
                     ≈ bsplinebasis(P, i, t)
                     ≈ bsplinebasis(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
            end
        end

        # Check the values of the derivative of B-spline basis funcitons
        for i in 1:n
            for t in rand(5)
                @test (bsplinebasis′₊₀(P, i, t)
                     ≈ bsplinebasis′₊₀(fP, i, t)
                     ≈ bsplinebasis′₋₀(P, i, t)
                     ≈ bsplinebasis′₋₀(fP, i, t)
                     ≈ bsplinebasis′(P, i, t)
                     ≈ bsplinebasis′(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
            end
        end

        # Check the derivative
        @testset "derivative" begin
            Random.seed!(57)
            for i in 1:n, t in rand(5)
                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
                b′   = bsplinebasis′(fP, i, t)
                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
                @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
            end
        end
    end

    @testset "4th degree basis" begin
        Random.seed!(42)

        p = 4
        v = rand(10)
        k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
        P = BSplineSpace(p, k)
        fP = FastBSplineSpace(p, k)
        n = dim(P)

        @test degree(P) == degree(fP) == 4

        # Check the values of B-spline basis funcitons
        for i in 1:n, t in rand(5)
            for t in rand(5)
                @test (bsplinebasis₊₀(P, i, t)
                     ≈ bsplinebasis₊₀(fP, i, t)
                     ≈ bsplinebasis₋₀(P, i, t)
                     ≈ bsplinebasis₋₀(fP, i, t)
                     ≈ bsplinebasis(P, i, t)
                     ≈ bsplinebasis(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
            end
        end

        # Check the values of the derivative of B-spline basis funcitons
        for i in 1:n
            for t in rand(5)
                @test (bsplinebasis′₊₀(P, i, t)
                     ≈ bsplinebasis′₊₀(fP, i, t)
                     ≈ bsplinebasis′₋₀(P, i, t)
                     ≈ bsplinebasis′₋₀(fP, i, t)
                     ≈ bsplinebasis′(P, i, t)
                     ≈ bsplinebasis′(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
            end
        end

        # Check the derivative
        @testset "derivative" begin
            Random.seed!(57)
            for i in 1:n, t in rand(5)
                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
                b′   = bsplinebasis′(fP, i, t)
                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
                @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
            end
        end
    end

    @testset "5th degree basis" begin
        Random.seed!(42)

        p = 5
        v = rand(10)
        k = Knots(v) + (p+1)*Knots(0, 1) + Knots(v[2], v[3])
        P = BSplineSpace(p, k)
        fP = FastBSplineSpace(p, k)
        n = dim(P)

        @test degree(P) == degree(fP) == 5

        # Check the values of B-spline basis funcitons
        for i in 1:n, t in rand(5)
            for t in rand(5)
                @test (bsplinebasis₊₀(P, i, t)
                     ≈ bsplinebasis₊₀(fP, i, t)
                     ≈ bsplinebasis₋₀(P, i, t)
                     ≈ bsplinebasis₋₀(fP, i, t)
                     ≈ bsplinebasis(P, i, t)
                     ≈ bsplinebasis(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis₊₀(P, i, t) ≈ bsplinebasis₊₀(fP, i, t) ≒ bsplinebasis₊₀(P, i, t₊₀) ≒ bsplinebasis₊₀(fP, i, t₊₀)
                @test bsplinebasis₋₀(P, i, t) ≈ bsplinebasis₋₀(fP, i, t) ≒ bsplinebasis₋₀(P, i, t₋₀) ≒ bsplinebasis₋₀(fP, i, t₋₀)
                @test bsplinebasis(P, i, t) ≈ bsplinebasis(fP, i, t)
            end
        end

        # Check the values of the derivative of B-spline basis funcitons
        for i in 1:n
            for t in rand(5)
                @test (bsplinebasis′₊₀(P, i, t)
                     ≈ bsplinebasis′₊₀(fP, i, t)
                     ≈ bsplinebasis′₋₀(P, i, t)
                     ≈ bsplinebasis′₋₀(fP, i, t)
                     ≈ bsplinebasis′(P, i, t)
                     ≈ bsplinebasis′(fP, i, t))
            end
            for t in k
                t₊₀ = nextfloat(t)
                t₋₀ = prevfloat(t)
                @test bsplinebasis′₊₀(P, i, t) ≈ bsplinebasis′₊₀(fP, i, t) ≒ bsplinebasis′₊₀(P, i, t₊₀) ≒ bsplinebasis′₊₀(fP, i, t₊₀)
                @test bsplinebasis′₋₀(P, i, t) ≈ bsplinebasis′₋₀(fP, i, t) ≒ bsplinebasis′₋₀(P, i, t₋₀) ≒ bsplinebasis′₋₀(fP, i, t₋₀)
                @test bsplinebasis′(P, i, t) ≈ bsplinebasis′(fP, i, t)
            end
        end

        # Check the derivative
        @testset "derivative" begin
            Random.seed!(57)
            for i in 1:n, t in rand(5)
                b′₊₀ = bsplinebasis′₊₀(fP, i, t)
                b′₋₀ = bsplinebasis′₋₀(fP, i, t)
                b′   = bsplinebasis′(fP, i, t)
                a′₊₀ = (bsplinebasis₊₀(fP, i, t+Δt) - bsplinebasis₊₀(fP, i, t-Δt)) / 2Δt
                a′₋₀ = (bsplinebasis₋₀(fP, i, t+Δt) - bsplinebasis₋₀(fP, i, t-Δt)) / 2Δt
                a′   = (bsplinebasis(fP, i, t+Δt) - bsplinebasis(fP, i, t-Δt)) / 2Δt
                @test b′₊₀ ≈ b′₋₀ ≈ b′ ≈ a′₊₀ ≈ a′₋₀ ≈ a′
            end
        end
    end

    @testset "_bsplinebasis" begin
        Random.seed!(42)

        k = Knots(rand(10).-1) + Knots(rand(10)) + Knots(rand(10).+1)
        ts = rand(10)

        for p in 0:BasicBSpline.MAX_DEGREE
            fP = FastBSplineSpace(p,k)
            P = BSplineSpace(p,k)
            for t in ts
                j = BasicBSpline._knotindex(fP,t)
                _B = BasicBSpline._bsplinebasis(fP,t,j)
                # _B′ = BasicBSpline._bsplinebasis′(fP,t,j)
                B = [bsplinebasis(fP,i,t) for i in j-p:j]
                B′ = [bsplinebasis′(fP,i,t) for i in j-p:j]
                @test norm(collect(_B) - B) < ε

                _B = bsplinebasis(P,t)[j-p:j]
                @test norm(_B - B) < ε

                _B = bsplinebasis₊₀(P,t)[j-p:j]
                @test norm(_B - B) < ε

                _B = bsplinebasis₋₀(P,t)[j-p:j]
                @test norm(_B - B) < ε

                _B′ = bsplinebasis′(P,t)[j-p:j]
                @test norm(_B′ - B′) < ε

                _B′ = bsplinebasis′₊₀(P,t)[j-p:j]
                @test norm(_B′ - B′) < ε

                _B′ = bsplinebasis′₋₀(P,t)[j-p:j]
                @test norm(_B′ - B′) < ε
            end
        end
    end
end


text:
# This file is a part of SimilaritySearch.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

import Base:
    push!, popfirst!, pop!, length, last, first, empty!

export Item, KnnResult, maxlength, covrad

struct Item
    objID::Int
    dist::Float64
end

mutable struct KnnResult #{T}
    k::Int32
    pool::Vector{Item}
    #    pool::Vector{Item{T}}

    function KnnResult(k::Integer)
        v = Vector{Item}()
        sizehint!(v, k)
        new(k, v)
    end

end

"""
    fix_order!(res::KnnResult)

Fixes the sorted state of the array. It implements a kind of insertion sort
It is efficient due to the expected distribution of the items being inserted
(few smaller than the ones already inside)
"""
@inline function fix_order!(res::KnnResult)
    arr = res.pool
    item = arr[end]
    i = length(arr)
    @inbounds while i > 1
        if item.dist < arr[i-1].dist
            arr[i] = arr[i-1]
        else
            arr[i] = item
            return
        end
        i -= 1
    end
    arr[1] = item
end

"""
    push!(p::KnnResult, objID::Integer, dist::AbstractFloat)

Appends an item into the result set
"""
#push!(p::KnnResult{Int64}, objID::I, dist::F) where {I <: Union{Int32,Int16}, F <: Real} = push!(p, convert(Int64, objID), convert(Float64, dist))
#push!(p::KnnResult{T}, objID::T, dist::F) where {T, F <: Union{Float16, Float32}} = push!(p, objID, convert(Float64, dist))

function push!(p::KnnResult, objID::Integer, dist::AbstractFloat)
    if length(p.pool) < p.k
        # fewer items than the maximum capacity
        push!(p.pool, Item(objID, dist))
        fix_order!(p)
        return true
    end

    if dist >= last(p).dist
        # p.k == length(p.pool) but item.dist doesn't improve the pool's radius
        return false
    end

    # p.k == length(p.pool) but item.dist improves the result set
    @inbounds p.pool[end] = Item(objID, dist)
    fix_order!(p)
    true
end

"""
    first(p::KnnResult)

Return the first item of the result set, the closest item
"""
function first(p::KnnResult)
    p.pool[1]
end

"""
    last(p::KnnResult)

Returns the last item of the result set
"""
function last(p::KnnResult)
    p.pool[end]
end

"""
    popfirst!(p::KnnResult)

Removes and returns the nearest neeighboor from the pool, an O(length(p.pool)) operation
"""
function popfirst!(p::KnnResult)
    popfirst!(p.pool)
end

"""
    pop!(p)

Removes and returns the last item in the pool, it is an O(1) operation
"""
function pop!(p::KnnResult)
    pop!(p.pool)
end

"""
    length(p::KnnResult)

length returns the number of items in the result set
"""
Base.length(p::KnnResult) = length(p.pool)

"""
    maxlength(p::KnnResult)

The maximum allowed cardinality (the k of knn)
"""
maxlength(p::KnnResult) = p.k

"""
    covrad(p::KnnResult)::Float64

Returns the coverage radius of the result set; if length(p) < K then typemax(Float32) is returned
"""
function covrad(p::KnnResult)::Float64
    return length(p.pool) < p.k ? typemax(Float32) : last(p).dist
end

"""
    empty!(p::KnnResult)

Clears the content of the result pool
"""
function empty!(p::KnnResult)
    resize!(p.pool, 0)
end

##### iterator interface
### KnnResult
function Base.iterate(p::KnnResult)
    return length(p) == 0 ? nothing : (first(p), 2)
end

function Base.iterate(p::KnnResult, state::Int)
    if state > length(p)
        return nothing
    end

    return p.pool[state], state + 1
end


text:
# Converted the code for airspace encounter model from MATLAB to Julia

# Originally written in MATLAB by Mykel Kochenderfer, mykel@stanford.edu
# Converted by Youngjun Kim, youngjun@stanford.edu
# Date: 4/23/2014

function em_sample(aem, num_transition_samples; initial_dist = nothing)

# EM_SAMPLE Outputs samples from an encounter model to files.
#   Outputs samples into two specified files from an encounter model
#   described in a file.
#
#   EM_SAMPLES takes as input the following arguments:
#   PARAMETERS_FILENAME: a string specifying the name of the parameters
#   file
#   INITIAL_OUTPUT_FILENAME: a string specifying the the name of the file
#   to store transition samples
#   NUM_INITIAL_SAMPLES: the number of samples to generate
#   NUM_TRANSITION_SAMPLES: the number of steps to sample from the
#   transition network

    p = aem.parameters

    # read parameters
    #p = Parameters()
    #em_read(p, parameters_filename)

    # create priors
    dirichlet_initial = bn_dirichlet_prior(p.N_initial)
    dirichlet_transition = bn_dirichlet_prior(p.N_transition)

    num_transition_samples  = Int(num_transition_samples)

    output = Array(AbstractFloat, (num_transition_samples, p.n_initial + 2))

    if initial_dist == nothing
        x = create_sample(p, dirichlet_initial, dirichlet_transition, num_transition_samples)

        for j = 1:p.n_initial
            if !isempty(p.boundaries[j])
                if x[j, 1] < p.boundaries[j][1] || x[j, 1] > p.boundaries[j][end]
                    print("sample: ", x[:, 1]')
                    println("i: ", j, ", value: ", x[j, 1], ", boundary: ", p.boundaries[j][end])
                    error("The value violates the boundary.")
                end
            elseif !(Int(x[j, 1]) in collect(1:p.r_initial[j]))
                print("sample: ", x[:, 1]')
                println("i: ", j, ", value: ", x[j, 1], ", boundary: ", p.r_initial[j])
                error("The value violates the boundary.")
            end
        end
    else
        x, IS = create_sample(p, dirichlet_initial, dirichlet_transition, num_transition_samples, initial_dist = initial_dist)
    end

    output[1:num_transition_samples, 1] = ones(num_transition_samples)
    output[1:num_transition_samples, 2] = collect(0:num_transition_samples - 1)
    output[1:num_transition_samples, 3:end] = x'

    if initial_dist == nothing
        return output
    else
        return output, IS
    end
end

function em_sample_n(aem, num_initial_samples, num_transition_samples)

# EM_SAMPLE Outputs samples from an encounter model to files.
#   Outputs samples into two specified files from an encounter model
#   described in a file.
#
#   EM_SAMPLES takes as input the following arguments:
#   PARAMETERS_FILENAME: a string specifying the name of the parameters
#   file
#   INITIAL_OUTPUT_FILENAME: a string specifying the the name of the file
#   to store transition samples
#   NUM_INITIAL_SAMPLES: the number of samples to generate
#   NUM_TRANSITION_SAMPLES: the number of steps to sample from the
#   transition network

    p = aem.parameters

    # read parameters
    #p = Parameters()
    #em_read(p, parameters_filename)

    # create priors
    dirichlet_initial = bn_dirichlet_prior(p.N_initial)
    dirichlet_transition = bn_dirichlet_prior(p.N_transition)

    num_initial_samples = Int(num_initial_samples)
    num_transition_samples  = Int(num_transition_samples)

    output_array = Array(AbstractFloat, (num_initial_samples * num_transition_samples, p.n_initial + 2))

    for i = 1:num_initial_samples
        x = create_sample(p, dirichlet_initial, dirichlet_transition, num_transition_samples)

        for j = 1:p.n_initial
            if !isempty(p.boundaries[j])
                if x[j, 1] < p.boundaries[j][1] || x[j, 1] > p.boundaries[j][end]
                    print("sample: ", i, " ", x[:, 1]')
                    println("i: ", j, ", value: ", x[j, 1], ", boundary: ", p.boundaries[j][end])
                    error("The value violates the boundary.")
                end
            elseif !(Int(x[j, 1]) in [1:p.r_initial[j]])
                print("sample: ", i, " ", x[:, 1]')
                println("i: ", j, ", value: ", x[j, 1], ", boundary: ", p.r_initial[j])
                error("The value violates the boundary.")
            end
        end

        output_array[((i - 1) * num_transition_samples + 1):(i * num_transition_samples), 1] = i * ones(num_transition_samples)
        output_array[((i - 1) * num_transition_samples + 1):(i * num_transition_samples), 2] = [0:num_transition_samples - 1]
        output_array[((i - 1) * num_transition_samples + 1):(i * num_transition_samples), 3:end] = x'
    end

    return output_array
end

function create_sample(p, dirichlet_initial, dirichlet_transition, sample_time; initial_dist = nothing)

    if initial_dist == nothing
        initial, events = dbn_sample(p.G_initial, p.G_transition, p.temporal_map, p.r_transition, p.N_initial, p.N_transition, dirichlet_initial, dirichlet_transition, p.boundaries, sample_time)
    else
        initial, events, IS = dbn_sample(p.G_initial, p.G_transition, p.temporal_map, p.r_transition, p.N_initial, p.N_transition, dirichlet_initial, dirichlet_transition, p.boundaries, sample_time, initial_dist = initial_dist)
    end

    if isempty(events)
        events = [sample_time 0 0]
    else
        events = [events; (sample_time - sum(events[:,1])) 0 0]
    end

    # Within-bin resampling
    events = resample_events(initial, events, p.resample_rates)

    # Dediscretize
    for i = 1:length(initial)
        if !isempty(p.boundaries[i])
            if initial_dist == nothing
                initial[i] = dediscretize(initial[i], p.boundaries[i], p.zero_bins[i])[1]
            else
                if IS[i, 1] == 0
                    initial[i] = dediscretize(initial[i], p.boundaries[i], p.zero_bins[i])[1]
                else
                    initial[i] = IS[i, 2]
                end
            end
        end
    end

    if !isempty(events)
        for i = 1:(size(events, 1) - 1)
            e2 = convert(Int, events[i, 2])
            events[i, 3] = dediscretize(events[i, 3], p.boundaries[e2], p.zero_bins[e2])[1]
        end
    end

    if initial_dist == nothing
        return events2samples(initial, events)
    else
        return events2samples(initial, events), IS
    end
end

function events2samples(initial, events)

    n = length(initial)
    D = zeros(n, Int(sum(events[:, 1])))

    x = vec(copy(initial))
    t = 0

    for i = 1:size(events, 1)
        event = events[i, :]

        delta_t = Int(event[1])

        if event[2] == 0
            if delta_t > 0
                t = t + 1
                D[:, t:(t + delta_t - 1)] = repmat(x, 1, delta_t)
            end
        else
            if delta_t > 0
                D[:, (t + 1):(t + delta_t)] = repmat(x, 1, delta_t)
                t = t + delta_t
            end

            e2 = convert(Int, event[2])
            x[e2] = event[3]
        end
    end

    return D
end

function dediscretize(d, parameters, zero_bins)

    if isempty(parameters)
        return d
    end

    x = zeros(length(d))
    n = length(d)

    for i = 1:n
        # generate uniform in interval [a, b]
        if d[i] == zero_bins
            x[i] = 0
        else
            dd = convert(Int, d[i])
            a = parameters[dd]
            b = parameters[dd + 1]
            x[i] = a + (b - a) * rand()
        end
    end

    return x
end

function resample_events(initial, events, rates)

    n = size(events, 1)

    newevents = Float64[]

    x = vec(copy(initial))

    for i = 1:n
        holdtime = events[i, 1]

        if holdtime == 0
            append!(newevents, vec(events[i, :]))
        else
            delta_t = 0

            for j = 1:holdtime
                changes = find(rand(size(rates)) .< rates)
                delta_t = delta_t + 1

                if !isempty(changes)
                    T = [[delta_t; zeros(length(changes) - 1)] changes x[changes]]
                    for k = 1:size(T, 1)
                        append!(newevents, vec(T[k, :]))
                    end

                    delta_t = 0
                end
            end

            append!(newevents, [delta_t, events[i,2], events[i,3]])
        end

        if events[i, 2] > 0
            e2 = convert(Int, events[i,2])
            x[e2] = events[i, 3]
        end
    end

    newevents = reshape(newevents, (3, Int(length(newevents)/3)))'

    return newevents
end

function dbn_sample(G_initial, G_transition, temporal_map, r, N_initial, N_transition, dirichlet_initial, dirichlet_transition, boundaries, t_max; initial_dist = nothing)

# DBN_SAMPLE Samples from a dynamic Bayesian network.
# Returns a sample of the intitial variables and a series of events from
# the dynamic Bayesian network with specified graphical structure,
# sufficient statistics, and prior.
#
#   DBN_SAMPLE takes as input the following
#   G_initial - initial distribution graph structure
#   G_transition - continuous transition model graph structure
#   r - a column vector specifying the number of bins for each variable
#   N_initial - sufficient statistics for initial distribution
#   N_transition - sufficient statistics for transition model
#   dirichlet_initial - Dirichlet prior
#   dirichlet_transition - Dirichlet prior
#   t_max - the maximum amount of time to run the simulation
#
#   DBN_SAMPLE returns the following
#   initial
#   events - matrix of (t, variable_index, new_value) (NOTE: t is the time
#   since the last event)

    if initial_dist == nothing
        initial = bn_sample(G_initial, r, N_initial, dirichlet_initial, boundaries, 1)
    else
        initial, IS = bn_sample(G_initial, r, N_initial, dirichlet_initial, boundaries, 1, initial_dist = initial_dist)
    end

    events = Float64[]
    n_initial = size(G_initial, 1)

    dynamic_variables = temporal_map[:, 2]

    order = bn_sort(G_transition)

    x = [initial zeros(1, length(dynamic_variables))]

    delta_t = 0

    for t = 2:t_max
        delta_t = delta_t + 1
        x_old = copy(x)

        for i = order
            if in(i, dynamic_variables)
                # only dynamic variables change
                parents = G_transition[:, i]

                j = 1
                if !isempty(find(parents))
                    j = asub2ind(r[parents], x[parents]')
                end

                x[i] = select_random(N_transition[i][:, j] + dirichlet_transition[i][:, j])
            end
        end

        # map back
        x[temporal_map[:, 1]] = x[temporal_map[:, 2]]

        if !isequal(x[1:n_initial], x_old[1:n_initial])
            # change (i.e. new event)
            for i = 1:n_initial
                if x[i] != x_old[i]
                    append!(events, vec([delta_t i x[i]]))
                    delta_t = 0
                end
            end
        end
    end

    events = reshape(events, (3, Int(length(events)/3)))'

    if initial_dist == nothing
        return initial, events
    else
        return initial, events, IS
    end
end

function bn_sample(G, r, N, alpha, boundaries, num_samples; initial_dist = nothing)

# BN_SAMPLE Produces a sample from a Bayesian network.
#   Returns a matrix whose rows consist of n-dimensional samples from the
#   specified Bayesian network.
#
#   S = BN_SAMPLE(G, R, N, ALPHA, NUM_SAMPLES) returns a matrix S whose
#   rows consist of n-dimensional samples from a Bayesian network with
#   graphical structure G, sufficient statistics N, and prior ALPHA. The
#   number of samples is specified by NUM_SAMPLES. The array R specifies
#   the number of bins associated with each variable.

    order = bn_sort(G)

    n = length(N)
    S = zeros((num_samples, n))

    if initial_dist != nothing
        IS = zeros(n, 4)
    end

    for sample_index = 1:num_samples
        # Generate each sample
        for i = order
            if initial_dist == nothing
                parents = G[:, i]

                j = 1
                if !isempty(find(parents))
                    j = asub2ind(r[[parents; falses(length(r) - size(G, 1))]], S[sample_index, parents])
                end

                S[sample_index, i] = select_random(N[i][:, j] + alpha[i][:, j])

            else
                f = initial_dist[i]

                if f == nothing
                    parents = G[:, i]

                    j = 1
                    if !isempty(find(parents))
                        j = asub2ind(r[[parents; falses(length(r) - size(G, 1))]], S[sample_index, parents])
                    end

                    S[sample_index, i] = select_random(N[i][:, j] + alpha[i][:, j])
                else
                    value, prob_new = f()

                    if !isempty(boundaries[i])
                        index = findfirst(x -> (x > value), boundaries[i]) - 1

                        if index == -1
                            index = r[i]
                        end
                    else
                        index = value
                    end

                    parents = G[:, i]

                    j = 1
                    if !isempty(find(parents))
                        j = asub2ind(r[[parents, falses(length(r) - size(G, 1))]], S[sample_index, parents])
                    end

                    A = N[i][:, j] + alpha[i][:, j]

                    if !isempty(boundaries[i])
                        prob = 1 / (boundaries[i][index + 1] - boundaries[i][index]) * A[index] / sum(A)
                    else
                        prob = A[index] / sum(A)
                    end

                    IS[i, :] = [index, value, prob_new, prob]

                    S[sample_index, i] = index
                end
            end
        end
    end

    if initial_dist == nothing
        return S
    else
        return S, IS
    end
end

function select_random(weights)
# SELECT_RANDOM Randomly selects an index according to specified weights.
#   Returns a randomly selected index according to the distribution
#   specified by a vector of weights.
#
#   INDEX = SELECT_RANDOM(WEIGHTS) returns a scalar index INDEX selected
#   randomly according to the specified weights WEIGHTS represented as an
#   array.

    s = cumsum(weights)
    r = s[end] * rand()
    index = findfirst(x -> (x >= r), s)

    return index
end

# Can't replace this with built-in sub2ind...
# NB: This needs to be checked for accuracy.
function asub2ind(siz, x)
# ASUB2IND Linear index from multiple subscripts.
#   Returns a linder index from multiple subscripts assuming a matrix of a
#   specified size.
#
#   NDX = ASUB2IND(SIZ,X) returns the linear index NDX of the element in a
#   matrix of dimension SIZ associated with subscripts specified in X.

     k = [1; cumprod(siz[1:end-1])]
     ndx = k' * (x' - 1) + 1
     lindex = convert(Int, ndx[1])

     # @assert lindex == sub2ind(siz, map(Int64, x)...) #different answer than sub2ind

     lindex
end

function bn_sort(G)
# BN_SORT Produces a topological sort of a Bayesian network.
#   Returns an array specifying the indices of variables in topological
#   order according to a graphical structure.
#
#   ORDER = BN_SORT(G) topologically sorts the variables in the specified
#   graph G and returns an array ORDER that contains the indices of the
#   variables in order. The matrix G is a square adjacency matrix.
#
# This code is based on the topological sort routine by Kevin Murphy. The
# original source code is available from:
# http://bnt.sourceforge.net/

# INPUT:
# G - a square adjacency matrix
#
# OUTPUT:
# order - an array of indices indicating order
#
# This code is based on the topological sort routine by Kevin Murphy. The
# original source code is available from:
# http://bnt.sourceforge.net/

    n = size(G)[1]
    indeg = zeros(Int, (1, n))
    zero_indeg = Int[]    # a stack of nodes with no parents

    for i = 1:n
        indeg[i] = sum(G[:, i])

        if indeg[i] == 0
            push!(zero_indeg, i)
        end
    end

    t = 1
    order = zeros(Int, (1, n))

    while !isempty(zero_indeg)
        v = pop!(zero_indeg)

        order[t] = v
        t = t + 1

        cs = find(G[v, :])

        for j = 1:length(cs)
            c = cs[j]
            indeg[c] = indeg[c] - 1

            if indeg[c] == 0
                push!(zero_indeg, c)
            end
        end
    end

    return order
end

function bn_dirichlet_prior(N)

    n = length(N)

    alpha = Array{Any}(n, 1)

    for i = 1:n
        r, q = size(N[i])
        alpha[i] = ones(r, q)
    end

    return alpha
end


text:
# ---
# title: 693. Binary Number with Alternating Bits
# id: problem693
# author: zhwang
# date: 2022-03-28
# difficulty: Easy
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/binary-number-with-alternating-bits/description/>
# hidden: true
# ---
#
# Given a positive integer, check whether it has alternating bits: namely, if
# two adjacent bits will always have different values.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 5
#     Output: true
#     Explanation: The binary representation of 5 is: 101
#
#
# **Example 2:**
#
#
#
#     Input: n = 7
#     Output: false
#     Explanation: The binary representation of 7 is: 111.
#
# **Example 3:**
#
#
#
#     Input: n = 11
#     Output: false
#     Explanation: The binary representation of 11 is: 1011.
#
# **Example 4:**
#
#
#
#     Input: n = 10
#     Output: true
#     Explanation: The binary representation of 10 is: 1010.
#
# **Example 5:**
#
#
#
#     Input: n = 3
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 231 - 1`
#
#
## @lc code=start
using LeetCode

function has_alternating_bits(n::Int)::Bool
    a = (n >> 1) ⊻ n
    return a & (a + 1) == 0
end

## @lc code=end


text:
using CUDA
using CUDAKernels
using KernelAbstractions: @kernel, @index
using Adapt: adapt_structure

using Oceananigans.Architectures: device, GPU, CPU, AbstractMultiArchitecture
using Oceananigans.Utils: work_layout

function set!(Φ::NamedTuple; kwargs...)
    for (fldname, value) in kwargs
        ϕ = getproperty(Φ, fldname)
        set!(ϕ, value)
    end
    return nothing
end

set!(u::AbstractField, v) = u .= v # fallback

# Niceties
const AbstractCPUField = AbstractField{X, Y, Z, <:CPU} where {X, Y, Z}
const AbstractReducedCPUField = AbstractReducedField{X, Y, Z, <:CPU} where {X, Y, Z}

""" Returns an AbstractReducedField on the CPU. """
function similar_cpu_field(u::AbstractReducedField)
    FieldType = typeof(u).name.wrapper
    return FieldType(location(u), CPU(), u.grid; dims=u.dims)
end

""" Set the CPU field `u` data to the function `f(x, y, z)`. """
function set!(u::AbstractCPUField, f::Function)
    f_field = FunctionField(location(u), f, u.grid)
    u .= f_field
    return nothing
end

#####
##### set! for fields on the GPU
#####

const AbstractGPUField = AbstractField{X, Y, Z, <:GPU} where {X, Y, Z}
const AbstractReducedGPUField = AbstractReducedField{X, Y, Z, <:GPU} where {X, Y, Z}

""" Returns a field on the CPU with `nothing` boundary conditions. """
function similar_cpu_field(u)
    FieldType = typeof(u).name.wrapper
    return FieldType(location(u), CPU(), adapt_structure(CPU(), u.grid), nothing)
end

""" Set the GPU field `u` to the array or function `v`. """
function set!(u::AbstractGPUField, v::Union{Array, Function})
    v_field = similar_cpu_field(u)
    set!(v_field, v)
    set!(u, v_field)
    return nothing
end

#####
##### Setting fields to fields
#####

""" Set the CPU field `u` data to the GPU field data of `v`. """
set!(u::AbstractCPUField, v::AbstractGPUField) = copyto!(parent(u), parent(v))

""" Set the GPU field `u` data to the CPU field data of `v`. """
set!(u::AbstractGPUField, v::AbstractCPUField) = copyto!(parent(u), parent(v))

set!(u::AbstractCPUField, v::AbstractCPUField) = parent(u) .= parent(v)
set!(u::AbstractGPUField, v::AbstractGPUField) = parent(u) .= parent(v)


text:
using GradDescent
using Test

include("constructors.jl")
include("quadatric.jl")
include("normal-mle.jl")
include("linear-regression.jl")
#include("optimizer-list.jl")


text:
# In this tutorial, we will learn
#  - How to implement multi-field PDEs
#  - How to build div-conforming FE spaces
#  - How to impose boundary conditions in multi-field problems
#
# ## Problem statement
#
# In this tutorial, we show how to solve a multi-field PDE in Gridap. As a model problem, we consider the Darcy equations with Dirichlet and Neumann boundary conditions. The PDE we want to solve is: find the flux vector $u$, and the pressure $p$ such that
#
# ```math
#    \left\lbrace
#    \begin{aligned}
#       \kappa^{-1} u + \nabla p = 0  \ &\text{in} \ \Omega,\\
#       \nabla \cdot u = f  \ &\text{in} \ \Omega,\\
#       u \cdot n = g \ &\text{on}\  \Gamma_{\rm D},\\
#       p = h \ &\text{on}\ \Gamma_{\rm N},\\
#    \end{aligned}
#    \right.
# ```
#
# being $n$ the outwards unit normal vector to the boundary $\partial\Omega$.  In this particular tutorial, we consider the unit square $\Omega \doteq (0,1)^2$ as the computational domain, the Neumann boundary $\Gamma_{\rm N}$ is the right and left sides of $\Omega$, and $\Gamma_{\rm D}$ is the bottom and top sides of $\Omega$. We consider $f = g \doteq 0$ and $h(x) \doteq x_1$, i.e., $h$ equal to 0 on the left side and 1 on the right side. The inverse of the permeability tensor, namely $\kappa^{-1}(x)$, is chosen equal to
#
# ```math
# \begin{pmatrix}
#   100 & 90 \\
#   90 & 100
# \end{pmatrix}
# \text{ for } \ x \in [0.4,0.6]^2, \text{ and }
# \begin{pmatrix}
#   1 & 0 \\
#   0 & 1
# \end{pmatrix}
# \ \text	{otherwise.}
# ```
#
# In order to state this problem in weak form, we introduce the following Sobolev spaces. $H(\mathrm{div};\Omega)$ is the space of vector fields in $\Omega$, whose components and divergence are in $L^2(\Omega)$. On the other hand, $H_g(\mathrm{div};\Omega)$ and $H_0(\mathrm{div};\Omega)$ are the subspaces of functions in $H(\mathrm{div};\Omega)$ such that their normal traces are equal to $g$ and $0$ respectively almost everywhere in $\Gamma_{\rm D}$. With these notations, the weak form reads: find $(u,p)\in H_g(\mathrm{div};\Omega)\times L^2(\Omega)$ such that $a((u,q),(v,q)) = b(v,q)$ for all $(v,q)\in H_0(\mathrm{div};\Omega)\times L^2(\Omega)$, where
#
# ```math
# \begin{aligned}
# a((u,q),(v,q)) &\doteq \int_{\Omega}  v \cdot \left(\kappa^{-1} u\right) \ {\rm d}\Omega - \int_{\Omega} (\nabla \cdot v)\ p \ {\rm d}\Omega + \int_{\Omega} q\ (\nabla \cdot u) \ {\rm d}\Omega,\\
# b(v,q) &\doteq \int_{\Omega} q\ f \ {\rm  d}\Omega - \int_{\Gamma_{\rm N}} (v\cdot n)\ h  \ {\rm  d}\Gamma.
# \end{aligned}
# ```
#
#
#  ## Numerical scheme
#
# In this tutorial, we use the Raviart-Thomas (RT)  space for the flux approximation [1]. On a reference square with sides aligned with the Cartesian axes, the \ac{rt} space of order $k$ is represented as $Q_{(k+1,k)} \times Q_{(k,k+1)}$, being the polynomial space defined as follows. The component  $w_\alpha$ of a vector field $w$ in $Q_{(k+1,k)} \times Q_{(k,k+1)}$ is obtained as the tensor product of univariate polynomials of order $k+1$ in direction $\alpha$ times univariate polynomials of order $k$ on the other directions. That is, $\nabla\cdot w \in Q_k$, where $Q_k$ is the multivariate polynomial space of degree at most $k$ in each of the spatial coordinates. Note that the definition of the \ac{rt} space also applies to arbitrary dimensions. The global FE space for the flux $V$ is obtained by mapping the cell-wise RT space into the physical space using the Piola transformation and enforcing continuity of normal traces across cells (see [1] for specific details).
#
#  We consider the subspace  $V_0$ of functions in $V$ with zero normal trace on $\Gamma_{\rm D}$, and the subspace $V_g$ of functions in $V$ with normal trace equal to the projection of $g$ onto the space of traces of $V$ on $\Gamma_{\rm D}$. With regard to the pressure, we consider the discontinuous space of cell-wise polynomials in $Q_k$.
#
# ## Discrete model
#
# We start the driver loading the Gridap package and constructing the geometrical model. We generate a $100\times100$ structured mesh for the domain $(0,1)^2$.

using Gridap
domain = (0,1,0,1)
partition = (100,100)
model = CartesianDiscreteModel(domain,partition)

# ## Multi-field FE spaces
#
# Next, we build the FE spaces. We consider the first order RT space for the flux and the discontinuous pressure space as described above.  This mixed FE pair satisfies the inf-sup condition and, thus, it is stable.

order = 1

V = TestFESpace(
  reffe=:RaviartThomas, order=order, valuetype=VectorValue{2,Float64},
  conformity=:HDiv, model=model, dirichlet_tags=[5,6])

Q = TestFESpace(
  reffe=:QLagrangian, order=order, valuetype=Float64,
  conformity=:L2, model=model)

# Note that the Dirichlet boundary for the flux are the bottom and top sides of the squared domain (identified with the boundary tags 5, and 6 respectively), whereas no Dirichlet data can be imposed on the pressure space. We select `conformity=:HDiv` for the flux (i.e., shape functions with $H^1(\mathrm{div};\Omega)$ regularity) and `conformity=:L2` for the pressure (i.e. discontinuous shape functions).
#
# From these objects, we construct the trial spaces. Note that we impose homogeneous boundary conditions for the flux.

uD = VectorValue(0.0,0.0)
U = TrialFESpace(V,uD)
P = TrialFESpace(Q)

# When the singe-field spaces have been designed, the multi-field test and trial spaces are expressed as arrays of single-field ones in a natural way.

Y = MultiFieldFESpace([V, Q])
X = MultiFieldFESpace([U, P])

# ## Numerical integration
#
# In this example we need to integrate in the interior of $\Omega$ and on the Neumann boundary $\Gamma_{\rm N}$. For the volume integrals, we extract the triangulation from the geometrical model and define the corresponding cell-wise quadrature of degree of exactness at least 2 as follows.

trian = Triangulation(model)
degree = 2
quad = CellQuadrature(trian,degree)

# In order to integrate the Neumann boundary condition, we only need to build an integration mesh for the right side of the domain (which is the only part of $\Gamma_{\rm N}$, where the Neumann function $h$ is different from zero). Within the model, the right side of $\Omega$ is identified with the boundary tag 8. Using this identifier, we extract the corresponding surface triangulation and create a quadrature with the desired degree of exactness.

neumanntags = [8,]
btrian = BoundaryTriangulation(model,neumanntags)
bquad = CellQuadrature(btrian,degree)

# ## Weak form
#
# We start by defining the permeability tensors inverses commented above using the `@law` macro.

const kinv1 = TensorValue(1.0,0.0,0.0,1.0)
const kinv2 = TensorValue(100.0,90.0,90.0,100.0)
@law function σ(x,u)
   if ((abs(x[1]-0.5) <= 0.1) && (abs(x[2]-0.5) <= 0.1))
      return kinv2⋅u
   else
      return kinv1⋅u
   end
end

# With this definition, we can express the integrand of the bilinear form as follows.

px = get_physical_coordinate(trian)

function a(x,y)
   v, q = y
   u, p = x
   v⋅σ(px,u) - (∇⋅v)*p + q*(∇⋅u)
end

# The arguments `x` and `y` of previous function represent a trial and a test function in the multi-field test and trial spaces `X` and `Y` respectively. In the first lines in the function definition, we unpack the single-field test and trial functions from the multi-field ones. E.g., `v` represents a test function for the flux and `q` for the pressure. These quantities can also be written as `y[1]` and `y[2]` respectively. From the single-field functions, we write the different terms of the bilinear form as we have done in previous tutorials.
#
# In a similar way, we can define the forcing term related to the Neumann boundary condition.

nb = get_normal_vector(btrian)
h = -1.0
function b_ΓN(y)
  v, q = y
  (v⋅nb)*h
end

# ## Multi-field FE problem
#
# Finally, we can assemble the FE problem and solve it. Note that we build the `AffineFEOperator` object using the multi-field trial and test spaces `Y` and `X`.

t_Ω = LinearFETerm(a,trian,quad)
t_ΓN = FESource(b_ΓN,btrian,bquad)
op = AffineFEOperator(X,Y,t_Ω,t_ΓN)
xh = solve(op)
uh, ph = xh

# Since this is a multi-field example, the `solve` function returns a multi-field solution `xh`, which can be unpacked in order to finally recover each field of the problem. The resulting single-field objects can be visualized as in previous tutorials (see next figure).

writevtk(trian,"darcyresults",cellfields=["uh"=>uh,"ph"=>ph])

# ![](../assets/darcy/darcy_results.png)
#
# ## References
#
# [1] F. Brezzi and M. Fortin. *Mixed and hybrid finite element methods*. Springer-Verlag, 1991.
#

max_stars_repo_name: andre-massing/Tutorials

max_issues_repo_name: andre-massing/Tutorials

max_forks_repo_name: andre-massing/Tutorials


text:
function trace_n_components!(data::BmmData, min_molecules_per_cell::Int)
    if !(:n_components in keys(data.tracer))
        data.tracer[:n_components] = Dict{Int, Int}[]
    end

    trace_nums = unique(max.(round.(Int, [0.0, 0.5, 1.0, 2.0, 5.0] .* min_molecules_per_cell), 1))
    n_mols_per_cell = num_of_molecules_per_cell(data)
    push!(data.tracer[:n_components], Dict(tn => sum(n_mols_per_cell .>= tn) for tn in trace_nums));
end

function trace_assignment_history!(data::BmmData, assignment_history_depth::Int; use_guids::Bool=true)
    if assignment_history_depth == 0
        return
    end

    if use_guids
        push!(data.tracer[:assignment_history], global_assignment_ids(data))
    else
        push!(data.tracer[:assignment_history], deepcopy(data.assignment))
    end
    if length(data.tracer[:assignment_history]) > assignment_history_depth
        data.tracer[:assignment_history] = data.tracer[:assignment_history][(end - assignment_history_depth + 1):end]
    end
end

function trace_component_history!(data::BmmData) # TODO: remove it
    push!(data.tracer[:component_history], deepcopy(data.components))
end

function merge_tracers(tracers::Array{Dict{Symbol, Any}, 1})::Dict{Symbol, Any}
    res = Dict{Symbol, Any}();

    tracers = filter(tr -> :n_components in keys(tr), tracers)
    if length(tracers) == 0
        return res
    end

    n_components_per_tracer = [tr[:n_components] for tr in tracers]
    n_components_res = deepcopy(n_components_per_tracer[1])

    for it in 1:length(n_components_res)
        for k in keys(n_components_res[it])
            for nc in n_components_per_tracer[2:end]
                n_components_res[it][k] += nc[it][k]
            end
        end
    end
    res[:n_components] = n_components_res

    assignments = get.(tracers, Ref(:assignment_history), Ref([]))
    if any(length.(assignments) .> 0)
        if !all(length.(assignments) .== length(assignments[1]))
            @warn "Tracers have different length of assignment_history. Can't merge it."
        else
            res[:assignment_history] = vcat.(assignments...)
        end
    end

    return res
end

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
using Flux, MLDatasets
using Flux: onehotbatch, onecold, DataLoader, Optimiser
using BSON:@save,@load

function ConvMixer(in_channels, kernel_size, patch_size, dim, depth, N_classes)
    f = Chain(
            Conv((patch_size, patch_size), in_channels=>dim, gelu; stride=patch_size),
            BatchNorm(dim),
            [
                Chain(
                    SkipConnection(Chain(Conv((kernel_size,kernel_size), dim=>dim, gelu; pad=SamePad(), groups=dim), BatchNorm(dim)), +),
                    Chain(Conv((1,1), dim=>dim, gelu), BatchNorm(dim))
                )
                for i in 1:depth
            ]...,
            AdaptiveMeanPool((1,1)),
            flatten,
            Dense(dim,N_classes)
        )
    return f
end

function get_data(batchsize; dataset = MLDatasets.CIFAR10, idxs = nothing)
    """
    idxs=nothing gives the full dataset, otherwise (for testing purposes) only the 1:idxs elements of the train set are given.
    """
    ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

    # Loading Dataset
    if idxs===nothing
        xtrain, ytrain = dataset(:train)[:]
        xtest, ytest = dataset(:test)[:]
	else
        xtrain, ytrain = dataset(:train)[1:idxs]
        xtest, ytest = dataset(:test)[1:Int(idxs/10)]
    end

    # Reshape Data to comply to Julia's (width, height, channels, batch_size) convention in case there are only 1 channel (eg MNIST)
    if ndims(xtrain)==3
        w = size(xtrain)[1]
        xtrain = reshape(xtrain, (w,w,1,:))
        xtest = reshape(xtest, (w,w,1,:))
    end

    ytrain, ytest = onehotbatch(ytrain, 0:9), onehotbatch(ytest, 0:9)

    train_loader = DataLoader((xtrain, ytrain), batchsize=batchsize, shuffle=true)
    test_loader = DataLoader((xtest, ytest), batchsize=batchsize)

    return train_loader, test_loader
end

function create_loss_function(dataloader, device)

    function loss(model)
        n = 0
        l = 0.0f0
        acc = 0.0f0

        for (x,y) in dataloader
            x,y = x |> device, y |> device
            z = model(x)
            l += Flux.logitcrossentropy(z, y, agg=sum)
            acc += sum(onecold(z).==onecold(y))
            n += size(x)[end]
        end
        l / n, acc / n
    end

    return loss

end

function train(n_epochs=100)

    #params: warning, the training can be long with these params
    train_loader, test_loader = get_data(128)
    η = 3e-4
    in_channel = 3
    patch_size = 2
    kernel_size = 7
    dim = 128
    dimPL = 2
    depth = 18
    use_cuda = true

    #logging the losses
    train_save = zeros(n_epochs, 2)
    test_save = zeros(n_epochs, 2)

    if use_cuda
        device = gpu
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    train_loss_fn = create_loss_function(train_loader, device)
    test_loss_fn = create_loss_function(test_loader, device)

    model = ConvMixer(in_channel, kernel_size, patch_size, dim, depth, 10) |> device

    ps = params(model)
    opt = Optimiser(
            WeightDecay(1f-3),
            ClipNorm(1.0),
            ADAM(η)
            )

    for epoch in 1:n_epochs
        for (x,y) in train_loader
            x,y = x|>device, y|>device
            gr = gradient(()->Flux.logitcrossentropy(model(x), y, agg=sum), ps)
            Flux.Optimise.update!(opt, ps, gr)
        end

        #logging
        train_loss, train_acc = train_loss_fn(model) |> cpu
        test_loss, test_acc = test_loss_fn(model) |> cpu
        train_save[epoch,:] = [train_loss, train_acc]
        test_save[epoch,:] = [test_loss, test_acc]

        if epoch%5==0
            @info "Epoch $epoch : Train loss = $train_loss || Validation accuracy = $test_acc."
        end

    end

    model = model |> cpu
    @save "model.bson" model
    @save "losses.bson" train_save test_save
end

if abspath(PROGRAM_FILE) == @__FILE__
    train()
end


text:
using Base.Test
using QuantumOptics

# Define parameters for spin coupled to electric field mode.
ωc = 1.2
ωa = 0.9
g = 1.0
γ = 0.5
κ = 1.1

Ntrajectories = 500
T = Float64[0.:0.1:10.;]

# Define operators
fockbasis = FockBasis(8)
spinbasis = SpinBasis(1//2)
basis = tensor(spinbasis, fockbasis)

sx = sigmax(spinbasis)
sy = sigmay(spinbasis)
sz = sigmaz(spinbasis)
sp = sigmap(spinbasis)
sm = sigmam(spinbasis)

# Hamiltonian
Ha = embed(basis, 1, 0.5*ωa*sz)
Hc = embed(basis, 2, ωc*number(fockbasis))
Hint = sm ⊗ create(fockbasis) + sp ⊗ destroy(fockbasis)
H = Ha + Hc + Hint
Hdense = full(H)

# Jump operators
Ja = embed(basis, 1, sqrt(γ)*sm)
Jc = embed(basis, 2, sqrt(κ)*destroy(fockbasis))
J = [Ja, Jc]
Jdense = map(full, J)

# Initial conditions
Ψ₀ = spinup(spinbasis) ⊗ fockstate(fockbasis, 5)
ρ₀ = Ψ₀ ⊗ dagger(Ψ₀)

# Test mcwf
tout, Ψt = timeevolution.mcwf(T, Ψ₀, Hdense, Jdense; seed=UInt(1), reltol=1e-7)
Ψ = Ψt[end]

tout, Ψt = timeevolution.mcwf(T, Ψ₀, H, J; seed=UInt(1), reltol=1e-6)
@test norm(Ψt[end]-Ψ) < 1e-5

tout, Ψt = timeevolution.mcwf(T, Ψ₀, H, J; seed=UInt(2), reltol=1e-6)
@test norm(Ψt[end]-Ψ) > 0.1

# Test mcwf_h
tout, Ψt = timeevolution.mcwf_h(T, Ψ₀, H, J; seed=UInt(1), reltol=1e-6)
@test norm(Ψt[end]-Ψ) < 1e-5

tout, Ψt = timeevolution.mcwf_h(T, Ψ₀, H, Jdense; seed=UInt(1), reltol=1e-6)
@test norm(Ψt[end]-Ψ) < 1e-5

tout, Ψt = timeevolution.mcwf_h(T, Ψ₀, Hdense, J; seed=UInt(1), reltol=1e-6)
@test norm(Ψt[end]-Ψ) < 1e-5

tout, Ψt = timeevolution.mcwf_h(T, Ψ₀, H, J; seed=UInt(2), reltol=1e-6)
@test norm(Ψt[end]-Ψ) > 0.1

# Test mcwf nh
Hnh = H - 0.5im*sum([dagger(J[i])*J[i] for i=1:length(J)])
Hnh_dense = full(Hnh)

tout, Ψt = timeevolution.mcwf_nh(T, Ψ₀, Hnh, J; seed=UInt(1), reltol=1e-6)
@test norm(Ψt[end]-Ψ) < 1e-5

tout, Ψt = timeevolution.mcwf_nh(T, Ψ₀, Hnh, Jdense; seed=UInt(1), reltol=1e-6)
@test norm(Ψt[end]-Ψ) < 1e-5

tout, Ψt = timeevolution.mcwf_nh(T, Ψ₀, Hnh_dense, J; seed=UInt(1), reltol=1e-6)
@test norm(Ψt[end]-Ψ) < 1e-5

tout, Ψt = timeevolution.mcwf_nh(T, Ψ₀, Hnh, J; seed=UInt(2), reltol=1e-6)
@test norm(Ψt[end]-Ψ) > 0.1

# Test convergence to master solution
tout_master, ρt_master = timeevolution.master(T, ρ₀, H, J)

ρ_average = DenseOperator[0 * ρ₀ for i=1:length(T)]
for i=1:Ntrajectories
    tout, Ψt = timeevolution.mcwf(T, Ψ₀, Hdense, Jdense; seed=UInt(i))
    for j=1:length(T)
        ρ_average[j] += (Ψt[j] ⊗ dagger(Ψt[j]))/Ntrajectories
    end
end
for i=1:length(T)
    err = tracedistance(ρt_master[i], ρ_average[i])
    @test err < 0.1
end

# Test single jump operator
J1 = [Ja]
J1_dense = map(full, J1)
J2 = [Ja, 0 * Jc]
J2_dense = map(full, J2)

tout_master, ρt_master = timeevolution.master(T, ρ₀, Hdense, J1_dense)

ρ_average_1 = DenseOperator[0 * ρ₀ for i=1:length(T)]
ρ_average_2 = DenseOperator[0 * ρ₀ for i=1:length(T)]
for i=1:Ntrajectories
    tout, Ψt_1 = timeevolution.mcwf(T, Ψ₀, Hdense, J1_dense; seed=UInt(i))
    tout, Ψt_2 = timeevolution.mcwf(T, Ψ₀, Hdense, J2_dense; seed=UInt(i))
    for j=1:length(T)
        ρ_average_1[j] += (Ψt_1[j] ⊗ dagger(Ψt_1[j]))/Ntrajectories
        ρ_average_2[j] += (Ψt_2[j] ⊗ dagger(Ψt_2[j]))/Ntrajectories
    end
end
for i=1:length(T)
    @test tracedistance(ρt_master[i], ρ_average_1[i]) < 0.1
    @test tracedistance(ρt_master[i], ρ_average_2[i]) < 0.1
end

# Test equivalence to schroedinger time evolution for no decay
J = DenseOperator[]
tout_schroedinger, Ψt_schroedinger = timeevolution.schroedinger(T, Ψ₀, Hdense)
tout_mcwf, Ψt_mcwf = timeevolution.mcwf(T, Ψ₀, Hdense, J)
tout_mcwf_h, Ψt_mcwf_h = timeevolution.mcwf_h(T, Ψ₀, Hdense, J)
tout_mcwf_nh, Ψt_mcwf_nh = timeevolution.mcwf_nh(T, Ψ₀, Hdense, J)

for i=1:length(T)
    @test norm(Ψt_mcwf[i] - Ψt_schroedinger[i]) < 1e-4
    @test norm(Ψt_mcwf_h[i] - Ψt_schroedinger[i]) < 1e-4
    @test norm(Ψt_mcwf_nh[i] - Ψt_schroedinger[i]) < 1e-4
end


text:
module ProgressView

export @monitor, monitor_flush, reset_tasks!, push_task!, pop_task!,
       maybe_display_tasks, display_tasks

const ESC = "\e"                # escape
const CR = "\r"                 # carriage return
const NL = "\n"                 # next line

const CSI = "$ESC["             # control sequence introducer
const CUU = "$(CSI)A"           # cursor up
const CUD = "$(CSI)B"           # cursor down
const ED = "$(CSI)J"            # erase in display
const EL = "$(CSI)K"            # erase in line

const update_every = 1.0       # update output every that many seconds

const stack = Vector{String}()  # current task stack

const lastoutput = Vector{String}() # currently displayed task stack
const lasttime = Ref(0.0)           # time at which stack was output

macro monitor(taskname::AbstractString, task)
    quote
        push_task!($taskname)
        local result = $(esc(task))
        pop_task!()
        result
    end
end

monitor_flush() = display_tasks()

function reset_tasks!()
    empty!(stack)
    display_tasks()
    return nothing
end

function push_task!(task::AbstractString)
    push!(stack, String(task))
    maybe_display_tasks()
    return nothing
end

function pop_task!()
    pop!(stack)
    maybe_display_tasks()
    return nothing
end

function maybe_display_tasks()
    currenttime = time_ns() / 1.0e+9
    if currenttime ≥ lasttime[] + update_every
        display_tasks()
        lasttime[] = currenttime
    end
    return nothing
end

function display_tasks()
    currentoutput = copy(stack)

    # Find first difference between current and last output
    tstart = 1
    while tstart ≤ length(lastoutput) &&
              tstart ≤ length(currentoutput) &&
              lastoutput[tstart] == currentoutput[tstart]
        tstart += 1
    end

    # Move cursor up to the first line that needs outputting
    tmove = length(tstart:length(lastoutput))
    print("$CUU"^tmove)
    # Output new lines
    for line in tstart:length(currentoutput)
        print("  "^(line - 1), "$line: ", currentoutput[line], "$EL$NL")
    end
    # Clear remaining lines
    print("$ED")

    # Save state
    copy!(lastoutput, currentoutput)

    return nothing
end

end


text:
# Masked Eps Greedy Policy
"""
Epsilon greedy policy that operates within a safety mask. Both actions from the greedy part and the random part are drawn from the safe actions returned by
the safety mask.
`MaskedEpsGreedyPolicy{S, A, M}(mdp::MDP{S, A}, epsilon::Float64, mask::M, rng::AbstractRNG)`
"""
struct MaskedEpsGreedyPolicy{M <: SafetyMask} <: Policy
    val::ValuePolicy # the greedy policy
    epsilon::Float64
    mask::M
    rng::AbstractRNG
end

MaskedEpsGreedyPolicy(mdp::MDP{S, A}, epsilon::Float64, mask::M, rng::AbstractRNG) where {S, A, M <: SafetyMask} = MaskedEpsGreedyPolicy{M}(ValuePolicy(mdp), epsilon, mask, rng)

function POMDPs.action(policy::MaskedEpsGreedyPolicy{M}, s) where M
    acts = safe_actions(policy.mask, s)
    if rand(policy.rng) < policy.epsilon
        return rand(policy.rng, acts)
    else
        return best_action(acts, policy.val, s)
    end
end

"""
A value policy that operates within a safety mask, it takes the action in the set of safe_actions that maximizes the given value function.
`MaskedValuePolicy{M <: SafetyMask}(val::ValuePolicy, mask::M`
"""
struct MaskedValuePolicy{M <: SafetyMask} <: Policy
    val::ValuePolicy
    mask::M
end

function POMDPs.action(policy::MaskedValuePolicy{M}, s) where M
    acts = safe_actions(policy.mask, s)
    return best_action(acts, policy.val, s)
end

function best_action(acts::Vector{A}, policy::ValuePolicy{A, M}, s) where {A,M}
    si = stateindex(policy.mdp, s)
    best_ai = 1
    best_val = policy.value_table[si, best_ai]
    for a in acts
        ai = actionindex(policy.mdp, a)
        val =  policy.value_table[si, ai]
        if val > best_val
            best_ai = ai
            best_val = val
        end
    end
    return policy.act[best_ai]
end


text:
import Base: ==, show, hash

abstract type AbstractFeature <: Function end
export AbstractFeature

"""
    𝑓 = Feature(method::Function, name=Symbol(method), keywords="", description="")

Construct a `Feature`, which is a function annotated with a `name`, `keywords` and short `description`.
Features can be called as functions while `getname(𝑓)`, `getkeywords(𝑓)` and `getdescription(𝑓)` can be used to access the annotations.
The function should have at minimum a method for `AbstractVector`.
The method on vectors will be applied column-wise to `Matrix` inputs, regardless of the function methods defined for `Matrix`.

# Examples
```julia-repl
𝑓 = Feature(sum, :sum, ["distribution"], "Sum of time-series values")
𝑓(1:10) # == sum(1:10) == 55
getdescription(𝑓) # "Sum of time-series values"
```
"""
struct Feature <: AbstractFeature
    method::Function
    name::Symbol
    keywords::Vector{String}
    description::String
    Feature(method::Function, name=Symbol(method), keywords=[""], description="") = new(method, name, keywords, description)
end
Feature(args...) = Feature{Float64}(args...)
export Feature

getmethod(𝑓::AbstractFeature)  = 𝑓.method
getname(𝑓::AbstractFeature)  = 𝑓.name
getkeywords(𝑓::AbstractFeature)  = 𝑓.keywords
getdescription(𝑓::AbstractFeature)  = 𝑓.description
export getmethod, getname, getkeywords, getdescription

(𝑓::AbstractFeature)(x::AbstractVector)  = getmethod(𝑓)(x)
(𝑓::AbstractFeature)(X::AbstractArray) = mapslices(getmethod(𝑓), X; dims=1)

# We assume that any features with the same name are the same feature
hash(𝑓::AbstractFeature, h::UInt) = hash(𝑓.name, h)
(==)(𝑓::AbstractFeature, 𝑓′::AbstractFeature) = hash(𝑓) == hash(𝑓′)

formatshort(𝑓::AbstractFeature) = ":"*string(getname(𝑓))*" "
show(𝑓::AbstractFeature) = print(formatshort(𝑓))
show(io::IO, 𝑓::AbstractFeature) = print(io, formatshort(𝑓))
show(io::IO, m::MIME"text/plain", 𝑓::AbstractFeature) = printstyled(io, formatshort(𝑓), color=:light_blue)


text:
module Trying

export @try

"""
    @try expression

Evaluate `expression` and return the result.  If it throws an error,
record the exception and backtrace in `Trying.ERRORS`.
"""
:(@try)

using Parameters

struct ErrorRecord
    expression
    exception
    backtrace
end

function Base.show(io::IO, rec::ErrorRecord)
    if get(io, :limit, false)
        printstyled(io, rec.expression; color=:blue)
        print(io, " threw ")
        err = sprint(
            summary,
            rec.exception,
            context = IOContext(io, :module => parentmodule(typeof(rec.exception))),
        )
        printstyled(io, err; color=:red)
    else
        invoke(show, Tuple{IO, Any}, io, rec)
    end
    return
end

function Base.show(io::IO, ::MIME"text/plain", rec::ErrorRecord)
    println(io, "ErrorRecord")
    printstyled(io, "Trying to evaluate expression:"; color=:blue)
    println(io)
    println(io, rec.expression)
    println(io)
    printstyled(io, "Error thrown:"; color=:blue)
    println(io)
    showerror(io, rec.exception, rec.backtrace)
end

"""
    Trying.ERRORS :: Vector{ErrorRecord}

Record of last `Trying.CONFIG.max_errors` errors.
"""
const ERRORS = ErrorRecord[]

@with_kw mutable struct TryingConfig
    max_errors::Int = 10
end

"""
    Trying.CONFIG

# Properties
- `max_errors::Int = $(TryingConfig().max_errors)`:
  Number of errors to be stored in `Trying.ERRORS`.
"""
const CONFIG = TryingConfig()

function resize_errors()
    inds = 1:length(ERRORS) - CONFIG.max_errors - 1
    if !isempty(inds)
        deleteat!(ERRORS, inds)
    end
    return
end

record_error(args...) = push!(ERRORS, ErrorRecord(args...))

function try_impl(ex)
    quote
        resize_errors()
        try
            $(esc(ex))
        catch exception
            record_error($(QuoteNode(ex)), exception, catch_backtrace())
            @info "Error thrown is available at `Trying.ERRORS[end]`."
            rethrow()
        end
    end
end

@eval macro $(:try)(ex)
    try_impl(ex)
end

end # module