text:
module SUN_generator
    using LinearAlgebra

    struct Generator
        NC::Int64
        generator::Array{Array{ComplexF64,2},1}
        function Generator(NC)
            return new(NC,make_generators(NC))
        end
    end

    function Base.getindex(g::Generator,i)
        return g.generator[i]
    end

    function Base.length(a::Generator)
        return length(a.generator)
    end

    function lie2matrix!(matrix,g::Generator,a)
        matrix .= 0
        NC = g.NC
        for (i,genmatrix) in enumerate(g.generator)
            matrix .+= a[i]*genmatrix
        end
        return
    end

    function matrix2lie!(a,g::Generator,A)
        for i=1:length(a)
            #println("i = $i")
            #display(g.generator[i]*g.generator[i])
            #println("\t")
            #println(tr(g.generator[i]*A)/2)
            a[i] = tr(g.generator[i]*A)/2
        end
        return
    end

    sigma = []
    s = ComplexF64[
        0 1
        1 0
    ]
    push!(sigma,s)
    s = ComplexF64[
        0 -im
        im 0
    ]
    push!(sigma,s)
    s = ComplexF64[
        1 0
        0 -1
    ]
    push!(sigma,s)

    function make_largematrix(σ,i,j,Nc)
        λ = zeros(ComplexF64,Nc,Nc)
        λ[i,i] = σ[1,1]
        λ[i,j] = σ[1,2]
        λ[j,i] = σ[2,1]
        λ[j,j] = σ[2,2]
        return λ
    end

    function normalization(A)
        Norm = sqrt(real(tr(A*A))/2)
        #print("    Normalization = √",real(tr(A*A))/2,"= ")
        #println("$Norm")
        return Norm
    end

    function make_generators(Nc)
        d=["x","y","z"]
        #
        lams = Array{Complex,2}[]
        #
        if Nc < 2
            error("Invalid Nc=$Nc")
        end
        #
        # off-diagonal part
        #println("Off diagonals")
        for i=1:Nc
            for j=i+1:Nc
                for a = 1:2
                    #println("$(d[a]): i j = $i $j")
                    A = make_largematrix(sigma[a],i,j,Nc)
                    #showmat(A)
                    N=normalization(A)
                    A = A/N
                    push!(lams,A)
                    #println(" ")
                end
            end
        end
        #
        # diagonal part
        #println("diagonals")
        for a=1:Nc-1
            A = zeros(ComplexF64,Nc,Nc)
            for i=1:a
                A[i,i] = 1
            end
            A[a+1,a+1] = -tr(A)
            #showmat(A)
            N=normalization(A)
            A = A/N
            push!(lams,A)
            #println(" ")
        end
        return lams
    end
end


text:
export
    inference, generate,
    loss, loss_pn, loss_pln,
    loss_f, callback_f,
    ICNFModel, ICNFDist

function inference(icnf::AbstractICNF{T}, mode::TestMode, xs::AbstractMatrix{T}, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractVector where {T <: AbstractFloat} end
function inference(icnf::AbstractICNF{T}, mode::TrainMode, xs::AbstractMatrix{T}, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractVector where {T <: AbstractFloat} end

function generate(icnf::AbstractICNF{T}, mode::TestMode, n::Integer, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractMatrix{T} where {T <: AbstractFloat} end
function generate(icnf::AbstractICNF{T}, mode::TrainMode, n::Integer, p::AbstractVector=icnf.p; rng::Union{AbstractRNG, Nothing}=nothing)::AbstractMatrix{T} where {T <: AbstractFloat} end

function loss(icnf::AbstractICNF{T}, xs::AbstractMatrix{T}, p::AbstractVector=icnf.p; agg::Function=mean) where {T <: AbstractFloat} end

function loss_pn(icnf::AbstractICNF{T}, xs::AbstractMatrix{T}, p::AbstractVector=icnf.p; agg::Function=mean, nλ::T=convert(T, 1e-4)) where {T <: AbstractFloat}
    lv = loss(icnf, xs, p; agg)
    prm_n = norm(p)
    lv + nλ*prm_n
end

function loss_pln(icnf::AbstractICNF{T}, xs::AbstractMatrix{T}, p::AbstractVector=icnf.p; agg::Function=mean, nλ::T=convert(T, 1e-4)) where {T <: AbstractFloat}
    lv = loss(icnf, xs, p; agg)
    prm_ln = log(norm(p))
    lv + nλ*prm_ln
end

# -- Flux interface

function (icnf::AbstractICNF{T})(xs::AbstractMatrix{T})::AbstractVector{T} where {T <: AbstractFloat}
    inference(icnf, TestMode(), xs)
end

function loss_f(icnf::AbstractICNF{T}, opt_app::FluxOptApp)::Function where {T <: AbstractFloat}
    function f(xs::AbstractMatrix{T})::T
        loss(icnf, xs)
    end
    f
end

function callback_f(icnf::AbstractICNF{T}, opt_app::FluxOptApp, loss::Function, data::DataLoader{T3})::Function where {T <: AbstractFloat, T2 <: AbstractMatrix{T}, T3 <: Tuple{T2}}
    xs, = first(data)
    function f()::Nothing
        vl = loss(icnf, xs)
        @info "Training" loss=vl
        nothing
    end
    f
end

# -- Optim interface

function loss_f(icnf::AbstractICNF{T}, opt_app::OptimOptApp, itrtr::AbstractVector)::Function where {T <: AbstractFloat}
    function f(p::AbstractVector{T})::T
        xs, = itrtr[1]
        loss(icnf, xs, p)
    end
    f
end

function callback_f(icnf::AbstractICNF{T}, opt_app::OptimOptApp, loss::Function, data::DataLoader{T3}, itrtr::AbstractVector)::Function where {T <: AbstractFloat, T2 <: AbstractMatrix{T}, T3 <: Tuple{T2}}
    xs, = first(data)
    function f(s::OptimizationState)::Bool
        vl = loss(icnf, xs, s.metadata["x"])
        @info "Training" loss=vl
        nxitr = iterate(data, itrtr[2])
        if isnothing(nxitr)
            true
        else
            itrtr .= nxitr
            false
        end
    end
    f
end

# -- SciML interface

function loss_f(icnf::AbstractICNF{T}, opt_app::SciMLOptApp)::Function where {T <: AbstractFloat}
    function f(p::AbstractVector, θ::SciMLBase.NullParameters, xs::AbstractMatrix{T})
        loss(icnf, xs, p)
    end
    f
end

function callback_f(icnf::AbstractICNF{T}, opt_app::SciMLOptApp, loss::Function, data::DataLoader{T3})::Function where {T <: AbstractFloat, T2 <: AbstractMatrix{T}, T3 <: Tuple{T2}}
    xs, = first(data)
    function f(p::AbstractVector{T}, l::T)::Bool
        vl = loss(icnf, xs, p)
        @info "Training" loss=vl
        false
    end
    f
end

# MLJ interface

mutable struct ICNFModel{T, T2} <: MLJICNF where {T <: AbstractFloat, T2 <: AbstractICNF{T}}
    m::T2
    loss::Function

    opt_app::OptApp
    optimizer::Any
    n_epochs::Integer
    adtype::SciMLBase.AbstractADType

    batch_size::Integer
end

function ICNFModel(
        m::T2,
        loss::Function=loss,
        ;
        opt_app::OptApp=FluxOptApp(),
        optimizer::Any=default_optimizer[typeof(opt_app)],
        n_epochs::Integer=128,
        adtype::SciMLBase.AbstractADType=Optimization.AutoZygote(),

        batch_size::Integer=128,
        ) where {T <: AbstractFloat, T2 <: AbstractICNF{T}}
    ICNFModel{T, T2}(m, loss, opt_app, optimizer, n_epochs, adtype, batch_size)
end

function MLJModelInterface.fit(model::ICNFModel, verbosity, X)
    x = collect(MLJModelInterface.matrix(X)')
    data = DataLoader((x,); batchsize=model.batch_size, shuffle=true, partial=true)
    ncdata = ncycle(data, model.n_epochs)
    initial_loss_value = model.loss(model.m, x)

    if model.opt_app isa FluxOptApp
        model.optimizer isa Flux.Optimise.AbstractOptimiser || error("model.optimizer must be a Flux optimizer")
        _loss = loss_f(model.m, model.opt_app)
        _callback = callback_f(model.m, model.opt_app, model.loss, data)
        _p = Flux.params(model.m)
        tst = @timed Flux.Optimise.train!(_loss, _p, ncdata, model.optimizer; cb=_callback)
    elseif model.opt_app isa OptimOptApp
        model.optimizer isa Optim.AbstractOptimizer || error("model.optimizer must be an Optim optimizer")
        itrtr = Any[nothing, nothing]
        itrtr .= iterate(ncdata)
        _loss = loss_f(model.m, model.opt_app, itrtr)
        _callback = callback_f(model.m, model.opt_app, model.loss, data, itrtr)
        ops = Optim.Options(
            x_abstol=-Inf, x_reltol=-Inf,
            f_abstol=-Inf, f_reltol=-Inf,
            g_abstol=-Inf, g_reltol=-Inf,
            outer_x_abstol=-Inf, outer_x_reltol=-Inf,
            outer_f_abstol=-Inf, outer_f_reltol=-Inf,
            outer_g_abstol=-Inf, outer_g_reltol=-Inf,
            f_calls_limit=0, g_calls_limit=0, h_calls_limit=0,
            allow_f_increases=true, allow_outer_f_increases=true,
            successive_f_tol=typemax(Int), iterations=typemax(Int), outer_iterations=typemax(Int),
            store_trace=false, trace_simplex=true, show_trace=false, extended_trace=true,
            show_every=1, callback=_callback, time_limit=Inf,
        )
        tst = @timed res = optimize(_loss, model.m.p, model.optimizer, ops)
        model.m.p .= res.minimizer
    elseif model.opt_app isa SciMLOptApp
        _loss = loss_f(model.m, model.opt_app)
        _callback = callback_f(model.m, model.opt_app, model.loss, data)
        optfunc = OptimizationFunction(_loss, model.adtype)
        optprob = OptimizationProblem(optfunc, model.m.p)
        tst = @timed res = solve(optprob, model.optimizer, ncdata; callback=_callback)
        model.m.p .= res.u
    end
    final_loss_value = model.loss(model.m, x)
    @info("Fitting",
        "elapsed time (seconds)"=tst.time,
        "garbage collection time (seconds)"=tst.gctime,
    )

    fitresult = nothing
    cache = nothing
    report = (
        stats=tst,
        initial_loss_value=initial_loss_value,
        final_loss_value=final_loss_value,
    )
    fitresult, cache, report
end

function MLJModelInterface.transform(model::ICNFModel, fitresult, Xnew)
    xnew = collect(MLJModelInterface.matrix(Xnew)')

    tst = @timed logp̂x = inference(model.m, TestMode(), xnew)
    @info("Transforming",
        "elapsed time (seconds)"=tst.time,
        "garbage collection time (seconds)"=tst.gctime,
    )

    DataFrame(px=exp.(logp̂x))
end

function MLJModelInterface.fitted_params(model::ICNFModel, fitresult)
    (
        learned_parameters=model.m.p,
    )
end

MLJBase.metadata_pkg(
    ICNFModel,
    package_name="ICNF",
    package_uuid="9bd0f7d2-bd29-441d-bcde-0d11364d2762",
    package_url="https://github.com/impICNF/ICNF.jl",
    is_pure_julia=true,
    package_license="MIT",
    is_wrapper=false,
)
MLJBase.metadata_model(
    ICNFModel,
    input_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
    target_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
    output_scitype=Table{AbstractVector{ScientificTypes.Continuous}},
    supports_weights=false,
    docstring="ICNFModel",
    load_path="ICNF.ICNFModel",
)

# Distributions interface

struct ICNFDist{T, T2} <: ICNFDistribution where {T <: AbstractFloat, T2 <: AbstractICNF{T}}
    m::T2
end

function ICNFDist(m::T2) where {T <: AbstractFloat, T2 <: AbstractICNF{T}}
    ICNFDist{T, T2}(m)
end

Base.length(d::ICNFDist) = d.m.nvars
Base.eltype(d::ICNFDist) = eltype(d.m.p)
Distributions._logpdf(d::ICNFDist, x::AbstractVector) = first(Distributions._logpdf(d, hcat(x)))
Distributions._logpdf(d::ICNFDist, A::AbstractMatrix) = inference(d.m, TestMode(), A)
Distributions._rand!(rng::AbstractRNG, d::ICNFDist, x::AbstractVector) = (x[:] = Distributions._rand!(rng, d, hcat(x)))
Distributions._rand!(rng::AbstractRNG, d::ICNFDist, A::AbstractMatrix) = (A[:] = generate(d.m, TestMode(), size(A, 2); rng))

ext: jl


text:
include("dynamics.jl")
using Plots
using ReferenceFrameRotations
#dynamics.jlのテスト
#Runge-Kutta法のテスト

len=63000
a=[zeros(len) zeros(len)]
a[1,2]=1.0
dt=0.01
#理論解
b=[sin.(range(0,(len-1)*dt,step=dt)) cos.(range(0,(len-1)*dt,step=dt))]
#微分する関数
function dif(v::Array{Float64,1},A)
    θ=atan(v[1],v[2])
    cos(θ),-sin(θ)
end
#積分
for  i in 2:len
    a[i,:]=RK4((a[i-1,:]),0,dt)
    a[i,:]./=norm(a[i,:])
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
COMMENT=r"^//.*"
SECTION=r"^//-+\s[a-zA-Z]*\s-*"
iscomment(s) = occursin(COMMENT, s)
issection(s) = iscomment(s) && occursin(SECTION, s)

function partition(p, s)
  idxs = findall(p, s)
  prepend!(idxs, 1)
  push!(idxs, length(s))
  ranges = [idxs[(i-1)]:idxs[i] for i in 2:length(idxs)]
  return [s[r] for r in ranges]
end


text:
# Move all NaNs to the end of the array A
function sortnans!(A, iₗ=firstindex(A), iᵤ=lastindex(A))
    # Count up NaNs
    Nₙₐₙ = 0
    @turbo for i = iₗ:iᵤ
        Nₙₐₙ += A[i] != A[i]
    end
    # If none, return early
    Nₙₐₙ == 0 && return A, iₗ, iᵤ

    # Otherwise, swap all NaNs
    i = iₗ
    j = iᵤ
    N = iᵤ - iₗ
    @inbounds for n = 0:N-Nₙₐₙ
        i = iₗ + n
        if A[i] != A[i]
            while A[j] != A[j]
                j -= 1
            end
            j <= i && break
            A[i], A[j] = A[j], A[i]
            j -= 1
        end
    end
    return A, iₗ, iᵤ - Nₙₐₙ
end
# For integers, don't need to check for NaNs
sortnans!(A::AbstractArray{<:Integer}, iₗ=firstindex(A), iᵤ=lastindex(A)) = A, iₗ, iᵤ

# Check for sortedness, assuming no NaNs
@inline function issortedrange(A, iₗ, iᵤ)
    @inbounds for i = iₗ+1:iᵤ
        if A[i-1] > A[i]
            return false
        end
    end
    return true
end

# Check for anti-sortedness, assuming no NaNs
@inline function isantisortedrange(A, iₗ, iᵤ)
    @inbounds for i = iₗ+1:iᵤ
        if A[i-1] < A[i]
            return false
        end
    end
    return true
end

# Reverse an array, faster than Base.reverse!
@inline function vreverse!(A, iₗ, iᵤ)
    N = (iᵤ - iₗ) + 1
    n = (N ÷ 2) - 1
    if N < 32
        @inbounds for i ∈ 0:n
            𝔦ₗ, 𝔦ᵤ = iₗ+i, iᵤ-i
            A[𝔦ₗ], A[𝔦ᵤ] = A[𝔦ᵤ], A[𝔦ₗ]
        end
    else
        @turbo for i ∈ 0:n
            𝔦ₗ = iₗ+i
            𝔦ᵤ = iᵤ-i
            l = A[𝔦ₗ]
            u = A[𝔦ᵤ]
            A[𝔦ₗ] = u
            A[𝔦ᵤ] = l
        end
    end
    return A
end

# Partially sort `A` around the `k`th sorted element and return that element
function quickselect!(A::AbstractArray, iₗ=firstindex(A), iᵤ=lastindex(A), k=(iₗ+iᵤ)÷2)
    # Pick a pivot for partitioning
    N = iᵤ - iₗ + 1
    A[iₗ], A[k] = A[k], A[iₗ]
    pivot = A[iₗ]

    # Count up elements that must be moved to upper partition
    Nᵤ = 0
    @turbo for i = (iₗ+1):iᵤ
        Nᵤ += A[i] >= pivot
    end
    Nₗ = N - Nᵤ

    # Swap elements between upper and lower partitions
    i = iₗ
    j = iᵤ
    @inbounds for n = 1:Nₗ-1
        i = iₗ + n
        if A[i] >= pivot
            while A[j] >= pivot
                j -= 1
            end
            j <= i && break
            A[i], A[j] = A[j], A[i]
            j -= 1
        end
    end
    # Move pivot to the top of the lower partition
    iₚ = iₗ + Nₗ - 1
    A[iₗ], A[iₚ] = A[iₚ], A[iₗ]
    # Recurse: select from partition containing k
    (iₗ <= k < iₚ) && quickselect!(A, iₗ, iₚ, k)
    (iₚ < k <= iᵤ) && quickselect!(A, iₚ+1, iᵤ, k)
    return A[k]
end

# Sort `A`, assuming no NaNs
function quicksort!(A, iₗ=firstindex(A), iᵤ=lastindex(A))
    if issortedrange(A, iₗ, iᵤ)
        # If already sorted, we're done here
        return A
    end
    # Otherwise, we have to sort
    N = iᵤ - iₗ + 1
    if isantisortedrange(A, iₗ, iᵤ)
        vreverse!(A, iₗ, iᵤ)
        return A
    elseif N == 3
        # We know we are neither sorted nor antisorted, so only four possibilities remain
        iₘ = iₗ + 1
        a,b,c = A[iₗ], A[iₘ], A[iᵤ]
        if a <= b
            if a <= c
                A[iₘ], A[iᵤ] = c, b             # a ≤ c ≤ b
            else
                A[iₗ], A[iₘ], A[iᵤ] = c, a, b   # c ≤ a ≤ b
            end
        else
            if a <= c
                A[iₗ], A[iₘ] = b, a             # b ≤ a ≤ c
            else
                A[iₗ], A[iₘ], A[iᵤ] = b, c, a   # b ≤ c ≤ a
            end
        end
        return A
    else
        # Pick a pivot for partitioning
        iₚ = iₗ + (N >> 2)
        A[iₗ], A[iₚ] = A[iₚ], A[iₗ]
        pivot = A[iₗ]

        # Count up elements that must be moved to upper partition
        Nᵤ = 0
        @turbo for i = (iₗ+1):iᵤ
            Nᵤ += A[i] >= pivot
        end
        Nₗ = N - Nᵤ

        # Swap elements between upper and lower partitions
        i = iₗ
        j = iᵤ
        @inbounds for n = 1:Nₗ-1
            i = iₗ + n
            if A[i] >= pivot
                while A[j] >= pivot
                    j -= 1
                end
                j <= i && break
                A[i], A[j] = A[j], A[i]
                j -= 1
            end
        end
        # Move pivot to the top of the lower partition
        iₚ = iₗ + Nₗ - 1
        A[iₗ], A[iₚ] = A[iₚ], A[iₗ]
        # Recurse: sort both upper and lower partitions
        quicksort!(A, iₗ, iₚ)
        quicksort!(A, iₚ+1, iᵤ)
    end
end

# Sort `A`, assuming no NaNs, multithreaded
function quicksortt!(A, iₗ=firstindex(A), iᵤ=lastindex(A), level=1)
    if issortedrange(A, iₗ, iᵤ)
        # If already sorted, we're done here
        return A
    end
    # Otherwise, we have to sort
    N = iᵤ - iₗ + 1
    if isantisortedrange(A, iₗ, iᵤ)
        vreverse!(A, iₗ, iᵤ)
        return A
    elseif N == 3
        # We know we are neither sorted nor antisorted, so only four possibilities remain
        iₘ = iₗ + 1
        a,b,c = A[iₗ], A[iₘ], A[iᵤ]
        if a <= b
            if a <= c
                A[iₘ], A[iᵤ] = c, b             # a ≤ c ≤ b
            else
                A[iₗ], A[iₘ], A[iᵤ] = c, a, b   # c ≤ a ≤ b
            end
        else
            if a <= c
                A[iₗ], A[iₘ] = b, a             # b ≤ a ≤ c
            else
                A[iₗ], A[iₘ], A[iᵤ] = b, c, a   # b ≤ c ≤ a
            end
        end
        return A
    else
        # Pick a pivot for partitioning
        iₚ = iₗ + (N >> 2)
        A[iₗ], A[iₚ] = A[iₚ], A[iₗ]
        pivot = A[iₗ]

        # Count up elements that must be moved to upper partition
        Nᵤ = 0
        @turbo for i = (iₗ+1):iᵤ
            Nᵤ += A[i] >= pivot
        end
        Nₗ = N - Nᵤ

        # Swap elements between upper and lower partitions
        i = iₗ
        j = iᵤ
        @inbounds for n = 1:Nₗ-1
            i = iₗ + n
            if A[i] >= pivot
                while A[j] >= pivot
                    j -= 1
                end
                j <= i && break
                A[i], A[j] = A[j], A[i]
                j -= 1
            end
        end
        # Move pivot to the top of the lower partition
        iₚ = iₗ + Nₗ - 1
        A[iₗ], A[iₚ] = A[iₚ], A[iₗ]
        # Recurse: sort both upper and lower partitions
        if level < 7
            @sync begin
                Threads.@spawn quicksortt!(A, iₗ, iₚ, level+1)
                Threads.@spawn quicksortt!(A, iₚ+1, iᵤ, level+1)
            end
        else
            quicksort!(A, iₗ, iₚ)
            quicksort!(A, iₚ+1, iᵤ)
        end
        return A
    end
end


text:
abstract type Node end

mutable struct ScalarNode <: Node
    tag::AbstractString
    value::AbstractString
    start_mark::Union{Mark, Nothing}
    end_mark::Union{Mark, Nothing}
    style::Union{Char, Nothing}
end

mutable struct SequenceNode <: Node
    tag::AbstractString
    value::Vector
    start_mark::Union{Mark, Nothing}
    end_mark::Union{Mark, Nothing}
    flow_style::Bool
end

mutable struct MappingNode <: Node
    tag::AbstractString
    value::Vector
    start_mark::Union{Mark, Nothing}
    end_mark::Union{Mark, Nothing}
    flow_style::Bool
end


text:
module ClimaCommsSA

using ClimaComms
import ClimaComms: Neighbor

using Distributed
using KernelAbstractions
using SharedArrays

const ctx_chnl = Ref(RemoteChannel())
get_ctx_chnl() = ctx_chnl[]
const nbr_chnls = Dict{Int, RemoteChannel}()
get_nbr_chnl(nbr) = nbr_chnls[nbr]

mutable struct CCBarrier
    sense::SharedArray{Bool, 1}
    worker_sense::SharedArray{Bool, 1}
    function CCBarrier()
        sense = SharedArray{Bool}((1,), init = (sa) -> sa[1] = false)
        worker_sense =
            SharedArray{Bool}((nprocs(),), init = (sa) -> fill!(sa, false))
        new(sense, worker_sense)
    end
end

mutable struct CCReducer{FT}
    sense::SharedArray{Bool, 1}
    worker_sense::SharedArray{Bool, 1}
    vals::SharedArray{FT, 1}
    function CCReducer{FT}() where {FT}
        sense = SharedArray{Bool}((1,), init = (sa) -> sa[1] = false)
        worker_sense =
            SharedArray{Bool}((nprocs(),), init = (sa) -> fill!(sa, false))
        vals = SharedArray{FT}((nprocs(),), init = (sa) -> fill!(sa, zero(FT)))
        new(sense, worker_sense, vals)
    end
end

struct SACommsContext <: ClimaComms.AbstractCommsContext
    neighbors::Vector{ClimaComms.Neighbor}
    barrier::CCBarrier
    reducer_f64::CCReducer
    reducer_f32::CCReducer

    function SACommsContext(neighbors::Vector{ClimaComms.Neighbor})
        if ClimaComms.iamroot(ClimaCommsSA.SACommsContext)
            barr = CCBarrier()
            red64 = CCReducer{Float64}()
            red32 = CCReducer{Float32}()
            if nprocs() > 1
                for i in 1:nprocs()
                    put!(ClimaCommsSA.get_ctx_chnl(), (barr, red64, red32))
                end
            end
        else
            (barr, red64, red32) = take!(ClimaCommsSA.get_ctx_chnl())
        end
        new(neighbors, barr, red64, red32)
    end
end

function Neighbor(
    ::Type{SACommsContext},
    pid,
    AT,
    FT,
    send_dims,
    recv_dims = send_dims,
)
    @assert AT <: Array # until support for CuArray is added
    if pid <= 0
        return ClimaComms.Neighbor{Missing}(pid, missing, missing)
    end

    if ClimaComms.mypid(SACommsContext) < pid
        nbr_chnls[pid] = RemoteChannel()
        send_buf = SharedArray{FT}(
            send_dims,
            init = (sa) -> fill!(sa, zero(FT)),
            pids = [ClimaComms.mypid(SACommsContext), pid],
        )
        recv_buf = SharedArray{FT}(
            recv_dims,
            init = (sa) -> fill!(sa, zero(FT)),
            pids = [ClimaComms.mypid(SACommsContext), pid],
        )
        put!(ClimaCommsSA.get_nbr_chnl(pid), (send_buf, recv_buf))
    else
        while true
            try
                nbr_chnls[pid] = remotecall_fetch(
                    ClimaCommsSA.get_nbr_chnl,
                    pid,
                    ClimaComms.mypid(SACommsContext),
                )
                break
            catch e
                if e isa RemoteException && e.captured.ex isa KeyError
                    yield()
                    continue
                else
                    rethrow()
                end
            end
        end
        (send_buf, recv_buf) = take!(ClimaCommsSA.get_nbr_chnl(pid))
    end
    return ClimaComms.Neighbor{SharedArray}(pid, send_buf, recv_buf)
end

ClimaComms.get_stage(sa::SharedArray) = sdata(sa)

function ClimaComms.init(::Type{SACommsContext})
    ctx_chnl[] = remotecall_fetch(ClimaCommsSA.get_ctx_chnl, 1)
    return myid(), nprocs()
end

ClimaComms.mypid(::Type{SACommsContext}) = Distributed.myid()
ClimaComms.mypid(ctx::SACommsContext) = Distributed.myid()
ClimaComms.iamroot(CC::Type{SACommsContext}) = ClimaComms.mypid(CC) == 1
ClimaComms.iamroot(ctx::SACommsContext) = ClimaComms.mypid(ctx) == 1
ClimaComms.nprocs(::Type{SACommsContext}) = Distributed.nprocs()
ClimaComms.nprocs(::SACommsContext) = Distributed.nprocs()
ClimaComms.singlebuffered(::Type{SACommsContext}) = true
ClimaComms.neighbors(ctx::SACommsContext) = ctx.neighbors

function ClimaComms.start(ctx::SACommsContext; dependencies = nothing) end

function ClimaComms.progress(ctx::SACommsContext) end

function ClimaComms.finish(ctx::SACommsContext; dependencies = nothing) end

function ClimaComms.barrier(ctx::SACommsContext)
    if nprocs() == 1
        return
    end

    barr = ctx.barrier
    if ClimaComms.iamroot(ctx)
        for pid in workers()
            while barr.worker_sense[pid] == barr.sense[1]
                yield()
            end
        end
        barr.sense[1] = !barr.sense[1]
    else
        barr.worker_sense[ClimaComms.mypid(ctx)] =
            !barr.worker_sense[ClimaComms.mypid(ctx)]
        while barr.sense[1] != barr.worker_sense[ClimaComms.mypid(ctx)]
            yield()
        end
    end
end

function ClimaComms.reduce(ctx::SACommsContext, val, op)
    if nprocs() == 1
        return val
    end

    reducer = if val isa Float64
        ctx.reducer_f64
    elseif val isa Float32
        ctx.reducer_f32
    else
        error("No reducer for type $(typeof(val))")
    end

    if ClimaComms.iamroot(ctx)
        for pid in workers()
            while reducer.worker_sense[pid] == reducer.sense[1]
                yield()
            end
            val = op(val, reducer.vals[pid])
        end
        reducer.sense[1] = !reducer.sense[1]
        return val
    else
        reducer.vals[ClimaComms.mypid(ctx)] = val
        reducer.worker_sense[ClimaComms.mypid(ctx)] =
            !reducer.worker_sense[ClimaComms.mypid(ctx)]
        while reducer.sense[1] != reducer.worker_sense[ClimaComms.mypid(ctx)]
            yield()
        end
    end
end

ClimaComms.abort(ctx::SACommsContext, status::Int) = exit(status)

end # module


text:
module StrideArrays

# Write your package code here.

using VectorizationBase, ArrayInterface,
    SLEEFPirates, VectorizedRNG,
    LoopVectorization, LinearAlgebra,
    Random#, StackPointers#,
    # SpecialFunctions # Perhaps there is a better way to support erf?

using VectorizationBase: align, gep, AbstractStridedPointer, AbstractSIMDVector, vnoaliasstore!, staticm1,
    static_sizeof, lazymul, vmul_fast, StridedPointer, gesp, zero_offsets, pause, zstridedpointer,
    val_dense_dims, val_stride_rank, preserve_buffer
using LoopVectorization: maybestaticsize, CloseOpen
using Static: StaticInt, Zero, One, StaticBool, True, False
using ArrayInterface: OptionallyStaticUnitRange, size, strides, offsets, indices,
    static_length, static_first, static_last, axes,
    dense_dims, stride_rank
using StrideArraysCore: AbstractStrideArray, AbstractStrideMatrix, AbstractStrideVector,
    AbstractPtrStrideArray, PtrArray, static_expr, rank_to_sortperm,
    StrideArray, StrideVector, StrideMatrix, similar_layout,
    @gc_preserve

using Octavian
using Octavian: MemoryBuffer

export @StrideArray, @gc_preserve, # @Constant,
    AbstractStrideArray, AbstractStrideVector, AbstractStrideMatrix,
    StrideArray, StrideVector, StrideMatrix,
    PtrArray,# PtrVector, PtrMatrix,
    # ConstantArray, ConstantVector, ConstantMatrix, allocarray,
    matmul!, matmul_serial!, mul!, *ˡ, StaticInt,
    matmul, matmul_serial
# LazyMap,

include("rand.jl")
include("blas.jl")
include("broadcast.jl")
include("miscellaneous.jl")

# include("precompile.jl")
# _precompile_()

end


text:
#=
This script creates the Artifact.toml that will lazily provide the proj-datumgrid resources.
For more information see https://proj.org/resource_files.html

The resources will be downloaded from GitHub releases from this URL:
    https://github.com/OSGeo/proj-datumgrid/releases
With as a fallback option the OSGeo website:
    https://download.osgeo.org/proj/

To update to newer versions of one of the proj-datumgrid bundles, update the entry in the
datumgrids vector below, updating the version and sha256_hash. Then run this script again,
move the newly generated Artifact.toml over the one in the root of this package, and make
a pull request with the updated Artifact.toml.

Note that as of writing, this script will download about 1 GB of data.
=#

using Pkg.Artifacts
using BinaryProvider

# this file will be either created, or appended to
artifacts_toml = joinpath(@__DIR__, "Artifacts.toml")

datumgrids = [
    (name="proj-datumgrid", version=v"1.8", sha256_hash="3ff6618a0acc9f0b9b4f6a62e7ff0f7bf538fb4f74de47ad04da1317408fcc15"),
    (name="proj-datumgrid-europe", version=v"1.6", sha256_hash="af2f4f364d84eb9d5ca23403f8540ae8daee54ec900903bc92d48afb1ed285f8"),
    (name="proj-datumgrid-north-america", version=v"1.4", sha256_hash="95a6cdcdc078caed412ea66d3fdd2f699286336d5ca6c98bffb2d7132585ce40"),
    (name="proj-datumgrid-oceania", version=v"1.2", sha256_hash="952c43886c9ed6a098c308de61161a00cbe83fbb0f6e9e43c0c25c323bb06a4c"),
    (name="proj-datumgrid-world", version=v"1.0", sha256_hash="a488a5a69d1af6ec2ee83a2c64c52adac52e6dbfafe883f0341340009a9f40ba"),
]

for datumgrid in datumgrids
    # get the elements we need to build the correct URL
    name = datumgrid.name
    version = string(datumgrid.version.major, '.', datumgrid.version.minor)
    if name == "proj-datumgrid"
        tag = version
    else
        region = replace(name, "proj-datumgrid-"=>"")
        tag = string(region, '-', version)
    end
    sha = datumgrid.sha256_hash
    url1 = "https://github.com/OSGeo/proj-datumgrid/releases/download/$tag/$name-$version.tar.gz"
    url2 = "https://download.osgeo.org/proj/$name-$version.tar.gz"
    download_info = [(url1, sha), (url2, sha)]

    # download the package and save the hash
    proj_hash = create_artifact() do artifact_dir
        download_verify_unpack(url1, sha, artifact_dir, verbose=true, ignore_existence=true)
    end

    # add the new entry to the Artifact.toml together with the download info
    bind_artifact!(artifacts_toml, name, proj_hash; download_info=download_info, lazy=true, force=true)

    @info "created artifact" name version proj_hash artifacts_toml
end


text:
struct Equation <: AbstractEquation
  R_0::Number
  B_0::Number
  I_P::Number

  G_T::AbstractSymbol
end


text:
# Julia wrapper for header: /opt/local/include/gdal_rat.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0


text:
struct SynthWorkspace{T<:AF}
    lwavgrid::AA{T,1}
    rwavgrid::AA{T,1}
    allwavs::AA{T,1}
    allints::AA{T,1}
    wavt::AA{T,1}
    bist::AA{T,1}
    dept::AA{T,1}
    widt::AA{T,1}
end

function SynthWorkspace(; ndepths::Integer=100)
    lwavgrid = zeros(ndepths)
    rwavgrid = zeros(ndepths)
    allwavs  = zeros(2 * ndepths)
    allints  = zeros(2 * ndepths)
    wavt     = zeros(ndepths)
    bist     = zeros(ndepths)
    dept     = zeros(ndepths)
    widt     = zeros(ndepths)
    return SynthWorkspace(lwavgrid, rwavgrid, allwavs, allints, wavt, bist, dept, widt)
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
########################################################################################################################
module SlurmParallelSubmitIEEE
########################################################################################################################
# imports and exports
using FilterExps, SmootherExps, EnsembleKalmanSchemes, DeSolvers, L96, JLD2, Debugger

########################################################################################################################
########################################################################################################################
## Time series data
########################################################################################################################
# observation time series to load into the experiment as truth twin
# time series are named by the model, seed to initialize, the integration scheme used to produce, number of analyses,
# the spinup length, and the time length between observation points

ts01 = "../data/time_series/IEEE39bus_time_series_seed_0000_diff_0.012_tanl_0.01_nanl_50000_spin_5000_h_0.010.jld2"
########################################################################################################################

########################################################################################################################
## Experiment parameter generation
########################################################################################################################
########################################################################################################################
# Filters
########################################################################################################################
#
#
########################################################################################################################
# Classic smoothers
########################################################################################################################
## arguments are
## time_series, method, seed, lag, shift, obs_un, obs_dim, γ, param_err, param_wlk, N_ens, state_infl, param_infl = args
#
## incomplete will determine if the script checks for exisiting data before submitting experiments
#incomplete = true
#
## note, nanl is hard coded in the experiment, and h is inferred from the time series
## data, this is only used for checking versus existing data with the incomplete parameter as above
#nanl = 25000
#h = 0.01
#tanl = 0.01
#diffusion = 0.012
#
## these values set parameters for the experiment when running from scratch, or can
## be tested versus existing data
#sys_dim = 40
#obs_dim = 20
#methods = ["etks"]
#seed = 0
#
## note MDA is only defined for shifts / lags where the lag is a multiple of shift
## MDA is never defined for the classic smoother, but we will use the same parameter
## discretizations for SDA for reference values
#lags = 1:3:52
#shifts = [1]
#
## observation parameters, gamma controls nonlinearity
##gammas = Array{Float64}(1:11)
#gammas = [1.0]
#obs_un = 0.1
#param_err = 0.03
#param_infl = 1.0
#param_wlk = round.(LinRange(0, 10, 11) * 0.001, digits=3)
#
## if varying nonlinearity in obs, typically take a single ensemble value
##N_ens = [21]
#N_ens = 11:2:41
#
## inflation values, finite size versions should only be 1.0 generally
##state_infl = [1.0]
#state_infl = LinRange(1.00, 1.10, 11)
#
## set the time series of observations for the truth-twin
#time_series = [ts01]
#
## time_series, method, seed, lag, shift, obs_un, obs_dim, γ, param_err, param_wlk, N_ens, state_infl, param_infl = args
## load the experiments as a tuple
#args = Tuple[]
#for ts in time_series
#    for method in methods
#        for γ in gammas
#            for l in 1:length(lags)
#                lag = lags[l]
#                #shifts = lags[1:l]
#                for shift in shifts
#                    for N in N_ens
#                        for s_infl in state_infl
#                            for wlk in param_wlk
#                                if incomplete == true
#                                    #tanl = parse(Float64,ts[71:74])
#                                    #diffusion = parse(Float64,ts[60:64])
#                                    name = method *
#                                                "-classic_IEEE39bus_param_seed_" * lpad(seed, 4, "0") *
#                                                "_diff_" * rpad(diffusion, 5, "0") *
#                                                "_sysD_" * lpad(sys_dim, 2, "0") *
#                                                "_obsD_" * lpad(obs_dim, 2, "0") *
#                                                "_obsU_" * rpad(obs_un, 4, "0") *
#                                                "_gamma_" * lpad(γ, 5, "0") *
#                                                "_paramE_" * rpad(param_err, 4, "0") *
#                                                "_paramW_" * rpad(wlk, 6, "0") *
#                                                "_nanl_" * lpad(nanl, 5, "0") *
#                                                "_tanl_" * rpad(tanl, 4, "0") *
#                                                "_h_" * rpad(h, 4, "0") *
#                                                "_lag_" * lpad(lag, 3, "0") *
#                                                "_shift_" * lpad(shift, 3, "0") *
#                                                "_mda_false" *
#                                                "_nens_" * lpad(N, 3,"0") *
#                                                "_stateInfl_" * rpad(round(s_infl, digits=2), 4, "0") *
#                                                "_paramInfl_" * rpad(round(param_infl, digits=2), 4, "0") *
#                                                ".jld2"
#
#                                    fpath = "/x/capa/scratch/cgrudzien/power_grid_data/" * method * "-classic/"
#                                    try
#                                        f = load(fpath*name)
#                                    catch
#                                        tmp = (ts, method, seed, lag, shift, obs_un, obs_dim, γ, param_err, wlk, N, s_infl, param_infl)
#                                        push!(args, tmp)
#                                    end
#                                else
#                                    tmp = (ts, method, seed, lag, shift, obs_un, obs_dim, γ, param_err, wlk, N, s_infl, param_infl)
#                                    push!(args, tmp)
#                                end
#                            end
#                        end
#                    end
#                end
#            end
#        end
#    end
#end
#
#name = "/home/cgrudzien/DataAssimilationBenchmarks/data/input_data/classic_param_smoother_input_args.jld2"
#save(name, "experiments", args)
#
#for j in 1:length(args)
#    f = open("./submit_job.sl", "w")
#    write(f,"#!/bin/bash\n")
#    write(f,"#SBATCH -n 1\n")
#    # slow partition is for Okapi, uncomment when necessary
#    #write(f,"#SBATCH -p slow\n")
#    write(f,"#SBATCH -o ensemble_run.out\n")
#    write(f,"#SBATCH -e ensemble_run.err\n")
#    write(f,"julia SlurmExperimentDriver.jl " * "\"" *string(j) * "\"" * " \"classic_smoother_param\"")
#    close(f)
#    my_command = `sbatch  submit_job.sl`
#    run(my_command)
#end
#
########################################################################################################################
# Single-iteration smoothers
########################################################################################################################
## arguments are
## time_series, method, seed, lag, shift, mda, obs_un, obs_dim, γ,
## param_err, param_wlk, N_ens, state_infl, param_infl = args
#
## incomplete will determine if the script checks for exisiting data before submitting experiments
#incomplete = false
#
## note, nanl is hard coded in the experiment, and h is inferred from the time series
## data, this is only used for checking versus existing data with the incomplete parameter as above
#nanl = 12000
#h = 0.01
#tanl = 0.01
#diffusion = 0.012
#
## these values set parameters for the experiment when running from scratch, or can
## be tested versus existing data
#sys_dim = 40
#obs_dim = 20
#methods = ["etks"]
#seed = 0
#
## note MDA is only defined for shifts / lags where the lag is a multiple of shift
## MDA is never defined for the classic smoother, but we will use the same parameter
## discretizations for SDA for reference values
#lags = 1:3:52
#shifts = [1]
#
## observation parameters, gamma controls nonlinearity
##gammas = Array{Float64}(1:11)
#gammas = [1.0]
#obs_un = 0.1
#param_err = 0.03
#param_infl = 1.0
##param_wlk = round.(LinRange(0, 10, 11) * 0.001, digits=3)
#param_wlk = [0.0000]
#
## if varying nonlinearity in obs, typically take a single ensemble value
##N_ens = [21]
#N_ens = 11:2:41
#
## inflation values, finite size versions should only be 1.0 generally
##state_infl = [1.0]
#state_infl = LinRange(1.00, 1.10, 11)
#
## set the time series of observations for the truth-twin
#time_series = [ts01]
#
## set MDA true or false
#mdas = [false, true]
#
## time_series, method, seed, lag, shift, mda, obs_un, obs_dim, γ, param_err, param_wlk, N_ens, state_infl, param_infl = args
## load the experiments as a tuple
#args = Tuple[]
#for ts in time_series
#    for mda in mdas
#        for method in methods
#            for γ in gammas
#                for l in 1:length(lags)
#                    lag = lags[l]
#                    #shifts = lags[1:l]
#                    for shift in shifts
#                        for N in N_ens
#                            for s_infl in state_infl
#                                for wlk in param_wlk
#                                    if incomplete == true
#                                        #tanl = parse(Float64,ts[71:74])
#                                        #diffusion = parse(Float64,ts[60:64])
#                                        name = method *
#                                                    "-single-iteration_IEEE39bus_param_seed_" * lpad(seed, 4, "0") *
#                                                    "_diff_" * rpad(diffusion, 5, "0") *
#                                                    "_sysD_" * lpad(sys_dim, 2, "0") *
#                                                    "_obsD_" * lpad(obs_dim, 2, "0") *
#                                                    "_obsU_" * rpad(obs_un, 4, "0") *
#                                                    "_gamma_" * lpad(γ, 5, "0") *
#                                                    "_paramE_" * rpad(param_err, 4, "0") *
#                                                    "_paramW_" * rpad(wlk, 6, "0") *
#                                                    "_nanl_" * lpad(nanl, 5, "0") *
#                                                    "_tanl_" * rpad(tanl, 4, "0") *
#                                                    "_h_" * rpad(h, 4, "0") *
#                                                    "_lag_" * lpad(lag, 3, "0") *
#                                                    "_shift_" * lpad(shift, 3, "0") *
#                                                    "_mda_" * string(mda) *
#                                                    "_nens_" * lpad(N, 3,"0") *
#                                                    "_stateInfl_" * rpad(round(s_infl, digits=2), 4, "0") *
#                                                    "_paramInfl_" * rpad(round(param_infl, digits=2), 4, "0") *
#                                                    ".jld2"
#
#                                        fpath = "/x/capa/scratch/cgrudzien/power_grid_data/" * method * "-single-iteration/"
#                                        try
#                                            f = load(fpath*name)
#                                        catch
#                                            tmp = (ts, method, seed, lag, shift, mda, obs_un, obs_dim, γ, param_err, wlk, N, s_infl, param_infl)
#                                            push!(args, tmp)
#                                        end
#                                    else
#                                        tmp = (ts, method, seed, lag, shift, mda, obs_un, obs_dim, γ, param_err, wlk, N, s_infl, param_infl)
#                                        push!(args, tmp)
#                                    end
#                                end
#                            end
#                        end
#                    end
#                end
#            end
#        end
#    end
#end
#
#name = "/home/cgrudzien/DataAssimilationBenchmarks/data/input_data/single_iteration_param_smoother_input_args.jld2"
#save(name, "experiments", args)
#
#for j in 1:length(args)
#    f = open("./submit_job.sl", "w")
#    write(f,"#!/bin/bash\n")
#    write(f,"#SBATCH -n 1\n")
#    # slow partition is for Okapi, uncomment when necessary
#    #write(f,"#SBATCH -p slow\n")
#    write(f,"#SBATCH -o ensemble_run.out\n")
#    write(f,"#SBATCH -e ensemble_run.err\n")
#    write(f,"julia SlurmExperimentDriver.jl " * "\"" *string(j) * "\"" * " \"single_iteration_smoother_param\"")
#    close(f)
#    my_command = `sbatch  submit_job.sl`
#    run(my_command)
#end
#
#
########################################################################################################################
# Iterative smoothers
########################################################################################################################
# arguments are
# time_series, method, seed, lag, shift, mda, obs_un, obs_dim, γ,
# param_err, param_wlk, N_ens, state_infl, param_infl = args

# incomplete will determine if the script checks for exisiting data before submitting experiments
incomplete = false

# note, nanl is hard coded in the experiment, and h is inferred from the time series
# data, this is only used for checking versus existing data with the incomplete parameter as above
nanl = 25000
h = 0.01
tanl = 0.01
diffusion = 0.012

# these values set parameters for the experiment when running from scratch, or can
# be tested versus existing data
sys_dim = 40
obs_dim = 20
methods = ["ienks-transform"]
seed = 0

# note MDA is only defined for shifts / lags where the lag is a multiple of shift
# MDA is never defined for the classic smoother, but we will use the same parameter
# discretizations for SDA for reference values
lags = 1:3:52
shifts = [1]

# observation parameters, gamma controls nonlinearity
#gammas = Array{Float64}(1:11)
gammas = [1.0]
obs_un = 0.1
param_err = 0.03
param_infl = 1.0
#param_wlk = round.(LinRange(0, 10, 11) * 0.001, digits=3)
param_wlk = [0.0001]

# if varying nonlinearity in obs, typically take a single ensemble value
#N_ens = [21]
N_ens = 11:2:41

# inflation values, finite size versions should only be 1.0 generally
#state_infl = [1.0]
state_infl = LinRange(1.00, 1.10, 11)

# set the time series of observations for the truth-twin
time_series = [ts01]

# set MDA true or false
mdas = [true]

# time_series, method, seed, lag, shift, mda, obs_un, obs_dim, γ, param_err, param_wlk, N_ens, state_infl, param_infl = args
# load the experiments as a tuple
args = Tuple[]
for ts in time_series
    for mda in mdas
        for method in methods
            for γ in gammas
                for l in 1:length(lags)
                    lag = lags[l]
                    #shifts = lags[1:l]
                    for shift in shifts
                        for N in N_ens
                            for s_infl in state_infl
                                for wlk in param_wlk
                                    if incomplete == true
                                        #tanl = parse(Float64,ts[71:74])
                                        #diffusion = parse(Float64,ts[60:64])
                                        name = method *
                                                    "_IEEE39bus_param_seed_" * lpad(seed, 4, "0") *
                                                    "_diff_" * rpad(diffusion, 5, "0") *
                                                    "_sysD_" * lpad(sys_dim, 2, "0") *
                                                    "_obsD_" * lpad(obs_dim, 2, "0") *
                                                    "_obsU_" * rpad(obs_un, 4, "0") *
                                                    "_gamma_" * lpad(γ, 5, "0") *
                                                    "_paramE_" * rpad(param_err, 4, "0") *
                                                    "_paramW_" * rpad(wlk, 6, "0") *
                                                    "_nanl_" * lpad(nanl, 5, "0") *
                                                    "_tanl_" * rpad(tanl, 4, "0") *
                                                    "_h_" * rpad(h, 4, "0") *
                                                    "_lag_" * lpad(lag, 3, "0") *
                                                    "_shift_" * lpad(shift, 3, "0") *
                                                    "_mda_" * string(mda) *
                                                    "_nens_" * lpad(N, 3,"0") *
                                                    "_stateInfl_" * rpad(round(s_infl, digits=2), 4, "0") *
                                                    "_paramInfl_" * rpad(round(param_infl, digits=2), 4, "0") *
                                                    ".jld2"

                                        fpath = "/x/capa/scratch/cgrudzien/power_grid_data/" * method * "/"
                                        try
                                            f = load(fpath*name)
                                        catch
                                            tmp = (ts, method, seed, lag, shift, mda, obs_un, obs_dim, γ, param_err, wlk, N, s_infl, param_infl)
                                            push!(args, tmp)
                                        end
                                    else
                                        tmp = (ts, method, seed, lag, shift, mda, obs_un, obs_dim, γ, param_err, wlk, N, s_infl, param_infl)
                                        push!(args, tmp)
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

name = "/home/cgrudzien/DataAssimilationBenchmarks/data/input_data/iterative_param_smoother_input_args.jld2"
save(name, "experiments", args)

for j in 1:length(args)
    f = open("./submit_job.sl", "w")
    write(f,"#!/bin/bash\n")
    write(f,"#SBATCH -n 1\n")
    # slow partition is for Okapi, uncomment when necessary
    #write(f,"#SBATCH -p slow\n")
    write(f,"#SBATCH -o ensemble_run.out\n")
    write(f,"#SBATCH -e ensemble_run.err\n")
    write(f,"julia SlurmExperimentDriver.jl " * "\"" *string(j) * "\"" * " \"iterative_smoother_param\"")
    close(f)
    my_command = `sbatch  submit_job.sl`
    run(my_command)
end

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

end

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks


text:
# If the numbers 1 to 5 are written out in words: one, two, three, four, five,
# then there are 3 + 3 + 5 + 4 + 4 = 19 letters used in total.
#
# If all the numbers from 1 to 1000 (one thousand) inclusive were written out
# in words, how many letters would be used?
#
# NOTE: Do not count spaces or hyphens. For example, 342 (three hundred and
# forty-two) contains 23 letters and 115 (one hundred and fifteen) contains 20
# letters. The use of "and" when writing out numbers is in compliance with
# British usage.

using ProjectEulerSolutions

# Convert integer to string, remove hyphens and spaces, then add up lengths.
function p017solution(n::Integer=2)::Integer

    total_length = 0
    for i = 1:n
        total_length += length(replace(integer_to_string(i), r"[- ]" => s""))
    end

    return total_length
end

p017 = Problems.Problem(p017solution)

Problems.benchmark(p017, 1_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using LinearAlgebra

const liblapack = Base.liblapack_name
const libblas = Base.libblas_name
import LinearAlgebra.BLAS.@blasfunc
import LinearAlgebra: BlasFloat, BlasInt

function ksvd!(A, U, S, VT, work, calc_work = false)
    m, n   = size(A)
    minmn  = min(m, n)
    lwork  = BlasInt(length(work))
    info   = Ref{BlasInt}()

    jobu = 'S'
    jobvt = 'S'

    # calculate work
    if calc_work
        lwork = BlasInt(-1)
    end

    ccall((@blasfunc(dgesvd_), liblapack), Cvoid,
                      (Ref{UInt8}, Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                       Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ptr{Float64},
                       Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{Float64},
                       Ref{BlasInt}, Ptr{BlasInt}),
                      jobu, jobvt, m, n, A, m, S, U, m, VT, n,
                      work, lwork, info)

    if calc_work
        lwork = BlasInt(real(work[1]))
        resize!(work, lwork)
    end
end

function kgemv!(tA, rank, alpha, A, x, beta, y)
    ccall((@blasfunc(dgemv_), libblas), Cvoid,
                (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ref{Float64},
                 Ptr{Float64}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
                 Ref{Float64}, Ptr{Float64}, Ref{BlasInt}),
                 tA, size(A,1), rank, alpha,
                 A, size(A,1), x, 1,
                 beta, y, 1)
end

# adapted from http://pi.math.cornell.edu/~web6140/TopTenAlgorithms/QRUpdate.html
# w is a work vector
function qrupdate!(Q, R, u, v, w)
    # Compute the QR factorization of Q*R + u*v':

    # Note that Q*R + u*v' = Q*(R + w*v') with w = Q'*u:
    mul!(w, Q', u)
    n = size(Q, 1)

    # Convert R+w*v' into upper-hessenberg form using n-1 Givens rotations:
    for k = n-1:-1:1
        c, s, r = LinearAlgebra.givensAlgorithm(w[k], w[k+1])
        w[k+1] = 0.; w[k] = r

        # Compute G*R[k:k+1,:] and Q[:,k:k+1]*G', where G = [c -s ; s c]
        for j = 1:n
            newrow = c*R[k,j] + s*R[k+1,j]
            R[k+1,j] = -s*R[k,j] + c*R[k+1,j]
            R[k,j] = newrow
            newcol = c*Q[j,k] + s*Q[j,k+1]
            Q[j,k+1] = -s*Q[j,k] + c*Q[j,k+1]
            Q[j,k] = newcol
        end
    end

    # R <- R + w*v' is now upper-hessenberg:
    R[1,:] .+= w[1].*v

    # Convert R from upper-hessenberg form to upper-triangular form using n-1 Givens rotations:
    for k = 1:n-1
        c, s, r = LinearAlgebra.givensAlgorithm(R[k,k], R[k+1,k])

        # Compute G*R[k:k+1,:] and Q[:,k:k+1]*G', where G = [c -s ; s c]
        for j = 1:n
            newrow = c*R[k,j] + s*R[k+1,j]
            R[k+1,j] = -s*R[k,j] + c*R[k+1,j]
            R[k,j] = newrow
            newcol = c*Q[j,k] + s*Q[j,k+1]
            Q[j,k+1] = -s*Q[j,k] + c*Q[j,k+1]
            Q[j,k] = newcol
        end
    end
end


text:
"""
    diagnostics(ss::StateSpace{T}) where T

Run diagnostics and print results
"""
function diagnostics(ss::StateSpace{T}) where T
    println("==============================================================")
    println("                    Running diagnostics...                    ")

    e = residuals(ss)
    pvalueJB = jarquebera(e)
    pvalueLB = ljungbox(e)
    pvalueH  = homoscedast(e)

    p = ss.model.dim.p
    for j = 1:p
        if p > 1
            println("--------------------------------------------------------------")
            println("                       Series $j of $p...                       ")
        end
        println("                Jarque-Bera: p-value = $(round(pvalueJB[j], digits = 5))                ")
        println("                  Ljung-Box: p-value = $(round(pvalueLB[j], digits = 5))                 ")
        println("           Homoscedasticity: p-value = $(round(pvalueH[j], digits = 5))              ")
    end

    println("==============================================================")

    return Diagnostics(pvalueJB, pvalueLB, pvalueH)
end

"""
    residuals(ss::StateSpace{T}) where T

Obtain standardized residuals from a state-space model
"""
function residuals(ss::StateSpace{T}) where T
    e = Matrix{T}(undef, ss.model.dim.n, ss.model.dim.p)
    for j = 1:ss.model.dim.p
        for t = 1:ss.model.dim.n
            e[t, j] = ss.filter.v[t, j] / sqrt(ss.filter.F[j, j, t])

            # Replace Inf with large numbers to avoid numerical issues
            e[t, j] == Inf && (e[t, j] = 1e6)
            e[t, j] == -Inf && (e[t, j] = -1e6)
        end
    end
    return e
end

"""
    jarquebera(e::Matrix{T}) where T

Run Jarque-Bera normality test and return p-values
"""
function jarquebera(e::Matrix{T}) where T
    n, p = size(e)
    pvalue = Vector{T}(undef, p)
    for j = 1:p
        m1 = sum(e[:, j])/n
        m2 = sum((e[:, j] .- m1).^2)/n
        m3 = sum((e[:, j] .- m1).^3)/n
        m4 = sum((e[:, j] .- m1).^4)/n
        b1 = (m3/m2^(3/2))^2
        b2 = (m4/m2^2)
        statistic = n * b1/6 + n*(b2 - 3)^2/24
        d = Chisq(2.)
        pvalue[j] = 1. - cdf(d, statistic)
    end
    return pvalue
end

"""
    ljungbox(e::Matrix{T}; maxlag::Int = 20) where T

Run Ljung-Box independence test and return p-values
"""
function ljungbox(e::Matrix{T}; maxlag::Int = 20) where T
    n, p = size(e)
    n <= maxlag && (maxlag = n-1)
    pvalue = Vector{T}(undef, p)
    for j = 1:p
        acor = autocor(e[:, j], 1:maxlag)
        Q = n*(n+2)*sum((acor[k]^2)/(n-k) for k = 1:maxlag)
        d = Chisq(maxlag)
        pvalue[j] = 1 - cdf(d, Q)
        isnan(pvalue[j]) && (pvalue[j] = 0.0)
    end
    return pvalue
end

"""
    homoscedast(e::Matrix{T}) where T

Run homoscedasticity test and return p-values
"""
function homoscedast(e::Matrix{T}) where T
    n, p = size(e)
    t1 = floor(Int, n/3)
    t2 = ceil(Int, 2*(n/3)+1)
    dist = FDist(n-t2+1, t1)
    pvalue = Vector{T}(undef, p)
    for j = 1:p
        H = sum(e[t, j]^2 for t = t2:n)/sum(e[t, j]^2 for t = 1:t1)
        Δ = abs(cdf(dist, H) - 0.5)
        H1 = quantile(dist, 0.5 - Δ)
        H2 = quantile(dist, 0.5 + Δ)
        pvalue[j] = cdf(dist, H1) + (1 - cdf(dist, H2))
    end
    return pvalue
end


text:
using Random

ENV["JULIA_EDITOR"] = "vim"

println("Setup successful")

max_stars_repo_path: Chapter01/06. How to customize Julia on startup/startup.jl


text:
module Transport
using JSON

using Mea.Periodize

# maybe define an options Dict, that would be taken by almost all of the functions here for example:
# options = Dict("fctper"=>"make_akw2green", "beta_"=>12.0, "cutoff"=>100.0, "maxevals"=>100000,
#                "libintegrator => "cubature")
# funtion(modelvector, options)
#or maybe just for coefstrans: coefstrans = Transport.coefstrans(modelvec, beta_, options=options)
#                                   or    = Transport.coefstrans(modelvec, beta_; options...)
# where options is now a Dict of Symbol=>value (options = Dict(Symbol("fctper")=>"make_akw2green ...))

const cutoffdefault = 20.0

function dfdw(beta_::Float64, ww::Float64)
    return (-beta_*exp(beta_*ww)/(1.0 + exp(beta_*ww))^2.0   )
end

function vz2_int(kk::Array{Float64 ,1}) #v_perp^2.0 integrated in z

    (kx, ky) = (kk[1], kk[2])
    result = 2.0
    #result = 2.0*(cos(kx) - cos(ky))^(4.0)
    return result
end

 function make_lab_kintegrand(model::Periodize.Model; fctper::String="make_akw2green")
    #fctper = "make_akw2green, make_akw2cum, make_akw2trace" as defined in Periodize.jl
    make_akw2 = getfield(Periodize, Symbol(fctper))
    akw2 = make_akw2(model)

    function lab_kintegrand(kk::Array{Float64, 1})
          return(akw2(kk)*2.0)#vz2_int(kk))
      end
       return lab_kintegrand
  end

  function make_lab_kintegrand_cuba(model::Periodize.Model; fctper::String="make_akw2green")

    make_akw2 = getfield(Periodize, Symbol(fctper))
    akw2 = make_akw2(model)

    function lab_kintegrand_cuba(kk::Array{Float64, 1}, ff::Array{Float64, 1})
        kk_scaled = (-pi + 2*pi*kk)
        ff[1] = (4.0*pi*pi*akw2(kk_scaled)*2.0)#vz2_int(kk_scaled))
    end
    return lab_kintegrand_cuba
   end

function calc_labk(modelvec::Periodize.ModelVector, beta_::Float64; cutoff::Float64=cutoffdefault, maxevals::Int64=100000,
                    libintegrator::String="cubature", fctper::String="make_akw2green") #calculate the k integral of the L_ab coeffiecient
    #returns: Array{Float64, 1} size of w_vec that will be integrated in frequency

  modelvec = deepcopy(modelvec)
  cutoff = cutoff
  cutoffidx=0
  for (ii, ww) in enumerate(modelvec.wvec_)
    if abs(ww*beta_) < cutoff
      cutoffidx = ii
      break
    end
  end

  modelvec.wvec_ = modelvec.wvec_[cutoffidx:end-cutoffidx]
  modelvec.sEvec_c_ = modelvec.sEvec_c_[cutoffidx:end-cutoffidx, :, :]

  kwargs = Dict(:fctper=>fctper)
  #println("kwargs... = ", kwargs...)
  if libintegrator == "cubature"
    result = Periodize.calcintegral(modelvec, make_lab_kintegrand, maxevals=maxevals; kwargs...)
  elseif libintegrator == "cuba"
    result = Periodize.calcintegral_cuba(modelvec, make_lab_kintegrand_cuba, maxevals=maxevals; kwargs...)
  else
    error("Invalid value for libintegrator")
  end

  integrand_w = zeros(Float64, size(result)[1])
  for ii in 1:size(result)[1]
      integrand_w[ii] = -dfdw(beta_, result[ii, 1])*result[ii, 2]
  end

  return (modelvec.wvec_, integrand_w)

end

function calc_l11(integrand_w::Array{Float64, 1}, wvec::Array{Float64, 1}, beta_::Float64)

    l11 = 0.0
    assert(size(wvec) == size(integrand_w))

   # frequency integration
   for ii in 1:size(wvec)[1]-1                       # y                                 x
       l11 += 0.5 * (integrand_w[ii] + integrand_w[ii+1]) * (wvec[ii+1] - wvec[ii]) / (2.0*pi)
   end

   l11 /= beta_
   #println(l11)
   return l11
end

function calc_l21(integrand_w::Array{Float64, 1}, wvec::Array{Float64, 1}, beta_::Float64)

    l21 = 0.0
    assert(size(wvec) == size(integrand_w))

    for jj in 1:size(integrand_w)[1]
        integrand_w[jj] *= wvec[jj]
    end

   # frequency integration
   for ii in 1:size(wvec)[1]-1                       # y                                 x
       l21 += 0.5 * (integrand_w[ii] + integrand_w[ii+1]) * (wvec[ii+1] - wvec[ii]) / (2.0*pi)
   end

   l21 /= beta_
   #println(l21)
    return l21
end

function calc_l22(integrand_w::Array{Float64, 1}, wvec::Array{Float64, 1}, beta_::Float64)

    l22 = 0.0
    assert(size(wvec) == size(integrand_w))

    for jj in 1:size(integrand_w)[1]
        integrand_w[jj] *= (wvec[jj])^(2.0)
    end

   # frequency integration
   for ii in 1:size(wvec)[1]-1                       # y                                 x
       l22 += 0.5 * (integrand_w[ii] + integrand_w[ii+1]) * (wvec[ii+1] - wvec[ii]) / (2.0*pi)
   end

   l22 /= beta_
   #println(l22)
   return l22
end

function calc_sigmadc(modelvec::Periodize.ModelVector, beta_::Float64; cutoff::Float64=cutoffdefault)

    (wvec, integrand_w) = calc_labk(modelvec, beta_, cutoff=cutoff)
    sigmadc = beta_*calc_l11(integrand_w, wvec,  beta_)
    #println("sigmadc = ", sigmadc)
    return sigmadc
end

# function calc_seebeck(modelvec::Periodize.ModelVector, beta_::Float64; cutoff=cutoffdefault)
#   seebeck = -beta_*calc_l21(modelvec, beta_, cutoff=cutoff)/calc_l11(modelvec, beta_, cutoff=cutoff)
#   #println("seebeck = ", seebeck)
#   return seebeck
# end

function calc_n(modelvec::Periodize.ModelVector, beta_::Float64; fout_name::String="dos_n.dat",   maxevals::Int64=100000, fctper::String="make_akwgreen")

    dos = Periodize.calcintegral(modelvec, getfield(Periodize, Symbol(fctper)), fout_name=fout_name, maxevals=maxevals)
    n = 0.0
    for jj in 1:size(dos)[1]
        ww = dos[jj, 1]
        dos[jj, 2] *= 1.0/(exp(beta_*ww) + 1.0)
    end

   # frequency integration
   for ii in 1:size(dos)[1]-1     # y                                 x
       n += 0.5 * (dos[ii, 2] + dos[ii+1, 2]) * (dos[ii+1, 1] - dos[ii, 1]) / (2.0*pi)
   end

   #println("n = ", n)
   return n
end

function coefstrans(modelvec::Periodize.ModelVector, beta_::Float64; cutoff::Float64=cutoffdefault, fout_name::String="dos.dat", maxevals::Int64=100000,
                    libintegrator::String="cubature", fctper::String="make_akw2green")

    dd = Dict{String, Float64}()
    (wvec, integrand_w) = calc_labk(modelvec, beta_, cutoff=cutoff, maxevals=maxevals, libintegrator=libintegrator, fctper=fctper)
    dd["l11"] =  calc_l11(integrand_w, wvec, beta_)
    dd["l21"] = calc_l21(integrand_w, wvec, beta_)
    dd["l22"] = calc_l22(integrand_w, wvec, beta_)
    dd["sigmadc"] = beta_*dd["l11"]
    dd["seebeck"] = -beta_*dd["l21"]/dd["l11"]

    ss = split(fctper, "2")
    fctdos = ss[1] * ss[2] # fctper = make_akw2green -> fctdos = make_akwgreen
    dd["n"] = calc_n(modelvec, beta_, fout_name=fout_name, maxevals=maxevals, fctper=fctdos)

    fout = "coefstrans" * fctper * ".json"
    if isfile(fout)
        mv(fout, fout * "." *  Dates.format(now(), "YYYY-uu-HHh:MM") )
    end

    open(fout, "w") do f
            write(f, JSON.json(dd))
    end

    return dd
end

end


text:
function OnBoardAntennaTurnAroundRangeModifier(arg0::Vector3D)
    return OnBoardAntennaTurnAroundRangeModifier((Vector3D,), arg0)
end

function get_parameters_drivers(obj::OnBoardAntennaTurnAroundRangeModifier)
    return jcall(obj, "getParametersDrivers", List, ())
end

function modify(obj::OnBoardAntennaTurnAroundRangeModifier, arg0::EstimatedMeasurement)
    return jcall(obj, "modify", void, (EstimatedMeasurement,), arg0)
end


text:
using Images
function MoG(img::Array{T,3}) where T <: Number
	(nr,nc,nch) = size(img)
    grad1, grad2 = imgradients(img, KernelFactors.sobel)
	mog = zeros(nr,nc,nch)
    for i=1:nch
        mog[:,:,i] = magnitude(grad1[:,:,1],grad2[:,:,2])
	end
	return mog
end
function MoG(img::Array{T,2}) where T<: Number
    (nr,nc) = size(img)
    grad1, grad2 = imgradients(img, KernelFactors.sobel)
    return magnitude(grad1,grad2)
end
function MoG(RGB::Array{T1,3}, D::Array{T2,2}) where {T1<:Number, T2<:Number}
	mogRGB   = MoG(RGB)
	mogD     = MoG(D)
	mog =  max.(mogRGB[:,:,1],mogRGB[:,:,2], mogRGB[:,:,3], mogD)
    mog = mog / maximum(mog)
    return mog
end
function MoG(I::Array{T1,2}, D::Array{T2,2}) where T1<:Number where T2<:Number
    mogI = MoG(I)
    mogD = MoG(D)
    mog  = max.(mogI, mogD)
    return mog./maximum(mog)
end
struct Edge{EdgeType}
    data::EdgeType
    node
end
struct Node{NodeDataType, EdgeType}
    index::Integer
    data::NodeDataType
    neighbors::Vector{Edge{EdgeType}}
end
mutable struct Graph{NodeDataType, EdgeType}
    nodes::Dict{Integer,Node{NodeDataType, EdgeType}}
    indices_count::Integer
    edge_weight::Function
    merge_nodes_data::Function
    function Graph{NodeDataType, EdgeType}(edge_weight::Function, merge_nodes_data::Function) where {NodeDataType, EdgeType}
        nodes = Dict{Integer,Node{NodeDataType, EdgeType}}()
        return new(nodes,0, edge_weight, merge_nodes_data)
    end
end

import Base.getindex
import Base.haskey
function getindex(r::Graph, k::Integer)
    return r.nodes[k].data
end
function haskey(r::Graph, k::Integer)
    return haskey(r.nodes,k)
end

function add_vertex!(r::Graph{D,E}, data::D,neighbors::Vector{Edge{E}}=Edge{E}[]) where  D where E
    r.indices_count+=1;
    r.nodes[r.indices_count] = Node{D,E}(r.indices_count, data, neighbors )
    return r.indices_count
end
function number_of_adjacents(r::Graph{D,E}, i::Integer)  where D where E
    return length(r.nodes[i].neighbors)
end
function adjacents(r::Graph{D,E}, i::Integer) where D where E
    return vec([(edge.node.index, edge.data) for edge in r.nodes[i].neighbors])
end
function remove_connection_to_node!(n::Node{D,E}, i::Integer) where D where E
    filter!(n->n.node.index !=i, n.neighbors)
end
function vertices(r::Graph{D,E}) where D where E
    return collect(keys(r.nodes))
end
function remove_vertex!(r::Graph{D,E}, i::Integer) where D where E
    for a in r.nodes[i].neighbors
        remove_connection_to_node!(a.node, i)
    end
    delete!(r.nodes, i)
end
function connect!(r::Graph{D,E}, i::Integer, j::Integer) where D where E
    edge = r.edge_weight(r[i],r[j])
    pos = findfirst(edge->edge.node.index==j,r.nodes[i].neighbors)
    if pos == 0
        push!(r.nodes[i].neighbors, Edge{E}(edge, r.nodes[j]))
    end
    pos = findfirst(edge->edge.node.index==i,r.nodes[j].neighbors)
    if pos==0
        push!(r.nodes[j].neighbors, Edge{E}(edge, r.nodes[i]))
    end
end
function merge_vertices!(r::Graph{D,E}, i::Integer, j::Integer) where D where E
    node_i = r.nodes[i]
    node_j = r.nodes[j]
    remove_vertex!(r, j)
    remove_vertex!(r, i)
    new_data        = r.merge_nodes_data(node_i.data, node_j.data)
    new_vertex_id   = add_vertex!(r, new_data)
    new_neighbors   = filter(x->(x.node.index!=i) && (x.node.index!=j),vcat(node_i.neighbors, node_j.neighbors))
    for edge in new_neighbors
        connect!(r, new_vertex_id, edge.node.index)
    end
end
abstract type AbstractRAGNode end
mutable struct RAGEdge
    wd::Float64
    wb::Float64
end
mutable struct RAGNode <: AbstractRAGNode
    w::MvWatsonExponential
    π::Float64
    mean_color::Vector
    boundary::Matrix{Bool}
    pixel_indices::Vector{Integer}
end

"""
```
function boundaries(segmented_image::Matrix; w::Integer=1)
```
Returns the boundary pixels of every region from the `segmented_image`
"""
function boundaries(segmented_image::Matrix; w::Integer=1)
    (nr, nc) = size(segmented_image)
    b = Dict()
    for j in unique(segmented_image)
        img_object = segmented_image.==j
        b[j] = dilate(img_object) - erode(img_object)
    end
    return b
end
function edge_weight(Ni::RAGNode, Nj::RAGNode,mog)
    wd = min(bregman_divergence(Ni.w, Nj.w), bregman_divergence(Nj.w, Ni.w))
    boundary_intersection = find(Ni.boundary .& Nj.boundary)
    wi = sum(mog[boundary_intersection]) / (length(boundary_intersection))
    return RAGEdge(wd,wi)
end

function create_node(::Type{RAGNode}, i, RGB, segmented_image, data, bounda)
    (nr,nc, nd) = size(RGB)
    pixel_indices = findall(segmented_image.==i)
    π = length(pixel_indices)/(nr*nc)
    wmean = vec(mean(data.normal[1][:,pixel_indices],2))
    N = MvWatsonExponential(wmean)
    (rows,cols) = ind2sub((size(RGB,1),size(RGB,2)),pixel_indices)
    mean_color = vec(mean(broadcast_getindex(RGB,rows,cols,collect(1:3)'),1))
   return RAGNode(N, π, mean_color, bounda[i], pixel_indices)
end
"""
```
function create_rag(m::Model, RGB, XYZ, segmented_image::Matrix{Int64}, data::RGBDData)
```
Create the region adjacency graph
"""

function create_rag(::Type{NodeType}, ::Type{EdgeType}, RGB, XYZ, segmented_image::Matrix{T}, data,mog) where NodeType where EdgeType where T<:Integer
    (nr,nc, nd) = size(XYZ)
    clusters    = unique(segmented_image)
    rag         = Graph{NodeType,EdgeType}((n1,n2)->edge_weight(n1,n2,mog), merge_data)
    seg_image   = SegmentedImage(segmented_image, collect(1:maximum(segmented_image)), Dict{Int,Real}(), Dict{Int,Int}() )
    G, vert_map = region_adjacency_graph(seg_image,(x,y)->1);

    bounda      = boundaries(segmented_image)
    vertex_map = Dict{Integer,Integer}()
    # Vertex Creation
    for i in clusters
        index = add_vertex!(rag, create_node(NodeType, i, RGB, segmented_image, data, bounda))
        vertex_map[i] = index
    end
    for edge in edges(G)
        i  = edge.src
        j  = edge.dst
        Ni = rag[vertex_map[i]]
        Nj = rag[vertex_map[j]]
        connect!(rag, vertex_map[i],vertex_map[j])
    end
    return rag
end

struct RegionMergingConfig
    NodeType::DataType
    EdgeType::DataType
    kp::Float64
    thd::Float64
    thb::Float64
    thr::Float64
    small_region_threshold::Integer
end
function RegionMergingConfig()
    return RegionMergingConfig(
        RAGNode,
        RAGEdge,
        5,
        3,
        0.2,
        0.9,
        5)
end
function candidancy(cfg::RegionMergingConfig, v::T) where T<:AbstractRAGNode
    return  v.w.concentration > cfg.kp
end
function elegibility(cfg::RegionMergingConfig, v::T, r::RAGEdge) where T<:AbstractRAGNode
    return ( r.wb < cfg.thb) && (r.wd < cfg.thd)
end
function should_merge(cfg::RegionMergingConfig, XYZ,  n1::T, n2::T, r::RAGEdge) where T<:AbstractRAGNode
    return candidancy(cfg,n2) && elegibility(cfg, n1, r) && consistency(cfg, XYZ, n1, n2)
end
function consistency(cfg::RegionMergingConfig, XYZ, v::T, r::T) where T<:AbstractRAGNode
    (indices_r, indices_c) = ind2sub((size(XYZ,1),size(XYZ,2)),vcat(v.pixel_indices, r.pixel_indices))
    points_3d      = broadcast_getindex(XYZ,indices_r,indices_c,collect(1:3)')
    if size(points_3d,1) > 3
        points_3d            = points_3d'
        (B, P, inliers)     = ransacfitplane(points_3d,Inf,false)
        plir                = length(inliers) / length(indices_r)
   return plir > cfg.thr
    else
        return false
    end
end
function merge_data(i::RAGNode, j::RAGNode)
    new_π         = i.π + j.π
    new_η         = (i.w.n*i.π + j.w.n*j.π)/(new_π)
    new_indices   = vcat(i.pixel_indices, j.pixel_indices)
    new_color     = (i.mean_color + j.mean_color)/2
    new_boundaries = i.boundary .| j.boundary
    new_node      = RAGNode(MvWatsonExponential(new_η), new_π, new_color, new_boundaries, new_indices)

end
"""
```
function remove_small_regions(rag,RGB; small_area::Integer=15)
```
Remove segments below a minimum size, assigning them to a neighbor segment with the most similar average color.
"""
function remove_small_regions(r::Graph{NodeType, EdgeType}, seg_image, mog, small_area::Integer=1000) where NodeType where EdgeType
    prev_number_of_segments = Inf
    number_of_segments = length(r.nodes)
    while number_of_segments < prev_number_of_segments
        prev_number_of_segments = number_of_segments
        areas = [(key, length(r[key].pixel_indices)) for key in vertices(r)]
        sort!(areas, by=x->x[2])
        for (k, area) in areas
            if haskey(r,k) && (length(r[k].pixel_indices) < small_area)
                adjacent_colors = zeros(size(r[k].mean_color,1),number_of_adjacents(r,k))
                j = 1
                for (i,w) in adjacents(r,k)
                    adjacent_colors[:,j] = r[i].mean_color
                    j=j+1
                end
                dist  = colwise(Euclidean(), r[k].mean_color,adjacent_colors)

                min_dist, min_neigh = findmin(dist)
                merge_vertices!(r, k, adjacents(r,k)[min_neigh][1])

            end
        end
        number_of_segments = length(r.nodes)
    end

end
"""
```
function merge_parts(cfg::RegionMergingConfig, m::Model, RGB, XYZ, segmented_image::Matrix{<:Integer}, data::RGBDData )
```

"""
function merge_parts(cfg::RegionMergingConfig, RGB, XYZ, segmented_image::Matrix{<:Integer}, data )
    (nr,nc)  = size(segmented_image)
    mog      = MoG(RGB,XYZ[:,:,3])
    rag      = create_rag(cfg.NodeType, cfg.EdgeType, RGB, XYZ, segmented_image, data, mog )
    remove_small_regions(rag, segmented_image, mog, cfg.small_region_threshold)

    last_number_of_nodes = 0
    number_of_nodes = length(rag.nodes)
    while number_of_nodes != last_number_of_nodes
        last_number_of_nodes = number_of_nodes
        for (i, v) in rag.nodes
            merged = false
            if candidancy(cfg,rag[i])
                regions  = sort(adjacents(rag,i), by=x->x[2].wd)
                for (r,w) in regions
                    if should_merge(cfg,XYZ, rag[i], rag[r], w)
                        merge_vertices!(rag,i, r)
                        merged = true
                        break
                    end
                end
                if merged
                    break
                end

            end
        end
        number_of_nodes = length(rag.nodes)
    end

    return rag
end

function process_labels(cfg::RegionMergingConfig, RGB, D, labels, xss)

#    labels  = round.(Integer,mapwindow(median!, labels, (3,3)))
    labels = Images.label_components(labels)
    #new_model = Model(maximum(labels))
    #n         = size(D,1)*size(D,2)
    #rag       = merge_parts(cfg, RGB, D, labels, xss)
    #cluster         = similar(labels)
    #j=1
    #for k in  vertices(rag)
    #    cluster[rag[k].pixel_indices]=j
    #    j=j+1
    #end

    #cluster  = round.(Integer,mapwindow(median!, cluster, (3,3)))
    #return cluster
    return labels
end


text:
using Random, Distributions, KernelDensity, Plots, LaTeXStrings; pyplot()
Random.seed!(1)

function tStat(mu0,mu,sig,n)
    sample = rand(Normal(mu,sig),n)
    xBar   = mean(sample)
    s      = std(sample)
    (xBar-mu0)/(s/sqrt(n))
end

mu0, mu1A, mu1B = 20, 22, 24
sig, n = 7, 5
N = 10^6
alpha = 0.05

dataH0  = [tStat(mu0,mu0,sig,n) for _ in 1:N]
dataH1A = [tStat(mu0,mu1A,sig,n) for _ in 1:N]
dataH1B = [tStat(mu0,mu1B,sig,n) for _ in 1:N]
dataH1C = [tStat(mu0,mu1B,sig,2*n) for _ in 1:N]
dataH1D = [tStat(mu0,mu1B,sig/2,2*n) for _ in 1:N]

tCrit = quantile(TDist(n-1),1-alpha)
estPwr(sample) = sum(sample .> tCrit)/N

println("Rejection boundary: ", tCrit)
println("Power under H0:  ", estPwr(dataH0))
println("Power under H1A: ", estPwr(dataH1A))
println("Power under H1B (mu's farther apart): ", estPwr(dataH1B))
println("Power under H1C (double sample size): ", estPwr(dataH1C))
println("Power under H1D (like H1C but std/2): ", estPwr(dataH1D))

kH0  = kde(dataH0)
kH1A = kde(dataH1A)
kH1D = kde(dataH1D)
xGrid = -10:0.1:15

plot(xGrid,pdf(kH0,xGrid),
	c=:blue, label="Distribution under H0")
plot!(xGrid,pdf(kH1A,xGrid),
	c=:red, label="Distribution under H1A")
plot!(xGrid,pdf(kH1D,xGrid),
	c=:green, label="Distribution under H1D")
plot!([tCrit,tCrit],[0,0.4],
	c=:black, ls=:dash, label="Critical value boundary",
	xlims=(-5,10), ylims=(0,0.4), xlabel=L"\Delta = \mu - \mu_0")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
struct BoolVariate <: AbstractBoolVariate
    name::String
    data::AbstractVector{Bool}
end

Base.length(boolvar::AbstractBoolVariate) = length(boolvar.data)

function slice(boolvar::BoolVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    slice(boolvar.data, fromobs, toobs, slicelength)
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
using BAI
using Test
using LinearAlgebra

tolerance = 1e-13

mu = [1.2; 3.5; 1.1; 0.4; 2.9]
omega = [0.1; 0.2; 0.4; 0.1; 0.2]
lambda = [2.3; 7.1; 2.4; 0.2; 7.1]
laplacian = zeros(Float64, 5, 5)
for i in 1:length(lambda), j in (i + 1):length(lambda)
    w = rand()
    laplacian[i, j], laplacian[j, i] = -w, -w
end
for i in 1:length(lambda)
    laplacian[i, i] -= sum(laplacian[:, i])
end
smoothness = 0
for i in 1:length(lambda), j in 1:length(lambda)
    if i != j
        global smoothness += -laplacian[i, j] * (lambda[i] - lambda[j])^2 / 2
    end
end
objective = 0
for i in 1:length(lambda)
    global objective += omega[i] * (mu[i] - lambda[i])^2 / 2
end
@test norm(BAI.smoothness(lambda, laplacian) - smoothness) < tolerance
@test norm(BAI.objective(mu, omega, lambda) - objective) < tolerance

max_stars_repo_name: tomaskocak/BAI

max_issues_repo_name: tomaskocak/BAI

max_forks_repo_name: tomaskocak/BAI


text:
using Transits.Orbits: relative_position

@testset "simple orbit" begin
    period = 3.456
    t0 = 1.45
    b = 0.5
    duration = 0.12
    t = t0 .+ range(-2.0*period, 2.0*period, length=5000)
    m0 = @. abs(mod(t - t0 + 0.5*period, period) - 0.5*period) < 0.5 * duration

    orbit = SimpleOrbit(;period, t0, b, duration)
    pos = @. relative_position(orbit, t)
    bs = map(v -> sqrt(v[1]^2.0 + v[2]^2.0), pos)
    zs = map(v->v[3], pos)
    m = @. (bs ≤ 1.0) & (zs > 0.0)
    @test orbit.ref_time == -0.278
    @test m ≈ m0

    #idxs = in_transit(orbit, t)
    #@test all(bs[idxs] .≤ 1.0)
    #@test all(zs[idxs] .> 0.0)

    # TODO: Do we want this?
    #pos  = @. star_position(orbit, t)
    #for vec in pos
    #    @test all(vec .≈ 0.0)
    #end
end

@testset "simple limbdark" begin
    period = 3.456
    t0 = 1.45
    b = 0.5
    duration = 0.12

    t = t0 .+ range(-2.0*period, 2.0*period, length=5000)
    m0 = @. abs(mod(t - t0 + 0.5*period, period) - 0.5*period) < 0.5 * duration
    orbit = SimpleOrbit(;period, t0, b, duration)
    ld = PolynomialLimbDark([0.2, 0.3])

    lc1 = @. ld.(orbit, t, 0.01)
    @test all(lc1[.!m0] .== 1.0)

    ldt = IntegratedLimbDark(ld)
    lc2 = @. ldt.(orbit, t, 0.01, 0.01)
    @assert lc1 ≉ lc2
end


text:
@testset "Schreier Vectors and Stabilizers" begin
    Random.seed!(1);
    SIZE=30
    G = PermutationGroups.AbstractAlgebra.SymmetricGroup(SIZE);
    S = [rand(G) for _ in 1:3];

    pt = 1
    tr = Orbit1(Transversal, S, pt)

    import PermutationGroups.Word

    for OrbT in OrbitTypes
        schr = Schreier(OrbT, S, pt, ^)
        @test schr[pt] == one(G)

        g,h,k = S

        @test representative(S, schr.orb, pt) == one(G)
        @test representative(S, schr.orb, pt^g) == g
        @test representative(S, schr.orb, (pt^g)^h) == g*h
        @test representative(S, schr.orb, pt^(g*h)) == g*h

        @test Word(inv.(S), schr.orb, pt) == Word(Int[])
        @test Word(inv.(S), schr.orb, pt^g) == Word(Int[1])
        @test Word(inv.(S), schr.orb, (pt^g)^h) == Word(Int[2,1])
        @test Word(inv.(S), schr.orb, pt^(g*h)) == Word(Int[2,1])

        z = g*h*k
        δ = pt^z
        @test Word(inv.(S), schr.orb, δ, ^) == Word([3,2,1])

        @test schr[δ] == z
        @test all([pt^schr[o] == o for o in schr.orb])
        @test all(schr[o] == tr[o] for o in tr)
    end
end

@testset "Schreier-Sims unit tests" begin
    S = [perm"(1,2,4,3)(5)", perm"(2,5,4)"]
    @test PermutationGroups.initial_bsgs(S)[1] == [1,2]
    @test PermutationGroups.initial_bsgs(S)[2] == [S, [S[2]]]

    a,b = [perm"(1,2,4,3)(5)", perm"(1,2,5,4)"]
    B, S = PermutationGroups.initial_bsgs([a,b])
    @test B == [1]
    @test S[1] == [a,b]

    trs = [Schreier(gens, pt) for (gens, pt) in zip(S, B)]
    @test collect(trs[1]) == [1,2,4,5,3]
    depth = 1
    Δ = trs[depth]

    # first point on the orbit is β = 1:
    u = Δ[1]
    # the first Schreier generator is trivial
    @test isone(u*a*getinv(Δ, 1^a))
    # but the second is not
    @test u*b*getinv(Δ, 1^b) == perm"(2,5)(3,4)"
    c = perm"(2,5)(3,4)"
    @test sift(c, B, trs) == (c, 2)

    # so we extend the basis
    push!(B, firstmoved(c))
    # add c to generators at depth 2
    push!(S, [c])
    # and compute the orbit of the added point at depth 2
    push!(trs, Schreier(S[depth+1], 2))
    depth = 2
    @test sift(c, B, trs) == (perm"(5)", 3)

    u = Δ[2]
    @test u == a
    # Schreier generators are trivial
    @test u*a == Δ[2^a]
    @test u*b == Δ[2^b]
    # so we move to the next point on the orbit
    u = Δ[4] # β = 4
    @test u*a == Δ[4^a]
    @test u*b*getinv(Δ, 4^b) == a^2*b
    d = a^2*b
    @test sift(d, B, trs) == (d, 2)

    push!(S[2], d)
    trs[2] = Schreier(S[2], B[2])
    @test collect(trs[2]) == [2,5,3,4]

    @test sift(a^2*b, B, trs) == (perm"(5)", 3)

    @test prod(length, trs) == 20

    B_new, S_new, trs_new = schreier_sims([a,b])
    @test B_new == B
    @test S_new == S
    @test trs_new == trs

end


text:
struct LinearSystem <: AbstractEnv
    A
    B
end

function State(env::LinearSystem)
    @unpack B = env
    n = size(B)[1]
    return function (x)
        @assert length(x) == n
        x
    end
end

function Params(env::LinearSystem)
    () -> nothing
end

function Dynamics!(env::LinearSystem)
    @unpack A, B = env
    @Loggable function dynamics!(dx, x, p, t; u)
        @log state = x
        @log input = u
        dx .= A*x + B*u
    end
end


text:
list_test = split(
"--- scanner 0 ---
404,-588,-901
528,-643,409
-838,591,734
390,-675,-793
-537,-823,-458
-485,-357,347
-345,-311,381
-661,-816,-575
-876,649,763
-618,-824,-621
553,345,-567
474,580,667
-447,-329,318
-584,868,-557
544,-627,-890
564,392,-477
455,729,728
-892,524,684
-689,845,-530
423,-701,434
7,-33,-71
630,319,-379
443,580,662
-789,900,-551
459,-707,401

--- scanner 1 ---
686,422,578
605,423,415
515,917,-361
-336,658,858
95,138,22
-476,619,847
-340,-569,-846
567,-361,727
-460,603,-452
669,-402,600
729,430,532
-500,-761,534
-322,571,750
-466,-666,-811
-429,-592,574
-355,545,-477
703,-491,-529
-328,-685,520
413,935,-424
-391,539,-444
586,-435,557
-364,-763,-893
807,-499,-711
755,-354,-619
553,889,-390

--- scanner 2 ---
649,640,665
682,-795,504
-784,533,-524
-644,584,-595
-588,-843,648
-30,6,44
-674,560,763
500,723,-460
609,671,-379
-555,-800,653
-675,-892,-343
697,-426,-610
578,704,681
493,664,-388
-671,-858,530
-667,343,800
571,-461,-707
-138,-166,112
-889,563,-600
646,-828,498
640,759,510
-630,509,768
-681,-892,-333
673,-379,-804
-742,-814,-386
577,-820,562

--- scanner 3 ---
-589,542,597
605,-692,669
-500,565,-823
-660,373,557
-458,-679,-417
-488,449,543
-626,468,-788
338,-750,-386
528,-832,-391
562,-778,733
-938,-730,414
543,643,-506
-524,371,-870
407,773,750
-104,29,83
378,-903,-323
-778,-728,485
426,699,580
-438,-605,-362
-469,-447,-387
509,732,623
647,635,-688
-868,-804,481
614,-800,639
595,780,-596

--- scanner 4 ---
727,592,562
-293,-554,779
441,611,-461
-714,465,-776
-743,427,-804
-660,-479,-426
832,-632,460
927,-485,-438
408,393,-506
466,436,-512
110,16,151
-258,-428,682
-393,719,612
-211,-452,876
808,-476,-593
-575,615,604
-485,667,467
-680,325,-822
-627,-443,-432
872,-547,-609
833,512,582
807,604,487
839,-516,451
891,-625,532
-652,-548,-490
30,-46,-14",
'\n'
)

list_challenge = let input = @__DIR__() * "/day19_input.txt"
    isfile(input) ? readlines(input) : list_test_1
end

# Test function

function testit(name, result, check::Function)
    if check(result)
        @info "$name successful, with result $result"
    else
        @warn "$name failed, with result $result"
    end
end

# Part 1

orientations = let
    rot90 = [
        [0 -1 0; 1 0 0; 0 0 1],
        [0 0 -1; 0 1 0; 1 0 0],
        [1 0 0; 0 0 -1; 0 1 0]
    ]
    Set(
        g1^n * g2^m for g1 in rot90 for g2 in rot90 for n in 0:3 for m in 0:3
    )
end

function get_scanners(list)
    breaks = [0; findall(==(""), list); length(list)+1]
    data = [list[breaks[i]+2:breaks[i+1]-1] for i in 1:length(breaks)-1]
    scanners = permutedims.(
        eval.(
            Meta.parse.(
                '[' .* [join(replace.(d, "," => " "), "; ") for d in data] .* ']'
            )
        )
    )
    return scanners
end

function part1(list)
    scanners = get_scanners(list)
    positioned = fill(false, length(scanners))
    translations = fill(zeros(Int, 3), length(scanners))
    transforms = fill([1 0 0; 0 1 0; 0 0 1], length(scanners))
    positioned[1] = true
    printstyled("positioned scanner $(k2-1)\n", color=:green)
    while any(.!positioned)
        printstyled("round\n", color=:yellow)
        for k1 in 1:length(scanners), k2 in 1:length(scanners)
            if positioned[k1] && !positioned[k2]
                sck1 = transforms[k1] * scanners[k1] .+ translations[k1]
                for g in orientations
                    sck2 =  g * scanners[k2]
                    us = [rk1 - rk2 for rk1 in eachcol(sck1) for rk2 in eachcol(sck2)]
                    uc = Dict(u => count(map(==(u), us)) for u in unique(us)) # could use StatsBase.countmap, but I want to avoid non-Standard packages
                    if any(≥(12), values(uc))
                        u = first(k for (k,v) in uc if v ≥ 12)
                        printstyled("positioned scanner $(k2-1)", color=:green)
                        printstyled(" ($(count(.!positioned)) scanners left)\n", color=:red)
                        positioned[k2] = true
                        translations[k2] = u
                        transforms[k2] = g
                    end
                end
            end
        end
    end

    positioned_scanners = [
        transforms[k] * scanners[k] .+ translations[k] for k in 1:length(scanners)
    ]
    return positioned, translations, transforms, scanners, positioned_scanners, length(Set(eachcol(reduce(hcat, positioned_scanners))))
end

testit("Test Part 1", part1(list_test)[end], ==(79))

testit("Challenge Part 1", part1(list_challenge)[end], ==(442))

# Part 2

function part2(list)
    _, translations, _ = part1(list)
    return maximum(sum(abs, r1 - r2) for r1 in translations for r2 in translations)
end

testit("Test Part 2", part2(list_test), ==(3621))
testit("Challenge Part 2", part2(list_challenge), ==(11079))

max_stars_repo_name: rmsrosa/adventofcode

max_issues_repo_name: rmsrosa/adventofcode

max_forks_repo_name: rmsrosa/adventofcode


text:
module DistanceAndSimilarity

include("Functions.jl")

export euclidianDistance
export jaccardDistance
export jaccardSimilarity
export minkowskiDistance
export manhattanDistance
export cosineDistance
export errorCheck
end

max_stars_repo_name: nicolasmagalhaes/DistanceAndSimilarity

max_issues_repo_name: nicolasmagalhaes/DistanceAndSimilarity

max_forks_repo_name: nicolasmagalhaes/DistanceAndSimilarity


text:
## SCALE TRANSFORMATIONS

#     Scale = SCALE()

# Object for dispatching on scales and functions when generating
# parameter ranges. We require different behaviour for scales and
# functions:

#      transform(Scale, scale(:log10), 100) = 2
#      inverse_transform(Scale, scale(:log10), 2) = 100

# but
#     transform(Scale, scale(log10), 100) = 100       # identity
#     inverse_transform(Scale, scale(log10), 100) = 2

struct SCALE end
Scale = SCALE()

scale(s::Symbol)   = Val(s)
scale(f::Function) = f

transform(::SCALE, ::Val{:linear}, x) = x
inverse_transform(::SCALE, ::Val{:linear}, x) = x

transform(::SCALE, ::Val{:log}, x) = log(x)
inverse_transform(::SCALE, ::Val{:log}, x) = exp(x)

transform(::SCALE, ::Val{:logminus}, x) = log(-x)
inverse_transform(::SCALE, ::Val{:logminus}, x) = -exp(x)

transform(::SCALE, ::Val{:log10}, x) = log10(x)
inverse_transform(::SCALE, ::Val{:log10}, x) = 10^x

transform(::SCALE, ::Val{:log2}, x) = log2(x)
inverse_transform(::SCALE, ::Val{:log2}, x) = 2^x

transform(::SCALE, f::Function, x) = x            # not a typo!
inverse_transform(::SCALE, f::Function, x) = f(x) # not a typo!

## SCALE INSPECTION (FOR EG PLOTTING)

"""
    scale(r::ParamRange)

Return the scale associated with a `ParamRange` object `r`. The possible
return values are: `:none` (for a `NominalRange`), `:linear`, `:log`, `:log10`,
`:log2`, or `:custom` (if `r.scale` is a callable object).
"""
scale(r::NominalRange) = :none
scale(r::NumericRange) = :custom
scale(r::NumericRange{B,T,Symbol}) where {B<:Boundedness,T} = r.scale

## ITERATORS

"""
    MLJTuning.iterator(r::NominalRange, [,n, rng])
    MLJTuning.iterator(r::NumericRange, n, [, rng])

Return an iterator (currently a vector) for a `ParamRange` object `r`.
In the first case iteration is over all `values` stored in the range
(or just the first `n`, if `n` is specified). In the second case, the iteration
is over approximately `n` ordered values, generated as follows:

First, exacltly `n` values are generated between `U` and `L`, with a
spacing determined by `r.scale`, where `U` and `L` are given by the
following table:

| `r.lower`   | `r.upper`  | `L`                 | `U`                 |
|-------------|------------|---------------------|---------------------|
| finite      | finite     | `r.lower`           | `r.upper`           |
| `-Inf`      | finite     | `r.upper - 2r.unit` | `r.upper`           |
| finite      | `Inf`      | `r.lower`           | `r.lower + 2r.unit  |
| `-Inf`      | `Inf`      | `r.origin - r.unit` | `r.origin + r.unit  |

If `r` isa a discrete range (`r isa
NumericRange{<:Any,<:Any,<:Integer}`) then the values are rounded,
with any duplicate values removed. Otherwise all the values are used
as is (and there are exacltly `n` of them).

If a random number generator `rng` is specified, then the values are
returned in random order (sampling without replacement), and otherwise
they are returned in numeric order, or in the order provided to the
range constructor, in the case of a `NominalRange`.
"""
iterator(r::ParamRange, n::Integer, rng::AbstractRNG) =
    StatsBase.sample(rng, iterator(r, n), n, replace=false)

iterator(r::NominalRange, n::Integer) =
    collect(r.values[1:min(n, length(r.values))])
iterator(r::NominalRange) = collect(r.values)
iterator(r::NominalRange, rng::AbstractRNG) =
    iterator(r, length(r.values), rng)

# nominal range, top level dispatch
function iterator(r::NumericRange{<:Bounded,T},
                  n::Int) where {T<:Real}
    L = r.lower
    U = r.upper
    return iterator(T, L, U, r.scale, n)
end

function iterator(r::NumericRange{<:LeftUnbounded,T},
                  n::Int) where {T<:Real}
    L = r.upper - 2r.unit
    U = r.upper
    return iterator(T, L, U, r.scale, n)
end

function iterator(r::NumericRange{<:RightUnbounded,T},
                  n::Int) where {T<:Real}
    L = r.lower
    U = r.lower + 2r.unit
    return iterator(T, L, U, r.scale, n)
end

function iterator(r::NumericRange{<:DoublyUnbounded,T},
                  n::Int) where {T<:Real}
    L = r.origin - r.unit
    U = r.origin + r.unit
    return iterator(T, L, U, r.scale, n)
end

# middle level

iterator(::Type{<:Real}, L, U, s, n) =
    iterator(L, U, s, n)

function iterator(I::Type{<:Integer}, L, U, s, n)
    raw = iterator(L, U, s, n)
    rounded = map(x -> round(I, x), raw)
    return unique(rounded)
end

# ground level

function iterator(L, U, s, n)
    transformed = range(transform(Scale, scale(s), L),
                stop=transform(Scale, scale(s), U),
                        length=n)
    inverse_transformed = map(transformed) do value
        inverse_transform(Scale, scale(s), value)
    end
    return inverse_transformed
end


text:
using PyPlot

L = 1.0
γ0 = 0.5
γL = -0.5
x = range(0, stop=L, length=201)
u(x, c) = ( ( L - x ) * γ0 + x * γL ) / L + (c/2) * x * ( L - x )

figure(1)
cvals = collect(6:-3:-6)
plot(x, u.(x,cvals'))
grid(true)
xlabel(L"x")
legend([latexstring("c=$(cvals[k])") for k = 1:length(cvals)])
savefig("bvp1d_const_rhs.pdf")

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
"""
Class for sparse variational Gaussian Processes

```julia
SVGP(X::AbstractArray{T1},y::AbstractArray{T2},
     kernel::Union{Kernel,AbstractVector{<:Kernel}},
     likelihood::LikelihoodType,inference::InferenceType,
     nInducingPoints::Integer;        verbose::Integer=0,Autotuning::Bool=true,
     atfrequency::Integer=1,IndependentPriors::Bool=true, OptimizeInducingPoints::Bool=false,ArrayType::UnionAll=Vector)
```

Argument list :

**Mandatory arguments**
 - `X` : input features, should be a matrix N×D where N is the number of observation and D the number of dimension
 - `y` : input labels, can be either a vector of labels for multiclass and single output or a matrix for multi-outputs (note that only one likelihood can be applied)
 - `kernel` : covariance function, can be either a single kernel or a collection of kernels for multiclass and multi-outputs models
 - `likelihood` : likelihood of the model, currently implemented : Gaussian, Student-T, Laplace, Bernoulli (with logistic link), Bayesian SVM, Multiclass (softmax or logistic-softmax) see [`Likelihood`](@ref likelihood_user)
 - `inference` : inference for the model, can be analytic, numerical or by sampling, check the model documentation to know what is available for your likelihood see the [`Compatibility table`](@ref compat_table)
 - `nInducingPoints` : number of inducing points
**Optional arguments**
 - `verbose` : How much does the model print (0:nothing, 1:very basic, 2:medium, 3:everything)
 - `Autotuning` : Flag for optimizing hyperparameters
 - `atfrequency` : Choose how many variational parameters iterations are between hyperparameters optimization
 - `IndependentPriors` : Flag for setting independent or shared parameters among latent GPs
 - `OptimizeInducingPoints` : Flag for optimizing the inducing points locations
 - `ArrayType` : Option for using different type of array for storage (allow for GPU usage)
"""
mutable struct SVGP{L<:Likelihood,I<:Inference,T<:Real,V<:AbstractVector{T}} <: AbstractGP{L,I,T,V}
    X::Matrix{T} #Feature vectors
    y::LatentArray #Output (-1,1 for classification, real for regression, matrix for multiclass)
    nSample::Int64 # Number of data points
    nDim::Int64 # Number of covariates per data point
    nFeature::Int64 # Number of features of the GP (equal to number of points)
    nLatent::Int64 # Number pf latent GPs
    IndependentPriors::Bool # Use of separate priors for each latent GP
    nPrior::Int64 # Equal to 1 or nLatent given IndependentPriors
    Z::LatentArray{Matrix{T}} #Inducing points locations
    μ::LatentArray{V}
    Σ::LatentArray{Symmetric{T,Matrix{T}}}
    η₁::LatentArray{V}
    η₂::LatentArray{Symmetric{T,Matrix{T}}}
    μ₀::LatentArray{MeanPrior{T}}
    Kmm::LatentArray{Symmetric{T,Matrix{T}}}
    invKmm::LatentArray{Symmetric{T,Matrix{T}}}
    Knm::LatentArray{Matrix{T}}
    κ::LatentArray{Matrix{T}}
    K̃::LatentArray{V}
    kernel::LatentArray{Kernel{T}}
    likelihood::Likelihood{T}
    inference::Inference{T}
    verbose::Int64
    Autotuning::Bool
    atfrequency::Int64
    OptimizeInducingPoints::Bool
    Trained::Bool
end

function SVGP(X::AbstractArray{T1},y::AbstractArray{T2},kernel::Union{Kernel,AbstractVector{<:Kernel}},
            likelihood::LikelihoodType,inference::InferenceType,
            nInducingPoints::Integer
            ;verbose::Integer=0,Autotuning::Bool=true,atfrequency::Integer=1,
            mean::Union{<:Real,AbstractVector{<:Real},MeanPrior}=ZeroMean(),
            IndependentPriors::Bool=true, OptimizeInducingPoints::Bool=false,ArrayType::UnionAll=Vector) where {T1<:Real,T2,LikelihoodType<:Likelihood,InferenceType<:Inference}

            X,y,nLatent,likelihood = check_data!(X,y,likelihood)
            @assert check_implementation(:SVGP,likelihood,inference) "The $likelihood is not compatible or implemented with the $inference"

            nPrior = IndependentPriors ? nLatent : 1
            nSample = size(X,1); nDim = size(X,2);
            kernel = [deepcopy(kernel) for _ in 1:nPrior]

            @assert nInducingPoints > 0 && nInducingPoints < nSample "The number of inducing points is incorrect (negative or bigger than number of samples)"
            Z = KMeansInducingPoints(X,nInducingPoints,nMarkov=10); Z=[deepcopy(Z) for _ in 1:nPrior]
            nFeature = nInducingPoints

            μ = LatentArray([zeros(T1,nFeature) for _ in 1:nLatent]); η₁ = deepcopy(μ);
            Σ = LatentArray([Symmetric(Matrix(Diagonal(one(T1)*I,nFeature))) for _ in 1:nLatent]);
            η₂ = -0.5*inv.(Σ);
            κ = LatentArray([zeros(T1,inference.Stochastic ? inference.nSamplesUsed : nSample, nFeature) for _ in 1:nPrior])
            Knm = deepcopy(κ)
            K̃ = LatentArray([zeros(T1,inference.Stochastic ? inference.nSamplesUsed : nSample) for _ in 1:nPrior])
            Kmm = LatentArray([similar(Σ[1]) for _ in 1:nPrior]); invKmm = similar.(Kmm)
            μ₀ = []
            if typeof(mean) <: Real
                μ₀ = [ConstantMean(mean) for _ in 1:nPrior]
            elseif typeof(mean) <: AbstractVector{<:Real}
                μ₀ = [EmpiricalMean(mean) for _ in 1:nPrior]
            else
                μ₀ = [mean for _ in 1:nPrior]
            end

            nSamplesUsed = nSample
            if inference.Stochastic
                @assert inference.nSamplesUsed > 0 && inference.nSamplesUsed < nSample "The size of mini-batch is incorrect (negative or bigger than number of samples), please set nMinibatch correctly in the inference object"
                nSamplesUsed = inference.nSamplesUsed
                opt = kernel[1].fields.variance.opt
                opt.α = opt.α*0.1
                setoptimizer!.(kernel,[copy(opt) for _ in 1:nLatent])
            end

            likelihood = init_likelihood(likelihood,inference,nLatent,nSamplesUsed)
            inference = init_inference(inference,nLatent,nFeature,nSample,nSamplesUsed)
            model = SVGP{LikelihoodType,InferenceType,T1,ArrayType{T1}}(X,y,
                    nSample, nDim, nFeature, nLatent,
                    IndependentPriors,nPrior,
                    Z,μ,Σ,η₁,η₂,
                    μ₀,Kmm,invKmm,Knm,κ,K̃,
                    kernel,likelihood,inference,
                    verbose,Autotuning,atfrequency,OptimizeInducingPoints,false)
            if isa(inference.optimizer_η₁[1],ALRSVI)
                init!(model.inference,model)
            end
            return model
end

function Base.show(io::IO,model::SVGP{<:Likelihood,<:Inference,T}) where T
    print(io,"Sparse Variational Gaussian Process with a $(model.likelihood) infered by $(model.inference) ")
end