text:
using LightGraphs
using MetaGraphs
using CGP
using TikzGraphs
using TikzPictures

CGP.Config.init("cfg/test.yaml")

function example_graph()
    fg = MetaDiGraph(SimpleDiGraph())
    nin = 3; nout = 2
    add_vertices!(fg, nin + nout + 3)
    types = [0*ones(nin); 2*ones(3); 1*ones(nout)]
    for i in eachindex(types)
        set_prop!(fg, i, :type, types[i])
    end
    set_prop!(fg, 4, :function, CGP.Config.f_sum)
    set_prop!(fg, 4, :type, 2)
    set_prop!(fg, 4, :name, "fsum")
    add_edge!(fg, Edge(1, 4))
    set_prop!(fg, 1, 4, :ci, 1)
    add_edge!(fg, Edge(2, 4))
    set_prop!(fg, 2, 4, :ci, 2)
    set_prop!(fg, 5, :function, CGP.Config.f_mult)
    set_prop!(fg, 5, :name, "fmult")
    add_edge!(fg, Edge(3, 5))
    set_prop!(fg, 3, 5, :ci, 1)
    add_edge!(fg, Edge(4, 5))
    set_prop!(fg, 4, 5, :ci, 2)
    set_prop!(fg, 6, :function, CGP.Config.f_aminus)
    set_prop!(fg, 6, :name, "aminus")
    add_edge!(fg, Edge(1, 6))
    set_prop!(fg, 1, 6, :ci, 1)
    add_edge!(fg, Edge(3, 6))
    set_prop!(fg, 3, 6, :ci, 2)
    add_edge!(fg, Edge(5, 7))
    set_prop!(fg, 5, 7, :ci, 0)
    add_edge!(fg, Edge(6, 8))
    set_prop!(fg, 6, 8, :ci, 0)

    set_prop!(fg, 1, :name, "in1")
    set_prop!(fg, 2, :name, "in2")
    set_prop!(fg, 3, :name, "in3")
    set_prop!(fg, 7, :name, "out1")
    set_prop!(fg, 8, :name, "out1")

    set_prop!(fg, 1, :colour, "blue!0")
    set_prop!(fg, 2, :colour, "blue!100")
    set_prop!(fg, 3, :colour, "blue!70")
    set_prop!(fg, 4, :colour, "blue!60")
    set_prop!(fg, 5, :colour, "blue!40")
    set_prop!(fg, 6, :colour, "blue!20")
    set_prop!(fg, 7, :colour, "blue!10")
    set_prop!(fg, 8, :colour, "red!10")

    # max(1, min(Int64(ceil(100 * abs(-1.5) / 2)), 100))

    fg
end

function example_pdf()
    mg = example_graph()
    names = map(x->get_prop(mg, x, :name), 1:nv(mg))
    styles =  Dict(map(x-> x => "fill=$(get_prop(mg, x, :colour))", 1:nv(mg)))
    t = TikzGraphs.plot(mg.graph, names, node_style="draw, rounded corners", node_styles=styles)
    TikzPictures.save(TikzPictures.PDF("example.pdf"), t)
end


text:
cufftfloat(x) = _cufftfloat(float(x))
_cufftfloat(::Type{T}) where {T<:cufftReals} = T
_cufftfloat(::Type{Float16}) = Float32
_cufftfloat(::Type{Complex{T}}) where {T} = Complex{_cufftfloat(T)}
_cufftfloat(::Type{T}) where {T} = error("type $T not supported")
_cufftfloat(x::T) where {T} = _cufftfloat(T)(x)

complexfloat(x::CuArray{Complex{<:cufftReals}}) = x
realfloat(x::CuArray{<:cufftReals}) = x

complexfloat(x::CuArray{T}) where {T<:Complex} = copy1(typeof(cufftfloat(zero(T))), x)
complexfloat(x::CuArray{T}) where {T<:Real} = copy1(typeof(complex(cufftfloat(zero(T)))), x)

realfloat(x::CuArray{T}) where {T<:Real} = copy1(typeof(cufftfloat(zero(T))), x)

function copy1(::Type{T}, x) where T
    y = CuArray{T}(undef, map(length, axes(x)))
    #copy!(y, x)
    y .= broadcast(xi->convert(T,xi),x)
end

# promote to a complex floating-point type (out-of-place only),
# so implementations only need Complex{Float} methods
for f in (:fft, :bfft, :ifft)
    pf = Symbol("plan_", f)
    @eval begin
        $f(x::CuArray{<:Real}, region=1:ndims(x)) = $f(complexfloat(x), region)
        $pf(x::CuArray{<:Real}, region) = $pf(complexfloat(x), region)
        $f(x::CuArray{<:Complex{<:Union{Integer,Rational}}}, region=1:ndims(x)) = $f(complexfloat(x), region)
        $pf(x::CuArray{<:Complex{<:Union{Integer,Rational}}}, region) = $pf(complexfloat(x), region)
    end
end
rfft(x::CuArray{<:Union{Integer,Rational}}, region=1:ndims(x)) = rfft(realfloat(x), region)
plan_rfft(x::CuArray{<:Real}, region) = plan_rfft(realfloat(x), region)

*(p::Plan{T}, x::CuArray) where {T} = p * copy1(T, x)
*(p::ScaledPlan, x::CuArray) = rmul!(p.p * x, p.scale)

ext: jl


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Convert a satellite state vector from an Earth-Centered Inertial (ECI)
#   reference frame to another ECI reference frame.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export svECItoECI

"""
    svECItoECI(sv::SatelliteStateVector, args...)

Convert the satellite state vector `sv` from an ECI frame to another ECI frame.
The arguments `args...` must match those of the function `rECItoECI` **wihtout**
the rotation representation.

"""
function svECItoECI(sv::SatelliteStateVector, args...)
    D = rECItoECI(DCM, args...)

    # Since both frames does not have a significant angular velocity between
    # them, then we just need to convert the representations.
    return SatelliteStateVector(t = sv.t, r = D*sv.r, v = D*sv.v, a = D*sv.a)
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mobile_analytics

using AWS.Compat
using AWS.UUIDs
"""
    PutEvents()

The PutEvents operation records one or more events. You can have up to 1,500 unique custom events per app, any combination of up to 40 attributes and metrics per custom event, and any number of attribute or metric values.

# Required Parameters
- `events`: An array of Event JSON objects
- `x-amz-Client-Context`: The client context including the client ID, app title, app version and package name.

# Optional Parameters
- `x-amz-Client-Context-Encoding`: The encoding used for the client context.
"""

put_events(events, x_amz_Client_Context; aws_config::AWSConfig=global_aws_config()) = mobile_analytics("POST", "/2014-06-05/events", Dict{String, Any}("events"=>events, "headers"=>Dict{String, Any}("x-amz-Client-Context"=>x_amz_Client_Context)); aws_config=aws_config)
put_events(events, x_amz_Client_Context, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = mobile_analytics("POST", "/2014-06-05/events", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("events"=>events, "headers"=>Dict{String, Any}("x-amz-Client-Context"=>x_amz_Client_Context)), args)); aws_config=aws_config)


text:
using PersistentCohomology
using PersistentCohomology: Span

using Test

using SparseArrays, StructArrays
using LinearAlgebra
using GaloisFields

@testset "vietorisrips" begin

    m = [0.0 0.2 0.9 0.0
         0.2 0.0 0.5 0.0
         0.9 0.5 0.0 0.0
         0.0 0.0 0.0 0.0]

    s = sparse(m)

    cplx = vietorisrips(s, 2)
    @test keys(cplx[1]) == StructArray([(1,), (2,), (3,), (4,)],)
    @test values(cplx[1]) == [0.0, 0.0, 0.0, 0.0]

    @test keys(cplx[2]) == StructArray([(2, 1), (3, 1), (3, 2)])
    @test values(cplx[2]) == [0.2, 0.9, 0.5]

    @test keys(cplx[3]) == StructArray([(3, 2, 1)])
    @test values(cplx[3]) == [0.9]

    @test length(cplx) == 3

    cplx = vietorisrips(s, 1)
    @test !issparse(cplx[1])
    @test keys(cplx[1]) == StructArray([(1,), (2,), (3,), (4,)],)
    @test values(cplx[1]) == [0.0, 0.0, 0.0, 0.0]

    @test keys(cplx[2]) == StructArray([(2, 1), (3, 1), (3, 2)])
    @test values(cplx[2]) == [0.2, 0.9, 0.5]

    @test length(cplx) == 2
end

@testset "persistentcocycles" begin
    GF = @GaloisField 3
    ts = [-pi, -pi/2, 0, pi/2]
    pts = map(t -> [cos(t), sin(t)], ts)

    m = round.([norm(pt1-pt2) for pt1 in pts, pt2 in pts], digits = 4)
    cplx = vietorisrips(sparse(m), 2)
    @test keys(cplx[1]) == StructArray(([1, 2, 3, 4],))
    @test keys(cplx[2]) == StructArray([(2, 1), (3, 1), (3, 2), (4, 1), (4, 2), (4, 3)])
    @test keys(cplx[3]) == StructArray([(3, 2, 1), (4, 2, 1), (4, 3, 1), (4, 3, 2)])

    short = m[1, 2]
    long = m[1, 3]

    cocycles, spans = persistentcocycles(GF, cplx, 1)

    @test spans[1][1] == Span{Float64}(0, Inf)
    @test findnz(values(cocycles[1][1])) == ([1, 2, 3, 4], GF.([1, 1, 1, 1]))
    @test spans[1][2] == Span{Float64}(0, short)
    @test findnz(values(cocycles[1][2])) == ([2], GF.([1]))
    @test spans[1][3] == Span{Float64}(0, short)
    @test findnz(values(cocycles[1][3])) == ([3], GF.([1]))
    @test spans[1][4] == Span{Float64}(0, short)
    @test findnz(values(cocycles[1][4])) == ([4], GF.([1]))
    @test length(cocycles[1]) == 4

    @test spans[2][1] == Span{Float64}(short, long)
    @test findnz(values(cocycles[2][1])) == ([6], GF.([1]))
    @test issparse(cocycles[2][1])
    @test findnz(cocycles[2][1]) == (StructVector([(4, 3)]), GF.([1]))
    @test length(cocycles[2]) == 1
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -12.331138821460632 <= q <= 11.287379993959185)

                     add_NL_constraint(m, :(sigmoid(0.16893740747119512 + 0.361142465903193*sigmoid(0.3759145594542668 + 0.23452173305395663*sigmoid(-0.8645963938709578 + -0.1467408759040536*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*sigmoid(-0.9609056399980664 + 0.06060702618594416*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*sigmoid(-0.9973666264873402 + 0.1126936245130743*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.21914525995597334*sigmoid(-0.017859443401365915 + -0.7776638835685663*sigmoid(-0.8645963938709578 + -0.1467408759040536*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*sigmoid(-0.9609056399980664 + 0.06060702618594416*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*sigmoid(-0.9973666264873402 + 0.1126936245130743*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.7287835676612944*sigmoid(0.6615910141517634 + 0.18731902906517428*sigmoid(-0.8645963938709578 + -0.1467408759040536*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*sigmoid(-0.9609056399980664 + 0.06060702618594416*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*sigmoid(-0.9973666264873402 + 0.1126936245130743*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + sigmoid(-0.9123133183860701 + 0.273513749523993*sigmoid(0.3759145594542668 + 0.23452173305395663*sigmoid(-0.8645963938709578 + -0.1467408759040536*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*sigmoid(-0.9609056399980664 + 0.06060702618594416*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*sigmoid(-0.9973666264873402 + 0.1126936245130743*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.052599029592083735*sigmoid(-0.017859443401365915 + -0.7776638835685663*sigmoid(-0.8645963938709578 + -0.1467408759040536*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*sigmoid(-0.9609056399980664 + 0.06060702618594416*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*sigmoid(-0.9973666264873402 + 0.1126936245130743*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.2642412570560455*sigmoid(0.6615910141517634 + 0.18731902906517428*sigmoid(-0.8645963938709578 + -0.1467408759040536*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*sigmoid(-0.9609056399980664 + 0.06060702618594416*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*sigmoid(-0.9973666264873402 + 0.1126936245130743*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + sigmoid(0.5514504616591211 + 0.27732859795053555*sigmoid(0.3759145594542668 + 0.23452173305395663*sigmoid(-0.8645963938709578 + -0.1467408759040536*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*sigmoid(-0.9609056399980664 + 0.06060702618594416*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*sigmoid(-0.9973666264873402 + 0.1126936245130743*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + -0.08702067567597149*sigmoid(-0.017859443401365915 + -0.7776638835685663*sigmoid(-0.8645963938709578 + -0.1467408759040536*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*sigmoid(-0.9609056399980664 + 0.06060702618594416*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*sigmoid(-0.9973666264873402 + 0.1126936245130743*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.40713382739201665*sigmoid(0.6615910141517634 + 0.18731902906517428*sigmoid(-0.8645963938709578 + -0.1467408759040536*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*sigmoid(-0.9609056399980664 + 0.06060702618594416*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*sigmoid(-0.9973666264873402 + 0.1126936245130743*sigmoid(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*sigmoid(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*sigmoid(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function CRT_SIA(m, Nx, ϕ, u, D, K_eq, nu_A, nu_B, nu_C, C_bc, C_0, N_comp, dt, T, tol)
  # boundary conditions
  BC = Array{Any}(N_comp)
  for i in 1:N_comp
      BC[i] = createBC(m)
      # Dirichlet BC on the left side
      BC[i].left.a[:] = 0.0
      BC[i].left.b[:] = 1.0
      BC[i].left.c[:] = C_bc[i]
  end

  # initial condition
  C_old = Array{Any}(N_comp)
  for i in 1:N_comp
      C_old[i] = createCellVariable(m, C_0[i], BC[i]) # initial value of concentration
  end

  # Advective term
  u_vec = createFaceVariable(m, u)
  M_adv = convectionUpwindTerm(u_vec) # matrix of coef for the advection term

  # Diffusive term
  M_diff = Array{Any}(N_comp)
  for i in 1:N_comp
    D_face = createFaceVariable(m, D[i]*ϕ)
    M_diff[i] = diffusionTerm(D_face)      # matrix of coef for the diffusion term
  end

  # Boundary terms
  M_bc = Array{Any}(N_comp)
  RHS_bc  = Array{Any}(N_comp)
  for i in 1:N_comp
      M_bc[i], RHS_bc[i] = boundaryConditionTerm(BC[i]) # matrix and RHS for the BC
  end

  # Constructing stationary part of each system matrix
  M_stat = Array{Any}(N_comp)
  for i in 1:N_comp
    M_stat[i] = M_adv - M_diff[i] + M_bc[i]
  end

  # Preparing transient terms
  M_t = Array{Any}(N_comp)
  RHS_t = Array{Any}(N_comp)

  maxit = 5000

  # Time-stepping
  for t in dt:dt:T
      # transport equation
      for i in 1:N_comp
          M_t[i], RHS_t[i] = transientTerm(C_old[i], dt, ϕ)
          M = M_t[i]+M_stat[i]
          RHS = RHS_t[i]+RHS_bc[i]
          C_new = solveLinearPDE(m, M, RHS)
          C_old[i] = copyCell(C_new)
      end

      # Iterative loop
      iter = 0
      error = 1000
      eps0 = 0
      eps1 = 0

      while (error > tol && iter < maxit)
        # reaction equation
        eps1 = epsilon.(internalCells(C_old[1]), internalCells(C_old[2]),
                internalCells(C_old[3]), nu_A, nu_B, nu_C, K_eq)
        C_old[1] = createCellVariable(m, internalCells(C_old[1])-nu_A*eps1, BC[1])
        C_old[2] = createCellVariable(m, internalCells(C_old[2])-nu_B*eps1, BC[2])
        C_old[3] = createCellVariable(m, internalCells(C_old[3])+nu_C*eps1, BC[3])

        error = sum(abs(eps1-eps0))
        eps0 = eps1
        iter += 1
    end

    #print(iter, ", ", t, ", ", error, "\n")
    #print(sum(internalCells(C_old[3])-K_eq*internalCells(C_old[1]).*internalCells(C_old[2])), "\n") #mass balance
  end

  return C_old
end

max_stars_repo_name: FrederikLehn/coupled-solvers

max_issues_repo_name: FrederikLehn/coupled-solvers

max_forks_repo_name: FrederikLehn/coupled-solvers


text:
function Base.copy(fut::AbstractFuture)
    res = Future(fut)

    # partitioned_using() do
    #     keep_all_sample_keys(res, fut)
    #     keep_sample_rate(res, fut)
    # end

    # TODO: Ensure that the CopyTo function can always be used to split anything
    # for which there is no specified PT. But make sure that other splitting
    # functions such as those that split arrays are expecting a name to be
    # provided. If we do that, we can un-comment the below and allow for
    # copying anything
    # pt(fut, PartitionType())
    # pt(res, PartitionType(), match=fut)
    pt(fut, Replicating())
    pt(res, Replicating(), match=fut)

    @partitioned fut res begin
        res = Base.copy(fut)
    end

    res
end

function Base.deepcopy(fut::AbstractFuture)
    res = Future(fut)

    # partitioned_using() do
    #     keep_all_sample_keys(res, fut)
    #     keep_sample_rate(res, fut)
    # end

    # TODO: Ensure that the CopyTo function can always be used to split anything
    # for which there is no specified PT. But make sure that other splitting
    # functions such as those that split arrays are expecting a name to be
    # provided. If we do that, we can un-comment the below and allow for
    # copying anything
    # pt(fut, PartitionType())
    # pt(res, PartitionType(), match=fut)
    pt(fut, Replicating())
    pt(res, Replicating(), match=fut)

    @partitioned fut res begin
        res = Base.deepcopy(fut)
    end

    res
end

# Array type

struct Array{T,N} <: AbstractFuture where {T,N}
    data::Future
    size::Future
    # TODO: Add offset for indexing
    # offset::Future

    Array{T,N}() where {T,N} = new(Future(), Future())
    Array{T,N}(A::Array{T,N}) where {T,N} = new(Future(), Future(A.size))
    Array{T,N}(data::Future, size::Future) where {T,N} = new(data, size)
end

const Vector{T} = Array{T,1}
const Matrix{T} = Array{T,2}

Banyan.convert(::Type{Future}, A::Array{T,N}) where {T,N} = A.data

# Array sample

Banyan.sample_axes(A::U) where U <: Base.AbstractArray{T,N} where {T,N} = [1:ndims(A)...]
Banyan.sample_keys(A::U) where U <: Base.AbstractArray{T,N} where {T,N} = sample_axes(A)

# `sample_divisions`, `sample_percentile`, and `sample_max_ngroups` should
# work with the `orderinghash` of values in the data they are used on

function Banyan.sample_divisions(A::U, key) where U <: Base.AbstractArray{T,N} where {T,N}
    if isempty(A)
        return []
    end

    max_ngroups = sample_max_ngroups(A, key)
    ngroups = min(max_ngroups, 512)
    data = sort([orderinghash(e) for e in eachslice(A, dims=key)])
    datalength = length(data)
    grouplength = div(datalength, ngroups)
    # We use `unique` here because if the divisions have duplicates, this could
    # result in different partitions getting the same divisions.
    # TODO: Ensure that `unique` doesn't change the order
    unique([
        # Each group has elements that are >= start and < end
        (
            data[(i-1)*grouplength + 1],
            data[i == ngroups ? datalength : i*grouplength + 1]
        )
        for i in 1:ngroups
    ])
end

function Banyan.sample_percentile(A::U, key, minvalue, maxvalue) where U <: Base.AbstractArray{T,N} where {T,N}
    if isempty(A) || isnothing(minvalue) || isnothing(maxvalue)
        return 0
    end

    count((begin oh = orderinghash(e); oh >= minvalue && oh <= maxvalue end for e in eachslice(A, dims=key))) / size(A, key)

    # TODO: Determine whether we need to assume a more coarse-grained percentile using the divisions
    # from `sample_divisions` and computing the percent of divisions that overlap with the range

    # # minvalue, maxvalue = orderinghash(minvalue), orderinghash(maxvalue)
    # divisions = sample_divisions(A, key)
    # percentile = 0
    # divpercentile = 1/length(divisions)
    # inminmax = false

    # # Iterate through divisions to compute percentile
    # for (i, (divminvalue, divmaxvalue)) in enumerate(divisions)
    #     # Check if we are between the minvalue and maxvalue
    #     if (i == 1 || minvalue >= divminvalue) && (i == length(divisions) || minvalue < divmaxvalue)
    #         inminmax = true
    #     end

    #     # Add to percentile
    #     if inminmax
    #         percentile += divpercentile
    #     end

    #     # Check if we are no longer between the minvalue and maxvalue
    #     if (i == 1 || maxvalue >= divminvalue) && (i == length(divisions) || maxvalue < divmaxvalue)
    #         inminmax = false
    #     end
    # end

    # percentile
end

# NOTE: The key used for arrays must be a single dimension - if you are doing
# something like sorting on multiple dimensions or grouping on multiple
# dimensions you can just use the first dimension as the key.

function Banyan.sample_max_ngroups(A::U, key) where U <: Base.AbstractArray{T,N} where {T,N}
    if isempty(A)
        return 0
    end

    data = sort([orderinghash(e) for e in eachslice(A, dims=key)])
    currgroupsize = 1
    maxgroupsize = 0
    prev = nothing
    prev_is_nothing = true # in case `prev` _can_ be nothing
    for curr in data
        if !prev_is_nothing && curr == prev
            currgroupsize += 1
        else
            maxgroupsize = max(maxgroupsize, currgroupsize)
            currgroupsize = 1
        end
        # TODO: Maybe use deepcopy here if eltype might be nested
        prev = copy(curr)
        prev_is_nothing = false
    end
    maxgroupsize = max(maxgroupsize, currgroupsize)
    div(size(A, key), maxgroupsize)
end
# TODO: Handle issue where mapslices requires dims to be a single dimension;
# probably need to vary mapslices on the dim itself and then use eachslices,
# get the orderinghash and then take minimum across that
# TODO: Change to use eachslice everywhere and ensure we use key not d
Banyan.sample_min(A::U, key) where U <: Base.AbstractArray{T,N} where {T,N} = isempty(A) ? nothing : minimum((orderinghash(e) for e in eachslice(A, dims=key)))
Banyan.sample_max(A::U, key) where U <: Base.AbstractArray{T,N} where {T,N} = isempty(A) ? nothing : maximum((orderinghash(e) for e in eachslice(A, dims=key)))

# Array creation

function read_hdf5(path; kwargs...)
    A_loc = RemoteSource(path; kwargs...)
    A_loc.src_name == "Remote" || error("$path does not exist")
    if is_debug_on()
        # @show A_loc.src_parameters
        # @show A_loc.size
    end
    A = Future(source=A_loc)
    Array{A_loc.eltype,A_loc.ndims}(A, Future(A_loc.size))
end

function write_hdf5(A, path; invalidate_source=true, invalidate_sample=true, kwargs...)
    # # A_loc = Remote(pathname, mount)
    # destined(A, Remote(path, delete_from_cache=true))
    # mutated(A)
    # # This doesn't rely on any sample properties so we don't need to wrap this
    # # in a `partitioned_with` to delay the PT construction to after sample
    # # properties are computed.
    # pt(A, Blocked(A) | Replicated())
    # # partition(A, Replicated()) # TODO: Use semicolon for keyword args
    # # Distributed will allow for balanced=true|false and any divisions or key
    # # but the PT library should be set up such that you can't split if
    # # divisions or key are not provided or balanced=false
    # # partition(A, Blocked())
    # # for axis in 1:min(4, ndims(A))
    # #     # Partition with distribution of either balanced, grouped, or unknown
    # #     partition(A, Blocked(key=a), mutated=true)
    # # end
    # @partitioned A begin end
    # compute(A)
    pt(A, Blocked(A) | Replicated())
    partitioned_computation(
        A,
        destination=RemoteDestination(path; invalidate_source=invalidate_source, invalidate_sample=invalidate_sample, kwargs...),
        new_source=_->RemoteSource(path)
    )
end

function Banyan.write_to_disk(A::Array{T,N}) where {T,N}
    pt(A, Blocked(A) | Replicated())
    partitioned_computation(A, destination=Disk())
end

function fill(v, dims::NTuple{N,Integer}) where {N}
    fillingdims = Future(source=Size(dims))
    A = Array{typeof(v),N}(Future(), Future(dims))
    v = Future(v)
    dims = Future(dims)

    # For futures that contains dims or nrows, we initialize them to store the
    # dimensions or # of rows of the whole value. Then if we need to mutate,
    # we reassign to the dimensions or # of rows of the part and assign a
    # reducing PT.

    # NOTE: If a value is being created for the first time in a code region, it
    # is being mutated. The only exception is a value that is already created
    # on the client and is merely being replicated in the code region.
    # partition(A, Replicated())
    # partition(dims, Replicated())

    # for axis in 1:min(4, ndims(A))
    #     # Partition with distribution of either balanced, grouped, or unknown
    #     partition(A, Blocked(key=a))
    # end
    # We use `partitioned_with` here to ensure that a sample of A is produced
    # first so that we can use the Blocked PT constructor which depends on A
    # having its sample taken
    partitioned_with() do
        # blocked
        # TODO: Ensure that we are properly creating new PAs
        pt(A, Blocked(A))
        pt(fillingdims, Divided(), match=A, on="key")

        # replicated
        pt(A, fillingdims, v, Replicated())
    end

    # TODO: Remove the println @macroexpand
    # println(@macroexpand begin @partitioned A v fillingdims begin
    #     println(A)
    #     println(v)
    #     A = fill(v, fillingdims)
    # end end)
    if is_debug_on()
        @partitioned A v fillingdims begin
            A = Base.fill(v, fillingdims)
        end
    else
        @partitioned A v fillingdims begin
            A = Base.fill(v, fillingdims)
        end
    end

    A
end

fill(v, dims::Integer...) = fill(v, Tuple(dims))

zeros(::Type{T}, args...; kwargs...) where {T} = fill(zero(T), args...; kwargs...)
zeros(args...; kwargs...) where {T} = zeros(Float64, args...; kwargs...)
ones(::Type{T}, args...; kwargs...) where {T} = fill(one(T), args...; kwargs...)
ones(args...; kwargs...) where {T} = ones(Float64, args...; kwargs...)
trues(args...; kwargs...) where {T} = fill(true, args...; kwargs...)
falses(args...; kwargs...) where {T} = fill(false, args...; kwargs...)

# Array properties

Base.ndims(A::Array{T,N}) where {T,N} = ndims(sample(A))
Base.size(A::Array{T,N}) where {T,N} = collect(A.size)
Base.length(V::Array{T,N}) where {T,N} = prod(collect(V.size))
Base.eltype(A::Array{T,N}) where {T,N} = eltype(sample(A))

function pts_for_copying(A, res)
    # balanced
    pt(A, Blocked(A, balanced=true))
    pt(res, Balanced(), match=A)

    # unbalanced
    pt(A, Blocked(A, balanced=false, scaled_by_same_as=res))
    pt(res, Unbalanced(scaled_by_same_as=A), match=A)

    # replicated
    pt(A, res, Replicated())
end

function Base.copy(A::Array{T,N})::Array{T,N} where {T,N}
    res = Future()

    partitioned_using() do
        keep_all_sample_keys(res, A)
        keep_sample_rate(res, A)
    end

    partitioned_with() do
        pts_for_copying(A, res)
    end

    @partitioned A res begin
        res = Base.copy(A)
    end

    Array{T,N}(res, deepcopy(A.size))
end

function Base.deepcopy(A::Array{T,N})::Array{T,N} where {T,N}
    res = Future()

    partitioned_using() do
        keep_all_sample_keys(res, A)
        keep_sample_rate(res, A)
    end

    partitioned_with() do
        pts_for_copying(A, res)
    end

    @partitioned A res begin
        res = Base.deepcopy(A)
    end

    Array{T,N}(res, deepcopy(A.size))
end

# Array operations

function Base.map(f, c::Array{T,N}...) where {T,N}
    f = Future(f)
    res = Future()

    partitioned_using() do
        # We shouldn't need to keep sample keys since we are only allowing data
        # to be blocked for now. The sample rate is kept because it might be
        # smaller if this is a column of the result of a join operation.
        # keep_all_sample_keys(res, fut)
        keep_sample_rate(res, first(c))
    end

    # TODO: Determine whether array operations need to use mutated_from or mutated_to

    partitioned_with() do
        # balanced
        pt(first(c), Blocked(first(c), balanced=true))
        pt(c[2:end]..., res, Blocked() & Balanced(), match=first(c), on=["key", "id"])

        # unbalanced
        pt(first(c), Blocked(first(c), balanced=false, scaled_by_same_as=res))
        pt(c[2:end]..., res, Unbalanced(scaled_by_same_as=first(c)), match=first(c))

        # replicated
        if !is_debug_on()
            pt(c..., res, f, Replicated())
        else
            pt(f, Replicated())
        end
    end

    # println(@macroexpand begin @partitioned f c res begin
    #     res = Base.map(f, c...)
    # end end)
    @partitioned f c res begin
        res = Base.map(f, c...)
        # @show res
        # @show typeof(res)
        # @show eltype(res)
    end

    # @show sample(res)
    # @show typeof(sample(res))
    # @show eltype(sample(res))

    Array{eltype(sample(res)),N}(res, deepcopy(first(c).size))
end

# NOTE: This function is shared between the client library and the PT library
function indexapply(op, objs...; index::Integer=1)
    lists = [obj for obj in objs if (obj isa AbstractVector || obj isa Tuple)]
    length(lists) > 0 || throw(ArgumentError("Expected at least one tuple as input"))
    index = index isa Colon ? length(first(lists)) : index
    operands = [((obj isa AbstractVector || obj isa Tuple) ? obj[index] : obj) for obj in objs]
    indexres = op(operands...)
    res = first(lists)
    if first(lists) isa Tuple
        res = [res...]
        res[index] = indexres
        Tuple(res)
    else
        res = copy(res)
        res[index] = indexres
        res
    end
end

function Base.mapslices(f, A::Array{T,N}; dims) where {T,N}
    if isempty(dims) return map(f, A) end

    f = Future(f)
    res_size = Future()
    # TODO: Ensure that this usage of Any is correct here and elsewhere
    res = Array{Any,Any}(Future(), res_size)
    dims = Future(dims)

    partitioned_using() do
        keep_sample_rate(res, A)
    end

    partitioned_with() do
        # Blocked PTs along dimensions _not_ being mapped along
        bpt = [bpt for bpt in Blocked(A) if !(dims isa Colon) && !(bpt.key in [dims...])]

        # balanced
        pt(A, bpt & Balanced())
        pt(res, Blocked() & Balanced(), match=A, on="key")

        # unbalanced
        pt(A, bpt & Unbalanced(scaled_by_same_as=res))
        pt(res, Unbalanced(scaled_by_same_as=A), match=A)

        # replicated
        # TODO: Determine why this MatchOn constraint is not propagating
        pt(res_size, ReducingWithKey(quote axis -> (a, b) -> indexapply(+, a, b, index=axis) end), match=A, on="key")
        pt(A, res, res_size, f, dims, Replicated())
    end

    @partitioned f A dims res begin
        res = Base.mapslices(f, A, dims=dims)
        res_size = Base.size(res)
    end

    res
end

# TODO: Implement reduce and sortslices

function Base.reduce(op, A::Array{T,N}; dims=:, kwargs...) where {T,N}
    if :init in keys(kwargs) throw(ArgumentError("Reducing with an initial value is not currently supported")) end

    op = Future(op)
    res_size = Future()
    res = dims isa Colon ? Future() : Array{Any,Any}(Future(), res_size)
    if is_debug_on()
        # @show dims # TODO: Ensure this isn't a function
    end
    dims = Future(dims)
    kwargs = Future(kwargs)

    partitioned_using() do
        keep_sample_rate(res, A)
    end

    partitioned_with() do
        # TODO: Duplicate annotations to handle the balanced and unbalanced cases
        # seperately
        # TODO: Have a better API where duplicating to handle balanced and unbalanced
        # isn't needed
        # TODO: Cascaade MatchOn constraints
        # TODO: Ensure that MatchOn value is being discovered

        for bpt in Blocked(A)
            pt(A, bpt)
            if collect(dims) isa Colon || bpt.key in [collect(dims)...]
                # NOTE: Be careful about trying to serialize things that would
                # require serializing the whole Banyan module. For example, if
                # this where Reducing(op) or if we tried Future(op) where op
                # could refer to a + function overloaded by BanyanArrays.
                pt(res, Reducing(collect(op)))
            else
                pt(res, bpt.balanced ? Balanced() : Unbalanced(scaled_by_same_as=A), match=A)
            end
        end
        pt(res_size, ReducingWithKey(quote axis -> (a, b) -> indexapply(+, a, b, index=axis) end), match=A, on="key")
        # TODO: Allow replication
        if !is_debug_on()
            pt(A, res, res_size, dims, kwargs, op, Replicated())
        else
            pt(dims, kwargs, op, Replicated())
        end
    end

    @partitioned op A dims kwargs res res_size begin
        if is_debug_on()
            # @show size(A)
            # @show dims # TODO: Figure out why dims is sometimes a function
        end
        res = Base.reduce(op, A; dims=dims, kwargs...)
        if res isa Base.Array
            res_size = Base.size(res)
        end
    end

    res
end

function Base.sortslices(A::Array{T,N}, dims; kwargs...) where {T,N}
    get(kwargs, :by, identity) == identity || throw(ArgumentError("Sorting by a function is not supported"))
    !haskey(kwargs, :order) || throw(ArgumentError("Sorting by an order is not supported"))

    # Determine what to sort by and whether to sort in reverse
    sortingdim = dims isa Colon ? 1 : first(dims)
    isreversed = get(kwargs, :rev, false)

    res = Future()
    dims = Future(dims)
    kwargs = Future(kwargs)

    partitioned_using() do
        keep_sample_rate(res, A)
    end

    partitioned_with() do
        # unbalanced -> unbalanced
        pt(A, Grouped(A, by=sortingdim, rev=isreversed, scaled_by_same_as=res, balanced=false))
        pt(res, Blocked() & Unbalanced(scaled_by_same_as=A), match=A, on=["key", "divisions", "id"])

        # balanced -> balanced
        pt(A, Grouped(A, by=sortingdim, rev=isreversed, balanced=true))
        pt(res, Blocked() & Balanced(), match=A, on=["key", "divisions", "id"])

        # replicated
        pt(A, res, dims, kwargs, Replicated())
    end

    @partitioned A dims kwargs res begin
        res = Base.sortslices(A, dims=dims, kwargs...)
    end

    Array{T,N}(res, deepcopy(A.size))
end

Base.sort(A::Array{T,N}; kwargs...) where {T,N} = sortslices(A, dims=:; kwargs...)

# Array aggregation

# TODO: Determine split between what happens in annotations vs PT constructors
# - Annotations
#   - Collecting axes/keys that are relevant to the computation by looking at arguments
#   - Iterating through keys that require special partitioning
#   - Registering a delayed PA and passing in job with info
# - PT Library
#   - Maintaining job_properties with keys and axes that are used

# for (op, agg) in [(:(sum), :(Base.:+)), (:(minimum), :(min)), (:(maximum), :(max))]
#     @eval begin
#         function $op(A::Array{T, N})
#             dims = dims isa Tuple ? dims : tuple(dims)
#             dimensions = Future(dims)
#             operator = Future($op)
#             res = Future()
#         end

#         function $op(A::Array{T, N}; dims=:) where {T, N}
#             dims = dims isa Tuple ? dims : tuple(dims)
#             dimensions = Future(dims)
#             operator = Future($op)
#             res_size = Future() # TODO: Compute this here using Future, mutation=
#             res = Array{T, N}(Future(), res_size)

#             partition(A, Replicated())
#             partition(res, Replicated())

#             # Partition by distributing across 1 of up to 4 axes
#             for axis in 1:min(4, ndims(A))
#                 # Partition with distribution of either balanced, grouped, or unknown
#                 partition(A, Blocked(key=axis))

#                 # Partition the result based on whether the reduction is across the dimension of partitioning
#                 if axis in dims
#                     partition(res, Reducing(reducer=$agg))
#                 else
#                     partition(res, Partitioned(), matches_with=A)
#                 end
#             end

#             # Partition all as replicated
#             partition(operator, Replicated())
#             partition(res_size, Replicated())

#             @partitioned A dimensions operator aggregaton res res_size begin
#                 res = $op(V, dims=dimensions)
#                 res_size = size(res)
#             end
#         end
#     end
# end

# function mean(V::Vector)
#     # TODO: Accept vector partitioned in different dimensions and return
#     # replicated with reducing or distributed with same id potentially grouped
#     # by an axis or balanced
#     sum(V) / length(V)
# end

# # Array element-wise operations

# # TODO: Implement a bunch of operations where we allow Replicated and Blocked
# # across a dimension where communication is not needed (e.g., accumulation
# # across dims=2 allows us to have Blocked along=1)

# Array unary operations

for op in [:-]
    @eval begin
        Base.$op(X::Array{T,N}) where {T,N} = map($op, X)
    end
end

for (op, agg) in [(:(sum), :(+)), (:(minimum), :(min)), (:(maximum), :(max))]
    # TODO: Maybe try ensuring that the Base.:+ here is not including the method from above
    @eval begin
        Base.$op(X::Array{T,N}; dims=:) where {T,N} = reduce($agg, X; dims=dims)
    end
end

# Array binary operations

# for op in [:+, :-, :>, :<, :(>=), :(<=), :(==), :!=]
# NOTE: Only addition and subtraction are supported
for op in [:+, :-]
    @eval begin
        Base.$op(A::Array{T,N}, B::Array{T,N}) where {T,N} = map($op, A, B)
    end
end

# TODO: Add broadcasting support

# # Binary operators
# for op in (:+, :-, :>, :<, :>=, :<=, :≥, :≤, :(==), :!=)
#     @eval begin
#         function Base.$op(A::Array{T, N}, B::Array{T, N}) where {T, N}
#             A_size = A.size
#             res_size = Future(A_size)
#             res = Array{T, N}(Future(), res_size)
#             op = Future(op)
#             # TODO:
#             # - make samples get computed instantly
#             # TODO: Put info in partitioned:
#             # - Future
#             # - PT (replicated, balanced, grouped, pseudo)
#             # - PCs (matching on, at most)
#             # - location (None) with sample, sample properties, total memory
#             # - mutating (constant)
#             # TODO: for columns: how to partition on grouped
#             # - maintain set of possible keys to group by
#             # - compute atmost whenever needed
#             # TODO: Implement new Distributed to support group-by views - how to partition on unknown?
#             # - distribution = balanced by axis|grouped by axis or key with divisions|unknown
#             # - identified by id = true|false
#             # - axis
#             # - key
#             # - divisions - ensure that a split function can't be used if it is required but not provided
#             # - ID (must be null for splitting or casting to)
#             # TODO: Auto-generate default random on casts and splits
#             # "!" for parameters that are to be randomly assigned, can be used both in PAs and in PT library

#             # TODO: Accept replicated or balanced or distributed with same ID for inputs and outputs while grouped by any valid key

#             partition(A, Replicated())
#             for axis in 1:min(4, ndims(A))
#                 partition(A, Blocked(key=axis))
#             end
#             partition(B, Partitioned(), matches_with=A)
#             partition(res, Partitioned(), matches_with=A, mutating=true)
#             partition(res_size, ReplicatedOrReducing(), match=A_size)
#             partition(op, Replicated())

#             # colnames = propertynames(sample(df))
#             # partition() do job
#             #     job.properties[:keys]
#             #     Grouped(propertynames(sample(df)))
#             #     Grouped(df, groupingkey, reverse=false)
#             # end
#             # partition(A, Balanced())
#             # partition(A, Distributed(like=B))
#             # partition(A_size, Replicated())
#             # partition(B, Distributed(), matches_with=A)
#             # partition(res, Distributed(), matches_with=A)
#             # partition(res_nrows, Replicated())
#             # partition(op, Replicated())

#             @partitioned res res_size A A_size B op begin
#                 res = op(A, B)
#                 res_size = A_size
#             end

#             res
#         end
#     end
# end

# # TODO: Implement unary operators

# # TODO: Implement element-wise operations

# # Array reshaping

# # TODO: Support sorting once we have a non-verbose way of specifying grouping for both dataframes and arrays

# # function sort(A::Array; kwargs...)
# #     # TODO: Accept replicated or grouped on key used for sotring and return same but with same ID but
# #     # newly generated ID if the input has id set to null because of casting or re-splitting
# #     # TODO: Accept A as either blocked on some other dimension or grouped on the dimension
# #     # used for sorting. Then have the result be blocked. And then use Blocked everywhere else and
# #     # ensure that getindex will only return Vector that is Blocked and not Grouped
# #     res = Future()
# #     partition(A, Grouped())
# # end

# function sortslices(A::Array)
#     # TODO: use the first of the dims to sort on to split the array
#     # TODO: use the first of the slice when splitting the array on the dim to get quantiles
# end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
using LinearAlgebra, PolynomialRoots
include("findlocalmaxima.jl")

function rootMUSIC(Y,K,SNR_db=20,m=0)
    # Estimation of signal parameters via rotational invariance techniques
    # (ESPRIT) algorithm for line spectral estimation from measurements.
    # Can handle multiple snapshots by constructing a low rank Hankel matrix
    # from the data.
    #
    # Inputs:
    # Y: Data matrix
    # K: number of sinusoids
    # SNR_db: SNR in each dimension of the data in dB
    # m: Hankel parameter of the Hankel matrix of the data
    #
    # Parameters:
    # N = size(Y,1): number of measurements per snapshot
    # L = size(Y,2): number of snapshots
    #
    # Output:
    # f_hat: estimate of the frequencies of the K sinusoids in [rad]

(N,L) = size(Y)

if m == 0
m = ceil(Int,N/2)
end

R = zeros(m,m)
for j = 1:L, i = m:N
    R = R + Y[i-m+1:i,j]*Y[i-m+1:i,j]'/(L*N)
end

    evdR = eigen(R)
    p = sortperm(real(evdR.values))[1:m-K]
    Un = evdR.vectors[:,p]

    M = Un*Un'

    # Construct polynomial
    polcoeffs = Array{Complex{Float64},1}(undef,2*m-1)
    for k = 1:2*m-1
        polcoeffs[k] = sum(diag(M,k-m))
    end

    # Find roots of polynomial
    polyroots = roots(polcoeffs)

    # Remove roots with multiplicity > 1
    idcs = []
    for i = 1:length(polyroots)-1, j = i+1:length(polyroots)
        if abs(polyroots[i]-polyroots[j]) < 0.5/SNR_db + 0.001
            push!(idcs,i)
        end
    end

    polyroots2 = polyroots[setdiff(1:length(polyroots),idcs)]

    # Find K roots closest to unit circle
    dist2UnitCircle = abs.(abs.(polyroots2) .- 1)
    p = sortperm(dist2UnitCircle)[1:K]
    closestRoots = polyroots2[p]

    # Calculate frequencies from roots
    f_hat = (angle.(closestRoots)/(2*pi) .+ 1).%1
end


text:
function Base.show(io::IO, result::TopsisResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::ARASResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Orderings: ")
    println(io, result.orderings)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::CODASResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::COPRASResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::CRITICResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::CoCoSoResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::EDASResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::ElectreResult)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::GreyResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ordering: ")
    println(io, result.ordering)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::MABACResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::MAIRCAResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::MARCOSResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::MooraResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::PrometheeResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::SawResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::VikorResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::WASPASResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

function Base.show(io::IO, result::WPMResult)
    println(io, "Scores:")
    println(io, result.scores)
    println(io, "Ranking: ")
    println(io, result.ranking)
    println(io, "Best indices:")
    println(io, result.bestIndex)
end

max_stars_repo_name: drvinceknight/JMcDM

max_issues_repo_name: drvinceknight/JMcDM

max_forks_repo_name: drvinceknight/JMcDM


text:
using MAT
using MultipleViewGeometry, Base.Test
using MultipleViewGeometry.ModuleCostFunction
using MultipleViewGeometry.ModuleTypes
using BenchmarkTools, Compat
using StaticArrays

#"C:\Users\Spock\Desktop\tempdata"

file = matopen("test/debug/F.mat")
F = read(file,"F")
close(file)

file = matopen("test/debug/mPts1List.mat")
mpts1list = read(file,"mPts1List")
close(file)

file = matopen("test/debug/mPts2List.mat")
mpts2list = read(file,"mPts2List")
close(file)

file = matopen("test/debug/covs.mat")
covs = read(file,"covs")
close(file)

file = matopen("test/debug/Ct.mat")
Ct = read(file,"Ct")
close(file)

npts, dim  = size(mpts1list[1])

s²  = 1e-7
s = sqrt(s² )
Λ₁ =  [SMatrix{3,3}(s² *diagm([1.0,1.0,0.0])) for i = 1:npts]
Λ₂ =  [SMatrix{3,3}(s² *diagm([1.0,1.0,0.0])) for i = 1:npts]

ntrial = 10000
ℱ = zeros(9,ntrial)
for itrial = 1:ntrial
  mpts1 = mpts1list[itrial]
  mpts2 = mpts2list[itrial]
  𝒪 = [Point2DH(vcat(mpts1[i,:],[1])) for i =1:npts]
  𝒪ʹ= [Point2DH(vcat(mpts2[i,:],[1])) for i =1:npts]

  𝐅₀ = estimate(FundamentalMatrix(),DirectLinearTransform(), (𝒪, 𝒪ʹ))
  𝐅 = estimate(FundamentalMatrix(),
                           FundamentalNumericalScheme(reshape(𝐅₀,9,1), 5, 1e-10),
                                                             (Λ₁,Λ₂), (𝒪, 𝒪ʹ))
  𝐟 = reshape(𝐅,9,1)
  𝐟 = 𝐟 / norm(𝐟)
  ℱ[:,itrial] = 𝐟
end

norm(F[:,1])
norm(ℱ[:,1])

F[:,1]
ℱ[:,1]

index = 2
mpts1 = mpts1list[index]
mpts2 = mpts2list[index]
𝒪 = [Point2DH(vcat(mpts1[i,:],[1])) for i =1:npts]
𝒪ʹ= [Point2DH(vcat(mpts2[i,:],[1])) for i =1:npts]
T1 = hartley_transformation(𝒪)
T2 = hartley_transformation(𝒪ʹ)

𝐅₀  = estimate(FundamentalMatrix(), DirectLinearTransform(), (𝒪, 𝒪ʹ))
𝐅 = estimate(FundamentalMatrix(),
                        FundamentalNumericalScheme(reshape(𝐅₀,9,1), 5, 1e-10),
                                                          (Λ₁,Λ₂), (𝒪, 𝒪ʹ))
𝐟 = reshape(𝐅,9,1)
𝐟 = 𝐟 / norm(𝐟)

𝐅 = 𝐅 / norm(𝐅)

𝛍 = mean(ℱ,2)
𝛍 = 𝛍 / norm(𝛍)
d = length(𝛍)
𝐏 = eye(d) - norm(𝛍)^-2 * (𝛍*𝛍')
Cₜ = zeros((d,d))
for itrial = 1:ntrial
    𝐟 =  ℱ[:,itrial]
    Cₜ = Cₜ + 𝐏*(𝐟-𝛍) *(𝐟-𝛍)'*𝐏'
end
Cₜ = Cₜ / ntrial

#
# 𝛍 = mean(F,2)
# 𝛍 = 𝛍 / norm(𝛍)
# d = length(𝛍)
# 𝐏 = eye(d) - norm(𝛍)^-2 * (𝛍*𝛍')
# Cₜ = zeros((d,d))
# for itrial = 1:ntrial
#     𝐟 =  F[:,itrial]
#     Cₜ = Cₜ + 𝐏*(𝐟-𝛍) *(𝐟-𝛍)'*𝐏'
# end
# Cₜ = Cₜ / ntrial

# Q  = ℱ;
# for itrial = 1:ntrial
#     q =  Q[:,itrial]
#     Q[:,itrial] = Q[:,itrial] / sign(q[end-1])
# end
#
# 𝛍 = mean(Q,2)
# 𝛍 = 𝛍 / norm(𝛍)
# d = length(𝛍)
# 𝐏 = eye(d) - norm(𝛍)^-2 * (𝛍*𝛍')
# Cₜ = zeros((d,d))
# for itrial = 1:ntrial
#     𝐟 =  Q[:,itrial]
#     Cₜ = Cₜ + 𝐏*(𝐟-𝛍) *(𝐟-𝛍)'*𝐏'
# end
# Cₜ = Cₜ / ntrial

#
# A = Cₜ
#
#
# # Validate the covariance matrix of an estimate based on the AML cost function.
C = covariance_matrix(AML(),FundamentalMatrix(), reshape(𝐅,9,1), (Λ₁,Λ₂), (𝒪 , 𝒪ʹ))
C1 = covariance_matrix(AML(),FundamentalMatrix(),  𝐟, (Λ₁,Λ₂), (𝒪 , 𝒪ʹ))
C2 = covariance_matrix_debug(AML(),FundamentalMatrix(), reshape(𝐅,9,1), (Λ₁,Λ₂), (𝒪 , 𝒪ʹ))

C3 = covariance_matrix_debug(AML(),FundamentalMatrix(), 𝐟, (Λ₁,Λ₂), (𝒪 , 𝒪ʹ))

Cn = covariance_matrix_normalised(AML(),FundamentalMatrix(), 𝐟, (Λ₁,Λ₂), (𝒪 , 𝒪ʹ))

# (diag(C3)./diag(Ct) .-1)*100
#
# norm(diag(C3)-diag(Ct)) / norm(diag(Ct)) * 100
#
norm(C3-Ct) / norm(Ct) * 100
norm(C2-Ct) / norm(Ct) * 100
norm(Cn-Ct) / norm(Ct) * 100
norm(C1-Ct) / norm(Ct) * 100

# Cx = covariance_matrix(AML(),FundamentalMatrix(),  𝐟 , (Λ₁,Λ₂), (𝒪 , 𝒪ʹ))
#
#
# P = (1/norm(𝐟)) * (eye(9) - ((𝐟*𝐟') / norm(𝐟)^2) )
#
# P2 = eye(9) - norm(𝐟)^-2 * (𝐟*𝐟')
#
# cost(AML(),FundamentalMatrix(), 𝐟, (Λ₁,Λ₂), (𝒪 , 𝒪ʹ))

# A = Float64.([3 0 2; 2 0 -2; 0 1 1])
# size(A)
# B = similar(A)
# rows, cols = size(A)
# for r = 1:rows
#     for c =1:cols
#         B[r,c] = det(A)
#     end
# end

Q = reshape(𝐟,(3,3))


text:
function matchregret_dbr(I::Int64, r::Vector{Vector{Float64}}, gs::GameSet,
                            pdbr::Vector{Float64}, σfix::Vector{Vector{Float64}})
    ni_stage, na_stage = gs.ni_stage, gs.na_stage
    na = numactions(I, ni_stage, na_stage)
    σ = Vector{Float64}(undef, na)
    if I <= ni_stage[3] || I > sum(ni_stage[1:5]) # attacker infoset
        for a in 1:na
            denom = sum(max(r[I][b], 0.0) for b in 1:na)
            if denom > 0.0
                σ[a] = pdbr[I] * σfix[I][a] + (1 - pdbr[I]) * (max(r[I][a], 0.0) / denom)
            else
                σ[a] = pdbr[I] * σfix[I][a] + (1 - pdbr[I]) * (1.0 / na)
            end
        end
    else # defender infoset
        for a in 1:na
            denom = sum(max(r[I][b], 0.0) for b in 1:na)
            if denom > 0.0
                σ[a] = max(r[I][a], 0.0) / denom
            else
                σ[a] = 1.0 / na
            end
        end
    end
    return σ
end

function updateutility_dbr!(player::Int64, I::Int64, h::SVector, σ::Vector{Float64},
        u::Float64, πi::Float64, πo::Float64,
        r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}},
        g::AirDefenseGame, gs::GameSet, depth::Int64, pdbr::Vector{Float64}, σfix::Vector{Vector{Float64}})
    A, ni_stage, na_stage = gs.A, gs.ni_stage, gs.na_stage
    na = numactions(I, ni_stage, na_stage)
    m = Vector{Float64}(undef, na)
    nextactions = actions(depth, A)
    # @show na[I]
    for a in 1:na
        ha = setindex(h, nextactions[a], depth)
        if getplayer(depth) == player
            πip = σ[a] * πi
            # up = updatetree_dbr!(ha, player, πip, πo, r, s, g, gs)
            up = updatetree_dbr!(ha, player, πip, πo, r, s, g, gs, depth + 1, pdbr, σfix)
            m[a] = up
            u = u + σ[a] * up
        else
            πop = σ[a] * πo
            # up = updatetree_dbr!(ha, player, πi, πop, r, s, g, gs)
            up = updatetree_dbr!(ha, player, πi, πop, r, s, g, gs, depth + 1, pdbr, σfix)
            u = u + up
        end
    end
    return m, u
end

function updateregret_dbr!(player::Int64, I::Int64, h::SVector, σ::Vector{Float64},
        u::Float64, m::Vector{Float64}, πi::Float64,
        r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}}, gs::GameSet, depth ::Int64)
    ni_stage, na_stage = gs.ni_stage, gs.na_stage
    if getplayer(depth) == player
        for a in 1:numactions(I, ni_stage, na_stage)
            r[I][a] = r[I][a] + m[a] - u
            s[I][a] = s[I][a] + πi * σ[a]
        end
    end
end

function updatetree_dbr!(h::SVector, player::Int64, πi::Float64, πo::Float64,
    r::Vector{Vector{Float64}}, s::Vector{Vector{Float64}}, g::AirDefenseGame,
    gs::GameSet, depth::Int64, pdbr::Vector{Float64}, σfix::Vector{Vector{Float64}})
    # πo is short for \vec{π}_{-i}
    A = gs.A
    (terminal(depth)) && (return leafutility(h, player, πi, πo, g, gs))
    if getplayer(depth) == 3
        as = actions(depth, A)
        chance_weights = pweights([chance(depth, a, g, gs) for a in as])
        a = as[sample(chance_weights)]
        ha = setindex(h, a, depth)
        # return updatetree_dbr!(ha, player, πi, πo, r, s, g, gs)
        return updatetree_dbr!(ha, player, πi, πo, r, s, g, gs, depth + 1, pdbr, σfix)
    end
    I = infoset(h, depth, g, gs)
    u = 0.0
    σ = matchregret_dbr(I, r, gs, pdbr, σfix)
    m, u = updateutility_dbr!(player, I, h, σ, u, πi, πo, r, s, g, gs, depth, pdbr, σfix)
    updateregret_dbr!(player, I, h, σ, u, m, πi, r, s, gs, depth)
    return u
end

function cfr_dbr(g::AirDefenseGame, gs::GameSet, pdbr::Vector{Float64}, σfix::Vector{Vector{Float64}};
        iterlimit::Int = 100_000, timelimit::Number = 600, tol::Float64 = 5e-5)
    ni, ni_stage, na_stage, players = gs.ni, gs.ni_stage, gs.na_stage, gs.players
    r = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # regret
    s = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # cumulative strategies
    u = 0.0
    u1 = Float64[]
    u1deriv = 0
    u1prev = 0
    converged = false
    iter = 1
    tstart = now()
    while !converged && tominute(now() - tstart) < timelimit && iter < iterlimit
        for player in players
            h0 = SVector(0, 0, 0, 0, 0, 0)
            u = updatetree_dbr!(h0, player, 1.0, 1.0, r, s, g, gs, 1, pdbr, σfix)
            (player == 1) && push!(u1, u)
        end
        u1deriv = u1deriv * (1 - 1 / iter) + (u1[iter] - u1prev) * (1 / iter)
        converged = abs(u1deriv) < tol
        u1prev = u1[iter]
        iter % 10 == 0 && print("\rIter: $iter")
        iter += 1
    end
    denoms = [sum(i) for i in s]
    σ = Vector{Vector{Float64}}(undef, gs.ni)
    for i in eachindex(s)
        if denoms[i] == 0
            σ[i] = fill(1 / length(s[i][:]), length(s[i][:]))
        else
            σ[i] = s[i][:] ./ denoms[i]
        end
    end
    return u1, r, s, σ, converged
end

function strat_cfr_dbr(s::Vector{Vector{Float64}})
    denoms = [sum(i) for i in s]
    σ = Vector{Vector{Float64}}(undef, length(s))
    for i in eachindex(s)
        if denoms[i] == 0
            σ[i] = fill(1 / length(s[i][:]), length(s[i][:]))
        else
            σ[i] = s[i][:] ./ denoms[i]
        end
    end
    return σ
end

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense