text:
function sinkhorn(backend::GPUBackend, state::WassersteinLossLayerState, inputs::Vector{Blob})
  pred  = inputs[1]
  label = inputs[2]

  data_type = eltype(pred)

  pred_size = get_fea_size(pred)
  pred_num  = get_num(pred)
  label_size= get_fea_size(label)

  # init as uniform distribution
  copy!(state.u, ones(data_type, pred_size, pred_num) / pred_size);
  u = state.u
  a = pred
  b = label
  K = state.K

  if isempty(state.tmps)
    state.tmps = Blob[
      make_blob(backend, data_type, size(b)),
      make_blob(backend, data_type, size(a)),
      make_blob(backend, ones(data_type, size(a))/pred_num)
    ]
  end

  for iter = 1:state.layer.sinkhorn_iter
    # u = a ./ (K * (b./(u'*K)'))

    # tmps[1] = K' * u
    CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_T, CuBLAS.OP_N,
        label_size, pred_num, pred_size, one(data_type), K.ptr, pred_size,
        u.ptr, pred_size, zero(data_type), state.tmps[1].ptr, label_size)

    # tmps[1] = b ./ tmps[1]
    CuVec.div2!(backend, b, state.tmps[1])

    # tmps[2] = K * tmps[1]
    CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_N,
        pred_size, pred_num, label_size, one(data_type), K.ptr, pred_size,
        state.tmps[1].ptr, label_size, zero(data_type), state.tmps[2].ptr, pred_size)

    # tmps[2] = a ./ tmps[2]
    CuVec.div2!(backend, a, state.tmps[2])

    # u = tmps[2]
    copy!(u, state.tmps[2])
  end

  # compute objective function
  #-------------------------------------
  # v = b ./ (K'*u)

  # tmps[1] = K' * u
  CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_T, CuBLAS.OP_N,
      label_size, pred_num, pred_size, one(data_type), K.ptr, pred_size,
      u.ptr, pred_size, zero(data_type), state.tmps[1].ptr, label_size)

  # tmps[1] = b ./ tmps[1]
  CuVec.div2!(backend, b, state.tmps[1])

  #-------------------------------------
  # loss = sum(u .* (KM * v)) / pred_num

  # tmps[2] = KM * tmp[1]
  CuBLAS.gemm(backend.cublas_ctx, CuBLAS.OP_N, CuBLAS.OP_N,
      pred_size, pred_num, label_size, one(data_type), state.KM.ptr, pred_size,
      state.tmps[1].ptr, label_size, zero(data_type), state.tmps[2].ptr, pred_size)

  # tmps[2] = u .* tmps[2]
  CuVec.mul!(backend, state.tmps[2], u)

  # tmps[3] == ones/pred_num
  # loss = sum(tmps[2]) / pred_num
  state.loss = CuBLAS.dot(backend.cublas_ctx, data_type, length(state.tmps[2]),
      state.tmps[2].ptr, 1, state.tmps[3].ptr, 1)

  # compute gradient
  copy!(state.alpha, u)
  CuVec.log!(backend, state.alpha)
  CuBLAS.scal(backend.cublas_ctx, length(state.alpha), convert(data_type, 1.0/state.layer.lambda/pred_num),
      state.alpha.ptr, 1)
end


text:
using ElectronProgressBars
using Juno
using Logging

with_logger(ElectronProgressBars.get_logger()) do
    n = 3
    m = 50
    Juno.progress(name="outer") do id
        for i in 1:n
            @debug "outer" progress=i/n _id=id
            Juno.progress(name="outer") do id
                for j in 1:m
                    @debug "inner" progress=j/m _id=id
                    sleep(0.01)
                end
            end
        end
    end
end

ext: jl


text:
module JsonGrinder
using Mill, JSON, Printf, Flux
using HierarchicalUtils

using Mill: ArrayNode, BagNode, ProductNode, catobs, NGramMatrix

const FloatType = Float32

include("switches.jl")
include("schema/schema.jl")
include("extractors/extractors.jl")
include("html_show_tools.jl")
include("hierarchical_utils.jl")
include("util.jl")

export ExtractScalar, ExtractCategorical, ExtractArray, ExtractDict, ExtractOneHot, ExtractVector, MultipleRepresentation, ExtractString, AuxiliaryExtractor, ExtractKeyAsField
export suggestextractor, schema, extractbatch, generate_html

Base.show(io::IO, ::T) where T <: Union{JSONEntry, AbstractExtractor} = print(io, nameof(T))
Base.show(io::IO, ::MIME"text/plain", @nospecialize(n::Union{JSONEntry, AbstractExtractor})) =
    HierarchicalUtils.printtree(io, n; trav=false, htrunc=3, vtrunc=3, breakline=false)
Base.getindex(n::Union{JSONEntry, AbstractExtractor}, i::AbstractString) = HierarchicalUtils.walk(n, i)

end # module


text:
using BenchmarkTools, SparseArrays, LinearAlgebra
using ExpmV

const SUITE = BenchmarkGroup()

dimensions = 2 .^ (5:10)
density = 10. .^ (-4 : -1)

for p in density
    SUITE[p] = BenchmarkGroup(["density"])
    for d in dimensions
        SUITE[p][d] = @benchmarkable ExpmV.expmv(rt,r,rv) setup=((rt,r,rv)=(rand(), sprandn($d,$d,$p/2)+1im*sprandn($d,$d,$p/2), normalize(randn($d)+1im*randn($d))))
    end
end


text:
"Parameters for right-hand side function of implicit partitioned Runge-Kutta methods."
mutable struct ParametersIPRK{DT, TT, D, S, ET <: NamedTuple} <: Parameters{DT,TT}
    equs::ET
    tab::PartitionedTableau{TT}
    Δt::TT

    t::TT
    q::Vector{DT}
    p::Vector{DT}

    function ParametersIPRK{DT,D}(equs::ET, tab::PartitionedTableau{TT}, Δt::TT) where {D, DT, TT, ET <: NamedTuple}
        new{DT, TT, D, tab.s, ET}(equs, tab, Δt, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

@doc raw"""
Implicit partitioned Runge-Kutta integrator cache.

### Fields

* `q̃`: initial guess of q
* `p̃`: initial guess of p
* `ṽ`: initial guess of v
* `f̃`: initial guess of f
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of q
* `P`: internal stages of p
* `V`: internal stages of v
* `F`: internal stages of f
* `Y`: vector field of internal stages of q
* `Z`: vector field of internal stages of p
"""
struct IntegratorCacheIPRK{ST,D,S} <: PODEIntegratorCache{ST,D}
    q̃::Vector{ST}
    p̃::Vector{ST}
    ṽ::Vector{ST}
    f̃::Vector{ST}
    s̃::Vector{ST}

    Q::Vector{Vector{ST}}
    P::Vector{Vector{ST}}
    V::Vector{Vector{ST}}
    F::Vector{Vector{ST}}
    Y::Vector{Vector{ST}}
    Z::Vector{Vector{ST}}

    function IntegratorCacheIPRK{ST,D,S}() where {ST,D,S}
        # create temporary vectors
        q̃ = zeros(ST,D)
        p̃ = zeros(ST,D)
        ṽ = zeros(ST,D)
        f̃ = zeros(ST,D)
        s̃ = zeros(ST,D)

        # create internal stage vectors
        Q = create_internal_stage_vector(ST, D, S)
        P = create_internal_stage_vector(ST, D, S)
        V = create_internal_stage_vector(ST, D, S)
        F = create_internal_stage_vector(ST, D, S)
        Y = create_internal_stage_vector(ST, D, S)
        Z = create_internal_stage_vector(ST, D, S)

        new(q̃, p̃, ṽ, f̃, s̃, Q, P, V, F, Y, Z)
    end
end

function IntegratorCache(params::ParametersIPRK{DT,TT,D,S}; kwargs...) where {DT,TT,D,S}
    IntegratorCacheIPRK{DT,D,S}(; kwargs...)
end

function IntegratorCache{ST}(params::ParametersIPRK{DT,TT,D,S}; kwargs...) where {ST,DT,TT,D,S}
    IntegratorCacheIPRK{ST,D,S}(; kwargs...)
end

@inline CacheType(ST, params::ParametersIPRK{DT,TT,D,S}) where {DT,TT,D,S} = IntegratorCacheIPRK{ST,D,S}

@doc raw"""
Implicit partitioned Runge-Kutta integrator solving the system
```math
\begin{aligned}
V_{n,i} &= v (Q_{n,i}, P_{n,i}) , &
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\
F_{n,i} &= f (Q_{n,i}, P_{n,i}) , &
P_{n,i} &= p_{n} + h  \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} , &
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} ,
\end{aligned}
```
Usually we are interested in Hamiltonian systems, where
```math
\begin{aligned}
V_{n,i} &= \dfrac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , &
F_{n,i} &= - \dfrac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) ,
\end{aligned}
```
and tableaus satisfying the symplecticity conditions
```math
\begin{aligned}
b_{i} \bar{a}_{ij} + \bar{b}_{j} a_{ji} &= b_{i} \bar{b}_{j} , &
\bar{b}_i &= b_i .
\end{aligned}
```
"""
struct IntegratorIPRK{DT, TT, D, S, PT <: ParametersIPRK{DT,TT},
                                    ST <: NonlinearSolver{DT},
                                    IT <: InitialGuessPODE{TT}} <: AbstractIntegratorPRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorIPRK(params::ParametersIPRK{DT,TT,D,S}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,ST,IT}
        new{DT, TT, D, S, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorIPRK{DT,D}(equations::NamedTuple, tableau::PartitionedTableau{TT}, Δt::TT) where {DT,TT,D}
        # get number of stages
        S = tableau.s

        # create params
        params = ParametersIPRK{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, 2*D*S, params, caches)

        # create initial guess
        iguess = InitialGuessPODE(get_config(:ig_extrapolation), equations[:v], equations[:f], Δt)

        # create integrator
        IntegratorIPRK(params, solver, iguess, caches)
    end

    function IntegratorIPRK{DT,D}(v::Function, f::Function, tableau::PartitionedTableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
        IntegratorIPRK{DT,D}(NamedTuple{(:v,:f)}((v,f)), tableau, Δt; kwargs...)
    end

    function IntegratorIPRK{DT,D}(v::Function, f::Function, h::Function, tableau::PartitionedTableau{TT}, Δt::TT; kwargs...) where {DT,TT,D}
        IntegratorIPRK{DT,D}(NamedTuple{(:v,:f,:h)}((v,f,h)), tableau, Δt; kwargs...)
    end

    function IntegratorIPRK(equation::Union{PODE{DT}, HODE{DT}}, tableau::PartitionedTableau{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorIPRK{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

@inline Base.ndims(::IntegratorIPRK{DT,TT,D,S}) where {DT,TT,D,S} = D

function initialize!(int::IntegratorIPRK, sol::AtomicSolutionPODE)
    sol.t̄ = sol.t - timestep(int)

    equations(int)[:v](sol.t, sol.q, sol.p, sol.v)
    equations(int)[:f](sol.t, sol.q, sol.p, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̄, sol.q̄, sol.p̄, sol.v̄, sol.f̄)
end

function update_params!(int::IntegratorIPRK, sol::AtomicSolutionPODE)
    # set time for nonlinear solver and copy previous solution
    int.params.t  = sol.t
    int.params.q .= sol.q
    int.params.p .= sol.p
end

function initial_guess!(int::IntegratorIPRK{DT}, sol::AtomicSolutionPODE{DT},
                        cache::IntegratorCacheIPRK{DT}=int.caches[DT]) where {DT}

    for i in eachstage(int)
        evaluate!(int.iguess, sol.q̄, sol.p̄, sol.v̄, sol.f̄,
                              sol.q, sol.p, sol.v, sol.f,
                              cache.Q[i], cache.P[i], cache.V[i], cache.F[i],
                              tableau(int).q.c[i], tableau(int).p.c[i])
    end
    for i in eachstage(int)
        for k in eachdim(int)
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+1] = 0
            int.solver.x[2*(ndims(int)*(i-1)+k-1)+2] = 0
            for j in eachstage(int)
                int.solver.x[2*(ndims(int)*(i-1)+k-1)+1] += tableau(int).q.a[i,j] * cache.V[j][k]
                int.solver.x[2*(ndims(int)*(i-1)+k-1)+2] += tableau(int).p.a[i,j] * cache.F[j][k]
            end
        end
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
                                        P::Vector{Vector{ST}}, F::Vector{Vector{ST}}, Z::Vector{Vector{ST}},
                                        params::ParametersIPRK{DT,TT,D,S}) where {ST,DT,TT,D,S}
    local tqᵢ::TT
    local tpᵢ::TT

    for i in 1:S
        for k in 1:D
            # copy y to Y and Z
            Y[i][k] = x[2*(D*(i-1)+k-1)+1]
            Z[i][k] = x[2*(D*(i-1)+k-1)+2]

            # compute Q and P
            Q[i][k] = params.q[k] + params.Δt * Y[i][k]
            P[i][k] = params.p[k] + params.Δt * Z[i][k]
        end

        # compute time of internal stage
        tqᵢ = params.t + params.Δt * params.tab.q.c[i]
        tpᵢ = params.t + params.Δt * params.tab.p.c[i]

        # compute v(Q,P) and f(Q,P)
        params.equs[:v](tqᵢ, Q[i], P[i], V[i])
        params.equs[:f](tpᵢ, Q[i], P[i], F[i])
    end
end

# Compute stages of implicit partitioned Runge-Kutta methods.
function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersIPRK{DT,TT,D,S},
                          caches::CacheDict) where {ST,DT,TT,D,S}

    # get cache for internal stages
    cache = caches[ST]

    # compute stages from nonlinear solver solution x
    compute_stages!(x, cache.Q, cache.V, cache.Y, cache.P, cache.F, cache.Z, params)

    # compute b = - [(Y-AV), (Z-AF)]
    for i in 1:S
        for k in 1:D
            b[2*(D*(i-1)+k-1)+1] = - cache.Y[i][k]
            b[2*(D*(i-1)+k-1)+2] = - cache.Z[i][k]
            for j in 1:S
                b[2*(D*(i-1)+k-1)+1] += params.tab.q.a[i,j] * cache.V[j][k]
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * cache.F[j][k]
            end
        end
    end
end

function integrate_step!(int::IntegratorIPRK{DT,TT}, sol::AtomicSolutionPODE{DT,TT},
                         cache::IntegratorCacheIPRK{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from atomic solution
    update_params!(int, sol)

    # compute initial guess
    initial_guess!(int, sol, cache)

    # reset atomic solution
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache.Q, cache.V, cache.Y,
                                  cache.P, cache.F, cache.Z, int.params)

    # compute final update
    update_solution!(sol.q, sol.q̃, cache.V, tableau(int).q.b, tableau(int).q.b̂, timestep(int))
    update_solution!(sol.p, sol.p̃, cache.F, tableau(int).p.b, tableau(int).p.b̂, timestep(int))

    # copy solution to initial guess
    update_vector_fields!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f)
end


text:
@testset "436.find-right-interval.jl" begin
    @test find_right_interval([(1, 2)]) == [-1]
    @test find_right_interval([(3, 4), (2, 3), (1, 2)]) == [-1, 1, 2]
    @test find_right_interval([(1, 4), (2, 3), (3, 4)]) == [-1, 3, -1]
end


text:
"""
    callSParseAim(hh, leads, lags)

Julia wrapper function that calls the c function callSparseAim.
"""
function callSparseAim( hh, leads, lags )

    # allocate space for the matrices and initialize inputs
    neq        = size(hh, 1)
    nstate     = 0
    hrows      = neq
    hcols      = size(hh, 2)
    qmax       = hrows*leads*(hrows*(lags+leads+1))
    retCodePtr = 0
    cofb       = zeros(neq, neq * lags)
    qmatrix    = zeros(neq*leads, hcols)

    # use the library libSPARSEAMA to call c function ...
    # libAndersonMoore is a shared library that combines sparseAMA
    # and LAPACK. LAPACK must be compiled with -fPIC.
    ccall(sym, Compat.Nothing,
          (  Ptr{Float64}, Int32, Int32, Int32, Int32,
          Int32, Int32, Int32,
          Ptr{Float64}, Ptr{Float64}, Ptr{Float64}  ),
          hh, hrows, hcols, neq, leads, lags, nstate,
          qmax, Ref{retCodePtr}[], cofb, qmatrix)

    return hh, cofb, qmatrix, retCodePtr
end


text:
#=
Given a list of numbers and a number k, return whether any two numbers from the list add up to k.

For example, given [10, 15, 3, 7] and k of 17, return true since 10 + 7 is 17.
=#

using Test
include("Solutions/problem1_addupto.jl")
@test addupto([10, 15, 3, 7], 17)

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
using DynamicProgrammingFibonacci.Dp

@test fibonacci(0) == 0
@test fibonacci(1) == 1
@test fibonacci(2) == 1
@test fibonacci(4) == 3
@test fibonacci(100) == 354224848179261915075
@test fibonacci(140) == 81055900096023504197206408605
@test fibonacci(300) == 222232244629420445529739893461909967206666939096499764990979600


text:
using Plots

function plot_rs(res, shp; nlayout = [4, 4])
    np = prod(nlayout)
    lons, lats = get_coords(shp)
    ps = []
    for i in 1:length(res)
        p = plot(res[i])
        scatter!(p, [lons[i]], [lats[i]], label = "outlet", title = shp.id[i])
        outfile = @sprintf("Figures/%02d.png", i)
        println(outfile)
        savefig(outfile)
        push!(ps, p)
    end
    ps
end


text:
################################################################################
## Copyright 2016. Juri Hinz and Jeremy Yee.
## Convex switching system model
################################################################################

## Define CSS type

immutable CSS

    ## Model
    grid::Array{Float64, 2} # (gnum, snum) grid
    weight::Array{Float64, 1} # (dnum) weights of sampled disturbances
    disturb::Array{Float64, 3} # (snum, snum, dnum) sampled distrubances
    control::Array # transition probabilities for positions
    Reward::Function # reward
    Scrap::Function # reward at last time point

    ## Parameters
    gnum::Int64 # grid size
    snum::Int64 # state size
    dnum::Int64 # disturbance sampling size
    pnum::Int64 # number of positions
    anum::Int64 # number of actions
    tnum::Int64 # number of time points
    fullcontrol::Bool # is the positions fully controllable

    ## Assigning the parameters
    function CSS(grid::Array{Float64, 2}, weight::Array{Float64, 1},
                 disturb::Array{Float64, 3}, control::Array,
                 Reward::Function, Scrap::Function, tnum::Int64)

        ## Extract parameters
        gnum = size(grid)[1]
        snum = size(grid)[2]
        dnum = length(weight)
        pnum = size(control)[1]
        anum = size(control)[2]

        if ndims(control) == 3
            fullcontrol = false
        elseif ndims(control) == 2
            fullcontrol = true
        else
            error("control not in correct format")
        end

        return new(grid, weight, disturb, control, Reward, Scrap, gnum, snum,
                   dnum, pnum, anum, tnum, fullcontrol)

    end

end


text:
merge_ext!(S::SeisData, Ω::Int64, rest::Array{Int64, 1}) = nothing

function get_δt(t::Int64, Δ::Int64)
  δts = rem(t, Δ)
  if δts > div(Δ,2)
    δts -= Δ
  end
  return δts
end

function check_alignment(Ti::Array{Int64,1}, Tj::Array{Int64,1}, Xi::Array{T,1}, Xj::Array{T,1}, Δ::Int64) where {T<:AbstractFloat}
  z = zero(Int64)

  # Rare case, but possible: they're the same data
  if first(Ti) == first(Tj) && last(Ti) == last(Tj) && Xi == Xj
    return Ti, Xi, false, z
  else

    # Find a time correction δt that gets applied to Ti, ts[i], te[i], etc.
    # Find overlapping region, if one exists
    L = length(Ti)
    for n = 1:4
      if n > L-1
        break
      end

      xi_f = view(Xi, n+1:L)
      xj_f = view(Xj, 1:L-n)
      if isapprox(xi_f, xj_f)
        return Tj[1:L-n], xj_f, false, 1*n
      end

      xi_b = view(Xi, 1:L-n)
      xj_b = view(Xj, n+1:L)
      if isapprox(xi_b, xj_b)
        return Ti[1:L-n], xi_f, false, -1*n
      end
    end
  end
  return vcat(Ti,Tj), vcat(Xi,Xj), true, z
end

function xtmerge!(t::Array{Int64,1}, x::Array{T,1}, d::Int64) where {T<:AbstractFloat}
  # Sanity check
  (length(t) == length(x)) || error(string("Badly set times (Nt=", length(t), ",Nx=", length(x), "); can't merge!"))

  # Sort
  i = sortperm(t)
  sort!(t)
  x[:] = x[i]

  # Check for duplicates
  J0 = findall((diff(t).==0).*(diff(x).==0))
  while !isempty(J0)
    deleteat!(x, J0)
    deleteat!(t, J0)
    J0 = findall(diff(t) .== 0)
  end

  J0 = findall(diff(t) .< d)
  while !isempty(J0)
    J1 = J0.+1
    K = [isnan.(x[J0]) isnan.(x[J1])]

    # Average nearly-overlapping x that are either both NaN or neither Nan
    ii = findall(K[:,1].==K[:,2])
    i0 = J0[ii]
    i1 = J1[ii]
    t[i0] = div.(t[i0].+t[i1], 2)
    x[i0] = 0.5.*(x[i0].+x[i1])

    # Delete nearly-overlapping x with only one NaN (and delete all x ∈ i1)
    i3 = findall(K[:,1].*(K[:,2].==false))
    i4 = findall((K[:,1].==false).*K[:,2])
    II = sort([J0[i4]; J1[i3]; i1])
    deleteat!(t, II)
    deleteat!(x, II)

    J0 = findall(diff(t) .< d)
  end
  return nothing
end

# get hash of each non-empty loc, resp, units, fs ( fs should never be empty )
function get_subgroups( LOC   ::Array{InstrumentPosition,1},
                        FS    ::Array{Float64,1},
                        RESP  ::Array{InstrumentResponse,1},
                        UNITS ::Array{String,1},
                        group ::Array{Int64,1} )
  zh = zero(UInt64)
  N_grp = length(group)
  N_grp == 1 && return([group])

  H = Array{UInt,2}(undef, N_grp, 4)
  for (n,g) in enumerate(group)
    H[n,1] = hash(getindex(FS, g))
    H[n,2] = isempty(getindex(LOC, g)) ? zh : hash(getindex(LOC, g))
    H[n,3] = isempty(getindex(RESP, g)) ? zh : hash(getindex(RESP, g))
    H[n,4] = isempty(getindex(UNITS, g)) ? zh : hash(getindex(UNITS, g))
  end

  # If an entire column is unset, we don't care
  H = H[:, findall([sum(H[:,i])>0 for i in 1:4])]
  (Nh,Nc) = size(H)

  # Find unique rows of H; sort
  H_filled = sum(H .!= zh, dims=2)[:]
  H_inds = sortperm(H_filled, rev=true)
  N_subgrp = length(H_inds)
  N_subgrp == 1 && return([group])

  subgrp_inds = Array{Array{Int64,1},1}(undef, N_subgrp)

  H = H[H_inds, :]
  group = group[H_inds]
  H_sub = deepcopy(H)

  for i = 1:N_subgrp
    subgrp_inds[i] = Array{Int64,1}(undef, 0)
    subgrp_hash = H_sub[i, :]
    for j = N_grp:-1:1
      m = prod([H[j,k] in (zh, subgrp_hash[k]) for k=1:Nc])
      if m
        push!(subgrp_inds[i], group[j])
        deleteat!(group, j)
        H = H[setdiff(1:end, j), :]
      end
      N_grp = length(group)
    end
    if N_grp == 0
      subgrp_inds = subgrp_inds[1:i]
      N_subgrp = length(subgrp_inds)
      break
    end
  end
  deleteat!(subgrp_inds, [isempty(subgrp_inds[j]) for j=1:length(subgrp_inds)])
  return subgrp_inds
end

function get_next_pair(W::Array{Int64,2})
  L = size(W,1)
  i = 1

  # dest loop
  while i < L
    si = getindex(W, i)
    ei = getindex(W, i + L)
    j = i + 1

    # src loop
    while j ≤ L
      sj = getindex(W, j)
      ej = getindex(W, j+L)
      if min(si ≤ ej, ei ≥ sj) == true
        # src    # dest
        return vcat(W[j,:], j), vcat(W[i,:], i)
      end
      j = j + 1
    end
    i = i + 1
  end
  return zeros(Int64, 7), zeros(Int64, 7)
end

# merges into the channel with the most recent data
function merge_non_ts!(S::GphysData, subgrp::Array{Int64,1})
  te = [maximum(t[:,2]) for t in S.t[subgrp]]
  Ω = subgrp[argmax(te)]
  T = vcat(S.t[subgrp]...)[:,2]
  X = vcat(S.x[subgrp]...)
  Z = unique(collect(zip(T,X)))
  T = first.(Z)
  X = last.(Z)
  ii = sortperm(T)
  S.t[Ω] = hcat(collect(1:1:length(T)), T[ii])
  S.x[Ω] = X[ii]
  return Ω
end


text:
# v0.6

using Combinatorics

function permsort(x::Array)
    for perm in permutations(x)
        if issorted(perm)
            return perm
        end
    end
end

x = randn(10)
@show x permsort(x)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module MLCourse

import Pkg

project_relative_path(xs...) = normpath(joinpath(dirname(dirname(pathof(MLCourse))), xs...))

include("notebooks.jl")
include_dependency(joinpath("..", "Project.toml"))
const _VERSION = VersionNumber(Pkg.TOML.parsefile(project_relative_path("Project.toml"))["version"])

using Requires, Git, Pluto, Distributed, PackageCompiler
using Markdown, Base64

function __init__()
@require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" begin
using Zygote
"""
    gradient_descent(f, x₀, η, T; callback = x -> nothing)

Run gradient descent on `f` with initial condition `x₀` and learning rate `η`
for `T` steps. The function `callback` is executed after every update of the
of the parameters and takes the current value of the parameters as input.
"""
function gradient_descent(f, x₀, η, T; callback = x -> nothing)
    x = copy(x₀) # use copy to not overwrite the input
    for t in 1:T
        x .-= η * gradient(f, x)[1] # update parameters in direction of -∇f
        callback(x) # the callback will be used to save intermediate values
    end
    x
end
end # Zygote
@require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
import Plots: plot!, scatter, scatter!, annotate!
export plot_residuals!, biplot
function plot_residuals!(x, y, f; kwargs...)
    for (xi, yi) in zip(x, y)
        plot!(fill(xi, 2), [yi, f(xi)]; color = :red, linewidth = 2,
              label = nothing, kwargs...)

    end
end
function plot_residuals!(x1, x2, y, f; kwargs...)
    for (x1i, x2i, yi) in zip(x1, x2, y)
        plot!(fill(x1i, 2), fill(x2i, 2), [yi, f(x1i, x2i)];
              color = :red, linewidth = 2, label = nothing, kwargs...)
    end
end
"""
    biplot(pca; pc = (1, 2), score_style = :text, loading = :all)

Show a biplot for a fitted machine `pca` and components `pc` (by default the
first two components are shown. If `score_style` is `:text` the index of each
data point (row in the input data of a `PCA` machine) is plotted at the
coordinates that indicate the scores of that data point. If `score_style` is an
integer the scores are plotted as points of size `score_style`. If `loadings` is
`:all`, all loadings are shown. If `loadings` is an integer ``n``, the ``n``
largest loadings along the PCs defined by `pc` are shown. If `loadings` is a
vector or integers, the loadings for the input dimensions given in this list of
integers is shown.

### Example
```
julia> using MLJ, MLJMultivariateStatsInterface, DataFrames, Plots
julia> mach = fit!(machine(PCA(), DataFrame(randn(10, 4) * randn(4, 6), :auto)))
julia> biplot(mach)
julia> biplot(mach, loadings = 3)
julia> biplot(mach, pc = (1, 3), score_style = 4, loadings = [3, 4])
```
"""
function biplot(m;
                pc = (1, 2),
                score_style = :text,
                loadings = :all)
    scores = MLJ.transform(m, m.data[1])
    p = scatter(getproperty(scores, Symbol(:x, pc[1])),
                getproperty(scores, Symbol(:x, pc[2])),
                label = nothing, aspect_ratio = 1, size = (600, 600),
                xlabel = "PC$(pc[1])", ylabel = "PC$(pc[2])",
                framestyle = :axis, markeralpha = score_style == :text ? 0 : 1,
                c = :gray,
                txt = score_style == :text ? Plots.text.(1:nrows(scores), 8, :gray) :
                                             nothing,
                markersize = score_style == :text ? 0 :
                             isa(score_style, Int) ? score_style : 3)
    plot!(p[1], inset = (1, Plots.bbox(0, 0, 1, 1)),
          right_margin = 10Plots.mm, top_margin = 10Plots.mm)
    p2 = p[1].plt.subplots[end]
    plot!(p2, aspect_ratio = 1, mirror = true, legend = false)
    params = fitted_params(m)
    ϕ = if hasproperty(params, :pca)
        params.pca.projection
    else
        params.projection
    end
    n = colnames(m.data[1])
    idxs = if loadings == :all
        1:length(n)
    elseif isa(loadings, Int)
        l = reshape(sum(abs2.(ϕ[:, [pc...]]), dims = 2), :)
        sortperm(l, rev = true)[1:loadings]
    else
        loadings
    end
    for i in idxs
        plot!(p2, [0, .9*ϕ[i, pc[1]]], [0, .9*ϕ[i, pc[2]]],
              c = :red, arrow = true)
        annotate!(p2, [(ϕ[i, pc[1]], ϕ[i, pc[2]],
                        (n[i], :red, :center, 8))])
    end
    scatter!(p2, 1.1*ϕ[:, pc[1]], 1.1*ϕ[:, pc[2]],
             markersize = 0, markeralpha = 0) # dummy for sensible xlims and xlims
    plot!(p2,
          background_color_inside = Plots.RGBA{Float64}(0, 0, 0, 0),
          tickfontcolor = Plots.RGB(1, 0, 0))
    p
end
end # Plots
@require DataFrames="a93c6f00-e57d-5684-b7b6-d8193f3e46c0" begin
using DataFrames
polysymbol(x, d) = Symbol(d == 0 ? "" : d == 1 ? "$x" : "$x^$d")

colnames(df::AbstractDataFrame) = names(df)
colnames(d::NamedTuple) = keys(d)
colname(names, predictor::Int) = names[predictor]
colname(names, predictor::Symbol) = predictor ∈ names || string(predictor) ∈ names ? Symbol(predictor) : error("Predictor $predictor not found in $names.")
function poly(data, degree, predictors::NTuple{1} = (1,))
    cn = colnames(data)
    col = colname(cn, predictors[1])
    res = DataFrame([getproperty(data, col) .^ k for k in 1:degree],
                    [polysymbol(col, k) for k in 1:degree])
    if hasproperty(data, :y)
        res.y = data.y
    end
    res
end
function poly(data, degree, predictors::NTuple{2})
    cn = colnames(data)
    col1 = colname(cn, predictors[1])
    col2 = colname(cn, predictors[2])
    res = DataFrame([getproperty(data, col1) .^ d1 .* getproperty(data, col2) .^ d2
                     for d1 in 0:degree, d2 in 0:degree if 0 < d1 + d2 ≤ degree],
                    [Symbol(polysymbol(col1, d1), polysymbol(col2, d2))
                     for d1 in 0:degree, d2 in 0:degree if 0 < d1 + d2 ≤ degree])
    if hasproperty(data, :y)
        res.y = data.y
    end
    res
end
poly(_, _, _) = error("Polynomials in more than 2 predictors are not implemented.")
end # DataFrames
@require MLJ="add582a8-e3ab-11e8-2d5e-e98b27df1bc7" begin
using MLJ
Base.@kwdef mutable struct Polynomial{T} <: Static
    degree::Int = 3
    predictors::T = (1,)
end
MLJ.transform(p::Polynomial, _, X) = poly(X, p.degree, p.predictors)
export Polynomial
PolynomialRegressor(; kwargs...) = error("`PolynomialRegressor(degree, regressor) is deprecated. Use e.g. `@pipeline(Polynomial(degree), LinearRegressor())` instead.")
end # MLJ
end # __init__

function fitted_linear_func(mach)
    θ̂ = fitted_params(mach)
    θ̂₀ = θ̂.intercept
    θ̂₁ = θ̂.coefs[1][2]
    x -> θ̂₀ + θ̂₁ * x
end;

_grid(x, y) = (repeat(x, length(y)), repeat(y, inner = length(x)))
function grid(x, y; names = (:X, :Y), output_format = NamedTuple)
    t = NamedTuple{names}(_grid(x, y))
    output_format(t)
end

function embed_figure(name)
    "![](data:img/png; base64,
         $(open(base64encode,
                project_relative_path("notebooks", "figures", name))))"
end

heaviside(x::T) where T = ifelse(x > 0, one(T), zero(T))

function start()
    nb = joinpath(project_relative_path(), "index.jl")
    exeflags = ["--project=$(project_relative_path())"]
    sysimg = project_relative_path("precompile", "mlcourse.so")
    isfile(sysimg) && push!(exeflags, "--sysimage=$sysimg")
    pid = Distributed.addprocs(1, exeflags = exeflags) |> first
    expr = :(using Pluto; Pluto.run(notebook = $nb,
                                    workspace_use_distributed = false,
                                    dismiss_update_notification = true))
    @async Distributed.remotecall_eval(Main, [pid], expr)
    println("Starting a Pluto notebook in your browswer.
Please use `MLCourse.stop()` to interrupt the Pluto notebook server.")
end

stop() = Distributed.interrupt()

function update(; kwargs...)
    @info "Performing an automatic update while keeping local changes.
    If this fails, please run manually `git pull` in the directory
    `$(project_relative_path())`."
    current_dir = pwd()
    cd(project_relative_path())
    if !isempty(readlines(`$(git()) diff --stat`))
        run(`$(git()) add -u`)
    end
    if !isempty(readlines(`$(git()) diff --cached`))
        run(`$(git()) -c user.name="student" -c user.email="student@mlcourse" commit -m "automatic commit of local changes"`)
    end
    run(`$(git()) pull -s recursive -X patience -X ours -X ignore-all-space --no-edit`)
    if !isempty(read(`$(git()) diff HEAD^ --name-only Manifest.toml`, String))
        Pkg.instantiate()
        MLCourse.create_sysimage(; kwargs...)
    end
    cd(current_dir)
end

function create_sysimage(; sysimage_path = project_relative_path("precompile", "mlcourse.so"))
    PackageCompiler.create_sysimage([:Pluto, :MLJ, :MLJLinearModels,
                                     :MLCourse, :StatsPlots, :DataFrames,
                                     :MLJFlux, :Distributions];
                                    sysimage_path,
                                    precompile_execution_file=project_relative_path("precompile", "warmup.jl"))
end

end # module

max_stars_repo_name: baltisberger/MLCourse

max_issues_repo_name: baltisberger/MLCourse

max_forks_repo_name: baltisberger/MLCourse


text:
using Test
using Dash
using Base64
@testset "Send String" begin
    data = "test string"
    res = dcc_send_string(data, "test_file.csv"; type = "text/csv")
    @test res[:content] == data
    @test res[:filename] == "test_file.csv"
    @test res[:type] == "text/csv"
    @test !res[:base64]

    data2 = "test 2 string"
    res = dcc_send_string(data2, "test_file.csv"; type = "text/csv") do io, data
        write(io, data)
    end
    @test res[:content] == data2
    @test res[:filename] == "test_file.csv"
    @test res[:type] == "text/csv"
    @test !res[:base64]
end
@testset "Send Bytes" begin
    data = "test string"
    res = dcc_send_bytes(Vector{UInt8}(data), "test_file.csv"; type = "text/csv")
    @test res[:content] == Base64.base64encode(data)
    @test res[:filename] == "test_file.csv"
    @test res[:type] == "text/csv"
    @test res[:base64]

    data2 = "test 2 string"
    res = dcc_send_bytes(write, data2, "test_file.csv"; type = "text/csv")
    @test res[:content] == Base64.base64encode(data2)
    @test res[:filename] == "test_file.csv"
    @test res[:type] == "text/csv"
    @test res[:base64]
end
@testset "Send File" begin

    file = "assets/test.png"
    res = dcc_send_file(file, nothing; type = "text/csv")
    @test res[:content] == Base64.base64encode(read(file))
    @test res[:filename] == "test.png"
    @test res[:type] == "text/csv"
    @test res[:base64]

    file = "assets/test.png"
    res = dcc_send_file(file, "ttt.jpeg"; type = "text/csv")
    @test res[:content] == Base64.base64encode(read(file))
    @test res[:filename] == "ttt.jpeg"
    @test res[:type] == "text/csv"
    @test res[:base64]
end


text:
using NewtonsMethod2
using Test

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


text:
"""
    assigndist_uniform(lower, upper)

Assign parameters to a uniform distribution with `lower` and `upper` uncertainty
bounds.
"""
function assigndist_uniform(lower, upper)
    Uniform(lower, upper)
end

export assigndist_uniform


text:
using Dates
const SRC_PATH = "../../../src"

include(joinpath(SRC_PATH, "MLPClassifier/MLPClassifier.jl"))
include(joinpath(SRC_PATH, "datasets/dataset_reader.jl"))
include(joinpath(SRC_PATH, "datasets/Drug-Classification/dataset1.jl"))
include(joinpath(SRC_PATH, "FileInterface/fileInterface.jl"))
include("../MLP_test_functions.jl")
include(joinpath(SRC_PATH, "MutualInformation/information-plane.jl"))
#include(joinpath(SRC_PATH, "MutualInformation/information-plane.jl"))

function show_results(dt, net)
    println("** Test set")
    println("Accuracy: ", accuracy(zip(dt.input_test, dt.output_test), net))
    println("Error rate: ", error_rate(zip(dt.input_test, dt.output_test), net))
    println("** Training set")
    println("Accuracy: ", accuracy(zip(dt.input_train, dt.output_train), net))
    println("Error rate: ", error_rate(zip(dt.input_train, dt.output_train), net))
end

t0 = time()
#topology = [15,5,10]
topology = [5]
#topology = [10,5]
#topology = [5,5]
Random.seed!(1234)
net1 = MLPClassifier(dataset1_nfeatures, dataset1_labels, topology)
#Random.seed!(1234)
#net2 = MLPClassifier(dataset1_nfeatures, dataset1_labels, topology)
net2 = deepcopy(net1)
#net = MLPClassifier(dataset1_nfeatures, dataset1_labels, [5])
path = joinpath(SRC_PATH, "datasets/Drug-Classification")
dt = get_dataset1_setup(path, normalise=true, separate=true)
dt_b = get_dataset1_setup(path, normalise=true, separate=true, balance=true)

#epochs = 10000
epochs = 5000

println("Topology: ", topology)
println("Epochs: ", epochs)

#Random.seed!(1234)
println("-> Not balanced")
println()
println("--> Results before training")
show_results(dt, net1)
fit!(net1, dt.input_train, dt.output_train; max_epochs = epochs)
println()
println("--> Results after training")
show_results(dt, net1)
println()
t1 = time()
#Random.seed!(1234)
println("-> Balanced")
println()
println("--> Results before training")
show_results(dt_b, net2)
fit!(net2, dt_b.input_train, dt_b.output_train; max_epochs = epochs)
println()
println("--> Results after training")
show_results(dt_b, net2)

println()

t2 = time()
println("Training time (not balanced):", t1 - t0," seconds.")
println("Training time (balanced):", t2 - t1," seconds.")

current_date = Dates.now()
time_str = Dates.format(current_date, "HH:MM:SS")
date_str = Dates.format(current_date, "dd-mm-yyyy")

#filename = "MLP-nets-dt1.jld2"
filename = string("MLP-nets-dt1-", date_str, "-", time_str, ".jld2")
addToFile(filename, "Net-notbalanced", net1)
addToFile(filename, "Net-balanced", net2)

#=
data = zip(dt.input_train, dt.output_train)
data_test = zip(dt.input_test, dt.output_test)
epochs = 10000
number_nets = 20
step_len = 1
nbins = 5
plot_IP_average(dt.n_features, dt.labels, topology, number_nets, data, epochs, step_len, "IP-dt1-MLP.png")
plot_MLP_error_rate(dt.n_features, dt.labels, topology, number_nets, data, data_test, epochs, step_len, "MLP-error_rate.png")
points_IP = calculate_IP_average_points(dt.n_features, dt.labels, topology, number_nets, data, epochs, step_len, nbins)
generate_plot_IP_average(points_IP, number_nets, epochs, step_len, "IP-dt1-MLP-$(nbins)bins.png")
=#
# Criar uma função para salvar as redes em um arquivo
# file = load("Resultados_30:09:2021/MLP-dt1-IP-points.jld2")

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
# %% show AC, MI for choosing optimal τ
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, Random

ττ = 0:90
ac = autocor(w, ττ)
smi = selfmutualinfo(w, ττ)
smi ./= maximum(smi)

fig = figure(figsize = (figx/2, figy))
ax = gca()
ax.plot(ττ, ac; label = "AC")
ax.plot(ττ, smi, label = "SMI", color = "C3")
τ = estimate_delay(w, "mi_min")
ax.scatter(τ, smi[τ]; color = "C2", s = 100, zorder = 99)
ax.set_ylim(0,1.05)
ax.set_xticks(0:30:90)
ax.set_xlabel("\$\\tau\$", labelpad = -20)
ax.legend()
fig.tight_layout(;pad = 0.25)
wsave(plotsdir("6", "mutualinfo_ac_tau"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
export
    Term,
    creator,
    annihilator,
    occupation,
    vacancy

"""
Represent product of fermionic creation/annihilation operators.

For a Fock space of flavours, this is an abstract representation of the product
of an arbitrary number of creation and annihilation operators, not necessarily
normal ordered.  One can write each such product in a canonical form:

    t(v,i,j,k) = v * prod(x -> c[x]', i) * prod(x -> c[x], reverse(j)) *
                 prod(x -> c[x]*c[x]', k)

where `i`, `j` are tuples of flavours, each one ordered ascendingly, and `k`
is a tuple of flavours disjoint from both `i` and `j`.  Note that above
expression is not a normal-ordered product either, but a "classification" of
each flavour with respect to the term's effect.

Instead of the `O(2^(2n))` dense or `O(2^n)` sparse storage required for the
elements of such a product, this class stores only an `O(1)` set of bitfields
allowing for fast on-the-fly computation.
"""
struct Term{T}
    mask::UInt64
    left::UInt64
    right::UInt64
    change::UInt64
    sign_mask::UInt64
    value::T
end

"""Conversion between terms types"""
Term{T}(term::Term{U}) where {T, U} = Term(term, U(term.value))

"""Convert term to itself"""
Term{T}(term::Term{T}) where {T} = term

"""Create a term from template, but assign new value to it"""
Term(term::Term, new_value) =
    Term(term.mask, term.left, term.right, term.change, term.sign_mask, new_value)

"""Constant term (energy shift)"""
Term(value::Number) =
    Term(UInt64(0), UInt64(0), UInt64(0), UInt64(0), UInt64(0), value)

"""Make creation operator"""
function creator(i::Integer; value=Int8(1), statmask::UInt64=~UInt64(0))
    @boundscheck (i >= 1 && i <= 64) || throw(BoundsError(1:64, i))
    mask = UInt64(1) << (i - 1)
    signmask = (mask - one(mask)) & statmask
    return Term(mask, mask, zero(mask), mask, signmask, value)
end

"""Make annihilation operator"""
function annihilator(i::Integer; value=Int8(1), statmask::UInt64=~UInt64(0))
    @boundscheck (i >= 1 && i <= 64) || throw(BoundsError(1:64, i))
    mask = UInt64(1) << (i - 1)
    signmask = (mask - one(mask)) & statmask
    return Term(mask, zero(mask), mask, mask, signmask, value)
end

"""Make occupation (number) operator"""
function occupation(i::Integer; value=Int8(1))
    @boundscheck (i >= 1 && i <= 64) || throw(BoundsError(1:64, i))
    mask = UInt64(1) << (i - 1)
    return Term(mask, mask, mask, zero(mask), zero(mask), value)
end

"""Make vacancy (1 - number) operator"""
function vacancy(i::Integer; value=Int8(1))
    @boundscheck (i >= 1 && i <= 64) || throw(BoundsError(1:64, i))
    mask = UInt64(1) << (i - 1)
    return Term(mask, zero(mask), zero(mask), zero(mask), zero(mask), value)
end

"""Multiplies two terms"""
function Base.:*(a::Term, b::Term)
    # Get value first, mainly to determine the type of result and ensure
    # type stability.  We have to make sure to get a signed type here, as
    # the sign may flip later (this is handled in a specialization).
    value = a.value * b.value

    # First, take care of the Pauli principle.  This is relatively easy, as
    # a.right and b.left are the "demands" of the two terms on their in-between
    # state, each confined to their mask.  This means we simply have to fulfill
    # both demands whenever the masks intersect.
    intersect = a.mask & b.mask
    if xor(a.right, b.left) & intersect != 0
        return Term(zero(value))
    end

    # Now that we know that the term is valid, let us determine the left and
    # right outer states: a.left and b.right are again demands on those states,
    # each confined to their mask.  Outside the mask, the other term may make
    # demands "through" the adjacent term.
    mask = a.mask | b.mask
    right = b.right | (a.right & ~b.mask)
    left = a.left | (b.left & ~a.mask)
    change = xor(left, right)

    # Each fundamental operator comes with a sign mask extending to the less
    # significant side of it.  Ignoring the flavours affected by the term for
    # the moment, multiplying terms simply means xor'ing the masks.
    sign_mask = xor(a.sign_mask, b.sign_mask)

    # Note that here, the order matters: for a normal-ordered term, we assume
    # that:
    #
    #       ψ_L = c[i[1]]' * ... * c[i[n]]' * c[j[n]] * ... * c[j[1]] ψ_R
    #
    # where both i and j are ordered ascendingly. This is because we can then
    # split this equation into two parts:
    #
    #       ψ_L = c[i[1]]' * ... * c[i[n]]' ψ_0
    #       ψ_R = c[j[1]]' * ... * c[j[n]]' ψ_0
    #
    # This means that for determining the sign, the intermediate state ψ_0 is
    # relevant, and so for a sign computation we must exclude any flavour in
    # ψ_L and ψ_R affected by the creators/annihilators.
    sign_mask &= ~mask

    # Finally we have to restore normal ordering by permuting operators.  We
    # know this may only yield a global sign.  So, we simply observe the
    # difference in effect of the product and normal-ordered product on a
    # single contributing "trial" state.  We choose the trial to be `right`,
    # as we know the normal-ordered product contributes and does not produce
    # a sign.
    trial = right
    permute = trial & b.sign_mask
    trial = xor(trial, b.change)
    permute = xor(permute, trial & a.sign_mask)
    value = copy_parity(permute, value)

    return Term(mask, left, right, change, sign_mask, value)
end

# This ensures type-stability is observed.
Base.:*(a::Term{A}, b::Term{B}) where {A <: Unsigned, B <: Unsigned} =
    Term(a, signed(a.value)) * Term(b, signed(b.value))

"""Check if state is mapped when multiplied by term from the left"""
is_mapped_right(t::Term, state) = state & t.mask == t.right

"""Check if state is mapped when multiplied by term from the right"""
is_mapped_left(t::Term, state) = state & t.mask == t.left

"""Change sign of number if parity of a bit pattern is odd"""
copy_parity(pattern::Integer, number) =
    (count_ones(pattern) & 1 != 0) ? -number : number

# We use this to ensure type stability: otherwise, passing an unsigned
# number would have to return a union, as -x converts it to a signed one.
# We could silently use signed() here, but it is better to handle this
# explicitly in application code.
copy_parity(::Integer, ::Unsigned) =
    throw(DomainError("Cannot copy parity sign to unsigned number"))

"""Map a state through term, assuming that it is indeed mapped"""
map_state_right(t::Term, state) =
    xor(state, t.change), copy_parity(state & t.sign_mask, t.value)

"""Map a state through term, assuming that it is indeed mapped"""
map_state_left(t::Term, state) = map_state_right(t, state)

"""
Return the state shift through the term.

Index of the side diagonal where the values of term are nonzero, or,
equivalently, a state with index `i` is mapped by `term` to either nothing
or a state with index `i + state_shift(term)`.
"""
state_shift(t::Term) = (t.left - t.right) % Int64

"""Check if terms can be added and give one term"""
addable(a::Term, b::Term) =
    a.mask == b.mask && a.right == b.right && a.left == b.left

"""Convert terms with values of one type to another"""
Base.convert(::Type{Term{T}}, term::Term) where T =
    Term(term, convert(T, term.value))

"""Promotion rules"""
Base.promote_rule(::Type{Term{T}}, ::Type{Term{U}}) where {T, U} =
    Term{promote_type(T, U)}

"""Promotion rules"""
Base.promote_rule(::Type{Term{T}}, ::Type{Term{U}}) where {T <: Unsigned, U <: Unsigned} =
    Term{promote_type(signed(T), signed(U))}

"""Explicit plus of term"""
Base.:+(term::Term) = Term(term, +term.value)

"""Negate term"""
Base.:-(term::Term) = Term(term, -term.value)

"""Scale term by a scalar factor"""
Base.:*(factor::Number, term::Term) = Term(term, factor * term.value)

"""Scale term by a scalar factor"""
Base.:*(term::Term, factor::Number) = Term(term, term.value * factor)

"""Add terms with the same operators"""
Base.:+(a::Term, b::Term) = Term(_check_addable(a, b), a.value + b.value)

"""Subtract terms with the same operators"""
Base.:-(a::Term, b::Term) = Term(_check_addable(a, b), a.value - b.value)

"""Check if terms are equal"""
Base.:(==)(lhs::Term, rhs::Term) = addable(lhs, rhs) && lhs.value == rhs.value

"""Check if terms are approximately equal"""
function Base.isapprox(lhs::Term, rhs::Term; atol::Real=0,
                       rtol::Real=Base.rtoldefault(lhs.value, rhs.value, atol),
                       nans::Bool=false)
    return addable(lhs, rhs) && isapprox(lhs.value, rhs.value; atol, rtol, nans)
end

"""Get the transpose of a term"""
Base.transpose(t::Term) =
    Term(t.mask, t.right, t.left, t.change, t.sign_mask, t.value)

"""Get the adjoint (conjugate transpose) of a term"""
Base.adjoint(t::Term) =
    Term(t.mask, t.right, t.left, t.change, t.sign_mask, conj(t.value))

"""Get term like template, but with zero value"""
Base.zero(t::Term) = Term(t, zero(t.value))

"""Check if the value of term is zero"""
Base.iszero(t::Term) = iszero(t.value)

function _check_addable(a::Term, b::Term)
    if !addable(a, b)
        throw(ArgumentError("can only add terms which have same operators"))
    end
    return a
end


text:
using Test

using CompScienceMeshes
using BEAST

"""
    isdivconforming(space)

Returns for every basis functions in the space the total of the absolute
values of non-compensated skeleton flux, i.e. flux flowing into from one
side of an edge and not flowing out on the other side. Should be all zeros
for divergence conforming functions, at least on a closed surface.
"""
function isdivconforming(space)

    geo = geometry(space)
    mesh = geo

    edges = skeleton(mesh,1)
    D = connectivity(edges, mesh, abs)
    rows = rowvals(D)
    vals = nonzeros(D)

    Flux = zeros(Float64, numcells(mesh),3)
    TotalFlux = zeros(Float64, numcells(edges), numfunctions(space))

    for i in 1 : numfunctions(space)

        fill!(Flux, 0)
        bfs = BEAST.basisfunction(space,i)
        for bf in bfs
            c = bf.cellid
            e = bf.refid
            x = bf.coeff
            Flux[c,e] += x
        end

        for E in 1 : numcells(edges)
            for j in nzrange(D,E)
                F = rows[j]
                e = vals[j]
                @assert 1 <= e <= 3
                TotalFlux[E,i] += Flux[F,e]
            end
        end
    end

    return TotalFlux
end

function touches_predicate(mesh)

    verts = skeleton(mesh,0)
    verts = sort(verts.faces, lt=isless)

    edges = skeleton(mesh,1)
    edges = sort(edges.faces, lt=isless)

    function pred(s)

        num_hits = 0
        for i in 1:length(s)
            if !isempty(searchsorted(verts, SVector(s[i]), lt=isless))
                num_hits += 1
            end
        end

        @assert num_hits < 3
        num_hits == 0 && return false
        num_hits == 1 && return true
        isempty(searchsorted(edges, s, lt=isless)) && return true
        return false
    end

    return pred
end

function interior(mesh::Mesh)
    D = dimension(mesh)
    edges = skeleton(mesh, D-1)
    pred = interior_tpredicate(mesh)
    submesh(pred, edges)
end

#meshfile = Pkg.dir("BEAST","test","sphere2.in")
meshfile = joinpath(dirname(@__FILE__),"assets","sphere316.in")
mesh = readmesh(meshfile)
@test numvertices(mesh) == 160
@test numcells(mesh) == 316

rt = raviartthomas(mesh)
@test numfunctions(rt) == 316 * 3 / 2

fine = barycentric_refinement(mesh)
edges = skeleton(mesh, 1)

bc = buffachristiansen(mesh)
@test numfunctions(bc) == 316 * 3 / 2

lc = isdivconforming(rt)
@test maximum(lc) < eps(Float64) * 1000
println("RT space is div-conforming")

lc = isdivconforming(bc);
@test maximum(lc) < eps(Float64) * 1000
println("BC space is div-conforming")

# Now repeat the exercise with an open mesh
mesh = meshrectangle(1.0, 1.0, 0.2);
fine = barycentric_refinement(mesh);

rt = raviartthomas(mesh)
bc = buffachristiansen(mesh)

@test numfunctions(rt) == 65
@test numfunctions(bc) == 65

int_pred = interior_tpredicate(mesh)
bnd_pred(s) = !int_pred(s)

leaky_edges = findall(sum(abs.(isdivconforming(rt)),dims=1) .!= 0)
@test length(leaky_edges) == 0

bnd = boundary(mesh)
bndtch_pred = touches_predicate(bnd)
edges = interior(mesh)
bndtch_edges = findall(bndtch_pred, edges.faces)
leaky_edges = findall(vec(sum(abs.(isdivconforming(bc)),dims=1)) .!= 0)
@test bndtch_edges == leaky_edges

## Test the charge of BC functions
#meshfile = Pkg.dir("BEAST","test","sphere2.in")
meshfile = joinpath(dirname(@__FILE__),"assets","sphere316.in")
mesh = readmesh(meshfile)
bc = buffachristiansen(mesh)
fine = geometry(bc)
charges = zeros(numcells(fine))

for fn in bc.fns
    abs_charge = 0.0
    net_charge = 0.0
    fill!(charges,0)
    for  _sh in fn
        cellid = _sh.cellid

        #cell = simplex(vertices(fine, fine.faces[cellid]))
        net_charge += _sh.coeff
        charges[cellid] += _sh.coeff
    end
    abs_charge = sum(abs.(charges))
    @test net_charge + 1 ≈ 1
    @test abs_charge ≈ 2
end

# THe BC construction function should throw for non-oriented surfaces
# width, height, h = 1.0, 0.5, 0.05
G1 = meshrectangle(1.0, 1.0, 0.25)
G2 = CompScienceMeshes.rotate(G1, 0.5π * x̂)
G = CompScienceMeshes.weld(G1,G2)
@test_throws AssertionError buffachristiansen(G)
# G3 = CompScienceMeshes.rotate(G1, 1.0π * x̂)


text:
HelloWorld = "{o}=} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {========{o} {o}=} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {========{o} {o}=} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {========{o} {o}=} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {========{o} {o}=} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{========{o} {o}=} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{========{o} {o}=} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {========{o} {o}=} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{========{o} {o}=} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {========{o}
{o}=} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {========{o} {o}=} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {========{o} {o}=} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{========{o} {o}=} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o}
{o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {o}=}
{={o} {o}=} {={o} {o}=} {={o} {o}=} {={o} {========{o} "


text:
using PyPlot

figure(figsize=(4.5,3))
vlinear=0.2
x=[-1,-vlinear,vlinear,1]
y=[-1,-1,1,1]
plot(x,y,color="black",linewidth=2)
xticks([-vlinear,0,vlinear],["\$-v_\\mathrm{linear}\$","","\$+v_\\mathrm{linear}\$"])
yticks([-1,0,1],["\$-f_\\mathrm{constant}\$","\$0\$","-\$+f_\\mathrm{constant}\$"])
xlabel("\$v\$")
ylabel("\$f\$")
ylim(-1.1,1.1)
grid(true)
subplots_adjust(left=0.25, right=0.95, top=0.95, bottom=0.15)
savefig("friction.png")
# savefig("friction.pdf")

max_forks_repo_name: christiankral/ElectroMechanicalDrives


text:
using Test, LinearAlgebra, Distributions, AdvancedHMC, Random, ForwardDiff
using AdvancedHMC.Adaptation: WelfordVar, NaiveVar, WelfordCov, NaiveCov, add_sample!, get_var, get_cov, reset!

# Check that the estimated variance is approximately correct.
@testset "Online v.s. naive v.s. true var/cov estimation" begin
    D = 10
    T = Float64
    sz = (D,)
    n_samples = 100_000

    var_welford = WelfordVar(T, sz)
    var_naive = NaiveVar(T, sz)
    var_estimators = [var_welford, var_naive]
    cov_welford = WelfordCov(T, sz)
    cov_naive = NaiveCov(T, sz)
    cov_estimators = [cov_welford, cov_naive]
    estimators = [var_estimators..., cov_estimators...]

    for dist in [MvNormal(zeros(D), ones(D)), Dirichlet(D, 1)]
        for _ = 1:n_samples
            s = rand(dist)
            for estimator in estimators
                add_sample!(estimator, s)
            end
        end

        @test get_var(var_welford) ≈ get_var(var_naive) atol=0.1D
        for estimator in var_estimators
            @test get_var(estimator) ≈ var(dist) atol=0.1D
        end

        @test get_cov(cov_welford) ≈ get_cov(cov_naive) atol=0.1D^2
        for estimator in cov_estimators
            @test get_cov(estimator) ≈ cov(dist) atol=0.1D^2
        end

        for estimator in estimators
            reset!(estimator)
        end
    end
end

@testset "Preconditioner constructors" begin
    θ = [0.0, 0.0, 0.0, 0.0]
    pc1 = Preconditioner(UnitEuclideanMetric) # default dim = 2
    pc2 = Preconditioner(DiagEuclideanMetric)
    pc3 = Preconditioner(DenseEuclideanMetric)

    # Var adaptor dimention should be increased to length(θ) from 2
    AdvancedHMC.adapt!(pc1, θ, 1.)
    AdvancedHMC.adapt!(pc2, θ, 1.)
    AdvancedHMC.adapt!(pc3, θ, 1.)
    @test AdvancedHMC.Adaptation.getM⁻¹(pc2) == ones(length(θ))
    @test AdvancedHMC.Adaptation.getM⁻¹(pc3) == LinearAlgebra.diagm(0 => ones(length(θ)))
end

@testset "Stan HMC adaptors" begin
    θ = [0.0, 0.0, 0.0, 0.0]

    adaptor1 = StanHMCAdaptor(
        Preconditioner(UnitEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )
    adaptor2 = StanHMCAdaptor(
        Preconditioner(DiagEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )
    adaptor3 = StanHMCAdaptor(
        Preconditioner(DenseEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )
    for a in [adaptor1, adaptor2, adaptor3]
        AdvancedHMC.initialize!(a, 1_000)
        @test a.state.window_start == 76
        @test a.state.window_end == 950
        @test a.state.window_splits == [100, 150, 250, 450, 950]
        AdvancedHMC.adapt!(a, θ, 1.)
    end
    @test AdvancedHMC.Adaptation.getM⁻¹(adaptor2) == ones(length(θ))
    @test AdvancedHMC.Adaptation.getM⁻¹(adaptor3) == LinearAlgebra.diagm(0 => ones(length(θ)))

    @test_deprecated StanHMCAdaptor(
        1_000,
        Preconditioner(DiagEuclideanMetric),
        NesterovDualAveraging(0.8, 0.5)
    )

    @testset "buffer > `n_adapts`" begin
        AdvancedHMC.initialize!(
            StanHMCAdaptor(
                Preconditioner(DenseEuclideanMetric),
                NesterovDualAveraging(0.8, 0.5)
            ),
            100
        )
    end
end

let D=10
    function runnuts(ℓπ, metric; n_samples=3_000)
        n_adapts = 1_500

        θ_init = rand(D)

        h = Hamiltonian(metric, ℓπ, ForwardDiff)
        prop = NUTS(Leapfrog(find_good_eps(h, θ_init)))
        adaptor = StanHMCAdaptor(
            Preconditioner(metric),
            NesterovDualAveraging(0.8, prop.integrator)
        )
        samples, stats = sample(h, prop, θ_init, n_samples, adaptor, n_adapts; verbose=false)
        return (samples=samples, stats=stats, adaptor=adaptor)
    end

    @testset "Adapted mass v.s. true variance" begin
        n_tests = 5

        @testset "DiagEuclideanMetric" begin
            for _ in 1:n_tests
                Random.seed!(1)

                # Random variance
                σ = ones(D) + abs.(randn(D))

                # Diagonal Gaussian
                target = MvNormal(zeros(D), σ)
                ℓπ = θ -> logpdf(target, θ)

                res = runnuts(ℓπ, DiagEuclideanMetric(D))
                @test res.adaptor.pc.var ≈ σ .^ 2 rtol=0.2

                res = runnuts(ℓπ, DenseEuclideanMetric(D))
                @test res.adaptor.pc.covar ≈ diagm(0 => σ .^ 2) rtol=0.25
            end
        end

        @testset "DenseEuclideanMetric" begin
            for _ in 1:n_tests
                # Random covariance
                m = randn(D, D)
                Σ = m' * m

                # Correlated Gaussian
                target = MvNormal(zeros(D), Σ)
                ℓπ = θ -> logpdf(target, θ)

                res = runnuts(ℓπ, DiagEuclideanMetric(D))
                @test res.adaptor.pc.var ≈ diag(Σ) rtol=0.2

                res = runnuts(ℓπ, DenseEuclideanMetric(D))
                @test res.adaptor.pc.covar ≈ Σ rtol=0.25
            end
        end

    end

    @testset "Initialisation adaptor by metric" begin
        target = MvNormal(zeros(D), 1)
        ℓπ = θ -> logpdf(target, θ)

        mass_init = fill(0.5, D)
        res = runnuts(ℓπ, DiagEuclideanMetric(mass_init); n_samples=1)
        @test res.adaptor.pc.var == mass_init

        mass_init = diagm(0 => fill(0.5, D))
        res = runnuts(ℓπ, DenseEuclideanMetric(mass_init); n_samples=1)
        @test res.adaptor.pc.covar == mass_init
    end
end

@testset "rand momentum variables from metric via vector of RNGs" begin
    D = 10
    n_chains = 5
    rng = [MersenneTwister(1) for _ in 1:n_chains]
    for metric in [
        UnitEuclideanMetric((D, n_chains)),
        DiagEuclideanMetric((D, n_chains)),
        # DenseEuclideanMetric((D, n_chains)) # not supported at the moment
    ]
        r = rand(rng, metric)
        all_same = true
        for i in 2:n_chains
            all_same = all_same && r[:,i] == r[:,1]
        end
        @test all_same
    end
end


text:
#=
definitions of real domains for polynomials
=#

"""
$(TYPEDEF)

Real domains for polynomial constraints.
"""
abstract type Domain{T <: Real} end

"""
$(SIGNATURES)

Dimension of a domain.
"""
function dimension(dom::Domain)::Int end

"""
$(SIGNATURES)

Sample `npts` points from the interior of the domain.
"""
function sample(dom::Domain, npts::Int)::Matrix end

"""
$(SIGNATURES)

Degree of a polynomial constraints defining a domain.
"""
function degree(dom::Domain)::Int end

"""
$(SIGNATURES)

Evaluations of the polynomial domain constraints at the points.
"""
function weights(dom::Domain, pts::AbstractMatrix)::Vector{Vector} end

"""
$(TYPEDEF)

Real vectors ``x ∈ ℝⁿ`` of dimension `n::Int`.
"""
mutable struct FreeDomain{T <: Real} <: Domain{T}
    n::Int
    function FreeDomain{T}(n::Int) where {T <: Real}
        @assert n >= 1
        dom = new{T}()
        dom.n = n
        return dom
    end
end

dimension(dom::FreeDomain) = dom.n

sample(dom::FreeDomain{T}, npts::Int) where {T <: Real} =
    sample(BoxDomain{T}(-ones(T, dom.n), ones(T, dom.n)), npts)::Matrix{T}

degree(::FreeDomain) = 0

weights(::FreeDomain{T}, ::AbstractMatrix{T}) where {T <: Real} = Vector{T}[]::Vector{Vector{T}}

"""
$(TYPEDEF)

Hyperbox ``x ∈ [l, u]`` with lower bounds `l::Vector{T}` and upper bounds
`u::Vector{T}`.
"""
mutable struct BoxDomain{T <: Real} <: Domain{T}
    l::Vector{T}
    u::Vector{T}
    function BoxDomain{T}(l::Vector{<:Real}, u::Vector{<:Real}) where {T <: Real}
        @assert length(l) == length(u)
        dom = new{T}()
        dom.l = l
        dom.u = u
        return dom
    end
end

dimension(dom::BoxDomain) = length(dom.l)

function sample(dom::BoxDomain{T}, npts::Int) where {T <: Real}
    dim = dimension(dom)
    pts = rand(T, npts, dim) .- T(0.5)
    shift = (dom.u + dom.l) .* T(0.5)
    for i in 1:npts
        @views pts[i, :] = pts[i, :] .* (dom.u - dom.l) + shift
    end
    return pts::Matrix{T}
end

degree(::BoxDomain) = 2

function weights(dom::BoxDomain{T}, pts::AbstractMatrix{T}) where {T <: Real}
    @views g = [(pts[:, i] .- dom.l[i]) .* (dom.u[i] .- pts[:, i]) for
        i in 1:size(pts, 2)]
    @assert all(all(gi .>= 0) for gi in g)
    return g::Vector{Vector{T}}
end

"""
$(TYPEDEF)

Hyperbox ``x ∈ [l, u]`` with lower bounds `l::Vector{T}` and upper bounds
`u::Vector{T}` and additional constraint functions.
"""
mutable struct CustomBoxDomain{T <: Real} <: Domain{T}
    l::Vector{T}
    u::Vector{T}
    gs::Vector{Function}
    deg::Int

    function CustomBoxDomain{T}(l::Vector{<:Real}, u::Vector{<:Real}, gs::Vector{Function}, deg::Int) where {T <: Real}
        @assert length(l) == length(u)
        dom = new{T}()
        dom.l = l
        dom.u = u
        dom.gs = gs
        dom.deg = max(deg, 2)
        return dom
    end
end

dimension(dom::CustomBoxDomain) = length(dom.l)

function sample(dom::CustomBoxDomain{T}, npts::Int) where {T <: Real}
    dim = dimension(dom)
    pts = Vector{Vector{T}}(undef, npts)
    shift = (dom.u + dom.l) .* T(0.5)
    scale = dom.u - dom.l
    k = 0
    randbox() = (rand(T, dim) .- T(0.5)) .* scale .+ shift
    while k < npts
        x = randbox()
        if all(g -> g(x...) > zero(T), dom.gs)
           k += 1
           pts[k] = x
        end
    end
    return reduce(vcat, transpose.(pts))::Matrix{T}
end

degree(dom::CustomBoxDomain) = dom.deg

function weights(dom::CustomBoxDomain{T}, pts::AbstractMatrix{T}) where {T <: Real}
   transpts = similar(transpose(pts))
   transpose!(transpts, pts)
   @views g = vcat([(pts[:, i] .- dom.l[i]) .* (dom.u[i] .- pts[:, i]) for i in 1:size(pts, 2)],
                   [[g(transpts[:, i]...) for i in 1:size(transpts, 2)] for g in dom.gs])
   @assert all(all(gi .>= 0) for gi in g)
   return g::Vector{Vector{T}}
end

# for hyperball and hyperellipse
function ball_sample(dom::Domain{T}, npts::Int) where {T <: Real}
    dim = dimension(dom)
    pts = T.(randn(npts, dim)) # randn doesn't work with all real types
    norms = sum(abs2, pts, dims = 2)
    pts ./= sqrt.(norms) # scale
    norms ./= 2
    gammainv = [gamma_inc(ai, dim / 2)[2] ^ inv(dim) for ai in norms]
    pts .*= gammainv
    return pts::Matrix{T}
end

"""
$(TYPEDEF)

Euclidean hyperball ``\\lVert (x-c) \\rVert_2 \\leq r`` with center
`c::Vector{T}` and positive radius `r::T`.
"""
mutable struct BallDomain{T <: Real} <: Domain{T}
    c::Vector{T}
    r::T
    function BallDomain{T}(c::Vector{T}, r::T) where {T <: Real}
        @assert r > 0
        dom = new{T}()
        dom.c = c
        dom.r = r
        return dom
    end
end

dimension(dom::BallDomain) = length(dom.c)

function sample(dom::BallDomain{T}, npts::Int) where {T <: Real}
    pts = ball_sample(dom, npts)
    pts .*= dom.r # scale
    pts .+= dom.c' # shift
    return pts::Matrix{T}
end

degree(::BallDomain) = 2

function weights(dom::BallDomain{T}, pts::AbstractMatrix{T}) where {T <: Real}
    g = [abs2(dom.r) - sum(abs2, pts[j, :] - dom.c) for j in 1:size(pts, 1)]
    @assert all(g .>= 0)
    return [g]::Vector{Vector{T}}
end

"""
$(TYPEDEF)

Hyperellipse ``(x-c)' Q (x-c) \\leq 1`` with center `c::Vector{T}` and symmetric
positive definite scaling/rotation matrix `Q::AbstractMatrix{T}`.
"""
mutable struct EllipsoidDomain{T <: Real} <: Domain{T}
    c::Vector{T}
    Q::AbstractMatrix{T}
    QU::UpperTriangular{T}
    function EllipsoidDomain{T}(
        c::Vector{T},
        Q::AbstractMatrix{T},
        ) where {T <: Real}
        @assert length(c) == size(Q, 1)
        @assert issymmetric(Q)
        F = cholesky(Q, check = false)
        @assert isposdef(F)
        dom = new{T}()
        dom.c = c
        dom.Q = Q
        dom.QU = F.U
        return dom
    end
end

dimension(dom::EllipsoidDomain) = length(dom.c)

function sample(dom::EllipsoidDomain{T}, npts::Int) where {T <: Real}
    pts = ball_sample(dom, npts)
    rdiv!(pts, dom.QU') # scale/rotate
    pts .+= dom.c' # shift
    return pts::Matrix{T}
end

degree(::EllipsoidDomain) = 2

function weights(
    dom::EllipsoidDomain{T},
    pts::AbstractMatrix{T},
    ) where {T <: Real}
    g = [1 - sum(abs2, dom.QU * (pts[j, :] - dom.c)) for j in 1:size(pts, 1)]
    @assert all(g .>= 0)
    return [g]::Vector{Vector{T}}
end


text:
using Oceananigans.Fields: ZFaceField, CenterField
using Oceananigans.TurbulenceClosures: AnisotropicMinimumDissipation

subfilter_diffusivity(::AnisotropicMinimumDissipation, diffusivity_fields, name) =
    getproperty(diffusivity_fields.κₑ, name)

"""
    subfilter_momentum_fluxes(model,
                              uz_scratch = Field{Face, Center, Face}(model.grid),
                              vz_scratch = Field{Center, Face, Face}(model.grid),
                              c_scratch = CenterField(model.grid))

Returns a dictionary of horizontally-averaged subfilter momentum fluxes.
"""
function subfilter_momentum_fluxes(model;
                                   uz_scratch = Field{Face, Center, Face}(model.grid),
                                   vz_scratch = Field{Center, Face, Face}(model.grid),
                                   c_scratch = CenterField(model.grid))

    u, v, w = model.velocities

    νₑ = model.diffusivity_fields.νₑ

    averages = Dict(
                    :νₑ_∂z_u => Field(Average(∂z(u) * νₑ, dims=(1, 2))),
                    :νₑ_∂z_v => Field(Average(∂z(v) * νₑ, dims=(1, 2))),
                    :νₑ_∂z_w => Field(Average(∂z(w) * νₑ, dims=(1, 2)))
                   )

    return averages
end

"""
    subfilter_tracer_fluxes(model, w_scratch=ZFaceField(model.grid))

Returns a dictionary of horizontally-averaged subfilter tracer fluxes.
"""
function subfilter_tracer_fluxes(model; w_scratch = ZFaceField(model.grid))

    averages = Dict()

    for tracer in propertynames(model.tracers)
        c = getproperty(model.tracers, tracer)
        κₑ = subfilter_diffusivity(model.closure, model.diffusivity_fields, tracer)

        name = Symbol(:κₑ_∂z_, tracer)

        averages[name] = Field(Average(∂z(c) * κₑ, dims=(1, 2)))
    end

    return averages
end

max_stars_repo_name: ali-ramadhan/LESbrary

max_issues_repo_name: ali-ramadhan/LESbrary

max_forks_repo_name: ali-ramadhan/LESbrary


text:
type QROperator{CO,MT,T} <: Operator{T}
    R::CO
    H::MT # Contains the Householder reflections
    ncols::Int   # number of cols already upper triangularized
end

QROperator(R::CachedOperator,H::AbstractArray,ncs::Int) =
    QROperator{typeof(R),typeof(H),eltype(H)}(R,H,ncs)

for OP in (:domainspace,:rangespace)
    @eval $OP(QR::QROperator) = $OP(QR.R)
end

getindex(QR::QROperator,k::Integer,j::Integer) = (QR[:Q]*QR[:R])[k,j]

immutable QROperatorR{QRT,T} <: Operator{T}
    QR::QRT
end

QROperatorR(QR) = QROperatorR{typeof(QR),eltype(QR)}(QR)
domainspace(R::QROperatorR) = domainspace(R.QR)
rangespace(R::QROperatorR) = ℓ⁰

function getindex(R::QROperatorR,k::Integer,j::Integer)
    if j < k
        zero(eltype(R))
    else
        resizedata!(R.QR,:,j)
        R.QR.R[k,j]
    end
end

bandinds(R::QROperatorR) = 0,bandinds(R.QR.R,2)

immutable QROperatorQ{QRT,T} <: Operator{T}
    QR::QRT
end

QROperatorQ(QR) = QROperatorQ{typeof(QR),eltype(QR)}(QR)

domainspace(Q::QROperatorQ) = ℓ⁰
rangespace(Q::QROperatorQ) = rangespace(Q.QR)

getindex(Q::QROperatorQ,k::Integer,j::Integer) = (Q'*eltype(Q)[zeros(k-1);1])[j]

function getindex(QR::QROperator,d::Symbol)
    d==:Q && return QROperatorQ(QR)
    d==:R && return QROperatorR(QR)

    error("Symbol not recognized")
end

function QROperator{T}(R::CachedOperator{T,AlmostBandedMatrix{T}})
    M = R.data.bands.l+1   # number of diag+subdiagonal bands
    H = Array(T,M,100)
    QROperator(R,H,0)
end

function QROperator{T}(R::CachedOperator{T,BandedMatrix{T}})
    M = R.data.l+1   # number of diag+subdiagonal bands
    H = Array(T,M,100)
    QROperator(R,H,0)
end

QROperator{T}(R::CachedOperator{T,RaggedMatrix{T}}) =
    QROperator(R,RaggedMatrix(T,0,Int[]),0)

function QROperator{T,AM<:AbstractMatrix}(R::CachedOperator{T,AM})
    error("Cannot create a QR factorization for $R")
end

function Base.qrfact(A::Operator)
    if isambiguous(domainspace(A)) || isambiguous(rangespace(A))
        throw(ArgumentError("Only non-ambiguous operators can be factorized."))
    end
    QROperator(cache(A;padding=true))
end

function Base.qr(A::Operator)
    QR = qrfact(A)
    QR[:Q],QR[:R]
end

for OP in (:(Base.qrfact),:(Base.qr))
    @eval $OP{OO<:Operator}(A::Array{OO}) = $OP(interlace(A))
end

function Base.det(R::QROperatorR;maxiterations::Int=10_000)
    QR = R.QR
    RD = R.data
    resizedata!(QR,:,1)
    ret = -RD[1,1]
    k = 2
    while abs(abs(RD[k-1,k-1])-1) > eps(eltype(A))
        resizedata!(QR,:,k)
        ret *= -RD[k,k]
        k+=1
        k > maxiterations && error("Determinant unlikely to converge after 10_000 iterations.")
    end

    ret
end
Base.det(QR::QROperatorQ) = 1

Base.det(QR::QROperator) = det(QR[:R])
Base.det(A::Operator) = det(qrfact(A))

## Multiplication routines

linsolve{CO,MT,T<:Real}(QR::QROperator{CO,MT,T},b::Vector{T};kwds...) =
    Fun(QR[:R]\Ac_mul_B(QR[:Q],b;kwds...),domainspace(QR))
linsolve{CO,MT,T<:Complex}(QR::QROperator{CO,MT,T},b::Vector{T};kwds...) =
    Fun(QR[:R]\Ac_mul_B(QR[:Q],b;kwds...),domainspace(QR))

linsolve{CO,MT,T,V<:Number}(QR::QROperator{CO,MT,T},b::Vector{V};kwds...) =
    linsolve(QR,Vector{T}(b);kwds...)

linsolve{CO,MT,T<:Real,V<:Complex}(QR::QROperator{CO,MT,T},b::Vector{V};kwds...) =
    linsolve(QR,real(b);kwds...)+im*linsolve(QR,imag(b);kwds...)
linsolve{CO,MT,T<:Complex,V<:Real}(QR::QROperator{CO,MT,T},b::Vector{V};kwds...) =
    linsolve(QR,Vector{T}(b);kwds...)

linsolve(QR::QROperator,b::Fun;kwds...) = linsolve(QR,coefficients(b,rangespace(QR));kwds...)
function linsolve(QR::QROperator,b::Vector{Any};kwds...)
    #TODO: PDEQR remove this is a hack
    if length(b) == 1 && isa(b[1],Fun)
        linsolve(QR,Fun(b[1],rangespace(QR));kwds...)
    else
        linsolve(QR,Fun(b,rangespace(QR));kwds...)
    end
end
linsolve{FF<:Fun}(QR::QROperator,b::Vector{FF};kwds...) = linsolve(QR,Fun(b,rangespace(QR));kwds...)
function linsolve(A::QROperator,B::Matrix;kwds...)
    ds=domainspace(A)
    ret=Array(Fun{typeof(ds),promote_type(mapreduce(eltype,promote_type,B),eltype(ds))},1,size(B,2))
    for j=1:size(B,2)
        ret[:,j]=linsolve(A,B[:,j];kwds...)
    end
    demat(ret)
end
linsolve(A::QROperator,b;kwds...) = linsolve(A,Fun(b);kwds...)

Base.Ac_mul_B(A::QROperatorQ,b::Vector{Any};kwds...) = Ac_mul_B(A,Fun(b,rangespace(A));kwds...)
Base.Ac_mul_B{FF<:Fun}(A::QROperatorQ,b::Vector{FF};kwds...) = Ac_mul_B(A,Fun(b,rangespace(A));kwds...)

Base.At_mul_B{T<:Real}(A::QROperatorQ{T},B::Union{Vector{T},Matrix{T}}) = Ac_mul_B(A,B)

Base.Ac_mul_B(A::QROperatorQ,B::Vector;tolerance=eps(eltype(A))/10,maxlength=1000000) =
        Ac_mul_Bpars(A,B,tolerance,maxlength)

Base.Ac_mul_B{QR,T,V<:Number}(A::QROperatorQ{QR,T},B::AbstractVector{V};opts...) =
    Ac_mul_B(A,Vector{T}(B))

Base.Ac_mul_B{QR,T}(A::QROperatorQ{QR,T},B::Fun;opts...) =
    Ac_mul_B(A,coefficients(B,rangespace(A)))

function linsolve(R::QROperatorR,b::Vector)
    if length(b) > R.QR.ncols
        # upper triangularize columns
        resizedata!(R.QR,:,length(b))
    end
    Fun(trtrs!(Val{'U'},R.QR.R,copy(b)),domainspace(R))
end

linsolve(R::QROperatorR,b::Fun{SequenceSpace};kwds...) = linsolve(R,b.coefficients;kwds...)
linsolve(A::QROperatorR,b::Fun;kwds...) = error("linsolve not implement for $(typeof(b)) right-hand sides")
linsolve(A::QROperatorR,b;kwds...) = linsolve(A,Fun(b);kwds...)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module RoundingTests

using Test
using Dates

@testset "conversion to and from the rounding epoch (ISO 8601 year 0000)" begin
    @test DatesPlus.epochdays2date(-1) == DatesPlus.Date(-1, 12, 31)
    @test DatesPlus.epochdays2date(0) == DatesPlus.Date(0, 1, 1)
    @test DatesPlus.epochdays2date(1) == DatesPlus.Date(0, 1, 2)
    @test DatesPlus.epochdays2date(736329) == DatesPlus.Date(2016, 1, 1)
    @test DatesPlus.epochms2datetime(-86400000) == DatesPlus.DateTime(-1, 12, 31)
    @test DatesPlus.epochms2datetime(0) == DatesPlus.DateTime(0, 1, 1)
    @test DatesPlus.epochms2datetime(86400000) == DatesPlus.DateTime(0, 1, 2)
    @test DatesPlus.epochms2datetime(Int64(736329) * 86400000) == DatesPlus.DateTime(2016, 1, 1)
    @test DatesPlus.date2epochdays(DatesPlus.Date(-1, 12, 31)) == -1
    @test DatesPlus.date2epochdays(DatesPlus.Date(0, 1, 1)) == 0
    @test DatesPlus.date2epochdays(DatesPlus.Date(2016, 1, 1)) == 736329
    @test DatesPlus.datetime2epochms(DatesPlus.DateTime(-1, 12, 31)) == -86400000
    @test DatesPlus.datetime2epochms(DatesPlus.DateTime(0, 1, 1)) == 0
    @test DatesPlus.datetime2epochms(DatesPlus.DateTime(2016, 1, 1)) == Int64(736329) * 86400000
end
@testset "Basic rounding" begin
    dt = DatesPlus.Date(2016, 2, 28)    # Sunday
    @test floor(dt, DatesPlus.Year) == DatesPlus.Date(2016)
    @test floor(dt, DatesPlus.Year(5)) == DatesPlus.Date(2015)
    @test floor(dt, DatesPlus.Year(10)) == DatesPlus.Date(2010)
    @test floor(dt, DatesPlus.Quarter) == DatesPlus.Date(2016, 1)
    @test floor(dt, DatesPlus.Quarter(6)) == DatesPlus.Date(2016, 1)
    @test floor(dt, DatesPlus.Month) == DatesPlus.Date(2016, 2)
    @test floor(dt, DatesPlus.Month(6)) == DatesPlus.Date(2016, 1)
    @test floor(dt, DatesPlus.Week) == DatesPlus.toprev(dt, DatesPlus.Monday)
    @test ceil(dt, DatesPlus.Year) == DatesPlus.Date(2017)
    @test ceil(dt, DatesPlus.Year(5)) == DatesPlus.Date(2020)
    @test ceil(dt, DatesPlus.Quarter) == DatesPlus.Date(2016, 4)
    @test ceil(dt, DatesPlus.Quarter(6)) == DatesPlus.Date(2017, 7)
    @test ceil(dt, DatesPlus.Month) == DatesPlus.Date(2016, 3)
    @test ceil(dt, DatesPlus.Month(6)) == DatesPlus.Date(2016, 7)
    @test ceil(dt, DatesPlus.Week) == DatesPlus.tonext(dt, DatesPlus.Monday)
    @test round(dt, DatesPlus.Year) == DatesPlus.Date(2016)
    @test round(dt, DatesPlus.Month) == DatesPlus.Date(2016, 3)
    @test round(dt, DatesPlus.Week) == DatesPlus.Date(2016, 2, 29)

    dt = DatesPlus.DateTime(2016, 2, 28, 15, 10, 50, 500)
    @test floor(dt, DatesPlus.Day) == DatesPlus.DateTime(2016, 2, 28)
    @test floor(dt, DatesPlus.Hour) == DatesPlus.DateTime(2016, 2, 28, 15)
    @test floor(dt, DatesPlus.Hour(2)) == DatesPlus.DateTime(2016, 2, 28, 14)
    @test floor(dt, DatesPlus.Hour(12)) == DatesPlus.DateTime(2016, 2, 28, 12)
    @test floor(dt, DatesPlus.Minute) == DatesPlus.DateTime(2016, 2, 28, 15, 10)
    @test floor(dt, DatesPlus.Minute(15)) == DatesPlus.DateTime(2016, 2, 28, 15, 0)
    @test floor(dt, DatesPlus.Second) == DatesPlus.DateTime(2016, 2, 28, 15, 10, 50)
    @test floor(dt, DatesPlus.Second(30)) == DatesPlus.DateTime(2016, 2, 28, 15, 10, 30)
    @test ceil(dt, DatesPlus.Day) == DatesPlus.DateTime(2016, 2, 29)
    @test ceil(dt, DatesPlus.Hour) == DatesPlus.DateTime(2016, 2, 28, 16)
    @test ceil(dt, DatesPlus.Hour(2)) == DatesPlus.DateTime(2016, 2, 28, 16)
    @test ceil(dt, DatesPlus.Hour(12)) == DatesPlus.DateTime(2016, 2, 29, 0)
    @test ceil(dt, DatesPlus.Minute) == DatesPlus.DateTime(2016, 2, 28, 15, 11)
    @test ceil(dt, DatesPlus.Minute(15)) == DatesPlus.DateTime(2016, 2, 28, 15, 15)
    @test ceil(dt, DatesPlus.Second) == DatesPlus.DateTime(2016, 2, 28, 15, 10, 51)
    @test ceil(dt, DatesPlus.Second(30)) == DatesPlus.DateTime(2016, 2, 28, 15, 11, 0)
    @test round(dt, DatesPlus.Day) == DatesPlus.DateTime(2016, 2, 29)
    @test round(dt, DatesPlus.Hour) == DatesPlus.DateTime(2016, 2, 28, 15)
    @test round(dt, DatesPlus.Hour(2)) == DatesPlus.DateTime(2016, 2, 28, 16)
    @test round(dt, DatesPlus.Hour(12)) == DatesPlus.DateTime(2016, 2, 28, 12)
    @test round(dt, DatesPlus.Minute) == DatesPlus.DateTime(2016, 2, 28, 15, 11)
    @test round(dt, DatesPlus.Minute(15)) == DatesPlus.DateTime(2016, 2, 28, 15, 15)
    @test round(dt, DatesPlus.Second) == DatesPlus.DateTime(2016, 2, 28, 15, 10, 51)
    @test round(dt, DatesPlus.Second(30)) == DatesPlus.DateTime(2016, 2, 28, 15, 11, 0)
end
@testset "Rounding for dates at the rounding epoch (year 0000)" begin
    dt = DatesPlus.DateTime(0)
    @test floor(dt, DatesPlus.Year) == dt
    @test floor(dt, DatesPlus.Month) == dt
    @test floor(dt, DatesPlus.Week) == DatesPlus.Date(-1, 12, 27)   # Monday prior to 0000-01-01
    @test floor(DatesPlus.Date(-1, 12, 27), DatesPlus.Week) == DatesPlus.Date(-1, 12, 27)
    @test floor(dt, DatesPlus.Day) == dt
    @test floor(dt, DatesPlus.Hour) == dt
    @test floor(dt, DatesPlus.Minute) == dt
    @test floor(dt, DatesPlus.Second) == dt
    @test ceil(dt, DatesPlus.Year) == dt
    @test ceil(dt, DatesPlus.Month) == dt
    @test ceil(dt, DatesPlus.Week) == DatesPlus.Date(0, 1, 3)       # Monday following 0000-01-01
    @test ceil(DatesPlus.Date(0, 1, 3), DatesPlus.Week) == DatesPlus.Date(0, 1, 3)
    @test ceil(dt, DatesPlus.Day) == dt
    @test ceil(dt, DatesPlus.Hour) == dt
    @test ceil(dt, DatesPlus.Minute) == dt
    @test ceil(dt, DatesPlus.Second) == dt
end
@testset "Rounding for multiples of a period" begin
    # easiest to test close to rounding epoch
    dt = DatesPlus.DateTime(0, 1, 19, 19, 19, 19, 19)
    @test floor(dt, DatesPlus.Year(2)) == DateTime(0)
    @test floor(dt, DatesPlus.Month(2)) == DateTime(0, 1)       # Odd number; months are 1-indexed
    @test floor(dt, DatesPlus.Week(2)) == DateTime(0, 1, 17)    # Third Monday of 0000
    @test floor(dt, DatesPlus.Day(2)) == DateTime(0, 1, 19)     # Odd number; days are 1-indexed
    @test floor(dt, DatesPlus.Hour(2)) == DateTime(0, 1, 19, 18)
    @test floor(dt, DatesPlus.Minute(2)) == DateTime(0, 1, 19, 19, 18)
    @test floor(dt, DatesPlus.Second(2)) == DateTime(0, 1, 19, 19, 19, 18)
    @test ceil(dt, DatesPlus.Year(2)) == DateTime(2)
    @test ceil(dt, DatesPlus.Month(2)) == DateTime(0, 3)        # Odd number; months are 1-indexed
    @test ceil(dt, DatesPlus.Week(2)) == DateTime(0, 1, 31)     # Fifth Monday of 0000
    @test ceil(dt, DatesPlus.Day(2)) == DateTime(0, 1, 21)      # Odd number; days are 1-indexed
    @test ceil(dt, DatesPlus.Hour(2)) == DateTime(0, 1, 19, 20)
    @test ceil(dt, DatesPlus.Minute(2)) == DateTime(0, 1, 19, 19, 20)
    @test ceil(dt, DatesPlus.Second(2)) == DateTime(0, 1, 19, 19, 19, 20)
end
@testset "Rounding for dates with negative years" begin
    dt = DatesPlus.DateTime(-1, 12, 29, 19, 19, 19, 19)
    @test floor(dt, DatesPlus.Year(2)) == DateTime(-2)
    @test floor(dt, DatesPlus.Month(2)) == DateTime(-1, 11)     # Odd number; months are 1-indexed
    @test floor(dt, DatesPlus.Week(2)) == DateTime(-1, 12, 20)  # 2 weeks prior to 0000-01-03
    @test floor(dt, DatesPlus.Day(2)) == DateTime(-1, 12, 28)   # Even; 4 days prior to 0000-01-01
    @test floor(dt, DatesPlus.Hour(2)) == DateTime(-1, 12, 29, 18)
    @test floor(dt, DatesPlus.Minute(2)) == DateTime(-1, 12, 29, 19, 18)
    @test floor(dt, DatesPlus.Second(2)) == DateTime(-1, 12, 29, 19, 19, 18)
    @test ceil(dt, DatesPlus.Year(2)) == DateTime(0)
    @test ceil(dt, DatesPlus.Month(2)) == DateTime(0, 1)        # Odd number; months are 1-indexed
    @test ceil(dt, DatesPlus.Week(2)) == DateTime(0, 1, 3)      # First Monday of 0000
    @test ceil(dt, DatesPlus.Day(2)) == DateTime(-1, 12, 30)    # Even; 2 days prior to 0000-01-01
    @test ceil(dt, DatesPlus.Hour(2)) == DateTime(-1, 12, 29, 20)
    @test ceil(dt, DatesPlus.Minute(2)) == DateTime(-1, 12, 29, 19, 20)
    @test ceil(dt, DatesPlus.Second(2)) == DateTime(-1, 12, 29, 19, 19, 20)
end
@testset "Rounding for dates that should not need rounding" begin
    for dt in [DatesPlus.DateTime(2016, 1, 1), DatesPlus.DateTime(-2016, 1, 1)]
        local dt
        for p in [DatesPlus.Year, DatesPlus.Month, DatesPlus.Day, DatesPlus.Hour, DatesPlus.Minute, DatesPlus.Second]
            local p
            @test floor(dt, p) == dt
            @test ceil(dt, p) == dt
        end
    end
end
@testset "Various available RoundingModes" begin
    dt = DatesPlus.DateTime(2016, 2, 28, 12)
    @test round(dt, DatesPlus.Day, RoundNearestTiesUp) == DatesPlus.DateTime(2016, 2, 29)
    @test round(dt, DatesPlus.Day, RoundUp) == DatesPlus.DateTime(2016, 2, 29)
    @test round(dt, DatesPlus.Day, RoundDown) == DatesPlus.DateTime(2016, 2, 28)
    @test_throws DomainError round(dt, DatesPlus.Day, RoundNearest)
    @test_throws DomainError round(dt, DatesPlus.Day, RoundNearestTiesAway)
    @test_throws DomainError round(dt, DatesPlus.Day, RoundToZero)
    @test round(dt, DatesPlus.Day) == round(dt, DatesPlus.Day, RoundNearestTiesUp)
end
@testset "Rounding datetimes to invalid resolutions" begin
    dt = DatesPlus.DateTime(2016, 2, 28, 12, 15)
    for p in [DatesPlus.Year, DatesPlus.Month, DatesPlus.Week, DatesPlus.Day, DatesPlus.Hour]
        local p
        for v in [-1, 0]
            @test_throws DomainError floor(dt, p(v))
            @test_throws DomainError ceil(dt, p(v))
            @test_throws DomainError round(dt, p(v))
        end
    end
end
@testset "Rounding for periods" begin
    x = DatesPlus.Second(172799)
    @test floor(x, DatesPlus.Week) == DatesPlus.Week(0)
    @test floor(x, DatesPlus.Day) == DatesPlus.Day(1)
    @test floor(x, DatesPlus.Hour) == DatesPlus.Hour(47)
    @test floor(x, DatesPlus.Minute) == DatesPlus.Minute(2879)
    @test floor(x, DatesPlus.Second) == DatesPlus.Second(172799)
    @test floor(x, DatesPlus.Millisecond) == DatesPlus.Millisecond(172799000)
    @test ceil(x, DatesPlus.Week) == DatesPlus.Week(1)
    @test ceil(x, DatesPlus.Day) == DatesPlus.Day(2)
    @test ceil(x, DatesPlus.Hour) == DatesPlus.Hour(48)
    @test ceil(x, DatesPlus.Minute) == DatesPlus.Minute(2880)
    @test ceil(x, DatesPlus.Second) == DatesPlus.Second(172799)
    @test ceil(x, DatesPlus.Millisecond) == DatesPlus.Millisecond(172799000)
    @test round(x, DatesPlus.Week) == DatesPlus.Week(0)
    @test round(x, DatesPlus.Day) == DatesPlus.Day(2)
    @test round(x, DatesPlus.Hour) == DatesPlus.Hour(48)
    @test round(x, DatesPlus.Minute) == DatesPlus.Minute(2880)
    @test round(x, DatesPlus.Second) == DatesPlus.Second(172799)
    @test round(x, DatesPlus.Millisecond) == DatesPlus.Millisecond(172799000)

    x = DatesPlus.Nanosecond(2000999999)
    @test floor(x, DatesPlus.Second) == DatesPlus.Second(2)
    @test floor(x, DatesPlus.Millisecond) == DatesPlus.Millisecond(2000)
    @test floor(x, DatesPlus.Microsecond) == DatesPlus.Microsecond(2000999)
    @test floor(x, DatesPlus.Nanosecond) == x
    @test ceil(x, DatesPlus.Second) == DatesPlus.Second(3)
    @test ceil(x, DatesPlus.Millisecond) == DatesPlus.Millisecond(2001)
    @test ceil(x, DatesPlus.Microsecond) == DatesPlus.Microsecond(2001000)
    @test ceil(x, DatesPlus.Nanosecond) == x
    @test round(x, DatesPlus.Second) == DatesPlus.Second(2)
    @test round(x, DatesPlus.Millisecond) == DatesPlus.Millisecond(2001)
    @test round(x, DatesPlus.Microsecond) == DatesPlus.Microsecond(2001000)
    @test round(x, DatesPlus.Nanosecond) == x
end

@testset "Rouding DateTime to Date" begin
    now_ = DateTime(2020, 9, 1, 13)
    for p in (Year, Month, Day)
        for r in (RoundUp, RoundDown)
            @test round(Date, now_, p, r) == round(Date(now_), p, r)
        end
        @test round(Date, now_, p) == round(Date, now_, p, RoundNearestTiesUp)
        @test floor(Date, now_, p) == round(Date, now_, p, RoundDown)
        @test ceil(Date, now_, p)  == round(Date, now_, p, RoundUp)
    end
end
@testset "Rounding for periods that should not need rounding" begin
    for x in [DatesPlus.Week(3), DatesPlus.Day(14), DatesPlus.Second(604800)]
        local x
        for p in [DatesPlus.Week, DatesPlus.Day, DatesPlus.Hour, DatesPlus.Second, DatesPlus.Millisecond, DatesPlus.Microsecond, DatesPlus.Nanosecond]
            local p
            @test floor(x, p) == p(x)
            @test ceil(x, p) == p(x)
        end
    end
end
@testset "Various available RoundingModes for periods" begin
    x = DatesPlus.Hour(36)
    @test round(x, DatesPlus.Day, RoundNearestTiesUp) == DatesPlus.Day(2)
    @test round(x, DatesPlus.Day, RoundUp) == DatesPlus.Day(2)
    @test round(x, DatesPlus.Day, RoundDown) == DatesPlus.Day(1)
    @test_throws DomainError round(x, DatesPlus.Day, RoundNearest)
    @test_throws DomainError round(x, DatesPlus.Day, RoundNearestTiesAway)
    @test_throws DomainError round(x, DatesPlus.Day, RoundToZero)
    @test round(x, DatesPlus.Day) == round(x, DatesPlus.Day, RoundNearestTiesUp)
end
@testset "Rounding periods to invalid resolutions" begin
    x = DatesPlus.Hour(86399)
    for p in [DatesPlus.Week, DatesPlus.Day, DatesPlus.Hour, DatesPlus.Second, DatesPlus.Millisecond, DatesPlus.Microsecond, DatesPlus.Nanosecond]
        local p
        for v in [-1, 0]
            @test_throws DomainError floor(x, p(v))
            @test_throws DomainError ceil(x, p(v))
            @test_throws DomainError round(x, p(v))
        end
    end
    for p in [DatesPlus.Year, DatesPlus.Month]
        local p
        for v in [-1, 0, 1]
            @test_throws MethodError floor(x, p(v))
            @test_throws MethodError ceil(x, p(v))
            @test_throws DomainError round(x, p(v))
        end
    end
end

end


text:
using Unitful, FiniteElementDiffEq

### Setup
dx = (1//2^(3))u"m"
mesh = notime_squaremesh([0u"m" 1u"m" 0u"m" 1u"m"],dx,:dirichlet) #Fails at meshgrid because unit ranges not supported

f = (x) -> sin(2π.*map((y)->y.val,x[:,1])).*cos(2π.*map((y)->y.val,x[:,2]))
prob = PoissonProblem(f,mesh,D=1u"m^2/s")

#=
sol = solve(mesh,prob)
=#

#Define a parabolic problem
T = 1u"s"
dt = (1//2^(7))u"s"
mesh = parabolic_squaremesh([0u"m" 1u"m" 0u"m" 1u"m"],dx,dt,T,:dirichlet)
f = (u,x,t)  -> (float(ones(size(x,1))))u"N" - .5u
u0_func = (x) -> map((x)->(x)u"N",zeros(size(x,1)))
u0 = u0_func(mesh.node)
prob = HeatProblem(u0,f,mesh)
#=
println("Euler")
sol = solve(prob,alg=:Euler)
=#

true