text:
struct OneShootSimulation{T, N, Ny, Nx, Nθ}
    f::Function
    g::Function
    x0::T
    y::Vector{T}
    time_span::Tuple{Int, Int}
    # Sensitivity equations
    Jf::Function
    Jg::Function
    # x
    x
    dxdθ
    dxdx0
    x_extended
    # x_buffer
    x_buffer
    dxdθ_buffer
    dxdx0_buffer
    x_buffer_extended
    # Simulated y
    ys
    dydθ
    dydx0
    ys_extended
    # Buffer y
    y_buffer
    # Loss function
    loss
end

function OneShootSimulation{T}(f::Function, g::Function,
                               x0::T, y::Vector{T}, k0::Int, θ,
                               loss=L2DistLoss())
    # Get sizes
    N = length(y)
    Ny = length(y[1])
    Nx = length(x0)
    Nθ = length(θ)
    time_span = (k0, k0+N-1)
    # Define x0_extended
    x = deepcopy(x0)
    dxdθ = zeros(Nx, Nθ)
    dxdx0 = eye(Nx, Nx)
    x_extended = (x, dxdθ, dxdx0)
    # Define x_buffer
    x_buffer = deepcopy(x0)
    dxdθ_buffer = zeros(Nx, Nθ)
    dxdx0_buffer = eye(Nx, Nx)
    x_buffer_extended = (x_buffer, dxdθ_buffer, dxdx0_buffer)
    # define ys extended
    ys_extended = Vector{Tuple{T, Matrix{Float64},Matrix{Float64}}}(N)
    ys = Vector{T}(N)
    dydθ = Vector{Matrix{Float64}}(N)
    dydx0 = Vector{Matrix{Float64}}(N)
    for i = 1:N
        ys[i] = zeros(Ny)
        dydθ[i] = zeros(Ny, Nθ)
        dydx0[i] = eye(Ny, Nx)
        ys_extended[i] = (ys[i], dydθ[i], dydx0[i])
    end
    # Define auxiliar buffers
    buffer_Jf = zeros(Nx, Nx)
    buffer_Jg = zeros(Ny, Nx)
    y_buffer = zeros(Ny)
    # Sensitivity equations
    Jf = sensitivity_equation(f, buffer_Jf)
    Jg = sensitivity_equation(g, buffer_Jg)
    # Simulate
    simulate_space_state!(ys_extended, Jf, Jg, x_extended,
                          time_span, (θ,); x0_2=x_buffer_extended)
    OneShootSimulation{T, N, Ny, Nx, Nθ}(
        f, g, deepcopy(x0), y, time_span, Jf, Jg, x, dxdθ, dxdx0, x_extended,
        x_buffer, dxdθ_buffer, dxdx0_buffer, x_buffer_extended,
        ys, dydθ, dydx0, ys_extended, y_buffer, loss)
end

function new_simulation!{T, N, Ny, Nx, Nθ}(
        oss::OneShootSimulation{T, N, Ny, Nx, Nθ}, x0::T, θ)
    # Define x0_extended
    copy!(oss.x0, x0)
    copy!(oss.x, x0)
    copy!(oss.dxdθ, zeros(Nx, Nθ))
    copy!(oss.dxdx0, eye(Nx, Nx))
    # Define x_buffer
    copy!(oss.x_buffer, x0)
    copy!(oss.dxdθ_buffer, zeros(Nx, Nθ))
    copy!(oss.dxdx0_buffer, eye(Nx, Nx))
    # Simulate
    simulate_space_state!(oss.ys_extended, oss.Jf, oss.Jg, oss.x_extended,
                          oss.time_span, (θ,);
                          x0_2=oss.x_buffer_extended)
    return oss
end

function cost_function{T, N, Ny, Nx, Nθ}(
       oss::OneShootSimulation{T, N, Ny, Nx, Nθ})
    cost = 0
    for i = 1:N
        cost += value(oss.loss, oss.y[i], oss.ys[i], AvgMode.Sum())
    end
    return cost
end

function gradient!{T, N, Ny, Nx, Nθ}(grad::Vector{Float64},
        oss::OneShootSimulation{T, N, Ny, Nx, Nθ},
        variable="θ", accumulat=false)
    if !accumulat
        fill!(grad, 0)
    end
    if variable == "θ"
        J = oss.dydθ
    else
        J = oss.dydx0
    end

    for i = 1:N
        deriv!(oss.y_buffer, oss.loss, oss.y[i], oss.ys[i])
        Base.LinAlg.BLAS.gemv!('T', 1.0, J[i], oss.y_buffer, 1.0, grad)
    end
    return grad
end

function hessian_approx!{T, N, Ny, Nx, Nθ}(hessp::Vector{Float64},
        oss::OneShootSimulation{T, N, Ny, Nx, Nθ}, p::Vector{Float64},
        variable="θ", accumulat=false)
    if !accumulat
        fill!(hessp, 0)
    end
    if variable == "θ"
        J = oss.dydθ
    else
        J = oss.dydx0
    end

    for i = 1:N
        A_mul_B!(oss.y_buffer, J[i], p)
        oss.y_buffer .*= deriv2.(oss.loss, oss.y[i], oss.ys[i])
        Base.LinAlg.BLAS.gemv!('T', 1.0, J[i], oss.y_buffer, 1.0, hessp)
    end
    return hessp
end

# Get functions (usefull when using julia native parallel functions)
get_x0(oss::OneShootSimulation) = oss.x0
get_x(oss::OneShootSimulation) = oss.x
get_dxdx0(oss::OneShootSimulation) = oss.dxdx0
get_dxdθ(oss::OneShootSimulation) = oss.dxdθ


text:
using FeedbackParticleFilters, Statistics, StatsBase

println("Testing WeightedParticleEnsemble.jl:")

@testset "WeightedParticleEnsemble.jl" begin

    print("  inner constructor for WeightedParticleEnsemble")
    N    = 6
    d    = 2
    x_pf = [-2.422480820086937 -0.592332203303167 -2.017301296096984 -1.5151245392598531 0.02565906919199346  0.15161614796874012;
            -1.4986273707574955 0.3549802657063476 0.3625127047102   -0.166475807288935  0.05739520762738308 -0.18730971131161062]
    w_pf = StatsBase.ProbabilityWeights([0.1, 0.1, 0.1, 0.1, 0.1, 0.5])
    ens  = WeightedParticleEnsemble(x_pf, w_pf)
    print("-")
    @test ens.positions == x_pf
    print("-")
    @test ens.weights   == w_pf
    println("DONE")

    print("  method dim")
    print("-")
    @test dim(ens) == (d+1) * N
    println("DONE")

    print("  method particle_dim")
    print("-")
    @test particle_dim(ens) == d
    println("DONE")

    print("  method no_of_particles")
    print("-")
    @test no_of_particles(ens) == N
    println("DONE")

    print("  method eff_no_of_particles")
    print("-")
    @test eff_no_of_particles(ens) == 3.3333333333333335
    println("DONE")

    print("  method get_pos")
    for i in 1:N
        print("-")
        @test get_pos(ens, i) == x_pf[:, i]
    end
    println("DONE")

    print("  method get_weight")
    for i in 1:N
        print("-")
        @test get_weight(ens, i) == w_pf[i]
    end
    println("DONE")

    print("  method sum_of_weights")
    print("-")
    @test sum_of_weights(ens) == 1.0
    println("DONE")

    print("  method mean")
    print("-")
    @test size(mean(ens)) == (d, 1)
    print("-")
    @test maximum(abs.(mean(ens)[:,1] - [-0.5763499049711248; -0.18267635565605528])) < 1e-15
    println("DONE")

    print("  method cov")
    print("-")
    @test size(cov(ens)) == (d, d)
    print("-")
    @test maximum(abs.(cov(ens) - [0.9378181191314968 0.17476858717998095; 0.17476858717998095 0.23760369813722895])) < 1e-15
    println("DONE")

    print("  method var")
    print("-")
    @test size(var(ens)) == (d, 1)
    print("-")
    @test maximum(abs.(var(ens)[:,1] - [0.9378181191314968; 0.23760369813722895])) < 1e-15
    println("DONE")

    print("  method resample!")
    x_pf_old = deepcopy(x_pf)
    resample!(ens)
    for i in 1:N
        print("-")
        @test get_weight(ens, i) == float(1/N)
    end
    for i in 1:N, j in 1:d
        print("-")
        @test get_pos(ens, i)[j] in x_pf_old
    end
    print("-")
    @test abs(eff_no_of_particles(ens) - N) < 1e-3
    println("DONE")

end; #WeightedParticleEnsemble.jl

ext: jl

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
function Base.:~(arg1::AbstractNumber)
    tmp = Base.:~(number(arg1))
    like(arg1, tmp)
end

function Base.:conj(arg1::AbstractNumber)
    tmp = Base.:conj(number(arg1))
    like(arg1, tmp)
end

function Base.:abs(arg1::AbstractNumber)
    tmp = Base.:abs(number(arg1))
    like(arg1, tmp)
end

function Base.:sin(arg1::AbstractNumber)
    tmp = Base.:sin(number(arg1))
    like(arg1, tmp)
end

function Base.:cos(arg1::AbstractNumber)
    tmp = Base.:cos(number(arg1))
    like(arg1, tmp)
end

function Base.:tan(arg1::AbstractNumber)
    tmp = Base.:tan(number(arg1))
    like(arg1, tmp)
end

function Base.:sinh(arg1::AbstractNumber)
    tmp = Base.:sinh(number(arg1))
    like(arg1, tmp)
end

function Base.:cosh(arg1::AbstractNumber)
    tmp = Base.:cosh(number(arg1))
    like(arg1, tmp)
end

function Base.:tanh(arg1::AbstractNumber)
    tmp = Base.:tanh(number(arg1))
    like(arg1, tmp)
end

function Base.:asin(arg1::AbstractNumber)
    tmp = Base.:asin(number(arg1))
    like(arg1, tmp)
end

function Base.:acos(arg1::AbstractNumber)
    tmp = Base.:acos(number(arg1))
    like(arg1, tmp)
end

function Base.:atan(arg1::AbstractNumber)
    tmp = Base.:atan(number(arg1))
    like(arg1, tmp)
end

function Base.:atan(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:atan(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:atan(a::AbstractNumber, b::Number)
    tmp = Base.:atan(number(a), b)
    like(a, tmp)
end

function Base.:atan(a::Number, b::AbstractNumber)
    tmp = Base.:atan(a, number(b))
    like(b, tmp)
end

function Base.:asinh(arg1::AbstractNumber)
    tmp = Base.:asinh(number(arg1))
    like(arg1, tmp)
end

function Base.:acosh(arg1::AbstractNumber)
    tmp = Base.:acosh(number(arg1))
    like(arg1, tmp)
end

function Base.:atanh(arg1::AbstractNumber)
    tmp = Base.:atanh(number(arg1))
    like(arg1, tmp)
end

function Base.:sec(arg1::AbstractNumber)
    tmp = Base.:sec(number(arg1))
    like(arg1, tmp)
end

function Base.:csc(arg1::AbstractNumber)
    tmp = Base.:csc(number(arg1))
    like(arg1, tmp)
end

function Base.:cot(arg1::AbstractNumber)
    tmp = Base.:cot(number(arg1))
    like(arg1, tmp)
end

function Base.:asec(arg1::AbstractNumber)
    tmp = Base.:asec(number(arg1))
    like(arg1, tmp)
end

function Base.:acsc(arg1::AbstractNumber)
    tmp = Base.:acsc(number(arg1))
    like(arg1, tmp)
end

function Base.:acot(arg1::AbstractNumber)
    tmp = Base.:acot(number(arg1))
    like(arg1, tmp)
end

function Base.:sech(arg1::AbstractNumber)
    tmp = Base.:sech(number(arg1))
    like(arg1, tmp)
end

function Base.:csch(arg1::AbstractNumber)
    tmp = Base.:csch(number(arg1))
    like(arg1, tmp)
end

function Base.:coth(arg1::AbstractNumber)
    tmp = Base.:coth(number(arg1))
    like(arg1, tmp)
end

function Base.:asech(arg1::AbstractNumber)
    tmp = Base.:asech(number(arg1))
    like(arg1, tmp)
end

function Base.:acsch(arg1::AbstractNumber)
    tmp = Base.:acsch(number(arg1))
    like(arg1, tmp)
end

function Base.:acoth(arg1::AbstractNumber)
    tmp = Base.:acoth(number(arg1))
    like(arg1, tmp)
end

function Base.:sinc(arg1::AbstractNumber)
    tmp = Base.:sinc(number(arg1))
    like(arg1, tmp)
end

function Base.:cosc(arg1::AbstractNumber)
    tmp = Base.:cosc(number(arg1))
    like(arg1, tmp)
end

function Base.:cosd(arg1::AbstractNumber)
    tmp = Base.:cosd(number(arg1))
    like(arg1, tmp)
end

function Base.:cotd(arg1::AbstractNumber)
    tmp = Base.:cotd(number(arg1))
    like(arg1, tmp)
end

function Base.:cscd(arg1::AbstractNumber)
    tmp = Base.:cscd(number(arg1))
    like(arg1, tmp)
end

function Base.:secd(arg1::AbstractNumber)
    tmp = Base.:secd(number(arg1))
    like(arg1, tmp)
end

function Base.:sind(arg1::AbstractNumber)
    tmp = Base.:sind(number(arg1))
    like(arg1, tmp)
end

function Base.:tand(arg1::AbstractNumber)
    tmp = Base.:tand(number(arg1))
    like(arg1, tmp)
end

function Base.:acosd(arg1::AbstractNumber)
    tmp = Base.:acosd(number(arg1))
    like(arg1, tmp)
end

function Base.:acotd(arg1::AbstractNumber)
    tmp = Base.:acotd(number(arg1))
    like(arg1, tmp)
end

function Base.:acscd(arg1::AbstractNumber)
    tmp = Base.:acscd(number(arg1))
    like(arg1, tmp)
end

function Base.:asecd(arg1::AbstractNumber)
    tmp = Base.:asecd(number(arg1))
    like(arg1, tmp)
end

function Base.:asind(arg1::AbstractNumber)
    tmp = Base.:asind(number(arg1))
    like(arg1, tmp)
end

function Base.:atand(arg1::AbstractNumber)
    tmp = Base.:atand(number(arg1))
    like(arg1, tmp)
end

function Base.:atand(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:atand(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:atand(a::AbstractNumber, b::Number)
    tmp = Base.:atand(number(a), b)
    like(a, tmp)
end

function Base.:atand(a::Number, b::AbstractNumber)
    tmp = Base.:atand(a, number(b))
    like(b, tmp)
end

function Base.:rad2deg(arg1::AbstractNumber)
    tmp = Base.:rad2deg(number(arg1))
    like(arg1, tmp)
end

function Base.:deg2rad(arg1::AbstractNumber)
    tmp = Base.:deg2rad(number(arg1))
    like(arg1, tmp)
end

function Base.:log(arg1::AbstractNumber)
    tmp = Base.:log(number(arg1))
    like(arg1, tmp)
end

function Base.:log(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:log(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:log(a::AbstractNumber, b::Number)
    tmp = Base.:log(number(a), b)
    like(a, tmp)
end

function Base.:log(a::Number, b::AbstractNumber)
    tmp = Base.:log(a, number(b))
    like(b, tmp)
end

function Base.:log2(arg1::AbstractNumber)
    tmp = Base.:log2(number(arg1))
    like(arg1, tmp)
end

function Base.:log10(arg1::AbstractNumber)
    tmp = Base.:log10(number(arg1))
    like(arg1, tmp)
end

function Base.:log1p(arg1::AbstractNumber)
    tmp = Base.:log1p(number(arg1))
    like(arg1, tmp)
end

function Base.:exponent(arg1::AbstractNumber)
    tmp = Base.:exponent(number(arg1))
    like(arg1, tmp)
end

function Base.:exp(arg1::AbstractNumber)
    tmp = Base.:exp(number(arg1))
    like(arg1, tmp)
end

function Base.:exp2(arg1::AbstractNumber)
    tmp = Base.:exp2(number(arg1))
    like(arg1, tmp)
end

function Base.:expm1(arg1::AbstractNumber)
    tmp = Base.:expm1(number(arg1))
    like(arg1, tmp)
end

function Base.:cbrt(arg1::AbstractNumber)
    tmp = Base.:cbrt(number(arg1))
    like(arg1, tmp)
end

function Base.:sqrt(arg1::AbstractNumber)
    tmp = Base.:sqrt(number(arg1))
    like(arg1, tmp)
end

function Base.:ceil(arg1::AbstractNumber)
    tmp = Base.:ceil(number(arg1))
    like(arg1, tmp)
end

function Base.:floor(arg1::AbstractNumber)
    tmp = Base.:floor(number(arg1))
    like(arg1, tmp)
end

function Base.:trunc(arg1::AbstractNumber)
    tmp = Base.:trunc(number(arg1))
    like(arg1, tmp)
end

function Base.:round(arg1::AbstractNumber)
    tmp = Base.:round(number(arg1))
    like(arg1, tmp)
end

function Base.:significand(arg1::AbstractNumber)
    tmp = Base.:significand(number(arg1))
    like(arg1, tmp)
end

function Base.:frexp(arg1::AbstractNumber)
    tmp = Base.:frexp(number(arg1))
    like(arg1, tmp)
end

function Base.:ldexp(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:ldexp(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:ldexp(a::AbstractNumber, b::Number)
    tmp = Base.:ldexp(number(a), b)
    like(a, tmp)
end

function Base.:ldexp(a::Number, b::AbstractNumber)
    tmp = Base.:ldexp(a, number(b))
    like(b, tmp)
end

function Base.:modf(arg1::AbstractNumber)
    tmp = Base.:modf(number(arg1))
    like(arg1, tmp)
end

function Base.:real(arg1::AbstractNumber)
    tmp = Base.:real(number(arg1))
    like(arg1, tmp)
end

function Base.:imag(arg1::AbstractNumber)
    tmp = Base.:imag(number(arg1))
    like(arg1, tmp)
end

function Base.:!(arg1::AbstractNumber)
    tmp = Base.:!(number(arg1))
    like(arg1, tmp)
end

function Base.:identity(arg1::AbstractNumber)
    tmp = Base.:identity(number(arg1))
    like(arg1, tmp)
end

function Base.:zero(arg1::AbstractNumber)
    tmp = Base.:zero(number(arg1))
    like(arg1, tmp)
end

function Base.:one(arg1::AbstractNumber)
    tmp = Base.:one(number(arg1))
    like(arg1, tmp)
end

function Base.:<<(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:<<(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:<<(a::AbstractNumber, b::Number)
    tmp = Base.:<<(number(a), b)
    like(a, tmp)
end

function Base.:<<(a::Number, b::AbstractNumber)
    tmp = Base.:<<(a, number(b))
    like(b, tmp)
end

function Base.:>>(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:>>(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:>>(a::AbstractNumber, b::Number)
    tmp = Base.:>>(number(a), b)
    like(a, tmp)
end

function Base.:>>(a::Number, b::AbstractNumber)
    tmp = Base.:>>(a, number(b))
    like(b, tmp)
end

function Base.:abs2(arg1::AbstractNumber)
    tmp = Base.:abs2(number(arg1))
    like(arg1, tmp)
end

function Base.:sign(arg1::AbstractNumber)
    tmp = Base.:sign(number(arg1))
    like(arg1, tmp)
end

function Base.:sinpi(arg1::AbstractNumber)
    tmp = Base.:sinpi(number(arg1))
    like(arg1, tmp)
end

function Base.:cospi(arg1::AbstractNumber)
    tmp = Base.:cospi(number(arg1))
    like(arg1, tmp)
end

function Base.:exp10(arg1::AbstractNumber)
    tmp = Base.:exp10(number(arg1))
    like(arg1, tmp)
end

function Base.:iseven(arg1::AbstractNumber)
    tmp = Base.:iseven(number(arg1))
    tmp
end

function Base.:ispow2(arg1::AbstractNumber)
    tmp = Base.:ispow2(number(arg1))
    tmp
end

function Base.:isfinite(arg1::AbstractNumber)
    tmp = Base.:isfinite(number(arg1))
    tmp
end

function Base.:isinf(arg1::AbstractNumber)
    tmp = Base.:isinf(number(arg1))
    tmp
end

function Base.:isodd(arg1::AbstractNumber)
    tmp = Base.:isodd(number(arg1))
    tmp
end

function Base.:isinteger(arg1::AbstractNumber)
    tmp = Base.:isinteger(number(arg1))
    tmp
end

function Base.:isreal(arg1::AbstractNumber)
    tmp = Base.:isreal(number(arg1))
    tmp
end

function Base.:isnan(arg1::AbstractNumber)
    tmp = Base.:isnan(number(arg1))
    tmp
end

function Base.:isempty(arg1::AbstractNumber)
    tmp = Base.:isempty(number(arg1))
    tmp
end

function Base.:iszero(arg1::AbstractNumber)
    tmp = Base.:iszero(number(arg1))
    tmp
end

function Base.:transpose(arg1::AbstractNumber)
    tmp = Base.:transpose(number(arg1))
    like(arg1, tmp)
end

function Base.:copysign(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:copysign(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:copysign(a::AbstractNumber, b::Number)
    tmp = Base.:copysign(number(a), b)
    like(a, tmp)
end

function Base.:copysign(a::Number, b::AbstractNumber)
    tmp = Base.:copysign(a, number(b))
    like(b, tmp)
end

function Base.:flipsign(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:flipsign(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:flipsign(a::AbstractNumber, b::Number)
    tmp = Base.:flipsign(number(a), b)
    like(a, tmp)
end

function Base.:flipsign(a::Number, b::AbstractNumber)
    tmp = Base.:flipsign(a, number(b))
    like(b, tmp)
end

function Base.:signbit(arg1::AbstractNumber)
    tmp = Base.:signbit(number(arg1))
    like(arg1, tmp)
end

function Base.:+(arg1::AbstractNumber)
    tmp = Base.:+(number(arg1))
    like(arg1, tmp)
end

function Base.:+(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:+(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:+(a::AbstractNumber, b::Number)
    tmp = Base.:+(number(a), b)
    like(a, tmp)
end

function Base.:+(a::Number, b::AbstractNumber)
    tmp = Base.:+(a, number(b))
    like(b, tmp)
end

function Base.:-(arg1::AbstractNumber)
    tmp = Base.:-(number(arg1))
    like(arg1, tmp)
end

function Base.:-(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:-(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:-(a::AbstractNumber, b::Number)
    tmp = Base.:-(number(a), b)
    like(a, tmp)
end

function Base.:-(a::Number, b::AbstractNumber)
    tmp = Base.:-(a, number(b))
    like(b, tmp)
end

function Base.:*(arg1::AbstractNumber)
    tmp = Base.:*(number(arg1))
    like(arg1, tmp)
end

function Base.:*(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:*(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:*(a::AbstractNumber, b::Number)
    tmp = Base.:*(number(a), b)
    like(a, tmp)
end

function Base.:*(a::Number, b::AbstractNumber)
    tmp = Base.:*(a, number(b))
    like(b, tmp)
end

function Base.:/(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:/(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:/(a::AbstractNumber, b::Number)
    tmp = Base.:/(number(a), b)
    like(a, tmp)
end

function Base.:/(a::Number, b::AbstractNumber)
    tmp = Base.:/(a, number(b))
    like(b, tmp)
end

function Base.:\(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:\(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:\(a::AbstractNumber, b::Number)
    tmp = Base.:\(number(a), b)
    like(a, tmp)
end

function Base.:\(a::Number, b::AbstractNumber)
    tmp = Base.:\(a, number(b))
    like(b, tmp)
end

function Base.:^(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:^(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:^(a::AbstractNumber, b::Number)
    tmp = Base.:^(number(a), b)
    like(a, tmp)
end

function Base.:^(a::Number, b::AbstractNumber)
    tmp = Base.:^(a, number(b))
    like(b, tmp)
end

function Base.:(==)(arg1::AbstractNumber)
    tmp = Base.:(==)(number(arg1))
    tmp
end

function Base.:(==)(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:(==)(number(arg1), number(arg2))
    tmp
end

function Base.:(==)(a::AbstractNumber, b::Number)
    tmp = Base.:(==)(number(a), b)
    tmp
end

function Base.:(==)(a::Number, b::AbstractNumber)
    tmp = Base.:(==)(a, number(b))
    tmp
end

function Base.:!=(arg1::AbstractNumber)
    tmp = Base.:!=(number(arg1))
    tmp
end

function Base.:!=(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:!=(number(arg1), number(arg2))
    tmp
end

function Base.:!=(a::AbstractNumber, b::Number)
    tmp = Base.:!=(number(a), b)
    tmp
end

function Base.:!=(a::Number, b::AbstractNumber)
    tmp = Base.:!=(a, number(b))
    tmp
end

function Base.:<(arg1::AbstractNumber)
    tmp = Base.:<(number(arg1))
    tmp
end

function Base.:<(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:<(number(arg1), number(arg2))
    tmp
end

function Base.:<(a::AbstractNumber, b::Number)
    tmp = Base.:<(number(a), b)
    tmp
end

function Base.:<(a::Number, b::AbstractNumber)
    tmp = Base.:<(a, number(b))
    tmp
end

function Base.:<=(arg1::AbstractNumber)
    tmp = Base.:<=(number(arg1))
    tmp
end

function Base.:<=(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:<=(number(arg1), number(arg2))
    tmp
end

function Base.:<=(a::AbstractNumber, b::Number)
    tmp = Base.:<=(number(a), b)
    tmp
end

function Base.:<=(a::Number, b::AbstractNumber)
    tmp = Base.:<=(a, number(b))
    tmp
end

function Base.:>(arg1::AbstractNumber)
    tmp = Base.:>(number(arg1))
    tmp
end

function Base.:>(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:>(number(arg1), number(arg2))
    tmp
end

function Base.:>(a::AbstractNumber, b::Number)
    tmp = Base.:>(number(a), b)
    tmp
end

function Base.:>(a::Number, b::AbstractNumber)
    tmp = Base.:>(a, number(b))
    tmp
end

function Base.:>=(arg1::AbstractNumber)
    tmp = Base.:>=(number(arg1))
    tmp
end

function Base.:>=(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:>=(number(arg1), number(arg2))
    tmp
end

function Base.:>=(a::AbstractNumber, b::Number)
    tmp = Base.:>=(number(a), b)
    tmp
end

function Base.:>=(a::Number, b::AbstractNumber)
    tmp = Base.:>=(a, number(b))
    tmp
end

function Base.:≈(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:≈(number(arg1), number(arg2))
    tmp
end

function Base.:≈(a::AbstractNumber, b::Number)
    tmp = Base.:≈(number(a), b)
    tmp
end

function Base.:≈(a::Number, b::AbstractNumber)
    tmp = Base.:≈(a, number(b))
    tmp
end

function Base.:min(arg1::AbstractNumber)
    tmp = Base.:min(number(arg1))
    like(arg1, tmp)
end

function Base.:min(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:min(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:min(a::AbstractNumber, b::Number)
    tmp = Base.:min(number(a), b)
    like(a, tmp)
end

function Base.:min(a::Number, b::AbstractNumber)
    tmp = Base.:min(a, number(b))
    like(b, tmp)
end

function Base.:max(arg1::AbstractNumber)
    tmp = Base.:max(number(arg1))
    like(arg1, tmp)
end

function Base.:max(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:max(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:max(a::AbstractNumber, b::Number)
    tmp = Base.:max(number(a), b)
    like(a, tmp)
end

function Base.:max(a::Number, b::AbstractNumber)
    tmp = Base.:max(a, number(b))
    like(b, tmp)
end

function Base.:div(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:div(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:div(a::AbstractNumber, b::Number)
    tmp = Base.:div(number(a), b)
    like(a, tmp)
end

function Base.:div(a::Number, b::AbstractNumber)
    tmp = Base.:div(a, number(b))
    like(b, tmp)
end

function Base.:fld(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:fld(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:fld(a::AbstractNumber, b::Number)
    tmp = Base.:fld(number(a), b)
    like(a, tmp)
end

function Base.:fld(a::Number, b::AbstractNumber)
    tmp = Base.:fld(a, number(b))
    like(b, tmp)
end

function Base.:rem(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:rem(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:rem(a::AbstractNumber, b::Number)
    tmp = Base.:rem(number(a), b)
    like(a, tmp)
end

function Base.:rem(a::Number, b::AbstractNumber)
    tmp = Base.:rem(a, number(b))
    like(b, tmp)
end

function Base.:mod(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:mod(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:mod(a::AbstractNumber, b::Number)
    tmp = Base.:mod(number(a), b)
    like(a, tmp)
end

function Base.:mod(a::Number, b::AbstractNumber)
    tmp = Base.:mod(a, number(b))
    like(b, tmp)
end

function Base.:mod1(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:mod1(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:mod1(a::AbstractNumber, b::Number)
    tmp = Base.:mod1(number(a), b)
    like(a, tmp)
end

function Base.:mod1(a::Number, b::AbstractNumber)
    tmp = Base.:mod1(a, number(b))
    like(b, tmp)
end

function Base.:cmp(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:cmp(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:cmp(a::AbstractNumber, b::Number)
    tmp = Base.:cmp(number(a), b)
    like(a, tmp)
end

function Base.:cmp(a::Number, b::AbstractNumber)
    tmp = Base.:cmp(a, number(b))
    like(b, tmp)
end

function Base.:&(arg1::AbstractNumber)
    tmp = Base.:&(number(arg1))
    like(arg1, tmp)
end

function Base.:&(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:&(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:&(a::AbstractNumber, b::Number)
    tmp = Base.:&(number(a), b)
    like(a, tmp)
end

function Base.:&(a::Number, b::AbstractNumber)
    tmp = Base.:&(a, number(b))
    like(b, tmp)
end

function Base.:|(arg1::AbstractNumber)
    tmp = Base.:|(number(arg1))
    like(arg1, tmp)
end

function Base.:|(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:|(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:|(a::AbstractNumber, b::Number)
    tmp = Base.:|(number(a), b)
    like(a, tmp)
end

function Base.:|(a::Number, b::AbstractNumber)
    tmp = Base.:|(a, number(b))
    like(b, tmp)
end

function Base.:xor(arg1::AbstractNumber)
    tmp = Base.:xor(number(arg1))
    like(arg1, tmp)
end

function Base.:xor(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = Base.:xor(number(arg1), number(arg2))
    like(arg1, tmp)
end

function Base.:xor(a::AbstractNumber, b::Number)
    tmp = Base.:xor(number(a), b)
    like(a, tmp)
end

function Base.:xor(a::Number, b::AbstractNumber)
    tmp = Base.:xor(a, number(b))
    like(b, tmp)
end

function Base.:clamp(arg1::AbstractNumber, arg2::AbstractNumber, arg3::AbstractNumber)
    tmp = Base.:clamp(number(arg1), number(arg2), number(arg3))
    like(arg1, tmp)
end

function SpecialFunctions.:airyai(arg1::AbstractNumber)
    tmp = SpecialFunctions.:airyai(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:airyaiprime(arg1::AbstractNumber)
    tmp = SpecialFunctions.:airyaiprime(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:airybi(arg1::AbstractNumber)
    tmp = SpecialFunctions.:airybi(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:airybiprime(arg1::AbstractNumber)
    tmp = SpecialFunctions.:airybiprime(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:airyaix(arg1::AbstractNumber)
    tmp = SpecialFunctions.:airyaix(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:airyaiprimex(arg1::AbstractNumber)
    tmp = SpecialFunctions.:airyaiprimex(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:airybix(arg1::AbstractNumber)
    tmp = SpecialFunctions.:airybix(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:airybiprimex(arg1::AbstractNumber)
    tmp = SpecialFunctions.:airybiprimex(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:besselh(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:besselh(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:besselh(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:besselh(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:besselh(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:besselh(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:besselh(arg1::AbstractNumber, arg2::AbstractNumber, arg3::AbstractNumber)
    tmp = SpecialFunctions.:besselh(number(arg1), number(arg2), number(arg3))
    like(arg1, tmp)
end

function SpecialFunctions.:besselhx(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:besselhx(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:besselhx(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:besselhx(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:besselhx(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:besselhx(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:besselhx(arg1::AbstractNumber, arg2::AbstractNumber, arg3::AbstractNumber)
    tmp = SpecialFunctions.:besselhx(number(arg1), number(arg2), number(arg3))
    like(arg1, tmp)
end

function SpecialFunctions.:besseli(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:besseli(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:besseli(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:besseli(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:besseli(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:besseli(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:besselix(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:besselix(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:besselix(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:besselix(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:besselix(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:besselix(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:besselj(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:besselj(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:besselj(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:besselj(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:besselj(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:besselj(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:besselj0(arg1::AbstractNumber)
    tmp = SpecialFunctions.:besselj0(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:besselj1(arg1::AbstractNumber)
    tmp = SpecialFunctions.:besselj1(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:besseljx(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:besseljx(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:besseljx(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:besseljx(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:besseljx(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:besseljx(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:besselk(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:besselk(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:besselk(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:besselk(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:besselk(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:besselk(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:besselkx(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:besselkx(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:besselkx(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:besselkx(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:besselkx(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:besselkx(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:bessely(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:bessely(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:bessely(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:bessely(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:bessely(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:bessely(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:bessely0(arg1::AbstractNumber)
    tmp = SpecialFunctions.:bessely0(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:bessely1(arg1::AbstractNumber)
    tmp = SpecialFunctions.:bessely1(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:besselyx(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:besselyx(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:besselyx(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:besselyx(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:besselyx(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:besselyx(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:dawson(arg1::AbstractNumber)
    tmp = SpecialFunctions.:dawson(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:erf(arg1::AbstractNumber)
    tmp = SpecialFunctions.:erf(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:erfc(arg1::AbstractNumber)
    tmp = SpecialFunctions.:erfc(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:erfcinv(arg1::AbstractNumber)
    tmp = SpecialFunctions.:erfcinv(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:erfcx(arg1::AbstractNumber)
    tmp = SpecialFunctions.:erfcx(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:erfi(arg1::AbstractNumber)
    tmp = SpecialFunctions.:erfi(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:erfinv(arg1::AbstractNumber)
    tmp = SpecialFunctions.:erfinv(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:eta(arg1::AbstractNumber)
    tmp = SpecialFunctions.:eta(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:digamma(arg1::AbstractNumber)
    tmp = SpecialFunctions.:digamma(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:invdigamma(arg1::AbstractNumber)
    tmp = SpecialFunctions.:invdigamma(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:polygamma(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:polygamma(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:polygamma(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:polygamma(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:polygamma(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:polygamma(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:polygamma(a::Integer, b::AbstractNumber)
    tmp = SpecialFunctions.:polygamma(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:trigamma(arg1::AbstractNumber)
    tmp = SpecialFunctions.:trigamma(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:hankelh1(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:hankelh1(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:hankelh1(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:hankelh1(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:hankelh1(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:hankelh1(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:hankelh1x(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:hankelh1x(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:hankelh1x(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:hankelh1x(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:hankelh1x(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:hankelh1x(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:hankelh2(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:hankelh2(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:hankelh2(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:hankelh2(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:hankelh2(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:hankelh2(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:hankelh2x(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:hankelh2x(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:hankelh2x(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:hankelh2x(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:hankelh2x(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:hankelh2x(a, number(b))
    like(b, tmp)
end

function SpecialFunctions.:zeta(arg1::AbstractNumber)
    tmp = SpecialFunctions.:zeta(number(arg1))
    like(arg1, tmp)
end

function SpecialFunctions.:zeta(arg1::AbstractNumber, arg2::AbstractNumber)
    tmp = SpecialFunctions.:zeta(number(arg1), number(arg2))
    like(arg1, tmp)
end

function SpecialFunctions.:zeta(a::AbstractNumber, b::Number)
    tmp = SpecialFunctions.:zeta(number(a), b)
    like(a, tmp)
end

function SpecialFunctions.:zeta(a::Number, b::AbstractNumber)
    tmp = SpecialFunctions.:zeta(a, number(b))
    like(b, tmp)
end


text:
using InterleavedImages, ImageMetadata, AxisArrays, Test

A1 = rand(3,3,3,5)
A2 = rand(3,3,3,5)
ss = InterleavedImage(A1,A2,4)

@test @inferred(ss[1,3,2,3]) == A1[1,3,2,2]

@test all(ss[:,:,:,1].==A1[:,:,:,1])
@test all(ss[:,:,:,2].==A2[:,:,:,1])

@test all(ss[:,:,:,1:2].==cat(A1[:,:,:,1], A2[:,:,:,1], dims=4))

#setindex!
ss[:,:,:,3] = zeros(3,3,3)
@test all(A1[:,:,:,2].==zeros(3,3,3))

#ImageMeta
meta1 = ImageMeta(A1)
meta2 = ImageMeta(A2)
ssm = InterleavedImage(meta1,meta2,4)

@test isa(ssm, ImageMeta)
@test all(ssm[:,:,:,1].==A1[:,:,:,1])
@test all(ssm[:,:,:,2].==A2[:,:,:,1])

#AxisArray
aa1 = AxisArray(A1)
aa2 = AxisArray(A2)
ssaa = InterleavedImage(aa1,aa2,4)
@test isa(ssaa, AxisArray)
@test all(ssaa[:,:,:,1].==A1[:,:,:,1])
@test all(ssaa[:,:,:,2].==A2[:,:,:,1])


text:
#****************************************************************************
#                     Molecular Dynamics Potentials (MDP)
#                            CESMIX-MIT Project
#
#  Contributing authors: Ngoc-Cuong Nguyen (cuongng@mit.edu, exapde@gmail.com)
#****************************************************************************

__precompile__()

module Integration

using Revise

export VelocityVerlet

include("latticestruct.jl");
include("regionstruct.jl");
include("domainstruct.jl");
include("groupstruct.jl");
include("integrationlib.jl");
include("velocityverlet.jl");

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Flase

world = World(
    v0 = 1.,
    n_dogs = 120,
    boxsize = 10.0,
    motion = BrownianMotion(
        noise = 0.5,
        friction = 1.0
        )
    )
simulation = InfiniteSimulation(;
    dt = 0.05,
    world = world,
    plotter = UnicodePlotter()
    )
Flase.runsim( simulation )


text:
decompose_trace(trace::Evolutionary.OptimizationTrace) = last(trace)

function Evolutionary.trace!(record::Dict{String,Any}, objfun, state, population, method::Evolutionary.AbstractOptimizer, options)
    record["x"] = population
end

function __map_optimizer_args(prob::OptimizationProblem, opt::Evolutionary.AbstractOptimizer;
    cb=nothing,
    maxiters::Union{Number, Nothing}=nothing,
    maxtime::Union{Number, Nothing}=nothing,
    abstol::Union{Number, Nothing}=nothing,
    reltol::Union{Number, Nothing}=nothing,
    kwargs...)

    mapped_args = (;)

    mapped_args = (; mapped_args..., kwargs...)

    if !isnothing(cb)
        mapped_args = (; mapped_args..., callback = cb)
    end

    if !isnothing(maxiters)
        mapped_args = (; mapped_args..., iterations=maxiters)
    end

    if !isnothing(maxtime)
        mapped_args = (; mapped_args..., time_limit=maxtime)
    end

    if !isnothing(abstol)
        mapped_args = (; mapped_args..., abstol=abstol)
    end

    if !isnothing(reltol)
        mapped_args = (; mapped_args..., reltol=reltol)
    end

    return Evolutionary.Options(;mapped_args...)
end

function __solve(prob::OptimizationProblem, opt::Evolutionary.AbstractOptimizer, data = DEFAULT_DATA;
                 cb = (args...) -> (false),
                 maxiters::Union{Number, Nothing} = nothing,
                 maxtime::Union{Number, Nothing} = nothing,
                 abstol::Union{Number, Nothing}=nothing,
                 reltol::Union{Number, Nothing}=nothing,
                 progress = false, kwargs...)
    local x, cur, state

    if data != DEFAULT_DATA
        maxiters = length(data)
    end

    cur, state = iterate(data)

    function _cb(trace)
        cb_call = cb(decompose_trace(trace).metadata["x"],trace.value...)
        if !(typeof(cb_call) <: Bool)
            error("The callback should return a boolean `halt` for whether to stop the optimization process.")
        end
        cur, state = iterate(data, state)
        cb_call
    end

    maxiters = _check_and_convert_maxiters(maxiters)
    maxtime = _check_and_convert_maxtime(maxtime)

    _loss = function(θ)
        x = prob.f(θ, prob.p, cur...)
        return first(x)
    end

    opt_args = _map_optimizer_args(prob, opt, cb=_cb, maxiters=maxiters, maxtime=maxtime,abstol=abstol, reltol=reltol; kwargs...)

    t0 = time()
    if isnothing(prob.ub) | isnothing(prob.ub)
        opt_res = Evolutionary.optimize(_loss, prob.u0, opt, opt_args)
    else
        cons = Evolutionary.BoxConstraints(prob.lb, prob.ub)
        opt_res = Evolutionary.optimize(_loss, cons, prob.u0, opt, opt_args)
    end
    t1 = time()
    opt_ret = Symbol(Evolutionary.converged(opt_res))

    SciMLBase.build_solution(prob, opt, Evolutionary.minimizer(opt_res), Evolutionary.minimum(opt_res); original=opt_res, retcode=opt_ret)
end


text:
using GLPlot
GLPlot.init()


text:
"""
`module  JAC.PhotoIonizationFluores`
    ... a submodel of JAC that contains all methods for computing photo-excitation-autoionization cross
        sections and rates.
"""
module PhotoIonizationFluores

    using ..AutoIonization, ..Basics, ..ManyElectron, ..Radial, ..Nuclear, ..PhotoEmission

    """
    `struct  PhotoIonizationFluores.Settings`
        ... defines a type for the details and parameters of computing photon-impact excitation-autoionization
            pathways |i(N)>  --> |m(N)>  --> |f(N-1)>.

        + multipoles              ::Array{EmMultipole,1}               ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge,1}                  ... Specifies the gauges to be included into the computations.
        + printBefore  ::Bool                               ... True, if all energies and lines are printed before their evaluation.
        + selectPathways          ::Bool                               ... True if particular pathways are selected for the computations.
        + selectedPathways        ::Array{Tuple{Int64,Int64,Int64},1}  ... List of list of pathways, given by tupels (inital, inmediate, final).
        + maxKappa                ::Int64                              ... Maximum kappa value of partial waves to be included.
    """
    struct Settings
        multipoles                ::Array{EmMultipole,1}
        gauges                    ::Array{UseGauge,1}
        printBefore    ::Bool
        selectPathways            ::Bool
        selectedPathways          ::Array{Tuple{Int64,Int64,Int64},1}
        maxKappa                  ::Int64
    end

    """
    `PhotoIonizationFluores.Settings()`  ... constructor for the default values of photon-impact excitation-autoionizaton settings.
    """
    function Settings()
        Settings( Basics.EmMultipole[], UseGauge[], false,  false, Tuple{Int64,Int64,Int64}[], 0)
    end

    # `Base.show(io::IO, settings::PhotoIonizationFluores.Settings)`
    #		... prepares a proper printout of the variable settings::PhotoIonizationFluores.Settings.
    function Base.show(io::IO, settings::PhotoIonizationFluores.Settings)
        println(io, "multipoles:              $(settings.multipoles)  ")
        println(io, "gauges:                  $(settings.gauges)  ")
        println(io, "printBefore:  $(settings.printBefore)  ")
        println(io, "selectPathways:          $(settings.selectPathways)  ")
        println(io, "selectedPathways:        $(settings.selectedPathways)  ")
        println(io, "maxKappa:                $(settings.maxKappa)  ")
    end

    """
    `struct  PhotoIonizationFluores.Channel`
        ... defines a type for a photon-impact excitaton & autoionization channel that specifies all quantum numbers, phases
            and amplitudes.

        + excitationChannel  ::PhotoEmission.Channel       ... Channel that describes the photon-impact excitation process.
        + augerChannel       ::AutoIonization.Channel      ... Channel that describes the subsequent Auger/autoionization process.
    """
    struct  Channel
        excitationChannel    ::PhotoEmission.Channel
        augerChannel         ::AutoIonization.Channel
    end

    """
    `struct  PhotoIonizationFluores.Pathway`
        ... defines a type for a photon-impact excitation pathway that may include the definition of different excitation
            and autoionization channels and their corresponding amplitudes.

        + initialLevel        ::Level                  ... initial-(state) level
        + intermediateLevel   ::Level                  ... intermediate-(state) level
        + finalLevel          ::Level                  ... final-(state) level
        + photonEnergy        ::Float64                 ... energy of the (incoming) electron
        + electronEnergy      ::Float64                 ... energy of the (finally outgoing, scattered) electron
        + crossSection        ::EmProperty              ... total cross section of this pathway
        + hasChannels         ::Bool                    ... Determines whether the individual excitation and autoionization
                                                            channels are defined in terms of their multipole, gauge,
                                                            free-electron kappa, phases and the total angular momentum/parity
                                                            as well as the amplitude, or not.
        + channels            ::Array{PhotoIonizationFluores.Channel,1}  ... List of channels of this pathway.
    """
    struct  Pathway
        initialLevel          ::Level
        intermediateLevel     ::Level
        finalLevel            ::Level
        photonEnergy          ::Float64
        electronEnergy        ::Float64
        crossSection          ::EmProperty
        hasChannels           ::Bool
        channels              ::Array{PhotoIonizationFluores.Channel,1}
    end

    """
    `PhotoIonizationFluores.Pathway()`
        ... 'empty' constructor for an photon-impact excitation-autoionization pathway between a specified
            initial, intermediate and final level.
    """
    function Pathway()
        Pathway(Level(), Level(), Level(), 0., 0., 0., false, PhotoIonizationFluores.Channel[] )
    end

    # `Base.show(io::IO, pathway::PhotoIonizationFluores.Pathway)`
    #		... prepares a proper printout of the variable pathway::PhotoIonizationFluores.Pathway.
    function Base.show(io::IO, pathway::PhotoIonizationFluores.Pathway)
        println(io, "initialLevel:               $(pathway.initialLevel)  ")
        println(io, "intermediateLevel:          $(pathway.intermediateLevel)  ")
        println(io, "finalLevel:                 $(pathway.finalLevel)  ")
        println(io, "photonEnergy                $(pathway.photonEnergy)  ")
        println(io, "electronEnergy              $(pathway.electronEnergy)  ")
        println(io, "crossSection:               $(pathway.crossSection)  ")
        println(io, "hasChannels:                $(pathway.hasChannels)  ")
        println(io, "channels:                   $(pathway.channels)  ")
    end

    """
    `PhotoIonizationFluores.computePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                            grid::Radial.Grid, settings::PhotoIonizationFluores.Settings; output=true)`
        ... to compute the photo-excitation-fluorescence amplitudes and all properties as requested by the given settings.
            A list of lines::Array{PhotoIonizationFluores.Lines} is returned.
    """
    function  computePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                              settings::PhotoIonizationFluores.Settings; output=true)
        println("")
        printstyled("PhotoIonizationFluores.computePathways(): The computation of photo-excitation-fluorescence amplitudes starts now ... \n", color=:light_green)
        printstyled("-------------------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        #
        println("Not yet implemented: Data structures and properties still need to be worked out.")
        #
        pathways = "Not yet implemented !"
        #
        if    output    return( pathways )
        else            return( nothing )
        end
    end

end # module


text:
# Part 1: Read image and reshape to 25 x 6 x N layers then report required stat
image = reshape(open(f -> parse.(Int8, collect(readline(f))), "day08_input.txt"), 25, 6, :)
layer = image[:, :, argmin(mapslices(sum, map(iszero, image), dims = [1, 2])[:])]
println("count(1) * count(2) = $(sum(layer .== 1) * sum(layer .== 2))")

# Part 2: Collapse layers, 0 is black, 1 is white, 2 is transparent, first layer is top
rendered = mapslices(x -> first(filter(!isequal(2), x)), image, dims=[3])[:, :, 1]
println.(join.(eachrow(map(x -> x == 1 ? "█" : " ", rendered'))))

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
__precompile__()
#
#@everywhere module KOS
module KOS
#
using Distributions
#
export init, fn_run_wos, fn_wos_core_estimation!, fn_run_wos_many, fn_wos_core_estimation,update_rhs, setup_run_wos, Params, chunk!,vec_fn_wos_single_samp,fn_wos_single_samp,chunkchunk!,chunkchunkchunk!,chunkchunk2!,chunkchunkchunk2!

#
function __init__()
end
#
type Params
  fn_ext::Function  # exterior data
  fn_rhs::Function  # source (right-hand side) function
  d::Function       # signed distance to the boundary
  alpha::Float64    # stability index (Laplace exponent is alpha/2)
end

#####################################################################################
# Begin of WOSRNG functions
######################################################################################
######################################################################################
type WOSRNG
  beta_samples::Array{Float64}  # samples
  angle_samples::Array{Float64,2}  #
  radius_samples::Array{Float64}  #
  cradius_samples::Array{Float64}  #
  # iterator (second index)
  it::Int64  #
  #
  rv_length::Int64 # number of samples in first index
  beta_dist # beta distribution
  alpha # alpha
end
#
function wosrng_reset!(w::WOSRNG)
  # set dimension and reset iterator
  w.it=1
end
#
function get_beta!(w::WOSRNG)
  return w.beta_samples[w.it]
end
#
function get_angle!(w::WOSRNG)
  # get uniform sample
  return w.angle_samples[w.it,1:2]
end
#
function get_angle2!(w::WOSRNG)
  # get uniform sample
  return w.angle_samples[w.it,3:4]
end
#
function get_radius!(w::WOSRNG) # first
  c     =w.cradius_samples[w.it]
  sample=w.radius_samples[w.it]
  return sample,c
end
#
#
function wosrng_load(beta_dist,alp)
  # construct WOSRNG object
  #print("wosrng_load with rv_length=",rv_length,"  ")
  rv_length=10
  th=2*pi*rand(rv_length,2);
  u1=[cos.(th[:,1]) sin.(th[:,1]) cos.(th[:,2]) sin.(th[:,2])]
  u2=rand(rv_length).^(1/alp)
  u3=cdf(beta_dist,1-u2.^2)
  b=sqrt.(rand(beta_dist,rv_length));
  #
  return WOSRNG(b,u1,u2,u3,0,
                rv_length,beta_dist,alp)
end
#
function wosrng_extend(w::WOSRNG, grow)
  w.it+=1
  if (w.rv_length<grow+w.it)
    w.rv_length+grow
    #
    th=2*pi*rand(2*grow,2)
    u1=[cos.(th[:,1]) sin.(th[:,1]) cos.(th[:,2]) sin.(th[:,2])]
    w.angle_samples =vcat(w.angle_samples, u1)
    #
    w.radius_samples=vcat(w.radius_samples, rand(grow).^(1/w.alpha))
    w.cradius_samples=cdf(w.beta_dist, 1-w.radius_samples.^2)
    #
    w.beta_samples  =vcat(w.beta_samples,sqrt.(rand(w.beta_dist,grow)) )
  end
end
#
function wosrng_reload!(w::WOSRNG)
  # New random random variables
  th=2*pi*rand(w.rv_length,2);
  u1=[cos.(th[:,1]) sin.(th[:,1]) cos.(th[:,2]) sin.(th[:,2])]
  u2=rand(w.rv_length).^(1/w.alpha)
  u3=cdf(w.beta_dist,1-u2.^2)
  b=sqrt.(rand(w.beta_dist,w.rv_length));
  w.angle_samples =u1
  w.radius_samples=u2
  w.cradius_samples=u3
  w.beta_samples  =b
  #
  w.it   =0
end
######################################
#####################################################################################
# End of WOSRNG functions
######################################################################################
######################################################################################
#  Begin distance/example PDE function
#
function d_ball_0_1( x::Array{Float64,1} )
  # out=d_ball_0_1(x)
  #
  # Distance of points to the domain boundary
  # Domain is a ball of radius 1 centred at (0,0)
  # x: point
  R=1.
  return R-sqrt(sum(x.^2))
end
#
function vec_d_ball_0_1( xs )
  # out=d_ball_0_1(x)
  # same as above
  # points are rows of xs
  R=1.
  return (R-sqrt.(sum(xs.^2,length(size(xs)))))[:]
end
#

function d_box( x )
  # out=d_box(x)
  #
  # Distance of points to the domain boundary
  # Domain is a box [a,b]x[a,b]
  # x
  a=-1; b=1.
  return minimum([x[1]-a,b-x[1],x[2]-a,b-x[2]])
end
#
function vec_fn_zero(x::Array{Float64,1})
  # Zero exterior function
  # x: point
  return 0.
end
#
function vec_fn_zero(x::Array{Float64,2})
  # Zero exterior function
  # x: point(s) as rows
  return zeros(size(x)[1])
end
#
function  vec_fn_ext1(x::Array{Float64,1}, alpha)
  # Exterior function/model solution: Free-space
  # Green's function centred at P
  # From D1.7 of Bucur
  # See Section 7, Example 1 of paper
  #
  # x: point(s) as columns for evaluation
  #
  P=[2. 0.] # any point outside ball
  tmp= sum((broadcast(-,x,P).^2)).^(-1+(alpha/2))
  return tmp
end
#
function  vec_fn_ext1(x::Array{Float64,2}, alpha)
  # Exterior function/model solution: Free-space
  # Green's function centred at P
  # From D1.7 of Bucur
  # See Section 7, Example 1 of paper
  #
  # x: point(s) as columns for evaluation
  #
  P=[2. 0.0] # any point outside ball
  tmp=sum((broadcast(-,x,P).^2),2).^(-1+(alpha/2))
  return tmp[:]
end
#
function vec_fn_ext2(x::Array{Float64,2})  #
  # Exterior function: Gaussian with P (vectorised)
  # See Section 7, Example 2 of paper
  # Solution (numerical) for this problem with g=0 is u_hom_check
  #
  # x: point(s) as columns for evaluation
  #
  P=[2 0]; # any point outside ball
  tmp=exp(-sum(broadcast(-,x,P).^2,2))
  return tmp[:]
end
#
function vec_fn_ext2(x::Array{Float64,1})
  P=[2 0]; # any point outside ball
  tmp=exp(-sum(broadcast(-,x,P).^2))
  return tmp
end
#
function fn_rhs1(x::Array{Float64,2}, alpha )
  # Source-term function: Constant example
  # From Dyda, page 549, Table 3, Example 2 (rescaled)
  # Exact solution for this problem with f=0 is fn_exact1
  # x: point(s) as columns for evaluation
  return ones(size(x)[1])::Array{Float64,1}
end
#
function fn_rhs1(x::Array{Float64,1}, alpha )
  # Source-term function: Constant example
  # From Dyda, page 549, Table 3, Example 2 (rescaled)
  # Exact solution for this problem with f=0 is fn_exact1
  # x: point(s) as columns for evaluation
  return 1::Float64
end
#
function fn_rhs2( x::Array{Float64,1}, alpha )
  # Source-term function: Non-constant example
  # From Dyda, page 549, Table 3, Example 2
  # Exact solution for this problem with f=0 is u4
  #
  # x: point(s) as columns for evaluation
  w=alpha/2
  c=(2^alpha)*gamma(2+w)*gamma(1+w)
  tmp=c*(1-(1+w)*sum(x.^2))
  return tmp::Float64
end
#
function fn_rhs2( x::Array{Float64,2}, alpha )
  # Source-term function: Non-constant example
  # From Dyda, page 549, Table 3, Example 2
  # Exact solution for this problem with f=0 is u4
  #
  # x: point(s) as columns for evaluation
  w=alpha/2
  c=(2^alpha)*gamma(2+w)*gamma(1+w)
  tmp=c*(1-(1+w)*sum(x.^2, 2))
  return tmp[:]::Array{Float64,1}
end
#
#
function fn_c_rhs( x::Array{Float64,1} )
  # Right-hand side designed to expose slow coupling rates.  #
  # x: point(s) as columns for evaluation
  tmp=2+sum(x.^2)
  return tmp::Float64
end
#
function fn_c_rhs( x::Array{Float64,2} )
  tmp=2+sum(x.^2, 2)
  return tmp[:]::Array{Float64,1}
end
#
#
#
function fn_c_ext( x::Array{Float64,1})
  # Exterior values
  # Designed to expose slow coupling rates
  tmp=sin((sum(x.^2)))
  return tmp::Float64
end
#
function fn_c_ext( x::Array{Float64,2})
  #
  tmp=sin.((sum(x.^2, 2)))
  return tmp[:]::Array{Float64,1}
end
#
function fn_exact1(x,alpha)
  # Model solution: Constant source example
  # From Dyda, page 549, Table 3, Example 1
  # x: point(s) as columns for evaluation
  w=alpha/2; z=1-alpha/2;
  c=((2^alpha)*gamma(1+w)*w*beta(w,z)/ gamma(z))
  if d_ball_0_1(x)>0
   return (1.-sum(x.^2)).^(alpha/2)/c
  else
   return 0.
  end
end
#
function fn_exact2(x::Array{Float64,1}, alpha )
  # Model solution: Non-constant source example
  # From Dyda, page 549, Table 3, Example 2
  # x: point(s) as columns for evaluation
  if d_ball_0_1(x)>0
    return (1-sum(x.^2)).^(1+(alpha/2))
  else
    return 0.
  end
end
#
function example0(alpha) # all zero
  p=KOS.Params(vec_fn_zero,vec_fn_zero,vec_d_ball_0_1,alpha)
  exact = x->0. # dummy
  return p, exact
end
#
function example2(alpha) # quadratic rhs
  println("Example 2: quadratic rhs")
  p=KOS.Params(vec_fn_zero,x->fn_rhs2(x,alpha),vec_d_ball_0_1,alpha)
  exact = x->fn_exact2(x,alpha) #yes
  return p, exact
end
#
function example1(alpha) # constant rhs
  println("Example 1: constant rhs")
  p=KOS.Params(vec_fn_zero,x->fn_rhs1(x,alpha), vec_d_ball_0_1, alpha)
  exact = x->fn_exact1(x,alpha) #yes
  return p, exact
end
#
function example3(alpha) # Green's fn
  println("Example 3: Green's fn exterior")
  p=KOS.Params(x->vec_fn_ext1(x,alpha),vec_fn_zero,vec_d_ball_0_1,alpha)
  return p, x->vec_fn_ext1(x,alpha)
end
#
function example4(alpha)
  # Green's fn
  println("Example 4: Green's fn exterior of box")
  p=KOS.Params(fn_ext2,fn_zero,d_box,alpha)
  return p, fn_ext2
end
#
function example5(alpha)
  # Difficult for coupling
  println("Example 5: difficult for coupling")
  p=KOS.Params(fn_c_ext,fn_c_rhs,vec_d_ball_0_1,alpha)
  return p, (x->0.0)
end
##########################################################################################
##########################################################################################
# end PDE functions
##########################################################################################
function setup_run_wos(p::KOS.Params,dict::Dict)
  alp=p.alpha; z=(2-alp)/2; w=alp/2
  # (in_c)= 2 pi* c(alp,2)*B(z,w) on page 1 of ideas.pdf
  in_c = alp^(-1) *2^(-alp+1) *gamma(w)^(-2) *beta(z,w)
  dict["c"]=in_c
  #
  beta_dist=Beta(w,z)
  dict["dist_beta"]=  beta_dist
  # in_val=E[(1-I(R^2,z,w))] for R=X^(1/alp) and X~u(0,1)
  # = int (1-I(x^(2/alp),z,w)) dx
  function integrand(X)
    return cdf(beta_dist,1-X.^(2/alp))
  end
  in_val, error=quadgk(integrand,0,1)
  dict["val"]=in_val
  function my_step(d::Array{Float64,1})
    return vec_wos_step(d,beta_dist)
  end
  if haskey(dict,"no_samp_MC")
    in_no       =dict["no_samp_MC"]
  else
    in_no=1; dict["no_samp_MC"]=1
  end
  my_rhs_ball =update_rhs(p,dict)
  dist_fn=y::Array{Float64,1}->p.d(y)[1]
  ext_fn=y::Array{Float64,1}->p.fn_ext(y)[1]
  #
  return my_step,my_rhs_ball,dist_fn,ext_fn
end
#
function update_rhs(p::Params,dict::Dict)
  in_val     = dict["val"]
  no_samp    = dict["no_samp_MC"]
  in_c       = dict["c"]
  beta_dist  = dict["dist_beta"]
  #
  my_rhs_ball=(rhon::Array{Float64,2},rn::Array{Float64,1}) ->vec_eval_MC0(rhon, rn, in_c, in_val, beta_dist,p.alpha,p.fn_rhs)::Array{Float64,1}
  #
   return my_rhs_ball
end

##########################################################################################
function eval_MC0(rhon::Array{Float64,1}, rn::Float64,
    in_c::Float64,in_value::Float64, beta_dist::Distributions.Beta{Float64}, alp::Float64, fn_rhs::Function)
# rhon=current position, rn=current radius,
# in_c, in_value useful constant
# beta_dist useful Beta distribution
# alp=alpha
# fn_rhs=function
  r   = rand()^(1/alp)
  tmp = polar_quad0(r*rn, rhon, fn_rhs) * cdf(beta_dist,1-r^2)
  return (in_c*(rn^alp)*(tmp+in_value*fn_rhs(rhon)))::Float64
end
##########################################################################################
function vec_eval_MC0(rhon::Array{Float64,2}, rn::Array{Float64,1},
    in_c::Float64,in_value::Float64, beta_dist::Distributions.Beta{Float64}, alp::Float64, fn_rhs::Function)
# rhon=current position, rn=current radius,
# in_c, in_value useful constant
# beta_dist useful Beta distribution
# alp=alpha
# fn_rhs=function
  r   = rand()^(1/alp)
  tmp = vec_polar_quad0(r*rn, rhon, fn_rhs) * cdf(beta_dist,1-r^2)
  return (in_c*(rn.^alp).*(tmp+in_value*fn_rhs(rhon)))::Array{Float64,1}
end
###########################################################################################
function eval_MC_coupled(rhon::Array{Float64,1}, rn::Float64,
    in_c::Float64,in_value::Float64, beta_dist::Distributions.Beta{Float64},
    alp::Float64, fn_rhs::Function, w::WOSRNG)
# rhon=current position, rn=current radius,
# in_c, in_value useful constant
# beta_dist useful Beta distribution
# alp=alpha
# fn_rhs=function
  r,c  = get_radius!(w)#^(1/alp)
  tmp= polar_quad_coupled(r*rn, rhon, fn_rhs, w) *c
  # cdf(beta_dist,1-r^2)
  return (in_c*(rn^alp)*(tmp+in_value*fn_rhs(rhon)))::Float64
end
#########################################################################################
function eval_MC_coupled_alt(rhon::Array{Float64,1}, rn::Float64,
    in_c::Float64,in_value::Float64, beta_dist::Distributions.Beta{Float64},
    alp::Float64, fn_rhs::Function, w::WOSRNG)
# rhon=current position, rn=current radius,
# in_c, in_value useful constant
# beta_dist useful Beta distribution
# alp=alpha
# fn_rhs=function
  r,c = get_radius!(w)#^(1/alp)
  tmp=polar_quad_coupled_alt(r*rn, rhon, fn_rhs, w) * c
  #cdf(beta_dist,1-r^2)
  return (in_c*(rn^alp)*(tmp+in_value*fn_rhs(rhon)))::Float64
end
##########################################################################################
function  polar_quad0(r::Float64, rhon::Array{Float64,1},
  fn_rhs::Function)
  # Function fn_rhs(rhon+r*x) in terms of polar coordinates of x
  theta = 2*pi*rand()
  X=r*[cos(theta),sin(theta)]
  return (fn_rhs(rhon+X)-fn_rhs(rhon))::Float64
end
##########################################################################################
function  vec_polar_quad0(r::Array{Float64,1}, rhon::Array{Float64,2},
  fn_rhs::Function)
  # Function fn_rhs(rhon+r*x) in terms of polar coordinates of x
  theta = 2*pi*rand()
  X=r*[cos(theta),sin(theta)]'
  return (fn_rhs(rhon+X)-fn_rhs(rhon))::Array{Float64,1}
end
########################################################################################
function  polar_quad_coupled_alt(r::Float64, rhon::Array{Float64,1},
  fn_rhs::Function, w::WOSRNG)
  #
  # Function fn_rhs(rhon+r*x) in terms of polar coordinates of x
  #
  X = get_angle!(w)
  X=r*X;  Y=r*[-X[2],X[1]]
  mysum=0.25*(fn_rhs(rhon+X)+fn_rhs(rhon-X)+
              fn_rhs(rhon+Y)+fn_rhs(rhon-Y))
  return (mysum-fn_rhs(rhon))::Float64
end
########################################################################################
function  polar_quad_coupled(r::Float64, rhon::Array{Float64,1},
  fn_rhs::Function, w::WOSRNG)
  #
  # Function fn_rhs(rhon+r*x) in terms of polar coordinates of x
  #
  X = r*get_angle!(w)
  return (fn_rhs(rhon+X)-fn_rhs(rhon))::Float64
end
#######################################################################################
function wos_step(d::Float64, dist::Distributions.Beta{Float64})
  r     = d/sqrt(rand(dist))
  theta = rand()*2*pi
  return r*[cos(theta); sin(theta)]::Array{Float64,1}
end
#######################################################################################
function vec_wos_step(d::Array{Float64,1}, dist::Distributions.Beta{Float64})
  r     = d/sqrt(rand(dist))
  theta = rand()*2*pi
  return (r.*[cos(theta), sin(theta)]')::Array{Float64,2}
end
########################################################################################
function wos_step_coupled(d::Float64, w::WOSRNG)
  r     = d/get_beta!(w)
  X = get_angle2!(w)
  return r*X::Array{Float64,1}
end
#########################################################################################
#########################################################################################
# end key WOS functions
#########################################################################################
function  fn_wos_core_estimation!(val::Array{Float64,2},
    x0::Array{Float64,1},my_step::Function,
    my_rhs_ball::Function, dist_fn::Function,
    ext_fn::Function,  no_its::Int64,  method::Bool)
  # Executes WoS random walks for initial data x0
  # Results returned in val (first col homogeneou part, second col inhom part)
  # Produce no_its number of setup_run_wos
  # method=TRUE for non-zero rhs (false otherwise)
  # must be "coupled"=false
  y=deepcopy(x0); d0=dist_fn(x0);
  #
  if (method) # Use when rhs non-zero
    for i=1:no_its
        d=d0;   y=x0  # initial position of WoS
        val[i,2]=0.   # inhom
        while (d>0.)  # walk
          # inhomog part for current sphere (y,d)
          val[i,2] +=my_rhs_ball(reshape(y,(1,2)),[d])[1]
          y+=my_step([d])[1] #  next sphere
          d=dist_fn(y)
        end # Exit domain
        # Calculate homogeneous part at the exit point
        val[i,1]=ext_fn(y)# homogenous
      end # End i loop
  else # Use when rhs is zero and method set to False
    for i=1:no_its
        d=d0; y=x0
        while (d>0.)
          y+=my_step(d)#  next sphere
          d=dist_fn(y)
        end # Exit domain
        # Calculate the homogeneous part at the exit point
        val[i,1]=ext_fn(y)
    end # End i loop
  end
end
########################################################################################################
function  fn_wos_core_estimation(x0::Array{Float64,1},
 my_step::Function, my_rhs_ball::Function,dist_fn::Function,ext_fn::Function,
  no_its::Int64,  method::Bool)
  # Executes WoS random walks for initial data x0
  # Results [mean,variance] returned as return
  # Produce no_its number of setup_run_wos
  # method=TRUE for non-zero rhs (false otherwise)
  # must be "coupled"=false
  y=deepcopy(x0); d0=dist_fn(x0); mom1=0.; mom2=0.
  #
  if (method) # Use when rhs non-zero
    for i=1:no_its
        tmp=0.; d=d0;   y=x0  # initial position of WoS
        while (d>0)
          tmp+=my_rhs_ball(y,d);
          y+=my_step(d) #  next sphere
          d=dist_fn(y)
        end # Exit domain
        # Calculate homogeneous part at the exit point
        tmp+=ext_fn(y)# homogenous
        mom1+=tmp; mom2+=tmp^2
      end # End i loop
  else # Use when rhs is zero and method set to False
    for i=1:no_its
        d=d0; y=x0
        while (d>0.)
          y+=my_step(d)#  next sphere
          d=dist_fn(y)
        end # Exit domain
        # Calculate the homogeneous part at the exit point
        tmp=ext_fn(y); mom1+=tmp; mom2+=tmp^2
    end # End i loop
  end
  mom1/=no_its; mom2/=no_its
  return mom1, (mom2-mom1^2)
end
##################################################################################
########################################################################################################
function  fn_wos_single_samp(x0s::Array{Float64,2},
 my_step::Function, my_rhs_ball::Function,dist_fn::Function,ext_fn::Function, index_fn::Function,
    method::Bool)
  # Executes WoS random walks for initial data x0
  # Results [mean,variance] returned as return
  # Produce no_its number of setup_run_wos
  # method=TRUE for non-zero rhs (false otherwise)
  no_pts=size(x0s)[1]; out=zeros(no_pts)
  for i=1:no_pts
    index_fn()#
    x0=x0s[i,:]; d0=dist_fn(x0)
    #
    if (method) # Use when rhs non-zero
      tmp=0.; d=d0;   y=x0  # initial position of WoS
      while (d>0)
        tmp+=my_rhs_ball(y,d);
        y+=my_step(d) #  next sphere
        d=dist_fn(y)
      end # Exit domain
      # Calculate homogeneous part at the exit point
      tmp+=ext_fn(y)# homogenous
    else # Use when rhs is zero and method set to False
      d=d0; y=x0
      while (d>0.)
        y+=my_step(d)#  next sphere
        d=dist_fn(y)
      end # Exit domain
      # Calculate the homogeneous part at the exit point
      tmp=ext_fn(y)
    end
    out[i]=tmp
  end
  return out
end
#
function  vec_fn_wos_single_samp(x0s::Array{Float64,2},
 my_step::Function, my_rhs_ball::Function, dist_fn::Function, ext_fn::Function)
  # Executes WoS random walks for vector initial data x0 (points in rows)
  # Gives a  single vector of WoS samples (for each IC)
  # method=TRUE for non-zero rhs (false otherwise)
  no_pts=size(x0s)[1]
  tmp=zeros(no_pts)
  init_dist=dist_fn(x0s)
  #
  alive_index=find(z->(z>0),init_dist)
  d_alive=init_dist[alive_index]
  y=deepcopy(x0s)  # initial positions of WoS
  while (length(alive_index)>0)
    tmp[alive_index]+=my_rhs_ball(y[alive_index,:],d_alive)
    y[alive_index,:]+=my_step(d_alive) #  next sphere
    d_alive    =dist_fn(y[alive_index,:])
    sub_alive  =find(z->(z>0),d_alive)
    alive_index=alive_index[sub_alive]
    d_alive    =d_alive[sub_alive]
  end # all points exit domain
  # Calculate homogeneous part at the exit point
  tmp+=ext_fn(y)# homogenous
  return tmp
end
##################################################################################
####################################
function chunkchunk!(delta_sum::SharedArray,
                    no_samps,my_step,
                    my_rhs_ball,dist_fn,ext_fn,
                    finest_level_int,correction_level,pa,pts)
  idx=max(1,myid()-1) # processor number
  tmp=zeros(correction_level[end])
  m=0
  while (m<no_samps )
    # run wos
    tmp[correction_level]=vec_fn_wos_single_samp(pts,
                   my_step, my_rhs_ball, dist_fn, ext_fn)
    m+=1
    # update as corrections
    for i=finest_level_int # update on finest interior
      delta  = tmp[i] - mean(tmp[pa[i,:]])
      delta_sum[i,idx]+= delta
    end
  end # loop over samples
end
####################################
function chunkchunk2!(delta_sum::SharedArray,
                    delta_sq_sum::SharedArray,
                    no_samps,my_step,
                    my_rhs_ball,dist_fn,ext_fn,
                    finest_level_int,correction_level,pa,pts)
  idx=max(1,myid()-1)
  tmp=zeros(correction_level[end])
  #
  m=0
  while (m<no_samps )
     # run wos
     tmp[correction_level]=vec_fn_wos_single_samp(pts,
                    my_step, my_rhs_ball, dist_fn, ext_fn)
     m+=1
     # update as corrections
     for i=finest_level_int # update on finest interior
        delta               = tmp[i] - mean(tmp[pa[i,:]])
        delta_sum[i,idx]   += delta
        delta_sq_sum[i,idx]+= delta^2
     end
   end # loop over samples
end
####################################
function chunkchunkchunk!(delta_sum::SharedArray, no_samps,
                          my_step, my_rhs_ball,dist_fn,ext_fn,
                          coarse_level, pts)
  #
  idx=max(1,myid()-1) # processor number
  m=0
  while (m<no_samps )
     # run wos
     delta_sum[coarse_level,idx]+=vec_fn_wos_single_samp(pts,
                   my_step, my_rhs_ball, dist_fn, ext_fn)
     m+=1
   end # loop over samples
end
####################################
function chunkchunkchunk2!(delta_sum::SharedArray, delta_sq_sum::SharedArray,
                    no_samps,
                    my_step, my_rhs_ball, dist_fn, ext_fn,
                    coarse_level, pts)
  idx=max(1,myid()-1) # processor number
  m=0
  while (m<no_samps )
     # run wos
     tmp=vec_fn_wos_single_samp(pts,
                   my_step, my_rhs_ball, dist_fn, ext_fn)
     delta_sum[coarse_level,idx]+=tmp
     delta_sq_sum[coarse_level,idx]+=tmp.^2
     m+=1
  end # loop over samples
end
######################
############################
function fn_run_wos_many(x0::Array{Float64,2},
  p::Params, dict::Dict, my_step::Function,
  my_rhs_ball::Function, dist_fn::Function,
  ext_fn::Function)
  #
  # Private driver routine for Walk-on-spheres algorithm
  #
  # x0: initial position
  # p: params for PDE
  # dict: dict of method parameters*
  # return mean_val, var_sample, err_est, no_samples
  #
  # initialise
  no_per_batch=dict["no_per_batch"]
  max_samples =dict["max_samples"]/dict["no_per_batch"]
  # two error estimates (for hom, inhom parts)
  err_est=dict["tol"]+ones(2); num=zeros(2); val=zeros(no_per_batch,2)
  # initialise
  no_particles=size(x0)[end]
  val=zeros(no_per_batch,2)
  mean_val=zeros(no_particles)
  #
  if dict["coupled"]==true
    wosrng_reload!(dict["WOSRNG"])
  end
  #
  assert(dict["coupled"]==false)
  index_fn=dict["index_fn"]
  #
  for i=1:(size(x0)[end])
    if dist_fn(x0[:,i])>0
      fn_wos_core_estimation!(val, x0[:,i], my_step, my_rhs_ball, dist_fn, ext_fn, index_fn, no_per_batch, true)
      mean_val[i]   +=2*mean(val)
    end
  end
  #
  return mean_val
end
##########################################################################################
#
function fn_run_wos(x0::Array{Float64,1}, p::Params, dict::Dict,
 my_step::Function, my_rhs_ball::Function, dist_fn::Function,
 ext_fn::Function)
  #
  # Private driver routine for Walk-on-spheres algorithm
  #
  # x0: initial position
  # p: params for PDE
  # dict: dict of method parameters*
  # return mean_val, var_sample, err_est, no_samples
  #
  # initialise
  no_per_batch=dict["no_per_batch"]
  max_samples =dict["max_samples"]/dict["no_per_batch"]
  # two error estimates (for hom, inhom parts)
  err_est=dict["tol"]+ones(2); num=zeros(2); val=zeros(no_per_batch,2)
  # initialise
  sum_val    =zeros(2); sum_sq_val=zeros(2);  mean_val=zeros(2)
  mean_sq_val=zeros(2); var_sample=zeros(2)
  #index_fn=dict["index_fn"]
  #
  while ( err_est[2]>(dict["tol"]/2) && num[1]<max_samples)
    # stopping criteria on inhomogeneous part
    num[1] += 1
    fn_wos_core_estimation!(val, x0, my_step, my_rhs_ball, dist_fn, ext_fn,  no_per_batch, true)
    sum_val   +=mean(val,1)[:]
    sum_sq_val+=mean(val.^2,1)[:]
    mean_val   =sum_val/num[1]
    mean_sq_val=sum_sq_val/num[1]
    var_sample =mean_sq_val-mean_val.^2
    err_est=sqrt.(var_sample/num[1])
  end
  num[2]=num[1]
  #
  while err_est[1]>(dict["tol"]/2) && num[1]<max_samples
    # homogeneous part only
    num[1]+=1
    fn_wos_core_estimation!(val,x0, my_step, my_rhs_ball,dist_fn,ext_fn,  no_per_batch,false)
    sum_val[1]   =sum_val[1]   +mean(val[:,1])
    sum_sq_val[1]=sum_sq_val[1]+mean(val[:,1].^2)
    mean_val[1]   =sum_val[1]/num[1];
    mean_sq_val[1]=sum_sq_val[1]/num[1];
    var_sample[1] =mean_sq_val[1]-mean_val[1]^2
    err_est[1]=sqrt(var_sample[1]/num[1]);
  end
  no_samples=no_per_batch * num
  #
  return mean_val, var_sample, err_est, no_samples
end
#
######################################################
end# end of module KOS


text:
function converter_cost(data)
    for(c,conv) in data["convdc_ne"]
        conv["cost"] = conv["Pacmax"]*0.083 *100+ 28
        #display(conv["cost"])
    end
end

function converter_parameters_rxb(data)
for (c,conv) in data["convdc_ne"]
    # display("converter")
    # xtf = conv["xtf"]
    # bf = conv["bf"]
    # xc = conv["xc"]
    # Pmax = conv["Pacmax"]
    # Qmax = conv["Qacmax"]
    # println("convter:$c", "xtf:$xtf", "bf:$bf", "xc:$xc", "P: $Pmax", "Q: $Qmax")

    bus = conv["busac_i"]
    #display(bus)
    base_kV = data["bus"]["$bus"]["base_kv"]
    base_S = sqrt((100*conv["Pacmax"])^2+(100*conv["Qacmax"])^2) #base MVA = 100
    base_Z = base_kV^2/base_S # L-L votlage/3 phase power
    base_Y= 1/base_Z
    #display("baseS:$base_S")
    conv["xtf"] = 0.10*100/base_S #new X =old X *(100MVA/old Sbase)
    # display(conv["xtf"])
    # display(base_Z)
    conv["rtf"] = conv["xtf"]/100
    conv["bf"] = 0.08*base_S/100
    conv["xc"] = 0.07*100/base_S #new X =old X *(100MVA/old Zbase)
    # display(conv["xc"])
    # display(base_Z)
    conv["rc"] = conv["xc"]/100 #new X =old X *(100MVA/old Zbase)
    rtf = conv["rtf"]
    xtf = conv["xtf"]
    bf = conv["bf"]
    Pmax = conv["Pacmax"]
    Pmin =  conv["Pacmin"]
    Qmax = conv["Qacmax"]
    Qmin =  conv["Qacmin"]

    conv["Imax"] = sqrt(Pmax^2+Qmax^2)
    xc = conv["xc"]
    rc = conv["rc"]
    Imax = conv["Imax"]

    # println("convter:$c", "xtf:$xtf", "bf:$bf", "xc:$xc","baseS:$base_S", "baseZ:$base_Z","Pmin: $Pmin", "Pmax: $Pmax", "Qmin: $Qmin", "Qmax: $Qmax")
    #println("rtf:$rtf","     ", "xtf:$xtf","     ", "bf:$bf", "     ","rc:$rc", "     ","xc:$xc", "     ","Imax:$Imax","     ","Pmin: $Pmin", "     ","Pmax: $Pmax" )
    if xtf > 0.1 || xc > 0.1
        # display("Casldfjew;afjz;lsdjf ;ljewbaf;ahjdfkj daslvncf;ahfdjasf newohf::$c")
    end
end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).
@recipe function f(
    parshapes::NamedTupleShape,
    samples::DensitySampleVector,
    parsel::NTuple{2,Symbol};
    intervals = standard_confidence_vals,
    colors = standard_colors,
    mean = false,
    std_dev = false,
    globalmode = false,
    localmode = true,
    diagonal = Dict(),
    upper = Dict(),
    right = Dict(),
    filter = true
)
    i = findfirst(x -> x == parsel[1], keys(parshapes))
    j = findfirst(x -> x == parsel[2], keys(parshapes))

    @series begin
        intervals --> intervals
        colors --> colors
        mean --> mean
        std_dev --> std_dev
        globalmode --> globalmode
        localmode --> localmode
        diagonal --> diagonal
        upper --> upper
        right --> right
        filter --> filter

        samples, (i,j)
    end
end

@recipe function f(
    posterior::PosteriorDensity,
    samples::DensitySampleVector,
    parsel::NTuple{2,Symbol};
    intervals = standard_confidence_vals,
    colors = standard_colors,
    mean = false,
    std_dev = false,
    globalmode = false,
    localmode = true,
    diagonal = Dict(),
    upper = Dict(),
    right = Dict(),
    filter = true
)

    i = findfirst(x -> x == parsel[1], keys(params_shape(posterior)))
    j = findfirst(x -> x == parsel[2], keys(params_shape(posterior)))

    @series begin
        intervals --> intervals
        colors --> colors
        mean --> mean
        std_dev --> std_dev
        globalmode --> globalmode
        localmode --> localmode
        diagonal --> diagonal
        upper --> upper
        right --> right
        filter --> filter

        samples, (i,j)
    end
end

@recipe function f(
    samples::DensitySampleVector,
    parsel::NTuple{2,Integer};
    intervals = standard_confidence_vals,
    colors = standard_colors,
    mean = false,
    std_dev = false,
    globalmode = false,
    localmode = true,
    diagonal = Dict(),
    upper = Dict(),
    right = Dict(),
    filter = true
)
    if filter
        samples = BAT.drop_low_weight_samples(samples)
    end

    pi_x, pi_y = parsel
    bins = get(plotattributes, :bins, 200)
    seriestype = get(plotattributes, :seriestype, :smallest_intervals)

    if seriestype == :histogram2d || seriestype == :histogram || seriestype == :hist
        colorbar = true
    else
        colorbar = false
    end

    xguide --> "\$\\theta_$(pi_x)\$"
    yguide --> "\$\\theta_$(pi_y)\$"

    h = fit(Histogram, (flatview(samples.params)[pi_x, :], flatview(samples.params)[pi_y, :]), FrequencyWeights(samples.weight), closed=:left, nbins=bins)

    if seriestype == :scatter
        base_markersize = get(plotattributes, :markersize, 1.5)

        acc = findall(x -> x > 0, samples.weight)
        rej = findall(x -> x <= 0, samples.weight)

        color = parse(RGBA{Float64}, get(plotattributes, :seriescolor, :green))
        label = get(plotattributes, :label, isempty(rej) ? "samples" : "accepted")

        @series begin
            seriestype := :scatter
            label := label
            markersize := [w < 1 ? base_markersize : base_markersize * sqrt(w) for w in samples.weight[acc]]
            markerstrokewidth := 0
            color := [w >= 1 ? color : RGBA(convert(RGB, color), color.alpha * w) for w in samples.weight[acc]]
            (flatview(samples.params)[pi_x, acc], flatview(samples.params)[pi_y, acc])
        end

        if !isempty(rej)
            @series begin
                seriestype := :scatter
                label := "rejected"
                markersize := base_markersize
                markerstrokewidth := 0
                color := :red
                (flatview(samples.params)[pi_x, rej], flatview(samples.params)[pi_y, rej])
            end
        end

    else

        @series begin

            seriestype --> seriestype
            intervals --> intervals
            colors --> colors
            diagonal --> diagonal
            upper --> upper
            right --> right

            h, (pi_x, pi_y)

        end

    end

    #------ stats -----------------------------
    stats = MCMCBasicStats(samples)

    mean_options = convert_to_options(mean)
    globalmode_options = convert_to_options(globalmode)
    localmode_options = convert_to_options(localmode)
    stddev_options = convert_to_options(std_dev)

    if mean_options != ()
        mx= stats.param_stats.mean[pi_x]
        my = stats.param_stats.mean[pi_y]

        Σ_all = stats.param_stats.cov
        Σx = Σ_all[pi_x, pi_x]
        Σy = Σ_all[pi_y, pi_y]

        @series begin
            seriestype := :scatter
            label := get(mean_options, "label", "mean") #: ($(@sprintf("%.2f", mx)), $(@sprintf("%.2f", my)))
            seriestype==:marginal ? subplot := 3 :
            markeralpha := get(mean_options, "markeralpha", 1)
            markercolor := get(mean_options, "markercolor", :black)
            markersize := get(mean_options, "markersize", 4)
            markershape := get(mean_options, "markershape", :circle)
            markerstrokealpha := get(mean_options, "markerstrokealpha", 1)
            markerstrokecolor := get(mean_options, "markerstrokecolor", :black)
            markerstrokestyle := get(mean_options, "markerstrokestyle", :solid)
            markerstrokewidth := get(mean_options, "markerstrokewidth", 1)
            colorbar := colorbar

            if(std_dev)
                xerror := sqrt(Σx)
                yerror := sqrt(Σy)
            end
           ([mx], [my])
        end
    end

   if globalmode_options != ()
        globalmode_x = stats.mode[pi_x]
        globalmode_y = stats.mode[pi_y]

        @series begin
            seriestype := :scatter
            label := get(globalmode_options, "label", "global mode") #: ($(@sprintf("%.2f", globalmode_x)), $(@sprintf("%.2f", globalmode_y)))
            seriestype==:marginal ? subplot := 3 :
            markeralpha := get(globalmode_options, "markeralpha", 1)
            markercolor := get(globalmode_options, "markercolor", :black)
            markersize := get(globalmode_options, "markersize", 4)
            markershape := get(globalmode_options, "markershape", :rect)
            markerstrokealpha := get(globalmode_options, "markerstrokealpha", 1)
            markerstrokecolor := get(globalmode_options, "markerstrokecolor", :black)
            markerstrokestyle := get(globalmode_options, "markerstrokestyle", :solid)
            markerstrokewidth := get(globalmode_options, "markerstrokewidth", 1)
            colorbar := colorbar

           ([globalmode_x], [globalmode_y])
        end
    end

    if localmode_options != ()
        localmode_values = calculate_localmode_2d(h)
        for (i, l) in enumerate(localmode_values)
         @series begin
            seriestype := :scatter
            if i==1 && length(localmode_values)==1
                label := get(localmode_options, "label", "local mode") #: ($(@sprintf("%.2f", l[1])), $(@sprintf("%.2f", l[2])))
            elseif i ==1
                label := get(localmode_options, "label", "local modes")
            else
                label :=""
            end

            seriestype == :marginal ? subplot := 3 :
            markeralpha := get(localmode_options, "markeralpha", 0)
            markercolor := get(localmode_options, "markercolor", :black)
            markersize := get(localmode_options, "markersize", 4)
            markershape := get(localmode_options, "markershape", :rect)
            markerstrokealpha := get(localmode_options, "markerstrokealpha", 1)
            markerstrokecolor := get(localmode_options, "markerstrokecolor", :black)
            markerstrokestyle := get(localmode_options, "markerstrokestyle", :solid)
            markerstrokewidth := get(localmode_options, "markerstrokewidth", 1)
            colorbar := colorbar

            ([localmode_values[i][1]], [localmode_values[i][2]])
            end
        end
    end
end