text:
macro params(struct_expr)
    header = struct_expr.args[2]
    fields = @view struct_expr.args[3].args[2:2:end]
    params = []
    for i in 1:length(fields)
        x = fields[i]
        T = gensym()
        if x isa Symbol
            push!(params, T)
            fields[i] = :($x::$T)
        elseif x.head == :(::)
            abstr = x.args[2]
            var = x.args[1]
            push!(params, :($T <: $abstr))
            fields[i] = :($var::$T)
        end
    end
    if header isa Symbol && length(params) > 0
        struct_expr.args[2] = :($header{$(params...)})
    elseif header.head == :curly
        append!(struct_expr.args[2].args, params)
    elseif header.head == :<:
        if struct_expr.args[2].args[1] isa Symbol
            name = struct_expr.args[2].args[1]
            struct_expr.args[2].args[1] = :($name{$(params...)})
        elseif header.head == :<: && struct_expr.args[2].args[1] isa Expr
            append!(struct_expr.args[2].args[1].args, params)
        else
            error("Unidentified type definition.")
        end
    else
        error("Unidentified type definition.")
    end
    esc(struct_expr)
end

@define_cu(IterativeSolvers.CGStateVariables, :u, :r, :c)

struct RaggedArray{TO, TV}
    offsets::TO
    values::TV
end
whichdevice(ra::RaggedArray) = whichdevice(ra.offsets)

function RaggedArray(vv::Vector{Vector{T}}) where T
    offsets = [1; 1 .+ accumulate(+, collect(length(v) for v in vv))]
    values = Vector{T}(undef, offsets[end]-1)
    for (i, v) in enumerate(vv)
        r = offsets[i]:offsets[i+1]-1
        values[r] .= v
    end
    RaggedArray(offsets, values)
end
@define_cu(RaggedArray, :offsets, :values)

function Base.getindex(ra::RaggedArray, i)
    @assert 1 <= i < length(ra.offsets)
    r = ra.offsets[i]:ra.offsets[i+1]-1
    @assert 1 <= r[1] && r[end] <= length(ra.values)
    return @view ra.values[r]
end
function Base.getindex(ra::RaggedArray, i, j)
    @assert 1 <= j < length(ra.offsets)
    r = ra.offsets[j]:ra.offsets[j+1]-1
    @assert 1 <= i <= length(r)
    return ra.values[r[i]]
end
function Base.setindex!(ra::RaggedArray, v, i, j)
    @assert 1 <= j < length(ra.offsets)
    r = ra.offsets[j]:ra.offsets[j+1]-1
    @assert 1 <= i <= length(r)
    ra.values[r[i]] = v
end

function find_varind(black, white, ::Type{TI}=Int) where TI
    nel = length(black)
    nel == length(white) || throw("Black and white vectors should be of the same length")
    varind = zeros(TI, nel)
    k = 1
    for i in 1:nel
        if !black[i] && !white[i]
            varind[i] = k
            k += 1
        end
    end
    return varind
end

function find_black_and_white(dh)
    black = falses(getncells(dh.grid))
    white = falses(getncells(dh.grid))
    if haskey(dh.grid.cellsets, "black")
        for c in grid.cellsets["black"]
            black[c] = true
        end
    end
    if haskey(dh.grid.cellsets, "white")
        for c in grid.cellsets["white"]
            white[c] = true
        end
    end

    return black, white
end

YoungsModulus(p) = getE(p)
PoissonRatio(p) = getν(p)

function compliance(Ke, u, dofs)
    comp = zero(eltype(u))
    for i in 1:length(dofs)
        for j in 1:length(dofs)
            comp += u[dofs[i]]*Ke[i,j]*u[dofs[j]]
        end
    end
    comp
end

function meandiag(K::AbstractMatrix)
    z = zero(eltype(K))
    for i in 1:size(K, 1)
        z += abs(K[i, i])
    end
    return z / size(K, 1)
end

density(var, xmin) = var*(1-xmin) + xmin

macro debug(expr)
    return quote
        if DEBUG[]
            $(esc(expr))
        end
    end
end


text:
using RDatasets;
iris = dataset("datasets", "iris");
iris[!, :is_setosa] = iris.Species .== "setosa";

using MLJ, JLBoostMLJ;
X, y = unpack(iris, x->!(x in [:is_setosa, :Species]), ==(:is_setosa));

using JLBoostMLJ:JLBoostClassifier;
model = JLBoostClassifier()

mljmachine  = machine(model, X, y)

fit!(mljmachine)

predict(mljmachine, X)

feature_importance(fitted_params(mljmachine).fitresult, X, y)

y_cate = categorical(y)

using JLBoost, JLBoostMLJ, MLJ
jlb = JLBoostClassifier()
r1 = range(jlb, :nrounds, lower=1, upper = 6)
r2 = range(jlb, :max_depth, lower=1, upper = 6)
r3 = range(jlb, :eta, lower=0.1, upper=1.0)

tm = TunedModel(model = jlb, ranges = [r1, r2, r3], measure = cross_entropy)
m = machine(tm, X, y_cate)

fit!(m)

fitted_params(m).best_model.max_depth
fitted_params(m).best_model.nrounds
fitted_params(m).best_model.eta

mljmodel = fit(model, 1, X, y)

predict(model, mljmodel.fitresult, X)

feature_importance(mljmodel.fitresult.treemodel, X, y)

ext: jl


text:
using NLsolve, Distributions

@doc raw"""
	CapitalAdjustModel(; kwargs...)

Construct capital investment model with one type of capital and one
productivity shock and convex adjustment costs. The recursive formulation is:

```math
\begin{equation*}
V(K,z) = \max_a K^\alpha e^z - i K - \frac{\gamma}{2} i^2 K + \beta \mathbb{E} V(K', z') \\
i \equiv \frac{a}{K} - (1-\delta) \\
\text{where } K' = a \\
z' = \rho z + \sigma \varepsilon, \quad \varepsilon \sim \mathcal{N}(0,1)
\end{equation*}
```

The firm solves the following problem to initialize at t=0:

```math
\begin{equation*}
V_0(z_0) = \max_{K_0} V(K_0,z_0) - (1 + C0) K_0 \\
\text{where } z_0 = \sqrt{\frac{\sigma^2}{1-\rho^2}} ε_0 \\
ε_0 \sim \mathcal{N}(0,1)
\end{equation*}
```

Example:

```julia
using StructuralDDP
prob = CapitalAdjustModel(nK=150, nz=15, ρ=0.5, σ=0.3, γ=2.)
```

"""
function CapitalAdjustModel(;
	α = 0.67,
	β = 0.9,
	ρ = 0.6,
	σ = 0.3,
	δ = 0.15,
	γ = 2.0,
	C0 = 0.6,
	F = 0.,

	nK = 100,
	nz = 20,

	# options for testing
	intdim = :Separable_ExogStates,
	initialize_exact = true)

    # calculate Amin, Amax
    stdz = sqrt(σ^2/(1-ρ^2))
    minz = -3*stdz
    maxz =  3*stdz
    vz = collect(LinRange(minz, maxz, nz))

    function steadystate(z)
        # check https://www3.nd.edu/~esims1/investment_notes.pdf

        function f2!(F,x)
            (V_K, K_log) = x

            F[1] = - 1 - γ + V_K
            F[2] = V_K - β*( exp(z)*α*exp(K_log)^(α-1) + (1-δ)*V_K)
        end

        return nlsolve(f2!, [1.1; 1.])
    end

    K_ss = exp(steadystate(stdz^2/2).zero[2])
    minK_log = steadystate(-2*stdz).zero[2]
    maxK_log = steadystate(3*stdz).zero[2]
    vK   = exp.(collect(LinRange(minK_log, maxK_log, nK)))

    tStateVectors = (vK, vz) # tuple of basis vectors

	function reward(vStates, vChoices)
	    K, z = vStates
	    Kprime = vChoices[1]
		i = Kprime/K - (1-δ)
		action = Kprime != K
	    return K^α * exp(z) - i*K - F*K*action - γ/2 * i^2 * K
	end

	# including output, i.e. for partial rewardfunc
	function reward(partial_reward, vStates, vChoices)
	    K = vStates[1]
		Kprime = vChoices[1]
	    i = Kprime/K - (1-δ)
		action = Kprime != K
	    return partial_reward - i*K - F*K*action - γ/2 * i^2 * K
	end

	reward_partial(vStates) = vStates[1]^α * exp(vStates[2])

	if intdim == :All

	    transfunc = function mytransfunc(vStates, K, Shock)
	        z = vStates[2]
	        zprime  = ρ*z + σ * Shock;
	        return K, zprime
	    end
		tChoiceVectors = (vK,)

	elseif intdim == :Separable

	    transfunc = function mytransfunc1(vStates, vChoices, Shock)
	        z = vStates[2]
	        zprime  = ρ*z + σ * Shock;
	        return zprime
	    end
		tChoiceVectors = (1,)

	elseif intdim == :Separable_States
	    transfunc = function mytransfunc2(vStates, Shock)
	        z = vStates[2]
	        zprime  = ρ*z + σ * Shock;
	        return zprime
	    end
		tChoiceVectors = (1,)

	elseif intdim == :Separable_ExogStates
	    transfunc = function mytransfunc3(vExogStates, shock)
	        z = vExogStates[1]
	        zprime  = ρ*z + σ * shock;
	        return zprime
	    end
		tChoiceVectors = (1,)
	else

		error("$intdim wrong intdim")
	end

	if initialize_exact
		tChoiceVectorsZero = (1,)

		initializationproblem(value, choice) = value - (1 + (1-β)/β + C0) * choice

		function initialize(shock)
		    z0 = shock * sqrt(σ^2 / (1-ρ^2))
			return z0
		end
		shockdist_initial = Normal()
	else
		tChoiceVectorsZero = nothing
		initializationproblem = nothing
        initialize = nothing
		shockdist_initial = nothing
	end

	DDP(tStateVectors,
        tChoiceVectors,
        reward,
        transfunc,
		Normal(), # give distribution of shocks: standard normal
		β;
		intdim = intdim,
        rewardfunc_partial = reward_partial,
        initializationproblem = initializationproblem,
        initializefunc = initialize,
		shockdist_initial = shockdist_initial,
		tChoiceVectorsZero = tChoiceVectorsZero,
        )
end


text:
"""
    JuliaSerializer <: Construct{Any}

Standard Julia serialization.
"""
struct JuliaSerializer <: Construct{Any} end

deserialize(::JuliaSerializer, s::IO; contextkw...) = Serialization.deserialize(s)
serialize(::JuliaSerializer, s::IO, obj; contextkw...) = Serialization.serialize(s, obj)


text:
module DiskArrays
export AbstractDiskArray, interpret_indices_disk

"""
Abstract DiskArray type that can be inherited by Array-like data structures that
have a significant random access overhead and whose access pattern follows
n-dimensional (hyper)-rectangles.
"""
abstract type AbstractDiskArray{T,N} <: AbstractArray{T,N} end

Base.ndims(::Type{<:AbstractDiskArray{T,N}}) where {T,N} = N

"""
    readblock!(A::AbstractDiskArray, A_ret, r::AbstractUnitRange...)

The only function that should be implemented by a `AbstractDiskArray`. This function
"""
function readblock!() end

"""
    writeblock!(A::AbstractDiskArray, A_in, r::AbstractUnitRange...)

Function that should be implemented by a `AbstractDiskArray` if write operations
should be supported as well.
"""
function writeblock!() end

function readblock!(A::AbstractDiskArray, A_ret, r::AbstractVector...)
  #Implement fallback method if DiskArray does not support strided reading
  #Currently this allocates an intermediate array. In the future, this
  #function could test how sparse the reading is and maybe be smarter here
  #by reading slices.
  mi,ma = map(minimum,r), map(maximum,r)
  A_temp = similar(A_ret,map((a,b)->b-a+1,mi,ma))
  readblock!(A,A_temp,map(:,mi,ma)...)
  A_ret .= view(A_temp,map(ir->ir.-(minimum(ir).-1),r)...)
  nothing
end

function writeblock!(A::AbstractDiskArray, A_ret, r::AbstractVector...)
  #Implement fallback method if DiskArray does not support strided reading
  #Currently this allocates an intermediate array. In the future, this
  #function could test how sparse the reading is and maybe be smarter here
  #by reading slices.
  mi,ma = map(minimum,r), map(maximum,r)
  A_temp = similar(A_ret,map((a,b)->b-a+1,mi,ma))
  A_temp[map(ir->ir.-(minimum(ir).-1),r)...] = A_ret
  writeblock!(A,A_temp,map(:,mi,ma)...)
  nothing
end

Base.ndims(::AbstractDiskArray{<:Any,N}) where N = N
Base.eltype(::AbstractDiskArray{T}) where T = T

function getindex_disk(a, i...)
  inds, trans = interpret_indices_disk(a,i)
  data = Array{eltype(a)}(undef, map(length,inds)...)
  readblock!(a,data,inds...)
  trans(data)
end

setindex_disk!(a::AbstractDiskArray{T},v::T,i...) where T<:AbstractArray =
  setindex_disk!(a,[v],i...)

function setindex_disk!(a::AbstractDiskArray,v::AbstractArray,i...)
  inds, trans = interpret_indices_disk(a,i)
  data = reshape(v,map(length,inds))
  writeblock!(a,data,inds...)
  v
end

"""
Function that translates a list of user-supplied indices into plain ranges and
integers for reading blocks of data. This function respects additional indexing
rules like omitting additional trailing indices.

The passed array handle A must implement methods for `Base.size` and `Base.ndims`
The function returns two values:

  1. a tuple whose length equals `ndims(A)` containing only unit
  ranges and integers. This contains the minimal "bounding box" of data that
  has to be read from disk.
  2. A callable object which transforms the hyperrectangle read from disk to
  the actual shape that represents the Base getindex behavior.
"""
function interpret_indices_disk(A, r::Tuple)
  throw(ArgumentError("Indices of type $(typeof(r)) are not yet "))
end

#Read the entire array and reshape to 1D in the end
function interpret_indices_disk(A, ::Tuple{Colon})
  return map(Base.OneTo, size(A)), Reshaper(prod(size(A)))
end

interpret_indices_disk(A, r::Tuple{<:CartesianIndex}) =
  interpret_indices_disk(A,r[1].I)

interpret_indices_disk(A, r::Tuple{<:CartesianIndices}) =
  interpret_indices_disk(A,r[1].indices)

function interpret_indices_disk(A, r::NTuple{N, Union{Integer, AbstractVector, Colon}}) where N
  if ndims(A)==N
    inds = map(_convert_index,r,size(A))
    resh = DimsDropper(findints(r))
    return inds, resh
  elseif ndims(A)<N
    foreach((ndims(A)+1):N) do i
      r[i]==1 || throw(BoundsError(A, r))
    end
    _, rshort = commonlength(size(A),r)
    return interpret_indices_disk(A, rshort)
  else
    size(A,N+1)==1 || throw(BoundsError(A, r))
    return interpret_indices_disk(A, (r...,1))
  end
end

function interpret_indices_disk(A, r::Tuple{<:AbstractArray{<:Bool}})
  ba = r[1]
  if ndims(A)==ndims(ba)
    inds = getbb(ba)
    resh = a -> a[view(ba,inds...)]
    return inds, resh
  elseif ndims(ba)==1
    interpret_indices_disk(A,(reshape(ba,size(A)),))
  else
    throw(BoundsError(A, r))
  end
end

function interpret_indices_disk(A, r::NTuple{1, AbstractVector})
  lininds = first(r)
  cartinds = CartesianIndices(A)
  mi, ma = extrema(view(cartinds, lininds))
  inds = map((i1,i2) -> i1:i2, mi.I,ma.I)
  resh = a -> map(lininds) do ii
      a[cartinds[ii] - mi + oneunit(mi)]
  end
  inds, resh
end

struct Reshaper{I}
  reshape_indices::I
end
(r::Reshaper)(a) = reshape(a,r.reshape_indices)
struct DimsDropper{D}
  d::D
end
(d::DimsDropper)(a) = length(d.d)==ndims(a) ? a[1] : dropdims(a,dims=d.d)
struct TransformStack{S}
  s::S
end
(s::TransformStack)(a) = ∘(s.s...)(a)

function getbb(ar::AbstractArray{Bool})
  maxval = CartesianIndex(size(ar))
  minval = CartesianIndex{ndims(ar)}()
  reduceop = (i1,i2)->begin i2===nothing ? i1 : (min(i1[1],i2),max(i1[2],i2)) end
  mi,ma = mapfoldl(reduceop,
    zip(CartesianIndices(ar),ar),
    init = (maxval,minval)) do ii
    ind,val = ii
    val ? ind : nothing
  end
  inds = map((i1,i2) -> i1:i2, mi.I,ma.I)
end

#Some helper functions
"For two given tuples return a truncated version of both so they have common length"
commonlength(a, b) = _commonlength((first(a),), (first(b),),Base.tail(a), Base.tail(b))
commonlength(::Tuple{}, b) = (),()
commonlength(a, ::Tuple{}) = (),()
_commonlength(a1, b1, a, b) = _commonlength((a1..., first(a)), (b1..., first(b)), Base.tail(a), Base.tail(b))
_commonlength(a1, b1, ::Tuple{}, b) = (a1,b1)
_commonlength(a1, b1, a, ::Tuple{}) = (a1, b1)

"Find the indices of elements containing integers in a Tuple"
findints(x) = _findints((),1,x...)
_findints(c, i, x::Integer, rest...) = _findints((c..., i), i+1, rest...)
_findints(c, i, x, rest...) = _findints(c, i+1, rest...)
_findints(c, i)  = c
#Normal indexing for a full subset of an array
_convert_index(i::Integer, s::Integer) = i:i
_convert_index(i::AbstractVector, s::Integer) = i
_convert_index(::Colon, s::Integer) = Base.OneTo(Int(s))

macro implement_getindex(t)
quote
function Base.getindex(a::$t,i...)
  getindex_disk(a,i...)
end
end
end

macro implement_setindex(t)
quote
Base.setindex!(a::$t,v::AbstractArray,i...) =
  setindex_disk!(a,v,i...)

# Add an extra method if a single number is given
Base.setindex!(a::$t{<:Any,N}, v, i...) where N =
  Base.setindex!(a,fill(v,ntuple(i->1,N)...), i...)

#Special care must be taken for logical indexing, we can not avoid reading the data
#before writing
function Base.setindex!(a::$t,v::AbstractArray,i::AbstractArray{<:Bool})
  inds, trans = interpret_indices_disk(a,(i,))
  data = Array{eltype(a)}(undef, map(length,inds)...)
  readblock!(a,data,inds...)

  data = reshape(v,map(length,inds))
  writeblock!(a,data,inds...)
  v
end
end
end

function Base.show(io::IO, ::MIME"text/plain", X::AbstractDiskArray)
  println(io, "Disk Array with size ", join(size(X)," x "))
end
function Base.show(io::IO, X::AbstractDiskArray)
  println(io, "Disk Array with size ", join(size(X)," x "))
end

include("chunks.jl")
include("ops.jl")
include("iterator.jl")
include("subarrays.jl")
include("permute_reshape.jl")

#The all-in-one macro
macro implement_diskarray(t)
quote
  @implement_getindex $t
  @implement_setindex $t
  if VERSION >= v"1.3.0"
    @implement_broadcast $t
    @implement_iteration $t
    @implement_mapreduce $t
    @implement_reshape $t
    @implement_permutedims $t
  end
end
end

@implement_diskarray AbstractDiskArray
end # module


text:
using Random

mutable struct Species
    # A subpopulation containing similar individiduals
    id::Int64                         # species's id
    age::Int64                        # species's age
    subpopulation::Vector{MLPChromosome} # species's individuals
    hasBest::Bool                     # Does this species has the best individual of the population?
    spawn_amount::Int64
    no_improvement_age::Int64         # the age species has shown no improvements on average
    last_avg_fitness::Float64
    representant::MLPChromosome
    function Species(g::Global, first_individual::MLPChromosome, previous_id=0)
        # A species requires at least one individual to come to existence
        id  = previous_id == 0 ? get_new_id(g) : previous_id
        first_individual.species_id = id
        new(id,
            0,     # species's age
            [first_individual],    # species' individuals
            false, # hasBest: Does this species has the best individual of the population?
            0,     # spawn_amount
            0,     # the age species has shown no improvements on average
            0.,    # last_avg_fitness
            first_individual)
    end
end

function get_new_id(g::Global)
    g.speciesCnt += 1
    return g.speciesCnt
end

function add(s::Species, individual::MLPChromosome)
    individual.species_id = s.id
    push!(s.subpopulation, individual)
    # choose a new random representant for the species
    s.representant = s.subpopulation[rand(1:length(s.subpopulation))]
end

# Returns the total number of individuals in this species
Base.length(s::Species) = length(s.subpopulation)

function Base.show(io::IO, s::Species)
    str = @sprintf("   Species: %2d   size: %3d   age: %3d   spawn: %3d  ",
                    s.id, length(s.subpopulation), s.age, s.spawn_amount)
    str = @sprintf("%1s  No improvement: %3d   avg. fitness: %1.8f",
                    str,s.no_improvement_age, s.last_avg_fitness)
    @printf(io,"%6s", str)
end

function tournamentSelection(s::Species, k=2)
    # Tournament selection with size k (default k=2).
    # randomly select k competitors
    selected = randperm(length(s.subpopulation))[1:k]
    chs = s.subpopulation[selected]
    best = chs[1]
    for ch in chs # choose best among randomly selected
        best = ch.fitness > best.fitness ? ch : best
    end
    return best
end

function average_fitness(s::Species)
    # Returns the raw average fitness for this species
    @assert length(s) > 0
    fitness_sum = 0.0
    fitness_sum = sum((ch)-> ch.fitness, s.subpopulation)
    current = fitness_sum / length(s)

    # if no_improvement_age > threshold, species will be removed
    if current > s.last_avg_fitness
        s.last_avg_fitness = current
        s.no_improvement_age = 0
    else
        s.no_improvement_age += 1
    end
    return current
end

function reproduce(g::Global, s::Species, fitness::Evaluation)
    # Returns a list of 'spawn_amount' new individuals

    offspring = Tuple{MLPChromosome, Bool}[] # new offspring for this species
    s.age += 1  # increment species age

#     @printf("Reproducing species %d with %d members\n", s.id, length(s))

    # this condition is useless since no species with spawn_amount < 0 will
    # reach this point - at least it shouldn't happen.
    @assert s.spawn_amount > 0

    sort!(s.subpopulation, by = ch-> ch.fitness, rev = true)

    if g.cf.elitism
        # TODO: Wouldn't it be better if we set elitism=2,3,4...
        # depending on the size of each species?
        push!(offspring, (s.subpopulation[1], true))
        s.spawn_amount -= 1
    end

    survivors = floor(Int, length(s) * g.cf.survival_threshold) # keep a % of the best individuals
#     println("survivors = $survivors $(length(s))")

    s.subpopulation = survivors > 0 ? s.subpopulation[1:survivors] : [s.subpopulation[1]]

    while(s.spawn_amount > 0)

        s.spawn_amount -= 1

        if length(s) > 1
            # Selects two parents from the remaining species and produces a single individual
            # Stanley selects at random, here we use tournament selection (although it is not
            # clear if has any advantages)
            parent1 = tournamentSelection(s)
            parent2 = tournamentSelection(s)

            @assert parent1.species_id == parent2.species_id
            child = crossover(g, parent1, parent2)
            child = mutate!(child,g)
            fitness(child)
            push!(offspring, (child, false))
        else
            # mutate only
            parent1 = s.subpopulation[1]
            # TODO: temporary hack - the child needs a new id (not the father's)
            child = crossover(g, parent1, parent1)
            child = mutate!(child,g)
            fitness(child)
            #push!(offspring, child)
            push!(offspring, (child, false))
        end
    end

    # reset species (new members will be added again when speciating)
    s.subpopulation = []

    # select a new random representant member
    representant_index = rand(1:length(offspring))
    s.representant = offspring[representant_index][1]
    #offspring[representant_index][2] = true
    offspring[representant_index] = (offspring[representant_index][1], true)

    return offspring
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
# NOTE: This file is for utilities used in tests — not for tests of utilities

"""
This module is copied from [this post by Brian Guenter](
https://discourse.julialang.org/t/treating-nan-as-error-helping-debugging/36933/9).
Essentially, you can wrap your `Real` type in `NaNCheck`, and pass it (or an
array) of such objects to a function you want to check to track down where
NaNs are coming into a calculation.  As soon as such an object is used in a
calculation, it raises an exception — essentially making this a
[signaling NaN](https://en.wikipedia.org/wiki/NaN#Signaling_NaN).

This can also be helpful for ensuring that no un-initialized memory is used in
a computation (which might silently give bad results).  Instead of creating an
`undef` array, just create the array from `NaNCheck`s.  The elements can be
set to some non-NaN value before use, and the calculation will proceed as
usual.  If they are not set to non-NaN, they will raise an exception whenever
they are used.

Such an array can be created for type `T` with something like

    arr = fill(NaNCheck{T}(NaN), arr_size)

"""
module NaNChecker

struct NaNCheck{T<:Real} <: Real
    val::T
    function NaNCheck{T}(a::S) where {T<:Real, S<:Real}
        @assert !(T <: NaNCheck)
        new{T}(T(a))
    end
end
export NaNCheck
Base.isnan(a::NaNCheck{T}) where{T} = isnan(a.val)
Base.isinf(a::NaNCheck{T}) where{T} = isinf(a.val)
Base.typemin(::Type{NaNCheck{T}}) where{T} = NaNCheck{T}(typemin(T))
Base.typemax(::Type{NaNCheck{T}}) where{T} = NaNCheck{T}(typemax(T))
Base.eps(::Type{NaNCheck{T}}) where {T} = NaNCheck{T}(eps(T))
Base.decompose(a::NaNCheck{T}) where {T} = Base.decompose(a.val)
Base.round(a::NaNCheck{T}, m::RoundingMode) where {T} = NaNCheck{T}(round(a.val, m))

struct NaNException <: Exception end

# (::Type{Float64})(a::NaNCheck{S}) where {S<:Real} = NaNCheck{Float64}(Float64(a.val))
(::Type{T})(a::NaNCheck{S}) where {T<:Integer,S<:Real} = T(a.val)
(::Type{NaNCheck{T}})(a::NaNCheck{S}) where {T<:Real,S<:Real} = NaNCheck{T}(T(a.val))
Base.promote_rule(::Type{NaNCheck{T}}, ::Type{T}) where {T<:Number} = NaNCheck{T}
Base.promote_rule(::Type{T}, ::Type{NaNCheck{T}}) where {T<:Number} = NaNCheck{T}
Base.promote_rule(::Type{S}, ::Type{NaNCheck{T}}) where {T<:Number, S<:Number} = NaNCheck{promote_type(T,S)}
Base.promote_rule(::Type{NaNCheck{T}}, ::Type{S}) where {T<:Number, S<:Number} = NaNCheck{promote_type(T,S)}
Base.promote_rule(::Type{NaNCheck{S}}, ::Type{NaNCheck{T}}) where {T<:Number, S<:Number} = NaNCheck{promote_type(T,S)}

for op = (:sin, :cos, :tan, :log, :exp, :sqrt, :abs, :-, :atan, :acos, :asin, :log1p, :floor, :ceil, :float)
    eval(quote
        function Base.$op(a::NaNCheck{T}) where{T}
            temp = NaNCheck{T}(Base.$op(a.val))
            if isnan(temp)
                throw(NaNException())
            end
            return temp
        end
    end)
end

for op = (:+, :-, :/, :*, :^, :atan)
    eval(quote
        function Base.$op(a::NaNCheck{T}, b::NaNCheck{T}) where{T}
            temp = NaNCheck{T}(Base.$op(a.val, b.val))
            if isnan(temp)
                throw(NaNException())
            end
            return temp
        end
    end)
end

for op =  (:<, :>, :<=, :>=, :(==), :isequal)
    eval(quote
        function Base.$op(a::NaNCheck{T}, b::NaNCheck{T}) where{T}
            temp = Base.$op(a.val, b.val)
            return temp
        end
    end)
end

end  # module NaNChecker
NaNCheck = NaNChecker.NaNCheck;

module ExplicitWignerMatrices

function d_explicit(n, m′, m, expiβ::Complex{T}) where T
    if abs(m′) < abs(m)
        return (-1)^(m-m′) * d_explicit(n, m, m′, expiβ)
    end
    if m′ < 0
        return (-1)^(m-m′) * d_explicit(n, -m′, -m, expiβ)
    end
    cosβ = expiβ.re
    sinβ = expiβ.im
    if (n,m′,m) == (0,0,0)
        T(1)
    elseif (n,m′,m) == (1,0,0)
        cosβ
    elseif (n,m′,m) == (1,1,-1)
        (1-cosβ) / 2
    elseif (n,m′,m) == (1,1,0)
        -sinβ / √T(2)
    elseif (n,m′,m) == (1,1,1)
        (1+cosβ) / 2
    elseif (n,m′,m) == (2,0,0)
        (3cosβ^2-1) / 2
    elseif (n,m′,m) == (2,1,-1)
        (1+cosβ-2cosβ^2) / 2
    elseif (n,m′,m) == (2,1,0)
        -√(T(3)/8) * 2 * sinβ * cosβ
    elseif (n,m′,m) == (2,1,1)
        (-1+cosβ+2cosβ^2) / 2
    elseif (n,m′,m) == (2,2,-2)
        (1-cosβ)^2 / 4
    elseif (n,m′,m) == (2,2,-1)
        -sinβ * (1-cosβ) / 2
    elseif (n,m′,m) == (2,2,0)
        √(T(3)/8) * sinβ^2
    elseif (n,m′,m) == (2,2,1)
        -sinβ * (1+cosβ) / 2
    elseif (n,m′,m) == (2,2,2)
        (1+cosβ)^2/4
    else
        T(NaN)
    end
end

function D_explicit(n, m′, m, expiα::Complex{T}, expiβ::Complex{T}, expiγ::Complex{T}) where T
    return expiα^(-m′) * d_explicit(n, m′, m, expiβ) * expiγ^(-m)
end

function d_formula(n, m′, m, expiβ::Complex{T}) where T
    # https://en.wikipedia.org/wiki/Wigner_D-matrix#Wigner_.28small.29_d-matrix
    cosβ = expiβ.re
    sin½β = √((1-cosβ)/2)
    cos½β = √((1+cosβ)/2)
    prefactor =√T(
        factorial(big(n + m′))
        * factorial(big(n - m′))
        * factorial(big(n + m))
        * factorial(big(n - m))
    )
    prefactor * sum(
        ifelse(iseven(m′ - m + s), 1, -1)
        * cos½β ^ (2n + m - m′ - 2s)
        * sin½β ^ (m′ - m + 2s)
        / (
            factorial(big(n + m - s))
            * factorial(big(s))
            * factorial(big(m′ - m + s))
            * factorial(big(n - m′ - s))
        )
        for s in max(0, m - m′):min(n + m, n - m′)
    )
end

function D_formula(n, m′, m, expiα::Complex{T}, expiβ::Complex{T}, expiγ::Complex{T}) where T
    # https://en.wikipedia.org/wiki/Wigner_D-matrix#Definition_of_the_Wigner_D-matrix
    # Note that the convention in this package is conjugated relative to the convention
    # used by Wikipedia, so we include that conjugation here.
    return expiα^(m′) * d_formula(n, m′, m, expiβ) * expiγ^(m)
end

end  # module ExplicitWignerMatrices

module NINJA

function Wigner_d(ι::T, ell, m, s) where {T<:Real}
    # Eq. II.8 of Ajith et al. (2007) 'Data formats...'
    k_min = max(0, m - s)
    k_max = min(ell + m, ell - s)
    prefactor = √T(
        factorial(big(ell + m))
        * factorial(big(ell - m))
        * factorial(big(ell + s))
        * factorial(big(ell - s))
    )
    sum(
        ifelse(iseven(k), 1, -1)
         * cos(ι / 2) ^ (2 * ell + m - s - 2 * k)
         * sin(ι / 2) ^ (2 * k + s - m)
         * prefactor
         / T(
             factorial(big(ell + m - k))
             * factorial(big(ell - s - k))
             * factorial(big(k))
             * factorial(big(k + s - m))
         )
        for k in k_min:k_max
    )
end

function sYlm(s, ell, m, ι::T, ϕ::T) where {T<:Real}
    # Eq. II.7 of Ajith et al. (2007) 'Data formats...'
    # Note the weird definition w.r.t. `-s`
    if abs(s) > ell || abs(m) > ell
        return zero(complex(T))
    end
    (
        ifelse(iseven(s), 1, -1)
        * √((2ell + 1) / (4T(π)))
        * Wigner_d(ι, ell, m, -s)
        * cis(m * ϕ)
    )
end

# Eqs. (II.9) through (II.13) of https://arxiv.org/abs/0709.0093v3
m2Y22(ι::T, ϕ::T) where {T<:Real} = √(5 / (64T(π))) * (1 + cos(ι))^2 * cis(2ϕ)
m2Y21(ι::T, ϕ::T) where {T<:Real} = √(5 / (16T(π))) * sin(ι) * (1 + cos(ι)) * cis(ϕ)
m2Y20(ι::T, ϕ::T) where {T<:Real} = √(15 / (32T(π))) * sin(ι)^2
m2Y2m1(ι::T, ϕ::T) where {T<:Real} = √(5 / (16T(π))) * sin(ι) * (1 - cos(ι)) * cis(-1ϕ)
m2Y2m2(ι::T, ϕ::T) where {T<:Real} = √(5 / (64T(π))) * (1 - cos(ι))^2 * cis(-2ϕ)

end  # module NINJA

αrange(::Type{T}, n=15) where T = T[
    0; nextfloat(T(0)); rand(T(0):eps(T(π)):T(π), n÷2); prevfloat(T(π)); T(π);
    nextfloat(T(π)); rand(T(π):eps(2T(π)):2T(π), n÷2); prevfloat(T(π)); 2T(π)
]
βrange(::Type{T}, n=15) where T = T[
    0; nextfloat(T(0)); rand(T(0):eps(T(π)):T(π), n); prevfloat(T(π)); T(π)
]
γrange(::Type{T}, n=15) where T = αrange(T, n)
v̂range(::Type{T}, n=15) where T = QuatVec{T}[𝐢; 𝐣; 𝐤; -𝐢; -𝐣; -𝐤; normalize.(randn(QuatVec{T}, n))]
function Rrange(::Type{T}, n=15) where T
    invsqrt2 = inv(√T(2))
    [
        [
            sign*R
            for R in [
                Rotor{T}(1);
                [Rotor{T}(𝐯) for 𝐯 in (𝐢,𝐣,𝐤)];
                [Rotor{T}(invsqrt2 + invsqrt2*𝐯) for 𝐯 in (𝐢,𝐣,𝐤)];
                [Rotor{T}(invsqrt2 - invsqrt2*𝐯) for 𝐯 in (𝐢,𝐣,𝐤)]
        ]
        for sign in (1,-1)
        ];
        randn(Rotor{T}, n)
    ]
end
epsilon(k) = ifelse(k>0 && isodd(k), -1, 1)

"""
    array_equal(a1, a2, equal_nan=false)

Ensure that arrays have same types and shapes, and all elements are the same.
If `equal_nan` is `true`, NaNs in the same place in each array will be
considered to be equal.

Note that this is slightly stricter than the numpy version of this function,
because arrays of different type will not be considered equal.

"""
function array_equal(a1::T1, a2::T2, equal_nan=false) where {T1, T2}
    if T1 !== T2 || size(a1) != size(a2)
        return false
    end
    all(e->e[1]==e[2] || (equal_nan && isnan(e1) && isnan(e2)), zip(a1, a2))
end


text:
@testset "$(rpad("Deterministic equations",80))" begin

    ################################################################################
    # Test ODE: Ordinary Differential Equation
    ################################################################################

    function f_ode(t, x, f)
        f[1] = x[1]
    end

    ode  = ODE(eltype(q₀), 1, 1, 1, f_ode, t₀, q₀)
    ode1 = ODE(f_ode, t₀, q₀)
    ode2 = ODE(f_ode, q₀)

    @test ode == ode1
    @test ode == ode2

    @test hash(ode1) == hash(ode2)

    @test ode == similar(ode, t₀, q₀)
    @test ode == similar(ode, q₀)

    ################################################################################
    # Test PODE: Partitioned Ordinary Differential Equation
    ################################################################################

    function v_pode(t, q, p, v)
        v[1] = q[1]
    end

    function f_pode(t, q, p, f)
        f[1] = 2p[1]
    end

    pode  = PODE(eltype(q₀), 1, 1, 1, v_pode, f_pode, t₀, q₀, p₀)
    pode1 = PODE(v_pode, f_pode, t₀, q₀, p₀)
    pode2 = PODE(v_pode, f_pode, q₀, p₀)

    @test pode == pode1
    @test pode == pode2

    @test hash(pode1) == hash(pode2)

    @test pode == similar(pode, t₀, q₀, p₀)
    @test pode == similar(pode, q₀, p₀)

    ################################################################################
    # Test IODE: Implicit Ordinary Differential Equation
    ################################################################################

    function iode_ϑ(t, q, v, p)
        p[1] = v[1]
    end

    function iode_f(t, q, v, f)
        f[1] = sin(q[1])
    end

    function iode_g(t, q, λ, g)
        g[1] = λ[1]
    end

    function iode_v(t, q, p, v)
        v[1] = p[1]
    end

    iode  = IODE(eltype(q₀), 1, 1, 1, iode_ϑ, iode_f, iode_g, t₀, q₀, p₀, λ₀; v=iode_v)
    iode1 = IODE(iode_ϑ, iode_f, iode_g, t₀, q₀, p₀, λ₀; v=iode_v)
    iode2 = IODE(iode_ϑ, iode_f, iode_g, t₀, q₀, p₀; v=iode_v)
    iode3 = IODE(iode_ϑ, iode_f, iode_g, q₀, p₀; v=iode_v)

    @test iode == iode1
    @test iode == iode2
    @test iode == iode3

    @test hash(iode1) == hash(iode2)

    @test iode == similar(iode, t₀, q₀, p₀, λ₀)
    @test iode == similar(iode, t₀, q₀, p₀)
    @test iode == similar(iode, q₀, p₀)

    ################################################################################
    # Test DAE: Differential Algebraic Equation
    ################################################################################

    function v_dae(t, x, v)
        v[1] = x[1]
        v[2] = x[2]
    end

    function u_dae(t, x, λ, u)
        u[1] = +λ[1]
        u[2] = -λ[1]
    end

    function ϕ_dae(t, x, λ, ϕ)
        ϕ[1] = x[2] - x[1]
    end

    dae  = DAE(eltype(q₀), 1, 2, 1, 1, v_dae, u_dae, ϕ_dae, t₀, x₀, λ₀)
    dae1 = DAE(v_dae, u_dae, ϕ_dae, t₀, x₀, λ₀)
    dae2 = DAE(v_dae, u_dae, ϕ_dae, x₀, λ₀)

    @test dae == dae1
    @test dae == dae2

    @test hash(dae1) == hash(dae2)

    @test dae == similar(dae, t₀, x₀, λ₀)
    @test dae == similar(dae, x₀, λ₀)

    ################################################################################
    # Test PDAE: Partitioned Differential Algebraic Equation
    ################################################################################

    function v_pdae(t, q, p, v)
        v[1] = q[1]
    end

    function f_pdae(t, q, p, f)
        f[1] = p[1]
    end

    function p_pdae(t, q, v, p)
        p[1] = v[1]
    end

    function u_pdae(t, q, p, λ, u)
        u[1] = λ[1]
    end

    function g_pdae(t, q, p, λ, g)
        g[1] = λ[1]
    end

    function ϕ_pdae(t, q, p, λ, ϕ)
        ϕ[1] = p[1] - q[1]
    end

    pdae  = PDAE(eltype(q₀), 1, 1, 1, 1, v_pdae, f_pdae, u_pdae, g_pdae, ϕ_pdae, t₀, q₀, p₀, λ₀)
    pdae1 = PDAE(v_pdae, f_pdae, u_pdae, g_pdae, ϕ_pdae, t₀, q₀, p₀, λ₀)
    pdae2 = PDAE(v_pdae, f_pdae, u_pdae, g_pdae, ϕ_pdae, q₀, p₀, λ₀)

    @test pdae == pdae1
    @test pdae == pdae2

    @test hash(pdae1) == hash(pdae2)

    @test pdae == similar(pdae, t₀, q₀, p₀, λ₀)
    @test pdae == similar(pdae, t₀, q₀, p₀)
    @test pdae == similar(pdae, q₀, p₀)

    idae  = IDAE(eltype(q₀), 1, 1, 1, 1, p_pdae, f_pdae, u_pdae, g_pdae, ϕ_pdae, t₀, q₀, p₀, λ₀)
    idae1 = IDAE(p_pdae, f_pdae, u_pdae, g_pdae, ϕ_pdae, t₀, q₀, p₀, λ₀)
    idae2 = IDAE(p_pdae, f_pdae, u_pdae, g_pdae, ϕ_pdae, q₀, p₀, λ₀)

    @test idae == idae1
    @test idae == idae2

    @test hash(idae1) == hash(idae2)

    @test idae == similar(idae, t₀, q₀, p₀, λ₀)
    @test idae == similar(idae, t₀, q₀, p₀)
    @test idae == similar(idae, q₀, p₀)

    ################################################################################
    # Test SODE: Split Ordinary Differential Equation
    ################################################################################

    function f_sode_1(t, x, f)
        f[1] = x[1]
    end

    function f_sode_2(t, x, f)
        f[1] = x[1]^2
    end

    f_sode = (f_sode_1, f_sode_2)

    sode  = SODE(eltype(q₀), 1, 1, 1, f_sode, t₀, q₀)
    sode1 = SODE(f_sode, t₀, q₀)
    sode2 = SODE(f_sode, q₀)

    @test sode == sode1
    @test sode == sode2

    @test hash(sode1) == hash(sode2)

    @test sode == similar(sode, t₀, q₀)
    @test sode == similar(sode, q₀)

end


text:
module ApproxFunRational
using Base, ApproxFun, ApproxFunBase, ApproxFunFourier, Reexport, FFTW, LinearAlgebra, ApproxFunOrthogonalPolynomials, ApproxFunRational#, MacroTools#, Reexport, AbstractFFTs, FFTW, InfiniteArrays, FillArrays, FastTransforms, IntervalSets,
            #DomainSets, SpecialFunctions

@reexport using ApproxFunBase

#import ApproxFunBase: Fun, AbstractTransformPlan, TransformPlan, ITransformPlan, ConcreteDerivative, ConcreteMultiplication, prectype, spacescompatible

import ApproxFunBase: normalize!, flipsign, FiniteRange, Fun, MatrixFun, UnsetSpace, VFun, RowVector,
                UnivariateSpace, AmbiguousSpace, SumSpace, SubSpace, WeightSpace, NoSpace, Space,
                HeavisideSpace, PointSpace, dimension, colstop, israggedbelow, rowstop,
                IntervalOrSegment, RaggedMatrix, AlmostBandedMatrix, chop, chop!,
                AnyDomain, ZeroSpace, ArraySpace, TrivialInterlacer, BlockInterlacer,
                AbstractTransformPlan, TransformPlan, ITransformPlan,
                ConcreteConversion, ConcreteMultiplication, ConcreteDerivative, ConcreteIntegral, CalculusOperator,
                ConcreteVolterra, Volterra, VolterraWrapper,
                MultiplicationWrapper, ConversionWrapper, DerivativeWrapper, Evaluation, EvaluationWrapper,
                Conversion, defaultConversion, defaultcoefficients, default_Fun, Multiplication, Derivative, Integral, bandwidths,
                ConcreteEvaluation, ConcreteDefiniteLineIntegral, ConcreteDefiniteIntegral, ConcreteIntegral,
                DefiniteLineIntegral, DefiniteIntegral, ConcreteDefiniteIntegral, ConcreteDefiniteLineIntegral, IntegralWrapper,
                ReverseOrientation, ReverseOrientationWrapper, ReverseWrapper, Reverse, NegateEven,
                Dirichlet, ConcreteDirichlet, DirichletWrapper,
                TridiagonalOperator, SubOperator, Space, @containsconstants, spacescompatible,
                hasfasttransform, canonicalspace, domain, setdomain, prectype, domainscompatible,
                plan_transform, plan_itransform, plan_transform!, plan_itransform!, transform, itransform, hasfasttransform,
                CanonicalTransformPlan, ICanonicalTransformPlan,
                Integral, ∞,
                domainspace, rangespace, boundary,
                union_rule, conversion_rule, maxspace_rule, conversion_type, maxspace, hasconversion, points,
                rdirichlet, ldirichlet, lneumann, rneumann, ivp, bvp,
                linesum, differentiate, integrate, linebilinearform, bilinearform,
                UnsetNumber, coefficienttimes, subspace_coefficients, sumspacecoefficients, specialfunctionnormalizationpoint,
                Segment, IntervalOrSegmentDomain, PiecewiseSegment, isambiguous, Vec, eps, isperiodic,
                arclength, complexlength,
                invfromcanonicalD, fromcanonical, tocanonical, fromcanonicalD, tocanonicalD, canonicaldomain, setcanonicaldomain, mappoint,
                reverseorientation, checkpoints, evaluate, mul_coefficients, coefficients, coefficientmatrix, isconvertible,
                clenshaw, ClenshawPlan, sineshaw,
                toeplitz_getindex, toeplitz_axpy!, sym_toeplitz_axpy!, hankel_axpy!, ToeplitzOperator, SymToeplitzOperator, hankel_getindex,
                SpaceOperator, ZeroOperator, InterlaceOperator,
                interlace!, reverseeven!, negateeven!, cfstype, pad!, alternatesign!, mobius,
                extremal_args, hesseneigvals, chebyshev_clenshaw, recA, recB, recC, roots,splitatroots,
                chebmult_getindex, intpow, alternatingsum,
                domaintype, diagindshift, rangetype, weight, isapproxinteger, default_Dirichlet, scal!, dotu,
                components, promoterangespace, promotedomainspace, choosedomainspace,
                block, blockstart, blockstop, blocklengths, isblockbanded, pointscompatible,
                AbstractProductSpace, MultivariateFun, BivariateSpace,
                @wrapperstructure, @wrapperspaces, @wrapper, @calculus_operator, resizedata!, slnorm, affine_setdiff,
                complexroots
import ApproxFunFourier: PeriodicDomain, PeriodicLine, Circle
import FFTW: plan_r2r!, fftwNumber, REDFT10, REDFT01, REDFT00, RODFT00, R2HC, HC2R,
                r2r!, r2r,  plan_fft, plan_ifft, plan_ifft!, plan_fft!, cFFTWPlan
import Base: values, convert, getindex, setindex!, *, +, -, ==, <, <=, >, |, !, !=, eltype, iterate,
                                >=, /, ^, \, ∪, transpose, size, tail, broadcast, broadcast!, copyto!, copy, to_index, (:),
                                similar, map, vcat, hcat, hvcat, show, summary, stride, sum, cumsum, sign, imag, conj, inv,
                                complex, reverse, exp, sqrt, abs, abs2, sign, issubset, values, in, first, last, rand, intersect, setdiff,
                                isless, union, angle, join, isnan, isapprox, isempty, sort, merge, promote_rule,
                                minimum, maximum, extrema, argmax, argmin, findmax, findmin, isfinite,
                                zeros, zero, one, promote_rule, repeat, length, resize!, isinf,
                                getproperty, findfirst, unsafe_getindex, fld, cld, div, real, imag,
                                @_inline_meta, eachindex, firstindex, lastindex, keys, isreal, OneTo,
                                Array, Vector, Matrix, view, ones, @propagate_inbounds, print_array,
                                split

import ApproxFunOrthogonalPolynomials: Laguerre

import LinearAlgebra: ⋅

export PeriodicLine, chop!, OscLaurent, SumFun, OscConstantSpace, zai, cai, condense, Cauchy, CauchyP, CauchyM, ⋅, fouriertransform, FourierTransform, combine!#, spacescompatible
#include("Domains/Domains.jl")

# Override the evaluation of Piecewise Funs
evaluate(f::Fun{T},x) where {T <: PiecewiseSpace{S}} where S = sum(map(F -> F(x),components(f)))

struct OscLaurent{D<:PeriodicLine,R} <: Space{D,R} # OscLaurent{D<:SPeriodicLine,R}?
    domain::D
    exp::Float64
    OscLaurent{D,R}(d,ex) where {D,R} = new{D,R}(d,ex)
    OscLaurent{D,R}(d) where {D,R} = new{D,R}(d,0.)
    OscLaurent{D,R}() where {D,R} = new{D,R}(D(),0.,0)
end
OscLaurent(d::PeriodicLine,exp::Float64) = OscLaurent{typeof(d),complex(prectype(d))}(d,exp)
OscLaurent(d::PeriodicLine) = OscLaurent(d,0.)
OscLaurent(α::Float64) = OscLaurent(PeriodicLine{false,Float64}(0.,1.),α)
OscLaurent(α::Float64,L::Float64) = OscLaurent(PeriodicLine{false,Float64}(0.,L),α)
OscLaurent() = OscLaurent(PeriodicLine())

include("SumFun.jl")

struct OscConstantSpace{D<:PeriodicLine,R} <: Space{D,R} # OscLaurent{D<:SPeriodicLine,R}?
    domain::D
    exp::Float64
    OscConstantSpace{D,R}(d,ex) where {D,R} = new{D,R}(d,ex)
    OscConstantSpace{D,R}(d) where {D,R} = new{D,R}(d,0.)
    OscConstantSpace{D,R}() where {D,R} = new{D,R}(D(),0.,0)
end
OscConstantSpace(d::PeriodicLine,exp::Float64) = OscConstantSpace{typeof(d),complex(prectype(d))}(d,exp)
OscConstantSpace(d::PeriodicLine) = OscConstantSpace(d,0.)
OscConstantSpace(α::Float64) = OscConstantSpace(PeriodicLine{false,Float64}(0.,1.),α)
OscConstantSpace(α::Float64,L::Float64) = OscConstantSpace(PeriodicLine{false,Float64}(0.,L),α)
OscConstantSpace() = OscConstantSpace(PeriodicLine())

maxspace(a::LaguerreWeight,b::LaguerreWeight) =  spacescompatible(a,b) ? a : PiecewiseSpace(a,b)
maxspace(a::OscLaurent,b::OscConstantSpace) = a
maxspace(b::OscConstantSpace,a::OscLaurent) = a

function getindex(C::ConcreteConversion{A,B,T},k::Integer,j::Integer) where {A<:OscConstantSpace,B<:OscLaurent,T}
    if j == k == 1
        return one(T)
    else
        return zero(T)
    end
end
Base.size(C::ConcreteConversion{A,B,T}) where {A<:OscConstantSpace,B<:OscLaurent,T} = (∞,1)
+(f::Fun{S},g::Fun{T}) where {S<:OscConstantSpace,T<:OscLaurent} = g+f
-(f::Fun{S},g::Fun{T}) where {S<:OscConstantSpace,T<:OscLaurent} = g-f
bandwidths(C::ConcreteConversion{A,B,T}) where {A<:OscConstantSpace,B<:OscLaurent,T} = (0,0)
Conversion(A::OscConstantSpace,B::OscLaurent) = ConcreteConversion(A::OscConstantSpace,B::OscLaurent)

function *(f::Fun{OscLaurent{D,R}},g::Fun{OscLaurent{D,R}})  where {D,R}
    sp = OscLaurent(f.space.domain,f.space.exp+g.space.exp)
    #m = maximum([length(f.coefficients),length(g.coefficients)])
    m = length(f.coefficients) + length(g.coefficients)
    Fun(sp,ApproxFun.transform(sp,values(pad(f,m)).*values(pad(g,m))))
end

evaluate(f::Fun{OscConstantSpace{D,R}},x::Float64) where {D,R} = f.coefficients[1]*exp(1im*f.space.exp*x)
function *(A::Array{T,2},b::Array{S,1})  where {T<:Fun,S<:Fun}
    c = [];
    for i = 1:size(A)[1]
        sum = A[i,1]*b[1]
        for j = 2:size(A)[2]
            sum = sum + A[i,j]*b[j]
        end
        append!(c,[sum])
    end
    c
end

function *(A::Array{T,2},b::Array{S,1})  where {T<:SumFun,S<:SumFun}
    c = [];
    for i = 1:size(A)[1]
        sum = A[i,1]*b[1]
        for j = 2:size(A)[2]
            sum = sum + A[i,j]*b[j]
        end
        append!(c,[sum])
    end
    c
end

spacescompatible(A::OscConstantSpace{D,R},B::OscLaurent{D,R})  where {D,R} = false
spacescompatible(B::OscLaurent{D,R},A::OscConstantSpace{D,R}) where {D,R} = A.exp ≈ B.exp
spacescompatible(B::OscConstantSpace{D,R},A::OscConstantSpace{D,R}) where {D,R} = A.exp ≈ B.exp
#hasconversion(A::OscConstantSpace,B::OscLaurent) = A.exp == B.exp
#hasconversion(B::OscLaurent,A::OscConstantSpace) = false

## A hack, definitely
## There is an issue because if F = Fun(f,OscLaurent(α)) then F != f, typically
## If checkpoints could depend on space and domain, this could be fixed easily
function Fun(f::Function,sp::OscLaurent{D,R}) where {D,R}
    g = Fun(f,Laurent(sp.domain))
    n = length(g.coefficients)
    n = round(Int64,n/length(g(points(g)[1])))+1
    return Fun(f,sp,n)
end

## Inner product for vectors.  Could extend to matrices using trace.
function ⋅(f::Fun{T},g::Fun{S}) where {S<:ApproxFun.ArraySpace{Q,1},T<:ApproxFun.ArraySpace{J,1}} where {Q,J}
    sum(map(⋅,f,g))
end

function ⋅(f::Fun{T},g::Fun{S},s::Function) where {S<:ApproxFun.ArraySpace{Q,1},T<:ApproxFun.ArraySpace{J,1}} where {Q,J}
    sum(s(transpose(conj(f))*g))[1]
end

function ⋅(f::Fun{T},g::Fun{S},s) where {S,T}
    sum(s(conj(f)*g))
end

function ⋅(f::Fun{T},g::Fun{S}) where {S,T}
    sum(conj(f)*g)
end

function ⋅(f::SumFun ,g::SumFun)
    sum(conj(f)*g)
end

#function ⋅(f::SumFun ,g::SumFun)
#    sum(conj(f)*g)
#end

function ⋅(f::Array{T,1},g::Array{S,1}) where {S<:SumFun,T<:SumFun}
    sum(sum(map(*,conj(f),g)))
end

function ⋅(f::Array{T,1},g::Array{S,1}) where {S<:Fun,T<:Fun}
    sum(map(⋅,f,g))
end

function condense(f::Fun)
    f
end

function condense(f::Array{T,1}) where T<:Fun
    map(condense,f)
end

function condense(f::Fun{T}) where {T <: SumSpace}
    sum(components(f))
end

function condense(f::Fun{T}) where {T <: ArraySpace{S}} where {S<:SumSpace}
    Fun(map( x -> sum(components(x)), Array(f)))
end

function condense(f::Fun{T}) where {T <: ArraySpace{S}} where {S<:Space}
    f
end

function chop(f::Array{T,1}) where T<:Fun
    map(chop,f)
end

function chop!(f::Array{T,1}) where T<:Fun
    map(chop!,f)
end

spacescompatible(a::OscLaurent{D,R},b::OscLaurent{D,R}) where {D,R} = a.exp ≈ b.exp
fourierpoints(n::Integer) = fourierpoints(Float64,n)
fourierpoints(::Type{T},n::Integer) where {T<:Number} = convert(T,π)*collect(0:2:2n-2)/n

##FFT That interlaces coefficients
zero_nan(x) = isnan(x) ? zero(typeof(x)) : x

function fsum(x::AbstractVector{T}) where T
    x[1] + sum(x[2:2:end].*((-1)^j for j in 1:length(x[2:2:end]))) + sum(x[3:2:end].*((-1)^j for j in 1:length(x[3:2:end])))
end

struct plan_mfft!
    P
end
plan_mfft!(x::AbstractVector{T}) where T = plan_mfft!(plan_fft!(x))

function *(tr::plan_mfft!,x::AbstractVector{T}) where T
    tr.P*x
    return x
end

struct plan_mifft!
    P
end
plan_mifft!(x::AbstractVector{T}) where T = plan_mifft!(plan_ifft!(x))

function *(tr::plan_mifft!,x::AbstractVector{T}) where T
    return tr.P*x
end

plan_transform!(sp::OscLaurent{D,R},x::AbstractVector{T}) where {D,R,T<:Complex} =
    TransformPlan(sp,plan_fft!(x),Val{true})
plan_itransform!(sp::OscLaurent{D,R},x::AbstractVector{T}) where {D,R,T<:Complex} =
    ITransformPlan(sp,plan_ifft!(x),Val{true})

plan_transform!(sp::OscLaurent{D,R},x::AbstractVector{T}) where {D,R,T} =
    error("In place variants not possible with real data.")
plan_itransform!(sp::OscLaurent{D,R},x::AbstractVector{T}) where {D,R,T} =
    error("In place variants not possible with real data.")

plan_transform(sp::OscLaurent{D,R},x::AbstractVector{T}) where {T<:Complex,D,R} =
    TransformPlan(sp,plan_transform!(sp,x),Val{false})
plan_itransform(sp::OscLaurent{D,R},x::AbstractVector{T}) where {T<:Complex,D,R} =
    ITransformPlan(sp,plan_itransform!(sp,x),Val{false})

function *(P::TransformPlan{T,OscLaurent{D,R},true},vals::AbstractVector{T}) where {T,D,R}
    n = length(vals)
    #vals = [zero_nan(j) for j in vals]
    vals = lmul!(inv(convert(T,n)),P.plan*vals)
    reverseeven!(interlace!(vals,1))
    vals
end

function *(P::ITransformPlan{T,OscLaurent{D,R},true},cfs::AbstractVector{T}) where {T,D,R}
    n = length(cfs)
    reverseeven!(cfs)
    cfs[:]=[cfs[1:2:end];cfs[2:2:end]]  # TODO: deinterlace!
    lmul!(n,cfs)
    P.plan*cfs
end

*(P::TransformPlan{T,OscLaurent{D,R},false},vals::AbstractVector{T}) where {T<:Complex,D,R} =
    P.plan*copy(vals)
*(P::TransformPlan{T,OscLaurent{D,R},false},vals::AbstractVector{T}) where {T,D,R} =
    P.plan*Vector{Complex{T}}(vals)
*(P::ITransformPlan{T,OscLaurent{D,R},false},vals::AbstractVector{T}) where {T<:Complex,D,R} =
    P.plan*copy(vals)
*(P::ITransformPlan{T,OscLaurent{D,R},false},vals::AbstractVector{T}) where {T,D,R} =
    P.plan*Vector{Complex{T}}(vals)

transform(::OscLaurent{D,R},vals,plan) where {D,R} = plan*vals
itransform(::OscLaurent{D,R},cfs,plan) where {D,R} = plan*cfs

transform(sp::OscLaurent{D,R},vals::AbstractVector) where {D,R} = plan_transform(sp,vals)*vals
itransform(sp::OscLaurent{D,R},cfs::AbstractVector) where {D,R} = plan_itransform(sp,cfs)*cfs

include("utils.jl")
include("calculus.jl")
include("operators.jl")
include("cauchy.jl")

function Base.conj(sp::OscLaurent{DD,RR}) where {DD,RR}
    OscLaurent(domain(sp),-sp.exp)
end

function Base.conj(f::Fun{OscLaurent{DD,RR}}) where {DD,RR}  ## Almost the same as Laurent
    ncoefficients(f) == 0 && return f

    cfs = Array{cfstype(f)}(undef, iseven(ncoefficients(f)) ? ncoefficients(f)+1 : ncoefficients(f))
    cfs[1] = conj(f.coefficients[1])
    cfs[ncoefficients(f)] = 0
    for k=2:2:ncoefficients(f)-1
        @inbounds cfs[k] = conj(f.coefficients[k+1])
    end
    for k=3:2:ncoefficients(f)+1
        @inbounds cfs[k] = conj(f.coefficients[k-1])
    end
    Fun(conj(space(f)),cfs)
end

end #module


text:
# This file is part of Kpax3. License is MIT.

# Pr(B1 = [1; 1; 1]) = 4 / 5
# Pr(B1 = [2; 1; 1]) = 1 / 5
#
# Pr(B2 = [1; 1; 1]) = (4 / 5) * (4 / 5) = 16 / 25
# Pr(B2 = [1; 2; 1]) = (4 / 5) * (1 / 5) = 4 / 25
# Pr(B2 = [2; 1; 1]) = (1 / 5) * (4 / 5) = 4 / 25
# Pr(B2 = [2; 2; 1]) = (1 / 5) * (1 / 5) * (1 / 2) = 1 / 50
# Pr(B2 = [2; 3; 1]) = (1 / 5) * (1 / 5) * (1 / 2) = 1 / 50
#
# Pr(B3 = [1; 1; 1]) = (4 / 5) * (4 / 5) * (4 / 5) = 64 / 125
# Pr(B3 = [1; 1; 2]) = (4 / 5) * (4 / 5) * (1 / 5) = 16 / 125
# Pr(B3 = [1; 2; 1]) = (4 / 5) * (1 / 5) * (4 / 5) = 16 / 125
# Pr(B3 = [1; 2; 2]) = (4 / 5) * (1 / 5) * (1 / 5) * (1 / 2) = 2 / 125
# Pr(B3 = [1; 2; 3]) = (4 / 5) * (1 / 5) * (1 / 5) * (1 / 2) = 2 / 125
# Pr(B3 = [2; 1; 1]) = (1 / 5) * (4 / 5) * (4 / 5) = 16 / 125
# Pr(B3 = [2; 1; 2]) = (1 / 5) * (4 / 5) * (1 / 5) * (1 / 2) = 2 / 125
# Pr(B3 = [2; 1; 3]) = (1 / 5) * (4 / 5) * (1 / 5) * (1 / 2) = 2 / 125
# Pr(B3 = [2; 2; 1]) = (1 / 5) * (1 / 5) * (1 / 2) * (4 / 5) = 2 / 125
# Pr(B3 = [2; 2; 2]) = (1 / 5) * (1 / 5) * (1 / 2) * (1 / 5) * (1 / 2) = 1/500
# Pr(B3 = [2; 2; 3]) = (1 / 5) * (1 / 5) * (1 / 2) * (1 / 5) * (1 / 2) = 1/500
# Pr(B3 = [2; 3; 1]) = (1 / 5) * (1 / 5) * (1 / 2) * (4 / 5) = 2 / 125
# Pr(B3 = [2; 3; 2]) = (1 / 5) * (1 / 5) * (1 / 2) * (1 / 5) * (1 / 2) = 1/500
# Pr(B3 = [2; 3; 3]) = (1 / 5) * (1 / 5) * (1 / 2) * (1 / 5) * (1 / 2) = 1/500

function test_mutation()
  R = [1; 1; 1]
  v = [3; 0; 0]

  o = Kpax3.KOffspring(copy(R), copy(v))

  mrate = 1 / 5

  N = 1000000

  pr = [64 / 125; 16 / 125; 16 / 125; 2 / 125; 2 / 125; 16 / 125; 2 / 125; 2 / 125;  2 / 125;  1 / 500; 1 / 500; 2 / 125;  1 / 500; 1 / 500]

  tmp = zeros(Float64, length(pr))

  for i in 1:N
    copyto!(o.R, R)
    copyto!(o.v, v)

    Kpax3.mutation!(o, mrate)

    if o.R == [1; 1; 1]
      @test o.v == [3; 0; 0]
      tmp[ 1] += 1
    elseif o.R == [1; 1; 2]
      @test o.v == [2; 1; 0]
      tmp[ 2] += 1
    elseif o.R == [1; 2; 1]
      @test o.v == [2; 1; 0]
      tmp[ 3] += 1
    elseif o.R == [1; 2; 2]
      @test o.v == [1; 2; 0]
      tmp[ 4] += 1
    elseif o.R == [1; 2; 3]
      @test o.v == [1; 1; 1]
      tmp[ 5] += 1
    elseif o.R == [2; 1; 1]
      @test o.v == [2; 1; 0]
      tmp[ 6] += 1
    elseif o.R == [2; 1; 2]
      @test o.v == [1; 2; 0]
      tmp[ 7] += 1
    elseif o.R == [2; 1; 3]
      @test o.v == [1; 1; 1]
      tmp[ 8] += 1
    elseif o.R == [2; 2; 1]
      @test o.v == [1; 2; 0]
      tmp[ 9] += 1
    elseif o.R == [2; 2; 2]
      @test o.v == [0; 3; 0]
      tmp[10] += 1
    elseif o.R == [2; 2; 3]
      @test o.v == [0; 2; 1]
      tmp[11] += 1
    elseif o.R == [2; 3; 1]
      @test o.v == [1; 1; 1]
      tmp[12] += 1
    elseif o.R == [2; 3; 2]
      @test o.v == [0; 2; 1]
      tmp[13] += 1
    elseif o.R == [2; 3; 3]
      @test o.v == [0; 1; 2]
      tmp[14] += 1
    end
  end

  tmp /= N

  @test isapprox(tmp, pr, rtol=0.01)

  nothing
end

test_mutation()


text:
# DS_DN.jl
# ========
#
# Parameterized procedure for identifying and calculating the number of
# synonymous and nonsynonymous mutations between codons, accounting for the
# different possible pathways of substitution.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/NaturalSelection.jl/blob/master/LICENSE.md

abstract type DS_DN_Counter end

const ONEPATH = Val{1}()
const TWOPATHS = Val{2}()
const THREEPATHS = Val{3}()

function splice_into(x::C, y::C, pos::Integer) where C <: Codon
    mask = UInt64(3) << bitindex(x, pos)
    return C((UInt64(x) & ~mask) | (UInt64(y) & mask))
end

function find_differences(x::Codon, y::Codon)
    diffs = 0x00
    @inbounds for pos in 1:3
        diffs = (x[pos] != y[pos]) | (diffs << 1)
    end
    return diffs, count_ones(diffs)
end

@inline function classify_mutation(x::Codon, y::Codon, code::GeneticCode, weight::N) where N <: Real
    @inbounds aresame = code[x] == code[y]
    DS = ifelse(aresame, weight, N(0))
    DN = ifelse(aresame, N(0), weight)
    return DS, DN
end

"""
    DS_DN_enumerator(::Type{T}, x::Codon, y::Codon, code::GeneticCode = DEFAULT_TRANS) where T<:DS_DN_Counter

This generic, parametric algorithm takes two codons, and computes the number of
synonymous and non-synonymous mutations between them.

The exact behaviour of the algorithm can be adjusted by defining different
types of T <: DS_DN_Counter.

Different counter types allows one to change how different types of
mutations are weighted, and how the results are accumulated.

For example, in this package there are currently two DS_DN_Counter types.

One is called `AllPaths`, and the other is called `ShortestPath`.

If one site between the two codons is different. It is obvious and simple
then to count whether it is a synonymous or nonsynonymous mutation.
In this scenario, both `AllPaths` and `ShortestPath` count the mutation(s), and
both counters weight them in the same way. Both counters act the same in this
scenario.

However, when there are two differences between the two codons, the number of
synonymous and non-synonymous differences per codon, sum to 2, and there are two
possible mutation pathways that could have happened, each possible pathway has
two steps.
For example, comparing CTA and GTT, the possible pathways are:

```
1: CTA (L) -> GTA (V) -> GTT (V) : 1 nonsynonymous change and 1 synonymous change.
2: CTA (L) -> CTT (L) -> GTT (V) : 1 nonsynonymous change and 1 synonymous change.
```

Furthuremore, for three changes, there are 6 possible mutation pathways that
could have happened, each with three steps.

For example, comparing CTA and GAT, the possible pathways are:

```
1: CTA (L) -> GTA (V) -> GAA (E) -> GAT (D) : 3 nonsynonymous changes.
2: CTA (L) -> GTA (V) -> GTT (V) -> GAT (D) : 2 nonsynonymous and 1 synonymous change.
3: CTA (L) -> CAA (Q) -> GAA (E) -> GAT (D) : 3 nonsynonymous changes.
4: CTA (L) -> CAA (Q) -> CAT (H) -> GAT (D) : 3 nonsynonymous changes.
5: CTA (L) -> CTT (L) -> GTT (V) -> GAT (D) : 2 nonsynonymous changes and 1 synonymous change.
6: CTA (L) -> CTT (L) -> CAT (H) -> GAT (D) : 2 nonsynonymous changes and 1 synonymous change.
```

In these two cases, the two counters differ. `AllPaths` considers all
paths equally likely, and weights them all equally. As it iterates through each
path it accumulates the weighted number of synonymous and non-synonymous subsitutions
from each path returning the accumulated total.

In contrast, the `ShortestPath` also counts the number of synonymous and
non-synonymous mutations for each path, but it only keeps, and finally returns
the number of synonymous and non-synonymous mutations for the most parsimonious
path.

The `AllPaths` counter Counts DS and DN according to the all-paths method
devised for use with Nei and Gojoborei's method of estimating dN and dS ratios.
It considers all paths equally likely and accumulates DN and DS as it considers
each path.

The `ShortestPath` method does not consider each path equally likely, it returns
the DN and DS for the most parsiomonious path out of all paths it considered.
"""
function DS_DN_enumerator(::Type{T}, x::Codon, y::Codon, code::GeneticCode = DEFAULT_TRANS) where T <: DS_DN_Counter
    if x == y # Early escape, codons are the same, no syn or nonsyn mutations.
        return 0.0, 0.0
    else
        diff_positions, n_diffs = find_differences(x, y) # Which positions are different.
        if n_diffs == 1
            DS, DN = classify_mutation(x, y, code, weighting(T, ONEPATH))
            # One site in the two codons is different. It is obvious and simple
            # then to count whether it is a synonymous or nonsynonymous mutation.
        elseif n_diffs == 2
            DS, DN = DS_DN_init(T)
            # For two changes, the number of synonymous and non-synonymous
            # differences per codon, sum to 2, there are two pathways,
            # each possible pathway having two steps.
            # For example, comparing CTA and GTT, the possible pathways are:
            # 1: CTA (L) -> GTA (V) -> GTT (V) : 1 nonsynonymous change and 1 synonymous change.
            # 2: CTA (L) -> CTT (L) -> GTT (V) : 1 nonsynonymous change and 1 synonymous change.
            @inbounds for pos in 1:3
                if ((diff_positions >> (3 - pos)) & 0x01) == 0x01
                    temp_cdn = splice_into(x, y, pos)

                    DS_a, DN_a = classify_mutation(x, temp_cdn, code, weighting(T, TWOPATHS))
                    DS_b, DN_b = classify_mutation(temp_cdn, y, code, weighting(T, TWOPATHS))
                    DS_i = DS_a + DS_b
                    DN_i = DN_a + DN_b

                    DS, DN = accumulate(T, DS, DN, DS_i, DN_i)

                end
            end
        elseif n_diffs == 3
            DS, DN = DS_DN_init(T)
            # For three changes, there are 6 pathways, each with three steps.
            # For example, comparing CTA and GAT, the possible pathways are:
            # 1: CTA (L) -> GTA (V) -> GAA (E) -> GAT (D) : 3 nonsynonymous changes.
            # 2: CTA (L) -> GTA (V) -> GTT (V) -> GAT (D) : 2 nonsynonymous and 1 synonymous change.
            # 3: CTA (L) -> CAA (Q) -> GAA (E) -> GAT (D) : 3 nonsynonymous changes.
            # 4: CTA (L) -> CAA (Q) -> CAT (H) -> GAT (D) : 3 nonsynonymous changes.
            # 5: CTA (L) -> CTT (L) -> GTT (V) -> GAT (D) : 2 nonsynonymous changes and 1 synonymous change.
            # 6: CTA (L) -> CTT (L) -> CAT (H) -> GAT (D) : 2 nonsynonymous changes and 1 synonymous change.
            @inbounds for path in SITE_PERMUTATIONS
                tmp_cdn_a = splice_into(x, y, path[1])
                tmp_cdn_b = splice_into(tmp_cdn_a, y, path[2])

                DS_a, DN_a = classify_mutation(x, tmp_cdn_a, code, weighting(T, THREEPATHS))
                DS_b, DN_b = classify_mutation(tmp_cdn_a, tmp_cdn_b, code, weighting(T, THREEPATHS))
                DS_c, DN_c = classify_mutation(tmp_cdn_b, y, code, weighting(T, THREEPATHS))
                DS_i = DS_a + DS_b + DS_c
                DN_i = DN_a + DN_b + DN_c

                DS, DN = accumulate(T, DS, DN, DS_i, DN_i)

            end
        end
        return DS, DN
    end
end


text:
const _BitInteger = (Int == Int64 ? Base.BitInteger64 : Base.BitInteger32)
const _BitSigned = (Int == Int64 ? Base.BitSigned64 : Base.BitSigned32)
const _BitUnsigned = (Int == Int64 ? Base.BitUnsigned64 : Base.BitUnsigned32)

### Mag
for (jf, af) in [(:+, :add!), (:-, :sub!), (:*, :mul!), (:/, :div!)]
    @eval Base.$jf(x::MagOrRef, y::MagOrRef) = $af(zero(x), x, y)
end
Base.:+(x::MagOrRef, y::Integer) = add!(zero(x), x, convert(UInt, y))
Base.:+(x::Integer, y::MagOrRef) = add!(zero(y), y, convert(UInt, x))
Base.:*(x::MagOrRef, y::Integer) = mul!(zero(x), x, convert(UInt, y))
Base.:*(x::Integer, y::MagOrRef) = mul!(zero(y), y, convert(UInt, x))
Base.:/(x::MagOrRef, y::Integer) = div!(zero(x), x, convert(UInt, y))

Base.:(^)(x::MagOrRef, e::Integer) = pow!(zero(x), x, convert(UInt, e))
rsqrt(x::MagOrRef) = rsqrt!(zero(x), x)
Base.hypot(x::MagOrRef, y::MagOrRef) = hypot!(zero(x), x, y)
root(x::MagOrRef, n::Integer) = root!(zero(x), x, convert(UInt, n))
neglog(x::MagOrRef) = neg_log!(zero(x), x)
expinv(x::MagOrRef) = expinv!(zero(x), x)
for f in [:inv, :sqrt, :log, :log1p, :exp, :expm1, :atan, :cosh, :sinh]
    @eval Base.$f(x::MagOrRef) = $(Symbol(f, :!))(zero(x), x)
end

### Arf
function Base.sign(x::ArfOrRef)
    isnan(x) && return Arf(NaN) # Follow Julia and return NaN
    return Arf(sgn(x))
end

Base.abs(x::ArfOrRef) = abs!(zero(x), x)
Base.:(-)(x::ArfOrRef) = neg!(zero(x), x)
for (jf, af) in [(:+, :add!), (:-, :sub!), (:*, :mul!), (:/, :div!)]
    @eval function Base.$jf(x::ArfOrRef, y::Union{ArfOrRef,_BitInteger})
        z = Arf(prec = _precision((x, y)))
        $af(z, x, y)
        return z
    end
end
function Base.:+(x::_BitInteger, y::ArfOrRef)
    z = zero(y)
    add!(z, y, x)
    return z
end
function Base.:*(x::_BitInteger, y::ArfOrRef)
    z = zero(y)
    mul!(z, y, x)
    return z
end
function Base.:/(x::_BitUnsigned, y::ArfOrRef)
    z = zero(y)
    ui_div!(z, x, y)
    return z
end
function Base.:/(x::_BitSigned, y::ArfOrRef)
    z = zero(y)
    si_div!(z, x, y)
    return z
end

function Base.sqrt(x::ArfOrRef)
    y = zero(x)
    sqrt!(y, x)
    return y
end
function rsqrt(x::ArfOrRef)
    y = zero(x)
    rsqrt!(y, x)
    return y
end
function root(x::ArfOrRef, k::Integer)
    y = zero(x)
    root!(y, x, convert(UInt, k))
    return y
end

### Arb and Acb
for (jf, af) in [(:+, :add!), (:-, :sub!), (:*, :mul!), (:/, :div!)]
    @eval Base.$jf(x::ArbOrRef, y::Union{ArbOrRef,ArfOrRef,_BitInteger}) =
        $af(Arb(prec = _precision((x, y))), x, y)
    @eval Base.$jf(x::AcbOrRef, y::Union{AcbOrRef,ArbOrRef,_BitInteger}) =
        $af(Acb(prec = _precision((x, y))), x, y)
    if jf == :(+) || jf == :(*)
        @eval Base.$jf(x::Union{ArfOrRef,_BitInteger}, y::ArbOrRef) =
            $af(Arb(prec = _precision((x, y))), y, x)
        @eval Base.$jf(x::Union{ArbOrRef,_BitInteger}, y::AcbOrRef) =
            $af(Acb(prec = _precision((x, y))), y, x)
    end
end

Base.:(-)(x::Union{ArbOrRef,AcbOrRef}) = neg!(zero(x), x)
Base.abs(x::ArbOrRef) = abs!(zero(x), x)
Base.:(/)(x::_BitUnsigned, y::ArbOrRef) = ui_div!(zero(y), x, y)

Base.:(^)(x::ArbOrRef, y::ArbOrRef) = pow!(Arb(prec = _precision((x, y))), x, y)
function Base.:(^)(x::ArbOrRef, y::_BitInteger)
    z = zero(x)
    x, n = (y >= 0 ? (x, y) : (inv!(z, x), -y))
    return pow!(z, x, convert(UInt, n))
end
Base.:(^)(x::AcbOrRef, y::Union{AcbOrRef,ArbOrRef,_BitInteger}) =
    pow!(Acb(prec = _precision((x, y))), x, y)

# We define the same special cases as Arb does, this avoids some
# overhead
Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{-2}) =
    (y = inv(x); sqr!(y, y))
#Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{-1}) - implemented in Base.intfuncs.jl
Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{0}) = one(x)
Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{1}) = copy(x)
Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{2}) = sqr(x)

Base.hypot(x::ArbOrRef, y::ArbOrRef) = hypot!(Arb(prec = _precision((x, y))), x, y)

root(x::Union{ArbOrRef,AcbOrRef}, k::Integer) = root!(zero(x), x, convert(UInt, k))

# Unary methods in Base
for f in [
    :inv,
    :sqrt,
    :log,
    :log1p,
    :exp,
    :expm1,
    :sin,
    :cos,
    :tan,
    :cot,
    :sec,
    :csc,
    :atan,
    :asin,
    :acos,
    :sinh,
    :cosh,
    :tanh,
    :coth,
    :sech,
    :csch,
    :atanh,
    :asinh,
    :acosh,
]
    @eval Base.$f(x::Union{ArbOrRef,AcbOrRef}) = $(Symbol(f, :!))(zero(x), x)
end

sqrtpos(x::ArbOrRef) = sqrtpos!(zero(x), x)
sqrt1pm1(x::ArbOrRef) = sqrt1pm1!(zero(x), x)
rsqrt(x::Union{ArbOrRef,AcbOrRef}) = rsqrt!(zero(x), x)
sqr(x::Union{ArbOrRef,AcbOrRef}) = sqr!(zero(x), x)

Base.sinpi(x::Union{ArbOrRef,AcbOrRef}) = sin_pi!(zero(x), x)
Base.cospi(x::Union{ArbOrRef,AcbOrRef}) = cos_pi!(zero(x), x)
tanpi(x::Union{ArbOrRef,AcbOrRef}) = tan_pi!(zero(x), x)
cotpi(x::Union{ArbOrRef,AcbOrRef}) = cot_pi!(zero(x), x)
cscpi(x::Union{ArbOrRef,AcbOrRef}) = csc_pi!(zero(x), x)
# Julias definition of sinc is equivalent to Arbs definition of sincpi
Base.sinc(x::Union{ArbOrRef,AcbOrRef}) = sinc_pi!(zero(x), x)
Base.atan(y::ArbOrRef, x::ArbOrRef) = atan2!(Arb(prec = _precision((y, x))), y, x)

function Base.sincos(x::Union{ArbOrRef,AcbOrRef})
    s, c = zero(x), zero(x)
    sin_cos!(s, c, x)
    return (s, c)
end
function sincospi(x::Union{ArbOrRef,AcbOrRef})
    s, c = zero(x), zero(x)
    sin_cos_pi!(s, c, x)
    return (s, c)
end
function sinhcosh(x::Union{ArbOrRef,AcbOrRef})
    s, c = zero(x), zero(x)
    sinh_cosh!(s, c, x)
    return (s, c)
end

### Acb
function Base.:(*)(x::AcbOrRef, y::Complex{Bool})
    if real(y)
        if imag(y)
            z = mul_onei!(zero(x), x)
            return add!(z, x, z)
        else
            return Acb(x)
        end
    end
    imag(y) && return mul_onei!(zero(x), x)
    return zero(x)
end
Base.:(*)(x::Complex{Bool}, y::AcbOrRef) = y * x

Base.real(z::AcbLike; prec = _precision(z)) = get_real!(Arb(prec = prec), z)
Base.imag(z::AcbLike; prec = _precision(z)) = get_imag!(Arb(prec = prec), z)
Base.conj(z::AcbLike) = conj!(Acb(prec = _precision(z)), z)
Base.abs(z::AcbLike) = abs!(Arb(prec = _precision(z)), z)

### min and max
for T in [MagOrRef, ArfOrRef, ArbOrRef]
    for op in [:min, :max]
        # min/max
        @eval Base.$op(x::$T, y::$T) = $(Symbol(op, :!))(zero(x), x, y)
        # minimum/maximum
        # The default implementation of minimum/maximum in Julia is
        # not correct for Arb, we implemented a correct version when
        # no specific dimension is given. The default give correct
        # results for Mag and Arf but wrong results for Arb in some
        # cases, be careful!
        @eval function Base.$(Symbol(op, :imum))(A::AbstractArray{<:$T})
            isempty(A) &&
                throw(ArgumentError("reducing over an empty collection is not allowed"))
            res = copy(first(A))
            for x in Iterators.drop(A, 1)
                $(Symbol(op, :!))(res, res, x)
            end
            return res
        end
    end
    @eval Base.minmax(x::$T, y::$T) = (min(x, y), max(x, y))

    @eval function Base.extrema(A::AbstractArray{<:$T})
        isempty(A) &&
            throw(ArgumentError("reducing over an empty collection is not allowed"))
        l = copy(first(A))
        u = copy(first(A))
        for x in Iterators.drop(A, 1)
            min!(l, l, x)
            max!(u, u, x)
        end
        return (l, u)
    end
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

mutable struct TypeDBClientException <: Exception
    error_message::Optional{AbstractGeneralError}
    params::Optional{Tuple}
    individual_message::Optional{String}
    cause::Optional{Exception}
end

function Base.show(io::IO, typedb_excption::TypeDBClientException)
    err_message = string(typedb_excption.error_message)
    joined_params = ""
    if typedb_excption.params !== nothing
        if !isempty(typedb_excption.params)
            replace_item = "_error_item"
            err_message = replace(err_message, replace_item=>string(typedb_excption.params[1]))
            joined_params = join(typedb_excption.params, "\n")
        end
    end
    println(io,"TypeDBClientException:")
    println(io,"message: ", err_message)
    println(io,"params: ",joined_params)
    println(io,"remark: ",typedb_excption.individual_message)
    print(io,"cause: ")
    println(io,typedb_excption.cause)

    return nothing
end

function TypeDBClientException(err::Type{<:AbstractGeneralError}, parameters...)
    def_err = _build_error_messages(err)
    _params = isempty(parameters) ? Tuple{}() : parameters
    return TypeDBClientException(def_err, _params, nothing, nothing)
end

function TypeDBClientException(err::gRPCServiceCallException, parameters...)
    def_err = _build_error_messages(GRPC_SERVER_ERROR)

    # reading the Number and Prefix from the Server message.
    bracket_range = findfirst(r"\[.+\]", err.message)
    if bracket_range !== nothing
        def_err.message_prefix = bracket_value = err.message[bracket_range]
        (range_code_prefix = findfirst(r"[A-Z]+", bracket_value)) !== nothing &&
	    (def_err.code_prefix = bracket_value[range_code_prefix])
        (range_code_number = findfirst(r"\d+", bracket_value)) !== nothing &&
            (def_err.code_number = parse(Int64, bracket_value[range_code_number]))
	    def_err.message_body = string(strip(err.message[bracket_range.stop+1:end]))
    else
        # determine the default values to detect the error message
        # if there is no error number from the server
	    def_err.message_prefix = ""
        def_err.message_body = string(strip(err.message))
    end

    return TypeDBClientException(def_err, isempty(parameters) ? Tuple{}() : parameters, nothing, nothing)
end

function TypeDBClientException(message::AbstractString, cause::Exception)
    err = _build_error_messages(GENERAL_UNKOWN_ERROR)
    return TypeDBClientException(err,nothing,message,cause)
end


text:
using GLVisualize, GeometryTypes, GLAbstraction, Colors

if !isdefined(:runtests)
	window = glscreen()
end

description = """
Example that will walk through the basic steps to create a mesh with colors from
Geometry primitives.
"""

baselen = 0.4f0
dirlen = 2f0
# create an array of differently colored boxes in the direction of the 3 axes
rectangles = [
    (HyperRectangle{3,Float32}(Vec3f0(baselen), Vec3f0(dirlen, baselen, baselen)), RGBA(1f0,0f0,0f0,1f0)),
    (HyperRectangle{3,Float32}(Vec3f0(baselen), Vec3f0(baselen, dirlen, baselen)), RGBA(0f0,1f0,0f0,1f0)),
    (HyperRectangle{3,Float32}(Vec3f0(baselen), Vec3f0(baselen, baselen, dirlen)), RGBA(0f0,0f0,1f0,1f0))
]
# convert to an array of normal meshes
# note, that the constructor is a bit weird. GLNormalMesh takes a tuple of
# a geometry and a color. This means, the geometry will be converted to a GLNormalMesh
# and the color will be added afterwards, so the resulting type is a GLNormalColorMesh
meshes = map(GLNormalMesh, rectangles)
# merge them into one big mesh
# the resulting type is a GLNormalAttributeMesh, since we merged meshes with different
# attributes (colors). An array of the colors will be created and each vertex in the
# mesh will be asigned to one of the colors found there.
colored_mesh = merge(meshes)
_view(visualize(colored_mesh), window)

# one could also create a GLNormalAttributeMesh manually:
sphere = Sphere{Float32}(Point3f0(0), 2f0)
# decompose decomposes a mesh or geometry into the primitive of the first argument
vertices = decompose(Point3f0, sphere, 50)
faces = decompose(GLTriangle, sphere, 50)

# create a few colors, can be N0f8 or Float32
colors = RGBA{N0f8}[RGBA{N0f8}(rand(), rand(), rand(), 1.) for i=1:5]
# assign every vertice a random index into the color array (0-based indexes)
attribute_id = rand(0f0:4f0, length(vertices))

sphere_mesh = GLNormalAttributeMesh(
    vertices=vertices, faces=faces,
    attributes=colors, attribute_id=attribute_id
)

# move the model a bit to the right
moveright = translationmatrix(Vec3f0(5,0,0))

_view(visualize(sphere_mesh, model=moveright), window)

if !isdefined(:runtests)
	renderloop(window)
end


text:
using Discourse

function convert_HTTP_Response_To_JSON(response)

    compressed = HTTP.payload(response);
    decompressed = transcode(GzipDecompressor, compressed);
    json = JSON.parse(IOBuffer(decompressed))
    return json

end

function _get(path, kwargs...)
    """
        HTTP Helper function
        ARGS:
            path:
                kwargs:
        RETURNS:
    """
    return _request(GET, path, params=kwargs)

end

function _request(verb, path, params::Dict, files::Dict, data::Dict, json::Dict)
    """
        Executes HTTP request to API and handles response
        ARGS:
            verb: HTTP verb as string; defined in Discourse.jl
            path: path on Discourse API
            params: dict of parameters to include in API
        RETURNS:
            response_body_data::Dictionary / None::None
    """
    url = DiscourseClient.host + path
    headers = Dict(
        "Accept" => "application/json; charset=utf-8",
        "Api-key" => DiscourseClient.api_key,
        "Api-username" => DiscourseClient.api_username
    )

    # times we should retry request if rate limited.
    retry_count = 4
    # Extra time (on top of that required by API) to wait on a retry
    retry_backoff = 1

    while retry_count>=0
        r = HTTP.request(
            verb,
            url,
            allow_redirects=false,
            params=params,
            files=files,
            data=data,
            json=json,
            headers=headers,
            timeout=DiscourseClient.timeout
        )
        if r.status == 200
            break
        end
    end
    json_content = "application/json; charset=utf-8"
    content_type = r.headers["content_type"]
    decoded = r.json()
    return decoded
end


text:
@info "Caesar.jl is loading tools related to Distributed.jl"

export plotToVideo_Distributed

function plotToVideo_Distributed( pltfnc::Function,
                                  args::Union{<:AbstractVector, <:Tuple};
                                  pool = WorkerPool(procs()[2:end]),
                                  videofile::AbstractString="/tmp/test.avi" )
  # prep locations for image results
  len = length(args)
  _PL = Vector{Any}(undef, len)
  #
  pool = WorkerPool(procs()[2:end])
  for (i,lb) in enumerate(args)
    # Threads.@spawn begin
      im = pltfnc(lb)
      _PL[i] = remotecall(convert, pool, Matrix{RGB}, im)
      @show i
      nothing
    # end
  end
  # fetch the remote results
  PL = fetch.(_PL)

  Caesar.writevideo(videofile, PL, fps=5)

  PL, videofile
end


text:
include("loadrealdata.jl")
include("run.jl")
using DataFrames

F = typeof(1.0)

activations = [
    Relu{F},
    Requ{F}
]

constraintsets = [
    L2Constraint,
    BasisL2Constraint
]

Cs = [1e-3, 1, 2, 4, 8, 12, 16, 24, 32, 48, 64, 100]

algconf(algo::Symbol, args::Tuple, kw::Dict) = Dict(:algo=>algo, :args=>args, :kwargs=>kw)
algconfigs = [algconf(a...) for a in (
    (:fw, (Cs,), Dict(:ad=>false)),
    (:fw, (Cs,), Dict(:ad=>true)),

    (:gb, (1.,), Dict(:ad=>false, :maxiter=>800)),
    (:gb, (1.,), Dict(:ad=>true, :maxiter=>800)),

    (:sw, (1/16.,), Dict(:ad=>false, :maxiter=>1000)),
    (:sw, (1/16.,), Dict(:ad=>true, :maxiter=>1000)),

    (:sw, (1/16.,), Dict(:lasso=>true, :maxiter=>1500)),
    (:sw, (1/16.,), Dict(:lasso=>true, :ad=>true, :maxiter=>1500)),

    (:nn, (Cs, 300), Dict(:stepsize=>3e-3, :maxiter=>5000)),
    (:xg, (), Dict(:numround=>500)),
)]

function getarg(idx, typ::Type{T}=Int, default::T=1) where T
    if length(ARGS) ≥ idx
        return parse(ARGS[idx])::T
    else
        return default
    end
end

dataidx = getarg(1, Symbol, :toy1)
algidx = getarg(2)
actidx  = getarg(4)
considx = getarg(3)
if length(ARGS) == 0
    println("Usage:\n\ttest.jl [dataset_idx=:toy1] [algorithm_idx=1] [activation_idx=1] [constraint_set_idx=1]")
    println("Using defaults...")
end

trdata, tedata = loadrealdata(dataidx)
l = all(tedata.y.^2 .== tedata.y) ? LOGIT : MSE
B = activations[actidx]
O = constraintsets[considx]
algconfig = algconfigs[algidx]

##
println("Loss: ", l)
println("Activation: ", B)
println("Constraint Set: ", O)
println("Algo: ", algconfig)
println("Data: ", dataidx)
##

r = onerun(trdata, tedata, l, B, O; algconfig...)
@show DataFrame(Iteration=1:length(r[1]),
               Cardinality=r[1], Norm=r[2],
               Trainerr=r[3], Testerr=r[4], Bayes=r[5],
               Altrainerr=r[6], Altesterr=r[7], AltBayes=r[8],
               Time=r[9])


text:
"""
Common exceptions
https://github.com/frictionlessdata/tableschema-jl#exceptions
"""
struct ConstraintError <: Exception
    message::String
    field::Field
    value
    expected

    ConstraintError(m::String, f::Field, v, e) = new(m, f, v, e)
    ConstraintError(m::String, f::Field, v) = new(m, f, v, nothing)
    ConstraintError(m::String, v, e) = new(m, nothing, v, e)
    ConstraintError(m::String, v) = new(m, nothing, v, nothing)
end

struct FieldError <: Exception
    message::String
    # key::String
    # line::Int16
end

struct SchemaError <: Exception
    message::String
    key

    SchemaError(m::String, k::String) = new(m, k)
    SchemaError(m::String) = new(m, nothing)
    SchemaError(f::FieldError) = new(f.message, nothing)
end

struct TableValidationException <: Exception
    var::String
end

struct CastError <: Exception
    message::String
    errors::Array{CastError}

    CastError(m::String, e::Array) = new(m, e)
    CastError(m::String) = new(m, [])
    CastError() = new("Cast error", [])
end


text:
function get_initial_population(
        model::Model,
        nth_param_set::Int64,
        n_population::Int64,
        n_gene::Int64,
        initial_threshold::Float64)::Matrix{Float64}
    open(
        joinpath(
            model.path,
            "logs",
            "$nth_param_set.log"
        ),
        "w",
    ) do f
        write(f, "Generating the initial population...\n\n")
    end
    population::Matrix{Float64} = fill(
        Inf, (n_population, n_gene + 1)
    )
    @inbounds @simd for i = 1:n_population
        while initial_threshold <= population[i, end]
            for j = 1:n_gene
                population[i, j] = rand()
            end
            population[i, end] = model.obj_func(population[i, 1:n_gene])
        end
        open(
            joinpath(
                model.path,
                "logs",
                "$nth_param_set.log"
            ),
            "a",
        ) do f
            write(f, @sprintf("%d / %d\n", i, n_population))
        end
    end
    open(
        joinpath(
            model.path,
            "logs",
            "$nth_param_set.log"
        ),
        "a"
    ) do f
        write(f, "\n----------------------------------------\n")
    end
    population = sortslices(population, dims=1, by=x -> x[end])

    return population
end

function get_initial_population_continue(
        model::Model,
        nth_param_set::Int64,
        n_population::Int64,
        n_gene::Int64,
        initial_threshold::Float64,
        p0_bounds::Vector{Float64})::Matrix{Float64}
    generation::Int64 = readdlm(
        joinpath(
            model.path,
            "fitparam",
            "$nth_param_set",
            "generation.dat"
        )
    )[1, 1]
    best_indiv::Vector{Float64} = readdlm(
        joinpath(
            model.path,
            "fitparam",
            "$nth_param_set",
            "fit_param$generation.dat"
        )
    )[:, 1]
    open(
        joinpath(
            model.path,
            "logs",
            "$nth_param_set.log"
        ),
        "a"
    ) do f
        write(f,
            "\n----------------------------------------\n" *
            "Generating the initial population...\n"
        )
    end
    population::Matrix{Float64} = fill(
        Inf, (n_population, n_gene + 1)
    )
    @inbounds @simd for i = 1:n_population
        while initial_threshold <= population[i, end]
            for gene_idx = 1:n_gene
                population[i, gene_idx] = model.bestIndivVal2randGene(
                    gene_idx, best_indiv, p0_bounds
                )
                if population[i, gene_idx] > 1.0
                    population[i, gene_idx] = 1.0
                elseif population[i, gene_idx] < 0.0
                    population[i, gene_idx] = 0.0
                end
            end
            population[i, end] = model.obj_func(population[i, 1:n_gene])
        end
        open(
            joinpath(
                model.path,
                "logs",
                "$nth_param_set.log"
            ),
            "a"
        ) do f
            write(f, @sprintf("%d / %d\n", i, n_population))
        end
    end
    open(
        joinpath(
            model.path,
            "logs",
            "$nth_param_set.log"
        ),
        "a"
    ) do f
        write(f, "\n----------------------------------------\n")
    end
    population = sortslices(population, dims=1, by=x -> x[end])

    return population
end


text:
# PolynomialApproximation.jl - polynomial approximations to FunctionalForm with ApproxFun package

using ApproxFun

# specify range of function. Too short + you end up with millions of polynomial coeffs
OHVssσFun=Fun(r->OHVssσ(r), 0.4..10)

# count number of coeffs
ncoefficients(OHVssσFun)

# plot the risidual
#plot( [x -> OHVssσFun(x) - OHVssσ(x)], 0.1,10 )

# derivative of the polynomial fit
OHVssσFunDeriv=OHVssσFun'

# Shouldn't really be doing computation in the src... Oh well!

println("Testing ApproxFun...")

orig=OOPP # OOPP is an EPL

for rcut in 0.05:0.1:0.5
    origFun=Fun(r->orig(r), rcut..10)
    coeff=ncoefficients(origFun)
    resid=origFun(1.0)-orig(1.0)
    println("rcut: $rcut ncoeff OHVssσFun=Fun: $coeff Residual@1.0: $resid")
    print("ApproxFun: ")
    @time map(origFun,0.1:0.01:10)
    print("Original: ")
    @time map(orig,0.1:0.01:10)
end


text:
function grab(from::T, p::Symbol, sample_idx::Symbol) where T
  return getfield(from, Symbol(p, sample_idx))
end

function update_etaC(m::Gtilde, s::T) where T
  anew = copy(collect(m.eta.alpha))

  foreach(lam -> anew[lam] +=1, s.lambdaC)
  m.beta || foreach(lam -> anew[lam] +=1, s.lambdaT)

  return rand(Dirichlet(anew))
end

function update_etaT(m::Gtilde, s::T) where T
  if m.beta
    anew = copy(collect(m.eta.alpha))
    foreach(lam -> anew[lam] +=1, s.lambdaT)
    return rand(Dirichlet(anew))
  else
    return s.etaC
  end
end

function update_lambda(m::Gtilde, s::T, i::Symbol) where T
  y = grab(m, :y, i)
  v = grab(s, :v, i)
  zeta = grab(s, :zeta, i)
  eta = grab(s, :eta, i)
  logeta = log.(eta)

  lambda = [let
              loc = s.mu + s.psi * zeta[n]
              scale = sqrt.(s.omega) / sqrt(v[n])
              logmix = normlogpdf.(loc, scale, y[n]) + logeta
              all(s.nu .< 10000) && (logmix += gammalogpdf.(s.nu/2, 2 ./ s.nu, v[n]))
              MCMC.wsample_logprob(logmix)
            end for n in eachindex(y)]

  return lambda
end

update_lambdaC(m::Gtilde, s::T) where T = update_lambda(m, s, :C)
update_lambdaT(m::Gtilde, s::T) where T = update_lambda(m, s, :T)

function update_v(m::Gtilde, s::T, i::Symbol) where T
  lami = grab(s, :lambda, i)
  zetai = grab(s, :zeta, i)
  yi = grab(m, :y, i)
  vi = grab(s, :v, i)

  nu = s.nu[lami]
  omega = s.omega[lami]
  psi = s.psi[lami]
  mu = s.mu[lami]

  shape = nu/2 .+ 1
  rate = (nu + zetai.^2 + ((yi - mu - psi .* zetai) .^ 2) ./ omega) / 2

  return @.rand(Gamma(shape, 1 / rate))
end

update_vC(m::Gtilde, s::T) where T = update_v(m, s, :C)
update_vT(m::Gtilde, s::T) where T = update_v(m, s, :T)

function update_zeta(m::Gtilde, s::T, i::Symbol) where T
  yi = grab(m, :y, i)
  zetai = grab(s, :zeta, i)
  lami = grab(s, :lambda, i)
  vi = grab(s, :v, i)

  psi = s.psi[lami]
  omega = s.omega[lami]
  mu = s.mu[lami]

  vnew = @. 1 / (vi + (psi ^ 2) * vi / omega)
  mnew = @. vnew * vi * psi * (yi - mu) / omega

  return @.rand(truncated(Normal(mnew, sqrt(vnew)), 0, Inf))
end

update_zetaC(m::Gtilde, s::T) where T = update_zeta(m, s, :C)
update_zetaT(m::Gtilde, s::T) where T = update_zeta(m, s, :T)

function update_mu(m::Gtilde, s::T) where T
  m_mu, s_mu = params(m.mu)
  vkernels = zero.(s.mu)
  mkernels = zero.(s.mu)

  for i in (:C, :T)
    y = grab(m, :y, i)
    lambda = grab(s, :lambda, i)
    zeta = grab(s, :zeta, i)
    v = grab(s, :v, i)
    for n in eachindex(lambda)
      k = lambda[n]
      vkernels[k] += v[n]
      mkernels[k] += (y[n] - s.psi[k] * zeta[n]) * v[n]
    end
  end

  vnew = 1 ./ (s_mu^-2 .+ vkernels ./ s.omega)
  mnew = vnew .* (m_mu/(s_mu^2) .+ mkernels ./ s.omega)

  return randn(m.K) .* sqrt.(vnew) + mnew
end

function update_tau(m::Gtilde, s::T) where T
  a, b = params(m.tau)
  new_shape = a + m.K * m.a_omega
  new_rate = b + sum(1 ./ s.omega)
  return rand(Gamma(new_shape, 1 / new_rate))
end

function update_sigma(m::Gtilde, s::T) where T
  return @.scalefromaltskewt(sqrt(s.omega), s.psi)
end

function update_phi(m::Gtilde, s::T) where T
  return @.skewfromaltskewt(sqrt(s.omega), s.psi)
end

function update_sigma_phi(m::Gtilde, s::T) where T
  return (sigma=update_sigma(m, s), phi=update_phi(m, s))
end

function update_omega(m::Gtilde, s::T) where T
  a, b = m.a_omega, s.tau
  akernel = zero.(s.omega)
  bkernel = zero.(s.omega)

  for i in (:C, :T)
    y = grab(m, :y, i)
    lam = grab(s, :lambda, i)
    zeta = grab(s, :zeta, i)
    v = grab(s, :v, i)
    for n in eachindex(y)
      k = lam[n]
      akernel[k] += 1
      bkernel[k] += v[n] * (y[n] - s.mu[k] - s.psi[k] * zeta[n])^2
    end
  end

  anew = a .+ akernel / 2
  bnew = b .+ bkernel / 2

  return @.rand(InverseGamma(anew, bnew))
end

function logprob_nu(m::Gtilde, s::T, nu::AbstractVector{<:Real}) where T
  logprior = sum(logpdf.(m.nu, nu))
  loglike = let
    _ll = 0.0
    for i in (:C, :T)
      lam = grab(s, :lambda, i)
      v = grab(s, :v, i)
      nu_i = nu[lam]
      _ll += sum(gammalogpdf.(nu_i / 2, 2 ./ nu_i, v))
    end
    _ll
  end
  return logprior + loglike
end

function update_psi(m::Gtilde, s::T) where T
  m_psi, s_psi = params(m.psi)
  vkernel = zero.(s.psi)
  mkernel = zero.(s.psi)

  # Update kernels
  for i in (:C, :T)
    y = grab(m, :y, i)
    lam = grab(s, :lambda, i)
    zeta = grab(s, :zeta, i)
    v = grab(s, :v, i)
    for n in eachindex(v)
      k = lam[n]
      vkernel[k] += zeta[n]^2 * v[n] / s.omega[k]
      mkernel[k] += zeta[n] * (y[n] - s.mu[k]) * v[n] / s.omega[k]
    end
  end

  vnew = 1 ./ (s_psi^-2 .+ vkernel)
  mnew = vnew .* (m_psi/s_psi^2 .+ mkernel)

  return randn(m.K) .* sqrt.(vnew) + mnew
end


text:
using Pkg;
Pkg.activate(joinpath(@__DIR__, "..")); Pkg.instantiate()
Pkg.activate(); Pkg.instantiate()

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Documenter, ArviZ

makedocs(
    modules = [ArviZ],
    sitename = "ArviZ.jl",
    pages = [
        "Home" => "index.md",
        "Quickstart" => "quickstart.md",
    ],
    checkdocs = :exports,
    format = Documenter.HTML(
        prettyurls = haskey(ENV, "CI"),
        sidebar_sitename = false,
        analytics = "UA-152384573-1",
    ),
)

deploydocs(repo = "github.com/arviz-devs/ArviZ.jl.git",)


text:
"""
Author: Kara I, 2021

The CLI entry point for Oboe. This file is responsible for parsing and validating
the CLI arguments and passing them to oboe-main.jl where the actual processing happens.

Run me as follows:
julia oboe-cli.jl --fips FIPS --agg AGG --name NAME [--force]

oboe-cli.jl
2021-07-10  v.0.0 Initial version
2021-07-14  v.0.1 Added support for --force
2021-07-15  v.0.1.1 YS: clarify help and validation error text
2021-08-12  v.0.2: updated to use FromFile-style module import
"""

using FromFile
using ArgParse
@from "oboe-main.jl" import OboeMain

#initialize CLI argument parsing
s = ArgParseSettings()

@add_arg_table! s begin
    "--fips"
        nargs = '+'
        help = "FIPS codes for the states that need to be processed, or ALL for all of contiguous US"
        arg_type = String
        required = true
    "--agg"
        nargs = 1
        help = "Pick aggregation type {tra|cty|ap|ste}"
        arg_type = String
        required = true
    "--name"
        nargs = 1
        help = "Output file names prefix, UPPERCASE. No hyphens - or underscores _"
        arg_type = String
        required = true
    "--force"
        action = :store_true
        help = "Overwrite existing output files"
end

function validate_args(agg, name)
    if agg ∉ ["tra", "cty", "ap", "ste"]
        throw(DomainError(agg, "Unknown aggregation type. Choose {tra|cty|ap|ste}"))
    end
    if any(char -> char ∈ ['-','_'], name)
        throw(DomainError(name, "Found a hyphen - or underscore _ in output prefix. These are not allowed"))
    end
end

useall = false # true if all contiguous US states should be used
useNW = false # true if the a~NW-tracts.dat dataset should be used instead of usa-tracts.dat

#Parse and process args
parsed_args = parse_args(ARGS, s)
#Get the target FIPS codes from the args, or use all states is ALL is given instead
fips = parsed_args["fips"]
if fips == ["ALL"]
    useall = true
elseif fips == ["NW"]
    useNW = true
    fips = ["41", "53"]
end

#Get the aggregation mode and the filename
agg = parsed_args["agg"][1]
name = uppercase(parsed_args["name"][1])
force = parsed_args["force"]

validate_args(agg, name)

#perform the data processing routine using the given args

if useall
    OboeMain.processOboe(name, agg, force=force)
else
    OboeMain.processOboe(name, agg, fips=fips, useNW=useNW, force=force)
end


text:
cd(@__DIR__)
using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()
using Documenter
using NNFEM
makedocs(sitename="NNFEM", modules=[NNFEM],
pages = Any[
    "index.md",
    "api.md"
],
authors = "Kailai Xu and Daniel (Zhengyu) Huang")

deploydocs(
    repo = "github.com/kailaix/NNFEM.jl.git",
)


text:
# This file was generated, do not modify it. # hide
train, test = partition(eachindex(y), 0.5, shuffle=true, rng=333)
fit!(mach, rows=train)
ypred = predict_mode(mach, rows=test)
misclassification_rate(ypred, y[test])


text:
module BidirectionalMaps

using Base: ImmutableDict

export ImmutableBimap

struct ImmutableBimap{S,T}
  left::ImmutableDict{S,T}
  right::ImmutableDict{T,S}

  function ImmutableBimap{S,T}(U, V) where {S,T}
      vec2pairs(X, Y) = [first(t) => last(t) for t in zip(X, Y)]
      left = ImmutableDict(vec2pairs(U, V)...)
      right = ImmutableDict(vec2pairs(V, U)...)
      new(left, right)
  end
end

end # module


text:
using JLD2, FileIO
using Test

include("bufferedreader.jl")
include("lookup3.jl")
include("internal.jl")
include("rw.jl")
include("append.jl")
include("groups.jl")
include("consistency.jl")
include("loadsave.jl")
include("modules.jl")
include("modules-nested.jl")
include("isreconstructed.jl")
include("backwards_compatibility.jl")
include("inlineunion.jl")
include("customserialization.jl")
include("compression.jl")


text:
import JuLIP

# TODO: get rid of this file altogether

using JuLIP.Potentials: @analytic,
                        AnalyticFunction,
                        F64fun

using Roots

import Base: Dict,
             ==

# -------------- Space Tranformations ---------------

# -------------- Cut-off Mechanisms ---------------

fcut_analyse(descr::String) = fcut_analyse(eval(Meta.parse(descr)))
fcut_analyse(args...) = fcut_analyse(args)
fcut_analyse(C::NBodyIPs.Cutoffs.NBCutoff) = C

function fcut_analyse(args::Tuple)
   sym = args[1]::Symbol
   args = args[2:end]
   if Symbol(sym) == :cos
      return CosCut(args...)
   elseif Symbol(sym) == :square
      return PolyCut(2, args[1])
   elseif Symbol(sym) == :cos2s
      return CosCut2s(args...)
   elseif Symbol(sym) == :penv
      return PolyCutSym(args...)
   elseif Symbol(sym) == :penv2s
      return PolyCut2s(args...)

   # elseif Symbol(sym) == :sw
   #    return let L=args[1], rcut=args[2]
   #       r -> cutsw(r, rcut, L), r -> cutsw_d(r, rcut, L), rcut
   #    end
   # elseif Symbol(sym) == :spline
   #    return let rc1=args[1], rc2=args[2]
   #       r -> cutsp(r, rc1, rc2), r -> cutsp_d(r, rc1, rc2), rc2
   #    end
   # elseif Symbol(sym) == :s2rat
   #    return let rnn=args[1], rin = args[2], rcut = args[3], p = args[4]
   #       f = @analytic r -> ( ((rnn/r)^p - (rnn/rcut)^p)^2 * ((rnn/r)^p - (rnn/rin)^p)^2 )
   #       r -> f.f(r) * (rin < r < rcut), r -> f.f_d(r) * (rin < r < rcut), rcut
   #    end

   else
      error("fcut_analyse: unknown symbol $(sym) for fcut.")
   end
end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ConnectionMonitorTestGroup <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    disable::Any # spec type: Union{ Nothing, Bool } # spec name: disable
    testConfigurations::Any # spec type: Union{ Nothing, Vector{String} } # spec name: testConfigurations
    sources::Any # spec type: Union{ Nothing, Vector{String} } # spec name: sources
    destinations::Any # spec type: Union{ Nothing, Vector{String} } # spec name: destinations

    function ConnectionMonitorTestGroup(;name=nothing, disable=nothing, testConfigurations=nothing, sources=nothing, destinations=nothing)
        o = new()
        validate_property(ConnectionMonitorTestGroup, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(ConnectionMonitorTestGroup, Symbol("disable"), disable)
        setfield!(o, Symbol("disable"), disable)
        validate_property(ConnectionMonitorTestGroup, Symbol("testConfigurations"), testConfigurations)
        setfield!(o, Symbol("testConfigurations"), testConfigurations)
        validate_property(ConnectionMonitorTestGroup, Symbol("sources"), sources)
        setfield!(o, Symbol("sources"), sources)
        validate_property(ConnectionMonitorTestGroup, Symbol("destinations"), destinations)
        setfield!(o, Symbol("destinations"), destinations)
        o
    end
end # type ConnectionMonitorTestGroup

const _property_map_ConnectionMonitorTestGroup = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("disable")=>Symbol("disable"), Symbol("testConfigurations")=>Symbol("testConfigurations"), Symbol("sources")=>Symbol("sources"), Symbol("destinations")=>Symbol("destinations"))
const _property_types_ConnectionMonitorTestGroup = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("disable")=>"Bool", Symbol("testConfigurations")=>"Vector{String}", Symbol("sources")=>"Vector{String}", Symbol("destinations")=>"Vector{String}")
Base.propertynames(::Type{ ConnectionMonitorTestGroup }) = collect(keys(_property_map_ConnectionMonitorTestGroup))
Swagger.property_type(::Type{ ConnectionMonitorTestGroup }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitorTestGroup[name]))}
Swagger.field_name(::Type{ ConnectionMonitorTestGroup }, property_name::Symbol) =  _property_map_ConnectionMonitorTestGroup[property_name]

function check_required(o::ConnectionMonitorTestGroup)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("testConfigurations")) === nothing) && (return false)
    (getproperty(o, Symbol("sources")) === nothing) && (return false)
    (getproperty(o, Symbol("destinations")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ConnectionMonitorTestGroup }, name::Symbol, val)
end


text:
import TimeZones: build

# ENV variable allows users to modify the default to be "latest". Do NOT use "latest"
# as the default here as can make it difficult to debug to past versions of working code.
# Note: Also allows us to only download a single tzdata version during CI jobs.
build(get(ENV, "JULIA_TZ_VERSION", "2019c"))