text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Test

using Statistics
using ArraysOfArrays, StatsBase

@testset "onlinestats" begin
    n = 10
    data1  = [(-1)^x*exp(x-n/2) for x in 1:n]
    data2 = reverse(data1, dims=1)

    w = [exp(x-n/2) for x in 1:n]

    data = VectorOfSimilarVectors(vcat(data1', data2'))
    m = innersize(data, 1)

    @testset "BAT.kbn_add" begin
        bg = +1E5
        sm = +1e-4
        res = bg + sm

        bgT = (bg, 0.0)
        smT = (sm, 0.0)

        @test collect(BAT.kbn_add((0.0,0.0),0.0)) ≈ [0.0,0.0]
        @test sum(BAT.kbn_add((0.0,0.0),eps())) ≈ eps()
        @test sum(BAT.kbn_add((bg,0.0),sm)) ≈ res
        @test sum(BAT.kbn_add((sm,0.0),bg)) ≈ res

        @test sum(BAT.kbn_add(bgT,smT)) ≈ res

        @test typeof(@inferred BAT.kbn_add(bgT,sm)) == NTuple{2, Float64}
        @test typeof(@inferred BAT.kbn_add(bgT,smT)) == NTuple{2, Float64}
        @test typeof(@inferred BAT.kbn_add((0.0f0, 0.0f0), 0.0f0)) == NTuple{2, Float32}
        @test typeof(@inferred BAT.kbn_add((0.0f0, 0.0f0), (0.0f0, 0.0f0))) == NTuple{2, Float32}
    end

    @testset "BAT.OnlineMvMean" begin
        @test typeof(@inferred BAT.OnlineMvMean(n)) <: AbstractVector{Float64}
        @test typeof(@inferred BAT.OnlineMvMean{Float32}(n)) <: AbstractVector{Float32}

        mvmean = BAT.OnlineMvMean(m)
        @test size(mvmean, 1) == m

        countM = 3
        mvmeans = Array{BAT.OnlineMvMean{Float64}}(undef, countM)

        for i in axes(mvmeans, 1)
            mvmeans[i] = BAT.OnlineMvMean(m)
        end

        for i in axes(data, 1)
            push!(mvmean, data[i], w[i]);
            push!(mvmeans[(i % countM) + 1], data[i], w[i]);
        end

        @test mvmean ≈ mean(data, Weights(w))

        res = merge(mvmeans...)
        @test res ≈ mean(data, Weights(w))

        mvmean = BAT.OnlineMvMean(m)
        res = append!(deepcopy(mvmean), data)
        @test res ≈ mean(data)
        res = append!(deepcopy(mvmean), data, w)
        @test res ≈ mean(data, Weights(w))
    end

    @testset "BAT.OnlineMvCov" begin
        @test typeof(@inferred BAT.OnlineMvCov(n)) <: AbstractMatrix{Float64}
        @test typeof(@inferred BAT.OnlineMvCov{Float32, ProbabilityWeights}(n)) <: AbstractMatrix{Float32}

        for wKind in [ProbabilityWeights, FrequencyWeights, AnalyticWeights, Weights]
            mvcov = BAT.OnlineMvCov{Float64, wKind}(m)
            for i in eachindex(data)
                push!(mvcov, data[i], w[i]);
            end

            @test mvcov ≈ cov(data, wKind(w); corrected = (wKind != Weights))
        end

        countMvCovs = 3
        mvcovs = Array{BAT.OnlineMvCov{Float64, ProbabilityWeights}}(undef, countMvCovs)
        for i in axes(mvcovs,1)
            mvcovs[i] = BAT.OnlineMvCov(m)
        end
        mvcovc = BAT.OnlineMvCov(m)

        for i in eachindex(data)
            push!(mvcovs[(i % countMvCovs) + 1], data[i], w[i]);
            push!(mvcovc, data[i], w[i]);
        end

        res = merge(mvcovs...)
        @test res ≈ cov(data, ProbabilityWeights(w); corrected = true)
        @test res ≈ mvcovc
        @test Float64(res.sum_w) ≈ Float64(mvcovc.sum_w)
        @test Float64(res.sum_w2) ≈ Float64(mvcovc.sum_w2)
        @test res.n ≈ mvcovc.n
        @test res.Mean_X ≈ mvcovc.Mean_X
        @test res.New_Mean_X ≈ mvcovc.New_Mean_X
        @test push!(res, data[1], zero(Float64)) ≈ mvcovc

        mvcov = BAT.OnlineMvCov(m)
        res = append!(deepcopy(mvcov), data)
        @test res ≈ cov(data)
        res = append!(deepcopy(mvcov), data, w)
        @test res ≈ cov(data, ProbabilityWeights(w); corrected = true)
    end

    @testset "BAT.BasicMvStatistic" begin
        bmvstats = BAT.BasicMvStatistics{Float64, ProbabilityWeights}(m)

        countBMS = 3
        bmvs = Array{BAT.BasicMvStatistics{Float64, ProbabilityWeights}}(undef, countBMS)
        for i in axes(bmvs,1)
            bmvs[i] = BAT.BasicMvStatistics{Float64, ProbabilityWeights}(m)
        end

        for i in eachindex(data)
            BAT.push!(bmvstats, data[i], w[i]);
            #BAT.push!(bmvs[(i % countBMS) + 1], data[i], w[i]);
        end
    end
end


text:
include("common.jl")
include("wrapped_rng.jl")
include("PCG/runtests.jl")
include("MersenneTwisters/runtests.jl")
include("Xorshifts/runtests.jl")
include("Random123/runtests.jl")

ext: jl


text:
module jump

    using JuMP
    using Ipopt
    using Distributions
    srand(12345)
    normal = Normal(0,0.01)

    function run()

        # create a model
        m = Model(solver=IpoptSolver())

        # define constants (N, price, etc)
        N = 100
        price = collect(linspace(0.05,0.95,N))
        beta0 = 0.1
        demand0 = (1.0-price) / (2*beta0)
        demand = demand0 - rand(normal,N)

        #define JuMP variables
        @variable(m,eps[1:N])
        @variable(m,0 <= beta <= 1)

        # define constraints and objective
        @objective(m,Min,sum(eps[i]^2 for i in 1:N))
        @constraint(m,constr[i=1:N], 1.0 - 2*beta*(demand[i]-eps[i])-price[i] == 0)

        # solve
        status = solve(m)
        Dict(:obj=>getobjectivevalue(m),:beta=>getvalue(beta),:eps=>getvalue(eps))
    end
end

module nlopt

    using Distributions
    using NLopt
    srand(12345)
    normal = Normal(0,0.01)
    N = 100
    price = collect(linspace(0.05,0.95,N))
    beta0 = 0.1
    demand0 = (1.0-price) / (2*beta0)
    demand = demand0 - rand(normal,N)

    # objective function
    function obj(x,g)
        # x = [e_1,...,e_N,beta] i.e. (N+1,1)

        # gradient = [grad w.r.t e_1,
        #            ...,grad w.r.t e_N,
        #                grad w.r.t beta] i.e. (N+1,1)
        if length(g)>0
           g[:] = vcat(2.0*x[1:(end-1)],0.0)
        end
        # value of objective
        r = sum(x[1:(end-1)].^2)
        return r
    end

    # function test_constr()
    #     x = rand(N+1)
    #     my_c = zeros(N+1)
    #     my_g = zeros(N+1,N)
    #     constr(my_c,x,my_g,N,demand,price);
    #     @test all(my_g .== Calculus.hessian()
    # end

    #
    function constr(r::Vector,x::Vector,g::Matrix,n,q,p)
        if length(g) > 0
            # g has to be n by m for nlopt
            g[:,:] = cat(1,diagm(2*x[end]*ones(n)),-2*(q'.-x[1:(end-1)]'))
        end
        # value of contraints
        r[:] = 1.0 - 2*x[end].*(q.-x[1:(end-1)]) .- p
    end
    constr_clos(r::Vector,x::Vector,g::Matrix) = constr(r::Vector,x::Vector,g::Matrix,N,demand,price)

    function run()
        opt = Opt(:LD_SLSQP,N+1)
        lower_bounds!(opt,[[-Inf for i in 1:N]...,0.0])
        upper_bounds!(opt,[[Inf for i in 1:N]...,1.0])
        min_objective!(opt,obj)
        equality_constraint!(opt,constr_clos,[1e-10 for i in 1:N])
        xtol_rel!(opt,1e-4)
        ftol_rel!(opt,1e-6)

        res = optimize(opt, vcat(rand(normal,N),0.9))
        println("beta0 = $beta0")
        println("beta  = $(res[2][end])")
        r = zeros(N)
        g =zeros(N+1,N)
        constr(r,res[2],g,N,demand,price);
        println("maximal error of constraint at solution = $(maxabs(r))")
        return res
    end

end # module

max_stars_repo_name: snowdj/CoursePack

max_issues_repo_name: snowdj/CoursePack

max_forks_repo_name: snowdj/CoursePack


text:
module SCP_Exact

include("../Misc.jl")
using .Misc

using Distributions, SpecialFunctions, Einsum, Combinatorics
import Base.Iterators: product

export generate, log_marginal

function generate(I::Ƶ, J::Ƶ, K::Ƶ, R::Ƶ; μ::ℜ=3.5, γ::ℜ=0.1) where {ℜ<:Real, Ƶ<:Int}
    a::ℜ, b::ℜ = I*J*K*γ, γ/μ
    α_R, α_KR, α_IKR, α_JKR = fill(a/R,R), fill(a/(K*R),K,R), fill(a/(I*K*R),I,K,R), fill(a/(J*K*R),J,K,R)

    λ = rand(Gamma(a,1.0/b))
    θ_R = rand(Dirichlet(α_R))
    θ_KR = reshape([θ_kr for r=1:R for θ_kr=rand(Dirichlet(α_KR[:,r]))],K,R)
    θ_IKR = reshape([θ_ikr for k=1:K, r=1:R for θ_ikr=rand(Dirichlet(α_IKR[:,k,r]))],I,K,R)
    θ_JKR = reshape([θ_jkr for k=1:K, r=1:R for θ_jkr=rand(Dirichlet(α_JKR[:,k,r]))],J,K,R)

    Λ_IJK::Array{ℜ} = [λ * sum(θ_R .* θ_KR[k,:] .* θ_IKR[i,k,:] .* θ_JKR[j,k,:]) for i=1:I, j=1:J, k=1:K]

    X::Array{ℜ} = map.(λ_ijk -> rand(Poisson(λ_ijk)), Λ_IJK)
    X⁺::Array{ℜ} = sum(X,dims=3)[:,:]
    return X, X⁺, λ, (θ_R, θ_KR, θ_IKR, θ_JKR)
end

function generate(I::Ƶ, J::Ƶ, K::Ƶ, R::Ƶ, T::Ƶ; μ::ℜ=3.5, γ::ℜ=0.1) where {ℜ<:Real, Ƶ<:Int}
    a::ℜ, b::ℜ = I*J*K*γ, γ/μ
    α_R, α_KR, α_IKR, α_JKR = fill(a/R,R), fill(a/(K*R),K,R), fill(a/(I*K*R),I,K,R), fill(a/(J*K*R),J,K,R)

    θ_R = rand(Dirichlet(α_R))
    θ_KR = reshape([θ_kr for r=1:R for θ_kr=rand(Dirichlet(α_KR[:,r]))],K,R)
    θ_IKR = reshape([θ_ikr for k=1:K, r=1:R for θ_ikr=rand(Dirichlet(α_IKR[:,k,r]))],I,K,R)
    θ_JKR = reshape([θ_jkr for k=1:K, r=1:R for θ_jkr=rand(Dirichlet(α_JKR[:,k,r]))],J,K,R)

    θ_IJK::Array{ℜ} = [sum(θ_R .* θ_KR[k,:] .* θ_IKR[i,k,:] .* θ_JKR[j,k,:]) for i=1:I, j=1:J, k=1:K]
    θ_flat = reshape(θ_IJK, I*J*K)

    X::Array{ℜ} = reshape(rand(Multinomial(T,θ_flat)),I,J,K)
    X⁺::Array{ℜ} = sum(X,dims=3)[:,:,1]
    return X, X⁺, (θ_R, θ_KR, θ_IKR, θ_JKR)
end

function comb(X_ij::Array{ℜ,1}, X⁺_ij::Ƶ) where {ℜ <: Real, Ƶ <: Int}
    K::Ƶ = length(X_ij)
    if isnan(X_ij[1])
        return combinations(1:(X⁺_ij+K-1), K-1)
    else
        return combinations(cumsum(Int.(X_ij) .+ 1)[1:K-1], K-1)
    end
end

function log_marginal(X::Array{ℜ,3}, X⁺::Array{ℜ,2}, R::Ƶ; μ::ℜ=3.5, γ::ℜ=0.1) where {ℜ<:Real, Ƶ<:Int}
    I::Ƶ, J::Ƶ, K::Ƶ = size(X)

    a::ℜ, b::ℜ = I*J*K*γ, γ/μ
    α_KR, α_IKR, α_JKR = fill(a/(K*R),K,R), fill(a/(I*K*R),I,K,R), fill(a/(J*K*R),J,K,R)

    s, S_KR, S_IKR, S_JKR, S⁺ = zeros(ℜ,R), zeros(ℜ,K,R), zeros(ℜ,I,K,R), zeros(ℜ,J,K,R), sum(X⁺)
    log_PX, log_PS = -Inf, -Inf
    i::Ƶ, j::Ƶ , k::Ƶ = 0, 0, 0

    log_C::ℜ = a*log(b) - (a + S⁺)*log(b + 1.0)
    log_C += sum(lgamma.(α_KR)) - sum(lgamma.(α_IKR)) - sum(lgamma.(α_JKR))

    X_full = zeros(ℜ,I,J,K)

    for X_div = product([comb(X[i,j,:], Int(X⁺[i,j])) for i=1:I,j=1:J]...)
        if K==1
            X_full[:,:,1] .= X⁺
        else
            X_full .= 0.0
            for (ij, part) = enumerate(X_div) ## fix X
                i, j = mod(ij-1,I)+1, div(ij-1,I)+1
                X_full[i,j,1] = part[1]-1
                for k=2:K-1
                    X_full[i,j,k] = part[k] - part[k-1] - 1
                end
                X_full[i,j,K] = X⁺[i,j] + K-1 - part[K-1]
            end
        end

        if R==1
            S_KR[:,1] .= sum(X_full,dims=(1,2))[1,1,:]
            S_IKR[:,:,1] .= sum(X_full,dims=2)[:,1,:]
            S_JKR[:,:,1] .= sum(X_full,dims=1)[1,:,:]

            log_PS = log_C
            log_PS += sum(lgamma.(S_IKR .+ α_IKR)) + sum(lgamma.(S_JKR .+ α_JKR)) - sum(lgamma.(S_KR .+ α_KR))
            log_PS -= sum(lgamma.(X_full .+ 1.0))
            log_PX = logsumexp([log_PX,log_PS])
        else
            for S_div = product(map(X_ijk -> combinations(1:(X_ijk+R-1),R-1), X_full)...)
                S_KR .= α_KR
                S_IKR .= α_IKR
                S_JKR .= α_JKR

                log_PS = 0.0
                for (ijk, part) = enumerate(S_div)
                    i, j, k = mod(ijk-1,I)+1, mod(div(ijk-1,I),J)+1, div(ijk-1,I*J)+1
                    s[1] = part[1]-1
                    for r=2:R-1
                        s[r] = part[r] - part[r-1] - 1
                    end
                    s[R] = X_full[i,j,k] + R-1 - part[R-1]

                    S_KR[k,:] .+= s
                    S_IKR[i,k,:] .+= s
                    S_JKR[j,k,:] .+= s
                    log_PS -= sum(lgamma.(s .+ 1.0))
                end
                log_PS += sum(lgamma.(S_IKR)) + sum(lgamma.(S_JKR)) - sum(lgamma.(S_KR)) + log_C
                log_PX = logsumexp([log_PX,log_PS])
            end
        end
    end
    return log_PX
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
module accesorios

using Unitful

"""
    entrada_recta(m="Rennels")

Devuelve el coeficiente de pérdida para una entrada afilada a una
tubería. Hay seis fuentes disponibles; cuatro de ellas recomiendan
K = 0.5, el más reciente método "Rennels" recomienda K = 0.57, y el
método "Miller" recomienda de K = 0.51, leído en un gráfico.

![sdsdsd](assets/entrada_recta.png#thumbnail)

---

**Los campos de entrada son:**
- `m :: String, opcional `

Es el metodo a usar, uno de Rennels, Swamee, Blevins, Crane, Idelchik y Miller

---

**Salida:**
- `K :: Float `
Coeficiente de perdida `K`

"""
function entrada_recta(m="Rennels")
    if m==="Rennels"
        0.57
    elseif m in ("Swamee", "Blevins", "Crane", "Idelchik")
        0.5
    elseif m==="Miller"
        0.5092676683721356
    else
        error("Los métodos validos son: Rennels, Swamee, Blevins, Crane, Idelchik y Miller")
    end
end

end


text:
function _precompile_()
	#ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
	@assert Base.precompile(Tuple{typeof(plot),Matrix{Float64}})   # time: 11.283913
	@assert Base.precompile(Tuple{typeof(helper_multi_cols),Dict{Symbol, Any},Matrix{Float64},Bool,String,String,String,String,Bool,Vector{Bool},Vector{String},String,Vector{String},Bool,Bool})   # time: 0.2952042
	@assert Base.precompile(Tuple{typeof(finish_PS_module),Dict{Symbol, Any},Vector{String},String,Bool,Bool,Bool,Matrix{Float64},Vararg{Any, N} where N})   # time: 0.1380821
	@assert Base.precompile(Tuple{typeof(make_color_column),Dict{Symbol, Any},String,String,Int64,Int64,Int64,Bool,Bool,Bool,Vector{String},Matrix{Float64},Nothing})   # time: 0.1258531
	@assert Base.precompile(Tuple{typeof(get_marker_name),Dict{Symbol, Any},Matrix{Float64},Vector{Symbol},Bool,Bool})   # time: 0.0417101
	@assert Base.precompile(Tuple{typeof(add_opt_cpt),Dict{Symbol, Any},String,Matrix{Symbol},Char,Int64,Matrix{Float64},Nothing,Bool,Bool,String,Bool})   # time: 0.0158795
	@assert Base.precompile(Tuple{typeof(put_in_legend_bag),Dict{Symbol, Any},Vector{String},Matrix{Float64}})   # time: 0.0101455

	@assert Base.precompile(Tuple{typeof(imshow),GMTgrid{Float32, 2}})   # time: 12.84012
	@assert Base.precompile(Tuple{typeof(finish_PS_module),Dict{Symbol, Any},Vector{String},String,Bool,Bool,Bool,GMTgrid{Float32, 2},Vararg{Any, N} where N})   # time: 0.1452611
	@assert Base.precompile(Tuple{Core.kwftype(typeof(grdimage)),NamedTuple{(:show,), Tuple{Bool}},typeof(grdimage),String,GMTgrid{Float32, 2}})   # time: 0.043773
	@assert Base.precompile(Tuple{typeof(common_get_R_cpt),Dict{Symbol, Any},String,String,String,Int64,GMTgrid{Float32, 2},Nothing,Nothing,String})   # time: 0.0355232
	@assert Base.precompile(Tuple{typeof(common_shade),Dict{Symbol, Any},String,GMTgrid{Float32, 2},GMTcpt,Nothing,Nothing,String})   # time: 0.0226077
end


text:
# Autogenerated wrapper script for Lua_jll for x86_64-linux-gnu
export liblua, lua, luac

using Readline_jll
JLLWrappers.@generate_wrapper_header("Lua")
JLLWrappers.@declare_library_product(liblua, "liblua.so")
JLLWrappers.@declare_executable_product(lua)
JLLWrappers.@declare_executable_product(luac)
function __init__()
    JLLWrappers.@generate_init_header(Readline_jll)
    JLLWrappers.@init_library_product(
        liblua,
        "lib/liblua.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        lua,
        "bin/lua",
    )

    JLLWrappers.@init_executable_product(
        luac,
        "bin/luac",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
## load packages
using JuMP, Clp

## Game Definition: Kuhn Poker
# N : player indices
# Two-player, nature not included
N = [1, 2]

# I : information sets
I1 = ["JX", "QX", "KX",
      "JX-check1-bet2", "QX-check1-bet2", "KX-check1-bet2"]
I11 = view(I1, 1:3)
I12 = view(I1, 4:6)
I2 = ["XJ-check1", "XJ-bet1", "XQ-check1", "XQ-bet1", "XK-check1", "XK-bet1"]
I = (I1, I2)

# Σ : sequences
Σ1 = ["∅",
      "JX-check1", "JX-bet1", "QX-check1", "QX-bet1", "KX-check1", "KX-bet1",
      "JX-check1-fold1", "JX-check1-call1", "QX-check1-fold1", "QX-check1-call1", "KX-check1-fold1", "KX-check1-call1"]
Σ2 = ["∅",
      "XJ-check2", "XJ-bet2", "XJ-fold2", "XJ-call2",
      "XQ-check2", "XQ-bet2", "XQ-fold2", "XQ-call2",
      "XK-check2", "XK-bet2", "XK-fold2", "XK-call2",]
Σ = (Σ1, Σ2)

# g : payoff function
g = [0 0 0 0 0 0 0 0 0 0 0 0 0;
     0 0 0 0 0 -1 0 0 0 -1 0 0 0;
     0 0 0 0 0 0 0 1 -2 0 0 1 -2;
     0 1 0 0 0 0 0 0 0 -1 0 0 0;
     0 0 0 1 2 0 0 0 0 0 0 1 -2;
     0 1 0 0 0 1 0 0 0 0 0 0 0;
     0 0 0 1 2 0 0 1 2 0 0 0 0;
     0 0 0 0 0 0 -1 0 0 0 -1 0 0;
     0 0 0 0 0 0 -2 0 0 0 -2 0 0;
     0 0 -1 0 0 0 0 0 0 0 -1 0 0;
     0 0 2 0 0 0 0 0 0 0 -2 0 0;
     0 0 -1 0 0 0 -1 0 0 0 0 0 0;
     0 0 2 0 0 0 2 0 0 0 0 0 0]

 # p: probability of each branch
 p = 1/6

## Nash Equilibrium Solver
ns2 = length(Σ2)
ni1 = length(I1)
mne = Model(Clp.Optimizer)
@variable(mne, r2[1:ns2] >= 0)
@variable(mne, v1[1:ni1])
@variable(mne, v0)
@objective(mne, Min, v0)
@constraint(mne, val1, v0 - (v1[1] + v1[2] + v1[3]) >= sum(g[1, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val2, v1[1] - v1[4] >= sum(g[2, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val3, v1[1] >= sum(g[3, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val4, v1[2] - v1[5] >= sum(g[4, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val5, v1[2] >= sum(g[5, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val6, v1[3] - v1[6] >= sum(g[6, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val7, v1[3] >= sum(g[7, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val8, v1[4] >= sum(g[8, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val9, v1[4] >= sum(g[9, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val10, v1[5] >= sum(g[10, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val11, v1[5] >= sum(g[11, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val12, v1[6] >= sum(g[12, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, val13, v1[6] >= sum(g[13, s2] * r2[s2] for s2 = 1:ns2))
@constraint(mne, r2[1] == 1.0)
@constraint(mne, r2[2] + r2[3] == r2[1])
@constraint(mne, r2[4] + r2[5] == r2[1])
@constraint(mne, r2[6] + r2[7] == r2[1])
@constraint(mne, r2[8] + r2[9] == r2[1])
@constraint(mne, r2[10] + r2[11] == r2[1])
@constraint(mne, r2[12] + r2[13] == r2[1])

## Solve
optimize!(mne)
val_constraints = [val1,
                    val2,
                    val3,
                    val4,
                    val5,
                    val6,
                    val7,
                    val8,
                    val9,
                    val10,
                    val11,
                    val12,
                    val13]

## Values and Duals
r1 = dual.(val_constraints)

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
using PyCall
TB = pyimport("textblob")

str =
"""Some people think that Star Wars The Last Jedi is an excellent movie,
with perfect, flawless storytelling and impeccable acting. Others
think that it was an average movie, with a simple storyline and basic
acting. However, the reality is almost everyone felt anger and
disappointment with its forced acting and bad storytelling."""

blob = TB.TextBlob(str)
[ i.sentiment for i in blob.sentences ]

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Printf

struct Target
    value::Float32
    reward::Float32
    policy::Vector{Float32}
end

value(x::Target) = x.value
reward(x::Target) = x.reward
policy(x::Target) = x.policy

tuple(x::Target) = x.value, x.reward, x.policy

function Base.show(io::IO, x::Target; hide_reward=false, encoded_label="")
    #println(stdout, "r($(x.reward)) --> v($(x.value)), π[$(join(x.policy,","))]")
    #println(stdout, " --> (R=$(x.reward)) --> (V=$(x.value), π=[$(join(x.policy, ", "))])")
    #@printf "%5s --> "p[1] p[2] p[3] p[4] p[5] p[6] p[7] p[8] p[9]
    if !hide_reward
        print("R: ")
        @printf "%s" Int(x.reward)
        print(",  ")
    end
    v = @sprintf "%.2f" x.value
    policy_components = map(pr->begin
        @sprintf "%.2f" pr
    end, x.policy) |> x -> join(x, ", ")
    print(io, "V: $v π: [$policy_components]")
    println()
end


text:
export raise23, raise45

"""
raise23(b[], x[], npts)

Subroutine that increases a second-order B-spline curve to a third-order B-spline curve.
Returns the new control polygon vertices array and the number of its points.

# Arguments

- `b::Array{Float64}`: 2-dimensional array containing the contol polygon vertices.
- `x::Array{Float64}`: knot vector.
- `npts::Int64`: number of control polygon veritces.

# Examples

- Example 3.10 pag 98
```jldoctest
julia> x = knot(4, 2)
6-element Array{Float64,1}:
 0.0
 0.0
 1.0
 2.0
 3.0
 3.0

julia> b = [1. 2 3 4; 1 2 2 1; 1 1 1 1]
3×4 Array{Float64,2}:
 1.0  2.0  3.0  4.0
 1.0  2.0  2.0  1.0
 1.0  1.0  1.0  1.0

julia> raise23(b, x, 4)[1]
3×7 Array{Float64,2}:
 1.0  1.5  2.0  2.5  3.0  3.5  4.0
 1.0  1.5  2.0  2.0  2.0  1.5  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0

julia> raise23(b, x, 4)[2]
 7
```
"""

function raise23(b::Array{Float64}, x::Array{Float64}, npts::Int64)::Tuple{Array{Float64, 2}, Int64}

    @assert length(b) == 3 * npts ("ERROR: number of control polygon vertices and basis vector length not compatible")
    @assert length(x) == npts + 2 ("ERROR: number of control polygon vertices and knot vector length not compatible")

    mpts::Int64 = 1                     #number of vertices in the new array b*
    bstar::Array{Float64} = Float64[]

    if x[1] < x[2]                      #first point of b*
        push!(bstar, b[1,1] / 2)
        push!(bstar, b[2,1] / 2)
        push!(bstar, b[3,1] / 2)
        mpts = mpts + 1
    end

    push!(bstar, b[1,1])                   #second point of b* (or first if x1 >= x2)
    push!(bstar, b[2,1])
    push!(bstar, b[3,1])

    for i = 2 : npts                    #filling b* array
        if x[i] < x[i+1]
            mpts = mpts + 1
            push!(bstar, (b[1,i] + b[1,i-1]) / 2)
            push!(bstar, (b[2,i] + b[2,i-1]) / 2)
            push!(bstar, (b[3,i] + b[3,i-1]) / 2)
        end

        mpts = mpts + 1
        push!(bstar, b[1,i])
        push!(bstar, b[2,i])
        push!(bstar, b[3,i])
    end

    if x[npts + 2] > x[npts + 1]
        mpts = mpts + 1
        push!(bstar, b[1,npts] / 2)
        push!(bstar, b[2,npts] / 2)
        push!(bstar, b[3,npts] / 2)
    end

    bstar = reshape(bstar, 3, :)       #transforms 1-dimensional bstar array in a 3 x mpts array

    return (bstar, mpts)

end

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

"""

raise45(b[], x[], npts)

Subroutine that increases a cubic B-spline curve to a quartic B-spline curve.
Returns the new control polygon vertices array and the number of its points.

# Arguments

- `b::Array{Float64}`: 2-dimensional array containing the contol polygon vertices.
- `x::Array{Float64}`: knot vector.
- `npts::Int64`: number of control polygon veritces.

# Examples

- Example 3.11 pag 100
```jldoctest
julia> x = knot(4, 4)
8-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 1.0
 1.0
 1.0
 1.0

julia> b = [1. 2 3 4; 1 2 2 1; 1 1 1 1]
3×4 Array{Float64,2}:
 1.0  2.0  3.0  4.0
 1.0  2.0  2.0  1.0
 1.0  1.0  1.0  1.0

julia> raise45(b, x, 4)[1]
3×5 Array{Float64,2}:
 1.0  1.75  2.5  3.25  4.0
 1.0  1.75  2.0  1.75  1.0
 1.0  1.0   1.0  1.0   1.0
```
- Example 3.12 pag 100

```jldoctest
julia> k = knot(7, 4)
11-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
 1.0
 2.0
 3.0
 4.0
 4.0
 4.0
 4.0

julia> b = [1. 2 3 4 5 6 7; 1 2 2 1 0 0 1; 1 1 1 1 1 1 1]
3×7 Array{Float64,2}:
 1.0  2.0  3.0  4.0  5.0  6.0  7.0
 1.0  2.0  2.0  1.0  0.0  0.0  1.0
 1.0  1.0  1.0  1.0  1.0  1.0  1.0

julia> raise45(b, k, 7)[1] #the result in position (2,8) is 2/24 and it is correct. On the book it is 20/24 (WRONG!)
3×11 Array{Float64,2}:
 1.0  1.75  2.25  2.95833  3.5  4.0  4.5  5.04167    5.75  6.25  7.0
 1.0  1.75  2.0   1.91667  1.5  1.0  0.5  0.0833333  0.0   0.25  1.0
 1.0  1.0   1.0   1.0      1.0  1.0  1.0  1.0        1.0   1.0   1.0
```
"""
function raise45(b::Array{Float64}, x::Array{Float64}, npts::Int64)::Tuple{Array{Float64,2},Int64}

    @assert length(b) == 3 * npts ("ERROR: number of control polygon vertices and basis vector length not compatible")
    @assert length(x) == npts + 4 ("ERROR: number of control polygon vertices and knot vector length not compatible")
    #mpts = npts + s + ord + 1 = npts + 5 + s ?????

    mpts::Int64 = 0
    bstar::Array{Float64} = Float64[]

    #coordinates of the first two points in b*
    for i = 1 : 3
        push!(bstar, b[i,1]) #b*[1]
    end
    for  i = 1 : 3
        push!(bstar, (b[i,1] + 3 * b[i,2]) / 4) #[b*[2]]
    end
    mpts = mpts + 2

    if npts == 4
        for i = 1 : 3
            push!(bstar, (b[i,2] + b[i,3]) / 2) #b*[3]
        end
        mpts = mpts + 1
    else
        for i = 1 : 3
            push!(bstar, (3 * b[i,2] + b[i,3]) / 4) #b*[3]
        end
        mpts = mpts + 1

        if npts == 5
            for i = 1 : 3
                push!(bstar, (b[i,2] + 6 * b[i,3] + b[i,4]) / 8) #b*[4]
            end
            mpts = mpts + 1
        else
            for i = 1 : 3
                push!(bstar, (3 * b[i,2] + 19 * b[i,3] + 2 * b[i,4]) / 24) #b*[4]
            end
            for i = 1 : 3
                push!(bstar, (b[i,3] + b[i,4]) / 2) #b*[5]
            end
            mpts = mpts + 2

            if npts > 6
                n4 = npts - 4
                for i = 3 : n4
                    for j = 1 : 3
                        push!(bstar, (b[j,i] + 10 * b[j,i+1] + b[j,i+2]) / 12) #b*[2i]
                    end
                    for j = 1 : 3
                        push!(bstar, (b[j,i+1] + b[j,i+2]) / 2) #b*[2(i+1)]
                    end
                    mpts = mpts + 2
                end
            end

            for i = 1 : 3
                push!(bstar, (3 * b[i,npts-1] + 19 * b[i,npts-2] + 2 * b[i,npts-3]) / 24) #b*[mpts-3]
            end
            mpts = mpts + 1
        end

        for i = 1 : 3
            push!(bstar, (3 * b[i,npts-1] + b[i,npts-2]) / 4) #b*[mpts-2]
        end
        mpts = mpts + 1
    end

    for i = 1 : 3
        push!(bstar, (b[i,npts] + 3 * b[i,npts-1]) / 4) #b*[mpts-1]
    end
    for i = 1 : 3
        push!(bstar, b[i,npts]) #b*[mpts]
    end
    mpts = mpts + 2

    bstar = reshape(bstar,3,:)  #transforms 1-dimensional bstar array in a 3 x mpts array

    return (bstar, mpts)

end


text:
@windows_only using RPMmd

function add_to_juliarc(line::String)
    f=open("$(ENV["HOME"])/juliarc.jl","a")
    println(f)
    println(f,line)
    close(f)
end

function hasbin{T<:String}(bin::String,extra::Vector{T}=[])
    path = split(ENV["PATH"],":")
    for p in path
        f = joinpath(p,bin)
        if isfile(f)
            return realpath(p)
        end
    end
    for p in extra
        if isdir(p)
            f = joinpath(p,bin)
            if isfile(f)
                #ENV["PATH"] = ENV["PATH"]*":"*p
                return realpath(p)
            end
        end
    end
    return ""
end

#example: SystemLibInstall("libexpat", [("brew","expat"), ("port","expat")])
typealias PackageInstallMap Dict{Symbol,Union(String,Cmd)}
type SystemLibInstall <: BuildStep
    filenames::Vector{String}
    known_map::PackageInstallMap
end
SystemLibInstall(s::Vector{String},known_map::PackageInstallMap) = SystemLibInstall(String[s,], PackageInstallMap())
SystemLibInstall(s::String) = SystemLibInstall(s, PackageInstallMap())
function SystemLibInstall(s,known_map::Vector)
    known = PackageInstallMap()
    for (manager,cmd) in known_map
        known[manager] = cmd
    end
    SystemLibInstall(s,known)
end

let brew_path=""
global brew
function brew(cmd::Symbol,arg=nothing)
    if cmd == :exists
        brew_path = hasbin("brew",["/usr/local/bin",])
        if brew_path != ""
            brewlib = abspath(brew_path,"..","lib")
            if !contains(LIBRARY_PATH, brewlib)
                push!(Sys.DL_LOAD_PATH, brewlib)
                add_to_juliarc("push!(Sys.DL_LOAD_PATH, \"$(escape(brewlib))\")")
            end
            return true
        end
        return false
    elseif cmd == :whatprovides
        return nothing
    elseif cmd == :install
        return run(`$brew_path install $arg`)
    end
end
end

let port_path = ""
global port
function port(cmd::Symbol,arg=nothing)
    if cmd == :exists
        port_path = hasbin("port",["/opt/local/bin",])
        if port_path != ""
            portlib = abspath(port_path,"..","lib")
            if !contains(LIBRARY_PATH, portlib)
                push!(Sys.DL_LOAD_PATH, portlib)
                add_to_juliarc("push!(Sys.DL_LOAD_PATH, \"$(escape(portlib))\")")
            end
            return true
        end
        return false
    elseif cmd == :whatprovides
        return nothing
    elseif cmd == :install
        print("System dependency install using `sudo`. You will be prompted to enter password.")
        return run(`sudo $port_path install $arg`)
    end
end
end

function rpmmd(cmd::Symbol,arg=nothing)
    if cmd == :exists
        try
            eval(Expr(:toplevel,:using,:RPMmd))
            return true
        catch
            return false
        end
    elseif cmd == :whatprovides
        return RPMmd.whatprovides("/$arg")
    elseif cmd == :install
        return RPMmd.install(arg)
    end
end

function yum(cmd::Symbol,arg=nothing)
    if cmd == :exists
        return hasbin("yum")
    elseif cmd == :whatprovides
        #pkgs = readall(`yum whatprovides */$arg`, STDIN)
        #TODO: scan for package name(s)
        return nothing
    elseif cmd == :install
        print("System dependency install using `sudo`. You will be prompted to enter password.")
        return run(`sudo yum install $arg`)
    end
end

function apt(cmd::Symbol,arg=nothing)
    if cmd == :exists
        return hasbin("apt-get")
    elseif cmd == :whatprovides
        return unique([split(pkg,":",2) for pkg in split(readall("apt-file search $arg"),'\n')])
    elseif cmd == :install
        print("System dependency install using `sudo`. You will be prompted to enter password.")
        return run("sudo apt-get install $arg")
    end
end

global system_man = nothing
function default_pkg_manager()
    global system_man
    if system_man !== nothing
        return system_man::Function
    end
    pkg_mangers = []
    @osx_only append!(pkg_managers,[port,brew])
    @windows_only push!(pkg_managers,[rpmmd,])
    @linux_only push!(pkg_managers,[apt,yum])
    for man in pkg_managers
        if man(:exists)
            system_man = man
            return man
        end
    end
    return nothing
end

function run(si::SystemLibInstall)
    man = default_pkg_manager()
    lib = Base.find_library(si.filenames)
    if lib != ""
        return true
    end
    if man === nothing
        return false
    end
    provides = get(si.known_map, string(man), nothing)
    if provides === nothing
        for file in si.filenames
            provides = man(:whatprovides, file)
            if provides !== nothing
                break
            end
        end
    end
    if provides === nothing
        return false
    end
    return man(:install, provides)
end


text:
# EKSPF uses FPFState
const EKSPFState = FPFState

struct EKSPF{TS<:HiddenStateModel, TO<:ObservationModel} <: AbstractFilteringAlgorithm{ContinuousTime, ContinuousTime}
    state_model::TS
    obs_model::TO
    N::Int
    function EKSPF(st_mod::HiddenStateModel{S1, ContinuousTime}, ob_mod::ObservationModel{S1, S2, ContinuousTime}, N::Int) where {S1, S2}
        if ob_mod isa CountingObservationModel
            return new{typeof(st_mod), typeof(ob_mod)}(st_mod, ob_mod, N)
        else
            error("Error: the EKSPF algorithm is only defined for counting process observations.")
        end
    end
end

#####################
###### METHODS ######
#####################

initial_condition(ekspf::EKSPF) = initial_condition(state_model(ekspf))
no_of_particles(ekspf::EKSPF) = ekspf.N
state_model(ekspf::EKSPF) = ekspf.state_model
obs_model(ekspf::EKSPF) = ekspf.obs_model
gain_estimation_method(ekspf::EKSPF) = ConstantGainEKSPF()

function Base.show(io::IO, ::MIME"text/plain", ekspf::EKSPF)
    print(io, "Ensemble Kushner-Stratonovich-Poisson filter algorithm
    hidden state model:                     ", state_model(ekspf),"
    observation model:                      ", obs_model(ekspf),"
    ensemble size:                          ", no_of_particles(ekspf))
end

function Base.show(io::IO, ekspf::EKSPF)
    print(io, "EKSPF with ", no_of_particles(ekspf)," particles")
end

######################
### MAIN ALGORITHM ###
######################

function initialize(ekspf::EKSPF)
    return EKSPFState(ekspf)
end

function update!(filter_state::EKSPFState, filt_algo::EKSPF, dY, dt)
    propagate!(filter_state, filt_algo, dt)
    assimilate!(dY, filter_state, filt_algo, dt)
end

function propagate!(filter_state::EKSPFState, filt_algo::EKSPF, dt)
    propagate!(filter_state.ensemble, filt_algo.state_model, dt)
end

function assimilate!(dN, filter_state::EKSPFState, filt_algo::EKSPF, dt)
    ensemble  = filter_state.ensemble
    eq        = filter_state.eq
    method    = ConstantGainEKSPF()

    error     = dN .- dt * eq.H
    solve!(eq, method)

    applygain!(ensemble, eq, error)
    update!(eq, ensemble)
    return ensemble.positions
end

################################
### CONVENIENCE CONSTRUCTORS ###
################################

function EKSPFState(ekspf::EKSPF)
    N   = no_of_particles(ekspf)
    ens = UnweightedParticleEnsemble(hcat([initialize(ekspf.state_model) for i in 1:N]...))
    eq  = PoissonEquation(observation_function(ekspf.obs_model), ens)
    return EKSPFState(ens, eq)
end

function EKSPF(filt_prob::AbstractFilteringProblem, N)
    st_mod = state_model(filt_prob)

    if obs_model(filt_prob) isa CountingObservationModel
        ob_mod = obs_model(filt_prob)
    else
        error("Error: the EKSPF algorithm is only defined for counting process observations.")
    end

    return EKSPF(st_mod, ob_mod, N)
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
# this does basis promotion, unlike Zygote's default for AbstractArrays
Zygote.accum(a::Field, b::Field) = a+b
# this may create a LazyBinaryOp, unlike Zygote's
Zygote.accum(a::FieldOp, b::FieldOp) = a+b

# constant functions, as far as AD is concerned
@nograd ProjLambert
@nograd fieldinfo
@nograd hasfield

# AD for Fourier Fields can be really subtle because such objects are
# still supposed to represent a pixel-space field, despite that
# they're stored as the half-plane real FFT coefficients. this leads
# to needing factors of Npix(=Ny*Nx) and rfft_degeneracy_fac (which
# gives weight of 2 to coefficient which would be there twice in the
# full-plane FFT) in a few cases below to make the logic all work. the
# Npix factor is handled by the Zfac function.
Zfac(::SpatialBasis{Map},     proj::FlatProj) = 1
Zfac(::SpatialBasis{Fourier}, proj::FlatProj) = proj.Ny * proj.Nx
Zfac(L::DiagOp{<:Field{B}}) where {B} = Zfac(B(), L.diag.metadata)

## constructors and getproperty
# these are needed to allow you to "bail out" to the underlyig arrays
# via f.arr and reconstruct them with FlatField(f.arr, f.metadata).
@adjoint function (::Type{F})(arr::A, metadata::M) where {B<:SpatialBasis{Map},M<:FieldMetadata,T,A<:AbstractArray{T},F<:BaseField{B}}
    F(arr, metadata), Δ -> (Δ.arr, nothing)
end
@adjoint function (::Type{F})(arr::A, metadata::M) where {B<:SpatialBasis{Fourier},M<:FieldMetadata,T,A<:AbstractArray{T},F<:BaseField{B}}
    F(arr, metadata), Δ -> (Δ.arr .* adapt(Δ.storage, T.(rfft_degeneracy_fac(metadata.Ny) ./ Zfac(B(), metadata))), nothing)
end
# the factors here need to cancel the ones in the corresponding constructors above
@adjoint function Zygote.literal_getproperty(f::BaseField{B}, ::Val{:arr}) where {B<:SpatialBasis{Map}}
    getfield(f,:arr), Δ -> (BaseField{B}(Δ, f.metadata),)
end
@adjoint function Zygote.literal_getproperty(f::BaseField{B,M,T}, ::Val{:arr}) where {B<:SpatialBasis{Fourier},M,T}
    getfield(f,:arr), Δ -> (BaseField{B}(Δ ./ adapt(typeof(Δ), T.(rfft_degeneracy_fac(f.Ny) ./ Zfac(B(), f.metadata))), f.metadata),)
end
Zygote.accum(f::BaseField, nt::NamedTuple{(:arr,:metadata)}) = (@assert(isnothing(nt.arr)); f)

# FieldTuple
@adjoint (::Type{FT})(fs) where {FT<:FieldTuple} = FT(fs), Δ -> (Δ.fs,)
@adjoint Zygote.literal_getproperty(f::FieldTuple, ::Val{:fs}) = getfield(f,:fs), Δ -> (FieldTuple(map((f,f̄) -> isnothing(f̄) ? zero(f) : f̄, getfield(f,:fs), Δ)),)

# BatchedReals
@adjoint Zygote.literal_getproperty(br::BatchedReal, ::Val{:vals}) = getfield(br,:vals), Δ -> (batch(real.(Δ)),)

## Field algebra

# Zygote has lots of rules for AbstractVectors / AbstractMatrices that
# don't quite work right due to the auto-basis conversions done for
# Fields, or which work right but trigger scalar indexing (thus don't
# work on GPU). this leads us to need a few more custom rules below
# than might be ideal, although its not too bad. see also:
# https://discourse.julialang.org/t/how-to-deal-with-zygote-sometimes-pirating-its-own-adjoints-with-worse-ones

# ℝᴺ -> ℝ¹
@adjoint sum(f::Field{B}) where {B} = sum(f), Δ -> (Δ*one(f),)
@adjoint norm(f::Field) = Zygote.pullback(f->sqrt(dot(f,f)), f)
@adjoint dot(f::Field{B1}, g::Field{B2}) where {B1,B2} = dot(f,g), Δ -> (Δ*B1(g), Δ*B2(f))
@adjoint (*)(f::Adjoint{<:Any,<:Field}, g::Field) = Zygote.pullback((f,g)->dot(f',g),f,g)
# ℝᴺˣᴺ -> ℝ¹
@adjoint logdet(L::ParamDependentOp, θ) = Zygote._pullback(θ->logdet(L(θ)), θ)
@adjoint logdet(L::DiagOp) = logdet(L), Δ -> (Δ * Zfac(L) * pinv(L)',)

# basis conversion
@adjoint (::Type{B})(f::Field{B′}) where {B<:Basis, B′} = B(f), Δ -> (B′(Δ),)

# algebra
@adjoint (+)(f::Field{B1}, g::Field{B2}) where {B1,B2} = f+g, Δ -> (B1(Δ), B2(Δ))

@adjoint (*)(a::Real, L::DiagOp) = a*L, Δ -> (tr(L'*Δ)/Zfac(L), a*Δ)
@adjoint (*)(L::DiagOp, a::Real) = a*L, Δ -> (a*Δ, tr(L'*Δ)/Zfac(L))

# operators
@adjoint function (*)(D::DiagOp{<:Field{B}}, v::Field{B′}) where {B,B′}
    D*v, Δ -> (Diagonal(B(Δ) .* conj.(B(v))), B′(D'*Δ))
end
@adjoint function (\)(D::DiagOp{<:Field{B}}, v::Field{B′}) where {B,B′}
    z = D \ v
    function back(Δ)
        v̄ = D' \ Δ
        -Diagonal(B(v̄) .* conj.(B(z))), B′(v̄)
    end
    z, back
end
# this makes it so we only have to define a L*f and L\f adjoint (like
# above) for any new operators, and we get f'*L and f'/L for free
# without necessarily needing to fully implement the AbstractMatrix
# interface for L
@adjoint *(f::Adjoint{<:Any,<:Field}, L::Union{DiagOp,ImplicitOp}) = Zygote.pullback((f,L)->(L'*f')', f, L)
@adjoint /(f::Adjoint{<:Any,<:Field}, L::Union{DiagOp,ImplicitOp}) = Zygote.pullback((f,L)->(L'\f')', f, L)
# special case for some ops which are constant by definition
@adjoint *(L::Union{FuncOp,DiagOp{<:∇diag}}, f::Field{B}) where {B} = L*f, Δ->(nothing, B(L'*Δ))

## FieldVectors

# following two definitions are almost definitely not totally right w.r.t
# putting stuff in the correct basis, although they're working for everything
# I've needed thus far

@adjoint function *(x::FieldOrOpRowVector, y::FieldVector)
    z = x * y
    # when x is a vector of Fields
    back(Δ) = (Δ * y', x' * Δ)
    # when x is a vector of Diagonals. in this case, Δ * basis(Δ)(y)' create an
    # OuterProdOp in the same basis as the Diagonals in x
    back(Δ::Field{B}) where {B} = (Δ * basis(Δ)(y)'), (x' * Δ)
    z, back
end

@adjoint function *(A::FieldOrOpMatrix, x::FieldOrOpVector)
    z = A * x
    back(Δ::FieldVector) = (B=basis(eltype(z)); (B(Δ) * B(x)', A' * Δ))
    back(Δ::FieldOrOpVector) = (Δ * x', A' * Δ) # not sure if anything special needed here?
    back(Δ) = (Δ * x', A' * Δ)
    z, back
end

# don't know why Zygote's default adjoint for this breaks in various ways but this is simple enough
@adjoint +(I::UniformScaling, L::Union{FieldOp, FieldOrOpMatrix}) = I+L, Δ->(nothing, Δ)

# Zygote/lib/array.jl:311 would suggest this should be:
#    M⁻¹, Δ->(-M⁻¹' * Δ * M⁻¹' + (- M * M⁻¹ * Δ' * M⁻¹ * M⁻¹' + Δ' * M⁻¹ * M⁻¹') + (M⁻¹' * M⁻¹ * Δ' - M⁻¹' * M⁻¹ * Δ' * M⁻¹ * M),)
# I haven't derived their version, but numerically the one gives the right answer where as their doesn't...
@adjoint function pinv(L::Union{FieldOp, FieldOrOpMatrix})
    L⁻¹ = pinv(L)
    L⁻¹, Δ->(-L⁻¹' * Δ * L⁻¹',)
end

@adjoint sqrt(L::DiagOp) = (z=sqrt(L);), Δ -> ((pinv(z)/2)'*Δ,)

# some stuff which arguably belongs in Zygote or ChainRules
# see also: https://github.com/FluxML/Zygote.jl/issues/316

@adjoint broadcasted(::typeof(\), x ::Numeric, y::Numeric) =
    broadcast(\, x, y), Δ -> (nothing, unbroadcast(x, @. -Δ*y/x^2), unbroadcast(y, @. Δ/x))

@adjoint (::Type{SA})(tup) where {SA<:SArray} = SA(tup), Δ->(tuple(Δ...),)

# workaround for https://github.com/FluxML/Zygote.jl/issues/686
@static if versionof(Zygote) > v"0.4.15"
    Zygote._zero(xs::StaticArray, T) = SizedArray{Tuple{size(xs)...},Union{T,Nothing}}(map(_->nothing, xs))
end

# workaround for Zygote not working through cat when dims is a Val
# adapted from solution by Seth Axen
# see https://github.com/FluxML/Zygote.jl/pull/881
@adjoint function cat(Xs::AbstractArray...; dims)
    cat(Xs...; dims = dims), Δ -> begin
        start = ntuple(_ -> 0, ndims(Δ))
        catdims = Base.dims2cat(dims)
        dXs = map(Xs) do x
            move = ntuple(d -> (d<=length(catdims) && catdims[d]) ? size(x,d) : 0, ndims(Δ))
            x_in_Δ = ntuple(d -> (d<=length(catdims) && catdims[d]) ? (start[d]+1:start[d]+move[d]) : Colon(), ndims(Δ))
            start = start .+ move
            dx = Δ[x_in_Δ...]
        end
    end
end

@adjoint adapt(to, x::A) where {A<:AbstractArray} = adapt(to, x), Δ -> (nothing, adapt(A, Δ))

# finite difference Hessian using Zygote gradients
# todo: delete, just use FiniteDifferences
function hessian(f, xs::Vector; ε=1f-3)
    hcat(finite_difference(xs->vcat(gradient(f,xs)[1]...),xs,ε=ε)...)
end

function finite_difference(f, xs::Vector; ε=1f-3, progress=false)
    @showprogress (progress ? 1 : Inf) map(1:length(xs)) do i
        xs₊ = copy(xs); xs₊[i] += ε
        xs₋ = copy(xs); xs₋[i] -= ε
        (f(xs₊) .- f(xs₋)) ./ (2ε)
    end
end


text:
using ModelingToolkit
using NeuralPDE
using Flux
using DiffEqFlux
using GalacticOptim

# using CUDA
# CUDA.allowscalar(false)

function NeuralPDE.PhysicsInformedNN(dx, chain::Chain, initθ; autodiff=false, strategy = GridTraining(), kwargs...)
   indim = size(nnet.layers[1].W, 2)
   outdim = size(chain.layers[end].W, 1)
   isuinplace = outdim == 1
   φ = get_phi(chain, isuinplace)
   derivative = get_derivative(indim, φ, autodiff, isuinplace)
   PhysicsInformedNN(dx, chain, initθ, φ, autodiff, derivative, strategy, kwargs)
end

function fitPDESystem!(K, α, µ, ρ₀, ρ₀ˢ, ρ₀ˡ,
                       F₁, F₂,
                       tspan, x₁span, x₂span,
                       Δt, Δx₁, Δx₂;
                       maxiters=1000, strategy = GridTraining())

   @parameters t, x₁, x₂, θ

    u₁ = Variable{Float32}(:u₁)
    u₂ = Variable{Float32}(:u₂)
    v₁ = Variable{Float32}(:v₁)
    v₂ = Variable{Float32}(:v₂)
   σ₁₁ = Variable{Float32}(:σ₁₁)
   σ₁₂ = Variable{Float32}(:σ₁₂)
   σ₂₂ = Variable{Float32}(:σ₂₂)
     p = Variable{Float32}(:p)

   @derivatives Dt'~t
   @derivatives Dx₁'~x₁
   @derivatives Dx₂'~x₂

   # System of PDEs
   equations = [
      Dt(u₁(t,x₁,x₂,θ)) + ( Dx₁(σ₁₁(t,x₁,x₂,θ)) + Dx₂(σ₁₂(t,x₁,x₂,θ)) )/ρ₀ˢ + Dx₁(p(t,x₁,x₂,θ))/ρ₀ - F₁(t,x₁,x₂) ~ 0,
      Dt(u₂(t,x₁,x₂,θ)) + ( Dx₁(σ₁₂(t,x₁,x₂,θ)) + Dx₂(σ₂₂(t,x₁,x₂,θ)) )/ρ₀ˢ + Dx₂(p(t,x₁,x₂,θ))/ρ₀ - F₂(t,x₁,x₂) ~ 0,
      Dt(v₁(t,x₁,x₂,θ)) + Dx₁(p(t,x₁,x₂,θ))/ρ₀ - F₁(t,x₁,x₂) ~ 0,
      Dt(v₂(t,x₁,x₂,θ)) + Dx₂(p(t,x₁,x₂,θ))/ρ₀ - F₂(t,x₁,x₂) ~ 0,
      Dt(σ₁₁(t,x₁,x₂,θ)) + (ρ₀ˡ/ρ₀ * K + 4µ/3)Dx₁(u₁(t,x₁,x₂,θ)) + (ρ₀ˡ/ρ₀ * K - 2µ/3)Dx₂(u₂(t,x₁,x₂,θ)) - ρ₀ˢ/ρ₀ * K * ( Dx₁(v₁(t,x₁,x₂,θ)) + Dx₂(v₂(t,x₁,x₂,θ)) ) ~ 0,
      Dt(σ₁₂(t,x₁,x₂,θ)) + µ * ( Dx₂(u₁(t,x₁,x₂,θ)) + Dx₁(u₂(t,x₁,x₂,θ)) ) ~ 0,
      Dt(σ₂₂(t,x₁,x₂,θ)) + (ρ₀ˡ/ρ₀ * K - 2µ/3)Dx₁(u₁(t,x₁,x₂,θ)) + (ρ₀ˡ/ρ₀ * K + 4µ/3)Dx₂(u₂(t,x₁,x₂,θ)) - ρ₀ˢ/ρ₀ * K * ( Dx₁(v₁(t,x₁,x₂,θ)) + Dx₂(v₂(t,x₁,x₂,θ)) ) ~ 0,
      Dt(p(t,x₁,x₂,θ)) - (K - α*ρ₀*ρ₀ˢ)*( Dx₁(u₁(t,x₁,x₂,θ)) + Dx₂(u₂(t,x₁,x₂,θ)) ) + α*ρ₀*ρ₀ˡ*( Dx₁(v₁(t,x₁,x₂,θ)) + Dx₂(v₂(t,x₁,x₂,θ)) ) ~ 0
   ]

   # Initial and boundary conditions
   boundconds = [
        # initial conditions at t = 0
       u₁(tspan.lower, x₁, x₂, θ) ~ 0,
       u₂(tspan.lower, x₁, x₂, θ) ~ 0,
       v₁(tspan.lower, x₁, x₂, θ) ~ 0,
       v₂(tspan.lower, x₁, x₂, θ) ~ 0,
      σ₁₁(tspan.lower, x₁, x₂, θ) ~ 0,
      σ₁₂(tspan.lower, x₁, x₂, θ) ~ 0,
      σ₂₂(tspan.lower, x₁, x₂, θ) ~ 0,
        p(tspan.lower, x₁, x₂, θ) ~ 0,
      # boundary conditions at x₁ = x₁⁻
       u₁(t, x₁span.lower, x₂, θ) ~ 0,
       u₂(t, x₁span.lower, x₂, θ) ~ 0,
       v₁(t, x₁span.lower, x₂, θ) ~ 0,
       v₂(t, x₁span.lower, x₂, θ) ~ 0,
      σ₁₁(t, x₁span.lower, x₂, θ) ~ 0,
      σ₁₂(t, x₁span.lower, x₂, θ) ~ 0,
      σ₂₂(t, x₁span.lower, x₂, θ) ~ 0,
        p(t, x₁span.lower, x₂, θ) ~ 0,
      # boundary conditions at x₁ = x₁⁺
       u₁(t, x₁span.upper, x₂, θ) ~ 0,
       u₂(t, x₁span.upper, x₂, θ) ~ 0,
       v₁(t, x₁span.upper, x₂, θ) ~ 0,
       v₂(t, x₁span.upper, x₂, θ) ~ 0,
      σ₁₁(t, x₁span.upper, x₂, θ) ~ 0,
      σ₁₂(t, x₁span.upper, x₂, θ) ~ 0,
      σ₂₂(t, x₁span.upper, x₂, θ) ~ 0,
        p(t, x₁span.upper, x₂, θ) ~ 0,
      # boundary conditions at x₂ = x₂⁻
      σ₁₂(t, x₁, x₂span.lower, θ) ~ 0,
      # boundary conditions at x₂ = x₂⁺
       u₁(t, x₁, x₂span.upper, θ) ~ 0,
       u₂(t, x₁, x₂span.upper, θ) ~ 0,
       v₁(t, x₁, x₂span.upper, θ) ~ 0,
       v₂(t, x₁, x₂span.upper, θ) ~ 0,
      σ₁₁(t, x₁, x₂span.upper, θ) ~ 0,
      σ₁₂(t, x₁, x₂span.upper, θ) ~ 0,
      σ₂₂(t, x₁, x₂span.upper, θ) ~ 0,
        p(t, x₁, x₂span.upper, θ) ~ 0
   ]
   if ρ₀ˡ == 0
      push!(boundconds, σ₂₂(t, x₁, x₂span.lower, θ) ~ -p(t, x₁, x₂span.lower, θ))
   else
      push!(boundconds, σ₂₂(t, x₁, x₂span.lower, θ) ~ 0, p(t, x₁, x₂span.lower, θ) ~ 0)
   end

   # Time and space domains
   domain = [t ∈ tspan, x₁ ∈ x₁span, x₂ ∈ x₂span]

   # Neural network
   nnet = FastChain(FastDense(3,  16, Flux.selu),
                    FastDense(16, 16, Flux.σ),
                    FastDense(16, 16, Flux.σ),
                    FastDense(16, 8)) #|> gpu
   # nnet = Chain(Dense(3,  16, Flux.selu),
   #              Dense(16, 16, Flux.σ),
   #              Dense(16, 16, Flux.σ),
   #              Dense(16, 8)) #|> gpu

   # Initial neural network training parameters
   initθ = initial_params(nnet) #|> gpu

   discretization = PhysicsInformedNN([Δt, Δx₁, Δx₂], nnet, initθ; strategy)

   pde_system = PDESystem(equations, boundconds, domain, [t, x₁, x₂], [u₁, u₂, v₁, v₂, σ₁₁, σ₁₂, σ₂₂, p])

   problem = NeuralPDE.discretize(pde_system, discretization)

   optimizer = Flux.ADAM(0.1f0)
   callback = function (p, l)
      println("Current loss is: $l")
      return false
   end

   @time res = GalacticOptim.solve(problem, optimizer; progress=false, cb=callback, maxiters)

   ϕ = discretization.phi
   θopt = res.minimizer

   return function sol(t, x₁, x₂)
      ϕ([t, x₁, x₂], θopt)
   end
end

# Specify constants
const K, α, µ, ρ₀, ρ₀ˢ, ρ₀ˡ = let
   ρ₀ᶠˢ = 1400
   ρ₀ᶠˡ = 997
   cp₁ = 2000
   cp₂ = 1300
   cs = 1300
   d₀ = 0.3

   ρ₀ˢ = (1 - d₀)ρ₀ᶠˢ
   ρ₀ˡ = d₀ * ρ₀ᶠˡ
   ρ₀ = ρ₀ˢ + ρ₀ˡ
   µ  = ρ₀ˢ * cs^2
   K  = (cp₁^2 + cp₂^2 - 8/3*ρ₀ˡ/ρ₀*cs^2 - √((cp₁^2 - cp₂^2)^2 - 64/9*cs^4*ρ₀ˡ*ρ₀ˢ/ρ₀^2) ) * ρ₀ * ρ₀ˢ/(2ρ₀ˡ)
   α₃ = (cp₁^2 + cp₂^2 - 8/3*ρ₀ˢ/ρ₀*cs^2 - √((cp₁^2 - cp₂^2)^2 - 64/9*cs^4*ρ₀ˡ*ρ₀ˢ/ρ₀^2) ) / (2ρ₀^2)
   α  = ρ₀ * α₃ + K/ρ₀^2

   Float32.((K, α, µ, ρ₀, ρ₀ˢ, ρ₀ˡ))
end

# Specify force functions F₁ and F₂
F₁(t, x₁, x₂) = 100exp(-t^2 - x₁^2 - x₂^2) * sin(t)
F₂(t, x₁, x₂) = 1000exp(-t^2 - x₁^2 - x₂^2) * cos(t)

# Specify variable intervals
 tspan = IntervalDomain(0.f0, 60.f0)
x₁span = IntervalDomain(-100.f0, 100.f0)
x₂span = IntervalDomain(0.f0, 100.f0)

# Specify discretization
Δt  = 10.f0
Δx₁ = 20.f0
Δx₂ = 10.f0

# Fit the neural network to the PDE system
w = fitPDESystem!(K, α, µ, ρ₀, ρ₀ˢ, ρ₀ˡ,
                  F₁, F₂,
                  tspan, x₁span, x₂span,
                  Δt, Δx₁, Δx₂; maxiters=1000)

w(tspan.lower, x₁span.lower, x₂span.lower)

w(1, 2, 3)

max_stars_repo_name: AzamatB/PoroelasticityPDESystem

max_issues_repo_name: AzamatB/PoroelasticityPDESystem

max_forks_repo_name: AzamatB/PoroelasticityPDESystem


text:
struct RecordHeader
    ts_sec::UInt32
    ts_usec::UInt32
    incl_len::UInt32
    orig_len::UInt32
end

function Base.bswap(x::RecordHeader)
    RecordHeader(
        bswap(x.ts_sec),
        bswap(x.ts_usec),
        bswap(x.incl_len),
        bswap(x.orig_len))
end


text:
using RiskMeasures
using Test

@testset "helper.jl" begin
    # Write your tests here.
    @test pnorm(2,[0.5 0.5],[1 2]) ≈ 1.5811 atol=0.01

    @test ontoSimplex([0.3 0.3 0.3 0.3]) == [0.25 0.25 0.25 0.25]
    @test ontoSimplex([-0.3 -0.3 0.3 0.3]) == [0.0 0.0 0.5 0.5]
    @test goldenSearch(x -> x^2, 1.0)[2] < 1E-16
    @test goldenSearch(x -> x^2, -1.0)[2] < 1E-16

    @test checkSpectral(x -> x) == false
    @test checkSpectral(x -> 2 .- 2 .*x) == false
    @test checkSpectral(x -> 0 .*x .+ 1.0) == true

    @test eval(math_expr(:+,1,2)) == 3
    @test eval(math_expr(:+,1)) == 1
    @test eval(add_expr(1, 2)) == 3
end

@testset "basic_RM.jl" begin
    # Write your tests here.
    states = [1.0, 2.0, 3.0, 4.0]; prob = [0.2, 0.4, 0.2, 0.2]
    @test Expectation([1.0, 2.0, 3.0, 4.0],[0.0, 0.0, 0.0, 0.0]) == sum([1 2 3 4])/4
    @test Expectation(states,prob) == 0.2*1 + 0.4*2 + 0.2*3 + 0.2*4

    @test mSD(states,prob,0.0f0,2.0f0) == Expectation(states,prob)
    @test mSD([1.0, 1.0], [0.0, 0.0],0.0f0,2.0f0) == 1
    @test mSD(states,prob,1.0f0,2.0f0) ≈ 3.16419 atol = 0.001
    @test mSD(states,prob,1.0f0,2.0f0) < mSD(states,prob,1.0f0,3.0f0)

    @test VaR(states,prob,0.0f0) == minimum(states)
    @test CTE(states,prob,0.0f0) == Expectation(states,prob)
end

@testset "valueRisk.jl" begin
    # Write your tests here.
    states = [1.0, 2.0, 3.0, 4.0]; prob = [0.2, 0.4, 0.2, 0.2]

    @test EVaR2(states, prob,0.0f0)[1] == Expectation(states,prob)
    @test EVaR2(states, prob,0.5f0)[1] ≈ 3.413183 atol = 0.001
    @test EVaR2([1.0, 1.0], [0.0, 0.0],0.5f0)[1] == 0
    @test EVaR(states, prob,0.0f0)[1] == Expectation(states,prob)
    @test EVaR(states, prob,0.5f0)[1] ≈ 3.413183 atol = 0.001
    @test EVaR([1.0, 1.0], [0.0, 0.0],0.5f0)[1] == 0.0
    @test EVaR(states, prob,5.5f0)[1] == maximum(states)
    @test AVaR(states,prob, 0.0f0)[1] == Expectation(states,prob)
    @test AVaR(states,prob, 0.5f0)[1] ≈ CTE(states,prob,0.5f0) atol = 0.0001
end

@testset "genRM.jl" begin
    # Write your tests here.
    states = [1.0, 2.0, 3.0, 4.0]; prob = [0.2, 0.4, 0.2, 0.2]

    @test entropic(states, prob,1.0f0) ≈ 2.91430 atol = 0.001
    @test entropic([1.0, 1.0], [0.0, 0.0],0.5f0) == 1.0
    @test entropic([1.0, 1.0], [0.0, 0.0],0.0f0) === nothing

    @test meanVariance([1.0, 1.0], [0.0, 0.0],-1.0f0) === nothing
    @test meanVariance([1.0, 1.0], [0.0, 0.0],1.0f0) == 1.0

    @test meanDeviation([1.0, 1.0], [0.0, 0.0],-1.0f0,2.0f0) === nothing
    @test meanDeviation([1.0, 1.0], [0.0, 0.0],1.0f0,2.0f0) == 1.0
    @test meanDeviation([1.0, 1.0], [0.0, 0.0],-1.0f0,0.5f0) === nothing
    @test meanDeviation([1.0, 1.0], [0.0, 0.0],1.0f0,0.5f0) === nothing

    @test meanSemiVariance([1.0, 1.0], [0.0, 0.0],1.0f0, 0.0f0) == 2.0
    @test meanSemiVariance([1.0, 1.0], [0.0, 0.0],-1.0f0, 0.0f0) === nothing

    @test meanSemiDevi([1.0, 1.0], [0.0, 0.0],1.0f0, 0.0f0, 2.0f0) == 2.0
    @test meanSemiDevi([1.0, 1.0], [0.0, 0.0],-1.0f0, 0.0f0, 2.0f0) === nothing
    @test meanSemiDevi([1.0, 1.0], [0.0, 0.0],1.0f0, 0.0f0, -2.0f0) === nothing

    @test spectral([0.0, 1.0], [0.2, 0.8], x -> 2.0*x) ≈ 0.96 atol = 0.0001
    @test spectral([0.0, 1.0], [0.2, 0.8], x -> x) === nothing

    @test distortion([0.0, 1.0], [1.0, 0.0], x -> x^2) == 0.0

    o1 = :( (Y)^2 *p); o2 = :sqrt; conds = [o1 :+ o2];
    states = [1.0, 2.0, 3.0]; prob = [0.3, 0.4, 0.3];
    @test GenCoherent(states, prob,conds)[1] ≈ 3.0 atol = 0.0001

    @test GenConvex(states, prob,conds,x->x) === nothing
end


text:
# Update move limits
function update_limits!(primal_data, m, k, x1, x2, s_init, s_incr, s_decr)
    @unpack L, x, U, = primal_data
    for j in 1:dim(m)
        if k == 1 || k == 2
            # Equation 11 in Svanberg
            L[j] = x[j] - s_init * (max(m,j) - min(m, j))
            U[j] = x[j] + s_init * (max(m,j) - min(m, j))
        else
            # Equation 12 in Svanberg
            if sign(x[j] - x1[j]) != sign(x1[j] - x2[j])
                L[j] = x[j] - (x1[j] - L[j]) * s_decr
                U[j] = x[j] + (U[j] - x1[j]) * s_decr
            # Equation 13 in Svanberg
            else
                L[j] = x[j] - (x1[j] - L[j]) * s_incr
                U[j] = x[j] + (U[j] - x1[j]) * s_incr
            end
        end
    end
end

function compute_mma!(primal_data, m)
    @unpack x, L, U, α, β, f_x, g, ∇f_x, ∇g, p, q, p0, q0, r = primal_data
    # Bound limits
    for j = 1:dim(m)
        μ = 0.1
        α[j] = max(L[j] + μ * (x[j] - L[j]), min(m, j))
        β[j] = min(U[j] - μ * (U[j] - x[j]), max(m, j))
    end

    r0 = 0.0
    for i in 0:length(constraints(m))
        if i == 0
            ri = f_x
            ∇fi = @view ∇f_x[:]
        else
             ri = g[i]
             ∇fi = @view ∇g[:,i]
        end
        for j in 1:dim(m)
            Ujxj = U[j] - x[j]
            xjLj = x[j] - L[j]
            if ∇fi[j] > 0
                p_ij = abs2(Ujxj) * ∇fi[j]
                q_ij = 0.0
            else
                p_ij = 0.0
                q_ij = -abs2(xjLj) * ∇fi[j]
            end
            ri -= p_ij / Ujxj + q_ij / xjLj
            if i == 0
                p0[j] = p_ij
                q0[j] = q_ij
            else
                p[j, i] = p_ij
                q[j, i] = q_ij
            end
        end
        if i == 0
            r0 = ri
        else
            r[i] = ri
        end
    end
    primal_data.r0 = r0
end

function compute_dual!(λ, primal_data)
    @unpack p0, q0, p, q, L, U, x, r0, r = primal_data
    update_x!(primal_data, λ)

    #Optimal value of Lagrangian at λ
    φ = r0 + dot(λ, r)
    @inbounds for j = 1:length(x)
        φ += (p0[j] + matdot(λ, p, j)) / (U[j] - x[j])
        φ += (q0[j] + matdot(λ, q, j)) / (x[j] - L[j])
    end
    return -φ
end

function compute_dual_grad!(∇φ, λ, primal_data)
    @unpack r, p, q, L, U, x = primal_data
    update_x!(primal_data, λ)
    for i = 1:length(λ)
        ∇φ[i] = r[i]
        for j = 1:length(x)
            ∇φ[i] += p[j,i] / (U[j] - x[j])
            ∇φ[i] += q[j,i] / (x[j] - L[j])
        end
    end
    # Negate since we have a maximization problem
    scale!(∇φ, -1.0)
    return
end

# Updates x to be the analytical optimal point in the dual
# problem for a given λ
function update_x!(primal_data, λ)
    @unpack x, p0, q0, p, q, L, U, α, β = primal_data
    @inbounds for j in 1:length(x)
        fpj = sqrt(p0[j] + matdot(λ, p, j))
        fqj = sqrt(q0[j] + matdot(λ, q, j))
        x[j] = (fpj * L[j] + fqj * U[j]) / (fpj + fqj)
        if x[j] > β[j]
            x[j] = β[j]
        elseif x[j] < α[j]
            x[j] = α[j]
        end
    end
    return
end


text:
module ForwardFootball

export leer_BasicInfo, leer_TrainInfo

include("leer_BasicInfo.jl")
include("leer_TrainInfo.jl")

end # module

max_stars_repo_name: vntr-CC/Forward-Football

max_issues_repo_name: vntr-CC/Forward-Football

max_forks_repo_name: vntr-CC/Forward-Football


text:
function LindbladJumpCallback()
    r1 = Ref{Float64}(rand(Float64))
    r2 = Ref{Float64}(rand(Float64))
    condition = (u, t, integrator) -> real(u' * u) - r1[]

    affect! = function (integrator)
        psi = copy(integrator.u)
        c_ops = copy(integrator.p)
        collaps_idx = 1
        if length(c_ops) == 1
            integrator.u = normalize(c_ops[1] * psi)
        else
            r2[] = rand()
            dp = 0
            for i in 1:length(c_ops)
                c_op = c_ops[i]
                dp += real(psi' * (c_op' * c_op) * psi)
            end
            prob = 0
            for i in 1:length(c_ops)
                c_op = c_ops[i]
                prob += real(psi' * (c_op' * c_op) * psi) / dp
                # println(prob)
                if prob >= r2[]
                    collaps_idx = i
                    # println(collaps_idx)
                    break
                end
            end
            integrator.u = normalize(c_ops[collaps_idx] * psi)
        end

        r1[] = rand()
    end

    initialize = function (c, u, t, integrator)
        r1[] = rand()
        u_modified!(integrator, false)
    end

    return ContinuousCallback(condition, affect!, initialize = initialize, save_positions = (false,false))
end

"""
    mcsolve(H::AbstractArray, psi0, t_l, c_ops; e_ops = [], n_traj = 100, ensemble_method = EnsembleSerial())

Time evolution of an open quantum system using quantum trajectories.
"""
function mcsolve(H::AbstractArray, psi0, t_l, c_ops; e_ops = [], n_traj = 100, ensemble_method = EnsembleSerial())
    N_c_ops = length(c_ops)
    N_t = length(t_l)
    ti, tf = t_l[1], t_l[end]
    tspan = (ti, tf)
    p = c_ops

    H_eff = H
    for i in 1:N_c_ops
        c_op = c_ops[i]
        H_eff += - 0.5im * c_op' * c_op
    end

    A = DiffEqArrayOperator(-1im * H_eff)
    prob = ODEProblem(A, psi0, tspan, p)
    cb1 = LindbladJumpCallback()
    cb2 = AutoAbstol(false; init_curmax=0.0)
    cb = CallbackSet(cb1, cb2)

    prob = ODEProblem{true}(A, psi0, tspan, p, callback = cb, abstol = 1e-7, reltol = 1e-5)
    ensemble_prob = EnsembleProblem(prob)
    # Tsit5() alg_hints=[:nonstiff] lsoda()
    sol = solve(ensemble_prob, alg_hints=[:nonstiff], ensemble_method, trajectories=n_traj, saveat = t_l)
    # sol = solve(ensemble_prob, LinearExponential(krylov=:adaptive),
    #     EnsembleThreads(), trajectories=n_traj, dt = (tf - ti) / N_t)

    if length(e_ops) != 0
        e_ops_expect = zeros(Float64, length(t_l), length(e_ops))
    end

    Threads.@threads for i in 1:length(t_l)
        for idx in 1:length(sol)
            normalize!(sol[idx].u[i])
        end
        for j in 1:length(e_ops)
            e_ops_expect[i, j] = mean([expect(e_ops[j], sol[idx2].u[i]) for idx2 in 1:length(sol)])
        end
    end

    # Threads.@threads for idx in 1:length(sol)
    #     for i in 1:length(sol[idx].t)
    #         normalize!(sol[idx].u[i])
    #     end
    # end

    if length(e_ops) == 0
        return sol
    else
        return sol, e_ops_expect
    end
end

function mesolve(H::AbstractArray, psi0, t_l, c_ops; e_ops = [], krylovdim = 10)
    N_c_ops = length(c_ops)
    N_t = length(t_l)
    ti, tf = t_l[1], t_l[end]
    tspan = (ti, tf)
    p = c_ops

    L = -1im * (spre(H) - spost(H))
    for i in 1:N_c_ops
        L += lindblad_dissipator( c_ops[i] )
    end

    rho0_vec = reshape(psi0 * psi0', length(psi0)^2)

    A = DiffEqArrayOperator(L)
    prob = ODEProblem(A, rho0_vec, tspan, abstol = 1e-7, reltol = 1e-5)
    sol = solve(prob, LinearExponential(krylov=:adaptive, m = krylovdim), dt = (tf - ti) / (N_t - 1))

    if length(e_ops) != 0
        e_ops_expect = zeros(Float64, length(t_l), length(e_ops))

        Threads.@threads for i in 1:length(t_l)
            for j in 1:length(e_ops)
                if ishermitian(e_ops[j])
                    e_ops_expect[i, j] = real(tr(e_ops[j] * reshape(sol.u[i], size(H, 1), size(H, 1))))
                else
                    e_ops_expect[i, j] = tr(e_ops[j] * reshape(sol.u[i], size(H, 1), size(H, 1)))
                end
            end
        end

        return sol, e_ops_expect
    else
        return sol
    end
end

function sesolve(H::AbstractArray, psi0, t_l; e_ops = [], krylovdim = 10)
    N_t = length(t_l)
    ti, tf = t_l[1], t_l[end]
    tspan = (ti, tf)

    A = DiffEqArrayOperator(-1im * H)
    prob = ODEProblem(A, psi0, tspan, abstol = 1e-7, reltol = 1e-5)
    sol = solve(prob, LinearExponential(krylov=:adaptive, m = krylovdim), dt = (tf - ti) / (N_t - 1))

    if length(e_ops) != 0
        e_ops_expect = zeros(Float64, length(t_l), length(e_ops))

        Threads.@threads for i in 1:length(t_l)
            for j in 1:length(e_ops)
                e_ops_expect[i, j] = expect(e_ops[j], sol.u[i])
            end
        end

        return sol, e_ops_expect
    else
        return sol
    end
end


text:
using .tensorboard: TextPluginData

"""
    log_text(logger::TBLogger, name::String, text::Any; step=step(logger))

Logs text with name `name` at step `step`
- text: If text is a 2-D or 3-D `Array`, it will be rendered as a table or a list. Any other data will be represented as string
"""
function log_text(logger::TBLogger, name::String, text::Any; step = nothing)
    summ = SummaryCollection(text_summary(name, text))
    write_event(logger.file, make_event(logger, summ, step=step))
end

"""
    construct_tensor_dims(data) -> TensorShapeProto_Dim

Constructs the `TensorShapeProto_Dim` protobuffer holding the size (and eventually
name) of a single dimension of a tensor/array.

If you define a custom array type, such as AxisArray or named array types,
you should define this method so that you can serialize the name to string.
"""
function construct_tensor_dims(data::AbstractArray)
    dims = Vector{TensorShapeProto_Dim}()
    for dim in size(data)
        push!(dims, TensorShapeProto_Dim(size = dim))
    end
    return dims
end

construct_tensor_dims(data) = [TensorShapeProto_Dim(size=1)]

function construct_tensor_string_val(data::AbstractArray)
    textstringval = Vector{Vector{UInt8}}()

    if ndims(data) == 2
        text = permutedims(data, (2, 1))
    end

    for string in data
        string = markdown_repr(string)
        push!(textstringval, Vector{UInt8}(string))
    end

    return textstringval
end

construct_tensor_string_val(data) = [Vector{UInt8}(markdown_repr(data))]

function text_summary(name::String, text::Any)
    # Create a string tensor
    textstringval = construct_tensor_string_val(text)

    # Structure holding the shape of the tensor
    dims = construct_tensor_dims(text)
    texttensorshape = TensorShapeProto(dim = dims)

    # Metadata for the text
    textcontent = serialize_proto(TextPluginData(version = 0))
    plugindata = SummaryMetadata_PluginData(plugin_name = "text", content = textcontent)
    smd = SummaryMetadata(plugin_data = plugindata)

    # Create the tensor
    texttensor = TensorProto(dtype = _DataType.DT_STRING, string_val = textstringval, tensor_shape = texttensorshape)
    Summary_Value(tag = name, metadata = smd, tensor = texttensor)
end

"""
    markdown_repr(x)

Returns the most suitable method for displaying `x` in a markdown document.
"""
function markdown_repr(x)
    for mime_type in ("text/markdown", "text/html", "text/plain")
        try
            if showable(mime_type, x)
                return repr(mime_type, x)
            end
        catch
            @debug "It said it was showable, but it lied" mime_type x exception=err
            # we blindly catch everything, because who knows what evil
            # people are doing in their misimplemented `repr` functions
        end
    end
    return repr(x)
end


text:
using ChebyshevQuantum
using Test

tol_val = 1e-11

@testset "test C interp" begin

    function f(x)
        return exp(-0.2*sin(3*x))
    end

    c = ChebyshevQuantum.Chebyshev.make_Cheb(f , a=2, b = 3)

    @test isapprox(c(2.3), f(2.3), atol=tol_val)

    function f2(x)
        if x < 0.1
            return exp(-0.2*sin(3*x))
        else
            return 1+x^2
        end
    end

    c = ChebyshevQuantum.Chebyshev.make_Cheb(f2 , a=-1, b = 1)

    x=0.0
    @test isapprox(c(x), f2(x), atol=tol_val)
    x=0.2
    @test isapprox(c(x), f2(x), atol=tol_val)

end


text:
## shortcut
typealias Rat Rational{Int}

"a type used to enter note lengths, e.g. 1//4 == quarter note length"
IR = Union{Int,Rat}

## midi conversion helper functions

"float [0,1] to 14bit midi int conversion"
ftoi14(val::Float64) = clamp(convert(Int,round(val*16383)),0,16383)

"float [0,1] to 7bit midi int conversion"
ftoi7(val::Float64) = clamp(convert(Int,round(val*127)),0,127)

"midi int 7bit to float [0,1] conversion"
itof14(val::Int) = ( val = clamp(val,0,16383) ; val/16383. )

"midi int 14bit to float [0,1] conversion"
itof7(val::Int) = ( val = clamp(val,0,127) ; val/127. )

"""
Abstract atomic audio expression. \\
Complex audio expressions (`Exp`) contain a list of these.
"""
abstract Atom

"when non-zero duration is required"
abstract Duratom <: Atom

"""
A complex audio expression. \\
It contains zero or more `Atom`(s). \\
It has its own duration seperate of its `Atom` members.
"""
immutable Exp
	dur::Rat			# duration
	as::Vector{Atom}	# order not important, atoms contain thier own offset values
end

Exp() = Exp(0//1,Atom[])

Base.convert(::Type{Exp}, a::Atom) = Exp(a.dur,Atom[a])

"""
	atomOrExp(z)

Tries to convert input to an `Atom`. \\
If that fails, tries to convert to `Exp`.
"""
atomOrExp(a::Atom) = a
atomOrExp(x::Exp) = x

function atomOrExp(z)
	try
		return convert(Atom,z)
	catch
		return convert(Exp,z)
	end
end

"""
	offsets(x::Exp)

Returns `(t1,t2)` where, \\
`t1` is the start of the earliest `Atom` \\
`t2` is the end of the latest `Atom`
"""
function offsets(x::Exp)
	t1 = foldr((a,v)->min(a.ofs,v),1//0,x.as)
	t2 = foldr((a,v)->max(a.ofs+a.dur,v),-1//0,x.as)
	t1, t2
end

## adds x to mt after its offset
## returns x's length in ticks for mt
## see SirenSeq.Midi for MidiTrack
toTrack!(mt::MidiTrack, x::Exp) = ( t1 = mt.ff ; for a in x.as ; t1 = max(toTrack!(mt,a),t1) ; end ; t1 )

## used by C(zs...) and atomsOp
atomsPush!(as::Vector{Atom}, a::Atom) = ( push!(as,a) ; a.dur )
atomsPush!(as::Vector{Atom}, x::Exp) = ( for a in x.as ; atomsPush!(as,a) ; end ; x.dur )
atomsPush!(as::Vector{Atom}, z) = atomsPush!(as,atomOrExp(z))

"""
	C(zs...)

Creates a non-atomic audio expression (`Exp`) with members of `zs` all played at the same time.
"""
function C(zs...)
	as = Atom[]
	dur = 0//1
	for z in zs
		dd = atomsPush!(as,z)
		dur = max(dur,dd)
	end
	Exp(dur,as)
end

atomsOp(op::Function, x::Exp) = ( as = Atom[] ; for a in x.as ; atomsPush!(as,op(a)) ; end ; as )

## return copy of `x` with midi channel of all Atoms changed to `v`
function channel(v::Int, x::Exp)
	as = atomsOp(z->channel(v,z),x)
	Exp(x.dur,as)
end

## return copy of `x` with all Atom durations and offsets multiplied by `v`
function dilate(v::Rational{Int}, x::Exp)
	as = atomsOp(z->dilate(v,z),x)
	Exp(x.dur*v,as)
end

## return copy of `x` with `v` added to all Atom offsets
function sshift(v::Rational{Int}, x::Exp)
 	as = atomsOp(z->sshift(v,z),x)
	Exp(x.dur,as)
end

## return copy of `x` with all Atom velocities mutiplied by `v`
function accel(v::Float64, x::Exp)
	as = atomsOp(z->accel(v,z),x)
	Exp(x.dur,as)
end

## return copy of `x` with `v` added to all note inteval values
function transl(v::Int, x::Exp)
	as = atomsOp(z->transl(v,z),x)
	Exp(x.dur,as)
end

## return copy of `x` with `v` added to all note octave values
function octave(v::Int, x::Exp)
	as = atomsOp(z->octave(v,z),x)
	Exp(x.dur,as)
end

## return copy of `x` with all note scales changed to `v`
function sscale(v::Scales.Scale, x::Exp)
	as = atomsOp(z->sscale(v,z),x)
	Exp(x.dur,as)
end

## used by S(zs...) and R(v,z)
function seqPush!(as::Vector{Atom}, a::Atom, dur::Rational{Int})
	a2 = sshift(dur,a)
	push!(as,a2)
	dur+a2.dur
end

function seqPush!(as::Vector{Atom}, x::Exp, dur::Rational{Int})
	for a in x.as
		a2 = sshift(dur,a)
		push!(as,a2)
	end
	dur+x.dur
end

seqPush!(as::Vector{Atom}, z, dur::Rational{Int}) = seqPush!(as,atomOrExp(z),dur)

"""
	S(zs...)

Creates a non-atomic audio expression (`Exp`) with members of `zs` all played in sequence.
"""
function S(zs...)
	as = Atom[]
	dur = 0//1
	for z in zs
		dur = seqPush!(as,z,dur)
	end
	Exp(dur,as)
end

"""
	R(v::Int, z)

Returns `z` repeated `v` times in sequence.
"""
function R(v::Int, z)
	@assert v > 0
	as = Atom[]
	dur = 0//1
	for i in 1:v
		dur = seqPush!(as,z,dur)
	end
	Exp(dur,as)
end

function Base.isequal(x::Exp, y::Exp)
	if x.dur != y.dur ; return false ; end
	as1 = sort!(x.as,lt=atomHashOrder)
	as2 = sort!(y.as,lt=atomHashOrder)
	if length(as1) != length(as2) ; return false ; end
	for i in 1:length(as1)
		if !isequal(as1[i],as2[i]) ; return false ; end
	end
	true
end

function Base.isequal(x::Atom, y::Exp)
	if x.dur != y.dur ; return false ; end
	if length(y.as) != 1 ; return false ; end
	isequal(x,t.as[1])
end

function atomHashOrder(x1::Atom, x2::Atom)
	if hash(typeof(x1)) < hash(typeof(x1)) ; return true ; end
	@assert typeof(x1) == typeof(x2)
	for t in fieldnames(typeof(x1))
		v1 = getfield(x1,t)
		v2 = getfield(x2,t)
		if v1 == v2 ; continue ; end
		return hash(v2) < hash(v1)
	end
	return false
end


text:
"""
    BrennanSchwartz()

[Brennan and Schwartz model](https://link.springer.com/chapter/10.1057/9781403946027_3).
"""
struct BrennanSchwartz <: Model end

"""
    evaluate(IRD, BrennanSchwartz(), n = 12)

Evaluate interest rate derivative `IRD` using `BrennanSchwartz` model.

# Arguments
- `IRD::InterestRateDerivative`: interest rate derivative
- `n`: number of paths to simulate
"""
function evaluate(IRD::InterestRateDerivative, m::BrennanSchwartz, n::Int64 = 12)
    Δt = IRD.t / n
    rates = [IRD.r]
    for i = 0:n
        dr = IRD.k * (IRD.θ - rates[end]) * Δt + IRD.σ * randn()
        append!(rates, rates[end] + dr)
    end
    return rates
end


text:
@testset "make_sin" begin
    n = 50
    xtmp, ytmp = make_sin(n; noise = 0.)

    @test length(xtmp) == length(ytmp) == n
    for i = 1:length(xtmp)
        @test sin.(xtmp[i]) ≈ ytmp[i]
    end
    # print(scatterplot(xtmp, ytmp; color = :blue, height = 5))
end

@testset "make_poly" begin
    coef = [.8, .5, 2]
    xtmp, ytmp = make_poly(coef, -10:.1:10; noise = 0)

    @test length(xtmp) == length(ytmp)
    for i = 1:length(xtmp)
        @test (coef[1] * xtmp[i]^2 + coef[2] * xtmp[i]^1 + coef[3]) ≈ ytmp[i]
    end
    # print(scatterplot(xtmp, ytmp; color = :blue, height = 5))
end

@testset "make_spiral" begin
    n = 97
    xtmp, ytmp = make_spiral(n; noise = 0.)

    @test length(xtmp[1, :]) == length(ytmp) == 2*n
    # test_plot = scatterplot(xtmp[1, 1:97], xtmp[2, 1:97], title="Spiral Function", color=:blue, name="pos")
    # print(scatterplot!(test_plot, xtmp[1, 98:194], xtmp[2, 98:194], color=:yellow, name="neg" ))
end


text:
module BSpline

export
  bspline_basismatrix!, bspline_basismatrix,
  derivative_bspline_basismatrix!, derivative_bspline_basismatrix,
  derivative2_bspline_basismatrix!, derivative2_bspline_basismatrix,
  bspline_basis, derivative_bspline_basis, derivative2_bspline_basis,
  num_basis,
  BSplineT

type BSplineT{T<:FloatingPoint}
  M::Int64                  # B-spline order (2 for linear, 3 for quadratic, etc.)
  knots::Vector{T}          # knot vector
end

BSplineT{T<:FloatingPoint}(M::Int64, knots::Vector{T}) = BSplineT{T}(M, knots)
# BSplineT{T<:FloatingPoint}(M::Int64, knots::LinSpace{T}) = BSplineT{T}(M, [knots;])

function BSplineT{T<:FloatingPoint}(M::Int64, lb::T, ub::T, num_knots::Int64)
  h = (ub - lb) / (num_knots - 1)

  BSplineT{T}(M,
             vcat(lb .- h * [M-1:-1:1;],
                  linspace(lb, ub, num_knots),
                  ub .+ h * [1:M-1;]))
end

num_basis{T<:FloatingPoint}(bs::BSplineT{T}) = length(bs.knots) - bs.M

# B-spline basis function value B(j,n) at x.
#
# Input arguments:
# j:   interval index, 0 =< j < num_basis
#
# x:   value where the basis function is to be evaluated
#
# Output arguments:
# y:
#    B-spline basis function value, nonzero for a knot span of n

function bspline_basis{T<:FloatingPoint}(j::Integer, bs::BSplineT{T}, x::T)
  n = bs.M
  t = bs.knots
  y = zero(T)
  if n > 1
    dn = x - t[j+1]
    dd = t[j+n] - t[j+1]
    if dd != zero(T)  # indeterminate forms 0/0 are deemed to be zero
      b = bspline_basis(j, BSplineT(n-1,t), x)
      y = y + b*(dn/dd)
    end

    dn = t[j+n+1] - x
    dd = t[j+n+1] - t[j+2]
    if dd != zero(T)
      b = bspline_basis(j+1, BSplineT(n-1,t), x)
      y = y + b*(dn/dd)
    end
  else
    y = t[j+1] <= x < t[j+2] ? one(T) : zero(T)
  end
  y
end

# B-spline basis function value matrix B(n) for x.
#
# Input arguments:
# n:
#    B-spline order (2 for linear, 3 for quadratic, etc.)
# t:
#    knot vector
# x :
#    an m-dimensional vector of values where the basis function is to be
#    evaluated
#
# Output arguments:
# out:
#    a matrix of m rows and numel(t)-n columns

function bspline_basismatrix!{T<:FloatingPoint}(
    out::StridedMatrix{T}, bs::BSplineT{T}, x::Vector{T}
    )
  for j = 1 : num_basis(bs)
    for k = 1 : length(x)
      out[k, j] = bspline_basis(j-1, bs, x[k])
    end
  end
  out
end

bspline_basismatrix{T<:FloatingPoint}(bs::BSplineT{T}, x::T) =
  bspline_basismatrix!(Array(T, (1, num_basis(bs))), bs, [x;])

bspline_basismatrix{T<:FloatingPoint}(bs::BSplineT{T}, x::Vector{T}) =
  bspline_basismatrix!(Array(T, (length(x), num_basis(bs))), bs, x)

# derivatives

function derivative_bspline_basis{T<:FloatingPoint}(j::Integer, bs::BSplineT{T}, x::T)
  n = bs.M
  t = bs.knots

  y = zero(T)
  if n > 1
    dn = n - 1
    dd = t[j+n] - t[j+1]
    if dd != zero(T)  # indeterminate forms 0/0 are deemed to be zero
      b = bspline_basis(j, BSplineT(n-1, t), x)
      y = y + b*(dn/dd)
    end

    dn = 1 - n
    dd = t[j+n+1] - t[j+2]
    if dd != zero(T)
      b = bspline_basis(j+1, BSplineT(n-1, t), x)
      y = y + b*(dn/dd)
    end
  else
    y = t[j+1] <= x < t[j+2] ? one(T) : zero(T)
  end
  y
end

function derivative_bspline_basismatrix!{T<:FloatingPoint}(
        out::StridedMatrix{T}, bs::BSplineT{T}, x::Vector{T}
        )
  for j = 1 : num_basis(bs)
    for k = 1 : length(x)
      out[k,j] = derivative_bspline_basis(j-1,bs,x[k])
    end
  end
  out
end

derivative_bspline_basismatrix{T<:FloatingPoint}(bs::BSplineT{T}, x::Vector{T}) =
  derivative_bspline_basismatrix!(Array(T, (length(x), num_basis(bs))), bs, x)

derivative_bspline_basismatrix{T<:FloatingPoint}(bs::BSplineT{T}, x::T) =
  derivative_bspline_basismatrix!(Array(T, (1, num_basis(bs))), bs, [x;])

# second derivative

function derivative2_bspline_basis{T<:FloatingPoint}(j::Integer, bs::BSplineT{T}, x::T)
  n = bs.M
  t = bs.knots
  y = zero(T)
  if n > 1
    dn = n - 1
    dd = t[j+n] - t[j+1]
    if dd != zero(T)  # indeterminate forms 0/0 are deemed to be zero
      b = derivative_bspline_basis(j, BSplineT(n-1, t), x)
      y = y + b*(dn/dd)
    end

    dn = 1 - n
    dd = t[j+n+1] - t[j+2]
    if dd != zero(T)
      b = derivative_bspline_basis(j+1, BSplineT(n-1, t), x)
      y = y + b*(dn/dd)
    end
  else
    y = t[j+1] <= x < t[j+2] ? one(T) : zero(T)
  end
  y
end

function derivative2_bspline_basismatrix!{T<:FloatingPoint}(
            out::StridedMatrix{T}, bs::BSplineT{T}, x::Vector{T}
            )
  for j = 1 : num_basis(bs)
    for k = 1 : length(x)
      out[k,j] = derivative2_bspline_basis(j-1, bs, x[k])
    end
  end
end

derivative_bspline_basismatrix{T<:FloatingPoint}(bs::BSplineT{T}, x::Vector{T}) =
  derivative2_bspline_basismatrix!(Array(T, (length(x), num_basis(bs))), bs, x)

derivative_bspline_basismatrix{T<:FloatingPoint}(bs::BSplineT{T}, x::T) =
  derivative2_bspline_basismatrix!(Array(T, (1, num_basis(bs))), bs, [x;])

end

max_stars_repo_name: mlakolar/npgm

max_issues_repo_name: mlakolar/npgm

max_forks_repo_name: mlakolar/npgm


text:
push!(LOAD_PATH, "./src/")

using AdaGram
using AdaGram.ArgParse
using Distributed

s = ArgParseSettings()

@add_arg_table s begin
  "model"
    help = "path to serialized model"
    arg_type = AbstractString
    required = true
  "text"
    help = "text to measure likelihood"
    arg_type = AbstractString
    required = true
  "--window"
    help = "window size"
    arg_type = Int
    default = 5
  "--workers"
    help = "number of workers"
    arg_type = Int
    default = 1
  "--minprob"
    help = "minimum probability of a prototype"
    arg_type = Float64
    default = 0.005
  "--batch"
    help = "size of buffer read into memory"
    arg_type = Int
    default = 16777216
  "--log"
    help = "save intermediate averages to the file"
    arg_type = AbstractString
end

args = parse_args(ARGS, s)

addprocs(args["workers"])
@everywhere using AdaGram

vm, dict = load_model(args["model"])
println(parallel_likelihood(vm, dict, args["text"], args["window"], args["minprob"];
  batch=args["batch"], log=args["log"]))


text:
using Test
using SASLib

using Dates
using Statistics: mean
using SharedArrays: SharedArray
using Tables

function getpath(dir, file)
    path = joinpath(dir, file)
    #println("================ $path ================")
    path
end
readfile(dir, file; kwargs...)  = readsas(getpath(dir, file); kwargs...)
openfile(dir, file; kwargs...)  = SASLib.open(getpath(dir, file), kwargs...)
getmetadata(dir, file; kwargs...) = metadata(getpath(dir, file), kwargs...)

# Struct used for column type conversion test case below
struct YearStr year::String end
Base.convert(::Type{YearStr}, v::Float64) = YearStr(string(round(Int, v)))

@testset "SASLib" begin

    @testset "object pool" begin
        println("Testing object pool...")

        # string pool
        default = ""
        x = SASLib.ObjectPool{String, UInt8}(default, 5)
        @test length(x) == 5
        @test size(x) == (5, )
        @test lastindex(x) == 5
        @test count(v -> v == default, x) == 5
        @test count(v -> v === default, x) == 5
        @test map(v -> "x$v", x) == [ "x", "x", "x", "x", "x" ]
        x[1] = "abc"
        @test x[1] == "abc"
        @test x.uniqueitemscount == 2
        x[2] = "abc"
        @test x[2] == "abc"
        @test x.uniqueitemscount == 2
        x[3] = "xyz"
        @test x.uniqueitemscount == 3
        @test x.itemscount == 5
        @test_throws BoundsError x[6] == ""

        # tuple pool
        y = SASLib.ObjectPool{Tuple, UInt8}((1,1,1), 100)
        y[1:100] = [(v, v, v) for v in 1:100]
        @test y[1] == (1,1,1)
        @test y[2] == (2,2,2)
        @test y[100] == (100,100,100)
        @test y.uniqueitemscount == 100   # first one is the same as the default

        # more error conditions
        z = SASLib.ObjectPool{Int, UInt8}(0, 1000)
        @test_throws BoundsError z[1:300] = 1:300
    end

    @testset "case insensitive dict" begin
        function testdict(lowercase_key, mixedcase_key, second_lowercase_key)

            T = typeof(lowercase_key)
            d = SASLib.CIDict{T,Int}()

            # getindex/setindex!
            d[lowercase_key] = 99
            @test d[lowercase_key] == 99
            @test d[mixedcase_key] == 99
            d[mixedcase_key] = 88           # should replace original value
            @test length(d) == 1            # still 1 element
            @test d[lowercase_key] == 88
            @test d[mixedcase_key] == 88

            # haskey
            @test haskey(d, lowercase_key) == true
            @test haskey(d, mixedcase_key) == true

            # iteration
            d[second_lowercase_key] = 77
            ks = T[]
            vs = Int[]
            for (k,v) in d
                push!(ks, k)
                push!(vs, v)
            end
            @test ks == [lowercase_key, second_lowercase_key]
            @test vs == [88, 77]

            # keys/values
            @test collect(keys(d)) == [lowercase_key, second_lowercase_key]
            @test collect(values(d)) == [88, 77]

            # show
            @test show(d) == nothing
        end
        testdict(:abc, :ABC, :def)
        testdict("abc", "ABC", "def")
    end

    @testset "open and close" begin
        handler = openfile("data_pandas", "test1.sas7bdat")
        @test typeof(handler) == SASLib.Handler
        @test handler.config.filename == getpath("data_pandas", "test1.sas7bdat")
        @test SASLib.close(handler) == nothing
    end

    @testset "read basic test files (test*.sas7bdat)" begin
        dir = "data_pandas"
        files = filter(x -> endswith(x, "sas7bdat") && startswith(x, "test"),
            Base.Filesystem.readdir("$dir"))
        for f in files
            result = readfile(dir, f)
            @test size(result) == (10, 100)
        end
    end

    @testset "incremental read" begin
        handler = openfile("data_pandas", "test1.sas7bdat")
        @test handler.config.filename == getpath("data_pandas", "test1.sas7bdat")
        result = SASLib.read(handler, 3)  # read 3 rows
        @test size(result, 1) == 3
        result = SASLib.read(handler, 4)  # read 4 rows
        @test size(result, 1) == 4
        result = SASLib.read(handler, 5)  # should read only 3 rows even though we ask for 5
        @test size(result, 1)  == 3
    end

    @testset "various data types" begin
        rs = readfile("data_pandas", "test1.sas7bdat")
        @test sum(rs[:Column1][1:5]) == 2.066
        @test count(isnan, rs[:Column1]) == 1
        @test rs[:Column98][1:3] == [ "apple", "dog", "pear" ]
        @test rs[:Column4][1:3] == [Date("1965-12-10"), Date("1977-03-07"), Date("1983-08-15")]
    end

    @testset "datetime with missing values" begin
        rs = readfile("data_pandas", "datetime.sas7bdat")
        @test size(rs) == (5, 4)
        @test rs[:mtg][1] == Date(2017, 11, 24)
        @test rs[:dt][5] == DateTime(2018, 3, 31, 14, 20, 33)
        @test count(ismissing, rs[:mtg]) == 1
        @test count(ismissing, rs[:dt]) == 3
    end

    @testset "include/exclude columns" begin
        fname = getpath("data_pandas", "productsales.sas7bdat")

        rs = readsas(fname, include_columns=[:MONTH, :YEAR])
        @test size(rs, 2) == 2
        @test sort(names(rs)) == sort([:MONTH, :YEAR])

        rs = readsas(fname, include_columns=[1, 2, 7])
        @test size(rs, 2) == 3
        @test sort(names(rs)) == sort([:ACTUAL, :PREDICT, :PRODUCT])

        rs = readsas(fname, exclude_columns=[:DIVISION])
        @test size(rs, 2) == 9
        @test !(:DIVISION in names(rs))

        rs = readsas(fname, exclude_columns=collect(2:10))
        @test size(rs, 2) == 1
        @test sort(names(rs)) == sort([:ACTUAL])

        # case insensitive include/exclude
        rs = readsas(fname, include_columns=[:month, :Year])
        @test size(rs, 2) == 2
        rs = readsas(fname, exclude_columns=[:diVisiON])
        @test size(rs, 2) == 9

        # test bad include/exclude param
        # see https://discourse.julialang.org/t/test-warn-doesnt-work-with-warn-in-0-7/9001
        @test_logs (:warn, "Unknown include column blah") (:warn,
            "Unknown include column Year") readsas(fname, include_columns=[:blah, :Year])
        @test_logs (:warn, "Unknown exclude column blah") (:warn,
            "Unknown exclude column Year") readsas(fname, exclude_columns=[:blah, :Year])
        # error handling
        @test_throws SASLib.ConfigError readsas(fname,
            include_columns=[1], exclude_columns=[1])
    end

    @testset "ResultSet" begin
        rs = readfile("data_pandas", "productsales.sas7bdat")

        # metadata for result set
        @test size(rs) == (1440, 10)
        @test size(rs,1) == 1440
        @test size(rs,2) == 10
        @test length(columns(rs)) == 10
        @test length(names(rs)) == 10

        # cell indexing
        @test rs[1][1] ≈ 925.0
        @test rs[1,1] ≈ 925.0
        @test rs[1,:ACTUAL] ≈ 925.0

        # row/column indexing
        @test typeof(rs[1]) == NamedTuple{
            (:ACTUAL, :PREDICT, :COUNTRY, :REGION, :DIVISION, :PRODTYPE, :PRODUCT, :QUARTER, :YEAR, :MONTH),
            Tuple{Float64,Float64,String,String,String,String,String,Float64,Float64,Dates.Date}}
        @test typeof(rs[:ACTUAL]) == Array{Float64,1}
        @test sum(rs[:ACTUAL]) ≈ 730337.0

        # iteration
        @test sum(r[1] for r in rs) ≈ 730337.0

        # portion of result set
        @test typeof(rs[1:2]) == SASLib.ResultSet
        @test typeof(rs[:ACTUAL, :PREDICT]) == SASLib.ResultSet
        @test rs[1:2][1][1] ≈ 925.0
        @test rs[:ACTUAL, :PREDICT][1][1] ≈ 925.0

        # setindex!
        rs[1,1] = 100.0
        @test rs[1,1] ≈ 100.0
        rs[1,:ACTUAL] = 200.0
        @test rs[1,:ACTUAL] ≈ 200.0

        # display related
        @test show(rs) == nothing
        @test SASLib.sizestr(rs) == "1440 rows x 10 columns"

        # Tables.jl interface - getproperty test
        @test rs.ACTUAL == rs[:ACTUAL]
        @test names(rs) == propertynames(rs)

        # Tables.jl coverage - indirect tests / usage
        @test Tables.schema(rs).names == Tuple(names(rs))
        @test Tables.schema(rs).types == Tuple(eltype.([rs[s] for s in names(rs)]))
        @test length(Tables.rowtable(rs)) == 1440
        @test length(Tables.columntable(rs)) == 10
        @test size(Tables.matrix(rs[:ACTUAL, :PREDICT])) == (1440,2)

        # Tables.jl coverage - direct tests
        @test Tables.istable(typeof(rs)) === true
        @test Tables.rowaccess(typeof(rs)) === true
        @test Tables.columnaccess(typeof(rs)) === true
        @test Tables.rows(rs) |> first |> propertynames |> Tuple == Tuple(names(rs))
        @test Tables.columns(rs) |> propertynames |> Tuple == Tuple(names(rs))
    end

    @testset "metadata" begin
        md = getmetadata("data_pandas", "test1.sas7bdat")
        @test md.filename == getpath("data_pandas", "test1.sas7bdat")
        @test md.encoding == "WINDOWS-1252"
        @test md.endianness == :LittleEndian
        @test md.compression == :none
        @test md.pagesize == 65536
        @test md.npages == 1
        @test md.nrows == 10
        @test md.ncols == 100
        @test length(md.columnsinfo) == 100
        @test md.columnsinfo[1] == Pair(:Column1, Float64)

        md = getmetadata("data_pandas", "productsales.sas7bdat")
        @test show(md) == nothing
        println()

        # convenient comparison routine since v0.6/v0.7 displays different order
        same(x,y) = sort(string.(collect(x))) == sort(string.(collect(y)))

        @test same(SASLib.typesof(Int64), (Int64,))
        @test same(SASLib.typesof(Union{Int64,Int32}), (Int64,Int32))
        @test same(SASLib.typesof(Union{Int64,Int32,Missing}),
            (Int64,Int32,Missing))

        @test SASLib.typesfmt((Int64,)) == "Int64"
        @test SASLib.typesfmt((Int64,Int32)) == "Int64/Int32"
        @test SASLib.typesfmt((Int64,Int32,Missing)) == "Int64/Int32/Missing"
        @test SASLib.typesfmt((Int64,Int32); excludemissing=true) == "Int64/Int32"
        @test SASLib.typesfmt((Int64,Int32,Missing); excludemissing=true) == "Int64/Int32"
        @test SASLib.colfmt(md)[1] == "ACTUAL(Float64)"
    end

    @testset "stat_transfer" begin
        rs = readfile("data_misc", "types.sas7bdat")
        @test sum(rs[:vbyte][1:2])   == 9
        @test sum(rs[:vint][1:2])    == 9
        @test sum(rs[:vlong][1:2])   == 9
        @test sum(rs[:vfloat][1:2])  ≈  10.14000010
        @test sum(rs[:vdouble][1:2]) ≈  10.14000000
    end

    # topical.sas7bdat contains columns labels which should be ignored anywas
    @testset "AHS2013" begin
        handler = openfile("data_AHS2013", "topical.sas7bdat")
        rs = SASLib.read(handler, 1000)
        @test size(rs) == (1000, 114)
        @test show(handler) == nothing
        SASLib.close(handler)
		# @test result[:page_count] == 10
        # @test result[:page_length] == 16384
        # @test result[:system_endianness] == :LittleEndian
        @test count(x -> x == "B", rs[:DPEVVEHIC]) == 648
        @test mean(filter(!isnan, rs[:PTCOSTGAS])) ≈ 255.51543209876544
    end

    @testset "file encodings" begin
        rs = readfile("data_reikoch", "extr.sas7bdat")
        # @test result[:file_encoding] == "CP932"
        @test rs[:AETXT][1] == "眠気"

        rs = readfile("data_pandas", "test1.sas7bdat", encoding = "US-ASCII")
        # @test result[:file_encoding] == "US-ASCII"
        @test rs[:Column42][3] == "dog"
    end

    @testset "handler object" begin
        handler = openfile("data_reikoch", "binary.sas7bdat")
        @test handler.U64 == true
        @test handler.byte_swap == true
        @test handler.column_data_lengths == [8,8,8,8,8,8,8,8,8,8,14]
        @test handler.column_data_offsets == [0,8,16,24,32,40,48,56,64,72,80]
        @test handler.column_names == ["I","I1","I2","I3","I4","I5","I6","I7","I8","I9","CHAR"]
        @test handler.column_symbols == [:I,:I1,:I2,:I3,:I4,:I5,:I6,:I7,:I8,:I9,:CHAR]
        @test handler.compression == 0x02
        @test handler.file_encoding == "ISO-8859-1"
        @test handler.file_endianness == :BigEndian
        @test handler.header_length == 8192
        @test handler.page_length == 8192
        @test handler.row_count == 100
        @test handler.vendor == 0x01
        @test handler.config.convert_dates == true
        @test handler.config.include_columns == []
        @test handler.config.exclude_columns == []
        @test handler.config.encoding == ""
    end

    @testset "array constructors" begin

        rs = readfile("data_AHS2013", "homimp.sas7bdat")
        @test typeof(rs[:RAS]) == SASLib.ObjectPool{String,UInt16}

        # string_array_fn test for specific string columns
        rs = readfile("data_AHS2013", "homimp.sas7bdat",
            string_array_fn = Dict(:RAS => REGULAR_STR_ARRAY))
        @test typeof(rs[:RAS]) == Array{String,1}
        @test typeof(rs[:RAH]) != Array{String,1}

        # string_array_fn test for all string columns
        rs = readfile("data_AHS2013", "homimp.sas7bdat",
            string_array_fn = Dict(:_all_ => REGULAR_STR_ARRAY))
        @test typeof(rs[:RAS])     == Array{String,1}
        @test typeof(rs[:RAH])     == Array{String,1}
        @test typeof(rs[:JRAS])    == Array{String,1}
        @test typeof(rs[:JRAD])    == Array{String,1}
        @test typeof(rs[:CONTROL]) == Array{String,1}

        # number_array_fn test by column name
        makesharedarray(n) = SharedArray{Float64}(n)
        rs = readfile("data_misc", "numeric_1000000_2.sas7bdat",
            number_array_fn = Dict(:f => makesharedarray))
        @test typeof(rs[:f]) == SharedArray{Float64,1}
        @test typeof(rs[:x]) == Array{Float64,1}

        # number_array_fn test for all numeric columns
        rs = readfile("data_misc", "numeric_1000000_2.sas7bdat",
        number_array_fn = Dict(:_all_ => makesharedarray))
        @test typeof(rs[:f]) == SharedArray{Float64,1}
        @test typeof(rs[:x]) == SharedArray{Float64,1}

    end

    # column type conversion
    @testset "user specified column types" begin

        # normal use case
        rs = readfile("data_pandas", "productsales.sas7bdat";
            verbose_level = 0, column_types = Dict(:YEAR => Int16, :QUARTER => Int8))
        @test eltype(rs[:YEAR]) == Int16
        @test eltype(rs[:QUARTER]) == Int8

        # error handling - warn() when a column cannot be converted
        rs = readfile("data_pandas", "productsales.sas7bdat";
            verbose_level = 0, column_types = Dict(:YEAR => Int8, :QUARTER => Int8))
        @test eltype(rs[:YEAR]) == Float64
        @test eltype(rs[:QUARTER]) == Int8
        #TODO expect warning for :YEAR conversion

        # case insensitive column symbol
        rs = readfile("data_pandas", "productsales.sas7bdat";
            verbose_level = 0, column_types = Dict(:Quarter => Int8))
        @test eltype(rs[:QUARTER]) == Int8

        # conversion to custom types
        rs = readfile("data_pandas", "productsales.sas7bdat";
            verbose_level = 0, column_types = Dict(:Year => YearStr))
        @test eltype(rs[:YEAR]) == YearStr

        # test Union type
        let T = Union{Int,Missing}
            rs = readfile("data_pandas", "productsales.sas7bdat";
                verbose_level = 0, column_types = Dict(:Year => T))
            @test eltype(rs[:YEAR]) == T
        end
    end

    # see output; keep this for coverage reason
    @testset "verbosity" begin
        rs = readfile("data_pandas", "test1.sas7bdat"; verbose_level = 2)
        @test size(rs, 1) > 0
    end

    @testset "just reads" begin
        for dir in ["data_pandas", "data_reikoch", "data_AHS2013", "data_misc"]
            for f in readdir(dir)
                if endswith(f, ".sas7bdat") &&
                        !(f in ["zero_variables.sas7bdat"])
                    rs = readfile(dir, f)
                    @test size(rs, 1) > 0
                end
            end
        end
    end

    @testset "exception" begin
        @test_throws SASLib.FileFormatError readsas("runtests.jl")
    end

end


text:
__precompile__()

module JuAFEM
using ForwardDiff
using Reexport
@reexport using Tensors
@reexport using WriteVTK

import Base: show, push!, @propagate_inbounds

include("exports.jl")

"""
Represents a reference shape which quadrature rules and interpolations are defined on.
Currently, the only concrete types that subtype this type are `RefCube` in 1, 2 and 3 dimensions,
and `RefTetrahedron` in 2 and 3 dimensions.
"""
abstract type AbstractRefShape end

struct RefTetrahedron <: AbstractRefShape end
struct RefCube <: AbstractRefShape end

"""
Abstract type which has `CellValues` and `FaceValues` as subtypes
"""
abstract type Values{dim,T,refshape} end
abstract type CellValues{dim,T,refshape} <: Values{dim,T,refshape} end
abstract type FaceValues{dim,T,refshape} <: Values{dim,T,refshape} end

include("utils.jl")

# Interpolations
include("interpolations.jl")

# Quadrature
include(joinpath("Quadrature", "quadrature.jl"))

# FEValues
include(joinpath("FEValues","cell_values.jl"))
include(joinpath("FEValues","face_values.jl"))
include(joinpath("FEValues","common_values.jl"))
include(joinpath("FEValues","face_integrals.jl"))

# Grid
include(joinpath("Grid", "grid.jl"))
include(joinpath("Grid", "grid_generators.jl"))

# Dofs
include(joinpath("Dofs", "DofHandler.jl"))
include(joinpath("Dofs", "ConstraintHandler.jl"))

include("iterators.jl")

# Assembly
include("assembler.jl")

# Export
include(joinpath("Export", "VTK.jl"))

# Other
include("deprecations.jl")

end # module


text:
using StructJuMP

numScens = 2
m = StructuredModel(num_scenarios=numScens)

@variable(m, x >= 1)
@variable(m, y <= 2)

@constraint(m,  x + y == 1)
@constraint(m, -x + y <= 1)
@objective(m, :Min, x*x + 0.5x*y + 0.25y*y - y)

rhs  = [5,4]
coef = [2,3]
qc   = [1,0.5]
ac   = [1,0.75]

for i = 1:numScens
    bl = StructuredModel(parent=m, id=i)
    @variable(bl, 0 <= w <= 1)
    @constraint(bl, coef[i]w - x - y <= rhs[i])
    @constraint(bl, coef[i]w + x     == rhs[i])
    @objective(bl, :Min, qc[i]*w*w + ac[i]*w)
end


text:
cd(@__DIR__); include("setups/grid7x3.jl");

## HAD
D = natural_eigdist(𝚽, 𝛌, Q; distance = :HAD)
E = transform(fit(MDS, D, maxoutdim=2, distances=true))
plt = grid7x3_mds_heatmaps(E, 𝚽; backend = :pyplot)
savefig(plt, "../figs/Grid7x3_MDS_HAD.png")
display(plt)