text:
module JSONStat

# TODO: Deal with 'Collection' class
# TODO: status field
# TODO: validate(jstat)

using JSON3, Tables, PrettyTables

# Fields to include in dataset metadata and dimension+category data and REPL print
const metadatafields = (:source, :updated, :note, :extension, :href, :type, :link, :error) # excluding :version, :class , :role (added to category)
const dimfields = (:note, :href)
const categoryfields = (:child, :coordinates, :unit)

# Output struct for defining specific show method and accessors
struct Datatable <: Tables.AbstractColumns
    name::String
    data::NamedTuple
    dimensions::Dict
    metadata::Dict
end

function Base.show(io::IO, dt::Datatable)
    println("\n ", length(data(dt)), "x", length(data(dt)[1]), " JSONStat.Datatable")
    printstyled(" ", name(dt), "\n"; bold=true)
    if haskey(metadata(dt),:source)
        println(" ", metadata(dt)[:source])
    end
    if haskey(metadata(dt),:updated)
         println(" Updated: ", metadata(dt)[:updated])
    end
    pretty_table(dt; alignment=:l)
end

# Tables.jl implementation
Tables.istable(::Type{<:Datatable}) = true
Tables.columnaccess(::Type{<:Datatable}) = true
Tables.columns(dt::Datatable) = dt

# getter methods to avoid getproperty clash
name(dt::Datatable) = getfield(dt, :name)
data(dt::Datatable) = getfield(dt, :data)
dimensions(dt::Datatable) = getfield(dt, :dimensions)
metadata(dt::Datatable) = getfield(dt, :metadata)

# Methods for Tables.jl acces to data
Tables.columnnames(dt::Datatable) = propertynames(data(dt))
Tables.getcolumn(dt::Datatable, i::Int) = getfield(data(dt), i)
Tables.getcolumn(dt::Datatable, nm::Symbol) = getproperty(data(dt), nm)

"""
    JSONStat.read(js::Union{Vector{UInt8},String})

Constructs JSONStat.datatable compatible with Tables.jl

# Example
```julia
    HTTP.get("https://json-stat.org/samples/canada.json").body |> JSONStat.read
```
"""
function read(json::Union{Vector{UInt8},String})

    jstat = JSON3.read(json)

    if jstat.class == "dataset"
        readdataset(jstat)
    elseif jstat.class == "collection"
        error("Collection class not supported yet")
    else
        error("Class not supported yet")
    end
end

function readdataset(jstat::JSON3.Object)

    l = prod(jstat.size)
    dim = parsedimensions(jstat)
    (typeof(jstat.value) == Dict{String,Any}) ? data = dicttovect(jstat, l) : data = jstat.value
    mask = @. !ismissing(data) # 1 for non missing values

    headers = Vector{Symbol}(undef, length(jstat.size) + 1)
    values = Vector{Any}(undef, length(jstat.size) + 1)

    # Dimension columns
    f = 1
    for (i, dᵢ) in enumerate(dim)
        headers[i] = dᵢ.label
        values[i] = repeat(dᵢ.categories; inner=div(l, dᵢ.size * f), outer=f)[mask]
        f = f * dᵢ.size
    end

    # Value column
    headers[end] = :Value
    values[end] = data[mask]

    return Datatable(jstat.label,
        (; zip(headers, values)...), # NamedTuple of values for easy Tables.jl compliance
        dimensiondata(jstat),
        metadata(jstat))
end

function validate(js::JSON3.Object) end

"""parsedimensions(jstat) gathers basic information for building dimension columns and their category values.
Information is organized in a Vector of NamedTuples where each tuple has column label, a vector of categories labels/keys and size"""
function parsedimensions(jstat::JSON3.Object)

    dims = Vector()

    for (id, sz) in zip(jstat.id, jstat.size) # Constructs a Vector of Tuples containing (dimension, labels, size)

        # 1. Column names are the label of each dimension if it has one, if not the id is used
        nm = haskey(jstat.dimension[id], :label) ? Symbol(jstat.dimension[id][:label]) : jstat.dimension[id]

        # 2. Gather categories in the correct order
        categories = jstat.dimension[id].category

        if !haskey(categories, :index) # if there is no index, categories are the labels
            push!(dims, (; id=id, label=nm, categories=collect(values(jstat.dimension[id]["category"]["label"])), size=sz))
        else # if there is an index
            if isa(categories.index, JSON3.Object) # if categories is an Object dictionary we sort it. If it is an array it is alreary sorted
                order = sortperm(collect(values(categories.index)))
                orderedkeys = collect(keys(categories.index))[order] # sorts categories by index
            else
                orderedkeys = categories.index
            end

            if haskey(categories, :label) # if there is also labels
                push!(dims,
                    (; id=id, label=nm, categories=[categories.label[idx] for idx in orderedkeys], size=sz))
            else
                push!(dims, (; id=id, label=nm, categories=orderedkeys, size=sz))
            end
        end
    end #for

    return dims
end

"""Stores dataset metadata"""
function metadata(jstat::JSON3.Object)
    mt = Dict()
    for fᵢ in metadatafields
        haskey(jstat, fᵢ) && (mt[fᵢ] = jstat[fᵢ])
    end
    mt
end

"""Stores dimension and category metadata"""
function dimensiondata(jstat::JSON3.Object) # NOT WORKING
      dimdata = Dict()
      for dᵢ in keys(jstat.dimension)
        lbl =  haskey(jstat.dimension[dᵢ], :label) ? jstat.dimension[dᵢ][:label] : string(jstat.dimension[dᵢ])
        auxd = Dict()

        auxd[:categories] = haskey(jstat.dimension[dᵢ].category,:label) ? # Adding categories labels/ids
                            collect(values(jstat.dimension[dᵢ].category.label)) : string.(keys(jstat.dimension[dᵢ].category.index))

        for i in keys(jstat.dimension[dᵢ]) # Adding dimension root info
            if i in dimfields
                auxd[i] = JSON3.copy(dᵢ[i])
            end
        end
        if haskey(jstat, :role) # Adding dimension role if any
            for (k,v) in jstat.role
                if string(dᵢ) in v
                auxd[:role] = string(k)
                end
            end
        end
        for i in keys(jstat.dimension[dᵢ].category) # Adding category info TODO: Handle :unit better (has category label...)
            if i in categoryfields
                auxd[i] = JSON3.copy(jstat.dimension[dᵢ].category[i])
            end
        end
        dimdata[lbl] = auxd
      end
      dimdata
end

"""Helper function. Creates an Array from JSONStat when sparse values are delivered as a dictionary"""
function dicttovect(jstat::JSON3.Object, l::Int)
    v = Vector{Any}(missing, l)
    for (k, j) in jstat.value
        v[parse(Int, k)+1] = j #JSONStat uses zero-indexing
    end
    return v
end

"""Pretty printing dimension data"""
pretty(d::Dict) = JSON3.pretty(d)

end # Module

max_stars_repo_name: lsaenzt/JSONStat

max_issues_repo_name: lsaenzt/JSONStat

max_forks_repo_name: lsaenzt/JSONStat


text:
module LoggerManager

import Base.print, Base.println, Base.info, Base.warn, Base.error

export @printf # cannot replace Base.@printf

using Suppressor
using TimeManager
using CoreExtended

export LOGGER_OUT
export LOGGER_ERROR

export print
export printf
export println
export info
export warn
export error

LOGGER_OUT = "out.log"
LOGGER_ERROR = "error.log"

function __init__()
	open(LOGGER_OUT, "w+")
	open(LOGGER_ERROR, "w+")
end

""" TODO """
macro printf(xs...) open(f -> (:(Base.@printf(f, $xs...)); :(Base.@printf(:STDOUT, $xs...))),LOGGER_OUT, "a+") end

@suppress begin print(xs...) = open(f -> (print(f, xs...); print(STDOUT, xs...)), LOGGER_OUT, "a+") end
@suppress begin println(xs...) = open(f -> (println(f, programTimeStr(), " ", xs...); println(STDOUT, xs...)), LOGGER_OUT, "a+") end
@suppress begin info(xs...) = open(f -> (info(f, programTimeStr(), " ", xs...); info(STDOUT, xs...)), LOGGER_OUT, "a+") end
@suppress begin error(xs...) = open(f -> (error(f, programTimeStr()," ", xs...); error(STDERR, xs...)), LOGGER_ERROR, "a+") end
@suppress begin warn(xs...) = open(f -> (warn(f, programTimeStr()," ", xs...); warn(STDERR, xs...)), LOGGER_ERROR, "a+") end

""" TODO """
function logException(ex::Exception)
	t=programTimeStr()

	println(STDERR, "[ ERROR ] EXCEPTION! See '$LOGGER_ERROR' for more info.")
	open(f -> println(f, t, " [ ERROR ] EXCEPTION! See '$LOGGER_ERROR' for more info."), LOGGER_OUT, "a+")
	open(function(f)
			println(f, t, " --- [ BACKTRACE ] ---")
			Base.showerror(f, ex, catch_backtrace())
			println(f, "\n---------------------")
	end, LOGGER_ERROR, "a+")
end

CoreExtended.linkToException(logException)

#messageToString(mode::String, time::String, name::String, msg::String) = string(time, " ", mode != "" ? string("[",mode,"] "): "" , name != "" ? string(name, " ") : "" , ": ", msg, "\n")
#logMsg(mode::Symbol, name::String, msg::String) = messageToString(string(mode), programTimeStr(), name, msg)

end #LoggerManager

ext: jl

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
# kd tree for fast searching of a list of shapes
export KDTree
using Statistics: median

"""
A K-D tree of shapes in `K` dimensions.

This is a binary tree.  Each node in the tree divides space along
coordinate `ix` into shapes that overlap the region `≤ x` and `≥ x`.
The leaves of the tree are just a list of shapes `s`.

This implementation is a little different from a standard K-D tree,
because a given shape may be in both branches of the tree.  (Normally,
K-D trees are used for nearest-neighbor searches for a list of *points*,
not shapes of nonzero size.)
"""
mutable struct KDTree{K,S<:Shape{K}}
    s::Vector{S}
    ix::Int
    x::Float64
    left::KDTree{K,S}  # shapes ≤ x in coordinate ix
    right::KDTree{K,S} # shapes > x in coordinate ix
    KDTree{K,S}(s::AbstractVector{S}) where {K,S<:Shape{K}} = new(s, 0)
    function KDTree{K,S}(ix::Integer, x::Real, left::KDTree{K,S}, right::KDTree{K,S}) where {K,S<:Shape{K}}
        1 ≤ ix ≤ K || throw(BoundsError())
        new(S[], ix, x, left, right)
    end
end

Base.ndims(::KDTree{K}) where {K} = K

"""
    KDTree(s::AbstractVector{<:Shape{K}})

Construct a K-D tree (`KDTree`) representation of a list of
`shapes` in order to enable rapid searching of an shape list.

When searching the tree, shapes that appear earlier in `s`
take precedence over shapes that appear later.
"""
function KDTree(s::AbstractVector{S}) where {K,S<:Shape{K}}
    (length(s) ≤ 4 || K == 0) && return KDTree{K,S}(s)

    # figure out the best dimension ix to divide over,
    # the dividing plane x, and the number (nl,nr) of
    # shapes that fall into the left and right subtrees
    b = map(bounds, s) # cartesian bounding boxes of all shapes
    ix = 0
    x = 0.0
    nl = nr = typemax(Int)
    for i = 1:K
        mx = median(map(b -> 0.5*(b[1][i] + b[2][i]), b))
        mnl = count(b -> b[1][i] ≤ mx, b) # lower bound ≤ mx
        mnr = count(b -> b[2][i] > mx, b) # upper bound > mx
        if max(mnl,mnr) < max(nl,nr)
            ix = i
            x = mx
            nl = mnl
            nr = mnr
        end
    end

    # don't bother subdividing if it doesn't reduce the # of shapes much
    4*max(nl,nr) > 3*length(s) && return KDTree{K,S}(s)

    # create the arrays of shapes in each subtree
    sl = Vector{S}(undef, nl)
    sr = Vector{S}(undef, nr)
    il = ir = 0
    for k in eachindex(s)
        if b[k][1][ix] ≤ x
            sl[il += 1] = s[k]
        end
        if b[k][2][ix] > x
            sr[ir += 1] = s[k]
        end
    end

    return KDTree{K,S}(ix, x, KDTree(sl), KDTree(sr))
end

depth(kd::KDTree) = kd.ix == 0 ? 0 : max(depth(kd.left), depth(kd.right)) + 1

Base.show(io::IO, kd::KDTree{K,S}) where {K,S} = print(io, "KDTree{$K,$S} of depth ", depth(kd))

function _show(io::IO, kd::KDTree, indent)
    indentstr = " "^indent
    if kd.ix == 0
        println(io, indentstr, length(kd.s), " shapes")
    else
        println(io, indentstr, "if x[", kd.ix, "] ≤ ", kd.x, ':')
        _show(io, kd.left, indent + 2)
        println(io, indentstr, "else:")
        _show(io, kd.right, indent + 2)
    end
end

function Base.show(io::IO, ::MIME"text/plain", kd::KDTree)
    println(io, kd, ':')
    _show(io, kd, 0)
end

function Base.findfirst(p::SVector{N}, s::Vector{S}) where {N,S<:Shape{N}}
    for i in eachindex(s)
        b = bounds(s[i])
        if all(b[1] .< p .< b[2]) && p ∈ s[i]  # check if p is within bounding box is faster
            return s[i]
        end
    end
    return nothing
end

function Base.findfirst(p::SVector{N}, kd::KDTree{N}) where {N}
    if isempty(kd.s)
        if p[kd.ix] ≤ kd.x
            return findfirst(p, kd.left)
        else
            return findfirst(p, kd.right)
        end
    else
        return findfirst(p, kd.s)
    end
end

"""
    findfirst(p::AbstractVector, kd::KDTree)

Return the first shape in `kd` that contains the point `p`; return `nothing` otherwise.
"""
Base.findfirst(p::AbstractVector{<:Real}, kd::KDTree{N}) where {N} = findfirst(SVector{N}(p), kd)
Base.findfirst(p::AbstractVector{<:Real}, s::Vector{<:Shape{N}}) where {N} = findfirst(SVector{N}(p), s)


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module BinaryQF
using Nemo, NumberTheory

export ModuleBinaryQF
export L002476, L008784, L031363, L034017, L035251, L038872, L038873
export L042965, L057126, L057127, L068228, L084916, L089270, L141158, L242660
export L243655, L244779, L244780, L244819, L243168, L244291, L007522
# not yet implemented: L033200

"""

A binary quadratic form over Z is a quadratic homogeneous polynomial in two variables with integer coefficients, ``q(x, y) = ax^2 + bxy + cy^2``.

A quadratic form ``q(x, y)`` represents an integer ``n`` if there exist integers ``x`` and ``y`` with ``q(x, y) = n``. We say that ``q`` primitively represents ``n`` if there exist relatively prime integers ``x`` and ``y`` such that ``q(x, y) = n``.

Ported from [BinaryQuadraticForms](http://oeis.org/wiki/User:Peter_Luschny/BinaryQuadraticForms) where you can find much more information on this subject.

* L002476, L008784, L031363, L034017, L035251, L038872, L038873, L042965, L057126, L057127, L068228, L084916, L089270, L141158, L242660, L243655, L244779, L244780, L244819, L243168, L244291, L007522
"""
const ModuleBinaryQF = ""

# export binaryQF
# export L002144, L056874, L047486, L244713, L033203, L035247, L045373
# export L001481, L002313, L002479, L003136, L007645, L028951, L033207

# Note that we list the *positive* numbers of a form  whereas in the OEIS often
# the '0' is included. The reason for our choice is explained in Andrew Sutherland's
# "Introduction to Arithmetic Geometry", MIT Open Course Ware, Fall 2013.

"""

Return integers that are represented by the binary quadratic form ``a x^2 + b xy + c y^2`` over Z. Parameter 'subset' is in {"positive", "primitively", "prime"}. Use it only as an intern function.
"""
function binaryQF(
    a::Int,
    b::Int,
    c::Int,
    bound = 100::Int,
    subset = "positive",
    verbose = false
)

    α, β, γ = a, b, c

    discriminant() = fmpz(β^2 - 4 * α * γ)
    isreduced() = (-α < β ≤ α < γ) || (fmpz(0) ≤ β ≤ α == γ)

    function roots(a::Int, b::Int, c::Int, n::Int, y::Int)
        #################################
        throw(ErrorException("not yet implemented"))
        ################################
        #    x = var('x')
        #    eq = a * x * x + b * x * y + c * y * y
        #    (eq - n).roots(multiplicities = false, ring = fmpz)
    end

    function sqr_disc(M, primitively = false)

        d = discriminant()
        d == 0 && throw(ValueError("discriminant must not be zero"))

        a, b, c = α, β, γ
        (a == 0 && c == 0) && return [b * n for n ∈ 1:div(M, abs(b))]

        D = isqrt(d)
    # a must be ≠ 0
        if a == 0
            a, c = c, 0
        end
        k = 2 * D
        m = 4 * a * D
        u = b + D
        v = b - D
        S = fmpz[]

    # Solvability in Z.
        for n ∈ 1:M
            h = fmpz(4 * a * n)  # a <> 0 and n <> 0
            for t ∈ Divisors(h) # returns fmpz
                g = fmpz(div(h, t))
            # if divides(k, g - t) && divides(m, g * u - t * v)
                if rem(g - t, k) == 0 && rem(g * u - t * v, m) == 0

                    if primitively
                        y = div(g - t, k)
                        R = roots(a, b, c, n, y)
                        if isPrimeTo(R[1], y)
                            push!(S, n)
                            break
                        end
                    else
                        push!(S, n)
                        break
                    end
                end
            end
        end
        sort([s for s ∈ Set(S)])
    end

    function imag_prime(M)
        #################################
        throw(ErrorException("not yet implemented"))
        ################################
        #    solve = pari('qfbsolve')
        #    Q = pari('Qfb')(α, β, γ)
        #    p = 1
        #    R = []
        #    while true
        #        p = next_prime(p)
        #        p > M && break
        #        solve(Q, p) && push!(R,p)
        #    end
        #    R
    end

    function imag_primitively(M)

        a, b, c = α, β, γ
        d = c - div(b * b, 4 * a)
        A = []

        for y ∈ 0:isqrt(div(M, d))
            r = y * b / (2 * a)
            s = sqrt((M - d * y * y) / a)
            for x ∈ Int(round(ceil(-s - r))):Int(round(floor(s - r)))
                isPrimeTo(x, y) && push!(A, a*x^2 + b*x*y + c*y^2)
            end
        end
        sort([s for s ∈ Set(A)])
    end

    function imag_positive(M)
        #################################
        throw(ErrorException("not yet implemented"))
        ################################
        #   L = [2*fmpz(α), fmpz(β), fmpz(β), 2*fmpz(γ)]
        #   G = Matrix(fmpz, 2, 2, L)
        #   A = pari('qfrep')(G, M, 1)
        #   [k+1 for k ∈ 0:M-1 if A[k] > 0]
    end

    function primitive_reps(a, h, b, M, S)

        if a ≤ M
            push!(S, a)
            if b ≤ M
                push!(S, b)
                if a ≤ (M - b) && h ≤ (M - a - b)
                    (a ≤ (M - a - h)) && primitive_reps(a, h + 2 * a, a + b + h, M, S)
                    (b ≤ (M - b - h)) && primitive_reps(a + b + h, h + 2 * b, b, M, S)
                end
            end
        end
    end

    function positive_primitives(bound, primitively)

        a, b, c = α, β, γ
        S = fmpz[]

        while true
            new_val = a + b + c
            if new_val > 0
                primitive_reps(a, b + 2 * a, new_val, bound, S)
                b += 2 * c
                a = new_val
            else
                if new_val < 0
                    b += 2 * a
                    c = new_val
                end
            end
            if a == α && b == β && c == γ
                break
            end
        end

        if !primitively
            X = fmpz[]
            for p ∈ S
                q = t = 1
                while q ≤ bound
                    push!(X, q)
                    q = t * t * p
                    t += 1
                end
            end
            S = X
        end

        sort([s for s ∈ Set(S)])
    end

    function reduce_real()

        d = discriminant()
        isSquare(d) && throw(ValueError("form must not have square discriminant"))
        droot = isqrt(d)
        a, b, c = α, β, γ

        while a ≤ 0 || c ≥ 0 || b ≤ abs(a + c)

            cAbs = c
            if cAbs < 0
                cAbs *= -1
            end

            # cAbs = 0 will not happen for a non square form
            delta = div(b + droot, 2 * cAbs)
            if c < 0
                delta *= -1
            end
            aa = c
            bb = 2 * c * delta - b
            cc = c * delta * delta - b * delta + a
            a, b, c = aa, bb, cc
        end

        return [a, b, c]
    end

    function reduce_imag()

        a, b, c = α, β, γ
        if a < 0
            a, b, c = -a, -b, -c
        end
        d = discriminant()

        while true
            A = (a == c && b < 0) || (c < a)
            L = (-a == b && a < c) || (a < abs(b))

            !(A || L) && break

            if A
                a, b, c = c, -b, a
            end

            if L
                b -= 2 * a * div(b, 2 * a)
                if abs(b) > a
                    b -= 2 * a
                end
                c = div(b * b - d, 4 * a)
            end
        end

        return [a, b, c]
    end

    # Return the unique reduced form equivalent to BQF(a,b,c)
    function reduced_form()

        isreduced() && return [α, β, γ]

        if discriminant() ≥ 0
            return reduce_real()
        else
            return reduce_imag()
        end
    end

    ### --- eval part of function starts here --- ###

    prime = subset == "prime"
    primitively = subset == "primitively"

    d = discriminant()
    d == 0 && throw(ValueError("discriminant must not be 0"))

    a, b, c = α, β, γ
    if verbose
        println(
            "Original form [",
            a, ", ",
            b, ", ",
            c, "] with discriminant ",
            d
        )
    end

    if isSquare(d)

        verbose && println("Square discriminant!")
        if prime # for efficiency
            primitively = false # for efficiency
        end # for efficiency
        pp = sqr_disc(bound, primitively)
        if prime
            pp = [m for m ∈ pp if isPrime(m)]
        end

    else

        α, β, γ = reduced_form()
        verbose && println("Reduced form  [", α, ", ", β, ", ", γ, "]")

        if d < 0
            if prime
                pp = imag_prime(bound)
            else
                if primitively
                    pp = imag_primitively(bound)
                else
                    pp = imag_positive(bound)
                end
            end
        ### real case, indefinite form ###
        else # d > 0 and not square
            if prime
                primitively = true
            end # for efficiency
            pp = positive_primitives(bound, primitively)
            if prime
                pp = [m for m ∈ pp if isPrime(m)]
            end
        end
    end

    if verbose
        msg0 = prime ? " primes " : " positive integers "
        msg1 = primitively ? "primitively " : ""
        msg2 = "represented up to "
        println("There are ", length(pp), msg0, msg1, msg2, bound)
    end
    pp
end # binary_QF

"""

Return positive numbers of the form ``n = x^2-3y^2`` of discriminant 12.
"""
L084916(bound::Int) = binaryQF(1, 0, -3, bound)

"""

Return positive numbers that are primitively represented by the indefinite quadratic form ``x^2 - 3y^2`` of discriminant 12.
"""
L243655(bound::Int) = binaryQF(1, 0, -3, bound, "primitively")

"""

Return primes congruent to 1 (mod 12).
"""
L068228(bound::Int) = binaryQF(1, 0, -3, bound, "prime")

"""

Return positive numbers of the form ``x^2+xy-2y^2``.
"""
L242660(bound::Int) = binaryQF(1, 1, -2, bound)

"""

Return positive numbers primitively represented by the binary quadratic form (1, 1, -2).
"""
L244713(bound::Int) = binaryQF(1, 1, -2, bound, "primitively")

"""

Return primes of the form ``6m + 1``.
"""
L002476(bound::Int) = binaryQF(1, 1, -2, bound, "prime")

"""

Return positive numbers of the form ``x^2 - 2y^2`` with integers ``x, y`` (discriminant is 8).
"""
L035251(bound::Int) = binaryQF(1, 0, -2, bound)

"""

Return mumbers n such that 2 is a square mod n.
"""
L057126(bound::Int) = binaryQF(1, 0, -2, bound, "primitively")

"""

Return primes p such that 2 is a square mod p; or, primes congruent to ``{1, 2, 7}`` mod ``8``.
"""
L038873(bound::Int) = binaryQF(1, 0, -2, bound, "prime")

#"""
#Return primes congruent to ±1 mod 8 apart from the first term.
#"""
#L001132(bound::Int) = binaryQF(1, 0, -2, bound, "prime")

"""

Return positive numbers of the form ``x^2+xy-y^2``; or, of the form ``5x^2-y^2``.
"""
L031363(bound::Int) = binaryQF(1, 1, -1, bound) # "positive"

"""

Return positive numbers represented by the integer binary quadratic form ``x^2+xy-y^2`` with ``x`` and ``y`` relatively prime.
"""
L089270(bound::Int) = binaryQF(1, 1, -1, bound, "primitively")

"""

Return primes represented by the integer binary quadratic form ``x^2+xy-y^2``.
"""
L141158(bound::Int) = binaryQF(1, 1, -1, bound, "prime")

"""

Return primes congruent to ``{0, 1, 4}`` mod 5 (cf. also `[A141158]`).
"""
L038872(bound::Int) = binaryQF(1, 1, -1, bound, "prime")

"""

Return positive integers not congruent to 2 mod 4; regular numbers modulo 4.
"""
L042965(bound::Int) = binaryQF(1, 0, -1, bound) # "positive"

"""

Return numbers that are congruent to ``{0, 1, 3, 5, 7}`` mod 8. Positive integers represented by the binary quadratic form ``x^2-y^2`` with ``x`` and ``y`` relatively prime.
"""
L047486(bound::Int) = binaryQF(1, 0, -1, bound, "primitively")

#"""
#Return odd primes; primes represented by the binary quadratic form ``x^2-y^2``.
#"""
#L065091(bound::Int) = binaryQF(1, 0, -1, bound, "prime")

"""

Return positive integers of the form ``x^2 + xy + y^2`` (Loeschian numbers).
"""
L003136(bound::Int) = binaryQF(1, 1, 1, bound) # "positive"

"""

Return positive integers that are primitively represented by ``x^2 + xy + y^2``.
"""
L034017(bound::Int) = binaryQF(1, 1, 1, bound, "primitively")

"""

Return primes that are represented by ``x^2 + xy + y^2`` (generalized cuban primes).
"""
L007645(bound::Int) = binaryQF(1, 1, 1, bound, "prime")

"""

Return positive integers that are the sum of 2 squares.
"""
L001481(bound::Int) = binaryQF(1, 0, 1, bound) # "positive"

"""

Return numbers ``n`` that are primitively represented by ``x^2 + y^2``. Also numbers ``n`` such that ``√(-1)`` mod ``n`` exists.
"""
L008784(bound::Int) = binaryQF(1, 0, 1, bound, "primitively")

"""

Return primes of form ``x^2 + y^2``; or primes congruent to 1 or 2 modulo 4.
"""
L002313(bound::Int) = binaryQF(1, 0, 1, bound, "prime")

"""

Return Pythagorean primes: primes of form ``4n + 1``.
"""
L002144(bound::Int) = binaryQF(1, 0, 1, bound, "prime")

"""

Return positive integers of the form ``x^2+xy+2y^2`` with ``x`` and ``y`` integers. See also A035248.
"""
L028951(bound::Int) = binaryQF(1, 1, 2, bound) # "positive"

"""

Return positive numbers primitively represented by the binary quadratic form (1, 1, 2).
"""
L244779(bound::Int) = binaryQF(1, 1, 2, bound, "primitively")

"""

Return primes represented by the binary quadratic form (1, 1, 2). Primes congruent to ``{0, 1, 2, 4}`` mod 7.
"""
L045373(bound::Int) = binaryQF(1, 1, 2, bound, "prime")

"""

Return primes of form ``x^2+7*y^2``.
"""
L033207(bound::Int) = binaryQF(1, 1, 2, bound, "prime")

"""

Return integers of form ``x^2 + 2y^2``.
"""
L002479(bound::Int) = binaryQF(1, 0, 2, bound)

"""

Return positive integers primitively represented by ``x^2 + 2y^2``.
"""
L057127(bound::Int) = binaryQF(1, 0, 2, bound, "primitively")

"""

Return primes of form ``x^2+2*y^2``. Primes congruent to ``{1, 2, 3}`` mod 8.
"""
L033203(bound::Int) = binaryQF(1, 0, 2, bound, "prime")

"""

Return odd primes of form ``x^2+2y^2``. Primes congruent to ``{1, 3}`` mod 8.
"""
L033200(bound::Int) = binaryQF(1, 0, 2, bound, "prime")

"""

Return positive numbers represented by the binary quadratic form (1, 1, 3). (See also A028954.)
"""
L035247(bound::Int) = binaryQF(1, 1, 3, bound) # "positive"

"""

Return positive numbers primitively represented by the binary quadratic form (1, 1, 3).
"""
L244780(bound::Int) = binaryQF(1, 1, 3, bound, "primitively")

"""

Return primes of form ``x^2+xy+3y^2``, discriminant -11.
"""
L056874(bound::Int) = binaryQF(1, 1, 3, bound, "prime")

"""

Return positive numbers primitively represented by the binary quadratic form (1, 0, 3).
"""
L244819(bound::Int) = binaryQF(1, 0, 3, bound, "primitively")

"""

Return positive integers of the form ``x^2+6xy-3y^2``.
"""
L243168(bound::Int) = binaryQF(1, 6, -3, bound)

"""

Return positive numbers primitively represented by the binary quadratic form (1, 6, -3).
"""
L244291(bound::Int) = binaryQF(1, 6, -3, bound, "primitively")

"""

Return primes of the form ``8n+7``, that is, primes congruent to -1 mod 8.
"""
L007522(bound::Int) = binaryQF(-1, 4, 4, bound, "prime")

#START-TEST-########################################################

using Test

function test()

    # Note that we list the *positive* numbers of a form  whereas in the OEIS
    # often the '0' is included.

    Data = Dict{Int,Array{fmpz}}(
        002476 => [7, 13, 19, 31, 37, 43, 61, 67, 73, 79, 97, 103],
        008784 => [1, 2, 5, 10, 13, 17, 25, 26, 29, 34, 37, 41],
        031363 => [1, 4, 5, 9, 11, 16, 19, 20, 25, 29, 31, 36],
        034017 => [1, 3, 7, 13, 19, 21, 31, 37, 39, 43, 49, 57],
        035251 => [1, 2, 4, 7, 8, 9, 14, 16, 17, 18, 23, 25],
        038872 => [5, 11, 19, 29, 31, 41, 59, 61, 71, 79, 89, 101],
        038873 => [2, 7, 17, 23, 31, 41, 47, 71, 73, 79, 89, 97],
        042965 => [1, 3, 4, 5, 7, 8, 9, 11, 12, 13, 15, 16],
        057126 => [1, 2, 7, 14, 17, 23, 31, 34, 41, 46, 47, 49],
        057127 => [1, 2, 3, 6, 9, 11, 17, 18, 19, 22, 27, 33],
        068228 => [13, 37, 61, 73, 97, 109, 157, 181, 193, 229, 241, 277],
        084916 => [1, 4, 6, 9, 13, 16, 22, 24, 25, 33, 36, 37],
        089270 => [1, 5, 11, 19, 29, 31, 41, 55, 59, 61, 71, 79],
        141158 => [5, 11, 19, 29, 31, 41, 59, 61, 71, 79, 89, 101],
        242660 => [1, 4, 7, 9, 10, 13, 16, 18, 19, 22, 25, 27],
        243655 => [1, 6, 13, 22, 33, 37, 46, 61, 69, 73, 78, 94],
        244779 => [1, 2, 4, 7, 8, 11, 14, 16, 22, 23, 28, 29],
        244780 => [1, 3, 5, 9, 11, 15, 23, 25, 27, 31, 33, 37],
        244819 => [1, 3, 4, 7, 12, 13, 19, 21, 28, 31, 37, 39],
        243168 => [1, 4, 9, 13, 16, 24, 25, 33, 36, 37, 49, 52],
        244291 => [1, 4, 13, 24, 33, 37, 52, 61, 69, 73, 88, 97],
        007522 => [7, 23, 31, 47, 71, 79, 103, 127, 151, 167, 191, 199]
    )

    L = [
        L002476,
        L008784,
        L031363,
        L034017,
        L035251,
        L038872,
        L038873,
        L042965,
        L057126,
        L057127,
        L068228,
        L084916,
        L089270,
        L141158,
        L242660,
        L243655,
        L244779,
        L244780,
        L244819,
        L243168,
        L244291,
        L007522
    ]

    SeqNum(seq) = parse(Int, string(seq)[end-5:end])

    @testset "BinaryQF" begin

        for seq ∈ L
            S = seq(300)
            anum = SeqNum(seq)
            data = Data[anum]
            @test all(S[1:11] .== data[1:11])
        end
    end
end

function demo()
    # binaryQF(1, 1, 1, bound)
    # binaryQF(1, 1, 1, bound, "primitively")
    # binaryQF(1, 1, 1, bound, "prime")
    println(L034017(50))
end

"""

L034017(10000)
    0.009646 seconds (133.97 k allocations: 2.701 MiB)
"""
function perf()
    GC.gc()
    @time L034017(10000)
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
# Object is mutable. Can we ensure the content never changes?

# modify the getproperty for :contents field by returning a copy.
function Base.getproperty(fc::FileContent, s::Symbol)
    if s in (:path,)
        return getfield(fc, s)
    elseif s === :contents
        getfield(fc, :loaded) || load_contents!(fc)
        return copy(getfield(fc, :contents))   # return copy
    else
        error("Unsupported property: $s")
    end
end

#=
julia> fc.contents
359-element Array{UInt8,1}:
 0x23
 0x23
 0x0a

julia> fc.contents[1] = 0x00
0x00

julia> fc.contents
359-element Array{UInt8,1}:
 0x23
 0x23
 0x0a

 =#


text:
##############################################################################
##
## IsNested
##
##############################################################################

function isnested(fe::FixedEffect, p::CategoricalVector)
    entries_in_p = Dict{eltype(fe.refs), eltype(p.refs)}()
    sizehint!(entries_in_p, fe.n)
    for (feref, pref) in zip(fe.refs, p.refs)
        x = get(entries_in_p, feref, 0)
        if x == 0
            # it's a new level, create entry
            entries_in_p[feref] = pref
        elseif x != pref
            # not nested: for the same level in a, two different levels in b
            return false
        end
    end
    return true
end

##############################################################################
##
## Remove Singletons
##
##############################################################################

function remove_singletons!(esample, fe::FixedEffect)
    cache = zeros(Int, fe.n)
    for i in 1:length(esample)
        if esample[i]
            cache[fe.refs[i]] += 1
        end
    end
    for i in 1:length(esample)
        if esample[i] && cache[fe.refs[i]] <= 1
            esample[i] = false
        end
    end
end

##############################################################################
##
## Subet and Make sure Interaction if Vector{Float64} (instead of missing)
##
##############################################################################

# index and convert interaction Vector{Float64,, Missing} to Vector{Missing}
function _subset(fe::FixedEffect, esample)
    FixedEffect{typeof(fe.refs), Vector{Float64}}(fe.refs[esample], convert(Vector{Float64}, view(fe.interaction, esample)), fe.n)
end


text:
# Full model objects
parameter_defaults_tests = @with_kw (ρ = 0.0215,
                                σ = 3.1725,
                                N = 10,
                                θ = 5.0018,
                                γ = 1.00,
                                κ = 0.0732,
                                ζ = 1.0,
                                η = 0.,
                                Theta = 1,
                                χ = 1/5.9577,
                                υ = 0.0484 ,
                                μ = -0.0189,
                                δ = 0.02,
                                d_0 = 3.0426,
                                d_T = 2.83834)

settings_defaults_tests = @with_kw (z_ex = unique([range(0., 0.1, length = 150); range(0.1, 1., length = 100); range(1., 5, length = 50)]),
                                T = 75.0,
                                tstops = unique([range(0.0, 10.0, length=41); range(10.0, 20.0, length=21); range(20.0, T, length=56)]),
                                interp = Interpolations.LinearInterpolation,
                                stationary_x0 = default_stationary_x0,
                                pre_shock_times = [-1, -5, -10, -15, -20],
                                fixedpoint_ftol = 1e-8,
                                fixedpoint_show_trace = false,
                                fixedpoint_m = 5,
                                fixedpoint_beta = 1.0,
                                fixedpoint_iterations = 500,
                                fixedpoint_x0 = default_fixedpoint_x0)

parameter_defaults = @with_kw (θ = 4.988976587938262,
                                κ = 0.104196324793307,
                                χ = 1/7.883537319864373,
                                μ = -0.031064624217571,
                                υ = 0.048301140601665,
                                σ = 3.166924135838110,
                                ζ = 1,
                                δ = 0.020,
                                ρ = 0.020338044668517,
                                N = 10,
                                γ = 1.00,
                                η = 0.,
                                Theta = 1,
                                d_0 = 3.022492825462601,
                                d_T = 2.820243542916341)

settings_defaults = @with_kw (z_ex = unique([range(0., 0.1, length = 90); range(0.1, 1., length = 120); range(1., 5, length = 60)]),
                                T = 75.0,
                                tstops = unique([range(0.0, 10.0, length=41); range(10.0, 20.0, length=21); range(20.0, T, length=56)]),
                                interp = Interpolations.LinearInterpolation,
                                stationary_x0 = default_stationary_x0,
                                pre_shock_times = [-1, -5, -10, -15, -20],
                                fixedpoint_ftol = 1e-8,
                                fixedpoint_show_trace = false,
                                fixedpoint_m = 5,
                                fixedpoint_beta = 1.0,
                                fixedpoint_iterations = 500,
                                fixedpoint_x0 = default_fixedpoint_x0)

function default_fixedpoint_x0(parameters, settings)
    if length(settings.tstops) == 116 # HARDCODE THE KNOWN NUMBER OF TSTOPS FOR THE KNOWN SOLUTION
        return [-1.0, -0.988616, -0.977299, -0.96605, -0.954868, -0.943753, -0.932706, -0.921726, -0.910813, -0.899968, -0.88919, -0.878479, -0.867836, -0.85726, -0.846752, -0.836311, -0.825937, -0.81563, -0.805391, -0.79522, -0.785115, -0.775079, -0.765109, -0.755207, -0.745372, -0.735605, -0.725905, -0.716272, -0.706707, -0.69721, -0.687779, -0.678416, -0.669121, -0.659893, -0.650732, -0.641639, -0.632613, -0.623655, -0.614764, -0.605941, -0.588427, -0.571182, -0.554207, -0.5375, -0.521062, -0.504894, -0.488996, -0.473363, -0.458001, -0.442913, -0.42809, -0.413535, -0.399253, -0.385243, -0.371497, -0.358024, -0.344837, -0.331924, -0.319282, -0.306909, -0.2827, -0.2596, -0.237648, -0.216737, -0.196736, -0.177814, -0.160075, -0.143496, -0.12804, -0.113637, -0.100236, -0.0877941, -0.0763549, -0.0659462, -0.0567017, -0.0484165, -0.0408617, -0.0342745, -0.0282847, -0.0234237, -0.0188387, -0.0150505, -0.0118514, -0.00928253, -0.00711861, -0.00539311, -0.00401417, -0.00291141, -0.00208732, -0.00149604, -0.00101103, -0.000673539, -0.000453234, -0.000278941, -0.000168804, -0.0001053, -5.716e-5, -3.03866e-5, -1.73678e-5, -8.03325e-6, -3.65108e-6, -1.99316e-6, -8.91957e-7, -4.8973e-7, -3.65493e-7, -2.81903e-7, -2.38864e-7, -2.06177e-7, -1.73548e-7, -1.40636e-7, -1.07106e-7, -7.23857e-8, -3.64788e-8, 0.0]
    elseif length(settings.tstops) == 88
        return [-1.0, -0.96668, -0.934115, -0.902297, -0.871215, -0.840863, -0.811232, -0.782311, -0.754093, -0.726569, -0.69973, -0.673566, -0.648068, -0.623228, -0.599038, -0.575489, -0.552577, -0.530443, -0.508997, -0.488145, -0.467877, -0.448185, -0.42906, -0.410492, -0.392473, -0.374994, -0.358046, -0.341621, -0.325708, -0.310301, -0.295389, -0.280964, -0.267019, -0.253544, -0.240532, -0.227973, -0.21586, -0.204186, -0.180241, -0.15818, -0.137926, -0.119404, -0.102549, -0.0872978, -0.0735807, -0.0613147, -0.0503832, -0.0407036, -0.0323833, -0.0253242, -0.0195078, -0.0147896, -0.010991, -0.00799515, -0.00567981, -0.00399821, -0.00273242, -0.00183129, -0.00116784, -0.000741846, -0.000442362, -0.000279069, -0.000153617, -8.31759e-5, -4.87385e-5, -2.34855e-5, -1.14183e-5, -6.49374e-6, -3.02169e-6, -1.62968e-6, -1.15147e-6, -8.22077e-7, -6.73772e-7, -5.97719e-7, -5.32128e-7, -4.79283e-7, -4.33382e-7, -3.88373e-7, -3.45597e-7, -3.04862e-7, -2.6164e-7, -2.16769e-7, -1.69966e-7, -1.18744e-7, -6.07212e-8, 0.0]
    else
        return range(-1, 0, length = length(settings.tstops))[2:end-1] #guess linear
    end
end

default_stationary_x0(parameters, settings) = [parameters.ρ * 0.75; 2.0; 0.57]

model_cachename(parameters, settings) = join(hash((
        parameters = parameters,
        settings = merge(settings, (interp = typeof(settings.interp),
                                    stationary_x0 = typeof(settings.stationary_x0),
                                    fixedpoint_x0 = typeof(settings.fixedpoint_x0)))
)))

function load_parameters(filepath::String)
    tmp = load(filepath) |> collect |> first
    return (ρ = tmp.rho, σ = tmp.sigma, θ = tmp.theta, κ = tmp.kappa, χ = tmp.chi, μ = tmp.mu, υ = tmp.upsilon,
            ζ = tmp.zeta, δ = tmp.delta, N = tmp.N, γ = tmp.gamma, η = tmp.eta, Theta = tmp.Theta,
            d_0 = tmp.d_0, d_T = tmp.d_T)
end


text:
###
### Hash
###
###
### MurmurHash3 function of BioSequence.
###
### The hash function defined here cares about the starting position of a
### character sequence in the underlying data. That means, even if the starting
### positions of two sequences (`s1` and `s2`) are different in their `data`
### field, their hash values are identical if `s1 == s1` is true.  The
### implementation is based on the 128bit MurmurHash3 function, which was written
### by Austin Appleby, and the source code is distributed under the public domain:
### https://github.com/aappleby/smhasher/blob/master/src/MurmurHash3.cpp
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

const c1 = 0x87c37b91114253d5
const c2 = 0x4cf5ad432745937f

@inline function rotl64(x::UInt64, r)
    return (x << (r & 63)) | (x >>> (-r & 63))
end

@inline function fmix64(k::UInt64)
    k = k ⊻ k >> 33
    k *= 0xff51afd7ed558ccd
    k = k ⊻ k >> 33
    k *= 0xc4ceb9fe1a85ec53
    k = k ⊻ k >> 33
    return k
end

@inline function murmur1(h1, k1)
    k1 *= c1
    k1 = rotl64(k1, 31)
    k1 *= c2
    h1 = h1 ⊻ k1
    return (h1, k1)
end

@inline function murmur2(h1, h2, k2)
    k2 *= c2
    k2 = rotl64(k2, 33)
    k2 *= c1
    h2 = h1 ⊻ k2
    return (h2, k2)
end

@inline function murmur(h1, h2, k1, k2)
    h1, k1 = murmur1(h1, k1)
    h1 = rotl64(h1, 27)
    h1 += h2
    h1 = h1 * 5 + 0x52dce729

    h2, k2 = murmur2(h1, h2, k2)
    h2 = rotl64(h2, 31)
    h2 += h1
    h2 = h2 * 5 + 0x38495ab5

    return (h1, h2, k1, k2)
end

function finalize(h1, h2, len)
    h1 = h1 ⊻ len
    h2 = h2 ⊻ len
    h1 += h2
    h2 += h1
    h1 = fmix64(h1)
    h2 = fmix64(h2)
    h1 += h2
    h2 += h1

    return h1
end

function tail(data, next, stop, h1, h2)
    k1 = k2 = zero(UInt64)

    # Use this to mask any noncoding bits in last chunk
    mask = bitmask(offset(stop - bits_per_symbol(stop)) + bits_per_symbol(stop))
    @inbounds if next < stop
        k1 = data[index(next)]
        next += 64
    end
    @inbounds if next < stop
        k2 = data[index(next)] & mask
    else
        k1 &= mask
    end

    h1, k1 = murmur1(h1, k1)
    h2, k2 = murmur2(h1, h2, k2)
    return (h1, h2)
end

# ref: MurmurHash3_x64_128
function Base.hash(seq::LongSequence, seed::UInt64)
    # Mix sequence length so that dna"A" and dna"AA"
    # return the different hash values.
    h1::UInt64 = h2::UInt64 = hash(length(seq), seed)
    next = bitindex(seq, 1)
    last = bitindex(seq, lastindex(seq) + 1)
    data = encoded_data(seq)

    @inbounds while last - next ≥ 128
        j = index(next)
        k1 = data[j]
        k2 = data[j+1]
        h1, h2, k1, k2 = murmur(h1, h2, k1, k2)
        next += 128
    end

    h1, h2 = tail(data, next, last, h1, h2)
    h1 = finalize(h1, h2, length(seq))

    return h1
end


text:
module CauerLowPassFilterModel

using TinyModia

using DifferentialEquations
using Test
@usingModiaPlot

include("../models/Electric.jl")

l1 = 1.304
l2 = 0.8586
c1 = 1.072*1u"F"
c2 = 1/(1.704992^2*l1)*1u"F"
c3 = 1.682*1u"F"
c4 = 1/(1.179945^2*l2)*1u"F"
c5 = 0.7262*1u"F"

CauerLowPassOPVWithoutNodes(i=1) = Model(
  # Cauer low pass filter with operational amplifiers

  C1 = Capacitor | Map(C=c1 + c2),
  C2 = Capacitor | Map(C=c2, v=Var(init=nothing)),
  C3 = Capacitor | Map(C=l1*1u"F"),
  C4 = Capacitor | Map(C=c4*i, v=Var(init=nothing)),
  C5 = Capacitor | Map(C=c2, v=Var(init=nothing)),
  R1 = Resistor | Map(R=1u"Ω"),
  R2 = Resistor | Map(R=1u"Ω"),
  R3 = Resistor | Map(R=1u"Ω"),
  Op1 = IdealOpAmp3Pin,
  G = Ground,
  R4 = Resistor | Map(R=-1u"Ω"),
  R5 = Resistor | Map(R=-1u"Ω"),
  Op2 = IdealOpAmp3Pin,
  Op3 = IdealOpAmp3Pin,
  G1 = Ground,
  R6 = Resistor | Map(R=1u"Ω"),
  R7 = Resistor | Map(R=1u"Ω"),
  C6 = Capacitor | Map(C=c2 + c3 + c4, v=Var(init=nothing)),
  R8 = Resistor | Map(R=-1u"Ω"),
  R9 = Resistor | Map(R=-1u"Ω"),
  R10 = Resistor | Map(R=1u"Ω"),
  Op4 = IdealOpAmp3Pin,
  Op5 = IdealOpAmp3Pin,
  C7 = Capacitor | Map(C=l2*1u"F"),
  C8 = Capacitor | Map(C=c4, i=Var(start=0u"A")),
  C9 = Capacitor | Map(C=c4 + c5),
  R11 = Resistor | Map(R=1u"Ω"),

  G2 = Ground,
  G3 = Ground,
  G4 = Ground,
  V = ConstantVoltage | Map(V=1.0u"V"),
  Ground1 = Ground,

connect = :[
  (Op1.in_p, G.p)
  (G1.p, Op2.in_p)
  (R1.n, Op1.in_n, C2.n, R2.n, C1.p, R3.p)
  (R3.n, C1.n, Op1.out, R4.p, C5.p)
  (R4.n, Op2.in_n, C3.p, R5.n)
  (C2.p, R5.p, Op3.out, C6.n, R9.p, C8.p)
  (C3.n, Op2.out, R2.p, R7.p)
  (R7.n, Op3.in_n, C5.n, R6.n, C6.p, C4.n)
  (C4.p, R8.p, R11.n, C9.n, Op5.out)
  (R9.n, R8.n, Op4.in_n, C7.p)
  (R6.p, C7.n, Op4.out, R10.p)
  (G2.p, Op3.in_p)
  (R11.p, C9.p, R10.n, Op5.in_n, C8.n)
  (Op4.in_p, G3.p)
  (Op5.in_p, G4.p)
  (V.p, Ground1.p)
  (V.n, R1.p)  ]
)

n = 10
println("Build array of $n Cauer low pass filters")
@time Filters = Model(
    filters = [CauerLowPassOPVWithoutNodes(0.1*i) for i in 1:n]
)

@test_skip begin
model = @instantiateModel(Filters, logDetails=false, logTiming=true, unitless=true)

println("Simulate")
@time simulate!(model, Tsit5(), stopTime = 60, requiredFinalStates =
    [-0.5000732186007855, -0.5002239998029879, 0.4996923410661849, -0.4996706636198064, -0.5000929741546876, -0.5001255383344897, -0.500147742207576, 0.49981113006573824, -0.4996196061432069, -0.5001392063603706, -0.5001815616147427, -0.5000459266587362, 0.49996685095395915, -0.49958392015463465, -0.5001886259176451, -0.5002389821958331, -0.4999173001411904, 0.5001611780149184, -0.49957016523210945, -0.5002393712119582, -0.5002944833191254, -0.49976183273648583, 0.5003940588092074, -0.49958591550449066, -0.5002887040151874, -0.500343478728503, -0.4995812991530469, 0.5006629055320556, -0.49963969356758026, -0.5003327561376948, -0.50037980690654, -0.49938004306850625, 0.5009615326279456, -0.4997407922705924, -0.5003662177501006, -0.500395375828166, -0.4991659724894459, 0.5012787730322307, -0.4998989295998676, -0.5003819572251011, -0.5003797586125351, -0.49895184523346914, 0.5015966895590688, -0.5001236602409048, -0.5003705619890064, -0.5003197458974357, -0.49875691633150787, 0.5018882791340892, -0.5004234377294563, -0.5003197904255381])
plot(model, ["filters_1.C9.v", "filters_2.C9.v", "filters_3.C9.v"])
end

end


text:
qadd(board::Array{Vector{Int}},
         qrow::Vector{Int}) = push!(copy(board), qrow)
qhits(a::Array{Int}, b::Array{Int}) =
           any(a .== b) || abs(a-b)[1] == abs(a-b)[2]
qhit(qrow::Vector{Int}, board::Array{Vector{Int}}) =
any(map(x->qhits(qrow, x), board))

function qsolve(m, n=m, board=Array(Vector{Int}, 0))
  if n == 0
    return board
  end
  for px = 1:m
    for py = 1:m
      if !qhit([px, py], board)
        s = qsolve(m, n-1, qadd (board, [px, py]))
        s != nothing && return s
      end
    end
  end
  return nothing
end

max_stars_repo_path: Module 3/Chapter03/queens.jl

max_issues_repo_path: Module 3/Chapter03/queens.jl

max_forks_repo_path: Module 3/Chapter03/queens.jl


text:
function model(du,u,p,t)
    ts = [50.0, 100.0, 200.0, 900.0]
    freq = 1.0 * 0.2
    amp = p["amp"]
    if t < ts[1]
        p["maxA"] = 0.0
    elseif t < ts[3]
        p["maxA"] = 2.3
    else
        p["maxA"] = 2.3*0.65
        p["gA"] = 0.005
    end
    if t < ts[2]
        p["maxF"] = 0.0
    elseif t < ts[3]
        p["maxF"] = 2.1
    else
        p["maxF"] = 2.1*0.65
        p["gF"] = 0.005
    end
    if t < ts[4]
        p["muA"] = 0.0
        p["muF"] = 0.0
    else
        Random.seed!(floor(Int, (t-ts[4])*freq)*p["y"]*floor(Int, p["x"]))
        p["muA"] =   sign(sin((t-ts[4])*pi*freq)) * amp * randn()*u[5]
        p["muF"] =  sign(sin((t-ts[4])*pi*freq)) * amp * randn()*u[6]
    end
    P5, P8, R5, R8, A, F, r2, c15, c25, c35, r9, c18, c28, c38, c48 = u
    x, gA, gF, maxA, maxF, s05, s08, K5A, K8F, n5A, n8F, K58, n58  = p["x"], p["gA"], p["gF"], p["maxA"], p["maxF"], p["s05"], p["s08"], p["K5A"], p["K8F"], p["n5A"], p["n8F"], p["K58"], p["n58"]
    K9F, K2A, n9F, n2A, K25, K98, n25, n98 = p["K9F"], p["K2A"], p["n9F"], p["n2A"], p["K25"], p["K98"], p["n25"], p["n98"]
    gc, K = p["gc"], p["K"]
    l05, l15, l25, l35, kP5 = p["l05"], p["l15"], p["l25"], p["l35"], p["kP5"]
    l08, l18, l28, l38, l48, kP8 = p["l08"], p["l18"], p["l28"], p["l38"], p["l48"], p["kP8"]
    kR5, kR8, a15, a25, a35, a18, a28, a38, a48, s2, s9 = p["kR5"], p["kR8"], p["a15"], p["a25"], p["a35"], p["a18"], p["a28"], p["a38"], p["a48"], p["s2"], p["s9"]
    kr2, kr9, b15, b25, b35, b18, b28, b38, b48 = p["kr2"], p["kr9"], p["b15"], p["b25"], p["b35"], p["b18"], p["b28"], p["b38"], p["b48"]
    s5, s8  = p["s5"], p["s8"]
    muA, muF = p["muA"], p["muF"]
    du[1] = dP5 = l05 * (R5 - 3 * c15 - 3 * c25 - c35) + l15 * 3 * c15 + l25 * 3 * c25 + l35 * c35 - kP5 * P5
    du[2] = dP8 = l08 * (R8 - 4 * c18 - 6 * c28 - 4 * c38 - c48) + l18 * 4 * c18 + l28 * 6 * c28 + 4 * l38 * c38 + l48 * c48 - kP8 * P8
    du[3] = dR5  = s5*((A/K5A)^n5A)/(1+(P8/K58)^n58+(A/K5A)^n5A) + s05 - kR5 * ((R5- 3 * c15 - 3 * c25 - 1 * c35 ) + a15 * 3 * c15 + a25 * 3 * c25 + a35 * 1 * c35)
    du[4] = dR8 = s8*((F/K8F)^n8F)/(1+(F/K8F)^n8F) + s08 - kR8 * ((R8 - 4 * c18 - 6 * c28 - 4 * c38 - c48) + a18 * 4 * c18 + a28 * 6 * c28 + a38 * 4 * c38 + a48 * c48)
    du[5] = dA = gA*(maxA * exp(-0.01 * x) - A) + muA
    du[6] = dF = gF*(maxF * exp(-0.01 * (99 - x)) - F) + muF
    du[7] = dr2 = s2 / (1 + (A/K2A)^n2A) - kr2 * (r2 - 1 * 3 * c15 - 2 * 3 * c25 - 3 * 1 * c35 ) - b15 * 1 * 3 * c15 - b25 * 2 * 3 * c25 - b35 * 3 * 1 * c35
    du[8] = dc15 = gc*(K*(R5-3*c15-3*c25-c35)*(r2-3*c15-3*2*c25-3*c35)-c15)
    du[9] = dc25 = gc * (K * c15 * (r2 - 3 * c15 - 3 * 2 * c25 - 3 * c35) - c25)
    du[10] = dc35 = gc * (K * c25 * (r2 - 3 * c15 - 3 * 2 * c25 - 3 * c35) - c35)
    du[11] = dr9 = s9 * (F/K9F)^n9F/(1+(F/K9F)^n9F) - kr9 * (r9 - 1 * 4 * c18 - 2 * 6 * c28 - 3 * 4 * c38 - 4 * 1 * c48) - b18 * 1 * 4 * c18  - b28 * 2 * 6 * c28 - b38 * 3 * 4 * c38 - b48 * 4 * 1 * c48
    du[12] = dc18 = gc * (K * (R8 - 4 * c18 - 6 * c28 - 4 * c38 - c48) * (r9 - 1 * 4 * c18 - 2 * 6 * c28 - 3 * 4 * c38 - 4 * 1 * c48) - c18)
    du[13] = dc28 = gc * (K * c18 * (r9 - 1 * 4 * c18 - 2 * 6 * c28 - 3 * 4 * c38 - 4 * 1 * c48) - c28)
    du[14] = dc38 = gc * (K * c28 * (r9 - 1 * 4 * c18 - 2 * 6 * c28 - 3 * 4 * c38 - 4 * 1 * c48) - c38)
    du[15] = dc48 = gc * (K * c38 * (r9 - 1 * 4 * c18 - 2 * 6 * c28 - 3 * 4 * c38 - 4 * 1 * c48) - c48)
end

d_ics1 = Dict("R8"=> 0.00900023576123889, "c15"=> 1.761337400332505e-06, "P8"=> 2.9560034347325604e-06, "P5"=> 3.4020485285993286e-07, "c25"=> 9.118941517979407e-05, "c35"=> 0.004721137474792935, "A"=> 1.5e-323, "c18"=> 8.844673354103184e-06, "R5"=> 0.000500023861019239, "c38"=> 0.001867111128611335, "F"=> 1.5e-323, "r2"=> 0.3524443000060298, "r9"=> 0.33001432105772416, "c28"=> 0.0002687664938900425, "c48"=>0.001)
svlist = ["P5", "P8", "R5", "R8", "A", "F", "r2", "c15", "c25", "c35", "r9", "c18","c28", "c38", "c48"]
u0 = [d_ics1[x] for x in svlist]
tspan = (0.0,1000.0)
p=Dict("x"=>50
    ,"amp"=>0.001
    ,"muA"=>0.0
    ,"muF"=>0.0
    ,"s2"=>0.16
    ,"kR5"=>1.0
    ,"kr2"=>1.0
    ,"b15"=>1.0
    ,"b25"=>1.0
    ,"b35"=>0.6
    ,"a15"=>0.3
    ,"a25"=>0.5
    ,"a35"=>1.5
    ,"b18"=>2.7
    ,"b28"=>2.1
    ,"b38"=>0.7
    ,"b48"=>0.1
    ,"a18"=>1.0
    ,"a28"=>1.0
    ,"a38"=>1.5
    ,"a48"=>2.7
    ,"l05"=>1.8
    ,"l15"=>0.5
    ,"l25"=>0.0
    ,"l35"=>0.0
    ,"kP5"=>1.0
    ,"K"=>1000.0
    ,"maxA"=>2.3
    ,"maxF"=>2.1
    ,"s9"=>0.81
    ,"kR8"=>1.0
    ,"kR8"=>1.0
    ,"kr9"=>1.0
    ,"l08"=>3.9
    ,"l18"=>0.0
    ,"l28"=>0.0
    ,"l38"=>0.0
    ,"l48"=>0.0
    ,"kP8"=>1.0
    ,"gA"=>1.0
    ,"gF"=>1.0
    ,"gc"=>10.0
    ,"s05"=>0.06
    ,"s5"=>2.0
    ,"s08"=>0.0
    ,"s8"=>1.19
    ,"K5A"=>7.2
    ,"n5A"=>2.0
    ,"K58"=>0.006
    ,"n58"=>2.0
    ,"K8F"=>0.86
    ,"n8F"=>6.0
    ,"K9F"=>0.1
    ,"n9F"=>6.0
    ,"K2A"=>3.2
    ,"n2A"=>6.0
    ,"K25"=>1000
    ,"n25"=>6.0
    ,"K98"=>1000
    ,"n98"=>6.0
    ,"y"=>1)


text:
# These constants determine the order in which the g coefficients are stored in the nullspace matrix, and the corresponding moving line matrix. Do not change these. Every other function depends on this.
# later will change code so changing these constants will automagically work.
const lineindex = 3
const spatialindex = 1
const orderindex = 2

"spatial dimension of curve represented as an array of coefficients `x[i] = ∑Bj(θ)*x[i,j]` where `Bj(θ)` is the curve basis"
curvedimension(x::Array) = size(x, spatialindex)

"highest polynomial power of the curve represented as an array of coefficients `x[i] = ∑Bj(θ)*x[i,j]` where `Bj(θ)` is the curve basis"
curveorder(x::Array) = size(x, orderindex) - 1

numcoefficients(x::Array) = size(x, orderindex)

"spatial dimension of the moving line represented as an array of coefficients `g[i] = ∑Bl(θ)*gl[i,j]` where `Bl(θ)` is the polynomial basis"
linedimension(x::Array) = size(x, 1) - 1

"number of lines in moving line array"
numberoflines(x::Array) = size(x, lineindex)

movingline(x::Array, linenum::Integer) = x[:, :, linenum]

"returns a matrix expressing the relationship `[x(θ) 1]⋅g(θ) = 0`. The vectors in the right nullspace of this matrix contain the coefficients of the moving lines `gᵢ(θ)`."
function orthogonalitymatrix(curve::Array{T,2}, movinglineorder) where {T}
    x = curve
    qx = curveorder(x)
    qg = movinglineorder
    d = curvedimension(x)

    C = zeros(T, qx + qg + 1, (d + 1) * (qg + 1))

    for i in 1:d, j in 1:(qx + 1), l in 1:(qg + 1)
        C[(j - 1) + (l - 1) + 1, (i - 1) * (qg + 1) + l] = x[i, j]
    end

    # now set the final term which does not include any components of the curve x(theta)
    for row in 1:(qg + 1)
        col = d * (qg + 1) + row
        C[row, col] = one(T)
    end

    return C
end

"returns 3D array indexed like this: `x[line curve order,spatial dimension, line number]``"
function extractmovinglines(Vt, nullspacesize, movinglineorder, dimension)
    rows, _ = size(Vt)
    # temp = Vt[(rows-nullspacesize + 1):end,:]
    # println("rows $rows null $nullspacesize")
    # show(IOContext(stdout), "text/plain", temp)

    lines = permutedims(Vt[(rows - nullspacesize + 1):end, :]) # use this instead of transpose because transpose is recursive so it tries to transpose the elements of the array and crashes for Expr element types.
    # show(IOContext(stdout), "text/plain", lines)
    linecoeffs = movinglineorder + 1
    return permutedims(reshape(lines, linecoeffs, (dimension + 1), :), (2, 1, 3))
end

"`movinglines[:,i]` is the ith moving line. For `li = movinglines[:,i] (dimension+1,lineorder) = size(li)`. `rline[:,1] = pt1` and `rline[:,2] = pt2`. The line equation is `pt1 + alpha*pt2`."
function matricesforeigen(ray::AbstractRay{T,N}, movinglines::Array{T,3}) where {T,N}
    pt1 = origin(ray)
    pt2 = direction(ray)

    dimension = size(pt1, 1)
    numcoeffs = numcoefficients(movinglines)
    numlines = numberoflines(movinglines)

    A = Array{T,2}(undef, numcoeffs, numlines)
    B = Array{T,2}(undef, numcoeffs, numlines)

    for i in 1:numlines
        li = movinglines[:, :, i]

        for l in 1:numcoeffs
            asum = zero(T)
            bsum = zero(T)

            for j in 1:dimension
                asum += pt1[j] * li[j, l]
                bsum += pt2[j] * li[j, l]
            end

            asum += li[dimension + 1, l] # add in the final term of the sum that is not dotted with rline

            bsum = -bsum # reverse sign to get matrices in A-alpha*B form
            A[l, i] = asum
            B[l, i] = bsum
        end
    end
    return (A, B)
end

function matrixsizes(dimension, movinglineorder, curveorder)
    cols, rows = ((dimension + 1) * (movinglineorder + 1), curveorder + movinglineorder + 1)  # cols is the number of coefficients for the moving line. The line has d+1 terms for the line(plane/hyperplane) equation and movinglineorder+1 coefficients for the curve describing the moving line.
    nullspacesize = cols - rows

    return cols, rows, nullspacesize
end

"Evaluates a curve defined in the power basis. Curves and moving lines accessed like this: `[xi,ci]` where `xi` is the dimension index, and `ci` is the coefficient index."
function evaluatecurve(x::Array{T,2}, theta::Real) where {T<:Real}
    # dim,coefficients = size(x)
    dim = curvedimension(x)
    numcoeffs = numcoefficients(x)
    result = Array{T,1}(undef, dim)

    for dimension in 1:dim
        power = one(T)
        sum = zero(T)
        for coefficient in 1:numcoeffs
            sum += x[dimension, coefficient] * power
            power *= theta
        end
        result[dimension] = sum
    end
    return result
end

function validintersection(rline::AbstractRay{T,N}, linealpha, curve::Array{T,2}, curvetheta) where {T<:Real,N}
    if linealpha < 0
        return false
    end
    linepoint = point(rline, linealpha)
    curvepoint = evaluatecurve(curve, curvetheta)
    return isapprox(linepoint, curvepoint, rtol = 1e-10)
end

function eigenresults(rline::AbstractRay{T,N}, curve::Array{T,2}) where {T<:Real,N} # force rline and curve to use same number type to avoid expensive runtime conversion.
    dim = curvedimension(curve)
    orderofcurve = curveorder(curve)
    movinglineorder = max(orderofcurve, Int64(ceil(orderofcurve * 2 / dim - 1))) # Hve two constraints for movingline order. Must be big enough to result in a nullspace >= orderofcurve for the first step of computing the moving lines. Also need need at least as many eigenvalues as possible intersections of the line with the curve. Eigenmatrix is square of size qgxqg so need this matrix to be at least of size qx by qx. Take the greater of the the two constraint values.

    cmat = orthogonalitymatrix(curve, movinglineorder)

    fact = svd(cmat, full = true)
    nullspace = fact.Vt

    _, _, nullspacesize = matrixsizes(dim, movinglineorder, orderofcurve)

    movinglines = extractmovinglines(nullspace, nullspacesize, movinglineorder, dim)
    A, B = matricesforeigen(rline, movinglines)

    eig = eigen(A[:, 1:4]', B[:, 1:4]') # transpose so [A-lambdaB]'*Pl = 0

    return eig
end

"returns an array of intersection points. Each element in the array is (`[x,y,...],alpha,theta)` where `[x,y,...]` is the n-dimensional intersection point, alpha is the line parameter value at the intersection point, and theta is the curve parameter value at the intersection point"
function intersections(rline::AbstractRay{T,N}, curve::Array{T,2})::Array{Tuple{Array{T,1},T,T}} where {T<:Real,N}
    eigenstuff = eigenresults(rline, curve)
    result = Array{Tuple{Array{T,1},T,T}}(undef, 0)

    for index in CartesianIndices(eigenstuff.values)
        eigenvector = view(eigenstuff.vectors, :, index)
        eigenvalue = eigenstuff.values[index]
        # only take the real part of the solution. Due to roundoff might have real roots with very small complex part.
        theta = real(eigenvector[2] / eigenvector[1]) # this is only correct for curves defined in power basis. Once new Curve types are defined this function will be determined by the curve basis of the curve. Something like curveParameter(::CurveBasis,eigenvector)
        println(theta)
        if validintersection(rline, eigenvalue, curve, theta)
            push!(result, (point(rline, eigenvalue), eigenvalue, theta)) # return point on curve and alpha and theta parametric values corresponding to this point.
        end
    end

    return result
end
export intersections


text:
"""
```julia
gradient(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
```

A helper function for computing the gradient of an expression with respect to
an array of variable expressions.
"""
function gradient(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
    [expand_derivatives(Differential(v)(O),simplify) for v in vars]
end

"""
```julia
jacobian(ops::AbstractVector{<:Expression}, vars::AbstractVector{<:Expression}; simplify = true)
```

A helper function for computing the Jacobian of an array of expressions with respect to
an array of variable expressions.
"""
function jacobian(ops::AbstractVector{<:Expression}, vars::AbstractVector{<:Expression}; simplify = true)
    [expand_derivatives(Differential(v)(O),simplify) for O in ops, v in vars]
end

"""
```julia
sparsejacobian(ops::AbstractVector{<:Expression}, vars::AbstractVector{<:Expression}; simplify = true)
```

A helper function for computing the sparse Jacobian of an array of expressions with respect to
an array of variable expressions.
"""
function sparsejacobian(ops::AbstractVector{<:Expression}, vars::AbstractVector{<:Expression}; simplify = true)
    I = Int[]
    J = Int[]
    du = Expression[]

    sp = jacobian_sparsity(ops, vars)
    I,J,_ = findnz(sp)

    exprs = Expression[]

    for (i,j) in zip(I, J)
        push!(exprs, expand_derivatives(Differential(vars[j])(ops[i]), simplify))
    end
    sparse(I,J, exprs, length(ops), length(vars))
end

"""
```julia
jacobian_sparsity(ops::AbstractVector{<:Expression}, vars::AbstractVector{<:Expression})
```

Return the sparsity pattern of the Jacobian of an array of expressions with respect to
an array of variable expressions.
"""
function jacobian_sparsity(du, u)
    dict = Dict(zip(to_symbolic.(u), 1:length(u)))

    i = Ref(1)
    I = Int[]
    J = Int[]

    # This rewriter notes down which u's appear in a
    # given du (whose index is stored in the `i` Ref)
    r = [@rule ~x::(x->haskey(dict, x)) => begin
        push!(I, i[])
        push!(J, dict[~x])
        nothing
    end] |> Rewriters.Chain |> Rewriters.Postwalk

    for ii = 1:length(du)
        i[] = ii
        r(to_symbolic(du[ii]))
    end

    sparse(I, J, true, length(du), length(u))
end

"""
```julia
hessian(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
```

A helper function for computing the Hessian of an expression with respect to
an array of variable expressions.
"""
function hessian(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
    first_derivs = vec(jacobian([O], vars, simplify=simplify))
    n = length(vars)
    H = Array{Expression, 2}(undef,(n, n))
    fill!(H, 0)
    for i=1:n
        for j=1:i
            H[j, i] = H[i, j] = expand_derivatives(Differential(vars[i])(first_derivs[j]))
        end
    end
    H
end

isidx(x) = x isa TermCombination

let
    _scalar = one(TermCombination)

    linearity_propagator = [
          @rule +(~~xs) => reduce(+, filter(isidx, ~~xs), init=_scalar)
          @rule *(~~xs) => reduce(*, filter(isidx, ~~xs), init=_scalar)
          @rule (~f)(~x::(!isidx)) => _scalar
          @rule (~f)(~x::isidx) => if haslinearity(~f, Val{1}())
              combine_terms(linearity(~f, Val{1}()), ~x)
          else
              error("Function of unknown linearity used: ", ~f)
          end
          @rule (^)(~x::isidx, ~y) => ~y isa Number && isone(~y) ? ~x : (~x) * (~x)
          @rule (~f)(~x, ~y) => begin
              if haslinearity(~f, Val{2}())
                  a = isidx(~x) ? ~x : _scalar
                  b = isidx(~y) ? ~y : _scalar
                  combine_terms(linearity(~f, Val{2}()), a, b)
              else
                  error("Function of unknown linearity used: ", ~f)
              end
          end] |> Rewriters.Chain |> Rewriters.Postwalk |> Rewriters.Fixpoint

    # we do this in a let block so that Revise works on the list of rules
    global hessian_sparsity

    """
    ```julia
    hessian_sparsity(ops::AbstractVector{<:Expression}, vars::AbstractVector{<:Expression})
    ```

    Return the sparsity pattern of the Hessian of an array of expressions with respect to
    an array of variable expressions.
    """
    function hessian_sparsity(f, u)
        idx(i) = TermCombination(Set([Dict(i=>1)]))
        dict = Dict(SymbolicUtils.to_symbolic.(u) .=> idx.(1:length(u)))
        found = []
        f = Rewriters.Prewalk(x->haskey(dict, x) ? dict[x] : x)(to_symbolic(f))
        _sparse(linearity_propagator(to_symbolic(f)), length(u))
    end
end

"""
```julia
islinear(ex::Expression, u)
```
Check if an expression is linear with respect to a list of variable expressions.
"""
function islinear(ex::Expression, u)
    isempty(hessian_sparsity(ex, u).nzval)
end

"""
```julia
sparsehessian(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
```

A helper function for computing the sparse Hessian of an expression with respect to
an array of variable expressions.
"""
function sparsehessian(O::Expression, vars::AbstractVector{<:Expression}; simplify = true)
    S = hessian_sparsity(O, vars)
    I, J, _ = findnz(S)
    exprs = Array{Expression}(undef, length(I))
    fill!(exprs, 0)
    prev_j = 0
    d = nothing
    for (k, (i, j)) in enumerate(zip(I, J))
        j > i && continue
        if j != prev_j
            d = expand_derivatives(Differential(vars[j])(O), false)
        end
        expr = expand_derivatives(Differential(vars[i])(d), simplify)
        exprs[k] = expr
        prev_j = j
    end
    H = sparse(I, J, exprs, length(vars), length(vars))
    for (i, j) in zip(I, J)
        j > i && (H[i, j] = H[j, i])
    end
    return H
end

function simplified_expr(O::Operation)
  if O isa Constant
    return O.value
  elseif isa(O.op, Differential)
    return :(derivative($(simplified_expr(O.args[1])),$(simplified_expr(O.op.x))))
  elseif isa(O.op, Variable)
    isempty(O.args) && return O.op.name
    return Expr(:call, Symbol(O.op), simplified_expr.(O.args)...)
  end
  if O.op === (^)
      if length(O.args) > 1  && O.args[2] isa Constant && O.args[2].value < 0
          return Expr(:call, :^, Expr(:call, :inv, simplified_expr(O.args[1])), -(O.args[2].value))
      end
  end
  return Expr(:call, Symbol(O.op), simplified_expr.(O.args)...)
end

simplified_expr(c::Constant) = c.value

function simplified_expr(eq::Equation)
    Expr(:(=), simplified_expr(eq.lhs), simplified_expr(eq.rhs))
end

simplified_expr(eq::AbstractArray) = simplified_expr.(eq)
simplified_expr(x::Integer) = x
simplified_expr(x::AbstractFloat) = x


text:
"""
Overview: functions used to manage dofs handling
"""

"""
> mutable struct Field
- neq : number of equations
- nen : total number of entries
- cpn : nodal uknowns per node
- dof : dofs ordering
- val : field values
"""
mutable struct Field
  neq::Int64
  nen::Int64
  cpn::Int64
  dof::Vector{Int64}
  val::Vector{Float64}
end

"""
> Field(mesh::Grid, cpn::Int64)
Field constructor
- mesh : Grid
- cpn : number of nodal unknowns per node
"""
function Field(mesh::Grid, cpn::Int64)
  #
  nen = mesh.nn*cpn
  dof = zeros(Int64,nen)
  val = zeros(Float64,nen)
  # filter nodes without domain
  for i = 1:mesh.nΩ
    iΩ = mesh.Ω[i]
    nn = size(iΩ.e2n)[1]
    for j = 1:nn
      for k = 1:cpn
        dof[k+(iΩ.e2n[j]-1)*cpn] = 1
      end
    end
  end
  # fix Dirichlet boundary conditions
  for i = 1:mesh.nΓ
    iΓ = mesh.Γ[i]
    nn = size(iΓ.e2n)[1]
    bcd = iΓ.bcdofs
    bcv = iΓ.bcvals
    for j = 1:cpn
      if (bcd[j]>0)
        for k = 1:nn
          dof[j+(iΓ.e2n[k]-1)*cpn] = -1
          val[j+(iΓ.e2n[k]-1)*cpn] = bcv[j]
        end
      end
    end
  end
  # count number of equations
  neq = 0
  for i = 1:nen
    if (dof[i]>0)
      neq += 1
      dof[i] = neq
    end
  end
  #
  return Field(neq,nen,cpn,dof,val)
  #
end

"""
> dofs!(ϕ::Field, nn::Int64, nodes::Vector{Int64}, dofs::Vector{nt64})
It copies the dofs associated to nodes list nodes into dofs.
- ϕ : field of interest
- nn : number of queries
- nodes : list of nodes
- dofs : number of degrees of freedom
"""
function dofs!(ϕ::Field, nn::Int64, nodes::Vector{Int64},
               dofs::Vector{Int64})
  #
  for i = 1:nn
    for j = 1:ϕ.cpn
      @inbounds dofs[j+(i-1)*ϕ.cpn] = ϕ.dof[j+(nodes[i]-1)*ϕ.cpn]
    end
  end
  #
  return nothing
end

"""
> dofs!(ϕ::Field, nn::Int64, nodes::SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, dofs::Vector{nt64})
It copies the dofs associated to nodes list nodes into dofs.
This version accepts SubArrays.
- ϕ : field of interest
- nn : number of queries
- nodes : list of nodes
- dofs : number of degrees of freedom
"""
function dofs!(ϕ::Field, nn::Int64, nodes::SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true},
               dofs::Vector{Int64})
  #
  for i = 1:nn
    for j = 1:ϕ.cpn
      @inbounds dofs[j+(i-1)*ϕ.cpn] = ϕ.dof[j+(nodes[i]-1)*ϕ.cpn]
    end
  end
  #
  return nothing
end

"""
> dofs_vals!(ϕ::Field, nn::Int64, nodes::Vector{Int64}, dofs::Vector{Int64}, vals::Vector{Float64})
It copies the dofs and values associated to nodes list nodes into dofs and vals.
- ϕ : field of interest
- nn : number of queries
- nodes : list of nodes
- dofs : number of degrees of freedom
- vals : values of the field
"""
function dofs_vals!(ϕ::Field, nn::Int64, nodes::Vector{Int64},
                    dofs::Vector{Int64}, vals::Vector{Float64})
  #
  for i = 1:nn
    for j = 1:ϕ.cpn
      @inbounds dofs[j+(i-1)*ϕ.cpn] = ϕ.dof[j+(nodes[i]-1)*ϕ.cpn]
      @inbounds vals[j+(i-1)*ϕ.cpn] = ϕ.val[j+(nodes[i]-1)*ϕ.cpn]
    end
  end
  #
  return nothing
end

"""
> update_field!(ϕ::Field,U::Vector{Float64})
it overwrites free dofs values using the solution U
- ϕ : field to update
- U : solution vector
"""
function update_field!(ϕ::Field,U::Vector{Float64})
  #
  for i = 1:ϕ.nen
    if (ϕ.dof[i]>0)
      @inbounds ϕ.val[i] = U[ϕ.dof[i]]
    end
  end
  #
  return nothing
end

"""
> increment_field!(ϕ::Field,U::Vector{Float64})
it increments field values. ϕ.val += U
- ϕ : field to update
- U : solution vector
"""
function increment_field!(ϕ::Field,U::Vector{Float64})
  #
  for i = 1:ϕ.nen
    if (ϕ.dof[i]>0)
      @inbounds ϕ.val[i] += U[ϕ.dof[i]]
    end
  end
  #
  return nothing
end


text:
export Model

mutable struct Model{D <: Distribution, T <: AbstractFloat}
    ω::Vector{T}
    A::Dict{Int, Matrix{T}}
    B::Dict{Int, Matrix{T}}
    scaling::Real
end

function deepcopy(gas::Model{D, T}) where {D, T}
    return Model{D, T}(deepcopy(gas.ω), deepcopy(gas.A), deepcopy(gas.B), deepcopy(gas.scaling))
end

function create_ω(num_params::Int)
    return fill(NaN, num_params)
end

function create_lagged_matrix(lags::Vector{Int}, time_varying_params::Vector{Int}, zeros_params::Vector{T}) where T
    mat = Dict{Int, Matrix{T}}()
    for i in lags
        mat[i] = Diagonal(zeros_params)
        for k in time_varying_params
            mat[i][k, k] = NaN
        end
    end
    return mat
end

function Model(p::Int, q::Int, D::Type{<:Distribution}, scaling::Real;
             time_varying_params::Vector{Int} = collect(1:num_params(D)))

    # Vector of unkowns
    zeros_params = fill(0.0, num_params(D))
    ω = create_ω(num_params(D))

    # Create A and B
    A = create_lagged_matrix(collect(1:p), time_varying_params, zeros_params)
    B = create_lagged_matrix(collect(1:q), time_varying_params, zeros_params)

    return Model{D, Float64}(ω, A, B, scaling)
end

function Model(ps::Vector{Int}, qs::Vector{Int}, D::Type{<:Distribution}, scaling::Real;
             time_varying_params::Vector{Int} = collect(1:num_params(D)))

    # Vector of unkowns
    zeros_params = fill(0.0, num_params(D))
    ω = create_ω(num_params(D))

    # Create A and B
    A = create_lagged_matrix(ps, time_varying_params, zeros_params)
    B = create_lagged_matrix(qs, time_varying_params, zeros_params)

    return Model{D, Float64}(ω, A, B, scaling)
end

"""
    Model

The constructor of a score-driven model. The model receives the lag structure, the
distribution and the scaling. You can define the lag structure in two different
ways, either by passing integers `p` and `q` to add all lags from `1` to `p` and `1` to `q` or by
passing vectors of integers `ps` and `qs` containing the desired lags. Once you build
the model all of the unknown parameters that must be estimated are represented as `NaN`.

```jldoctest
# Passing p and q
julia> Model(2, 2, LogNormal, 0.5)
Model{LogNormal,Float64}([NaN, NaN], Dict(2=>[NaN 0.0; 0.0 NaN],1=>[NaN 0.0; 0.0 NaN]), Dict(2=>[NaN 0.0; 0.0 NaN],1=>[NaN 0.0; 0.0 NaN]), 0.5)

# Passing ps and qs
julia> Model([1, 12], [1, 12], Gamma, 0.0)
Model{Gamma,Float64}([NaN, NaN], Dict(12=>[NaN 0.0; 0.0 NaN],1=>[NaN 0.0; 0.0 NaN]), Dict(12=>[NaN 0.0; 0.0 NaN],1=>[NaN 0.0; 0.0 NaN]), 0.0)
```

If you don't want all the parameters to be considered time-varying you can express it
through the keyword argument `time_varying_params`, there you should pass a vector
containing a number that represents which parameter should be time-varying. As an example
in the Normal distribution `time_varying_params = [1]` indicates that only ``\\mu`` is
time-varying. You can find the table with the dictionary (number => parameter) in the
section [ScoreDrivenModels distributions](@ref).

```jldoctest
julia> Model([1, 12], [1, 12], Normal, 1.0; time_varying_params = [1])
Model{Normal,Float64}([NaN, NaN], Dict(12=>[NaN 0.0; 0.0 0.0],1=>[NaN 0.0; 0.0 0.0]), Dict(12=>[NaN 0.0; 0.0 0.0],1=>[NaN 0.0; 0.0 0.0]), 1.0)
```
"""
function Model end

function number_of_lags(gas::Model)
    return max(maximum(keys(gas.A)), maximum(keys(gas.B)))
end

"""
    Unknowns

Structure that stores the positions of the parameters to be estimated in the [`Model`](@ref).
"""
mutable struct Unknowns
    ω::Vector{Int}
    A::Dict{Int, Vector{Int}}
    B::Dict{Int, Vector{Int}}
end

function fill_psitilde!(gas::Model, psitilde::Vector{T}, unknowns::Unknowns) where T
    offset = 0
    # fill ω
    for i in unknowns.ω
        offset += 1
        @inbounds gas.ω[i] = psitilde[offset]
    end
    # fill A
    for (k, v) in unknowns.A
        for i in v
            offset += 1
            @inbounds gas.A[k][i] = psitilde[offset]
        end
    end
    # fill B
    for (k, v) in unknowns.B
        for i in v
            offset += 1
            @inbounds gas.B[k][i] = psitilde[offset]
        end
    end
    return
end

function find_unknowns(gas::Model)
    unknowns_A = Dict{Int, Vector{Int}}()
    unknowns_B = Dict{Int, Vector{Int}}()

    unknowns_ω = find_unknowns(gas.ω)

    for (k, v) in gas.A
        unknowns_A[k] = find_unknowns(v)
    end
    for (k, v) in gas.B
        unknowns_B[k] = find_unknowns(v)
    end
    return Unknowns(unknowns_ω, unknowns_A, unknowns_B)
end

function dim_unknowns(gas::Model)
    return length(find_unknowns(gas))
end

function length(unknowns::Unknowns)
    len = length(values(unknowns.ω))
    for (k, v) in unknowns.A
        len += length(v)
    end
    for (k, v) in unknowns.B
        len += length(v)
    end
    return len
end

function log_lik(psitilde::Vector{T}, y::Vector{T}, gas::Model{D, T},
                 initial_params::Matrix{T}, unknowns::Unknowns, n::Int) where {D, T}

    # Use the unkowns vectors to fill the right positions
    fill_psitilde!(gas, psitilde, unknowns)

    if isnan(initial_params[1]) # Means default stationary initialization
        params = score_driven_recursion(gas, y)
    else
        params = score_driven_recursion(gas, y; initial_params = initial_params)
    end

    return log_likelihood(D, y, params, n) / length(y)
end


text:
#=
Constants needed for various methods
=#

#Detection Power

const power_win=20
const power_win1=power_win+1
const power_win0=power_win-1

#Wavelets
#=
const coiflets_five=[0.00021208083980379827
0.00035858968789573785
-0.0021782363581090178
-0.004159358781386048
0.010131117519849788
0.023408156785839195
-0.02816802897093635
-0.09192001055969624
0.05204316317624377
0.4215662066908515
-0.7742896036529562
0.4379916261718371
0.06203596396290357
-0.10557420870333893
-0.0412892087501817
0.03268357426711183
0.01976177894257264
-0.009164231162481846
-0.006764185448053083
0.0024333732126576722
0.0016628637020130838
-0.0006381313430451114
-0.00030225958181306315
0.00014054114970203437
4.134043227251251e-05
-2.1315026809955787e-05
-3.7346551751414047e-06
2.0637618513646814e-06
1.6744288576823017e-07
-9.517657273819165e-08
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
];

#=
Multiscale Correlation of Wavelet Coefficients

Yang et al 2011
=#

const bigJ=20

const wave_a=collect(0.5:.1:1.5);

#This is REAL ugly, and edges aren't very good. Also code isn't very Julian

const coiflets_scaled=zeros(Float64,length(wave_a),bigJ)

coif_intp=interpolate(coiflets_five,BSpline(Quadratic(Line())),OnCell())
interpvec=collect(1:.1:40)
coif_up=Float64[coif_intp[i] for i in interpvec]

for i=1:length(wave_a)
    t=collect(1:bigJ)./wave_a[i]
    for j=1:bigJ
        coiflets_scaled[i,j]=coif_up[indmin(abs(interpvec-t[j]))]
    end
end

const onesquarea=Float64[1/sqrt(wave_a[i]) for i=1:11]
=#
#=
CLASSIT clustering algorithm
=#

const acuity=1.0 / sqrt(2.0 * pi) #Probably adjust this based on variance

const cob_cutoff=0.2

#=
Discrete Derivatives - Maximum Difference Test
=#

const DD_inds=[1; 3; 7]


text:
module PlasmoTest
#保留五位小数的函数
digit5(x) = round(x*100000)/100000
export digit5
end # module

max_stars_repo_name: BoyangSheng/Skill-studing

max_issues_repo_name: BoyangSheng/Skill-studying

max_forks_repo_name: BoyangSheng/Skill-studying


text:
# The List implementation is inspired by https://github.com/JuliaCollections/DataStructures.jl
export List, Nil, Cons, list, cons, nil, head, tail
abstract type List{T} end
struct Nil{T} <: List{T} end
struct Cons{T} <: List{T}
  head::T
  tail::List{T}
end

nil(T) = Nil{T}()
nil() = nil(Any)
cons(h, t::List{T}) where {T} = Cons{T}(h, t)
list() = nil()
Base.isempty(lst::Nil) = true
Base.isempty(lst::Cons) = false
Base.getindex(lst::Nil, i) = throw(BoundsError(lst, i))
Base.getindex(lst::Cons, i) = i == 1 ? lst.head : getindex(lst.tail, i-1)

list(elts...) =
  let l = nil()
    for i = length(elts):-1:1
      l = cons(elts[i], l)
    end
    l
  end

list(elts::T...) where {T} =
  let l = nil(T)
    for i = length(elts):-1:1
      l = cons(elts[i], l)
    end
    l
  end
head(x::Cons) = x.head
tail(x::Cons) = x.tail

import Base.==
==(x::Nil, y::Nil) = true
==(x::Cons, y::Cons) = (x.head == y.head) && (x.tail == y.tail)

Base.show(io::IO, lst::Nil) = print(io, "list()")
Base.show(io::IO, lst::Cons) =
  begin
    print(io, "list(")
    show(io, head(lst))
    for e in tail(lst)
      print(io, ", ")
      show(io, e)
    end
    print(io, ")")
  end

Base.length(l::Nil) = 0
Base.length(l::Cons) =
  let n = 0
    for i in l
      n += 1
    end
    n
  end

Base.map(f::Base.Callable, lst::List) = list((f(e) for e in lst)...)
Base.filter(f::Function, lst::List) = list((e for e in lst if f(e))...)

Base.cat() = list()
Base.cat(lst::List, lsts::List...) =
  let T = typeof(lst).parameters[1]
    n = length(lst)
    for l in lsts
      T2 = typeof(l).parameters[1]
      T = typejoin(T, T2)
      n += length(l)
    end
    elems = Vector{T}(undef, n)
    i = 1
    for e in lst
      elems[i] = e
      i += 1
    end
    for lst in lsts
      for e in lst
      elems[i] = e
      i += 1
      end
    end
    let l = nil(T)
      for i = i-1:-1:1
        l = cons(elems[i], l)
      end
      l
    end
  end

Base.iterate(l::List, ::Nil) = nothing
Base.iterate(l::List, state::Cons = l) = state.head, state.tail

# Lists can be converted to Arrays

Base.convert(::Type{Array{S,1}}, l::List{T}) where {S, T <: S} = collect(T, l)

export random, random_range, set_random_seed,
       division, map_division

#previous_random::Int = 12345
previous_random = 12345

set_random_seed(v::Int) =
  global previous_random = v

next_random(previous_random::Int) =
  let test = 16807*rem(previous_random,127773) - 2836*div(previous_random,127773)
    if test > 0
      if test > 2147483647
        test - 2147483647
      else
        test
      end
    else
      test + 2147483647
    end
  end

next_random!() =
  begin
    global previous_random = next_random(previous_random)
    previous_random
  end

random(x::Int) = rem(next_random!(), x)

random(x::Real) = x*next_random!()/2147483647.0

random_range(x0, x1) =
  if x0 == x1
    x0
  else
    x0 + random(x1 - x0)
  end

export RGB, rgb, rgb_radiance

struct RGB
  r::Real
  g::Real
  b::Real
end

rgb(r::Real=0, g::Real=0, b::Real=0) =
    RGB(convert(UInt8, r), convert(UInt8, g), convert(UInt8, b))

rgb_radiance(c::RGB) = 0.265*c.r+0.67*c.g+0.065*c.b

required() = error("Required parameter")

division(t0, t1, n::Real, include_last::Bool=true) =
  let n = convert(Int, n), iter = range(t0, stop=t1, length=n + 1)
    collect(include_last ? iter : take(iter, n))
  end

# Generic objects are processed using map_division
division(obj::Any, n::Real) = map_division(identity, obj, n)

map_division(f, t0, t1, n::Real, include_last::Bool=true) =
  let n = convert(Int, n), iter = range(t0, stop=t1, length=n + 1)
    map(f, include_last ? iter : take(iter, n))
  end

map_division(f, u0, u1, nu::Real, include_last_u::Bool, v0, v1, nv::Real) =
  map_division(u -> map_division(v -> f(u, v), v0, v1, nv),
               u0, u1, nu, include_last_u)

map_division(f, u0, u1, nu::Real, v0, v1, nv::Real, include_last_v::Bool=true) =
  map_division(u -> map_division(v -> f(u, v), v0, v1, nv, include_last_v),
               u0, u1, nu)

map_division(f, u0, u1, nu::Real, include_last_u::Bool, v0, v1, nv::Real, include_last_v::Bool) =
  map_division(u -> map_division(v -> f(u, v), v0, v1, nv, include_last_v),
               u0, u1, nu, include_last_u)

# Grasshopper compatibility

export series, crossref, remap, cull

series(start::Real, step::Real, count::Int) =
  range(start, step=step, length=count)

export crossref_holistic
crossref_holistic(arr1, arr2) =
  vcat([arr1[i] for i in range(1, stop=length(arr1)) for j in arr2]...),
  vcat([arr2 for i in arr1]...)

crossref(as, bs) = [(a, b) for a in as, b in bs]

remap(in, (min_in, max_in), (min_out, max_out)) =
  min_out + (max_out-min_out)/(max_in-min_in)*(in-min_in)

cull(template, as) =
  [a for (a, t) in zip(as, cycle(template)) if t]

## To present errors in the backends that call back to Julia
export errormsg

errormsg(e) =
  let b = backtrace()
    sprint((io, e) -> showerror(io, e, b, backtrace=true), e)
  end


text:
function StreamingAemWriter(arg0::Appendable, arg1::TimeScale, arg2::Map)
    return StreamingAemWriter((Appendable, TimeScale, Map), arg0, arg1, arg2)
end

function new_segment(obj::StreamingAemWriter, arg0::Map)
    return jcall(obj, "newSegment", StreamingAemWriter_AEMSegment, (Map,), arg0)
end

function write_header(obj::StreamingAemWriter)
    return jcall(obj, "writeHeader", void, ())
end


text:
using PartialSvdStoch, Test

using LowRankApprox

function testvrpca_withLowRankApprox(mat, k)
    #
    normMat = norm(mat)
    # get initial spectrum from psvd
    opts = LRAOptions(rank=k, rtol=0.001)
    @info "testvrpca_withLowRankApprox : initialization"
    @time Uguess, Sguess, Vguess = psvd(mat,opts)
    matGuess = Uguess * diagm(Sguess) * transpose(Vguess)
    deltaMatGuess = norm(matGuess - mat)
    @info "rank obtained : " size(Uguess)[2]
    @info " residual Error , initial L2 norm " deltaMatGuess normMat
    @info " relative error at initialization " deltaMatGuess/normMat
    #
    eta = 0.01
    batchSize = 20
    @info "doing vrpca svd  with rank initialization..." eta batchSize
    pcaPb = VrPCA(mat, k)
    @time iterateVrPCA(pcaPb, eta, batchSize)
    U = PartialSvdStoch.getindex(pcaPb,:U)
    S = PartialSvdStoch.getindex(pcaPb,:S)
    Vt = PartialSvdStoch.getindex(pcaPb,:Vt)
    # it happens that the obtained rank is less than asked for
    obtainedRank = pcaPb.k
    matVrpca = U.value * diagm(S.value) * Vt.value
    deltaMatIter = norm(matVrpca - mat)
    normMat = norm(mat)
    @info "deltaS after vr"  norm(matVrpca - mat) normMat
    relErr = deltaMatIter/normMat
    @info " relative error after vrpca iterations : " deltaMatIter/normMat
    relErr < 0.05 ? true : false
end

"""
# function testvrpca_withepsil(mat, epsil)

This test extracts the rank to get approximation at given precision epsil
see test : vrpca_epsil
"""
function testvrpca_withepsil(mat, epsil)
    #
    normMat = norm(mat)
    @info "testvrpca_withepsil : doing vrpca initialization with epsil" epsil
    pcaPb = VrPCA(mat, epsil)
    #
    eta = 0.01
    batchSize = 20
    @info "doing vrpca svd iteration" eta batchSize
    @time iterateVrPCA(pcaPb, eta, batchSize)
    #
    U = PartialSvdStoch.getindex(pcaPb,:U)
    S = PartialSvdStoch.getindex(pcaPb,:S)
    Vt = PartialSvdStoch.getindex(pcaPb,:Vt)
    matVrpca = U.value * diagm(S.value) * Vt.value
    relErr = norm(matVrpca - mat)/normMat
    @info "residual error after vr"  norm(matVrpca - mat) normMat relErr
    relErr < 0.05 ? true : false
end


text:
using SimpleDrawing, Plots, ProgressMeter

"""
    mechanism(a,b,offset,t)

Draw a picture of the spirograph mechanism. Draw with `spirograph` first, then follow
with `mechanism` using the same parameters.
"""
function mechanism(a::Int, b::Int, offset::Real, t::Real = 0)
    draw_circle(0, 0, a, linecolor = :black)
    θ = Spirograph.inner_angle(a, b, t)
    x = a + b
    z = (a + b) * exp(im * θ) * im
    draw_circle(z, abs(b), linecolor = :black)
    pen = Spirograph.spot(a, b, offset, t)
    draw_segment(z,pen,color=:black)
    draw_point(pen, color = :red, linecolor = :red, marker = 3, linewidth = 0)
    finish()
end

function spiro_movie(a, b, offset, step = Spirograph._DEFAULT_STEP)
    zvals = Spirograph.spiro_points(a, b, offset, step)
    T = Spirograph.stop_t(a, b)
    tvals = collect(0:step:T)
    npts = length(zvals)
    x = real(zvals)
    y = imag(zvals)

    bound = 0
    if b>0
        bound = a+b+max(abs(offset), b)
    else
        bound = a + max(abs(offset)-abs(b), 0)
    end

    PM = Progress(npts)

    movie = @animate for j = 1:npts
        newdraw()
        draw_rectangle(-bound, -bound, bound, bound, linecolor = :white)
        xx = x[1:j]
        yy = y[1:j]
        plot!(xx, yy, linecolor = :blue)
        mechanism(a, b, offset, tvals[j])
        next!(PM)
        finish()
    end
    return movie

end


text:
using Random
using Statistics

mutable struct ReplayBuffers
    writable_uuid::String
    d::Dict{String,ReplayBuffer}
end

function ReplayBuffers(writeable::ReplayBuffer)
    ReplayBuffers(
        string(writeable.uuid),
        Dict(
            string(writeable.uuid) => writeable
        )
    )
end

function getgamestepquality(x::ReplayBuffers, total_steps)
    avg_step = getavgstep(x)
    if ismissing(avg_step) || iszero(total_steps)
        return missing
    end
    round(avg_step / total_steps, digits=2)
end

function calcnumreanalyzed(x::ReplayBuffers)
    addends = gettotalreanalyses.(values(x.d))
    isempty(addends) ? 0 : sum(addends)
end

function getavgstep(x::ReplayBuffers)
    steps = getfield.(values(x.d), :step_mean)
    isempty(steps) ? 0 : round(mean(steps), digits=1)
end

function Base.length(x::ReplayBuffers)
    sum = 0
    for key in keys(x.d)
        sum += length(x.d[key])
    end
    sum
end

Base.lastindex(x::ReplayBuffers) = length(x)
Base.iterate(x::ReplayBuffers, state=1) = state > length(x) ? nothing : (x[state], state+1)
getbuffer(x::ReplayBuffers, uuid_str::String) = valornothing(x.d, uuid_str)
getwritablebuffer(x::ReplayBuffers) = getbuffer(x, x.writable_uuid)

function Base.getindex(x::ReplayBuffers, idx::Integer)
    for key in keys(x.d)
        buffer = x.d[key]
        len = length(buffer)
        if idx <= len
            return buffer[idx]
        end
        idx -= len
    end
    @assert false "This should not be accessible"
end

function Base.rand(x::ReplayBuffers)
    if isempty(x)
        return nothing
    end
    x[rand(1:length(x))]
end

function randtrainingbatch(x::ReplayBuffers, config::MuZeroConfig)
    while isempty(x)
    end
    batch_size = config.batch_size
    games = [rand(x) for _ in 1:batch_size]
    TrainingBatch(
        config.game_type,
        config.action_space_size,
        [TrainingSample(config, g) for g in games])
end

function setbuffer!(x::ReplayBuffers, buffer::ReplayBuffer)
    x.d[string(buffer.uuid)] = buffer
    nothing
end

function setbuffer!(x::ReplayBuffers, buffer_path::String)
    try
        setbuffer!(x, deserializeexpanduser(buffer_path))
    catch e
        @error "setbuffer! errored out: ($(buffer_path))" exception=(e, catch_backtrace())
    end
end

hasbuffer(x::ReplayBuffers, uuid_str::String) = haskey(x.d, uuid_str)


text:
# prediction from single tree - assign each observation to its final leaf
function predict!(pred, tree::Tree, X::AbstractArray{T, 2}) where {T<:Real}
    @inbounds @threads for i in 1:size(X, 1)
        id = 1
        x = view(X, i, :)
        while tree.nodes[id].split
            if x[tree.nodes[id].feat] < tree.nodes[id].cond
                id = tree.nodes[id].left
            else
                id = tree.nodes[id].right
            end
        end
        pred[i] += tree.nodes[id].pred
    end
    return pred
end

# prediction from single tree - assign each observation to its final leaf
function predict(tree::Tree, X::AbstractArray{T, 2}, K) where T<:Real
    pred = zeros(SVector{K,T}, size(X, 1))
    predict!(pred, tree, X)
    return pred
end

# prediction from single tree - assign each observation to its final leaf
function predict(model::GBTree, X::AbstractArray{T, 2}) where T<:Real
    pred = zeros(SVector{model.params.K,T}, size(X, 1))
    for tree in model.trees
        predict!(pred, tree, X)
    end
    pred = reinterpret(T, pred)
    if typeof(model.params.loss) == Logistic
        @. pred = sigmoid(pred)
    elseif typeof(model.params.loss) == Poisson
        @. pred = exp(pred)
    elseif typeof(model.params.loss) == Gaussian
        pred = transpose(reshape(pred, model.params.K, :))
        pred[:,2] = exp.(pred[:,2])
    elseif typeof(model.params.loss) == Softmax
        pred = transpose(reshape(pred, model.params.K, :))
        for row in eachrow(pred)
            row .= softmax(row)
        end
    end
    return pred
end

# prediction in Leaf - GradientRegression
function pred_leaf(loss::S, node::TrainNode{L,T}, params::EvoTreeRegressor, δ²) where {S<:GradientRegression,L,T}
    - params.η .* node.∑δ ./ (node.∑δ² .+ params.λ .* node.∑𝑤)
end

# prediction in Leaf - MultiClassRegression
function pred_leaf(loss::S, node::TrainNode{L,T}, params::EvoTreeRegressor, δ²) where {S<:MultiClassRegression,L,T}
    - params.η * node.∑δ ./ (node.∑δ² + params.λ * node.∑𝑤[1])
end

# prediction in Leaf - L1Regression
function pred_leaf(loss::S, node::TrainNode{L,T}, params::EvoTreeRegressor, δ²) where {S<:L1Regression,L,T}
    params.η .* node.∑δ ./ (node.∑𝑤 .* (1 .+ params.λ))
end

# prediction in Leaf - QuantileRegression
function pred_leaf(loss::S, node::TrainNode{L,T}, params::EvoTreeRegressor, δ²) where {S<:QuantileRegression,L,T}
    SVector{1,Float64}(params.η * quantile(reinterpret(Float64, δ²[collect(node.𝑖)]), params.α) / (1 + params.λ))
    # pred = params.η * quantile(δ²[collect(node.𝑖)], params.α) / (1 + params.λ)
end

# prediction in Leaf - GaussianRegression
function pred_leaf(loss::S, node::TrainNode{L,T}, params::EvoTreeRegressor, δ²) where {S<:GaussianRegression,L,T}
    - params.η * node.∑δ ./ (node.∑δ² + params.λ * node.∑𝑤[1])
end


text:
using UnicodePlots
using Base.Test
using Compat

# write your own tests here
include("tst_plots.jl")


text:
using Test
using EzXML
using GraphIO.GraphML

@testset "GraphML" begin
    for g in values(allgraphs)
        readback_test(GraphMLFormat(), g, testfail=true)
    end
    fname = joinpath(testdir, "testdata", "warngraph.graphml")

    @test_logs (:warn, "Skipping unknown node 'warnnode' - further warnings will be suppressed") match_mode=:any loadgraphs(fname, GraphMLFormat())
    @test_logs (:warn, "Skipping unknown XML element 'warnelement' - further warnings will be suppressed") match_mode=:any loadgraph(fname, "graph", GraphMLFormat())
    d = loadgraphs(fname, GraphMLFormat())
    write_test(GraphMLFormat(), d)
end


text:
module TowerTests

include("nodes-test.jl")

end  # module TowerTests