text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -3.3773826343471955 <= q <= 0.8255956851677421)

                     add_NL_constraint(m, :(softplus(-0.5817410001457377 + -0.38051045291676777*softplus(-0.7752766980375418 + 0.20887750352529233*softplus(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + -0.43558326218422394*softplus(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + -0.5061279411570236*softplus(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.828517562717189*softplus(-0.31567568822764036 + 0.26872838410166544*softplus(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.5533650825489542*softplus(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.11778701597368801*softplus(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.13996198025628548*softplus(0.6085123540494801 + -0.012952465394535206*softplus(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.6228552648202648*softplus(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.7276708970518393*softplus(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2])))) + softplus(-0.08252400134510385 + 0.3442818429504264*softplus(-0.7752766980375418 + 0.20887750352529233*softplus(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + -0.43558326218422394*softplus(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + -0.5061279411570236*softplus(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.7813134096358496*softplus(-0.31567568822764036 + 0.26872838410166544*softplus(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.5533650825489542*softplus(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.11778701597368801*softplus(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.5888247957791601*softplus(0.6085123540494801 + -0.012952465394535206*softplus(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.6228552648202648*softplus(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.7276708970518393*softplus(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2])))) + softplus(-0.2140493400207517 + 0.6401317985056014*softplus(-0.7752766980375418 + 0.20887750352529233*softplus(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + -0.43558326218422394*softplus(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + -0.5061279411570236*softplus(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + -0.4453350709432602*softplus(-0.31567568822764036 + 0.26872838410166544*softplus(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.5533650825489542*softplus(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.11778701597368801*softplus(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2]))) + 0.9148591417614278*softplus(0.6085123540494801 + -0.012952465394535206*softplus(0.3913052313874794 + -0.22863414084682043*$(x[1]) + 0.3783772947965627*$(x[2])) + 0.6228552648202648*softplus(-0.4746875961429966 + 0.3096610657098613*$(x[1]) + -0.16020248217350597*$(x[2])) + 0.7276708970518393*softplus(0.7539608271466713 + -0.29418097465562054*$(x[1]) + -0.3042211456088699*$(x[2])))) - $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 validate_available_devices(device::Type{D}, devices) where {D <: PSY.Device}
    if isempty(devices)
        @warn("The data doesn't include devices of type $(device), consider changing the device models")
        return false
    end
    return true
end

function validate_available_services(service::Type{S}, services) where {S <: PSY.Service}
    if isempty(services)
        @warn("The data doesn't include services of type $(service), consider changing the service models")
        return false
    end
    return true
end

ext: jl


text:
@testset "decisionTree.jl" begin

    data_surv = CSO.dataset("cso_raw")
    data_trudeau = CSO.dataset("Montreal-Trudeau_15min_precipitations")
    data = innerjoin(data_surv, data_trudeau, on=:Date)
    dropmissing!(data, :Duration)
    train, test = train_test_year(data, 2020)

    varnames = [:d15min, :d1h, :d6h]
    Xtrain = Matrix(train[:,varnames])
    ytrain = train.Duration
    Xtest = Matrix(test[:,varnames])
    model = DecisionTreeClassifier(max_depth=2)

    tree = DecisionTree.fit!(model, Xtrain, ytrain)
    predictions = DecisionTree.predict(tree, Xtest)

    RMSD = rmsd(predictions, test.Duration)

    #The RMSD is well derived
    @test RMSD ≈ DTree(model, "Duration", varnames, train, test)

    #Predictions are well derived
    @test predictions ≈ DTree(model, "Duration", varnames, train, test, false, true, false)

    #The tree is well built
    @test tree == DTree(model, "Duration", varnames, train, test, true, false, false)

end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
using RecipesBase

function Base.real(h::History)
    steps, vals = get(h)
    _h = History(real(eltype(vals)), eltype(steps))
    _h.iterations = steps
    _h.values = real(vals)
    return _h
end

function Base.imag(h::History)
    steps, vals = get(h)
    _h = History(real(eltype(vals)), eltype(steps))
    _h.iterations = steps
    _h.values = imag(vals)
    return _h
end

function Base.getindex(h::History, step::Number)
    steps, vals = get(h)
    id = findfirst(steps .== step)
    return vals[id]
end

function Base.getindex(h::History, ids)
    steps, vals = get(h)

    id_start = findfirst(steps .>= first(ids))
    id_end = findlast(steps .<= last(ids))

    _h = History(eltype(vals), eltype(steps))
    _h.iterations = steps[id_start:id_end]
    _h.values = vals[id_start:id_end]
    return _h
end

Base.firstindex(h::History) = first(h.iterations)
Base.lastindex(h::History) = last(h.iterations)

@recipe function plot(h::History)
    markershape --> :cross
    title       --> "Value History"
    get(h)
end


text:
# module wrap_cpp
module wt

using Clang.cindex

export
    CustomArg,
    ArgProxy,
    WrMethod
    #emitc_init,
    #emitc_post

################################################################################
# Method entry point
################################################################################

immutable WrappedMethod
    name::ASCIIString
    method::CXXMethod
    parent::ClassDecl
    args::Array{Any,1}
    hasretval::Bool
end

###
# Parse a C++ method and build a collection of
# argument types, possibly with proxy implementation.
###
function analyze_method(method::CXXMethod)
    arg_list = get_args(method)

    hasreturn = returns_value(method)

    args = Any[]
    for arg in arg_list
        proxy = get_proxy(spelling(cu_type(arg)))
        (proxy == None) ? push!(args, arg) : push!(args, proxy)
    end
    return WrappedMethod(spelling(method),
                         method,
                         cindex.getCursorLexicalParent(method),
                         args,
                         hasreturn)
end

################################################################################
# Extensible type-proxy interface
################################################################################
# Usage:
#   - create custom handler ArgType <: ArgProxy
#   - write emitc_init and emitc_post for that type
#   - set TypeProxies[type_spelling] = ArgType

abstract TypeProxy <: CLType

typealias __TypeProxiesType Dict{String, Type}
const TypeProxies = Dict{String, Type}()

function get_proxy(typename::String)
    get(TypeProxies::__TypeProxiesType, typename, None)
end

################################################################################
# Generation of C wrapper
################################################################################

#emitc(out::IO, arg::ArgProxy) = emitc(out, arg.cltype)

function emit(out::IO, t::CLType)
    print(out, cl_to_c[typeof(t)])
end

function emit(out::IO, t::Union(Record, Typedef))
    print(out, spelling(cindex.getTypeDeclaration(t)))
end

function emit(out::IO, t::cindex.ConstantArray)
    arrtype = cindex.getArrayElementType(t)
    arrsize = cindex.getArraySize(t)
    emit(out, arrtype)
    opensquare(out); print(out, arrsize); closesquare(out)
end

function emit(out::IO, t::cindex.Pointer)
    pointee = cindex.pointee_type(t)
    emit(out, pointee)
    print(out, "*")
end

function emit(out::IO, parm::cindex.ParmDecl)
    emit(out, cu_type(parm))
end

function check_args(args)
    check = true
    for arg in args
        argt = cu_type(arg)
        check &= isa(argt, cindex.LValueReference)
        check &= isa(argt, cindex.FirstBuiltin)
    end
    return check
end

function emit_args(out::IO, args)
    emit_arg(out,arg,docomma=false) = begin
        emit(out,arg)
        space(out)
        print(out,spelling(arg))
        docomma && (comma(out); space(out))
    end
    for i = 1:length(args)
        emit_arg(out,args[i], i < length(args))
    end
end

function wrap(out::IO, method::cindex.CXXMethod, nameid::Int)
    buf = IOBuffer()

    # bail out if any argument is not supported
    args = get_args(method)
    if (!check_args(args)) return end

    methodname = spelling(method)
    parentdecl = cindex.getCursorLexicalParent(method)
    parentname = spelling(parentdecl)

    # emit return type and name
    rettype = return_type(method)
    emit(buf, rettype)
    space(buf)
    print(buf, parentname, "_", methodname, nameid)

    # check for a return value
    hasreturn = returns_value(rettype)

    # emit arguments
    openparen(buf)
        print(buf, parentname, "* this")
        length(args) > 0 && comma(buf)
        space(buf)
        emit_args(buf, args)
    closeparen(buf)
    space(buf)

    opencurly(buf)
    newline(buf)
    print(buf, "  ")

    # emit method call
    if(hasreturn)
        print(buf, "return")
        space(buf)
    end
    print(buf, "this->", methodname)
    openparen(buf)
        emit_args(buf, args)
    closeparen(buf)

    # close
    newline(buf)
    closecurly(buf)
    newline(buf)

    # print the buffer to the output
    print(out, takebuf_string(buf))
    # TODO: return information about the thing we wrapped?
end

function wrap(out::IO, top::cindex.ClassDecl)
    println("Wrapping class: ", name(top))
    MethodCount = Dict{ASCIIString, Int}()

    for decl in [c for c in children(top)]
        declname = spelling(decl)
        if isa(decl, cindex.CXXMethod)
            id = (MethodCount[declname] = get(MethodCount, declname, 0) + 1)
            wrap(out, decl, id)
        end
    end
end

cl_to_c = {
    cindex.VoidType       => "void",
    cindex.FirstBuiltin   => "void",
    cindex.BoolType       => "bool",
    cindex.Char_U         => "unsigned char",
    cindex.UChar          => "unsigned char",
    cindex.Char16         => "char16_t",
    cindex.Char32         => "char32_t",
    cindex.UShort         => "unsigned short",
    cindex.UInt           => "unsigned int",
    cindex.ULong          => "unsigned long",
    cindex.ULongLong      => "unsigned long long",
    cindex.Char_S         => "char",
    cindex.SChar          => "char",
    cindex.WChar          => "wchar_t",
    cindex.Short          => "short",
    cindex.IntType        => "int",
    cindex.Long           => "long",
    cindex.LongLong       => "long long",
    cindex.Float          => "float",
    cindex.Double         => "double",
    cindex.LongDouble     => "long double",
    cindex.Enum           => "enum",
    cindex.NullPtr        => "NULL",
    cindex.UInt128        => "uint128_t"
    }

################################################################################
# Construction of Julia wrapper
################################################################################

# TODO: FIXME
function wrapjl(out::IO, t::cindex.ConstantArray)
    print(out, "FIXEDARRAY")
end

function wrapjl(out::IO, t::Union(cindex.Record, cindex.Typedef))
    print(out, spelling(cindex.getTypeDeclaration(t)))
end

function wrapjl(out::IO, arg::cindex.CLType)
    print(out, string(cl_to_jl[typeof(arg)]))
end

function wrapjl(out::IO, ptr::cindex.Pointer)
    pointee = pointee_type(ptr)

    print(out, "Ptr")
    opencurly(out)
    wrapjl(out, pointee)
    closecurly(out)
end

function wrapjl(out::IO, parm::cindex.ParmDecl)
    wrapjl(out, cu_type(parm))
end

function wrapjl_args(out::IO, args)
    emit_arg(out,arg,docomma=false) = begin
        print(out,spelling(arg))
        print(out,"::")
        wrapjl(out,arg)
        docomma && (comma(out); space(out))
    end
    for i = 1:length(args)
        emit_arg(out,args[i], i < length(args))
    end
end

function wrapjl(out::IO, method::cindex.CXXMethod, id::Int)
    buf = IOBuffer()

    methodname = spelling(method)
    parentdecl = cindex.getCursorLexicalParent(method)
    parentname = spelling(parentdecl)

    args = get_args(method)
    if (!check_args(args)) return end                       # TODO: warning?

    print(buf, "@method")                                   # emit Julia call
    space(buf)
    print(buf, parentname)                                  # parent class
    space(buf)
    print(buf, methodname)
    space(buf)
    wrapjl(buf, return_type(method))                        # return type
    space(buf)
    openparen(buf)
    wrapjl_args(buf, args)                                  # arguments
    closeparen(buf)
    space(buf)
    print(buf, methodname, id)             # C name
    newline(buf)
    print(out, takebuf_string(buf))
end

function wrapjl(out::IO, class::cindex.ClassDecl)
    MethodCount = Dict{ASCIIString, Int}()

    for decl in [c for c in children(class)]
        declname = spelling(decl)
        if isa(decl, cindex.CXXMethod)
            id = (MethodCount[declname] = get(MethodCount, declname, 0) + 1)
            wrapjl(out, decl, id)
        end
    end
end

cl_to_jl = {
    cindex.VoidType         => Void,
    cindex.BoolType         => Bool,
    cindex.Char_U           => Uint8,
    cindex.UChar            => :Cuchar,
    cindex.Char16           => Uint16,
    cindex.Char32           => Uint32,
    cindex.UShort           => Uint16,
    cindex.UInt             => Uint32,
    cindex.ULong            => :Culong,
    cindex.ULongLong        => :Culonglong,
    cindex.Char_S           => Uint8,
    cindex.SChar            => Uint8,
    cindex.WChar            => Char,
    cindex.Short            => Int16,
    cindex.IntType          => :Cint,
    cindex.Long             => :Clong,
    cindex.LongLong         => :Clonglong,
    cindex.Float            => :Cfloat,
    cindex.Double           => :Cdouble,
    cindex.LongDouble       => Float64,
    cindex.Enum             => :Cint,
    cindex.NullPtr          => C_NULL,
    cindex.UInt128          => Uint128,
    cindex.FirstBuiltin     => Void,
    "size_t"                => :Csize_t,
    "ptrdiff_t"             => :Cptrdiff_t
    }

################################################################################
# Utility functions
################################################################################

space(io::IO)      = print(io, " ")
openparen(io::IO)  = print(io, "(")
closeparen(io::IO) = print(io, ")")
opencurly(io::IO)  = print(io, "{")
closecurly(io::IO) = print(io, "}")
opensquare(io::IO) = print(io, "[")
closesquare(io::IO)= print(io, "]")
comma(io::IO)      = print(io, ",")
newline(io::IO)    = print(io, "\n")

function base_classes(class::cindex.ClassDecl)
    return [cindex.getCursorReferenced(c)
            for c in cindex.search(class, cindex.CXXBaseSpecifier)]
end

returns_value(m::CXXMethod) = returns_value(cu_type(m))
function returns_value(rtype::CLType)
    hasret = true
    hasret &= ~isa(rtype, VoidType)
    return hasret
end

function get_args(method::cindex.CXXMethod)
    args = CLCursor[]
    for c in children(method)
        isa(c,cindex.ParmDecl) ? push!(args, c) : break
    end
    return args
end

end # module wrap_cpp


text:
using FeedbackParticleFilters, Statistics

println("Testing ConstantGainApproximation.jl:")

@testset "ConstantGainApproximation.jl" begin

    h(x)      = x[1:2].^2
    pos       = randn(3, 10)
    eq        = PoissonEquation(h, pos)
    solve!(eq, ConstantGainApproximation())

    print("  solver")
    for i in 2:10, j in 1:2
        print("-")
        @test eq.gain[:,1,j] == eq.gain[:,i,j]
    end
    for i in 1:3, j in 1:2
        print("-")
        @test eq.gain[i,1,j] == Statistics.mean(eq.positions[i,:] .* (eq.H[j,:] .- eq.mean_H[j,:]))
    end
    println("DONE")

end; #ConstantGainApproximation.jl

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
struct state
    x :: Vector{<:Real}
    log_likelihood :: Vector{<:Real}
    grad_log_likelihood :: Vector{<:Real}
end

function state(x::Vector{<:Real}, log_likelihood :: Real)
    return state(x, [log_likelihood], zeros(length(x)))
end

function state(x::Vector{<:Real}, log_likelihood :: Real, grad_log_likelihood :: Vector{<:Real})
    return state(x, [log_likelihood], grad_log_likelihood)
end

struct chain_diagnostics
    log_likelihood  :: Vector{<:Real}
    acceptance_ratio :: Vector{<:Real}
end

function chain_diagnostics(num_samples :: Int)
    return chain_diagnostics(zeros(num_samples), zeros(num_samples))
end


text:
using DimensionalData, Test, LinearAlgebra, Statistics

A = [1.0 2.0 3.0;
     4.0 5.0 6.0]
x, y, z = X([:a, :b]), Y(10.0:10.0:30.0), Z()
dimz = x, y
da1 = DimArray(A, (x, y); name=:one)
da2 = DimArray(Float32.(2A), (x, y); name=:two)
da3 = DimArray(Int.(3A), (x, y); name=:three)
da4 = DimArray(cat(4A, 5A, 6A, 7A; dims=3), (x, y, z); name=:extradim)

s = DimStack((da1, da2, da3))
mixed = DimStack((da1, da2, da4))

@testset "Constructors" begin
    @test DimStack((one=A, two=2A, three=3A), dimz) == s
    @test DimStack(da1, da2, da3) == s
    @test DimStack((one=da1, two=da2, three=da3), dimz) == s
end

@testset "Properties" begin
    @test DimensionalData.dims(s) == dims(da1)
    @test keys(data(s)) == (:one, :two, :three)
    @test keys(data(mixed)) == (:one, :two, :extradim)
    da1x = s[:one]
    @test parent(da1x) === parent(da1)
    @test dims(da1x) === dims(da1)
    @test size(da1x) === (2, 3)
    @test size(mixed) === (2, 3, 4)
    @test size(da1x, X) === 2
    @test size(mixed, 3) === 4
    @test axes(da1x) === (Base.OneTo(2), Base.OneTo(3))
    @test axes(mixed) === (Base.OneTo(2), Base.OneTo(3), Base.OneTo(4))
    @test axes(da1x, X) === Base.OneTo(2)
    @test axes(mixed, 2) === Base.OneTo(3)
end

@testset "map" begin
    @test values(map(a -> a .* 2, s)) == values(DimStack(2da1, 2da2, 2da3))
    @test dims(map(a -> a .* 2, s)) == dims(DimStack(2da1, 2da2, 2da3))
    @test map(a -> a[1], s) == (one=1.0, two=2.0, three=3.0)
    @test values(map(a -> a .* 2, s)) == values(DimStack(2da1, 2da2, 2da3))
    @test map(+, s, s, s) == map(a -> a .* 3, s)
    @test_throws ArgumentError map(+, s, mixed)
end

@testset "Methods with no arguments" begin
    @testset "permuting methods" begin
        @test data(permutedims(s)) ==
            (one=[1.0 4.0; 2.0 5.0; 3.0 6.0],
             two=[2.0 8.0; 4.0 10.0; 6.0 12.0],
             three=[3.0 12.0; 6.0 15.0; 9.0 18.0])
        @test adjoint(s) == DimStack(adjoint(da1), adjoint(da2), adjoint(da3))
        @test transpose(s) == DimStack(transpose(da1), transpose(da2), transpose(da3))
        @test Transpose(s) == DimStack(Transpose(da1), Transpose(da2), Transpose(da3))
        @test data(rotl90(s)) ==
            (one=[3.0 6.0;  2.0 5.0;  1.0 4.0],
             two=[6.0f0 12.0f0; 4.0f0 10.0f0; 2.0f0 8.0],
           three=[9 18; 6 15; 3 12])
        @test rotl90(s) == DimStack(rotl90(da1), rotl90(da2), rotl90(da3))
        @test rotr90(s) == DimStack(rotr90(da1), rotr90(da2), rotr90(da3))
        @test rot180(s) == DimStack(rot180(da1), rot180(da2), rot180(da3))
    end

    @test cor(s) isa DimStack
    @test cov(s) isa DimStack

    @test inv(s[1:2, 1:2]) isa DimStack

    @testset "reducing methods" begin
        @test sum(s) === (one=21.0, two=42.0f0, three=63)
        @test prod(s) === (one=720.0, two=46080.0f0, three=524880)
        @test Base.minimum(s) === (one=1.0, two=2.0f0, three=3)
        @test maximum(s) === (one=6.0, two=12.0f0, three=18)
        @test extrema(s) === (one=(1.0, 6.0), two=(2.0f0, 12.0f0), three=(3, 18))
        @test mean(s) === (one=3.5, two=7.0f0, three=10.5)
        @test std(s) === (one=1.8708286933869707, two=3.7416575f0, three=5.612486080160912)
        @test var(s) === (one=3.5, two=14.0f0, three=31.5)
        @test median(s) === (one=3.5, two=7.0f0, three=10.5)
    end
end

@testset "Methods with a dims keyword argument" begin
    @testset "reducing methods" begin
        sum_ = sum(s; dims=X)
        @test sum_  == DimStack(sum(da1; dims=X), sum(da2; dims=X), sum(da3; dims=X))
        @test dims(sum_) == (X([:combined], Categorical()), dims(da1, Y))
        @test prod(s; dims=X) == DimStack(prod(da1; dims=X), prod(da2; dims=X), prod(da3; dims=X))
        @test dims(maximum(s; dims=X)) ==
            dims(DimStack(maximum(da1; dims=X), maximum(da2; dims=X), maximum(da3; dims=X)))
        @test minimum(s; dims=X) == DimStack(minimum(da1; dims=X), minimum(da2; dims=X), minimum(da3; dims=X))
        @test mean(s; dims=X()) == DimStack(mean(da1; dims=X), mean(da2; dims=X), mean(da3; dims=X))
        @test std(s; dims=X) == DimStack(std(da1; dims=X), std(da2; dims=X), std(da3; dims=X))
        @test var(s; dims=X) == DimStack(var(da1; dims=X), var(da2; dims=X), var(da3; dims=X))
        @test median(s; dims=X) == DimStack(median(da1; dims=X), median(da2; dims=X), median(da3; dims=X))
        mean(mixed; dims=X)
        mean(mixed; dims=Y)
        @test_broken mean(mixed; dims=Z)
    end

    @testset "dim duplicating methods" begin
        @test cor(s; dims=X) isa DimStack
        @test cov(s; dims=X) isa DimStack
    end

    @testset "dim dropping methods" begin
        @test DimensionalData.layers(dropdims(s[X([1])]; dims=X)) ==
            (one=DimArray([1.0, 2.0, 3.0], dims(da1, Y)),
             two=DimArray([2.0, 4.0, 6.0], dims(da1, Y)),
             three=DimArray([3.0, 6.0, 9.0], dims(da1, Y)))
    end

end

@testset "permuting methods with an argument" begin
    @test permutedims(s, (Y, X)) ==
        DimStack(permutedims(da1, (Y, X)), permutedims(da2, (Y, X)), permutedims(da3, (Y, X)))
    @test PermutedDimsArray(s, (Y(), X())) ==
        DimStack(PermutedDimsArray(da1, (Y, X)), PermutedDimsArray(da2, (Y, X)), PermutedDimsArray(da3, (Y, X)))
    rot = rotl90(s, 1)
    @test rot isa DimStack
    @test typeof(data(rot)) == NamedTuple{(:one, :two, :three),Tuple{Matrix{Float64},Matrix{Float32},Matrix{Int64}}}
    @test data(rot) ==
        (one=[3.0 6.0;  2.0 5.0;  1.0 4.0],
         two=[6.0f0 12.0f0; 4.0f0 10.0f0; 2.0f0 8.0f0],
       three=[9 18; 6 15; 3 12])
    @test rot[:one][X(At(:a)), Y(At(10.0))] == da1[X(At(:a)), Y(At(10.0))]
    @test rotr90(s, 2) == DimStack(rotr90(da1, 2), rotr90(da2, 2), rotr90(da3, 2))
    @test rot180(s, 1) == DimStack(rot180(da1, 1), rot180(da2, 1), rot180(da3, 1))
end

@testset "reducing methods that take a function" begin
    f = x -> 2x
    @test data(minimum(f, s; dims=X)) ==
        (one=[2.0 4.0 6.0], two=[4.0 8.0 12.0], three=[6.0 12.0 18.0])
    @test data(mean(x -> 2x, s; dims=X())) ==
        (one=mean(2da1; dims=X) , two=mean(2da2; dims=X), three=mean(2da3; dims=X))
    @test mean(f, s; dims=X) == DimStack(mean(f, da1; dims=X), mean(f, da2; dims=X), mean(f, da3; dims=X))
    @test reduce(+, s; dims=X) == DimStack(reduce(+, da1; dims=X), reduce(+, da2; dims=X), reduce(+, da3; dims=X))
    @test sum(f, s; dims=X) == DimStack(sum(f, da1; dims=X), sum(f, da2; dims=X), sum(f, da3; dims=X))
    @test prod(f, s; dims=X) == DimStack(prod(f, da1; dims=X), prod(f, da2; dims=X), prod(f, da3; dims=X))
    @test maximum(f, s; dims=X) == DimStack(maximum(f, da1; dims=X), maximum(f, da2; dims=X), maximum(f, da3; dims=X))
    @test extrema(f, s; dims=X) == DimStack(extrema(f, da1; dims=X), extrema(f, da2; dims=X), extrema(f, da3; dims=X))
    @test reduce(+, s) == (one=21.0, two=42.0, three=63.0)
    @test sum(f, s) == (one=42.0, two=84.0, three=126.0)
    @test prod(f, s) == (one=46080., two=2.94912e6, three=3.359232e7)
    @test Base.minimum(f, s) == (one=2.0, two=4.0, three=6.0)
    @test maximum(f, s) == (one=12.0, two=24.0, three=36.0)
    @test extrema(f, s) == (one=(2.0, 12.0), two=(4.0, 24.0), three=(6.0, 36.0))
    @test mean(f, s) == (one=7.0, two=14.0, three=21)
end


text:
function loadtexttree(filename::String)
    branches = CSV.read(joinpath(abspath(dirname(@__FILE__)),filename));
    return branches
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
include("../../MaximinOPF/src/MaximinOPF.jl")
using PowerModels
using JuMP
using SCS
using Ipopt
using ProxSDP
PowerModels.silence()

testcase = Dict(
	"file" => "data/case9.m",
	"PMOption" => SparseSDPWRMPowerModel,
 	"name" => "case9K3",
 	"attack_budget" => 3,
 	"inactive_indices" => [],
 	"protected_indices" => []
	)

pm_data = PowerModels.parse_file(testcase["file"])
pm_data["attacker_budget"] = testcase["attack_budget"] ###Adding another key and entry
pm_data["inactive_branches"] = testcase["inactive_indices"] ###Adding another key and entry
pm_data["protected_branches"] = testcase["protected_indices"] ###Adding another key and entry

sdp_solver=with_optimizer(ProxSDP.Optimizer, log_verbose=false, tol_primal=1e-6, tol_dual=1e-6 )
conic_solver=with_optimizer(SCS.Optimizer,verbose=0)
ip_solver=with_optimizer(Ipopt.Optimizer,print_level=0)

# nonconvex AC forms
nonconvex_ac=[ACPPowerModel, ACRPowerModel, ACTPowerModel]
for pm_form in nonconvex_ac
    println("Formulating and solving the form ",pm_form)
    pm=MaximinOPF.MinimaxOPFModel(pm_data, pm_form)
    JuMP.set_optimizer(pm.model,ip_solver)
    JuMP.optimize!(pm.model)
    println("\tOptimal value using powerform ", pm_form, " is: ",JuMP.objective_value(pm.model), " with status ",JuMP.termination_status(pm.model))
end

# linear approximations
linear_approx=[DCPPowerModel, DCMPPowerModel, NFAPowerModel]
for pm_form in linear_approx
    println("Formulating and solving the form ",pm_form)
    pm=MaximinOPF.MinimaxOPFModel(pm_data, pm_form)
    JuMP.set_optimizer(pm.model,conic_solver)
    JuMP.optimize!(pm.model)
    println("\tOptimal value using powerform ", pm_form, " is: ",JuMP.objective_value(pm.model), " with status ",JuMP.termination_status(pm.model))
end

# quadratic approximations
quadratic_approx=[DCPLLPowerModel, LPACCPowerModel]
for pm_form in quadratic_approx
    println("Formulating and solving the form ",pm_form)
    pm=MaximinOPF.MinimaxOPFModel(pm_data, pm_form)
    JuMP.set_optimizer(pm.model,ip_solver)
    JuMP.optimize!(pm.model)
    println("\tOptimal value using powerform ", pm_form, " is: ",JuMP.objective_value(pm.model), " with status ",JuMP.termination_status(pm.model))
end
# quadratic relaxations
quadratic_relax=[SOCWRPowerModel, SOCBFPowerModel, QCRMPowerModel, QCLSPowerModel]
for pm_form in quadratic_relax
    println("Formulating and solving the form ",pm_form)
    pm=MaximinOPF.MinimaxOPFModel(pm_data, pm_form)
    JuMP.set_optimizer(pm.model,ip_solver)
    JuMP.optimize!(pm.model)
    println("\tOptimal value using powerform ", pm_form, " is: ",JuMP.objective_value(pm.model), " with status ",JuMP.termination_status(pm.model))
end
quad_conic_relax=[SOCWRConicPowerModel, SOCBFConicPowerModel]
for pm_form in quad_conic_relax
    println("Formulating and solving the form ",pm_form)
    pm=MaximinOPF.MinimaxOPFModel(pm_data, pm_form)
    JuMP.set_optimizer(pm.model,conic_solver)
    JuMP.optimize!(pm.model)
    println("\tOptimal value using powerform ", pm_form, " is: ",JuMP.objective_value(pm.model), " with status ",JuMP.termination_status(pm.model))
end
# sdp relaxations
sdp_relax=[SDPWRMPowerModel, SparseSDPWRMPowerModel]
for pm_form in sdp_relax
    println("Formulating and solving the form ",pm_form)
    pm=MaximinOPF.MinimaxOPFModel(pm_data, pm_form)
    JuMP.set_optimizer(pm.model,sdp_solver)
    JuMP.optimize!(pm.model)
    println("\tOptimal value using powerform ", pm_form, " is: ",JuMP.objective_value(pm.model), " with status ",JuMP.termination_status(pm.model))
end


text:
# LeZi-Update predictor

"""
    LeZiUpdate{SymbolType}()

Creates a LeZiUpdate predictor for `SymbolType`. `SymbolType` can be any valid type including
`Char`, `Int` etc.,

## Examples
```julia-repl
julia> p = LeZiUpdate{Char}()
DiscretePredictors.LeZiUpdate{Char}(Array{Char,1}[], [*] (0)
, Char[], Char[])

julia> p = LeZiUpdate{Int}()
DiscretePredictors.LeZiUpdate{Int}(Array{Int,1}[], [*] (0)
, Int[], Int[])

julia> p = LeZiUpdate{String}()
DiscretePredictors.LeZiUpdate{String}(Array{String,1}[], [*] (0)
, String[], String[])
```

Reference:
Bhattacharya, Amiya, and Sajal K. Das. "LeZi-update: An information-theoretic framework for personal mobility tracking in PCS networks." Wireless Networks 8.2/3 (2002): 121-135.
"""
mutable struct LeZiUpdate{T} <: BasePredictor{T}
    dictionary::Vector{Vector{T}}
    model::Trie{T,Int}
    phrase::Vector{T}
    context::Vector{T}

    LeZiUpdate{T}() where {T} = new( Vector{Vector{}}(), Trie{T,Int}(), Vector{T}(), Vector{T}() )
    # (::Type{LeZiUpdate{T}}){T}() = new{T}( Vector{Vector{}}(), Trie{T,Int}(), Vector{T}(), Vector{T}() );
end

function add!( p::LeZiUpdate{T}, sym::T ) where {T}
    p.model.value += 1   # Update number of symbols seen by model

    push!( p.phrase, sym )
    if p.phrase ∉ p.dictionary
        push!( p.dictionary, p.phrase )

        suffix = p.phrase[1:end];
        while !isempty( suffix )
            prefix  = suffix[1:end];
            while !isempty( prefix )
                if haskey( p.model, prefix )
                    p.model[prefix] += 1
                else
                    p.model[prefix] = 1
                end
                popfirst!( prefix )
            end
            pop!( suffix )
        end
        p.context   = p.phrase[2:end]
        p.phrase    = Vector{T}()
    else
        p.context   = p.phrase[1:end]
    end

    nothing
end

function predict( p::LeZiUpdate{T} ) where {T}
    # Create a dictionary with symbols
    symbols = Dict( k => (p.model[[k]]/p.model.value) for k ∈ keys(children(p.model,Vector{T}())) )

    buffer  = Vector{T}()
    for i = length(p.context):-1:1
        pushfirst!( buffer, p.context[i] )
        # println( "    Buffer: ", buffer );
        list_of_children    = children( p.model, buffer )
        # Get sum of all children values
        s   = isempty(list_of_children) ?
                0 :
                mapreduce( k->list_of_children[k].value, +, keys(list_of_children) );
        # Apply escape probability
        esc_prob = (p.model[buffer]-s)/p.model[buffer]
        for symbol in keys(symbols)
            symbols[symbol] *= esc_prob
        end
        # Get each child probability
        for k ∈ keys(list_of_children)
            # println( "        ", k, " -> ", list_of_children[k].value )
            symbols[k] += list_of_children[k].value/p.model[buffer]
        end
    end

    return symbols;
end

function info_string( p::LeZiUpdate{T} ) where {T}
    return @sprintf( "LeZiUpdate" )
end

function unique_string( p::LeZiUpdate{T} ) where {T}
    return @sprintf( "LZUP" )
end

"""
seq: a|aa|b|ab|bb|bba|abc|c|d|dc|ba|aaa

dic{}, phrase[]

a :>
    dic{a}, phrase[]
    model{a:1}
a :>
    dic{a}, phrase[a]
    model{a:1}
a :>
    dic{a,aa}, phrase[]
    model{a:3,aa:1}
b :>
    dic{a,aa,b}, phrase[]
    model{a:3,aa:1,b:1}
a :>
    dic{a,aa,b}, phrase[a]
    model{a:3,aa:1,b:1}
b :>
    dic{a,aa,b,ab}, phrase[]
    model{a:4,aa:1,ab:1,b:2}
b :>
    dic{a,aa,b,ab}, phrase[b]
    model{a:4,aa:1,ab:1,b:2}
"""


text:
# Top-level source for NodalDG
#
#
# Dependencies:
#  FastGaussQuadrature  ( https://github.com/ajt60gaibb/FastGaussQuadrature.jl )
#  SpecialFunctions (for gamma function)
#

module NodalDG

include("meshgen1d.jl")

include("grid1d.jl")

include("lserk.jl")

include("jacobi-basis.jl")

include("grid2d.jl")

end


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

mutable struct NetworkConfigurationDiagnosticResponse <: SwaggerModel
    results::Any # spec type: Union{ Nothing, Vector{NetworkConfigurationDiagnosticResult} } # spec name: results

    function NetworkConfigurationDiagnosticResponse(;results=nothing)
        o = new()
        validate_property(NetworkConfigurationDiagnosticResponse, Symbol("results"), results)
        setfield!(o, Symbol("results"), results)
        o
    end
end # type NetworkConfigurationDiagnosticResponse

const _property_map_NetworkConfigurationDiagnosticResponse = Dict{Symbol,Symbol}(Symbol("results")=>Symbol("results"))
const _property_types_NetworkConfigurationDiagnosticResponse = Dict{Symbol,String}(Symbol("results")=>"Vector{NetworkConfigurationDiagnosticResult}")
Base.propertynames(::Type{ NetworkConfigurationDiagnosticResponse }) = collect(keys(_property_map_NetworkConfigurationDiagnosticResponse))
Swagger.property_type(::Type{ NetworkConfigurationDiagnosticResponse }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_NetworkConfigurationDiagnosticResponse[name]))}
Swagger.field_name(::Type{ NetworkConfigurationDiagnosticResponse }, property_name::Symbol) =  _property_map_NetworkConfigurationDiagnosticResponse[property_name]

function check_required(o::NetworkConfigurationDiagnosticResponse)
    true
end

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


text:
using DrWatson
@quickactivate "StatReth"

# %%
# Sampling from grid approximate posterior
using Plots
using StatsBase
using Distributions
using StatsPlots

# %%
n = 1_000
p_grid = range(0, 1, length = n)
col_
prior = ones(n)

# %%
"""
Fair warning!
I went completely overboard with this one, I got so bad I even started
taking notes in German. Please try to keep your sanity.
"""

# function pi_dist(samples)
#     l, u = quantile(samples, (0.005, 0.995))
#     return u - l
# end

function auswertung(n, k, p_grid, prior)
    likelihood = @. pdf(Binomial(n, p_grid), k)
    # posterior = likelihood .* prior
    # posterior ./= sum(posterior)
    # posterior
    likelihood ./ sum(likelihood)
end

function posterior_über_k(n, p, p_grid, prior)
    gesamt_p_für_p = zeros(length(p_grid))
    for k in 0:n
        p_dass_dieses_k_auftritt = pdf(Binomial(n, p), k)
        p_für_p = auswertung(n, k, p_grid, prior)
        gesamt_p_für_p .+= p_dass_dieses_k_auftritt .* p_für_p
    end
    gesamt_p_für_p
end

function pi_dist(n, p, p_grid, prior)
    post_über_k = posterior_über_k(n, p, p_grid, prior)
    weights = pweights(post_über_k)
    l, u = quantile(p_grid, weights, [0.005, 0.995])
    return u - l
    # ns = 10_000
    # samples = sample(p_grid, weights, ns)
    # samples
end

function über_p(n, p_grid, prior)
    maximum(p_grid) do p
        println(n, " ", p)
        posterior_samples(n, p, p_grid, prior)
        # pi_dist(predpos_samples)
    end
end

# %%
posterior_über_k(5000, 0.7, p_grid, prior) |> plot

let
    likelihood = @. pdf(Binomial(50, p_grid), 35)
    posterior = likelihood .* prior
    posterior ./= sum(posterior)

    plot!(posterior)
end

# %%
let
    predpos_samples = predictiv_posterior(50, 0.7, p_grid, prior)
    density(predpos_samples)

    pi_dist(predpos_samples)
end

über_p(5500, p_grid, prior)

# %%
ns = 100:100:3000
ns = 3100:100:5000
res = zeros(length(ns))
Threads.@threads for i in 1:length(ns)
    n = ns[i]
    res[i] = über_p(n, p_grid, prior)
end
plot(ns, res)

# %%
ns = 100:100:5000
res = [0.36436936936936937, 0.26126126126126126, 0.21421421421421427, 0.18719219219219219, 0.1661711711711712, 0.15415415415415412, 0.1411461461461454, 0.13213713713713715, 0.12612612612612617, 0.11811811811811818, 0.11311311311311306, 0.10710710710710714, 0.10410410410410414, 0.10010010010010012, 0.09709709709709713, 0.09310310310310232, 0.09109109109109109, 0.0880880880880881, 0.08609109109109109, 0.08308308308308315, 0.08208208208208201, 0.08008008008008005, 0.07907907907907907, 0.07707707707707712, 0.07507507507507505, 0.07407407407407407, 0.07207207207207211, 0.07107107107107108, 0.07007007007007005, 0.06906906906906907, 0.06536298120084288, 0.06435885201026825, 0.06335555117095265, 0.062440202049018456, 0.06153473199637327, 0.06067345230968496, 0.05986667635189841, 0.059026361745720124, 0.058315847683361066, 0.05757238632468448, 0.05684749802817202, 0.05620289679737939, 0.05553017612660588, 0.05487561453559159, 0.05430071679577558, 0.053702781658199195, 0.05307944329882003, 0.05256906573675041, 0.05204639141882944, 0.051503807272551605]
scatter(ns, res)


text:
import HoTRG_lemon.SpinModule:
		initialize!,
		isZeroTemperature,
		isClassical,
		getModelname,
		getStates,
		isSymmetricFactorization,
		isPottsModel,
		getHamiltonian,
		getTemperature,
		getExternalfield,
		getEnvParameters,
		getMeasureOperator,
		getFactorW,
		getFactorWp,
		setTemperature!,
		setExternalfield!,
		setEnvParameters!

abstract Lattice{T}

#---
# buildFunction
function buildLattice(latticeinfo::LatticeInfo, spinmodel::SpinModel)
	if isClassical(latticeinfo) != isClassical(spinmodel)
		error("info does not match")
	end
	if isClassical(latticeinfo)
		return buildClassicalLattice(latticeinfo,spinmodel)
	else
		return buildQuantumLattice(latticeinfo,spinmodel)
	end
end

function buildLattice(latticecode::AbstractString, spinmodel::SpinModel)
	latticeinfo = LatticeInfo(latticecode)
	return buildLattice(latticeinfo, spinmodel)
end

#---
# functions

#---
# inherit methods lattice.spinmodel
single_args_non_preserving = [
							:isZeroTemperature,
							:isClassical,
							:isSymmetricFactorization,
							:isPottsModel,
							:getModelname,
							:getStates,
							:getHamiltonian,
							:getTemperature,
							:getExternalfield,
							:getEnvParameters,
							:getFactorW,
							:getFactorWp,
							:getMeasureOperator
							]

for f in single_args_non_preserving
	eval(quote

		function $(f)(lattice::Lattice)
		  return $(f)(lattice.spinmodel)
		end

	end)
end

set_functions = [
				:setTemperature!,
				:setExternalfield!,
				]

for f in set_functions
	eval(quote

		  function $(f){T}(lattice::Lattice{T}, parameter::T)
		  $(f)(lattice.spinmodel, parameter)
			initialize!(lattice)
		end

	end)
end

function setEnvParameters!{T}(lattice::Lattice, temp::T, externalfield::T)
	setEnvParameters!(lattice.spinmodel, temp, externalfield)
	initialize!(lattice)
end

#---
# functions for Lattice
function initialize!{T}(lattice::Lattice{T})
	println("specify lattice!: return random tensor")
	lattice.tensorT = rand(2,2,2,2)
	lattice.tensorTtilde = rand(2,2,2,2)
end

function getTensorT(lattice::Lattice)
	return lattice.tensorT, lattice.tensorTtilde
end

function getTensorT(lattice::Lattice, whichIsTensorT::Int)
	tensorT, tensorTtilde = getTensorT(lattice)
	if whichIsTensorT == 1
		return tensorT, tensorTtilde
	else
		return tensorTtilde, tensorT
	end
end

function setTensorT!(lattice::Lattice, tensorT::Array; tilde = false)
	if tilde
		lattice.tensorTtilde = tensorT
	else
		lattice.tensorT = tensorT
	end
end


text:
"""
    Neurify(max_iter::Int64, tree_search::Symbol)

Neurify combines symbolic reachability analysis with constraint refinement to
minimize over-approximation of the reachable set.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: AbstractPolytope
3. Output: LazySet

# Return
`CounterExampleResult`

# Method
Symbolic reachability analysis and iterative interval refinement (search).
- `max_iter` default `100`.
- `tree_search` default `:DFS` - depth first search.
- `optimizer` default `GLPK.Optimizer`

# Property
Sound but not complete.

# Reference
[S. Wang, K. Pei, J. Whitehouse, J. Yang, and S. Jana,
"Efficient Formal Safety Analysis of Neural Networks,"
*CoRR*, vol. abs/1809.08098, 2018. arXiv: 1809.08098.](https://arxiv.org/pdf/1809.08098.pdf)
[https://github.com/tcwangshiqi-columbia/Neurify](https://github.com/tcwangshiqi-columbia/Neurify)
"""

@with_kw struct Neurify <: Solver
    max_iter::Int64     = 100
    tree_search::Symbol = :DFS
    optimizer = GLPK.Optimizer
end

function solve(solver::Neurify, problem::Problem)
    isbounded(problem.input) || throw(UnboundedInputError("Neurify can only handle bounded input sets."))

    # Because of over-approximation, a split may not bisect the input set.
    # Therefore, the gradient remains unchanged (since input didn't change).
    # And this node will be chosen to split forever.
    # To prevent this, we split each node only once if the gradient of this node hasn't changed.
    # Each element in splits is a tuple (layer_index, node_index, node_gradient).

    nnet, output = problem.network, problem.output
    reach_list = []
    domain = init_symbolic_grad(problem.input)
    splits = Set()
    for i in 1:solver.max_iter
        if i > 1
            domain, splits = select!(reach_list, solver.tree_search)
        end

        reach = forward_network(solver, nnet, domain, collect=true)
        # The first entry is the input set
        popfirst!(reach)
        result, max_violation_con = check_inclusion(solver, nnet, last(reach).sym, output)

        if result.status === :violated
            return result
        elseif result.status === :unknown
            subdomains = constraint_refinement(nnet, reach, max_violation_con, splits)
            for domain in subdomains
                push!(reach_list, (init_symbolic_grad(domain), copy(splits)))
            end

        end
        isempty(reach_list) && return CounterExampleResult(:holds)
    end
    return CounterExampleResult(:unknown)
end

function check_inclusion(solver::Neurify, nnet::Network,
                         reach::SymbolicInterval, output)
    # The output constraint is in the form A*x < b
    # We try to maximize output constraint to find a violated case, or to verify the inclusion.
    # Suppose the output is [1, 0, -1] * x < 2, Then we are maximizing reach.Up[1] * 1 + reach.Low[3] * (-1)

    input_domain = domain(reach)

    model = Model(solver); set_silent(model)
    x = @variable(model, [1:dim(input_domain)])
    add_set_constraint!(model, input_domain, x)

    max_violation = 0.0
    max_violation_con = nothing
    for (i, cons) in enumerate(constraints_list(output))
        # NOTE can be taken out of the loop, but maybe there's no advantage
        # NOTE max.(M, 0) * U  + ... is a common operation, and maybe should get a name. It's also an "interval map".
        a, b = cons.a, cons.b
        c = max.(a, 0)'*reach.Up + min.(a, 0)'*reach.Low

        @objective(model, Max, c * [x; 1] - b)
        optimize!(model)

        if termination_status(model) == OPTIMAL
            if compute_output(nnet, value(x)) ∉ output
                return CounterExampleResult(:violated, value(x)), nothing
            end

            viol = objective_value(model)
            if viol > max_violation
                max_violation = viol
                max_violation_con = a
            end
        else
            # TODO can we be more descriptive?
            error("No solution, please check the problem definition.")
        end

    end

    if max_violation > 0.0
        return CounterExampleResult(:unknown), max_violation_con
    else
        return CounterExampleResult(:holds), nothing
    end
end

function constraint_refinement(nnet::Network,
                               reach::Vector{<:SymbolicIntervalGradient},
                               max_violation_con::AbstractVector{Float64},
                               splits)

    i, j, influence = get_max_nodewise_influence(nnet, reach, max_violation_con, splits)
    # We can generate three more constraints
    # Symbolic representation of node i j is Low[i][j,:] and Up[i][j,:]
    aL, bL = reach[i].sym.Low[j, 1:end-1], reach[i].sym.Low[j, end]
    aU, bU = reach[i].sym.Up[j, 1:end-1], reach[i].sym.Up[j, end]

    # custom intersection function that doesn't do constraint pruning
    ∩ = (set, lc) -> HPolytope([constraints_list(set); lc])

    subsets = [domain(reach[1])] # all the reaches have the same domain, so we can pick [1]

    # If either of the normal vectors is the 0-vector, we must skip it.
    # It cannot be used to create a halfspace constraint.
    # NOTE: how can this come about, and does it mean anything?
    if !iszero(aL)
        subsets = subsets .∩ [HalfSpace(aL, -bL), HalfSpace(aL, -bL), HalfSpace(-aL, bL)]
    end
    if !iszero(aU)
        subsets = subsets .∩ [HalfSpace(aU, -bU), HalfSpace(-aU, bU), HalfSpace(-aU, bU)]
    end
    return filter(!isempty, subsets)
end

function get_max_nodewise_influence(nnet::Network,
                                    reach::Vector{<:SymbolicIntervalGradient},
                                    max_violation_con::AbstractVector{Float64},
                                    splits)

    LΛ, UΛ = reach[end].LΛ, reach[end].UΛ
    is_ambiguous_activation(i, j) = (0 < LΛ[i][j] < 1) || (0 < UΛ[i][j] < 1)

    # We want to find the node with the largest influence
    # Influence is defined as gradient * interval width
    # The gradient is with respect to a loss defined by the most violated constraint.
    LG = UG = max_violation_con
    i_max, j_max, influence_max = 0, 0, -Inf

    # Backpropagation to calculate the node-wise gradient
    for i in reverse(1:length(nnet.layers))
        layer = nnet.layers[i]
        sym = reach[i].sym
        if layer.activation isa ReLU
            for j in 1:n_nodes(layer)
                if is_ambiguous_activation(i, j)
                    # taking `influence = max_gradient * reach.r[i][j]*k` would be
                    # different from original paper, but can improve the split efficiency.
                    # where `k = n-i+1`, i.e. counts up from 1 as you go back in layers.

                    # radius wrt to the jth node/hidden dimension
                    r = radius(sym, j)
                    influence = max(abs(LG[j]), abs(UG[j])) * r
                    if influence >= influence_max && (i, j, influence) ∉ splits
                        i_max, j_max, influence_max = i, j, influence
                    end
                end
            end
        end

        LG_hat = max.(LG, 0.0) .* LΛ[i] .+ min.(LG, 0.0) .* UΛ[i]
        UG_hat = min.(UG, 0.0) .* LΛ[i] .+ max.(UG, 0.0) .* UΛ[i]

        LG, UG = interval_map(layer.weights', LG_hat, UG_hat)
    end

    # NOTE can omit this line in the paper version
    (i_max == 0 || j_max == 0) && error("Can not find valid node to split")

    push!(splits, (i_max, j_max, influence_max))

    return (i_max, j_max, influence_max)
end

# Symbolic forward_linear
function forward_linear(solver::Neurify, L::Layer, input::SymbolicIntervalGradient)
    output_Low, output_Up = interval_map(L.weights, input.sym.Low, input.sym.Up)
    output_Up[:, end] += L.bias
    output_Low[:, end] += L.bias
    sym = SymbolicInterval(output_Low, output_Up, domain(input))
    return SymbolicIntervalGradient(sym, input.LΛ, input.UΛ)
end

# Symbolic forward_act
function forward_act(solver::Neurify, L::Layer{ReLU}, input::SymbolicIntervalGradient)
    n_node = n_nodes(L)
    output_Low, output_Up = copy(input.sym.Low), copy(input.sym.Up)
    LΛᵢ, UΛᵢ = zeros(n_node), ones(n_node)
    # Symbolic linear relaxation
    # This is different from ReluVal
    for j in 1:n_node
        up_low, up_up = bounds(upper(input), j)
        low_low, low_up = bounds(lower(input), j)

        up_slope = relaxed_relu_gradient(up_low, up_up)
        low_slope = relaxed_relu_gradient(low_low, low_up)

        output_Up[j, :] .*= up_slope
        output_Up[j, end] += up_slope * max(-up_low, 0)

        output_Low[j, :] .*= low_slope

        LΛᵢ[j], UΛᵢ[j] = low_slope, up_slope
    end
    sym = SymbolicInterval(output_Low, output_Up, domain(input))
    LΛ = push!(input.LΛ, LΛᵢ)
    UΛ = push!(input.UΛ, UΛᵢ)
    return SymbolicIntervalGradient(sym, LΛ, UΛ)
end

function forward_act(solver::Neurify, L::Layer{Id}, input::SymbolicIntervalGradient)
    n_node = n_nodes(L)
    LΛ = push!(input.LΛ, ones(n_node))
    UΛ = push!(input.UΛ, ones(n_node))
    return SymbolicIntervalGradient(input.sym, LΛ, UΛ)
end


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

name = "LCIO"
version = v"02.16.1"

# Collection of sources required to build LCIO
sources = [
    GitSource("https://github.com/iLCSoft/LCIO.git", "8f9e86b93b7d5d83221fabb872ed7e82f1638476"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/LCIO*/
mkdir build && cd build

cmake .. -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release
cmake --build . --target install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
filter!(!Sys.isfreebsd, platforms)
filter!(!Sys.iswindows, platforms)
filter!(p -> arch(p) != "armv7l", platforms)
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("liblcio", :liblcio),
    LibraryProduct("libsio", :libsio)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Zlib_jll"),
]

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

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
@doc Markdown.doc"""
    lift(a::padic) -> fmpz

Returns the positive canonical representative in $\mathbb{Z}$. $a$ needs
to be integral.
"""
function lift(a::padic)
  b = fmpz()
  R = parent(a)
  ccall((:padic_get_fmpz, libflint), Nothing, (Ref{fmpz}, Ref{padic}, Ref{FlintPadicField}), b, a, R)
  return b
end

function _lift(a::padic)
  R = parent(a)
  v = valuation(a)
  if v >= 0
    return fmpq(lift(a))
  else
    m = prime(R)^-v
    return fmpq(lift(m * a))//m
  end
end

function Base.setprecision(f::Generic.Poly{padic}, N::Int)
  f = map_coefficients(x->setprecision(x, N), f, parent = parent(f))
  return f
end

function setprecision!(f::Generic.Poly{padic}, N::Int)
  for i=1:length(f)
    f.coeffs[i] = setprecision!(f.coeffs[i], N)
  end
  return f
end

"""
The log of `1-x`, x needs to have a valuation >1
"""
function my_log_one_minus_inner(x::fmpz, pr::Int, v::Int, p)
  #need N s.th. Nv-log_p(N) >= pr
  #as a function of N, this has a min at log(p)/v
  #the N needs to > pr/v + s.th. small
  lp = log(p)
  N = div(pr, v) + 1
  while N*v-log(N)/lp < pr
    @show N += 1
  end
  return -x*my_eval(map(i->fmpq(1, i), 1:N), fmpq(x))
end

function my_log_one_minus(x::padic)
  v = valuation(x)
  lg = parent(x)(0)
  le = 1
  #write x as (1-py_1)(1-p^2y_2)(1-p^4y_3)...
  #with y_i small , i.e. |y_i| <= p^(2^(i-1))
  pp = prime(parent(x), 2)
  X = 1-x
  while true
    y = lift(1-X) % pp
    lg += parent(x)(my_log_one_minus_inner(y, precision(x), le, prime(parent(x))))
    X = X*inv(parent(x)(1-y))
    pp *= pp
    le *= 2
    if isone(X) || le >= precision(x)
      return lg
    end
    @assert valuation(X-1) >= le
  end

end

function my_log_one_minus_inner(x::Generic.Res{fmpq_poly}, pr::Int, v::Int, p)
  #need N s.th. Nv-log_p(N) >= pr
  #as a function of N, this has a min at log(p)/v
  #the N needs to > pr/v + s.th. small
  lp = log(p)
  N = div(pr, v) + 1
  while N*v-log(N)/lp < pr
    @show N += 1
  end
  R = parent(x)
  return -x*my_eval(map(i->R(fmpq(1, i)), 1:N), x)
end

function my_log_one_minus(x::qadic)
  v = valuation(x)
  lg = parent(x)(0)
  le = 1
  #write x as (1-py_1)(1-p^2y_2)(1-p^4y_3)...
  #with y_i small , i.e. |y_i| <= p^(2^(i-1))
  pp = prime(parent(x))^2
  X = 1-x
  R, _ = PolynomialRing(QQ, cached = false)
  S = ResidueRing(R, map_coefficients(x->QQ(lift(x)), defining_polynomial(parent(x)), parent = R))
  while true
    Y = 1-X
    y = S(R([lift(coeff(Y, i)) % pp for i=0:length(Y)]))
    lg += parent(x)(my_log_one_minus_inner(y, precision(x), le, prime(parent(x))).data)
    X = X*inv(parent(x)(1-y.data))
    pp *= pp
    le *= 2
    if isone(X) || le >= precision(x)
      return lg
    end
    @assert valuation(X-1) >= le
  end

end


text:
####################################################################
####################### alpha reconstruction #######################
####################################################################
using PyPlot

nxy = 256
alpharec = zeros(Float64,nxy,nxy,3)
alpharecsky = zeros(Float64,nxy,nxy,3)
alpharecdata = zeros(Float64,nxy,nxy,3)
nu0 = (nu[end]+nu[1])/2

N = 15
ν = nu./nu0
A = [ones(N) log10(ν) log10(ν).^2]
for i = 1:nxy, j = 1:nxy

    if minimum(squeeze(squeeze(x[i,j,:],1),1)) <= 0
    alpharec[i,j,:] = 0

    elseif minimum(squeeze(squeeze(x[i,j,:],1),1)) > 0
    b = squeeze(squeeze(log10(x[i,j,:]),1),1)
    c = A\b
    alpharec[i,j,1] = c[1]
    alpharec[i,j,2] = c[2]
    alpharec[i,j,3] = c[3]
    end
end

for i = 1:nxy, j = 1:nxy

    if minimum(squeeze(squeeze(sky[i,j,:],1),1)) <= 0
    alpharecsky[i,j,:] = 0

    elseif minimum(squeeze(squeeze(sky[i,j,:],1),1)) > 0
    b = squeeze(squeeze(log10(sky[i,j,:]),1),1)
    c = A\b
    alpharecsky[i,j,1] = c[1]
    alpharecsky[i,j,2] = c[2]
    alpharecsky[i,j,3] = c[3]
    end
end

for i = 1:nxy, j = 1:nxy

    if minimum(squeeze(squeeze(mydata[i,j,:],1),1)) <= 0
    alpharecdata[i,j,:] = 0

elseif minimum(squeeze(squeeze(mydata[i,j,:],1),1)) > 0
    b = squeeze(squeeze(log10(mydata[i,j,:]),1),1)
    c = A\b
    alpharecdata[i,j,1] = c[1]
    alpharecdata[i,j,2] = c[2]
    alpharecdata[i,j,3] = c[3]
    end
end


text:
function test_execute()
    execute(TEST_CONFIG, `ls -al`; name="myproject")
    execute(TEST_CONFIG, `julia -v`; dir=joinpath(TEST_TMP_DIR, "test-playground"))
end

test_execute()


text:
module NeuralFieldEq

import Base: show

using FFTW
using Distributions: mean
using LinearAlgebra: mul!, ldiv!

# Module containing auxiliary functions
include("AuxFunctions.jl")
using .AuxFunctions

# Define the structures types for probNFE and solveNFE
include("NFEStruct.jl")

export Input1D, Input2D  # Public structures to wrap up the inputs to probNFE
export probNFE, solveNFE # Public functions to compute NFE
export SolveOutSto1D, SolveOutSto2D # Structures to print output of solveNFE

# Interface functions
include("probNFE.jl")
include("solveNFE1D.jl")
include("solveNFE2D.jl")

# Format outputs for probSNFE and solveSNFE
include("printOut.jl")

end


text:
int_rules_7_2_3 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 7.2.3*(d+e*x^2)^p*(a+b*arccosh(c*x))^n =#
    @apply_utils Antiderivative((~d1 + ~(e1') * ~x) ^ ~(p') * (~d2 + ~(e2') * ~x) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative((~d1 * ~d2 + ~e1 * ~e2 * (~x) ^ 2) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~n], ~x) && (EqQ(~d2 * ~e1 + ~d1 * ~e2, 0) && IntegerQ(~p))
    #= Antiderivative(((~a')+(~b')*acosh((~c')*(~x)))^(~n')/sqrt((~d)+(~e')*(~x)^2),~x) :=  1/c*Simp(sqrt(1+c*x)*sqrt(-1+c*x)/sqrt(d+e*x^2))*Subst(Antiderivative((a+b*x)^n, x),x,acosh(c*x)) <-- FreeQ([a,b,c,d,e,n],x) && EqQ(c^2*d+e,0) =#
    @apply_utils Antiderivative(1 / (sqrt(~d + ~(e') * (~x) ^ 2) * (~(a') + ~(b') * acosh(~(c') * ~x))), ~x) => (1 / (~b * ~c)) * Simp((sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x)) / sqrt(~d + ~e * (~x) ^ 2)) * log(~a + ~b * acosh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && EqQ((~c) ^ 2 * ~d + ~e, 0)
    @apply_utils Antiderivative(1 / (sqrt(~d1 + ~(e1') * ~x) * sqrt(~d2 + ~(e2') * ~x) * (~(a') + ~(b') * acosh(~(c') * ~x))), ~x) => (1 / (~b * ~c)) * Simp(sqrt(1 + ~c * ~x) / sqrt(~d1 + ~e1 * ~x)) * Simp(sqrt(-1 + ~c * ~x) / sqrt(~d2 + ~e2 * ~x)) * log(~a + ~b * acosh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2], ~x) && (EqQ(~e1, ~c * ~d1) && EqQ(~e2, -(~c) * ~d2))
    @apply_utils Antiderivative((~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n') / sqrt(~d + ~(e') * (~x) ^ 2), ~x) => (1 / (~b * ~c * (~n + 1))) * Simp((sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x)) / sqrt(~d + ~e * (~x) ^ 2)) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && NeQ(~n, -1))
    @apply_utils Antiderivative((~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n') / (sqrt(~d1 + ~(e1') * ~x) * sqrt(~d2 + ~(e2') * ~x)), ~x) => (1 / (~b * ~c * (~n + 1))) * Simp(sqrt(1 + ~c * ~x) / sqrt(~d1 + ~e1 * ~x)) * Simp(sqrt(-1 + ~c * ~x) / sqrt(~d2 + ~e2 * ~x)) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~n], ~x) && (EqQ(~e1, ~c * ~d1) && (EqQ(~e2, -(~c) * ~d2) && NeQ(~n, -1)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * acosh(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(SimplifyIntegrand(u / (sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IGtQ(~p, 0))
    @apply_utils Antiderivative(sqrt(~d + ~(e') * (~x) ^ 2) * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => ((~x * sqrt(~d + ~e * (~x) ^ 2) * (~a + ~b * acosh(~c * ~x)) ^ ~n) / 2 - ((~b * ~c * ~n) / 2) * Simp(sqrt(~d + ~e * (~x) ^ 2) / (sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x))) * Antiderivative(~x * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1), ~x)) - (1 / 2) * Simp(sqrt(~d + ~e * (~x) ^ 2) / (sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x))) * Antiderivative((~a + ~b * acosh(~c * ~x)) ^ ~n / (sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && GtQ(~n, 0))
    @apply_utils Antiderivative(sqrt(~d1 + ~(e1') * ~x) * sqrt(~d2 + ~(e2') * ~x) * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => ((~x * sqrt(~d1 + ~e1 * ~x) * sqrt(~d2 + ~e2 * ~x) * (~a + ~b * acosh(~c * ~x)) ^ ~n) / 2 - ((~b * ~c * ~n) / 2) * Simp(sqrt(~d1 + ~e1 * ~x) / sqrt(1 + ~c * ~x)) * Simp(sqrt(~d2 + ~e2 * ~x) / sqrt(-1 + ~c * ~x)) * Antiderivative(~x * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1), ~x)) - (1 / 2) * Simp(sqrt(~d1 + ~e1 * ~x) / sqrt(1 + ~c * ~x)) * Simp(sqrt(~d2 + ~e2 * ~x) / sqrt(-1 + ~c * ~x)) * Antiderivative((~a + ~b * acosh(~c * ~x)) ^ ~n / (sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x)), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2], ~x) && (EqQ(~e1, ~c * ~d1) && (EqQ(~e2, -(~c) * ~d2) && GtQ(~n, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => ((~x * (~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n) / (2 * ~p + 1) + ((2 * ~d * ~p) / (2 * ~p + 1)) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p - 1) * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x)) - ((~b * ~c * ~n) / (2 * ~p + 1)) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / ((1 + ~c * ~x) ^ ~p * (-1 + ~c * ~x) ^ ~p)) * Antiderivative(~x * (1 + ~c * ~x) ^ (~p - 1 / 2) * (-1 + ~c * ~x) ^ (~p - 1 / 2) * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (GtQ(~n, 0) && GtQ(~p, 0)))
    @apply_utils Antiderivative((~d1 + ~(e1') * ~x) ^ ~(p') * (~d2 + ~(e2') * ~x) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => ((~x * (~d1 + ~e1 * ~x) ^ ~p * (~d2 + ~e2 * ~x) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n) / (2 * ~p + 1) + ((2 * ~d1 * ~d2 * ~p) / (2 * ~p + 1)) * Antiderivative((~d1 + ~e1 * ~x) ^ (~p - 1) * (~d2 + ~e2 * ~x) ^ (~p - 1) * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x)) - ((~b * ~c * ~n) / (2 * ~p + 1)) * Simp((~d1 + ~e1 * ~x) ^ ~p / (1 + ~c * ~x) ^ ~p) * Simp((~d2 + ~e2 * ~x) ^ ~p / (-1 + ~c * ~x) ^ ~p) * Antiderivative(~x * (1 + ~c * ~x) ^ (~p - 1 / 2) * (-1 + ~c * ~x) ^ (~p - 1 / 2) * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2], ~x) && (EqQ(~e1, ~c * ~d1) && (EqQ(~e2, -(~c) * ~d2) && (GtQ(~n, 0) && GtQ(~p, 0))))
    @apply_utils Antiderivative((~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n') / (~d + ~(e') * (~x) ^ 2) ^ (3 / 2), ~x) => (~x * (~a + ~b * acosh(~c * ~x)) ^ ~n) / (~d * sqrt(~d + ~e * (~x) ^ 2)) + ((~b * ~c * ~n) / ~d) * Simp((sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x)) / sqrt(~d + ~e * (~x) ^ 2)) * Antiderivative((~x * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1)) / (1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && GtQ(~n, 0))
    @apply_utils Antiderivative((~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n') / ((~d1 + ~(e1') * ~x) ^ (3 / 2) * (~d2 + ~(e2') * ~x) ^ (3 / 2)), ~x) => (~x * (~a + ~b * acosh(~c * ~x)) ^ ~n) / (~d1 * ~d2 * sqrt(~d1 + ~e1 * ~x) * sqrt(~d2 + ~e2 * ~x)) + ((~b * ~c * ~n) / (~d1 * ~d2)) * Simp(sqrt(1 + ~c * ~x) / sqrt(~d1 + ~e1 * ~x)) * Simp(sqrt(-1 + ~c * ~x) / sqrt(~d2 + ~e2 * ~x)) * Antiderivative((~x * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1)) / (1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2], ~x) && (EqQ(~e1, ~c * ~d1) && (EqQ(~e2, -(~c) * ~d2) && GtQ(~n, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => ((-(~x) * (~d + ~e * (~x) ^ 2) ^ (~p + 1) * (~a + ~b * acosh(~c * ~x)) ^ ~n) / (2 * ~d * (~p + 1)) + ((2 * ~p + 3) / (2 * ~d * (~p + 1))) * Antiderivative((~d + ~e * (~x) ^ 2) ^ (~p + 1) * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x)) - ((~b * ~c * ~n) / (2 * (~p + 1))) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / ((1 + ~c * ~x) ^ ~p * (-1 + ~c * ~x) ^ ~p)) * Antiderivative(~x * (1 + ~c * ~x) ^ (~p + 1 / 2) * (-1 + ~c * ~x) ^ (~p + 1 / 2) * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (GtQ(~n, 0) && (LtQ(~p, -1) && NeQ(~p, -3 / 2))))
    @apply_utils Antiderivative((~d1 + ~(e1') * ~x) ^ ~p * (~d2 + ~(e2') * ~x) ^ ~p * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => ((-(~x) * (~d1 + ~e1 * ~x) ^ (~p + 1) * (~d2 + ~e2 * ~x) ^ (~p + 1) * (~a + ~b * acosh(~c * ~x)) ^ ~n) / (2 * ~d1 * ~d2 * (~p + 1)) + ((2 * ~p + 3) / (2 * ~d1 * ~d2 * (~p + 1))) * Antiderivative((~d1 + ~e1 * ~x) ^ (~p + 1) * (~d2 + ~e2 * ~x) ^ (~p + 1) * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x)) - ((~b * ~c * ~n) / (2 * (~p + 1))) * Simp((~d1 + ~e1 * ~x) ^ ~p / (1 + ~c * ~x) ^ ~p) * Simp((~d2 + ~e2 * ~x) ^ ~p / (-1 + ~c * ~x) ^ ~p) * Antiderivative(~x * (1 + ~c * ~x) ^ (~p + 1 / 2) * (-1 + ~c * ~x) ^ (~p + 1 / 2) * (~a + ~b * acosh(~c * ~x)) ^ (~n - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2], ~x) && (EqQ(~e1, ~c * ~d1) && (EqQ(~e2, -(~c) * ~d2) && (GtQ(~n, 0) && (LtQ(~p, -1) && NeQ(~p, -3 / 2)))))
    @apply_utils Antiderivative((~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n') / (~d + ~(e') * (~x) ^ 2), ~x) => (-1 / (~c * ~d)) * Subst(Antiderivative((~a + ~b * ~x) ^ ~n * csch(~x), ~x), ~x, acosh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IGtQ(~n, 0))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~n, ~x) => (Simp(sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x) * (~d + ~e * (~x) ^ 2) ^ ~p) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) - ((~c * (2 * ~p + 1)) / (~b * (~n + 1))) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / ((1 + ~c * ~x) ^ ~p * (-1 + ~c * ~x) ^ ~p)) * Antiderivative(~x * (1 + ~c * ~x) ^ (~p - 1 / 2) * (-1 + ~c * ~x) ^ (~p - 1 / 2) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~p], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && (LtQ(~n, -1) && IntegerQ(2 * ~p)))
    @apply_utils Antiderivative((~d1 + ~(e1') * ~x) ^ ~(p') * (~d2 + ~(e2') * ~x) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~n, ~x) => (sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x) * (~d1 + ~e1 * ~x) ^ ~p * (~d2 + ~e2 * ~x) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) - ((~c * (2 * ~p + 1)) / (~b * (~n + 1))) * Simp((~d1 + ~e1 * ~x) ^ ~p / (1 + ~c * ~x) ^ ~p) * Simp((~d2 + ~e2 * ~x) ^ ~p / (-1 + ~c * ~x) ^ ~p) * Antiderivative(~x * (-1 + (~c) ^ 2 * (~x) ^ 2) ^ (~p - 1 / 2) * (~a + ~b * acosh(~c * ~x)) ^ (~n + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~p], ~x) && (EqQ(~e1, ~c * ~d1) && (EqQ(~e2, -(~c) * ~d2) && (LtQ(~n, -1) && IntegerQ(~p + 1 / 2))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => (1 / (~b * ~c)) * Simp((~d + ~e * (~x) ^ 2) ^ ~p / ((1 + ~c * ~x) ^ ~p * (-1 + ~c * ~x) ^ ~p)) * Subst(Antiderivative((~x) ^ ~n * sinh(-(~a) / ~b + ~x / ~b) ^ (2 * ~p + 1), ~x), ~x, ~a + ~b * acosh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (EqQ((~c) ^ 2 * ~d + ~e, 0) && IGtQ(2 * ~p, 0))
    @apply_utils Antiderivative((~d1 + ~(e1') * ~x) ^ ~(p') * (~d2 + ~(e2') * ~x) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => (1 / (~b * ~c)) * Simp((~d1 + ~e1 * ~x) ^ ~p / (1 + ~c * ~x) ^ ~p) * Simp((~d2 + ~e2 * ~x) ^ ~p / (-1 + ~c * ~x) ^ ~p) * Subst(Antiderivative((~x) ^ ~n * sinh(-(~a) / ~b + ~x / ~b) ^ (2 * ~p + 1), ~x), ~x, ~a + ~b * acosh(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~n], ~x) && (EqQ(~e1, ~c * ~d1) && (EqQ(~e2, -(~c) * ~d2) && IGtQ(2 * ~p, 0)))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)), ~x) => With([u = IntHide((~d + ~e * (~x) ^ 2) ^ ~p, ~x)], Dist(~a + ~b * acosh(~c * ~x), u, ~x) - ~b * ~c * Antiderivative(SimplifyIntegrand(u / (sqrt(1 + ~c * ~x) * sqrt(-1 + ~c * ~x)), ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~d, ~e], ~x) && (NeQ((~c) ^ 2 * ~d + ~e, 0) && (IGtQ(~p, 0) || ILtQ(~p + 1 / 2, 0)))
    #= Antiderivative(((~d)+(~e')*(~x)^2)^(~p')*((~a')+(~b')*acosh((~c')*(~x)))^(~n),~x) :=  1/(b*c^(2*p+1))*Subst(Antiderivative(x^n*(c^2*d+e*cosh(-a/b+x/b)^2)^p*sinh(-a/b+ x/b),x),x,a+b*acosh(c*x)) <-- FreeQ([a,b,c,d,e,n],x) && IGtQ(p,0) =#
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Antiderivative(ExpandIntegrand((~a + ~b * acosh(~c * ~x)) ^ ~n, (~d + ~e * (~x) ^ 2) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n], ~x) && (NeQ((~c) ^ 2 * ~d + ~e, 0) && (IntegerQ(~p) && (~p > 0 || IGtQ(~n, 0))))
    @apply_utils Antiderivative((~d + ~(e') * (~x) ^ 2) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable((~d + ~e * (~x) ^ 2) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d, ~e, ~n, ~p], ~x)
    @apply_utils Antiderivative((~d1 + ~(e1') * ~x) ^ ~(p') * (~d2 + ~(e2') * ~x) ^ ~(p') * (~(a') + ~(b') * acosh(~(c') * ~x)) ^ ~(n'), ~x) => Unintegrable((~d1 + ~e1 * ~x) ^ ~p * (~d2 + ~e2 * ~x) ^ ~p * (~a + ~b * acosh(~c * ~x)) ^ ~n, ~x) <-- FreeQ([~a, ~b, ~c, ~d1, ~e1, ~d2, ~e2, ~n, ~p], ~x)
end

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.2 Inverse hyperbolic cosine/7.2.3 (d+e x^2)^p (a+b arccosh(c x))^n.jl


text:
#This file was used to generate the .JLD file used in the generate plots code

using MAT

using Plots
using Measures
using Flux
using DifferentialEquations
using DiffEqFlux
using LaTeXStrings
using Random
using Polynomials

vars = matread("C:/Users/Raj/Desktop/2.168/COVID-19-master_latest/COVID-19-master/csse_covid_19_data/csse_covid_19_daily_reports/Rise_Italy_Track.mat")

Infected = vars["Italy_Infected_All"]
Recovered = vars["Italy_Recovered_All"]
Dead= vars["Italy_Dead_All"]
Time = vars["Italy_Time"]

Infected = Infected - Recovered - Dead

#Mobility
#Add apple mobility index data here
vars = matread("C:/Users/Raj/Desktop/Julia Lab Internship/COVID Modelling 1/Italy_Mobility_20June.mat")
Mobility = vars["Mobility"]
Mobility = hcat(Mobility, Mobility[end])
Mobility = Mobility ./ 100

Random.seed!(150)

#Define neural network architecture for contact and recovery rate
ann1 = Chain(Dense(1,10,relu), Dense(10, 10), Dense(10,1))
p1,re = Flux.destructure(ann1)

ann2 = Chain(Dense(3, 10,relu), Dense(10,1))
p1n,re2 = Flux.destructure(ann2)
p3 = [p1; p1n]
ps = Flux.params(p3)

##Define model. Beta (contact rate) is parametrized by the mobility data while recovery rate is not.
function QSIR(du, u, p, t)
    index = Int(floor(t))
    NN1 = abs(re(p[1:141])([Mobility[index]])[1])
    un = [u[1]; u[2]; u[3]]
    NN2 = abs(re2(p[142:192])(un)[1])
    du[1]=  - 0.0003*NN1*u[1]*(u[2])/u0[1]
    du[2] = 0.0003*NN1*u[1]*(u[2])/u0[1] - NN2*u[2]/u0[1]
    du[3] = NN2*u[2]/u0[1]
end

u0 = Float64[60000000.0, 593, 62]
tspan = (1, 95.0)
datasize = 95;
#p = param([0.18])
#p = ([0.2, 0.013, 2])

prob = ODEProblem(QSIR, u0, tspan, p3)
t = range(tspan[1],tspan[2],length=datasize)

###Check if solution is calculated
sol = Array(concrete_solve(prob, Tsit5(),u0, p3, saveat=t))

function predict_adjoint() # Our 1-layer neural network
  Array(concrete_solve(prob,Tsit5(),u0,p3,saveat=t))
end

function loss_adjoint()
 prediction = predict_adjoint()
 #loss = sum(abs2, log.(abs.(Infected[1:end])) .- log.(abs.(prediction[2, :]))) + 1*(sum(abs2, log.(abs.(Recovered[1:end] + Dead[1:end]) ) .- log.(abs.(prediction[3, :] ))))
 loss = sum(abs2, log.(abs.(Infected[1:end])) .- log.(abs.(prediction[2, :]))) + (sum(abs2, log.(abs.(Recovered[1:end] + Dead[1:end]) ) .- log.(abs.(prediction[3, :] ))))
end

Loss = []
P3  =[]
#anim = Animation()
datan = Iterators.repeated((), 20000)
opt = ADAM(0.01)
cb = function()
  display(loss_adjoint())
  global Loss = append!(Loss, loss_adjoint())
  global P3 = append!(P3, p3)

end

cb()

#Number of iterations required = 40000. Can be more depending on model architecture

Flux.train!(loss_adjoint, ps, datan, opt, cb = cb)

####Loss is stagnating around 50! Need to reduce it a lot. How?
L = findmin(Loss)
idx = L[2]
idx1 = (idx-1)*192 +1
idx2 = idx*192
p3n = P3[idx1: idx2]

prediction = Array(concrete_solve(prob,Tsit5(),u0,p3n,saveat=t))

#Compare preeiction with data
bar(Time[1:end], Infected[1:end], alpha = 0.4, xaxis = "Days post 500 infected", yaxis = "Italy: Number of cases", label = "Data: Infected", legend = :topleft, framestyle = :box, left_margin = 5mm, color = :red)
plot!(t, prediction[2, :], xaxis = "Days post 500 infected", yaxis = "Italy: Number of cases", label = "Prediction", legend = :topright, framestyle = :box, left_margin = 5mm, bottom_margin = 5mm, top_margin = 5mm,  grid = :off, color = :red, linewidth  = 3, ylims = (0, 20000), foreground_color_legend = nothing, background_color_legend = nothing, yguidefontsize = 14, xguidefontsize = 14,  xtickfont = font(12, "TimesNewRoman"), ytickfont = font(12, "TimesNewRoman"), legendfontsize = 12)
bar!(Time[1:end], Recovered[1:end] + Dead[1:end], alpha = 0.4, xaxis = "Days post 500 infected", yaxis = "Italy: Number of cases", label = "Data: Recovered", legend = :topleft, framestyle = :box, left_margin = 5mm, color = :blue)
plot!(t, prediction[3, :], ylims = (-0.05*maximum(Recovered + Dead),2.5*maximum(Recovered + Dead)), right_margin = 5mm, xaxis = "Days post 500 infected", yaxis = "Italy: Number of cases", label = "Prediction ", legend = :topleft, framestyle = :box, left_margin = 5mm, bottom_margin =5mm, top_margin = 5mm, grid = :off, color = :blue, linewidth  = 3, foreground_color_legend = nothing, background_color_legend = nothing,  yguidefontsize = 14, xguidefontsize = 14,  xtickfont = font(12, "TimesNewRoman"), ytickfont = font(12, "TimesNewRoman"), legendfontsize = 12)

##DO FORECASTING FOR MORE TIME POINTS (MAYBE 10 -14 days)


text:
import RData
import RDatasets


text:
using Printf
using Dates
using JLD2
using ClimateMachine.VTK: writevtk, writepvtu
using ClimateMachine.Callbacks
using ClimateMachine.GenericCallbacks

abstract type AbstractCallback end

struct Default <: AbstractCallback end
struct Info <: AbstractCallback end
struct StateCheck{T} <: AbstractCallback
    number_of_checks::T
end

Base.@kwdef struct JLD2State{T, V, B} <: AbstractCallback
    iteration::T
    filepath::V
    overwrite::B = true
end

Base.@kwdef struct VTKOutput{T, V, B, S} <: AbstractCallback
    iteration::T
    outputdir::V
    overwrite::B = true
    number_sample_points::S
end

struct ExponentialFiltering <: AbstractCallback end

function create_callbacks(simulation, odesolver)
    callbacks = simulation.callbacks

    if isempty(callbacks)
        return ()
    else
        cbvector = [create_callback(callback, simulation, odesolver) for callback in callbacks]
        return tuple(cbvector...)
    end
end

function create_callback(::Default, simulation::Simulation, odesolver)
    cb_info = create_callback(Info(), simulation, odesolver)
    cb_state_check = create_callback(StateCheck(10), simulation, odesolver)

    return (cb_info, cb_state_check)
end

function create_callback(::Info, simulation::Simulation, odesolver)
    Q = simulation.state
    timeend = simulation.time.finish
    mpicomm = MPI.COMM_WORLD

    starttime = Ref(now())
    cbinfo = ClimateMachine.GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (s = false)
        if s
            starttime[] = now()
        else
            energy = norm(Q)
            @info @sprintf(
                """Update
                simtime = %8.2f / %8.2f
                runtime = %s
                norm(Q) = %.16e""",
                ClimateMachine.ODESolvers.gettime(odesolver),
                timeend,
                Dates.format(convert(Dates.DateTime, Dates.now() - starttime[]), Dates.dateformat"HH:MM:SS"),
                energy
            )

            if isnan(energy)
                error("NaNs")
            end
        end
    end

    return cbinfo
end

function create_callback(callback::StateCheck, simulation::Simulation, _...)
    sim_length = simulation.time.finish - simulation.time.start
    timestep = simulation.timestepper.timestep
    nChecks = callback.number_of_checks

    nt_freq = floor(Int, sim_length / timestep / nChecks)

    cbcs_dg = ClimateMachine.StateCheck.sccreate([(simulation.state, "state")], nt_freq)

    return cbcs_dg
end

function create_callback(output::JLD2State, simulation::Simulation, odesolver)
    # Initialize output
    output.overwrite && isfile(output.filepath) && rm(output.filepath; force = output.overwrite)

    Q = simulation.state
    mpicomm = MPI.COMM_WORLD
    iteration = output.iteration

    steps = ClimateMachine.ODESolvers.getsteps(odesolver)
    time = ClimateMachine.ODESolvers.gettime(odesolver)

    file = jldopen(output.filepath, "a+")
    JLD2.Group(file, "state")
    JLD2.Group(file, "time")
    file["state"][string(steps)] = Array(Q)
    file["time"][string(steps)] = time
    close(file)

    jldcallback = ClimateMachine.GenericCallbacks.EveryXSimulationSteps(iteration) do (s = false)
        steps = ClimateMachine.ODESolvers.getsteps(odesolver)
        time = ClimateMachine.ODESolvers.gettime(odesolver)
        @info steps, time
        file = jldopen(output.filepath, "a+")
        file["state"][string(steps)] = Array(Q)
        file["time"][string(steps)] = time
        close(file)
        return nothing
    end
end

function create_callback(output::JLD2State, simulation::Simulation, odesolver)
    # Initialize output
    output.overwrite && isfile(output.filepath) && rm(output.filepath; force = output.overwrite)

    Q = simulation.state
    mpicomm = MPI.COMM_WORLD
    iteration = output.iteration

    steps = ClimateMachine.ODESolvers.getsteps(odesolver)
    time = ClimateMachine.ODESolvers.gettime(odesolver)

    file = jldopen(output.filepath, "a+")
    JLD2.Group(file, "state")
    JLD2.Group(file, "time")
    file["state"][string(steps)] = Array(Q)
    file["time"][string(steps)] = time
    close(file)

    jldcallback = ClimateMachine.GenericCallbacks.EveryXSimulationSteps(iteration) do (s = false)
        steps = ClimateMachine.ODESolvers.getsteps(odesolver)
        time = ClimateMachine.ODESolvers.gettime(odesolver)
        @info steps, time
        file = jldopen(output.filepath, "a+")
        file["state"][string(steps)] = Array(Q)
        file["time"][string(steps)] = time
        close(file)
        return nothing
    end
end

struct SolverConfig{SO, MP, QQ, DG, NA}
    solver::SO
    mpicomm::MP
    Q::QQ
    dg::DG
    name::NA
end

function create_callback(output::VTKOutput, simulation, odesolver)

    mpicomm = MPI.COMM_WORLD
    iteration = output.iteration

    solver_config = SolverConfig(simulation.odesolver, mpicomm, simulation.state, simulation.dgmodel, simulation.name)
    vtkcallback = ClimateMachine.Callbacks.vtk(iteration, solver_config, output.outputdir, output.number_sample_points)
end

function create_callback(output::ExponentialFiltering, simulation, odesolver)
    filterorder = 20
    filter = ExponentialFilter(simulation.dgmodel.grid, 0, filterorder)
    cbfilter = GenericCallbacks.EveryXSimulationSteps(1) do
        Filters.apply!(
            simulation.state,
            (:θ,),
            simulation.dgmodel.grid,
            filter,
            state_auxiliary = simulation.dgmodel.state_auxiliary,
        )
        nothing
    end
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
# ---
# title: 664. Strange Printer
# id: problem664
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Depth-first Search
# link: <https://leetcode.com/problems/strange-printer/description/>
# hidden: true
# ---
#
# There is a strange printer with the following two special requirements:
#
#   1. The printer can only print a sequence of the same character each time.
#   2. At each turn, the printer can print new characters starting from and ending at any places, and will cover the original existing characters.
#
# Given a string consists of lower English letters only, your job is to count
# the minimum number of turns the printer needed in order to print it.
#
# **Example 1:**
#
#
#
#     Input: "aaabbb"
#     Output: 2
#     Explanation: Print "aaa" first and then print "bbb".
#
#
# **Example 2:**
#
#
#
#     Input: "aba"
#     Output: 2
#     Explanation: Print "aaa" first and then print "b" from the second place of the string, which will cover the existing character 'a'.
#
#
# **Hint** : Length of the given string will not exceed 100.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using LinearElasticity
using Test

@testset "LinearElasticity.jl" begin
    include("tensors.jl")
    include("arithmetic.jl")
    include("invariants.jl")
end


text:
# test CommonConvWrapper

using Test
using NLsolve
using IncrementalInference
using Manifolds
using Statistics
using TensorCast

import IncrementalInference: getSample

mutable struct FunctorWorks
  a::Array{Float64,2}
end

@testset "FunctorWorks" begin

function (fw::FunctorWorks)(x)
  fw.a[1,1] = -1.0
  nothing
end

A = rand(2,3)
At = deepcopy(A)
At[1,1] = -1.0

fvar = FunctorWorks(A)
fvar(0.0)
@test At == A

end

mutable struct FunctorArray{T}
  # This was a tuple, but array will like work better in the long term
  fnc!::Function
  a::Array{T, 1}
end

##

@testset "FunctorArray" begin

##

function testarray!(a1::Array{Float64, 2}, a2::Array{Float64,2})
  a1[1,1] = -1.0
  @show a1
  nothing
end

function (fi::FunctorArray)(x)
  fi.fnc!(fi.a...)
end

A = rand(2,3)
B = rand(2,3)
t = Array{Array{Float64,2},1}()
push!(t,A)
push!(t,B)
At = deepcopy(A)
At[1,1] = -1.0

fvar = FunctorArray(testarray!, t)
fvar([0.0])
@test At == A

##

end

##

# abstract Nonparametric <: Function
# This is what the internmediate user would be contributing
mutable struct Pose1Pose1Test{T} <: AbstractRelativeRoots
  Dx::T
end

function getSample(cf::CalcFactor{<:Pose1Pose1Test})
  return rand(cf.factor.Dx, 1)
end

#proposed standardized parameter list, does not have to be functor
function (cf::CalcFactor{<:Pose1Pose1Test})(Dx,
                                            p1,
                                            p2 )
  #
  return Dx[1] - (p2[1] - p1[1])
end

##

@testset "Test with CommonConvWrapper for un-permuted root finding..." begin

##

N = 110
p1 = rand(1,N)
p2 = rand(1,N)
t = Array{Array{Float64,2},1}()
push!(t,p1)
push!(t,p2)

odo = Pose1Pose1Test(Normal(100.0,1.0))

fg = initfg()
X0 = addVariable!(fg, :x0, ContinuousEuclid{1})
initManual!(fg, :x0, [zeros(1) for _ in 1:100])
X1 = addVariable!(fg, :x1, ContinuousEuclid{1})
addFactor!(fg, [:x0;:x1], odo, graphinit=false)

pts = approxConv(fg, getFactor(fg, :x0x1f1), :x1)

## Now check the contents of internal CCW

ccw = IIF._getCCW(fg, :x0x1f1)

ptr_ = ccw.params[ccw.varidx]
@cast tp1[i,j] := ptr_[j][i]
@test 90.0 < Statistics.mean(tp1) < 110.0
ptr_ = ccw.params[1]
@cast tp2[i,j] := ptr_[j][i]
@test -10.0 < Statistics.mean(tp2) < 10.0

##

println("and in the reverse direction")

initManual!(fg, :x1, [100*ones(1) for _ in 1:100])

pts = approxConv(fg, getFactor(fg, :x0x1f1), :x0)

ptr_ = ccw.params[1]
@cast tp1[i,j] := ptr_[j][i]
@test -10.0 < Statistics.mean(tp1) < 10.0
ptr_ = ccw.params[2]
@cast tp2[i,j] := ptr_[j][i]
@test 90.0 < Statistics.mean(tp2) < 110.0

##

end

# use the range only example, should give a circle with nothing in the middle

@testset "Generic convolution testing in factor graph context..." begin

##

N=100
p1 = [randn(1) for _ in 1:N]
d1 = manikde!(p1, Euclidean(1))
p2 = [randn(1) for _ in 1:N]
t = Vector{Vector{Vector{Float64}}}()
push!(t,p1)
push!(t,p2)

fg = initfg()

v1=addVariable!(fg, :x1, ContinuousScalar, N=N)
v2=addVariable!(fg, :x2, ContinuousScalar, N=N)
bws = getBW(d1)[:,1]
f1 = addFactor!(fg, [v1], Prior(manikde!(p1, bws, TranslationGroup(1))) )

odo = Pose1Pose1Test(Normal(100.0,1.0))
f2 = addFactor!(fg, [v1;v2], odo)

tree = buildTreeReset!(fg)

pts_ = getBelief(fg,:x1) |> getPoints
@cast pts[i,j] := pts_[j][i]
@test abs(Statistics.mean(pts)-0.0) < 10.0
pts_ = getBelief(fg,:x2) |> getPoints
@cast pts[i,j] := pts_[j][i]
@test abs(Statistics.mean(pts)-0.0) < 10.0

##

tree = solveTree!(fg)

##

pts_ = getBelief(fg,:x1) |> getPoints
@cast pts[i,j] := pts_[j][i]
@test abs(Statistics.mean(pts)-0.0) < 10.0

pts_ = getBelief(fg,:x2) |> getPoints
@cast pts[i,j] := pts_[j][i]
@test abs(Statistics.mean(pts)-100.0) < 10.0

##

end


text:
using NumericalMethodsforEngineers

ab = SymmetricBandedMatrix(2, Float64[0 0 16; 0 4 5; 8 -4 22])
bb = SymmetricBandedMatrix(3, Float64[0 0 0 16; 0 0 4 2; 0 8 -4 22; 13 5 9 1])
cb = SymmetricBandedMatrix(4, Float64[0 0 0 0 1; 0 0 0 2 3; 0 0 4 5 6; 0 7 8 9 10; 11 12 13 14 15])
db = SymmetricBandedMatrix(2, Float64[0 0 16; 0 4 5; 8 -4 22; 1 2 3; 0 7 6 ; 11 0 15])
eb = SymmetricBandedMatrix(3, Float32[0 0 0 16; 0 0 4 5; 0 8 -4 22; 1 2 3 4; 0 0 7 6 ; 0 11 0 15; 0 0 0 20])
ebi = SymmetricBandedMatrix(3, [0 0 0 16; 0 0 4 5; 0 8 -4 22; 1 2 3 4; 0 0 7 6; 0 11 0 15; 0 0 0 20])
ebf = SymmetricBandedMatrix(3, [0 0 0 16+0.1im; 0 0 4 5; 0 8 -4+9im 22; 1 2 3 4; 0 0 7 6 ; 0 11 0 15; 0 0 0 20])
ebif = SymmetricBandedMatrix(3, [0 0 0 16; 0 0 4+im 5; 0 8 -4 22; 1 2 3 4; 0 0 7 6; 0 11 0 15; 0 0 0 20])
ebif32 = SymmetricBandedMatrix(3, Complex{Float32}[0 0 0 16; 0 0 4+im 5; 0 8 -4 22; 1 2 3 4; 0 0 7 6; 0 11 0 15; 0 0 0 20])

a = full(ab)
res = tosymmetricbandedmatrix(2, a)
@test res == ab

b = full(bb)
res = tosymmetricbandedmatrix(3, b,)
@test res == bb

c = full(cb)
res = tosymmetricbandedmatrix(4, c)
@test res == cb

d = full(db)
res = tosymmetricbandedmatrix(2, d)
@test res == db

e = full(eb)
res = tosymmetricbandedmatrix(3, e)
@test res == eb
res2 = tosymmetricbandedmatrix(e)
@test res2 == res

ei = full(ebi)
res = tosymmetricbandedmatrix(3, ei)
@test res == ebi
res2 = tosymmetricbandedmatrix(ei)
@test res2 == res

ef = full(ebf)
res = tosymmetricbandedmatrix(3, ef)
@test res == ebf

eif = full(ebif)
res = tosymmetricbandedmatrix(3, eif)
@test res == ebif

eif32 = full(ebif32)
res = tosymmetricbandedmatrix(3, eif32)
@test res == ebif32

res2 = copy(res)
@test res == res2
@test res !== res2

res2.bmat[4, 3] = 8.0+9im
@test res != res2
@test res !== res2

res3 = res2
@test res3 == res2
@test res3 === res2

f=sprand(6, 6, 0.3)
ff = Matrix(f)
ffs=ff.+ff'
bm = tosymmetricbandedmatrix(ffs)
@test full(bm) == ffs

ms = sprand(5, 5, 0.1)
# Make sure it is symmetric
sms = issymmetric(ms) ? ms : ms .+ ms'
bms = issymmetric(sms) ? convert(SymmetricBandedMatrix, sms) : convert(SymmetricBandedMatrix, round(sms, 8))
@test full(bms) == sms