text:
"""
    ndreducible(A => (:i, :j, ...), B => (...), ..., C => (...))

Given pairs of an array (or a `Broadcasted`) and indices,
`ndreducible` "couples" the arrays with shared index and plan the best
nested `for`-loop along all the indices.  `ndreducible` returns a
fold-able object which can be used with `foldl` and `foreach` defined
in Transducers.jl.  Use `referenceable` from Referenceables.jl to
mutate arrays.

See also: [`plan`](@ref).

# Examples

```math
C_{ij} = ∑_k A_{ik} B_{kj} + C_{ij}
```

```jldoctest
julia> using NDReducibles
       using Referenceables: referenceable

julia> A = rand(-10:10, 5, 4)
       B = rand(-10:10, 4, 3)
       C = zeros(Int, 5, 3);

julia> foreach(
           ndreducible(
               referenceable(C) => (:i, :j),
               A => (:i, :k),
               B => (:k, :j),
           ),
           simd = :ivdep,  # optional
       ) do (c, a, b)
           c[] += a * b
           return  # not required but useful for performance
       end;

julia> C == A * B
true
```

```math
∑_i A_i B_i C_i
```

```jldoctest
julia> using NDReducibles
       using Transducers: MapSplat

julia> A = rand(-10:10, 5)
       B = rand(-10:10, 5)
       C = rand(-10:10, 5);

julia> foldl(
           +,
           MapSplat(*),
           ndreducible(
               A => (:i,),
               B => (:i,),
               C => (:i,),
           ),
           simd = true,  # optional
       ) == sum(A .* B .* C)
true
```

```math
∑ \\{ A_i B_i C_i \\;|\\; A_i > 0 \\}
```

```jldoctest
julia> using NDReducibles
       using Transducers: MapSplat, Filter

julia> A = rand(-10:10, 5)
       B = rand(-10:10, 5)
       C = rand(-10:10, 5);

julia> foldl(
           +,
           Filter(((a, b, c),) -> a > 0) |> MapSplat(*),
           ndreducible(
               A => (:i,),
               B => (:i,),
               C => (:i,),
           ),
           simd = true,  # optional
       ) == sum(A .* (A .> 0) .* B .* C)
true
```
"""
function ndreducible(pairs...)
    patterns = map(AccessPattern, pairs)
    # TODO: check consistency of patterns
    return NDReducible(patterns, plan(patterns...))
end

struct NDReducible{T <: NTuple{<:Any, AccessPattern}, P} <: Foldable
    # TODO: NDReducible <: Reducible
    patterns::T
    plan::P
end

struct NAY end  # not accessed yet

@inline function Transducers.__foldl__(rf, acc, coll::NDReducible)
    accessed = ntuple(_ -> NAY(), length(coll.patterns))
    result = _ndreduce(rf, acc, accessed, (), coll.patterns, coll.plan.indices...)
    @return_if_reduced result
    return complete(rf, result)
end

@inline function _ndreduce(rf, acc, accessed, indexmap0, patterns, i, indices...)
    for j in axisfor(patterns, i)
        indexmap = (indexmap0..., i => j)
        acc = @return_if_reduced _ndreduce(
            rf,
            acc,
            tryaccess(patterns, indexmap, accessed),
            indexmap,
            patterns,
            indices...,
        )
    end
    return acc
end

@inline function _ndreduce(rf, acc, accessed, indexmap0, patterns, i)
    @simd_if rf for j in axisfor(patterns, i)
        indexmap = (indexmap0..., i => j)
        acc = @next(rf, acc, tryaccess(patterns, indexmap, accessed))
    end
    return acc
end

axisfor(patterns, i) = unreduced(_axisfor(patterns, i) :: Reduced)

@inline _axisfor(patterns, i) =
    foldlargs(nothing, patterns...) do _, p::AccessPattern
        foldlargs(1, p.indices...) do n, j
            if j === i
                reduced(axes(p.indexable)[n])
            else
                n + 1
            end
        end
    end

@inline tryaccess(patterns, indexmap, accessed) =
    map(patterns, accessed) do (p::AccessPattern), value
        value isa NAY || return value
        idxvalue = map(p.indices) do (i::Index)
            foldlargs(NAY(), indexmap...) do found, (iname, ivalue)
                found isa NAY || return found
                (iname::Index) === i ? (ivalue::Int) : NAY()
            end
        end
        if idxvalue isa Tuple{Vararg{Int}}
            @inbounds p.indexable[idxvalue...]
        else
            NAY()
        end
    end


text:
module TotalVariation

using  DSP, SparseArrays, LinearAlgebra
export gstv, tv
#See ``Total Variation Denoising With Overlapping Group Sparsity'' by
# Ivan Selesnick and Po-Yu Chen (2013)

#Group Sparse Total Variation Denoising
function gstv(y::Vector{Float64}, k::Int, λ::Float64;
    show_cost::Bool=false, iter::Int=100)
    #Initialize Solution
    N  = length(y)
    x  = copy(y)

    #Differential of input
    b = diff(y)

    #Precalculate D D' where D is the first-order difference matrix
    DD::SparseMatrixCSC{Float64,Int} = spdiagm(-1=>-ones(N-2), 0=>2*ones(N-1), 1=>-ones(N-2))

    #Value To Prevent Singular Matrices
    epsilon = 1e-15

    #Convolution Mask - spreads D x over a larger area
    #This regularizes the problem with applying a gradient to a larger area.
    #at k=1 the normal total variational sparse solution (for D x) is found.
    h = ones(k)

    for i=1:iter
        u::Vector{Float64}              = diff(x)
        r::Vector{Float64}              = sqrt.(max.(epsilon, DSP.conv(u.^2,h)))
        Λ::Vector{Float64}              = DSP.conv(1 ./ r, h)[k:end-(k-1)]
        F::SparseMatrixCSC{Float64,Int} = sparse(Diagonal(1 ./ Λ))/λ + DD
        if(show_cost)
            #1/2||y-x||_2^2 + λΦ(Dx)
            #Where Φ(.) is the group sparse regularizer
            println("Cost at iter ",i," is ", 0.5*sum(abs2.(x.-y)) + λ*sum(r))
        end

        tmp::Vector{Float64} = F\b
        dfb::Vector{Float64} = diff(tmp)

        x[1]       = y[1]       + tmp[1]
        x[2:end-1] = y[2:end-1] + dfb[:]
        x[end]     = y[end]     - tmp[end]
    end

    return x
end

#Normal Total Variation Problem
function tv(y::Vector{Float64}, λ::Float64;
    show_cost::Bool=false, iter::Int=100)
    gstv(y, 1, λ, show_cost=show_cost, iter=iter)
end

# package code goes here
end # module

ext: jl


text:
module OpenSpielWrapper

using ..AlphaZero
using OpenSpiel

struct Spec{G} <: GI.AbstractGameSpec
  spiel_game :: G # loaded with OpenSpiel.load_game
  suppress_warnings :: Bool
  heuristic_value :: Function
  vectorize_state :: Function
  symmetries :: Function
  action_string :: Function
  parse_action :: Function
  read_state :: Function
end

mutable struct Env{G, S, P} <: GI.AbstractGameEnv
  spec :: Spec{G}
  state :: S
  curplayer :: P
end

#####
##### Default functions and constructors
#####

default_heurisic_value(env) = 0.0
# not handling terminal states, not flip board for two-player games
default_vectorize_state(spec, state) = convert(Array{Float32}, observation_tensor(state))
default_symmetries(spec, state) = GI.symmetries(spec, state)
default_action_string(spec, action) = GI.action_string(spec, action)
default_parse_action(spec, str::String) = GI.parse_action(spec, str)
default_read_state(spec) = GI.read_state(spec)

"""
    Spec(spiel_game; kwargs...) <: AbstractGameSpec

Wrap an OpenSpiel game object (`CxxWrap.StdLib.SharedPtrAllocated{OpenSpiel.Game}`)
into an AlphaZero game spec.

# Keyword arguments

The following optional functions from `GameInterface` are not present in
OpenSpiel.jl and can be provided as keyword arguments:

  - `heuristic_value`
  - `vectorize_state`
  - `symmetries`
  - `action_string`
  - `parse_action`
  - `read_state`

You can silence warnings by setting `suppress_warnings=true`.
"""

function Spec(
  spiel_game;
  suppress_warnings = false,
  heuristic_value = default_heurisic_value,
  vectorize_state = default_vectorize_state,
  symmetries = default_symmetries,
  action_string = default_action_string,
  parse_action = default_parse_action,
  read_state = default_read_state
)
  nplayers = num_players(spiel_game)
  @assert nplayers == 1 || nplayers == 2 "
  AlphaZero.jl only supports games with one or two players."
  if heuristic_value == default_heurisic_value && !suppress_warnings
    @warn "The `default_heuristic_value` function was not provided and so " *
      "algorithms such as MinMax may not work with this environment."
  end
  return Spec(
    spiel_game,
    suppress_warnings,
    heuristic_value,
    vectorize_state,
    symmetries,
    action_string,
    parse_action,
    read_state
  )
end

#####
##### GameInterface API
#####

function GI.init(spec::Spec)
  state = new_initial_state(spec.spiel_game)
  curplayer = current_player(state)
  return Env(spec, state, curplayer)
end

GI.spec(env::Env) = env.spec

# Queries on specifications

GI.two_players(spec::Spec) = num_players(spec.spiel_game) == 2

# TODO kinda ugly, may be game relevant
GI.actions(spec::Spec) = collect(0:num_distinct_actions(spec.spiel_game)-1)

GI.vectorize_state(spec::Spec, state) = spec.vectorize_state(spec, state)

# Operations on environments

GI.current_state(env::Env) = deepcopy(env.state)

function GI.set_state!(env::Env, state)
  env.state = deepcopy(state)
  env.curplayer = current_player(state)
end

GI.game_terminated(env::Env) = is_terminal(env.state)

function GI.white_playing(env::Env)
  wp = env.curplayer==0
  return is_terminal(env.state) && GI.two_players(env.spec) ? !wp : wp
end

# not all games implement legal_actions_mask
function GI.actions_mask(env::Env)
  try
    return convert(Array{Bool}, legal_actions_mask(env.state))
  catch e
    @error "AlphaZero.jl only works with OpenSpiel environments that implement " *
      "the `legal_actions_mask` function."
    rethrow(e)
  end
end

#TODO is_chance_node, is_simultaneous_node
function GI.play!(env::Env, action)
  apply_action(env.state, action)
  if !is_terminal(env.state)
    env.curplayer = current_player(env.state)
  end
end

GI.white_reward(env::Env) = returns(env.state)[1]

# available_actions dispatched
GI.available_actions(env::Env) = legal_actions(env.state)

GI.heuristic_value(env::Env) = env.spec.heuristic_value(env)

# Symmetries

GI.symmetries(spec, state) = spec.symmetries(spec.spiel_game, state)

# Interactive utilities

function GI.render(env::Env)
  show(env.state)
  print("\n\n") # GameInterface expects the rendering to end with an empty line
end

# action_string, parse action require Env instead of Spec for OpenSpiel functions
GI.action_string(env::Env, action) = action_to_string(env.state, action)
GI.action_string(spec::Spec, action) = spec.action_string(spec, action)

GI.parse_action(env::Env, str::String) = string_to_action(env.state, str)
GI.parse_action(spec::Spec, str::String) = spec.parse_action(spec, str)

GI.read_state(spec::Spec) = spec.read_state(spec.spiel_game)

end # module


text:
"""
for defining discrete Cartesian and Polar lattice grids used in Iris.
"""
module Lattices

export Lattice
export Bravais
export latticeindex

files = (
    "1D/Lattices1D.jl",
    "2D/Lattices2D.jl",
    # "3D/Lattices3D.jl"
)

using ..Points
using Formatting
using LinearAlgebra
using RecipesBase
using StaticArrays
using Statistics

import ..Points.Cartesian
import ..Points.Polar
import ..Points.Spherical
AbstractLatticeType = Points.AbstractCoordinateType
struct Bravais <: AbstractLatticeType end

struct Lattice{N,TYPE}
    constants::NTuple{N,Float64} # lattice spacings
    e::NTuple{N,Point{N,Cartesian}} # unit vectors
    origin::Point{N,Cartesian}

    function Lattice(::TYPE,
                primitives::NTuple{N,Point{N}},
                origin::Point{N},
                ) where {N,CE,CO} where TYPE<:AbstractLatticeType

        1≤N≤3 || throw(ErrorException("lattice only defined for dimensions ≤ 3, d=$N"))
        constants = _lattice_constants(primitives)
        vectors = _lattice_vectors(primitives)
        return new{N,TYPE}(constants, vectors, origin)
    end
end

_lattice_constants(primitives::Tuple) = map(p->norm(p),primitives)
_lattice_vectors(primitives::Tuple) = map(p->p/norm(p),primitives)

"""
    Bravais(primitives...; [origin]) -> bravaislattice
"""
Bravais(args...; kwargs...) = Lattice(Bravais(), args...; kwargs...)

"""
    Cartesian(dx,...; kwargs) -> clat

    Cartesian((dx,...); kwargs...) -> clat
"""
Cartesian(args...; kwargs...) = Lattice(Cartesian(), args...; kwargs...)

"""
    Polar(dx,...; [origin]) -> polarlat
"""
Polar(args...; kwargs...) = Lattice(Polar(),args...; kwargs...)

"""
    Spherical(dx,...; [origin]) -> sphlat
"""
Spherical(dx...; kwargs...) = Lattice(Spherical(), dx...; kwargs...)

Lattice(t::AbstractLatticeType, dx::Tuple; kwargs...) = Lattice(t,dx...; kwargs...)
Lattice(t::AbstractLatticeType, dx::Number, dy...; kwargs...) = Lattice(t,(float(dx),float.(dy)...); kwargs...)

"""
    Lattice(::Bravais, primitives...; [origin]) -> bravaislattice
"""
Lattice(::Bravais, primitive::Point, primitives...; kwargs... ) = Lattice(Bravais(),(primitive, primitives...); kwargs...)
function Lattice(::Bravais,
            primitives::NTuple{N,Point};
            origin = Point(ntuple(i->0.0,N)),
            ) where N
    constants = _lattice_constants(primitives)
    return Lattice(Bravais(), Cartesian.(primitives), _lattice_origin(origin,constants))
end
"""
    Lattice(::Union{Cartesian,Polar,Spherical}, dx...; [origin]) -> lattice
"""
function Lattice(type::Union{Cartesian,Polar,Spherical},
            constants::NTuple{N,Float64};
            origin = Point(ntuple(i->0.0,N)),
            kwargs...
            ) where N

    if haskey(kwargs,:angles)
        angles = kwargs[:angles]
    elseif haskey(kwargs,:angle)
        angles = kwargs[:angle]
    else
        angles = zeros(Float64,(N-1)N÷2)#SVector{M,Float64}(ntuple(i->0.0,(N-1)N÷2)),
    end

    M = length(angles)
    M==(N-1)N÷2 || throw(ArgumentError("incorrect number of angles for dimension, should be $((N-1)N÷2), but is $M"))
    return Lattice(type, _lattice_primitives(type,angles,constants), _lattice_origin(origin,constants))
end

_lattice_origin(origin::Point{N},constants::NTuple{N}) where N = origin
_lattice_origin(origin,constants::NTuple{N}) where N = Point(origin...)

foreach(include,files)

"""
    latticeindex(lattice, point) -> ind::Float64
"""
latticeindex(lat::Lattice{N,C},x::Real,y...) where {N,C} = latticeindex(lat,Point{C}(x,y...))

end # module


text:
# constants.jl
#
# author: Wade Shen
# swade@ll.mit.edu
# Copyright © 2014 Massachusetts Institute of Technology, Lincoln Laboratory
# version 0.1
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
import Base.getindex
export default_space, english_stoplist, news_email_header, unk_token, html_entity_table

const default_space = r"[\s\p{Zs}]+"

# PHP's table
const html_entity_table = Dict(
                               " " => " ",
                               "¡" => "¡",
                               "¢" => "¢",
                               "£" => "£",
                               "¤" => "¤",
                               "¥" => "¥",
                               "¦" => "¦",
                               "§" => "§",
                               "¨" => "¨",
                               "©" => "©",
                               "ª" => "ª",
                               "«" => "«",
                               "¬" => "¬",
                               "" => "",
                               "®" => "®",
                               "¯" => "¯",
                               "°" => "°",
                               "±" => "±",
                               "²" => "²",
                               "³" => "³",
                               "´" => "´",
                               "µ" => "µ",
                               "¶" => "¶",
                               "·" => "·",
                               "¸" => "¸",
                               "¹" => "¹",
                               "º" => "º",
                               "»" => "»",
                               "¼" => "¼",
                               "½" => "½",
                               "¾" => "¾",
                               "¿" => "¿",
                               "À" => "À",
                               "Á" => "Á",
                               "Â" => "Â",
                               "Ã" => "Ã",
                               "Ä" => "Ä",
                               "Å" => "Å",
                               "Æ" => "Æ",
                               "Ç" => "Ç",
                               "È" => "È",
                               "É" => "É",
                               "Ê" => "Ê",
                               "Ë" => "Ë",
                               "Ì" => "Ì",
                               "Í" => "Í",
                               "Î" => "Î",
                               "Ï" => "Ï",
                               "Ð" => "Ð",
                               "Ñ" => "Ñ",
                               "Ò" => "Ò",
                               "Ó" => "Ó",
                               "Ô" => "Ô",
                               "Õ" => "Õ",
                               "Ö" => "Ö",
                               "×" => "×",
                               "Ø" => "Ø",
                               "Ù" => "Ù",
                               "Ú" => "Ú",
                               "Û" => "Û",
                               "Ü" => "Ü",
                               "Ý" => "Ý",
                               "Þ" => "Þ",
                               "ß" => "ß",
                               "à" => "à",
                               "á" => "á",
                               "â" => "â",
                               "ã" => "ã",
                               "ä" => "ä",
                               "å" => "å",
                               "æ" => "æ",
                               "ç" => "ç",
                               "è" => "è",
                               "é" => "é",
                               "ê" => "ê",
                               "ë" => "ë",
                               "ì" => "ì",
                               "í" => "í",
                               "î" => "î",
                               "ï" => "ï",
                               "ð" => "ð",
                               "ñ" => "ñ",
                               "ò" => "ò",
                               "ó" => "ó",
                               "ô" => "ô",
                               "õ" => "õ",
                               "ö" => "ö",
                               "÷" => "÷",
                               "ø" => "ø",
                               "ù" => "ù",
                               "ú" => "ú",
                               "û" => "û",
                               "ü" => "ü",
                               "ý" => "ý",
                               "þ" => "þ",
                               "ÿ" => "ÿ",
                               "&" => "&",
                               """ => "\"",
                               "<" => "<",
                               ">" => ">",
                               )

# -------------------------------------------------------------------------------------------------------------------------
# From Lewis04's SMART stop list
# -------------------------------------------------------------------------------------------------------------------------
const english_stoplist = Set(AbstractString["a",
                                            "a's",
                                            "able",
                                            "about",
                                            "above",
                                            "according",
                                            "accordingly",
                                            "across",
                                            "actually",
                                            "after",
                                            "afterwards",
                                            "again",
                                            "against",
                                            "ain't",
                                            "all",
                                            "allow",
                                            "allows",
                                            "almost",
                                            "alone",
                                            "along",
                                            "already",
                                            "also",
                                            "although",
                                            "always",
                                            "am",
                                            "among",
                                            "amongst",
                                            "an",
                                            "and",
                                            "another",
                                            "any",
                                            "anybody",
                                            "anyhow",
                                            "anyone",
                                            "anything",
                                            "anyway",
                                            "anyways",
                                            "anywhere",
                                            "apart",
                                            "appear",
                                            "appreciate",
                                            "appropriate",
                                            "are",
                                            "aren't",
                                            "around",
                                            "as",
                                            "aside",
                                            "ask",
                                            "asking",
                                            "associated",
                                            "at",
                                            "available",
                                            "away",
                                            "awfully",
                                            "b",
                                            "be",
                                            "became",
                                            "because",
                                            "become",
                                            "becomes",
                                            "becoming",
                                            "been",
                                            "before",
                                            "beforehand",
                                            "behind",
                                            "being",
                                            "believe",
                                            "below",
                                            "beside",
                                            "besides",
                                            "best",
                                            "better",
                                            "between",
                                            "beyond",
                                            "both",
                                            "brief",
                                            "but",
                                            "by",
                                            "c",
                                            "c'mon",
                                            "c's",
                                            "came",
                                            "can",
                                            "can't",
                                            "cannot",
                                            "cant",
                                            "cause",
                                            "causes",
                                            "certain",
                                            "certainly",
                                            "changes",
                                            "clearly",
                                            "co",
                                            "com",
                                            "come",
                                            "comes",
                                            "concerning",
                                            "consequently",
                                            "consider",
                                            "considering",
                                            "contain",
                                            "containing",
                                            "contains",
                                            "corresponding",
                                            "could",
                                            "couldn't",
                                            "course",
                                            "currently",
                                            "d",
                                            "definitely",
                                            "described",
                                            "despite",
                                            "did",
                                            "didn't",
                                            "different",
                                            "do",
                                            "does",
                                            "doesn't",
                                            "doing",
                                            "don't",
                                            "done",
                                            "down",
                                            "downwards",
                                            "during",
                                            "e",
                                            "each",
                                            "edu",
                                            "eg",
                                            "eight",
                                            "either",
                                            "else",
                                            "elsewhere",
                                            "enough",
                                            "entirely",
                                            "especially",
                                            "et",
                                            "etc",
                                            "even",
                                            "ever",
                                            "every",
                                            "everybody",
                                            "everyone",
                                            "everything",
                                            "everywhere",
                                            "ex",
                                            "exactly",
                                            "example",
                                            "except",
                                            "f",
                                            "far",
                                            "few",
                                            "fifth",
                                            "first",
                                            "five",
                                            "followed",
                                            "following",
                                            "follows",
                                            "for",
                                            "former",
                                            "formerly",
                                            "forth",
                                            "four",
                                            "from",
                                            "further",
                                            "furthermore",
                                            "g",
                                            "get",
                                            "gets",
                                            "getting",
                                            "given",
                                            "gives",
                                            "go",
                                            "goes",
                                            "going",
                                            "gone",
                                            "got",
                                            "gotten",
                                            "greetings",
                                            "h",
                                            "had",
                                            "hadn't",
                                            "happens",
                                            "hardly",
                                            "has",
                                            "hasn't",
                                            "have",
                                            "haven't",
                                            "having",
                                            "he",
                                            "he's",
                                            "hello",
                                            "help",
                                            "hence",
                                            "her",
                                            "here",
                                            "here's",
                                            "hereafter",
                                            "hereby",
                                            "herein",
                                            "hereupon",
                                            "hers",
                                            "herself",
                                            "hi",
                                            "him",
                                            "himself",
                                            "his",
                                            "hither",
                                            "hopefully",
                                            "how",
                                            "howbeit",
                                            "however",
                                            "i",
                                            "i'd",
                                            "i'll",
                                            "i'm",
                                            "i've",
                                            "ie",
                                            "if",
                                            "ignored",
                                            "immediate",
                                            "in",
                                            "inasmuch",
                                            "inc",
                                            "indeed",
                                            "indicate",
                                            "indicated",
                                            "indicates",
                                            "inner",
                                            "insofar",
                                            "instead",
                                            "into",
                                            "inward",
                                            "is",
                                            "isn't",
                                            "it",
                                            "it'd",
                                            "it'll",
                                            "it's",
                                            "its",
                                            "itself",
                                            "j",
                                            "just",
                                            "k",
                                            "keep",
                                            "keeps",
                                            "kept",
                                            "know",
                                            "knows",
                                            "known",
                                            "l",
                                            "last",
                                            "lately",
                                            "later",
                                            "latter",
                                            "latterly",
                                            "least",
                                            "less",
                                            "lest",
                                            "let",
                                            "let's",
                                            "like",
                                            "liked",
                                            "likely",
                                            "little",
                                            "look",
                                            "looking",
                                            "looks",
                                            "ltd",
                                            "m",
                                            "mainly",
                                            "many",
                                            "may",
                                            "maybe",
                                            "me",
                                            "mean",
                                            "meanwhile",
                                            "merely",
                                            "might",
                                            "more",
                                            "moreover",
                                            "most",
                                            "mostly",
                                            "much",
                                            "must",
                                            "my",
                                            "myself",
                                            "n",
                                            "name",
                                            "namely",
                                            "nd",
                                            "near",
                                            "nearly",
                                            "necessary",
                                            "need",
                                            "needs",
                                            "neither",
                                            "never",
                                            "nevertheless",
                                            "new",
                                            "next",
                                            "nine",
                                            "no",
                                            "nobody",
                                            "non",
                                            "none",
                                            "noone",
                                            "nor",
                                            "normally",
                                            "not",
                                            "nothing",
                                            "novel",
                                            "now",
                                            "nowhere",
                                            "o",
                                            "obviously",
                                            "of",
                                            "off",
                                            "often",
                                            "oh",
                                            "ok",
                                            "okay",
                                            "old",
                                            "on",
                                            "once",
                                            "one",
                                            "ones",
                                            "only",
                                            "onto",
                                            "or",
                                            "other",
                                            "others",
                                            "otherwise",
                                            "ought",
                                            "our",
                                            "ours",
                                            "ourselves",
                                            "out",
                                            "outside",
                                            "over",
                                            "overall",
                                            "own",
                                            "p",
                                            "particular",
                                            "particularly",
                                            "per",
                                            "perhaps",
                                            "placed",
                                            "please",
                                            "plus",
                                            "possible",
                                            "presumably",
                                            "probably",
                                            "provides",
                                            "q",
                                            "que",
                                            "quite",
                                            "qv",
                                            "r",
                                            "rather",
                                            "rd",
                                            "re",
                                            "really",
                                            "reasonably",
                                            "regarding",
                                            "regardless",
                                            "regards",
                                            "relatively",
                                            "respectively",
                                            "right",
                                            "s",
                                            "said",
                                            "same",
                                            "saw",
                                            "say",
                                            "saying",
                                            "says",
                                            "second",
                                            "secondly",
                                            "see",
                                            "seeing",
                                            "seem",
                                            "seemed",
                                            "seeming",
                                            "seems",
                                            "seen",
                                            "self",
                                            "selves",
                                            "sensible",
                                            "sent",
                                            "serious",
                                            "seriously",
                                            "seven",
                                            "several",
                                            "shall",
                                            "she",
                                            "should",
                                            "shouldn't",
                                            "since",
                                            "six",
                                            "so",
                                            "some",
                                            "somebody",
                                            "somehow",
                                            "someone",
                                            "something",
                                            "sometime",
                                            "sometimes",
                                            "somewhat",
                                            "somewhere",
                                            "soon",
                                            "sorry",
                                            "specified",
                                            "specify",
                                            "specifying",
                                            "still",
                                            "sub",
                                            "such",
                                            "sup",
                                            "sure",
                                            "t",
                                            "t's",
                                            "take",
                                            "taken",
                                            "tell",
                                            "tends",
                                            "th",
                                            "than",
                                            "thank",
                                            "thanks",
                                            "thanx",
                                            "that",
                                            "that's",
                                            "thats",
                                            "the",
                                            "their",
                                            "theirs",
                                            "them",
                                            "themselves",
                                            "then",
                                            "thence",
                                            "there",
                                            "there's",
                                            "thereafter",
                                            "thereby",
                                            "therefore",
                                            "therein",
                                            "theres",
                                            "thereupon",
                                            "these",
                                            "they",
                                            "they'd",
                                            "they'll",
                                            "they're",
                                            "they've",
                                            "think",
                                            "third",
                                            "this",
                                            "thorough",
                                            "thoroughly",
                                            "those",
                                            "though",
                                            "three",
                                            "through",
                                            "throughout",
                                            "thru",
                                            "thus",
                                            "to",
                                            "together",
                                            "too",
                                            "took",
                                            "toward",
                                            "towards",
                                            "tried",
                                            "tries",
                                            "truly",
                                            "try",
                                            "trying",
                                            "twice",
                                            "two",
                                            "u",
                                            "un",
                                            "under",
                                            "unfortunately",
                                            "unless",
                                            "unlikely",
                                            "until",
                                            "unto",
                                            "up",
                                            "upon",
                                            "us",
                                            "use",
                                            "used",
                                            "useful",
                                            "uses",
                                            "using",
                                            "usually",
                                            "uucp",
                                            "v",
                                            "value",
                                            "various",
                                            "very",
                                            "via",
                                            "viz",
                                            "vs",
                                            "w",
                                            "want",
                                            "wants",
                                            "was",
                                            "wasn't",
                                            "way",
                                            "we",
                                            "we'd",
                                            "we'll",
                                            "we're",
                                            "we've",
                                            "welcome",
                                            "well",
                                            "went",
                                            "were",
                                            "weren't",
                                            "what",
                                            "what's",
                                            "whatever",
                                            "when",
                                            "whence",
                                            "whenever",
                                            "where",
                                            "where's",
                                            "whereafter",
                                            "whereas",
                                            "whereby",
                                            "wherein",
                                            "whereupon",
                                            "wherever",
                                            "whether",
                                            "which",
                                            "while",
                                            "whither",
                                            "who",
                                            "who's",
                                            "whoever",
                                            "whole",
                                            "whom",
                                            "whose",
                                            "why",
                                            "will",
                                            "willing",
                                            "wish",
                                            "with",
                                            "within",
                                            "without",
                                            "won't",
                                            "wonder",
                                            "would",
                                            "would",
                                            "wouldn't",
                                            "x",
                                            "y",
                                            "yes",
                                            "yet",
                                            "you",
                                            "you'd",
                                            "you'll",
                                            "you're",
                                            "you've",
                                            "your",
                                            "yours",
                                            "yourself",
                                            "yourselves",
                                            "z",
                                            "zero",
                                            "--number--",
                                            ])

const news_email_header = r"^((to|from|subject|organization|keywords|summary|lines|distribution|reply-to|nntp-posting-host|disclaimer|x-newsreader|article-i.d.):|(-|\|).*$|-+$|in article.*writes:$|in article.*$|.* writes:$)"

const unk_token = "--UNK--"


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ e44af6da-9045-11eb-15ab-ad339f35ca21
using ImageFiltering,FileIO,ImageShow,TestImages,Images

# ╔═╡ 002e63da-9046-11eb-20b5-4f1bf8ff7260
img = load("../hello-world.png");

# ╔═╡ 03a23ec6-9046-11eb-035c-ffb12919833f
imgg = imfilter(img, Kernel.gaussian(3));

# ╔═╡ 09e7aea6-9046-11eb-3d06-c3980898c7bb
imgl = imfilter(img, Kernel.Laplacian());

# ╔═╡ 3ffbb60e-9048-11eb-1b5e-57ef5a52c22d
# imgs = imfilter(img, Kernel.sobel())

# ╔═╡ 54499d7c-9048-11eb-0b1b-13f6373e402c
img_coins = load("../1.jpg")

# ╔═╡ 9d38abe4-905d-11eb-1ba2-958ed4d9e9af
img_gray = Gray.(img_coins)

# ╔═╡ da50f694-905d-11eb-29dd-d160c78bc260
kernel = [1 0 -1; 2 0 -2;1 0 -1];

# ╔═╡ e082aeb8-905d-11eb-047c-5bb3725816ab
sobel_x = imfilter(img_gray, kernel);

# ╔═╡ e42883f6-905d-11eb-07bf-a9f6acecf1a1
grad = imfilter(sobel_x, kernel')

# ╔═╡ acceb0e8-905e-11eb-027a-a1c08f28641c

# ╔═╡ 02069d1a-905e-11eb-3471-9f3545ac23f4
grad .^ 4

# ╔═╡ 10460d34-905e-11eb-15ab-57d2dd4b66ef
begin
	kernel_1 = [1 2 1]';
kernel_2 = [1 0 -1];
grad_sep = imfilter(img_gray, (kernel_1, kernel_2, kernel_2', kernel_1'))
end

# ╔═╡ Cell order:
# ╠═e44af6da-9045-11eb-15ab-ad339f35ca21
# ╠═002e63da-9046-11eb-20b5-4f1bf8ff7260
# ╠═03a23ec6-9046-11eb-035c-ffb12919833f
# ╠═09e7aea6-9046-11eb-3d06-c3980898c7bb
# ╠═3ffbb60e-9048-11eb-1b5e-57ef5a52c22d
# ╠═54499d7c-9048-11eb-0b1b-13f6373e402c
# ╠═9d38abe4-905d-11eb-1ba2-958ed4d9e9af
# ╠═da50f694-905d-11eb-29dd-d160c78bc260
# ╠═e082aeb8-905d-11eb-047c-5bb3725816ab
# ╠═e42883f6-905d-11eb-07bf-a9f6acecf1a1
# ╠═acceb0e8-905e-11eb-027a-a1c08f28641c
# ╠═02069d1a-905e-11eb-3471-9f3545ac23f4
# ╠═10460d34-905e-11eb-15ab-57d2dd4b66ef

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
module RestClient

using HTTP

macro GET_str(data)
    return prepare_request(:GET, data)
end

macro POST_str(data)
    return prepare_request(:POST, data)
end

macro PUT_str(data)
    return prepare_request(:PUT, data)
end

macro DELETE_str(data)
    return prepare_request(:DELETE, data)
end

macro OPTIONS_str(data)
    return prepare_request(:OPTIONS, data)
end

function make_request(method, data)
    #=
    transforms " Key1 : Value 233 44 " to "Key1" => "Value 233 44"
    (split by : and trims)
    =#
    parse_header(header) = Pair(strip.(split(header, ":"))...)

    lines = split(data, "\n")
    url = lines[1]
    empty_index = findfirst(isequal(""), lines)
    body = HTTP.nobody
    if empty_index === nothing
        headers = []
    else
        headers = parse_header.(lines[2:empty_index - 1])
        if empty_index < length(lines)
            body = join(lines[empty_index + 1:end], "\n")
        end
    end

    return HTTP.request(method, url, headers, body)
end

function prepare_request(method, data)
    data = replace(data, "\"" => "\\\"")
    m = Base.string(method)
    return esc(
            :(RestClient.make_request(Symbol($m), $(Meta.parse("\"$data\""))))
        )
end

export @GET_str, @POST_str,  @PUT_str, @DELETE_str, @OPTIONS_str

end


text:
@testset "Testing model calibration" begin

    path = joinpath(Dolo.pkg_path, "examples", "models")
    model = Dolo.yaml_import(joinpath(path, "rbc_iid.yaml"))

    calib_0 = copy( model.calibration.flat )

    # check one value
    beta = model.calibration.flat[:beta]
    delta = model.calibration.flat[:delta]
    rk = model.calibration.flat[:rk]
    k0 = model.calibration[:states][2]
    @test beta==0.99
    @test rk==1/beta-1+delta
    @test abs(k0-9.35498)<1e-5

    # check calibrations are still correct after modifications
    Dolo.set_calibration!(model, :beta, 0.98)
    beta = model.calibration.flat[:beta]
    delta = model.calibration.flat[:delta]
    rk = model.calibration.flat[:rk]
    k0 = model.calibration[:states][2]
    @test beta==0.98
    @test rk==1/beta-1+delta
    @test abs(k0-6.37024)<1e-4

    # can we change many parameters at the same time?
    Dolo.set_calibration!(model; beta=0.97, i=:(delta*k-0.01))
    rk = model.calibration.flat[:rk]
    beta = model.calibration.flat[:beta]
    delta = model.calibration.flat[:delta]
    k0 = model.calibration[:states][2]
    i0 = model.calibration[:controls][2]
    @test beta==0.97
    @test rk==1/beta-1+delta
    @test i0==delta*k0-0.01

    # put everything back to original values
    set_calibration!(model; beta=0.99, i=:(delta*k))
    # check one value
    beta = model.calibration.flat[:beta]
    delta = model.calibration.flat[:delta]
    rk = model.calibration.flat[:rk]
    k0 = model.calibration[:states][2]
    i0 = model.calibration[:controls][2]
    @test beta==0.99
    @test rk==1/beta-1+delta
    @test abs(k0-9.35498)<1e-5
    @test i0==delta*k0

end


text:
import JSON

function DashBase.to_dash(p::PlotlyJS.SyncPlot)
    data = JSON.lower(p.plot)
    pop!(data, :config, nothing)
    return data
end


text:
module NKLandscapes

include("aliases.jl")

include("landscape.jl")
include("genotype.jl")
include("population.jl")
include("metapopulation.jl")
include("neighbors.jl")
include("walks.jl")
include("enumeration.jl")

include("selection/moran.jl")
include("selection/proportional.jl")
include("selection/tournament.jl")
include("selection/elite.jl")

include("mutation/bitwise.jl")
include("mutation/bitstring.jl")

include("migration.jl")

end


text:
using PackageCompiler
create_sysimage(:GalewskyShallowWaterThetaMethod,
  sysimage_path=joinpath(@__DIR__,"..","GalewskyShallowWaterThetaMethod.so"),
  precompile_execution_file=joinpath(@__DIR__,"warmup.jl"))

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences


text:
@info "Running tests"
using SpectralDistances# Distributions
using Test, LinearAlgebra, Statistics, Random, ControlSystems, InteractiveUtils, SparseArrays # For subtypes
using DSP, Distances, DoubleFloats
import GLPK, Convex, JuMP, SCS

using SpectralDistances: ngradient, nhessian, njacobian, polyconv, hproots, rev

@testset "SpectralDistances.jl" begin

    @testset "complete distmat" begin
        @info "Testing complete distmat"

        P = randn(2,40)
        D = pairwise(SqEuclidean(), P, dims=2)
        W = rand(size(D)...) .> 0.3 # Create a random mask
        W = (W + W') .> 0           # It makes sense for the mask to be symmetric
        W[diagind(W)] .= true
        D0 = W .* D                 # Remove missing entries

        D2,S = complete_distmat(D0, W)

        @test (norm(D-D2)/norm(D)) < 1e-5
        @test (norm(W .* (D-D2))/norm(D)) < 1e-5

        X  = Diagonal(sqrt.(S.S))*S.Vt

        # Verify that reconstructed `X` is correct up to rotation and translation
        A = [X[1:2,:]' ones(size(D,1))]
        P2 = (A*(A \ P'))'
        @test norm(P-P2)/norm(P) < 1e-3

    end

    @testset "Convolutional" begin
        @info "Testing Convolutional"
        include("test_convolutional.jl")
    end

    @testset "DTW" begin
        @info "Testing DTW"
        include("test_dtw.jl")
    end

    @testset "Model creation" begin
        @info "Testing Model creation"

        @testset "AR" begin
            @info "Testing AR"

            a = randn(3)
            m = AR(a)
            @test_nowarn show(m)
            @test_nowarn show([m, m])
            @test m.a == a
            @test length(m.p) == 2

            m = AR(Continuous(), a)
            @test m.ac == a
            @test length(m.pc) == 2

            r = ContinuousRoots(hproots(rev(a)))
            m = AR(r)
            @test m.pc == r

            r = DiscreteRoots(hproots(rev(a)))
            m = AR(r)
            @test m.p == r
        end

        @testset "ARMA" begin
            @info "Testing ARMA"

            a = randn(3)
            b = [1, 0.1]
            m = ARMA(b,a)
            @test s1(m.a) == s1(a)
            @test length(m.p) == 2
            @test m.b == b
            @test length(m.z) == 1

            m = ARMA(Continuous(), b,a)
            @test s1(m.ac) == s1(a)
            @test length(m.pc) == 2
            @test m.bc == b
            @test length(m.zc) == 1

            r = ContinuousRoots(hproots(rev(a)))
            z = ContinuousRoots(hproots(rev(b)))
            m = ARMA(z,r)
            @test m.pc == r
            @test m.zc == z

            r = DiscreteRoots(hproots(rev(a)))
            z = DiscreteRoots(hproots(rev(b)))
            m = ARMA(z,r)
            @test m.p == r
            @test m.z == z
        end

    end

    @testset "time" begin
        @info "Testing time"
        include("test_time.jl")
    end

    @testset "plot" begin
        @info "Testing plot"
        include("test_plotting.jl")
    end

    @testset "Solvers" begin
        @info "Testing Solvers"

        C = Float64.(Matrix(I(2)))
        a = [1., 0]
        b = [0, 1.]
        Γ,u,v = @inferred sinkhorn(C,a,b,β=0.01)
        @test Γ ≈ [0 1; 0 0]

        a = [0.5, 0.5]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = @inferred sinkhorn(C,a,b,β=0.1)
        @test Γ ≈ [0 0.5; 0 0.5]

        a = [1., 0]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = @inferred sinkhorn_log(C,a,b,β=0.01)
        @test Γ ≈ [0 1; 0 0]

        a = [0.5, 0.5]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = @inferred sinkhorn_log(C,a,b,β=0.01)
        @test Γ ≈ [0 0.5; 0 0.5]

        a = [0.5, 0.5]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = @inferred sinkhorn_log!(C,a,b,β=0.01)
        @test Γ ≈ [0 0.5; 0 0.5]

        a = [1., 0]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = @inferred IPOT(C,a,b,β=1)
        @test Γ ≈ [0 1; 0 0]

        a = [0.5, 0.5]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = @inferred IPOT(C,a,b)
        @test Γ ≈ [0 0.5; 0 0.5]

        a = [1., 0]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = @inferred sinkhorn_unbalanced(C,a,b,Balanced(),β=0.01)
        @test Γ ≈ [0 1; 0 0]

        a = [0.5, 0.5]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = @inferred sinkhorn_unbalanced(C,a,b,Balanced())
        @test Γ ≈ [0 0.5; 0 0.5]

        a = [1., 0]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = ot_jump(C,a,b)
        @test Γ ≈ [0 1; 0 0]

        a = [0.5, 0.5]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = ot_jump(C,a,b)
        @test Γ ≈ [0 0.5; 0 0.5]

        a = [1., 0]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = ot_convex(C,a,b)
        @test Γ ≈ [0 1; 0 0]

        a = [0.5, 0.5]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = ot_convex(C,a,b)
        @test Γ ≈ [0 0.5; 0 0.5]

        a = [0.5, 0.5]
        b = [0, 1.]
        C = Float64.(Matrix(I(2)))
        Γ,u,v = ot_convex(C,a,b)
        @test Γ ≈ [0 0.5; 0 0.5]

        models = examplemodels(2)
        Γ = transport_plan(OptimalTransportRootDistance(domain=Continuous(), weight=unitweight), models...; tol=1e-4)
        @test sum(Γ) ≈ 1 atol=0.001

    end

    @testset "gradients" begin
        @info "Testing gradients"
        include("test_diff.jl")
    end

    @testset "Energy" begin
        @info "testing Energy"
        for σ² = [0.1, 1., 2., 3]
            m = AR(ContinuousRoots([-1.]), σ²)
            e = @inferred spectralenergy(Continuous(),m)
            @test e ≈ σ²
            x = sqrt(σ²)randn(10000)
            m = LS(na=10)(x)
            @test spectralenergy(Continuous(),m) ≈ σ² atol=0.15
            m = TLS(na=10)(x)
            @test spectralenergy(Continuous(),m) ≈ σ² atol=0.15
        end
        # y = filt(numvec(Discrete(),m), denvec(Discrete(),m), x)
        # @test var(y) ≈ var(x)
    end

    @testset "Model estimation" begin
        y = sin.(0:0.1:100)
        fm = @inferred LS(na=2, λ=0)
        m = @inferred fm(y)
        @test fm(m) == m
        @test imag.(m.pc) ≈ [-0.1, 0.1] rtol=1e-4
        mc = @inferred SpectralDistances.change_precision(Float32, m)
        @test m ≈ mc
        @test eltype(mc.a) == Float32

        fm = @inferred TLS(na=2)
        m  = @inferred fm(y)
        @test imag.(m.pc) ≈ [-0.1, 0.1] rtol=1e-4

        fm = IRLS(na=2)
        m = fm(y)
        @test_broken imag.(m.pc) ≈ [-0.1, 0.1] rtol=1e-4

        y = sin.(0:0.1:1000) .+ 0.01 .*randn.()
        fm = PLR(na=2, nc=1)
        m = fm(y)
        @test imag.(log(m.p)) ≈ [-0.1, 0.1] rtol=1e-3

        y = sin.(0:0.1:100) .+ sin.(2 .* (0:0.1:100) .+ 0.3)
        fm = @inferred LS(na=4, λ=0)
        m = @inferred fm(y)
        @test imag.(m.pc) ≈ [-0.2, -0.1, 0.1, 0.2] rtol=1e-1

        fm = TLS(na=4)
        m = fm(y)
        @test spectralenergy(Continuous(), m) ≈ var(y) rtol=1e-3
        @test imag.(m.pc) ≈ [-0.2, -0.1, 0.1, 0.2] rtol=1e-4

        fm = IRLS(na=4)
        m = fm(y)
        @test @inferred(spectralenergy(Continuous(), m)) ≈ var(y) rtol=1e-3
        @test_broken imag.(m.pc) ≈ [-0.2, -0.1, 0.1, 0.2] rtol=1e-4

        y = sin.(0:0.1:1000) .+ sin.(2 .* (0:0.1:1000)) .+ 0.001 .*randn.()
        fm = @inferred PLR(na=4, nc=1, λ=0.0)
        m = fm(y)
        @test imag.(m.pc) ≈ [-0.2, -0.1, 0.1, 0.2] rtol=0.6

    end

@testset "modeldistance" begin
    t = 1:300
    ϵ = 1e-7
    fitmethod = LS
    for fitmethod in [LS, TLS]
        @info "Testing fitmethod $(string(fitmethod))"
        ls_loss = ModelDistance(fitmethod(na=2), EuclideanRootDistance(domain=Discrete()))
        @test @inferred(ls_loss(sin.(t), sin.(t))) < ϵ
        @test ls_loss(sin.(t), -sin.(t)) < ϵ # phase flip invariance
        @test ls_loss(sin.(t), sin.(t .+ 1)) < ϵ # phase shift invariance
        @test ls_loss(sin.(t), sin.(t .+ 0.1)) < ϵ # phase shift invariance
        @test ls_loss(10sin.(t), sin.(t .+ 1)) < ϵ # amplitude invariance
        @test ls_loss(sin.(t), sin.(1.1 .* t)) < 0.2 # small frequency shifts gives small errors
        @test ls_loss(sin.(0.1t), sin.(1.1 .* 0.1t)) < 0.1 # small frequency shifts gives small errors

        @test ls_loss(sin.(t), sin.(1.1 .* t)) ≈ ls_loss(sin.(0.1t), sin.(0.2t)) rtol=1e-2  # frequency shifts of relative size should result in the same error, probably only true for p=1
        ls_loss = ModelDistance(fitmethod(na=10), OptimalTransportRootDistance(domain=Discrete()))
        @test @inferred(ls_loss(filtfilt(ones(10),[10], randn(1000)), filtfilt(ones(10),[10], randn(1000)))) < 0.1 # Filtered through same filter, this test is very non-robust for TLS
        @test ls_loss(filtfilt(ones(10),[10], randn(1000)), filtfilt(ones(10),[10], randn(1000))) < ls_loss(filtfilt(ones(4),[4], randn(1000)), filtfilt(ones(10),[10], randn(1000))) # Filtered through different filters, this test is not robust
    end
    @testset "PLR" begin
        fitmethod = PLR
        @info "Testing fitmethod $(string(fitmethod))"
        t = 1:1000 .+ 0.01 .* randn.()
        ϵ = 0.01
        @info "Creating model distance"
        ls_loss = ModelDistance(fitmethod(na=2,nc=1), EuclideanRootDistance(domain=Discrete()))
        @info "Calc loss"
        @test ls_loss(sin.(t), sin.(t)) < ϵ
        @test ls_loss(sin.(t), -sin.(t)) < ϵ # phase flip invariance
        @test ls_loss(sin.(t), sin.(t .+ 1)) < ϵ # phase shift invariance
        @test ls_loss(sin.(t), sin.(t .+ 0.1)) < ϵ # phase shift invariance
        @info "Halfway"
        @test ls_loss(10sin.(t), sin.(t .+ 1)) < ϵ # amplitude invariance
        @test ls_loss(sin.(t), sin.(1.1 .* t)) < 0.2 # small frequency shifts gives small errors
        @test ls_loss(sin.(0.1t), sin.(1.1 .* 0.1t)) < 0.1 # small frequency shifts gives small errors

        @test ls_loss(sin.(t), sin.(1.1 .* t)) ≈ ls_loss(sin.(0.1t), sin.(0.2t)) rtol=1e-3  # frequency shifts of relative size should result in the same error, probably only true for p=1
        @info "Almost there"
        # ls_loss = ModelDistance(fitmethod(na=10,nc=2), OptimalTransportRootDistance(domain=Discrete())) # This test segfaulted sometimes, couldn't figure out the problem
        # @test ls_loss(randn(1000), randn(1000)) > 0.05
        # @info "Done"
        # @test ls_loss(filtfilt(ones(10),[10], randn(1000)), filtfilt(ones(10),[10], randn(1000))) < 1 # Filtered through same filter, this test is very non-robust for TLS
        # @test ls_loss(filtfilt(ones(10),[10], randn(1000)), filtfilt(ones(10),[10], randn(1000))) < ls_loss(filtfilt(ones(4),[4], randn(1000)), filtfilt(ones(10),[10], randn(1000))) # Filtered through different filters, this test is not robust
    end
end

@testset "discrete_grid_transportplan" begin
    x = [1.,0,0]
    y = [0,0.5,0.5]

    g = @inferred discrete_grid_transportplan(x,y)
    @test sum(g,dims=1)[:] == y
    @test sum(g,dims=2)[:] == x

    # test robustness for long vectors
    x = s1(rand(1000))
    y = s1(rand(1000))
    D = [abs2(abs(i-j)/1000) for i in eachindex(x), j in eachindex(y)]
    g = discrete_grid_transportplan(x,y)
    @test sum(g,dims=1)[:] ≈ y
    @test sum(g,dims=2)[:] ≈ x
    @test dot(g, D) ≈ discrete_grid_transportcost(x,y) rtol=0.01

    g = discrete_grid_transportplan(y,x)
    @test sum(g,dims=1)[:] ≈ x
    @test sum(g,dims=2)[:] ≈ y
    @test dot(g, D) ≈ discrete_grid_transportcost(y,x) rtol=0.01

    # test exception for unequal masses
    x = s1(rand(Float32,1000))
    y = rand(Float32,1000)
    @test_throws ErrorException discrete_grid_transportplan(x,y)
    @test_throws ErrorException discrete_grid_transportplan(y,x)

    @test_throws ErrorException discrete_grid_transportcost(x,y)
    @test_throws ErrorException discrete_grid_transportcost(y,x)

end

# p = randn(20)
# @btime Zygote.gradient(p) do p
#     r = riroots(p)
#     sum([r[1]; r[2]])
# end
#
#
# jacobian(fd, p->begin
#     r = riroots(p)
#     sum([r[1]; r[2]])
# end, p)

@testset "eigenvalue manipulations" begin
    @test @inferred(SpectralDistances.reflectd(2)) ≈ 0.5 + 0.0im
    @test @inferred(SpectralDistances.reflectd(complex(0,2))) ≈ 0 + 0.5im
    @test @inferred(SpectralDistances.reflectd(complex(0,-2))) ≈ 0 - 0.5im

    e = @inferred SpectralDistances.hproots(randn(7))
    ed = @inferred DiscreteRoots(e)
    @test real(e) == real.(e)
    @test real(ed) == real.(ed)
    @test issorted(ed.r, by=angle)
    @test all(<(1) ∘ abs, ed)
    ec = log(ed)
    @test issorted(ec, by=SpectralDistances.imageigsortby)
    @test all(<(0) ∘ real, ec)
    @test @inferred(domain_transform(Continuous(), ed)) isa ContinuousRoots
    @test domain_transform(Continuous(), ed) ≈ SpectralDistances.eigsort(SpectralDistances.reflectc.(log.(ed))) ≈ ContinuousRoots(ed)
    @test @inferred(domain_transform(Discrete(), ed)) == ed
    @test domain(ed) isa Discrete
    @test domain(ContinuousRoots(ed)) isa Continuous
    @test log(ed) isa ContinuousRoots

    @test all(<(1) ∘ abs, reflect(ed))
    @test all(<(0) ∘ real, reflect(ec))

    E = embedding(ec)
    @test length(E) == 12
    @test E == embedding(Vector, ec)
    E = embedding(Matrix, ec)
    @test size(E) == (6,2)
    E = embedding(Matrix, ec, false)
    @test size(E) == (3,2)

    @test_logs (:error, r"real poles") embedding(ContinuousRoots([0,1]), false)

    Random.seed!(0)
    @test SpectralDistances.determine_domain(0.1randn(10)) isa Discrete
    @test SpectralDistances.determine_domain(randn(10).-4) isa Continuous
    @test_throws Exception SpectralDistances.determine_domain(0.1randn(10).-0.3)

    @testset "weight functions" begin
        @info "Testing weight functions"
        r = @inferred ContinuousRoots(randn(ComplexF64, 4))
        m = @inferred AR(r, 2.0)
        @test sum(unitweight(r)) == 1
        @test sum(unitweight(m)) == 1
        @test length(unitweight(r)) == 4
        @test length(unitweight(m)) == 4

        @test m.b^2*sum(residueweight(r)) ≈ sum(residueweight(m)) # Calling residueweight with a model should mutiply energy by b²
        @test length(residueweight(r)) == 4
        @test length(residueweight(m)) == 4
    end

end

@testset "ControlSystems interoperability" begin
    @info "Testing ControlSystems interoperability"
    m = @inferred AR(ContinuousRoots([-1]))
    g = tf(1,[1.,1])
    @test tf(m) == g
    @test_broken tf(m,1) == c2d(g,1)
    @test m*g == g*g
    @test ARMA(g).ac == [1,1]
    @test ARMA(g).bc == [1]
    @test @inferred(denvec(Continuous(), m)) == denvec(g)[1]
    @test numvec(Continuous(), m) == numvec(g)[1]
    @test @inferred(pole(Continuous(), m)) == pole(g)
    @test all(ControlSystems.bode(m) .≈ bode(g))
    @test all(ControlSystems.nyquist(m) .≈ nyquist(g))
    @test ControlSystems.freqresp(m, exp10.(LinRange(-1, 1, 10))) ≈ freqresp(g, exp10.(LinRange(-1, 1, 10)))
    @test all(ControlSystems.step(m, 10) .≈ step(g, 10))

    bodeplot(m)
    pzmap(m)
    nyquistplot(m)

end

@testset "polynomial acrobatics" begin
    a = randn(5)
    b = randn(5)
    @test @inferred(polyconv(a,b)) ≈ DSP.conv(a,b)

    a = randn(5)
    b = randn(10)
    @test polyconv(a,b) ≈ DSP.conv(a,b)
    @test polyconv(b,a) ≈ DSP.conv(b,a)

    a[1] = 1
    @test @inferred(roots2poly(roots(reverse(a)))) ≈ a

end

@testset "preprocess roots with residue weight" begin
    rd  = @inferred EuclideanRootDistance(domain=Continuous(), weight=residueweight)
    m = @inferred AR([1., -0.1])
    @test @inferred(rd(m,m)) < eps()
    @test @inferred(SpectralDistances.preprocess_roots(rd, m))[] ≈ -2.3025850929940455
end

@testset "residues and roots" begin

    # @test inv(1:3) == 1:3
    # @test inv([1,1,2]) == [1,1,3]
    # @test inv([1,2,2]) == [1,1,2]

    a = randn(5); a[1]=1;a
    r = roots(reverse(a))
    @test roots2poly(r) ≈ a

    b,a = 1,0.1randn(5)
    a[1] = 1
    r   = @inferred SpectralDistances.hproots(reverse(a))
    G   = tf(b,a)
    w   = im
    F = evalfr(G,w)[]
    a2 = roots2poly(r)
    @assert a ≈ a2
    roots(reverse(a2))
    @test sum(r) do r
        powers = length(a)-1:-1:1
        ap     = powers.*a[1:end-1]
        G      = tf(b,ap)
        evalfr(G,r)[] *1/(w - r)
    end ≈ F

    res = @inferred residues(a,1)
    @test sum(eachindex(r)) do i
        res[i]/(w-r[i])
    end ≈ F

    @test abs2.(res) ≈ SpectralDistances.abs2residues!(zeros(length(a)-1),a,1)

    @test prod(eachindex(r)) do i
        1/(w-r[i])
    end ≈ F

    n = 4
    a1 = randn(n+1); a1[1] = 1
    r1 = SpectralDistances.reflectc.(hproots(reverse(a1)))
    r1 = complex.(0.01real.(r1), imag.(r1))
    r1 = @inferred SpectralDistances.normalize_energy(ContinuousRoots(r1))

    a2 = randn(n+1); a2[1] = 1
    r2 = SpectralDistances.reflectc.(hproots(reverse(a2)))
    r2 = complex.(0.01real.(r2), imag.(r2))
    r2o = @inferred SpectralDistances.normalize_energy(ContinuousRoots(r2))

    m1,m2 = AR(a1), AR(a2)

    dist2 = @inferred RationalOptimalTransportDistance(domain=Continuous(), p=2, interval=(-20.,20.))
    f    = w -> evalfr(SpectralDistances.domain(dist2), SpectralDistances.magnitude(dist2), w, m1)
    @test SpectralDistances.c∫(f,dist2.interval...)[end] ≈ spectralenergy(Continuous(), m1) rtol=1e-3
    f    = w -> evalfr(SpectralDistances.domain(dist2), SpectralDistances.magnitude(dist2), w, m2)
    @test SpectralDistances.c∫(f,dist2.interval...)[end] ≈ spectralenergy(Continuous(), m2) rtol=1e-3

    @test spectralenergy(Continuous(), AR(ContinuousRoots([-1.]))) ≈ spectralenergy(tf(1,[1,1])) rtol=1e-3

    m1,m2 = ARMA([1], a1), ARMA([1], a2)

    f    = w -> evalfr(SpectralDistances.domain(dist2), SpectralDistances.magnitude(dist2), w, m1)
    @test_broken SpectralDistances.c∫(f,dist2.interval...)[end] ≈ spectralenergy(Continuous(), m1) rtol=1e-3 # The error is in the integration, spectgralenergy produces same result as for AR above
    f    = w -> evalfr(SpectralDistances.domain(dist2), SpectralDistances.magnitude(dist2), w, m2)
    @test_broken SpectralDistances.c∫(f,dist2.interval...)[end] ≈ spectralenergy(Continuous(), m2) rtol=1e-3

end

@testset "d(m,m)" begin
    a = roots2poly([0.9 + 0.1im, 0.9 - 0.1im])
    m = AR(a)
    for D in [  subtypes(SpectralDistances.AbstractDistance);
                subtypes(SpectralDistances.AbstractRootDistance);
                subtypes(SpectralDistances.AbstractCoefficientDistance)]
        (!isempty(methods(D)) && (:domain ∈ fieldnames(D))) || continue
        d = @inferred D(domain=Continuous())
        println(D)
        @test @inferred(d(m,m)) < eps() + 0.001*(d isa OptimalTransportRootDistance)
        d isa Union{RationalOptimalTransportDistance, RationalCramerDistance} && continue
        d = D(domain=Discrete())
        println(D)
        @test d(m,m) < eps() + 0.001*(d isa OptimalTransportRootDistance)
    end
end

@testset "d(m,m̃)" begin
    a1 = [1,-0.1,0.8]
    m1 = @inferred AR(a1)
    a2 = [1,-0.1,0.801]
    m2 = AR(a2)
    for D in [  subtypes(SpectralDistances.AbstractDistance);
                subtypes(SpectralDistances.AbstractRootDistance);
                subtypes(SpectralDistances.AbstractCoefficientDistance)]
        (!isempty(methods(D)) && (:domain ∈ fieldnames(D))) || continue
        d = D(domain=Continuous())
        println(D)
        # @show d(m1,m2)
        @test @inferred(d(m1,m2)) > 1e-10
        d isa Union{RationalOptimalTransportDistance, RationalCramerDistance} && continue
        d = D(domain=Discrete())
        println(D)
        # @show d(m1,m2)
        @test d(m1,m2) > 1e-10
    end
end

@testset "distmat" begin
    e = complex.(randn(3), randn(3))
    D = @inferred SpectralDistances.distmat(SqEuclidean(), e)
    @test issymmetric(D)
    @test tr(D) == 0
    @test @inferred(SpectralDistances.distmat_euclidean(e,e)) ≈ D

    A = randn(10,10)
    A[diagind(A)] .= rand(10);
    A = A + A';
    B = SparseMatrixCSC(A)

    symmetrize!(A)
    symmetrize!(B)
    @test A == B
    @test issymmetric(A)
end

@testset "Welch" begin
    x1 = @inferred SpectralDistances.bp_filter(randn(3000), (0.01,0.1))
    x2 = SpectralDistances.bp_filter(randn(3000), (0.01,0.12))
    x3 = SpectralDistances.bp_filter(randn(3000), (0.01,0.3))
    dist = @inferred WelchOptimalTransportDistance(p=1)
    @test @inferred(dist(x1,x2)) < dist(x1,x3)
    dist = WelchOptimalTransportDistance(p=2)
    @test dist(x1,x2) < dist(x1,x3)

    dist = @inferred WelchLPDistance(p=1)
    @test @inferred(dist(x1,x2)) < dist(x1,x3)
    dist = WelchLPDistance(p=2)
    @test dist(x1,x2) < dist(x1,x3)

    dist = @inferred EnergyDistance()
    @test @inferred(dist(x1,x2)) < dist(x1,x3)
end

@testset "Bures" begin
    fm = LS(na=4)
    x1 = SpectralDistances.bp_filter(randn(3000), (0.01,0.1))  |> fm
    x2 = SpectralDistances.bp_filter(randn(3000), (0.01,0.12)) |> fm
    x3 = SpectralDistances.bp_filter(randn(3000), (0.01,0.3))  |> fm
    dist = BuresDistance()
    @test dist(x1,x2) < dist(x1,x3)
end

@testset "KernelWassersteinRootDistance" begin
    fm = LS(na=4)
    x1 = SpectralDistances.bp_filter(randn(3000), (0.01,0.1))  |> fm
    x2 = SpectralDistances.bp_filter(randn(3000), (0.01,0.12)) |> fm
    x3 = SpectralDistances.bp_filter(randn(3000), (0.01,0.3))  |> fm
    dist = @inferred KernelWassersteinRootDistance(domain=Continuous())
    @test @inferred(dist(x1,x2)) < dist(x1,x3)
end

@testset "ClosedForm" begin
    x1 = SpectralDistances.bp_filter(randn(3000), (0.01,0.1))
    x2 = SpectralDistances.bp_filter(randn(3000), (0.01,0.12))
    x3 = SpectralDistances.bp_filter(randn(3000), (0.01,0.3))
    fm = LS(na=4)
    dist = ModelDistance(fm, RationalOptimalTransportDistance(domain=Continuous(), p=1, interval=(-15., 15)))
    @test dist(x1,x2) < dist(x1,x3)
    dist = RationalOptimalTransportDistance(domain=Continuous(), p=1, interval=(0., 15.))
    @test dist(fm(x1),welch_pgram(x2)) < dist(fm(x1),welch_pgram(x3))
    @test dist(fm(x1),welch_pgram(x2)) < dist(fm(x1),welch_pgram(x3))

    dist = ModelDistance(fm, RationalOptimalTransportDistance(domain=Continuous(), p=2, interval=(-15., 15)))
    @test dist(x1,x2) < dist(x1,x3)
    dist = RationalOptimalTransportDistance(domain=Continuous(), p=2, interval=(0., 15.))
    @test dist(fm(x1),welch_pgram(x2)) < dist(fm(x1),welch_pgram(x3))

    w = LinRange(0, 2pi*15, 3000)
    dist = ModelDistance(fm, RationalOptimalTransportDistance(domain=Continuous(), p=1, interval=(0., 15)))
    ddist = ModelDistance(fm, DiscretizedRationalDistance(w, sqrt.(SpectralDistances.distmat_euclidean(w,w))))
    @test ddist(x1, x1) == 0
    @test_broken ddist(x1, x2) ≈ dist(x1, x2)

    w = LinRange(0, 2pi*15, 3000)
    dist = ModelDistance(fm, RationalOptimalTransportDistance(domain=Continuous(), p=2, interval=(0., 15)))
    ddist = ModelDistance(fm, DiscretizedRationalDistance(w, (SpectralDistances.distmat_euclidean(w,w))))
    @test ddist(x1, x1) == 0
    @test_broken ddist(x1, x2) ≈ dist(x1, x2)

end

@testset "Histogram" begin
    x1 = randn(100)
    x2 = randn(100) .+ 1
    x3 = randn(100) .+ 2
    dist = SpectralDistances.OptimalTransportHistogramDistance(p=1)
    @test dist(x1,x2) < dist(x1,x3)
    dist = SpectralDistances.OptimalTransportHistogramDistance(p=2)
    @test dist(x1,x2) < dist(x1,x3)
end

@testset "Theoretical scaling from paper" begin
    n = 4
    a1 = randn(n+1); a1[1] = 1
    r1 = SpectralDistances.reflectc.(roots(reverse(a1)))
    r1 = ContinuousRoots(complex.(0.01real.(r1), imag.(r1)))
    # r1 = SpectralDistances.normalize_energy(ContinuousRoots(r1))

    a2 = randn(n+1); a2[1] = 1
    r2 = SpectralDistances.reflectc.(roots(reverse(a2)))
    r2 = ContinuousRoots(complex.(0.01real.(r2), imag.(r2)))
    # r2 = SpectralDistances.normalize_energy(ContinuousRoots(r2))

    function scaleddist(alpha,p)
        r1a,r2a = ContinuousRoots(alpha.*r1), ContinuousRoots(alpha.*r2)
        m1,m2 = AR(r1a), AR(r2a)
        # dist = RationalOptimalTransportDistance(domain=Continuous(), p=p, interval=(-20.,20.))
        rdist = EuclideanRootDistance(domain=Continuous(), weight=residueweight, p=p)
        rdist(m1,m2)
    end
    for α ∈ 1:0.5:3, p ∈ 1:3
        @test α^(1. -2n +p) * scaleddist(1,p) ≈ scaleddist(α,p)

        @test α^(1. -n)*residues(ContinuousRoots(r1)) ≈ residues(ContinuousRoots(α*r1))
        @test α^(1. -2n)*residueweight(r1) ≈ residueweight(α*r1)
    end

end

@testset "Interpolations" begin
    # NOTE: Only testing that it runs, not that it does the right thing. Have a look at the figure in the docs to deteermine if it works well
    n = 4
    r1 = complex.(-0.01 .+ 0.001randn(3), 2randn(3))
    r1 = ContinuousRoots([r1; conj.(r1)])
    r2 = complex.(-0.01 .+ 0.001randn(3), 2randn(3))
    r2 = ContinuousRoots([r2; conj.(r2)])
    r1,r2 = normalize_energy.((r1, r2))
    A1 = AR(r1)
    A2 = AR(r2)
    t = 0.1
    dist = RationalOptimalTransportDistance(domain=Continuous(), p=2, interval=(0., exp10(1.01)))
    interp = SpectralDistances.interpolator(dist, A1, A2)
    w = exp10.(LinRange(-1.5, 1, 300))
    for t = LinRange(0, 1, 7)
        Φ = clamp.(interp(w,t), 1e-10, 100)
    end

    rdist = EuclideanRootDistance(domain=Continuous(), p=2)
    interp = SpectralDistances.interpolator(rdist, A1, A2, normalize=false)
    w = exp10.(LinRange(-1.5, 1, 300))
    for t = LinRange(0, 1, 7)
        Φ = interp(w,t)
    end
end

@testset "barycenter" begin
    @info "Testing barycenter"
    include("test_barycenter.jl")
end

@testset "Slerp" begin
    a = n1([1,2,3])
    b = n1([4,4,3])
    c = SpectralDistances.slerp(a,b,0.5)
    @test c ≈ n1((a+b)/2)
end

@testset "Utils" begin
    X = randn(1000,10)
    x1,x2 = twoD(X)
    @test length(x1) == length(x2) == size(X,1)
    x1,x2,x3 = threeD(X)
    @test length(x1) == length(x3) == size(X,1)

    @test norm(n1(x1)) ≈ 1
    @test sum(s1(abs.(x1))) ≈ 1
    @test var(v1(x1)) ≈ 1
    @test mean(v1(x1)) ≈ 0 atol=sqrt(eps())
    @test mean(SpectralDistances.m1(x1)) ≈ 0 atol=0.08
    @test var(SpectralDistances.m1(x1)) ≈ 1 atol=0.2

    # @test norm(n1(X,1)) ≈ 1
    @test  all(sum(s1(abs.(X),1), dims=1) .≈ 1)
    @test  all(var(v1(X,1), dims=1) .≈ 1)
    @test  mean(v1(X,1)) ≈ 0 atol = sqrt(eps())
    @test  mean(SpectralDistances.m1(X,1)) ≈ 0 atol=0.08
    @test  var(SpectralDistances.m1(X,1)) ≈ 1 atol=0.2
end

@testset "Unbalanced transport" begin
    @info "Testing Unbalanced transport"

    fm = LS(na = 10)
    m1 = fm(filtfilt(ones(10), [10], randn(1000)))
    m2 = fm(filtfilt(ones(5), [5], randn(1000)))
    dist = OptimalTransportRootDistance(domain = Continuous(), p=1, divergence=Balanced())
    d1 = evaluate(dist,m1,m2)
    dist = OptimalTransportRootDistance(domain = Continuous(), p=1, divergence=KL(1.0))
    d2 = evaluate(dist,m1,m2)
    dist = OptimalTransportRootDistance(domain = Continuous(), p=1, divergence=KL(10.0))
    d3 = evaluate(dist,m1,m2)
    dist = OptimalTransportRootDistance(domain = Continuous(), p=1, divergence=KL(0.01))
    d4 = evaluate(dist,m1,m2)

    @test d1 > d3 > d2 > d4
end

end

##


text:
@testset "distribution_free.jl" begin
    @testset "bounds" begin
        # default bounds for base kernels
        CalibrationTests.uniformbound(ExponentialKernel()) == 1
        CalibrationTests.uniformbound(SqExponentialKernel()) == 1
        CalibrationTests.uniformbound(WhiteKernel()) == 1

        # default bounds for kernels with input transformations
        CalibrationTests.uniformbound(SqExponentialKernel() ∘ ScaleTransform(rand())) == 1
        CalibrationTests.uniformbound(ExponentialKernel() ∘ ScaleTransform(rand(10))) == 1

        # default bounds for scaled kernels
        CalibrationTests.uniformbound(42 * ExponentialKernel()) == 42

        # default bounds for tensor product kernels
        kernel = (3.2 * SqExponentialKernel()) ⊗ (2.7 * WhiteKernel())
        CalibrationTests.uniformbound(kernel) == 3.2 * 2.7

        # default bounds for kernel terms
        CalibrationTests.uniformbound(BlockUnbiasedSKCE(kernel)) == 2 * 3.2 * 2.7
    end

    @testset "estimator and estimates" begin
        kernel = (ExponentialKernel() ∘ ScaleTransform(0.1)) ⊗ WhiteKernel()

        for skce in (BiasedSKCE(kernel), UnbiasedSKCE(kernel), BlockUnbiasedSKCE(kernel))
            for nclasses in (2, 10, 100), nsamples in (10, 50, 100)
                # sample predictions and targets
                dist = Dirichlet(nclasses, 1)
                predictions = [rand(dist) for _ in 1:nsamples]
                targets_consistent = [
                    rand(Categorical(prediction)) for prediction in predictions
                ]
                targets_onlyone = ones(Int, length(predictions))

                # for both sets of targets
                for targets in (targets_consistent, targets_onlyone)
                    test = DistributionFreeSKCETest(skce, predictions, targets)

                    @test test.estimator == skce
                    @test test.n == nsamples
                    @test test.estimate ≈ skce(predictions, targets)
                    @test test.bound == CalibrationTests.uniformbound(skce)
                end
            end
        end
    end

    @testset "consistency" begin
        kernel = (ExponentialKernel() ∘ ScaleTransform(0.1)) ⊗ WhiteKernel()
        αs = 0.05:0.1:0.95
        nsamples = 100

        pvalues_consistent = Vector{Float64}(undef, 100)

        for skce in (BiasedSKCE(kernel), UnbiasedSKCE(kernel), BlockUnbiasedSKCE(kernel))
            for nclasses in (2, 10)
                dist = Dirichlet(nclasses, 1)
                predictions = [Vector{Float64}(undef, nclasses) for _ in 1:nsamples]
                targets_consistent = Vector{Int}(undef, nsamples)

                for i in eachindex(pvalues_consistent)
                    # sample predictions and targets
                    for j in 1:nsamples
                        rand!(dist, predictions[j])
                        targets_consistent[j] = rand(Categorical(predictions[j]))
                    end

                    # define test
                    test_consistent = DistributionFreeSKCETest(
                        skce, predictions, targets_consistent
                    )

                    # estimate pvalue
                    pvalues_consistent[i] = pvalue(test_consistent)
                end

                # compute empirical test errors
                errors = ecdf(pvalues_consistent).(αs)
                @test all(((α, p),) -> p < α, zip(αs, errors))
            end
        end
    end
end


text:
export iauBpn2xy
"""
Extract from the bias-precession-nutation matrix the X,Y coordinates
of the Celestial Intermediate Pole.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
    rbpn      double[3][3]  celestial-to-true matrix (Note 1)

Returned:
    x,y       double        Celestial Intermediate Pole (Note 2)

Notes:

    1. The matrix rbpn transforms vectors from GCRS to true equator (and
    CIO or equinox) of date, and therefore the Celestial Intermediate
    Pole unit vector is the bottom row of the matrix.

    2. The arguments x,y are components of the Celestial Intermediate
    Pole unit vector in the Geocentric Celestial Reference System.

Reference:

    "Expressions for the Celestial Intermediate Pole and Celestial
    Ephemeris Origin consistent with the IAU 2000A precession-
    nutation model", Astron.Astrophys. 400, 1145-1154
    (2003)

    n.b. The celestial ephemeris origin (CEO) was renamed "celestial
        intermediate origin" (CIO) by IAU 2006 Resolution 2.

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# void iauBpn2xy(double rbpn[3][3], double *x, double *y)

function iauBpn2xy(rbpn::AbstractMatrix{<:Real})
    # Preallocate return values
    ref_x = Ref{Float64}(0.0)
    ref_y = Ref{Float64}(0.0)

    ccall((:iauBpn2xy, libsofa_c), Cvoid,
        (Ptr{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
        convert(Array{Float64, 2}, rbpn'), ref_x, ref_y)

    return ref_x[], ref_y[]
end


text:
module DCAbuild
using DCAlign
using PlmDCA, FastaIO, GaussDCA
using Statistics, Printf, DelimitedFiles, Logging, StatsBase

include("insertions.jl")      # all insertions related functions
include("utils.jl")	          # i/o functions,
include("build_model.jl")     # main script for learning a Potts model, insertion penalties and gap penalties
include("fasta_utils.jl")	  # fasta utilities
include("seed_align.jl")	  # pipeline to an aligned seed from an unaligned one

function __init__()
	if Sys.which("mafft") === nothing
		 error("Please install mafft binaries at https://mafft.cbrc.jp/alignment/software/mafft-7.453-with-extensions-src.tgz and make the executables available to the system")
	end
	if Sys.which("hmmbuild") === nothing
		error("Please install hmmer suite at http://eddylab.org/software/hmmer/hmmer.tar.gz and make the executables available to the system")
	end
end

end #end module


text:
module TestUnitAsString

using  Modia
using  Test
using  Modia.Measurements
import Modia.MonteCarloMeasurements
import Modia.MonteCarloMeasurements.StaticParticles

v1 = 2.0u"m/s"
unit_v1 = unit(v1)
v1_unit = Modia.unitAsString( unit_v1  )   # = "m*s^-1"
v2_withoutUnit = 2.0
code = :( $v2_withoutUnit*@u_str($v1_unit) )  # = 2.0u"m*s^-1"
v2 = eval(code)
#@show v1
#@show v1_unit
#@show v2

@test v1==v2

mutable struct Data{FloatType <: AbstractFloat}
    velocity::Modia.quantity(FloatType, u"m/s")
end

data = Data{Float64}(2.0u"mm/s")
v3   = data.velocity
#@show v3   # v3 = 0.002 m s^-1

@test Modia.unitAsString( unit(v3) ) == "m*s^-1"

v = 2.0
@test typeof(v) <: AbstractFloat
@test Modia.baseType(typeof(v)) == Float64
@test Modia.isQuantity(             typeof(v)) == false
@test Modia.isMeasurements(          typeof(v)) == false
@test Modia.isMonteCarloMeasurements(typeof(v)) == false
Base.floatmax(typeof(v))

v = 2.0u"m/s"
@test typeof(v) <: Unitful.Quantity
@test Modia.baseType(typeof(v)) == Float64
@test Modia.isQuantity(             typeof(v)) == true
@test Modia.isMeasurements(          typeof(v)) == false
@test Modia.isMonteCarloMeasurements(typeof(v)) == false
Base.floatmax(typeof(v))

### Measurements
v = 2.0 ± 0.1
@test typeof(v) <: Measurements.Measurement
@test Modia.baseType(typeof(v)) == Float64
@test Modia.isQuantity(              typeof(v)) == false
@test Modia.isMeasurements(          typeof(v)) == true
@test Modia.isMonteCarloMeasurements(typeof(v)) == false
Base.floatmax(typeof(v))

v = (2.0 ± 0.1)u"m/s"
@test typeof(v) <: Unitful.Quantity
@test Modia.baseType(typeof(v)) <: Measurements.Measurement
@test Modia.baseType(Modia.baseType(typeof(v))) == Float64
@test Modia.isQuantity(              typeof(v)) == true
@test Modia.isMeasurements(          typeof(v)) == true
@test Modia.isMonteCarloMeasurements(typeof(v)) == false
Base.floatmax(typeof(v))

### MonteCarloMeasurements
const nparticles = 100
normal(mue,sigma) = StaticParticles(nparticles,MonteCarloMeasurements.Distributions.Normal(mue,sigma))

v = normal(2.0, 0.1)
@test typeof(v) <: Modia.MonteCarloMeasurements.AbstractParticles
@test Modia.baseType(typeof(v)) == Float64
@test Modia.isQuantity(              typeof(v)) == false
@test Modia.isMeasurements(          typeof(v)) == false
@test Modia.isMonteCarloMeasurements(typeof(v)) == true
Base.floatmax(typeof(v))

v = normal(2.0, 0.1)u"m/s"
@test typeof(v) <: Modia.MonteCarloMeasurements.AbstractParticles
@test Modia.baseType(typeof(v)) <: Unitful.Quantity
@test Modia.baseType(Modia.baseType(typeof(v))) == Float64
@test Modia.isQuantity(              typeof(v)) == true
@test Modia.isMeasurements(          typeof(v)) == false
@test Modia.isMonteCarloMeasurements(typeof(v)) == true
Base.floatmax(typeof(v))

logCode = false

FirstOrder1 = Model(
    T = 0.2u"s",
    k = 2.0,
    x = Var(init=0.3),
    equations = :[T * der(x) + x = k]
)
FirstOrder2 = FirstOrder1 | Map(k = 2.0u"m/s", x = Var(init = 0.3u"m/s"))
FirstOrder3 = FirstOrder1 | Map(k = 2.0 ± 0.1, x = Var(init = 0.3 ± 0.1))
FirstOrder4 = FirstOrder1 | Map(k = (2.0 ± 0.1)u"m/s" , x = Var(init = (0.3 ± 0.1)u"m/s"))
FirstOrder5 = FirstOrder1 | Map(k = normal(2.0, 0.1/3), x = Var(init = normal(0.3, 0.1/3)))
FirstOrder6 = FirstOrder1 | Map(k = normal(2.0, 0.1/3)u"m/s", x = Var(init = normal(0.3, 0.1/3)u"m/s"))

firstOrder1 = @instantiateModel(FirstOrder1, logCode=logCode)
firstOrder2 = @instantiateModel(FirstOrder2, logCode=logCode)
firstOrder3 = @instantiateModel(FirstOrder3, logCode=logCode, FloatType=Measurements.Measurement{Float64})
firstOrder4 = @instantiateModel(FirstOrder4, logCode=logCode, FloatType=Measurements.Measurement{Float64})
firstOrder5 = @instantiateModel(FirstOrder5, logCode=logCode, FloatType=StaticParticles{Float64,nparticles})
firstOrder6 = @instantiateModel(FirstOrder6, logCode=logCode, FloatType=StaticParticles{Float64,nparticles})

simulate!(firstOrder1)
simulate!(firstOrder2)
simulate!(firstOrder3)
simulate!(firstOrder4)
simulate!(firstOrder5)
simulate!(firstOrder6)

end


text:
using LinearExpressions

using Test

include("tests.jl")