text:
using DailyCheckIn, Test, Documenter

@testset "Documentation" begin
    using Documenter, DailyCheckIn

    DocMeta.setdocmeta!(DailyCheckIn,
                       :DocTestSetup,
                       :(using DailyCheckIn;
                         givenname = "José Bayoán";
                         surname = "Santiago Calderón";
                         division = "SDAD";
                         health = "good";
                         email = "jbs3hp@virginia.edu";
                         wd = RemoteWebDriver(
                             Capabilities("chrome"),
                             host = get(ENV, "WEBDRIVER_HOST", "localhost"),
                             port = parse(Int, get(ENV, "WEBDRIVER_PORT", "4444")),
                             path = ""
                             );
                         session = DailyCheckIn.WebDriver.Session(wd);
                         ENV["COLUMNS"] = 120; ENV["LINES"] = 30;),
                       recursive = true)
    # doctest(DailyCheckIn, fix = true)
    makedocs(sitename = "DailyCheckIn",
             modules = [DailyCheckIn],
             pages = [
                 "Home" => "index.md",
                 "API" => "api.md",
                ],
             source = joinpath("..", "docs", "src"),
             build = joinpath("..", "docs", "build"),
            #  source = joinpath("docs", "src"),
            #  build = joinpath("docs", "build"),
             )
    @test true
end

max_stars_repo_name: uva-bi-sdad/COVID-Daily-Check-In

max_issues_repo_name: uva-bi-sdad/COVID-Daily-Check-In

max_forks_repo_name: uva-bi-sdad/COVID-Daily-Check-In


text:
int_rules_6_1_13 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 6.1.13*(d+e*x)^m*sinh(a+b*x+c*x^2)^n =#
    @apply_utils Antiderivative(sinh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~a + ~b * ~x + ~c * (~x) ^ 2), ~x) - (1 / 2) * Antiderivative(ℯ ^ ((-(~a) - ~b * ~x) - ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(cosh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (1 / 2) * Antiderivative(ℯ ^ (~a + ~b * ~x + ~c * (~x) ^ 2), ~x) + (1 / 2) * Antiderivative(ℯ ^ ((-(~a) - ~b * ~x) - ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x)
    @apply_utils Antiderivative(sinh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrigReduce(sinh(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative(cosh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrigReduce(cosh(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative(sinh(~v) ^ ~(n'), ~x) => Antiderivative(sinh(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- IGtQ(~n, 0) && (QuadraticQ(~v, ~x) && Not(QuadraticMatchQ(~v, ~x)))
    @apply_utils Antiderivative(cosh(~v) ^ ~(n'), ~x) => Antiderivative(cosh(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- IGtQ(~n, 0) && (QuadraticQ(~v, ~x) && Not(QuadraticMatchQ(~v, ~x)))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) * sinh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * cosh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && EqQ(~b * ~e - 2 * ~c * ~d, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) * cosh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * sinh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && EqQ(~b * ~e - 2 * ~c * ~d, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) * sinh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * cosh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) - ((~b * ~e - 2 * ~c * ~d) / (2 * ~c)) * Antiderivative(sinh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ(~b * ~e - 2 * ~c * ~d, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) * cosh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * sinh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) - ((~b * ~e - 2 * ~c * ~d) / (2 * ~c)) * Antiderivative(cosh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && NeQ(~b * ~e - 2 * ~c * ~d, 0)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * sinh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * (~d + ~e * ~x) ^ (~m - 1) * cosh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) - (((~e) ^ 2 * (~m - 1)) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 2) * cosh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ(~m, 1) && EqQ(~b * ~e - 2 * ~c * ~d, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * cosh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (~e * (~d + ~e * ~x) ^ (~m - 1) * sinh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) - (((~e) ^ 2 * (~m - 1)) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 2) * sinh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ(~m, 1) && EqQ(~b * ~e - 2 * ~c * ~d, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * sinh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => ((~e * (~d + ~e * ~x) ^ (~m - 1) * cosh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) - ((~b * ~e - 2 * ~c * ~d) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 1) * sinh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x)) - (((~e) ^ 2 * (~m - 1)) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 2) * cosh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ(~m, 1) && NeQ(~b * ~e - 2 * ~c * ~d, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * cosh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => ((~e * (~d + ~e * ~x) ^ (~m - 1) * sinh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (2 * ~c) - ((~b * ~e - 2 * ~c * ~d) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 1) * cosh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x)) - (((~e) ^ 2 * (~m - 1)) / (2 * ~c)) * Antiderivative((~d + ~e * ~x) ^ (~m - 2) * sinh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (GtQ(~m, 1) && NeQ(~b * ~e - 2 * ~c * ~d, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * sinh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * sinh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (~e * (~m + 1)) - ((2 * ~c) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 2) * cosh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (LtQ(~m, -1) && EqQ(~b * ~e - 2 * ~c * ~d, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * cosh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => ((~d + ~e * ~x) ^ (~m + 1) * cosh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (~e * (~m + 1)) - ((2 * ~c) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 2) * sinh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (LtQ(~m, -1) && EqQ(~b * ~e - 2 * ~c * ~d, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * sinh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (((~d + ~e * ~x) ^ (~m + 1) * sinh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (~e * (~m + 1)) - ((~b * ~e - 2 * ~c * ~d) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 1) * cosh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x)) - ((2 * ~c) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 2) * cosh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (LtQ(~m, -1) && NeQ(~b * ~e - 2 * ~c * ~d, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~m * cosh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => (((~d + ~e * ~x) ^ (~m + 1) * cosh(~a + ~b * ~x + ~c * (~x) ^ 2)) / (~e * (~m + 1)) - ((~b * ~e - 2 * ~c * ~d) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 1) * sinh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x)) - ((2 * ~c) / ((~e) ^ 2 * (~m + 1))) * Antiderivative((~d + ~e * ~x) ^ (~m + 2) * sinh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (LtQ(~m, -1) && NeQ(~b * ~e - 2 * ~c * ~d, 0))
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * sinh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => Unintegrable((~d + ~e * ~x) ^ ~m * sinh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * cosh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2), ~x) => Unintegrable((~d + ~e * ~x) ^ ~m * cosh(~a + ~b * ~x + ~c * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * sinh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrigReduce((~d + ~e * ~x) ^ ~m, sinh(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative((~(d') + ~(e') * ~x) ^ ~(m') * cosh(~(a') + ~(b') * ~x + ~(c') * (~x) ^ 2) ^ ~n, ~x) => Antiderivative(ExpandTrigReduce((~d + ~e * ~x) ^ ~m, cosh(~a + ~b * ~x + ~c * (~x) ^ 2) ^ ~n, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~m], ~x) && IGtQ(~n, 1)
    @apply_utils Antiderivative((~u) ^ ~(m') * sinh(~v) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * sinh(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- FreeQ(~m, ~x) && (IGtQ(~n, 0) && (LinearQ(~u, ~x) && (QuadraticQ(~v, ~x) && Not(LinearMatchQ(~u, ~x) && QuadraticMatchQ(~v, ~x)))))
    @apply_utils Antiderivative((~u) ^ ~(m') * cosh(~v) ^ ~(n'), ~x) => Antiderivative(ExpandToSum(~u, ~x) ^ ~m * cosh(ExpandToSum(~v, ~x)) ^ ~n, ~x) <-- FreeQ(~m, ~x) && (IGtQ(~n, 0) && (LinearQ(~u, ~x) && (QuadraticQ(~v, ~x) && Not(LinearMatchQ(~u, ~x) && QuadraticMatchQ(~v, ~x)))))
end

ext: jl

max_stars_repo_path: src/rules/6 Hyperbolic functions/6.1 Hyperbolic sine/6.1.13 (d+e x)^m sinh(a+b x+c x^2)^n.jl


text:
# Modified from my (tbreloff) implementation of OnlinePCA in the OnlineStats package.
# That implementation had this note:
#       Roughly based on Weng et al (2003): "Candid covariance-free incremental principal component analysis"
#       used https://github.com/kevinhughes27/pyIPCA/blob/master/pyIPCA/ccipca.py as a reference

abstract PreprocessStep <: Learnable
immutable NoPreprocessing <: PreprocessStep end
transform!(prep::NoPreprocessing, x) = x

#=
Solving for a dimension-reduced Y = VX, where X (nin x nobs) is the original data,
and Y (nout x nobs) is projected. V (nout x nin) is a matrix where the columns are the first
nout eigenvectors of Σ = XX'/n (the covariance of X).
e (nout x 1) is a vector of eigenvectors of the covariance matrix.

We compute e and V by incrementally updating e and U, where U is a (nout x nin) matrix with the properties:
		Uᵢ = iᵗʰ column of U
		Vᵢ = iᵗʰ column of V
			 = iᵗʰ eigenvector of X-covariance
			 = Uᵢ / ‖Uᵢ‖
		eᵢ = iᵗʰ eigenvalue of X-covariance
		   = ‖Uᵢ‖
note: V is the normalized version of U
=#

type Whiten{T,METHOD,W<:Weight} <: PreprocessStep
    nin::Int        # number of input vars
    nout::Int       # number of principal components
    wgt::W
    U::Matrix{T}    # (nout x nin)
    V::Matrix{T}    # (nout x nin)
    e::Vector{T}    # (nout x 1)
    x̄::Vector{T}    # (nin x 1) mean(x)
    input::SumNode{T,1}
    output::OutputNode{T,1}
 end
 Whiten(nin::Int, nout::Int = nin; kw...) = Whiten(Float64, nin, nout; kw...)

function Whiten{T}(::Type{T},
                    nin::Int,
                    nout::Int = nin;
                    method::Symbol = :zca,  # choose from pca, whitened_pca, zca
                    lookback::Int = 100,
                    α::Float64 = NaN,
                    wgt::Weight =  BoundedEqualWeight(isnan(α) ? lookback : α)
                   )
    @assert method in (:pca, :whitened_pca, :zca)
    if method == :zca && nin != nout
        error("ZCA whitening requires nin==nout.  Got: nin=$nin, nout=$nout")
    end
    Whiten{T,method,typeof(wgt)}(
        nin,
        nout,
        wgt,
        zeros(T,nout,nin),
        zeros(T,nout,nin),
        zeros(T,nout),
        zeros(T,nin),
        InputNode(T, nin),
        OutputNode(T, nout)
    )
end

function Base.empty!{T}(o::Whiten{T})
    for a in (o.U,o.V,o.e,o.x̄)
        fill!(a, zero(T))
    end
    o.wgt.nobs = 0
    o
end

function center!{T<:AbstractArray}(x::T, x̄::T)
    for i=1:length(x)
        x[i] -= x̄[i]
    end
end

function learn!{T}(o::Whiten{T}, x::AbstractVector = input_value(o))
    # get the smoothing param
    OnlineStats.updatecounter!(o.wgt)
	λ = weight(o.wgt)

    # update the mean
    smooth!(o.x̄, x, λ)

    # copy and center x
    x = x - o.x̄
    # center!(x, o.x̄)

    #=
    TODO:
        - remove temporaries
        - handle ZCA whitening: for V'V = Σ⁻¹, D = diag(e)
            y_pca = D⁻¹/² * V'
            y_zca = V * y_pca
    =#

    n = min(o.nout, o.wgt.nobs)
	@inbounds for i in 1:n
        # TODO this could be more robust?
		if o.e[i] == zero(T)
			# initialize ith principal component
			Uᵢ = x
		else
			# update the ith eigvec/eigval
			Uᵢ = view(o.U, i, :)
			smooth!(Uᵢ, x * (dot(x, Uᵢ) / o.e[i]), λ)
        end
		eᵢ = norm(Uᵢ)
        # if eᵢ ≈ zero(T)
        #     warn("Got small e[$i] ($(eᵢ)) while computing iteration $(o.wgt.nobs) for:\n  $o")
        # end
		Vᵢ = Uᵢ / max(eᵢ, T(1e-8))

		# store these updates
		o.U[i,:] = Uᵢ
		o.V[i,:] = Vᵢ
		o.e[i] = eᵢ

        if i < n
    		# subtract projection on ith PC for the next step
    		x -= dot(x, Vᵢ) * Vᵢ
        end
	end
end

grad!(o::Whiten) = error("Can't backprop through Whiten")

# -------------------------------------------------------------------
# transform defs

# TODO: cut down on temporaries

function normal_pca{T}(o::Whiten{T})
    x = input_value(o) - o.x̄
    y = output_value(o)
    y[:] = o.V * x
    y
end

function whitened_pca{T}(o::Whiten{T})
    y = normal_pca(o)
    y[:] = y ./ sqrt.(o.e .+ T(1e-8))
    y
end

function zca{T}(o::Whiten{T})
    y = whitened_pca(o)
    y[:] = o.V' * y
    y
end

# "normal" pca... should match MultivariateStats.fit(PCA, x, maxdimout=k)
#   y = V (x - mean(x))
transform!{T}(o::Whiten{T,:pca}) = normal_pca(o)

# whitened pca... multiply by D⁻¹/², where D = diag(e) is a diagonal matrix
# with the eigenvalues e on the diagonal
#   yᵢ = D⁻¹/² V (x - mean(x))
transform!{T}(o::Whiten{T,:whitened_pca}) = whitened_pca(o)

# whitened zca... multiply by (V' * D⁻¹/²), where D = diag(e) is a diagonal matrix
# with the eigenvalues e on the diagonal.
# Note: ZCA is similar to whitened PCA, except that it rotates the whitened projection
#       back into "input space".  If you want whitening while maintaining the directions
#       of the inputs, this is the one to use.
#   yᵢ = D⁻¹/² V (x - mean(x))
transform!{T}(o::Whiten{T,:zca}) = zca(o)


text:
struct L0
    a::Int
end

struct L1 <: Layout
    a::Int
end

struct L2 <: Layout
    a::LForwardReference{L2, :L3}
end

struct L3 <: Layout
    a::L2
end

struct L4 <: Layout
    a::LForwardReference{L4, :L_not_defined}
end

struct L5 <: Layout
    a::LVarVector{Int,1}
end

struct L6 <: Layout
    a::LVarVector{Int,1}
    b::Int
end
struct L7 <: Layout
    a::Ptr{Ptr{Int}}
end

@testset "simple layout templates" begin
    @test_throws ArgumentError is_template_fixed(L0)
    @test_throws ArgumentError is_template_fixed(Vector{Int})
    @test is_template_fixed(Float64)
    @test is_template_fixed(L1)
    @test is_template_fixed(Ptr{L1})
    @test is_template_fixed(Ptr{L5})
    @test is_template_fixed(LFixedVector{Int,10})
    @test_throws ArgumentError is_template_fixed(LFixedVector{Int})
    @test is_template_variable(LVarVector{Int,1})
    @test is_template_fixed(L2)
    @test_throws UndefVarError is_template_fixed(L4)
    @test !is_template_fixed(L5)
    @test_throws ArgumentError is_template_fixed(L6)
    @test_throws ArgumentError is_template_fixed(L7)
    @test is_template_fixed(LFixedVector{L1,1})
    @test is_template_fixed(LFixedVector{L5,0})
    @test !is_template_fixed(LFixedVector{L5,1})
    @test_throws ArgumentError !is_template_fixed(LFixedVector{L5,2})
    @test !is_template_fixed(LVarVector{L1,1})
    @test !is_template_fixed(LVarVector{L5,0})
end


text:
# Jasper den Ouden 02-08-2012
# Robert J Ennis   03-13-2013
# Placed in the public domain

module GetC

export @getCFun

const __ofTypeSym = symbol("::")

macro getCFun(lib, jlFun, cFun)
    cName = cFun.args[1].args[1]

    arguments = map(function (arg)
                        if isa(arg, Symbol)
                            arg = Expr(__ofTypeSym, arg)
                        end
                        return arg
                    end, cFun.args[1].args[2:end])

    # Get info out of arguments of `cFun`
    argumentNames = map(arg->arg.args[1], arguments)
    returnType    = cFun.args[2]
    inputTypes    = map(arg->arg.args[2], arguments)

    # Construct the result.
    cSym       = Expr(:quote,cName)
    symAndLib = :($cSym, $lib)

    body = Expr(:ccall, symAndLib, returnType, Expr(:tuple,inputTypes...), argumentNames...)
    ret  = Expr(:function, Expr(:call, jlFun, argumentNames...), body)
    return esc(ret)

end

end #module GetC


text:
module EmptyDomainsTests

using Gridap
import Gridap: ∇
using Test

# Analytical functions
u(x) = VectorValue( x[1]^2 + 2*x[2]^2, -x[1]^2 )
∇u(x) = TensorValue( 2*x[1], 4*x[2], -2*x[1], zero(x[1]) )
Δu(x) = VectorValue( 6, -2 )
p(x) = x[1] + 3*x[2]
∇p(x) = VectorValue(1,3)
s(x) = -Δu(x)
f(x) = -Δu(x) + ∇p(x)
g(x) = tr(∇u(x))
∇(::typeof(u)) = ∇u
∇(::typeof(p)) = ∇p

# Mesh
cells = (10,10)
domain = (0,1,0,1)
model = CartesianDiscreteModel(domain,cells)
labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet",[1,2,5])
add_tag_from_tags!(labels,"neumann",[6,7,8])

# Domains and measures
Ω = Interior(model,Int[])
Ωs = Interior(model,Int[])
Ωf = Interior(model,Int[])
Λ = Boundary(model,Int[])
Γ = Interface(Ωf,Ωs)
n_Λ = get_normal_vector(Λ)
n_Γ = get_normal_vector(Γ)
k = 2
degree = 2*k
dΩ = Measure(Ω,degree)
dΩs = Measure(Ωs,degree)
dΩf = Measure(Ωf,degree)
dΛ = Measure(Λ,degree)
dΓ = Measure(Γ,degree)

# FE Spaces
reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},k)
reffe_p = ReferenceFE(lagrangian,Float64,k-1,space=:P)
V = TestFESpace(Ω,reffe_u,dirichlet_tags="dirichlet")
Q = TestFESpace(Ωf,reffe_p)
U = TrialFESpace(V,u)
P = Q
Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

# Weak form
a((u,p),(v,q)) =
  ∫( ∇(v)⊙∇(u) )*dΩs +
  ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )*dΩf

l((v,q)) =
  ∫( v⋅s )*dΩs +
  ∫( v⋅f + q*g )*dΩf +
  ∫( v⋅(n_Λ⋅∇u) - (n_Λ⋅v)*p )*dΛ +
  ∫( - (n_Γ.⁺⋅v.⁺)*p )*dΓ

# FE problem
op = AffineFEOperator(a,l,X,Y)
A = get_matrix(op)
b = get_vector(op)
@test size(A) == (0,0)
@test size(b) == (0,)
uh, ph = FEFunction(X,Float64[])

# Errors
eu = u - uh
ep = p - ph
eu_l2 = sqrt(sum(∫( eu⋅eu )*dΩ))
eu_h1 = sqrt(sum(∫( eu⋅eu + ∇(eu)⊙∇(eu) )*dΩ))
ep_l2 = sqrt(sum(∫( ep*ep )*dΩf))
tol = 1.0e-9
@test eu_l2 < tol
@test eu_h1 < tol
@test ep_l2 < tol

end # module


text:
include("common.jl")

function IDT_algorithm{T <: Real}(time::Array{T, 1}, coords::Array{T, 2}, durthresh, disthresh)
    i = 1
    n = length(time) - 1
    fixcluster = zeros(n)
    D(x, y) = (maximum(x) - minimum(x)) + (maximum(y) - minimum(y))

    while i < n
        window = getwindow(time, coords, durthresh, i)
        dispersion = D(window[:, 1], window[:, 2])
        winsize = size(window, 1)

        if (winsize + i) < n && dispersion <= disthresh
            start = i + winsize

            while dispersion <= disthresh
                start += 1
                window = vcat(window, coords[start, :])
                dispersion = D(window[:, 1], window[:, 2])
            end

            fixcluster[i] = 1
            i += size(window, 1)
        else
            i += 1
        end
    end

    hcat(time, coords, assign(convert(BitArray{1}, fixcluster)))
end

function getwindow{T <: Real}(time::Array{T, 1}, coords::Array{T, 2}, durthresh, start::Int64)

    next = start + 1
    maxn = size(coords, 1)

    while next < maxn && ((time[next] - time[start]) < durthresh)
        next += 1
    end

    coords[start:next, :]
end

function IDT{T <: Real}(time::Array{T, 1}, x::Array{T, 1}, y::Array{T, 1}; durthresh = 0, disthresh = 0)
    IDT_algorithm(time, hcat(x, y), durthresh, disthresh)
end

function IDT{T <: Real}(time::Array{T, 1}, coords::Array{T, 2}; durthresh = 0, disthresh = 0)
    IDT_algorithm(time, coords, durthresh, disthresh)
end

function IDT{T <: Real}(dat::Array{T, 2}; durthresh = 0, disthresh = 0)
    @assert size(dat, 2) == 3
    IDT_algorithm(dat[:, 1], dat[:, 2:3], durthresh, disthresh)
end


text:
function shuffle!(a::Vector)::typeof(a)
    l = length(a);
    for i in 1:l
        ri = rand(1:l)
        t = a[ri];
        a[ri] = a[i]
        a[i] = t;
    end
    return a;
end


text:
export LinearApproximator, LinearVApproximator, LinearQApproximator

using LinearAlgebra: dot

struct LinearApproximator{N,O} <: AbstractApproximator
    weights::Array{Float64,N}
    optimizer::O
end

#####
# LinearVApproximator
#####

const LinearVApproximator = LinearApproximator{1}

LinearVApproximator(;n, init=0., opt=Descent(1.0)) = LinearApproximator(fill(init, n), opt)

(V::LinearVApproximator)(s) = dot(s, V.weights)

function RLBase.update!(V::LinearVApproximator, correction::Pair)
    w = V.weights
    s, Δ = correction
    w̄ = s .* Δ
    Flux.Optimise.update!(V.optimizer, w, w̄)
end

#####
# LinearQApproximator
#####

const LinearQApproximator = LinearApproximator{2}

LinearQApproximator(;n_state, n_action, init=0., opt=Descent(1.0)) = LinearApproximator(fill(init, n_state, n_action), opt)

(Q::LinearQApproximator)(s) = [dot(s, c) for c in eachcol(Q.weights)]
(Q::LinearQApproximator)(s, a) = dot(s, @view(Q.weights[:, a]))

function RLBase.update!(Q::LinearQApproximator, correction::Pair)
    (s, a), Δ = correction
    @views w = Q.weights[:, a]
    w̄ = s .* Δ
    Flux.Optimise.update!(Q.optimizer, w, w̄)
end


text:
import FastaIO.writefasta
export writefasta, readfastastrains

# Closure for assigning labels to strains
let strain_number = 0
	global new_strain_number() = strain_number+=1
	global reset_strain_number() = strain_number = 0
end
function new_strain_label()
	return "strain_$(new_strain_number())"
end

"""
	writefasta(s::AbstractStrain, fields; fillvals = false)
	writefasta(f::IO, s::AbstractStrain, fields; fillvals=false)
	writefasta(S::Array{<:AbstractStrain}, fields; fillvals = false)
	writefasta(f::String, S::Array{<:AbstractStrain}, fields; fillvals=false, mode="w")

Write strain `s` to a fasta format. Header is built using `fields` from `s.data`.
"""
function writefasta(s::AbstractStrain, fields; fillvals = false)
	# Make the header
	header = mapreduce(*, fields, init = "") do x
		if haskey(s.data, String(x))
			"$(get(s.data, String(x), 0))|"
		elseif haskey(s.data, Symbol(x))
			"$(get(s.data, Symbol(x), 0))|"
		elseif fillvals
			"?|"
		else
			@error "Field `$x` not in strain"
		end
		# "$(get(s.data, x) do ; fillvals ? "?|" : @error "Field $x not in strain"; end)|"  -- One liner, but does not handle the symbol/string
	end[1:end-1] # Removing trailing '|'

	return (header, String(s.seq))
end
writefasta(f::IO, s::AbstractStrain, fields; fillvals=false) = writefasta(f, [writefasta(s, fields, fillvals=fillvals)])
function writefasta(f::String, S::Array{<:AbstractStrain}, fields; fillvals=false, mode="w")
	open(f, mode) do io
		for s in S
			writefasta(io, s, fields, fillvals=fillvals)
		end
	end
end
function writefasta(S::Array{<:AbstractStrain}, fields; fillvals = false)
	out = Array{Tuple{String, String},1}(undef, length(S))
	for (i,s) in enumerate(S)
		out[i] = writefasta(s, fields, fillvals = fillvals)
	end
	return out
end

function readfastastrains(f::Union{AbstractString,IO}, ::Val{A}, headerfields; separator = '|', strainfilters=[x->true], ignore_read_errors=false) where A <: BioSequences.Alphabet
	strains = Array{Strain{A},1}(undef, 0)
	nfiltered = 0
	nunread = 0
	ntot = 0
	typeof(f) <: AbstractString ? println("Reading $f...") : println("Reading alignment...")
	for (n,s) in FastaReader(f)
		dat = parse_header(n, headerfields, separator)
		st = try
			Strain(LongSequence{A}(s), dat)
		catch err
			if ignore_read_errors
				st = Strain(:aa)
			elseif occursin("X",s)
				st = Strain(LongSequence{A}(replace(s, "X"=>"N")), dat)
			elseif occursin("x",s)
				st = Strain(LongSequence{A}(replace(s, "x"=>"N")), dat)
			else
				println(s)
				error(err)
			end
		end
		if BioTools.isempty(st)
			nunread += 1
		elseif mapreduce(f->f(st), *, strainfilters, init=true)
			push!(strains, st)
		else
			nfiltered += 1
		end
		ntot += 1
	end
	println("Read $(length(strains)) strains out of $ntot. Filtered $nfiltered. Could not read $nunread")
	return strains
end
function readfastastrains(f::Union{AbstractString,IO}, ::Val{A}, headerfields; separator = '|', strainfilters=[x->true], ignore_read_errors=false) where A <: Integer
	strains = Array{ArtificialStrain{A},1}(undef, 0)
	nfiltered = 0
	nunread = 0
	ntot = 0
	typeof(f) <: AbstractString ? println("Reading $f...") : println("Reading alignment...")
	for (n,s) in FastaReader(f)
		dat = parse_header(n, headerfields, separator)
		st = ArtificialStrain([parse(A,x) for x in s], dat)
		if BioTools.isempty(st)
			nunread += 1
		elseif mapreduce(f->f(st), *, strainfilters, init=true)
			push!(strains, st)
		else
			nfiltered += 1
		end
		ntot += 1
	end
	println("Read $(length(strains)) strains out of $ntot. Filtered $nfiltered. Could not read $nunread")
	return strains
end
"""
	readfastastrains(f::Union{AbstractString,IO}, sequence_type::Symbol, headerfields; separator = '|', strainfilters=[x->true])
	readfastastrains(f::Union{AbstractString,IO}, sequence_type::DataType, headerfields; separator = '|', strainfilters=[x->true])

Possible symbols are `$(BioTools.sequenceymbols)`. Type of output will depend on the symbol used.

Implemented filters are `hasdate` and `gapfilter`
"""
readfastastrains(f::Union{AbstractString,IO}, sequence_type::Symbol, headerfields; separator = '|', strainfilters=[x->true], ignore_read_errors=false) = readfastastrains(f, Val(BioTools.type(sequence_type)), headerfields, separator=separator, strainfilters=strainfilters, ignore_read_errors=ignore_read_errors)
readfastastrains(f::Union{AbstractString,IO}, sequence_type::DataType, headerfields; separator = '|', strainfilters=[x->true], ignore_read_errors=false) = readfastastrains(f, Val(sequence_type), headerfields, separator=separator, strainfilters=strainfilters, ignore_read_errors=ignore_read_errors)

"""
- `headerfields` is the list of fields that should be parsed. `?` are ignored. If the header `h` is longer than `headerfields`, the end of it is ignored.
- `specialfields` are fields which require a special treatment. An example is the date (by default with `special_fields`), which calls `parse_date`.
"""
function parse_header(h, headerfields, separator; specialfields=special_fields)
	sh = split(h, separator)
	dat = Dict()
	for (i,f) in enumerate(headerfields)
		if in(f, specialfields)
			dat[f] = parse_special_field[f](sh[i])
		elseif !in(f, ignored_header_fields)
			dat[f] = sh[i]
		end
	end
	return dat
end

function parse_date(s::AbstractString)
	date = missing
	try
		date = Date(s)
	catch
		if !isnothing(match(r"-[0-9][0-9]*-XX$",s))
			# Day is missing, which does not matter
			date = Date(s[1:end-3])
		elseif !isnothing(match(r"[1-2][0-9][0-9][0-9]",s))
			# Month is missing, year is not --> date is missing
			date = missing
			# y = parse(Int64, match(r"[1-2][0-9][0-9][0-9]",hp["date"]).match)
		else
			date = missing
		end
	end
	return date
end

parse_float(s::AbstractString) = parse(Float64, s)

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
"""
# script test_noise

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-08-20

Script to add noise to input image, run LAMINART and plot layers and activation.
"""

using DrWatson
@quickactivate "MEngProject"
using MEngProject,
    CUDA,
    DifferentialEquations,
    PyPlot,
    NNlib,
    ImageFiltering,
    Images,
    ImageIO,
    MEngProject.LaminartKernels,
    MEngProject.LaminartInitFunc,
    MEngProject.Utils,
    BenchmarkTools,
    Test

using OrdinaryDiffEq,
    ParameterizedFunctions, LSODA, Sundials, DiffEqDevTools, Noise

batch = 1001

files = ["stairs_100gs.png"]
# files = ["kan_sq_cont_l.png"]
# files = ["kan_sq_cont_l.png", "stairs_100gs.png"]

tspan = (0.0f0, 800f0)

batch_ = string(batch, "_", rand(1000:9999))
mkdir(plotsdir(string("noise", batch_)))
for file in files[1:end]
    for noise in [0.01f0, 0.05f0, 0.1f0, 0.2f0, 0.5f0, 1f0]
        img = convert(Array{Float32,2}, load(datadir("img", files[1])))
        img = mult_gauss(img, noise)
        fig, ax = plt.subplots()
        im = ax.imshow(img, cmap = matplotlib.cm.gray, vmax = 1, vmin = 0)
        plt.title("Input image, noise \$σ=$noise\$")
        plt.axis("off")
        fig.tight_layout()
        plt.savefig(plotsdir(
            string("noise", batch_),
            string(file, "_input_noise_", noise, ".png"),
        ))
        close("all")
        img = nothing

        try
            p = LaminartInitFunc.parameterInit_conv_gpu_noise(
                datadir("img", file),
                Parameters.parameters_f32,
                noise,
            )

            u0 = cu(reshape(
                zeros(Float32, p.dim_i, p.dim_j * (5 * p.K + 2)),
                p.dim_i,
                p.dim_j,
                5 * p.K + 2,
                1,
            ))

            arr1 = similar(@view u0[:, :, 1:2, :])
            arr2 = similar(@view u0[:, :, 1:1, :])

            f = LaminartFunc.LamFunction(
                arr1, #x
                similar(arr1), #m
                similar(arr1), #s
                arr2, #x_lgn,
                similar(arr1), #C,
                similar(arr1), #H_z,
                similar(arr1), # dy_temp,
                similar(arr1), # dm_temp,
                similar(arr1), # dz_temp,
                similar(arr1), # ds_temp,
                similar(arr2), # dv_temp,
                similar(arr1), # H_z_temp,
                similar(arr2), #  V_temp_1,
                similar(arr2), #  V_temp_2,
                similar(arr1), #  Q_temp,
                similar(arr1), #   P_temp
            )

            prob = ODEProblem(f, u0, tspan, p)
            # 	@benchmark sol = solve(prob)
            sol = solve(prob)

            for t ∈ [25, 50, 100, 200, 400, 800]
                # 				for t ∈ [25,50,100]
                @inbounds begin
                    v0 = @view sol(t)[:, :, :, 1]
                    axMax = findmax(v0)[1]

                    for k ∈ 1:2:10
                        fig, ax = plt.subplots()

                        v1 = @view sol(t)[:, :, k, 1]
                        v2 = @view sol(t)[:, :, k+1, 1]
                        im = ax.imshow(
                            v1,
                            cmap = matplotlib.cm.PRGn,
                            vmax = axMax,
                            vmin = -axMax,
                        )
                        im2 = ax.imshow(
                            v2,
                            cmap = matplotlib.cm.RdBu_r,
                            vmax = axMax,
                            vmin = -axMax,
                            alpha = 0.5,
                        )

                        cbar = fig.colorbar(im2, shrink = 0.9, ax = ax)
                        cbar.ax.set_xlabel("\$k=2\$")
                        cbar = fig.colorbar(im, shrink = 0.9, ax = ax)
                        cbar.set_alpha(0.5)
                        cbar.draw_all()
                        cbar.ax.set_xlabel("\$k=1\$")
                        layer = Utils.layers[k]
                        plt.title("Layer: $layer, \$t=$t\$, noise \$σ=$noise\$")
                        plt.axis("off")
                        fig.tight_layout()
                        plt.savefig(plotsdir(
                            string("noise", batch_),
                            string(
                                file,
                                "_noise_",
                                noise,
                                "_t",
                                t,
                                "_",
                                Utils.la[k],
                                ".png",
                            ),
                        ))
                        close("all")
                    end

                    k = 11
                    fig, ax = plt.subplots()
                    v1 = @view sol[:, :, k, 1, t]
                    v2 = @view sol[:, :, k+1, 1, t]
                    im = ax.imshow(
                        v1,
                        cmap = matplotlib.cm.PRGn,
                        vmax = axMax,
                        vmin = -axMax,
                    )
                    im2 = ax.imshow(
                        v2,
                        cmap = matplotlib.cm.RdBu_r,
                        vmax = axMax,
                        vmin = -axMax,
                        alpha = 0.5,
                    )

                    cbar = fig.colorbar(im2, shrink = 0.9, ax = ax)
                    cbar.ax.set_xlabel("\$v^-\$")
                    cbar = fig.colorbar(im, shrink = 0.9, ax = ax)
                    cbar.set_alpha(0.5)
                    cbar.draw_all()
                    cbar.ax.set_xlabel("\$v^+\$")

                    layer = Utils.layers[k]
                    plt.title("Layer: $layer, \$t=$t\$, noise \$σ=$noise\$")
                    plt.axis("off")
                    fig.tight_layout()

                    plt.savefig(plotsdir(
                        string("noise", batch_),
                        string(
                            file,
                            "_noise_",
                            noise,
                            "_t",
                            t,
                            "_",
                            Utils.la[k],
                            ".png",
                        ),
                    ))
                    close("all")
                end
            end
        finally
            nothing
        end

        # time plot
        fig, axs = plt.subplots()

        for k ∈ 1:12
            v3 = @view sol[:, :, k, 1, end]
            v4 = @view sol[findmax(v3)[2][1], findmax(v3)[2][2], k, 1, :]
            layer = Utils.layers_1[k]
            axs.plot(sol.t, v4, Utils.lines[k], label = "$layer", alpha = 0.8)
        end
        axs.set_xlabel("Time")
        axs.set_ylabel("Activation")
        plt.title("Noise \$σ=$noise\$")
        plt.legend()
        fig.tight_layout()
        plt.savefig(plotsdir(
            string("noise", batch_),
            string(file, "_noise_", noise, "_time.png"),
        ))
        u0 = nothing
        p = nothing
        arr1 = nothing
        arr2 = nothing
        f = nothing
        prob = nothing
        sol = nothing
        close("all")
        CUDA.reclaim()
    end
end


text:
module TestGraft
export @syms
using Debug, Debug.AST, Debug.Eval
using Debug.Graft.instrument, Debug.is_trap, Debug.Meta.quot
import Debug.AST.is_emittable
export @test_graft

macro assert_fails(ex)
    quote
        if (try; $(esc(ex)); true; catch e; false; end)
            error($("@assert_fails: $ex didn't fail!"))
        end
    end
end

type Graft; ex; end
is_emittable(::Node{Graft}) = false

trap(::Any, ::Scope) = nothing
function trap(g::Node{Plain}, scope::Scope)
    g = exof(g)
    if isa(g, QuoteNode)
        g = g.value
        if isa(g, Node{Graft})
            debug_eval(scope, valueof(g).ex)
        end
    end
end

macro instrument_assume_global(trap_ex, ex)
    @gensym trap_var
    esc(quote
        const $trap_var = $trap_ex
        $(instrument(is_trap, trap_var, ex))
    end)
end

macro graft(ex)
    quot(Node(Graft(ex)))
end
macro test_graft(ex)
    :(@instrument_assume_global trap $(esc(ex)))
end

type T
    x
    y
end

@test_graft let
    # test read and write, both ways
    local x = 11
    @graft (@assert x == 11; x = 2)
    @assert x == 2

    # test tuple assignment in graft
    local y, z
    @graft y, z = 12, 23
    @assert (y, z) == (12, 23)

    # test updating operator in graft
    q = 3
    @graft q += 5
    @assert q == 8

    # test assignment to ref in graft (don't rewrite)
    a = [1:5;]
    @graft a[2]  = 45
    @graft a[3] += 1
    @assert isequal(a, [1,45,4,4,5])

    # test assignment to field in graft (don't rewrite)
    t = T(7,83)
    @graft t.x  = 71
    @graft t.y -= 2
    @assert (t.x == 71) && (t.y == 81)

    # test that k inside the grafted let block is separated from k outside
    k = 21
    @graft let k=5
        @assert k == 5
    end
    @assert k == 21

    # test that assignments to local vars don't leak out
    @graft let
        l = 3
    end

    # test that assignements to shared vars do leak out
    s = 1
    @graft let
        s = 6
    end
    @assert s == 6
end

@test_graft begin
    type TT
        x::Int
        function TT(x)
            local z
            @graft z=2x
            new(z)
        end
    end

    t = TT(42)
    @assert t.x == 84
end

@test_graft begin
    type TTT
        x::Int
        function TTT(x=11; m=1)
            local z
            @graft z=m*x
            new(z)
        end
    end

    t = TTT(m=3)
    @assert t.x == 33
end

@assert_fails @test_graft let
    @graft a = 3
end
@assert_fails @test_graft let
    @graft local d1 = 3
end
@assert_fails @test_graft let
    @graft local d2
end

end # module


text:
function quarter_date_to_number(date::Date)
    y = Dates.year(date)
    q = Dates.quarterofyear(date)
    return y + (q-1)*0.25
end

function quarter_number_to_date(datenum::Real)
    if datenum % 0.25 != 0
        throw(DomainError())
    end

    y = convert(Int, floor(datenum))
    q = convert(Int, (datenum % 1) / 0.25) + 1
    return quartertodate("$y-Q$q")
end

function plot_extension()
    be = typeof(Plots.backend())
    if be == Plots.GRBackend
        :pdf
    elseif be in [Plots.PlotlyBackend, Plots.PlotlyJSBackend]
        :html
    else
        :pdf
    end
end

function describe_series(m::AbstractModel, var::Symbol, class::Symbol;
                         detexify::Bool = false)
    res = if class in [:obs, :pseudo]
        dict = if class == :obs
            m.observable_mappings
        elseif class == :pseudo
            m.pseudo_observable_mappings
        end
        dict[var].name
    elseif class == :states
        string(var)
    elseif class in [:shocks, :stdshocks]
        replace(string(var), r"_sh$", "")
    else
        error("Invalid class: " * string(class))
    end

    detexify ? DSGE.detexify(res) : res
end

function series_ylabel(m::AbstractModel, var::Symbol, class::Symbol;
                       untrans::Bool = false,
                       fourquarter::Bool = false)
    if untrans && fourquarter
        error("Only one of untrans or fourquarter can be true")
    end

    if class in [:obs, :pseudo]
        dict = if class == :obs
            m.observable_mappings
        elseif class == :pseudo
            m.pseudo_observable_mappings
        end
        transform = dict[var].rev_transform

        if transform in [loggrowthtopct_annualized_percapita, loggrowthtopct_annualized, loggrowthtopct, loggrowthtopct_percapita]
            if untrans
                return "Q/Q Log Growth Rate"
            elseif fourquarter
                return "Percent 4Q Growth"
            else
                return "Percent Q/Q Annualized"
            end
        elseif transform in [logleveltopct_annualized_percapita, logleveltopct_annualized]
            if untrans
                return "Log Level"
            elseif fourquarter
                return "Percent 4Q Growth"
            else
                return "Percent Q/Q Annualized"
            end
        elseif transform == quartertoannual
            if untrans
                return "Percent Q/Q"
            else
                return "Percent Annualized"
            end
        elseif transform == identity
            ""
        else
            error("series_ylabel not implemented for transform: $transform")
        end
    elseif class == :stdshocks
        return "Standard Deviations"
    elseif class in [:states, :shocks]
        return ""
    else
        error("Invalid class: " * string(class))
    end
end

function save_plot(p::Plots.Plot, output_file::String = ""; verbose::Symbol = :low)
    if !isempty(output_file)
        output_dir = dirname(output_file)
        !isdir(output_dir) && mkpath(output_dir)
        Plots.savefig(output_file)

        if VERBOSITY[verbose] >= VERBOSITY[:low]
            println("Saved $output_file")
        end
    end
end


text:
# small functions supporting fciqmc!()
# versions without dependence on MPI.jl
using Base.Threads: nthreads

"""
    MultiScalar

Wrapper over a tuple that supports `+`, `-`, `min`, and `max`. Used with MPI communication
because `SVector`s are treated as arrays by `MPI.Allreduce` and `Tuples` do not support
scalar operations.

# Example

Suppose you want to compute the sum of a vector `dv` and also get the number of positive
elements it has in a single pass. You can use `MultiScalar`:

```julia
julia> dv = DVec(:a => 1, :b => -2, :c => 1);

julia> s, p = mapreduce(+, values(dv)) do v
    Rimu.MultiScalar(v, Int(sign(v) == 1))
end;

julia> s, p
(0, 2)
```

This will work with `MPIData`.

Note that only `MultiScalar`s with the same types can be operated on. This is a feature, as it
forces type stability.
"""
struct MultiScalar{T<:Tuple}
    tuple::T
end
MultiScalar(args...) = MultiScalar(args)
MultiScalar(v::SVector) = MultiScalar(Tuple(v))
MultiScalar(m::MultiScalar) = m
MultiScalar(arg) = MultiScalar((arg,))

const SVecOrTuple = Union{SVector,Tuple}

for op in (:+, :*, :max, :min)
    @eval function Base.$op(a::MultiScalar{T}, b::MultiScalar{T}) where {T}
        return MultiScalar($op.(a.tuple, b.tuple))
    end
    @eval function Base.$op(a::MultiScalar, b::SVecOrTuple)
        return $op(a, MultiScalar(b))
    end
    @eval function Base.$op(a::SVecOrTuple, b::MultiScalar)
        return $op(MultiScalar(a), b)
    end
end

Base.iterate(m::MultiScalar, args...) = iterate(m.tuple, args...)
Base.length(m::MultiScalar) = length(m.tuple)

# three-argument version
"""
    sort_into_targets!(target, source, stats) -> agg, wm, agg_stats
Aggregate coefficients from `source` to `agg` and from `stats` to `agg_stats`
according to thread- or MPI-level parallelism. `wm` passes back a reference to
working memory.
"""
sort_into_targets!(target, w, stats) =  w, target, stats
# default serial (single thread, no MPI) version: don't copy just swap

combine_stats(stats) = sum(stats)
combine_stats(stats::MultiScalar) = stats # special case for fciqmc_step!() using ThreadsX

function sort_into_targets!(target, ws::NTuple{NT,W}, statss) where {NT,W}
    # multi-threaded non-MPI version
    zero!(target)
    for w in ws # combine new walkers generated from different threads
        add!(target, w)
    end
    return target, ws, combine_stats(statss)
end
# three argument version for MPIData to be found in RMPI.jl


text:
module BenchHeap

using DataStructures
using BenchmarkTools
using Random

function push_heap(h::AbstractHeap, xs::Vector)
    n = length(xs)

    for i = 1 : n
        push!(h, xs[i])
    end
end

function pop_heap(h::AbstractHeap)
    n = length(h)

    for i = 1 : n
        pop!(h)
    end
end

suite = BenchmarkGroup()

heaptypes = [BinaryHeap, MutableBinaryHeap]
aexps = [1,3]
datatypes = [Int, Float64]
baseorderings = Dict(
    "Min" => Base.ForwardOrdering,
    #"Max" => Base.ReverseOrdering,
    )
fastfloatorderings = Dict(
    "Min" => DataStructures.FasterForward,
    "Max" => DataStructures.FasterReverse,
    )

for heap in heaptypes
    for aexp in aexps
        for dt in datatypes
            Random.seed!(0)
            a = rand(dt, 10^aexp)

            # Dict types to force use of abstract type if containing single value
            orderings = Dict{String, DataType}(baseorderings)
            if dt == Float64
                # swap to faster ordering operation
                for (k,v) in orderings
                    if haskey(fastfloatorderings, k)
                        orderings["Slow"*k] = v
                        orderings[k] = fastfloatorderings[k]
                    end
                end
            end

            for (ord_str, ord) in orderings
                prepath = [string(heap)]
                postpath = [string(dt), "10^"*string(aexp), ord_str]
                suite[vcat(prepath, ["make"], postpath)] =
                    @benchmarkable $(heap){$dt,$ord}($a)
                suite[vcat(prepath, ["push"], postpath)] =
                    @benchmarkable push_heap(h, $a) setup=(h=$(heap){$dt,$ord}())
                suite[vcat(prepath, ["pop"], postpath)] =
                    @benchmarkable pop_heap(h) setup=(h=$(heap){$dt,$ord}($a))
            end
        end
    end
end

fast_extreme_orderings = Dict(
    nsmallest => DataStructures.FasterForward(),
    nlargest => DataStructures.FasterReverse(),
    )

for func in [nlargest, nsmallest]
    fastord = fast_extreme_orderings[func]
    for aexp in [4]
        Random.seed!(0);
        a = rand(10^aexp);
        for nexp in [2]
            n = 10^nexp
            suite[["Slow " * string(func), "a=rand(10^"*string(aexp)*")", "n=10^"*string(nexp)]] =
                @benchmarkable $(func)($n, $a)
            suite[[string(func), "a=rand(10^"*string(aexp)*")", "n=10^"*string(nexp)]] =
                @benchmarkable DataStructures.nextreme($fastord, $n, $a)
        end
    end
end

end  # module

BenchHeap.suite

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Stheno: AbstractGP

"""
    GPARegressor(gps::Vector{<:AbstractGP}; replaces=false, imputes=false)

Construct a GPARegressor object with the vector of Gaussian processes `gps`.
"""
struct GPARegressor{Tgps<:Vector{<:AbstractGP}}
    gps::Tgps
    replaces::Bool
    imputes::Bool
end
function GPARegressor(gps::Vector{<:AbstractGP}; replaces=false, imputes=false)
    return GPARegressor(gps, replaces, imputes)
end

replaces(gpar::GPARegressor) = gpar.replaces
imputes(gpar::GPARegressor) = gpar.imputes

"""
    FiniteGPARegressor{Tgpar<:GPARegressor, TX}

Finite-dimensional GPAR object. Lightweight wrapper around input data.
"""
struct FiniteGPARegressor{Tgpar<:GPARegressor, TX}
    gpar::Tgpar
    X::TX
end

"""
    (gpar::GPARegressor)(x::AbstractVector)

Index into `gpar` at `x` to construct a `FiniteGPARegressor`.
"""
(gpar::GPARegressor)(x::AbstractVector) = FiniteGPARegressor(gpar, x)

"""
    rand(rng::AbstractRNG, gpar::FiniteGPARegressor, S::Integer)

Generate `S` i.i.d. draws from `gpar` using `rng`. If `gpar` isn't finite-dimensional, or
equivalently each element of `gpar.gps` isn't a finite-dimensional GP, and error is thrown.
"""
function Stheno.rand(rng::AbstractRNG, finite_gpar::FiniteGPARegressor, S::Integer)
    gpar, X = finite_gpar.gpar, finite_gpar.X
    for d in 1:length(gpar.gps)

    end
    return
end


text:
# Single Models ---------------------------------------------------------------------------------

# splice model and loss functions
start_val(model::Union{Sem, SemFiniteDiff, SemForwardDiff}; kwargs...) =
    start_val(model, model.observed, model.imply, model.diff, model.loss.functions...; kwargs...)

# Fabin 3 starting values for RAM(Symbolic)
start_val(model, observed::Union{SemObsCommon, SemObsMissing}, imply::Union{RAM, RAMSymbolic}, diff, args...; kwargs...) =
    start_fabin3(model; kwargs...)

# Ensemble Models --------------------------------------------------------------------------------
start_val(model::SemEnsemble; kwargs...) = start_simple(model; kwargs...)


text:
using EmpiricalBayes
using Test

using Distributions

marginal_grid = collect(range(-6, stop=6, length=1001));
true_dist = Normal(0,1)
d_true = NormalConvolutionProblem(true_dist, marginal_grid);

m = 7_000
Xs = rand(d_true, m)

estimate(Xs, BrownGreenshtein, PosteriorTarget(PosteriorMeanNumerator(2.0)))
posterior_stats(d_true, PosteriorTarget(PosteriorMeanNumerator(2.0)))

estimate(Xs, BrownGreenshtein, PosteriorTarget(PosteriorMeanNumerator(4.0)))
posterior_stats(d_true, PosteriorTarget(PosteriorMeanNumerator(4.0)))

estimate(Xs, BrownGreenshtein, PosteriorTarget(PosteriorMeanNumerator(3.0)))
posterior_stats(d_true, PosteriorTarget(PosteriorMeanNumerator(3.0)))


text:
mutable struct SurrealFinite <: Surreal
    # mutable means I can set the has value when I need it, instead of
    #   (i) always having to calculate it recursively
    #   (ii) calculating it for even temporary surreals
    shorthand::String
    L::Array{SurrealFinite,1}
    R::Array{SurrealFinite,1}
    h::UInt64 # this is only set the first time the hash function is called
    # constructor should check that L < R
    function SurrealFinite(shorthand::String, L::Array{SurrealFinite}, R::Array{SurrealFinite}, h::UInt64)
        global Count
        Count['c'] += 1
        if length(L) > 1
            L = sort( unique(L), by=x->(x,hash(x)) ) # use hash as tie break in sort so that order is deterministic
        end
        if length(R) > 1
            R = sort( unique(R), by=x->(x,hash(x)) ) # use hash as tie break in sort so that order is deterministic
        end
        # println("L = $L, R = $R")
        # use the fact they are sorted to not do a complete comparison
        # also means that == is much easier than if they were unsorted
        if isempty(L) || isempty(R) || L[end] < R[1]
            return new(shorthand, L, R, h)
        else
            error("Surreal number must have L < R (currently they are $(L) and $(R) )")
        end
    end
end
SurrealFinite(shorthand::String, L::Array, R::Array ) =
    SurrealFinite( shorthand, convert(Array{SurrealFinite},L), convert(Array{SurrealFinite},R), zero(UInt64))
SurrealFinite(L::Array, R::Array ) =
    SurrealFinite( "", convert(Array{SurrealFinite},L), convert(Array{SurrealFinite},R), zero(UInt64))
≀(L::Array, R::Array) = SurrealFinite(L::Array, R::Array )

const SurrealShort  = SurrealFinite
const SurrealDyadic = SurrealFinite

function hash(x::SurrealFinite, h::UInt)
    if x.h == 0
        x.h = hash(x.L, convert(UInt64, 0) ) * hash(x.R, convert(UInt64, 1))
    end
    return hash(x.h, h)
end
function hash(X::Array{SurrealFinite}, h::UInt)
    if isempty(X)
        hash(0, h)
    elseif length(X) == 1
        hash(X[1], h)
    else
        hash(X[1], h) * hash( X[2:end], h) # note, order is technically not important
    end
end
# really slow version: hash(X::Array{SurrealFinite}, h::UInt) = isempty(X) ? hash(0,h) : hash( prod(hash.(X, h )), h)
hash(x::SurrealFinite) = hash(x, zero(UInt64) )
hash(x::SurrealFinite, h::Integer) = hash(x, convert(UInt64, h) )
hash(X::Array{SurrealFinite}) = hash(X, zero(UInt64) )

# global dictionaries to avoid repeated calculations of the same things
#   in particular to speed up calculations, but also to ensure that resulting
#   structures are actually DAGs, i.e., pointers to the same "number" point at the same bit of memory
# i guess eventually these should be replaced by let-closures: https://docs.julialang.org/en/v0.6.2/manual/variables-and-scoping/
const ExistingSurreals = Dict{UInt64, SurrealFinite}()
const ExistingConversions = Dict{UInt64,Rational}()
const ExistingCanonicals = Dict{Rational,UInt64}()
# ExistingProducts = Dict{SurrealFinite, Dict{SurrealFinite,SurrealFinite}}()
const ExistingProducts   = Dict{UInt64, Dict{UInt64,UInt64}}()
const ExistingSums       = Dict{UInt64, Dict{UInt64,UInt64}}()
const ExistingNegations  = Dict{UInt64, UInt64}()
const Count = Dict{Char, Integer}('+'=>0, '*'=>0, '-'=>0, 'c'=>0, '='=>0, '≦'=>0)

# function size(d::Dict)
#     [length(d[k]) for k in sort(collect(keys(d)))]
# end
function reset!(d::Dict)
    for k in keys(d)
        d[k] = 0
    end
end

function clearcache()
    global ExistingSurreals
    global ExistingConversions
    global ExistingCanonicals
    global ExistingProducts
    global ExistingSums
    global ExistingNegations
    global Count
    empty!(ExistingSurreals)
    empty!(ExistingConversions)
    empty!(ExistingCanonicals)
    empty!(ExistingProducts)
    empty!(ExistingSums)
    empty!(ExistingNegations)
    reset!(Count)
    return 1
end

function convert(::Type{SurrealFinite}, n::Int )
    global ExistingSurreals
    global ExistingCanonicals
    global SurrealZero
    global SurrealOne
    if n==0
        result = SurrealZero
    elseif n==1
        result = SurrealOne
    elseif haskey(ExistingCanonicals, Rational(n))
        return ExistingSurreals[ ExistingCanonicals[Rational(n)] ]
    elseif n>1
        result = SurrealFinite(string(n), [convert(SurrealFinite, n-1)], ϕ )
    else
        result = SurrealFinite(string(n), ϕ, [convert(SurrealFinite, n+1)] )
    end
    hr = hash(result)
    if haskey(ExistingSurreals, hr)
       result = ExistingSurreals[hr] # don't double up on memory
    else
       ExistingSurreals[hr] = result
    end
    ExistingCanonicals[Rational(n)] = hr
    return result
    # should check to make sure abs(n) is not too big, i.e., causes too much recursion
end
function convert(::Type{SurrealFinite}, r::Rational )
    global ExistingSurreals
    global ExistingCanonicals
    if haskey(ExistingCanonicals, r)
        return ExistingSurreals[ ExistingCanonicals[r] ]
    elseif isinteger(r)
        result = convert(SurrealFinite, convert(Int64,r) )
    elseif ispow2(r.den)
        # Non-integer dyadic numbers
        n = convert(Int64, round( log2(r.den) ))
        if abs(r) + n < 60 # 60 is a bit arbitrary -- could experiment more on real limits
            result = SurrealFinite("$r", [convert(SurrealFinite, r - 1//2^n)], [convert(SurrealFinite, r + 1//2^n)] )
        else
            error("Generation too large")
        end
    else
        error("we can't do these yet as they require infinite sets")
        # return convert(SurrealFinite, r.num) // convert(SurrealFinite, r.den)
    end
    hr = hash(result)
    if haskey(ExistingSurreals, hr)
       result = ExistingSurreals[hr] # don't double up on memory
    else
       ExistingSurreals[hr] = result
    end
    ExistingCanonicals[r] = hr
    return result
end
function convert(::Type{SurrealFinite}, f::Float64 )
    if isinteger(f)
        return convert(SurrealFinite, convert(Int64,f) )
    elseif isfinite(f)
        # WARNING this could be very slow and memory hungry as exact representations could
        #    be up to 2^52 generations out, so can even cause seg fault, e.g., for 1122342342.23422522
        s = sign(f)
        c = significand(f)
        q = exponent(f)
        @static if VERSION < v"0.7.0"
            t = bits(f)[13:end]
        else
            t = bitstring(f)[13:end]
        end
        p1 = parse(Int, "0b1" * t)
        p2 = p1 * 2.0^(-52 + q)
        r = convert(Rational,s) * p1 // 2^(52 - q)  # rational representation of the Float
        if abs(r.num/r.den) + log2(abs(r.den)) < 60
            # fractional case
            return convert(SurrealFinite, r)
        else
            error("generation too large")
        end
    else
        @static if VERSION < v"0.7.0"
            throw(DomainError())
        else
            throw(DomainError(f))
        end
    end
end
dali(x) = convert(SurrealFinite, x)
SurrealFinite(x) = convert(SurrealFinite, x)

function convert(::Type{Rational}, s::SurrealFinite )
    global ExistingConversions
    h = hash(s)
    if haskey(ExistingConversions, h)
        return ExistingConversions[h]
    else
        if s ≅ zero(s)
            result = 0 // 1
        elseif s ≅ one(s)
            result = 1 // 1
        elseif s < zero(s)
            result = -convert(Rational, -s)
        elseif (sf = floor(Integer, s)) ≅ s
            result = Rational( sf )
        else # 0 < x < 1
            # do a binary search from top down, first valid is the simplest
            xl = isempty(s.L) ? convert(SurrealFinite, sf)   : maximum(s.L)
            xr = isempty(s.R) ? convert(SurrealFinite, sf+1) : minimum(s.R)
            not_end = true
            k = 0
            a = sf
            b = sf+1
            while not_end && k < 24 # 24 is arbitrary, but it would be painful to go lower
                k += 1
                d = (b+a) // 2
                # print("k=$k, a=$a, b=$b, d=$d \n")

                c = convert(SurrealFinite,  d)
                if xl < c < xr
                    result = Rational( d )
                    not_end = false
                elseif c <= xl
                    a = d
                elseif c >= xr
                    b = d
                else
                    error("this case should not happen")
                end
            end
        end
    end
    ExistingConversions[h] = result
    return result
end

# some catch alls
convert(::Type{T}, s::SurrealFinite ) where {T <: AbstractFloat} = convert(T, convert(Rational, s) )
convert(::Type{T}, s::SurrealFinite ) where {T <: Integer} = convert(T, convert(Rational, s) )
convert(::Type{Rational{T}}, s::SurrealFinite ) where {T <: Integer} = convert(Rational{T}, convert(Rational, s) )
AbstractFloat(x::SurrealFinite) = convert(AbstractFloat, x)
# Int64(x::SurrealFinite) = convert(Int64, x)

function convert(::Type{String}, s::SurrealFinite )
    # try to work out a nice way to print it
    if isinteger(s)
        return string(convert(Int, s))
    else
        r = convert(Rational, s)
        return string(r.num) * "/" * string(r.den)
    end
end

# promote all numbers to surreals for calculations
promote_rule(::Type{T}, ::Type{SurrealFinite}) where {T<:Real} = SurrealFinite

@static if VERSION < v"0.7.0"
    const ϕ = Array{SurrealFinite,1}(0) # empty array of SurrealFinites
    const ∅ = ϕ
else
    const ϕ = Array{SurrealFinite,1}(undef,0) # empty array of SurrealFinites
    const ∅ = ϕ
end
const SurrealZero = SurrealFinite("0", ϕ, ϕ )
const SurrealOne  = SurrealFinite("1", [ zero(SurrealFinite) ], ϕ )
const SurrealMinusOne  = SurrealFinite("-1", ϕ, [ zero(SurrealFinite) ] )
const SurrealTwo  = SurrealFinite("2", [SurrealOne], ϕ )
const SurrealThree  = SurrealFinite("3", [SurrealTwo], ϕ )
const SurrealMinusTwo  = SurrealFinite("-2", ϕ, [SurrealOne] )
zero(::SurrealFinite) = SurrealZero # always use the same zero
one(::SurrealFinite)  = SurrealOne  # always use the same one
# ↑ = one(SurrealFinite)  # this causes an error???
# ↓ = -one(SurrealFinite)

# relations
#   these are written in terms of the definition, but could
#   rewrite in terms of max/min to make marginally faster
#   or in terms of set operations to make more succinct
function leq(x::SurrealFinite, y::SurrealFinite, processed_list::Dict{UInt64,Bool})
#    for t in x.L
#        if y <= t
#            return false
#        end
#    end
#    for t in y.R
#        if t <= x
#            return false
#        end
#    end
    # global Count
    # Count['≦'] += 1

    k = hash(x,0) * hash(y,1)
    if !haskey(processed_list,k)
        # if !isempty(x.L) && y <= x.L[end]
        if !isempty(x.L) && leq(y, x.L[end],  processed_list)
            processed_list[k] = false
        # elseif !isempty(y.R) && x >= y.R[1]
        elseif !isempty(y.R) && leq(y.R[1], x,  processed_list)
            processed_list[k] = false
        else
            processed_list[k] = true
        end
    end
    return processed_list[k]
end
<=(x::SurrealFinite, y::SurrealFinite) = leq(x, y, Dict{UInt64,Bool}())
<(x::SurrealFinite, y::SurrealFinite) = x<=y && !(y<=x)
# ===(x::SurrealFinite, y::SurrealFinite) = x<=y && y<x # causes an error
#   === is 'egal', and hardcoded for mutables to test they are same object in memory
≅(x::SurrealFinite, y::SurrealFinite) = x<=y && y<=x
≅(x::Real, y::Real) = ≅(promote(x,y)...)
≇(x::SurrealFinite, y::SurrealFinite) = !( x ≅ y )
≇(x::Real, y::Real) = ≇(promote(x,y)...)
# ==(x::SurrealFinite, y::SurrealFinite) = size(x.L) == size(y.L) &&
#                                         size(x.R) == size(y.R) &&
#                                         all(x.L .== y.L) &&
#                                         all(x.R .== y.R)

# much faster version because doesn't evaluate the whole array every time
#   potential bug because sort order of equivalent values is not defined
function equals(x::SurrealFinite, y::SurrealFinite, processed_list::Dict{UInt64,Bool})
    # global Count
    # Count['='] += 1
    k = hash(x,0) * hash(y,1)
    if !haskey(processed_list,k)
        if length(x.L) != length(y.L) || length(x.R) != length(y.R)
            processed_list[k] = false
            return processed_list[k]
        end
        for i=1:length(x.L)
#             if x.L[i] != y.L[i]
            if !equals(x.L[i], y.L[i], processed_list)
                processed_list[k] = false
                return processed_list[k]
            end
        end
        for i=1:length(x.R)
#            if x.R[i] != y.R[i]
            if !equals(x.R[i], y.R[i], processed_list)
                processed_list[k] = false
                return processed_list[k]
            end
        end
        processed_list[k] = true
    end
    return processed_list[k]
end
==(x::SurrealFinite, y::SurrealFinite) = x===y || equals(x, y, Dict{UInt64,Bool}())
iszero(x::SurrealFinite) = isempty(x.L) && isempty(x.R)

# comparisons between sets (i.e., arrays) are all-to-all, so
#   (1) don't have to be the same size
#   (2) the < is not exactly the same as the < defined above
function <=(X::Array{SurrealFinite}, Y::Array{SurrealFinite} )
    if isempty(X) || isempty(Y)
        return true
    else
        for x in X
            for y in Y
                if !(x <= y)
                    return false
                end
            end
        end
        return true
    end
end
function <(X::Array{SurrealFinite}, Y::Array{SurrealFinite} )
    if isempty(X) || isempty(Y)
        return true
    else
        for x in X
            for y in Y
                if !( x < y )
                    return false
                end
            end
        end
        return true
    end
end

# unary operators
function -(x::SurrealFinite)
    global ExistingSurreals
    global ExistingNegations
    global Count
    hx = hash(x) # build the dictionary in terms of hashs, because it is used quite a bit
    if iszero(x)
        result = x
    elseif haskey(ExistingNegations, hx)
        return ExistingSurreals[ ExistingNegations[hx] ]
    elseif isempty(x.shorthand)
        result = SurrealFinite("", -x.R, -x.L )
    elseif x.shorthand[1] == '-'
        result = SurrealFinite(x.shorthand[2:end], -x.R, -x.L )
    else
        result = SurrealFinite("-"*x.shorthand, -x.R, -x.L )
    end
    hr = hash(result)
    if haskey(ExistingSurreals, hr)
        result = ExistingSurreals[hr] # don't double up on memory
    else
        ExistingSurreals[hr] = result
    end
    if !haskey(ExistingSurreals, hx)
        ExistingSurreals[hx] = x
    end
    ExistingNegations[hx] = hr
    ExistingNegations[hr] = hx
    Count['-'] += 1
    return result
end

# somewhat incomplete (for reasons described), slightly a cheat, and potentially very slow
function /(x::SurrealFinite, y::SurrealFinite)
    xr = convert(Rational, x)
    yr = convert(Rational, y)
    if y ≅ zero(y)
        error(InexactError)
    elseif y ≅ 2
        return x * convert(SurrealFinite,  1 // 2)
    elseif isinteger(y) && ispow2(yr.num)
        return x * convert(SurrealFinite,  1 // yr)
    elseif isinteger(xr.num/yr)
        return convert(SurrealFinite,  (xr.num/yr) // xr.den  )
    else
        error(InexactError)
    end
end

# binary operators
function +(x::SurrealFinite, y::SurrealFinite)
    global ExistingSurreals
    global ExistingSums
    global Count
    hx = hash(x) # build the dictionary in terms of hashs, because it is used quite a bit
    hy = hash(y) #   and these amortise the cost of initial calculation of hashs
    if haskey(ExistingSums, hx) && haskey(ExistingSums[hx], hy)
        return ExistingSurreals[ ExistingSums[hx][hy] ]
    elseif iszero(x)
        result = y
    elseif iszero(y)
        result = x
    else
#       result = SurrealFinite( [x.L .+ y; x .+ y.L],
        #                               [x.R .+ y; x .+ y.R] )
        shorthand = "($(x.shorthand) + $(y.shorthand))"
        result = SurrealFinite(shorthand,
                               [[x + y for x in x.L]; [x + y for y in y.L]],
                               [[x + y for x in x.R]; [x + y for y in y.R]] )
    end
    hr = hash(result)
    # println(" $x+$y:   $hr   $result")
    if haskey(ExistingSurreals, hr)
       result = ExistingSurreals[hr] # don't double up on memory
    else
       ExistingSurreals[hr] = result
    end
    # ExistingSurreals[hr] = result  # swap to just using this if we want to see all the ways an addition can given the same result

    if !haskey(ExistingSums, hx)
        ExistingSums[hx] = Dict{UInt64,UInt64}()
    end
    if !haskey(ExistingSums, hy)
        ExistingSums[hy] = Dict{UInt64,UInt64}()
    end
    ExistingSums[hx][hy] = hr
    ExistingSums[hy][hx] = hr
    Count['+'] += 1
    return result
end

# +(x::SurrealFinite, y::SurrealFinite) = SurrealFinite([x.L .+ y; x .+ y.L],
#                                                       [x.R .+ y; x .+ y.R] )

# can't do like this because of empty arrays I think
#+(x::SurrealFinite, y::SurrealFinite) = SurrealFinite([x.L + y.L],
#                                                      [x.R + y.R] )
+(X::Array{SurrealFinite}, Y::Array{SurrealFinite}) = vec([s+t for s in X, t in Y])

-(x::SurrealFinite, y::SurrealFinite) = x + -y
-(X::Array{SurrealFinite}, Y::Array{SurrealFinite}) = X + -Y

function *(x::SurrealFinite, y::SurrealFinite)
    global ExistingSurreals
    global ExistingProducts
    global Count
    hx = hash(x)
    hy = hash(y)
    if haskey(ExistingProducts, hx) && haskey(ExistingProducts[hx], hy)
        return ExistingSurreals[ ExistingProducts[hx][hy] ]
    elseif iszero(x) || iszero(y)
        result = zero(x)
#     elseif x == 1
#        result = y
#    elseif y == 1
#        result = x
    else
        # tmp1 = vec([s*y + x*t - s*t for s in x.L, t in y.L])
        # tmp2 = vec([s*y + x*t - s*t for s in x.R, t in y.R])
        # tmp3 = vec([s*y + x*t - s*t for s in x.L, t in y.R])
        # tmp4 = vec([s*y + x*t - s*t for s in x.R, t in y.L])
        # L = [ tmp1; tmp2 ]
        # R = [ tmp3; tmp4 ]
        shorthand = "($(x.shorthand) * $(y.shorthand))"
        L = [vec([s*y + x*t - s*t for s in x.L, t in y.L]);
             vec([s*y + x*t - s*t for s in x.R, t in y.R])]
        R = [vec([s*y + x*t - s*t for s in x.L, t in y.R]);
             vec([s*y + x*t - s*t for s in x.R, t in y.L])]
        # L = [subtimes( x, y, x.L, y.L ); subtimes( x, y, x.R, y.R )]
        # R = [subtimes( x, y, x.L, y.R ); subtimes( x, y, x.R, y.L )]
        result = SurrealFinite(shorthand, L, R)
    end
    if !haskey(ExistingProducts, hx)
        ExistingProducts[hx] = Dict{UInt64,UInt64}()
    end
    if !haskey(ExistingProducts, hy)
        ExistingProducts[hy] = Dict{UInt64,UInt64}()
    end
    hr = hash(result)
    if haskey(ExistingSurreals, hr)
        result = ExistingSurreals[hr] # don't double up on memory
    else
        ExistingSurreals[hr] = result
    end
    ExistingProducts[hx][hy] = hr
    ExistingProducts[hy][hx] = hr
    Count['*'] += 1
    return result
end
#function subtimes( x, y, X::Array{SurrealFinite}, Y::Array{SurrealFinite} )
#    n = length(X)
#    m = length(Y)
#    result = Array{SurrealFinite}(m*n)
#    for i=1:n
#        for j=1:m
#            result[(i-1)*m + j] = X[i]*y + x*Y[j] - X[i]*Y[j]
#        end
#    end
#    return result
#end
*(x::SurrealFinite, Y::Array{SurrealFinite}) = [ x*s for s in Y ]
*(X::Array{SurrealFinite}, y::SurrealFinite) = y*X

#####################################################3
# read/write routines for surreals

# print the first level in full (ignoring the top level shorthand if present)
#   these should be replaced using "expand"
pf(io::IO, x::SurrealFinite) = println(io, "{ ", x.L, " | ", x.R, " }")
pf(x::SurrealFinite) = pf(stdout, x)

"""
    expand(x::SurrealFinite; level=0)

 Writes a surreal as a string with varying levels of expansion.

## Arguments
* `x::SurrealFinite`: the number of elements to expand
* `level=0`: the amount of expansion

    + 0 : write shorthand if it exists, or ``\\{ X_L \\| X_R \\}`` if not
    + 1 : ``\\{ X_L \\| X_R \\}``
    + 2 : expand out ``X_L`` and ``X_R`` recursively

## Examples
```jldoctest
julia> expand( convert(SurrealFinite, 2))
"2"
julia> expand( convert(SurrealFinite, 2); level=1)
"{ 1 | ∅ }"
julia> expand( convert(SurrealFinite, 2); level=2)
"{ { { ∅ | ∅ } | ∅ } | ∅ }"
```
"""
function expand(x::SurrealFinite; level=0)
    if level==0
        s = x.shorthand != "" ? x.shorthand : expand(x; level=1)
    elseif level==1
        tmpL = isempty(x.L) ? "∅" : join(convert.(String, x.L), ',')
        tmpR = isempty(x.R) ? "∅" : join(convert.(String, x.R), ',')
        s = "{ " * tmpL * " | " * tmpR  * " }"
        return s
    elseif level>=1
        return "{ " * expand(x.L;level=level) * " | " * expand(x.R;level=level) * " }"
    end
end
expand(X::Array{SurrealFinite}; level=0) = isempty(X) ? "∅" : join(expand.(X; level=level), ',')

# isnumeric(s::AbstractString) = ismatch(r"^-?\d*.\d*$", s) || ismatch(r"^-?\d*//-?\d*$", s)

# this has to parse a surreal written into a string
function convert(::Type{SurrealFinite}, s::AbstractString )
    # interpret, (i) numbers as canonical, (ii) phi, \phi, ϕ correctly, (iii) structure
    @static if VERSION < v"0.7.0"
        s = replace(s, r"\s+", "") # remove white space
        s = replace(s, r"phi|\\phi|ϕ|∅|\{\}", "") # replace phi or \phi with empty set
    else
        s = replace(s, r"\s+" => "") # remove white space
        s = replace(s, r"phi|\\phi|ϕ|∅|\{\}" => "") # replace phi or \phi with empty set
    end
    # println("     s1 = $s")
    not_end = true
    basic_number = r"\{([^{}|]*)\|([^{|}]*)\}"
    while not_end
        @static if VERSION < v"0.7.0"
            if ismatch(basic_number, s)
                s = replace(s, basic_number, s"SurrealFinite([\1],[\2])") # remove comments
            else
                not_end = false
            end
        else
            if occursin(basic_number, s)
                s = replace(s, basic_number => s"SurrealFinite([\1],[\2])") # remove comments
            else
                not_end = false
            end
        end
    end
    # println("     s2 = $s")
    return eval(Meta.parse(s))
end

# read in the full format
function read(io::IO, ::Type{SurrealFinite}, n::Int=1)
    @static if VERSION < v"0.7.0"
        X = Array{SurrealFinite,1}(n)
    else
        X = Array{SurrealFinite,1}(undef,n)
    end
    k = 1
    while !eof(io) && k<=n
        # println(" $k:   $line")
        @static if VERSION < v"0.7.0"
            line = replace(readline(io), r"#.*", s"") # remove comments
            if ismatch(r"\S", line)
                X[k] = convert(SurrealFinite, line)
                k += 1
            end
        else
            line = replace(readline(io), r"#.*" => s"") # remove comments
            if occursin(r"\S", line)
                X[k] = convert(SurrealFinite, line)
                k += 1
            end
        end
    end
    return X
end
# read(filename::AbstractString, args...) = open(io->read(io, args...), filename)

# write out a string suitable for inclusion into latex docs
function surreal2tex(io::IO, x::SurrealFinite; level=0)
    s = expand(x; level=level)
    @static if VERSION < v"0.7.0"
        s = replace(s, r"\{", " \\{ ")
        s = replace(s, r"ϕ", " \\emptyset ")
        s = replace(s, r"∅", " \\emptyset ")
        s = replace(s, r"\|", " \\mid ")
        s = replace(s, r"\}", " \\} ")
        s = replace(s, r"(\d+)//(\d+)", s"\\frac{\1}{\2}")
    else
        s = replace(s, r"\{" => " \\{ ")
        s = replace(s, r"ϕ" => " \\emptyset ")
        s = replace(s, r"∅" => " \\emptyset ")
        s = replace(s, r"\|" => " \\mid ")
        s = replace(s, r"\}" => " \\} ")
        s = replace(s, r"(\d+)//(\d+)" => s"\\frac{\1}{\2}")
    end
    println(io,s)
end
surreal2tex(x::SurrealFinite; level=0) = surreal2tex(stdout, x; level=level)

# standard show will use shorthand when available
function show(io::IO, x::SurrealFinite)
    if io==stdout && x.shorthand != ""
        @static if VERSION < v"0.7.0"
            print_with_color(:bold, io, x.shorthand ) # could be :red
        else
            printstyled(io, x.shorthand; bold=true)
        end
    elseif x.shorthand != ""
        print(io, x.shorthand )
    else
        # print( io, "<", x.L, ":", x.R, ">")
        print( io, "{ ", x.L, " | ", x.R, " }")
    end
end
# show(io::IO, X::Array{SurrealFinite}) = print(io, "{", join(X, ", "), "}")
function show(io::IO, X::Array{SurrealFinite})
    if isempty(X)
        print(io, "∅")
    else
        print(io, join(X, ", "))
    end
end
# special "canonicalised" output
spf(x::SurrealFinite) = print("{ ", canonicalise.(x.L), " | ", canonicalise.(x.R), " }")

"""
    surreal2dag(x::SurrealFinite)
    surreal2dag(io::IO, x::SurrealFinite)

 Writes a surreal representation as a DAG out in DOT format for drawing using GraphVis,
 and returns the number of nodes in the graph. Returns the number of nodes in the graph.

## Arguments
* `io::IO`: output stream, default is stdout
* `x::SurrealFinite`: the number to write out

## Examples
```jldoctest
julia> surreal2dag(convert(SurrealFinite, 0))
digraph "0.0" {
   node_1 [shape=none,margin=0,label=
         <<TABLE BORDER="0" CELLBORDER="1" CELLSPACING="0" CELLPADDING="4">
         <TR><TD COLSPAN="2">0</TD></TR>
         <TR><TD PORT="L"> ∅ </TD><TD PORT="R"> ∅ </TD></TR>
         </TABLE>>,
         ];
}
1
```
"""
function surreal2dag(io::IO, x::SurrealFinite)
    println(io, "digraph \"", float(x), "\" {")
    k = 1
    SurrealsinDAG = Dict{SurrealFinite,Int}()
    m = surreal2dag_f(io, x, k, SurrealsinDAG)
    println(io, "}")
    return m
end
function surreal2dag_f(io::IO, x::SurrealFinite, k::Integer, SurrealsinDAG)
    m = k
    if !haskey(SurrealsinDAG, x)
        SurrealsinDAG[x] = m
        surreal2node(io, x, k)
        c = 1
        for s in x.L
            if !haskey(SurrealsinDAG, s)
                m += 1
                # println(io, "   node_$k:L -> node_$m;")
                println(io, "   node_$k:\"" *  convert(String, s) * "," * string(c) * "\" -> node_$m [color=\"red3\"];")
                m = surreal2dag_f(io, s, m, SurrealsinDAG)
            else
                println(io, "   node_$k:\"" *  convert(String, s) * "," * string(c)  * "\" -> node_$(SurrealsinDAG[s]) [color=\"red3\"];")
            end
            c += 1
        end
        c = 1
        for s in x.R
            if !haskey(SurrealsinDAG, s)
                m += 1
                # println(io, "   node_$k:R -> node_$m;")
                println(io, "   node_$k:\"" *  convert(String, s) * "," * string(c) * "\" -> node_$m [color=\"blue3\"];")
                m = surreal2dag_f(io, s, m, SurrealsinDAG)
            else
                println(io, "   node_$k:\"" *  convert(String, s) * "," * string(c) * "\" -> node_$(SurrealsinDAG[s]) [color=\"blue3\"];")
            end
            c += 1
        end
    end
    return m
end
surreal2dag(x::SurrealFinite) = surreal2dag(stdout, x)

"""
    surreal2dot(x::SurrealFinite)
    surreal2dot(io::IO, x::SurrealFinite)

 Writes a surreal representation as a tree out in DOT format for drawing using GraphVis,
 and returns the number of nodes in the graph.

## Arguments
* `io::IO`: output stream, default is stdout
* `x::SurrealFinite`: the number to write out

## Examples
```jldoctest
julia> surreal2dot(convert(SurrealFinite, 1))
digraph "1.0" {
   node_1 [shape=none,margin=0,label=
         <<TABLE BORDER="0" CELLBORDER="1" CELLSPACING="0" CELLPADDING="4">
         <TR><TD COLSPAN="2">1</TD></TR>
         <TR><TD PORT="L"> <TABLE BORDER="0" CELLBORDER="0" CELLPADDING="0"><TR><TD PORT="0,1"> 0 </TD>   </TR></TABLE> </TD><TD PORT="R"> ϕ </TD></TR>
         </TABLE>>,
         ];
   node_1:"0,1" -> node_2;
   node_2 [shape=none,margin=0,label=<<B>0</B>>]
}
2
```
"""
function surreal2dot(io::IO, x::SurrealFinite)
    println(io, "digraph \"", float(x), "\" {")
    k = 1
    m = surreal2dot_f(io, x, k)
    println(io, "}")
    return m
end
function surreal2dot_f(io::IO, x::SurrealFinite, k::Integer)
    m = k
    if x == zero(x)
        println(io, "   node_$k [shape=none,margin=0,label=<<B>0</B>>]")
    else
        #if x.shorthand==""
        #    S = convert(String, x)
        #else
        #    S = x.shorthand
        #end
        surreal2node(io, x, k)
        c = 1
        for s in x.L
            m += 1
            # println(io, "   node_$k:L -> node_$m;")
            println(io, "   node_$k:\"" *  convert(String, s) * "," * string(c) * "\" -> node_$m;")
            m = surreal2dot_f(io, s, m)
            c += 1
        end
        c = 1
        for s in x.R
            m += 1
            # println(io, "   node_$k:R -> node_$m;")
            println(io, "   node_$k:\"" *  convert(String, s) * "," * string(c)  * "\" -> node_$m;")
            m = surreal2dot_f(io, s, m)
           c += 1
         end
    end
    return m
end
function surreal2node(io::IO, x::SurrealFinite, k::Integer; extra_args::String="")
    S = convert(String, x)
    if isempty(x.L)
        L = "∅"
    else
        L = "<TABLE BORDER=\"0\" CELLBORDER=\"0\" CELLPADDING=\"0\"><TR>"
        c = 1
        for s in x.L
            tmp = convert(String, s)
            L *= "<TD PORT=\"$tmp," * string(c) * "\"> " * tmp  * " </TD>   "
            c += 1
        end
        L *= "</TR></TABLE>"
    end
    if isempty(x.R)
        R = "∅"
    else
        R = "<TABLE BORDER=\"0\" CELLBORDER=\"0\" CELLPADDING=\"0\"><TR>"
        c = 1
        for s in x.R
            tmp = convert(String, s)
            R *= "<TD PORT=\"$tmp," * string(c) * "\"> " * tmp * " </TD>   "
            c += 1
        end
        R *= "</TR></TABLE>"
    end
    print(io, "   ")
    if k>=0
        label = "$k"
    else
        label = "m$k"
    end
    println(io, """
                node_$label [shape=none,margin=0,label=
                         <<TABLE BORDER=\"0\" CELLBORDER=\"1\" CELLSPACING=\"0\" CELLPADDING=\"4\">
                         <TR><TD COLSPAN=\"2\">$S</TD></TR>
                         <TR><TD PORT=\"L\"> $L </TD><TD PORT=\"R\"> $R </TD></TR>
                         </TABLE>>,$extra_args
                         ];""")
end
surreal2dot(x::SurrealFinite) = surreal2dot(stdout, x)

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

"""
    generation(x::SurrealFinite)

 Finds the birthday of a surreal number, which is 1 + the max of any of its components.

## Arguments
* `x::SurrealFinite`: the number to operate on

## Examples
```jldoctest
julia> generation( convert(SurrealFinite, 1) )
1
```
"""
function generation(x::SurrealFinite)
    if x==zero(x)
        return 0
    else
        return max( maximum( generation.( [x.L; 0]) ),
                    maximum( generation.( [x.R; 0]) )) + 1
    end
end

"""
    canonicalise(s::SurrealFinite)

 Convert a surreal number form into its equivalent canonical form. It performs the
 conversion by converting to a Rational and then back to a Surreal.

## Examples
```jldoctest
julia> convert(SurrealFinite, 1) - convert(SurrealFinite, 1)
{ { ϕ | { ϕ | ϕ } } | { { ϕ | ϕ } | ϕ } }
julia> pf( canonicalise( convert(SurrealFinite, 1) - convert(SurrealFinite, 1) ) )
{ ϕ | ϕ }
```
"""
canonicalise(s::SurrealFinite) = convert(SurrealFinite, convert(Rational, s))
iscanonical(s::SurrealFinite) = canonicalise(s) == s

function unique2!( X::Array{SurrealFinite} )
    # our own unique that is based on ≅
    sort!(X)
    i = 1
    while i < length(X)
        if X[i] ≅ X[i+1]
            splice!( X, i )
        end
        i += 1
    end
end

###### standard math routines ##############################

sign(x::SurrealFinite) = x<zero(x) ? -one(x) : x>zero(x) ? one(x) : zero(x)
# abs(x::SurrealFinite) = x<zero(x) ? -x : x

# subtraction is much slower than comparison, so got rid of old versions
# these aren't purely surreal arithmetic, but everything could be, just would be slower
function floor(T::Type, s::SurrealFinite)
    if s < zero(s)
        if s >= -one(s)
            return -one(T)
        elseif s >= -2
            return convert(T, -2)
        else
            k = 1
            while s < -2^(k+1) && k < 12 # 12 is arbitrary, but it would be painful to go higher
                k += 1
            end
            if k==12
                error("s is too large for the current floor function")
            end
            a = -2^(k+1)
            b = -2^k
            for i=1:k
                d = (a+b) / 2
                if s < d
                    b = d
                else
                    a = d
                end
            end
            return convert(T, a) # N.B. returns canonical form of floor
        end
    elseif s < one(s)
        return zero(T)
    elseif s < 2
        return one(T)
    else
        # start with geometric search to bound the number
        k = 1
        while s >= 2^(k+1) && k < 12 # 12 is arbitrary, but it would be painful to go higher
            k += 1
        end
        if k==12
            error("s is too large for the current floor function")
        end
        # now a binary search to narrow it down
        a = 2^k
        b = 2^(k+1)
        for i=1:k
            d = (a+b) / 2
            if s < d
                b = d
            else
                a = d
            end
         end
        return convert(T, a) # N.B. returns canonical form of floor
    end
end
function ceil(T::Type, s::SurrealFinite)
    if isinteger(s)
        return floor(T,s)
    else
        return floor(T,s) + 1
    end
end

# implicit rounding mode is 'RoundNearestTiesUp'
#   to bbe consistent, should do the other rounding modes, and a precision, but the latter is hard
function round(T::Type, s::SurrealFinite)
    return floor(T, s + 1//2)
end

trunc(T::Type, s::SurrealFinite) = s>=0 ? floor(T,s) : -floor(T,-s)

# simple versions
floor(s::SurrealFinite) = floor(SurrealFinite, s)
ceil(s::SurrealFinite) = ceil(SurrealFinite, s)
round(s::SurrealFinite) = round(SurrealFinite, s)
trunc(s::SurrealFinite) = trunc(SurrealFinite, s)

# this should still be rewritten in terms of searches
function mod(s::SurrealFinite, n::SurrealFinite)
    if !isinteger(n)
        error("n should be an integer")
    end
    if s ≅ zero(s)
        return s
    elseif s < zero(s)
        return mod(s + n, n)
    elseif s >= n
        return mod(s - n, n)
    else
        return s
    end
end

function isinteger(s::SurrealFinite)
    if s ≅ floor(s)
        return true
    else
        return false
    end
end

isdivisible(s::SurrealFinite, n::SurrealFinite) = isinteger(s) ? mod(s,n) ≅ zero(s) : false
isodd(s::SurrealFinite)  = isinteger(s) ? !isdivisible(s, convert(SurrealFinite,2) ) : false
iseven(s::SurrealFinite) = isinteger(s) ?  isdivisible(s, convert(SurrealFinite,2) ) : false
# ispow2

isinf(s::SurrealFinite) = false
isnan(s::SurrealFinite) = false
isfinite(s::SurrealFinite) = true

parents(s::SurrealFinite) = [s.L; s.R]
function ⪯(s::SurrealFinite, t::SurrealFinite)# is s an ancestor of t, or equal to it
    if s == t
        return true
    elseif iszero(s)
        return true
    else
        for p in parents(t)
            if s⪯p
                return true
            end
        end
    end
    return false
end
isancestor(s::SurrealFinite, t::SurrealFinite) = s ⪯ t && s!=t
≺(s::SurrealFinite, t::SurrealFinite) = isancestor(s, t)
≻(s::SurrealFinite, t::SurrealFinite) = t ≺ s
⪰(s::SurrealFinite, t::SurrealFinite) = t ⪯ s

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

# extra analysis functions

# structure to hold stats on a surreal form
struct SurrealDAGstats
    nodes::Integer # nodes in its DAG
    tree_nodes::Integer # nodes in a tree-based representation of the graph
    edges::Integer # edges in its DAG
    generation::Integer # generation/birthday of the surreal
    longest_path::Array{SurrealFinite,1}
    paths::Int128 # number of paths from source to sink
    value::Rational
    minval::Rational
    maxval::Rational
end

#    LP=true means keep track of a longest path
#    V =true means keep track of values
# have these as switches, because they are expensive to calculate on large DAGs
function dag_stats(s::SurrealFinite, processed_list::Dict{SurrealFinite,SurrealDAGstats}; LP=false, V=true)
    if s == zero(s)
        nodes = 1
        tree_nodes = 1
        edges = 0
        generation = 0
        longest_path = LP ? [zero(s)] : ϕ # important this be \phi and not an arbitrary empty array
        paths = 1
        value = 0
        minval = 0
        maxval = 0
    else
        P = parents(s)
        nodes = 1
        tree_nodes = 1
        edges = length(P)
        generation = 0
        longest_path = LP ? [zero(s)] : ϕ
        paths = 0
        value = V ? convert(Rational, s) : 0
        minval = value
        maxval = value
        for p in P
            if !haskey(processed_list, p)
                (stats,l) = dag_stats(p, processed_list; LP=LP, V=V)
                nodes += stats.nodes
                edges += stats.edges
                minval = min(minval, stats.minval)
                maxval = max(maxval, stats.maxval)
            else
                stats = processed_list[p]
            end
            if stats.generation > generation
                generation = stats.generation
                longest_path = LP ? copy(stats.longest_path) : ϕ
            end
            tree_nodes += stats.tree_nodes
            paths += stats.paths
        end
        generation += 1
        longest_path = LP ? [longest_path; s] : ϕ
    end
    stats = SurrealDAGstats(nodes, tree_nodes, edges, generation, longest_path, paths, value, minval, maxval)
    processed_list[s] = stats # could make this a reverse list, ...
    return ( stats, processed_list )
end
dag_stats(s::SurrealFinite; LP=false, V=true) = dag_stats(s, Dict{SurrealFinite,SurrealDAGstats}(); LP=LP, V=V)[1]
nodes(s::SurrealFinite) = dag_stats(s).nodes
edges(s::SurrealFinite) = dag_stats(s).edges
paths(s::SurrealFinite) = dag_stats(s).paths
tree_nodes(s::SurrealFinite) = dag_stats(s).tree_nodes
function breadth(s::SurrealFinite)
    d = dag_stats(s)
    return d.maxval - d.minval
end
function width(s::SurrealFinite)
    # not yet defined exactly what I mean here --
    #    maybe maximum number of nodes for a fixed generation?
    0
end

# for debugging uniqueness of nodes in memory
#    check whether an entry in processed list is unique according to hash tables
#    the values in U should all be 1, but if we do say 2*4, without the caching in +,
#       we'll get a '2' because identical forms are being allocated multiple times in memory
function increment( U::Dict{UInt64, Int64}, h::UInt64 )
    if haskey(U, h)
        U[h] += 1
    else
        U[h] = 1
    end
end
function uniqueness(s::SurrealFinite,
                    P::Dict{UInt64, UInt64},
                    Q::Dict{UInt64, SurrealFinite},
                    U::Dict{UInt64, Int64},
                    V::Dict{UInt64, SurrealFinite} )
    h = hash(s)
    increment(U, h)
    @static if VERSION < v"0.7.0"
        os = object_id(s)
    else
        os = objectid(s)
    end
    P[ os ] = h
    Q[ os ] = s
    V[h] = s
    for p in parents(s)
        @static if VERSION < v"0.7.0"
            op = object_id(p)
        else
            op = objectid(p)
        end
        if !haskey( P, op )
            uniqueness(p, P, Q, U, V )
        end
    end
    return (P, Q, U, V )
end
uniqueness(s::SurrealFinite) = uniqueness(s,
                                          Dict{UInt64, UInt64}(),
                                          Dict{UInt64, SurrealFinite}(),
                                          Dict{UInt64, Int64}(),
                                          Dict{UInt64, SurrealFinite}() )
uniqueness_max( U::Dict{UInt64, Int64} ) = maximum( values(U) )
function uniqueness_failure(  U::Dict{UInt64, Int64}, V::Dict{UInt64, SurrealFinite} )
    U2 = filter( (k,v) -> v>1 , U) # not Julia 1.0 compliant
    V2 = [V[i] for i in keys(U2) ]
    return V2
end

###########################################################
### work towards a formal DAG representation and tools

function surrealDAG(s::SurrealFinite)
    dag_dict = dag_stats(s, Dict{SurrealFinite,SurrealDAGstats}())[2]
    # dag = Array{SurrealDAGstats}(length(dag_dict)) # sort from simplest to least simple
    dag = Array{SurrealFinite}(length(dag_dict)) # sort from simplest to least simple
    i = 1
    for x in sort( collect(keys(dag_dict)); by = x->(dag_dict[x].generation,dag_dict[x].value,dag_dict[x].nodes) )
        dag[i] = x
        i += 1
    end
    return dag # an array containing sorted list of all nodes of the surreal form
end
function dag_add(x::SurrealFinite, y::SurrealFinite)
    # bit of a cheat -- not reworking links, just relying on sub-sums
    # just intended to check the sums and cartesian products
    dag_x = surrealDAG(x)
    n_x = length(dag_x)
    dag_y = surrealDAG(y)
    n_y = length(dag_y)
    dag_z = Array{SurrealFinite}( n_x * n_y )
    for i=1:n_x
        for j=1:n_y
            dag_z[(i-1)*n_y + j] = dag_x[i] + dag_y[j]
        end
    end
    return unique(dag_z)
end


text:
"""
    Module for computing the Kmeans approximation for finding inducing points, it is based on the k-means++ algorithm
"""
module KMeansModule

using Distributions
using StatsBase
using LinearAlgebra
using Clustering
using AugmentedGaussianProcesses.KernelModule

export KMeansInducingPoints

#Return K inducing points from X, m being the number of Markov iterations for the seeding
function KMeansInducingPoints(X::AbstractArray{T,N},nC::Integer;nMarkov::Integer=10,kweights::Vector{<:Real}=[0.0]) where {T<:Real,N}
    C = copy(transpose(KmeansSeed(X,nC,nMarkov)))
    if kweights!=[0.0]
        Clustering.kmeans!(X',C,weights=kweights,tol=1e-3)
    else
        Clustering.kmeans!(X',C,tol=1e-3)
    end
    return copy(transpose(C))
end

"""Find good approximates seeds for the k-means++ algorithm, `X` is the data, `nC` is the number of centers and `nMarkov` is the number of Markov Chain samples before stopping. Implementation of [Fast and efficient seeding for KMeans based on [Fast and Provably Good Seeding for k-Means](https://las.inf.ethz.ch/files/bachem16fast.pdf)"""
function KmeansSeed(X::AbstractMatrix{T},nC::Integer,nMarkov::Integer) where {T<:Real} #X is the data, nC the number of centers wanted, m the number of Markov iterations
  NSamples = size(X,1)
  #Preprocessing, sample first random center
  init = rand(1:NSamples)
  C = zeros(T,nC,size(X,2))
  C[1,:] .= X[init,:]
  q = zeros(T,NSamples)
  @inbounds for i in 1:NSamples
    @views q[i] = 0.5*norm(X[i,:].-C[1])^2
  end
  sumq = sum(q)
  q = Weights(q/sumq .+ 1.0/(2*NSamples),1)
  uniform = Distributions.Uniform(0,1)
  x = zeros(T,size(X,2))
  y = similar(x)
  @inbounds for i in 2:nC
    x .= vec(X[StatsBase.sample(1:NSamples,q,1),:]) #weighted sampling,
    mindist = mindistance(x,C,i-1)
    for j in 2:nMarkov
      y .= vec(X[StatsBase.sample(q),:]) #weighted sampling
      dist = mindistance(y,C,i-1)
      if (dist/mindist > rand(uniform))
        x .= y;  mindist = dist
      end
    end
    C[i,:].=x
  end
  return C
end

"""Compute the minimum distance between a vector `x` and a collection of centers `C` given the number of computed centers `nC`"""
function mindistance(x::AbstractVector{T},C::AbstractMatrix{T},nC::Integer) where {T}#Point to look for, collection of centers, number of centers computed
  mindist = T(Inf)
  @inbounds for i in 1:nC
    mindist = min(norm(x-C[i,:])^2,mindist)
  end
  return mindist
end

end #module


text:
baremodule TaskDAGRecorders

export @sync, Threads

macro sync end
macro spawn end

function sample end

const NTHREADS = Core.Ref(8)
nthreads() = NTHREADS[]

"""
    TaskDAGRecorders.Internal

This module contains the internal details.
"""
module Internal

import ..TaskDAGRecorders: @sync, @spawn
using ..TaskDAGRecorders: TaskDAGRecorders

include("utils.jl")
include("recording.jl")

end  # module Internal

"""
    TaskDAGRecorders.Threads

A shim of `Base.Threads` that provides a subset of `Threads` for instrumenting
`@spawn` and `@sync`.
"""
baremodule Threads
using ..TaskDAGRecorders: @sync, @spawn, nthreads
end  # baremodule Threads

"""
    TaskDAGRecorders.NTHREADS::Ref{Int}

The number of threads that `TaskDAGRecorders.Threads` pretends to have.
"""
NTHREADS

end  # baremodule TaskDAGRecorders


text:
function codegen(::Type{Val{:iterate}}, ::Type{NUTS}, job::BasicMCJob)
  result::Expr
  whilebody = []
  ifwhilebody = []
  burninbody = []
  ifburninbody = []
  body = []
  dindex::Integer

  vform = variate_form(job.pstate)
  if vform != Univariate && vform != Multivariate
    error("Only univariate or multivariate parameter states allowed in NUTS code generation")
  end

  # push!(body, :(_job.sstate.count += 1))

  if job.tuner.verbose
    push!(body, :(_job.sstate.tune.proposed += 1))
  end

  if vform == Univariate
    push!(body, :(_job.sstate.momentum = randn()))
  elseif vform == Multivariate
    push!(body, :(_job.sstate.momentum[:] = randn(_job.pstate.size)))
  end

  push!(body, :(_job.sstate.oldhamiltonian = hamiltonian(_job.pstate.logtarget, _job.sstate.momentum)))

  if vform == Univariate
    push!(body, :(_job.sstate.pstateplus.value = _job.pstate.value))
    push!(body, :(_job.sstate.pstateplus.gradlogtarget = _job.pstate.gradlogtarget))
    push!(body, :(_job.sstate.momentumplus = _job.sstate.momentum))
    push!(body, :(_job.sstate.pstateminus.value = _job.pstate.value))
    push!(body, :(_job.sstate.pstateminus.gradlogtarget = _job.pstate.gradlogtarget))
    push!(body, :(_job.sstate.momentumminus = _job.sstate.momentum))
  elseif vform == Multivariate
    push!(body, :(_job.sstate.pstateplus.value = copy(_job.pstate.value)))
    push!(body, :(_job.sstate.pstateplus.gradlogtarget = copy(_job.pstate.gradlogtarget)))
    push!(body, :(_job.sstate.momentumplus = copy(_job.sstate.momentum)))
    push!(body, :(_job.sstate.pstateminus.value = copy(_job.pstate.value)))
    push!(body, :(_job.sstate.pstateminus.gradlogtarget = copy(_job.pstate.gradlogtarget)))
    push!(body, :(_job.sstate.momentumminus = copy(_job.sstate.momentum)))
  end

  push!(body, :(_job.sstate.j = 0))
  push!(body, :(_job.sstate.n = 1))
  push!(body, :(_job.sstate.s = true))
  if in(:accept, job.outopts[:diagnostics])
    push!(body, :(_job.sstate.update = false))
  end

  push!(body, :(_job.sstate.u = log(rand())+_job.sstate.oldhamiltonian))

  push!(whilebody, :(_job.sstate.v = rand(Bool) ? 1 : -1))

  push!(whilebody, :(
    if _job.sstate.v == -1
      _job.sstate.pstateminus,
      _job.sstate.momentumminus,
      _,
      _,
      _job.sstate.pstateprime,
      _job.sstate.nprime,
      _job.sstate.sprime =
        build_tree!(
          _job.sstate,
          _job.sstate.pstateminus,
          _job.sstate.momentumminus,
          _job.sstate.u,
          _job.sstate.v,
          _job.sstate.j,
          _job.sstate.tune.step,
          _job.parameter.logtarget!,
          _job.parameter.gradlogtarget!,
          _job.sampler
        )
    else
      _,
      _,
      _job.sstate.pstateplus,
      _job.sstate.momentumplus,
      _job.sstate.pstateprime,
      _job.sstate.nprime,
      _job.sstate.sprime =
        build_tree!(
          _job.sstate,
          _job.sstate.pstateplus,
          _job.sstate.momentumplus,
          _job.sstate.u,
          _job.sstate.v,
          _job.sstate.j,
          _job.sstate.tune.step,
          _job.parameter.logtarget!,
          _job.parameter.gradlogtarget!,
          _job.sampler
        )
    end
  ))

  if vform == Univariate
    push!(ifwhilebody, :(_job.pstate.value = _job.sstate.pstateprime.value))
    push!(ifwhilebody, :(_job.pstate.gradlogtarget = _job.sstate.pstateprime.gradlogtarget))
    if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
      push!(ifwhilebody, :(_job.pstate.gradloglikelihood = _job.sstate.pstateprime.gradloglikelihood))
    end
    if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
      push!(ifwhilebody, :(_job.pstate.gradlogprior = _job.sstate.pstateprime.gradlogprior))
    end
  elseif vform == Multivariate
    push!(ifwhilebody, :(_job.pstate.value = copy(_job.sstate.pstateprime.value)))
    push!(ifwhilebody, :(_job.pstate.gradlogtarget = copy(_job.sstate.pstateprime.gradlogtarget)))
    if in(:gradloglikelihood, job.outopts[:monitor]) && job.parameter.gradloglikelihood! != nothing
      push!(ifwhilebody, :(_job.pstate.gradloglikelihood = copy(_job.sstate.pstateprime.gradloglikelihood)))
    end
    if in(:gradlogprior, job.outopts[:monitor]) && job.parameter.gradlogprior! != nothing
      push!(ifwhilebody, :(_job.pstate.gradlogprior = copy(_job.sstate.pstateprime.gradlogprior)))
    end
  end
  push!(ifwhilebody, :(_job.pstate.logtarget = _job.sstate.pstateprime.logtarget))
  if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
    push!(ifwhilebody, :(_job.pstate.loglikelihood = _job.sstate.pstateprime.loglikelihood))
  end
  if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
    push!(ifwhilebody, :(_job.pstate.logprior = _job.sstate.pstateprime.logprior))
  end
  if in(:accept, job.outopts[:diagnostics])
    push!(ifwhilebody, :(_job.sstate.update = true))
  end

  push!(
    whilebody,
    Expr(:if, :(_job.sstate.sprime && (rand() < _job.sstate.nprime/_job.sstate.n)), Expr(:block, ifwhilebody...))
  )

  push!(whilebody, :(_job.sstate.j += 1))
  push!(whilebody, :(_job.sstate.n += _job.sstate.nprime))
  push!(
    whilebody,
    :(
      _job.sstate.s =
      _job.sstate.sprime &&
      !uturn(
        _job.sstate.pstateplus.value,
        _job.sstate.pstateminus.value,
        _job.sstate.momentumplus,
        _job.sstate.momentumminus
      )
    )
  )

  push!(body, Expr(:while, :(_job.sstate.s && (_job.sstate.j < _job.sampler.maxndoublings)), Expr(:block, whilebody...)))

  dindex = findfirst(job.outopts[:diagnostics], :accept)
  if dindex != 0
    push!(body, :( _job.pstate.diagnosticvalues[$dindex] = _job.sstate.update))
  end

  dindex = findfirst(job.outopts[:diagnostics], :ndoublings)
  if dindex != 0
    push!(body, :( _job.pstate.diagnosticvalues[$dindex] = _job.sstate.j))
  end

  if !job.plain
    push!(body, :(produce()))
  end

  @gensym _iterate

  result = quote
    function $_iterate(_job::BasicMCJob)
      $(body...)
    end
  end

  result
end


text:
# Makie recipe

using Batsrus
using Makie

include("utility.jl")

# 1D example
filename = "1d__raw_2_t25.60000_n00000258.out"
data = readdata(filename, dir="test/data")

var = "p"
VarIndex_ = findindex(data, var)

fig = Figure(resolution = (800, 600))
ax = Axis(fig[1,1], title="1D shocktube", xlabel="x", ylabel="Pressure [nPa]")
line1 = lines!(ax, data.x[:], data.w[:,VarIndex_], color = :red,
   label="Pressure")
axislegend()

fig
#=
# 2D example
#filename = "z=0_raw_1_t25.60000_n00000258.out"
#data = readdata(filename, dir="test/data")

var = "p"
plotrange = [-Inf,Inf,-Inf,Inf]
plotinterval = 0.1

fig = Figure(resolution = (800, 600))

xi, yi, wi = getdata(data, var, plotrange, plotinterval)

fig[1,1] = Axis(fig, title="2D shocktube")
c = heatmap!(fig[1,1], xi, yi, wi)
cbar = Colorbar(fig[1,2], c, label="Pressure [nPa]", width=20)

fig
=#