text:
# Julia wrapper for header: /usr/include/scip/scip_conflict.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeConflicthdlr(scip, name, desc, priority, conflictcopy, conflictfree, conflictinit, conflictexit, conflictinitsol, conflictexitsol, conflictexec, conflicthdlrdata)
    ccall((:SCIPincludeConflicthdlr, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_CONFLICTHDLRDATA}), scip, name, desc, priority, conflictcopy, conflictfree, conflictinit, conflictexit, conflictinitsol, conflictexitsol, conflictexec, conflicthdlrdata)
end

function SCIPincludeConflicthdlrBasic(scip, conflicthdlrptr, name, desc, priority, conflictexec, conflicthdlrdata)
    ccall((:SCIPincludeConflicthdlrBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_CONFLICTHDLR}}, Cstring, Cstring, Cint, Ptr{Cvoid}, Ptr{SCIP_CONFLICTHDLRDATA}), scip, conflicthdlrptr, name, desc, priority, conflictexec, conflicthdlrdata)
end

function SCIPsetConflicthdlrCopy(scip, conflicthdlr, conflictcopy)
    ccall((:SCIPsetConflicthdlrCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONFLICTHDLR}, Ptr{Cvoid}), scip, conflicthdlr, conflictcopy)
end

function SCIPsetConflicthdlrFree(scip, conflicthdlr, conflictfree)
    ccall((:SCIPsetConflicthdlrFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONFLICTHDLR}, Ptr{Cvoid}), scip, conflicthdlr, conflictfree)
end

function SCIPsetConflicthdlrInit(scip, conflicthdlr, conflictinit)
    ccall((:SCIPsetConflicthdlrInit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONFLICTHDLR}, Ptr{Cvoid}), scip, conflicthdlr, conflictinit)
end

function SCIPsetConflicthdlrExit(scip, conflicthdlr, conflictexit)
    ccall((:SCIPsetConflicthdlrExit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONFLICTHDLR}, Ptr{Cvoid}), scip, conflicthdlr, conflictexit)
end

function SCIPsetConflicthdlrInitsol(scip, conflicthdlr, conflictinitsol)
    ccall((:SCIPsetConflicthdlrInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONFLICTHDLR}, Ptr{Cvoid}), scip, conflicthdlr, conflictinitsol)
end

function SCIPsetConflicthdlrExitsol(scip, conflicthdlr, conflictexitsol)
    ccall((:SCIPsetConflicthdlrExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONFLICTHDLR}, Ptr{Cvoid}), scip, conflicthdlr, conflictexitsol)
end

function SCIPfindConflicthdlr(scip, name)
    ccall((:SCIPfindConflicthdlr, libscip), Ptr{SCIP_CONFLICTHDLR}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetConflicthdlrs(scip)
    ccall((:SCIPgetConflicthdlrs, libscip), Ptr{Ptr{SCIP_CONFLICTHDLR}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNConflicthdlrs(scip)
    ccall((:SCIPgetNConflicthdlrs, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPsetConflicthdlrPriority(scip, conflicthdlr, priority)
    ccall((:SCIPsetConflicthdlrPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONFLICTHDLR}, Cint), scip, conflicthdlr, priority)
end

function SCIPisConflictAnalysisApplicable(scip)
    ccall((:SCIPisConflictAnalysisApplicable, libscip), UInt32, (Ptr{SCIP_},), scip)
end

function SCIPinitConflictAnalysis(scip, conftype, iscutoffinvolved)
    ccall((:SCIPinitConflictAnalysis, libscip), SCIP_RETCODE, (Ptr{SCIP_}, SCIP_CONFTYPE, UInt32), scip, conftype, iscutoffinvolved)
end

function SCIPaddConflictLb(scip, var, bdchgidx)
    ccall((:SCIPaddConflictLb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}), scip, var, bdchgidx)
end

function SCIPaddConflictRelaxedLb(scip, var, bdchgidx, relaxedlb)
    ccall((:SCIPaddConflictRelaxedLb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}, Cdouble), scip, var, bdchgidx, relaxedlb)
end

function SCIPaddConflictUb(scip, var, bdchgidx)
    ccall((:SCIPaddConflictUb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}), scip, var, bdchgidx)
end

function SCIPaddConflictRelaxedUb(scip, var, bdchgidx, relaxedub)
    ccall((:SCIPaddConflictRelaxedUb, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, Ptr{SCIP_BDCHGIDX}, Cdouble), scip, var, bdchgidx, relaxedub)
end

function SCIPaddConflictBd(scip, var, boundtype, bdchgidx)
    ccall((:SCIPaddConflictBd, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BOUNDTYPE, Ptr{SCIP_BDCHGIDX}), scip, var, boundtype, bdchgidx)
end

function SCIPaddConflictRelaxedBd(scip, var, boundtype, bdchgidx, relaxedbd)
    ccall((:SCIPaddConflictRelaxedBd, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BOUNDTYPE, Ptr{SCIP_BDCHGIDX}, Cdouble), scip, var, boundtype, bdchgidx, relaxedbd)
end

function SCIPaddConflictBinvar(scip, var)
    ccall((:SCIPaddConflictBinvar, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPisConflictVarUsed(scip, var, boundtype, bdchgidx, used)
    ccall((:SCIPisConflictVarUsed, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_VAR}, SCIP_BOUNDTYPE, Ptr{SCIP_BDCHGIDX}, Ptr{UInt32}), scip, var, boundtype, bdchgidx, used)
end

function SCIPgetConflictVarLb(scip, var)
    ccall((:SCIPgetConflictVarLb, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPgetConflictVarUb(scip, var)
    ccall((:SCIPgetConflictVarUb, libscip), Cdouble, (Ptr{SCIP_}, Ptr{SCIP_VAR}), scip, var)
end

function SCIPanalyzeConflict(scip, validdepth, success)
    ccall((:SCIPanalyzeConflict, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cint, Ptr{UInt32}), scip, validdepth, success)
end

function SCIPanalyzeConflictCons(scip, cons, success)
    ccall((:SCIPanalyzeConflictCons, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_CONS}, Ptr{UInt32}), scip, cons, success)
end


text:
@testset "generalwhittle" begin
    @testset "Gradient" begin
        @testset "univariate" begin
            @test_throws ArgumentError WhittleLikelihoodInference.grad_generalwhittle!(ones(2), ones(10), ones(2,10), ones(11))
            @test_throws ArgumentError WhittleLikelihoodInference.grad_generalwhittle!(ones(2), ones(10), ones(2,11), ones(10))
        end
        @testset "multivariate" begin
            S,I = fill([@SMatrix ones(2,2) for i in 1:10],2)
            ∇S = [@SMatrix ones(2,2) for j in 1:2, i in 1:10]
            @test_throws ArgumentError WhittleLikelihoodInference.grad_generalwhittle!(ones(2), S[1:9], ∇S[1:9], I)
            @test_throws ArgumentError WhittleLikelihoodInference.grad_generalwhittle!(ones(2), S[1:9], ∇S, I[1:9])
        end
    end
    @testset "Hessian" begin
        @testset "univariate" begin
            @test_throws ArgumentError WhittleLikelihoodInference.hess_generalwhittle!(ones(2), ones(10), ones(2,10), ones(3,10), ones(11))
            @test_throws ArgumentError WhittleLikelihoodInference.hess_generalwhittle!(ones(2), ones(10), ones(2,11), ones(3,10), ones(10))
            @test_throws ArgumentError WhittleLikelihoodInference.hess_generalwhittle!(ones(2), ones(10), ones(2,10), ones(3,11), ones(10))
        end
        @testset "multivariate" begin
            S,I = fill([@SMatrix ones(2,2) for i in 1:10],2)
            ∇S = [@SMatrix ones(2,2) for j in 1:2, i in 1:10]
            ∇²S = [@SMatrix ones(2,2) for j in 1:3, i in 1:10]
            @test_throws ArgumentError WhittleLikelihoodInference.hess_generalwhittle!(ones(2), S[1:9], ∇S[1:9], ∇²S[1:9], I)
            @test_throws ArgumentError WhittleLikelihoodInference.hess_generalwhittle!(ones(2), S[1:9], ∇S, ∇²S[1:9], I[1:9])
            @test_throws ArgumentError WhittleLikelihoodInference.hess_generalwhittle!(ones(2), S[1:9], ∇S[1:9], ∇²S, I[1:9])
        end
    end
    @testset "Hessian offdiag" begin
        @testset "univariate" begin
            @test_throws ArgumentError WhittleLikelihoodInference.offdiag_add_hess_generalwhittle!(ones(2), ones(10), ones(2,10), ones(2,10), ones(11))
            @test_throws ArgumentError WhittleLikelihoodInference.offdiag_add_hess_generalwhittle!(ones(2), ones(10), ones(2,11), ones(2,10), ones(10))
            @test_throws ArgumentError WhittleLikelihoodInference.offdiag_add_hess_generalwhittle!(ones(2), ones(10), ones(2,10), ones(2,11), ones(10))
        end
        @testset "multivariate" begin
            S,I = fill([@SMatrix ones(2,2) for i in 1:10],2)
            ∇S1,∇S2 = fill([@SMatrix ones(2,2) for j in 1:2, i in 1:10],2)
            @test_throws ArgumentError WhittleLikelihoodInference.offdiag_add_hess_generalwhittle!(ones(2), S[1:9], ∇S1[1:9], ∇S2[1:9], I)
            @test_throws ArgumentError WhittleLikelihoodInference.offdiag_add_hess_generalwhittle!(ones(2), S[1:9], ∇S1, ∇S2[1:9], I[1:9])
            @test_throws ArgumentError WhittleLikelihoodInference.offdiag_add_hess_generalwhittle!(ones(2), S[1:9], ∇S1[1:9], ∇S2, I[1:9])
        end
    end
end

ext: jl


text:
export HurwitzMatrix,
    HurwitzDeterminants,
    StabilityMatrix

struct HurwitzDeterminants{T}
    M::MatElem{T}
end

Base.convert(::Type{AbstractMatrix},M::HurwitzDeterminants) = Array(M.M)
Base.show(io::IO, M::HurwitzDeterminants) = print(io, M.M)
Base.show(io::IO, ::MIME"text/plain", z::HurwitzDeterminants{T}) where {T} =
    print(io, "HurwitzDeterminants{$T} iterator:\n", z)

# function Base.iterate(M::HurwitzDeterminants)
#     state = 1
#     return Base.iterate(M, state)
# end

function Base.iterate(M::HurwitzDeterminants, state=1)
    if state > size(M.M, 1)
        return nothing
    else
        return Nemo.det(M.M[1:state,1:state]), state + 1
    end
end

# https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration-1
Base.length(M::HurwitzDeterminants) = size(M.M, 1)
Base.eltype(::Type{HurwitzDeterminants{T}}) where {T} = T

# Input: polynomial p in one variable
function HurwitzMatrix(R::AbstractAlgebra.Ring, p)
    n = length(p) - 1
    M = zero_matrix(R, n, n)
    for i in 1:n
        for j in 1:n
            if (n - 2 * i + j) in 0:n
                M[i,j] = coeff(p, n - 2 * i + j)
            end
        end
    end
    return M
end

function evenodd_coeffs(p)
    l = length(p); ## degree +1
    ## For the zero polynomial, 'odd' would be a 0-element array
    if l < 1
        even = [0]
        odd  = [0]
    else
        even = [coeff(p, i) for i in 0:2:l]
        odd  = [coeff(p, i) for i in 1:2:l]
    end
    return even, odd
end

function evenodd_polys(p)
    even, odd = evenodd_coeffs(p)
    return parent(p)(evalxsqrt(even)), parent(p)(evalxsqrt(odd))
end

function evalxsqrt(poly)
    n = size(poly, 1)
    R = parent(poly[1])
    out = fill(R(0), 2 * n - 1)
    out[1] = poly[1]
    for i in 2:n
        out[2 * i - 1] = poly[i]
    end
    return out
end

function HurwitzMatrix1(R::AbstractAlgebra.Ring, p)
    sylvester_matrix(evenodd_polys(p)...)
end

# ## k's are parameters
# R, (k1, k2) = PolynomialRing(Nemo.QQ, ["k1", "k2"]);
# RQ = FractionField(R)

# ## The ring for the polynomial system
# S, (x, y, z, w, t) = PolynomialRing(RQ, ["x", "y", "z", "w", "t"])
# SQ = FractionField(S)

# ## The ring for the characteristic polynomial.
# T, l = PolynomialRing(SQ, "A")

# P = T([x * y * k1, k1 * k2^2 + 2, k1, //(k1, k2 + k1 * k2) * x])
# Q = T([x, y, z, w, t])
# p, q = evenodd_polys(Q)

# sylvester_matrix(evenodd_polys(p)...)
# -x * resultant_subresultant(p, q) == det(HurwitzMatrix(S, Q))

# # add_row(a::MatrixElem, s::RingElement, i::Int, j::Int, cols = 1:ncols(a))

# R, (k1, k2, k3) = PolynomialRing(Nemo.QQ, ["k1", "k2", "k3"]);
# RQ = FractionField(R)

# ## The ring for the polynomial system
# S, (x, y, z) = PolynomialRing(RQ, ["x", "y", "z"])

# ## The ring for the characteristic polynomial.
# T, l = PolynomialRing(S, "A")

# chp = T([k1 * k2 * y + k2 * k3 * z + k1 * k3, k2 * y + k2 * z + k1 + k3, S(1)])

# p, q = evenodd_polys(chp)
# resultant_subresultant(p,q)
# M = sylvester_matrix(p, q)

# det(M) == resultant_subresultant(p, q)

# det(M) == det(HurwitzMatrix(S, chp))

function ISRO!(M, j)
    c = ncols(M)
    p = M[j,j]
    for i in (j + 1):c
        if !(iszero(M[i,j]))
            add_row!(M, M[i,j] // p, j, i, c)
        end
    end
end

function ISRO!(M)
    return j -> ISRO!(M, j)
end

function HurwitzDeterminants(R::AbstractAlgebra.Ring, p)
    return HurwitzDeterminants(HurwitzMatrix(R, p))
end

# mindeg = size(N,1)-rank(matrix(R,N))
function StabilityMatrix(R::AbstractAlgebra.Ring, J, mindeg::Integer=size(J, 1) - rank(J))
    P, x = PolynomialRing(R, "x")
    p = divexact(charpoly(P, J), x^mindeg)
    H = []
    iter = HurwitzDeterminants(R, p)
    print("==============================================\n")
    print("==============================================\n")
    print("$(length(iter)) expressions will be studied\n")
    for (i, D) in enumerate(iter)
        print("\n-- $i th determinant --\n")
        signs = unique(sign.(coeffs(D)))
        print("$(signs)\n")
        if -1 in signs
            H = [H; D]
        end
    end
    print("==============================================\n")
    print("==============================================\n")
    return H
end

function StabilityMatrix(io, R::AbstractAlgebra.Ring, J, mindeg::Integer=size(J, 1) - rank(J))
    P, x = PolynomialRing(R, "x")
    p = divexact(charpoly(P, J), x^mindeg)
    H = []
    iter = HurwitzDeterminants(R, p)
    write(io, "==============================================\n")
    write(io, "==============================================\n")
    write(io, "$(length(iter)) expressions will be studied\n")
    for (i, D) in enumerate(iter)
        write(io, "\n-- $i th determinant --\n")
        # signs = unique(sign.(coeffs(D)))
        signs = unique(sign.(reshape(hcat, coeffs.(coeffs(D)))))
        write(io, "$(signs)\n")
        if -1 in signs
            push!(H, D)
        end
    end
    write(io, "==============================================\n")
    write(io, "==============================================\n")
    return H
end


text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module NumberTheory
using  Nemo, Products

export ModuleNumberTheory
export τ, σ, σ2, ϕ, ω, Ω, ⊥, ⍊, Divisors, PrimeDivisors, Factors, Radical, mods
export V000005, V000010, V000203, V001222, V001221, V008683, V181830, V034444
export I003277, L003277, V061142, V034386, V002110, I050384, L050384, V001157
export Divides, isPrime, isCyclic, isStrongCyclic, isOdd, PrimeList
export isPrimeTo, isStrongPrimeTo, isNonnegative, isPositive, isEven, isSquare
export isComposite, isSquareFree, isPrimePower, isPowerOfPrimes, isPerfectPower

"""
* τ, σ, σ2, ϕ, ω, Ω, ⊥, ⍊
* Divisors, PrimeDivisors, Factors, Radical, mods, Divides, isPrime, isCyclic, isStrongCyclic, isOdd, PrimeList, isPrimeTo, isStrongPrimeTo, isNonnegative, isPositive, isEven, isSquare, isComposite, isSquareFree, isPrimePower, isPowerOfPrimes, isPerfectPower
* V000005, V000010, V000203, V001222, V001221, V008683, V181830, V034444, I003277, L003277, V061142, V034386, V002110, I050384, L050384, V001157
"""
const ModuleNumberTheory = ""

"""
Return true if n is prime false otherwise.
"""
isPrime(n) = Nemo.isprime(fmpz(n))
#function isPrime(n)::Bool Nemo.isprime(fmpz(n)) end

"""
Return factors of ``n``.
"""
Factors(n) = n == 0 ? [] : Nemo.factor(fmpz(n))

"""
Return the positive integers dividing ``n``.
"""
function Divisors(m, dosort=false)
    n = ZZ(m)
    n == ZZ(0) && return fmpz[]
    Nemo.isprime(n) && return [ZZ(1), n]
    d = [ZZ(1)]
    for (p, e) in Nemo.factor(n)
        d *= permutedims([p^i for i in 0:e])
        d = reshape(d, length(d))
    end
    dosort && sort!(d)
    d
end

"""
Return the prime numbers dividing ``n``.
"""
function PrimeDivisors(n)
    n == 0 && return []
    isPrime(n) && return [fmpz(n)]
    f = Factors(n)
    sort!([p for (p, e) in f])
end

"""
Return the radical of ``n`` which is the product of the prime numbers dividing ``n`` (also called the squarefree kernel of ``n``).
"""
Radical(n) = ∏(PrimeDivisors(n))

"""
Return ``Ω(n)``, the number of prime divisors of ``n`` counted with multiplicity (cf. A001222).
"""
function Ω(n)
    n == fmpz(0) && return 0
    isPrime(n) && return fmpz(1)
    f = Factors(n)
    sum([e for (p, e) in f])
end

"""
Return the number of prime divisors of ``n`` counted with multiplicity.
"""
V001222(n) = Ω(n)

"""
Return the result of replacing each prime factor of n with 2.
"""
V061142(n) = 1 << Int(Ω(n))

"""
Return ``ω(n)``,  the number of distinct prime divisors of ``n`` (cf. A001221).
"""
ω(n) = fmpz(length(PrimeDivisors(n)))

"""
Return the number of distinct prime divisors of ``n``.
"""
V001221(n) = ω(n)

"""
Return the number of unitary divisors of ``n``, ``d`` such that ``d`` divides ``n`` and ``d ⊥ n/d``.
"""
V034444(n::Int) = 1 << Int(ω(n))

"""
Return ``τ(n)`` (a.k.a. ``σ_0(n)``), the number of divisors of ``n`` (cf A000005).
"""
τ(n) = Nemo.sigma(fmpz(n), 0)

"""
Return the number of divisors of ``n``.
"""
V000005(n) = τ(n)

"""
Return ``σ(n)`` (a.k.a. ``σ_1(n)``), the sum of the divisors of ``n`` (cf. A000203).
"""
σ(n) = Nemo.sigma(fmpz(n), 1)

"""
Return ``σ2(n)`` (a.k.a. ``σ_2(n)``), the sum of squares of the divisors of ``n`` (cf. A001157).
"""
σ2(n) = Nemo.sigma(fmpz(n), 2)

# σ2(n) = sum(d^2 for d in Divisors(n))

"""
Return ``σ2(n)`` (a.k.a. ``σ_2(n)``), the sum of squares of the divisors of ``n``.
"""
V001157(n) = σ2(n)

"""
Return the Euler totient ``ϕ(n)``, numbers which are ``≤ n`` and prime to ``n``.
"""
ϕ(n) = Nemo.eulerphi(fmpz(n))

"""
Return the number of integers ``≤ n`` and prime to ``n``.
"""
V000010(n) = ϕ(n)

"""
Return the value of the Möbius function ``μ(n)`` which is the sum of the primitive n-th roots of unity.
"""
μ(n) = Nemo.moebiusmu(fmpz(n))

"""
Return the value of the Möbius function ``μ(n)`` which is the sum of the primitive n-th roots of unity.
"""
V008683(n) = μ(n)

"""
Return the sum of the divisors of ``n``.
"""
V000203(n) = σ(n)

"""
Query if ``m`` is prime to ``n``.
"""
isPrimeTo(m, n) = Nemo.gcd(fmpz(m), fmpz(n)) == fmpz(1)

"""
Query if ``m`` is prime to ``n``. Knuth, Graham and Patashnik write in "Concrete Mathematics": "Hear us, O mathematicians of the world! Let us not wait any longer! We can make many formulas clearer by defining a new notation now! Let us agree to write m ⊥ n, and to say "m is prime to n", if m and n are relatively prime."
"""
⊥(m, n) = isPrimeTo(m, n)

"""
Query if ``m`` is strong prime to ``n``. ``m`` is strong prime to ``n`` iff ``m`` is prime to ``n`` and ``m`` does not divide ``n-1``.
"""
isStrongPrimeTo(m, n) = isPrimeTo(m, n) && n ∉ Divisors(m - 1)

"""
Query if ``m`` is strong prime to ``n``. ``m`` is strong prime to ``n`` iff ``m`` is prime to ``n`` and ``m`` does not divide ``n-1``.
"""
⍊(m, n) = isStrongPrimeTo(m, n)

function NumbersStronglyPrimeTo(n::Int)
    P = fmpz[m for m in 1:n if ⊥(m, n)]
    D = Divisors(n - 1)
    return setdiff(P, D)
end

"""
Return the number of integers ``≤ n`` which are strong prime to ``n``.
"""
V181830(n) = n == 0 ? 0 : ϕ(n) - τ(n - 1)

"""
Is ``n`` a cyclic number? ``n`` such that there is just one group of order ``n``.
"""
isCyclic(n) = n == 0 ? false : ⊥(n, ϕ(n))

"""
Iterate over the first ``n`` cyclic numbers.
"""
I003277(n::Int) = Iterators.take(Iterators.filter(isCyclic, Iterators.countfrom(1)), n)

"""
Return the first ``n`` cyclic numbers in an array.
"""
L003277(n::Int) = collect(I003277(n))

"""
Is ``n`` a strong cyclic number?
"""
isStrongCyclic(n) = n == 0 ? false : ⍊(n, ϕ(n))

"""
Iterate over the first ``n`` strong cyclic numbers.
"""
I050384(n::Int) = Iterators.take(Iterators.filter(isStrongCyclic, Iterators.countfrom(1)), n)

"""
Return the first ``n`` strong cyclic numbers in an array.
"""
L050384(n::Int) = collect(I050384(n))

"""
Return the least absolute remainder. mods uses the symmetric representation for integers modulo m, i.e. remainders will be reduced to integers in the range ``[-``div``(|m| - 1, 2), ``div``(|m|, 2)]``.
"""
function mods(b, a)
    b == 0 && return a
    h = a >> 1
    (q, r) = divrem(b, a)
    if h <  r  r -= a end
    if h < -r  r += a end
    r
end

"""
Is the integer ``n`` nonnegative?
"""
isNonnegative(n) = n ≥ 0

"""
Is the integer ``n`` positive?
"""
isPositive(n) = n > 0

"""
Is the integer ``n`` a square number?
"""
isSquare(n) = Nemo.issquare(fmpz(n))

"""
Is the integer ``n`` a composite number?
"""
isComposite(n) = (n > 0) && 1 < Ω(n)

"""
Is the integer ``n`` a squarefree number?
"""
isSquareFree(n) = (n > 0) && ω(n) == Ω(n)

"""
Is the integer ``n`` a prime power?
"""
isPrimePower(n) = ω(n) == 1

"""
Is the integer ``n`` a power of primes?
"""
isPowerOfPrimes(n) = (n == 1) || (ω(n) == 1)

"""
Is the integer ``n`` a perfect powers?
"""
isPerfectPower(n) = ω(n) == 1 && Ω(n) ≠ 1

"""
Return `true` if b is divisible by a, otherwise return `false`.
"""
Divides(a, b) = a ≠ 0 && rem(fmpz(b), fmpz(a)) == fmpz(0)

# Defined in the module Base.
"""
Is ``n`` divisble by 2?
"""
isEven(n) = Base.iseven(n)
# n % 2 == 0 ? true : false

"""
Is ``n`` indivisble by 2?
"""
isOdd(n) = Base.isodd(n)
# n % 2 != 0 ? true : false

"""
Return the primorial of ``n``, the product of the primes ``≤ n``.
"""
V034386(n) = Nemo.primorial(n)

"""
Return a list of the first ``n`` primes.
"""
PrimeList(len::Int) = Iterators.take(Iterators.filter(isPrime, Iterators.countfrom(1)), len)

"""
Return the product of first ``n`` primes.
"""
V002110(n) = ∏(PrimeList(n))

# In the module GaussFactorial are the definitions of
# HasPrimitiveRoot
# HasNoPrimitiveRoot

# In the module Abundant
# isAbundant

# Further indicators, less suited for computations.

# isA008578(n::Int) = all(⊥(k, n)     for k in 1:n-1)
# isA002182(n::Int) = all(τ(k) < τ(n) for k in 1:n-1)
# isA002110(n::Int) = all(ω(k) < ω(n) for k in 1:n-1)
# isA131577(n::Int) = all(Ω(k) < Ω(n) for k in 1:n-1)

#START-TEST-########################################################

using Test, SeqTests, SeqUtils

function test()

    # 0-based version of sequences
    Data = Dict{Int,Array{fmpz}}(034386 => [1, 1, 2, 6, 6, 30, 30, 210, 210, 210],
    061142 => [1, 1, 2, 2, 4, 2, 4, 2, 8, 4],
    002110 => [1, 2, 6, 30, 210, 2310, 30030, 510510, 9699690, 223092870],
    000005 => [0, 1, 2, 2, 3, 2, 4, 2, 4, 3],
    000010 => [0, 1, 1, 2, 2, 4, 2, 6, 4, 6],
    000203 => [0, 1, 3, 4, 7, 6, 12, 8, 15, 13],
    001222 => [0, 0, 1, 1, 2, 1, 2, 1, 3, 2],
    001221 => [0, 0, 1, 1, 1, 1, 2, 1, 1, 1],
    008683 => [0, 1, -1, -1, 0, -1, 1, -1, 0, 0],
    181830 => [0, 1, 0, 0, 0, 1, 0, 2, 2, 2],
    034444 => [1, 1, 2, 2, 2, 2, 4, 2, 2, 2])

    @testset "NumTheory" begin
        @test τ(7560) == 64
        @test τ(46080) == 66
        @test τ(25920) == 70

        @test σ(7560) == 28800
        @test σ(46080) == 159666
        @test σ(25919) == 25920
        @test σ(25920) == 92202

        @test ω(7560) == 4
        @test ω(46080) == 3
        @test ω(25919) == 1
        @test ω(25920) == 3

        @test Ω(7560) == 8
        @test Ω(46080) == 13
        @test Ω(25919) == 1
        @test Ω(25920) == 11

        @test Radical(58564) == 22
        @test Radical(58565) == 58565

        FB(n::Int) = (r = 1; for k in 1:n ⊥(n, k) && (r = mods(r * k, n)) end; r)
        FA(n::Int) = mods(∏([j for j in 1:n if ⊥(j, n)]), n)

        for n in 1:20
            @test FA(n) == FB(n)
        end

        if is_oeis_installed()
            V = [V061142, V000005, V000010, V000203, V001222, V001221, V008683]
            W = [V034386, V002110, V181830, V034444]
            for v in V SeqTest(v, 'V', 1) end

            L = [L003277]
            SeqTest(L, 'L')
        end
    end

    composita = [false, false, false, false, true, false, true, false]
    @testset "Queries" begin
        for n in 0:7
            @test isComposite(n) == composita[n + 1]
        end
    end
end

function demo()
    for n in 390:400
        println(n, " ---")
        println(Factors(n))
        println(Divisors(n))
        println(τ(n), ", ", σ(n))
        println(PrimeDivisors(n))
        println(Radical(n))
    end

    println()
    println([n for n in 1:200 if isCyclic(n)])
    println([n for n in 1:200 if isStrongCyclic(n)])
    println(L050384(24))

end

"""
[Divisors(n) for n in 1:10000]
    0.278807 seconds (1.32 M allocations: 39.099 MiB, 46.04% gc time)
[Radical(n)  for n in 1:10000]
    0.070448 seconds (257.87 k allocations: 16.681 MiB)
"""
function perf()
    @time [Divisors(n) for n in 1:10000]
    @time [Radical(n)  for n in 1:10000]
end

function main()
    test()
    demo()
    perf()
end

main()

 end # module

#=

390 ---
1 * 5 * 13 * 2 * 3
Nemo.fmpz[1, 5, 13, 65, 2, 10, 26, 130, 3, 15, 39, 195, 6, 30, 78, 390]
16, 1008
Nemo.fmpz[2, 3, 5, 13]
390
391 ---
1 * 17 * 23
Nemo.fmpz[1, 17, 23, 391]
4, 432
Nemo.fmpz[17, 23]
391
392 ---
1 * 7^2 * 2^3
Nemo.fmpz[1, 7, 49, 2, 14, 98, 4, 28, 196, 8, 56, 392]
12, 855
Nemo.fmpz[2, 7]
14
393 ---
1 * 131 * 3
Nemo.fmpz[1, 131, 3, 393]
4, 528
Nemo.fmpz[3, 131]
393
394 ---
1 * 197 * 2
Nemo.fmpz[1, 197, 2, 394]
4, 594
Nemo.fmpz[2, 197]
394
395 ---
1 * 5 * 79
Nemo.fmpz[1, 5, 79, 395]
4, 480
Nemo.fmpz[5, 79]
395
396 ---
1 * 2^2 * 11 * 3^2
Nemo.fmpz[1, 2, 4, 11, 22, 44, 3, 6, 12, 33, 66, 132, 9, 18, 36, 99, 198, 396]
18, 1092
Nemo.fmpz[2, 3, 11]
66
397 ---
1 * 397
Nemo.fmpz[1, 397]
2, 398
Nemo.fmpz[397]
397
398 ---
1 * 2 * 199
Nemo.fmpz[1, 2, 199, 398]
4, 600
Nemo.fmpz[2, 199]
398
399 ---
1 * 7 * 3 * 19
Nemo.fmpz[1, 7, 3, 21, 19, 133, 57, 399]
8, 640
Nemo.fmpz[3, 7, 19]
399
400 ---
1 * 5^2 * 2^4
Nemo.fmpz[1, 5, 25, 2, 10, 50, 4, 20, 100, 8, 40, 200, 16, 80, 400]
15, 961
Nemo.fmpz[2, 5]
10

[1, 2, 3, 5, 7, 11, 13, 15, 17, 19, 23, 29, 31, 33, 35, 37, 41, 43,
47, 51, 53, 59, 61, 65, 67, 69, 71, 73, 77, 79, 83, 85, 87, 89, 91,
95, 97, 101, 103, 107, 109, 113, 115, 119, 123, 127, 131, 133, 137,
139, 141, 143, 145, 149, 151, 157, 159, 161, 163, 167, 173, 177, 179,
181, 185, 187, 191, 193, 197, 199]
[1, 15, 33, 35, 51, 65, 69, 77, 85, 87, 91, 95, 115, 119, 123, 133,
141, 143, 145, 159, 161, 177, 185, 187]
[1, 15, 33, 35, 51, 65, 69, 77, 85, 87, 91, 95, 115, 119, 123, 133,
141, 143, 145, 159, 161, 177, 185, 187]

  0.275991 seconds (1.32 M allocations: 39.054 MiB, 51.98% gc time)
  0.070448 seconds (257.87 k allocations: 16.681 MiB)

=#


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

if false
    import Base: Base, @show
else
    macro show(s)
        return :(println(stdout, $(QuoteNode(s)), " = ", $(esc(s))))
    end
end

include("compiler/ssair/slot2ssa.jl")
include("compiler/ssair/inlining.jl")
include("compiler/ssair/verify.jl")
include("compiler/ssair/legacy.jl")
include("compiler/ssair/EscapeAnalysis/EscapeAnalysis.jl")
include("compiler/ssair/passes.jl")

max_stars_repo_name: Physolia/julia

max_issues_repo_name: Eurus-Holmes/julia

max_forks_repo_name: Nexuscompute/julia


text:
using Test
using LinearAlgebra

include("setupTests.jl")

@testset "jInv" begin
	include("Mesh/runtests.jl")
	include("ForwardShare/runtests.jl")
	include("Utils/runtests.jl")
	include("InverseSolve/runtests.jl")
	include("LinearSolvers/runtests.jl")
end


text:
using CSV, DataFrames, Loess, PyCall, PyPlot
import DataStructures
PyDict(matplotlib."rcParams")["axes.formatter.use_locale"] = true
sexes = ["män", "kvinnor"]
sc9717 = CSV.File("data/secirc9717eurshort.csv") |> DataFrame
eshort = DataStructures.OrderedDict(:ami => "akut hjärtinfarkt", :othihd => "annan IHD",
	:othhd => "annan hjärtsjukdom", :str => "slaganfall",
	:othcirc => "annan cirkulation")

function lpred(yrs, col)
	yrs_f = Vector{Float64}(yrs)
	col_f = Vector{Float64}(col)
	lmod = loess(yrs_f, col_f)
	predict(lmod, yrs_f)
end

function sc9717sex(s)
	df_0 = sc9717[sc9717[:sex].==s, :]
	yrs = df_0[:year]
	df = df_0[:, 3:end]
	df_l = DataFrame()
	for es in keys(eshort)
		df_l[es] = lpred(yrs, df[es])
	end
	Dict(:sex => s, :yrs => yrs, :df => df, :df_l => df_l)
end

function sc9717plot(scdict, col)
	df = scdict[:df]
	df_l = scdict[:df_l]
	yrs = scdict[:yrs]
	pl = plot(yrs, log.(df[col]./10^5), "*", label = eshort[col])
	plot(yrs, log.(df_l[col]./10^5), "-", color = pl[1].get_color())
	grid(1)
	xticks(yrs, rotation = 45)
	legend()
	title("Åldersstandardiserade dödstal $(sexes[scdict[:sex]]) Sverige")
	ylabel("log(dödstal)")
end

sc9717plall(scdict) = map(c -> sc9717plot(scdict, c), collect(keys(eshort)))

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
module GMExtensions

    using Statistics
    using DrWatson
    using ProgressMeter
    using ValueHistories
    using BSON

    using Plots
    using StatsPlots

    using Flux
    using Zygote
    using ConditionalDists
    using IPMeasures
    using GenerativeModels

    using GenerativeModels: AbstractVAE, AbstractGM, ACD

    export save_checkpoint, load_checkpoint

    include("saveload.jl")
    include("layers.jl")
    include("plotrecipes.jl")
    include("callbacks.jl")
    # include("train.jl")
    include("reidentify.jl")
    include("conv.jl")
end


text:
include("test_gds.jl")
include("test_laplacian_rectangle.jl")
include("test_laplacian_sphere.jl")
# electrostatics
include("test_capacitance_spheres.jl")
include("test_capacitance_two_pads.jl")
# magnetostatics
include("test_inductance_wires.jl")
include("test_magnetostatics_superconductor.jl")
# electrodynamics
include("test_modes_box.jl")
include("test_modes_superconductors.jl")
include("test_microstrip_resonator.jl")
include("test_cpw_resonator.jl")


text:
module HodgkinHuxleySingle

using Test
using Conductor, ModelingToolkit, OrdinaryDiffEq, Unitful, IfElse
import ModelingToolkit: isparameter
import Conductor: Na, K
using Unitful: mV, mS, cm, µm, µA, ms

@testset "Hodgkin Huxley Single Compartment" begin

# Symbolic model construction

Vₘ = MembranePotential()

nav_kinetics = [
    Gate(AlphaBetaRates,
         αₘ = IfElse.ifelse(Vₘ == -40.0, 1.0, (0.1*(Vₘ + 40.0))/(1.0 - exp(-(Vₘ + 40.0)/10.0))),
         βₘ = 4.0*exp(-(Vₘ + 65.0)/18.0), p = 3)
    Gate(AlphaBetaRates,
         αₕ = 0.07*exp(-(Vₘ+65.0)/20.0),
         βₕ = 1.0/(1.0 + exp(-(Vₘ + 35.0)/10.0)))]

kdr_kinetics = [
    Gate(AlphaBetaRates,
         αₙ = IfElse.ifelse(Vₘ == -55.0, 0.1, (0.01*(Vₘ + 55.0))/(1.0 - exp(-(Vₘ + 55.0)/10.0))),
         βₙ = 0.125 * exp(-(Vₘ + 65.0)/80.0), p = 4)]

alphamss = ((4.0exp(-3.6111111111111107 - (0.05555555555555555Vₘ)) +
            IfElse.ifelse(Vₘ == -40.0, 1.0, (4.0 + 0.1Vₘ)*((1.0 - exp(-4.0 - (0.1Vₘ)))^-1)))^-1) *
            IfElse.ifelse(Vₘ == -40.0, 1.0, (4.0 + 0.1Vₘ)*((1.0 - exp(-4.0 - (0.1Vₘ)))^-1))

@test isequal(nav_kinetics[1].ss, alphamss)

@named NaV = IonChannel(Sodium, nav_kinetics, 120mS/cm^2)
@named Kdr = IonChannel(Potassium, kdr_kinetics, 36mS/cm^2)
@named leak = PassiveChannel(Leak, 0.3mS/cm^2)

@test isparameter(leak.sys.g) && !isparameter(NaV.sys.g)
channels = [NaV, Kdr, leak]
@test [length(equations(x.sys)) for x in channels] == [3,2,0]

gradients = Equilibria([Na => 50.0mV, K => -77.0mV, Leak => -54.4mV])
area = 4*pi*(20µm)^2
pulse(t, current) = IfElse.ifelse(t > 100.0, IfElse.ifelse(t < 200.0, 0.0004, 0.0), 0.0)

@named neuron = Soma(channels, gradients, stimulus = pulse, area = ustrip(Float64, cm^2, area));

@test length.([equations(neuron.sys),
               states(neuron.sys),
               parameters(neuron.sys)]) == [12,13,8]

t = 300.
sim_sys = Simulation(neuron, time = t*ms, system = true)

@test length.([equations(sim_sys),
               states(sim_sys),
               parameters(sim_sys)]) == [5,5,8]

expect_syms = [:Vₘ, :NaV₊m, :NaV₊h, :Kdr₊n, :Isyn]
@test all(x -> hasproperty(sim_sys, x), expect_syms)

### Hand-written model setup ###
# Steady-state functions
minf(V0) = inv((4.0 + 0.1V0)*inv(1.0 - exp(-4.0 - 0.1V0)) + 4.0*exp(-3.6111111111111107 - 0.05555555555555555V0)) *
           (4.0 + 0.1V0) * inv(1.0 - exp(-4.0 - 0.1V0))
ninf(V0) = (0.55 + 0.01V0) * inv(1.0 - exp(-5.5 - 0.1V0)) * inv(0.125exp(-0.8125 - 0.0125V0) +
                                                                (0.55 + 0.01V0) * inv(1.0 - exp(-5.5 - 0.1V0)))
hinf(V0) = 0.07*exp(-3.25 - 0.05V0) * inv(0.07*exp(-3.25 - 0.05V0) + inv(1.0 + exp(-3.5 - 0.1V0)) )
hand_pulse(t) = ifelse(100. < t < 200. , 0.0004 , 0.0)

# Initial conditions
V0 = -65.0
u0 = zeros(Float64, 4)
p = zeros(Float64, 8)

# Serialized MTK function initial conditions
u0[1] = V0       # Vₘ
u0[2] = minf(V0) # NaV₊m
u0[3] = hinf(V0) # NaV₊h
u0[4] = ninf(V0) # Kdr₊n

# Parameters
p[1]  = 0.001      # cₘ
p[2]  = 5.02655e-5 # aₘ
p[3]  = 50.0       # ENa
p[4]  = -77.0      # EK
p[5]  = -54.4      # El
p[6]  = 120.0      # NaV₊gbar
p[7]  = 36.0       # Kdr₊gbar
p[8]  = 0.3        # leak₊g

# Cleaned up version of what default MTK output looks like
function hodgkin_huxley!(du, u, p, t)
    @inbounds let Vₘ         = u[1],
                  NaV₊m      = u[2],
                  NaV₊h      = u[3],
                  Kdr₊n      = u[4],

                  cₘ          = p[1],
                  aₘ          = p[2],
                  ENa         = p[3],
                  EK          = p[4],
                  El          = p[5],
                  NaV₊gbar    = p[6],
                  Kdr₊gbar    = p[7],
                  leak₊g      = p[8],
                  Iapp        = hand_pulse(t)
    du[1] = (Iapp - aₘ*leak₊g*(Vₘ - El) - Kdr₊gbar*aₘ*Kdr₊n^4*(Vₘ - EK) -
             NaV₊gbar*aₘ*NaV₊h*NaV₊m^3*(Vₘ - ENa))*inv(aₘ)*inv(cₘ)
    du[2] = ifelse(Vₘ == -40.0 , 1.0 , (4.0 + 0.1Vₘ) * inv(1.0 - exp(-4.0 - 0.1Vₘ))) *
            (1 - NaV₊m) - 4.0*NaV₊m*exp(-3.6111111111111107 + -0.05555555555555555Vₘ)

    du[3] = 0.07*exp(-3.25 - 0.05Vₘ)*(1 - NaV₊h) - NaV₊h*inv(1.0 + exp(-3.5 - 0.1Vₘ))
    du[4] = ifelse(Vₘ == -55.0 , 0.1 , (0.55 + 0.01Vₘ) * inv(1.0 - exp(-5.5 - 0.1Vₘ)))*(1 - Kdr₊n) -
            0.125 * Kdr₊n * exp(-0.8125 - 0.0125Vₘ)
    end
    return nothing
end

byhand_prob = ODEProblem{true}(hodgkin_huxley!, u0, (0.,300.), p)
mtk_prob = ODAEProblem(sim_sys, [], (0., t), [])

byhand_sol = solve(byhand_prob, Rosenbrock23(), reltol=1e-9, abstol=1e-9, saveat=0.025);
current_mtk_sol = solve(mtk_prob, Rosenbrock23(), reltol=1e-9, abstol=1e-9, saveat=0.025);

tsteps = 0.0:0.025:300.0
byhand_out = Array(byhand_sol(tsteps, idxs=1))
current_mtk_out = current_mtk_sol(tsteps)[Vₘ]

@test isapprox(byhand_out, current_mtk_out, rtol=0.0001)

end #testset

end # module


text:
using Flux
println("NNPDE_tests")
using DiffEqFlux
println("Starting Soon!")
using ModelingToolkit
using DiffEqBase
using Test, NeuralPDE
using GalacticOptim
using Optim
using CUDA
using Quadrature

using Random
Random.seed!(100)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end
CUDA.allowscalar(false)
const gpuones = cu(ones(1))

## ODE
@parameters θ
@variables u(..)
Dθ = Differential(θ)

# 1D ODE
eq = Dθ(u(θ)) ~ θ^3 + 2*θ + (θ^2)*((1+3*(θ^2))/(1+θ+(θ^3))) - u(θ)*(θ + ((1+3*(θ^2))/(1+θ+θ^3)))

# Initial and boundary conditions
bcs = [u(0.) ~ 1.0]

# Space and time domains
domains = [θ ∈ IntervalDomain(0.0,1.0)]
# Discretization
dt = 0.1
# Neural network
inner = 12
chain = FastChain(FastDense(1,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,inner,Flux.σ),
                  FastDense(inner,1),(u,p)->gpuones .* u)

initθ = DiffEqFlux.initial_params(chain) |> gpu

strategy = NeuralPDE.GridTraining(dt)
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             strategy;
                                             init_params = initθ
                                             )

pde_system = PDESystem(eq,bcs,domains,[θ],[u])
prob = NeuralPDE.discretize(pde_system,discretization)
symprob = NeuralPDE.symbolic_discretize(pde_system,discretization)
res = GalacticOptim.solve(prob, ADAM(1e-1); cb = cb, maxiters=1000)
phi = discretization.phi

analytic_sol_func(t) = exp(-(t^2)/2)/(1+t+t^3) + t^2
ts = [domain.domain.lower:dt/10:domain.domain.upper for domain in domains][1]
u_real  = [analytic_sol_func(t) for t in ts]
u_predict  = [first(Array(phi(t,res.minimizer))) for t in ts]

@test u_predict ≈ u_real atol = 0.2

# t_plot = collect(ts)
# plot(t_plot ,u_real)
# plot!(t_plot ,u_predict)

## 1D PDE

@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2

eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
bcs = [u(0,x) ~ cos(x),
        u(t,0) ~ exp(-t),
        u(t,1) ~ exp(-t) * cos(1)]

domains = [t ∈ IntervalDomain(0.0,1.0),
          x ∈ IntervalDomain(0.0,1.0)]

pdesys = PDESystem(eq,bcs,domains,[t,x],[u])

inner = 12
chain = FastChain(FastDense(2,inner,Flux.σ),
                  FastDense(12,12,Flux.σ),
                  FastDense(inner,1),(u,p)->gpuones .* u)

quadrature_strategy = NeuralPDE.QuadratureTraining(quadrature_alg=HCubatureJL(),
                                                   reltol = 1e-2, abstol = 1e-2,
                                                   maxiters = 50)
initθ = initial_params(chain) |>gpu
discretization = NeuralPDE.PhysicsInformedNN(chain,
                                             quadrature_strategy;
                                             init_params = initθ)
prob = NeuralPDE.discretize(pdesys,discretization)
symprob = NeuralPDE.symbolic_discretize(pdesys,discretization)

res = GalacticOptim.solve(prob, ADAM(0.1); cb = cb, maxiters=1000)
phi = discretization.phi

u_exact = (t,x) -> exp.(-t) * cos.(x)
ts,xs = [domain.domain.lower:0.01:domain.domain.upper for domain in domains]
u_predict = reshape([first(Array(phi([t,x],res.minimizer))) for t in ts for x in xs],(length(ts),length(xs)))
u_real = reshape([u_exact(t,x) for t in ts  for x in xs ], (length(ts),length(xs)))
diff_u = abs.(u_predict .- u_real)

@test u_predict ≈ u_real atol = 0.3

# p1 = plot(ts, xs, u_real, linetype=:contourf,title = "analytic");
# p2 = plot(ts, xs, u_predict, linetype=:contourf,title = "predict");
# p3 = plot(ts, xs, diff_u,linetype=:contourf,title = "error");
# plot(p1,p2,p3)


text:
module Plugins

import Base.length, Base.iterate, Base.get, Base.getindex

export PluginStack, Plugin,
    hooks, hooklist, hook_cache,
    call_optional,
    customtype

"""
    abstract type Plugin

Provides default implementations of lifecycle hooks.
"""
abstract type Plugin end

"""
    PluginStack(plugins, hookfns = [])

Manages the plugins loaded into an application.

It provides fast access to the plugins by symbol, e.g. `pluginstack[:logger]`. Collection methods and iteration interface are implemented.

The pluginstack is created from a list of plugins, and optionally a list of hook functions. If hook functions are provided,
the `hooks()`` function can be called to

"""
mutable struct PluginStack
    plugins::Array{Plugin}
    hookfns
    symbolcache::Dict{Symbol, Plugin}
    hookcache::Union{NamedTuple, Nothing}
    PluginStack(plugintypes, hookfns = []; options...) = begin
        plugins = instantiate(plugintypes; options...)
        stack = new(plugins, hookfns, symbolcache(plugins), nothing)
        rebuild_cache!(stack)
        return stack
    end
end

symbolcache(plugins) = Dict([(symbol(plugin), plugin) for plugin in plugins])

Base.push!(stack::PluginStack, items...) = update!(stack, Base.push!, items...)
Base.pushfirst!(stack::PluginStack, items...) = update!(stack, Base.pushfirst!, items...)
Base.pop!(stack::PluginStack, items...) = update!(stack, Base.pop!, items...)
Base.popfirst!(stack::PluginStack, items...) = update!(stack, Base.popfirst!, items...)
Base.empty!(stack::PluginStack, items...) = update!(stack, Base.empty!, items...)

function update!(stack::PluginStack, op, args...)
    retval = op(stack.plugins, args...)
    rebuild_cache!(stack)
    return retval
end

Base.isempty(stack::PluginStack) = Base.isempty(stack.plugins)
Base.length(stack::PluginStack) = length(stack.plugins)
Base.iterate(stack::PluginStack) = iterate(stack.plugins)
Base.iterate(stack::PluginStack, state) = iterate(stack.plugins, state)

Base.get(stack::PluginStack, key::Symbol, default=nothing) = get(stack.symbolcache, key, default)
Base.getindex(stack::PluginStack, idx) = getindex(stack.plugins, idx)
Base.getindex(stack::PluginStack, key::Symbol) = get(stack, key)
Base.setindex!(stack::PluginStack, args...) = update!(stack, setindex!, args...)

include("hooks.jl")
include("deps.jl")
include("assembled_types.jl")
include("stages.jl")
include("lifecycle.jl")

end # module


text:
import Core.Intrinsics.ccall
ccall(:jl_new_main_module, Any, ())

baremodule Base

eval(x) = Core.eval(Base,x)
eval(m,x) = Core.eval(m,x)

include = Core.include

include("exports.jl")

if false
    # simple print definitions for debugging. enable these if something
    # goes wrong during bootstrap before printing code is available.
    show(x::ANY) = ccall(:jl_static_show, Void, (Ptr{Void}, Any),
                         Intrinsics.pointerref(Intrinsics.cglobal(:jl_uv_stdout,Ptr{Void}),1), x)
    print(x::ANY) = show(x)
    println(x::ANY) = ccall(:jl_, Void, (Any,), x)
    print(a::ANY...) = for x=a; print(x); end
end

## Load essential files and libraries

include("base.jl")
include("reflection.jl")
include("build_h.jl")
include("version_git.jl")
include("c.jl")

# core operations & types
include("promotion.jl")
include("tuple.jl")
include("range.jl")
include("expr.jl")
include("error.jl")

# core numeric operations & types
include("bool.jl")
include("number.jl")
include("int.jl")
include("operators.jl")
include("pointer.jl")

include("float.jl")
include("complex.jl")
include("rational.jl")

# core data structures (used by type inference)
include("abstractarray.jl")
include("reduce.jl")

include("subarray.jl")
include("array.jl")
include("bitarray.jl")
include("intset.jl")
include("dict.jl")
include("set.jl")
include("hashing.jl")
include("iterator.jl")

# compiler
import Core.Undef  # used internally by compiler
include("inference.jl")

# For OS specific stuff in I/O
include("osutils.jl")

const DL_LOAD_PATH = ByteString[]
@osx_only push!(DL_LOAD_PATH, "@executable_path/../lib/julia")
@osx_only push!(DL_LOAD_PATH, "@executable_path/../lib")

# strings & printing
include("char.jl")
include("ascii.jl")
include("utf8.jl")
include("utf16.jl")
include("utf32.jl")
include("iobuffer.jl")
include("string.jl")
include("utf8proc.jl")
importall .UTF8proc
include("regex.jl")
include("base64.jl")
importall .Base64

# system & environment
include("libc.jl")
include("env.jl")
include("errno.jl")
using .Errno
include("path.jl")
include("intfuncs.jl")

# I/O
include("task.jl")
include("io.jl")
include("show.jl")
include("stream.jl")
include("socket.jl")
include("stat.jl")
include("fs.jl")
importall .FS
include("process.jl")
include("multimedia.jl")
importall .Multimedia
ccall(:jl_get_uv_hooks, Void, ()) # TODO: should put this in _init
include("grisu.jl")
import .Grisu.print_shortest
include("file.jl")
include("methodshow.jl")

# core math functions
include("floatfuncs.jl")
include("math.jl")
importall .Math
const (√)=sqrt
const (∛)=cbrt
include("float16.jl")

# multidimensional arrays
include("cartesian.jl")
using .Cartesian
include("multidimensional.jl")

include("primes.jl")

# concurrency and parallelism
include("serialize.jl")
include("multi.jl")

# Polling (requires multi.jl)
include("poll.jl")

# code loading
include("loading.jl")

begin
    SOURCE_PATH = ""
    include = function(path)
        prev = SOURCE_PATH
        path = joinpath(dirname(prev),path)
        SOURCE_PATH = path
        Core.include(path)
        SOURCE_PATH = prev
    end
end

# reduction along dims
include("reducedim.jl")  # macros in this file relies on string.jl

# basic data structures
include("ordering.jl")
importall .Order
include("collections.jl")

# Combinatorics
include("sort.jl")
importall .Sort
include("combinatorics.jl")

# rounding utilities
include("rounding.jl")
importall .Rounding

# BigInts and BigFloats
include("gmp.jl")
importall .GMP
include("mpfr.jl")
importall .MPFR
big(n::Integer) = convert(BigInt,n)
big(x::FloatingPoint) = convert(BigFloat,x)
big(q::Rational) = big(num(q))//big(den(q))
big(z::Complex) = complex(big(real(z)),big(imag(z)))
@vectorize_1arg Number big

# more hashing definitions
include("hashing2.jl")

# random number generation
include("librandom.jl")
include("random.jl")
importall .Random

# distributed arrays and memory-mapped arrays
include("darray.jl")
include("mmap.jl")
include("sharedarray.jl")

# utilities - version, timing, help, edit, metaprogramming
include("version.jl")
include("datafmt.jl")
importall .DataFmt
include("deepcopy.jl")
include("util.jl")
include("interactiveutil.jl")
include("replutil.jl")
include("test.jl")
include("meta.jl")
include("i18n.jl")
include("help.jl")
using .I18n
using .Help
push!(I18n.CALLBACKS, Help.clear_cache)

# SIMD loops
include("simdloop.jl")
importall .SimdLoop

# frontend
include("Terminals.jl")
include("LineEdit.jl")
include("REPLCompletions.jl")
include("REPL.jl")
include("client.jl")

# sparse matrices and linear algebra
include("sparse.jl")
importall .SparseMatrix
include("linalg.jl")
importall .LinAlg
const ⋅ = dot
const × = cross
include("broadcast.jl")
importall .Broadcast

# statistics
include("statistics.jl")

# signal processing
include("fftw.jl")
include("dsp.jl")
importall .DSP

# (s)printf macros
include("printf.jl")
importall .Printf

# system information
include("sysinfo.jl")
import .Sys.CPU_CORES

# mathematical constants
include("constants.jl")

# Numerical integration
include("quadgk.jl")
importall .QuadGK

# deprecated functions
include("deprecated.jl")

# package manager
include("pkg.jl")
const Git = Pkg.Git

# base graphics API
include("graphics.jl")

# profiler
include("profile.jl")
importall .Profile

function __init__()
    # Base library init
    reinit_stdio()
    Multimedia.reinit_displays() # since Multimedia.displays uses STDOUT as fallback
    fdwatcher_init()
end

include("precompile.jl")

include = include_from_node1

end # baremodule Base

using Base
importall Base.Operators

Base.isfile("userimg.jl") && Base.include("userimg.jl")

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# Each new term in the Fibonacci sequence is generated by adding the previous two terms.
# By starting with 1 and 2, the first 10 terms will be:

# 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, ...

# By considering the terms in the Fibonacci sequence whose values do not exceed
# four million, find the sum of the even-valued terms.

is_even(x) = ((x % 2) == 0)
sum = 0
fibonacci(n) = n < 2 ? n : fibonacci(n-1) + fibonacci(n-2)
i = 2
current = fibonacci(i)

while current < 4e6

    println(current)
    global current, i, sum
    sum += is_even(current) ? current : 0

    i += 1
    current = fibonacci(i)
end
println()
println(sum)
println()

max_stars_repo_name: pedvide/project-euler

max_issues_repo_name: pedvide/project-euler

max_forks_repo_name: pedvide/project-euler


text:
if Base.VERSION >= v"1.6.0-DEV.843"
    depwarn(msg, funcsym) = Base.depwarn(msg, funcsym; force=true)
else
    depwarn(msg, funcsym) = Base.depwarn(msg, funcsym)
end

function ReadNextElement( S::String , List::Array{ String, 1 } )
    depwarn("`ReadNextElement(s, elements)` is deprecated, please use `tryparseelement(s, idx, elements)` instead.", :ReadNextElement)
    symbol, idxnext = tryparseelement(S, 1, List)
    return symbol, S[idxnext:end]
end

function ReadNextNumeric( S::String)
    depwarn("`ReadNextNumeric(s)` is deprecated, please use `tryparseint16(s, idx)` instead.", :ReadNextNumeric)
    i, idxnext = tryparseint16(S, 1)
    return i, S[idxnext:end]
end

function ReadNextCharge( S::String)
    depwarn("`ReadNextCharge(s)` is deprecated, please use `tryparsecharge(s, idx)` instead.", :ReadNextCharge)
    c, idxnext = tryparsecharge(S, 1)
    return c, S[idxnext:end]
end


text:
using KitBase
using LinearAlgebra, Plots, JLD2
using Optim, DiffEqFlux
using CSV,DataFrames
using PyCall

import KitML
cd(@__DIR__)
include("math.jl")

# model initialization
# Force CPU Learning
py"""
import os
os.environ['CUDA_VISIBLE_DEVICES'] = '-1'
"""

model = KitML.load_model("best_model.h5"; mode = :tf)

begin
    #single cell case
end
# Variable initialization
begin
    # space
    x0 = -1.5
    x1 = 1.5
    y0 = -1.5
    y1 = 1.5
    nx = 100#100
    ny = 100#100
    dx = (x1 - x0) / nx
    dy = (y1 - y0) / ny

    pspace = KitBase.PSpace2D(x0, x1, nx, y0, y1, ny)

    # time
    tEnd = 1.0
    cfl = 0.8
    dt = cfl / 2 * (dx * dy) / (dx + dy)

    # quadrature
    quadratureorder = 20
    points, triangulation = KitBase.octa_quadrature(quadratureorder)
    weights = KitBase.quadrature_weights(points, triangulation)
    nq = size(points, 1)

    # particle
    SigmaS = 0 * ones(ny + 4, nx + 4)
    SigmaA = 0 * ones(ny + 4, nx + 4)
    SigmaT = SigmaS + SigmaA

    # moments
    L = 1
    ne = GetBasisSize(L,3)
    phi = zeros(Float32, ne, nx, ny)
    α = zeros(Float32, ne, nx, ny)
    m = ComputeSphericalBasisKarth(points,1,3)
end

global t = 0.0
flux1 = zeros(ne, nx + 1, ny)
flux2 = zeros(ne, nx, ny + 1)

αT = zeros(nx*ny, ne)
phiT = zeros(nx*ny, ne)

#metaIter = 2
for metaIter in 1:500
# csv logging
df = DataFrame(Iter = Int64[], RelErrCells = Float64[], MaxRelror = Float64[])

# initial distribution
begin
    init_field(x, y) = (1.5+cos(2*π*x)*cos(2*π*y))
    for j = 1:nx
        for i = 1:ny
            y = y0 + (i - 0.5) * dy
            x = x0 + (j - 0.5) * dx
            # only zeroth order moment is non-zero
            phi[1, i, j] = init_field(x, y)
            phi[2, i, j] = 0.9/3.0 * phi[1, i, j]
            phi[3, i, j] = 0.9/3.0 * phi[1, i, j]
            phi[4, i, j] = 0.9/3.0 * phi[1, i, j]
        end
    end
end

println("metaIteration $(metaIter)")

# mechanical solver

for iter in 1:200
    println("iteration $(iter) , meta-iter $(metaIter)")

    #restructuring of the matrix
    @inbounds for i = 1:ne
        phiT[:, i] .= phi[i, :, :][:]
    end
    #Call the neural network
    αT = KitML.neural_closure(model, phiT)
    @inbounds for i = 1:ne
        α[i, :, :] .= reshape(αT[:, i], nx, ny)
    end

    #### Introduce helper variables
    phi_old = copy(phi)
    phi_train = Vector()
    h_train = Vector()
    errorsFound = false
    #####

    # realizability reconstruction
    error = 0.0
    count = 0
    maxError = 0.0
    @inbounds for j = 1:ny
        for i = 1:nx
            phi[:, i, j] .= KitBase.realizable_reconstruct(α[:, i, j], m, weights, KitBase.maxwell_boltzmann_dual_prime)
            # Check the prediction error
            error = norm( (phi[:, i, j]-phi_old[:, i, j] )/phi_old[:, i, j] ,2)^2 # take care of the criterion
            if error > maxError
                maxError = error
            end
            #println(error)
            if error > 5e-4 || phi[:, i, j] == NaN
                #println(error)
                errorsFound = true
                count = count +1
            end
        end
    end

    println("Max error: $(maxError)")
    #retrain network on cells with too much error, if there is any
    if errorsFound
        for j = 1:ny
            for i = 1:nx
                #add point to retraining batch. compute alpha and h w.r.t phi_old
                res = KitBase.optimize_closure(α[:, i, j], m, weights, phi_old[:, i, j], KitBase.maxwell_boltzmann_dual)
                α[:, i, j] .= res.minimizer
                phi[:, i, j] .= KitBase.realizable_reconstruct(α[:, i, j], m, weights, KitBase.maxwell_boltzmann_dual_prime)

                #Setup trianing data
                h = computeEntropyH(α[:, i, j], m, weights,maxwell_boltzmann_primal, KitBase.maxwell_boltzmann_dual_prime)
                h_train =vcat(h_train, h)
                phi_train =vcat(phi_train, phi[:, i, j])
            end
        end

        print("Current amount of  errorous datapoints: ")
        print(size(h_train)[1]/(nx*ny)*100)
        print(" % ")
        print(" in iter ")
        println(iter)

        # Transform into matrix format
        phiTrainMat = zeros(size(h_train)[1],4)
        hTrainMat = zeros(size(h_train)[1],1)
        for i in 1:size(h_train)[1]
            phiTrainMat[i,1] = phi_train[4*(i-1)+1]
            phiTrainMat[i,2] = phi_train[4*(i-1)+2]
            phiTrainMat[i,3] = phi_train[4*(i-1)+3]
            phiTrainMat[i,4] = phi_train[4*(i-1)+4]
            hTrainMat[i,1]= h_train[i]
        end

        print(size(h_train))
        model.fit(phiTrainMat,hTrainMat,epochs =200,batch_size = 32,validation_split=0.0,verbose= 1)
        println("model retrained")
        model.save("best_model.h5")
    else
        println("Error smaller than 5e-4")
    end

    push!(df,[iter , count/(nx*ny)*100, maxError])

    # flux
    fη1 = zeros(nq)
    @inbounds for j = 1:ny
        for i = 1:nx
            im1 = i-1
            if i == 1  #periodic boundaries
                im1 = nx
            end

            KitBase.flux_kfvs!(fη1, KitBase.maxwell_boltzmann_dual.(α[:, im1, j]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:], points[:, 1], dt)

            for k in axes(flux1, 1)
                flux1[k, i, j] = sum(m[k, :] .* weights .* fη1)
            end
        end
    end

    fη2 = zeros(nq)
    @inbounds for i = 1:nx
        for j = 1:ny
            jm1 = j-1
            if j == 1  #periodic boundaries
                jm1 = ny
            end
            KitBase.flux_kfvs!(fη2, KitBase.maxwell_boltzmann_dual.(α[:, i, jm1]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:], points[:, 2], dt)

            for k in axes(flux2, 1)
                flux2[k, i, j] = sum(m[k, :] .* (weights .* fη2))
            end
        end
    end

    # update
    @inbounds for j = 1:ny
        for i = 1:nx
            ip1 = i + 1
            jp1 = j + 1

            #periodic boundaries
            if i==nx
                ip1 = 1
            end
            if j==ny
                jp1 = 1
            end

            for q = 1:1
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, ip1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, jp1]) / dy +
                    (SigmaS[i, j] * phi[q, i, j] - SigmaT[i, j] * phi[q, i, j]) * dt
            end

            for q = 2:ne
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, ip1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, jp1]) / dy +
                    (- SigmaT[i, j] * phi[q, i, j]) * dt
            end
        end
    end

    global t += dt
    #contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], phi[1, :, :], clims=(0.5, 2.5))
end

    CSV.write("Hist_metaIter_$(metaIter).csv",df)
end

#contourf(pspace.x[1:nx, 1], pspace.y[1, 1:ny], phi[1, :, :])
#cd(@__DIR__)
#gif(anim, "periodic_reinforcment_$(metaIter).gif")


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

include("../../src/Sparso.jl")
include("../../src/simple-show.jl")
using Sparso

set_options(SA_ENABLE, SA_VERBOSE, SA_USE_SPMP, SA_CONTEXT, SA_REPLACE_CALLS)#, SA_DISABLE_COLLECTIVE_STRUCTURE_PREDICTION)

include("./ipm-ref.jl")
include("utils.jl")

if length(ARGS) == 0
  A, b, p = load_ipm_input("ipm/mps/osa-14")
else
  A, b, p = load_ipm_input(ARGS[1])
end

m = size(A, 1)
n = size(A, 2)
println("Problem size = [$m $n]")
#println("\tsum of A=", sum(A))
#println("\tsum of b=", sum(b))
#println("\tsum of p=", sum(p))

println("Original: ")
x, ref_total_time, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time,
    iter, relResidual, objval = ipm_ref(A, b, p)
x, ref_total_time, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time,
    iter, relResidual, objval = ipm_ref(A, b, p)
@printf "\nOriginal ref_total_time = %f\n" ref_total_time
@printf "Original spgemm = %f fact = %f blas1 = %f trslv = %f spmv = %f\n" spgemm_time fact_time blas1_time trslv_time spmv_time
@printf "Original iter %2i, resid = %9.2e, objval = %e\n" iter relResidual objval

println("\n\nWithout context-sensitive optimization: ")
Sparso.set_reuse_inspection(false)
@acc x, ref_total_time, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time,
    iter, relResidual, objval = ipm_ref(A, b, p)
Sparso.set_knob_log_level(1)
@acc x, ref_total_time, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time,
    iter, relResidual, objval = ipm_ref(A, b, p)
@printf "\nAccelerated acc_total_time = %f\n" ref_total_time
@printf "Accelerated spgemm = %f fact = %f blas1 = %f trslv = %f spmv = %f\n" spgemm_time fact_time blas1_time trslv_time spmv_time
@printf "Accelerated iter %2i, resid = %9.2e, objval = %e\n" iter relResidual objval
Sparso.set_reuse_inspection(true)

println("\n\nAccelerated: ")
@acc x, ref_total_time, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time,
    iter, relResidual, objval = ipm_ref(A, b, p)
Sparso.set_knob_log_level(1)
@acc x, ref_total_time, spgemm_time, fact_time, blas1_time, trslv_time, spmv_time,
    iter, relResidual, objval = ipm_ref(A, b, p)
@printf "\nAccelerated acc_total_time = %f\n" ref_total_time
@printf "Accelerated spgemm = %f fact = %f blas1 = %f trslv = %f spmv = %f\n" spgemm_time fact_time blas1_time trslv_time spmv_time
@printf "Accelerated iter %2i, resid = %9.2e, objval = %e\n" iter relResidual objval

#println("\n\nWith manual context-sensitive optimization: ")
#ipm_ref_simplified_with_context_opt(A, b, p)
#println("\tsum of x=", sum(x))

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
# ---
# title: 295. Find Median from Data Stream
# id: problem295
# author: zhwang
# date: 2022-02-26
# difficulty: Hard
# categories: Heap, Design
# link: <https://leetcode.com/problems/find-median-from-data-stream/description/>
# hidden: true
# ---
#
# Median is the middle value in an ordered integer list. If the size of the list
# is even, there is no middle value. So the median is the mean of the two middle
# value.
#
# For example,
#
# `[2,3,4]`, the median is `3`
#
# `[2,3]`, the median is `(2 + 3) / 2 = 2.5`
#
# Design a data structure that supports the following two operations:
#
#   * void addNum(int num) - Add a integer number from the data stream to the data structure.
#   * double findMedian() - Return the median of all elements so far.
#
#
#
# **Example:**
#
#
#
#     addNum(1)
#     addNum(2)
#     findMedian() -> 1.5
#     addNum(3)
#     findMedian() -> 2
#
#
#
#
# **Follow up:**
#
#   1. If all integer numbers from the stream are between 0 and 100, how would you optimize it?
#   2. If 99% of all integer numbers from the stream are between 0 and 100, how would you optimize it?
#
#
## @lc code=start
using LeetCode

using DataStructures

Base.@kwdef struct MedianFinder
    maxheap::BinaryMaxHeap{Int} = BinaryMaxHeap{Int}()
    minheap::BinaryMinHeap{Int} = BinaryMinHeap{Int}()
end

## Base.show(io::IO, heap::MedianFinder) = print(io, "$(heap.maxheap.valtree):max -- min:$(reverse(heap.minheap.valtree))")

function add_num!(heap::MedianFinder, num::Int)
    ## maximal element of maxheap <= minimal element of minheap
    ## length of maxheap <= length of minheap
    hi, ha = heap.minheap, heap.maxheap
    if length(ha) < length(hi)
        ## e.g. [2]--[4,5] =6> [2]--[4,5,6] => [2,4]--[5,6]
        push!(hi, num)
        push!(ha, pop!(hi))
    else
        ## e.g. [4]--[5] =2> [2,4]--[5] => [2]--[4,5]
        push!(ha, num)
        push!(hi, pop!(ha))
    end
    nothing
end

function find_median(heap::MedianFinder)::Float64
    hi, ha = heap.minheap, heap.maxheap
    length(hi) == length(ha) ? (first(hi) + first(ha)) / 2 : first(hi)
end
## @lc code=end

## @lc test=start
@testset "295.find-median-from-data-stream.jl" begin
    @testset "295-case 1" begin
        heap = MedianFinder()
        add_num!(heap, 1)
        add_num!(heap, 2)
        @test find_median(heap) == 1.5
        add_num!(heap, 3)
        @test find_median(heap) == 2
    end
    @testset "295-case 2" begin
        heap = MedianFinder()
        add_num!(heap, 5)
        add_num!(heap, 4)
        @test find_median(heap) == 4.5
        add_num!(heap, 2)
        @test find_median(heap) == 4
        add_num!(heap, 6)
        @test find_median(heap) == 4.5
    end
end
## @lc test=end

## @lc info=start
# link: [solution to 295](https://leetcode-cn.com/problems/find-median-from-data-stream/solution/pythonjulia-er-fen-cha-zhao-dui-pai-xu-b-5bw1/)
## @lc info=end


text:
# Struct that defines a constrained timeseries distribution
mutable struct ConstrainedTimeseriesDistribution
    vartype::Symbol # :continuous, :discrete
    lb::Array{Float64,1} # Lowerbound for continuous variables
    ub::Array{Float64,1} # Upperbound for continous variables
    feasible::Array{Array{Int64,1},1} # Feasible set for discrete variables
    timeseries_distribution::TimeseriesDistribution
end

# Construct a ConstrainedTimeseriesDisribution from a IID{Uniform} distribution
function ConstrainedTimeseriesDistribution(d::IID{Uniform{Float64}})
    lb = d.distribution.a*ones(N_pts(d))
    ub = d.distribution.b*ones(N_pts(d))
    ConstrainedTimeseriesDistribution(:continuous, lb, ub, Array{Int64}[], d)
end

# Construct a ConstrainedTimeseriesDisribution from a IID{Normal} distribution
function ConstrainedTimeseriesDistribution(d::IID{Normal{Float64}})
    lb = -Inf*ones(N_pts(d))
    ub = Inf*ones(N_pts(d))
    ConstrainedTimeseriesDistribution(:continuous, lb, ub, Array{Int64}[], d)
end

# Construct a ConstrainedTimeseriesDisribution from a IID{Categorical} distribution
function ConstrainedTimeseriesDistribution(d::IID{Categorical{Float64, Array{Float64, 1}}})
    feasible = [collect(1:length(d.distribution.p)) for i=1:N_pts(d)]
    ConstrainedTimeseriesDistribution(:discrete, Float64[], Float64[], feasible, d)
end

# Construct a ConstrainedTimeseriesDisribution from a IID{Bernoulli} distribution
function ConstrainedTimeseriesDistribution(d::IID{Bernoulli{Float64}})
    feasible = [[0,1] for i=1:N_pts(d)]
    ConstrainedTimeseriesDistribution(:discrete, Float64[], Float64[], feasible, d)
end

# Construct a ConstrainedTimeseriesDisribution from a GaussianProcess distribution
function ConstrainedTimeseriesDistribution(d::GaussianProcess)
    lb = -Inf*ones(length(d.x))
    ub = Inf*ones(length(d.x))
    ConstrainedTimeseriesDistribution(:continuous, lb, ub, Array{Int64}[], d)
end

# Convert function for ease of construction
Base.convert(::Type{ConstrainedTimeseriesDistribution}, x::TimeseriesDistribution) =  ConstrainedTimeseriesDistribution(x)

# Number of points in the timeseries
N_pts(t::ConstrainedTimeseriesDistribution) = N_pts(t.timeseries_distribution)

# Check if the ConstrainedTimeseriesDistribution is over a discrete random variable
isdiscrete(d::ConstrainedTimeseriesDistribution) = d.vartype == :discrete

# Check if the ConstrainedTimeseriesDistribution is over a continuous random variable
iscontinuous(d::ConstrainedTimeseriesDistribution) = d.vartype == :continuous

# Checks the validity of the constraints
function isfeasible(d::ConstrainedTimeseriesDistribution)
    isdiscrete(d) ? all(length.(d.feasible) .> 0) : all(d.ub .>= d.lb)
end

# Sampler for ConstrainedTimeseriesDistribution
function Base.rand(rng::AbstractRNG, d::ConstrainedTimeseriesDistribution)
    isdiscrete(d) ? rand(rng, d.timeseries_distribution, d.feasible) : rand(rng, d.timeseries_distribution, d.lb, d.ub)
end

# Defines a Multivariate, constrained time series distribution
MvTimeseriesDistribution = Dict{Symbol, ConstrainedTimeseriesDistribution}

# Logpdf of a MvTimeseries sample according to the provided distribution
Distributions.logpdf(t::MvTimeseriesDistribution, y::Dict{Symbol, Array}) = sum([logpdf(t[sym].timeseries_distribution, y[sym]) for sym in keys(t)])

# Check the validity of the constraints
isfeasible(t::MvTimeseriesDistribution) = all(isfeasible.(values(t)))

# Get the number of points in the MVTimeseriesDistribution. Checks for consistency across variables
function N_pts(t::MvTimeseriesDistribution)
    Ns = N_pts.(values(t))
    @assert all([Ns[1] == Ns[i] for i=2:length(Ns)])
    Ns[1]
end

# Adjusts the lower and upper constraints according to a greater than comparison
function greaterthan!(d::ConstrainedTimeseriesDistribution, val::Float64, truthvals::AbstractArray)
    @assert iscontinuous(d)
    truthvals
    d.lb[truthvals .== true] .= max.(d.lb[truthvals .== true], val)
    d.ub[truthvals .== false] .= min.(d.ub[truthvals .== false], val)
end

# Adjusts the lower and upper constraints according to a less than comparison
lessthan!(d::ConstrainedTimeseriesDistribution, val::Float64, truthvals::AbstractArray) = greaterthan!(d, val, flex_not.(truthvals))

# Adjusts the lower and upper constraints according to an equality comparison
function continuous_equality!(d::ConstrainedTimeseriesDistribution, val::Float64, truthvals::AbstractArray)
    truthvals[truthvals .!= true] .= :anybool
    lessthan!(d, val, truthvals)
    greaterthan!(d, val, truthvals)
end

# Adjusts the feasible set according to an equality comparison
function discrete_equality!(d::ConstrainedTimeseriesDistribution, val::Int64, truthvals::AbstractArray)
    @assert isdiscrete(d)
    for i = 1:length(truthvals)
        val_index = findfirst(d.feasible[i] .== val)
        if truthvals[i] == true
            d.feasible[i] = !isnothing(val_index) ? [val] : Int64[]
        elseif truthvals[i] == false
            !isnothing(val_index) && deleteat!(d.feasible[i], val_index)
        end
    end
end

# Constrains the MvTimeseriesDistribution according to the expression and the corresponding truth
function constrain_timeseries!(t::MvTimeseriesDistribution, expr::Expr, truthvals::AbstractArray)
    head, sym, val = (expr.head == :call) ? expr.args : [expr.head, expr.args...]
    val, d = Meta.eval(val), t[sym]
    if isdiscrete(d) && head == Symbol(".==")
        discrete_equality!(d, val, truthvals)
    elseif iscontinuous(d) && head == Symbol(".==")
        continuous_equality!(d, val, truthvals)
    elseif iscontinuous(d) &&  head == Symbol(".<=")
        lessthan!(d, val, truthvals)
    elseif iscontinuous(d) && head == Symbol(".>=")
        greaterthan!(d, val, truthvals)
    else
        @error string("Unknown expression ", expr)
    end
end

# Constrains the time series for each constraint in the array of constraints
function constrain_timeseries!(t::MvTimeseriesDistribution, constraints)
    for (expr, truthvals) in constraints
        constrain_timeseries!(t, expr, truthvals)
    end
end

# Samples a time series fromt the provided MvTimeseriesDistribution
function Base.rand(rng::AbstractRNG, d::MvTimeseriesDistribution)
    results = Dict{Symbol, Array}()
    for (sym, const_dist) in d
        results[sym] = rand(rng, const_dist)
    end
    results
end

# Error type for infeasible constraint
struct InfeasibleConstraint <: Exception
    msg::String
end

# Printing for the error
Base.showerror(io::IO, e::InfeasibleConstraint) = print(io, "Infeasible Constraint:  ", e.msg)

# Tries to sample a time series that satisfies the expression from the provided MvTimeseriesDistribution
function Base.rand(rng::AbstractRNG, expr::Expr, d::MvTimeseriesDistribution; validity_trials = 2)
    N = N_pts(d)
    for i=1:validity_trials
        d2 = MvTimeseriesDistribution()
        for (k, v) in d
            d2[k] = ConstrainedTimeseriesDistribution(v.vartype, copy(v.lb), copy(v.ub), deepcopy(v.feasible), v.timeseries_distribution)
        end
        constraints = sample_constraints(expr, N, rng::AbstractRNG)
        constrain_timeseries!(d2, constraints)
        isfeasible(d2) && return rand(rng, d2)
    end
    throw(InfeasibleConstraint(string("Couldn't find feasible constraints for ", expr)))
end


text:
export compute_fem_viscoelasticity_strain_energy_term

@doc raw"""
    compute_fem_viscoelasticity_strain_energy_term(ε0, σ0, ε, A, B, m, n, h)

Given the constitutive relation
```math
\sigma^{n+1} = S \sigma^n + H (\varepsilon^{n+1}-\varepsilon^n),
```
this function computes
```math
\int_A {\sigma:\delta \varepsilon}\mathrm{d} x = \underbrace{\int_A { B \varepsilon^{n+1}:\delta \varepsilon}\mathrm{d} x}  + \underbrace{ \int_A { A \sigma^{n+1}:\delta \varepsilon}\mathrm{d} x - \int_A { B \varepsilon^{n+1}:\delta \varepsilon}\mathrm{d} x }_f
```
and returns $f$
"""
function compute_fem_viscoelasticity_strain_energy_term(ε0::Union{Array{Float64, 2}, PyObject},
    σ0::Union{Array{Float64, 2}, PyObject},
    S::Union{Array{Float64, 2}, PyObject}, H::Union{Array{Float64, 2}, PyObject},
    m::Int64, n::Int64, h::Float64)
    ε0, σ0, S, H = convert_to_tensor([ε0, σ0, S, H], [Float64, Float64, Float64, Float64])
    f1 = -compute_strain_energy_term(ε0 * H, m, n, h)
    f2 = compute_strain_energy_term(σ0 * S, m, n, h)
    return f1 + f2
end


text:
using ConstructionBase
using Test
using LinearAlgebra

struct Empty end
struct AB{A,B}
    a::A
    b::B
end

@testset "constructorof" begin
    @test constructorof(Empty)() === Empty()
    @inferred constructorof(AB{Int, Int})
    @test constructorof(AB{Int, Int})(1, 2) === AB(1,2)
    @test constructorof(AB{Int, Int})(1.0, 2) === AB(1.0,2)
    @test constructorof(typeof((a=1, b=2)))(1.0, 2) === (a=1.0, b=2)
    @test constructorof(NamedTuple{(:a, :b)})(1.0, 2) === (a=1.0, b=2)
    @test constructorof(Tuple)(1.0, 2) === (1.0, 2)
    @test constructorof(Tuple{Nothing, Missing})(1.0, 2) === (1.0, 2)
end

@testset "getproperties" begin
    o = AB(1, 2)
    @test getproperties(o) === (a=1, b=2)
    @inferred getproperties(o)
    @test getproperties(Empty()) === NamedTuple()

    t = (1,2,3.0)
    @inferred getproperties(t)
    @test t === getproperties(t)
    @test () === getproperties(())
end

@testset "setproperties" begin

    @test setproperties(NamedTuple(), NamedTuple()) === NamedTuple()
    @test setproperties((), NamedTuple()) === ()
    @test setproperties(NamedTuple(), ()) === NamedTuple()
    @test setproperties((), ()) === ()
    @test setproperties(1, ()) === 1
    @test setproperties(1, NamedTuple()) === 1

    @test setproperties((1,), ()) === (1,)
    @test setproperties((1,), NamedTuple()) === (1,)
    @test setproperties((a=1,), ()) === (a=1,)
    @test setproperties((a=1,), NamedTuple()) === (a=1,)
    @test setproperties(AB(1,2), ()) === AB(1,2)
    @test setproperties(AB(1,2), NamedTuple()) === AB(1,2)

    @test setproperties(AB(1,2), (a=2, b=3))   === AB(2,3)
    @test setproperties(AB(1,2), (a=2, b=3.0)) === AB(2,3.0)
    @test setproperties(AB(1,2), a=2, b=3.0) === AB(2,3.0)

    res = @test_throws ArgumentError setproperties(AB(1,2), (a=2, this_field_does_not_exist=3.0))
    msg = sprint(showerror, res.value)
    @test occursin("this_field_does_not_exist", msg)
    @test occursin("overload", msg)
    @test occursin("ConstructionBase.setproperties", msg)

    res = @test_throws ArgumentError setproperties(AB(1,2), a=2, this_field_does_not_exist=3.0)
    msg = sprint(showerror, res.value)
    @test occursin("this_field_does_not_exist", msg)
    @test occursin("overload", msg)
    @test occursin("ConstructionBase.setproperties", msg)

    @test setproperties(42, NamedTuple()) === 42
    @test setproperties(42) === 42

    @test setproperties(Empty(), NamedTuple()) === Empty()
    @test setproperties(Empty()) === Empty()

    @test setproperties((a=1, b=2), (a=1.0,)) === (a=1.0, b=2)
    @test setproperties((a=1, b=2), a=1.0) === (a=1.0, b=2)

    @inferred setproperties(AB(1,2), a=2, b=3.0)
    @inferred setproperties(Empty(), NamedTuple())
    @inferred setproperties((a=1, b=2), a=1.0)
    @inferred setproperties((a=1, b=2), (a=1.0,))

    @test setproperties((1,), ()) === (1,)
    @test setproperties((1,), (10,)) === (10,)
    @test_throws ArgumentError setproperties((1,), (10,20)) === (10,)

    @inferred setproperties((1,2,3), (1,2,3))
    @test setproperties((1,2,3), ()) === (1,2,3)
    @test setproperties((1,2,3), (10.0,)) === (10.0,2,3)
    @test setproperties((1,2,3), (10.0,20)) === (10.0,20,3)
    @test setproperties((1,2,3), (10.0,20,30)) === (10.0,20,30)
    @test_throws ArgumentError setproperties((1,2,3), (10.0,20,30,40))

    @test_throws ArgumentError setproperties((a=1,b=2), (10,20))
    @test_throws ArgumentError setproperties((), (10,))
    @test_throws ArgumentError setproperties((1,2), (a=10,b=20))
end

struct CustomSetproperties
    _a::Int
end
function ConstructionBase.setproperties(o::CustomSetproperties, patch::NamedTuple)
    if isempty(patch)
        o
    elseif propertynames(patch) == (:a,)
        CustomSetproperties(patch.a)
    else
        error()
    end
end

@testset "custom setproperties unambiguous on empty" begin
    o = CustomSetproperties(1)
    @test o === setproperties(o)
    @test o === setproperties(o, NamedTuple())
    @test CustomSetproperties(2) === setproperties(o, a=2)
    @test CustomSetproperties(2) === setproperties(o, (a=2,))
end

@testset "constructors for non-standadard Base and LinearAlgebra etc objects" begin
    A1 = zeros(5, 6)
    A2 = ones(Float32, 5, 6)

    @testset "SubArray" begin
        subarray = view(A1, 1:2, 3:4)
        @test constructorof(typeof(subarray))(getproperties(subarray)...) === subarray
        @test all(constructorof(typeof(subarray))(A2, (Base.OneTo(2), 3:4), 0, 0) .== Float32[1 1; 1 1])
        @inferred constructorof(typeof(subarray))(getproperties(subarray)...)
        @inferred constructorof(typeof(subarray))(A2, (Base.OneTo(2), 3:4), 0, 0)
    end

    @testset "ReinterpretArray" begin
        ra1 = reinterpret(Float16, A1)
        @test constructorof(typeof(ra1))(A1) === ra1
        @test constructorof(typeof(ra1))(getproperties(ra1)...) === ra1
        ra2 = constructorof(typeof(ra1))(A2)
        @test size(ra2) == (10, 6)
        @test eltype(ra2) == Float16
        @inferred constructorof(typeof(ra1))(getproperties(ra1)...)
        @inferred constructorof(typeof(ra1))(A2)
    end

    @testset "PermutedDimsArray" begin
        pda1 = PermutedDimsArray(A1, (2, 1))
        @test constructorof(typeof(pda1))(A1) === pda1
        @test constructorof(typeof(pda1))(getproperties(pda1)...) === pda1
        @test eltype(constructorof(typeof(pda1))(A2)) == Float32
        @inferred constructorof(typeof(pda1))(getproperties(pda1)...)
        @inferred constructorof(typeof(pda1))(A2)
    end

    @testset "Tridiagonal" begin
        d = randn(12)
        dl = randn(11)
        du = randn(11)
        tda = Tridiagonal(dl, d, du)
        @test isdefined(tda, :du2) == false
        @test constructorof(typeof(tda))(dl, d, du) === tda
        @test constructorof(typeof(tda))(getproperties(tda)...) === tda
        # lu factorization defines du2
        tda_lu = lu!(tda).factors
        @test isdefined(tda_lu, :du2) == true
        @test constructorof(typeof(tda_lu))(getproperties(tda_lu)...) === tda_lu
        @test constructorof(typeof(tda_lu))(getproperties(tda)...) !== tda_lu
        @test constructorof(typeof(tda_lu))(getproperties(tda)...) === tda
        @inferred constructorof(typeof(tda))(getproperties(tda)...)
        @inferred constructorof(typeof(tda))(getproperties(tda_lu)...)
    end

    @testset "LinRange" begin
        lr1 = LinRange(1, 7, 10)
        lr2 = LinRange(1.0f0, 7.0f0, 10)
        @test constructorof(typeof(lr1))(1, 7, 10, nothing) === lr1
        @test constructorof(typeof(lr1))(getproperties(lr2)...) === lr2
        @inferred constructorof(typeof(lr1))(getproperties(lr1)...)
        @inferred constructorof(typeof(lr1))(getproperties(lr2)...)
    end

end

@testset "Anonymous function constructors" begin
    function multiplyer(a, b)
        x -> x * a * b
    end

    mult11 = multiplyer(1, 1)
    @test mult11(1) === 1
    mult23 = @inferred constructorof(typeof(mult11))(2.0, 3.0)
    @inferred mult23(1)
    @test mult23(1) === 6.0
    multbc = @inferred constructorof(typeof(mult23))("b", "c")
    @inferred multbc("a")
    @test multbc("a") == "abc"
end

struct Adder{V} <: Function
    value::V
end
(o::Adder)(x) = o.value + x

struct Adder2{V} <: Function
    value::V
    int::Int
end
(o::Adder2)(x) = o.value + o.int + x

struct AddTuple{T} <: Function
    tuple::Tuple{T,T,T}
end
(o::AddTuple)(x) = sum(o.tuple) + x

# A function with an inner constructor with checks
struct Rotation{M} <: Function
    matrix::M
    function Rotation(m)
        @assert isapprox(det(m), 1)
        @assert isapprox(m*m', I)
        new{typeof(m)}(m)
    end
end

@testset "Custom function object constructors still work" begin
    add1 = Adder(1)
    @test add1(1) === 2
    add2 = @inferred ConstructionBase.constructorof(typeof(add1))(2.0)
    @inferred add2(1)
    @test add2(1) == 3.0
    add12 = Adder2(1, 2)
    @test @inferred add12(3) ==  6
    add22 = @inferred ConstructionBase.constructorof(typeof(add12))(2.0, 2)
    @test @inferred add22(3) ==  7.0

    addtuple123 = AddTuple((1, 2, 3))
    @test addtuple123(1) === 7
    addtuple234 = @inferred ConstructionBase.constructorof(typeof(addtuple123))((2.0, 3.0, 4.0))
    @inferred addtuple234(1)
    @test addtuple234(1) === 10.0

    @testset "inner constructor without type parameters is still called" begin
        @test_throws AssertionError constructorof(Rotation{Matrix{Float64}})(zeros(3,3))
    end
end

function funny_numbers(n)::Tuple
    types = [
        Int128, Int16, Int32, Int64, Int8,
        UInt128, UInt16, UInt32, UInt64, UInt8,
        Float16, Float32, Float64,
    ]
    Tuple([T(true) for T in rand(types, n)])
end

function funny_numbers(::Type{NamedTuple}, n)::NamedTuple
    t = funny_numbers(n)
    pairs = map(1:n) do i
        Symbol("a$i") => t[i]
    end
    (;pairs...)
end

for n in [0,1,20,40]
    Sn = Symbol("S$n")
    types = [Symbol("A$i") for i in 1:n]
    fields = [Symbol("a$i") for i in 1:n]
    typed_fields = [:($ai::$Ai) for (ai,Ai) in zip(fields, types)]
    @eval struct $(Sn){$(types...)}
        $(typed_fields...)
    end
    @eval funny_numbers(::Type{$Sn}) = ($Sn)(funny_numbers($n)...)
end

@testset "inference" begin
    @testset "Tuple n=$n" for n in [0,1,2,3,4,5,10,20,30,40]
        t = funny_numbers(n)
        @test length(t) == n
        @test getproperties(t) === t
        @inferred getproperties(t)
        for k in 0:n
            t2 = funny_numbers(k)
            @inferred setproperties(t, t2)
            @test setproperties(t, t2)[1:k] === t2
            @test setproperties(t, t2) isa Tuple
            @test length(setproperties(t, t2)) == n
            @test setproperties(t, t2)[k+1:n] === t[k+1:n]
        end
    end
    @inferred getproperties(funny_numbers(100))
    @inferred setproperties(funny_numbers(100), funny_numbers(90))
    @testset "NamedTuple n=$n" for n in [0,1,2,3,4,5,10,20,30,40]
        nt = funny_numbers(NamedTuple, n)
        @test nt isa NamedTuple
        @test length(nt) == n
        @test getproperties(nt) === nt
        @inferred getproperties(nt)
        for k in 0:n
            nt2 = funny_numbers(NamedTuple, k)
            @inferred setproperties(nt, nt2)
            @test Tuple(setproperties(nt, nt2))[1:k] === Tuple(nt2)
            @test setproperties(nt, nt2) isa NamedTuple
            @test length(setproperties(nt, nt2)) == n
            @test Tuple(setproperties(nt, nt2))[k+1:n] === Tuple(nt)[k+1:n]
        end
    end
    @inferred getproperties(funny_numbers(NamedTuple, 100))
    @inferred setproperties(funny_numbers(NamedTuple, 100), funny_numbers(NamedTuple, 90))

    @inferred setproperties(funny_numbers(S0), funny_numbers(NamedTuple, 0))
    @inferred setproperties(funny_numbers(S1), funny_numbers(NamedTuple, 1))
    @inferred setproperties(funny_numbers(S20), funny_numbers(NamedTuple, 18))
    @inferred setproperties(funny_numbers(S40), funny_numbers(NamedTuple, 38))
    @inferred getproperties(funny_numbers(S0))
    @inferred getproperties(funny_numbers(S1))
    @inferred getproperties(funny_numbers(S20))
    @inferred getproperties(funny_numbers(S40))
end


text:
"""
Struct for DomainwallFermion
"""
struct DomainwallFermion_3D_wing{NC,WilsonFermion} <: AbstractFermionfields_3D{NC}
    w::Array{WilsonFermion,1}
    NC::Int64
    NX::Int64
    #NY::Int64
    #NZ::Int64
    NT::Int64
    L5::Int64
    Dirac_operator::String
    NWilson::Int64

    function DomainwallFermion_3D_wing(L5,NC::T,NX::T,NT::T) where T<: Integer
        x = WilsonFermion_2D_wing(NC,NX,NT)
        xtype = typeof(x)
        w = Array{xtype,1}(undef,L5)
        w[1] = x
        for i=2:L5
            w[i] = similar(x)
        end
        #println(w[2][1,1,1,1,1,1])
        NWilson = length(x)
        Dirac_operator = "Domainwall"
        return new{NC,xtype}(w,NC,NX,NT,L5,Dirac_operator,NWilson)
    end

end

function Base.similar(x::DomainwallFermion_3D_wing{NC,WilsonFermion} ) where {NC,WilsonFermion}
    return DomainwallFermion_3D_wing(x.L5,NC,x.NX,x.NT)
end

function D5DWx!(xout::DomainwallFermion_3D_wing{NC,WilsonFermion} ,U::Array{G,1},
    x::DomainwallFermion_3D_wing{NC,WilsonFermion} ,m,A,L5) where  {NC,WilsonFermion,G <: AbstractGaugefields}

    #temp = temps[4]
    #temp1 = temps[1]
    #temp2 = temps[2]
    clear_fermion!(xout)
    ratio = 1
    #ratio  =-1
    #ratio = xout.L5/L5

    for i5=1:L5
        j5=i5
        D4x!(xout.w[i5],U,x.w[j5],A,2) #Dw*x
        #Dx!(xout.w[i5],U,x.w[j5],A) #Dw*x
        #Wx!(xout.w[i5],U,x.w[j5],temps) #Dw*x
        #4+m
        #massfac4d = 1/(2*A.κ) #4+M
        #massfactor = (massfac4d - 2)+1
        massfactor = -(1/(2*A.κ) + 1)

        set_wing_fermion!(xout.w[i5])
        add!(ratio,xout.w[i5],ratio*massfactor,x.w[j5])
        #add!(ratio,xout.w[i5],ratio,x.w[j5]) #D = x + Dw*x
        set_wing_fermion!(xout.w[i5])

        j5=i5+1
        if 1 <= j5 <= L5
            #-P_-
            if L5 != 2
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        j5=i5-1
        if 1 <= j5 <= L5
            #-P_+
            if L5 != 2
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        if L5 != 1
            if i5==1
                j5 = L5
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end

            if i5== L5
                j5 = 1
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

    end
    set_wing_fermion!(xout)

    #=
    if L5 != xout.L5
        for i5=L5+1:xout.L5
            axpy!(1,x.w[i5],xout.w[i5])
        end
    end
    =#

    return
end

function D5DWdagx!(xout::DomainwallFermion_3D_wing{NC,WilsonFermion} ,U::Array{G,1},
    x::DomainwallFermion_3D_wing{NC,WilsonFermion} ,m,A,L5) where  {NC,WilsonFermion,G <: AbstractGaugefields}

    #temp = temps[4]
    #temp1 = temps[1]
    #temp2 = temps[2]
    clear_fermion!(xout)

    ratio = 1
    #ratio  =-1
    #ratio = xout.L5/L5

    for i5=1:L5
        j5=i5

        #Ddagx!(xout.w[i5],U,x.w[j5],A) #Ddagw*x
        D4dagx!(xout.w[i5],U,x.w[j5],A,2) #Ddagw*x
        #Wdagx!(xout.w[i5],U,x.w[j5],temps) #Ddagw*x
        #4+M
        massfactor = -(1/(2*A.κ) + 1)

        set_wing_fermion!(xout.w[i5])
        add!(ratio,xout.w[i5],ratio*massfactor,x.w[j5]) #D = x + Dw*x

        #add!(ratio,xout.w[i5],ratio,x.w[j5]) #D = x + Ddagw*x
        set_wing_fermion!(xout.w[i5])

        j5=i5+1
        if 1 <= j5 <= L5
            #-P_-
            if L5 != 2
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        j5=i5-1
        if 1 <= j5 <= L5
            #-P_+
            if L5 != 2
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        if L5 != 1
            if i5==1
                j5 = L5
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end

            if i5==L5
                j5 = 1
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

    end

    #=
    if L5 != xout.L5
        for i5=L5+1:xout.L5
            axpy!(1,x.w[i5],xout.w[i5])
        end
    end
    =#

    set_wing_fermion!(xout)

    return
end


text:
# struct to model a hierarchy

@tm1def mutable struct Hierarchy
  Name::Union{String,Nothing}
  Edges::Union{Vector,Nothing}
  Elements::Union{Vector{Element},Nothing}
  ElementAttributes::Union{Vector,Nothing}
  Subset::Union{Vector,Nothing}
  DefaultMember::Union{Vector,Nothing}
  Attributes::Union{Dict,Nothing}
  # There are other fields returned but not sure they're interesting
end

Hierarchy(name::AbstractString) = Hierarchy(Dict("Name" => name))

namefield(hierarchy::Hierarchy) = hierarchy.Name

# functions for endpoints

@api_default function get_hierarchy(
  api::TM1API,
  dimension_name::AbstractString,
  hierarchy_name::AbstractString;
  options...,
)
  # need to check whether these are the right params
  params = Dict("\$expand" => "Edges,Elements,ElementAttributes,Subsets,DefaultMember")
  result = tm1_get_json(
    api,
    "Dimensions('" * dimension_name * "')/Hierarchies('" * hierarchy_name * "')";
    params = params,
    options...,
  )
  Hierarchy(result)
end

@api_default function delete_hierarchy(
  api::TM1API,
  dimension_name::AbstractString,
  hierarchy_name::AbstractString;
  options...,
)

  tm1_delete(
    api,
    "Dimensions('" * dimension_name * "')/Hierarchies('" * hierarchy_name * "')";
    options...,
  )
end


text:
using ShaderAbstractions, LinearAlgebra
using ShaderAbstractions: VertexArray
using Test, GeometryTypes
import GeometryBasics

struct WebGL <: ShaderAbstractions.AbstractContext end

m = GLNormalMesh(Sphere(Point3f0(0), 1f0))

mvao = VertexArray(m)
instances = VertexArray(positions = rand(GeometryBasics.Point{3, Float32}, 100))

x = ShaderAbstractions.InstancedProgram(
    WebGL(),
    "void main(){}\n", "void main(){}\n",
    mvao,
    instances,
    model = Mat4f0(I),
    view = Mat4f0(I),
    projection = Mat4f0(I),
)

@test x.program.fragment_source == read(joinpath(@__DIR__, "test.frag"), String)
@test x.program.vertex_source == read(joinpath(@__DIR__, "test.vert"), String)


text:
import ProximalOperators
import RecursiveArrayTools

@inline function ProximalOperators.prox(
    h::ProximalOperators.ProximableFunction,
    x::RecursiveArrayTools.ArrayPartition,
    gamma...,
)
    # unwrap
    y, fy = ProximalOperators.prox(h, x.x, gamma...)
    # wrap
    return RecursiveArrayTools.ArrayPartition(y), fy
end

@inline function ProximalOperators.gradient(
    h::ProximalOperators.ProximableFunction,
    x::RecursiveArrayTools.ArrayPartition,
)
    # unwrap
    grad, fx = ProximalOperators.gradient(h, x.x)
    # wrap
    return RecursiveArrayTools.ArrayPartition(grad), fx
end

@inline ProximalOperators.prox!(
    y::RecursiveArrayTools.ArrayPartition,
    h::ProximalOperators.ProximableFunction,
    x::RecursiveArrayTools.ArrayPartition,
    gamma...,
) = ProximalOperators.prox!(y.x, h, x.x, gamma...)

@inline ProximalOperators.gradient!(
    y::RecursiveArrayTools.ArrayPartition,
    h::ProximalOperators.ProximableFunction,
    x::RecursiveArrayTools.ArrayPartition,
) = ProximalOperators.gradient!(y.x, h, x.x)


text:
function point_sampling_triangle(v1, v2, v3)
    r1, r2 = rand(2)
    sqrt_r1 = sqrt(r1)
    point = (1 - sqrt_r1) * v1
    point += sqrt_r1 * (1 - r2) * v2
    point += sqrt_r1 * r2 * v3
    return point
end

function point_sampling(
    verts::AbstractVector{<:AbstractPoint{3,T1}},
    triangles::AbstractVector{<:TriangleFace},
    bins::AbstractVector{<:Real}) where T1 <: AbstractFloat
    # @assert issorted(bins) "Last argument must be sorted in increasing order"
    # @assert isapprox(bins[end], 1.0) "The last element must be 1.0"
    random_number = rand()
    idx = searchsortedfirst(bins, random_number)
    face = triangles[idx]
    v1, v2, v3 = verts[face[1]], verts[face[2]], verts[face[3]]
    point = point_sampling_triangle(v1, v2, v3)
    return point
end

function points_sampling(
    verts::AbstractVector{<:AbstractPoint{3,T1}},
    triangles::AbstractVector{<:TriangleFace},
    n_sample::Int) where T1 <: AbstractFloat
    @assert n_sample > 0 "n must be a positive integer"
    areas_array = area.(Ref(verts), triangles)
    cumsum_area = similar(areas_array)
    cumsum_area = cumsum!(cumsum_area, areas_array)
    total_area = cumsum_area[end]
    bins = cumsum_area / total_area
    random_numbers = rand(n_sample)
    idx = searchsortedfirst.(Ref(bins), random_numbers)
    selected_triangles = triangles[idx]
    points = Vector{Point{3,T1}}(undef, n_sample)
    for i in 1:n_sample
        face = selected_triangles[i]
        v1, v2, v3 = verts[face[1]], verts[face[2]], verts[face[3]]
        points[i] = point_sampling_triangle(v1, v2, v3)
    end
    return points
end

function points_sampling(mesh::Mesh, n_sample)
    verts = coordinates(mesh)
    triangles = faces(mesh)
    return points_sampling(verts, triangles, n_sample)
end