text:
# The content in this file is for test purpose only
"""
    cpvagk(f, t, a, b, tol=256*eps())

Calculate the Cauchy principle value integration of the form ``𝒫∫_a^b f(x)/(x-t) dx``. The algorithm is adapted from [P. Keller, 02.01.2015](https://www.sciencedirect.com/science/article/pii/S0377042715004422)
"""
function cpvagk(f, t, a, b, tol=256*eps())
    #   Adapted from P. Keller, 02.01.2015

    #   CPVAGK(F,T,A,B,TOL) = CPVInt[a,b]F(x)/(x-T)dx.
    #
    #   Based on QUADGK - the built-in Matlab adaptive integrator.
    #
    #   Q = CPVAGK(F,T,A,B,TOL) approximates the above integral
    #   and tries to make the absolute error less than TOL.
    #   The default value of TOL is 256*eps.
    #
    #   [Q,ERR] = CPVAGK(F,T,A,B,TOL) also provides the estimate of
    #   the absolute error. We suggest to always check the value of ERR.
    #   In case a QUADGK warning message appear, the approximate error bound
    #   provided inside the warning message should be ingnored!
    #   Consider only the value of ERR as error estimate.
    #
    #   Test (research) version.

    #   P. Keller, 02.01.2015

    if b<=t || t <= a
        error("CPVAGK: Invalid arguments")
    end
    s = min(1.0, t-a, b-t); #computing the split point(s)
    a = float(a)
    b = float(b)
    s = float(s)
    t = float(t)
    epsab = epsf = eps()
    ft = f(t)
    abt = 0.25*(abs(a) + abs(b) + 2*abs(t))
    tcond = 1 + 0.5*(max(1,abt)-1)^2/max(1,abt);
    #Estimating possible accuracy loss...
    # a) related to rounding t and interval endpoints:
    v = [a+2*epsf*abs(a), b-2*epsf*abs(b)];
    fd = max(max(1.0, abs(t)), abs.(v)...).*abs.(f.(v))./[t-a,b-t];
    erra = 0.75*epsab * sum(fd);
    # b) related to the magnitude of the 1st and 2nd derivative:
    sts = 255*s/256  # divided differences five steps...
    sth = min(sts, min(2,(b-a))./[82,70,32,22]...)
    ste = min(sts, exp(3/8*log(max(1,abs(t))*epsf)))
    dfs = max(abs(f(t+sts)-ft), abs(f(t-sts)-ft)) / sts  # <|f'[t-s,t+s]|
    dfh = max(abs(f(t+sth)-ft), abs(f(t-sth)-ft))./ sth./ [3/2, 7/4, 2, 3]
    df1 = abs(f(t+ste^2)-f(t-ste^2))/ste^2;  # ~~ 2*|f'(t)|
    df2 = sqrt(abs(f(t+ste)-2*ft+f(t-ste)))/ste;  # ~~ |f"(t)|^(1/2)
    errb1 = 17*tcond*epsf*(max(df1,dfh...) + abs(ft));
    errb2 = 10*tcond*epsf*df2;
    # c) related to the distance from interval endpoints:
    errc = 0.75*epsf*abt*abs(ft)/min(b-t,t-a)
    # Setting (safe) error tolerances, absolute and relative...
    atol = max(tol/2, 8*dfs*epsf, erra, errb1, errb2, errc);
    rtol = max(100*epsf, errb2, errc/max(sqrt(epsf), abs(ft)));
    # Setting parameters for quadgk()
    params = (atol=atol, rtol=rtol, maxevals=163840)
    # Computing the integrals...
    q1 = 0
    e1 = 0
    q2 = 0
    e2 = 0
    if s < t-a
        q1, e1 = quadgk((x)->(f(x)-ft)./(x-t), a, t-s; params...)
    end
    if s < b-t
        q2, e2 = quadgk((x)->(f(x)-ft)./(x-t), t+s, b; params...)
    end
    q3, e3 = quadgk((x)->(f(t+x)-f(t-x))./x, 0, s; params...)
    # Computing the final result and error estimation...
    q = q1 + q2 + q3 + ft*log((b-t)/(t-a))
    err = e1 + e2 + e3 + max(erra,errc) + errb1 + errb2 + 10*tcond*epsf*max(abs.([q1,q2,q3,q])...)
    q, err
end

"""
$(SIGNATURES)

Calculate the Lamb shift of spectrum `γ`. `atol` is the absolute tolerance for Cauchy principal value integral.
"""
function lambshift_cpvagk(w, γ; atol = 1e-7)
    g(x) = γ(x) / (x - w)
    cpv, cperr = cpvagk(γ, w, w - 1.0, w + 1.0)
    negv, negerr = quadgk(g, -Inf, w - 1.0)
    posv, poserr = quadgk(g, w + 1.0, Inf)
    v = cpv + negv + posv
    err = cperr + negerr + poserr
    if (err > atol) || (isnan(err))
        @warn "Absolute error of integration is larger than the tolerance."
    end
    -v / 2 / pi
end


text:
module LowRankDPP

include("DPPDataPreparation.jl")
include("DPPLearning.jl")
include("DPPLearningBayesian.jl")
include("DPPPrediction.jl")
include("DPPEvaluation.jl")
include("DPPExamples.jl")

end

ext: jl


text:
using Test
using YaoSym, YaoArrayRegister, YaoBlocks, SymEngine

@test ket"111" + 2ket"111" == 3ket"111"
@test bra"111" * ket"111" == 1
@test (bra"111" + bra"101") * ket"111" == 1

@test ket"111"^2 == ket"111111"
@test ket"101" * ket"111" == ket"101111"
@test bra"110"^2 == bra"110110"
@test bra"110" * bra"111" == bra"110111" == bra"110_111"

# printing does not error
@test print(ket"001") == nothing
@test print(bra"001") == nothing
@test print(bra"111" + bra"101")  == nothing


text:
using RandomExtensions, Random, SparseArrays
using Test

@testset "Distributions" begin
    # Normal/Exponential
    @test rand(Normal()) isa Float64
    @test rand(Normal(0.0, 1.0)) isa Float64
    @test rand(Exponential()) isa Float64
    @test rand(Exponential(1.0)) isa Float64
    @test rand(Normal(Float32)) isa Float32
    @test rand(Exponential(Float32)) isa Float32
    @test rand(Normal(ComplexF64)) isa ComplexF64

    # pairs/complexes
    @test rand(Combine(Pair, 1:3, Float64)) isa Pair{Int,Float64}
    z = rand(Combine(Complex, 1:3, 6:9))
    @test z.re ∈ 1:3
    @test z.im ∈ 6:9
    @test z isa Complex{Int}
    z = rand(Combine(ComplexF64, 1:3, 6:9))
    @test z.re ∈ 1:3
    @test z.im ∈ 6:9
    @test z isa ComplexF64

    # Uniform
    @test rand(Uniform(Float64)) isa Float64
    @test rand(Uniform(1:10)) isa Int
    @test rand(Uniform(1:10)) ∈ 1:10
    @test rand(Uniform(Int)) isa Int
end

@testset "Containers" for rng in ([], [MersenneTwister(0)], [RandomDevice()])
    # Set
    for s = (rand(rng..., 1:99, Set{Int}, 10),
             rand(rng..., 1:99, Set, 10))
        @test s isa Set{Int}
        @test length(s) == 10
        @test rand(s) ∈ 1:99
    end
    s = Set([1, 2])
    @test s === rand!(s)
    @test first(s) ∉ (1, 2) # extremely unlikely
    @test length(s) == 2
    @test s === rand!(s, 3:9) <= Set(3:9)
    @test length(s) == 2

    # Dict
    for s = (rand(rng..., Combine(Pair, 1:99, 1:99), Dict, 10),
             rand(rng..., Combine(Pair, 1:99, 1:99), Dict{Int,Int}, 10))
        @test s isa Dict{Int,Int}
        @test length(s) == 10
        p = rand(s)
        @test p.first ∈ 1:99
        @test p.second ∈ 1:99
    end
    s = Dict(1=>2, 2=>1)
    @test s === rand!(s)
    @test length(s) == 2
    @test first(s).first ∉ (1, 2) # extremely unlikely
    rand!(s, Combine(Pair, 3:9, Int))
    @test length(s) == 2
    @test first(s).first ∈ 3:9

    # sparse
    @test rand(rng..., Float64, .5, 10) isa SparseVector{Float64}
    @test rand(rng..., .5, 10) isa SparseVector{Float64}
    @test rand(rng..., Int, .5, 10) isa SparseVector{Int}
    @test rand(rng..., Float64, .5, 10, 3) isa SparseMatrixCSC{Float64}
    @test rand(rng..., .5, 10, 3) isa SparseMatrixCSC{Float64}
    @test rand(rng..., Int, .5, 10, 3) isa SparseMatrixCSC{Int}

    # BitArray
    @test rand(rng..., BitArray, 10) isa BitVector
    @test rand(rng..., BitVector, 10) isa BitVector
    @test_throws MethodError rand(rng..., BitVector, 10, 20) isa BitVector
    @test rand(rng..., BitArray, 10, 3) isa BitMatrix
    @test rand(rng..., BitMatrix, 10, 3) isa BitMatrix
    @test_throws MethodError rand(rng..., BitVector, 10, 3) isa BitMatrix

    # String
    s = rand(rng..., String)
    @test s isa String
    @test length(s) == 8
    s = rand(rng..., String, 10)
    @test s isa String
    @test length(s) == 10
    s = rand(rng..., "asd", String)
    @test length(s) == 8
    @test Set(s) <= Set("asd")
end

@testset "Rand" for rng in ([], [MersenneTwister(0)], [RandomDevice()])
    for XT = zip(([Int], [1:3], []), (Int, Int, Float64))
        X, T = XT
        r = Rand(rng..., X...)
        @test collect(Iterators.take(r, 10)) isa Vector{T}
        @test r() isa T
        @test r(2, 3) isa Matrix{T}
        @test r(.3, 2, 3) isa SparseMatrixCSC{T}
    end
    for d = (Uniform(1:10), Uniform(Int))
        @test collect(Iterators.take(d, 10)) isa Vector{Int}
    end
end


text:
include("events.jl")

const DEFAULT_TAGS = Dict{AbstractString,AbstractString}("!" => "!", "!!" => "tag:yaml.org,2002:")

struct ParserError
    context::Union{AbstractString, Nothing}
    context_mark::Union{Mark, Nothing}
    problem::Union{AbstractString, Nothing}
    problem_mark::Union{Mark, Nothing}
    note::Union{AbstractString, Nothing}

    function ParserError(context=nothing, context_mark=nothing,
                         problem=nothing, problem_mark=nothing,
                         note=nothing)
        new(context, context_mark, problem, problem_mark, note)
    end
end

function show(io::IO, error::ParserError)
    if error.context != nothing
        print(io, error.context, " at ", error.context_mark, ": ")
    end
    print(io, error.problem, " at ", error.problem_mark)
end

mutable struct EventStream
    input::TokenStream
    next_event::Union{Event, Nothing}
    state::Union{Function, Nothing}
    states::Vector{Function}
    marks::Vector{Mark}
    yaml_version::Union{Tuple, Nothing}
    tag_handles::Dict{AbstractString, AbstractString}
    end_of_stream::Union{StreamEndEvent, Nothing}

    function EventStream(input::TokenStream)
        new(input, nothing, parse_stream_start, Function[], Mark[],
            nothing, Dict{AbstractString, AbstractString}(), nothing)
    end
end

function peek(stream::EventStream)
    if stream.next_event === nothing
        if stream.state === nothing
            return nothing
        elseif stream.end_of_stream !== nothing
            stream.state = nothing
            return stream.end_of_stream
        else
            x = stream.state(stream)
            #@show x
            stream.next_event = x
            #stream.next_event = stream.state(stream)
        end
    end

    return stream.next_event
end

function forward!(stream::EventStream)
    if stream.next_event === nothing
        if stream.state === nothing
            nothing
        elseif stream.end_of_stream !== nothing
            stream.state = nothing
            return stream.end_of_stream
        else
            stream.next_event = stream.state(stream)
        end
    end

    e = stream.next_event
    stream.next_event = nothing
    return e
end

function process_directives(stream::EventStream)
    stream.yaml_version = nothing
    stream.tag_handles = Dict{AbstractString, AbstractString}()
    while typeof(peek(stream.input)) == DirectiveToken
        token = forward!(stream.input)
        if token.name == "YAML"
            if stream.yaml_version !== nothing
                throw(ParserError(nothing, nothing,
                                  "found duplicate YAML directive",
                                  token.start_mark))
            end
            major, minor = token.value
            if major != 1
                throw(ParserError(nothing, nothing,
                    "found incompatible YAML document (version 1.* is required)",
                    token.start_mark))
            end
            stream.yaml_version = token.value
        elseif token.name == "TAG"
            handle, prefix = token.value
            if haskey(stream.tag_handles, handle)
                throw(ParserError(nothing, nothing,
                    "duplicate tag handle $(handle)", token.start_mark))
            end
            stream.tag_handles[handle] = prefix
        end
    end

    if stream.tag_handles !== nothing
        value = stream.yaml_version, copy(stream.tag_handles)
    else
        value = stream.yaml_version, nothing
    end

    for key in keys(DEFAULT_TAGS)
        if !haskey(stream.tag_handles, key)
            stream.tag_handles[key] = DEFAULT_TAGS[key]
        end
    end

    value
end

# Parser state functions

function parse_stream_start(stream::EventStream)
    token = forward!(stream.input) :: StreamStartToken
    event = StreamStartEvent(token.span.start_mark, token.span.end_mark,
                             token.encoding)
    stream.state = parse_implicit_document_start
    event
end

function parse_implicit_document_start(stream::EventStream)
    token = peek(stream.input)
    if !in(typeof(token), [DirectiveToken, DocumentStartToken, StreamEndToken])
        stream.tag_handles = DEFAULT_TAGS
        event = DocumentStartEvent(token.span.start_mark, token.span.start_mark,
                                   false)

        push!(stream.states, parse_document_end)
        stream.state = parse_block_node

        event
    else
        parse_document_start(stream)
    end
end

function parse_document_start(stream::EventStream)
    # Parse any extra document end indicators.
    while typeof(peek(stream.input)) == DocumentEndToken
        stream.input = Iterators.rest(stream.input)
    end

    # Parse explicit document.
    token = peek(stream.input)
    if typeof(token) != StreamEndToken
        start_mark = token.span.start_mark
        version, tags = process_directives(stream)
        if typeof(peek(stream.input)) != DocumentStartToken
            throw(ParserError(nothing, nothing,
                "expected '<document start>' but found $(typeof(token))"))
        end
        token = forward!(stream.input)
        event = DocumentStartEvent(start_mark, token.span.end_mark,
                                   true, version, tags)
        push!(stream.states, parse_document_end)
        stream.state = parse_document_content
        event
    else
        # Parse the end of the stream
        token = forward!(stream.input)
        event = StreamEndEvent(token.span.start_mark, token.span.end_mark)
        @assert isempty(stream.states)
        @assert isempty(stream.marks)
        stream.state = nothing
        event
    end
end

function parse_document_end(stream::EventStream)
    token = peek(stream.input)
    start_mark = end_mark = token.span.start_mark
    explicit = false
    if typeof(token) == DocumentEndToken
        forward!(stream.input)
        end_mark = token.span.end_mark
        explicit = true
        stream.end_of_stream = StreamEndEvent(token.span.start_mark,
                                              token.span.end_mark)
    end
    event = DocumentEndEvent(start_mark, end_mark, explicit)
    stream.state = parse_document_start
    event
end

function parse_document_content(stream::EventStream)
    if in(peek(stream.input), [DirectiveToken, DocumentStartToken, DocumentEndToken,StreamEndToken])
        event = process_empty_scalar(stream, peek(stream.input).span.start_mark)
        stream.state = pop!(stream.states)
        event
    else
        parse_block_node(stream)
    end
end

function parse_block_node(stream::EventStream)
    parse_node(stream, block=true)
end

function parse_flow_node(stream::EventStream)
    parse_node(stream)
end

function parse_block_node_or_indentless_sequence(stream::EventStream)
    parse_node(stream, block=true, indentless_sequence=true)
end

function parse_node(stream::EventStream; block=false, indentless_sequence=false)
    token = peek(stream.input)
    if typeof(token) == AliasToken
        forward!(stream.input)
        stream.state = pop!(stream.states)
        return AliasEvent(token.span.start_mark, token.span.end_mark, token.value)
    end

    anchor = nothing
    tag = nothing
    start_mark = end_mark = tag_mark = nothing
    if typeof(token) == AnchorToken
        forward!(stream.input)
        start_mark = token.span.start_mark
        end_mark = token.span.end_mark
        anchor = token.value
        token = peek(stream.input)
        if typeof(token) == TagToken
            forward!(stream.input)
            tag_mark = token.span.start_mark
            end_mark = token.span.end_mark
            tag = token.value
        end
    elseif typeof(token) == TagToken
        forward!(stream.input)
        start_mark = token.span.start_mark
        end_mark = token.span.end_mark
        tag = token.value
        token = peek(stream.input)
        if typeof(token) == AnchorToken
            forward!(stream.input)
            end_mark = token.end_mark
            anchor = token.value
        end
    end

    if tag !== nothing
        handle, suffix = tag
        if handle !== nothing
            if !haskey(stream.tag_handles, handle)
                throw(ParserError("while parsing a node", start_mark,
                                  "found undefined tag handle $(handle)",
                                  tag_mark))
            end
            tag = string(stream.tag_handles[handle], suffix)
        else
            tag = suffix
        end
    end

    token = peek(stream.input)
    if start_mark == nothing
        start_mark = end_mark = token.span.start_mark
    end

    event = nothing
    implicit = tag === nothing || tag == "!"
    if indentless_sequence && typeof(token) == BlockEntryToken
        end_mark = token.span.end_mark
        event = SequenceStartEvent(start_mark, end_mark, anchor, tag, implicit,
                                   false)
        stream.state = parse_indentless_sequence_entry
    else
        if typeof(token) == ScalarToken
            forward!(stream.input)
            end_mark = token.span.end_mark
            if (token.plain && tag == nothing) || tag == "!"
                implicit = true, false
            elseif tag == nothing
                implicit = false, true
            else
                implicit = false, false
            end
            event = ScalarEvent(start_mark, end_mark, anchor, tag, implicit,
                                token.value, token.style)
            stream.state = pop!(stream.states)
        elseif typeof(token) == FlowSequenceStartToken
            end_mark = token.span.end_mark
            event = SequenceStartEvent(start_mark, end_mark, anchor, tag,
                                       implicit, true)
            stream.state = parse_flow_sequence_first_entry
        elseif typeof(token) == FlowMappingStartToken
            end_mark = token.span.end_mark
            event = MappingStartEvent(start_mark, end_mark, anchor, tag,
                                      implicit, true)
            stream.state = parse_flow_mapping_first_key
        elseif block && typeof(token) == BlockSequenceStartToken
            end_mark = token.span.start_mark
            event = SequenceStartEvent(start_mark, end_mark, anchor, tag,
                                       implicit, false)
            stream.state = parse_block_sequence_first_entry
        elseif block && typeof(token) == BlockMappingStartToken
            end_mark = token.span.start_mark
            event = MappingStartEvent(start_mark, end_mark, anchor, tag,
                                      implicit, false)
            stream.state = parse_block_mapping_first_key
        elseif anchor !== nothing || tag !== nothing
            event = ScalarEvent(start_mark, end_mark, anchor, tag,
                                (implicit, false), "")
            stream.state = pop!(stream.states)
        else
            node = block ? "block" : "flow"
            throw(ParserError("while parsing a $(node) node", start_mark,
                    "expected the node content, but found $(typeof(token))",
                    token.span.start_mark))
        end
    end

    event
end

function parse_block_sequence_first_entry(stream::EventStream)
    token = forward!(stream.input)
    push!(stream.marks, token.span.start_mark)
    parse_block_sequence_entry(stream)
end

function parse_block_sequence_entry(stream::EventStream)
    token = peek(stream.input)
    if typeof(token) == BlockEntryToken
        forward!(stream.input)
        if !in(typeof(peek(stream.input)), [BlockEntryToken, BlockEndToken])
            push!(stream.states, parse_block_sequence_entry)
            return parse_block_node(stream)
        else
            stream.state = parse_block_sequence_entry
            return process_empty_scalar(stream, token.span.end_mark)
        end
    end

    if typeof(token) != BlockEndToken
        throw(ParserError("while parsing a block collection", stream.marks[end],
                          "expected <block end>, but found $(typeof(token))",
                          token.span.start_mark))
    end

    forward!(stream.input)
    pop!(stream.marks)
    stream.state = pop!(stream.states)
    SequenceEndEvent(token.span.start_mark, token.span.end_mark)
end

function parse_indentless_sequence_entry(stream::EventStream)
    token = peek(stream.input)
    if typeof(token) == BlockEntryToken
        forward!(stream.input)
        if !in(typeof(peek(stream.input)), [BlockEntryToken, KeyToken, ValueToken, BlockEndToken])
            push!(stream.states, parse_indentless_sequence_entry)
            return parse_block_node(stream)
        else
            stream.state = parse_indentless_sequence_entry
            return process_empty_scalar(stream, token.end_mark)
        end
    end

    stream.state = pop!(stream.states)
    SequenceEndEvent(token.span.start_mark, token.span.end_mark)
end

function parse_block_mapping_first_key(stream::EventStream)
    token = forward!(stream.input)
    push!(stream.marks, token.span.start_mark)
    parse_block_mapping_key(stream)
end

function parse_block_mapping_key(stream::EventStream)
    token = peek(stream.input)
    if typeof(token) == KeyToken
        forward!(stream.input)
        if !in(typeof(peek(stream.input)), [KeyToken, ValueToken, BlockEndToken])
            push!(stream.states, parse_block_mapping_value)
            return parse_block_node_or_indentless_sequence(stream)
        else
            stream.state = parse_block_mapping_value
            return process_empty_scalar(stream, token.span.end_mark)
        end
    end

    if typeof(token) != BlockEndToken
        throw(ParserError("while parsing a block mapping", stream.marks[end],
                          "expected <block end>, but found $(typeof(token))",
                          token.span.start_mark))
    end

    forward!(stream.input)
    pop!(stream.marks)
    stream.state = pop!(stream.states)
    MappingEndEvent(token.span.start_mark, token.span.end_mark)
end

function parse_block_mapping_value(stream::EventStream)
    token = peek(stream.input)
    if typeof(token) == ValueToken
        forward!(stream.input)
        if !in(typeof(peek(stream.input)), [KeyToken, ValueToken, BlockEndToken])
            push!(stream.states, parse_block_mapping_key)
            parse_block_node_or_indentless_sequence(stream)
        else
            stream.state = parse_block_mapping_key
            process_empty_scalar(stream, token.span.end_mark)
        end
    else
        stream.state = parse_block_mapping_key
        process_empty_scalar(stream, token.span.start_mark)
    end
end

function parse_flow_sequence_first_entry(stream::EventStream)
    token = forward!(stream.input)
    push!(stream.marks, token.span.start_mark)
    parse_flow_sequence_entry(stream, first_entry=true)
end

function parse_flow_sequence_entry(stream::EventStream; first_entry=false)
    token = peek(stream.input)
    if typeof(token) != FlowSequenceEndToken
        if !first_entry
            if typeof(token) == FlowEntryToken
                forward!(stream.input)
            else
                throw(ParserError("while parsing a flow sequence",
                                  stream.mark[end],
                                  "expected ',' or ']', but got $(typeof(token))",
                                  token.span.start_mark))
            end
        end

        token = peek(stream.input)
        if typeof(token) == KeyToken
            stream.state = parse_flow_sequence_entry_mapping_key
            return MappingStartEvent(token.span.start_mark, token.span.end_mark,
                                     nothing, nothing, true, true)
        elseif typeof(token) != FlowSequenceEndToken
            push!(stream.states, parse_flow_sequence_entry)
            return parse_flow_node(stream)
        end
    end

    forward!(stream.input)
    pop!(stream.marks)
    stream.state = pop!(stream.states)
    SequenceEndEvent(token.span.start_mark, token.span.end_mark)
end

function parse_flow_sequence_entry_mapping_key(stream::EventStream)
    token = forward!(stream.input)
    if !in(typeof(token), [ValueToken, FlowEntryToken, FlowSequenceEndToken])
        push!(stream.states, parse_flow_sequence_entry_mapping_value)
        parse_flow_node(stream)
    else
        stream.state = parse_flow_sequence_entry_mapping_value
        process_empty_scalar(stream, token.span.end_mark)
    end
end

function parse_flow_sequence_entry_mapping_value(stream::EventStream)
    token = peek(stream.input)
    if typeof(token) == ValueToken
        forward!(stream.input)
        if !in(typeof(peek(stream.input)), [FlowEntryToken, FlowSequenceEndToken])
            push!(stream.states, parse_flow_sequence_entry_mapping_end)
            parse_flow_node(stream)
        else
            stream.state = parse_flow_sequence_entry_mapping_end
            process_empty_scalar(stream, token.span.end_mark)
        end
    else
        stream.state = parse_flow_sequence_entry_mapping_end
        process_empty_scalar(stream, token.span.start_mark)
    end
end

function parse_flow_sequence_entry_mapping_end(stream::EventStream)
    stream.state = parse_flow_sequence_entry
    token = peek(stream.input)
    MappingEndEvent(token.span.start_mark, token.span.end_mark)
end

function parse_flow_mapping_first_key(stream::EventStream)
    token = forward!(stream.input)
    push!(stream.marks, token.span.start_mark)
    parse_flow_mapping_key(stream, first_entry=true)
end

function parse_flow_mapping_key(stream::EventStream; first_entry=false)
    token = peek(stream.input)
    if typeof(token) != FlowMappingEndToken
        if !first_entry
            if typeof(token) == FlowEntryToken
                forward!(stream.input)
            else
                throw(ParserError("while parsing a flow mapping",
                                  stream.marks[end],
                                  "expected ',' or '}', but got $(typeof(token))",
                                  token.span.start_mark))
            end
        end

        token = peek(stream.input)
        if typeof(token) == KeyToken
            forward!(stream.input)
            if !in(typeof(peek(stream.input)), [ValueToken, FlowEntryToken, FlowMappingEndToken])
                push!(stream.states, parse_flow_mapping_value)
                return parse_flow_node(stream)
            else
                stream.state = parse_flow_mapping_value
                return process_empty_scalar(stream, token.span.end_mark)
            end
        elseif typeof(token) != FlowMappingEndToken
            push!(stream.states, parse_flow_mapping_empty_value)
            return parse_flow_node(stream)
        end
    end

    forward!(stream.input)
    pop!(stream.marks)
    stream.state = pop!(stream.states)
    MappingEndEvent(token.span.start_mark, token.span.end_mark)
end

function parse_flow_mapping_value(stream::EventStream)
    token = peek(stream.input)
    if typeof(token) == ValueToken
        forward!(stream.input)
        if !in(typeof(peek(stream.input)), [FlowEntryToken, FlowMappingEndToken])
            push!(stream.states, parse_flow_mapping_key)
            parse_flow_node(stream)
        else
            stream.state = parse_flow_mapping_key
            process_empty_scalar(stream, token.span.end_mark)
        end
    else
        stream.state = parse_flow_mapping_key
        process_empty_scalar(stream, token.span.start_mark)
    end
end

function parse_flow_mapping_empty_value(stream::EventStream)
    stream.state = parse_flow_mapping_key
    process_empty_scalar(stream, peek(stream.input).span.start_mark)
end

function process_empty_scalar(stream::EventStream, mark::Mark)
    ScalarEvent(mark, mark, nothing, nothing, (true, false), "", nothing)
end


text:
# Fallback constructor for diffusivity types without precomputed diffusivities.
DiffusivityFields(arch, grid, tracer_names, bcs, closure) = nothing
DiffusivityFields(arch, grid, tracer_names, closure) = nothing

# For a tuple of closures
DiffusivityFields(arch, grid, tracer_names, bcs, closure_tuple::Tuple) =
    Tuple(DiffusivityFields(arch, grid, tracer_names, bcs, closure) for closure in closure_tuple)

#####
##### For closures that only require an eddy viscosity νₑ field.
#####

const ViscosityClosures = Union{AbstractSmagorinsky, AbstractLeith}

DiffusivityFields(arch, grid, tracer_names, ::ViscosityClosures;
                  νₑ = CellField(arch, grid, DiffusivityBoundaryConditions(grid))) = (νₑ=νₑ,)

function DiffusivityFields(arch, grid, tracer_names, bcs, ::ViscosityClosures)

    νₑ_bcs = :νₑ ∈ keys(bcs) ? bcs[:νₑ] : DiffusivityBoundaryConditions(grid)

    νₑ = CellField(arch, grid, νₑ_bcs)

    return (νₑ=νₑ,)
end

#####
##### For closures that also require tracer diffusivity fields κₑ on each tracer.
#####

const ViscosityDiffusivityClosures = Union{VAMD, RAMD}

function TracerDiffusivityFields(arch, grid, tracer_names; kwargs...)

    κ_fields = Tuple(c ∈ keys(kwargs) ?
                     kwargs[c] :
                     CellField(arch, grid, DiffusivityBoundaryConditions(grid))
                     for c in tracer_names)

    return NamedTuple{tracer_names}(κ_fields)
end

function TracerDiffusivityFields(arch, grid, tracer_names, bcs)

    κ_fields = Tuple(c ∈ keys(bcs) ? CellField(arch, grid, bcs[c]) :
                                     CellField(arch, grid, DiffusivityBoundaryConditions(grid))
                     for c in tracer_names)

    return NamedTuple{tracer_names}(κ_fields)
end

function DiffusivityFields(arch, grid, tracer_names, ::ViscosityDiffusivityClosures;
                           νₑ = CellField(arch, grid, DiffusivityBoundaryConditions(grid)),
                           kwargs...)

    κₑ = TracerDiffusivityFields(arch, grid, tracer_names; kwargs...)

    return (νₑ=νₑ, κₑ=κₑ)
end

function DiffusivityFields(arch, grid, tracer_names, bcs, ::ViscosityDiffusivityClosures)

    νₑ_bcs = :νₑ ∈ keys(bcs) ? bcs[:νₑ] : DiffusivityBoundaryConditions(grid)

    νₑ = CellField(arch, grid, νₑ_bcs)

    κₑ = :κₑ ∈ keys(bcs) ? TracerDiffusivityFields(arch, grid, tracer_names, bcs[:κₑ]) :
                           TracerDiffusivityFields(arch, grid, tracer_names)

    return (νₑ=νₑ, κₑ=κₑ)
end


text:
using DiffEqBase

####################################################################################################
# this function takes into accound a parameter passed to the vector field
# Putting the options `save_start=false` seems to give bugs with Sundials
function flowTimeSol(x, p, tm, pb::ODEProblem; alg = Euler(), kwargs...)
	_prob = DiffEqBase.remake(pb; u0 = x, tspan = (zero(tm), tm), p = p)
	# the use of concrete_solve makes it compatible with Zygote
	sol = DiffEqBase.concrete_solve(_prob, alg; save_everystep = false, kwargs...)
	return (t = sol.t[end], u = sol[end])
end

# this function is a bit different from the previous one as it is geared towards parallel computing of the flow.
function flowTimeSol(x::AbstractMatrix, p, tm, epb::EnsembleProblem; alg = Euler(), kwargs...)
	# modify the function which asigns new initial conditions
	# see docs at https://docs.sciml.ai/dev/features/ensemble/#Performing-an-Ensemble-Simulation-1
	_prob_func = (prob, ii, repeat) -> prob = remake(prob, u0 = x[:, ii], tspan = (zero(tm[ii]), tm[ii]), p = p)
	_epb = setproperties(epb, output_func = (sol, i) -> ((t = sol.t[end], u = sol[end]), false), prob_func = _prob_func)
	sol = DiffEqBase.solve(_epb, alg, EnsembleThreads(); trajectories = size(x, 2), save_everystep = false, kwargs...)
	# sol.u contains a vector of tuples (sol_i.t[end], sol_i[end])
	return sol.u
end

flow(x, p, tm, pb::ODEProblem; alg = Euler(), kwargs...) = flowTimeSol(x, p, tm, pb; alg = alg, kwargs...).u
flow(x, p, tm, pb::EnsembleProblem; alg = Euler(), kwargs...) = flowTimeSol(x, p, tm, pb; alg = alg, kwargs...)
flow(x, tm, pb::Union{ODEProblem, EnsembleProblem}; alg = Euler(), kwargs...) =  flow(x, nothing, tm, pb; alg = alg, kwargs...)
################
# function used to compute the derivative of the flow, so pb encodes the variational equation
function dflow(x::AbstractVector, p, dx, tm, pb::ODEProblem; alg = Euler(), kwargs...)
	n = length(x)
	_prob = DiffEqBase.remake(pb; u0 = vcat(x, dx), tspan = (zero(tm), tm), p = p)
	# the use of concrete_solve makes it compatible with Zygote
	sol = DiffEqBase.concrete_solve(_prob, alg, save_everystep = false; kwargs...)[end]
	return (t = tm, u = sol[1:n], du = sol[n+1:end])
end

# same for Parallel computing
function dflow(x::AbstractMatrix, p, dx, tm, epb::EnsembleProblem; alg = Euler(), kwargs...)
	N = size(x,1)
	_prob_func = (prob, ii, repeat) -> prob = remake(prob, u0 = vcat(x[:, ii], dx[:, ii]), tspan = (zero(tm[ii]), tm[ii]), p = p)
	_epb = setproperties(epb, output_func = (sol,i) -> ((t = sol.t[end], u = sol[end][1:N], du = sol[end][N+1:end]), false), prob_func = _prob_func)
	sol = DiffEqBase.solve(_epb, alg, EnsembleThreads(); trajectories = size(x, 2), save_everystep = false, kwargs...)
	return sol.u
end

dflow(x, dx, tspan, pb::Union{ODEProblem, EnsembleProblem}; alg = Euler(), kwargs...) = dflow(x, nothing, dx, tspan, pb; alg = alg, kwargs...)
################
# this function takes into accound a parameter passed to the vector field
function dflow_fd(x, p, dx, tm, pb::ODEProblem; alg = Euler(), δ = 1e-9, kwargs...)
	sol1 = flow(x .+ δ .* dx, p, tm, pb; alg = alg, kwargs...)
	sol2 = flow(x 			, p, tm, pb; alg = alg, kwargs...)
	return (t = tm, u = sol2, du = (sol1 .- sol2) ./ δ)
end

function dflow_fd(x, p, dx, tm, pb::EnsembleProblem; alg = Euler(), δ = 1e-9, kwargs...)
	sol1 = flow(x .+ δ .* dx, p, tm, pb; alg = alg, kwargs...)
	sol2 = flow(x 			, p, tm, pb; alg = alg, kwargs...)
	return [(t = sol1[ii][1], u = sol2[ii][2], du = (sol1[ii][2] .- sol2[ii][2]) ./ δ) for ii = 1:size(x,2) ]
end
dflow_fd(x, dx, tm, pb::Union{ODEProblem, EnsembleProblem}; alg = Euler(), δ = 1e-9, kwargs...) = dflow_fd(x, nothing, dx, tm, pb; alg = alg, δ = δ, kwargs...)
################
# this gives access to the full solution, convenient for Poincaré shooting
# this function takes into accound a parameter passed to the vector field and returns the full solution from the ODE solver. This is useful in Poincare Shooting to extract the period.
function flowFull(x, p, tm, pb::ODEProblem; alg = Euler(), kwargs...)
	_prob = DiffEqBase.remake(pb; u0 = x, tspan = (zero(tm), tm), p = p)
	sol = DiffEqBase.solve(_prob, alg; kwargs...)
end

function flowFull(x, p, tm, epb::EnsembleProblem; alg = Euler(), kwargs...)
	_prob_func = (prob, ii, repeat) -> prob = remake(prob, u0 = x[:, ii], tspan = (zero(tm[ii]), tm[ii]), p = p)
	_epb = setproperties(epb, output_func = (sol,i) -> ((t = sol.t, u = sol), false), prob_func = _prob_func)
	sol = DiffEqBase.solve(_epb, alg, EnsembleThreads(); trajectories = size(x, 2), kwargs...)
end
flowFull(x, tm, pb::Union{ODEProblem, EnsembleProblem}; alg = Euler(), kwargs...) = flowFull(x, nothing, tm, pb; alg = alg, kwargs...)
####################################################################################################
# Structures related to computing ODE/PDE Flows
"""
$(TYPEDEF)
$(TYPEDFIELDS)

## Simplified constructors
There are some simple constructors for which you only have to pass a `prob::ODEProblem` or `prob::EnsembleProblem` (for parallel computation) from `DifferentialEquations.jl` and an ODE time stepper like `Tsit5()`. Hence, you can do for example

	fl = Flow(F, prob, Tsit5(); kwargs...)

If your vector field depends on parameters `p`, you can define a `Flow` using

	fl = Flow(F, p, prob, Tsit5(); kwargs...)

Finally, you can pass two `ODEProblem` where the second one is used to compute the variational equation:

	fl = Flow(F, p, prob1::ODEProblem, alg1, prob2::ODEProblem, alg2; kwargs...)

"""
struct Flow{TF, Tf, Tts, Tff, Td, Tse}
	"The vector field `(x, p) -> F(x, p)` associated to a Cauchy problem,"
	F::TF

	"The flow (or semigroup) associated to the Cauchy problem `(x, p, t) -> flow(x, p, t)`. Only the last time point must be returned."
	flow::Tf

	"Flow which returns the tuple (t, u(t))"
	flowTimeSol::Tts

	"The flow (or semigroup) associated to the Cauchy problem `(x, p, t) -> flow(x, p, t)`. The whole solution on the time interval (0,t) must be returned. It is not strictly necessary to provide this."
	flowFull::Tff

	"The differential `dflow` of the flow w.r.t. `x`, `(x, p, dx, t) -> dflow(x, p, dx, t)`. One important thing is that we require `dflow(x, dx, t)` to return a Named Tuple: `(t = t, u = flow(x, p, t), du = dflow(x, p, dx, t))`, the last composant being the value of the derivative of the flow."
	dflow::Td

	"Serial version of dflow"
	dfserial::Tse
end

# constructors
Flow() = Flow(nothing, nothing, nothing, nothing, nothing, nothing)
Flow(F, fl, df) = Flow(F, fl, nothing, nothing, df, df)

# callable struct
(fl::Flow)(x, p, t)     			  = fl.flow(x, p, t)
(fl::Flow)(x, p, dx, t; kw2...) 	  = fl.dflow(x, p, dx, t; kw2...)
(fl::Flow)(::Val{:Full}, x, p, t) 	  = fl.flowFull(x, p, t)
(fl::Flow)(::Val{:TimeSol}, x, p, t)  = fl.flowTimeSol(x, p, t)
(fl::Flow)(::Val{:SerialdFlow}, x, p, dx, t) = fl.dfserial(x, p, dx, t)

"""
Creates a Flow variable based on a `prob::ODEProblem` and ODE solver `alg`. The vector field `F` has to be passed, this will be resolved in the future as it can be recovered from `prob`. Also, the derivative of the flow is estimated with finite differences.
"""
# this constructor takes into accound a parameter passed to the vector field
function Flow(F, p, prob::Union{ODEProblem, EnsembleProblem}, alg; kwargs...)
	probserial = prob isa EnsembleProblem ? prob.prob : prob
	return Flow(F,
		(x, p, t) ->		   flow(x, p, t, prob; alg = alg, kwargs...),
		(x, p, t) ->	flowTimeSol(x, p, t, prob; alg = alg, kwargs...),
		(x, p, t) ->       flowFull(x, p, t, prob; alg = alg, kwargs...),
		# we remove the callback in order to use this for the Jacobian in Poincare Shooting
		(x, p, dx, t; kw2...) -> dflow_fd(x, p, dx, t, prob; alg = alg, kwargs..., kw2...),
		# serial version of dflow. Used for the computation of Floquet coefficients
		(x, p, dx, t; kw2...) -> dflow_fd(x, p, dx, t, probserial; alg = alg, kwargs..., kw2...)
		)
end

function Flow(F, p, prob1::Union{ODEProblem, EnsembleProblem}, alg1, prob2::Union{ODEProblem, EnsembleProblem}, alg2; kwargs...)
	probserial1 = prob1 isa EnsembleProblem ? prob1.prob : prob1
	probserial2 = prob2 isa EnsembleProblem ? prob2.prob : prob2
	return Flow(F,
		(x, p, t) ->			  flow(x, p, t, prob1, alg = alg1; kwargs...),
		(x, p, t) ->	   flowTimeSol(x, p, t, prob1; alg = alg1, kwargs...),
		(x, p, t) -> 		  flowFull(x, p, t, prob1, alg = alg1; kwargs...),
		(x, p, dx, t; kw2...) -> dflow(x, p, dx, t, prob2; alg = alg2, kwargs..., kw2...),
		# serial version of dflow. Used for the computation of Floquet coefficients
		(x, p, dx, t; kw2...) -> dflow(x, p, dx, t, probserial2; alg = alg2, kwargs..., kw2...),
		)
end


text:
#=

Draw radar map in resulting world frame (after solve)

=#

##

using Distributed
# addprocs(10)

using Colors
using Caesar
@everywhere using Caesar

import Rotations as _Rot

# using DistributedFactorGraphs.DocStringExtensions
using Dates
using JSON2
# using BSON
using Serialization
using FixedPointNumbers
using Manifolds
using StaticArrays
using ImageMagick, FileIO
using Images

using ImageDraw
using ImageView
using ImageTransformations

using RoMEPlotting
Gadfly.set_default_plot_size(35cm,20cm)

## load the solved graph and reattach the datastores

dfg_datafolder = "/tmp/caesar/philos"

fg = loadDFG("$dfg_datafolder/results_4/x0_5_x280")

##

ds = FolderStore{Vector{UInt8}}(:radar, "$dfg_datafolder/data/radar")
addBlobStore!(fg, ds)

ds = FolderStore{Vector{UInt8}}(:gps_fix, "$dfg_datafolder/data/gps")
addBlobStore!(fg, ds)

# add if you want lidar also
ds = FolderStore{Vector{UInt8}}(:lidar, "$dfg_datafolder/data/lidar")
addBlobStore!(fg, ds)

# add if you want lidar also
ds = FolderStore{Vector{UInt8}}(:camera, "$dfg_datafolder/data/camera")
addBlobStore!(fg, ds)

## visualize the radar data

if false
  imgs = map(x->Gray{N0f8}.(x), fetchDataImage.(fg, sortDFG(ls(fg)), :RADAR_IMG));
  writevideo("/tmp/caesar/philos/radar280.mp4", imgs; fps=5, player="totem")
end

##

M = SpecialEuclidean(2)
e0 = identity_element(M)

##

# get two radar sweep point clouds and rotate them to the world frame
function _fetchDataPointCloudWorld(fg::AbstractDFG, lb::Symbol)
  de,db = getData(fg, lb, :RADAR_SWEEP)
  # deserialize is not legacy safe, rather use BSON.@save(io, Caesar._PCL.PCLPointCloud2(pc))
  bPC = deserialize(PipeBuffer(db))
  wPb = exp(M, e0, hat(M, e0, getPPE(fg, lb).suggested))
  apply(M, wPb, bPC)
end

##

wPC0 = _fetchDataPointCloudWorld(fg, :x0)
img0 = makeImage!(wPC0, (-1500,1500); color=Gray{N0f8}(0.01))
imshow(img0)

wPC50 = _fetchDataPointCloudWorld(fg, :x50)
# makeImage!(wPC50, (-1500,1500)) |> imshow

wPC100 = _fetchDataPointCloudWorld(fg, :x100)
# makeImage!(wPC100, (-2000,2000)) |> imshow

##

_X_ = 3000
_Y_ = 1250

wPC__ = _fetchDataPointCloudWorld.(fg, [Symbol(:x,i) for i in 0:5:280]);
img_wPC_ = map(pc->makeImage!(pc, (-1250,_X_),(-_Y_,_Y_); color=Gray{Float64}(0.05)), wPC__);
img_wPC = +(img_wPC_...);

imshow(img_wPC)

## movie sequence in world frame

IMGS = Vector{Matrix{RGB{Float64}}}()

img_wPC = RGB.(img_wPC_[1])
for gim in img_wPC_[2:end]
  # green overlay of latest
  imgG = (px->RGB(0,10*px,0)).(gim);
  imgL = img_wPC + imgG
  push!(IMGS, imgL)

  # grow full world map
  img_wPC += RGB.(gim)
end

writevideo("/tmp/caesar/philos/results_4/radarmap.mp4", IMGS; fps=3, player="totem")

## ================================================================================
## Keeping a few code fragments below
## ================================================================================

# for l in [Symbol(:x,i) for i in 10:10:280]
#   wPC = cat(wPC, _fetchDataPointCloudWorld(fg, l); reuse=true )
# end
# img_wPC = makeImage!(wPC, (-1250,_X_),(-_Y_,_Y_))
# imshow(img_wPC)

## manually north align

# nPC = apply(M, ProductRepr([0.0;0.0], _Rot.RotMatrix(135*pi/180.0)), wPC)

# img_nPC = makeImage!(nPC, (-2500,1000),(-2000,1500), color=Gray{N0f8}(0.1))
# img_nPC[1,1] = Gray(1)
# imshow(img_nPC)

# ##

# using ImageMagick, FileIO
# pngimg = toFormat(format"PNG", img_wPC)
# save(dfg_datafolder*"/results_3/wPC.png", img_wPC)

# ##

#


text:
abstract type StochasticProcess end

# process

# length of state X
function stateSize(self::StochasticProcess)
    throw(ArgumentError("Implementation of method stateSize() required."))
end

# lengh of stochastic factors, dW
function factors(self::StochasticProcess)
    throw(ArgumentError("Implementation of method factors() required."))
end

# initial values for simulation
function initialValues(self::StochasticProcess)
    throw(ArgumentError("Implementation of method initialValues() required."))
end

# evolve X(t0) -> X(t0+dt) using independent Brownian increments dW
# t0, dt are assumed float, X0, X1, dW are array
function evolve(self::StochasticProcess, t0, X0, dt, dW, X1)
    throw(ArgumentError("Implementation of method evolve() required."))
end

# hybrid evolution

# the short rate over an integration time period
# this is required for drift calculation in multi-asset and hybrid models
function shortRateOverPeriod(self::StochasticProcess, t0, dt, X0, X1)
    throw(ArgumentError("Implementation of method shortRateOverPeriod() required."))
end

function zeroBondVolatility(self::StochasticProcess, t, T)
    throw(ArgumentError("Implementation of method zeroBondVolatility() required."))
end

function zeroBondVolatilityPrime(self::StochasticProcess, t, T)
    throw(ArgumentError("Implementation of method zeroBondVolatilityPrime() required."))
end

# model

# the numeraire in the domestic currency used for discounting future payoffs
function numeraire(self::StochasticProcess, t, X)
    throw(ArgumentError("Implementation of method numeraire() required."))
end

# an asset price for a given currency alias
function asset(self::StochasticProcess, t, X, alias)
    throw(ArgumentError("Implementation of method asset() required."))
end

# a domestic/foreign currency zero coupon bond
function zeroBond(self::StochasticProcess, t, T, X, alias)
    throw(ArgumentError("Implementation of method zeroBond() required."))
end

# additional simulated credit quantities, we use a Cox process
# for references, see
#   - Brigo/Mercurio, 2007, Sec. 22.2.3
#   - Brigo/Vrins, Disentangling wrong-way risk, 2016, Sec. 3.3

# Cumulated intensity; probability of tau > t, conditional on F_t
function hazardProcess(self::StochasticProcess, t, X, alias)
    throw(ArgumentError("Implementation of method hazardProcess() required."))
end

# instantanous probability of default
function hazardRate(self::StochasticProcess, t, X, alias)
    throw(ArgumentError("Implementation of method hazardRate() required."))
end

# probavility of survival consitional on information at t
function survivalProb(self::StochasticProcess, t, T, X, alias)
    throw(ArgumentError("Implementation of method survivalProb() required."))
end

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
using Base.Test
# fourier1dfilter.jl
import WaveletScattering: AbstractFilter, Analytic1DFilter,
    Coanalytic1DFilter, FourierSymmetric1DFilter,
    FullResolution1DFilter,
    Vanishing1DFilter, VanishingWithMidpoint1DFilter,
    critical_log2_sampling, getindex, littlewoodpaleyadd!,
    nextpow2_exponent, renormalize!, spin
# meta.jl
import WaveletScattering: ΨMeta1D, ΦMeta, get_bandwidth,
    get_centerfrequency, get_j, get_n_orientations,
    get_qualityfactor, get_scale, get_γ, get_θ, get_χ
# spec1d.jl
import WaveletScattering: Spec1D

# AbstractFourierFilter 1D constructor
spec = Spec1D(log2_size = 4)
N = 16
y = zeros(Float32, N); y[16] = 1.0
@test isa(AbstractFilter(y, spec), Analytic1DFilter{Float32})
y = zeros(Float32, N); y[2] = 1.0
@test isa(AbstractFilter(y, spec), Coanalytic1DFilter{Float32})
y = ones(Float32, N)
@test isa(AbstractFilter(y, spec), FullResolution1DFilter{Float32})
y = zeros(Float32, N); y[2] = 1.0; y[16] = 1.0
@test isa(AbstractFilter(y, spec), Vanishing1DFilter{Float32})
y = zeros(Float32, N); y[2] = 1.0; y[9] = 1.0; y[1] = 1.0
@test isa(AbstractFilter(y, spec), VanishingWithMidpoint1DFilter{Float32})

# multiplication operator with scalar
# Base.:*{T<:Number}(ψ::Analytic1DFilter{T}, b::Number)
ψ = Analytic1DFilter(Float32[0.1, 0.3], 2)
ψ2 = ψ * Float32(2.0)
@test isa(ψ2, Analytic1DFilter{Float32})
@test ψ2.pos ≈ Float32[0.2, 0.6]
@test ψ2.posfirst == 2
# Base.:*{T<:Number}(ψ::Coanalytic1DFilter{T}, b::Number)
ψ = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
ψ2 = ψ * Float32(2.0)
@test isa(ψ2, Coanalytic1DFilter{Float32})
@test ψ2.neg ≈ Float32[0.2, 0.6, 0.8]
@test ψ2.neglast == -3
# Base.:*{T<:Number}(ψ::FullResolution1DFilter{T}, b::Number)
ψ = FullResolution1DFilter(Float32[0.01, 0.1, 0.2, 0.3])
ψ2 = ψ * Float32(2.0)
@test isa(ψ2, FullResolution1DFilter{Float32})
@test ψ2.coeff ≈ Float32[0.02, 0.2, 0.4, 0.6]
# Base.:*{T<:Number}(ψ::FourierSymmetric1DFilter{T}, b::Number)
ψ = FourierSymmetric1DFilter(Float32[0.8, 0.3, 0.1], one(Float32))
ψ2 = ψ * Float32(2.0)
@test isa(ψ2, FourierSymmetric1DFilter{Float32})
@test ψ2.leg ≈ Float32[1.6, 0.6, 0.2]
@test ψ2.zero ≈ Float32(2.0)
# Base.:*{T<:Number}(ψ::Vanishing1DFilter{T}, b::Number)
an = Analytic1DFilter(Float32[0.1, 0.3], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
ψ = Vanishing1DFilter(an, coan)
ψ2 = ψ * Float32(2.0)
@test isa(ψ2, Vanishing1DFilter{Float32})
@test ψ2.an.pos ≈ Float32[0.2, 0.6]
@test ψ2.an.posfirst == 2
@test ψ2.coan.neg ≈ Float32[0.2, 0.6, 0.8]
@test ψ2.coan.neglast == -3
# Base.:*{T<:Number}(ψ::VanishingWithMidpoint1DFilter{T}, b::Number))
an = Analytic1DFilter(Float32[0.1, 0.3], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
midpoint = Float32(0.5)
ψ = VanishingWithMidpoint1DFilter(an, coan, midpoint)
ψ2 = ψ * Float32(2.0)
@test isa(ψ2, VanishingWithMidpoint1DFilter{Float32})
@test ψ2.an.pos ≈ Float32[0.2, 0.6]
@test ψ2.an.posfirst == 2
@test ψ2.coan.neg ≈ Float32[0.2, 0.6, 0.8]
@test ψ2.coan.neglast == -3
@test ψ2.midpoint ≈ Float32(1.0)

# getindex
# getindex{T}(ψ::Analytic1DFilter{T}, i::Integer)
ψ = Analytic1DFilter(Float32[0.1, 0.3], 2)
@test Float32[ ψ[ω] for ω in 1:4 ] == Float32[0.0, 0.1, 0.3, 0.0]
# getindex{T}(ψ::Analytic1DFilter{T}, I::UnitRange{Int64})
@test ψ[-1:1] == Float32[0.0 ; 0.0 ; 0.0]
@test ψ[1:4] == Float32[0.0 ; 0.1 ; 0.3; 0.0]
@test ψ[2:3] == Float32[0.1 ; 0.3]
@test ψ[1:3] == Float32[0.0 ; 0.1 ; 0.3]
@test ψ[2:4] == Float32[0.1 ; 0.3 ; 0.0]
@test ψ[5:7] == Float32[0.0 ; 0.0 ; 0.0]
# getindex{T}(ψ::Coanalytic1DFilter{T}, i::Integer)
ψ = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
@test Float32[ψ[ω] for ω in -6:-2] == Float32[0.0, 0.1, 0.3, 0.4, 0.0]
# getindex{T}(ψ::Coanalytic1DFilter{T}, I::UnitRange{Int64})
@test ψ[-3:-1] == Float32[0.4 ; 0.0 ; 0.0]
@test ψ[-4:-3] == Float32[0.3 ; 0.4]
@test ψ[-5:-4] == Float32[0.1 ; 0.3]
@test ψ[-7:-5] == Float32[0.0 ; 0.0 ; 0.1]
@test ψ[-10:-8] == Float32[0.0 ; 0.0 ; 0.0]
@test ψ[2:3] == Float32[0.0 ; 0.0]
@test ψ[-6:-2] == Float32[0.0 ; 0.1 ; 0.3 ; 0.4 ; 0.0]
# getindex{T}(ψ::FullResolution1DFilter{T}, i::Integer)
ψ = FullResolution1DFilter(Float32[0.01, 0.1, 0.2, 0.3])
@test Float32[ψ[ω] for ω in -3:2] == Float32[0.0, 0.2, 0.3, 0.01, 0.1, 0.0]
# getindex{T}(ψ::FullResolution1DFilter{T}, I::UnitRange{Int64})
@test ψ[-3:2] == Float32[0.0, 0.2, 0.3, 0.01, 0.1, 0.0]
@test ψ[-5:-3] == Float32[0.0, 0.0, 0.0]
@test ψ[3:5] == Float32[0.0, 0.0, 0.0]
@test ψ[5:3] == Float32[]
# getindex{T}(ψ::Vanishing1DFilter{T}, i::Integer)
an = Analytic1DFilter(Float32[0.1, 0.3], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
ψ = Vanishing1DFilter(an, coan)
@test [ψ[ω] for ω in -6:4] ==
    Float32[0.0, 0.1, 0.3, 0.4, 0.0, 0.0, 0.0, 0.0, 0.1, 0.3, 0.0]
# getindex{T}(ψ::Vanishing1DFilter{T}, I::UnitRange{Int64})
@test ψ[-6:4] == Float32[0.0, 0.1, 0.3, 0.4, 0.0, 0.0, 0.0, 0.0, 0.1, 0.3, 0.0]
# getindex{T}(ψ::VanishingWithMidpoint1DFilter{T}, I::Integer)
an = Analytic1DFilter(Float32[0.1, 0.3, 0.4], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -2)
midpoint = Float32(0.5)
ψ = VanishingWithMidpoint1DFilter(an, coan, midpoint)
@test Float32[ψ[ω] for ω in -6:5] ==
    Float32[0.0, 0.5, 0.1, 0.3, 0.4, 0.0, 0.0, 0.0, 0.1, 0.3, 0.4, 0.0]
# getindex{T}(ψ::VanishingWithMidpoint1DFilter{T}, I::UnitRange{Int64})
@test ψ[-6:5] ==
    Float32[0.0, 0.5, 0.1, 0.3, 0.4, 0.0, 0.0, 0.0, 0.1, 0.3, 0.4, 0.0]

# critical_log2_sampling
an = Analytic1DFilter(Float32[0.1, 0.3, 0.4], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -2)
midpoint = Float32(0.5)
spec = Spec1D()
ψ = VanishingWithMidpoint1DFilter(an, coan, midpoint)
@test critical_log2_sampling(an, spec) == -12
@test critical_log2_sampling(coan, spec) == -16
@test critical_log2_sampling(ψ, spec) == 0

# littlewoodpaleyadd!
# littlewoodpaleyadd!(lp::Vector, ψ::Analytic1DFilter)
lp = zeros(Float32, 8)
ψ = Analytic1DFilter(Float32[0.1, 0.3], 2)
littlewoodpaleyadd!(lp, ψ); lp = zeros(Float32, 8) # warmup
allocatedmemory = @allocated littlewoodpaleyadd!(lp, ψ)
@test allocatedmemory <= 1e3 # on some machines (e.g. Travis's Linux) it is >0
@test lp ≈ Float32[0.0, 0.0, 0.01, 0.09, 0.0, 0.0, 0.0, 0.0]
# littlewoodpaleyadd!(lp::Vector, ψ::Coanalytic1DFilter)
lp = zeros(Float32, 8)
ψ = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
littlewoodpaleyadd!(lp, ψ); lp = zeros(Float32, 8) # warmup
allocatedmemory = @allocated littlewoodpaleyadd!(lp, ψ)
@test allocatedmemory <= 1e3 # on some machines (e.g. Travis's Linux) it is >0
@test lp ≈ Float32[0.0, 0.0, 0.0, 0.01, 0.09, 0.16, 0.0, 0.0]
# littlewoodpaleyadd!(lp::Vector, ψ::FourierSymmetric1DFilter)
lp = zeros(Float32, 8)
ψ = FourierSymmetric1DFilter(Float32[0.8, 0.3, 0.1], one(Float32))
littlewoodpaleyadd!(lp, ψ); lp = zeros(Float32, 8) # warmup
allocatedmemory = @allocated littlewoodpaleyadd!(lp, ψ)
@test allocatedmemory <= 1e3 # on some machines (e.g. Travis's Linux) it is >0
@test lp ≈ Float32[1.0, 0.64, 0.09, 0.01, 0.0, 0.01, 0.09, 0.64]
# littlewoodpaleyadd!(lp::Vector, ψ::Vanishing1DFilter)
lp = zeros(Float32, 8)
an = Analytic1DFilter(Float32[0.1, 0.3], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -1)
ψ = Vanishing1DFilter(an, coan)
littlewoodpaleyadd!(lp, ψ); lp = zeros(Float32, 8) # warmup
allocatedmemory = @allocated littlewoodpaleyadd!(lp, ψ)
@test allocatedmemory <= 1e3 # on some machines (e.g. Travis's Linux) it is >0
@test lp ≈ Float32[0.0, 0.0, 0.01, 0.09, 0.0, 0.01, 0.09, 0.16]
# littlewoodpaleyadd!(lp::Vector, ψ::VanishingWithMidpoint1DFilter)
lp = zeros(Float32, 8)
an = Analytic1DFilter(Float32[0.1, 0.3], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -1)
midpoint = Float32(0.5)
ψ = VanishingWithMidpoint1DFilter(an, coan, midpoint)
littlewoodpaleyadd!(lp, ψ); lp = zeros(Float32, 8) # warmup
allocatedmemory = @allocated littlewoodpaleyadd!(lp, ψ)
@test allocatedmemory <= 1e3 # on some machines (e.g. Travis's Linux) it is >0
@test lp ≈ Float32[0.0, 0.0, 0.01, 0.09, 0.25, 0.01, 0.09, 0.16]

# Base.maximum
# Base.maximum(ψ::Analytic1DFilter)
@test maximum(Analytic1DFilter([0.1, -0.3], 2)) ≈ 0.3
# Base.maximum(ψ::Coanalytic1DFilter)
@test maximum(Coanalytic1DFilter([0.1, 0.3, 0.4*im], -3)) ≈ 0.4
# Base.maximum(ψ::FullResolution1DFilter)
@test maximum(FullResolution1DFilter([0.01, 0.1*im, 0.2, -0.3])) ≈ 0.3
# Base.maximum(ψ::FourierSymmetric1DFilter)
@test maximum(FourierSymmetric1DFilter([0.2, 0.1], 1.0)) ≈ 1.0
# Base.maximum(ψ::Vanishing1DFilter)
an = Analytic1DFilter(Float32[0.1, 0.3], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
@test maximum(Vanishing1DFilter(an, coan)) ≈ Float32(0.4)
# Base.maximum(ψ::VanishingWithMidpoint1DFilter)
an = Analytic1DFilter(Float32[0.1, 0.3], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
midpoint = Float32(0.5)
ψ = VanishingWithMidpoint1DFilter(an, coan, midpoint)
@test maximum(ψ) ≈ Float32(0.5)

# nextpow2_exponent
@test nextpow2_exponent(5) == 3
@test nextpow2_exponent(4) == 2
@test nextpow2_exponent(3) == 2
@test nextpow2_exponent(2) == 1
@test nextpow2_exponent(1) == 0
@test nextpow2_exponent(0) == 0
@test nextpow2_exponent(-1) == 0
@test nextpow2_exponent(-2) == -1

# renormalize!
# case Q=1
spec = Spec1D()
(nΘs, nΧs, nJs) = size(spec.ψmetas)
D = typeof(spec.domain)
ψs = Array(AbstractFilter{spec.signaltype,D}, (nΘs, nΧs, nJs))
ψs[1, :, :] =
    pmap(AbstractFilter, spec.ψmetas[1, :, :], fill(spec, nΧs, nJs))
(nΘs > 1) && spin!(ψs)
ϕ = AbstractFilter(spec.ϕmeta, spec)
lp = renormalize!(ϕ, ψs, spec)
@test all(lp.< 1.01)
N = 1 << spec.log2_size
ξs = map(get_centerfrequency, spec.ψmetas)
firstω = round(Int, N * ξs[end])
lastω = round(Int, N * ξs[1])
@test all(lp[1+(firstω:lastω)] .> 0.5)
# case Q>1, max_s = Inf
spec = Spec1D(n_filters_per_octave = 8)
(nΘs, nΧs, nJs) = size(spec.ψmetas)
D = typeof(spec.domain)
ψs = Array(AbstractFilter{spec.signaltype,D}, (nΘs, nΧs, nJs))
ψs[1, :, :] =
    pmap(AbstractFilter, spec.ψmetas[1, :, :], fill(spec, nΧs, nJs))
(nΘs > 1) && spin!(ψs)
ϕ = AbstractFilter(spec.ϕmeta, spec)
lp = renormalize!(ϕ, ψs, spec)
@test all(lp.< 1.01)
N = 1 << spec.log2_size
ξs = map(get_centerfrequency, spec.ψmetas)
firstω = round(Int, N * ξs[end])
lastω = round(Int, N * ξs[1])
@test all(lp[1+(firstω:lastω)] .> 0.5)
# # case Q>1, max_s < Inf
spec = Spec1D(max_scale = 4410, n_filters_per_octave = 8)
(nΘs, nΧs, nJs) = size(spec.ψmetas)
D = typeof(spec.domain)
ψs = Array(AbstractFilter{spec.signaltype,D}, (nΘs, nΧs, nJs))
ψs[1, :, :] =
    pmap(AbstractFilter, spec.ψmetas[1, :, :], fill(spec, nΧs, nJs))
(nΘs > 1) && spin!(ψs)
ϕ = AbstractFilter(spec.ϕmeta, spec)
lp = renormalize!(ϕ, ψs, spec)
#@test all(lp.< 1.01)
N = 1 << spec.log2_size
ξs = map(get_centerfrequency, spec.ψmetas)
firstω = round(Int, N * ξs[end])
lastω = round(Int, N * ξs[1])
# @test all(lp[1+(firstω:lastω)] .> 0.5)

# spin
# spin(::Analytic1DFilter)
ψ = Analytic1DFilter(Float32[0.1, 0.3], 2)
ψspinned = spin(ψ)
@test isa(ψspinned, Coanalytic1DFilter{Float32})
@test ψspinned.neg ≈ Float32[0.3, 0.1]
@test ψspinned.neglast == -2
# spin(::Coanalytic1DFilter)
ψ = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
ψspinned = spin(ψ)
@test isa(ψspinned, Analytic1DFilter{Float32})
@test ψspinned.pos ≈ Float32[0.4, 0.3, 0.1]
@test ψspinned.posfirst == 3
# spin(::FullResolution1DFilter)
ψ = FullResolution1DFilter(Float32[0.1, 0.2, 0.3, 0.4])
ψspinned = spin(ψ)
@test isa(ψspinned, FullResolution1DFilter{Float32})
@test ψspinned.coeff ≈ Float32[0.4, 0.3, 0.2, 0.1]
# spin(::Vanishing1DFilter)
an = Analytic1DFilter(Float32[0.1, 0.3], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
ψ = Vanishing1DFilter(an, coan)
ψspinned = spin(ψ)
@test isa(ψspinned, Vanishing1DFilter{Float32})
@test ψspinned.coan.neg ≈ Float32[0.3, 0.1]
@test ψspinned.coan.neglast == -2
@test ψspinned.an.pos ≈ Float32[0.4, 0.3, 0.1]
@test ψspinned.an.posfirst == 3
# spin(::VanishingWithMidpoint1DFilter)
an = Analytic1DFilter(Float32[0.1, 0.3], 2)
coan = Coanalytic1DFilter(Float32[0.1, 0.3, 0.4], -3)
midpoint = Float32(0.5)
ψ = VanishingWithMidpoint1DFilter(an, coan, midpoint)
ψspinned = spin(ψ)
@test isa(ψspinned, VanishingWithMidpoint1DFilter{Float32})
@test ψspinned.coan.neg ≈ Float32[0.3, 0.1]
@test ψspinned.coan.neglast == -2
@test ψspinned.an.pos ≈ Float32[0.4, 0.3, 0.1]
@test ψspinned.an.posfirst == 3
@test ψspinned.midpoint ≈ Float32(0.5)


text:
module PredictMDFuture

end # module


text:
using Base.Test
import Images

@test Images.minclamp(-7, 10.0) === 10
@test Images.minclamp(-7, -10.0) === -7
@test Images.maxclamp(-7, 10.0) === -7
@test Images.maxclamp(-7, -10.0) === -10

@test Images.clamp(Uint8, -5) === 0x00
@test Images.clamp(Uint8, 5) === 0x05
@test Images.clamp(Uint8, 2000) === 0xff
@test Images.clamp(Float32, -3.2) === -3.2f0

@mytest_throws InexactError Images.clamp(Uint16, 3.2)
@test Images.truncround(Uint16, 3.2) === uint16(3)
@test Images.truncround(Uint16, -3.2) == uint16(0)
@test Images.truncround(Uint16, 1e20) == typemax(Uint16)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "UASM"
# Note: upstream doesn't tag new versions, but creates multiple versioned
# branches they keep pushing all the time, so it isn't even clear when a version
# is "stable".  Yes, you read it right.
version = v"2.55"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/Terraspace/UASM",
              "43bf08acc641166ee4766a33a5accc6e84535b3a"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/UASM*/
make -f gccLinux64.mak
install -Dvm 0755 "GccUnixR/uasm" "${bindir}/uasm${exeext}"
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    # The build system is... uhm... complicated.  It's pointless to fight it for
    # multiple platforms, we need this mainly as a native tool when building.
    Platform("x86_64", "linux"; libc="musl"),
]

# The products that we will ensure are always built
products = [
    ExecutableProduct("uasm", :uasm),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"5")

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
export DiracDelta

for TYP in (:DiracSpace,:PointSpace)
    @eval begin
        immutable $TYP{T}<:RealUnivariateSpace{AnyDomain}
          points::Vector{T}
          $TYP(pts::Vector{T})=new(sort(pts))
        end

        $TYP(points::AbstractVector) = $TYP{eltype(points)}(points)
        $TYP() = $TYP(Float64[])
        $TYP(point::Number) = $TYP([point])
        $TYP(p::Point)=$TYP(p.x)
        dimension(d::$TYP)=length(d.points)

        domain(DS::$TYP)=mapreduce(Point,union,DS.points)
        setdomain(DS::$TYP,d::UnionDomain)=$TYP(map(d->d.x,d.domains))
        points(sp::$TYP,n::Integer)=sp.points[1:n]

        spacescompatible(a::$TYP,b::$TYP)=a.points==b.points
        canonicalspace(a::$TYP)=a

        union_rule(a::$TYP,b::$TYP)=$TYP(sort(union(a.points,b.points)))

        function coefficients(cfs::Vector,fromspace::$TYP,tospace::$TYP)
            if spacescompatible(fromspace,tospace)
                return cfs
            end

            @assert length(cfs) ≤ length(fromspace.points)

            # this first for-loop removes coefficients of Dirac points that are zero
            nonzerofromspacepoints = eltype(fromspace.points)[]
            nonzeroDiraccfs = eltype(cfs)[]
            for i = 1:length(cfs)
                if cfs[i] != 0
                    push!(nonzerofromspacepoints, fromspace.points[i])
                    push!(nonzeroDiraccfs, cfs[i])
                end
            end

            # if the points that remain can be represented in the tospace
            if issubset(nonzerofromspacepoints,tospace.points)
                finalDiraccfs = eltype(cfs)[]
                j=1 #counter for the nonzerofromspacepoints
                for i = 1:length(tospace.points)
                    if j > length(nonzerofromspacepoints)
                        break
                    elseif nonzerofromspacepoints[j]==tospace.points[i]
                        push!(finalDiraccfs,nonzeroDiraccfs[j])
                        j += 1
                    else
                        push!(finalDiraccfs,0)
                    end
                end
                finalDiraccfs
            else
                error("The domain of the space being converted from has points that cannot be represented in the space you are converting to.")
            end
        end
    end
end

Space(d::Point)=PointSpace(d)

identity_fun(S::PointSpace)=Fun(S.points,S)
identity_fun(S::DiracSpace)=Fun(S.points,PointSpace(S.points))
transform(S::PointSpace,v::Vector,plan...)=v

function evaluate(f::AbstractVector,PS::PointSpace,x::Number)
    p = findfirst(y->isapprox(x,y),PS.points)
    if p==0
        zero(eltype(f))
    else
        f[p]
    end
end

evaluate(f::AbstractVector,PS::PointSpace,x::AbstractVector)=
    map(y->evaluate(f,PS,y),x)

Base.sum{DS<:DiracSpace}(f::Fun{DS})=sum(f.coefficients[1:dimension(space(f))])

DiracDelta(x::Number)=Fun([1.],DiracSpace(x))
DiracDelta()=DiracDelta(0.)

function Base.cumsum{S<:DiracSpace,T<:Real}(f::Fun{S},d::Interval{T})
    pts=space(f).points
    @assert pts ==sort(pts)
    cfs=cumsum(f.coefficients)
    if first(d) < first(pts) && last(d) > last(pts)
        Fun([0;cfs],HeavisideSpace([first(d);pts;last(d)]))
    elseif first(d) == first(pts) && last(d) > last(pts)
        Fun(cfs,HeavisideSpace([pts;last(d)]))
    elseif first(d) < first(pts) && last(d) == last(pts)
        Fun([0;cfs],HeavisideSpace([first(d);pts]))
    elseif first(d) == first(pts) && last(d) == last(pts)
        Fun(cfs[1:end-1],HeavisideSpace(pts))
    else
        error("Implement")
    end
end

# for TYP in (:ReSpace,:Space)
#   @eval begin
#     function coefficients(cfs::Vector,fromspace::$TYP,tospace::DiracSpace)
#       [0*tospace.points;coefficients(cfs,fromspace,tospace.space)]
#     end
#   end
# end

# function coefficients(cfs::Vector,fromspace::DiracSpace,tospace::Space)
#   n = length(fromspace.points)
#   if n == 0 || cfs[1:n] == 0*cfs[1:n]
#       coefficients(fromspace.space,tospace)
#   else
#     error("The space you are converting from has Dirac deltas that cannot be represented in the space you are converting to.")
#   end
# end

# function Multiplication{PS<:PointSpace}(f::Fun{PS},PS2::PointSpace)
#     @assert space(f).points==PS2.points
#     FiniteOperator(diagm(f.coefficients),PS2,PS2)
# end

# function Multiplication{PS<:PointSpace}(f::Fun{PS},DS::DiracSpace)
#     @assert space(f).points==DS.points
#     FiniteOperator(diagm(f.coefficients),DS,DS)
# end

# function Multiplication{DS<:DiracSpace}(f::Fun{DS},PS::PointSpace)
#     @assert space(f).points==PS.points
#     FiniteOperator(diagm(f.coefficients),PS,space(f))
# end

function coefficienttimes{PS<:PointSpace,DS<:DiracSpace}(f::Fun{PS},g::Fun{DS})
    @assert space(f).points==space(g).points
    Fun(f.coefficients.*g.coefficients,space(g))
end

function coefficienttimes{PS<:PointSpace,DS<:DiracSpace}(f::Fun{DS},g::Fun{PS})
    @assert space(f).points==space(g).points
    Fun(f.coefficients.*g.coefficients,space(f))
end

function coefficienttimes{PS<:PointSpace,PS2<:PointSpace}(f::Fun{PS},g::Fun{PS2})
    @assert space(f).points==space(g).points
    Fun(f.coefficients.*g.coefficients,space(g))
end

./{PS<:PointSpace}(f::Fun,g::Fun{PS}) = f.*inv(g)
Base.inv{PS<:PointSpace}(f::Fun{PS}) = Fun(1./f.coefficients,space(f))


text:
using GeomDAE
using Base.Test

function fx(t, x, fx)
    fx .= x
end

function fq(t, q, p, fq)
    fq .= q
end

function fp(t, q, p, fp)
    fp .= q.^2
end

include("equations/equations_tests.jl")
include("integrators/timeseries_tests.jl")
include("integrators/solutions_tests.jl")
include("integrators/integrators_test.jl")
include("integrators/tableaus_tests.jl")
include("solvers/linear_solvers_tests.jl")
include("solvers/nonlinear_solvers_tests.jl")
#include("utils/hdf5_tests.jl")


text:
include("/home/alkan/c/nese/jtb/yeast_dyn.jl")

# --- switch edges, keeping incoming and outgoing degrees ---
function switchnet (intr)
    n = length(intr[:,1])
    el = (Int64,Int64,Int64)[] # empty edge list
    for i=1:n for j=1:n
        if intr[i,j] != 0
            el = [el,(i,j,intr[i,j])]
        end
    end end

    ## shuffle edgelist
    ne = length(el)
    for rep=1:ne
        switching = true
        while switching
            r = rand(1:ne)              # pick a random edge
            rr = rand(1:ne)
            while r==rr rr=rand(1:ne) end # make sure edges are different
            if r<rr
                (r,rr) = (rr,r) # critical (otherwise idxs shift during splice)
            end
            if el[r][3]==el[rr][3]    # forget unless both + or both -
                if el[r][1]!=el[rr][2] && el[r][2]!=el[rr][1] &&
                    el[r][1]!=el[rr][1] && el[r][2]!=el[rr][2] # 4 different nodes
                    e1=(el[r][1],el[rr][2],el[r][3]) # new edges
                    e2=(el[rr][1],el[r][2],el[r][3])
                    e3=(e1[1],e1[2],-e1[3]) # new edges
                    e4=(e2[1],e2[2],-e2[3]) # new edges
                    if !in(e1,el) && !in(e2,el) && !in(e3,el) && !in(e4,el)
                        splice!(el,r)
                        splice!(el,rr)
                        append!(el,[e1,e2])
                        switching = false
                    end end end         # ifs
        end                             # while
    end                                 # for
    nw = zeros(Int64,n,n)
    for i=1:ne                          # construct new nw
        nw[el[i][1],el[i][2]] = el[i][3]
    end
return nw
end

# --- return a (n,ne,np) preserved shuffled version of a network
function shufflenet (intr,slf)
    n = length(slf)
    nw = similar(intr)                  # this will be the randomized version
    shf = shuffle(slf)                  # keep number of self degredating genes
    tl = tril(intr)                     # remove diag and join upper/lower triang
    tu = triu(intr)
    rect = tu[:,2:end]+tl[:,1:end-1]    # shuffle this
    rect = reshape(shuffle(reshape(rect,length(rect),1)[:,1]),n,n-1)
    rrect = [rect zeros(Int64,n)]       # now insert back the shuffled diagonal
    lrect = [zeros(Int64,n) rect]
    nw = tril(rrect - diagm(diag(rrect))) + triu(lrect-diagm(diag(lrect)))
    return (nw,shf)
end

# --- generate a random (n,ne,np) network.
# --- d = matrix diagonal
# --- ne = number of edges
# --- np = number of positive edges
function randomnkp (d,ne,np)
    n = length(d)                       # int matrix diagonal - never change
    nn = ne-np                          # negative edges
    nz = n*(n-1)-ne                     # number of zero-edges
#    print("in randomnkp() with $ne, $np")
    intlist = [ones(Int64,np);-ones(Int64,nn);zeros(Int64,nz)]
    rect = reshape(shuffle(intlist),n,n-1)
    rrect = [rect zeros(Int64,n)]
    lrect = [zeros(Int64,n) rect]
    interrand = tril(rrect - diagm(diag(rrect))) + triu(lrect-diagm(diag(lrect)))
#    println("done reshaping.")
    return interrand
end

max_stars_repo_name: akabakcioglu/network-dynamics

max_issues_repo_name: akabakcioglu/network-dynamics

max_forks_repo_name: akabakcioglu/network-dynamics


text:
module GodRays

using ..Ahorn, Maple

placements = Maple.GodRays

function Ahorn.canFgBg(effect::Maple.GodRays)
    return true, true
end

end


text:
using Test
using Random
using Distributions
using RandomTree

Random.seed!(0)

function treetest(tree)
    @test typeof(tree) <: FiniteTree
    @test isa(tree, FiniteTree)

    mdeg = maxdegree(tree)

    for i in 1:20
        degseq = RandomTree.degrees(tree)
        walk_seq = cumsum(broadcast(-, degseq, 1))
        @test all(x-> x>=0 && x <= mdeg, degseq)
        @test all(x-> x>=0, walk_seq[1:end-1])
        @test walk_seq[end] == -1
    end
end

@testset "RandomTree.jl" begin
@testset "tree.jl" begin
    @testset "CayleyTree" begin
        tree = CayleyTree(50)
        treetest(tree)
    end

    @testset "BinaryTree" begin
        tree = BinaryTree(50)
        treetest(tree)
    end

    @testset "CatalanTree" begin
        tree = CatalanTree(50)
        treetest(tree)
    end

    @testset "MotzkinTree" begin
        tree = MotzkinTree(50)
        treetest(tree)
    end

    @testset "DAryTree" begin
        for d in 2:5
            tree = DAryTree(50, d)
            treetest(tree)
        end
    end

    @testset "RandomRecursiveTree" begin
        tree = RandomRecursiveTree(50)
        treetest(tree)
    end

    @testset "FullDAryTree" begin
        tree = FullDAryTree(5, 1)
        @test length(tree) == 6
        @test collect(degrees(tree)) == [1, 1, 1, 1, 1, 0]
        tree = FullDAryTree(3, 2)
        @test length(tree) == 2^4-1
        @test collect(degrees(tree)) == [2, 2, 2, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, 0]
        treetest(tree)
    end

    @testset "Degree sequence" begin
        testseq = [
                    [[0, 0, 1, 2, 1], 3, [1, 2, 1, 0, 0]],
                    [[0, 0, 0, 2, 2], 4, [2, 2, 0, 0, 0]],
                    [[1, 1, 1, 1, 0], 1, [1, 1, 1, 1, 0]],
                    ]

        for case in testseq
            @test RandomTree.find_startindex(case[1]) == case[2]
        end

        for case in testseq
            @test RandomTree.fix_degseq(case[1]) == case[3]
        end

        dist = Geometric(1/2)
        n = 100
        for i in 1:200
            degcount = RandomTree.degrees_count(n, dist, n)
            @test all(x-> x>= 0, degcount)
            @test sum((i-1)*degcount[i] for i in 1:length(degcount)) == n-1
        end

        for d in 2:5
            dist = Binomial(d, 1/d)
            n = 100
            for i in 1:200
                degcount = RandomTree.degrees_count(n, dist, d)
                @test length(degcount) <= d+1
                @test all(x-> x>= 0, degcount)
                @test sum((i-1)*degcount[i] for i in 1:length(degcount)) == n-1
            end
        end

        @test RandomTree.degrees(10, [4,4,1,1]) == vec([0 0 0 0 1 1 1 1 2 3])
        @test RandomTree.degrees(10, [4,3,3]) == vec([0 0 0 0 1 1 1 2 2 2])
    end
end

@testset "walker.jl" begin
    tree = CayleyTree(5)

    @testset "DepthWalker" begin
        walker = DepthWalker(tree)
        @test walk(tree, walker) == [2, 1, 1, 0, 0]
        @test result(walker) ==  [0, 1, 2, 3, 1]

        testseq= [
                    [[1, 2, 1, 0, 0], [0, 1, 2, 3, 2]],
                    [[2, 1, 0, 1, 0], [0, 1, 2, 1, 2]],
                    [[2, 1, 1, 0, 0], [0, 1, 2, 3, 1]],
                    ]

        for (degseq, depth_seq) in testseq
            walker = DepthWalker(tree)
            @test degseq == walk(degseq, walker)
            @test result(walker) == depth_seq
        end
    end

    @testset "SubtreeSizeWalker" begin
        walker = SubtreeSizeWalker(tree)
        @test walk(tree, walker) == [2, 1, 1, 0, 0]
        @test result(walker) ==  [5, 3, 2, 1, 1]

        testseq= [
                    [[1, 2, 1, 0, 0], [5, 4, 2, 1, 1]],
                    [[2, 1, 0, 1, 0], [5, 2, 1, 2, 1]],
                    [[2, 1, 1, 0, 0], [5, 3, 2, 1, 1]],
                    ]

        for (degseq, subtree_size_seq) in testseq
            walker = SubtreeSizeWalker(tree)
            @test degseq == walk(degseq, walker)
            @test result(walker) == subtree_size_seq
        end
    end

    @testset "GraphWalker" begin
        testseq= [
                    [[1, 1, 2, 0, 0], [(1, 2), (2, 3), (3, 4), (3, 5)]],
                    [[2, 1, 1, 0, 0], [(1, 2), (2, 3), (3, 4), (1, 5)]],
                    ]
        for (degseq, ret_edges) in testseq
            walker = GraphWalker(tree)
            @test size(walker.tree_digraph) == size(tree)
            @test degseq == walk(degseq, walker)
            @test RandomTree.edges(result(walker)) == ret_edges
        end
    end
end

@testset "simulator" begin
    SIZE = 101

    trees = [CayleyTree(SIZE), BinaryTree(SIZE)]

    @testset "KcutSimulator" begin
        for tree = trees
            k = 2
            sim = KcutSimulator(tree, k)

            @test repr(sim) != ""

            @test sim.tree === tree
            @test sim.k == k

            r1, r2 = simulation(sim)
            @test  SIZE >= r1 >= r2 >0

            for ret = simulation(sim, 20)
                r1, r2 = ret
                @test  SIZE >= r1 >= r2 >0
            end
        end
    end

    @testset "LogProductSimulator" begin
        # For binary tree of size 5, this is always log(15)
        sim = LogProductSimulator(BinaryTree(5))
        @test simulation(sim) ≈ log(5)+log(3)

        for tree = trees
            for pow = 1:3
                sim = LogProductSimulator(tree, pow)

                @test repr(sim) != ""

                @test sim.tree === tree

                ret = simulation(sim)
                @test ret isa Real
                @test ret > log(SIZE)^pow
            end
        end
    end

    @testset "SizePowerSimulator" begin
        sim = SizePowerSimulator(BinaryTree(5))
        @test simulation(sim) ≈ sqrt(5)+sqrt(3)+3

        for tree = trees
            for pow = 1:3
                sim = SizePowerSimulator(tree, pow)

                @test repr(sim) != ""

                @test sim.tree === tree

                ret = simulation(sim)
                @test ret isa Real
                @test ret > SIZE^pow
            end
        end
    end

    @testset "HeightSimulator" begin
        for tree = trees
            sim = HeightSimulator(tree)

            @test repr(sim) != ""

            ret = simulation(sim)
            @test ret >= 1 && ret <= size(tree)

            for ret = simulation(sim, 20)
                ret = simulation(sim)
                @test ret >= 1 && ret <= size(tree)
            end
        end
    end

    @testset "TotalPathSimulator" begin
        for tree = trees
            sim = TotalPathSimulator(tree)

            @test repr(sim) != ""

            ret = simulation(sim)
            @test ret >= size(tree)-1 && ret <= size(tree)^2

            for ret = simulation(sim, 20)
                ret = simulation(sim)
                @test ret >= 1 && ret <= size(tree)^2
            end
        end
    end
end

@testset "FixedGraph" begin
    tgraph = RandomTree.FixedTreeGraph(6, [(1, 2); (2, 3); (3, 4); (1, 5); (4, 6)])

    sortedge!(tgraph)
    @test edges(tgraph) == [(1, 5), (1, 2), (2, 3), (3, 4), (4, 6)]

    position_array = RandomTree.positions(tgraph)

    @test position_array == [1 2; 3 1; 4 1; 5 1; 0 0; 0 0]
    @test collect(degrees(tgraph)) == [2, 1, 1, 1, 0, 0]
end
end


text:
#
# Provides default frames, and the means for users to create their own.
#

"""
$(TYPEDEF)

A `supertype` for all coordinate frames used in space.
"""
abstract type OrbitalFrame end

"""
$(TYPEDEF)

A `supertype` for all inertial coordinate frames in space.
"""
abstract type Inertial <: OrbitalFrame end

"""
$(TYPEDEF)

A `supertype` for all bodycentric-inertial coordinate frames in space.
"""
abstract type BodycentricInertial <: Inertial end

"""
$(TYPEDEF)

A `supertype` for all barycentric-ineretial coordinate frames in space.
"""
abstract type BarycentricInertial <: Inertial end

"""
$(TYPEDEF)

A `supertype` for all rotating coordinate frames in space.
"""
abstract type Rotating <: OrbitalFrame end

"""
$(TYPEDEF)

A `supertype` for all bodycentric-rotating coordinate frames in space.
"""
abstract type BodycentricRotating <: Rotating end

"""
$(TYPEDEF)

A `supertype` for all barycentric-rotating coordinate frames in space.
"""
abstract type BarycentricRotating <: Rotating end

"""
$(SIGNATURES)

Returns `true` or `false` to indicate whether the provided frame is `Inertial`.
"""
isinertial(::Type{T}) where T <: OrbitalFrame = T <: Inertial

"""
$(SIGNATURES)

Returns `true` or `false` to indicate whether the provided frame is `Rotating`.
"""
isrotating(::Type{T}) where T <: OrbitalFrame = T <: Rotating

"""
$(SIGNATURES)

Returns `true` or `false` to indicate whether the provided frame is bodycentric.
"""
isbodycentric(::Type{T}) where T <: OrbitalFrame = T <: Union{<:BodycentricInertial, BodycentricRotating}

"""
$(SIGNATURES)

Returns `true` or `false` to indicate whether the provided frame is barycentric.
"""
isbarycentric(::Type{T}) where T <: OrbitalFrame = T <: Union{<:BarycentricInertial, BarycentricRotating}

"""
$(SIGNATURES)

A simple convenience macro for defining new coordinate frames in space.
Creates a new `abstract` sub-type of the final argument, `property`.

# Extended Help

**Usage**

`@frame <name>`
`@frame <name> is <property>`

Here, `<property>` is previously defined `OrbitalFrame`.
which must all be subtypes of `OrbitalFrames.OrbitalFrame`.
By using the shortened syntax, `@frame <name>`, your new `abstract`
type `<name>` will only have `OrbitalFrames.OrbitalFrame` as
a `supertype`.

**Examples**
```julia
@frame AbstractFrame # second argument defaults to the top-level `OrbitalFrame`
@frame EarthCenteredInertial is BodycentricInertial
@frame SunEarthSynodic is BarycentricInertial
```
"""
macro frame(name)
    quote
        abstract type $name <: $OrbitalFrame end
    end
end
macro frame(name, _is, super)
    @assert :($_is) == :is "Incorrect Usage! Use `is` as the second argument. For example: `@frame ECI is BodycentricInertialFrame`."
    quote
        abstract type $name <: $super end
    end
end


text:
struct PointForce{T}
  x :: Vector{Float64}
  f0 :: Union{Float64,Vector{Float64}}
  t0 :: Float64
  σt :: Float64
  pulse :: ModulatedField
end

function PointForce(u::GridData,x0::Vector{Float64},f0,t0,σt,sys::NavierStokes)
  forcefield = GeneratedField(deepcopy(u),_createfields(f0,x0,sys),sys.grid)
  pulse = PulseField(forcefield,t0,σt)
  return PointForce{typeof(u)}(x0,f0,t0,σt,pulse)
end

PointForce(u::GridData,x0::Tuple{Float64,Float64},f0,t0,σt,sys::NavierStokes) =
        PointForce(u,[x0...],f0,t0,σt,sys)

function _createfields(f0::Vector{T},x0,sys) where {T<:Real}
  fieldvec = AbstractSpatialField[]
  for f in f0
    push!(fieldvec,SpatialGaussian(cellsize(sys),x0...,f*cellsize(sys)^2))
  end
  return fieldvec
end

function _createfields(f0::T,x0,sys) where {T<:Real}
  return SpatialGaussian(cellsize(sys),x0...,f0*cellsize(sys)^2)
end

(p::PointForce)(t) = p.pulse(t)

function Base.show(io::IO,f::PointForce{T}) where {T}
    println(io,"Transient point force applied on the $(T) field.")
    println(io,"   strength = $(f.f0)")
    println(io,"   location = $(f.x)")
    println(io,"   central time = $(f.t0)")
    println(io,"   half-interval = $(f.σt)")
end


text:
module MyTest2

function sayhello()
    println("Hello world")
end

function Add_two_numbers(n1::Int,n2::Int)
    return n1+n2
end

export Add_two_numbers
export sayhello
end


text:
include("./artificial_data.jl")

using CellwiseEditingDifferentiation
using Color
using Gadfly

for i in 1:nsites
  colors = distinguishable_colors(ncells[i])

  layers = Layer[]

  pchain = readdlm(joinpath(OUTDIR, @sprintf("pchain_%s_site%02d.txt", string(simulationid), i)))

  for j in 1:ncells[i]
    push!(layers, layer(
      x=pchain[:, j],
      Stat.histogram(bincount=50, density=true),
      Geom.line,
      Theme(default_color=colors[j])
    )[1])
  end

  vplot = plot(
    layers,
    Guide.xlabel("p$i"),
    Guide.title("Histograms of conditional p$i")
    # Guide.manual_color_key("Distribution", ["cell "*string(k) for k in cells[i]], [c for c in colors])
  )

  draw(PDF(joinpath(OUTDIR, @sprintf("conditional_phist_%s_site%02d.pdf", string(simulationid), i)), 4inch, 3inch), vplot)
end


text:
# Weeks–Chandler–Anderson potential function, a repulsive potential.

@inline lj_potential_wca_uij(
    r; ϵ = 1.0, σ = 1.0
) = norm_sqr(r) ≥ (R_MIN * σ) ^ 2 ? 0.0 :
    (lj_potential_uij(r, ϵ = ϵ, σ = σ) + ϵ)

@inline lj_potential_wca_uij(
    r₁, r₂;
    ϵ = 1.0, σ = 1.0,
    dist = default_distance
) = lj_potential_wca_uij(dist(r₁, r₂), ϵ = ϵ, σ = σ)

@inline lj_potential_wca_fij(
    r; ϵ = 1.0, σ = 1.0
) = norm_sqr(r) ≥ (R_MIN * σ) ^ 2 ? zero(Vector3) :
    lj_potential_fij(r, ϵ = ϵ, σ = σ)

@inline lj_potential_wca_fij(
    r₁, r₂;
    ϵ = 1.0, σ = 1.0,
    dist = default_distance
) = lj_potential_wca_fij(dist(r₁, r₂), ϵ = ϵ, σ = σ)


text:
using Statistics

###Input Layer

@doc raw"""
    layerForProp(
        cLayer::Input,
        X::AbstractArray = Array{Any,1}(undef,0);
        FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
        kwargs...,
    )

Perform forward propagation for `Input` `Layer`

# Arguments

- `cLayer` := the layer to perform for prop on

- `X` := is the input data of the `Input` `Layer`

- `FCache` := a cache holder of the for prop

# Return

- A `Dict{Symbol, AbstractArray}(:A => Ao)`

"""
function layerForProp(
    cLayer::Input,
    X::AbstractArray = Array{Any,1}(undef,0);
    FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
    kwargs...,
)
    if length(X) != 0
        # cLayer.A = X
        cLayer.inputS = cLayer.outputS = size(X)
    end
    cLayer.forwCount += 1
    # Base.GC.gc()
    return Dict(:A=>X)
end

###FCLayer forprop

@doc raw"""
    layerForProp(
        cLayer::FCLayer,
        Ai::AbstractArray = Array{Any,1}(undef,0);
        FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
        kwargs...,
    )

Perform forward propagation for `FCLayer` `Layer`

# Arguments

- `cLayer` := the layer to perform for prop on

- `Ai` := is the input activation of the `FCLayer` `Layer`

- `FCache` := a cache holder of the for prop

# Return

- A `Dict{Symbol, AbstractArray}(:A => Ao, :Z => Z)`
"""
function layerForProp(
    cLayer::FCLayer,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
    kwargs...,
)
    prevLayer = cLayer.prevLayer
    if length(Ai) == 0
        Ai = FCache[prevLayer][:A]
    end
    cLayer.inputS = cLayer.prevLayer.outputS
    Z = cLayer.W * Ai .+ cLayer.B
    actFun = cLayer.actFun
    # Z = cLayer.Z
    cLayer.outputS = size(Z)
    A = eval(:($actFun($Z)))
    cLayer.forwCount += 1
    # Base.GC.gc()
    return Dict(:Z=>Z, :A=>A)
end #function layerForProp(cLayer::FCLayer)

###AddLayer forprop

@doc raw"""
    layerForProp(
        cLayer::AddLayer;
        FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
        kwargs...,
    )

Perform forward propagation for `AddLayer` `Layer`

# Arguments

- `cLayer` := the layer to perform for prop on

- `FCache` := a cache holder of the for prop

# Return

- A `Dict{Symbol, AbstractArray}(:A => Ao)`
"""
function layerForProp(
    cLayer::AddLayer;
    FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
    kwargs...,
)

    A = similar(FCache[cLayer.prevLayer[1]][:A]) .= 0
    # if all(
    #     i -> (i.forwCount == cLayer.prevLayer[1].forwCount),
    #     cLayer.prevLayer,
    # )
        for prevLayer in cLayer.prevLayer
            A .+= FCache[prevLayer][:A]
        end
    # end #if all
    cLayer.forwCount += 1
    # Base.GC.gc()
    return Dict(:A=>A)
end #function layerForProp(cLayer::AddLayer)

### ConcatLayer forprop

function layerForProp(
    cLayer::ConcatLayer;
    FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
    kwargs...,
)

    N = ndims(FCache[cLayer.prevLayer[1]][:A])
    A = cat([FCache[prevLayer][:A] for prevLayer in cLayer.prevLayer]...; dims=N-1)

    cLayer.forwCount += 1
    # Base.GC.gc()
    return Dict(:A=>A)
end #function layerForProp(cLayer::AddLayer)

###Activation forprop

@doc raw"""
    layerForProp(
        cLayer::Activation,
        Ai::AbstractArray = Array{Any,1}(undef,0);
        FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
        kwargs...,
    )

Perform forward propagation for `Activation` `Layer`

# Arguments

- `cLayer` := the layer to perform for prop on

- `Ai` := is the input activation of the `Activation` `Layer`

- `FCache` := a cache holder of the for prop

# Return

- A `Dict{Symbol, AbstractArray}(:A => Ao)`
"""
function layerForProp(
    cLayer::Activation,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
    kwargs...,
)
    prevLayer = cLayer.prevLayer
    if length(Ai) == 0
        Ai = FCache[prevLayer][:A]
    end
    actFun = cLayer.actFun
    A = eval(:($actFun($Ai)))
    cLayer.inputS = cLayer.outputS = size(Ai)
    cLayer.forwCount += 1
    # Ai = nothing
    # Base.GC.gc()
    return Dict(:A=>A)
end #function layerForProp(cLayer::Activation)

### Flatten

@doc raw"""
    layerForProp(
        cLayer::Flatten,
        Ai::AbstractArray = Array{Any,1}(undef,0);
        FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
        kwargs...,
    )

Perform forward propagation for `Flatten` `Layer`

# Arguments

- `cLayer` := the layer to perform for prop on

- `Ai` := is the input activation of the `Flatten` `Layer`

- `FCache` := a cache holder of the for prop

# Return

- A `Dict{Symbol, AbstractArray}(:A => Ao)`
"""
function layerForProp(
    cLayer::Flatten,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
    kwargs...,
)
    prevLayer = cLayer.prevLayer
    if length(Ai) == 0
        Ai = FCache[prevLayer][:A]
    end

    cLayer.inputS = size(Ai)

    cLayer.outputS = (prod(cLayer.inputS[1:end-1]), cLayer.inputS[end])

    A = reshape(Ai,cLayer.outputS)

    cLayer.forwCount += 1
    # Ai = nothing
    # Base.GC.gc()
    return Dict(:A=>A)
end #function layerForProp(cLayer::Activation)

###BatchNorm

@doc raw"""
    layerForProp(
        cLayer::BatchNorm,
        Ai::AbstractArray = Array{Any,1}(undef,0);
        FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
        kwargs...,
    )

Perform forward propagation for `BatchNorm` `Layer` and trainable parameters W and B

# Arguments

- `cLayer` := the layer to perform for prop on

- `Ai` := is the input activation of the `BatchNorm` `Layer`

- `FCache` := a cache holder of the for prop

# Return

- `Dict(
        :μ => μ,
        :Ai_μ => Ai_μ,
        :Ai_μ_s => Ai_μ_s,
        :var => var,
        :Z => Z,
        :A => Ao,
        :Ap => Ap,
        )`

"""
function layerForProp(
    cLayer::BatchNorm,
    Ai::AbstractArray = Array{Any,1}(undef,0);
    FCache::Dict{Layer,Dict{Symbol, AbstractArray}},
    kwargs...,
)

    prediction = getindex(kwargs, :prediction; default=false)
    prevLayer = cLayer.prevLayer
    if length(Ai) == 0
        Ai = FCache[prevLayer][:A]
    end

    cLayer.ϵ = eltype(cLayer.W)(cLayer.ϵ)

    if prediction
        cLayer.forwCount += 1
        NDim = cLayer.dim
        μ = mean(Ai, dims = 1:NDim)
        Ai_μ = Ai .- μ
        Num = prod(size(Ai)[1:NDim])
        Ai_μ_s = Ai_μ .^ 2
        var = sum(Ai_μ_s, dims = 1:NDim) ./ Num
        Z = Ai_μ ./ sqrt.(var .+ cLayer.ϵ)
        Ao = cLayer.W .* Z .+ cLayer.B
        return Dict(:A => Ao)
    end #prediction

    # initWB!(cLayer)
    # initVS!(cLayer, model.optimizer)

    tA = eltype(Ai)
    N = ndims(Ai)
    Ai = permutedims(Ai, [N,(1:N-1)...])
    if tA==Float16
        Ai = Float64.(Ai)
    end
    NDim = cLayer.dim + 1

    μ = (mean(Ai, dims = 1:NDim))
    Ai_μ = (Ai .- μ)
    Num = prod(size(Ai)[1:NDim])
    Ai_μ_s = (Ai_μ .^ 2)
    var = tA.(sum(Ai_μ_s, dims = 1:NDim) ./ Num)
    Z = tA.(Ai_μ ./ sqrt.(var .+ cLayer.ϵ))
    Ap = cLayer.W .* Z .+ cLayer.B
    Ao = permutedims(Ap, [(2:N)..., 1])
    cLayer.forwCount += 1
    # Ai_μ = nothing
    cLayer.inputS = cLayer.outputS = size(Ao)
    # Base.GC.gc()
    return Dict(
        :μ => tA.(μ),
        :Ai_μ => tA.(Ai_μ),
        :Ai_μ_s => tA.(Ai_μ_s),
        :var => tA.(var),
        :Z => tA.(Z),
        :A => tA.(Ao),
        :Ap => tA.(Ap),
    )
end #function layerForProp(cLayer::BatchNorm)

export layerForProp


text:
using Documenter
using Wordlegames

DocMeta.setdocmeta!(Wordlegames, :DocTestSetup, :(using Wordlegames); recursive=true)

makedocs(
    sitename = "Wordlegames.jl",
    authors="Douglas Bates <dmbates@gmail.com> and contributors",
    repo="https://github.com/dmbates/Wordlegames.jl/blob/{commit}{path}#{line}",
    doctest = true,
    format = Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://dmbates.github.io/Wordlegames.jl",
        assets=String[],
    ),
    modules = [Wordlegames],
    pages=[
        "Home" => "index.md",
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(;
    repo = "github.com/dmbates/Wordlegames.jl",
    devbranch="main",
    push_preview=true,
)


text:
module StackOverflow

using HTTP
using CodecZlib
using JSON

include("questions.jl")
export
    getrecentquestionsfortag

include("privileges.jl")
export
    checkmyprivileges,
    getprivileges

include("utils.jl")
export
    convert_HTTP_Response_To_JSON

include("answers.jl")
export
    getquestions

include("badges.jl")
export
    getrecentbadgerecipients

end # module


text:
#
#  BoundarForces.jl
#  CollagenPolymer
#
#  Created by Christopher Revell on 04/08/2020.
#
#

module BoundaryForces

# Import Julia packages
using LinearAlgebra
using StaticArrays

# Import local program modules

@views function boundaryForces!(boundaryList,pos,F,ϵ,nGrid,boxSize,dxMatrix,rₘ)

	# Boundary forces
	for (particle,dimension,edge) in boundaryList
		if edge == 1
			dxMag = boxSize/2.0 + pos[particle][dimension]
			if dxMag < rₘ
				Fmag = (12.0*ϵ/rₘ)*(exp(6.0*(1.0-dxMag/rₘ))-exp(12.0*(1.0-dxMag/rₘ)))
				F[particle] -= (Fmag/dxMag)*dxMatrix[dimension,:]
			end
		end
		if edge == nGrid
			dxMag = boxSize/2.0 - pos[particle][dimension]
			if dxMag < rₘ
				Fmag = (12.0*ϵ/rₘ)*(exp(6.0*(1.0-dxMag/rₘ))-exp(12.0*(1.0-dxMag/rₘ)))
				F[particle] += (Fmag/dxMag)*dxMatrix[dimension,:]
			end
		end
	end

	return nothing

end

export boundaryForces!

end

max_stars_repo_name: chris-revell/collagen-model-jl

max_issues_repo_name: chris-revell/collagen-model-jl

max_forks_repo_name: chris-revell/collagen-model-jl


text:
module Aeconomica

import HTTP, JSON, DataFrames, Dates

export fetch_series, set_apikey, fetch_series_as_at, fetch_dataset

mutable struct APIKey
    key::String
end

const stored_apikey = APIKey("")

function apikey()
    if stored_apikey.key == ""
        throw(ErrorException("You have not set an `apikey`. Please run `set_apikey(\"YOURKEY\")` first."))
    else
        stored_apikey.key
    end
end

function set_apikey(key)
    stored_apikey.key = key
end

include("utils.jl")
include("series.jl")
include("series_as_at.jl")
include("dataset.jl")

end # module


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(kalman.deltaT, (ASCIIString,))
    precompile(kalman.recordSampleTime, (ASCIIString, Int32,))
    precompile(kalman.fThunk, (ASCIIString,))
    precompile(kalman.initialGuess, (Array{Float64, 1}, Array{Float64, 1},))
    precompile(kalman.qThunk, (ASCIIString, Float64,))
    precompile(kalman.newModel, (ASCIIString, Float64, Float64, Int32,))
end

max_stars_repo_name: waTeim/gpsFlow

max_issues_repo_name: waTeim/flying-monkey

max_forks_repo_name: waTeim/flying-monkey


text:
"""
    train!(model::AbstractGP; iterations::Integer=100, callback, convergence)

Function to train the given GP `model`.
## Arguments
- `model` : AbstractGP model with either an `Analytic`, `AnalyticVI` or `NumericalVI` type of inference

## Keyword Arguments
- `iterations::Int` : Number of iterations (not necessarily epochs!)for training
- `callback::Function=nothing` : Callback function called at every iteration. Should be of type `function(model,iter) ...  end`
- `convergence::Function=nothing` : Convergence function to be called every iteration, should return a scalar and take the same arguments as `callback`
"""
function train!(
    model::AbstractGP{T},
    iterations::Int=100;
    callback=nothing,
    convergence=nothing,
) where {T}
    if model.verbose > 0
        @info "Starting training $model with $(nSamples(model)) samples, $(nFeatures(model)) features and $(nLatent(model)) latent GP" *
              (nLatent(model) > 1 ? "s" : "")
    end

    iterations > 0 || error("Number of iterations should be positive")
    # model.evol_conv = [] #Array to check on the evolution of convergence
    local_iter = 1
    conv = Inf
    p = Progress(iterations; dt=0.2, desc="Training Progress: ")
    prev_elbo = -Inf
    while true #loop until one condition is matched
        try #Allow for keyboard interruption without losing the model
            update_parameters!(model) #Update all the variational parameters
            set_trained!(model, true)
            if !isnothing(callback)
                callback(model, nIter(model)) #Use a callback method if set by user
            end
            if (nIter(model) % model.atfrequency == 0) &&
               (nIter(model) >= 3) &&
               (local_iter != iterations)
                update_hyperparameters!(model) #Update the hyperparameters
            end
            # Print out informations about the convergence
            if verbose(model) > 2 || (verbose(model) > 1 && local_iter % 10 == 0)
                if inference(model) isa GibbsSampling
                    next!(p; showvalues=[(:samples, local_iter)])
                else
                    if (verbose(model) == 2 && local_iter % 10 == 0)
                        elbo = objective(model)
                        prev_elbo = elbo
                        ProgressMeter.update!(p, local_iter - 1)
                        ProgressMeter.next!(
                            p; showvalues=[(:iter, local_iter), (:ELBO, elbo)]
                        )
                    elseif verbose(model) > 2
                        elbo = objective(model)
                        prev_elbo = elbo
                        ProgressMeter.next!(
                            p; showvalues=[(:iter, local_iter), (:ELBO, prev_elbo)]
                        )
                    end
                end
            end
            local_iter += 1
            model.inference.nIter += 1
            (local_iter <= iterations) || break #Verify if the number of maximum iterations has been reached
        # (iter < model.nEpochs && conv > model.ϵ) || break; #Verify if any condition has been broken
        catch e
            if isa(e, InterruptException)
                @warn "Training interrupted by user at iteration $local_iter"
                break
            else
                rethrow(e)
            end
        end
    end
    if verbose(model) > 0
        @info "Training ended after $(local_iter - 1) iterations. Total number of iterations $(nIter(model))"
    end
    compute_kernel_matrices!(model, true) # Compute final version of the matrices for predictions
    post_step!(model)
    set_trained!(model, true)
    return nothing
end

function update_parameters!(model::GP)
    compute_kernel_matrices!(model) #Recompute the matrices if necessary (when hyperparameters have been updated)
    analytic_updates!(model)
    return nothing
end

function update_parameters!(model::VGP)
    compute_kernel_matrices!(model) #Recompute the matrices if necessary (always for the stochastic case, or when hyperparameters have been updated)
    variational_updates!(model)
    return nothing
end

## Update all variational parameters of the sparse variational GP Model ##
function update_parameters!(m::SVGP)
    if isStochastic(inference(m))
        setMBIndices!(
            inference(m),
            StatsBase.sample(1:nSamples(m), nMinibatch(inference(m)); replace=false),
        )
        setxview!(inference(m), view_x(data(m), MBIndices(inference(m))))
        setyview!(inference(m), view_y(likelihood(m), data(m), MBIndices(inference(m))))
    end
    compute_kernel_matrices!(m) #Recompute the matrices if necessary (always for the stochastic case, or when hyperparameters have been updated)
    variational_updates!(m)
    return nothing
end

function update_parameters!(m::MOVGP)
    compute_kernel_matrices!(m) #Recompute the matrices if necessary (always for the stochastic case, or when hyperparameters have been updated)
    update_A!(m)
    variational_updates!(m)
    return nothing
end

function update_parameters!(m::MOSVGP)
    if isStochastic(m.inference)
        setMBIndices!(
            inference(m),
            StatsBase.sample(1:nSamples(m), nMinibatch(inference(m)); replace=false),
        )
        setxview!(inference(m), view_x(data(m), MBIndices(inference(m))))
        setyview!(m.inference, view_y(likelihood(m), data(m), MBIndices(m.inference)))
    end
    compute_kernel_matrices!(m) #Recompute the matrices if necessary (always for the stochastic case, or when hyperparameters have been updated)
    update_A!(m)
    variational_updates!(m)
    return nothing
end

function update_parameters!(m::VStP)
    compute_kernel_matrices!(m) #Recompute the matrices if necessary (always for the stochastic case, or when hyperparameters have been updated)
    local_prior_updates!(m, input(m))
    variational_updates!(m)
    return nothing
end

function compute_kernel_matrices!(m::GP{T}, ::Bool) where {T}
    compute_K!(getf(m), input(m), T(jitt))
    setHPupdated!(inference(m), false)
    return nothing
end

@traitfn function compute_kernel_matrices!(
    m::TGP, update::Bool=false
) where {T,TGP<:AbstractGP{T};IsFull{TGP}}
    if isHPupdated(inference(m)) || update
        compute_K!.(getf(m), Ref(input(m)), T(jitt))
    end
    setHPupdated!(inference(m), false)
    return nothing
end

@traitfn function compute_kernel_matrices!(
    m::TGP, update::Bool=false
) where {T,TGP<:AbstractGP{T};!IsFull{TGP}}
    if isHPupdated(inference(m)) || update
        compute_K!.(getf(m), T(jitt))
    end
    #If change of hyperparameters or if stochatic
    if isHPupdated(inference(m)) || isStochastic(inference(m)) || update
        compute_κ!.(getf(m), Ref(xview(m)), T(jitt))
    end
    setHPupdated!(inference(m), false)
    return nothing
end


text:
using JuMP
using SCIP
using GLPK

solver = GndStr.IterGBD(
    JuMP.optimizer_with_attributes(
        SCIP.Optimizer,
        "display/width" => 139,
        "limits/memory" => 5000.0
    ),
    GLPK.Optimizer,
    addnogoods=true,
    addcritpath=false,
    maxiter=Int(1e9),
    timelimit=3600.0)