text:
module PSAPlots
    using RecipesBase
    using Pseudospectra
    using LinearAlgebra
    using PlotUtils
    using Printf

    import Pseudospectra: vec2ax, psa_compute
    import Pseudospectra: psmode_inv_lanczos, oneeigcond

    include("utils.jl")
    include("defaults.jl")
    #include("colormap.jl")

    include("recipes/spectralportrait.jl")
    include("recipes/plotmode.jl")
    include("recipes/mtxpowersplot.jl")
    include("recipes/mtxexpsplot.jl")
    include("recipes/surfplot.jl")

    export spectralportrait, spectralportrait!
    export plotmode, plotmode!
    export mtxpowersplot, mtxpowersplot!
    export mtxexpsplot, mtxexpsplot!
    export surfplot, surfplot!
end


text:
macro test_take(input, n)
    x = gensym()
    v = esc(input)
    quote
    br = Any[]
    for $x in take($v, $n)
        push!(br, $x)
    end
    mr = Any[]
    @itr for $x in take($v, $n)
        push!(mr, $x)
    end
    @test br == mr
    end
end

ext: jl


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

struct DomTreeNode
    level::Int
    children::Vector{Int}
end
DomTreeNode() = DomTreeNode(1, Vector{Int}())

struct DomTree
    idoms::Vector{Int}
    nodes::Vector{DomTreeNode}
end

"""
    Checks if bb1 dominates bb2
"""
function dominates(domtree::DomTree, bb1::Int, bb2::Int)
    bb1 == bb2 && return true
    target_level = domtree.nodes[bb1].level
    source_level = domtree.nodes[bb2].level
    source_level < target_level && return false
    for _ in (source_level - 1):-1:target_level
        bb2 = domtree.idoms[bb2]
    end
    return bb1 == bb2
end

bb_unreachable(domtree::DomTree, bb::Int) = bb != 1 && domtree.nodes[bb].level == 1

function update_level!(domtree::Vector{DomTreeNode}, node::Int, level::Int)
    domtree[node] = DomTreeNode(level, domtree[node].children)
    foreach(domtree[node].children) do child
        update_level!(domtree, child, level+1)
    end
end

struct DominatedBlocks
    domtree::DomTree
    worklist::Vector{Int}
end

function dominated(domtree::DomTree, root::Int)
    doms = DominatedBlocks(domtree, Vector{Int}())
    push!(doms.worklist, root)
    doms
end

function iterate(doms::DominatedBlocks, state::Nothing=nothing)
    isempty(doms.worklist) && return nothing
    bb = pop!(doms.worklist)
    for dominated in doms.domtree.nodes[bb].children
        push!(doms.worklist, dominated)
    end
    return (bb, nothing)
end

# Construct Dom Tree
# Simple algorithm - TODO: Switch to the fast version (e.g. https://tanujkhattar.wordpress.com/2016/01/11/dominator-tree-of-a-directed-graph/)
function construct_domtree(cfg::CFG)
    nblocks = length(cfg.blocks)
    dom_all = BitSet(1:nblocks)
    dominators = BitSet[n == 1 ? BitSet(1) : copy(dom_all) for n = 1:nblocks]
    changed = true
    while changed
        changed = false
        for n = 2:nblocks
            isempty(cfg.blocks[n].preds) && continue
            firstp, rest = Iterators.peel(Iterators.filter(p->p != 0, cfg.blocks[n].preds))
            new_doms = copy(dominators[firstp])
            for p in rest
                intersect!(new_doms, dominators[p])
            end
            push!(new_doms, n)
            changed = changed || (new_doms != dominators[n])
            dominators[n] = new_doms
        end
    end
    # Compute idoms
    idoms = fill(0, nblocks)
    for i = 2:nblocks
        doms = collect(dominators[i])
        for dom in doms
            i == dom && continue
            hasany = false
            for p in doms
                if p !== i && p !== dom && dom in dominators[p]
                    hasany = true; break
                end
            end
            hasany && continue
            idoms[i] = dom
        end
    end
    # Compute children
    domtree = DomTreeNode[DomTreeNode() for _ = 1:nblocks]
    for (idx, idom) in Iterators.enumerate(idoms)
        (idx == 1 || idom == 0) && continue
        push!(domtree[idom].children, idx)
    end
    # Recursively set level
    update_level!(domtree, 1, 1)
    DomTree(idoms, domtree)
end

max_stars_repo_name: DragonXDoom/julia

max_issues_repo_name: DragonXDoom/julia

max_forks_repo_name: DragonXDoom/julia


text:
import .Tree: ChildrenIndex, root_node

function _hierarchy(
    io::IO,
    cidx::ChildrenIndex,
    v::Int,
    indent::String,
    last::Bool,
    id::F,
) where {F}
    println(io, indent, last ? "└─" : "├─", id(v))
    indent = indent * (last ? "  " : "│ ")
    last_i = length(cidx[v])
    for (i, c) in enumerate(cidx[v])
        _hierarchy(io, cidx, c, indent, i == last_i, id)
    end
end

"""
    hierarchy([io,] cidx, [indent])

Print a tree given by a `ChildrenIndex` hierarchically onto `io` (default `stdout`).

# Example

```jldoctest
julia> import GraphIdx.Tree: ChildrenIndex, hierarchy

julia> hierarchy(ChildrenIndex([1, 1, 1, 3]))
1
├─2
└─3
  └─4

```
"""
function hierarchy(io::IO, cidx::ChildrenIndex, indent::String = "",
                   id::F = identity) where {F}
    local r = root_node(cidx)
    println(io, indent, id(r))
    last_i = length(cidx[r])
    for (i, c) in enumerate(cidx[r])
        _hierarchy(io, cidx, c, indent, i == last_i, id)
    end
end

hierarchy(cidx::ChildrenIndex, id::F = identity) where {F} =
    hierarchy(stdout, cidx, "", id)

"""
    hierarchy_string(cidx)

Same as `hierarchy` but output as String.
More precisely, print to a buffer and return as String.
"""
function hierarchy_string(cidx::ChildrenIndex, id::F = identity)::String where {F}
    buf = IOBuffer()
    hierarchy(buf, cidx, "", id)
    return String(take!(buf))
end

hierarchy_with(f::F, cidx::ChildrenIndex) where {F<:Function} =
    hierarchy(cidx, f)

"""
    parenthesis([io,] cidx [, stack])

Print the tree in parenthesis notation: recursively print `<node id>(subtree)`.
Per default output on `stdout`.

# Example

```jldoctest
julia> import GraphIdx.Tree: ChildrenIndex, parenthesis

julia> parenthesis(ChildrenIndex([1, 1, 1, 3]))
1(2()3(4()))

```
"""
function parenthesis(io::IO, cidx::ChildrenIndex, stack::Vector{Int} = Int[])
    dfs_walk(cidx, stack) do v::Int
        if v < 0    # just discovered
            print(io, -v, "(")
        else
            print(io, ")")
        end
    end
end

parenthesis(cidx::ChildrenIndex, stack::Vector{Int} = Int[]) =
    parenthesis(stdout, cidx, stack)

"""
    parenthesis_string(cidx [, stack])

Like [`parenthesis`](@ref) but return output as `String`.
"""
function parenthesis_string(cidx::ChildrenIndex, stack::Vector{Int} = Int[])::String
    buf = IOBuffer()
    parenthesis(buf, cidx, stack)
    return String(take!(buf))
end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
module TestGraphicalLinearAlgebra
using Test

using Random
using IterativeSolvers

using Catlab, Catlab.Theories, Catlab.WiringDiagrams, Catlab.Programs
using Catlab.LinearAlgebra.GraphicalLinearAlgebra
import LinearAlgebra: norm, svd

import LinearMaps: BlockDiagonalMap, UniformScaling
using LinearOperators

# Theories
###########

A, B = Ob(FreeLinearFunctions, :A, :B)
f, g = Hom(:f, A, B), Hom(:g, A, B)

# Domains and codomains
@test dom(plus(f,g)) == A
@test codom(plus(f,g)) == B

# Unicode syntax
@test f+g == plus(f,g)

# Evaluation
############

# LinearMaps instance
#--------------------

f, g = LinearMap([0 1; 1 0]), LinearMap([1 2 3; 4 5 6; 7 8 9])
h, k = LinearMap([1 -1; -1 1; 0 1]), LinearMap([1 1; 2 2; 3 3])

x, y = [2, 1], [7, 3, 5]
@test (f⋅f)*x == f*(f*x)
@test (f⋅h)*x == h*(f*x)
@test (f⊕g)*[x;y] == [f*x; g*y]
@test swap(dom(f),dom(g)) * [x;y] == [y;x]

@test mcopy(dom(f))*x == [x;x]
@test delete(dom(f))*x == []
@test plus(dom(f))*[x;x] == 2x
@test zero(dom(f))*Float64[] == zero(x)

@test (h+k)*x == h*x + k*x
@test scalar(dom(f),3)*x == 3*x
@test antipode(dom(f))*x == -1*x
@test adjoint(g)*y == g'*y

# GMRES with LinearMaps
#----------------------

V, W = Ob(FreeLinearFunctions, :V, :W)
f, g, h = Hom(:f, V, W), Hom(:g, V, W), Hom(:h, W, W)

M = plus(f,g)⋅h
# M = plus(f,g)

d = Dict(f=>LinearMap([1 1 1 1; 2 2 2 3; 3 3 3 3.0]),
         g=>LinearMap([3 1 1 1; -1 -1 -1 -1; 0 2 0 0.0]),
         h=>LinearMap([1 4 -1; -1 1 0; 1 1 -1.0]),
         V=>LinearMapDom(4),
         W=>LinearMapDom(3))
F(ex) = functor((LinearMapDom, LinearMap), ex, generators=d)

Random.seed!(1) # Set seed to avoid sporadic test failures.
n = 4
m = 3
x = ones(Float64, n)
A = F(M)'*F(M)
b = A*x
x̂ᴰ = Matrix(A)\b

@test IterativeSolvers.zerox(F(M), b) == IterativeSolvers.zerox(Matrix(F(M)), b)
x̂ = gmres(A, b; maxiter=20)

b̂ = A*x̂
@test norm(b̂-b) < 1e-8

n = 32
M = plus(f,g)⋅h
M₄ = foldl(⊕, [M, M, M, M, M, M, M, M])
x = ones(Float64, n)
# A needs to be square and nonsingular
A = ( F(M₄)'*F(M₄) ) + UniformScaling(n)
b = A*x
x̂ᴰ = Matrix(A)\b

x̂, h = gmres(A, b; tol=1e-12, maxiter=160, log=true, verbose=false)

b̂ = A*x̂
@test norm(b̂-b) < 1e-8
@test norm(x̂-x)/norm(x) < 1e-8

Σ = svdl(A; k=10)[1]

# Eigensystem tests

λ, Q = powm(F(M)'*F(M))
@test norm((F(M)'*F(M))*Q - λ*Q) < 1e-3
λ, Q = powm(F(M⋅adjoint(M)))
@test norm(F(M⋅adjoint(M))*Q - λ*Q) < 1e-3

Σ, L = svdl(F(M))
@test all(Σ .>= 0)

λ, Q = powm(A)
@test norm(A*Q - λ*Q) < 1e-8

Σ, L = svdl(A)
@test all(Σ .>= 0)

M′ = ( M⊕id(V) ) + ( id(V)⊕M )
M₂ = ( M′⋅adjoint(M′) ) + id(V⊕V)
@test all(svd(Matrix(F(M₂))).S .> 1e-10)
b = ones(size(F(M₂), 1))
x̂ = gmres(F(M₂), b)
b̂ = F(M₂)*x̂
@test norm(b̂ .- 1) < 1e-12

#x₂ = vcat(x,-x)
#b₂= F(M₂)*x₂
#x̂₂ = gmres(F(M₂), b₂)
#b̂₂ = F(M₂)*x̂₂
#@test norm(b̂₂-b₂) < 1e-4
#@test norm(x̂₂-x₂) < 1e-4

# LinearOps instance
#-------------------

f, g = LinearOperator([0 1; 1 0]), LinearOperator([1 2 3; 4 5 6; 7 8 9])
h, k = LinearOperator([1 -1; -1 1; 0 1]), LinearOperator([1 1; 2 2; 3 3])

x, y = [2, 1], [7, 3, 5]
@test (f⋅f)*x == f*(f*x)
@test (f⋅h)*x == h*(f*x)
@test (f⊕g)*[x;y] == [f*x; g*y]
@test swap(dom(f),dom(g)) * [x;y] == [y;x]

@test mcopy(dom(f))*x == [x;x]
@test delete(dom(f))*x == []
@test plus(dom(f))*[x;x] == 2x
@test zero(dom(f))*Float64[] == zero(x)

@test (h+k)*x == h*x + k*x
@test scalar(dom(f),3)*x == 3*x
@test antipode(dom(f))*x == -1*x
@test adjoint(g)*y == g'*y

# LinearOps functor
#------------------

A, B = Ob(FreeLinearFunctions, :A, :B)

f, g  = Hom(:f, A, B), Hom(:g, A, B)
h, k  = Hom(:f, B, A), Hom(:g, B, A)

fop, gop = LinearOperator([0 1 -1; -1 1 0]), LinearOperator([1 2 3; 4 5 6])
hop, kop = LinearOperator([1 -1; -1 1; 0 1]), LinearOperator([1 1; 2 2; 3 3])

d = Dict(f=>fop, h=>hop, g => gop, k => kop)

F(ex)=functor((LinearOpDom, LinearOperator), ex, generators=d)

@test Matrix(F(compose(f,h))*[1 1; 2 3; 4 5]) == Matrix(hop*fop*[1 1; 2 3; 4 5])
@test Matrix(F(compose(g,k))*[1 1; 2 3; 4 5]) == Matrix(kop*gop*[1 1; 2 3; 4 5])

# Catlab evaluate
#----------------

A, B = Ob(FreeLinearFunctions, :A, :B)
f, g, h, k = Hom(:f, A, A), Hom(:g, B, B), Hom(:h, A, B), Hom(:k, A, B)
generators = Dict(
  A => 2, B => 3,
  f => [0 1; 1 0],
  g => [1 2 3; 4 5 6; 7 8 9],
  h => [1 -1; -1 1; 0 1],
  k => [1 1; 2 2; 3 3],
)
ev = (args...) -> evaluate(args...; generators=generators)

val = arg -> generators[arg]
x, y = [2, 1], [7, 3, 5]
@test ev(f, x) == val(f)*x

@test ev(f⋅f, x) == val(f)*(val(f)*x)
@test ev(f⋅h, x) == val(h)*(val(f)*x)
@test ev(f⊕g, x, y) == (val(f)*x, val(g)*y)

@test ev(plus(A), [1,3], [2,4]) == [3,7]
@test ev(zero(A)) == [0,0]
@test ev(h+k, x) == val(h)*x + val(k)*x
@test ev(scalar(A,3),x) == 3*x
@test ev(antipode(A),x) == -1*x

# Wiring diagrams
#################

A, B = Ob(FreeLinearFunctions, :A, :B)
f, g, h, k = Hom(:f, A, A), Hom(:g, B, B), Hom(:h, A, B), Hom(:k, A, B)

d = to_wiring_diagram(f⋅h⋅g)
@test nboxes(d) == 3
@test dom(d) == Ports([:A])
@test codom(d) == Ports([:B])

@present Mat(FreeLinearFunctions) begin
    X::Ob
    Y::Ob
    Z::Ob
    V::Ob
    W::Ob

    f::Hom(X,Y)
    g::Hom(Y,Z)
    h::Hom(Y⊕Z, V)
    fˣ::Hom(X,X)
    fʸ::Hom(Y,Y)
    fᶻ::Hom(Z,Z)
    fᵛ::Hom(V,V)
    fʷ::Hom(W,W)
end

@test_skip A = @program Mat (x::X,y::Y) begin
    v = h(f(x),g(y))
    return v
end

@test_skip codom(A) == generator(Mat, :V)

@test_skip A = @program Mat (x::X,y::Y) begin
    v = h(f(x),g(y))
    return fᵛ(v)
end

@test_skip codom(A) == generator(Mat, :V)

end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    Arrow.BoolVector

A bit-packed array type, similar to [`ValidityBitmap`](@ref), but which
holds boolean values, `true` or `false`.
"""
struct BoolVector{T} <: ArrowVector{T}
    arrow::Vector{UInt8} # need to hold a reference to arrow memory blob
    pos::Int
    validity::ValidityBitmap
    ℓ::Int64
    metadata::Union{Nothing, Dict{String, String}}
end

Base.size(p::BoolVector) = (p.ℓ,)

@propagate_inbounds function Base.getindex(p::BoolVector{T}, i::Integer) where {T}
    @boundscheck checkbounds(p, i)
    if T >: Missing
        @inbounds !p.validity[i] && return missing
    end
    a, b = fldmod1(i, 8)
    @inbounds byte = p.arrow[p.pos + a - 1]
    # check individual bit of byte
    return ArrowTypes.fromarrow(T, getbit(byte, b))
end

@propagate_inbounds function Base.setindex!(p::BoolVector, v, i::Integer)
    @boundscheck checkbounds(p, i)
    x = convert(Bool, v)
    a, b = fldmod1(i, 8)
    @inbounds byte = p.arrow[p.pos + a - 1]
    @inbounds p.arrow[p.pos + a - 1] = setbit(byte, x, b)
    return v
end

arrowvector(::BoolKind, x::BoolVector, i, nl, fi, de, ded, meta; kw...) = x

function arrowvector(::BoolKind, x, i, nl, fi, de, ded, meta; kw...)
    validity = ValidityBitmap(x)
    len = length(x)
    blen = cld(len, 8)
    bytes = Vector{UInt8}(undef, blen)
    b = 0xff
    j = k = 1
    for y in x
        if y === false
            b = setbit(b, false, j)
        end
        j += 1
        if j == 9
            @inbounds bytes[k] = b
            b = 0xff
            j = 1
            k += 1
        end
    end
    if j > 1
        bytes[k] = b
    end
    return BoolVector{eltype(x)}(bytes, 1, validity, len, meta)
end

function compress(Z::Meta.CompressionType, comp, p::P) where {P <: BoolVector}
    len = length(p)
    nc = nullcount(p)
    validity = compress(Z, comp, p.validity)
    data = compress(Z, comp, view(p.arrow, p.pos:(p.pos + cld(p.ℓ, 8) - 1)))
    return Compressed{Z, P}(p, [validity, data], len, nc, Compressed[])
end

function makenodesbuffers!(col::BoolVector, fieldnodes, fieldbuffers, bufferoffset, alignment)
    len = length(col)
    nc = nullcount(col)
    push!(fieldnodes, FieldNode(len, nc))
    @debug 1 "made field node: nodeidx = $(length(fieldnodes)), col = $(typeof(col)), len = $(fieldnodes[end].length), nc = $(fieldnodes[end].null_count)"
    # validity bitmap
    blen = nc == 0 ? 0 : bitpackedbytes(len, alignment)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    # adjust buffer offset, make primitive array buffer
    bufferoffset += blen
    blen = bitpackedbytes(len, alignment)
    push!(fieldbuffers, Buffer(bufferoffset, blen))
    @debug 1 "made field buffer: bufferidx = $(length(fieldbuffers)), offset = $(fieldbuffers[end].offset), len = $(fieldbuffers[end].length), padded = $(padding(fieldbuffers[end].length, alignment))"
    return bufferoffset + blen
end

function writebuffer(io, col::BoolVector, alignment)
    @debug 1 "writebuffer: col = $(typeof(col))"
    @debug 2 col
    writebitmap(io, col, alignment)
    n = Base.write(io, view(col.arrow, col.pos:(col.pos + cld(col.ℓ, 8) - 1)))
    return n + writezeros(io, paddinglength(n, alignment))
end


text:
@testset "binary_trend" begin
    # sample data
    function generate_binary_data(nsamples)
        # generate predictions
        dist = Dirichlet(2, 1)
        predictions = [rand(dist) for _ in 1:nsamples]

        # generate targets
        targets_consistent = [rand() < predictions[i][1] ? 1 : 2 for i in 1:nsamples]
        targets_onlytwo = fill(2, nsamples)

        return (predictions, targets_consistent), (predictions, targets_onlytwo)
    end
    data_consistent, data_only_two = generate_binary_data(500)

    # define tensor product kernel (using the mean total variation distance as bandwidth)
    kernel = (ExponentialKernel() ∘ ScaleTransform(3)) ⊗ WhiteKernel()

    @testset "Consistency test" begin
        # define estimators
        estimators = (
            BiasedSKCE(kernel),
            UnbiasedSKCE(kernel),
            (BlockUnbiasedSKCE(kernel, b) for b in (2, 10, 50, 100))...,
        )

        for estimator in estimators
            test_consistent = @inferred(ConsistencyTest(estimator, data_consistent...))
            @test @inferred(pvalue(test_consistent)) > 0.6
            print(test_consistent)

            test_only_two = @inferred(ConsistencyTest(estimator, data_only_two...))
            @test @inferred(pvalue(test_only_two)) < 1e-6
            print(test_only_two)
        end
    end

    @testset "Distribution-free tests" begin
        # define estimators
        estimators = (
            BiasedSKCE(kernel),
            UnbiasedSKCE(kernel),
            (BlockUnbiasedSKCE(kernel, b) for b in (2, 10, 50, 100))...,
        )

        for estimator in estimators
            test_consistent = @inferred(
                DistributionFreeSKCETest(estimator, data_consistent...)
            )
            @test @inferred(pvalue(test_consistent)) > 0.7
            println(test_consistent)

            test_only_two = @inferred(DistributionFreeSKCETest(estimator, data_only_two...))
            if estimator isa Union{UnbiasedSKCE,BlockUnbiasedSKCE}
                @test @inferred(pvalue(test_only_two)) < 0.4
            else
                @test @inferred(pvalue(test_only_two)) < 1e-6
            end
            println(test_only_two)
        end
    end

    @testset "Asymptotic block SKCE test" begin
        for blocksize in (2, 10, 50, 100)
            test_consistent = @inferred(
                AsymptoticBlockSKCETest(kernel, blocksize, data_consistent...)
            )
            @test @inferred(pvalue(test_consistent)) > 0.7
            println(test_consistent)

            test_only_two = @inferred(
                AsymptoticBlockSKCETest(kernel, blocksize, data_only_two...)
            )
            @test @inferred(pvalue(test_only_two)) < 1e-6
            println(test_only_two)
        end
    end

    @testset "Asymptotic SKCE test" begin
        test_consistent = @inferred(AsymptoticSKCETest(kernel, data_consistent...))
        @test @inferred(pvalue(test_consistent)) > 0.7
        println(test_consistent)

        test_only_two = @inferred(AsymptoticSKCETest(kernel, data_only_two...))
        @test @inferred(pvalue(test_only_two)) < 1e-6
        println(test_only_two)
    end

    @testset "Asymptotic CME test" begin
        # define estimator (uniformly distributed test locations)
        testpredictions = [rand(Dirichlet(2, 1)) for _ in 1:5]
        testtargets = rand(1:2, 5)
        estimator = @inferred(UCME(kernel, testpredictions, testtargets))

        test_consistent = @inferred(AsymptoticCMETest(estimator, data_consistent...))
        @test @inferred(pvalue(test_consistent)) > 0.8
        println(test_consistent)

        test_only_two = @inferred(AsymptoticCMETest(estimator, data_only_two...))
        @test @inferred(pvalue(test_only_two)) < 1e-6
        println(test_only_two)
    end
end


text:
struct BasicFootTrajectory{T}
    p0::SVector{3, T}
    xaxis::SVector{2, T}
    Δx::Polynomial{4, T}
    Δz::Polynomial{3, T}
    s::Polynomial{6, T}
    tf::T
end

function BasicFootTrajectory(t0::Number, tf::Number, p0::SVector{3}, Δzmid::Number, pf::SVector{3}, zdf::Number)
    # Coordinates of final position relative to initial, projected onto swing plane (x-z)
    Δpf = pf - p0
    Δxf = sqrt(Δpf[1]^2 + Δpf[2]^2)
    T = eltype(Δxf)
    Δzf = convert(T, Δpf[3])
    xaxis = SVector(Δpf[1], Δpf[2]) / Δxf

    # Trajectories as function of interpolation variable `s`
    s0 = zero(T)
    sf = oneunit(T)
    smid = (s0 + sf) / 2
    Δx = fit_cubic(x0=s0, xf=sf, y0=zero(T), yd0=zero(T), yf=Δxf, ydf=zero(T))
    Δz = let
        num_coeffs = Val(3)
        A = hcat(
            SVector(SUP.coefficient_gradient(s0, num_coeffs, Val(0))),
            SVector(SUP.coefficient_gradient(smid, num_coeffs, Val(0))),
            SVector(SUP.coefficient_gradient(sf, num_coeffs, Val(0))),
        ) |> transpose
        b = SVector(tuple(zero(T), max(0, Δzf) + Δzmid, Δzf))
        Polynomial(Tuple(A \ b))
    end

    # Final value of derivative of height as a function of `s`
    dzds_f = derivative(Δz)(sf)

    # Final value of derivative of s as a function of t
    sdf = zdf / dzds_f

    # Interpolation trajectory
    s = fit_quintic(x0=t0, xf=tf, y0=s0, yd0=zero(s0), ydd0=zero(s0), yf=sf, ydf=sdf, yddf=zero(sf))

    BasicFootTrajectory(p0, xaxis, Δx, Δz, s, tf)
end

function Base.convert(::Type{BasicFootTrajectory{T}}, p::SVector{3}, tf::Number) where {T}
    xaxis = SVector(zero(T), zero(T))
    Δx = Polynomial(ntuple(i -> zero(T), Val(4)))
    Δz = Polynomial(ntuple(i -> zero(T), Val(3)))
    s = Polynomial(ntuple(i -> zero(T), Val(6)))
    BasicFootTrajectory(p, xaxis, Δx, Δz, s, tf)
end

function (traj::BasicFootTrajectory)(t::Number, ::Val{2})
    s = traj.s(t)
    sd = derivative(traj.s)(t)
    sdd = derivative(derivative(traj.s))(t)
    Δx = traj.Δx(s)
    Δxd = derivative(traj.Δx)(s) * sd
    Δxdd = derivative(traj.Δx)(s) * sdd + derivative(derivative(traj.Δx))(s) * sd^2
    Δz = traj.Δz(s)
    Δzd = derivative(traj.Δz)(s) * sd
    Δzdd = derivative(traj.Δz)(s) * sdd + derivative(derivative(traj.Δz))(s) * sd^2
    p = traj.p0 + SVector(Δx * traj.xaxis[1], Δx * traj.xaxis[2], Δz)
    pd = SVector(Δxd * traj.xaxis[1], Δxd * traj.xaxis[2], Δzd)
    pdd = SVector(Δxdd * traj.xaxis[1], Δxdd * traj.xaxis[2], Δzdd)
    p, pd, pdd
end

(traj::BasicFootTrajectory)(t::Number) = first(traj(t, Val(2)))


text:
function nlp_cvx_103_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - linear objective
    # - intersection convex quadratic constraints
    # Variants
    #   010 - one binding constraint (inflection point)
    #   011 - one binding constraint (inflection point)
    #   012 - one binding constraint (non-inflection point)
    #   013 - one binding constraint (non-inflection point)
    #   014 - intersection point

    model = Model(optimizer)

    @variable(model, x)
    @variable(model, y)

    @objective(model, Min, y)
    @NLconstraint(model, x^2 <= y)
    @NLconstraint(model, -x^2 + 1 >= y)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, 0, tol = objective_tol)
    return check_solution([x, y], [0, 0], tol = primal_tol)
end


text:
"""
    report(res::IS.Results, out_path::String, design_template::String)

This function uses weave to either generate a LaTeX or HTML
file based on the report_design.jmd (julia markdown) file
that it reads. Out_path in the weave function dictates
where the created file gets exported.

# Arguments
- `results::IS.Results`: The results to be plotted
- `out_path::String`: folder path to the location the report should be generated
- `design_template::String = "file_path"`: directs the function to the julia markdown report design, the default

# Example
```julia
results = solve_op_problem!(OpModel)
out_path = "/Users/downloads"
report(results, out_path, template)
```

# Accepted Key Words
- `doctype::String = "md2html"`: create an HTML, default is PDF via latex
- `backend::Plots.backend() = plotlyjs()`: sets the plots backend, default is gr()
"""
function report(res::IS.Results, out_path::String, design_template::String; kwargs...)
    doctype = get(kwargs, :doctype, "md2pdf")
    backend = get(kwargs, :backend, Plots.gr())
    initial_time = get(kwargs, :initial_time, nothing)
    len = get(kwargs, :horizon, nothing)

    !isfile(design_template) &&
        throw(ArgumentError("The provided template file is invalid"))
    args = Dict("results" => res, "backend" => backend)
    Weave.weave(
        design_template,
        out_path = out_path,
        latex_cmd = ["xelatex"],
        doctype = doctype,
        args = args,
    )
end


text:
function test_curtailableLoad()

    T = 2
    load = [2.0, 3.0]

    l = DR.DER.CurtailableLoad(index=0, T=T, load=load, binflag=false)

    model = GLPK.Optimizer()

    lmin = [1.0, 0.0]
    lmax = [10.0, 10.0]
    price = [1.0, -0.01]

    # Instantiate initial model
    h = DR.DER.House(
        0, T,
        lmin, lmax, price,
        [l]
    )

    model, var2idx, con2idx = DR.DER.build_model!(h, model)

    MOI.optimize!(model)

    @test MOI.get(model, MOI.ObjectiveValue()) ≈ 1.0-0.03

    for t in 1:T
        pnet = var2idx[(:house, h.index, :pnet, t)]
        pres = var2idx[(:curt, l.index, :pnet, t)]
        if price[t] >= 0.0
            @test MOI.get(model, MOI.VariablePrimal(), pnet) ≈ lmin[t]
            @test MOI.get(model, MOI.VariablePrimal(), pres) ≈ lmin[t]
        else
            @test MOI.get(model, MOI.VariablePrimal(), pnet) ≈ min(lmax[t], load[t])
            @test MOI.get(model, MOI.VariablePrimal(), pres) ≈ min(lmax[t], load[t])
        end
    end

    return nothing

end

@testset "Curt." begin
    test_curtailableLoad()
end


text:
"""
Based on
https://medium.com/vooban-ai/hyperopt-tutorial-for-optimizing-neural-networks-hyperparameters-e3102814b919
"""
module TestStatusFailSkip

using Test
using TreeParzen

function objective_fail(space)::Float64
    x = space[:x]
    y = space[:y]

    if y > 1
        # Make use of status fail as an example of skipping on error
        return Inf
    end
    return x^2 + y^2
end

space = Dict(
    :x => HP.Uniform(:x, -5.0, 5.0),
    :y => HP.Uniform(:y, -5.0, 5.0),
)
best = fmin(objective_fail, space, 100)
@test isapprox(best[:x], 0.051738577119312464, rtol = 1e2)
@test isapprox(best[:y], -0.09365785427863083, rtol = 1e2)

space = Dict(
    :x => HP.Uniform(:x, -5.0, 5.0),
)

# test to catch ArgumentError in evaluate hyperparams
@test_throws ArgumentError fmin(objective_fail, space, 10)

end
true


text:
using JsonGrinder, JSON, Test, SparseArrays, Mill
using HierarchicalUtils
import HierarchicalUtils: printtree

@testset "ExtractMultipleRepresentation" begin
	ex = MultipleRepresentation((ExtractCategorical(["Olda", "Tonda", "Milda"]),
		JsonGrinder.ExtractString()))
	e = ex("Olda")

	@test length(e.data) == 2
	@test size(e.data[1].data) == (4, 1)
	@test e.data[1].data == MaybeHotMatrix([2],4)
	@test size(e.data[2].data) == (2053, 1)
	@test findall(x->x > 0, SparseMatrixCSC(e.data[2].data)) .|> Tuple == [(206, 1), (272, 1), (624, 1), (738, 1), (1536, 1), (1676, 1)]

	e = ex(extractempty)
	@test nobs(e) == 0
	@test nobs(e[:e1]) == 0
	@test e[:e1].data isa MaybeHotMatrix{Int64,Array{Int64,1},Int64,Bool}
	@test nobs(e[:e2]) == 0
	@test e[:e2].data isa NGramMatrix{String,Array{String,1},Int64}

	ex2 = MultipleRepresentation((ExtractCategorical(["Olda", "Tonda", "Milda"]),
		JsonGrinder.ExtractString()))
	@test hash(ex) === hash(ex2)
	@test ex == ex2

	ex3 = MultipleRepresentation((ExtractCategorical(["Polda", "Tonda", "Milada"]),
		JsonGrinder.ExtractString()))
	@test hash(ex) !== hash(ex3)
	@test ex != ex3
end

@testset "show" begin
	ex = MultipleRepresentation((ExtractCategorical(["Olda", "Tonda", "Milda"]),
		JsonGrinder.ExtractString()))
	e = ex("Olda")

	@test buf_printtree(ex, trav=true) ==
	"""
	MultiRepresentation [""]
	  ├── e1: Categorical d = 4 ["E"]
	  └── e2: String ["U"]"""

	@test buf_printtree(e, trav=true) ==
	"""
	ProductNode with 1 obs [""]
	  ├── e1: ArrayNode(4×1 MaybeHotMatrix with Bool elements) with 1 obs ["E"]
	  └── e2: ArrayNode(2053×1 NGramMatrix with Int64 elements) with 1 obs ["U"]"""
end


text:
AbstractMeasurement = MonteCarlo.AbstractMeasurement
ConfigurationMeasurement = MonteCarlo.ConfigurationMeasurement

IsingMeasurement = MonteCarlo.IsingMeasurement
IsingEnergyMeasurement = MonteCarlo.IsingEnergyMeasurement
IsingMagnetizationMeasurement = MonteCarlo.IsingMagnetizationMeasurement

GreensMeasurement = MonteCarlo.GreensMeasurement
BosonEnergyMeasurement = MonteCarlo.BosonEnergyMeasurement

AbstractObservable = MonteCarloObservable.AbstractObservable

struct DummyModel <: MonteCarlo.Model end
struct DummyMeasurement <: AbstractMeasurement end

@testset "Type Hierarchy" begin
    @test IsingEnergyMeasurement <: IsingMeasurement
    @test IsingMagnetizationMeasurement <: IsingMeasurement
    @test IsingMeasurement <: AbstractMeasurement
    @test GreensMeasurement <: AbstractMeasurement
    @test BosonEnergyMeasurement <: AbstractMeasurement
end

@testset "Interface" begin
    m = DummyMeasurement()
    dummy_model = DummyModel()
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=1.0)

    @test nothing == MonteCarlo.prepare!(m, mc, model)
    @test_throws MethodError MonteCarlo.measure!(m, mc, model, 1)
    @test nothing == MonteCarlo.finish!(m, mc, model)
    @test MonteCarlo.default_measurements(mc, dummy_model) == Dict{Symbol, AbstractMeasurement}()
end

@testset "Checking defaults" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=1.0)

    defaults = MonteCarlo.default_measurements(mc, model)
    @test !isempty(defaults)

    @test haskey(defaults, :conf) && defaults[:conf] isa ConfigurationMeasurement
    obs = observables(defaults[:conf])
    @test haskey(obs, "Configurations") && obs["Configurations"] isa AbstractObservable

    @test haskey(defaults, :Magn) && defaults[:Magn] isa IsingMagnetizationMeasurement
    obs = observables(defaults[:Magn])
    @test haskey(obs, "Magnetization per site") && obs["Magnetization per site"] isa AbstractObservable
    @test haskey(obs, "Magnetic susceptibility") && obs["Magnetic susceptibility"] isa AbstractObservable
    @test haskey(obs, "Total magnetization") && obs["Total magnetization"] isa AbstractObservable
    @test haskey(obs, "Total magnetization squared") && obs["Total magnetization squared"] isa AbstractObservable

    @test haskey(defaults, :Energy) && defaults[:Energy] isa IsingEnergyMeasurement
    obs = observables(defaults[:Energy])
    @test haskey(obs, "Energy per site") && obs["Energy per site"] isa AbstractObservable
    @test haskey(obs, "Specific heat") && obs["Specific heat"] isa AbstractObservable
    @test haskey(obs, "Total energy") && obs["Total energy"] isa AbstractObservable
    @test haskey(obs, "Total energy squared") && obs["Total energy squared"] isa AbstractObservable

    @test isempty(mc.thermalization_measurements)
    @test !isempty(mc.measurements)
    @test haskey(mc.measurements, :conf) && mc.measurements[:conf] isa ConfigurationMeasurement
    @test haskey(mc.measurements, :Magn) && mc.measurements[:Magn] isa IsingMagnetizationMeasurement
    @test haskey(mc.measurements, :Energy) && mc.measurements[:Energy] isa IsingEnergyMeasurement

    model = HubbardModelAttractive(dims=2, L=2)
    mc = DQMC(model, beta=1.0)

    defaults = MonteCarlo.default_measurements(mc, model)
    @test !isempty(defaults)
    @test haskey(defaults, :conf) && defaults[:conf] isa ConfigurationMeasurement

    @test haskey(defaults, :Greens) && defaults[:Greens] isa GreensMeasurement
    obs = observables(defaults[:Greens])
    @test haskey(obs, "Equal-times Green's function") && obs["Equal-times Green's function"] isa AbstractObservable

    @test haskey(defaults, :BosonEnergy) && defaults[:BosonEnergy] isa BosonEnergyMeasurement
    obs = observables(defaults[:BosonEnergy])
    @test haskey(obs, "Bosonic Energy") && obs["Bosonic Energy"] isa AbstractObservable

    @test isempty(mc.thermalization_measurements)
    @test !isempty(mc.measurements)
    @test haskey(mc.measurements, :conf) && mc.measurements[:conf] isa ConfigurationMeasurement
    @test haskey(mc.measurements, :Greens) && mc.measurements[:Greens] isa GreensMeasurement
    @test haskey(mc.measurements, :BosonEnergy) && mc.measurements[:BosonEnergy] isa BosonEnergyMeasurement
end

@testset "Retrieving Measurements" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=1.0)

    ms = MonteCarlo.measurements(mc)
    @test mc.thermalization_measurements == ms[:TH]
    @test mc.measurements == ms[:ME]

    obs = MonteCarlo.observables(mc)
    @test keys(obs[:TH]) == keys(ms[:TH])
    @test keys(obs[:ME]) == keys(ms[:ME])

    @test haskey(obs[:ME][:conf], "Configurations")
    @test typeof(obs[:ME][:conf]["Configurations"]) <: AbstractObservable

    @test haskey(obs[:ME][:Energy], "Total energy")
    @test typeof(obs[:ME][:Energy]["Total energy"]) <: AbstractObservable
    @test haskey(obs[:ME][:Energy], "Total energy squared")
    @test typeof(obs[:ME][:Energy]["Total energy squared"]) <: AbstractObservable
    @test haskey(obs[:ME][:Energy], "Energy per site")
    @test typeof(obs[:ME][:Energy]["Energy per site"]) <: AbstractObservable
    @test haskey(obs[:ME][:Energy], "Specific heat")
    @test typeof(obs[:ME][:Energy]["Specific heat"]) <: AbstractObservable

    @test haskey(obs[:ME][:Magn], "Total magnetization")
    @test typeof(obs[:ME][:Magn]["Total magnetization"]) <: AbstractObservable
    @test haskey(obs[:ME][:Magn], "Total magnetization squared")
    @test typeof(obs[:ME][:Magn]["Total magnetization squared"]) <: AbstractObservable
    @test haskey(obs[:ME][:Magn], "Magnetization per site")
    @test typeof(obs[:ME][:Magn]["Magnetization per site"]) <: AbstractObservable
    @test haskey(obs[:ME][:Magn], "Magnetic susceptibility")
    @test typeof(obs[:ME][:Magn]["Magnetic susceptibility"]) <: AbstractObservable
end

@testset "Interacting with Measurements" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=1.0)

    @test length(mc.measurements) == 3
    delete!(mc, AbstractMeasurement)
    @test isempty(mc.measurements)

    @test_throws ErrorException push!(mc, :conf => Int64)
    @test_throws MethodError MonteCarlo.unsafe_push!(mc, :conf => 1.0)
    @test_throws ErrorException push!(mc, :conf => ConfigurationMeasurement, :bad_stage)
    @test_throws ErrorException MonteCarlo.unsafe_push!(mc, :conf => ConfigurationMeasurement(mc, model), :bad_stage)
    @test_throws ErrorException delete!(mc, :conf, :bad_stage)
    @test_throws ErrorException delete!(mc, ConfigurationMeasurement, :bad_stage)

    push!(mc, :conf => ConfigurationMeasurement)
    @test haskey(mc.measurements, :conf) && mc.measurements[:conf] isa ConfigurationMeasurement
    delete!(mc, :conf)
    @test !haskey(mc.measurements, :conf)

    MonteCarlo.unsafe_push!(mc, :Magn => IsingMagnetizationMeasurement(mc, model))
    @test haskey(mc.measurements, :Magn) && mc.measurements[:Magn] isa IsingMagnetizationMeasurement
    delete!(mc, IsingMagnetizationMeasurement)
    @test !haskey(mc.measurements, :Magn)

    @test isempty(mc.thermalization_measurements)
    push!(mc, :conf => ConfigurationMeasurement, :TH)
    @test haskey(mc.thermalization_measurements, :conf) &&
        mc.thermalization_measurements[:conf] isa ConfigurationMeasurement
    delete!(mc, :conf, :TH)

    MonteCarlo.unsafe_push!(mc, :conf => ConfigurationMeasurement(mc, model), :TH)
    @test haskey(mc.thermalization_measurements, :conf) &&
        mc.thermalization_measurements[:conf] isa ConfigurationMeasurement
    delete!(mc, ConfigurationMeasurement, :TH)
    @test !haskey(mc.thermalization_measurements, :TH)
end

@testset "Saving and Loading" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=1.0)
    run!(mc, thermalization=10, sweeps=10, verbose=false)
    push!(mc, :conf => ConfigurationMeasurement, :TH)

    obs = observables(mc)
    save_measurements!(mc, "test.jld", force_overwrite=true)
    _obs = load_measurements("test.jld")
    @test obs == _obs
    rm("test.jld")
end


text:
include("../src/Groebner.jl")

if !isdefined(Main, :Groebner)
    import Groebner
end

import AbstractAlgebra
using BenchmarkTools
using Logging
global_logger(ConsoleLogger(stderr, Logging.Error))

# BenchmarkTools.DEFAULT_PARAMETERS.seconds = 6
BenchmarkTools.DEFAULT_PARAMETERS.samples = 4

function benchmark_system_my(system)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]])

    # gb = Groebner.groebner(system)
    # println("length = $(length(gb))")
    # println("degree = $(maximum(AbstractAlgebra.total_degree, gb))")

    @btime Groebner.groebner($system, reduced=false)

end

function run_f4_ff_degrevlex_benchmarks(ground)
    systems = [
    ("cyclic 7", reverse(Groebner.cyclicn(7, ground=ground))),
    ("cyclic 8", reverse(Groebner.cyclicn(8, ground=ground))),
    ("root 12", reverse(Groebner.rootn(12, ground=ground))),
    ("root 13", reverse(Groebner.rootn(13, ground=ground))),
    #("root 14", Groebner.rootn(14, ground=ground)),
        ("katsura 9", reverse(Groebner.katsuran(9, ground=ground))),
        ("noon 7"    ,Groebner.noonn(7, ground=ground)),
        ("noon 8"    ,Groebner.noonn(8, ground=ground)),
        ("eco 10"    ,Groebner.eco10(ground=ground))
    ]

    for (name, system) in systems
        println("$name")
        benchmark_system_my(system)
    end
end

println()
ground = AbstractAlgebra.GF(2^31 - 1)
run_f4_ff_degrevlex_benchmarks(ground)

#=
cyclic 7
  209.316 ms (229506 allocations: 79.33 MiB)
cyclic 8
  4.351 s (1362066 allocations: 455.37 MiB)
root 12
  61.705 ms (179930 allocations: 32.39 MiB)
root 13
  211.091 ms (451378 allocations: 67.86 MiB)
root 14
  809.867 ms (1144701 allocations: 176.53 MiB)
katsura 9
  801.246 ms (149540 allocations: 59.57 MiB)
noon 7
  214.115 ms (419411 allocations: 91.51 MiB)
noon 8
  1.842 s (2074539 allocations: 492.01 MiB)
eco 10
  179.854 ms (137529 allocations: 41.41 MiB)
=#


text:
using DataStructures

function letterfreq(file::AbstractString; fltr::Function=(_) -> true)
    sort(Dict(counter(filter(fltr, read(file, String)))))
end

display(letterfreq("src/Letter_frequency.jl"; fltr=isletter))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("array_utils.jl")
include("fdiff_utils.jl")
include("poly_utils.jl")
include("misc_utils.jl")


text:
function dropdown(::Bulma, options::Observable; multiple = false, className = "", kwargs...)
    className = InteractBase.mergeclasses(className, multiple ? "select is-multiple" : "select")
    div_select = Node(:div, className=className)
    dropdown(NativeHTML(), options; multiple = multiple, div_select = div_select, kwargs...)
end


text:
export AUC, gini#, CrossEntropy

using DataFrames: groupby, combine, DataFrame, sort!
using CategoricalArrays: CategoricalVector

# import MLJBase: CrossEntropy
#
# CrossEntropy(x, y::CategoricalVector) = CrossEntropy(x, y.refs .- 1)

function AUC_plot_data(score, target::CategoricalVector;  kwargs...)
    @assert length(levels(target)) == 2
    _AUC_plot_data(score, 2 .- target.refs; kwargs...)
end

AUC_plot_data(score, target;  kwargs...) = _AUC_plot_data(score, target;  kwargs...)

function _AUC_plot_data(score, target;  plotauc = false)
    grp_df = groupby(DataFrame(score = score, target = target), :score)
    tmpdf = combine(grp_df,
        :target => sum => :ntarget,
        nrow => :n
    )
    sort!(tmpdf,:score)
    nrows = length(score)
    cutarget = accumulate(+, tmpdf[!,:ntarget]) ./ sum(tmpdf[!,:ntarget])
    cu = accumulate(+, tmpdf[!,:n]) ./ sum(tmpdf[!,:n])

    if plotauc
        plot(vcat(0,cu), vcat(0,cutarget))
        plot!([0,1],[0,1])
    end

    cu, cutarget
end

"""
    AUC(score, target; plotauc = false)

Return the AUC. To generate a plot set `plotauc=true`.
"""
function AUC(score, target; kwargs...)
    cu, cutarget = AUC_plot_data(score, target; kwargs...)
    sum((cutarget[2:end] .+ cutarget[1:end-1]).*(cu[2:end].-cu[1:end-1])./2)
end

"""
    gini(score, target; plotauc = false)

Return the `gini = (AUC - 0.5)/0.05 = 2AUC - 1`. AUC is the area under the curve while gini
is the ratio of (AUC minus the area of the bottom triangle) vs (Area of upper triangle).

For AUC a random model has AUC = 0.5 but for gini a random model has gini = 0.0

To generate a plot set `plotauc=true`.
"""
function gini(score, target; plotauc = false)
    auc = AUC(score,target; plotauc = plotauc)
    2*auc-1
end


text:
# Author: Jeremy Werner

using MultivariateStats

euclidean_distance(vector1::Array, vector2::Array) = √(sumabs2(vector1 - vector2))

function voting_function(vector_in::Array)
	counts_map = countmap(vector_in)
	k = [item[1] for item in counts_map]
	v = [item[2] for item in counts_map]
	indicators = sortperm(v, rev = true)
	if(length(indicators) > 1)
		if(k[indicators[1]] == k[indicators[2]])
			warn("Voting Tie! Please select k that is not a multiple of the number of output classes")
		end
	end
	k[indicators[1]]
end

# Takes matrix (with rows as observations) and converts it to an array of observations as arrays
function data_converter(input::Matrix)
	(rows,cols) = size(input)
	output = Array(Any, rows)
	for k in 1:rows
		output[k] = input[k,:]
	end
	output
end

# Returns array of categorical predictions for the predictors input
function kNN(X::Matrix, y::Array, k::Integer, predictors::Matrix)
	N = size(predictors, 1)
	predictions = Array(ASCIIString, N)
	train_size = size(X,1)
	distances = Array(Float64, train_size)

	if ndims(X) > 1
		if size(X,2) > 2
		 	pca_model = fit(PCA, transpose(X))
		 	X_array = data_converter(transpose(transform(pca_model, transpose(X))))
		 	predict_array = data_converter(transpose(transform(pca_model, transpose(predictors))))
		else
			X_array = data_converter(X)
			predict_array = data_converter(predictors)
		end
	else
		X_array = data_converter(X)
		predict_array = data_converter(predictors)
	end

	for i = 1:N
		# Use a distance metric and voting scheme of your choice here (as defined above)
		for j in 1:train_size
			distances[j] = euclidean_distance(X_array[j], predict_array[i])
		end
		indexes = sortperm(distances)
		predictions[i] = voting_function(y[indexes[1:k]])
	end
	predictions
end # kNN


text:
@testset "unused file" begin
    @test true
 end


text:
baremodule PathnameWrapper

using ...Pathnames: pathname
import Base

pwd() = pathname(Base.pwd())

readdir(dir::AbstractString=Base.pwd(); join::Bool=false, sort::Bool=true) = pathname.(Base.readdir(dir, join=join, sort=sort))

tempname(parent = Base.tempdir(); cleanup=true) = pathname(Base.tempname(parent, cleanup=cleanup))

tempdir() = pathname(Base.tempdir())

mktempdir(parent=Base.tempdir(); prefix="jl_", cleanup=true) = pathname(Base.mktempdir(parent, prefix=prefix, cleanup=cleanup))

homedir() = pathname(Base.homedir())

end # Wrappers


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

struct RUsage
    ru_utime_sec::Clong         #  user CPU time used
    ru_utime_usec::Clong        #  user CPU time used
    ru_stime_sec::Clong         #  system CPU time used
    ru_stime_usec::Clong        #  system CPU time used
    ru_maxrss::Clong            #  maximum resident set size
    ru_ixrss::Clong             #  integral shared memory size
    ru_idrss::Clong             #  integral unshared data size
    ru_isrss::Clong             #  integral unshared stack size
    ru_minflt::Clong            #  page reclaims (soft page faults)
    ru_majflt::Clong            #  page faults (hard page faults)
    ru_nswap::Clong             #  swaps
    ru_inblock::Clong           #  block input operations
    ru_oublock::Clong           #  block output operations
    ru_msgsnd::Clong            #  IPC messages sent
    ru_msgrcv::Clong            #  IPC messages received
    ru_nsignals::Clong          #  signals received
    ru_nvcsw::Clong             #  voluntary context switches
    ru_nivcsw::Clong            #  involuntary context switches
end

function get_vmsize()
    ru = Vector{RUsage}(uninitialized, 1)
    ccall(:getrusage, Cint, (Cint, Ptr{Cvoid}), 0, ru)
    return ru[1].ru_maxrss
end

function run_mtest(name, testf)
    print("Testing $name...")
    for i in 1:2
        print("priming process...")
        testf()
    end
    vm1 = get_vmsize()
    println("monitored run...")
    testf()
    vm2 = get_vmsize()

    diff = vm2 - vm1
    WARN = (diff > 1000) ? "<===================================== WARNING" : ""
    println("Memory Test ($name) : VMSize difference : $diff KB $WARN")
end

function mtest_create_strings()
    for i in 1:10^8
        string("$i")
    end
    GC.gc()
end

function mtest_remotecall_fetch()
    for i in 1:10^5
        remotecall_fetch(myid, 1)
    end
    GC.gc()
end

run_mtest("create_strings", () -> mtest_create_strings())
run_mtest("remotecall_fetch", () -> mtest_remotecall_fetch())

max_stars_repo_name: blegat/julia

max_issues_repo_name: blegat/julia

max_forks_repo_name: blegat/julia


text:
using ParamLevelSet
using jInv.Mesh

println("Testing the rotation itself:")

n = [64,64,64];
Mesh = getRegularMesh([0.0;3.0;0.0;3.0;0.0;3.0],n);
alpha = [1.5;2.5;-2.0;-1.0];
beta = [2.5;2.0;-1.5;2.5];
Xs = [1.5 1.5 1.5 ; 2.0 2.0 2.0; 1.2 2.3 1.5; 2.2 1.5 2.0];
m = wrapTheta(alpha,beta,Xs);

u, = ParamLevelSetModelFunc(Mesh,m;computeJacobian=0);

theta_phi = deg2rad.([37.0 24.0 ; 27.0 14.0]);
b = [5.0 5.0 7.0 ; -2.3 -5.2 3.5]*(sum(Mesh.h)/3.0);

Xc = convert(Array{Float32,2},getCellCenteredGrid(Mesh));

N = prod(size(u));
ur1 = zeros(N ,2);
ur2 = copy(ur1);

u = reshape(u,tuple(n...));

@elapsed ur1[:,1] = reshape(rotateAndMove3D(u,theta_phi[1,:],b[1,:]./Mesh.h)[1],N,1);
@elapsed ur1[:,2] = reshape(rotateAndMove3D(u,theta_phi[2,:],b[2,:]./Mesh.h)[1],N,1);

mr2,Jrot = rotateAndMoveRBFsimple(m,Mesh,theta_phi,b;computeJacobian = 1);
@elapsed ur2[:,1] .= ParamLevelSetModelFunc(Mesh,mr2[:,1];computeJacobian=0,Xc = Xc)[1];
@elapsed ur2[:,2] .= ParamLevelSetModelFunc(Mesh,mr2[:,2];computeJacobian=0,Xc = Xc)[1];

err = norm(ur1-ur2)./(norm(ur1));
if err < 1.5e-1
	println("RBF Rotation test succeeded: ",err);
else
	error("RBF Rotation test not succeeded: ",err);
end

println("Testing the sensitivity of the rotation (simple version - fixed theta/phi and b):");
println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");

m0 = m;
dm = 0.01*randn(length(m));

hh = 1.0;
for k=0:5
	hhh = (0.5^k)*hh;
	mt = rotateAndMoveRBFsimple(m + hhh*dm,Mesh,theta_phi,b;computeJacobian = 0)[1];
	println("norm(mt-mr2): ",norm(mt[:]-mr2[:]),", norm(mt - mr2 - Jrot*dm): ",norm(mt[:] - mr2[:] - Jrot*hhh*dm));
end

println("Testing the sensitivity of the rotation (with varying theta/phi and b):");
println("~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~");
nRBF = div(length(m),5);
nMoves = size(theta_phi,1);
m_wrapped = wrapRBFparamAndRotationsTranslations(m,theta_phi,b);
m_t,theta_phi_t,b_t = splitRBFparamAndRotationsTranslations(m_wrapped,nRBF,nMoves);

if norm(m_t - m) > 1e-5 || norm(theta_phi_t - theta_phi)>1e-5 || norm(b_t - b)>1e-5
	error("wrap/split RBFparamAndRotationsTranslations failed");
end

m0 = wrapRBFparamAndRotationsTranslations(m,theta_phi,b);
dm = 0.01*randn(length(m0));

mr2,Jrot = rotateAndMoveRBF(m,Mesh,theta_phi,b;computeJacobian = 1);

hh = 1.0;
for k=0:10
	hhh = (0.5^k)*hh;
	(mt,theta_phi_t,b_t) = splitRBFparamAndRotationsTranslations(m0+hhh*dm,nRBF,nMoves);
	mt2 = rotateAndMoveRBF(mt,Mesh,theta_phi_t,b_t;computeJacobian = 0)[1];
	println("norm(mt2-mr2): ",norm(mt2[:]-mr2[:]),", norm(mt2 - mr2 - Jrot*dm): ",norm(mt2[:] - mr2[:] - Jrot*hhh*dm));
end


text:
env = ENV["AHORN_ENV"]

redirect_stderr(stdout)

using Pkg
Pkg.activate(env)

ahornPath = Base.find_package("Ahorn")

infofilePath = joinpath(dirname(env), "info.txt")
open(infofilePath, "a") do info
    println(info, "")
    println(info, "Additional info from prepare-for-redistribution.jl:")
    println(info, "Ahorn install path: $(something(ahornPath, "?"))")
    try
        local ctx = Pkg.Types.Context()
        println(info, "Ahorn version: $(string(ctx.env.manifest[ctx.env.project.deps["Ahorn"]].tree_hash)[1:7])")
    catch e
        println(info, "Ahorn version: unknown")
    end
end

if ahornPath === nothing then
    return
end

depot = joinpath(dirname(env), "julia-depot")

rm(joinpath(depot, "clones"), force = true, recursive = true)
rm(joinpath(depot, "compiled"), force = true, recursive = true)

packages = joinpath(depot, "packages")
rm(joinpath(packages, "Maple"), force = true, recursive = true)
rm(joinpath(packages, "Ahorn"), force = true, recursive = true)

max_stars_repo_name: DemoJameson/Ahorn-VHD

max_issues_repo_name: DemoJameson/Ahorn-VHD

max_forks_repo_name: DemoJameson/Ahorn-VHD


text:
"""
    CountEdgesTask

This module includes the composite type CountEdgesTask which includes both the
        generic DaemonTask.Info and an AbstractString as a payload
"""
module CountEdgesTask

using JSON
using PyCall
@pyimport neuroglancer.pipeline as pl

using ...SimpleTasks.Types
using ContactAnalysis
using ContactAnalysis.Precomputed.PrecomputedWrapper

import SimpleTasks.Tasks.DaemonTask
import SimpleTasks.Tasks.BasicTask
import SimpleTasks.Services.Datasource

export CountEdgesTaskDetails, NAME, execute

function str_to_slices(s)
    m=match(r"(\d+)-(\d+)_(\d+)-(\d+)_(\d+)-(\d+)",s)
    bounds = map(x->parse(Int,x), m.captures)
    return (bounds[1]:bounds[2], bounds[3]:bounds[4], bounds[5]:bounds[6])
end

function slices_to_str(s)
    t=map(slice_to_str,s)
    return "$(t[1])_$(t[2])_$(t[3])"
end

function slice_to_str(s)
    return "$(s.start)-$(s.stop)"
end

type CountEdgesPayloadInfo
    segmentation_storage::AbstractString
    slices::AbstractString
    scale_idx::Int64
    output_storage::AbstractString
end

function CountEdgesPayloadInfo(d::Dict{String, Any})
    return CountEdgesPayloadInfo(d["segmentation_storage"],
                                 d["slices"],
                                 d["scale_idx"],
                                 d["output_storage"])
end

type CountEdgesTaskDetails <: DaemonTaskDetails
    basic_info::BasicTask.Info
    payload_info::CountEdgesPayloadInfo
end

"""
Parse task from JSON
"""
function CountEdgesTaskDetails(basic_info::BasicTask.Info, d::Dict{String, Any})
    payload_info = CountEdgesPayloadInfo(d)
    return CountEdgesTaskDetails(basic_info, payload_info)
end

const NAME = "CountEdgesTask"

function DaemonTask.prepare(task::CountEdgesTaskDetails,
                                datasource::DatasourceService)
    return true
end

function DaemonTask.execute(task::CountEdgesTaskDetails,
                                datasource::DatasourceService)

    if length(task.payload_info.segmentation_storage) == ""
        return DaemonTask.Result(true, [])
    end

    segmentation_storage = task.payload_info.segmentation_storage
    str_slices = task.payload_info.slices
    slices = str_to_slices(str_slices)
    scale_idx = task.payload_info.scale_idx
    output_storage = task.payload_info.output_storage

    seg = PrecomputedWrapper(segmentation_storage, scale_idx)
    res = seg.val[:_scale]["resolution"]/1000 # convert nm to um
    weights = []
    for (i,r) in enumerate(res) # calculate surface areas of faces
        push!(weights, prod(res[filter(a->!(i in a), 1:length(res))]))
    end
    d = ContactAnalysis.count_edges(seg[slices...], weights)

    out = pl.Storage(output_storage)
    out[:put_file](file_path="$(str_slices)_edge_counts.json",
                        content="$(JSON.json(d))")
    out[:wait]()

    return DaemonTask.Result(true, [str_slices])
end

function DaemonTask.finalize(task::CountEdgesTaskDetails,
                datasource::DatasourceService, result::DaemonTask.Result)
    return true
end

end # module CountEdgesTask


text:
function writePolyMeshPoints(nodes)
    open("src/testCases/writeMesh/constant/polyMesh/points", "w") do f
    write(f, "FoamFile\n")
    write(f, "{\n")
    write(f, "    version     2.0;\n")
    write(f, "    format      ascii;\n")
    write(f, "    arch        \"LSB;label=32;scalar=64\";\n")
    write(f, "    class       vectorField;\n")
    write(f, "    location    \"constant/polyMesh\";\n")
    write(f, "    object      points;\n")
    write(f, "}\n")
    write(f, "$(length(nodes))\n")
    write(f, "(\n")
    for node in nodes
        write(f, "($(node.centroid[1]) $(node.centroid[2]) $(node.centroid[3]))\n")
    end
    write(f, ")\n")
    end
end

function writePolyMeshFaces(faces)
    open("src/testCases/writeMesh/constant/polyMesh/faces", "w") do f
    write(f, "FoamFile\n")
    write(f, "{\n")
    write(f, "    version     2.0;\n")
    write(f, "    format      ascii;\n")
    write(f, "    arch        \"LSB;label=32;scalar=64\";\n")
    write(f, "    class       faceList;\n")
    write(f, "    location    \"constant/polyMesh\";\n")
    write(f, "    object      faces;\n")
    write(f, "}\n")
    write(f, "$(length(faces))\n")
    write(f, "(\n")
    for face in faces
        write(f, "$(length(face.iNodes))")
        write(f, "(")
        for iNode in face.iNodes
            write(f, " $(iNode - 1)")
        end
        write(f, ")\n")
    end
    write(f, ")\n")
    end
end

function writePolyMeshOwner(faces)
    open("src/testCases/writeMesh/constant/polyMesh/owner", "w") do f
    write(f, "FoamFile\n")
    write(f, "{\n")
    write(f, "    version     2.0;\n")
    write(f, "    format      ascii;\n")
    write(f, "    arch        \"LSB;label=32;scalar=64\";\n")
    write(f, "    class       labelList;\n")
    write(f, "    location    \"constant/polyMesh\";\n")
    write(f, "    object      owner;\n")
    write(f, "}\n")
    write(f, "$(length(faces))\n")
    write(f, "(\n")
    for face in faces
        write(f, "$(face.iOwner - 1)\n")
    end
    write(f, ")\n")
    end
end

function writePolyMeshNeighbour(faces, boundaries)
    open("src/testCases/writeMesh/constant/polyMesh/neighbour", "w") do f
    write(f, "FoamFile\n")
    write(f, "{\n")
    write(f, "    version     2.0;\n")
    write(f, "    format      ascii;\n")
    write(f, "    arch        \"LSB;label=32;scalar=64\";\n")
    write(f, "    class       labelList;\n")
    write(f, "    location    \"constant/polyMesh\";\n")
    write(f, "    object      neighbour;\n")
    write(f, "}\n")
    write(f, "$(boundaries[1].startFace - 1)\n")
    write(f, "(\n")
    for iFace in 1:boundaries[1].startFace - 1
        write(f, "$(faces[iFace].iNeighbour - 1)\n")
    end
    write(f, ")\n")
    end
end

function writePolyMeshBoundaries(boundaries)
    open("src/testCases/writeMesh/constant/polyMesh/boundary", "w") do f
    write(f, "FoamFile\n")
    write(f, "{\n")
    write(f, "    version     2.0;\n")
    write(f, "    format      ascii;\n")
    write(f, "    arch        \"LSB;label=32;scalar=64\";\n")
    write(f, "    class       polyBoundaryMesh;\n")
    write(f, "    location    \"constant/polyMesh\";\n")
    write(f, "    object      boundary;\n")
    write(f, "}\n")
    write(f, "$(length(boundaries))\n")
    write(f, "(\n")
    for boundary in boundaries
        write(f, "$(boundary.name)\n")
        write(f, "{\n")
        write(f, "type       $(boundary.type);\n")
        write(f, "nFaces     $(boundary.nFaces);\n")
        write(f, "startFace  $(boundary.startFace - 1);\n")
        write(f, "}\n")
    end
    write(f, ")\n")
    end
end

function writeOpenFoamMesh(mesh)
    writePolyMeshPoints(mesh.nodes)
    writePolyMeshFaces(mesh.faces)
    writePolyMeshOwner(mesh.faces)
    writePolyMeshNeighbour(mesh.faces, mesh.boundaries)
    writePolyMeshBoundaries(mesh.boundaries)
end


text:
"""
    encode(::Type{SimpleEncoder}, s::String)

Encode the input `String` object as `SimpleEncoder`.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> basmala = arabic(verses(crpsdata[1][1])[1])
julia> encode(SimpleEncoder, basmala)
"Ba+Kasra | Seen+Sukun | Meem+Kasra | <space> | HamzatWasl | Lam | Lam+Shadda+Fatha | Ha+Kasra | <space> | HamzatWasl | Lam | Ra+Shadda+Fatha | HHa+Sukun | Meem+Fatha | AlifKhanjareeya | Noon+Kasra | <space> | HamzatWasl | Lam | Ra+Shadda+Fatha | HHa+Kasra | Ya | Meem+Kasra"
```
"""
function encode(::Type{SimpleEncoder}, s::String)
    sencoder = SimpleEncoder()
    words = ""; i = 1
    for c in s
        if c === ' '
            words *= " | <space>"
        else
            words *= encode(sencoder, c, i)
        end
        i += 1
    end
    return words
end

function encode(encoder::SimpleEncoder, c::Union{Char,String}, i::Int64)
    if in(Symbol(c), AR_DIACS)
        return string("+", encoder.encode[Symbol(c)])
    else
        if i > 1
            return string(" | ", encoder.encode[Symbol(c)])
        else
            return string(encoder.encode[Symbol(c)])
        end
    end
end

"""
    encode(s::Union{Char,String}, encoder::AbstractEncoder)

Transliterate the input `String` object using a custom `encoder`. Custom `encoder` is
generated using the `@transliterator`.
"""
function encode(s::Union{Char,String}, encoder::AbstractEncoder)
    words = ""
    for c in s
        if c === ' '
            words *= " "
        else
            if c === '\U0622'
                words *= string(encoder.encode[Symbol('\U0627')], encoder.encode[Symbol('\U0653')])
            else
                words *= string(encoder.encode[Symbol(c)])
            end
        end
    end
    return words
end

"""
    encode(s::String)

Transliterate the input `String` object using `Buckwalter`.

# Examples
```julia-repl
julia> data = QuranData()
julia> crps, tnzl = load(data)
julia> crpsdata = table(crps)
julia> tnzldata = table(tnzl)
julia> basmala = arabic(verses(crpsdata[1][1])[1])
julia> encode(basmala)
"bisomi {ll~ahi {lr~aHoma`ni {lr~aHiymi"
```
"""
function encode(s::String)
    trans = Transliterator()
    return encode(s, trans)
end


text:
using QuantLib
using Printf
using Dates

function print_calibration_basket(basket::Vector{S}) where {S <: SwaptionHelper}
  println("Expiry         Maturity       Nominal        Rate               Pay/Rec    Market ivol")
  println("======================================================================================")
  for i in eachindex(basket)
    helper = basket[i]
    underlying = underlying_swap!(helper)
    endDate = underlying.fixedSchedule.dates[end]
    nominal = underlying.nominal
    vol = helper.volatility.value
    rate = underlying.fixedRate
    expiry = helper.swaption.exercise.dates[1]
    swapT = underlying.swapT == QuantLib.Receiver() ? "Receiver" : "Payer"
    println("$expiry    $endDate      $nominal       $rate   $swapT      $vol")
  end
end

function print_model_calibration(basket::Vector{S}, volatilities::Vector{Float64}) where {S <: SwaptionHelper}
  println("Expiry   Model Sigma  ModelPx  MktPx     ModelVol  MktVol")
  println("=========================================================")
  for i in eachindex(basket)
    helper = basket[i]
    expiry = helper.swaption.exercise.dates[1]
    modelSig = volatilities[i]
    modelPx = model_value!(helper)
    mktPx = helper.calibCommon.marketValue
    modelVol = implied_volatility!(helper, modelPx, 1e-6, 1000, 0.0, 2.0)
    marketVol = helper.volatility.value

    println(@sprintf("%s  %.6f  %.6f  %.6f  %.2f     %.2f", expiry, modelSig, modelPx, mktPx, modelVol, marketVol))
  end
end

function main()
  refDate = Dates.Date(2014, 4, 30)
  set_eval_date!(settings, refDate)

  forward6mLevel = 0.025
  oisLevel = 0.02

  forward6mQuote = Quote(forward6mLevel)
  oisQuote = Quote(oisLevel)

  yts6m = FlatForwardTermStructure(0, QuantLib.Time.TargetCalendar(), forward6mQuote, QuantLib.Time.Actual365())
  ytsOis = FlatForwardTermStructure(0, QuantLib.Time.TargetCalendar(), oisQuote, QuantLib.Time.Actual365())

  euribor6m = euribor_index(QuantLib.Time.TenorPeriod(Dates.Month(6)), yts6m)

  volLevel = 0.2
  volQuote = Quote(volLevel)

  swaptionVol = ConstantSwaptionVolatility(0, QuantLib.Time.TargetCalendar(), QuantLib.Time.ModifiedFollowing(), volQuote, QuantLib.Time.Actual365())

  strike = 0.04

  effectiveDate = QuantLib.Time.advance(Dates.Day(2), QuantLib.Time.TargetCalendar(), refDate)
  maturityDate = QuantLib.Time.advance(Dates.Year(10), QuantLib.Time.TargetCalendar(), effectiveDate)

  schedBizConvention = QuantLib.Time.ModifiedFollowing()

  fixedSchedule = QuantLib.Time.Schedule(effectiveDate, maturityDate, QuantLib.Time.TenorPeriod(Dates.Year(1)), schedBizConvention, schedBizConvention,
                  QuantLib.Time.DateGenerationForwards(), false, QuantLib.Time.TargetCalendar())
  floatingSchedule = QuantLib.Time.Schedule(effectiveDate, maturityDate, QuantLib.Time.TenorPeriod(Dates.Month(6)), schedBizConvention, schedBizConvention,
                  QuantLib.Time.DateGenerationForwards(), false, QuantLib.Time.TargetCalendar())

  exerciseDates = Vector{Date}(undef, 9)

  for i = 2:10
    exerciseDates[i - 1] = QuantLib.Time.advance(Dates.Day(-2), QuantLib.Time.TargetCalendar(), fixedSchedule.dates[i])
  end

  stepDates = exerciseDates[1:end - 1]
  sigmas = fill(0.01, length(exerciseDates))
  reversion = 0.01

  gsr = GSR(yts6m, stepDates, sigmas, reversion)

  swaptionEngine = Gaussian1DSwaptionEngine(gsr, 64, 7.0, true, false, ytsOis)
  nonstandardSwaptionEngine = Gaussian1DNonstandardSwaptionEngine(gsr, 64, 7.0, true, false, Quote(-1.0), ytsOis)

  underlying = NonstandardSwap(VanillaSwap(Payer(), 1.0, fixedSchedule, strike, QuantLib.Time.BondThirty360(), euribor6m, 0.0, floatingSchedule, QuantLib.Time.Actual360()))

  exercise = BermudanExercise(exerciseDates)

  swaption = NonstandardSwaption(underlying, exercise, nonstandardSwaptionEngine)

  swapBase = EuriborSwapIsdaFixA(QuantLib.Time.TenorPeriod(Dates.Year(10)), yts6m, ytsOis)

  basket = calibration_basket(swaption, nonstandardSwaptionEngine, swapBase, swaptionVol, NaiveBasketType())

  print_calibration_basket(basket)

  for i in eachindex(basket)
    basket[i] = update_pricing_engine(basket[i], swaptionEngine)
  end

  method = QuantLib.Math.LevenbergMarquardt()

  ec = QuantLib.Math.EndCriteria(1000, 10, 1.0e-8, 1.0e-8, 1.0e-8)
  calibrate_volatilities_iterative!(gsr, basket, method, ec)

  print_model_calibration(basket, get_volatilities(gsr))

  println(npv(swaption))
end


text:
const MASK_RANKS = (0x00000000000000ff,
                    0x000000000000ff00,
                    0x0000000000ff0000,
                    0x00000000ff000000,
                    0x000000ff00000000,
                    0x0000ff0000000000,
                    0x00ff000000000000,
                    0xff00000000000000)

const INSIDE = reduce(|, MASK_RANKS)

const MASK_FILES = (0x8080808080808080,
                    0x4040404040404040,
                    0x2020202020202020,
                    0x1010101010101010,
                    0x0808080808080808,
                    0x0404040404040404,
                    0x0202020202020202,
                    0x0101010101010101)

const CLEAR_FILE_A = 0x7f7f7f7f7f7f7f7f
const CLEAR_FILE_H = 0xfefefefefefefefe


text:
using Flux
using Flux.Data: DataLoader
using Flux.Optimise: Optimiser, WeightDecay
using Flux: onehotbatch, onecold, normalise, Dropout
using Flux.Losses: logitcrossentropy
using Statistics, Random
using Logging: with_logger
using TensorBoardLogger: TBLogger, tb_overwrite, set_step!, set_step_increment!
using ProgressMeter: @showprogress
import MLDatasets
import BSON
using CUDA
using Glob
using SignalAnalysis
using WAV
using Images
using Random:shuffle
println("Import Done")

labeldict = Dict("left" => 0, "right" => 1, "up" => 2,  "down" => 3)

function ConvNet5(; imgsize=(28,28,1), nclasses=4)
    return Chain(
            		Flux.normalise,
			Conv((3, 3), imgsize[end]=>80 ,elu),
			MaxPool((2, 2)),
			Dropout(0.25),
		    	flatten,
			Dense(13520, 128, relu),
			Dropout(0.5),
			Dense(128, nclasses)
          )
end

function read_audio(filepath)
	val, fs = wavread(filepath)
	lengthneed= 16000 - length(val)
	created = zeros(Float32, lengthneed)
	result = vcat(val[:],created)
	spec= spectrogram(result, 255, 128; fs=fs);
        y = imresize(spec.power, (28, 28));
        y = reshape(convert(Array{Float32},y ), (28,28,1));
	return y
end
function get_label(filepath)
	label = split(filepath, "/")[end-1]
	return Float32.(labeldict[label])
end
function loaddata(paths)
	data = zeros(Float32, ( 28,28, 1, length(paths)))
	for (num ,i) in enumerate(paths)
		data[:,:,:, num] = read_audio(i)
	end
	return data
end

function get_data(args)
    DATASET_PATH = "./mini_speech_commands"
    left = glob("left/*.wav", "./mini_speech_commands")
    right = glob("right/*.wav", "./mini_speech_commands")
    up = glob("up/*.wav", "./mini_speech_commands")
    down = glob("down/*.wav", "./mini_speech_commands")
    datalist = vcat(left, right, up, down)
    filenames = shuffle(datalist)
    train_files = filenames[1:2500]
    val_files = filenames[2500:3300]
    test_files = filenames[3300:4000]

    xtrain = loaddata(train_files)
    ytrain =  onehotbatch(map(x->get_label(x), train_files), 0:3)

    xtest = loaddata(val_files)
    ytest =   onehotbatch(map(x->get_label(x), val_files), 0:3)

    train_loader = DataLoader((xtrain, ytrain), batchsize=args.batchsize, shuffle=true)
    test_loader = DataLoader((xtest, ytest),  batchsize=args.batchsize, shuffle=true)

    return train_loader, test_loader
end

loss(ŷ, y) = logitcrossentropy(ŷ, y)

function eval_loss_accuracy(loader, model, device)
    l = 0f0
    acc = 0
    ntot = 0
    for (x, y) in loader
        x, y = x |> device, y |> device
        ŷ = model(x)
        l += loss(ŷ, y) * size(x)[end]
        acc += sum(onecold(ŷ |> cpu) .== onecold(y |> cpu))
        ntot += size(x)[end]
    end
    return (loss = l/ntot |> round4, acc = acc/ntot*100 |> round4)
end

println("Functions Done")

## utility functions
num_params(model) = sum(length, Flux.params(model))
round4(x) = round(x, digits=4)

# arguments for the `train` function
Base.@kwdef mutable struct Args
    η = 3e-4             # learning rate
    λ = 0                # L2 regularizer param, implemented as weight decay
    batchsize = 128      # batch size
    epochs = 50          # number of epochs
    seed = 0             # set seed > 0 for reproducibility
    use_cuda = true      # if true use cuda (if available)
    infotime = 1 	     # report every `infotime` epochs
    checktime = 5        # Save the model every `checktime` epochs. Set to 0 for no checkpoints.
    tblogger = true      # log training with tensorboard
    savepath = "runs/"    # results path
end

function train(; kws...)
    args = Args(; kws...)
    args.seed > 0 && Random.seed!(args.seed)
    use_cuda = args.use_cuda && CUDA.functional()

    if use_cuda
        device = gpu
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    ## DATA
    train_loader, test_loader = get_data(args)
    @info "Dataset Speech Commands: $(train_loader.nobs) train and $(test_loader.nobs) test examples"

    ## MODEL AND OPTIMIZER
    model = ConvNet5() |> device
    @info "ConvNet5 model: $(num_params(model)) trainable params"

    ps = Flux.params(model)

    opt = ADAM(args.η)
    if args.λ > 0 # add weight decay, equivalent to L2 regularization
        opt = Optimiser(WeightDecay(args.λ), opt)
    end

    ## LOGGING UTILITIES
    if args.tblogger
        tblogger = TBLogger(args.savepath, tb_overwrite)
        set_step_increment!(tblogger, 0) # 0 auto increment since we manually set_step!
        @info "TensorBoard logging at \"$(args.savepath)\""
    end

    function report(epoch)
        train = eval_loss_accuracy(train_loader, model, device)
        test = eval_loss_accuracy(test_loader, model, device)
        println("Epoch: $epoch   Train: $(train)   Test: $(test)")
        if args.tblogger
            set_step!(tblogger, epoch)
            with_logger(tblogger) do
                @info "train" loss=train.loss  acc=train.acc
                @info "test"  loss=test.loss   acc=test.acc
            end
        end
    end

    ## TRAINING
    @info "Start Training"
    report(0)
    for epoch in 1:args.epochs
        @showprogress for (x, y) in train_loader
            x, y = x |> device, y |> device
            gs = Flux.gradient(ps) do
                    ŷ = model(x)
                    loss(ŷ, y)
                end

            Flux.Optimise.update!(opt, ps, gs)
        end

        ## Printing and logging
        epoch % args.infotime == 0 && report(epoch)
        if args.checktime > 0 && epoch % args.checktime == 0
            !ispath(args.savepath) && mkpath(args.savepath)
            modelpath = joinpath(args.savepath, "model_$(epoch/5).bson")
            let model = cpu(model) #return model to cpu before serialization
                BSON.@save modelpath model epoch
            end
            @info "Model saved in \"$(modelpath)\""
        end
    end
end
println("Training Coming")

train()

max_stars_repo_name: ashwani-rathee/Sardor

max_issues_repo_name: ashwani-rathee/Sardor

max_forks_repo_name: ashwani-rathee/Sardor


text:
push!(LOAD_PATH, "../src/")

using Documenter
using HarmonicBalance

makedocs(
	sitename="HarmonicBalance.jl",
	format = Documenter.HTML(mathengine=MathJax()),
	pages = [
		"Background" => Any[
			"background/harmonic_balance.md"
			"background/stability_response.md"
		],
		"Examples" => Any[
			"examples/single_Duffing.md"
			"examples/linear_response.md"
			"examples/time_dependent.md"
			"examples/single_parametron_1D.md"
			"examples/single_parametron_2D.md"
			"examples/single_parametron_time_dep.md"
			],
		"Manual" => Any[
			"manual/entering_eom.md"
			"manual/extracting_harmonics.md"
			"manual/solving_harmonics.md"
			"manual/plotting.md"
			"manual/time_dependent.md"
			"manual/linear_response.md"
			"manual/saving.md"
			]
		]
)

deploydocs(
    repo = "github.com/NonlinearOscillations/HarmonicBalance.jl.git",
    push_preview = false
)


text:
struct SimpleSequentialInhibition{T <: Real}
    minx::T
    maxx::T
    miny::T
    maxy::T
    R::T
end
""" SimpleSequentialInhibition(window, R)

Define a simple sequential inhibition (a.k.a random sequential adsorbtion) point
process in `window` with minimum interpoint distance `R`.

The rand function has keyword argument `T` specifying the maximum number of
uniform proposal points as `T*area(window)/R^2`.

The algorithm used for sampling from the process is based on (Wang, 1994) with a
few modifications. For the skipped uniform proposals only the expected count is
used instead of a random variable.

J.-S. Wang, A fast algorithm for random sequential adsorption of discs,
International Journal of Modern Physics C 1994 05:04, 707-715.

```julia
rand(SimpleSequentialInhibition((x=(0.0, 100.0), y=(0.0, 100.0)), 10.0), T=10.0^20)
```
"""
SimpleSequentialInhibition(window, R) =
    SimpleSequentialInhibition{Float64}(window.x[1], window.x[2], window.y[1], window.y[2], R)
const RandomSequentialAdsorbtion = SimpleSequentialInhibition

# Algorithm 1
function tooclose(p, p0, R)
    x0, y0 = p0
    R2 = Float64(R^2)
    for i1 in eachindex(p)
        x1, y1 = p[i1]
        d2 = (x1-x0)^2 + (y1-y0)^2
        if d2 < R2
            return true
        end
    end
    false
end
function randSSIfillsimple!(p, (maxx, maxy), R, T)
    t = 0.0
    dt = R^2 / (maxx*maxy)
    while t <= T
        t += dt
        p1 = (rand()*maxx, rand()*maxy)
        if !tooclose(p, p1, R)
            push!(p, p1)
        end
    end
    p
end
function randSSIsimple(window, R, T)
    p = randp(window, 1)
    randSSIfillsimple!(p, window, R, T)
end

# Algorithm 2 with a grid based search for nearby points.
function tooclose!(closepoints::Array{<:Integer, 3}, pp, (x0,y0), R)
    xi = unsafe_trunc(Int, x0/R) + 2
    yi = unsafe_trunc(Int, y0/R) + 2
    R2 = Float64(R^2)
    for vyi in yi-1:yi+1
        for vxi in xi-1:xi+1
            ii = 1
            while ii <= 3 && (pi = closepoints[ii, vxi, vyi]) > 0
                @inbounds x1, y1 = pp[pi]
                d2 = (x1-x0)^2 + (y1-y0)^2
                if d2 < R2
                    return true
                end
                ii += 1
            end
        end
    end
    i = 1
    while closepoints[i, xi, yi] > 0
        i += 1
    end
    closepoints[i, xi, yi] = length(pp) + 1
    return false
end
function initializeclosepoints(pp, (maxx,maxy), R)
    closepoints = fill(-1, (3, Int(cld(maxx, R)+2), Int(cld(maxy, R))+2))
    for (pi, (x,y)) in enumerate(pp)
        xi = Int(fld(x, R)) + 2
        yi = Int(fld(y, R)) + 2
        i = 1
        while closepoints[i, xi, yi] > 0
            i += 1
        end
        closepoints[i, xi, yi] = pi
    end
    closepoints
end
# Part of alg 3
function iscovered(x0, x1, y0, y1, discs, closepoints::Array{Int,3}, R, tmp)
    relevantdiscs, tmpintervals = tmp
    xi = unsafe_trunc(Int, (x0 + x1)/(2R)) + 2
    yi = unsafe_trunc(Int, (y0 + y1)/(2R)) + 2
    empty!(relevantdiscs)
    for vyi in yi-1:yi+1
        for vxi in xi-1:xi+1
            ii = 1
            while ii <= 3 && (pi = closepoints[ii, vxi, vyi]) > 0
                x = discs[pi]
                if isrelevant(x0, x1, y0, y1, x, R)
                    push!(relevantdiscs, x)
                end
                ii += 1
            end
        end
    end

    iscovered(x0, x1, y0, y1, relevantdiscs, R, tmpintervals)
end

function randSSIfill!(pp, (maxx, maxy), R, T)
    window = (x=(0, maxx), y=(0, maxy))
    closepoints = initializeclosepoints(pp, (maxx, maxy), R)
    t = 0
    dt = R^2 / (maxx*maxy)
    while t <= T
        t += dt
        p1 = (rand()*maxx, rand()*maxy)

        if !tooclose!(closepoints, pp, p1, R)
            push!(pp, p1)
        end
    end
    pp
end

function randSSI((maxx, maxy), R, T)
    p = [(rand()*maxx, rand()*maxy)]
    randSSIfill!(p, (maxx, maxy), R, T)
end

# Algorithm 3 with occupied cells
function isrelevant(x0, x1, y0, y1, (x, y), R)
    x0 <= x <= x1 && y0-R <= y <= y0+R && return true
    y0 <= y <= y1 && x0-R <= x <= x0+R && return true
    R2 = R^2
    (x0-x)^2 + (y0-y)^2 <= R2 && return true
    (x0-x)^2 + (y1-y)^2 <= R2 && return true
    (x1-x)^2 + (y1-y)^2 <= R2 && return true
    (x1-x)^2 + (y0-y)^2 <= R2 && return true
    return false
end
function iscovered_tmp()
    [Vector{NTuple{2, Float64}}() for _ in 1:4]
end
function iscovered(x0, x1, y0, y1, discs, R)
    iscovered(x0, x1, y0, y1, discs, R, iscovered_tmp())
end
function iscovered(x0, x1, y0, y1, discs, R, tmp)
    intervalss = tmp
    # 1.
    #discs = filter(x -> isrelevant(x0, x1, y0, y1, x, R), discs)
    length(discs) == 0 && return false

    # 2. Covered by a single disc
    R2 = R^2
    for (x,y) in discs
        (x0-x)^2 + (y0-y)^2 <= R2 &&
        (x0-x)^2 + (y1-y)^2 <= R2 &&
        (x1-x)^2 + (y1-y)^2 <= R2 &&
        (x1-x)^2 + (y0-y)^2 <= R2 && return true
    end

    # 3. Not all corners covered
    any(discs) do (x,y)
        (x0-x)^2 + (y0-y)^2 <= R2  end || return false
    any(discs) do (x,y)
        (x0-x)^2 + (y1-y)^2 <= R2  end || return false
    any(discs) do (x,y)
        (x1-x)^2 + (y1-y)^2 <= R2  end || return false
    any(discs) do (x,y)
        (x1-x)^2 + (y0-y)^2 <= R2  end || return false

    # 4.
    # first edge x0,y0 - x1,y0
    # y0 < y1
    sx = x1-x0
    sy = y1-y0
    foreach(empty!, intervalss)

    addinterval(i,a,b) = push!(intervalss[i], (a,b))
    for (x,y) in discs
        dx2 = R^2 - (y-y0)^2
        if dx2 >= 0.0
            dx = sqrt(dx2)
            a0 = x-x0-dx
            a1 = x-x0+dx
            if a0 <= sx && a1 >= 0
                addinterval(1, a0, a1)
            end
        end
        dx2 = R^2 - (y-y1)^2
        if dx2 >= 0.0
            dx = sqrt(dx2)
            a0 = x-x0-dx
            a1 = x-x0+dx
            if a0 <= sx && a1 >= 0
                addinterval(2, a0, a1)
            end
        end

        dx2 = R^2 - (x-x0)^2
        if dx2 >= 0.0
            dx = sqrt(dx2)
            a0 = y-y0-dx
            a1 = y-y0+dx
            if a0 <= sy && a1 >= 0
                addinterval(3, a0, a1)
            end
        end
        dx2 = R^2 - (x-x1)^2
        if dx2 >= 0.0
            dx = sqrt(dx2)
            a0 = y-y0-dx
            a1 = y-y0+dx
            if a0 <= sy && a1 >= 0
                addinterval(4, a0, a1)
            end
        end
    end
    for i in eachindex(intervalss)
        intervals = intervalss[i]
        # Use in-place sorting
        sort!(intervals, alg=Base.Sort.InsertionSort, by=first)
        b = 0.0
        @assert intervals[1][1] <= 0.0
        for (a0, a1) in intervals
            a0 <= b || return false
            b = max(a1, b)
        end
        @assert b >= ifelse(i <= 2, sx, sy)
    end
    length(discs) == 2 && return true

    # 5.
    for i in 1:length(discs)
        for j in i+1:length(discs)
            ax, ay = discs[i]
            bx, by = discs[j]
            d2 = (ax-bx)^2 + (ay-by)^2
            d2 <= 4R2 || continue
            d = sqrt(d2)
            cx = 0.5*(ax + bx)
            cy = 0.5*(ay + by)
            l = sqrt(R^2 - d2/4)
            fx = cx + l*(ay-by)/d
            fy = cy - l*(ax-bx)/d
            gx = cx - l*(ay-by)/d
            gy = cy + l*(ax-bx)/d
            # if f or g is not in the box it doesn't need to be checked
            fcovered = !(x0 <= fx <= x1 && y0 <= fy <= y1)
            gcovered = !(x0 <= gx <= x1 && y0 <= gy <= y1)
            for k in 1:length(discs)
                if k != i && k != j
                    x,y = discs[k]
                    fcovered |= (fx-x)^2 + (fy-y)^2 <= R2
                    gcovered |= (gx-x)^2 + (gy-y)^2 <= R2
                end
            end
            (fcovered && gcovered) || return false
        end
    end
    return true
end

function randSSIfill2!(pp, (maxx, maxy), R, T)
    window = (x=(0, maxx), y=(0, maxy))
    closepoints = initializeclosepoints(pp, (maxx, maxy), R)
    tmp = (Vector{NTuple{2, Float64}}(), iscovered_tmp())
    covered = [iscovered(x, x+R, y, y+R, pp, closepoints, R, tmp) for x in 0:R:maxx, y in 0:R:maxy]
    t = 0.0
    dt = R^2 / (maxx*maxy)
    while t <= T
        t += dt
        (x,y) = (rand()*maxx, rand()*maxy)
        p1 = (x,y)
        xi = unsafe_trunc(Int, x/R) + 1
        yi = unsafe_trunc(Int, y/R) + 1
        if !covered[xi, yi]
            tc = tooclose!(closepoints, pp, p1, R)
            if !tc
                push!(pp, p1)
            end
        end
    end
    pp
end
function build_notcovered(closepoints, pp, (maxx, maxy), R, tmp)
    notcovered = Vector{NTuple{2, Float64}}()
    for x in 0:R:maxx, y in 0:R:maxy
        if x < maxx && y < maxy && !iscovered(x, min(x+R, maxx), y, min(y+R, maxy), pp, closepoints, R, tmp)
            push!(notcovered, (x,y))
        end
    end
    notcovered
end
function rebuild_notcovered!(notcovered, closepoints, pp, (maxx, maxy), R, cellsize, tmp)
    filter!(notcovered) do (x,y)
        x < maxx && y < maxy &&
        !iscovered(x, min(x+cellsize, maxx), y, min(y+cellsize, maxy), pp, closepoints, R, tmp)
    end
    # Here is maybe a possibility of infinite loop if cellsize is too small
    maxs = max(maxx, maxy)
    if maxs == maxs + 0.5cellsize
        throw(("Cellsize got too small", cellsize))
    end
    if  0 < length(notcovered) <= 100
        l = length(notcovered)
        resize!(notcovered, 4*l)
        for i in l:-1:1
            x0, y0 = notcovered[i]
            x1 = x0+0.5cellsize
            y1 = y0+0.5cellsize
            notcovered[4i] = (x1, y1)
            notcovered[4i-1] = (x1, y0)
            notcovered[4i-2] = (x0, y0)
            notcovered[4i-3] = (x0, y1)
        end
        rebuild_notcovered!(notcovered, closepoints, pp, (maxx, maxy), R, 0.5cellsize, tmp)
    else
        cellsize
    end
end
function randSSIfill3!(pp, (maxx, maxy), R, T)
    window = (x=(0, maxx), y=(0, maxy))
    closepoints = initializeclosepoints(pp, (maxx, maxy), R)
    tmp = (Vector{NTuple{2, Float64}}(), iscovered_tmp())
    cellsize = R
    notcovered = build_notcovered(closepoints, pp, (maxx, maxy), cellsize, tmp)

    # proposal_area = (length(notcovered)*cellsize^2)
    # The expected time used before getting a proposal in the not covered squares is
    # total area / proposal_area
    # Normally the time increment is R^2 / total area
    dt = 1/length(notcovered)
    t = 0.0
    proposals = 0
    total = 0
    acc = 0

    isempty(notcovered) && return (status=:jamm, cellsize=cellsize, time=t, totalproposals=total)
    while t <= T
        t += dt
        proposals += 1
        total += 1
        if total > 1_000_000
            println(("Something is wrong", length(pp), length(notcovered), (cellsize, t, dt)))
            throw("asdf")
        end
        if 3acc >= length(notcovered) || proposals >= 10*length(notcovered)
            cellsize = rebuild_notcovered!(notcovered, closepoints, pp, (maxx, maxy), R, cellsize, tmp)
            dt = R^2 /(length(notcovered) * cellsize^2)
            if isempty(notcovered)
                return (status=:jamm, cellsize=cellsize, time=t, totalproposals=total)
            end
            proposals = 0
            acc = 0
        end
        (x0, y0) = rand(notcovered)
        (x,y) = (x0+rand()*cellsize, y0+rand()*cellsize)
        p1 = (x,y)

        if x <= maxx && y <= maxy && !tooclose!(closepoints, pp, p1, R)
            acc += 1
            push!(pp, p1)
        end
    end
    (status=:ok, cellsize=cellsize, time=t, totalproposals=total)
end

function Base.rand(d::SimpleSequentialInhibition; T=10.0^20, alg=:gridrand, info=(info, p, d) -> nothing)
    window = (x=(d.minx, d.maxx), y=(d.miny, d.maxy))
    w = d.maxx-d.minx
    h = d.maxy-d.miny
    status = :ok
    if alg == :simple
        return randSSIsimple(window, d.R, T)
    elseif alg == :gridsearch
        p = randSSI((w, h), d.R, T)
    elseif alg == :gridrand
        # https://arxiv.org/pdf/cond-mat/9402066.pdf
        p = randSSI((w, h), d.R, min(16.0, T))
        p = [(rand()*w, rand()*h)]
        status = randSSIfill3!(p, (w, h), d.R, T)
    else
        throw(ArgumentError("$alg is not a valid algorithm for SimpleSequentialInhibition."))
    end
    for i in eachindex(p)
        x,y = p[i]
        p[i] = (x+d.minx, y+d.miny)
    end
    info(status, p, d)
    PointPattern(p, window)
end


text:
function test_triangle_istriangular()
    @test find([istriangular(i) for i in 1:400]) == [1, 3, 6, 10, 15, 21, 28, 36, 45, 55, 66, 78, 91, 105, 120, 136, 153, 171, 190, 210, 231, 253, 276, 300, 325, 351, 378]
end

function test_triangle_nthtriangular()
    @test nthtriangular(1) == 1
    @test nthtriangular(2) == 3
    @test nthtriangular(10) == 55
end

function test_triangle_ntriangular()
    @test ntriangular(10) == [1, 3, 6, 10, 15, 21, 28, 36, 45, 55]
    @test eltype(ntriangular(10)) == Int
    @test eltype(ntriangular(10, BigInt)) == BigInt
end

function test_triangle_sometriangular()
    @test collect(sometriangular(35)) == [1, 3, 6, 10, 15, 21, 28]
    @test collect(sometriangular(36)) == [1, 3, 6, 10, 15, 21, 28, 36]
    @test eltype(sometriangular(36)) == Int
    @test eltype(sometriangular(big(36))) == BigInt
end

function test_triangle_exacttriangular()
    @test collect(exacttriangular(10)) == [1, 3, 6, 10, 15, 21, 28, 36, 45, 55]
    @test eltype(exacttriangular(10)) == Int
    @test eltype(exacttriangular(10, BigInt)) == BigInt
end

function test_triangle_primitive_pythagorean_triples()
    @test collect(take(primitive_pythagorean_triples(), 13)) == Array{Int,1}[
        [3,4,5],
        [5,12,13],
        [21,20,29],
        [15,8,17],
        [7,24,25],
        [55,48,73],
        [45,28,53],
        [39,80,89],
        [119,120,169],
        [77,36,85],
        [33,56,65],
        [65,72,97],
        [35,12,37]
    ]
    @test eltype(primitive_pythagorean_triples()) == Array{Int,1}
    @test eltype(primitive_pythagorean_triples(BigInt)) == Array{BigInt,1}

    @test collect(primitive_pythagorean_triples(11)) == Array{Int,1}[]
    @test collect(primitive_pythagorean_triples(12)) == Array{Int,1}[[3,4,5]]
    @test collect(primitive_pythagorean_triples(100)) == Array{Int,1}[
        [3,4,5],
        [5,12,13],
        [21,20,29],
        [15,8,17],
        [7,24,25],
        [35,12,37],
        [9,40,41]
    ]
    @test eltype(primitive_pythagorean_triples(100)) == Array{Int,1}
    @test eltype(primitive_pythagorean_triples(big(100))) == Array{BigInt,1}

    @test count(_->true, primitive_pythagorean_triples(10)) == 0
    @test count(_->true, primitive_pythagorean_triples(100)) == 7
    @test count(_->true, primitive_pythagorean_triples(1000)) == 70
    @test count(_->true, primitive_pythagorean_triples(10000)) == 703
    @test count(_->true, primitive_pythagorean_triples(100000)) == 7026
    @test count(_->true, primitive_pythagorean_triples(1000000)) == 70229
end

function test_triangle_all()
    print(rpad("Math.NumberTheory.Triangle...", 50, ' '))

    test_triangle_istriangular()
    test_triangle_nthtriangular()
    test_triangle_ntriangular()
    test_triangle_sometriangular()
    test_triangle_exacttriangular()

    test_triangle_primitive_pythagorean_triples()

    println("PASS")
end


text:
## 12.jl
# 1列目をcol1.txtに，2列目をcol2.txtに保存

f = open("./dat/hightemp.txt")
o1 = open("./dat/col1.txt", "w")
o2 = open("./dat/col2.txt", "w")
for line in readlines(f)
 col1,col2 = split(chomp(line), "\t")[1:2]
 println(o1, col1)
 println(o2, col2)
end
close(f) #省略可
close(o1)
close(o2)


text:
# Context management

export
    CuContext, destroy, push, pop

const CTX_SCHED_AUTO            = 0x00
const CTX_SCHED_SPIN            = 0x01
const CTX_SCHED_YIELD           = 0x02
const CTX_SCHED_BLOCKING_SYNC   = 0x04
const CTX_MAP_HOST              = 0x08
const CTX_LMEM_RESIZE_TO_MAX    = 0x10

immutable CuContext
    handle::Ptr{Void}

    function CuContext(dev::CuDevice, flags::Integer)
            pctx_box = ptrbox(Ptr{Void})
            @cucall(:cuCtxCreate, (Ptr{Ptr{Void}}, Cuint, Cint),
                                  pctx_box, flags, dev.handle)
            new(ptrunbox(pctx_box))
    end

    CuContext(dev::CuDevice) = CuContext(dev, 0)

    function CuContext()
            pctx_box = ptrbox(Ptr{Void})
            @cucall(:cuCtxGetCurrent, (Ptr{Ptr{Void}},), pctx_box)
            new(ptrunbox(pctx_box))
    end
end

function destroy(ctx::CuContext)
    @cucall(:cuCtxDestroy, (Ptr{Void},), ctx.handle)
end

function push(ctx::CuContext)
    @cucall(:cuCtxPushCurrent, (Ptr{Void},), ctx.handle)
end

function pop(ctx::CuContext)
    pctx_box = ptrbox(Ptr{Void})
    @cucall(:cuCtxPopCurrent, (Ptr{Ptr{Void}},), pctx_box)
    return CuContext(ptrunbox(pctx_box))
end

function device(ctx::CuContext)
    @assert CuContext() == ctx
    # TODO: cuCtxGetDevice returns the device ordinal, but as a CUDevice*?
    #       This can't be right...
    device_box = ptrbox(Cint)
    @cucall(:cuCtxGetDevice, (Ptr{Cint},), device_box)
    return CuDevice(ptrunbox(device_box))
end

synchronize(ctx::CuContext) = @cucall(:cuCtxSynchronize, (Ptr{Void},), ctx.handle)


text:
################################################################################
#
#  AlgAssRelOrd
#
################################################################################

# S is the element type of the base field of the algebra, T the fractional ideal
# type of this field
mutable struct AlgAssRelOrd{S, T, U} <: Ring
  algebra::U
  dim::Int
  pseudo_basis#::Vector{Tuple{AbsAlgAssElem{S}, T}}
  basis_matrix::Generic.MatSpaceElem{S}
  basis_mat_inv::Generic.MatSpaceElem{S}
  basis_pmatrix::PMat{S, T}

  disc # an integral ideal in the base field

  ismaximal::Int                   # 0 Not known
                                   # 1 Known to be maximal
                                   # 2 Known to not be maximal

  trred_matrix::Generic.MatSpaceElem{S}

  inv_coeff_ideals::Vector{T}

  isnice::Bool
  nice_order#Tuple{AlgAssAbsOrd, T}
  nice_order_ideal::T

  function AlgAssRelOrd{S, T, U}(A::AbsAlgAss{S}) where {S, T, U}
    z = new{S, T, U}()
    z.algebra = A
    z.dim = dim(A)
    z.ismaximal = 0
    z.isnice = false
    return z
  end

  function AlgAssRelOrd{S, T, U}(A::U, M::PMat{S, T}) where {S, T, U}
    z = AlgAssRelOrd{S, T, U}(A)
    z.basis_pmatrix = M
    z.basis_matrix = M.matrix
    return z
  end

  function AlgAssRelOrd{S, T, U}(A::U, M::Generic.MatSpaceElem{S}) where {S, T, U}
    z = AlgAssRelOrd{S, T, U}(A)
    z.basis_matrix = M
    z.basis_pmatrix = pseudo_matrix(M)
    return z
  end
end

################################################################################
#
#  AlgAssRelOrdElem
#
################################################################################

mutable struct AlgAssRelOrdElem{S, T, U} <: RingElem
  parent::AlgAssRelOrd{S, T, U}
  elem_in_algebra::AbsAlgAssElem{S}
  coordinates::Vector{S}
  has_coord::Bool

  function AlgAssRelOrdElem{S, T, U}(O::AlgAssRelOrd{S, T, U}) where {S, T, U}
    z = new{S, T, U}()
    z.parent = O
    z.elem_in_algebra = zero(algebra(O))
    z.coordinates = Vector{S}(undef, degree(O))
    z.has_coord = false
    return z
  end

  function AlgAssRelOrdElem{S, T, U}(O::AlgAssRelOrd{S, T, U}, a::AbsAlgAssElem{S}) where {S, T, U}
    z = new{S, T, U}()
    z.parent = O
    z.elem_in_algebra = a
    z.coordinates = Vector{S}(undef, degree(O))
    z.has_coord = false
    return z
  end

  function AlgAssRelOrdElem{S, T, U}(O::AlgAssRelOrd{S, T, U}, a::AbsAlgAssElem{S}, arr::Vector{S}) where {S, T, U}
    z = new{S, T, U}()
    z.parent = O
    z.elem_in_algebra = a
    z.coordinates = arr
    z.has_coord = true
    return z
  end
end

################################################################################
#
#  AlgAssRelOrdIdl
#
################################################################################

mutable struct AlgAssRelOrdIdl{S, T, U}
  algebra::U

  pseudo_basis::Vector{Tuple{AbsAlgAssElem{S}, T}}
  # The basis matrices are in the BASIS of the ALGEBRA!
  basis_pmatrix::PMat{S, T}
  basis_matrix::Generic.MatSpaceElem{S}
  basis_mat_inv::Generic.MatSpaceElem{S}

  # Basis pseudo-matrices with respect to orders
  basis_pmatrix_wrt::Dict{AlgAssRelOrd{S, T}, PMat{S, T}}

  # Left and right order:
  # The largest orders of which the ideal is a left resp. right ideal.
  left_order::AlgAssRelOrd{S, T, U}
  right_order::AlgAssRelOrd{S, T, U}

  # Any order contained in the left or right order, that is, an order of which
  # the ideal is a (possibly fractional) ideal.
  order::AlgAssRelOrd{S, T, U}

  # isleft and isright with respect to `order`
  isleft::Int                      # 0 Not known
                                   # 1 Known to be a left ideal
                                   # 2 Known not to be a left ideal
  isright::Int                     # as for isleft

  iszero::Int                      # 0: don't know, 1: known to be zero, 2: known to be not zero

  norm::Dict{AlgAssRelOrd{S, T, U}, T} # The ideal has different norms with respect
                                       # to different orders
  normred::Dict{AlgAssRelOrd{S, T, U}, T}

  function AlgAssRelOrdIdl{S, T, U}(A::AbsAlgAss{S}) where {S, T, U}
    z = new{S, T, U}()
    z.algebra = A
    z.isleft = 0
    z.isright = 0
    z.iszero = 0
    z.basis_pmatrix_wrt = Dict{AlgAssRelOrd{S, T, U}, PMat{S, T}}()
    z.norm = Dict{AlgAssRelOrd{S, T, U}, T}()
    z.normred = Dict{AlgAssRelOrd{S, T, U}, T}()
    return z
  end

  function AlgAssRelOrdIdl{S, T, U}(A::AbsAlgAss{S}, M::PMat{S, T}) where {S, T, U}
    z = AlgAssRelOrdIdl{S, T, U}(A)
    z.basis_pmatrix = M
    z.basis_matrix = M.matrix
    return z
  end
end


text:
using SparseDiffTools
using Graphs
using Graphs: SimpleGraph
using Random
#=
            Graph g0
vertex nuumber followed by color in parentheses

            6(3)
           /    \
          /      \
         /        \
       5(1)--------3(2)
       |  \        |
       |   \       |
       |    \      |
       |   1(2)----2(1)
       |   /
       |  /
      4(3)
=#

g0 = SimpleGraph(6)
add_edge!(g0, 1,2)
add_edge!(g0, 1,4)
add_edge!(g0, 1,5)
add_edge!(g0, 3,2)
add_edge!(g0, 3,5)
add_edge!(g0, 3,6)
add_edge!(g0, 4,5)
add_edge!(g0, 5,6)

#=
        Graph g1

          1
          |
          |
     5----2----3
         / \
        /   \
      6      4

=#
g1 = SimpleGraph(6)
add_edge!(g1, 2,1)
add_edge!(g1, 3,2)
add_edge!(g1, 4,2)
add_edge!(g1, 5,2)
add_edge!(g1, 6,2)

#=
        Graph g2

            1------2
          / |    / |
         /  |   /  |
        /   | /    |
       3----4------5

=#
g2 = SimpleGraph(5)
add_edge!(g2, 1,2)
add_edge!(g2, 1,3)
add_edge!(g2, 1,4)
add_edge!(g2, 4,2)
add_edge!(g2, 5,2)
add_edge!(g2, 3,4)
add_edge!(g2, 4,5)

#test custom graphs first
coloring0 = SparseDiffTools.color_graph(g0, SparseDiffTools.BacktrackingColor())
coloring1 = SparseDiffTools.color_graph(g1, SparseDiffTools.BacktrackingColor())
coloring2 = SparseDiffTools.color_graph(g2, SparseDiffTools.BacktrackingColor())

for v = 1:nv(g0)
    color = coloring0[v]
    for j in inneighbors(g0, v)
        if coloring0[j] == color
             @test false
        end
    end
end

for v = 1:nv(g1)
    color = coloring1[v]
    for j in inneighbors(g1, v)
        if coloring1[j] == color
             @test false
        end
    end
end

for v = 1:nv(g2)
    color = coloring2[v]
    for j in inneighbors(g2, v)
        if coloring2[j] == color
             @test false
        end
    end
end

test_graphs = Array{SimpleGraph, 1}(undef, 0)
for _ in 1:5
    nv = rand(5:20)
    ne = rand(1:100)
    graph = SimpleGraph(nv)
    for e in 1:ne
        v1 = rand(1:nv)
        v2 = rand(1:nv)
        while v1 == v2
            v2 = rand(1:nv)
        end
        add_edge!(graph, v1, v2)
    end
    push!(test_graphs, copy(graph))
end

for i in 1:5
    g = test_graphs[i]
    out_colors = SparseDiffTools.color_graph(g,SparseDiffTools.BacktrackingColor())

    for v = vertices(g)
        color = out_colors[v]
        for j in inneighbors(g, v)
            @test out_colors[j] != color
        end
    end
end