text:
using Test
using HMMBase
using Distributions
using LinearAlgebra
using Random

import HMMBase: issquare, warn_logl

Random.seed!(2019)

@testset "Base" begin
    hmm1 = HMM([0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(10, 1)])
    hmm2 = HMM([0.9 0.1; 0.1 0.9], [MvNormal([0.0,0.0], [1.0,1.0]), MvNormal([10.0,10.0], [1.0,1.0])])

    @test size(hmm1) == (2, 1)
    @test size(hmm2) == (2, 2)

    @test nparams(hmm1) == 3
    @test nparams(hmm2) == 3
end

@testset "Base (2)" begin
    hmm1 = HMM([0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(10, 1)])
    hmm2 = HMM([0.9 0.1; 0.1 0.9], [MvNormal([0.0,0.0], [1.0,1.0]), MvNormal([10.0,10.0], [1.0,1.0])])

    @test issquare(hmm1.A)
    @test istransmat(hmm1.A)

    @test hmm1 != hmm2
    @test hmm1 == copy(hmm1)
    @test hmm1 !== copy(hmm1) # !== (identity) not != (equality)
end

@testset "Base (3)" begin
    hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(10, 1)])

    perm = [1, 2]
    hmmp = permute(hmm, perm)
    @test hmmp == hmm

    perm = [2, 1]
    hmmp = permute(hmm, perm)
    @test hmmp != hmm
    @test hmmp.B == hmm.B[perm]
    @test diag(hmmp.A) == diag(hmm.A)[perm]
end

@testset "Base (4)" begin
    hmm1 = HMM([0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(10, 1)])
    hmm2 = HMM([0.9 0.1; 0.1 0.9], [MvNormal([0.0,0.0], [1.0,1.0]), MvNormal([10.0,10.0], [1.0,1.0])])
    hmm3 = HMM([0.9 0.1; 0.1 0.9], [MvNormal(ones(3)), MvNormal(ones(3))])

    # Univariate HMMs should return observations vectors
    # (consistent with Distributions.jl)
    z1, y1 = rand(hmm1, 1000, seq = true)
    y11 = rand(hmm1, z1)

    @test size(z1) == (1000,)
    @test size(y1) == (1000,)
    @test size(y11) == size(y1)

    # Multivariate HMMs should return a `TxK` matrix
    # (different from Distributions.jl which returns `KxT`)
    z2, y2 = rand(hmm2, 1000, seq = true)
    y22 = rand(hmm2, z2)

    @test size(z2) == (1000,)
    @test size(y2) == (1000,2)
    @test size(y22) == size(y2)

    # Rand called with T < 1 should return empty arrays
    z3, y3 = rand(hmm2, 0, seq = true)
    y33 = rand(hmm2, z3)
    @test size(z3) == (0,)
    @test size(y3) == (0,2)
    @test size(y33) == size(y3)

    # Multivariate HMM should work with more observations than states:
    # related to the issue: https://github.com/maxmouchet/HMMBase.jl/issues/12
    y = rand(hmm3, 1000)
    @test size(y) == (1000, 3)
end

@testset "Constructors" begin
    # Test that errors are raised
    # Wrong trans. matrix
    @test_throws ArgumentError HMM(ones(2, 2), [Normal();Normal()])
    # Wrong trans. matrix dimensions
    @test_throws ArgumentError HMM([0.8 0.1 0.1; 0.1 0.1 0.8], [Normal(0, 1), Normal(10, 1)])
    # Wrong number of distributions
    @test_throws ArgumentError HMM([0.8 0.2; 0.1 0.9], [Normal(0, 1), Normal(10, 1), Normal()])
    # Wrong distributions size
    @test_throws ArgumentError HMM([0.8 0.2; 0.1 0.9], [MvNormal(randn(3)), MvNormal(randn(10))])
    # Wrong initial state
    @test_throws ArgumentError HMM([0.1;0.1], [0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(10, 1)])
    # Wrong initial state length
    @test_throws ArgumentError HMM([0.1;0.1;0.8], [0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(10, 1)])
end

@testset "Stationnary Distributions" begin
    hmm1 = HMM([0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(10, 1)])
    hmm2 = HMM([1.0 0.0; 0.0 1.0], [Normal(0, 1), Normal(10, 1)])

    dists1 = statdists(hmm1)
    dists2 = statdists(hmm2)

    @test length(dists1) == 1
    @test length(dists2) == 2

    @test permutedims(dists2[1]) ≈ [1.0 0.0] * (hmm2.A^1000)
    @test permutedims(dists2[2]) ≈ [0.0 1.0] * (hmm2.A^1000)
end

@testset "Messages (1)" begin
    # Example from https://en.wikipedia.org/wiki/Forward%E2%80%93backward_algorithm
    A = [0.7 0.3; 0.3 0.7]
    B = [Categorical([0.9, 0.1]), Categorical([0.2, 0.8])]
    hmm = HMM(A, B)

    O = [1,1,2,1,1]

    α, logtot1 = forward(hmm, O)
    α = round.(α, digits = 4)

    β, logtot2 = backward(hmm, O)
    β = round.(β, digits = 4)

    γ = posteriors(hmm, O)
    γ = round.(γ, digits = 4)

    @test α == [
        0.8182 0.1818;
        0.8834 0.1166;
        0.1907 0.8093;
        0.7308 0.2692;
        0.8673 0.1327;
    ]

    @test β == [
        0.5923 0.4077;
        0.3763 0.6237;
        0.6533 0.3467;
        0.6273 0.3727;
        1.0    1.0;
    ]

    @test γ == [
        0.8673 0.1327;
        0.8204 0.1796;
        0.3075 0.6925;
        0.8204 0.1796;
        0.8673 0.1327;
    ]

    @test logtot1 ≈ logtot2
end

@testset "Messages (2)" begin
    hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(10, 1)])
    y = rand(hmm, 1000)

    α1, logtot1 = forward(hmm, y)
    α2, logtot2 = forward(hmm, y, logl = true)

    @test logtot1 ≈ logtot2
    @test α1 ≈ α2

    β1, logtot3 = backward(hmm, y)
    β2, logtot4 = backward(hmm, y, logl = true)

    @test logtot3 ≈ logtot4
    @test β1 ≈ β2

    @test size(α1) == size(α2) == size(β1) == size(β2)
    @test logtot1 ≈ logtot2 ≈ logtot3 ≈ logtot4

    γ1 = posteriors(hmm, y)
    γ2 = posteriors(hmm, y, logl = true)

    @test γ1 ≈ γ2
end

@testset "Messages (3)" begin
    hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0, 0), Normal(10, 0)])
    y = rand(hmm, 1000)

    # The likelihood of a Normal distribution with std. = 0
    # equals either 0 or +Inf (-Inf, +Inf in log domain).
    # This cause the forward/backward algorithms to return NaNs,
    # We mark the tests as broken, since I don't know if this can be fixed.
    # The workaround is to set `robust = true`.
    _, logtot1 = forward(hmm, y)
    _, logtot2 = forward(hmm, y, logl = true)
    _, logtot3 = backward(hmm, y)
    _, logtot4 = backward(hmm, y, logl = true)

    @test_broken !isnan(logtot1)
    @test_broken !isnan(logtot2)
    @test_broken !isnan(logtot3)
    @test_broken !isnan(logtot4)

    # See comment in viterbi.jl
    @test_broken @test_throws BoundsError viterbi(hmm, y)
    # Works in log implementation
    @test_nowarn viterbi(hmm, y, logl = true)

    _, logtot5 = forward(hmm, y, robust = true)
    _, logtot6 = forward(hmm, y, logl = true, robust = true)
    _, logtot7 = backward(hmm, y, robust = true)
    _, logtot8 = backward(hmm, y, logl = true, robust = true)

    @test !isnan(logtot5)
    @test !isnan(logtot6)
    @test !isnan(logtot7)
    @test !isnan(logtot8)

    @test logtot5 ≈ logtot6 ≈ logtot7 ≈ logtot8

    @test viterbi(hmm, y, robust = true) == viterbi(hmm, y, logl = true, robust = true)
end

@testset "Viterbi" begin
    hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(100, 1)])
    z, y = rand(hmm, 1000, seq = true)

    zv1 = viterbi(hmm, y)
    zv2 = viterbi(hmm, y, logl = true)

    @test zv1 == z
    @test zv1 == zv2
end

@testset "MLE" begin
    hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(10, 1)])
    y = rand(hmm, 1000)

    # Likelihood should not decrease
    _, hist = fit_mle(hmm, y)
    @test issorted(round.(hist.logtots, digits=9))

    _, hist = fit_mle(hmm, y, robust = true)
    @test issorted(round.(hist.logtots, digits=9))

    _, hist = fit_mle(hmm, y, maxiter = 0)
    @test hist.iterations == 0
    @test !hist.converged
end

@testset "Utilities (1)" begin
    # Make sure that we do not relabel the states if they are in 1...K
    mapping, _  = gettransmat([3,3,1,1,2,2], relabel = true)
    for (k, v) in mapping
        @test mapping[k] == v
    end

    mapping, transmat = gettransmat([3,3,8,8,3,3], relabel = true)
    @test mapping[3] == 1
    @test mapping[8] == 2
    @test transmat == [2 / 3 1 / 3; 1 / 2 1 / 2]

    transmat = randtransmat(10)
    @test issquare(transmat)
    @test istransmat(transmat)

    @test_logs (:warn, r".+") warn_logl(rand(10, 2) .- 10)
    @test_nowarn warn_logl(rand(10, 2))
end

@testset "Utilities (2)" begin
    ref  = [1,1,2,2,3,3]
    seq1 = [2,2,3,3,1,1]
    seq2 = [1,1,1,1,2,2]

    @test remapseq(seq1, ref) == ref
    @test remapseq(seq2, ref) == [1,1,1,1,3,3]
end

@testset "Reproducibility" begin
    hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0, 1), Normal(10, 1)])

    z1, y1 = rand(MersenneTwister(0), hmm, 1000, seq = true)
    z2, y2 = rand(MersenneTwister(0), hmm, 1000, seq = true)
    z3, y3 = rand(hmm, 1000, seq = true)
    @test z1 == z2 != z3
    @test y1 == y2 != y3

    A1 = randtransmat(MersenneTwister(0), Dirichlet(4, 1.0))
    A2 = randtransmat(MersenneTwister(0), Dirichlet(4, 1.0))
    A3 = randtransmat(Dirichlet(4, 1.0))
    @test A1 == A2 != A3

    A4 = randtransmat(MersenneTwister(0), 4)
    A5 = randtransmat(MersenneTwister(0), 4)
    A6 = randtransmat(4)
    @test A4 == A5 != A6
end

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase


text:
@graphs begin
	@node phase_1 = Start() -> p1_strat
	@node p1_strat = CheckStrategy("military") -> [n_true=p1_mili_1, n_false=p1_corr_1]

	# Military
	@node p1_mili_1 = PerformAction("Recover action dice.") -> p1_mili_2
	@node p1_mili_2 = PerformAction("Draw event cards.") -> p1_mili_3
	@node p1_mili_3 = BinaryCondition("Holding more than 6 cards.") -> [n_true = p1_mili_discard, n_false = p1_mili_end]
	@node p1_mili_discard = PerformAction("""
										  Discard event cards down to 6.

										  Priority:
										  1. Doesn't use the term "Fellowship revealed"
										  2. Character card
										  3. Strategy card
										  4. Descending order of initiative
										  5. Doesn't place a tile
										  6. Random
										  """) -> p1_mili_end
	@node p1_mili_end = End("End of Phase") -> []

	# Corruption
	@node p1_corr_1 = PerformAction("Recover action dice.") -> p1_corr_2
	@node p1_corr_2 = PerformAction("Draw event cards.") -> p1_corr_3
	@node p1_corr_3 = BinaryCondition("Holding more than 6 cards.") -> [n_true = p1_corr_discard, n_false = p1_corr_end_1]
	@node p1_corr_discard = BinaryCondition("Holding more than 1 strategy cards.") -> [n_true = p1_corr_discard_1, n_false = p1_corr_discard_2]
	@node p1_corr_discard_1 = PerformAction("""
											Discard event cards down to 6.

											Priority:
											1. Doesn't use the term "Fellowship revealed"
											2. Doesn't place a tile
											3. Strategy card
											4. Character card
											5. Descending order of initiative
											6. Random
											""") -> p1_corr_end_2
	@node p1_corr_discard_2 = PerformAction("""
											Discard event cards down to 6.

											Priority:
											1. Doesn't use the term "Fellowship revealed"
											2. Doesn't place a tile
											3. Character card
											4. Strategy card
											5. Descending order of initiative
											6. Random
											""") -> p1_corr_end_2
	@node p1_corr_end_1 = End("End of Phase") -> []
	@node p1_corr_end_2 = End("End of Phase") -> []
end

ext: jl

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
# * N-body operators

abstract type QuantumOperator end

"""
    ishermitian(op::QuantumOperator)

By default, all `QuantumOperator`s are Hermitian; this can be
overridden for subtypes to explicitly declare an operator
non-Hermitian.
"""
LinearAlgebra.ishermitian(::QuantumOperator) = true

"""
    NBodyOperator{N}

Abstract N-body operator coupling `N` bodies each between two Slater determinants
"""
abstract type NBodyOperator{N} <: QuantumOperator end

"""
    numbodies(::NBodyOperator{N})

Returns the number of bodies coupled by the N-body operator, i.e. `N`.
"""
numbodies(::NBodyOperator{N}) where N = N

# # Examples

# ```jldoctest
# julia> ZeroBodyOperator
# Ω₀

# julia> OneBodyOperator
# Ω₁

# julia> TwoBodyOperator
# Ω₂

# julia> NBodyOperator{6}
# Ω₆
# ```

# Base.show(io::IO, ::Type{NBodyOperator{N}}) where N =
#     write(io, "Ω", to_subscript(N))

const ZeroBodyOperator = NBodyOperator{0}
const OneBodyOperator = NBodyOperator{1}
const TwoBodyOperator = NBodyOperator{2}

"""
    IdentityOperator{N}

The N-body identity operator. Leaves the orbital(s) acted upon
unchanged.
"""
struct IdentityOperator{N} <: NBodyOperator{N} end

Base.show(io::IO, ::IdentityOperator{N}) where N =
    write(io, to_boldface("I"), to_subscript(N))

# * Linear combination

"""
    LinearCombinationOperator(operators)

Represents a linear combination of [`NBodyOperator`](@ref)s.
"""
struct LinearCombinationOperator <: QuantumOperator
    operators::Vector{<:Pair{<:NBodyOperator,<:Number}}
end

Base.zero(::LinearCombinationOperator) =
    LinearCombinationOperator(Pair{<:NBodyOperator,<:Number}[])
Base.zero(::Type{LinearCombinationOperator}) =
    LinearCombinationOperator(Pair{<:NBodyOperator,<:Number}[])

Base.iszero(op::LinearCombinationOperator) =
    isempty(op.operators)

Base.:(==)(a::LinearCombinationOperator,b::LinearCombinationOperator) =
    a.operators == b.operators

function Base.hash(lco::LinearCombinationOperator, h::UInt)
    for (op,coeff) in lco.operators
        h = hash(op, hash(coeff, h))
    end
    h
end

"""
    numbodies(lco::LinearCombinationOperator)

Returns the maximum number of bodies coupled by any of the N-body
operators in the [`LinearCombinationOperator`](@ref).
"""
numbodies(lco::LinearCombinationOperator) =
    iszero(lco) ? 0 : maximum(numbodies.(first.(lco.operators)))

function Base.show(io::IO, op::LinearCombinationOperator)
    if iszero(op)
        write(io, "0")
        return
    end
    for (i,(o,n)) in enumerate(op.operators)
        i > 1 && write(io, " ")
        showcoeff(io, n, i > 1)
        write(io, string(o))
    end
end

Base.:(+)(a::NBodyOperator, b::NBodyOperator) =
    LinearCombinationOperator([a=>1,b=>1])
Base.:(-)(a::NBodyOperator, b::NBodyOperator) =
    LinearCombinationOperator([a=>1,b=>-1])

Base.:(+)(a::LinearCombinationOperator, b::NBodyOperator) =
    LinearCombinationOperator(vcat(a.operators, b=>1))
Base.:(-)(a::LinearCombinationOperator, b::NBodyOperator) =
    LinearCombinationOperator(vcat(a.operators, b=>-1))

Base.:(+)(a::NBodyOperator, b::LinearCombinationOperator) =
    LinearCombinationOperator(vcat(a=>1, b.operators))
Base.:(-)(a::NBodyOperator, b::LinearCombinationOperator) =
    LinearCombinationOperator(vcat(a=>1, Pair{<:NBodyOperator,<:Number}[o[1]=>-o[2]
                                                                       for o in b.operators]))

Base.:(+)(a::LinearCombinationOperator, b::LinearCombinationOperator) =
    LinearCombinationOperator(vcat(a.operators, b.operators))
Base.:(-)(a::LinearCombinationOperator, b::LinearCombinationOperator) =
    LinearCombinationOperator(vcat(a.operators, Pair{<:NBodyOperator,<:Number}[o[1]=>-o[2]
                                                                               for o in b.operators]))

Base.:(*)(a::Number, b::NBodyOperator) =
    LinearCombinationOperator([b=>a])
Base.:(*)(a::NBodyOperator, b::Number) =
    LinearCombinationOperator([a=>b])

Base.:(*)(a::Number, b::LinearCombinationOperator) =
    LinearCombinationOperator(Pair{<:NBodyOperator,<:Number}[op=>a*c for (op,c) in b.operators])
Base.:(*)(a::LinearCombinationOperator, b::Number) =
    LinearCombinationOperator(Pair{<:NBodyOperator,<:Number}[op=>b*c for (op,c) in a.operators])

Base.filter(fun::Function, op::LinearCombinationOperator) =
    LinearCombinationOperator(filter(oc -> fun(oc[1]), op.operators))

# * Contracted operators

"""
    ContractedOperator(a, o, b)

An [`NBodyOperator`](@ref) representing the contraction of the
operator `o` over the orbital sets `a` and `b`. The lengths of `a` and
`b` have to equal, and they cannot exceed the dimension of `o`.
"""
struct ContractedOperator{N,P,Q,A,O<:NBodyOperator{P},B} <: NBodyOperator{N}
    a::NTuple{Q,A}
    o::O
    b::NTuple{Q,B}
    function ContractedOperator(a::NTuple{Q,A}, o::O, b::NTuple{Q,B}) where {P,Q,A,O<:NBodyOperator{P},B}
        P ≥ Q || throw(ArgumentError("Cannot contract $(P)-body operator over $(Q) coordinates"))
        N = P - Q
        new{N,P,Q,A,O,B}(a, o, b)
    end
end

Base.:(==)(a::ContractedOperator, b::ContractedOperator) =
    a.a == b.a && a.o == b.o && a.b == b.b

"""
    in(orbital, co::ContractedOperator)

Test if `orbital` is among the right set of orbitals of the
[`ContractedOperator`](@ref) `co`. Useful to test if `co` is an
integral operator with respect to `orbital`.
"""
Base.in(orbital::O, co::ContractedOperator) where O = orbital ∈ co.b
"""
    in(corbital::Conjugate, co::ContractedOperator)

Test if `corbital` is among the left set of orbitals of the
[`ContractedOperator`](@ref) `co`. Useful to test if `co` is an
integral operator with respect to `corbital`.
"""
Base.in(corbital::Conjugate{O}, co::ContractedOperator) where O = corbital.orbital ∈ co.a

"""
    contract(orbital_matrix_element, i...)

Contract the `orbital_matrix_element` over all coordinates `i...`.
"""
contract(ome, i::Integer...) =
    throw(ArgumentError("`contract` not implemented for $(ome)"))

"""
    complement(N, i...)

Generate the complement to `i...` in the set `1:N`. Useful for
contracting [`OrbitalMatrixElement`](@ref)s over all coordinates
_except_ `i...`.
"""
complement(N::Integer, i::Integer...) = setdiff(1:N, i)
complement(N::Integer, i::AbstractVector{<:Integer}) = setdiff(1:N, i)

function Base.show(io::IO, o::ContractedOperator{N}) where N
    write(io, "⟨", join(string.(o.a), " "))
    write(io, "|")
    show(io, o.o)
    write(io, "|")
    write(io, join(string.(o.b), " "), "⟩")
end

export QuantumOperator, NBodyOperator, numbodies,
    ZeroBodyOperator, OneBodyOperator, TwoBodyOperator,
    IdentityOperator, ContractedOperator


text:
function BNB_control(user, balances_all;
        settings_path=joinpath(Pinguu.base_dir, "BNB_control/BNB_settings.json"))
    acc = JSON.parsefile(settings_path)
    prices = Dict()
    for (user_id, v) in user
        if haskey(acc, user_id)
            client = user[user_id]["client"]
            if !haskey(acc[user_id], "order")
                acc[user_id]["order"] = Dict()
                open(settings_path, "w") do f
                    JSON.print(f, acc, 4)
                end
            end
            if acc[user_id]["active"] == 1
                @info("---")
                @info("User: " * user_id)
                sleep(0.25)
                if !haskey(prices, acc[user_id]["symbol"])
                    ticker = client.get_ticker(symbol=acc[user_id]["symbol"])["lastPrice"]
                    price = parse(Float64, ticker)
                    prices[acc[user_id]["symbol"]] = Dict()
                    prices[acc[user_id]["symbol"]]["price"] = price
                    prices[acc[user_id]["symbol"]]["ticker"] = ticker
                else
                    price  = prices[acc[user_id]["symbol"]]["price"]
                    ticker = prices[acc[user_id]["symbol"]]["ticker"]
                end

                if balances_all[user_id]["BNB"]["estimated_free"] * price < acc[user_id]["low_limit"] && length(acc[user_id]["order"]) == 0
                    #
                    ##
                    ###
                    # If BNB level is low & no order ongoing => Check if enough base money is available
                    if balances_all[user_id][acc[user_id]["base_coin"]]["estimated_free"] <= acc[user_id]["val"]
                        @info(string("Base volume not sufficient: ", round(balances_all[user_id][acc[user_id]["base_coin"]]["estimated_free"], digits=2), " ", acc[user_id]["base_coin"]))
                        continue
                    end
                    qty = ceil(acc[user_id]["val"]/price, digits=3)
                    sleep(0.5)
                    result = try client.order_market_buy(
                        symbol=acc[user_id]["symbol"],
                        quantity=qty);
                        catch
                        @info("order_market_buy failed")
                        nothing
                    end
                    if result !== nothing
                        acc[user_id]["order"] = result
                        open(settings_path, "w") do f
                            JSON.print(f, acc, 4)
                        end
                    end
                    if user_id == "1"
                        telegram = JSON.parsefile(joinpath(Pinguu.base_dir, "telegram.json"))
                        msg = "[BNB Control] Kauf getätigt!"
                        TelegramClient(Pinguu.telegram; chat_id=telegram[user_id])
                        sendMessage(text = msg)
                    end
                elseif length(acc[user_id]["order"]) != 0
                    #
                    ##
                    ###
                    # Check ongoing order
                    order = try client.get_order(
                        symbol=acc[user_id]["symbol"],
                        orderId=acc[user_id]["order"]["orderId"])
                    catch
                        @info("get_order failed")
                        nothing
                    end
                    if order !== nothing
                        if order["status"] == client.ORDER_STATUS_FILLED
                            yr,m,d = split(string(today()), "-")
                            filename = joinpath(Pinguu.base_dir * "BNB_control/users/" * user_id * "/closed_orders/" * yr * "/" * m * "/", "closed_" * string(now()) * ".json")
                            !isdir(dirname(filename)) ? mkpath(dirname(filename)) : ""
                            open(filename, "w") do f
                                JSON.print(f, order, 4)
                            end
                            acc[user_id]["order"] = Dict()
                            open(settings_path, "w") do f
                                JSON.print(f, acc, 4)
                            end
                            if user_id == "1"
                                telegram = JSON.parsefile(joinpath(Pinguu.base_dir, "telegram.json"))
                                msg = "[BNB Control] Kauf abgeschlossen!"
                                TelegramClient(Pinguu.telegram; chat_id=telegram[user_id])
                                sendMessage(text = msg)
                            end
                        else
                            @info("Order status: " * order["status"])
                        end
                    end
                else
                    #
                    ##
                    ###
                    # Nothing to do
                    @info(string("BNB level: ", round(balances_all[user_id]["BNB"]["estimated_free"] * price, digits=2), " ", acc[user_id]["base_coin"]))
                    @info(string("BNB low limit: ", round(acc[user_id]["low_limit"], digits=2), " ", acc[user_id]["base_coin"]))
                end
            end
        end
    end
end
export BNB_control


text:
function InitCoords(configuration::Configuration,balancingMPI::BalancingMPI)
    n_total=sum(configuration.moles_length)
    coords=zeros(n_total,3)
    index=1
    for comm_mole_i in balancingMPI.mole_index_comms
        for mole_i in comm_mole_i
            molecule=configuration.moles[mole_i]
            mole_length=size(molecule.coords)[1]
            coords[index:mole_length+index-1,1:3]=molecule.coords
            index=index+mole_length
        end
    end
    @info "Initialize Coordinates Size $(n_total)"
    return coords
end

function InitMasses(configuration::Configuration,balancingMPI::BalancingMPI)
    n_total=sum(configuration.moles_length)
    masses=zeros(n_total)
    index=1
    for comm_mole_i in balancingMPI.mole_index_comms
        for mole_i in comm_mole_i
            molecule=configuration.moles[mole_i]
            mole_length=size(molecule.coords)[1]
            masses[index:mole_length+index-1]=molecule.mass
            index=index+mole_length
        end
    end
    @info "Initialize Masses Size $(n_total)"
    return masses
end

function InitTypes(configuration::Configuration,balancingMPI::BalancingMPI)
    n_total=sum(configuration.moles_length)
    types=zeros(Int64,n_total)
    index=1
    for comm_mole_i in balancingMPI.mole_index_comms
        for mole_i in comm_mole_i
            molecule=configuration.moles[mole_i]
            mole_length=size(molecule.coords)[1]
            types[index:mole_length+index-1]=molecule.type
            index=index+mole_length
        end
    end
    @info "Initialize Types Size $(n_total)"
    return types
end

function InitBonds(configuration::Configuration,balancingMPI::BalancingMPI)
    n_total=0
    for molei in configuration.moles
        n_total+=size(molei.bonds)[1]
    end
    bonds=zeros(Int64,n_total,3)
    index=1
    for comm_mole_i in balancingMPI.mole_index_comms
        total_atom_index=0
        for mole_i in comm_mole_i
            molecule=configuration.moles[mole_i]
            bond_length=size(molecule.bonds)[1]
            mole_length=size(molecule.coords)[1]
            if bond_length!=0
                bonds[index:bond_length+index-1,1]=molecule.bonds[:,1]
                bonds[index:bond_length+index-1,2:3]=molecule.bonds[:,2:3].+total_atom_index
            end
            total_atom_index+=mole_length
            index=index+bond_length
        end
    end
    @info "Initialize Bonds Size $(n_total)"
    return bonds
end

function InitVels(configuration::Configuration,balancingMPI::BalancingMPI,velocity::ZeroVelocity)
    n_total=sum(configuration.moles_length)
    vels=zeros(n_total,3)
    @info "Initialize Velocities Size $(n_total)"
    return vels
end

function InitVels(configuration::Configuration,balancingMPI::BalancingMPI,velocity::NoVelocity)
    n_total=sum(configuration.moles_length)
    vels=zeros(n_total,3)
    index=1
    for comm_mole_i in balancingMPI.mole_index_comms
        for mole_i in comm_mole_i
            molecule=configuration.moles[mole_i]
            mole_length=size(molecule.vels)[1]
            vels[index:mole_length+index-1,1:3]=molecule.vels
            index=index+mole_length
        end
    end
    @info "Initialize Velocities Size $(n_total)"
    return vels
end

function InitForces(configuration::Configuration)
    n_total=sum(configuration.moles_length)
    forces=zeros(n_total,3)
    @info "Initialize Forces Size $(n_total)"
    return forces
end

function InitEnergies(configuration::Configuration)
    n_total=sum(configuration.moles_length)
    energies=zeros(n_total)
    @info "Initialize Energy Size $(n_total)"
    return energies
end


text:
#= ======= FUNCTIONS =============================================== =#
function printBoard(board :: Array{Int, 2}) :: Void
    show(IOContext(STDOUT, limit=true), "text/plain", map((x -> (x == 0) ? '_'  :  ((x == 1) ? 'X' : 'O')), board))
    println()
end

# Board and Last Played Column -> isGameActive and Winner
function winState(board :: Array{Int, 2}, lastCol :: Int)
    # Look to see if board is filled
    ! (0 in board) && return true, 0
    (board == zeros(board)) && return false, 0

    lastRow = findlast(board[:,lastCol], 0) + 1  # Last row used
    player = board[lastRow,lastCol]  # Player in lastCol,lastRow

    # Search down for verticle win
    if lastRow <= 4 && all( x -> x == player, board[lastRow:(lastRow + 4)])
        return true, player
    end

    isWin = false

    # Fluff Matrix time
    # Extend matrix by two on all sides
    # Embed board in large board
    fluff = zeros(Int,(12,13))
    fluff[4:9,4:10] = board[:]

    # Shift location
    lastRow += 3
    lastCol += 3

    # Slice possible win paths
    downDiagLists = [[fluff[lastRow - 3 + i + j, lastCol - 3 + i + j] for i in 0:3] for j in 0:3]
    horizLists = [fluff[lastRow,(lastCol-3+i):(lastCol + i)] for i in 0:3]
    upDiagLists = [[fluff[lastRow + 3 - i - j, lastCol - 3 + i + j] for i in 0:3] for j in 0:3]
    lists = vcat(vcat(downDiagLists, horizLists), upDiagLists)

    if any(x -> all(y -> y == player, x), lists)
        return true, player
    end

    return false, 0
end

function place!(board :: Array{Int, 2}, player :: Int, colNum :: Int) :: Void
    column = board[:,colNum]
    @assert 0 in column "Error: Tried to place outside the board."
    board[findlast(column, 0), colNum] = player
    return nothing
end

#=
randomChoice( board )
Takes a board as a matrix,
returns a randod valid move
=#
function randomChoice(board :: Array{Int, 2}) :: Int
    return rand([col for col in 1:size(board,2) if 0 in board[:,col]])
end

#= ================================================================= =#

# MAIN
# Initialize the board
board = Array{Int, 2}(6, 7)
board = zeros(board)

#Initialize the game
player1 = 1
player2 = -1
currentPlayer = player1
lastMove = 0
printBoard(board)

while !winState(board, lastMove)[1]
    lastMove = randomChoice(board)
    place!(board, currentPlayer, lastMove)
    printBoard(board)
    currentPlayer *= -1
end

println(currentPlayer != 1 ? 'X' : 'O', " wins")

max_stars_repo_name: TheGluck/Connect-Four-AI

max_issues_repo_name: TheGluck/Connect-Four-AI

max_forks_repo_name: TheGluck/Connect-Four-AI


text:
push!(LOAD_PATH, dirname(@__DIR__))
using EnergySystemModeling, Parameters, JSON

output = "output"

parameters = load_json(EnergySystemModeling.Params, joinpath(output, "parameters.json"))
variables = load_json(EnergySystemModeling.Variables, joinpath(output, "variables.json"))
objectives = load_json(EnergySystemModeling.Objectives, joinpath(output, "objectives.json"))
expressions = Expressions(parameters, variables)

using Plots
using StatsPlots
pyplot()

savefig(plot_objective_values(objectives),
        joinpath(output, "objectives.svg"))

for n in parameters.N
    savefig(plot_generation_dispatch(parameters, variables, expressions, n),
            joinpath(output, "generation_dispatch_n$n.svg"))
    savefig(plot_generation_capacities(parameters, variables, expressions, n),
            joinpath(output, "generation_capacities_n$n.svg"))
    savefig(plot_storage_level(parameters, variables, expressions, n),
            joinpath(output, "storage_n$n.svg"))
    savefig(plot_box(parameters, variables, expressions, n),
           joinpath(output, "boxplot$n.svg"))

end
savefig(plot_storage_capacities(parameters, variables, expressions),
        joinpath(output, "storage_capacities.svg"))

savefig(plot_box_all(parameters, variables, expressions),
        joinpath(output, "boxplotall.svg"))

savefig(plot_dispatch_bars(parameters, variables, expressions),
        joinpath(output, "dispatchbars.svg"))

for l in 1:length(parameters.L)
    savefig(plot_transmission_flow(parameters, variables, expressions, l),
            joinpath(output, "transmission_flow_l$l.svg"))
end

savefig(plot_transmission_capacities(parameters, variables, expressions),
        joinpath(output, "transmission_capacities.svg"))

savefig(plot_transmission_bars(parameters, variables, expressions),
        joinpath(output, "transmission_bars.svg"))

savefig(plot_loss_of_load(parameters, variables, expressions),
        joinpath(output, "loss_of_load.svg"))


text:
"""
    SSIM([kernel], [(α, β, γ)]; crop=false) <: FullReferenceIQI
    assess(iqi::SSIM, img, ref)
    assess_ssim(img, ref; crop=false)

Structural similarity (SSIM) index is an image quality assessment method based
on degradation of structural information.

The SSIM index is composed of three components: luminance, contrast, and
structure; `ssim = 𝐿ᵅ * 𝐶ᵝ * 𝑆ᵞ`, where `W := (α, β, γ)` controls relative
importance of each components. By default `W = (1.0, 1.0, 1.0)`.

In practice, a mean version SSIM is used. At each pixel, SSIM is calculated
locally with neighborhoods weighted by `kernel`, returning a ssim map;
`ssim` is actaully `mean(ssim_map)`.
By default `kernel = KernelFactors.gaussian(1.5, 11)`.

!!! tip

    The default parameters comes from [1]. For benchmark usage, it is recommended to not
    change the parameters, because most other SSIM implementations follows the same settings.
    Keyword `crop` controls whether the boundary of ssim map should be dropped; you need
    to set it `true` to reproduce the result in [1].

# Example

`assess_ssim(img, ref)` should be sufficient to get a benchmark for algorithms. One
could also instantiate a customed SSIM, then pass it to `assess` or use it as a
function. For example:

```julia
iqi = SSIM(KernelFactors.gaussian(2.5, 17), (1.0, 1.0, 2.0))
assess(iqi, img, ref)
iqi(img, ref) # both usages are equivalent
```

!!! info

    SSIM is defined only for gray images. How `RGB` and other `Color3` images are handled may vary
    in different implementations. For `RGB` images, channels are handled seperately when calculating
    𝐿, 𝐶, and 𝑆. Generic `Color3` images are converted to `RGB` first before calculation.

# Reference

[1] Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P. (2004). Image quality assessment: from error visibility to structural similarity. _IEEE transactions on image processing_, 13(4), 600-612.

[2] Wang, Z., Bovik, A. C., Sheikh, H. R., & Simoncelli, E. P. (2003). The SSIM Index for Image Quality Assessment. Retrived May 30, 2019, from http://www.cns.nyu.edu/~lcv/ssim/
"""
struct SSIM{A<:AbstractVector} <: FullReferenceIQI
    kernel::A
    W::NTuple{3, Float64}
    crop::Bool
    function SSIM(kernel::AbstractVector=SSIM_KERNEL, W::NTuple=SSIM_W; crop=false)
        ndims(kernel) == 1 || throw(ArgumentError("only 1-d kernel is valid"))
        issymetric(kernel) || @warn "SSIM kernel is assumed to be symmetric"
        all(W .>= 0) || throw(ArgumentError("(α, β, γ) should be non-negative, instead it's $(W)"))
        kernel = centered(kernel)
        new{typeof(kernel)}(kernel, W, crop)
    end
end

# default values from [1]
const SSIM_KERNEL = KernelFactors.gaussian(1.5, 11) # kernel
const SSIM_W = (1.0, 1.0, 1.0) # (α, β, γ)

Base.:(==)(ia::SSIM, ib::SSIM) = ia.kernel == ib.kernel && ia.W == ib.W && ia.crop == ib.crop

# api
# By default we don't crop the padding boundary to match the ssim result from
# MATLAB Image Processing Toolbox, which is used more broadly than the original
# implementaion [2] (written in MATLAB as well).
# -- Johnny Chen <johnnychen94@hotmail.com>
(iqi::SSIM)(x, ref) = mean(_ssim_map(iqi, x, ref))

@doc (@doc SSIM)
assess_ssim(x, ref; crop=false) = SSIM(;crop=crop)(x, ref)

# Parameters `(K₁, K₂)` are used to avoid instability when denominator is very
# close to zero. Different from origianl implementation [2], we don't make it
# public since the ssim result is insensitive to these parameters according to
# [1].
# -- Johnny Chen <johnnychen94@hotmail.com>
const SSIM_K = (0.01, 0.03)

# SSIM is defined only for gray images,
# RGB images are treated as 3d gray images,
# other Color3 images are converted to RGB first.
function _ssim_map(iqi::SSIM,
                   x::Union{GenericGrayImage, AbstractArray{<:AbstractRGB}},
                   ref::Union{GenericGrayImage, AbstractArray{<:AbstractRGB}},
                   peakval = 1.0,
                   K = SSIM_K)
    if size(x) ≠ size(ref)
        err = ArgumentError("images should be the same size, instead they're $(size(x))-$(size(ref))")
        throw(err)
    end
    if axes(x) ≠ axes(ref)
        x = OffsetArrays.no_offset_view(x)
        ref = OffsetArrays.no_offset_view(ref)
    end
    α, β, γ = iqi.W
    C₁, C₂ = @. (peakval * K)^2
    C₃ = C₂/2

    T = promote_type(float(eltype(ref)), float(eltype(x)))
    x = of_eltype(T, x)
    ref = of_eltype(T, ref)

    if [α, β, γ] ≈ [1.0, 1.0, 1.0]
        ssim_map = __ssim_map_fast(x, ref, iqi.kernel, C₁, C₂; crop=iqi.crop)
    else
        l, c, s = __ssim_map_general(x, ref, iqi.kernel, C₁, C₂, C₃; crop=iqi.crop)
        # ensure that negative numbers in s are not being raised to powers less than 1
        # this is a non-standard implementation and it could vary from implementaion to implementaion
        γ < 1.0 && (s = max.(s, zero(eltype(s))))
        ssim_map = @. l^α * c^β * s^γ # equation (12) in [1]
    end
    return ssim_map
end

function _ssim_map(iqi::SSIM,
          x::AbstractArray{<:Color3},
          ref::AbstractArray{<:Color3},
          peakval = 1.0, K = SSIM_K)
    C = ccolor(RGB, floattype(eltype(x)))
    _ssim_map(iqi, of_eltype(C, x), of_eltype(C, ref), peakval, K)
end

# helpers
function issymetric(kernel)
    origin = first(axes(kernel, 1))
    center = (length(kernel)-1) ÷ 2 + origin
    kernel[origin:center] ≈ kernel[end:-1:center]
end

function __ssim_map_fast(x, ref, window, C₁, C₂; crop)
    μx², μxy, μy², σx², σxy, σy² = _ssim_statistics(x, ref, window; crop=crop)
    # this is a special case when [α, β, γ] ≈ [1.0, 1.0, 1.0]
    # equation (13) in [1]
    return @. ((2μxy + C₁)*(2σxy + C₂))/((μx²+μy²+C₁)*(σx² + σy² + C₂))
end

function __ssim_map_general(x, ref, window, C₁, C₂, C₃; crop)
    μx², μxy, μy², σx², σxy, σy² = _ssim_statistics(x, ref, window; crop = crop)

    T = eltype(σx²)
    σx_σy = @. sqrt( max(σx²*σy², zero(T) ))
    l = @. (2μxy + C₁)/(μx² + μy²) # equation (6) in [1]
    c = @. (2σx_σy + C₂)/(σx² + σy² + C₂) # equation (9) in [1]
    s = @. (σxy + C₃)/(σx_σy + C₃) # equation (10) in [1]

    # MS-SSIM needs these, so we don't multiply them together here
    return l, c, s
end

function _ssim_statistics(x::GenericImage, ref::GenericImage, window; crop)
    # For RGB and other Color3 images, we don't slide the window at the color channel.
    # In other words, these characters will be calculated channelwisely
    window = kernelfactors(Tuple(repeated(window, ndims(ref))))

    region = map(window, axes(x)) do w, a
        o = length(w) ÷ 2
        # Even if crop=true, it crops only when image is larger than window
        length(a) > length(w) ? (first(a)+o:last(a)-o) : a
    end
    R = crop ? CartesianIndices(region) : CartesianIndices(x)

    # don't slide the window in the channel dimension
    μx = view(imfilter(x,   window, "symmetric"), R) # equation (14) in [1]
    μy = view(imfilter(ref, window, "symmetric"), R) # equation (14) in [1]
    μx² = _mul.(μx, μx)
    μy² = _mul.(μy, μy)
    μxy = _mul.(μx, μy)
    σx² = view(imfilter(_mul.(x,   x  ), window, "symmetric"), R) .- μx² # equation (15) in [1]
    σy² = view(imfilter(_mul.(ref, ref), window, "symmetric"), R) .- μy² # equation (15) in [1]
    σxy = view(imfilter(_mul.(x,   ref), window, "symmetric"), R) .- μxy # equation (16) in [1]

    # after that, channel dimension can be treated generically so we expand them here
    return channelview.((μx², μxy, μy², σx², σxy, σy²))
end

# _ssim_statistics does not expected to work with other Color3 types
@inline _mul(x::Number, y::Number) = x * y
@inline _mul(x::C, y::C) where C <: Union{AbstractGray, AbstractRGB} = mapc(*, x, y)


text:
function FieldLUDecomposition(arg0::FieldMatrix)
    return FieldLUDecomposition((FieldMatrix,), arg0)
end

function get_determinant(obj::FieldLUDecomposition)
    return jcall(obj, "getDeterminant", FieldElement, ())
end

function get_l(obj::FieldLUDecomposition)
    return jcall(obj, "getL", FieldMatrix, ())
end

function get_p(obj::FieldLUDecomposition)
    return jcall(obj, "getP", FieldMatrix, ())
end

function get_pivot(obj::FieldLUDecomposition)
    return jcall(obj, "getPivot", Vector{jint}, ())
end

function get_solver(obj::FieldLUDecomposition)
    return jcall(obj, "getSolver", FieldDecompositionSolver, ())
end

function get_u(obj::FieldLUDecomposition)
    return jcall(obj, "getU", FieldMatrix, ())
end


text:
using KnetLayers

"""
julia> predict(model,[3 2 1 4 5 9 3 5 6 6 1 2 5;])
1×14 Array{Int64,2}:
 1  2  2  2  3  4  4  5  5  5  6  7  9  11
"""
struct S2S; encoder; decoder; output; loss; end

# Initialize model
model = S2S(LSTM(input=11,hidden=128,embed=9),
            LSTM(input=11,hidden=128,embed=9),
            Multiply(input=128,output=11),
            CrossEntropyLoss())

# Forward functions
(m::S2S)(x) = m.output(m.decoder(pad(10,sort(x,dims=2)), m.encoder(x; hy=true).hidden).y)
(m::S2S)(x,ygold) = m.loss(m(x),ygold)
predict(m::S2S,x) = getindex.(argmax(Array(m(x)), dims=1)[1,:,:], 1)

# Helper functions for padding
pad(p::Int,x::Array; dims=2) = cat(fill(p,size(x,1)),x;dims=dims)
pad(x::Array,p::Int; dims=2) = cat(x,fill(p,size(x,1));dims=dims)

# Create sorting data: 10 is used as start token and 11 is stop token.
dataxy(x) = (x,pad(sort(x, dims=2), 11))
batchSize, maxLength, dataLength= 64, 15, 2000; # Batch size and maximum sequence length for training
data = [dataxy([rand(1:9) for j=1:batchSize, k=1:rand(1:maxLength)]) for i=1:dataLength]

# Train your model
progress!(adam(model,data))
# Test the model
@show predict(model,[3 2 1 4 5 9 3 5 6 6 1 2 5;])


text:
AVOA = RecursiveArrayTools.AbstractVectorOfArray

function Base.Array(VA::AVOA{T,N,A}) where {T,N,A<:AbstractVector{<:ComponentVector}}
    return ComponentArray(reduce(hcat, VA.u), only(getaxes(VA.u[1])), FlatAxis())
end


text:
module SetFloatParameterTrigger

using ..Ahorn, Maple

@mapdef Trigger "ShaderHelper/SetFloatParameterTrigger" SetFloatParameterTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, shader::String="shaderhelper/testshader", value::Float=0, key::String="key")

const placements = Ahorn.PlacementDict(
    "Set Float Parameter Trigger (Shader Helper)" => Ahorn.EntityPlacement(
        SetFloatParameterTrigger,
        "rectangle"
    )
)

end

max_stars_repo_name: lordseanington/ShaderHelper

max_issues_repo_name: lordseanington/ShaderHelper

max_forks_repo_name: lordseanington/ShaderHelper


text:
ENV["JUPYTER"]="/home/idies/miniconda3/bin/jupyter"
using Pkg
Pkg.add("IJulia")
Pkg.add("Plots")
Pkg.add("PyPlot")
Pkg.add("PyCall")

max_stars_repo_name: sciserver/sciserver-base-images

max_issues_repo_name: sciserver/sciserver-base-images

max_forks_repo_name: sciserver/sciserver-base-images


text:
# using Base.Test
function testdata(T_)
    # srand(1)

    n,m      = 2,1
    At_      = [0.95 0.1; 0 0.95]
    Bt_      = reshape([0.2; 1],2,1)
    u        = randn(1,T_)
    x        = zeros(n,T_)
    for t = 1:T_-1
        if t == T_÷2
            At_ = [0.5 0.05; 0 0.5]
        end
        x[:,t+1] = At_*x[:,t] + Bt_*u[:,t] + 0.2randn(n)
    end
    xm = x + 0.2randn(size(x));
    x',xm',u',n,m
end

T = 1000
x,xm,u,n,m = testdata(T)

nc = 2

model = LPVSS(x, u, nc; normalize=true, λ = 1e-3)
xh = model(x,u)

eRMS = √(mean((xh[1:end-1,:]-x[2:end,:]).^2))
println("RMS error: ", eRMS)
@test eRMS <= 0.37

# using Plots
# plotlyjs(show=true)
# plot(xh[1:end-1,:], lab="Prediction", c=:red, layout=(2,1))
# plot!(x[2:end,:], lab="True", c=:blue)

nc = 4
v = 1:T
model = LPVSS(x, u, v, nc; normalize=true, λ = 1e-3)
xh = model(x,u,v)

eRMS = √(mean((xh[1:end-1,:]-x[2:end,:]).^2))
println("RMS error: ", eRMS)
@test eRMS <= 0.28

# using Plots
# plotlyjs(show=true)
# plot(xh[1:end-1,:], lab="Prediction", c=:red, layout=(2,1))
# plot!(x[2:end,:], lab="True", c=:blue)


text:
module PartialBase

include("./../../constraints/geometric/constants.jl")
include("./1_if.jl")
include("./2_xy_cover.jl")
include("./3_multiplication.jl")
using JuMP
using .Constants
using .IfPartialBase
using .XYCoverPartialBase
using .MultiplicationPartialBase

export cons_ss_partial_base

function cons_ss_partial_base(model)

    model = if_pb(model)
    model = xy_cover(model)
    model = multiplication(model)

    # partial base stability constraint
    @constraint(model, c_pbs[i=1:(length(boxes)-1), j=1:length(boxes)-1; i!=j],
        model[:pbs_box_i][i] >=
        pb_ss_var*area_of_box[i]
    )

    return model
end

end


text:
const CLIENT_NAME = "SCP"

const VIRES_MAGIC_NUMBER = @compat UInt16(40108)
const VIRES_VERSION_HEADER = @compat UInt16(1)
const VIRES_HEADER_STRING_LEN = 64
const VIRES_HEADER_LEN = 2*2 + 2*VIRES_HEADER_STRING_LEN + 4

const NULL = zero(UInt8)

const TIMEOUT_DEFAULT = 2.0 # [s]
const TIMEOUT_LONG_DEFAULT = 10.0 # [s]

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive


text:
using SpecialFunctions
import Base.Broadcast

const monadic_linear = [deg2rad, +, rad2deg, transpose, -, conj]

const monadic_nonlinear = [asind, log1p, acsch, erfc, digamma, acos, asec, acosh, airybiprime, acsc, cscd, log, tand, log10, csch, asinh, airyai, abs2, gamma, lgamma, erfcx, bessely0, cosh, sin, cos, atan, cospi, cbrt, acosd, bessely1, acoth, erfcinv, erf, dawson, inv, acotd, airyaiprime, erfinv, trigamma, asecd, besselj1, exp, acot, sqrt, sind, sinpi, asech, log2, tan, invdigamma, airybi, exp10, sech, erfi, coth, asin, cotd, cosd, sinh, abs, besselj0, csc, tanh, secd, atand, sec, acscd, cot, exp2, expm1, atanh]

diadic_of_linearity(::Val{(true, true, true)}) = [+, rem2pi, -, >, isless, <, isequal, max, min, convert]
diadic_of_linearity(::Val{(true, true, false)}) = [*]
#diadic_of_linearit(::(Val{(true, false, true)}) = [besselk, hankelh2, bessely, besselj, besseli, polygamma, hankelh1]
diadic_of_linearity(::Val{(true, false, false)}) = [/]
diadic_of_linearity(::Val{(false, true, false)}) = [\]
diadic_of_linearity(::Val{(false, false, false)}) = [hypot, atan, mod, rem, lbeta, ^, beta]
diadic_of_linearity(::Val) = []

haslinearity(f, nargs) = false

# linearity of a single input function is either
# Val{true}() or Val{false}()
#
for f in monadic_linear
    @eval begin
        haslinearity(::typeof($f), ::Val{1}) = true
        linearity(::typeof($f), ::Val{1}) = Val{true}()
    end
end
# linearity of a 2-arg function is:
# Val{(linear11, linear22, linear12)}()
#
# linearIJ refers to the zeroness of d^2/dxIxJ
for f in monadic_nonlinear
    @eval begin
        haslinearity(::typeof($f), ::Val{1}) = true
        linearity(::typeof($f), ::Val{1}) = Val{false}()
    end
end

for linearity_mask = 0:2^3-1
    lin = Val{map(x->x!=0, (linearity_mask & 4,
                            linearity_mask & 2,
                            linearity_mask & 1))}()

    for f in diadic_of_linearity(lin)
        @eval begin
            haslinearity(::typeof($f), ::Val{2}) = true
            linearity(::typeof($f), ::Val{2}) = $lin
        end
    end
end

# TermCombination datastructure

struct TermCombination
    terms::Set{Dict{Int, Int}} # idx => pow
end

@eval Base.one(::Type{TermCombination}) = $(TermCombination(Set([Dict{Int,Int}()])))
@eval Base.zero(::Type{TermCombination}) = $(TermCombination(Set{Dict{Int,Int}}()))

#=
function Base.:(==)(comb1::TermCombination, comb2::TermCombination)
    comb1.terms == comb2.terms && return true

    n1 = reduce(max, (k for (k,_) in Iterators.flatten(comb1.terms)), init=0)
    n2 = reduce(max, (k for (k,_) in Iterators.flatten(comb2.terms)), init=0)
    n = max(n1, n2)

    _sparse(comb1, n) == _sparse(comb2, n)
end
=#

function Base.:+(comb1::TermCombination, comb2::TermCombination)
    if isone(comb1) && !iszero(comb2)
        return comb2
    elseif isone(comb2) && !iszero(comb1)
        return comb1
    elseif comb1 === comb2
        return comb1
    end
    TermCombination(union(comb1.terms, comb2.terms))
end

Base.:+(comb1::TermCombination) = comb1

function _merge(dict1, dict2)
    d = copy(dict1)
    for (k, v) in dict2
        d[k] = min(2, get(dict1, k, 0) + v)
    end
    d
end

function Base.:*(comb1::TermCombination, comb2::TermCombination)
    if isone(comb1)
        return comb2
    elseif isone(comb2)
        return comb1
    elseif comb1 === comb2 # squaring optimization
        terms = comb1.terms
        # turns out it's enough to track
        # a^2*b^2
        # and a^2 + b^2 + ab
        # have the same hessian sparsity
        t = Dict(k=>2 for (k,_) in
                 Iterators.flatten(terms))
        TermCombination(Set([t]))
        #=
        # square each term
        t1 = [Dict(k=>2 for (k,_) in dict)
              for dict in comb1.terms]
        # multiply each term
        t2 = Dict{Int,Int}[]
        for i in 1:length(terms)
            for j in i+1:length(terms)
                push!(t2, _merge(terms[i], terms[j]))
            end
        end
        TermCombination(union(t1, t2))
        =#
    else
        Set([_merge(dict1, dict2)
             for dict1 in comb1.terms,
             dict2 in comb2.terms]) |> TermCombination
    end
end
Base.:*(comb1::TermCombination) = comb1
Base.iszero(c::TermCombination) = isempty(c.terms)
Base.isone(c::TermCombination) = all(isempty, c.terms)

function _sparse(t::TermCombination, n)
    I = Int[]
    J = Int[]
    for dict in t.terms
        kv = collect(pairs(dict))
        for i in 1:length(kv)
            k, v = kv[i]
            if v>=2
                push!(I, k)
                push!(J, k)
            end
            for j in i+1:length(kv)
                if v >= 1 && kv[j][2] >= 1
                    push!(I, k)
                    push!(J, kv[j][1])
                end
            end
        end
    end
    s1 = sparse(I,J,fill!(BitVector(undef, length(I)), true),n,n)
    s1 .| s1'
end

# 1-arg functions
combine_terms(::Val{true}, term) = term
combine_terms(::Val{false}, term) = term * term

# 2-arg functions
function combine_terms(::Val{linearity}, term1, term2) where linearity

    linear11, linear22, linear12 = linearity
    term = zero(TermCombination)
    if linear11
        if !linear12
            term += term1
        end
    else
        term += term1 * term1
    end

    if linear22
        if !linear12
            term += term2
        end
    else
        term += term2 * term2
    end

    if linear12
        term += term1 + term2
    else
        term += term1 * term2
    end
    term
end


text:
const DefaultRefType = UInt32

## Pools

# Type params:
# * `T` type of categorized values
# * `R` integer type for referencing category levels
# * `V` categorical value type
mutable struct CategoricalPool{T, R <: Integer, V}
    index::Vector{T}        # category levels ordered by their reference codes
    invindex::Dict{T, R}    # map from category levels to their reference codes
    order::Vector{R}        # 1-to-1 map from `index` to `level` (position of i-th category in `levels`)
    levels::Vector{T}       # category levels ordered by externally specified order
    valindex::Vector{V}     # "category value" objects 1-to-1 matching `index`
    ordered::Bool

    function CategoricalPool{T, R, V}(index::Vector{T},
                                      invindex::Dict{T, R},
                                      order::Vector{R},
                                      levels::Vector{T},
                                      valindex::Vector{V},
                                      ordered::Bool) where {T, R, V}
        if iscatvalue(T)
            throw(ArgumentError("Level type $T cannot be a categorical value type"))
        end
        if !iscatvalue(V)
            throw(ArgumentError("Type $V is not a categorical value type"))
        end
        if leveltype(V) !== T
            throw(ArgumentError("Level type of the categorical value ($(leveltype(V))) and of the pool ($T) do not match"))
        end
        if reftype(V) !== R
            throw(ArgumentError("Reference type of the categorical value ($(reftype(V))) and of the pool ($R) do not match"))
        end
        new(index, invindex, order, levels, valindex, ordered)
    end
end

struct LevelsException{T, R} <: Exception
    levels::Vector{T}
end

struct OrderedLevelsException{T, S} <: Exception
    newlevel::S
    levels::Vector{T}
end

## Values

"""
Default categorical value type for
referencing values of type `T`.
"""
struct CategoricalValue{T, R <: Integer}
    level::R
    pool::CategoricalPool{T, R, CategoricalValue{T, R}}
end

"""
`String` categorical value.
Provides `AbstractString` interoperability.
"""
struct CategoricalString{R <: Integer} <: AbstractString
    level::R
    pool::CategoricalPool{String, R, CategoricalString{R}}
end

# union of all categorical value types
const CatValue{R} = Union{CategoricalValue{T, R} where T,
                          CategoricalString{R}}

## Arrays

# Type params:
# * `T` original type of elements before categorization, could be Union{T, Missing}
# * `N` array dimension
# * `R` integer type for referencing category levels
# * `V` original type of elements (excluding Missing) before categorization
# * `C` categorical value type
# * `U` type of missing value, `Union{}` if missing values are not accepted
abstract type AbstractCategoricalArray{T, N, R, V, C, U} <: AbstractArray{Union{C, U}, N} end
const AbstractCategoricalVector{T, R, V, C, U} = AbstractCategoricalArray{T, 1, R, V, C, U}
const AbstractCategoricalMatrix{T, R, V, C, U} = AbstractCategoricalArray{T, 2, R, V, C, U}

struct CategoricalArray{T, N, R <: Integer, V, C, U} <: AbstractCategoricalArray{T, N, R, V, C, U}
    refs::Array{R, N}
    pool::CategoricalPool{V, R, C}

    function CategoricalArray{T, N}(refs::Array{R, N},
                                    pool::CategoricalPool{V, R, C}) where
                                                 {T, N, R <: Integer, V, C}
        T === V || T == Union{V, Missing} || throw(ArgumentError("T ($T) must be equal to $V or Union{$V, Missing}"))
        U = T >: Missing ? Missing : Union{}
        new{T, N, R, V, C, U}(refs, pool)
    end
end
const CategoricalVector{T, R, V, C, U} = CategoricalArray{T, 1, V, C, U}
const CategoricalMatrix{T, R, V, C, U} = CategoricalArray{T, 2, V, C, U}


text:
module scordelis_lo_dsg3_verification

using Test
using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!

function _execute(n = 8, visualize = true)
    # analytical solution for the vertical deflection and the midpoint of the
    # free edge
    analyt_sol=-0.3024;
    # Parameters:
    E=4.32e8;
    nu=0.0;
    thickness = 0.25; # geometrical dimensions are in feet
    R = 25.0;
    L = 50.0;
    formul = FEMMShellT3FFModule

    tolerance = R/n/1000
    fens, fes = T3block(40/360*2*pi,L/2,n,n);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i, 1]; y=fens.xyz[i, 2];
        fens.xyz[i, :] .= (R*sin(a), y, R*(cos(a)-1))
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    femm.mult_el_size = 0.2
    femm.drilling_stiffness_scale = 1.0
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # rigid diaphragm
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [1,3,5]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y
    l1 = selectnode(fens; box = Float64[-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,5,6]
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Midpoint of the free edge
    nl = selectnode(fens; box = Float64[sin(40/360*2*pi)*25 sin(40/360*2*pi)*25 L/2 L/2 -Inf Inf], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt[0, 0, -90, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    resultpercent = dchi.values[nl, 3][1]/analyt_sol*100

    return resultpercent
end

function test_convergence()
    formul = FEMMShellT3FFModule
    results = [
    66.54771615057949,
    85.54615143134853,
    89.85075281481419,
    92.50616661644985,
    95.40469210310079,
    97.65376880486126
    ]
    for (n, res) in  zip([4, 8, 10, 12, 16, 24], results)
        v = _execute(n, false)
        @test isapprox(res, v, rtol = 1.0e-4)
        # @show v
    end
    return true
end

test_convergence()

end # module

"""
The geometry consists of a “hook” in the form of a curved strip rigidly clamped
at one end and loaded with a unit in-plane shear along the width at the other
end. It has two circular segments that are connected at the tangent point. The
smaller segment has a mean radius of 0.3556 m (14 inches) and spans 60° from
the clamped end to the tangent point. The larger segment spans 150° from the
tangent point to the free end and has a mean radius of 1.1684 m (46 inches).
The hook is 0.0508 m (2 inches) thick and 0.508 m (20 inches) wide, modeled as
linear elastic with an elastic modulus of 22.77 MPa (3300 psi) and a Poisson's
ratio of 0.35. In most tests the shear force is applied through the use of a
distributing coupling constraint. The coupling constraint provides coupling
between a reference node on which the load is prescribed and the nodes located
on the free end. The distributed nodal loads on the free end are equivalent to
a uniformly distributed load of 8.7563 N/m (0.05 lb/in). In two of the tests an
equivalent shear force is applied as a distributed shear traction instead.

"""
module raasch_dsg3_verification

using Test
using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!

function _execute(input = "raasch_s4_1x9.inp", visualize = true)
    E = 3300.0;
    nu = 0.35;
    thickness  =  2.0;
    tolerance = thickness/2
    # analytical solution for the vertical deflection under the load
    analyt_sol = 5.02;
    R = 46.0;
    formul = FEMMShellT3FFModule

    output = import_ABAQUS(joinpath(input))
    fens = output["fens"]
    fes = output["fesets"][1]

    connected = findunconnnodes(fens, fes);
    fens, new_numbering = compactnodes(fens, connected);
    fes = renumberconn!(fes, new_numbering);

    fens, fes = Q4toT3(fens, fes)

    # plots = cat(plot_space_box([[0 0 -R/2]; [R/2 R/2 R/2]]),
    #     plot_nodes(fens),
    #     plot_midsurface(fens, fes);
    # dims = 1)
    # pl = render(plots)

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

        # Report
    # @info "Raasch hook, formulation=$(formul)"
    # @info "Mesh: $input"

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    femm.mult_el_size = 0.2
    femm.drilling_stiffness_scale = 1.0
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Clamped end
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end

    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    bfes = meshboundary(fes)
    l1 = selectelem(fens, bfes, box = [97.9615 97.9615 -16 -16 0 20], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, 0, 0.05, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # @infiltrate
    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    nl = selectnode(fens; box = Float64[97.9615 97.9615 -16 -16 0 0], inflate = tolerance)
    targetu =  dchi.values[nl, 3][1]
    # @info "Target: $(round(targetu, digits=8)),  $(round(targetu/analyt_sol, digits = 4)*100)%"

    return targetu/analyt_sol*100
end

function test_convergence()
    formul = FEMMShellT3FFModule
    results = [
    91.7059961843231,
    95.9355786538892,
    97.19276899988246,
    98.38896641657374
    ]
    for (m, res) in zip(["1x9", "3x18", "5x36", "10x72"], results)
        v = _execute("raasch_s4_" * m * ".inp", false)
        @test isapprox(res, v, rtol = 1.0e-4)
        # @show v
    end
    return true
end

test_convergence()

end # module

"""
The initially twisted cantilever beam is one of the standard test
problems for verifying the finite-element accuracy [1]. The beam is
clamped at one end and loaded either with unit in-plane or
unit out-of-plane force at the other. The centroidal axis of the beam is
straight at the undeformed  configuration, while its cross-sections are
twisted about the centroidal axis from 0 at the clamped end to pi/2 at
the free end.

Reference:
Zupan D, Saje M (2004) On "A proposed standard set of problems to test
finite element accuracy": the twisted beam. Finite Elements in Analysis
and Design 40: 1445-1451.

Thicker cross section (t = 0.32)
#     Loading in the Z direction
dir = 3; uex = 0.005424534868469; # Harder: 5.424e-3;
#     Loading in the Y direction
dir = 2; uex = 0.001753248285256; # Harder: 1.754e-3;

Thinner cross section (t = 0.0032)
#     Loading in the Z direction
dir = 3; uex = 0.005256;
#     Loading in the Y direction
dir = 2; uex = 0.001294;

MacNeal,  R. H., and R. L. Harder, “A Proposed Standard Set of Problems to Test Finite Element Accuracy,” Finite Elements in Analysis Design, vol. 11, pp. 3–20, 1985.

Simo,  J. C., D. D. Fox, and M. S. Rifai, “On a Stress Resultant Geometrically Exact Shell Model. Part II: The Linear Theory; Computational Aspects,” Computational Methods in Applied Mechanical Engineering, vol. 73, pp. 53–92, 1989.

"""
module twisted_beam_dsg3_verification

using Test
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!

params_thicker_dir_3 = (t =  0.32, force = 1.0, dir = 3, uex = 0.005424534868469);
params_thicker_dir_2 = (t =  0.32, force = 1.0, dir = 2, uex = 0.001753248285256);

params_thinner_dir_3 = (t =  0.0032, force = 1.0e-6, dir = 3, uex = 0.005256);
params_thinner_dir_2 = (t =  0.0032, force = 1.0e-6, dir = 2, uex = 0.001294);

function _execute(t = 0.32, force = 1.0, dir = 3, uex = 0.005424534868469, nL = 8, nW = 2, visualize = true)
    E = 0.29e8;
    nu = 0.22;
    W = 1.1;
    L = 12.0;
    formul = FEMMShellT3FFModule

    tolerance = W/nW/100
    fens, fes = T3block(L,W,nL,nW,:a);
    fens.xyz = xyz3(fens)
    for i in 1:count(fens)
        a=fens.xyz[i,1]/L*(pi/2); y=fens.xyz[i,2]-(W/2); z=fens.xyz[i,3];
        fens.xyz[i,:]=[fens.xyz[i,1],y*cos(a)-z*sin(a),y*sin(a)+z*cos(a)];
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), t), mater)
    # femm.mult_el_size = 0.2
    # femm.drilling_stiffness_scale = 1.0
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Clamped end
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in 1:6
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    nl = selectnode(fens; box = Float64[L L 0 0 0 0], tolerance = tolerance)
    loadbdry = FESetP1(reshape(nl, 1, 1))
    lfemm = FEMMBase(IntegDomain(loadbdry, PointRule()))
    v = FFlt[0, 0, 0, 0, 0, 0]
    v[dir] = force
    fi = ForceIntensity(v);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    result =  dchi.values[nl, dir][1]/uex*100
    return result
end

function test_convergence()
    formul = FEMMShellT3FFModule
    results = [
    39.709921740907355,
    68.87306876326497,
    86.01944734315117,
    95.04101960524827,
    53.10262177376127,
    83.8593790803426,
    94.91359387874728,
    98.21549248655576,
    48.16757753755567,
    79.43420077873479,
    92.54464819755955,
    96.85008269135751,
    50.577029703967334,
    80.34160167730624,
    92.48675665271801,
    96.7096641005938,
    ]
    for n in [2, 4, 8, 16, ]
        v = _execute(params_thicker_dir_2..., 2*n, n, false)
        # println(v, ",")
        @test isapprox(v, popat!(results, 1), rtol = 1.0e-3)
    end
    for n in [2, 4, 8, 16, ]
        v = _execute(params_thicker_dir_3..., 2*n, n, false)
        # println(v, ",")
        @test isapprox(v, popat!(results, 1), rtol = 1.0e-3)
    end
    for n in [2, 4, 8, 16, ]
        v = _execute(params_thinner_dir_2..., 2*n, n, false)
        # println(v, ",")
        @test isapprox(v, popat!(results, 1), rtol = 1.0e-3)
    end
    for n in [2, 4, 8, 16, ]
        v = _execute(params_thinner_dir_3..., 2*n, n, false)
        # println(v, ",")
        @test isapprox(v, popat!(results, 1), rtol = 1.0e-3)
    end
    return true
end

test_convergence()

end # module

"""
MODEL DESCRIPTION

Z-section cantilever under torsional loading.

Linear elastic analysis, Young's modulus = 210 GPa, Poisson's ratio = 0.3.

All displacements are fixed at X=0.

Torque of 1.2 MN-m applied at X=10. The torque is applied by two
uniformly distributed shear loads of 0.6 MN at each flange surface.

Objective of the analysis is to compute the axial stress at X = 2.5 from fixed end.

NAFEMS REFERENCE SOLUTION

Axial stress at X = 2.5 from fixed end (point A) at the midsurface is -108 MPa.
"""
module LE5_Z_cantilever_dsg3_verification

using Test
using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!

function _execute(input = "nle5xf3c.inp", nrefs = 0, visualize = true)
    E = 210e9;
    nu = 0.3;
    L = 10.0;
    thickness = 0.1
    formul = FEMMShellT3FFModule

    tolerance = thickness/1000
    output = import_ABAQUS(joinpath(dirname(@__FILE__()), input))
    fens = output["fens"]
    fes = output["fesets"][1]

    connected = findunconnnodes(fens, fes);
    fens, new_numbering = compactnodes(fens, connected);
    fes = renumberconn!(fes, new_numbering);

    for r in 1:nrefs
        fens, fes = T3refine(fens, fes)
    end

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    # Report
        # @info "Mesh: $input, nrefs = $nrefs"

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    femm.mult_el_size = 0.2
    femm.drilling_stiffness_scale = 1.0
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # plots = cat(plot_space_box([[0 0 -L/2]; [L/2 L/2 L/2]]),
    #     plot_nodes(fens),
    #     plot_midsurface(fens, fes);
    # dims = 1)
    # pl = render(plots)

    # Apply EBC's
    # plane of symmetry perpendicular to X
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,2,3,]
        setebc!(dchi, l1, true, i)
    end

    applyebc!(dchi)
    numberdofs!(dchi);
    # @show dchi.nfreedofs

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    nl = selectnode(fens; box = Float64[10.0 10.0 1.0 1.0 0 0], tolerance = tolerance)
    loadbdry1 = FESetP1(reshape(nl, 1, 1))
    lfemm1 = FEMMBase(IntegDomain(loadbdry1, PointRule()))
    fi1 = ForceIntensity(FFlt[0, 0, +0.6e6, 0, 0, 0]);
    nl = selectnode(fens; box = Float64[10.0 10.0 -1.0 -1.0 0 0], tolerance = tolerance)
    loadbdry2 = FESetP1(reshape(nl, 1, 1))
    lfemm2 = FEMMBase(IntegDomain(loadbdry2, PointRule()))
    fi2 = ForceIntensity(FFlt[0, 0, -0.6e6, 0, 0, 0]);
    F = distribloads(lfemm1, geom0, dchi, fi1, 3) + distribloads(lfemm2, geom0, dchi, fi2, 3);

    # @infiltrate
    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    return minimum(dchi.values[:, 3]), maximum(dchi.values[:, 3])
end

function test_convergence()
    formul = FEMMShellT3FFModule
    # @info "LE5 Z-cantilever, formulation=$(formul)"
    for n in [0,  ]
        res = _execute("nle5xf3c.inp", n, false)
        @test res[1] ≈ (-0.01568028415401719, 0.015401663810490641)[1]
        @test res[2] ≈ (-0.01568028415401719, 0.015401663810490641)[2]
    end
    return true
end

# res = test_st3dsg("nle5xf3c.inp", 0, false)
# @test res[1] ≈ (-0.01568028415401719, 0.015401663810490641)[1]
# @test res[2] ≈ (-0.01568028415401719, 0.015401663810490641)[2]

test_convergence()

end # module

"""
The barrel vault (Scordelis-Lo) roof is one of the benchmarks for linear elastic
analysis of shells.

The candidate element's usefulness in irregular geometries (and most practical
cases involve a high degree of geometric irregularity) is tested. As would be
expected,the irregular mesh results are not as good as those provided by a
regular mesh with the same number of variables.

Problem description

The physical basis of the problem is a deeply arched roof supported only
by diaphragms at its curved edges (an aircraft hanger), deforming under its own
weight. It is interesting to observe that the geometry is such that the
centerpoint of the roof moves upward under the self-weight (downwardly directed)
load. Perhaps this is one reason why the problem is not straightforward
numerically.
"""
module barrel_vault_test

using Test
using LinearAlgebra
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

function _execute(input = "barrelvault_s3r_fineirreg.inp", visualize = true)
    E = 3.0e6;
    nu = 0.0;
    thickness = 3.0;
    tolerance = thickness/20
    # analytical solution for the vertical deflection under the load
    analyt_sol = -0.3024*12;
    R = 25.0*12;
    L = 50.0*12;
    formul = FEMMShellT3FFModule

    output = import_ABAQUS(joinpath(dirname(@__FILE__()), input))
    fens = output["fens"]
    fes = output["fesets"][1]

     # boundingbox(fens.xyz)
    connected = findunconnnodes(fens, fes);
    fens, new_numbering = compactnodes(fens, connected);
    fes = renumberconn!(fes, new_numbering);

    fens, fes = mergenodes(fens, fes, thickness/10)

    # plots = cat(plot_space_box([[0 0 -R/2]; [R/2 R/2 R/2]]),
    #     plot_nodes(fens),
    #     plot_midsurface(fens, fes);
    # dims = 1)
    # pl = render(plots)

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

    # @info "Mesh: $input"

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    femm.mult_el_size = 0.2
    femm.drilling_stiffness_scale = 0.1
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # vtkwrite("barrel_vault_test-mesh.vtu", fens, fes)
    # bfes = meshboundary(fes)
    # vtkwrite("barrel_vault_test-boundary-mesh.vtu", fens, bfes)

    # Apply EBC's
    # Rigid diaphragm end
    l1 = selectnode(fens; box = Float64[-Inf Inf -Inf Inf 0 0], inflate = tolerance)
    for i in [1,2]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Z
    l1 = selectnode(fens; box = Float64[-Inf Inf -Inf Inf L/2 L/2], inflate = tolerance)
    for i in [3,4,5]
        setebc!(dchi, l1, true, i)
    end
    # plane of symmetry perpendicular to Y (apex)
    l1 = selectnode(fens; box = Float64[-Inf Inf 0 0 -Inf Inf], inflate = tolerance)
    for i in [2,4,6]
        setebc!(dchi, l1, true, i)
    end

    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    # Midpoint of the free edge
    nl = selectnode(fens; box = Float64[-Inf Inf sin(40/360*2*pi)*R sin(40/360*2*pi)*R L/2 L/2], inflate = tolerance)
    lfemm = FEMMBase(IntegDomain(fes, TriRule(3)))
    fi = ForceIntensity(FFlt[-0.625, 0, 0, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);

    # @infiltrate
    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])

    targetu =  dchi.values[nl, 1][1]
    # @info "Solution: $(round(targetu, digits=8)),  $(round(targetu/analyt_sol, digits = 4)*100)%"

    # Generate a graphical display of resultants
    cylindrical!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
        r = -vec(XYZ); r[3] = 0.0
        csmatout[:, 3] .= vec(r)/norm(vec(r))
        csmatout[:, 2] .= (0.0, 0.0, 1.0)
        cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
        return csmatout
    end
    ocsys = CSys(3, 3, cylindrical!)
    test_results = [
    (-1520.6449167366522, 14.067403309095397)
    (-73.8262145426215,   425.93651541819503)
    (-0.005341121492547284, 953.0929383629322)
    ]
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("m$nc", fld.values))
        @test isapprox(minimum(fld.values), test_results[nc][1], rtol = 0.01)
        @test isapprox(maximum(fld.values), test_results[nc][2], rtol = 0.01)
    end
    # vtkwrite("$(input)-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    test_results = [
    (-306.83173146926623,  309.5860993742647),
    (-1011.4832977998705,  2167.0403478574167),
    (-687.3500043290137,   69.38703021678862),
    ]
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("n$nc", fld.values))
        @test isapprox(minimum(fld.values), test_results[nc][1], rtol = 0.01)
        @test isapprox(maximum(fld.values), test_results[nc][2], rtol = 0.01)
    end
    # vtkwrite("$(input)-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    test_results = [
    (-13.784764125688811, 21.165312065421237),
    (-7.4216963152916255, 25.679801383392967)
    ]
    scalars = []
    for nc in 1:2
        fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("q$nc", fld.values))
        @test isapprox(minimum(fld.values), test_results[nc][1], rtol = 0.01)
        @test isapprox(maximum(fld.values), test_results[nc][2], rtol = 0.01)
    end
    # vtkwrite("$(input)-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

    # Visualization
    if !visualize
        return true
    end
    scattersysvec!(dchi, (R/2)/maximum(abs.(U)).*U)
    update_rotation_field!(Rfield0, dchi)
    plots = cat(plot_space_box([[0 0 -R]; [R R R]]),
        plot_nodes(fens),
        plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
    dims = 1)
    pl = render(plots)
    return true
end

function test_convergence()

    # @info "Scordelis-Lo Abaqus model"

    _execute("barrelvault_stri3_irreg.inp", false)

    # _execute("barrelvault_s3r_fineirreg.inp", false)

    return true
end

test_convergence()

end # module


text:
mutable struct Duration
  start_day::Date
  end_day::Date
end

mutable struct Params
  year::Int64
  month::Int64
  day::Int64
  stime::String # Dates.Time
  etime::String # Dates.Time

  multiday::Bool
  duration::Duration
end

mutable struct Header
  samples_per_record::Int64
  seconds_per_record::Int64
  records_per_packet::Int64
  trial_size::Int64
  eggs_reporting::Int64
  start_time::DateTime
  end_time::DateTime
  seconds_of_data::Int64

  Header() = new()
end

mutable struct Result
  header::Header
  data::DataFrame
  raw::String
end


text:
# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Unit prices
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# compute UNIT_PRICE
transform!(df, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] => ByRow((v,q) -> v/q) => :UNIT_PRICE)

# compute UNIT_PRICE_CHANGE
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "PRODUCT_NC", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, :UNIT_PRICE => pct_change => :UNIT_PRICE_CHANGE)

# compute MOM_PRICE_CHANGE
cols_grouping = ["DECLARANT_ISO", "PARTNER_ISO", "PRODUCT_NC", "FLOW"]
gdf = groupby(df, cols_grouping)
df = transform(gdf, [:PERIOD, :UNIT_PRICE] => mom_change => :MOM)

# export data
df_export = subset(df, :PERIOD => ByRow( x-> string(x) == year*sub_month))
CSV.write(dir_io * "clean/" * "df_prices_" * year * sub_month * ".csv", df_export)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------
# Share of Belgium
# -------------------------------------------------------------------------------------------------------------------------------------------------------------

# compute EU_WORLD (EU imports/exports to the entire world)
cols_grouping = ["PRODUCT_NC", "FLOW", "PERIOD"]
prod_class_other = ["PRODUCT_SITC", "PRODUCT_CPA2002", "PRODUCT_CPA2008", "PRODUCT_CPA2_1", "PRODUCT_BEC", "PRODUCT_BEC5"]
cols_grouping = [cols_grouping; prod_class_other]

gdf = groupby(df, cols_grouping)
df_sum_EU_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)
df_sum_EU_WORLD.TRADE_TYPE .= "total"

# compute EU_EXTRA and EU_INTRA (EU imports/exports to extra/intra EU countries)
cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
cols_grouping = [cols_grouping; prod_class_other]
gdf = groupby(df, cols_grouping)
df_sum_EU_EXTRA_INTRA = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

# compute BE_WORLD (BE imports/exports to the entire world)
df_BE = subset(df, :DECLARANT_ISO => ByRow(x -> x == "BE"))
cols_grouping = ["PRODUCT_NC", "FLOW", "PERIOD"]
cols_grouping = [cols_grouping; prod_class_other]

gdf = groupby(df_BE, cols_grouping)
df_sum_BE_WORLD = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)
df_sum_BE_WORLD.TRADE_TYPE .= "total"

# compute BE_EXTRA and BE_INTRA (BE imports/exports to extra/intra EU countries)
cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
cols_grouping = [cols_grouping; prod_class_other]
gdf = groupby(df_BE, cols_grouping)
df_sum_BE_EXTRA_INTRA = combine(gdf, [:VALUE_IN_EUROS, :QUANTITY_IN_KG] .=> (x -> sum(skipmissing(x))), renamecols=false)

# merge tables
df_sum_EU = vcat(df_sum_EU_EXTRA_INTRA, df_sum_EU_WORLD)
df_sum_BE = vcat(df_sum_BE_EXTRA_INTRA, df_sum_BE_WORLD)

# join, rename and compute shares in percentages (?)
cols_join = ["TRADE_TYPE", "PRODUCT_NC", "FLOW", "PERIOD"]
cols_join = [cols_join; prod_class_other]
df_join = leftjoin(df_sum_EU, df_sum_BE, on=cols_join, makeunique=true)
rename!(df_join, ["VALUE_IN_EUROS", "QUANTITY_IN_KG", "VALUE_IN_EUROS_1", "QUANTITY_IN_KG_1"] .=> ["VALUE_EU", "QUANITY_EU", "VALUE_BE", "QUANITY_BE"])
df_join.VALUE_SHARE = df_join.VALUE_BE ./ df_join.VALUE_EU .* 100
df_join.QUANTITY_SHARE = df_join.QUANITY_BE ./ df_join.QUANITY_EU .* 100

# possibly need to sort before grouping, seems that grouping sorts incorrectly sometimes?
sort!(df_join, [:PRODUCT_NC, :FLOW, :TRADE_TYPE, :PERIOD]) # dont understand why I need the sorting here? otherwise some MOM are not taken?

# compute SHARE_CHANGE
cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW"]
gdf = groupby(df_join, cols_grouping)
df_share_BE = transform(gdf, [:VALUE_SHARE, :QUANTITY_SHARE] .=> pct_change .=> [:VALUE_SHARE_CHANGE, :QUANTITY_SHARE_CHANGE])

# compute MOM_CHANGE
cols_grouping = ["TRADE_TYPE", "PRODUCT_NC", "FLOW"]
gdf = groupby(df_share_BE, cols_grouping)
df_share_BE = transform(gdf, [:PERIOD, :VALUE_SHARE] => mom_change => :MOM_VALUE,
                             [:PERIOD, :QUANTITY_SHARE] => mom_change => :MOM_QUANTITY)

# export data
df_export = subset(df_share_BE, :PERIOD => ByRow( x-> string(x) == year*sub_month))
CSV.write(dir_io * "clean/" * "df_share_BE_" * year * sub_month * ".csv", df_export)

# -------------------------------------------------------------------------------------------------------------------------------------------------------------


text:
using SpecialFunctions; # need to evaluate the Hankel function (Green's function)

function matricesAlpert(k,prt,order)

    (nskip, chi, weights) = Alpert_params(order)
    # nskip = a
    ninterp = order + 3;

    τ = prt.τ
    ν = prt.ν
    s = prt.s
    ds= prt.ds

    N = prt.N
    Δs = prt.Δs

    x1 = prt.x[:,1];     x2 = prt.x[:,2]
    dx1 = prt.dx[:,1];   dx2 = prt.dx[:,2]

    S  = zeros(N,N)+im*zeros(N,N)
    K  = zeros(N,N)+im*zeros(N,N)
    Kp = zeros(N,N)+im*zeros(N,N)
    T  = zeros(N,N)+im*zeros(N,N)

    for i=1:N

        for j=1:N

            r = sqrt((x1[i]-x1[j])^2+(x2[i]-x2[j])^2);

            H0 = hankelh1(0,k*r); H1 = hankelh1(1,k*r)

            # if i!=j

                S2  = im/4*H0*τ[j]*Δs*ds[j]
                K2  = im*k/4*H1/r*((x1[i]-x1[j])*dx2[j]-(x2[i]-x2[j])*dx1[j])*Δs*ds[j]
                Kp2 = im*k/4*H1/r*((x1[j]-x1[i])*dx2[i]-(x2[j]-x2[i])*dx1[i])*(τ[j]/τ[i])*Δs*ds[j]

                T2 = ((x1[i]-x1[j])*dx1[j]+(x2[i]-x2[j])*dx2[j])*((x1[i]-x1[j])*dx1[i]+(x2[i]-x2[j])*dx2[i])/r^2
                T2 = im/2*T2*(k^2*H0-2*k*H1/r)+im*k/2*(dx1[i]*dx1[j]+dx2[i]*dx2[j])*H1/r
                T2 = -0.5*T2*Δs*ds[j]/τ[j]

            # end

            S[i,j]  = S2
            K[i,j]  = K2
            Kp[i,j] = Kp2
            T[i,j]  = T2

        end

        # @. S[i, mod(i+(-nskip+1:nskip-1)-1, N)+1] = 0

        S[i, mod.(i.+(-nskip+1:nskip-1).-1, N).+1] .= 0
        K[i, mod.(i.+(-nskip+1:nskip-1).-1, N).+1] .= 0
        Kp[i,mod.(i.+(-nskip+1:nskip-1).-1, N).+1] .= 0
        T[i, mod.(i.+(-nskip+1:nskip-1).-1, N).+1] .= 0

    end

    for p = 1:size(chi, 1)

        if prt.PD == 0
            t = Δs*collect(0:N-1)
            w = mod.(t .+ Δs*chi[p], 2*pi)
            dw = ones(size(w))
        else
            t = Δs*collect(1:N)
            (w,dw) = polySigmoid(mod.(t .+ Δs*chi[p], 2*pi),prt.PD)
        end

        (aux_x, aux_dx) = prt.param(w)
        # @show aux_x
        # return

        aux_x1  = aux_x[:,1]
        aux_x2  = aux_x[:,2]

        aux_dx1 = aux_dx[:,1]
        aux_dx2 = aux_dx[:,2]

        tau_aux = sqrt.(aux_dx1.^2+aux_dx2.^2)

        aux_n1 = aux_dx2./tau_aux
        aux_n2 =-aux_dx1./tau_aux

        aux_r = sqrt.((x1-aux_x1).^2+(x2-aux_x2).^2)

        aux_H0 = hankelh1.(0,k*aux_r)
        aux_H1 = hankelh1.(1,k*aux_r)

        ind = floor.(chi[p] - ninterp/2 + 1) .+ (0 : ninterp - 1)

        L   = lagrangeInterp(ind', chi[p])

        kerS = im/4*aux_H0.*tau_aux.*dw
        kerK = im*k/4*aux_H1./aux_r.*((x1-aux_x1).*aux_dx2-(x2-aux_x2).*aux_dx1).*dw
        kerKp= im*k/4*aux_H1./aux_r.*((aux_x1-x1).*ν[:,1]+(aux_x2-x2).*ν[:,2]).*(tau_aux.*dw)

        kerT = ((x1-aux_x1).*aux_dx1+(x2-aux_x2).*aux_dx2).*((x1-aux_x1).*dx1+(x2-aux_x2).*dx2)./aux_r.^2
        kerT = im/2*kerT.*(k^2*aux_H0-2*k*aux_H1./aux_r)+
                im*k/2*(dx1.*aux_dx1+dx2.*aux_dx2)./aux_r.*aux_H1
        kerT = -0.5*kerT.*dw./tau_aux

        for i0=1:N
           j0  = i0 .+ ind
           j0  = floor.(Int,mod.(j0.-1, N) .+ 1)
           S[i0,j0] .= S[i0,j0]  + Δs*weights[p]*kerS[i0]*L[:,1]
           K[i0,j0] .= K[i0,j0]  + Δs*weights[p]*kerK[i0]*L[:,1]
           Kp[i0,j0].= Kp[i0,j0] + Δs*weights[p]*kerKp[i0]*L[:,1]
           T[i0,j0] .= T[i0,j0]  + Δs*weights[p]*kerT[i0]*L[:,1]
        end

    end

    T = k^2*S.*(ν*ν') + T

    mat = (S=S,K=K,Kp=Kp,T=T);

end

########################################################################

function  Alpert_params(order)

    if order == 2
        a = 1
        chi_p = [.591549430918953e-01]
        w_p = [5.000000000000000e-01]

    elseif order == 6
        a = 3
        chi_p= [4.004884194926570e-03;
                7.745655373336686e-02;
                3.972849993523248e-01;
                1.075673352915104e+00;
                2.003796927111872e+00]

        w_p = [ 1.671879691147102e-02;
                1.636958371447360e-01;
                4.981856569770637e-01;
                8.372266245578912e+00;
                9.841730844088381e+00]

    elseif order == 10
        a = 6
        chi_p = [1.175089381227308e-03;
             1.877034129831289e-02;
             9.686468391426860e-02;
             3.004818668002884e-01;
             6.901331557173356e-01;
             1.293695738083659e+00;
             2.090187729798780e+00;
             3.016719313149212e+00;
             4.001369747872486e+00;
             5.000025661793423e+00]

        w_p = [4.560746882084207e-03;
           3.810606322384757e-02;
           1.293864997289512e-01;
           2.884360381408835e-01;
           4.958111914344961e-01;
           7.077154600594529e-01;
           8.741924365285083e-01;
           9.661361986515218e-01;
           9.957887866078700e-01;
           9.998665787423845e-01]

    elseif order == 16
        a = 10;
        chi_p = [8.371529832014113e-04
                  1.239382725542637e-02
                  6.009290785739468e-02
                  1.805991249601928e-01
                  4.142832599028031e-01
                  7.964747731112430e-01
                  1.348993882467059e+00
                  2.073471660264395e+00
                  2.947904939031494e+00
                  3.928129252248612e+00
                  4.957203086563112e+00
                  5.986360113977494e+00
                  6.997957704791519e+00
                  7.999888757524622e+00
                  8.999998754306120e+00]

        w_p = [3.190919086626234e-03
               2.423621380426338e-02
               7.740135521653088e-02
               1.704889420286369e-01
               3.029123478511309e-01
               4.652220834914617e-01
               6.401489637096768e-01
               8.051212946181061e-01
               9.362411945698647e-01
               1.014359775369075e+00
               1.035167721053657e+00
               1.020308624984610e+00
               1.004798397441514e+00
               1.000395017352309e+00
               1.000007149422537e+00]

    end

    chi_p = [chi_p' -chi_p[end:-1:1]']';
    w_p = [w_p' w_p[end:-1:1]']';

    return (a, chi_p, w_p)

end

##############################################################################

function lagrangeInterp(xx, t)
    #  L is a row vector calculating the coefficients of lagrange interpolation
    #  at t xx are data
    n = length(xx) # here xx = index = [-2, -1, 0, 1, 2]
    denom = xx'*ones(1, n) - ones(n, 1)*xx
    numer = t*ones(n, n) - ones(n, 1)*xx
    for i = 1:n
      denom[i, i] = 1.0
      numer[i, i] = 1.0
    end
    temp = numer./denom
    L = prod(temp,dims=2)
    return L
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp