text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module ModiaMath.Result

Organize and plot simulation result data (time series).
The result data of [`ModiaMath.simulate!`](@ref) is returned in one of the
formats supported by this module. The [`ModiaMath.plot`](@ref) function of this
module allows to plot the result data by giving the signal names.
The legends/labels of the plots are automatically constructed by the
signal names and their unit. Example:

```julia
ModiaMath.plot(result, [ (:phi,:r)      (:phi,:phi2,:w);
                         (:w,:w2,:phi2) (:phi,:w)      ],
               heading="Matrix of plots")
```

generates the following plot:

![Matrix-of-Plots](../../resources/images/matrix-of-plots.svg)

# Main developer

[Martin Otter](https://rmc.dlr.de/sr/de/staff/martin.otter/),
[DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module Result

export getStringDictResult, SymbolDictResult, StringDictResult, ResultWithVariables
export resultHeading, resultTimeSeries, resultTable, plot
export RawResult, nResults, storeRawResult!
export closeFigure, closeAllFigures
export variablesDependingOnStruct

# using/imports
import ..ModiaMath
import DataFrames
using  Unitful
using  Requires

# Constants
const headingSize = 10

# include code
include("rawResult.jl")
include("result.jl")
include("plotResult.jl")

# Dummy plot functions if PyPlot is not in the current environment
include("plotResult_with_Nothing.jl")

# If PyPlot is in the current environment, import it in the REPL and install plot functions based on PyPlot
function __init__()
    if !Requires.isprecompiling()
        @eval Main begin
            pyplotAvailable = false
            try
                import PyPlot
                global pyplotAvailable = true
            catch
                println("    PyPlot not available (plot commands will be ignored).\n",
                        "    Try to install PyPlot. See hints here:\n",
                        "    https://github.com/ModiaSim/ModiaMath.jl/wiki/Installing-PyPlot-in-a-robust-way.")
            end

            if pyplotAvailable
                try
                    import PyCall
                catch
                    println("... From ModiaMath: PyCall is added, since needed to set default options for PyPlot")
                    import Pkg
                    Pkg.add("PyCall")
                    import PyCall
                end

                pyplot_rc = PyCall.PyDict(PyPlot.matplotlib["rcParams"])
                pyplot_rc["axes.formatter.limits"] = [-3,4]
                pyplot_rc["font.size"]        = 8.0
                pyplot_rc["lines.linewidth"]  = 1.0
                pyplot_rc["grid.linewidth"]   = 0.5
                pyplot_rc["axes.grid"]        = true
                pyplot_rc["axes.titlesize"]   = "medium"
                pyplot_rc["figure.titlesize"] = "medium"

                # Use separate plot windows (no inline plots)
                PyPlot.pygui(true)
            end
        end
    end

    @require PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" include("plotResult_with_PyPlot.jl")
end

end


text:
using AstroTime
using Dates
using Test

AstroTime.load_test_eop()

const speed_of_light = 299792458.0 # m/s
const astronomical_unit = 149597870700.0 # m

# TODO: Change me
@timescale GMT TAI
AstroTime.Epochs.getoffset(::GMTScale, ::InternationalAtomicTime, _, _) = 0.0
AstroTime.Epochs.getoffset(::InternationalAtomicTime, ::GMTScale, _, _) = 0.0

@timescale SCET TAI
function AstroTime.Epochs.getoffset(::SCETScale, ::InternationalAtomicTime,
                                    _, _, distance)
    return distance / speed_of_light
end
function AstroTime.Epochs.getoffset(::InternationalAtomicTime, ::SCETScale,
                                    _, _, distance)
    return -distance / speed_of_light
end

@timescale Dummy TDB

@timescale Lonely
@timescale Together Lonely
function AstroTime.Epochs.getoffset(::LonelyScale, ::TogetherScale, _, _)
    return 5.0
end
function AstroTime.Epochs.getoffset(::TogetherScale, ::LonelyScale, _, _)
    return -5.0
end

@timescale OneWay TDB true

@testset "AstroTime" begin
    include("periods.jl")
    include("astrodates.jl")
    include("epochs.jl")
    @testset "Custom Time Scales" begin
        # TODO: Change me
        tai = TAIEpoch(2000, 1, 1, 0, 0, 0.1)
        gmt = GMTEpoch(2000, 1, 1, 0, 0, 0.1)
        @test tai == TAIEpoch(gmt)
        @test gmt == GMTEpoch(tai)

        @test string(GMT) == "GMT"
        @test typeof(GMT) == GMTScale
        @test string(gmt) == "2000-01-01T00:00:00.100 GMT"
        @test find_path(GMT, TAI) == [GMT, TAI]

        @test string(SCET) == "SCET"
        @test typeof(SCET) == SCETScale
        @test find_path(SCET, TAI) == [SCET, TAI]
        scet = SCETEpoch(2000, 1, 1, 0, 0, 0.1)
        tai_exp = TAIEpoch(2000, 1, 1, 0, 8, 19.10478383615643)
        @test TAIEpoch(scet, astronomical_unit) ≈ tai_exp
        @test SCETEpoch(tai_exp, astronomical_unit) ≈ scet

        @test Dates.default_format(SCETEpoch) == AstroTime.EPOCH_ISO_FORMAT[]

        @test string(Dummy) == "Dummy"
        @test typeof(Dummy) == DummyScale
        @test find_path(Dummy, TAI) == [Dummy, TDB, TT, TAI]
        dummy = DummyEpoch(2000, 1, 1)
        @test_throws NoOffsetError TAIEpoch(dummy)
        @test_throws NoOffsetError TDBEpoch(dummy)

        lonely = LonelyEpoch(2000, 1, 1)
        together = TogetherEpoch(lonely)
        @test together == TogetherEpoch(2000, 1, 1, 0, 0, 5.0)
        @test lonely == LonelyEpoch(together)
        @test find_path(Lonely, TDB) == []

        @test find_path(TDB, OneWay) == [TDB, OneWay]
        @test find_path(OneWay, TDB) == []
        @test_throws NoPathError TDBEpoch(OneWayEpoch(2000, 1, 1))
    end
end

ext: jl

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
module ZOOclient

using Base.Dates.now

export Dimension, dim_print
export Objective, get_history_bestsofar, obj_clean_history
export zoo_min
export Parameter, print_population, get_positive_data, get_negative_data
export Solution, sol_print, sol_equal, sol_write, write_population

export zoolog
export rng, my_precision, set_seed, set_precision
export construct_init_sample

include("ZOOclient/utils/tool_function.jl")
include("ZOOclient/utils/zoo_global.jl")
include("ZOOclient/dimension.jl")
include("ZOOclient/objective.jl")
include("ZOOclient/parameter.jl")
include("ZOOclient/solution.jl")
include("ZOOclient/algos/racos/racos_classification.jl")
include("ZOOclient/algos/racos/racos_common.jl")
include("ZOOclient/algos/racos/racos.jl")
include("ZOOclient/algos/racos/sracos.jl")
include("ZOOclient/algos/racos/racos_optimization.jl")
include("ZOOclient/algos/asracos/asracos.jl")
include("ZOOclient/algos/asracos/asracos_opt.jl")
include("ZOOclient/algos/pracos/pracos.jl")
include("ZOOclient/algos/pracos/pracos_opt.jl")
include("ZOOclient/algos/psracos/psracos.jl")
include("ZOOclient/algos/psracos/psracos_opt.jl")
include("ZOOclient/algos/pposs/pposs_opt.jl")
include("ZOOclient/optimize.jl")

end


text:
# by Lazaro Alonso
# port from the original code used in:
# L Alonso, et. al. https://doi.org/10.1093/comnet/cnx053

# using CairoMakie # HIDE
using LinearAlgebra, Random, GLMakie
GLMakie.activate!() # hide
function RRGAdjacencyM3D(;radius = 0.17, nodes = 500, rseed = 123)
    Random.seed!(rseed)
    xy = rand(nodes, 3)
    x = xy[:,1]
    y = xy[:,2]
    z = xy[:,3]

    matrixAdjDiag = Diagonal(√2 * randn(nodes))
    matrixAdj = zeros(nodes, nodes)
    for point in 1:nodes-1
        xseps = (x[point + 1:end] .- x[point]).^2
        yseps = (y[point + 1:end] .- y[point]).^2
        zseps = (z[point + 1:end] .- z[point]).^2

        distance = sqrt.(xseps.+ yseps .+ zseps)
        dindx = findall(distance .<= radius) .+ point
        if length(dindx) > 0
            rnd = randn(length(dindx))
            matrixAdj[point, dindx] = rnd
            matrixAdj[dindx, point] = rnd
        end
    end
    return (matrixAdj .+ matrixAdjDiag, x, y, z)
end
adjacencyM3D, x, y, z = RRGAdjacencyM3D()

function getGraphEdges3D(adjMatrix3D, x, y, z)
    xyzos = []
    weights = []
    for i in 1:length(x), j in i+1:length(x)
        if adjMatrix3D[i,j] != 0.0
            push!(xyzos, [x[i], y[i], z[i]])
            push!(xyzos, [x[j], y[j], z[j]])
            push!(weights, adjMatrix3D[i,j])
            push!(weights, adjMatrix3D[i,j])
        end
    end
    return (Point3f0.(xyzos), Float32.(weights))
end

function plotGraph3D(adjacencyM3D, x, y, z)
    cmap = (:Spectral_11, 0.75)
    adjmin = minimum(adjacencyM3D)
    adjmax = maximum(adjacencyM3D)
    diagValues = diag(adjacencyM3D)
    segm, weights = getGraphEdges3D(adjacencyM3D, x, y, z)

    fig, ax, pltobj = linesegments(segm, color = weights, colormap = cmap,
        linewidth = abs.(weights), colorrange = (adjmin, adjmax),
        figure = (;resolution = (1000,800)),
        axis = (; type = Axis3, aspect =(1,1,1), perspectiveness = 0.5))
    meshscatter!(ax, x, y, z, color = diagValues, markersize = abs.(diagValues)./90,
        colorrange = (adjmin, adjmax), colormap = cmap)
    #limits!(ax, -0.02,1.02,-0.02,1.02)
    Colorbar(fig[1,2], pltobj, label = "weights", height = Relative(0.5))
    fig
end
#plotGraph(adjacencyM, x, y)
rrgraph3D = with_theme(theme_black()) do
    plotGraph3D(adjacencyM3D, x, y, z)
end
save(joinpath(@__DIR__, "output", "RRGraph3D.png"), rrgraph3D, px_per_unit = 2.0) # HIDE
using Pkg # HIDE
Pkg.status(["GLMakie", "LinearAlgebra", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
module FunctionEnsembles

export Ensemble, @ensemble, @functionfirst

using Lazy
struct Ensemble{T}
    funs::Vector{T}
end
Ensemble(e::Ensemble) = e

Base.getindex(e::Ensemble, inds...) = e.funs[inds...]
@forward Ensemble.funs Base.length, Base.size, Base.isempty, Base.getindex, Base.setindex!, Base.push!, Base.deleteat!, Base.insert!, Base.prepend!, Base.iterate

function (e::Ensemble)(args...; kwargs...)
    results = map(e.funs) do f
        f(args...; kwargs...)
    end
end

function (e::Ensemble)(op::Function, v0, args...; kwargs...)
    reduce(op, v0, f(args...; kwargs...) for f in e.funs)
end

macro ensemble(e)
    @assert e isa Expr "Call syntax: @ensemble v(args...) where v::Vector{Function}"
    if e.head == :call
        quote
            Ensemble($(esc(e.args[1])))($(esc(e.args[2:end]...)))
        end
    elseif e.head == :(=)
        quote
            $(esc(e.args[1])) = Ensemble($(esc(e.args[2].args[1])))($(esc(e.args[2].args[2:end]...)))
        end
    else
        error("Call syntax: @ensemble v(args...) where v::Vector{Function}")
    end
end

macro functionfirst(e)
    @assert e isa Expr "Call syntax: @functionfirst fun1(fun2, iterable)"
    @assert e.head == :call  "Call syntax: @functionfirst fun1(fun2, iterable)"
    quote
        outerfun = $(esc(e.args[1]))
        innerfun = $(esc(e.args[2]))
        arg = $(esc(e.args[3]))
        outerfun(innerfun.(arg))
    end
end

end # module


text:
using Dates

# Rata Die milliseconds for 1970-01-01T00:00:00
const RATAEPOCH = Dates.value(DateTime(1970))

"""
    is_odd(x::Int) -> Bool
"""
is_odd(x::Int) = x & 0x1 == 1

"""
    list_files(dir::AbstractString=".")
Returns the files in the directory dir
"""
function list_files(dir::AbstractString=".")
    rs = Vector{String}()
    for (root, dirs, files) in walkdir(dir)
        for file in files
            push!(rs, joinpath(root, file))
        end
    end
    rs
end

"""
    time_millis(dt::DateTime=now())
Returns millisecond since 1970-01-01T00:00:00
"""
time_millis(dt::DateTime=now()) = Dates.value(dt) - RATAEPOCH


text:
function temperatureasintersection(a::Array{Float64,1}, b::Array{Float64,1}, σ2::Array{Float64,1}, T::Float64, s::Float64)

    x = s*unitgalaxyvector(T)

    p = dot(a, (x-b)) / dot(a,a)

    logpdf(MvNormal(a*p + b, Diagonal(σ2)), x)

end

function temperatureasintersection(a::Array{Float64,1}, b::Array{Float64,1}, σ2::Array{Float64,1}, Trange::AbstractRange{Float64}, Srange::AbstractRange{Float64})

    f(T, s) = temperatureasintersection(a,b,σ2,T,s)

    @showprogress map(P -> f(P[1], P[2]), Iterators.product(Trange, Srange))

end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    AbstractGraph,
    UndirectedGraph, DirectedGraph, HyperGraph,
    OrderedGraph, OrderedDiGraph, OrderedHyperGraph,
    AbstractNode, AbstractEdge, UndirectedEdge, DirectedEdge,
    @cachefirst, @cache, hascache, clearcache!, clone,
    neighbors, outneighbors, inneighbors,
    findedgekey, findalledgekeys, getedge, hasedge, nodeattr, edgeattr,
    adjacencies, successors, predecessors,
    incidences, outincidences, in_incidences,
    nodeset, edgeset, edgesiter, neighborsiter, nodeattrs, edgeattrs,
    nodecount, edgecount, neighborcount, degree, indegree, outdegree,
    addnode!, addedge!, setnodeattr!, setedgeattr!,
    unlinknodes, unlinkedges,
    nodeattrtype, edgeattrtype

# Graph types

abstract type AbstractGraph end
abstract type UndirectedGraph <: AbstractGraph end
abstract type DirectedGraph <: AbstractGraph end
abstract type HyperGraph <: AbstractGraph end
abstract type OrderedGraph <: UndirectedGraph end
abstract type OrderedDiGraph <: DirectedGraph end
abstract type OrderedHyperGraph <: HyperGraph end

# TODO: use traits
# https://github.com/JuliaLang/julia/issues/2345

# Node and edge attributes

abstract type AbstractNode end
abstract type AbstractEdge end
abstract type UndirectedEdge <: AbstractEdge end # TODO: unnecessary?
abstract type DirectedEdge <: AbstractEdge end # TODO: unnecessary?

"""
    @cachefirst expression

A macro for molecular descriptor array cache mechanism.
"""
macro cachefirst(ex)
    func = ex.args[1].args[1]
    dummy = gensym()
    ex.args[1].args[1] = dummy
    return esc(quote
        $ex
        Core.@__doc__ function $func(graph; kwargs...)
            if isdefined(graph, :cache) && isempty(kwargs)
                # Return cache
                symf = nameof($func)
                symf in keys(graph.cache) && return graph.cache[symf]
            end
            return $dummy(graph; kwargs...)
        end
    end)
end

"""
    @cache expression

A macro that enables to force recalculate and set cache of descriptor arrays.
"""
macro cache(ex)
    func = ex.args[1]
    graph = ex.args[2]
    return esc(quote
        $graph.cache[nameof($func)] = $func($graph)
    end)
end

hascache(graph, key) = isdefined(graph, :cache) && haskey(graph.cache, key)

"""
    clearcache!(graph::AbstractGraph) -> nothing

Clear calculated property caches.

Calling `clearcache!` is recommended when the graph nodes/edges are added, removed or reindexed. You can [`Graph.clone`](@ref) the graph before destructive operation instead.
"""
function clearcache!(graph::AbstractGraph)
    empty!(graph.cache)
    return
end

"""
    clone(graph::AbstractGraph) -> AbstractGraph

Return deep copy of the graph.
"""
function clone end

# Lookup

"""
    neighbors(graph, i) -> Dict{Int,Int}

Return the mapping of incident edges and adjacent nodes of node `i`.
If the graph is directed graph, both outneighbors and inneighbors are mapped.
"""
neighbors(graph::UndirectedGraph, i::Int) = graph.neighbormap[i]
neighbors(graph::DirectedGraph, i::Int
    ) = merge(outneighbors(graph, i), inneighbors(graph, i))

neighbors(graph::HyperGraph, i::Int
    ) = Dict(inc => graph.edges[inc] for inc in graph.incidences[i])

"""
    outneighbors(graph::DirectedGraph, i::Int) -> Dict{Int,Int}

Return the mapping of successor node keys and out edge keys connected to
the given node.
"""
outneighbors(graph::DirectedGraph, i::Int) = graph.outneighbormap[i]

"""
    inneighbors(graph::DirectedGraph, i::Int) -> Dict{Int,Int}

Return the mapping of predecessor node keys and in edge keys connected to
the given node.
"""
inneighbors(graph::DirectedGraph, i::Int) = graph.inneighbormap[i]

function findedgekey(graph::UndirectedGraph, u::Int, v::Int)
    for (inc, adj) in neighbors(graph, u)
        adj == v && return inc
    end
end

function findedgekey(graph::DirectedGraph, source::Int, target::Int)
    for (inc, adj) in outneighbors(graph, source)
        adj == target && return inc
    end
end

function findedgekey(graph::HyperGraph, u::Int, v::Int)
    for (inc, adj) in neighbors(graph, u)
        v in adj && return inc
    end
end

function findalledgekeys(graph::UndirectedGraph, u::Int, v::Int)
    hasmultiedge(graph) || return [findedgekey(graph, u, v)]
    keys = Int[]
    for (inc, adj) in neighbors(graph, u)
        adj == v && push!(keys, inc)
    end
    return keys
end

function findalledgekeys(graph::DirectedGraph, source::Int, target::Int)
    hasmultiedge(graph) || return [findedgekey(graph, source, target)]
    keys = Int[]
    for (inc, adj) in outneighbors(graph, source)
        adj == target && push!(keys, inc)
    end
    return keys
end

function findalledgekeys(graph::HyperGraph, u::Int, v::Int)
    hasmultiedge(graph) || return [findedgekey(graph, u, v)]
    keys = Int[]
    for (inc, adj) in neighbors(graph, u)
        v in adj && push!(keys, inc)
    end
    return keys
end

"""
    getedge(graph::AbstractGraph, i::Int) -> Tuple{Int,Int}
    getedge(graph::HyperGraph, i::Int) -> Set{Int}

Return an edge.
"""
getedge(graph::AbstractGraph, i::Int) = graph.edges[i]

"""
    hasedge(graph::UndirectedGraph, u::Int, v::Int) -> Bool
    hasedge(graph::DirectedGraph, source::Int, target::Int) -> Bool

Return whether the given two nodes are connected by at least one edge.
"""
hasedge(graph::AbstractGraph, u::Int, v::Int
    ) = findedgekey(graph, u, v) !== nothing

"""
    nodeattr(graph::AbstractGraph, i::Int) -> AbstractNode

Return the attribute object of node `i`.
"""
nodeattr(graph::AbstractGraph, i::Int) = graph.nodeattrs[i]

"""
    edgeattr(graph::AbstractGraph, i::Int) -> AbstractEdge

Return the attribute object of edge `i`.
"""
edgeattr(graph::AbstractGraph, i::Int) = graph.edgeattrs[i]

"""
    edgeattr(graph::AbstractGraph, u::Int, v::Int
        ) -> Union{AbstractEdge,Nothing}

Return the attribute object of an edge that connects `u` and `v`. If not found,
return nothing.
"""
function edgeattr(graph::AbstractGraph, u::Int, v::Int)
    k = findedgekey(graph, u, v)
    return k === nothing ? nothing : graph.edgeattrs[k]
end

# Collection

adjacencies(g::AbstractGraph, i) = Set(values(neighbors(g, i)))
successors(g::DirectedGraph, i) = Set(values(outneighbors(g, i)))
predecessors(g::DirectedGraph, i) = Set(values(inneighbors(g, i)))

incidences(g::AbstractGraph, i) = Set(keys(neighbors(g, i)))
outincidences(g::DirectedGraph, i) = Set(keys(outneighbors(g, i)))
in_incidences(g::DirectedGraph, i) = Set(keys(inneighbors(g, i)))

"""
    nodeset(graph::Union{OrderedGraph,OrderedDiGraph}) -> Set{Int}

Return the set of node keys.
"""
nodeset(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}
    ) = Set(1:nodecount(graph))

"""
    edgeset(graph::Union{OrderedGraph,OrderedDiGraph}) -> Set{Int}

Return the set of edge keys.
"""
edgeset(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}
    ) = Set(1:edgecount(graph))

edgesiter(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}
    ) = graph.edges

neighborsiter(graph::Union{OrderedGraph}) = graph.neighbormap

"""
    nodeattrs(graph::Union{OrderedGraph,OrderedDiGraph}) -> Vector{AbstractNode}

Return graph node attributes.
"""
nodeattrs(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}
    ) = graph.nodeattrs

"""
    edgeattrs(graph::Union{OrderedGraph,OrderedDiGraph}) -> Vector{AbstractEdge}

Return graph edge attributes.
"""
edgeattrs(graph::Union{OrderedGraph,OrderedDiGraph,OrderedHyperGraph}
    ) = graph.edgeattrs

# Counting

"""
    nodecount(graph::AbstractGraph) -> Int

Return the number of graph nodes.
"""
nodecount(graph::UndirectedGraph) = length(graph.neighbormap)
nodecount(graph::DirectedGraph) = length(graph.outneighbormap)
nodecount(graph::HyperGraph) = length(graph.incidences)

"""
    edgecount(graph::AbstractGraph) -> Int

Return the number of graph edges.
"""
edgecount(graph::AbstractGraph) = length(graph.edges)

"""
    neighborcount(graph::AbstractGraph, n::Int) -> Int
    degree(graph::AbstractGraph, n::Int) -> Int

Return the number of adjacent nodes of the node 'n'.
"""
neighborcount(graph::AbstractGraph, n::Int) = length(neighbors(graph, n))
degree = neighborcount

"""
    outdegree(graph::DirectedGraph, n::Int) -> Int

Return the number of outneighbors of the node 'n'.
"""
outdegree(graph::DirectedGraph, n::Int) = length(outneighbors(graph, n))

"""
    indegree(graph::DirectedGraph, n::Int) -> Int

Return the number of inneighbors of the node 'n'.
"""
indegree(graph::DirectedGraph, n::Int) = length(inneighbors(graph, n))

# Editing

"""
    addnode!(graph) -> Int
    addnode!(graph, attr) -> Int

Add new node and return the node index. If the node attribute type is required,
specify the node attribute object by `node` keyword.
"""
function addnode!(graph::UndirectedGraph, attr::AbstractNode)
    push!(graph.neighbormap, Dict())
    push!(graph.nodeattrs, attr)
    return length(graph.neighbormap)
end

function addnode!(graph::UndirectedGraph)
    isdefined(graph, :nodeattrs) && throw(ErrorException("nodeattr required"))
    push!(graph.neighbormap, Dict())
    return length(graph.neighbormap)
end

function addnode!(graph::DirectedGraph, attr::AbstractNode)
    push!(graph.outneighbormap, Dict())
    push!(graph.inneighbormap, Dict())
    push!(graph.nodeattrs, attr)
    return length(graph.outneighbormap)
end

function addnode!(graph::DirectedGraph)
    isdefined(graph, :nodeattrs) && throw(ErrorException("nodeattr required"))
    push!(graph.outneighbormap, Dict())
    push!(graph.inneighbormap, Dict())
    return length(graph.outneighbormap)
end

function addnode!(graph::HyperGraph, attr::AbstractNode)
    push!(graph.incidences, Set())
    push!(graph.nodeattrs, attr)
    return length(graph.incidences)
end

function addnode!(graph::HyperGraph)
    isdefined(graph, :nodeattrs) && throw(ErrorException("nodeattr required"))
    push!(graph.incidences, Set())
    return length(graph.incidences)
end

"""
    addedge!(graph, u, v) -> Int
    addedge!(graph, u, v, attr) -> Int

Add new edge and return the edge index. If the edge attribute type is required,
specify the edge attribute object by `attr` keyword.
"""
function _addedge!(graph::UndirectedGraph, u::Int, v::Int)
    push!(graph.edges, (u, v))
    i = edgecount(graph)
    graph.neighbormap[u][i] = v
    graph.neighbormap[v][i] = u
    return i
end

function addedge!(graph::UndirectedGraph, u::Int, v::Int, attr::AbstractEdge)
    push!(graph.edgeattrs, attr)
    return _addedge!(graph, u, v)
end

function addedge!(graph::UndirectedGraph, u::Int, v::Int)
    isdefined(graph, :edgeattrs) && throw(ErrorException("edgeattr required"))
    return _addedge!(graph, u, v)
end

function _addedge!(graph::DirectedGraph, s::Int, t::Int)
    push!(graph.edges, (s, t))
    i = edgecount(graph)
    graph.outneighbormap[s][i] = t
    graph.inneighbormap[t][i] = s
    return i
end

function addedge!(graph::DirectedGraph, s::Int, t::Int, attr::AbstractEdge)
    push!(graph.edgeattrs, attr)
    return _addedge!(graph, s, t)
end

function addedge!(graph::DirectedGraph, s::Int, t::Int)
    isdefined(graph, :edgeattrs) && throw(ErrorException("edgeattr required"))
    return _addedge!(graph, s, t)
end

function _addedge!(graph::HyperGraph, edge::Set{Int})
    push!(graph.edges, edge)
    i = edgecount(graph)
    for n in edge
        push!(graph.incidences[n], i)
    end
    return i
end

function addedge!(graph::HyperGraph, edge::Set{Int}, attr::AbstractEdge)
    push!(graph.edgeattrs, attr)
    return _addedge!(graph, edge)
end

function addedge!(graph::HyperGraph, edge::Set{Int})
    isdefined(graph, :edgeattrs) && throw(ErrorException("edgeattr required"))
    return _addedge!(graph, edge)
end

"""
    setnodeattr!(graph::AbstractGraph, i::Int, attr::AbstractNode)

Update the node attribute.
"""
function setnodeattr!(graph::AbstractGraph, i::Int, attr::AbstractNode)
    graph.nodeattrs[i] = attr
    return
end

"""
    setedgeattr!(graph::AbstractGraph, i::Int, attr::AbstractNode)

Update the edge attribute.
"""
function setedgeattr!(graph::AbstractGraph, i::Int, attr::AbstractEdge)
    graph.edgeattrs[i] = attr
    return
end

"""
    unlinknodes(graph, nodes)

Delete given nodes and its incident edges from the graph.
"""
function unlinknodes end

"""
    unlinkedges(graph, edges)

Delete given edges from the graph.
"""
function unlinkedges end

# Types

"""
    nodeattrtype(graph) -> Type

Return node attribute type.
"""
nodeattrtype(mol::AbstractGraph) = eltype(mol.nodeattrs)

"""
    edgeattrtype(graph) -> Type

Return edge attribute type.
"""
edgeattrtype(mol::AbstractGraph) = eltype(mol.edgeattrs)


text:
using DataStreams, Test

df = (a=[1,2,3,1], b=["hey", "ho", "neighbor", "hi"], c=[4.0, 5.0, 6.6, 6.0], d=[0,0,0,0])

# sink = Data.Table
cell(x::Data.Table, row, col) = x[col][row]
cell(x::Data.RowTable, row, col) = x[row][col]

@testset "Data.query" begin
for sink in (Data.Table, Data.RowTable)
    # select *
    res = Data.query(df, [(col=1,), (col=2,), (col=3,), (col=4,)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String, Float64, Int)
    @test Data.header(sch) == ["a", "b", "c", "d"]
    @test size(sch) == (4, 4)
    @test cell(res, 1, 1) == 1
    @test cell(res, 1, 4) == 0
    @test cell(res, 4, 1) == 1
    @test cell(res, 4, 4) == 0

    # select *
    res = Data.query(df)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String, Float64, Int)
    @test Data.header(sch) == ["a", "b", "c", "d"]
    @test size(sch) == (4, 4)
    @test cell(res, 1, 1) == 1
    @test cell(res, 1, 4) == 0
    @test cell(res, 4, 1) == 1
    @test cell(res, 4, 4) == 0

    # select * offset 1
    res = Data.query(df, [(col=1,), (col=2,), (col=3,), (col=4,)], sink, offset=1)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String, Float64, Int)
    @test Data.header(sch) == ["a", "b", "c", "d"]
    @test size(sch) == (3, 4)

    # select * offset 1 limit 1
    res = Data.query(df, [(col=1,), (col=2,), (col=3,), (col=4,)], sink, offset=1, limit=1)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String, Float64, Int)
    @test Data.header(sch) == ["a", "b", "c", "d"]
    @test size(sch) == (1, 4)

    # select * offset 1 limit 5
    res = Data.query(df, [(col=1,), (col=2,), (col=3,), (col=4,)], sink, offset=1, limit=5)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String, Float64, Int)
    @test Data.header(sch) == ["a", "b", "c", "d"]
    @test size(sch) == (3, 4)

    # select * offset 4 limit 5
    res = Data.query(df, [(col=1,), (col=2,), (col=3,), (col=4,)], sink, offset=4, limit=5)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String, Float64, Int)
    @test Data.header(sch) == ["a", "b", "c", "d"]
    @test size(sch) == (0, 4)

    # select * where a = 2
    res = Data.query(df, [(col=1, filter=x->x==2), (col=2,), (col=3,)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String, Float64)
    @test Data.header(sch) == ["a", "b", "c"]
    @test size(sch) == (1, 3)

    # select a, b where a = 2
    res = Data.query(df, [(col=1, filter=x->x==2), (col=2,)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String)
    @test Data.header(sch) == ["a", "b"]
    @test size(sch) == (1, 2)

    # select b, a where a = 2
    res = Data.query(df, [(col=2,), (col=1, filter=x->x==2)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (String, Int)
    @test Data.header(sch) == ["b", "a"]
    @test size(sch) == (1, 2)

    # select c
    res = Data.query(df, [(col=3,)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Float64,)
    @test Data.header(sch) == ["c"]
    @test size(sch) == (4, 1)

    # select c where c < 6.0
    res = Data.query(df, [(col=3, filter=x->x < 6.0)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Float64,)
    @test Data.header(sch) == ["c"]
    @test size(sch) == (2, 1)

    # select a + 1
    res = Data.query(df, [(compute=x->x +1, computeargs=(1,))], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int,)
    @test Data.header(sch) == ["Column1"]
    @test size(sch) == (4, 1)

    # select c * 2.0 as c where c < 6.0
    res = Data.query(df, [(col=3, hide=true, filter=x->x < 6.0), (name=:c, compute=x->x * 2.0, computeargs=(3,))], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Float64,)
    @test Data.header(sch) == ["c"]
    @test size(sch) == (2, 1)

    # select c * 2.0 as c, c / 2.0 as c2 where c < 6.0
    res = Data.query(df, [(col=3, hide=true, filter=x->x < 6.0),
                    (name=:c, compute=x->x * 2.0, computeargs=(3,)),
                    (name=:c2, compute=x->x / 2.0, computeargs=(3,))], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Float64, Float64)
    @test Data.header(sch) == ["c", "c2"]
    @test size(sch) == (2, 2)

    # select a, sum(c) group by a
    res = Data.query(df, [(col=1, group=true), (col=3, aggregate=sum)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, Float64)
    @test Data.header(sch) == ["a", "c"]
    @test size(sch) == (3, 2)

    # select a, sum(c - d) group by a
    res = Data.query(df, [(col=1, group=true), (name=:c, compute=(x,y)->x - y, computeargs=(3,4), aggregate=sum)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, Float64)
    @test Data.header(sch) == ["a", "c"]
    @test size(sch) == (3, 2)

    # select a group by a
    res = Data.query(df, [(col=1, group=true)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int,)
    @test Data.header(sch) == ["a"]
    @test size(sch) == (3, 1)

    # select a, SUM(c) + SUM(d) as f group by a
    res = Data.query(df, [(col=1, group=true), (name=:f, computeaggregate=(c, d)->sum(c) + sum(d), computeargs=(3, 4))], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, Float64)
    @test Data.header(sch) == ["a", "f"]
    @test size(sch) == (3, 2)

    # select a, SUM(c) + SUM(d) as f where b like 'h%' group by a
    res = Data.query(df, [(col=1, group=true), (col=2, filter=x->x[1] == 'h', hide=true), (name=:f, computeaggregate=(c, d)->sum(c) + sum(d), computeargs=(3, 4))], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, Float64)
    @test Data.header(sch) == ["a", "f"]
    @test size(sch) == (2, 2)

    # select a, sum(c) group by a having sum(c) < 7.0
    res = Data.query(df, [(col=1, group=true), (col=3, aggregate=sum, having=x->x < 7.0)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, Float64)
    @test Data.header(sch) == ["a", "c"]
    @test size(sch) == (2, 2)

    # select a order by a
    res = Data.query(df, [(col=1, sort=true)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int,)
    @test Data.header(sch) == ["a"]
    @test size(sch) == (4, 1)

    # select a order by a desc
    res = Data.query(df, [(col=1, sort=true, sortasc=false)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int,)
    @test Data.header(sch) == ["a"]
    @test size(sch) == (4, 1)

    # select a, b order by a, b
    res = Data.query(df, [(col=1, sort=true), (col=2, sort=true)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String)
    @test Data.header(sch) == ["a", "b"]
    @test size(sch) == (4, 2)

    # select a, b order by a, b desc
    res = Data.query(df, [(col=1, sort=true), (col=2, sort=true, sortasc=false)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String)
    @test Data.header(sch) == ["a", "b"]
    @test size(sch) == (4, 2)

    # select a, b order by a desc, b desc
    res = Data.query(df, [(col=1, sort=true, sortasc=false), (col=2, sort=true, sortasc=false)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String)
    @test Data.header(sch) == ["a", "b"]
    @test size(sch) == (4, 2)

    # select a, b order by a desc, b asc
    res = Data.query(df, [(col=1, sort=true, sortasc=false), (col=2, sort=true)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String)
    @test Data.header(sch) == ["a", "b"]
    @test size(sch) == (4, 2)

    # select a, b order by b, a
    res = Data.query(df, [(col=1, sort=true, sortindex=2, sortasc=false), (col=2, sort=true, sortindex=1)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String)
    @test Data.header(sch) == ["a", "b"]
    @test size(sch) == (4, 2)

    # select a, b where a = 1 order by b desc, a
    res = Data.query(df, [(col=1, filter=x->x==1, sort=true, sortindex=2), (col=2, sort=true, sortindex=1, sortasc=false)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, String)
    @test Data.header(sch) == ["a", "b"]
    @test size(sch) == (2, 2)

    # select a, sum(c) group by a
    res = Data.query(df, [(col=1, group=true, sort=true), (col=3, aggregate=sum)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, Float64)
    @test Data.header(sch) == ["a", "c"]
    @test size(sch) == (3, 2)

    # select a, sum(c) group by a order by a having sum(c) < 7.0
    res = Data.query(df, [(col=1, group=true, sort=true), (col=3, aggregate=sum, having=x->x < 7.0)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, Float64)
    @test Data.header(sch) == ["a", "c"]
    @test size(sch) == (2, 2)

    # select a, sum(c) group by a order by b desc having sum(c) < 7.0
    res = Data.query(df, [(col=1, group=true), (col=3, aggregate=sum, sort=true, sortasc=false, having=x->x < 7.0)], sink)
    sch = Data.schema(res)
    @test Data.types(sch) == (Int, Float64)
    @test Data.header(sch) == ["a", "c"]
    @test size(sch) == (2, 2)
end
end # testset "Data.query"


text:
module EscapeAnalysis

import Core:
    CodeInfo,
    CodeInstance,
    MethodInstance,
    Argument,
    SSAValue,
    PiNode,
    PhiNode,
    UpsilonNode,
    PhiCNode,
    ReturnNode,
    GotoNode,
    GotoIfNot

const CC = Core.Compiler

import .CC:
    AbstractInterpreter,
    NativeInterpreter,
    WorldView,
    WorldRange,
    InferenceParams,
    OptimizationParams,
    get_world_counter,
    get_inference_cache,
    lock_mi_inference,
    unlock_mi_inference,
    add_remark!,
    may_optimize,
    may_compress,
    may_discard_trees,
    verbose_stmt_info,
    code_cache,
    get_inference_cache,
    OptimizationState,
    IRCode,
    optimize,
    widenconst,
    argextype

import Base.Meta:
    isexpr

import Base:
    destructure_callex

using InteractiveUtils

let __init_hooks__ = []
    global __init__, register_init_hook!
    __init__() = foreach(f->f(), __init_hooks__)
    register_init_hook!(@nospecialize(f)) = push!(__init_hooks__, f)
end

mutable struct EscapeAnalyzer{Info} <: AbstractInterpreter
    native::NativeInterpreter
    source::Union{Nothing,IRCode}
    info::Union{Nothing,Info}
end

CC.InferenceParams(interp::EscapeAnalyzer)    = InferenceParams(interp.native)
CC.OptimizationParams(interp::EscapeAnalyzer) = OptimizationParams(interp.native)
CC.get_world_counter(interp::EscapeAnalyzer)  = get_world_counter(interp.native)

CC.lock_mi_inference(::EscapeAnalyzer,   ::MethodInstance) = nothing
CC.unlock_mi_inference(::EscapeAnalyzer, ::MethodInstance) = nothing

CC.add_remark!(interp::EscapeAnalyzer, sv, s) = add_remark!(interp.native, sv, s)

CC.may_optimize(interp::EscapeAnalyzer)      = may_optimize(interp.native)
CC.may_compress(interp::EscapeAnalyzer)      = may_compress(interp.native)
CC.may_discard_trees(interp::EscapeAnalyzer) = may_discard_trees(interp.native)
CC.verbose_stmt_info(interp::EscapeAnalyzer) = verbose_stmt_info(interp.native)

CC.get_inference_cache(interp::EscapeAnalyzer) = get_inference_cache(interp.native)

# CC.code_cache(interp::EscapeAnalyzer) = code_cache(interp.native)

const GLOBAL_CODE_CACHE = IdDict{MethodInstance,CodeInstance}()
__clear_code_cache!() = empty!(GLOBAL_CODE_CACHE)

function CC.code_cache(interp::EscapeAnalyzer)
    worlds = WorldRange(get_world_counter(interp))
    return WorldView(GlobalCache(), worlds)
end

struct GlobalCache end

CC.haskey(wvc::WorldView{GlobalCache}, mi::MethodInstance) = haskey(GLOBAL_CODE_CACHE, mi)

CC.get(wvc::WorldView{GlobalCache}, mi::MethodInstance, default) = get(GLOBAL_CODE_CACHE, mi, default)

CC.getindex(wvc::WorldView{GlobalCache}, mi::MethodInstance) = getindex(GLOBAL_CODE_CACHE, mi)

function CC.setindex!(wvc::WorldView{GlobalCache}, ci::CodeInstance, mi::MethodInstance)
    GLOBAL_CODE_CACHE[mi] = ci
    add_callback!(mi) # register the callback on invalidation
    return nothing
end

function add_callback!(linfo)
    if !isdefined(linfo, :callbacks)
        linfo.callbacks = Any[invalidate_cache!]
    else
        if !any(@nospecialize(cb)->cb===invalidate_cache!, linfo.callbacks)
            push!(linfo.callbacks, invalidate_cache!)
        end
    end
    return nothing
end

function invalidate_cache!(replaced, max_world, depth = 0)
    delete!(GLOBAL_CODE_CACHE, replaced)

    if isdefined(replaced, :backedges)
        for mi in replaced.backedges
            mi = mi::MethodInstance
            if !haskey(GLOBAL_CODE_CACHE, mi)
                continue # otherwise fall into infinite loop
            end
            invalidate_cache!(mi, max_world, depth+1)
        end
    end
    return nothing
end

# analysis
# ========

"""
    abstract type EscapeInformation end

A lattice for escape information, which has the following elements:
- `NoInformation`: the top element of this lattice, meaning no information is derived
- `NoEscape`: the second topmost element of this lattice, meaning it will not escape from this local frame
- `ReturnEscape`: a lattice that is lower than `NoEscape`, meaning it will escape to the caller
- `Escape`: the bottom element of this lattice, meaning it will escape to somewhere

An abstract state will be initialized with the top(most) elements, and an escape analysis
will transition these elements from the top to the bottom.
"""
abstract type EscapeInformation end

struct NoInformation <: EscapeInformation end
struct Escape        <: EscapeInformation end
struct NoEscape      <: EscapeInformation end
struct ReturnEscape  <: EscapeInformation end

⊑(x::EscapeInformation, y::EscapeInformation) = x == y
⊑(::Escape,             ::EscapeInformation)  = true
⊑(::EscapeInformation,  ::NoInformation)      = true
⊑(::Escape,             ::NoInformation)      = true # avoid ambiguity
⊑(::ReturnEscape,       ::NoEscape)           = true

⊔(x::EscapeInformation, y::EscapeInformation) = x⊑y ? y : y⊑x ? x : NoInformation()
⊓(x::EscapeInformation, y::EscapeInformation) = x⊑y ? x : y⊑x ? y : Escape()

# extend lattices of escape information to lattices of mappings of arguments and SSA stmts to escape information
# ⊓ and ⊔ operate pair-wise, and from there we can just rely on the Base implementation for dictionary equality comparison
struct EscapeState
    arguments::Vector{EscapeInformation}
    ssavalues::Vector{EscapeInformation}
end
function EscapeState(nslots::Int, nargs::Int, nstmts::Int)
    arguments = EscapeInformation[
        i ≤ nargs ? ReturnEscape() : NoEscape() for i in 1:nslots]
    ssavalues = EscapeInformation[NoInformation() for _ in 1:nstmts]
    return EscapeState(arguments, ssavalues)
end
Base.copy(s::EscapeState) = EscapeState(copy(s.arguments), copy(s.ssavalues))

⊔(X::EscapeState, Y::EscapeState) = EscapeState(
    EscapeInformation[x ⊔ y for (x, y) in zip(X.arguments, Y.arguments)],
    EscapeInformation[x ⊔ y for (x, y) in zip(X.ssavalues, Y.ssavalues)])
⊓(X::EscapeState, Y::EscapeState) = EscapeState(
    EscapeInformation[x ⊓ y for (x, y) in zip(X.arguments, Y.arguments)],
    EscapeInformation[x ⊓ y for (x, y) in zip(X.ssavalues, Y.ssavalues)])
Base.:(==)(X::EscapeState, Y::EscapeState) = X.arguments == Y.arguments && X.ssavalues == Y.ssavalues

const GLOBAL_ESCAPE_CACHE = IdDict{MethodInstance,EscapeState}()
__clear_escape_cache!() = empty!(GLOBAL_ESCAPE_CACHE)

# An escape analysis implementation based on the algorithm described in the paper [MM02].
# The analysis works on the lattice of `EscapeInformation` and transitions lattice elements
# from the top to the bottom in a backward way, i.e. data flows from usage cites to definitions.
#
# [MM02] A Graph-Free approach to Data-Flow Analysis.
#        Markas Mohnen, 2002, April.
#        https://api.semanticscholar.org/CorpusID:28519618

# TODO
# - implement more builtin handling
# - (related to above) do alias analysis to some extent
# - maybe flow-sensitivity (with sparse analysis state)
function find_escapes(ir::IRCode, nargs::Int)
    (; stmts, sptypes, argtypes) = ir
    nstmts = length(stmts)
    state = EscapeState(length(ir.argtypes), nargs, nstmts) # flow-insensitive, only manage a single state

    while true
        local anyupdate = false
        local changes = Pair{Any,EscapeInformation}[]

        for pc in nstmts:-1:1
            stmt = stmts.inst[pc]

            # collect escape information
            if isa(stmt, Expr)
                head = stmt.head
                if head === :call
                    ft = widenconst(argextype(first(stmt.args), ir, sptypes, argtypes))
                    # TODO implement more builtins, make them more accurate
                    if ft === Core.IntrinsicFunction # XXX we may break soundness here, e.g. `pointerref`
                        continue
                    elseif ft === typeof(isa) || ft === typeof(typeof) || ft === typeof(Core.sizeof)
                        continue
                    elseif ft === typeof(getfield) || ft === typeof(tuple)
                        info = state.ssavalues[pc]
                        info === NoInformation() && (info = NoEscape())
                        for arg in stmt.args[2:end]
                            push!(changes, arg => info)
                        end
                    else
                        for arg in stmt.args[2:end]
                            push!(changes, arg => Escape())
                        end
                    end
                elseif head === :invoke
                    linfo = first(stmt.args)
                    escapes_for_call = get(GLOBAL_ESCAPE_CACHE, linfo, nothing)
                    if isnothing(escapes_for_call)
                        for arg in stmt.args[3:end]
                            push!(changes, arg => Escape())
                        end
                    else
                        for (arg, info) in zip(stmt.args[2:end], escapes_for_call.arguments)
                            push!(changes, arg => info)
                        end
                    end
                elseif head === :new
                    info = state.ssavalues[pc]
                    info === NoInformation() && (info = NoEscape())
                    for arg in stmt.args[2:end]
                        push!(changes, arg => info)
                    end
                elseif head === :(=)
                    lhs, rhs = stmt.args
                    if isa(lhs, GlobalRef)
                        push!(changes, rhs => Escape())
                    end
                elseif head === :enter || head === :leave || head === :pop_exception
                    continue
                else # TODO: this is too conservative
                    for arg in stmt.args
                        push!(changes, arg => Escape())
                    end
                end
            elseif isa(stmt, PiNode)
                if isdefined(stmt, :val)
                    info = state.ssavalues[pc]
                    push!(changes, stmt.val => info)
                end
            elseif isa(stmt, PhiNode)
                info = state.ssavalues[pc]
                values = stmt.values
                for i in 1:length(values)
                    if isassigned(values, i)
                        push!(changes, values[i] => info)
                    end
                end
            elseif isa(stmt, PhiCNode)
                info = state.ssavalues[pc]
                values = stmt.values
                for i in 1:length(values)
                    if isassigned(values, i)
                        push!(changes, values[i] => info)
                    end
                end
            elseif isa(stmt, UpsilonNode)
                if isdefined(stmt, :val)
                    info = state.ssavalues[pc]
                    push!(changes, stmt.val => info)
                end
            elseif isa(stmt, ReturnNode)
                if isdefined(stmt, :val)
                    push!(changes, stmt.val => ReturnEscape())
                end
            else
                @assert stmt isa GotoNode || stmt isa GotoIfNot || stmt isa GlobalRef || stmt === nothing # TODO remove me
                continue
            end

            isempty(changes) && continue

            # propagate changes
            new = copy(state)
            for (x, info) in changes
                if isa(x, Argument)
                    new.arguments[x.n] = new.arguments[x.n] ⊓ info
                elseif isa(x, SSAValue)
                    new.ssavalues[x.id] = new.ssavalues[x.id] ⊓ info
                end
            end
            empty!(changes)

            # convergence check and worklist update
            if new != state
                state = new

                anyupdate |= true
            end
        end

        anyupdate || break
    end

    return state
end

# entries
# =======

function CC.optimize(interp::EscapeAnalyzer, opt::OptimizationState, params::OptimizationParams, @nospecialize(result))
    nargs = Int(opt.nargs) - 1
    ir = run_passes_with_escape_analysis(interp, opt.src, nargs, opt)
    return CC.finish(interp, opt, params, ir, result)
end

# HACK enable copy and paste from Core.Compiler
function run_passes_with_escape_analysis end
register_init_hook!() do
@eval CC begin
    function $EscapeAnalysis.run_passes_with_escape_analysis(interp::$EscapeAnalyzer, ci::CodeInfo, nargs::Int, sv::OptimizationState)
        preserve_coverage = coverage_enabled(sv.mod)
        ir = convert_to_ircode(ci, copy_exprargs(ci.code), preserve_coverage, nargs, sv)
        ir = slot2reg(ir, ci, nargs, sv)
        #@Base.show ("after_construct", ir)
        # TODO: Domsorting can produce an updated domtree - no need to recompute here
        @timeit "compact 1" ir = compact!(ir)
        @timeit "Inlining" ir = ssa_inlining_pass!(ir, ir.linetable, sv.inlining, ci.propagate_inbounds)
        #@timeit "verify 2" verify_ir(ir)
        ir = compact!(ir)
        @timeit "collect escape information" escapes = $find_escapes(ir, nargs+1)
        $setindex!($GLOBAL_ESCAPE_CACHE, escapes, sv.linfo)
        interp.source = copy(ir)
        interp.info = escapes
        #@Base.show ("before_sroa", ir)
        @timeit "SROA" ir = getfield_elim_pass!(ir)
        #@Base.show ir.new_nodes
        #@Base.show ("after_sroa", ir)
        ir = adce_pass!(ir)
        #@Base.show ("after_adce", ir)
        @timeit "type lift" ir = type_lift_pass!(ir)
        @timeit "compact 3" ir = compact!(ir)
        #@Base.show ir
        if JLOptions().debug_level == 2
            @timeit "verify 3" (verify_ir(ir); verify_linetable(ir.linetable))
        end
        return ir
    end
end
end # register_init_hook!() do

macro analyze_escapes(ex0...)
    return InteractiveUtils.gen_call_with_extracted_types_and_kwargs(__module__, :analyze_escapes, ex0)
end

function analyze_escapes(@nospecialize(f), @nospecialize(types=Tuple{});
                         world = get_world_counter(),
                         interp = Core.Compiler.NativeInterpreter(world))
    interp = EscapeAnalyzer{EscapeState}(interp, nothing, nothing)

    code_typed(f, types; optimize=true, world, interp)
    return interp.source, interp.info
end

# in order to run a whole analysis from ground zero (e.g. for benchmarking, etc.)
__clear_caches!() = (__clear_code_cache!(); __clear_escape_cache!())

export
    analyze_escapes,
    @analyze_escapes

end # module EscapeAnalysis


text:
function Integrate(nx, xx, yy, adj)

    t = 0.0f0
    if (adj)
	for i = nx : -1 : 1
	    t += yy[i]
	    xx[i] += t;
	end
    else
	for i = 1:nx
	    t += xx[i]
	    yy[i] += t;
	end
    end

end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
abstract Model # this is either an LSTM or RNN

type NNMatrix # Neural net layer's weights & gradients
    n::Int
    d::Int
    w::Matrix{Float64} # matrix of weights
    dw::Matrix{Float64} # matrix of gradients
    NNMatrix(n::Int) = new(n, 1, zeros(n,1), zeros(n,1))
    NNMatrix(n::Int, d::Int) = new(n, d, zeros(n,d), zeros(n,d))
    NNMatrix(n::Int, d::Int, w::Matrix{Float64}, dw::Matrix{Float64}) = new(n, d, w, dw)
    NNMatrix(w::Array) = new(size(w,1), size(w,2), w, zeros(size(w,1),size(w,2)))
    function NNMatrix(w::Matrix{Float64}, dw::Matrix{Float64})
      if size(w) != size(dw)
        error("Matrix and Gradient Matrix sizes do not match");
      end
      new(size(w,1), size(w,2), w, dw)
    end
end

randNNMat(n::Int, d::Int, std::FloatingPoint=1.) = NNMatrix(n, d, randn(n,d)*std, zeros(n,d))
onesNNMat(n::Int, d::Int) = NNMatrix(n, d, ones(n, d), zeros(n, d))

function softmax(m::NNMatrix)
    out = NNMatrix(m.n,m.d)
    maxval = maximum(m.w)
    out.w[:] = exp(m.w - maxval)
    out.w /= sum(out.w)
    return out
end

max_forks_repo_name: haroldsoh/RecurrentNN


text:
# evaluate and understand the 4 types of linear phase FIR filter responses.

# from type 1, degree L to type 2, degree L.
# 𝑎 here is b_tilde from Flilip's thesis Section 2.2.3.
# a here is b from thesis.
function converttotype2(𝑎::Vector{T}) where T
    L = length(𝑎)

    a = Vector{T}(undef,L)

    a[1] = 𝑎[1] + 𝑎[2]/2
    for k = 2:L-1
        a[k] = (𝑎[k-1+1]+𝑎[k+1])/2
    end

    a[L] = 𝑎[L]/2

    return a
end

# 𝑎 is c_tilde from Flilip's thesis Section 2.2.3.
# a here is c from thesis.
function converttotype3(𝑎::Vector{T}) where T
    L = length(𝑎)
    @assert L >= 3

    a = Vector{T}(undef,L)

    a[1] = 𝑎[1] - 𝑎[3]/2
    for k = 2:L-2
        a[k] = (𝑎[k-1+1]-𝑎[k+1+1])/2
    end

    a[L-1] = 𝑎[L-1]/2
    a[L] = 𝑎[L]/2

    return a
end

# 𝑎 is d_tilde from Flilip's thesis Section 2.2.3.
# a here is d from thesis.
function converttotype4(𝑎::Vector{T}) where T
    L = length(𝑎)

    a = Vector{T}(undef,L)

    a[1] = 𝑎[1] - 𝑎[2]/2
    for k = 2:L-1
        a[k] = (𝑎[k-1+1]-𝑎[k+1])/2
    end

    a[L] = 𝑎[L]/2

    return a
end

function getfilterfuncs(type_num::Int)
    @assert 1 <= type_num <= 4

    getafunc = getatype1
    gethfunc = gethtype1

    if type_num == 2
        getafunc = getatype2
        gethfunc = gethtype2

    elseif type_num == 3
        getafunc = getatype3
        gethfunc = gethtype3

    elseif type_num == 4
        getafunc = getatype4
        gethfunc = gethtype4
    end

    return getafunc, gethfunc
end

#### type 1
function gethtype1(a::Vector{T}) where T
    L = length(a) -1
    N = 2*L+1

    h = Vector{T}(undef, N)
    fill!(h,-1.23)

    #tmp = Vector{T}(undef, L+1)
    h[L+1] = a[1]
    for k = 1:L
        h[L-k+1] = a[k+1]/2
    end

    h[L+2:end] = reverse(h[1:L])

    return h
end

function getatype1(h::Vector{T}) where T
    N = length(h)
    @assert isodd(length(h))

    L = div(N-1,2)

    a = Vector{T}(undef,L+1)
    fill!(a,-1.23)

    a[1] = h[L+1]
    for k = 1:L
        a[k+1] = 2*h[L-k+1]
    end

    return a
end

function evalHeqn211(h::Vector{T}, ω ) where T

    a = getatype1(h)
    L = length(a)-1

    return sum( a[k+1]*cos(k*ω) for k = 0:L)
end

#### type 2

function evalHeqn212(h::Vector{T}, ω ) where T

    a = getatype2(h)
    L = length(a)

    return sum( a[k]*cos((k-0.5)*ω) for k = 1:L)
end

function getatype2(h::Vector{T}) where T

    N = length(h)
    @assert iseven(N)

    L = div(N,2)

    a = Vector{T}(undef,L)
    fill!(a,-1.23)

    for k = 1:L
        a[k] = 2*h[L-k+1]
    end

    return a
end

function gethtype2(a::Vector{T}) where T
    L = length(a)
    N = 2*L

    h = Vector{T}(undef, N)
    fill!(h,-1.23)

    for k = 1:L
        h[L-k+1] = a[k]/2
    end

    h[L+1:end] = reverse(h[1:L])

    return h
end

#### type 3

function evalHeqn213(h::Vector{T}, ω ) where T

    a = getatype3(h)
    L = length(a)

    return sum( a[k]*sin(k*ω) for k = 1:L)
end

function gethtype3(a::Vector{T}) where T
    L = length(a)
    N = 2*L+1

    h = Vector{T}(undef, N)
    fill!(h,-1.23)

    #tmp = Vector{T}(undef, L+1)
    h[L+1] = zero(T)
    for k = 1:L
        h[L-k+1] = a[k]/2
    end

    h[L+2:end] = reverse(-h[1:L])

    return h
end

function getatype3(h::Vector{T}) where T
    N = length(h)
    @assert isodd(length(h))

    L = div(N-1,2)
    @assert isnumericallyclose(h[L+1],zero(T))

    a = Vector{T}(undef,L)
    fill!(a,-1.23)

    a[1] = zero(T)
    for k = 1:L
        a[k] = 2*h[L-k+1]
    end

    return a
end

#### type 4

function evalHeqn214(h::Vector{T}, ω ) where T

    a = getatype4(h)
    L = length(a)

    return sum( a[k]*sin((k-0.5)*ω) for k = 1:L )
end

function getatype4(h::Vector{T}) where T

    N = length(h)
    @assert iseven(N)

    L = div(N,2)

    a = Vector{T}(undef,L)
    fill!(a,-1.23)

    for k = 1:L
        a[k] = 2*h[L-k+1]
    end

    return a
end

function gethtype4(a::Vector{T}) where T
    L = length(a)
    N = 2*L

    h = Vector{T}(undef, N)
    fill!(h,-1.25)

    for k = 1:L
        h[L-k+1] = a[k]/2
    end

    h[L+1:end] = reverse(-h[1:L])

    return h
end

#### generic for all 4 types.

function evalHeqn210(h,ω,type_num)
    L = 0
    if type_num == 3 || type_num == 4
        # case negative symmetric.
        L = 1
    end

    N = length(h)
    G = naiveimpulsersp(h,ω)

    return G*exp( -im*(L*π - (N-1)*ω)/2 )
end

function naiveimpulsersp(h,ω)
    N = length(h)
    return sum( h[n+1]*exp(-im*ω*n) for n = 0:N-1 )
end

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
# Autogenerated wrapper script for CURL_jll for armv7l-linux-gnueabihf
export curl

using LibSSH2_jll
using Zlib_jll
using nghttp2_jll
using MbedTLS_jll
using LibCURL_jll
JLLWrappers.@generate_wrapper_header("CURL")
JLLWrappers.@declare_executable_product(curl)
function __init__()
    JLLWrappers.@generate_init_header(LibSSH2_jll, Zlib_jll, nghttp2_jll, MbedTLS_jll, LibCURL_jll)
    JLLWrappers.@init_executable_product(
        curl,
        "bin/curl",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using MyPackage
using Test

@testset "MyPackage.jl" begin
    # Write your tests here.
end


text:
using Requires
using AutoInteract
using Interact
using Plots
using Base.Test
using Interact.update_view

# TODO: properly test AxisArrays and other complex implementations.

function mocktoggle(value, label)
    widget = togglebutton(label, value = value)
    widget.signal.name = label * "_toggle"
    return widget
end

function mockslider(value, label, range)
    widget = slider(range, value = value, label = label)
    widget.signal.name = label * "_slider"
    return widget
end

# Override update method to do nothing on push!
Interact.update_view(::Any) = nothing

# Signal numbers will be different every time, so we compare properties
function compare_slider(widget1, widget2)
    @test widget1.label == widget2.label
    @test widget1.value == widget2.value
    @test widget1.range == widget2.range
    @test widget1.orientation == widget2.orientation
    @test widget1.signal.name == widget2.signal.name
end

function compare_toggle(widget1, widget2)
    @test widget1.label == widget2.label
    @test widget1.value == widget2.value
    @test widget1.signal.name == widget2.signal.name
end

# Test components to escape circular references
function compare_layouts(l1, l2)
    isequal(l1.layout.box.children, l2.layout.box.children)
    isequal(l1.layout.box.vert, l2.layout.box.vert)
    isequal(l1.vert, l2.vert)
end

mutable struct teststruct
    a::Array{Float64}
    b::Bool
    i::Int
    f::Float64
end

@testset "interface" begin

    @require Unitful begin
        @testset "interface from Unitful Unit" begin
            slider = make_widgets(1.0u"g", "Unit")
            @test slider.label == "Unit"
            @test slider.value == 1.0u"g"
            @test slider.range == 0.0u"g":0.02u"g":2.0u"g"
            @test slider.signal.name == "Unit_slider"
        end
    end

    @testset "interface from Float64" begin
        ref_slider = mockslider(1.0, "Float64", 0.0:0.02:2.0)
        slider = make_widgets(1.0, "Float64")
        compare_slider(ref_slider, slider)
        # Interface
        iface = make_interface(slider)
        # No vboxes are added for single sliders
        ref_iface = slider
        compare_slider(ref_iface, iface)
    end

    @testset "interface from Bool" begin
        ref_toggle = mocktoggle(true, "Bool")
        toggle = make_widgets(true, "Bool")
        compare_toggle(ref_toggle, toggle)
        # Interface
        iface = make_interface(toggle)
        # No vboxes are added for single sliders
        ref_iface = toggle
        compare_toggle(ref_iface, iface)
    end

    @testset "interface from Int" begin
        # Positive
        ref_slider = mockslider(1, "Int", 0:1:2)
        slider = make_widgets(1, "Int")
        compare_slider(ref_slider, slider)
        # Negative
        ref_slider = mockslider(-2, "Int", -4:1:0)
        slider = make_widgets(-2, "Int")
        compare_slider(ref_slider, slider)
        # Custom range
        ref_slider = mockslider(9, "Int", -99:9:99)
        slider = make_widgets(9, "Int", -99:9:99)
        compare_slider(ref_slider, slider)

        # Interface
        iface = make_interface(slider)
        # No vboxes are added for single sliders
        ref_iface = slider
        compare_slider(ref_iface, iface)
    end

    @testset "interface from tuple" begin
        ref_sliders = mockslider.((1.0, 2.0, 3.0), ("unknown1", "unknown2", "unknown3"), (0.0:0.02:2.0, 0.0:0.04:4.0, 0.0:0.06:6.0))
        sliders = make_widgets((1.0, 2.0, 3.0))
        compare_slider.(ref_sliders, sliders)
        iface = make_interface(sliders)
        ref_iface = vbox(sliders...)
    end

    @testset "Arrays < maxlen are ignored" begin
        set_defaults(:maxlen_widget_array => 2)
        ref_sliders = nothing
        sliders = make_widgets([1.0, 2.0, 3.0])
        @test ref_sliders == sliders
    end

    @testset "Arrays > minlen make sliders" begin
        set_defaults(:maxlen_widget_array => 3)
        ref_sliders = mockslider.((1.0, 2.0, 3.0), ("s1", "s2", "s3"), (0.0:0.02:2.0, 0.0:0.04:4.0, 0.0:0.06:6.0))
        sliders = make_widgets([1.0, 2.0, 3.0], "s")
        compare_slider.(ref_sliders, sliders)
    end

    @testset "interface from struct" begin
        set_defaults(:maxlen_widget_array => 0)
        data = teststruct([1.0,2.0,3.0], true, 10, 10.0)
        b = mocktoggle(true, "b")
        i = mockslider(10, "i", 0:1:20)
        f = mockslider(10.0, "f", 0.0:0.2:20.0)
        widgets = make_widgets(data)
        @test collect(keys(widgets)) == [:b, :i, :f]
        compare_toggle(b, widgets[:b])
        compare_slider(i, widgets[:i])
        compare_slider(f, widgets[:f])

        iface = make_interface(widgets)
        ref_iface = vbox(widgets[:b], widgets[:i], widgets[:f])

        # Test components to escape circular references
        isequal(iface.layout.box.children, ref_iface.layout.box.children)
        isequal(iface.layout.box.vert, ref_iface.layout.box.vert)
        isequal(iface.vert, ref_iface.vert)
    end
end

@testset "plottables" begin
    plotly()
    array = collect(1.0:100.0)
    data = teststruct(array, true, 10, 10.0)
    b = mocktoggle(true, "b")
    i = mockslider(10, "i", 0:1:21)
    f = mockslider(10.0, "f", 0.0:0.2:20.0)
    plottables = make_plottables(data)
    signals = get_signals(plottables)

    @testset "Plottable heriarchy only contains plottables objects" begin
        @test collect(keys(plottables)) == [:a]
    end

    @testset "Signals are initially false" begin
        @test signals.value == Dict(:a => false)
    end

    @testset "No plots if plottables signals are false" begin
        @test typeof(plot_all(data, signals.value)) == Void
    end

    @testset "Plots an array if its signal is true" begin
        signals.value[:a] = true
        p = plot_all(data, signals.value)
        @test length(p.subplots) == 1
        @test p.subplots[1].series_list[1].d[:y] == array
    end

    @testset "Makes widget for arrays larger than > minlen" begin
        set_defaults(:minlen_plottable_array => 3)
        @test make_plottables([1, 2], "too small") == nothing
        compare_toggle(make_plottables([1, 2, 3], "ok"), mocktoggle(false, "ok"))
    end
end

@testset "can apply modified signals back to data" begin
    data = teststruct([1.0,2.0,3.0], true, 10, 10.0)
    widgets = make_widgets(data)
    signals = get_signals(widgets)

    @testset "initial values are correct" begin
        @test data.b == true
        @test data.i == 10
        @test data.f == 10.0
    end

    push!(widgets[:b].signal, false)
    push!(widgets[:i].signal, 99)
    push!(widgets[:f].signal, 99.9)
    sleep(0.01) # Wait for async updates
    apply_signals!(data, signals.value)

    @testset "updates have been aplied" begin
        @test data.b == false
        @test data.i == 99
        @test data.f == 99.9
    end
end


text:
# SPDX-License-Identifier: MIT

using CUDA, CUDA.CUSPARSE, SparseArrays
import LogExpFunctions: logsumexp
using MarkovModels
import MarkovModels: logaddexp
using Test

const D = 5 # vector/matrix dimension
const B = 2 # batch size
const S = 3 # number of states
const N = 5 # number of frames

const T = Float64
const SF = LogSemifield{T}

function makefsm(SF, S)
    fsm = FSM{SF}()

    prev = addstate!(fsm, pdfindex = 1)
    setinit!(prev)
    link!(fsm, prev, prev)
    for s in 2:S
        state = addstate!(fsm, pdfindex = s)
        link!(fsm, prev, state)
        link!(fsm, state, state)
        prev = state
    end
    setfinal!(prev)
    renormalize!(fsm)
    fsm
end

function forward(A, init, lhs::AbstractMatrix{T}) where T
    N = size(lhs, 2)
    log_α = similar(lhs)
    log_α[:, 1] = lhs[:, 1] + init
    for n in 2:N
        log_α[:, n] = lhs[:,n] + dropdims(logsumexp(A .+ log_α[:,n-1], dims = 1),
                                          dims = 1)
    end
    log_α
end

function backward(Aᵀ, final, lhs::AbstractMatrix{T}) where T
    N = size(lhs, 2)
    log_β = similar(lhs)
    log_β[:, end] = final
    for n in N-1:-1:1
        log_β[:, n] = dropdims(logsumexp(Aᵀ .+ log_β[:,n+1] .+ lhs[:,n+1], dims = 1),
                               dims = 1)
    end
    log_β
end

function forward_backward(A, Aᵀ, init, final, lhs)
    log_α = forward(A, init, lhs)
    log_β = backward(Aᵀ, final, lhs)
    log_γ = log_α .+ log_β
    sums = logsumexp(log_γ, dims = 1)
    exp.(log_γ .- sums), minimum(sums)
end

@testset "logaddexp" begin
    @test (@inferred logaddexp(2.0, 3.0)) ≈ log(exp(2.0) + exp(3.0))
    @test (@inferred logaddexp(10002.0, 10003.0)) ≈ 10000 + logaddexp(2.0, 3.0)
end

if CUDA.functional()
    @testset "linalg" begin
        D = 5

        x = ones(T, D)
        y = spzeros(T, D)
        y[2] = T(3)
        xg = CuArray(x)
        yg = CuSparseVector(y)
        out = similar(xg)
        MarkovModels.elmul_svdv!(out, yg, xg)
        @test all(convert(Vector{T}, out) .≈ x .* y)

        x = ones(T, D)
        y = spzeros(T, D)
        xg = CuArray(x)
        yg = CuSparseVector(y)
        out = similar(xg)
        MarkovModels.elmul_svdv!(out, yg, xg)
        @test all(convert(Vector{T}, out) .≈ x .* y)

        for Yᵀ in [spdiagm(y), spzeros(T, length(y))]
            Yᵀ = spdiagm(y)
            Yᵀg = CuSparseMatrixCSC(Yᵀ)
            Yg = CuSparseMatrixCSR( Yᵀg.colPtr, Yᵀg.rowVal, Yᵀg.nzVal, (D,D))
            out = similar(xg)
            MarkovModels.mul_smdv!(out, Yg, xg)
            @test all(convert(Vector{T}, out) .≈ transpose(Yᵀ) * x)
        end

        x = ones(T, D)
        y = spzeros(T, D)
        X = ones(T, D, D)
        Y = spzeros(T, D, D)
        xg = CuArray(x)
        yg = CuSparseVector(y)
        Xg = CuArray(X)
        Yg = CuSparseMatrixCSC(Y)
        Yg = CuSparseMatrixCSR(Yg.colPtr, Yg.rowVal, Yg.nzVal, Yg.dims)
        out = similar(Xg)
        MarkovModels.elmul_svdm!(out, yg, Xg)
        @test all(convert(Matrix{T}, out) .≈ X .* reshape(y, :, 1))

        x = ones(T, D)
        y = spzeros(T, D)
        X = ones(T, D, D)
        Y = spzeros(T, D, D)
        Y[1,1] = T(1)
        Y[1,2] = T(2)
        Y[2,1] = T(3)
        xg = CuArray(x)
        yg = CuSparseVector(y)
        Xg = CuArray(X)
        Yg = CuSparseMatrixCSC(Y)
        Yg = CuSparseMatrixCSR(Yg.colPtr, Yg.rowVal, Yg.nzVal, Yg.dims)
        out = similar(Xg)
        MarkovModels.elmul_svdm!(out, yg, Xg)
        @test all(convert(Matrix{T}, out) .≈ X .* reshape(y, :, 1))

        A = spzeros(3, 3)
        Ag = CuSparseMatrixCSC(A)
        Ag = CuSparseMatrixCSR(Ag.colPtr, Ag.rowVal, Ag.nzVal, Ag.dims)
        Y = ones(3, 3)
        Yg = CuArray(Y)
        out = similar(Yg)
        MarkovModels.mul_smdm!(out, Ag, Yg)
        @test all(convert(Matrix{T}, out) .≈ A' * Y)

        A = spzeros(3, 3)
        A[1, 1] = 1/2
        A[1, 2] = 1/2
        A[2, 2] = 1/2
        A[2, 3] = 1/2
        A[3, 3] = 1/2
        Ag = CuSparseMatrixCSC(A)
        Ag = CuSparseMatrixCSR(Ag.colPtr, Ag.rowVal, Ag.nzVal, Ag.dims)
        Y = ones(3, 3)
        Yg = CuArray(Y)
        out = similar(Yg)
        MarkovModels.mul_smdm!(out, Ag, Yg)
        @test all(convert(Matrix{T}, out) .≈ A' * Y)

        A = ones(T, 3, 2)
        Ag = CuArray(A)
        B = ones(T, 2, 3)
        Bg = CuArray(B)
        out = similar(Ag, 3, 3)
        MarkovModels.mul_dmdm!(out, Ag, Bg)
        @test all(convert(Matrix{T}, out) .≈ A * B)
    end
end

@testset "forward_backward" begin
    lhs = ones(T, S, N)

    fsm = makefsm(SF, S)
    cfsm = compile(fsm, allocator = zeros)

    γ_ref, ttl_ref = forward_backward(
        convert(Matrix{T}, cfsm.A),
        convert(Matrix{T}, cfsm.Aᵀ),
        convert(Vector{T}, cfsm.π),
        convert(Vector{T}, cfsm.ω),
        convert(Matrix{T}, lhs)
    )

    γ_dcpu, ttl_dcpu = @inferred stateposteriors(cfsm, lhs)
    @test all(γ_ref .≈ γ_dcpu)
    @test ttl_ref ≈ ttl_dcpu

    cfsm = compile(fsm, allocator = spzeros)
    γ_scpu, ttl_scpu = @inferred stateposteriors(cfsm, lhs)
    @test all(γ_ref .≈ γ_scpu)
    @test ttl_ref ≈ ttl_scpu

    if CUDA.functional()
        cfsm = compile(fsm, allocator = zeros) |> gpu
        γ_dgpu, ttl_dgpu = @inferred stateposteriors(cfsm, CuArray(lhs))
        @test all(γ_ref .≈ convert(Matrix{T}, γ_dgpu))
        @test ttl_ref ≈ ttl_dgpu

        cfsm = compile(fsm, allocator = spzeros) |> gpu
        γ_sgpu, ttl_sgpu = @inferred stateposteriors(cfsm, CuArray(lhs))
        @test all(γ_ref .≈ convert(Matrix{T}, γ_sgpu))
        @test ttl_ref ≈ ttl_sgpu
    end
end

@testset "batch forward_backward" begin
    lhs = ones(T, S, N, B)

    fsm = makefsm(SF, S)
    cfsm = compile(fsm, allocator = zeros)

    γ_ref, ttl_ref = forward_backward(
        convert(Matrix{T}, cfsm.A),
        convert(Matrix{T}, cfsm.Aᵀ),
        convert(Vector{T}, cfsm.π),
        convert(Vector{T}, cfsm.ω),
        convert(Matrix{T}, lhs[:, :, 1])
    )

    γ_dcpu, ttl_dcpu = @inferred stateposteriors(cfsm, lhs)
    for b in 1:B
        @test all(γ_ref .≈ γ_dcpu[:,:,b])
        @test ttl_ref ≈ ttl_dcpu[b]
    end

    cfsm = compile(fsm, allocator = spzeros)
    γ_scpu, ttl_scpu = @inferred stateposteriors(cfsm, lhs)
    for b in 1:B
        @test all(γ_ref .≈ γ_scpu[:,:,b])
        @test ttl_ref ≈ ttl_scpu[b]
    end

    if CUDA.functional()
        cfsm = compile(fsm, allocator = zeros) |> gpu
        γ_dgpu, ttl_dgpu = @inferred stateposteriors(cfsm, CuArray(lhs))
        for b in 1:B
            @test all(γ_ref .≈ convert(Array{T,3}, γ_dgpu)[:,:,b])
            @test ttl_ref ≈ convert(Vector{T}, ttl_dgpu)[b]
        end

        cfsm = compile(fsm, allocator = spzeros) |> gpu
        γ_sgpu, ttl_sgpu = @inferred stateposteriors(cfsm, CuArray(lhs))
        for b in 1:B
            @test all(γ_ref .≈ convert(Array{T,3}, γ_sgpu)[:,:,b])
            @test ttl_ref ≈ convert(Vector{T}, ttl_sgpu)[b]
        end
    end
end

max_stars_repo_name: MartinKocour/MarkovModels

max_issues_repo_name: MartinKocour/MarkovModels

max_forks_repo_name: MartinKocour/MarkovModels