text:
using Test, QuakeML
using Dates: DateTime

@testset "Accessors" begin
    @testset "has_" begin
        let event = QuakeML.Event(public_id="smi:QuakeML.jl/event/a")
            @test !has_focal_mechanism(event)
            @test !has_magnitude(event)
            @test !has_origin(event)
            push!(event.focal_mechanism,
                QuakeML.FocalMechanism(public_id="smi:QuakeML.jl/focal_mechanism/a"))
            push!(event.magnitude,
                    QuakeML.Magnitude(public_id="smi:QuakeML.jl/magnitude/a", mag=1))
            push!(event.origin, QuakeML.Origin(public_id="smi:QuakeML.jl/origin/a",
                latitude=1, longitude=1, time=Dates.now()))
            @test has_focal_mechanism(event)
            @test has_magnitude(event)
            @test has_origin(event)
        end
    end

    @testset "preferred_origin" begin
        let qml = QuakeML.read(joinpath(@__DIR__, "data", "nepal_mw7.2.qml"))
            @test preferred_origin(qml.event[1]).public_id ==
                qml.event[1].preferred_origin_id
            @test preferred_origin(qml.event[1]) == qml.event[1].origin[2]
            empty!(qml.event[1].origin)
            @test_throws ArgumentError preferred_origin(qml.event[1])
        end
        origin1 = QuakeML.Origin(time=Dates.now(), longitude=0, latitude=0,
            public_id="smi:QuakeML.jl/origin/a")
        origin2 = QuakeML.Origin(time=Dates.now(), longitude=1, latitude=1,
            public_id="smi:QuakeML.jl/origin/b")
        origin3 = QuakeML.Origin(time=Dates.now(), longitude=2, latitude=2,
            public_id="smi:QuakeML.jl/origin/c")
        event = QuakeML.Event(public_id="smi:QuakeML.jl/event/a",
            preferred_origin_id="smi:QuakeML.jl/origin/b")
        # No origins
        @test_throws ArgumentError preferred_origin(event)
        # Returns first regardless
        push!(event.origin, origin1) # Add one origin with wrong ID
        @test preferred_origin(event) == origin1
        @test_nowarn preferred_origin(event, verbose=true)
        # Finds the preferred one
        push!(event.origin, origin2) # Add the right origin
        @test preferred_origin(event) == origin2
        @test preferred_origin(event).longitude.value == 1
        # Returns the first if no match
        empty!(event.origin)
        append!(event.origin, [origin1, origin3]) # Two wrong origins
        @test preferred_origin(event) == origin1
        @test_nowarn preferred_origin(event)
        @test_logs((:warn, "no origin with preferred id; returning the first origin"),
            preferred_origin(event, verbose=true))
        # Returns the first if no preferred origin at all
        event = QuakeML.Event(public_id="smi:QuakeML.jl/event/a",
            origin=[origin1, origin2, origin3])
        @test preferred_origin(event) == origin1
    end

    @testset "preferred_origins" begin
        let qml = QuakeML.read(joinpath(@__DIR__, "data", "nepal_mw7.2.qml"))
            @test preferred_origins(qml) == [preferred_origin(qml.event[1])]
            empty!(qml.event[1].origin)
            @test_throws ArgumentError preferred_origins(qml)
        end
    end

    @testset "preferred_magnitude" begin
        mag1 = QuakeML.Magnitude(mag=1, public_id="smi:QuakeML.jl/magnitude/a")
        mag2 = QuakeML.Magnitude(mag=2, public_id="smi:QuakeML.jl/magnitude/b")
        mag3 = QuakeML.Magnitude(mag=3, public_id="smi:QuakeML.jl/magnitude/c")
        event = QuakeML.Event(public_id="smi:QuakeML.jl/event/a",
            preferred_magnitude_id="smi:QuakeML.jl/magnitude/b")
        # No magnitudes
        @test_throws ArgumentError preferred_magnitude(event)
        # Returns first regardless
        push!(event.magnitude, mag1)
        @test preferred_magnitude(event) == mag1
        @test_nowarn preferred_magnitude(event, verbose=true)
        # Finds the preferred one
        push!(event.magnitude, mag2)
        @test preferred_magnitude(event) == mag2
        @test preferred_magnitude(event).mag.value == 2
        # Return the first if no match
        empty!(event.magnitude)
        append!(event.magnitude, [mag1, mag3])
        @test preferred_magnitude(event) == mag1
        @test_nowarn preferred_magnitude(event)
        @test_logs((:warn, "no magnitude with preferred id; returning the first magnitude"),
            preferred_magnitude(event, verbose=true))
        # Returns the first if no preferred magnitude at all
        event = QuakeML.Event(public_id="smi:QuakeML.jl/event/a",
            magnitude=[mag1, mag2, mag3])
        @test preferred_magnitude(event) == mag1

    end

    @testset "preferred_magnitudes" begin
        let qml = QuakeML.read(joinpath(@__DIR__, "data", "nepal_mw7.2.qml"))
            @test preferred_magnitudes(qml) == [preferred_magnitude(qml.event[1])]
            empty!(qml.event[1].magnitude)
            @test_throws ArgumentError preferred_magnitudes(qml)
        end
    end

    @testset "preferred_focal_mechanism" begin
        fm1 = QuakeML.FocalMechanism(public_id="smi:QuakeML.jl/focmech/a")
        fm2 = QuakeML.FocalMechanism(public_id="smi:QuakeML.jl/focmech/b")
        fm3 = QuakeML.FocalMechanism(public_id="smi:QuakeML.jl/focmech/c")
        event = QuakeML.Event(public_id="smi:QuakeML.jl/event/a",
            preferred_focal_mechanism_id="smi:QuakeML.jl/focmech/b")
        # No focal mechanisms
        @test_throws ArgumentError preferred_focal_mechanism(event)
        # Returns first regardless
        push!(event.focal_mechanism, fm1)
        @test preferred_focal_mechanism(event) == fm1
        @test_nowarn preferred_focal_mechanism(event, verbose=true)
        # Finds the preferred one
        push!(event.focal_mechanism, fm2)
        @test preferred_focal_mechanism(event) == fm2
        # Return the first if no match
        empty!(event.focal_mechanism)
        append!(event.focal_mechanism, [fm1, fm3])
        @test preferred_focal_mechanism(event) == fm1
        @test_nowarn preferred_focal_mechanism(event)
        @test_logs((:warn, "no focal mechanism with preferred id; returning the first focal mechanism"),
            preferred_focal_mechanism(event, verbose=true))
        # Returns the first if no preferred focal mechanism at all
        event = QuakeML.Event(public_id="smi:QuakeML.jl/event/a",
            focal_mechanism=[fm1, fm2, fm3])
        @test preferred_focal_mechanism(event) == fm1

    end

    @testset "preferred_focal_mechanisms" begin
        let qml = QuakeML.read(joinpath(@__DIR__, "data", "nepal_mw7.2.qml"))
            @test preferred_focal_mechanisms(qml) == [preferred_focal_mechanism(qml.event[1])]
            empty!(qml.event[1].focal_mechanism)
            @test_throws ArgumentError preferred_focal_mechanisms(qml)
        end
    end

end

max_stars_repo_name: anowacki/QuakeML


text:
using CompScienceMeshes

m1 = meshrectangle(1.0, 0.5, 0.1, 3)
m1 = CompScienceMeshes.translate(m1, [0.0, 0.5, 0.0])
m2 = meshrectangle(1.0, 0.5, 0.1, 3)

nc1 = numcells(m1)
nc2 = numcells(m2)

nv1 = numvertices(m1)
nv2 = numvertices(m2)

M1 = weld(m1,m2)
@assert numcells(M1) == nc1 + nc2
@assert numvertices(M1) == nv1 + nv2 - 11

M2 = weld(m1,-m2)
@assert numcells(M2) == nc1 + nc2
@assert numvertices(M2) == nv1 + nv2 - 11

@assert CompScienceMeshes.isoriented(M1)
@assert !CompScienceMeshes.isoriented(M2)

e1 = skeleton(m1,1)
e2 = skeleton(m2,1)
E1 = skeleton(M1,1)

S11 = CompScienceMeshes.embedding(e1,E1)
S21 = CompScienceMeshes.embedding(e2,E1)

ext: jl


text:
using Test, Dates, Random, ODBC, Tables, WeakRefStrings, DataFrames, DecFP, CSV, SQLite

# You can also specify which database you want to test for as an environment variable:
# ENV["ODBC_TEST_DRIVERS"] = "mysql"
const TEST_DRIVERS = split(get(ENV, "ODBC_TEST_DRIVERS", "mysql, postgresql, mssql"), r"\s*,\s*")
const TEST_MYSQL = "mysql" in TEST_DRIVERS
const TEST_POSTGRESQL = "postgresql" in TEST_DRIVERS
const TEST_MSSQL = "mssql" in TEST_DRIVERS

# To run these tests first run these statements from the command line in the ODBC root directory:
# odbcinst -i -s -h -f test/setup/mysqltest.odbc.ini
# odbcinst -i -s -h -f test/setup/postgresqltest.odbc.ini
# odbcinst -i -s -h -f test/setup/mssqltest.odbc.ini

# You can modify the above files as needed to point to your mysql, postgres, or mssql
# database and drivers, or to specify your username and password

# To setup docker containers to test against you can run:
# docker run --name mysql -p 3306:3306  -e MYSQL_ALLOW_EMPTY_PASSWORD=true -d mysql
# docker run --name postgres -p 5432:5432 -e POSTGRES_PASSWORD="" -d postgres
# docker run --name mssql -e ACCEPT_EULA=Y -e 'MSSQL_SA_PASSWORD=YourStrong!Passw0rd' -p 1433:1433 -d microsoft/mssql-server-linux:2017-GDR

# You will also have needed to install drivers for mysql, postgres, and mssql and add those
# to your local odbcinst.ini file. For example:
# [MySQL]
# Description = MySQL driver
# Driver = /usr/local/lib/libmyodbc5a.so
# FileUsage = 1

# [PostgreSQL ANSI]
# Description = ODBC for PostgreSQL 9.6 ANSI
# Driver = /usr/local/lib/psqlodbca.so
# FileUsage = 1

# [ODBC Driver 13 for SQL Server]
# Description=Microsoft ODBC Driver 13 for SQL Server
# Driver=/usr/local/lib/libmsodbcsql.13.dylib
# UsageCount=1

@testset "ODBC.jl" begin

    @show TEST_DRIVERS
    @show ODBC.drivers()
    @show ODBC.dsns()
    @show ODBC.API.odbc_dm
    @show run(`odbcinst -q -d`)
    run(`uname -a`)

    if TEST_MYSQL
        include("mysql.jl")
    else
        @warn "Skipping mysql tests..."
    end

    if TEST_POSTGRESQL
        include("postgresql.jl")
    else
        @warn "Skipping postgresql tests..."
    end

    if TEST_MSSQL
        include("mssql.jl")
    else
        @warn "Skipping mssql tests..."
    end
end


text:
using Pkg
Pkg.activate(".")
using PackageCompiler
using Libdl
dir = @__DIR__
pwd = @__DIR__

new_project = dir

sysimage_path = joinpath(pwd, "sys." * Libdl.dlext)
script = joinpath(@__DIR__, "gen_ThetisSupport_script.jl")
create_sysimage(:JSON; sysimage_path=sysimage_path,
    project=new_project,
    precompile_execution_file=joinpath(@__DIR__, "gen_ThetisSupport_execution.jl"),
    precompile_statements_file=joinpath(@__DIR__, "gen_ThetisSupport_precompile.jl"),
    script=script)

max_stars_repo_name: AkillesAILimited/Thetis

max_issues_repo_name: AkillesAILimited/Thetis

max_forks_repo_name: AkillesAILimited/Thetis


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct ReactivePowerDroop <: ReactivePowerControl
        kq::Float64
        ωf::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of a Reactive Power droop controller

# Arguments
- `kq::Float64`: frequency droop gain, validation range: (0, nothing)
- `ωf::Float64`: frequency droop gain, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct ReactivePowerDroop <: ReactivePowerControl
    "frequency droop gain"
    kq::Float64
    "frequency droop gain"
    ωf::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function ReactivePowerDroop(kq, ωf, ext=Dict{String, Any}(), )
    ReactivePowerDroop(kq, ωf, ext, [:qm], 1, InfrastructureSystemsInternal(), )
end

function ReactivePowerDroop(; kq, ωf, ext=Dict{String, Any}(), )
    ReactivePowerDroop(kq, ωf, ext, )
end

# Constructor for demo purposes; non-functional.
function ReactivePowerDroop(::Nothing)
    ReactivePowerDroop(;
        kq=0,
        ωf=0,
        ext=Dict{String, Any}(),
    )
end

"""Get ReactivePowerDroop kq."""
get_kq(value::ReactivePowerDroop) = value.kq
"""Get ReactivePowerDroop ωf."""
get_ωf(value::ReactivePowerDroop) = value.ωf
"""Get ReactivePowerDroop ext."""
get_ext(value::ReactivePowerDroop) = value.ext
"""Get ReactivePowerDroop states."""
get_states(value::ReactivePowerDroop) = value.states
"""Get ReactivePowerDroop n_states."""
get_n_states(value::ReactivePowerDroop) = value.n_states
"""Get ReactivePowerDroop internal."""
get_internal(value::ReactivePowerDroop) = value.internal


text:
module CairoMakie

using AbstractPlotting, LinearAlgebra
using Colors, GeometryBasics, FileIO, StaticArrays
import Cairo

using AbstractPlotting: Scene, Lines, Text, Image, Heatmap, Scatter, @key_str, broadcast_foreach
using AbstractPlotting: convert_attribute, @extractvalue, LineSegments, to_ndim, NativeFont
using AbstractPlotting: @info, @get_attribute, Combined
using AbstractPlotting: to_value, to_colormap, extrema_nan

const LIB_CAIRO = if isdefined(Cairo, :libcairo)
    Cairo.libcairo
else
    Cairo._jl_libcairo
end

const OneOrVec{T} = Union{
    T,
    Vec{N1, T} where N1,
    NTuple{N2, T} where N2,
}

# re-export AbstractPlotting
for name in names(AbstractPlotting)
    @eval import AbstractPlotting: $(name)
    @eval export $(name)
end

include("infrastructure.jl")
include("utils.jl")
include("fonts.jl")
include("primitives.jl")
include("overrides.jl")

function __init__()
    activate!()
    AbstractPlotting.register_backend!(AbstractPlotting.current_backend[])
end

function display_path(type::String)
    if !(type in ("svg", "png", "pdf", "eps"))
        error("Only \"svg\", \"png\", \"eps\" and \"pdf\" are allowed for `type`. Found: $(type)")
    end
    return joinpath(@__DIR__, "display." * type)
end

function activate!(; inline = true, type = "svg")
    AbstractPlotting.current_backend[] = CairoBackend(display_path(type))
    AbstractPlotting.use_display[] = !inline
    return
end

end


text:
@testset "LightGraphs functions" begin

# This rather messy file is more of a TODO list than a proper testset at the moment.

#=
@testset "a_star" begin
    g = testgraphA()
    @info a_star(g, 1, 4)
    @info a_star(g, 1, 4, ones(nv(g), nv(g)) )
end

@testset "add_vertex!" begin
    g = testgraphA()
    add_vertex!(g)
    Base.print_matrix(stdout, weights(g))
end
@testset "add_vertices!" begin
    g = testgraphA()
    @test add_vertices!(g, 2) == 2
    Base.print_matrix(stdout, weights(g))
end
@testset "adjacency_matrix" begin
    g = testgraphA()
    @info adjacency_matrix(g)
end
@testset "adjacency_spectrum" begin
    g = testgraphA()
    @info adjacency_spectrum(g)
end
@testset "articulation" begin
    g = SimpleValueGraph(PathGraph(5))
    add_edge!(g, 2, 3, 10.0)
    add_edge!(g, 3, 4, 0.0)
    @test sort(articulation(g)) == [2,3,4]
end
@testset "attracting_components" begin
    # TODO only defined for directed graphs
end

@testset "barabasi_albert!" begin
    g = testgraphA()
    nv_old = nv(g)
    ne_old = ne(g)
    n = 10
    k = 3
    barabasi_albert!(g, n, k)
    @test nv(g) == n
    @test ne(g) == (n - nv_old) * k + ne_old
end
@testset "bellman_ford_shortest_paths" begin
    g = testgraphA()
    @test bellman_ford_shortest_paths(g, 1).dists == bellman_ford_shortest_paths(g, 1, weights(g)).dists
    @info bellman_ford_shortest_paths(g, [2, 4])
end
@testset "bfs_parents" begin
    g = testgraphA()
    @test length(bfs_parents(g, 1)) == nv(g)
    @test length(bfs_parents(g, 1, dir=:in)) == nv(g)
end
@testset "bfs_tree" begin
    g = testgraphA()
    @test nv(bfs_tree(g, 1)) == nv(g)
    @test ne(bfs_tree(g, 1, dir=:in)) == nv(g) - 1
end
@testset "bipartite_map" begin
    g = SimpleValueGraph(CompleteBipartiteGraph(3,5))
    @info bipartite_map(g)
    @test length(bipartite_map(SimpleValueGraph(CycleGraph(5)))) == 0
end
@testset "blockcounts" begin
    # TODO not sure how that works
end
@testset "blockfractions" begin
    # TODO not sure how that works
end
@testset "cartesian_product" begin
    g = testgraphA()
    g2 = cartesian_product(g, g)
    @test nv(g2) == nv(g)^2
    Base.print_matrix(stdout, weights(g2))
    # TODO cartesian product does not yet consider weights
end
=#

@testset "common_neighbors" begin
    g = testgraphA()
    @test common_neighbors(g, 1, 2) == [4,5,6,7]
end

#=
@testset "connected_components" begin
    g = testgraphA()
    @test length(connected_components(blockdiag(g,g))) == 2
end

@testset "core_number" begin
    g = SimpleValueGraph(StarGraph(7))
    @test unique(core_number(g)) == [1]
end

@testset "crosspath" begin
    # TODO this does not work as intended
end

@testset "cycle_basis" begin
    # TODO only works for AbstractSimpleGraphs
end

@testset "degree" begin
    @info "degree"
    g = SimpleValueGraph(CompleteBipartiteGraph(2,3))
    @test degree(g) == [3, 3, 2, 2, 2]
    @test degree(g, 1) == 3
end

@testset "degree_centrality" begin
    @info "degree_centrality"
    g = erdos_renyi(10, 20)
    h = SimpleValueGraph(g)
    @test degree_centrality(g) == degree_centrality(h)
end

@testset "dfs_tree" begin
    @info "dfs_tree"
    g = erdos_renyi(10, 20)
    h = graph_with_randvals(g)
    @test dfs_tree(g, 2) == dfs_tree(h, 2)
end

@testset "diameter" begin
    @info "diameter"
    #= seems to crash for some reasons
    g = erdos_renyi(10, 30)
    h = graph_with_randvals(g)
    @test diameter(g) == diameter(h)
    =#
end

@testset "difference" begin
    # TODO this operators does not work correctly yet
end

@testset "diffusion" begin
    @info "diffusion"
    # TODO this seems to work, but diffusion should maybe consider the weights
    g = testgraphA()
    @info diffusion(g, 0.5, 4)
    @info diffusion_rate(g, 0.5, 4)
end

@testset "dijkstra_shortest_paths" begin
    @info "diijkstar_shortest_paths"
    # TODO test with different weights
    # TODO crashes with negative weights or zero weights
    g = SimpleValueGraph(5)
    add_edge!(g, 1, 2, 2.5)
    add_edge!(g, 1, 3, 0.1)
    add_edge!(g, 1, 4, 1.5)
    add_edge!(g, 1, 5)
    add_edge!(g, 2, 5)

    @test dijkstra_shortest_paths(g, 1).dists[2:end] == [2.0, 0.1, 1.5, 1.0]
end

@testset "eccentricity" begin
    # TODO not working, relies on dijkstra
end

@testset "edit_distance" begin
    # TODO should edit_distance use the weights?
end

@testset "egonet" begin
    # TODO not working, could be because it relies on induced subgraph
end

@testset "eigenvector_centrality" begin
    g = testgraphA()
    @info eigenvector_centrality(g)
end

@testset "floyd_warshall_shortest_paths" begin
    # TODO verify that this really works with negative values
    g = graph_with_randvals(erdos_renyi(10, 30))
    @info floyd_warshall_shortest_paths(g)
end

@testset "gdistances" begin
    g = erdos_renyi(10, 20)
    h = SimpleValueGraph(g)
    @test gdistances(h, 2) == gdistances(g, 2)
    dists = fill(typemax(eltype(h)), nv(h))
    gdistances!(h, 3, dists)
    @test dists == gdistances(g, 3)
end

@testset "global_clustering_coefficient" begin
    g = erdos_renyi(10, 20)
    h = SimpleValueGraph(g)
    @test global_clustering_coefficient(h) == global_clustering_coefficient(g)
end

@testset "greedy_color" begin
    g = graph_with_randvals(CompleteMultipartiteGraph([2,3,4,5]))
    @test greedy_color(g, sort_degree=false, reps=5).num_colors >= 4
    @test greedy_color(g, sort_degree=true, reps=5).num_colors >= 4
end

@testset "has_negative_edge_cycle" begin
    # TODO this functions ignores edge weights at the moment
end

@testset "has_path" begin
    g = erdos_renyi(20, 0.1)
    h = graph_with_randvals(g)
    @test all(v -> has_path(g, 2, v, exclude_vertices=[3,4]) == has_path(h, 2, v, exclude_vertices=[3,4]), vertices(g))
end

@testset "has_self_loops" begin
    g = graph_with_randvals(testgraphA())
    @test !has_self_loops(g)
    add_edge!(g, 5, 5)
    @test has_self_loops(g)
end

@testset "label_propagation" begin
    # TODO function does not consider weights at the moment
    g = CompleteGraph(4)
    h = graph_with_randvals(blockdiag(g, g))
    add_edge!(h, 1, nv(g) + 1)
    @test length(unique(label_propagation(h)[1])) == 2
end

@testset "laplacian_matrix" begin
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @test laplacian_matrix(g) == laplacian_matrix(h)
    # TODO does not consider weights at the moment
end

@testset "laplacian_spectrum" begin
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @test laplacian_spectrum(g) ≈ laplacian_spectrum(h)
    # TODO does not consider weights at the moment
end

=#

@testset "laplacian_spectrum" begin
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @test local_clustering(g, 2) == local_clustering(h, 2)
end

@testset "laplacian_spectrum" begin
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @test local_clustering_coefficient(g, 2) == local_clustering_coefficient(h, 2)
end

@testset "maximal_cliques" begin
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @test maximal_cliques(g) == maximal_cliques(h)
end

#=

@testset "maximum_adjacency_visit" begin
    # TODO ducmentation is wrong, default weights are only 1.0 for simple graphs
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @test maximum_adjacency_visit(g, weights(h)) == maximum_adjacency_visit(h)
    @test maximum_adjacency_visit(g, weights(h)) == maximum_adjacency_visit(h, weights(h))
end

@testset "merge_vertices" begin
    # TODO merge_vertices and merge_vertices! do not work correctly yet
end

@testset "mincut" begin
    # TODO check what happens for negative weights
    g = SimpleValueGraph(3)
    add_edge!(g, 1, 2, 0.1)
    add_edge!(g, 1, 3, 1.0)
    add_edge!(g, 2, 3, 0.1)
    @test mincut(g)[2] ≈ 0.2
end

@testset "neighborhood" begin
    # TODO what about negative values?
    g = SimpleValueGraph(4)
    add_edge!(g, 1, 2, 0.5)
    add_edge!(g, 2, 3, 1.5)
    add_edge!(g, 3, 4, 1.5)
    @test sort(neighborhood(g, 2, 0.2)) == [2]
    @test sort(neighborhood(g, 2, 1.2)) == [1, 2]
    @test sort(neighborhood(g, 2, 2.2)) == [1, 2, 3]

    dists = sort(neighborhood_dists(g, 2, 10.0))
    compare_to = [(1, 0.5), (2, 0.0), (3, 1.5), (4, 3.0)]
    @test all(i -> dists[i][2] ≈ compare_to[i][2], 1:length(compare_to))
end

@testset "non_backtracking_matrix" begin
    # TODO this functions and non_backtracking_randomwalk do not work yet
end
=#

@testset "num_self_loops" begin
    g = graph_with_randvals(testgraphA())
    @test num_self_loops(g) == 0
    add_edge!(g, 3, 3)
    @test num_self_loops(g) == 1
    add_edge!(g, 4, 4)
    @test num_self_loops(g) == 2
end

#=
@testset "outdegree" begin
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @test outdegree(g) == outdegree(h)
    @test all(v -> outdegree(g, v) == outdegree(h, v), vertices(g))
end

@testset "outdegree_centrality" begin
    # TODO outdegree_centrality is not documented
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @test outdegree_centrality(g) == outdegree_centrality(h)
end

@testset "pagerank" begin
    # TODO pagerank seems not to consider weights
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @test pagerank(g) ≈ pagerank(h)
end

@testset "period" begin
    # TODO works only with directed graphs
end

@testset "periphery" begin
    # TODO not working, relies on dijkstra
end

@testset "periphery" begin
    # TODO not working, relies on dijkstra
end

@testset "prim_mst" begin
    # TODO seems, to work but returns SimpleGraph Edges, verify that this is the right approach
end

@testset "radiality_centrality" begin
    # TODO not working, relies on dijkstra
end

@testset "radius" begin
    # TODO not working, relies on dijkstra
end

@testset "randomwalk" begin
    # TODO maybe should consider edge_weights
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @info randomwalk(h, 2, 20)
end

@testset "savegraph" begin
    # TODO not yet working
end
=#

@testset "saw" begin
    g = erdos_renyi(20, 30)
    h = graph_with_randvals(g)
    @info saw(h, 2, 100)
end

#=
@testset "simplecycles_hawick_james" begin
    # TODO seems only to work with SimpleGraphs
end

@testset "spectral_distance" begin
    # TODO should to involve weights?
    g1 = erdos_renyi(20, 30)
    g2 = erdos_renyi(20, 35)
    h1 = graph_with_randvals(g1)
    h2 = graph_with_randvals(g2)
    @test spectral_distance(g1, g2) == spectral_distance(h1, h2)
end

@testset "squash" begin
    # TODO does not work yet, maybe implement using convert
end

@testset "stress_centrality" begin
    # TODO does use weights?
    # TODo does not work, relies on dijkstra
end

@testset "strongly_connected_components" begin
    # TODO only works for directed graphs
end

@testset "symmetric difference" begin
    # TODO not correctly implemented yet
end
=#

@testset "tensor_product" begin
    # TODO does not work correctly yet
end

@testset "topological_sort_by_dfs" begin
    # TODO only works for directed graphs
end

#=
@testset "randomwalk" begin
    # TODO maybe should consider edge_weights
    g = erdos_renyi(10, 40)
    h = graph_with_randvals(g)
    @test triangles(g) == triangles(h)
end
=#

@testset "tensor_product" begin
    # TODO does not work correctly yet
end

@testset "union" begin
    # TODO does not work correctly yet
end

@testset "vertex_cover" begin
    g = erdos_renyi(10, 40)
    h = graph_with_randvals(g)
    @test vertex_cover(g, DegreeVertexCover()) == vertex_cover(h, DegreeVertexCover())
end

@testset "yen_k_shortest_paths" begin
    # TODO seems to freeze for negtive cycles
end

end


text:
#=
common code for JuMP examples
=#

import JuMP
const MOI = JuMP.MOI
const MOIU = MOI.Utilities

MOIU.@model(SOCExpPSD,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan,),
    (MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
    MOI.SecondOrderCone, MOI.RotatedSecondOrderCone,
    MOI.ExponentialCone, MOI.PositiveSemidefiniteConeTriangle,),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction,),
    true,
    )

MOIU.@model(ExpPSD,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan,),
    (MOI.Reals, MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives,
    MOI.ExponentialCone, MOI.PositiveSemidefiniteConeTriangle,),
    (),
    (),
    (MOI.ScalarAffineFunction,),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction,),
    true,
    )

abstract type ExampleInstanceJuMP{T <: Real} <: ExampleInstance{T} end

# fallback: just check optimal status
function test_extra(inst::ExampleInstanceJuMP, model::JuMP.Model)
    @test JuMP.termination_status(model) == MOI.OPTIMAL
end

function run_instance(
    ex_type::Type{<:ExampleInstanceJuMP},
    inst_data::Tuple,
    extender::Union{Symbol, Nothing} = nothing,
    inst_options::NamedTuple = NamedTuple(),
    solver_type = Hypatia.Optimizer;
    default_options::NamedTuple = NamedTuple(),
    test::Bool = true,
    verbose::Bool = true,
    )
    new_options = merge(default_options, inst_options)

    verbose && println("setup model")
    setup_time = @elapsed (model, model_stats) =
        setup_model(ex_type, inst_data, extender, new_options, solver_type)

    verbose && println("solve and check")
    check_time = @elapsed solve_stats = solve_check(model, test = test)

    return (; model_stats..., solve_stats..., setup_time,
        check_time, :script_status => "Success")
end

function setup_model(
    ex_type::Type{<:ExampleInstanceJuMP},
    inst_data::Tuple,
    extender::Union{Symbol, Nothing},
    solver_options::NamedTuple,
    solver_type;
    rseed::Int = 1,
    )
    # setup example instance and JuMP model
    Random.seed!(rseed)
    inst = ex_type(inst_data...)
    model = build(inst)

    is_hypatia_opt = (solver_type == Hypatia.Optimizer)
    opt = hyp_opt = (is_hypatia_opt ? Hypatia.Optimizer(; solver_options...) :
        Hypatia.Optimizer(use_dense_model = false))
    if !isnothing(extender)
        # use MOI automated extended formulation
        extT = @eval $extender{Float64}()
        opt = MOI.Bridges.full_bridge_optimizer(
            MOIU.CachingOptimizer(extT, opt), Float64)
        # for PolyJuMP/SumOfSquares models
        for B in model.bridge_types
            MOI.Bridges.add_bridge(opt, B{Float64})
        end
    end
    backend = JuMP.backend(model)
    MOIU.reset_optimizer(backend, opt)
    MOIU.attach_optimizer(backend)
    isnothing(extender) || MOIU.attach_optimizer(backend.optimizer.model)
    flush(stdout); flush(stderr)

    hyp_model = hyp_opt.model
    if is_hypatia_opt
        model.ext[:inst] = inst
    else
        # not using Hypatia to solve, so setup new JuMP model with Hypatia data
        (A, b, c, G, h) =
            (hyp_model.A, hyp_model.b, hyp_model.c, hyp_model.G, hyp_model.h)
        (cones, cone_idxs) = (hyp_model.cones, hyp_model.cone_idxs)
        new_model = JuMP.Model()
        new_model.ext[:hyp_data] = hyp_model
        JuMP.@variable(new_model, x[1:length(c)])
        JuMP.@objective(new_model, Min, dot(c, x))
        eq_refs = JuMP.@constraint(new_model, A * x .== b)
        cone_refs = Vector{JuMP.ConstraintRef}(undef, length(cones))
        for (k, cone_k) in enumerate(cones)
            idxs = cone_idxs[k]
            h_k = h[idxs]
            G_k = G[idxs, :]
            moi_set = cone_from_hyp(cone_k)
            if Hypatia.needs_untransform(moi_set)
                Hypatia.untransform_affine(moi_set, h_k)
                @inbounds @views for j in 1:size(G_k, 2)
                    Hypatia.untransform_affine(moi_set, G_k[:, j])
                end
            end
            cone_refs[k] = JuMP.@constraint(new_model, h_k - G_k * x in moi_set)
        end
        new_model.ext[:x_var] = x
        new_model.ext[:eq_refs] = eq_refs
        new_model.ext[:cone_refs] = cone_refs

        opt = solver_type(; solver_options...)
        JuMP.set_optimizer(new_model, () -> opt)
        model = new_model
        flush(stdout); flush(stderr)
    end

    model_stats = get_model_stats(hyp_model)
    return (model, model_stats)
end

function solve_check(
    model::JuMP.Model;
    test::Bool = true,
    )
    JuMP.optimize!(model) # TODO make sure it doesn't copy again
    flush(stdout); flush(stderr)

    opt = JuMP.backend(model).optimizer
    if !isa(opt, Hypatia.Optimizer)
        backend_model = JuMP.backend(model).optimizer.model
        if backend_model isa MOIU.CachingOptimizer
            opt = backend_model.optimizer
        end
    end

    if opt isa Hypatia.Optimizer
        test && test_extra(model.ext[:inst], model)
        flush(stdout); flush(stderr)
        (solve_stats, _) = process_result(opt.model, opt.solver)
        return solve_stats
    end

    test && @info("cannot run example tests if solver is not Hypatia")

    solve_time = JuMP.solve_time(model)
    iters = MOI.get(model, MOI.BarrierIterations())
    primal_obj = JuMP.objective_value(model)
    dual_obj = JuMP.dual_objective_value(model)
    moi_status = MOI.get(model, MOI.TerminationStatus())
    if haskey(moi_hyp_status_map, moi_status)
        hyp_status = moi_hyp_status_map[moi_status]
    else
        @warn("MOI status $moi_status not handled")
        hyp_status = Solvers.UnknownStatus
    end

    hyp_data = model.ext[:hyp_data]
    eq_refs = model.ext[:eq_refs]
    cone_refs = model.ext[:cone_refs]
    x = JuMP.value.(model.ext[:x_var])
    y = (isempty(eq_refs) ? Float64[] : -JuMP.dual.(eq_refs))
    s_cones = Vector{Vector{Float64}}(undef, length(cone_refs))
    z_cones = Vector{Vector{Float64}}(undef, length(cone_refs))
    for (k, cr) in enumerate(cone_refs)
        moi_set = MOI.get(cr.model, MOI.ConstraintSet(), cr)
        idxs = Hypatia.permute_affine(moi_set, 1:length(hyp_data.cone_idxs[k]))
        s_k = Hypatia.rescale_affine(moi_set, JuMP.value.(cr))
        z_k = Hypatia.rescale_affine(moi_set, JuMP.dual.(cr))
        s_cones[k] = s_k[idxs]
        z_cones[k] = z_k[idxs]
    end
    s = vcat(s_cones...)
    z = vcat(z_cones...)

    rel_obj_diff = (primal_obj - dual_obj) / (1 + abs(dual_obj))
    compl = dot(s, z)
    (x_viol, y_viol, z_viol) =
        certificate_violations(hyp_status, hyp_data, x, y, z, s)
    flush(stdout); flush(stderr)

    solve_stats = (;
        :status => hyp_status, solve_time, iters, primal_obj, dual_obj,
        rel_obj_diff, compl, x_viol, y_viol, z_viol,
        )
    return solve_stats
end

# get Hypatia status from MOI status
moi_hyp_status_map = Dict(
    MOI.OPTIMAL => Solvers.Optimal,
    MOI.INFEASIBLE => Solvers.PrimalInfeasible,
    MOI.DUAL_INFEASIBLE => Solvers.DualInfeasible,
    MOI.SLOW_PROGRESS => Solvers.SlowProgress,
    MOI.ITERATION_LIMIT => Solvers.IterationLimit,
    MOI.TIME_LIMIT => Solvers.TimeLimit,
    MOI.OTHER_ERROR => Solvers.UnknownStatus,
    )

# get MOI cones (defined in MOI itself, not Hypatia) from some Hypatia cones
cone_from_hyp(cone::Cones.Cone) =
    error("cannot transform Hypatia cone of type $(typeof(cone)) to MOI cone")

cone_from_hyp(cone::Cones.Nonnegative) = MOI.Nonnegatives(Cones.dimension(cone))

cone_from_hyp(cone::Cones.EpiNormInf{T, T}) where {T <: Real} =
    (Cones.use_dual_barrier(cone) ? MOI.NormOneCone :
    MOI.NormInfinityCone)(Cones.dimension(cone))

cone_from_hyp(cone::Cones.EpiNormEucl) =
    MOI.SecondOrderCone(Cones.dimension(cone))

cone_from_hyp(cone::Cones.EpiPerSquare) =
    MOI.RotatedSecondOrderCone(Cones.dimension(cone))

cone_from_hyp(cone::Cones.HypoPerLog) =
    (@assert Cones.dimension(cone) == 3; MOI.ExponentialCone())

cone_from_hyp(cone::Cones.EpiRelEntropy) =
    MOI.RelativeEntropyCone(Cones.dimension(cone))

cone_from_hyp(cone::Cones.HypoGeoMean) =
    MOI.GeometricMeanCone(Cones.dimension(cone))

cone_from_hyp(cone::Cones.GeneralizedPower) =
    (@assert Cones.dimension(cone) == 3; MOI.PowerCone{Float64}(cone.α[1]))

cone_from_hyp(cone::Cones.EpiNormSpectral{T, T}) where {T <: Real} =
    (Cones.use_dual_barrier(cone) ? MOI.NormNuclearCone :
    MOI.NormSpectralCone)(cone.n, cone.m)

cone_from_hyp(cone::Cones.PosSemidefTri{T, T}) where {T <: Real} =
    MOI.PositiveSemidefiniteConeTriangle(cone.side)

cone_from_hyp(cone::Cones.HypoPerLogdetTri{T, T}) where {T <: Real} =
    MOI.LogDetConeTriangle(cone.side)

cone_from_hyp(cone::Cones.HypoRootdetTri{T, T}) where {T <: Real} =
    MOI.RootDetConeTriangle(cone.side)


text:
function poordump(@nospecialize(o), limit::Bool=false) # richdump minus annotations
    dom = Vector{Firelight.Node}()
    io = IOBuffer()
    Firelight.richdump(IOContext(io, :limit => limit), dom, o)
    @test !isempty(dom)
    return String(take!(io)), dom
end

@test poordump("Hello\nWorld")[1] == "String \"Hello\\nWorld\""
@test poordump(Core.svec())[1] == "empty SimpleVector"
@test poordump(Core.svec(1, "Hello\nWorld", []))[1] ==
    "SimpleVector\n    1: 1\n    2: \"Hello\\nWorld\"\n    3: Any[]"
@test poordump(:hello)[1] == "Symbol hello"
@test poordump(Firelight)[1] == "Module Firelight"
@test poordump(Main)[1] == "Module Main"
@test poordump(Core.Compiler)[1] == "Module Core.Compiler"
@test poordump(C_NULL)[1] == "Ptr{Nothing} @0x0000000000000000"
@test poordump(1)[1] == "Int64    1"
@test poordump(0x1)[1] == "UInt8    0x01"
@test poordump(3 + 4im)[1] == "Complex{Int64}\n    re: 3\n    im: 4"
@test poordump(fill(0x1, 100))[1] ==
    "Array{UInt8}((100,)) UInt8[0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01]"
@test poordump(fill(0x1, 100), true)[1] ==
    "Array{UInt8}((100,)) UInt8[0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01  …  0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01, 0x01]"
@test poordump(fill(0x1, 5), true)[1] ==
    "Array{UInt8}((5,)) UInt8[0x01, 0x01, 0x01, 0x01, 0x01]"
@test poordump(fill("1", 100))[1] ==
    "Array{String}((100,))
    1: \"1\"\n    2: \"1\"\n    3: \"1\"\n    4: \"1\"\n    5: \"1\"\n    6: \"1\"\n    7: \"1\"\n    8: \"1\"\n    9: \"1\"\n    10: \"1\"
    11: \"1\"\n    12: \"1\"\n    13: \"1\"\n    14: \"1\"\n    15: \"1\"\n    16: \"1\"\n    17: \"1\"\n    18: \"1\"\n    19: \"1\"\n    20: \"1\"
    21: \"1\"\n    22: \"1\"\n    23: \"1\"\n    24: \"1\"\n    25: \"1\"\n    26: \"1\"\n    27: \"1\"\n    28: \"1\"\n    29: \"1\"\n    30: \"1\"
    31: \"1\"\n    32: \"1\"\n    33: \"1\"\n    34: \"1\"\n    35: \"1\"\n    36: \"1\"\n    37: \"1\"\n    38: \"1\"\n    39: \"1\"\n    40: \"1\"
    41: \"1\"\n    42: \"1\"\n    43: \"1\"\n    44: \"1\"\n    45: \"1\"\n    46: \"1\"\n    47: \"1\"\n    48: \"1\"\n    49: \"1\"\n    50: \"1\"
    51: \"1\"\n    52: \"1\"\n    53: \"1\"\n    54: \"1\"\n    55: \"1\"\n    56: \"1\"\n    57: \"1\"\n    58: \"1\"\n    59: \"1\"\n    60: \"1\"
    61: \"1\"\n    62: \"1\"\n    63: \"1\"\n    64: \"1\"\n    65: \"1\"\n    66: \"1\"\n    67: \"1\"\n    68: \"1\"\n    69: \"1\"\n    70: \"1\"
    71: \"1\"\n    72: \"1\"\n    73: \"1\"\n    74: \"1\"\n    75: \"1\"\n    76: \"1\"\n    77: \"1\"\n    78: \"1\"\n    79: \"1\"\n    80: \"1\"
    81: \"1\"\n    82: \"1\"\n    83: \"1\"\n    84: \"1\"\n    85: \"1\"\n    86: \"1\"\n    87: \"1\"\n    88: \"1\"\n    89: \"1\"\n    90: \"1\"
    91: \"1\"\n    92: \"1\"\n    93: \"1\"\n    94: \"1\"\n    95: \"1\"\n    96: \"1\"\n    97: \"1\"\n    98: \"1\"\n    99: \"1\"\n    100: \"1\""
@test poordump(fill("1", 100), true)[1] ==
    "Array{String}((100,))\n    1: \"1\"\n    2: \"1\"\n    3: \"1\"\n    4: \"1\"\n    5: \"1\"\n  ...\n    96: \"1\"\n    97: \"1\"\n    98: \"1\"\n    99: \"1\"\n    100: \"1\""
@test poordump(fill("1", 5), true)[1] ==
    "Array{String}((5,))\n    1: \"1\"\n    2: \"1\"\n    3: \"1\"\n    4: \"1\"\n    5: \"1\""
@test poordump(Union{Float32, Float64})[1] ==  "Union{Float32, Float64}"
@test poordump(Float32)[1] == "Float32 <: AbstractFloat"
@test poordump(Complex)[1] == "UnionAll where var: {\n    T<:Real\n}\nbody: Complex{T} <: Number\n    re::T\n    im::T"
@test poordump(ComplexF32)[1] == "Complex{Float32} <: Number\n    re::Float32\n    im::Float32"
struct LinearAlgebra_Diagonal{T, V<:AbstractVector{T}} <: AbstractMatrix{T}
    diag::V
end
let LinearAlgebra_Diagonal_name = "LinearAlgebra_Diagonal"
    @test poordump(LinearAlgebra_Diagonal)[1] ==
        "UnionAll where var: {\n    T,\n    V<:AbstractArray{T,1}\n}\nbody: $LinearAlgebra_Diagonal_name{T,V} <: AbstractArray{T,2}\n    diag::V"
    @test poordump(LinearAlgebra_Diagonal.body.body)[1] ==
        "$LinearAlgebra_Diagonal{T,V<:AbstractArray{T,1}} <: AbstractArray{T,2}\n    diag::V"
end
@test poordump(Union{})[1] == "Core.TypeofBottom    Union{}"
@test poordump(typeof(Union{}))[1] == "Core.TypeofBottom <: Type{Union{}}"
@test poordump(Type{Union{}})[1] == "Type{Union{}} <: Any"

max_stars_repo_name: vtjnash/Firelight


text:
export getConvDiffFEMConstraintsAMPL
"""
function getConvDiffFEMConstraintsAMPL

Builds matrices to model constraints to be used in AMPL.
To this end, we take the FEM discretization on all nodes of
the mesh and incorporate the Dirichlet boundary condition.
This yields the constraints

Stiffness * u = Mass*w     (on non-Dirichlet nodes)
            u = gd         (on Dirichled nodes)

We represent these constraints as

A * u = B*w + b, where

A = [Stiffness; I[dir,:]], B = [Mass; zeros], b = [zeros; gd(dir)]

!! Note: The linear system arising from this approach will not be symmetric.
!!       It can be made so by eliminating the Dirichlet nodes. This is done in
!!       getConvDiffFEMParam, which is preferable in terms of efficiency

Required Input:
	M         - describes computational mesh

Keyword Arguments:
	sig       - scalar, >0
	v         - velocity (assumed to be spatially homogeneous)
	gd        - function for Dirichlet boundary conditions
	bc        - Tuple, defining the boundary conditions
	P         - receiver matrix
	Fields    - stores fields, i.e., u
	Ainv      - factorization of PDE

Outputs:
    A = [Stiffness; I[dir,:]] - FEM discretization of constraints on non-Dirichlet nodes
    B = [Mass; zeros]         - mass matrix to model source contribution to non-Dirichlet nodes
    b = [0; gd(dir)]          - vector to encode Dirichlet conditions.

"""
function getConvDiffFEMConstraintsAMPL(M::RegularMesh; v=Matrix{Int}(I,M.dim,1), gd::Function=X->zeros(size(X,1)),
	                                           sig::Number=0.01,bc=(:dir,:neu,:neu,:neu,:neu,:neu))

# get mass and stiffness matrices for all nodes
if M.dim==2
	Mass, Massc, Afull = getFEMMatrices2D(M,sig=sig,v=v)
else
	Mass, Massc, Afull = getFEMMatrices3D(M,sig=sig,v=v)
end

# get indices of nodes on Dirichlet boundary and indices of interior nodes
iddir = ConvDiffMIPDECO.getBoundaryIndicesFEM(M, bc)  # Dirichlet nodes
idint = setdiff(1:size(Afull,1),iddir) # interior nodes

# take constraints from interior nodes and add rows for Dirichlet condition.
# note that rows correspond to test functions, which are not needed on Dirichlet
# boundary
Id = sparse(I, size(Afull,2), size(Afull,2))
#A  = [Afull[idint,:];  Id[iddir,:]]
#B  = [Massc[idint,:]; spzeros(length(iddir),size(Massc,2))]
A  = [Id[iddir,:];Afull[idint,:]]
B  = [spzeros(length(iddir),size(Massc,2)); Massc[idint,:]]
xn = getNodalGrid(M)
#b  = [zeros(length(idint)); gd(xn[iddir,:])]
b  = [ gd(xn[iddir,:]);zeros(length(idint))]
return A,B,b
end

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
module Params
using Random

# Agent and object numbers, length of time
const bigm = 200 # Number of stocks
const bign = 500 # Number of investors
const bigt = 300 # Number of time periods
const bigk = 200 # Number of funds

# Market parameters
const marketstartval = 100 # Market index starting value
const drift = 0.000307888 # Daily market index drift, 8%pa
const marketvol = 1 # Market volatility

# Stock parameters
const betamean = 1 # Average stock beta
const betastd = 0.3 # Stock beta dispersion
const stockstartval = 100 # Price of stocks at the beginning
const stockvolrange = 0.001:0.001:0.01 # Range of stock volatilitiess

# Investor parameters
const reviewprobability = 1/63 # Daily likelihood of performance review
const perfwindow = 1:100 # Performance window for investor
const invcaprange = (50,150) # Investors' range of initial capital
const thresholdmean = 0 # Average investor return threshold for her fund
const thresholdstd = 0.05 # Standard deviation of investor return thresholds

# Fund parameters
const portfsizerange = 10:100 # Range of number of stocks in funds' portfolio
const impactrange = 0.0001:0.0001:0.001 # Stock price impact per currency unit

const plotpath = "./plots/"

# Stylised fact-check parameters
# Number of lags?
# Size of confidence interval?

end #module


text:
module Writers
export Writer, getaccumulator

using StructEquality
using TypeClasses
using DataTypesBasic

"""
like Pair, however assumes that `combine` is defined for the accumulator `acc`

Note that `neutral` may indeed be undefined
"""
struct Writer{Acc, Value}
  acc::Acc
  value::Value
  function Writer(acc, value=nothing)
    new{typeof(acc), typeof(value)}(acc, value)
  end
end

@def_structequal Writer

"""
    getaccumulator(writer::Writer)

Returns the accumulator of the Writer value.

Examples
--------

```jldoctest
julia> using TypeClasses

julia> getaccumulator(Writer("example-accumulator"))
"example-accumulator"
```
"""
getaccumulator(writer::Writer) = writer.acc
Base.get(writer::Writer) = writer.value
Base.getindex(writer::Writer) = writer.value

Base.eltype(::Type{Writer{Acc, T}}) where {Acc, T} = T
Base.eltype(::Type{<:Writer}) = Any

end  # module


text:
printstyled("  Nodal types\n", color=:light_green)
printstyled("    method extensions\n", color=:light_green)
fs = 100.0
j = 7
j1 = 3
j2 = 5
j3 = 6
j4 = 8
nc = 10
nx = 2^12
v1 = 3.0
v2 = 4.0
J = j1:j2

printstyled("      append!\n", color=:light_green)
U = read_nodal("silixa", fstr)
S = deepcopy(U)
S2 = S[j1:j2]
n = S.n
append!(S, S2)
@test S.n == n + S2.n
for f in SeisIO.Nodal.nodalfields
  F1 = getfield(S, f)
  F2 = getfield(S2, f)
  for (i,j) in enumerate(n+1:n+S2.n)
    @test F1[j] == F2[i]
  end
end
@test S.data[:, n+1:n+S2.n] == S2.data

printstyled("      convert\n", color=:light_green)
printstyled("        NodalChannel ⟺ SeisChannel\n", color=:light_green)
C = randSeisChannel()
C1 = convert(NodalChannel, C)
C2 = convert(SeisChannel, C1)
@test C == C2

printstyled("        NodalChannel ⟺ EventChannel\n", color=:light_green)
C1 = EventChannel(C)
C2 = convert(NodalChannel, C1)
C3 = convert(SeisChannel, C2)
@test typeof(C1) == EventChannel
@test typeof(C2) == NodalChannel
@test C == C3

printstyled("        NodalData ⟺ SeisData\n", color=:light_green)
S = randSeisData(nc, nx=nx, s=1.0)
t = mk_t(nx, S.t[1][1,2])
for i in 1:S.n
  S.fs[i] = fs
  S.t[i] = copy(t)
end
D = convert(NodalData, S)
for f in SeisIO.Nodal.nodalfields
  @test length(getfield(D, f)) == nc
end
@test size(D.data) == (nx, nc)
@test NodalData(S) == D

S = convert(SeisData, U)
for i in 1:S.n
  @test isapprox(S.x[i], U.x[i])
end
@test SeisData(U) == S

printstyled("        NodalData ⟺ EventTraceData\n", color=:light_green)
Ev = convert(EventTraceData, S)
D = convert(NodalData, Ev)
@test typeof(Ev) == EventTraceData
for f in SeisIO.Nodal.nodalfields
  @test length(getfield(D, f)) == S.n
end
@test size(D.data) == size(U.data)
@test NodalData(Ev) == D
@test EventTraceData(D) == Ev

printstyled("      deleteat!\n", color=:light_green)
S = deepcopy(U)
deleteat!(S, 3)
@test S.data[:, 1:2] == U.data[:, 1:2]
@test S.data[:, 3:S.n] == U.data[:, 4:U.n]

S = deepcopy(U)
deleteat!(S, J)
@test S.data[:, 1:2] == U.data[:, 1:2]
@test S.data[:, 3:S.n] == U.data[:, 6:U.n]

printstyled("      getindex\n", color=:light_green)
S = deepcopy(U)

C = S[j]
@test S.data[:,j] == C.x
for f in SeisIO.Nodal.nodalfields
  F = getfield(S, f)
  @test getindex(F, j) == getfield(C, f)
end

S2 = S[j1:j2]
for f in SeisIO.Nodal.nodalfields
  F1 = getfield(S, f)
  F2 = getfield(S2, f)
  for (i,j) in enumerate(j1:j2)
    @test F1[j] == F2[i]
  end
end
@test S.data[:, j1:j2] == S2.data

printstyled("      isempty\n", color=:light_green)
S = NodalData()
@test isempty(S) == true
C = NodalChannel()
@test isempty(C) == true
L = NodalLoc()
@test isempty(L) == true

printstyled("      isequal\n", color=:light_green)
S = deepcopy(U)
@test S == U
S = NodalData()
T = NodalData()
@test S == T

C = NodalChannel()
D = NodalChannel()
@test C == D

L1 = NodalLoc()
L2 = NodalLoc()
@test L1 == L2

printstyled("      push!\n", color=:light_green)
U = read_nodal("silixa", fstr)
S = deepcopy(U)
n = S.n
C = S[S.n]
push!(S, C)
@test S.n == n + 1
for f in SeisIO.Nodal.nodalfields
  F1 = getindex(getfield(S, f), n)
  F2 = getfield(C, f)
  @test F1 == F2
end
@test S.data[:, n+1] == S.data[:, n]
@test ===(S.data[:, n+1], S.data[:, n]) == false

D = convert(SeisChannel, C)
push!(S, D)

printstyled("      setindex!\n", color=:light_green)
S = deepcopy(U)
S2 = getindex(S, j1:j2)
setindex!(S, S2, j3:j4)
@test S.data[:,j1:j2] == S.data[:, j3:j4]

C = S[j2]
setindex!(S, C, j1)
@test S.data[:, j1] == S.data[:, j2]

printstyled("      show\n", color=:light_green)
redirect_stdout(out) do

  for i = 1:10
    for T in (NodalLoc, NodalChannel, NodalData)
      repr(T(), context=:compact=>true)
      repr(T(), context=:compact=>false)
      show(T())
    end
  end
end

printstyled("      sizeof\n", color=:light_green)
S = deepcopy(U)

# Add some notes
for j in 1:3
  j_min = rand(1:3)
  j_max = rand(j_min:S.n)
  note!(S, j_min:j_max, randstring(rand(1:100)))
end

# Add some things to :misc
for i in 1:4
  j_min = rand(1:3)
  j_max = rand(j_min:S.n)

  # add 1-6 entries each to the :misc field of each channel
  for j in j_min:j_max
    D = S.misc[j]
    pop_nodal_dict!(D, n=rand(1:6))
  end
end
@test sizeof(S) > 168

C = S[1]
note!(C, randstring(rand(1:100)))
note!(C, randstring(rand(1:100)))
pop_nodal_dict!(C.misc, n=rand(3:18))
@test sizeof(C) > 120

printstyled("      sort!\n", color=:light_green)
S = deepcopy(U)
id1 = S.id[1]
id2 = S.id[2]
S.id[1] = id2
S.id[2] = id1
sort!(S)
@test S.data[:,1] == U.data[:,2]
@test S.data[:,2] == U.data[:,1]
S.x[1][1] = v1
S.x[2][1] = v2
@test S.data[1,1] == v1
@test S.data[1,2] == v2


text:
# Ported from PRCcs/src/mcmc_gtb.py

function mcmc(a, b, phi, sst_df, n, ld_blk, blk_size, n_iter, n_burnin, thin, chrom, beta_std, seed; verbose=false)
    # seed
    if seed !== nothing
        Random.seed!(seed)
    end

    # derived stats
    beta_mrg = copy(sst_df.BETA)
    maf = copy(sst_df.MAF)
    n_pst = (n_iter-n_burnin)/thin
    p = length(sst_df.SNP)
    n_blk = length(ld_blk)

    # initialization
    beta = zeros(p)
    psi = ones(p)
    sigma = 1.0
    if phi === nothing
        phi = 1.0
        phi_updt = true
    else
        phi_updt = false
    end

    beta_est = zeros(p)
    psi_est = zeros(p)
    sigma_est = 0.0
    phi_est = 0.0

    # MCMC
    for itr in 1:n_iter
        if itr % 100 == 0
            verbose && @info "(Chromosome $chrom) MCMC iteration $itr"
        end

        mm = 1; quad = 0.0
        for kk in 1:n_blk
            if blk_size[kk] == 0
                continue
            else
                idx_blk = mm:(mm+blk_size[kk]-1)
                dinvt = ld_blk[kk] .+ Diagonal(1.0 ./ psi[idx_blk])
                dinvt_chol = cholesky(dinvt).U
                beta_tmp = (transpose(dinvt_chol) \ beta_mrg[idx_blk]) .+ sqrt(sigma/n) .* randn(length(idx_blk))
                beta[idx_blk] = dinvt_chol \ beta_tmp
                quad += dot(transpose(beta[idx_blk]) * dinvt, beta[idx_blk])
                mm += blk_size[kk]
            end
        end

        err = max(n/2.0*(1.0-2.0*sum(beta.*beta_mrg)+quad), n/2.0*sum(beta .^ 2 ./ psi))
        sigma = 1.0/rand(Gamma((n+p)/2.0, 1.0/err))

        delta = rand.(Gamma.(a+b, 1.0 ./ (psi .+ phi)))

        for jj in 1:p
            psi[jj] = gigrnd(a-0.5, 2.0*delta[jj], n*beta[jj]^2/sigma)
        end
        psi[psi .> 1] .= 1.0

        if phi_updt
            w = rand(Gamma(1.0, 1.0/(phi+1.0)))
            phi = rand(Gamma(p*b+0.5, 1.0/(sum(delta)+w)))
        end

        # posterior
        if (itr>n_burnin) && (itr % thin == 0)
            beta_est = beta_est + beta/n_pst
            psi_est = psi_est + psi/n_pst
            sigma_est = sigma_est + sigma/n_pst
            phi_est = phi_est + phi/n_pst
        end
    end

    # convert standardized beta to per-allele beta
    if !beta_std
        beta_est ./= sqrt.(2.0 .* maf .* (1.0 .- maf))
    end

    # print estimated phi
    if phi_updt && verbose
        @info @sprintf("Estimated global shrinkage parameter: %1.2e", phi_est)
    end

    return beta_est
end


text:
using KUparser, KUnet, HDF5, JLD, ArgParse, Compat, CUDArt

function fscore{T<:Parser}(pt::Type{T}, data::Vector{Corpus}, ndeps::Integer, feats::Fvec, args::Dict)
    # Assuming this is running on a worker, do not use ncpu when parsing.
    # Initialize training set using oparse on first corpus
    @date (p,x,y) = oparse(pt, data[1], ndeps, feats); p=nothing
    net = initnet(pt, data, ndeps, feats, args)
    @date for i=1:args["epochs"]
        train(net, x, y; batch=args["tbatch"], loss=KUnet.logploss, shuffle=true)
    end
    x = y = nothing
    if (args["parser"] == "gparser")
        @date p = gparse(pt, data[2], ndeps, feats, net, args["pbatch"])
    elseif (args["parser"] == "bparser")
        @date p = bparse(pt, data[2], ndeps, feats, net, args["nbeam"], args["pbatch"])
    else
        error("Unknown parser: "*args["parser"])
    end
    @show e = evalparse(p, data[2])
    net=p=nothing; gc()
    return e[1]  # e[1] is UAS including punct
end

function initnet{T<:Parser}(pt::Type{T}, data::Vector{Corpus}, ndeps::Integer, feats::Fvec, args::Dict)
    s1 = data[1][1]
    p1 = pt(wcnt(s1),ndeps)
    xrows=KUparser.flen(p1, s1, feats)
    yrows=p1.nmove
    # Initialize net:
    args["nogpu"] && KUnet.gpu(false)
    srand(42)
    net=newnet(relu, [xrows; args["hidden"]; yrows]...)
    net[end].f=KUnet.logp
    for k in [fieldnames(UpdateParam); :dropout]
        haskey(args, string(k)) || continue
        v = args[string(k)]
        if isempty(v)
            continue
        elseif length(v)==1
            setparam!(net, k, v[1])
        else
            @assert length(v)==length(net) "$k should have 1 or $(length(net)) elements"
            for i=1:length(v)
                setparam!(net[i], k, v[i])
            end
        end
    end
    return net
end

function flip(a, x)
    b = copy(a)
    i = findfirst(b, x)
    i == 0 ? push!(b, x) : deleteat!(b, i)
    return b
end

function flip!(a, x)
    i = findfirst(a, x)
    i == 0 ? push!(a, x) : deleteat!(a, i)
    return a
end


text:
# module Bukdu

export Changeset

"""
    Changeset

Used with `HTML5.Form.change`.
"""
mutable struct Changeset
    model::Type
    changes::NT where {NT<:NamedTuple}
    Changeset(model, changes = NamedTuple()) = new(model, changes)
end

# module Bukdu


text:
# Unit Tests of functions within src/sw_core.jl

# using c_sw

# using Test

function function_tests()
    @testset "c_sw function tests" begin
    # @test c_sw!()

    # @test divergence_corner()

    # @test d2a2c_vect()

    # @test edge_interpolate4()

    # @test fill2_4corners()

    # @test fill_4corners()

    end
end

@testset "All tests" begin
    function_tests()
end


text:
export parselog

function parselog(logfile::String, verbose::Bool)

    verbose && println("CPUMonitor: Begin Parsing Log...")
    # Read log file as array with " " delim
    log = readdlm(logfile)

    # Find unique process IDs
    pids = unique(log[:,1])
    count = 1

    stat = Array{CPUStat,1}(length(pids))

    # Collect stats from each PID into an instance of CPUStat
    for pid in pids

        inds = find(x -> x == pid, log[:,1])

        stat[count] = CPUStat(pid, log[inds, 9], log[inds, 10])

        count += 1

    end

    verbose && println("CPUMonitor: End Parsing Log...")

    return stat
end


text:
module Magnetic

using ModelingToolkit, Symbolics, IfElse, OrdinaryDiffEq

@parameters t
D = Differential(t)

include("utils.jl")

include("FluxTubes/sensors.jl")

include("QuasiStatic/FluxTubes/basic.jl")
include("QuasiStatic/FluxTubes/sources.jl")
include("QuasiStatic/FluxTubes/sensors.jl")

export MagneticPort,
    PositiveMagneticPort,
    NegativeMagneticPort,
    Ground,
    TwoPortElementary,
    TwoPortExtended,
    TwoPort,
    Idle,
    Short,
    Crossing,
    ConstantPermeance,
    ConstantReluctance,
    ConstantMagneticPotentialDifference,
    ConstantMagneticFlux,
    AbsoluteSensor,
    # FluxTubes sensors
    MagneticFluxSensor, MagneticPotentialDifferenceSensor

end #module


text:
"""
    SeisMute(d; <keyword arguments>) -> Array{Real,2}

Mute trace noise before the first arrival.
The muting function depends on the offset between source and receiver and the velocity of the first layer as
t_mute = t_0 + offset/v1.

# Arguments
- `d`: Array of input traces.

# Keyword arguments
- `offset=[0.]`: Vector of distances between source and receiver.
- `tmute=0.`: initial muting time
- `vmute=1500.`: Velocity of the first layer.
- `taper=0.1`: Taper of the muting function
- `dt::Real=0.004`: sampling interval in secs.

"""
function SeisMute(in;offset=[0.],tmute=0.,vmute=1500.,taper=0.1,dt=0.004)

	vec= size(in)
	in = reshape(in,vec[1],:)
	nt,nx = size(in)
	out = copy(in)
	for it = 1:nt
		for ix = 1:nx
			t = sqrt(tmute^2 + (offset[ix]/vmute).^2)
			if ((it-1)*dt <= t )
					if ((it-1)*dt > t - taper)
						out[it,ix] *= 1 - (t - (it-1)*dt)/taper
					else
						out[it,ix] = 0.
					end
			end
		end
	end
	out = reshape(out,vec)
	return out
end

"""
    SeisMute(in,out,parameters; <keyword arguments>)

Mutes noise before first arrival of traces in a file. Saves the muted traces to an output file.
The muting function depends on the offset between source and receiver and the velocity of the first layer as
t_mute = t_0 + offset/v1.

# Arguments
- `in::String`: Input file - Seis format.
- `out::String`: Output file - Seis format.
- `parameters` : list of the keyword arguments for the function SeisMute.

# Keyword arguments
- `group="gather"` : Options are all, some or gather
- `key=["imx","imy"]` : Defines type of gather
- `itrace=1` : Initial trace number
- `ntrace=10000` : Total number of traces to process at once

# Example
```
julia> param = Dict(:tmute=>0.0, :vmute=>10000, :taper=>0.05,:dt=>0.01)
julia> SeisMute(filein,fileout, param,group="some")
```

"""
function SeisMute(in::String,out::String,parameters;group="gather",key=["imx","imy"],itrace=1,ntrace=10000)

	tmute = get(parameters,:tmute,0.)
	vmute = get(parameters,:vmute,1500.)
	taper = get(parameters,:taper,0.1)

	if (group=="all")
		headers = SeisMain.SeisReadHeaders(in);
		offset = SeisMain.ExtractHeader(headers,"h")
		dt = headers[1].d1
		parameters = Dict(:offset=>offset,:tmute=>tmute,:vmute=>vmute,:taper=>taper,:dt=>dt)

		SeisProcessFile(in,out,[SeisMute],[parameters];group=group)
	else
		itrace_in = itrace #1
		#itrace_out = 1
		nx = SeisMain.GetNumTraces(in)
		while itrace_in <= nx
			h = SeisMain.SeisReadHeaders(in,group=group,key=key,itrace=itrace_in,ntrace=ntrace)
			offset = SeisMain.ExtractHeader(h,"h")
			dt = h[1].d1
			#num_traces = size(h,1)
			parameters = Dict(:offset=>offset,:tmute=>tmute,:vmute=>vmute,:taper=>taper,:dt=>dt)
			#SeisProcessFile(in,out,[SeisMute],[parameters];group=group,key=key,itrace=itrace_in,ntrace=ntrace)
			d1,h1,e1 = SeisMain.SeisRead(in,group=group,key=key,itrace=itrace_in,ntrace=ntrace)

			d2 = SeisMute(d1;parameters...)
			d1 = copy(d2)

			SeisMain.SeisWrite(out,d1,h1,e1,itrace=itrace_in)

			#println(itrace_in," ",ntrace)

			itrace_in += ntrace

			#itrace_in += num_traces
			#itrace_out += num_traces
		end
	end
end


text:
global const k_B        = 1.38066e-16
global const c_light    = 2.9979e10
global const σ_T        = 6.65245e-25
global const m_e        = 9.10953e-28
global const m_p        = 1.6726e-24
global const h_planck   = 6.6261e-27
global const q_e        = 1.602176487e-20 * c_light
global const yPrefac    = σ_T * k_B / (m_e * c_light^2)
global const eV2cgs     = 1.60218e-12
global const cgs2eV     = 1.0/eV2cgs
global const C_crit     = 3q_e / ( 4π * m_e * c_light ) # Donnert+16, MNRAS 462, 2014–2032 (2016), Eg. 20
                                                        #  -> converted to dimensionless momentum
global const γ          = 5.0/3.0
global const mJy_factor = 1.e26       # conversion factor from [erg/cm^3/Hz/s] to mJy/cm.


text:
using Plots

include("QuadraturePoints.jl")
include("element_matrices.jl")

#Some Constants
DFloat = Float64
machine_zero=eps(DFloat)

#ipoints=1: Lobatto
#ipoints=2: Legendre
#ipoints=3: Chebyshev
#ipoints=4: Equi-spaced

#{{{ Main
function main()

    #Fix Interpolation and Integration Order
    N=1
    Q=N+1
    Np=N+1
    Nq=Q+1
    ipoints=1
    qpoints=1

    #Allocate Arrays
    M=zeros(DFloat,Np,Np)
    D=zeros(DFloat,Np,Np)

    #Select Interpolation Points
    if ipoints == 1
        (ξ,ω) = QuadraturePoints.lobatto_gauss(Np,DFloat) #Lobatto
    elseif ipoints == 2
        (ξ,ω) = QuadraturePoints.legendre_gauss(Np,DFloat) #Legendre
    elseif ipoints == 3
        (ξ,ω) = QuadraturePoints.chebyshev_gauss(Np,DFloat) #Chebyshev
    elseif ipoints == 4
        (ξ,ω) = QuadraturePoints.equispaced_points(Np,DFloat) #Equi-spaced
    end

    #Select Quadrature Points
    if qpoints == 1
        (ξq,ωq) = QuadraturePoints.lobatto_gauss(Nq,DFloat) #Lobatto
    elseif qpoints == 2
        (ξq,ωq) = QuadraturePoints.legendre_gauss(Nq,DFloat) #Legendre
    end

    #Construct Lagrange Polynomials
    (ψ,dψ) = QuadraturePoints.lagrange_basis(Np,Nq,ξ,ξq)

    (M,D) = element_matrices(ψ,dψ,Np,Nq,ωq,DFloat)

    #Check correctness
    for j=1:Np, i=1:Np
        @show(i,j,M[i,j])
    end
    row_sum=zeros(DFloat,Np)
    for i=1:Np
        row_sum[i]=sum(D[i,:])
    end
    @show(row_sum)
    @show(sum(M))

    #Plot Interpolation
    println("Done") #output

end
#}}} Main

#----------------------------------#
# Run the main function
#----------------------------------#
main()


text:
#######################
# static DSL function #
#######################

using FunctionalCollections: PersistentVector

struct Params
    prob_outlier::Float64
    slope::Float64
    intercept::Float64
    inlier_std::Float64
    outlier_std::Float64
end

@gen (static) function datum(x::Float64, (grad)(params::Params))
    is_outlier = @trace(bernoulli(params.prob_outlier), :z)
    std = is_outlier ? params.outlier_std : params.inlier_std
    y = @trace(normal(x * params.slope + params.intercept, std), :y)
    return y
end

data_fn = Map(datum)

"""
my documentation
"""
@gen (static) function model(xs::Vector{Float64})
    n = length(xs)
    inlier_std = @trace(gamma(1, 1), :inlier_std)
    outlier_std = @trace(gamma(1, 1), :outlier_std)
    slope = @trace(normal(0, 2), :slope)
    intercept = @trace(normal(0, 2), :intercept)
    params::Params = Params(0.5, slope, intercept, inlier_std, outlier_std)
    ys = @trace(data_fn(xs, fill(params, n)), :data)
    return ys
end

@gen (static) function at_choice_example_1(i::Int)
    ret = @trace(bernoulli(0.5), :x => i)
end

# @trace(choice_at(bernoulli)(0.5, i), :x)

@gen (static) function at_choice_example_2(i::Int)
    ret = @trace(bernoulli(0.5), :x => i => :y)
end

# @trace(call_at(choice_at(bernoulli))(0.5, i, :y), :x)

@gen function foo(mu)
    @trace(normal(mu, 1), :y)
end

@gen (static) function at_call_example_1(i::Int)
    mu = 1.123
    ret = @trace(foo(mu), :x => i)
end

@gen (static) function at_call_example_2(i::Int)
    mu = 1.123
    ret = @trace(foo(mu), :x => i => :y)
end

# Modules to test load_generated_functions
module MyModuleA
using Gen
@gen (static) function foo(mu)
    y = @trace(normal(mu, 1), :y)
    return y
end
load_generated_functions(@__MODULE__)
end

module MyModuleB
using Gen
@gen (static) function foo(mu)
    y = @trace(normal(mu, 1), :y)
    return y
end
@load_generated_functions
end

@testset "static DSL" begin

function get_node_by_name(ir, name::Symbol)
    nodes = filter((node) -> (node.name == name), ir.nodes)
    @assert length(nodes) == 1
    nodes[1]
end

function get_node_by_addr(ir, addr::Symbol)
    nodes = filter(n -> hasfield(typeof(n), :addr) && (n.addr==addr), ir.nodes)
    @assert length(nodes) == 1
    nodes[1]
end

@testset "check IR" begin

#####################
# check IR of datum #
#####################

ir = Gen.get_ir(typeof(datum))
# argument nodes
@test length(ir.arg_nodes) == 2
x = ir.arg_nodes[1]
params = ir.arg_nodes[2]
@test x.name == :x
@test x.typ == :Float64
@test !x.compute_grad
@test params.name == :params
@test params.typ == :Params
@test params.compute_grad

# choice nodes and call nodes
@test length(ir.choice_nodes) == 2
@test length(ir.call_nodes) == 0

# is_outlier
is_outlier = ir.choice_nodes[1]
@test is_outlier.addr == :z
@test is_outlier.typ == QuoteNode(Bool)
@test is_outlier.dist == bernoulli
@test length(is_outlier.inputs) == 1

# std
std = get_node_by_name(ir, :std)
@test isa(std, Gen.JuliaNode)
@test std.name == :std
@test std.typ == QuoteNode(Any)
@test length(std.inputs) == 2
in1 = std.inputs[1]
in2 = std.inputs[2]
@test (in1 === is_outlier && in2 === params) || (in2 === is_outlier && in1 === params)

# y
y = ir.choice_nodes[2]
@test y.addr == :y
@test y.typ == QuoteNode(Float64)
@test y.dist == normal
@test length(y.inputs) == 2
@test y.inputs[2] === std

# y_mean
y_mean = y.inputs[1]
@test isa(y_mean, Gen.JuliaNode)
@test y_mean.typ == QuoteNode(Any)
@test length(y_mean.inputs) == 2
in1 = y_mean.inputs[1]
in2 = y_mean.inputs[2]
@test (in1 === x && in2 === params) || (in2 === x && in1 === params)

# prob outlier
prob_outlier = is_outlier.inputs[1]
@test isa(prob_outlier, Gen.JuliaNode)
@test length(prob_outlier.inputs) == 1
@test prob_outlier.inputs[1] === params
@test prob_outlier.typ == QuoteNode(Any)

@test ir.return_node === y

#####################
# check IR of model #
#####################

ir = Gen.get_ir(typeof(model))
@test length(ir.arg_nodes) == 1
xs = ir.arg_nodes[1]
@test xs.name == :xs
@test xs.typ == :(Vector{Float64})
@test !xs.compute_grad

# choice nodes and call nodes
@test length(ir.choice_nodes) == 4
@test length(ir.call_nodes) == 1

# inlier_std
inlier_std = ir.choice_nodes[1]
@test inlier_std.addr == :inlier_std
@test inlier_std.typ == QuoteNode(Float64)
@test inlier_std.dist == gamma
@test length(inlier_std.inputs) == 2

# outlier_std
outlier_std = ir.choice_nodes[2]
@test outlier_std.addr == :outlier_std
@test outlier_std.typ == QuoteNode(Float64)
@test outlier_std.dist == gamma
@test length(outlier_std.inputs) == 2

# slope
slope = ir.choice_nodes[3]
@test slope.addr == :slope
@test slope.typ == QuoteNode(Float64)
@test slope.dist == normal
@test length(slope.inputs) == 2

# intercept
intercept = ir.choice_nodes[4]
@test intercept.addr == :intercept
@test intercept.typ == QuoteNode(Float64)
@test intercept.dist == normal
@test length(intercept.inputs) == 2

# data
ys = ir.call_nodes[1]
@test ys.addr == :data
@test ys.typ == QuoteNode(PersistentVector{Float64})
@test ys.generative_function == data_fn
@test length(ys.inputs) == 2
@test ys.inputs[1] == xs

# params
params = get_node_by_name(ir, :params)
@test isa(params, Gen.JuliaNode)
@test params.name == :params
@test params.typ == :Params
@test length(params.inputs) == 4
@test slope in params.inputs
@test intercept in params.inputs
@test inlier_std in params.inputs
@test outlier_std in params.inputs

# n
n = get_node_by_name(ir, :n)
@test isa(n, Gen.JuliaNode)
@test n.name == :n
@test n.typ == QuoteNode(Any)
@test length(n.inputs) == 1
@test n.inputs[1] === xs

# params_vec
params_vec = ys.inputs[2]
@test isa(params_vec, Gen.JuliaNode)
@test params_vec.typ == QuoteNode(Any)
@test length(params_vec.inputs) == 2
in1 = params_vec.inputs[1]
in2 = params_vec.inputs[2]
@test (in1 === params && in2 === n) || (in2 === params && in1 === n)

@test ir.return_node === ys
end

@testset "at_choice" begin

# at_choice_example_1
ir = Gen.get_ir(typeof(at_choice_example_1))
i = get_node_by_name(ir, :i)
ret = get_node_by_addr(ir, :x)
@test isa(ret, Gen.GenerativeFunctionCallNode)
@test ret.addr == :x
@test length(ret.inputs) == 2
@test isa(ret.inputs[1], Gen.JuliaNode) # () -> 0.5
@test ret.inputs[2] === i
at = ret.generative_function
@test isa(at, Gen.ChoiceAtCombinator)
@test at.dist == bernoulli

# at_choice_example_2
ir = Gen.get_ir(typeof(at_choice_example_2))
i = get_node_by_name(ir, :i)
ret = get_node_by_addr(ir, :x)
@test isa(ret, Gen.GenerativeFunctionCallNode)
@test ret.addr == :x
@test length(ret.inputs) == 3
@test isa(ret.inputs[1], Gen.JuliaNode) # () -> 0.5
@test isa(ret.inputs[2], Gen.JuliaNode) # () -> :y
@test ret.inputs[3] === i
at = ret.generative_function
@test isa(at, Gen.CallAtCombinator)
at2 = at.kernel
@test isa(at2, Gen.ChoiceAtCombinator)
@test at2.dist == bernoulli
end

@testset "at_call" begin

# at_call_example_1
ir = Gen.get_ir(typeof(at_call_example_1))
i = get_node_by_name(ir, :i)
ret = get_node_by_addr(ir, :x)
@test isa(ret, Gen.GenerativeFunctionCallNode)
@test ret.addr == :x
@test length(ret.inputs) == 2
@test isa(ret.inputs[1], Gen.JuliaNode) # () -> 0.5
@test ret.inputs[2] === i
at = ret.generative_function
@test isa(at, Gen.CallAtCombinator)
@test at.kernel == foo

#at_call_example_2
ir = Gen.get_ir(typeof(at_call_example_2))
i = get_node_by_name(ir, :i)
ret = get_node_by_addr(ir, :x)
@test isa(ret, Gen.GenerativeFunctionCallNode)
@test ret.addr == :x
@test length(ret.inputs) == 3
@test isa(ret.inputs[1], Gen.JuliaNode) # () -> 0.5
@test isa(ret.inputs[1], Gen.JuliaNode) # () -> 0.5
@test isa(ret.inputs[2], Gen.JuliaNode) # () -> :y
@test ret.inputs[3] === i
at = ret.generative_function
@test isa(at, Gen.CallAtCombinator)
at2 = at.kernel
@test isa(at2, Gen.CallAtCombinator)
@test at2.kernel == foo
end

@testset "tracked diff annotation" begin

@gen (static, diffs) function bar(x)
    return x
end

@test Gen.get_options(typeof(bar)).track_diffs

@gen (static) function bar(x)
    return x
end

@test !Gen.get_options(typeof(bar)).track_diffs

end

@testset "no julia cache annotation" begin

@gen (static, nojuliacache) function bar(x)
    return x
end

@test !Gen.get_options(typeof(bar)).cache_julia_nodes

@gen (static) function bar(x)
    return x
end

@test Gen.get_options(typeof(bar)).cache_julia_nodes

end

@testset "trainable parameters" begin

@gen (static) function foo()
    @param theta::Float64
    return theta
end

ir = Gen.get_ir(typeof(foo))
theta = get_node_by_name(ir, :theta)
@test isa(theta, Gen.TrainableParameterNode)
@test ir.return_node === theta

end

@testset "returning a trace directly" begin

@gen (static) function f1()
    x = @trace(normal(0, 1), :foo)
    return x
end

@gen (static) function f2()
    return @trace(normal(0, 1), :foo)
end

ir1 = Gen.get_ir(typeof(f1))
ir2 = Gen.get_ir(typeof(f2))
return_node1 = ir1.return_node
return_node2 = ir2.return_node
@test isa(return_node2, typeof(return_node1))
@test return_node2.dist == return_node1.dist

inputs1 = return_node1.inputs
inputs2 = return_node2.inputs
@test 0 == inputs1[1].fn() == inputs2[1].fn()
@test 1 == inputs1[2].fn() == inputs2[2].fn()

@test return_node2.name != return_node1.name
@test return_node2.addr == return_node1.addr
@test return_node2.typ === return_node1.typ

end

@testset "@trace within Julia expressions" begin

@gen (static) function circle(r, t, std)
    pt = (@trace(normal(r*cos(t), std), :x), @trace(normal(r*sin(t), std), :y))
    return pt
end

@gen (static) function beta_bernoulli()
    return @trace(bernoulli(@trace(beta(1, 1), :z)), :x)
end

load_generated_functions()

x, y = circle(2, pi/2, 0)
@test isapprox(x, 0.0, atol=1e-8) && isapprox(y, 2.0, atol=1e-8)
tr, w = generate(circle, (1, 0, 1), choicemap(:x => 1, :y => 1))
@test isapprox(w, logpdf(normal, 1, 1, 1) + logpdf(normal, 1, 0, 1))

tr, w = generate(beta_bernoulli, (), choicemap(:z => 1.0))
@test tr[:x] == get_retval(tr) == true
tr, w = generate(beta_bernoulli, (), choicemap(:z => 0.0))
@test tr[:x] == get_retval(tr) == false
tr, w = generate(beta_bernoulli, (), choicemap(:z => 0.5, :x => true))
@test isapprox(w, logpdf(beta, 0.5, 1, 1) + logpdf(bernoulli, true, 0.5))

end

@testset "repeated assignments" begin

@gen (static) function foo(x)
    x = x + 1
    return x
end

@gen (static) function bar(x)
    y = x + 1
    y = y + 1
    return y + 1
end

load_generated_functions()

ir = Gen.get_ir(typeof(foo))
x1 = get_node_by_name(ir, :x)
x2 = ir.return_node
@test isa(x1, Gen.ArgumentNode)
@test isa(x2, Gen.JuliaNode)
@test x1 !== x2
@test x2.inputs[1] === x1
@test foo(0) == 1

ir = Gen.get_ir(typeof(bar))
y1 = get_node_by_name(ir, :y)
y2 = ir.return_node.inputs[1]
@test y1 !== y2
@test y2.inputs[1] === y1
@test bar(0) == 3

end

@testset "tuple assignments" begin

@gen (static) function foo(params::Tuple)
    p, (mu1, sigma1), (mu2, sigma2) = params
    mu, sigma = @trace(bernoulli(p), :x) ? (mu1, sigma1) : (mu2, sigma2)
    return @trace(normal(mu, sigma), :y)
end

ir = Gen.get_ir(typeof(foo))
params = get_node_by_name(ir, :params)
p = get_node_by_name(ir, :p)
mu1 = get_node_by_name(ir, :mu1)
sigma2 = get_node_by_name(ir, :sigma2)
@test p.inputs[1] === params
@test mu1.inputs[1].inputs[1] === params
@test sigma2.inputs[1].inputs[1] === params

return_node = ir.return_node
inputs = return_node.inputs
@test inputs[1] === get_node_by_name(ir, :mu)
@test inputs[2] === get_node_by_name(ir, :sigma)
@test inputs[1].fn((:fst, :snd)) == :fst
@test inputs[2].fn((:fst, :snd)) == :snd

load_generated_functions()

params = (0.5, (0, 1), (1, 2))
tr, w = generate(foo, (params,), choicemap(:x => true, :y => 1))
@test get_retval(tr) == 1
@test isapprox(w, logpdf(bernoulli, true, 0.5) + logpdf(normal, 1, 0, 1))

tr, w = generate(foo, (params,), choicemap(:x => false, :y => 1))
@test get_retval(tr) == 1
@test isapprox(w, logpdf(bernoulli, false, 0.5) + logpdf(normal, 1, 1, 2))

end

@testset "use of 'end'" begin

@gen (static) function foo()
    x = [1, 2, 3, 4, 5, 6]
    y = x[3:end]
    return y
end

load_generated_functions()

@test foo() == [3, 4, 5, 6]

end

@testset "'return nothing'" begin

@gen (static) function foo()
    return nothing
end

load_generated_functions()

@test foo() == nothing

end

@testset "getindex(trace)" begin

@gen (static) function bar(r)
    a = @trace(normal(0, 1), :a)
    return r
end
@gen (static) function foo()
    x = @trace(bar(1), :x)
    yz = @trace(bar(2), :y => :z)
    u = @trace(normal(0, 1), :u)
    vw = @trace(normal(0, 1), :v => :w)
    ret = 7
    return ret
end

Gen.load_generated_functions()

constraints = choicemap()
constraints[:u] = 1.1
constraints[:v => :w] = 1.2
constraints[:x => :a] = 1.3
constraints[:y => :z => :a] = 1.4
trace, = generate(foo, (), constraints)

# random choices
@test trace[:u] == 1.1
@test trace[:v => :w] == 1.2
@test trace[:x => :a] == 1.3
@test trace[:y => :z => :a] == 1.4

# auxiliary state
@test trace[:x] == 1
@test trace[:y => :z] == 2

# return value
@test trace[] == 7

end

@testset "docstrings" begin
    io = IOBuffer()
    print(io, @doc model)
    @test String(take!(io)) == "my documentation\n"
end

@testset "one-line definitions" begin

@gen (static) foo(x) = (y = @trace(normal(x, 1), :y); return y)

load_generated_functions()

tr, w = generate(foo, (0,), choicemap(:y => 1))
@test get_retval(tr) == 1
@test isapprox(w, logpdf(normal, 1, 0, 1))

end

@testset "module-defined static functions" begin

tr, w = generate(MyModuleA.foo, (0,), choicemap(:y => 1))
@test get_retval(tr) == 1
@test isapprox(w, logpdf(normal, 1, 0, 1))

tr, w = generate(MyModuleB.foo, (0,), choicemap(:y => 1))
@test get_retval(tr) == 1
@test isapprox(w, logpdf(normal, 1, 0, 1))

end

@testset "static gen function choicemaps" begin
@gen (static) function bar2()
    b ~ normal(0, 1)
    return b
end
@gen (static) function bar1()
    a ~ bar2()
    x ~ normal(0, 1)
    return x
end
Gen.load_generated_functions()
tr = simulate(bar1, ())
ch = get_choices(tr)
@test has_value(ch, :x)
@test !has_value(ch, :y)
@test_throws KeyError get_submap(ch, :x)
@test has_value(get_submap(ch, :a), :b)
@test get_submap(ch, :y) == EmptyChoiceMap()
@test length(get_values_shallow(ch)) == 1
@test length(get_submaps_shallow(ch)) == 1
end

end # @testset "static DSL"

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
using BenchmarkTools
using CSV
using ProgressLogging
using PSFModels: Gaussian
using PythonCall
using Random
using Statistics
using SubpixelRegistration

skr = pyimport("skimage.registration")

rng = Xoshiro(1134)

# image dimension sizes
SIZES = 50:100:1000
UPSAMPLE_FACTORS = (1, 100, 1000)

rows = []
@progress "sizes" for N in SIZES
    center = (N + 1) / 2
    img = Gaussian(x=center, y=center, fwhm=10)[1:N, 1:N]
    shift = 10 .* randn(rng, 2)
    img_shifted = fourier_shift(img, shift)
    @progress "upsample factors" for ups in UPSAMPLE_FACTORS
        bench_jl = @benchmark phase_offset($img, $img_shifted, upsample_factor=$ups)
        bench_py = @benchmark skr.phase_cross_correlation($img, $img_shifted, upsample_factor=$ups, return_error=true)
        bundle = (;
            size=N, upsample_factor=ups,
            time_jl=median(bench_jl).time/1e9, time_py=median(bench_py).time/1e9
        )
        @info "" bundle...
        push!(rows, bundle)
    end
end

outpath = joinpath(dirname(pathof(SubpixelRegistration)), "..", "bench", "benchmark_results.csv")
CSV.write(outpath, rows)


text:
abstract type AbstractGame end

abstract type AbstractContainer end

abstract type AbstractPlayer <: AbstractContainer end
abstract type AbstractBag <: AbstractContainer end

abstract type AbstractRound end
abstract type AbstractState end
abstract type AbstractClock end

abstract type AbstractLog end
abstract type AbstractGrid end

abstract type AbstractBoard end
abstract type AbstractEntry end

abstract type AbstractPiece end
abstract type AbstractBlock end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""CrossVersionObjectReference contains enough information to let you identify the referred resource.

    IoK8sApiAutoscalingV2beta1CrossVersionObjectReference(;
        apiVersion=nothing,
        kind=nothing,
        name=nothing,
    )

    - apiVersion::String : API version of the referent
    - kind::String : Kind of the referent; More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds\"
    - name::String : Name of the referent; More info: http://kubernetes.io/docs/user-guide/identifiers#names
"""
mutable struct IoK8sApiAutoscalingV2beta1CrossVersionObjectReference <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    name::Any # spec type: Union{ Nothing, String } # spec name: name

    function IoK8sApiAutoscalingV2beta1CrossVersionObjectReference(;apiVersion=nothing, kind=nothing, name=nothing)
        o = new()
        validate_property(IoK8sApiAutoscalingV2beta1CrossVersionObjectReference, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiAutoscalingV2beta1CrossVersionObjectReference, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiAutoscalingV2beta1CrossVersionObjectReference, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        o
    end
end # type IoK8sApiAutoscalingV2beta1CrossVersionObjectReference

const _property_map_IoK8sApiAutoscalingV2beta1CrossVersionObjectReference = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("name")=>Symbol("name"))
const _property_types_IoK8sApiAutoscalingV2beta1CrossVersionObjectReference = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("name")=>"String")
Base.propertynames(::Type{ IoK8sApiAutoscalingV2beta1CrossVersionObjectReference }) = collect(keys(_property_map_IoK8sApiAutoscalingV2beta1CrossVersionObjectReference))
Swagger.property_type(::Type{ IoK8sApiAutoscalingV2beta1CrossVersionObjectReference }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAutoscalingV2beta1CrossVersionObjectReference[name]))}
Swagger.field_name(::Type{ IoK8sApiAutoscalingV2beta1CrossVersionObjectReference }, property_name::Symbol) =  _property_map_IoK8sApiAutoscalingV2beta1CrossVersionObjectReference[property_name]

function check_required(o::IoK8sApiAutoscalingV2beta1CrossVersionObjectReference)
    (getproperty(o, Symbol("kind")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAutoscalingV2beta1CrossVersionObjectReference }, name::Symbol, val)
end


text:
"""
    minimum_weight_perfect_matching(g, weights::AEdgeMap{T}; cutoff=typemax{T}) where {T<:Real}

Given a graph `g` and an edgemap `weights` containing non-negative weights associated
to edges, returns a matching with the mimimum total weight among the ones containing
exactly `nv(g)/2` edges.

Edges in `g` not present in `weights` will not be considered for the matching.
To reduce computational time, a `cutoff` argument can be given. Only edges
with weight lower than `cutoff` will be considered for the matching.

This function relies on the BlossomV.jl package, a julia wrapper
around Kolmogorov's BlossomV algorithm.

The **return value** is a `NamedTuple` object `match` with fields:

    match.weight: total weight of the matching

    match.mate:    `mate[i] = j` if vertex `i` is matched to vertex `j`.
                   `mate[i] = -1` for unmatched vertices.

**Example**
```juliarepl
g = random_regular_graph(30, 4)
w = EdgeMap(g, e -> rand())
match = minimum_weight_perfect_matching(g, w)
print(match.weight)
print(match.mate)
```
"""
function minimum_weight_perfect_matching(
        g::AGraph,
        w::AEdgeMap{T};
        cutoff = typemax(T)) where {T<:Real}

    m = BlossomV.Matching(T, nv(g))
    for e in edges(g)
        c = get(w, e, cutoff)
        if c < cutoff
            BlossomV.add_edge(m, src(e)-1, dst(e)-1, c)
        end
    end
    BlossomV.solve(m)

    mate = fill(-1, nv(g))
    totweight = T(0)
    for i=1:nv(g)
        j = BlossomV.get_match(m, i-1) + 1
        mate[i] = j <= 0 ? -1 : j
        if i < j
            totweight += w[i, j]
        end
    end
    return (weight=totweight, mate=mate)
end


text:
using SymbolicTensors
using SymPy
using Test
using LinearAlgebra

@testset "SymbolicTensors.jl" begin
    @testset "types" begin
        field = TensorIndexType("field","f")
        @test typeof(field) == TensorIndexType
        @indices field i j
        @test typeof(i) == typeof(j) == TensorIndex
        A = TensorHead("A",[field])
        B = TensorHead("B",[field])
        @test typeof(A) == TensorHead
        @test typeof(A(i)) == IndexedTensor
        @test typeof(A(i)*A(-i)) == SymPy.Sym
        @test typeof(A(i) + B(i)) == TensAdd
        @test typeof(A(i) * B(j)) == TensMul
    end
    @testset "arithmetic" begin
        field = TensorIndexType("field","f")
        @indices field i j k l m n
        A = TensorHead("A",[field])
        B = TensorHead("B",[field])
        C = TensorHead("C",[field,field],get_tsymmetry.fully_symmetric(-2))
        @test A(i) + A(i) == 2*A(i)
        @test A(i) - A(i) == 0
        @test canon_bp(A(i)*A(j) + A(j)*A(i)) == 2*A(i)*A(j)
        @test canon_bp(C(-i,-j)*A(i)*A(j)) == 0
        # -- odd numbers of antisymmetric tensors contracted give zero
        @test canon_bp(C(i,-j)*C(j,-k)*C(k,-i)) == 0
        @test scalarIsEqual(A(i)*A(j)*field.metric(-i,-j),A(i)*A(-i),field.metric)
        @test scalarIsEqual(A(i)*B(j)*field.metric(-i,-j),A(i)*B(-i),field.metric)
        # should check properties of deltas, epsilons here too
        # -- should add some from xTensor
    end
    # -- should calculate some known Riemann scalars
    @testset "derivatives" begin
        field = TensorIndexType("field","f")
        @indices field i j k
        A = TensorHead("A",[field])
        @test diff(A(i),A(j)) == field.metric(i,-j)
        @test diff(A(-i),A(j)) == field.metric(-i,-j)
        @test diff(A(i),A(-j)) == field.metric(i,j)
        g = 1/(1+A(i)*A(-i))^2 * field.delta(-i,-j)
        gg = diff(g(-i,-j),A(k))
        @test diff(field.delta(-i,-j),A(k)) == 0
        # product rule w/ scalar chain rule
        @test scalarIsEqual( diff( A(i)*A(-i)*A(j), A(k)),
            A(i)*A(-i)*field.metric(j,-k) + 2*A(j)*A(-k) ,field.metric)
        # distribute over TensAdd
        @test canon_bp(diff( (A(i)*A(j) + A(j)*A(i))/2, A(k))) ==
            canon_bp(A(i)*field.metric(j,-k) + A(j)*field.metric(i,-k))
        #@test gg(-i,-j,-k) ==
        #h = canon_bp(1//2*(gg(-i,-j,-k) + gg(-i,-k,-j) - gg(-j,-k,-i) ))
        #GG = TensorHead("GG",[field,field,field])
        #h = canon_bp(1//2*(GG(-i,-j,-k) + GG(-i,-k,-j) - GG(-j,-k,-i) ))
        #@test h(i,-j,-k) == h(i,-k,-j)
    end

    @testset "errors" begin
    end

    @testset "replacements" begin
        field = TensorIndexType("field","f")
        @indices field i j k
        A = TensorHead("A",[field])
        field.set_metric(field.delta)
        @test replace_with_arrays(A(i)*A(j)*field.metric(-i,-j), Dict(field.delta(-i,-j)=>[1 0 ; 0 1], A(i)=>[0,4] )) == 16
    end

    @testset "quoting" begin
        vars = symbols("a b c d e f")
        a,b,c,d,e,f = vars
        function equalQuotes(q1,q2)
            if q1.head == :function
                q1 = q1.args[2]
                q2 = q2.args[2]
            end
            global a,b,c,d,e,f
            a,b,c,d,e,f = randn(6)
            return eval(q1) == eval(q2)
        end
        #Quote a scalar
        @test Quote(vars[1]) == convert(Expr,vars[1])
        cse_ex = sympy.cse(a^2 + 1/(1+a^2))
        @test equalQuotes(Quote(cse_ex),quote
              x0 = a ^ 2
              x0 + (1 + x0) ^ -1
        end)
        @test equalQuotes(Quote("hub",cse_ex,[a]),:(function hub(a)
              x0 = a ^ 2
              x0 + (1 + x0) ^ -1
        end))
        arr = [a^2,b^2+a^2,3c^2]
        cse_arr = sympy.cse(arr)
        @test equalQuotes(Quote(arr),quote
            reshape([a ^ 2, a ^ 2 + b ^ 2, __prod__(3, c ^ 2)], (3,))
        end)
        @test equalQuotes(Quote(cse_arr),quote
            x0 = a ^ 2
            reshape([x0, x0 + b ^ 2, __prod__(3, c ^ 2)], (3, 1))
        end)
        @test equalQuotes(Quote("_A",cse_arr,[a,b,c]),:(function _A(a, b, c)
            x0 = a ^ 2
            reshape([x0, x0 + b ^ 2, __prod__(3, c ^ 2)], (3, 1))
        end))
    end
    @testset "Riemann scalar computations" begin
        field = TensorIndexType("field","f")
        @indices field i j k l m n o
        F = TensorHead("F",[field])
        δ = field.delta

        simp(x) = contract_metric(canon_bp(x),field.metric)

        α = symbols("a",positive=true)
        g = 6 * α * δ(-i,-j) / (1-F(i)*F(j)*δ(-i,-j))^2
        gg = diff(g(-i,-j),F(k))
        gg = simp(gg)

        h = (gg(-i,-j,-k) - gg(-j,-k,-i) + gg(-k,-i,-j))/2
        h =  h(-i,-j,-k)
        h = simp(h)

        dh = diff(h(-i,-j,-k),F(l))
        dh = simp(dh)

        hh = (field.metric(m,n)*h(-m,-i,-j))*h(-n,-k,-l)
        hh = simp(hh)

        Riemann = dh(-i,-l,-j,-k) - dh(-i,-k,-l,-j) + hh(-j,-k,-i,-l) - hh(-j,-l,-i,-k)
        Riemann = simp(Riemann)

        Ric = Riemann(-i,-j,-k,-l)*field.metric(i,k)
        Ric = simp(Ric)

        R = Ric(-i,-j)*field.metric(j,i)
        R = simp(R)

        ## turn into arrays
        vars = symbols("x y z",real=true)
        x,y,z = vars
        d = length(vars)
        repl = Dict{Any,Any}(F(m) => vars, field.delta(-m,-n) => Matrix{Int64}(I,d,d))
        garr = simplify(replace_with_arrays(g(-i,-j),repl))

        repl[field] = garr
        #repl[field.metric(m,-n)] = Matrix{Int64}(I,d,d)
        harr = factor.(replace_with_arrays(h,repl))
        dharr = factor.(replace_with_arrays(dh,repl))
        hharr = factor.(replace_with_arrays(hh,repl))
        Rarr = factor(replace_with_arrays(R,repl))

        testpt = Dict(x=>0.1,y=>-0.8,z=>0.05,α=>1)
        @test float(subs(harr[1,1,1],testpt)) ≈ float(subs(-12*x / (-1+x^2+y^2+z^2)^3,testpt))
        @test float(subs(dharr[1,1,1,1],testpt)) ≈ float(subs(12*(1+5*x^2 - y^2 - z^2) / (-1+x^2+y^2+z^2)^4,testpt))
        @test float(subs(Rarr,testpt)) ≈ -4

        h = get_christoffel(F,field,g)
        Riemann = get_riemann(F,h,field,g)
        Ric = Riemann(-i,-j,-k,-l)*field.metric(i,k)
        Ric = simp(Ric)

        R = Ric(-i,-j)*field.metric(j,i)
        R = simp(R)
        vars = symbols("x y z",real=true)
        x,y,z = vars
        d = length(vars)
        repl = Dict{Any,Any}(F(m) => vars, field.delta(-m,-n) => Matrix{Int64}(I,d,d))
        garr = simplify(replace_with_arrays(g(-i,-j),repl))

        repl[field] = garr
        #repl[field.metric(m,-n)] = Matrix{Int64}(I,d,d)
        Rarr = factor(replace_with_arrays(R,repl))
        @test float(subs(Rarr,testpt)) ≈ -4
    end
end


text:
type Solution
    x
    value
    attach::Nullable
    function Solution(; x=[], value=0, attach = Nullable())
        new(x, value, attach)
    end
end

function sol_equal(sol1, sol2)
    if sol1.value != 0 && sol2.value!=0
        if abs(sol1.value - sol2.value) > my_precision
            return false
        end
    end
    if length(sol1.x) != length(sol2.x)
        return false
    end
    for i = 1:length(sol1.x)
        if abs(sol1.x[i] - sol2.x[i]) > my_precision
            return false
        end
    end
    return true
end

# find minimum solution in iset
function find_min(iset)
    min = Inf
    min_index = 0
    res = Solution()
    i = 0
    for sol in iset
        i += 1
        if sol.value < min
            min = sol.value
            min_index = i
            res = sol
        end
    end
    return res, min_index
end

# find maximum solution in iset
function find_max(iset)
    max = -Inf
    max_index = 0
    res = Solution()
    i = 0
    for sol in iset
        i += 1
        if sol.value > max
            max = sol.value
            max_index = i
            res = sol
        end
    end
    return res, max_index
end

function sol_print(sol)
    zoolog("value: $(sol.value)")
    zoolog("x: $(sol.x)")
end

function sol_write(f, sol)
    write(f, "value=$(sol.value)\n")
    write(f, "x=$(sol.x)\n")
end

function write_population(filename, positive_data, negative_data)
	f = open(filename, "w")
	write(f, "########################################\n")
	write(f, "positive_data: \n")
	for sol in positive_data
		sol_write(f, sol)
	end
	write(f, "########################################\n")
	write(f, "negative_data: \n")
	for sol in negative_data
		sol_write(f, sol)
	end
	close(f)
end


text:
using Printf
using StringEncodings
mutable struct  CSVinfo
     fieldname ::String
     signalStringName::String
     keyFieldNames::String
     fieldNameDict::Dict
     Area2FiberLen
     # PosNames
     function CSVinfo(;fieldname::String ="定位界标1,定位界标2,xxxx主机IP,通道号,定位界标1光程,定位界标2光程,子防区编号,防区起点GPS经度坐标,防区起点GPS纬度坐标,防区终点GPS经度坐标,防区终点GPS纬度坐标,左右岸,地名信息,界标之间围栏长度,管理处,分局",
          signalStringName::String="定位界标1光程",
          keyFieldNames::String ="定位界标1,定位界标1光程,定位界标2,定位界标2光程"
           )
          Area2FiberLen=[ ]
          new(fieldname,signalStringName,keyFieldNames ,Dict(),Area2FiberLen)
     end

end

#
function initDict!(CSVinfo1::CSVinfo)
     CSVinfo1.fieldNameDict = Dict()
     ts = split(CSVinfo1.fieldname,"," )
     for t = 1: length( ts )
          CSVinfo1.fieldNameDict[ ts[t] ] = 0
     end
     # return fieldNameDict
end

function readinCSVFile1!(CSVinfo1::CSVinfo;csvfile="../xxxx分局防区信息.csv",coding1=enc"GBK")
     ret=false;
     initDict!(CSVinfo1 )
     fields= split( CSVinfo1.keyFieldNames,",")
     fieldlen =length( fields)
     # csvfileAbsname = csvfile

     csvfileAbsname =joinpath(pwd(), csvfile)
     fin =nothing
     # println( string(coding1));println()
     if string(coding1 ) =="UTF8"
          fin = open(csvfileAbsname,"r" )
     elseif string(coding1 ) =="UTF8-BOM"
          println( csvfileAbsname)
          fin = open(csvfileAbsname,"r" )
         # @info "here ",fin
          i1 = read(fin,UInt8); #@info "here ", i1
          i2= read(fin,UInt8);#@info "here ", i2
          i3 = read(fin,UInt8);#@info "here ", i3
          if i1 ==  UInt8(239) &&  i2 == UInt8(187) && i3 == UInt8(191)
               println( "$csvfile: UTF8-BOM file ")
          else
               return false;
               # break
          end
     else
          fin = open(csvfileAbsname, coding1, "r");
          println( "$csvfile:$coding1 file ")
          # break
     end
     try
               CSVinfo1.Area2FiberLen =[]
               for line in readlines(fin)
                    ts = split(line,"," )
                    if occursin(CSVinfo1.signalStringName , line) ## 首行
                         for t = 1: length( ts )
                              CSVinfo1.fieldNameDict[ ts[t] ] = t
                         end
                    else
                         a=[]

                         for i =1:fieldlen
                              field = fields[i]
                              idx = CSVinfo1.fieldNameDict[ field]
                              val  = ts[ idx ]
                              append!(a,[val])
                         end
                         append!(CSVinfo1.Area2FiberLen,[ a ])
                    end
               end
               # println("open and read ok !!!!!!!!!!!!!");println()
               ret=true;
     catch
          # println(e)
          # println("open and read Error !!!!!!!!!!!!!");println()
     finally
         # close(fin )
     end
     # println("open and read result: $ret ----");println()
     return ret
end

function readinCSVFile!(CSVinfo1::CSVinfo;csvfile="../xxxx分局防区信息.csv")
     ret= false;
     while ret ==false
          for bianma=[enc"GBK",enc"UTF-8"]
               println()
               ret = readinCSVFile1!(CSVinfo1,csvfile=csvfile,coding1=bianma)
               if ret == true
                    break
               end
          end
     end
end
function show(CSVinfo1::CSVinfo )
     fields= split( CSVinfo1.keyFieldNames,",")
     fieldlen =length( fields)
     for  row  in CSVinfo1.Area2FiberLen
          for i =1:fieldlen
               fieldname = fields[ i ]
               print("$fieldname:$(row[i]) ,")
          end
          println()
     end
end
function write2CSV(;CSVinfo1::CSVinfo,filename::String )
     fileobj = open( filename,"w")
     # EF BB BF
     write( fileobj,0xEF)
     write( fileobj,0xBB)
     write( fileobj,0xBF)

     write( fileobj, CSVinfo1.keyFieldNames )
     write(fileobj,"\r\n")

     # fields= split( CSVinfo1.keyFieldNames,",")
     # fieldlen =length( fields)
     total = length( CSVinfo1.Area2FiberLen)
     for  idx  = 1: total
          row = CSVinfo1.Area2FiberLen[idx]

          for i =1:length( row )-1

               write( fileobj, "$(row[i])," )
               # fieldname = fields[ i ]
               # print("$fieldname:$(row[i]) ,")
          end
          write( fileobj, "$(row[end])\r\n" )
          # println()
     end
end

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
export MultiPPOManager

mutable struct MultiPPOManager <: AbstractPolicy
    agents::Dict{Int, Agent}
    agents_inv::Dict{Agent, Vector{Int}}
    traces
end

Base.getindex(A::MultiPPOManager, x) = getindex(A.agents, x)
# used for simultaneous environments.
function (π::MultiPPOManager)(env::AbstractEnv)
    actions = Dict()
    for (agent::Agent, players) in π.agents_inv

        states = [state(env, p) for p in players]
        states = cat(states..., dims=ndims(states[1])+1)
        # no mask rn
        dists = prob(agent.policy, states, nothing)

        for (p, dist) in zip(players, dists)
            action = rand.(agent.policy.rng, dist)
            if ndims(action) == 2
                action_log_prob = sum(logpdf.(dist, action), dims = 1)
            else
                action_log_prob = logpdf.(dist, action)
            end

            actions[p] = EnrichedAction(action; action_log_prob = [action_log_prob])
        end
    end
    actions
end

function (π::MultiPPOManager)(stage::PreActStage, env::AbstractEnv, actions)
    # update each agent's trajectory.
    for (agent, players) in π.agents_inv
        update_trajectory!(agent, env, stage, actions, players)
        RLBase.update!(agent.policy, agent.trajectory, env, stage)
    end
end

function (π::MultiPPOManager)(stage::PostActStage, env::AbstractEnv)
    # only need to update trajectory.
    for (agent, players) in π.agents_inv
        update_trajectory!(agent, env, stage, players)
    end
end

function (π::MultiPPOManager)(stage::PreEpisodeStage, env::AbstractEnv)
    # collect state and a dummy action to each agent's trajectory here.
    for (agent, players) in π.agents_inv
        update_trajectory!(agent, env, stage, players)
    end
end

function (π::MultiPPOManager)(stage::PostEpisodeStage, env::AbstractEnv)
    # collect state and a dummy action to each agent's trajectory here.
    for (agent, players) in π.agents_inv
        update_trajectory!(agent, env, stage, players)
    end
end

function update_trajectory!(
    agent::Agent,
    env::AbstractEnv,
    ::PreEpisodeStage,
    players::Vector{Int}
)
    t = agent.trajectory
    if length(t) > 0
        pop!(t[:state])
        pop!(t[:action])
        pop!(t[:action_log_prob])
    end
end

function update_trajectory!(
    agent::Agent,
    env::AbstractEnv,
    ::PostEpisodeStage,
    players::Vector{Int}
)
    # Note that for trajectories like `CircularArraySARTTrajectory`, data are
    # stored in a SARSA format, which means we still need to generate a dummy
    # action at the end of an episode.

    state_dict      = state(env, players)
    states          = ecat([state_dict[p] for p in players]...)
    actions         = [RLCore.get_dummy_action(action_space(env, p)) for p in players]
    action_log_prob = [0.0 for p in players]

    push!(agent.trajectory;
          state           = states,
          action          = actions,
          action_log_prob = action_log_prob,)
end

function update_trajectory!(
    agent::Agent,
    env::AbstractEnv,
    ::PreActStage,
    actions::Dict,
    players::Vector{Int}
)
    state_dict = state(env, players)

    push!(agent.trajectory;
          state           = ecat([state_dict[p] for p in players]...),
          action          = ecat([actions[p].action for p in players]...),
          action_log_prob = ecat([actions[p].meta.action_log_prob[1] for p in players]...))
end

function update_trajectory!(
    agent::Agent,
    env::AbstractEnv,
    ::PostActStage,
    players::Vector{Int}
)
    reward_dict = reward(env, players)
    term_dict = is_terminated(env, players)
    push!(agent.trajectory[:reward], ecat([reward_dict[p] for p in players]...))
    push!(agent.trajectory[:terminal], ecat([term_dict[p]   for p in players]...))
end


text:
function read_input(input_filename)
    lines = readlines(input_filename)
    grid = zeros(Int, length(lines), length(lines[1]))
    for iy in 1:length(lines)
        for ix in 1:length(lines[iy])
            grid[iy, ix] = parse(Int, lines[iy][ix])
        end
    end
    return grid
end

function find_neighbours(coords, grid)
    neighbours = []
    for xdiff in -1:1
        for ydiff in -1:1
            if ((xdiff == 0) && (ydiff ==0) ||
                (((coords[1] + xdiff) in 1:size(grid)[1]) == false) ||
                (((coords[2] + ydiff) in 1:size(grid)[2]) == false))
                continue
            end
            push!(neighbours, CartesianIndex(coords[1]+xdiff, coords[2]+ydiff))
        end
    end
    return neighbours
end

function flash(coords, grid, flash_count)
    flash_count += 1
    grid[coords] = 0
    neighbours = find_neighbours(coords, grid)
    for neighbour in neighbours
        if grid[neighbour] == 0
            continue
        elseif grid[neighbour] == 9
            grid, flash_count = flash(neighbour, grid, flash_count)
        else
            grid[neighbour] += 1
        end
    end
    return grid, flash_count
end

function process_step(grid, flash_count)
    # initially increment every element by one
    grid = map(x->(x+1), grid)
    flasher_coords = findall(x->x>9, grid)
    for flasher in flasher_coords
        grid, flash_count = flash(flasher, grid, flash_count)
    end
    return grid, flash_count
end

function process_many_steps(grid)
    flash_count = 0
    i=0
    while true
        i += 1
        grid, flash_count = process_step(grid, flash_count)
        if (count(x->x==0, grid) == length(grid))
            return i
        end
    end
end

grid = read_input("input.txt")
all_flash_step = process_many_steps(grid)

println("answer is ", all_flash_step)