text:
export LM_HO_Param, LM_AHO_Param, AHO_Param, NonEquil_AHO_Param, Imag_Sig_Param, Imag_Sig_Lmd_Param,
       get_action, get_D_action, get_sde_funcs, getParamArray

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

####### MODEL ###########
abstract type AbstractModel end
abstract type AbstractOscillator <: AbstractModel end
abstract type Abstract_LM_Oscillator <: AbstractOscillator end

"""
    Imagenary sigma model parameters
"""
struct Imag_Sig_Param <: Abstract_LM_Oscillator
    σ::Complex
end

Base.string(p::Imag_Sig_Param) = string("σ: ", p.σ)

"""
    Imagenary sigma and lambda model parameters
"""
struct Imag_Sig_Lmd_Param <: Abstract_LM_Oscillator
    σ::Complex
    λ::Complex
end

Base.string(p::Imag_Sig_Lmd_Param) = string("σ: ", p.σ, "\n", "λ: ", p.λ)

"""
    Larg-Mass Harmonic oscillator parameters
"""
struct LM_HO_Param <: Abstract_LM_Oscillator
    m::Float64
end

Base.string(p::LM_HO_Param) = string("m: ", p.m)

"""
    Larg-Mass Anharmonic oscillator parameters
"""
struct LM_AHO_Param <: Abstract_LM_Oscillator
    m::Float64
    λ::Float64
end

Base.string(p::LM_AHO_Param) = string("m: ", p.m, "\n", "λ: ", p.λ)

"""
    Anharmonic oscillator parameters
"""
struct AHO_Param <: AbstractOscillator
    g::Float64
    m::Float64
    λ::Float64
end

Base.string(p::AHO_Param) = string("g: ", p.g, "\n", "m: ", p.m, "\n", "λ: ", p.λ)

"""
    Non-Equilibrium Anharmonic oscillator parameters
"""
struct NonEquil_AHO_Param <: AbstractOscillator
    g::Float64
    m::Float64
    λ::Float64

    # Source term parameters
    ϕ::Float64
    ξ::Float64
    σ::Float64
    ϕ_dot::Float64
    η::Float64
end

NonEquil_AHO_Param(m::Float64,λ::Float64,ϕ::Float64,ξ::Float64,σ::Float64,η::Float64) = NonEquil_AHO_Param(1,m,λ,ϕ,ξ,σ,0,η)
NonEquil_AHO_Param(m::Float64,λ::Float64,ϕ::Float64,ξ::Float64,σ::Float64) = NonEquil_AHO_Param(1,m,λ,ϕ,ξ,σ,0,0)
NonEquil_AHO_Param(m::Float64,λ::Float64) = NonEquil_AHO_Param(1,m,λ,1,1,1,0,0)

function Base.string(p::NonEquil_AHO_Param)
    string("g: ", p.g, "\n", "m: ", p.m, "\n", "λ: ", p.λ, "\n",
          "Source_term: Gaussian\n",
          "ϕ: ", p.ϕ, "\n", "ξ: ", p.ξ, "\n", "ϕ_dot: ", p.ϕ_dot, "\n",
          "η: ", p.η, "\n", "σ: ", p.σ)
end

###################################
# CONSTRUCTING THE PARAMETER ARRAY
###################################
"""
    Constructing the parameter array

    Arguments:
        p: Parameters of model
        a: Contour points spacing

    Return:
        LVector of arguments
"""
function getParamArray(p::T,a::AbstractArray) where T<:AbstractModel
    LVector(p=p,a=a,
            vals=fill(zeros(length(a)), 6),
            as = (circshift(a,1),circshift(a,-1),circshift(a,2))
    )
end
#########################
####### IMAG SIGMA ######
#########################
"""
(σ1+iσ2)(x1+ix2) = σ1*x1-σ2*x2 + i(σ1*x2+σ2*x2)
"""
function a_Imag_Sig_real(du,u,param,t)

    p = param.p

    du[1] = -(real(p.σ)*u[1] - imag(p.σ)*u[2])
    du[2] = -(real(p.σ)*u[2] + imag(p.σ)*u[1])
end

function b_Imag_Sig_real(du,u,param,t)
    du[1] = sqrt(2)
end

function get_sde_funcs(::Imag_Sig_Param)
    return a_Imag_Sig_real, b_Imag_Sig_real, nothing
end

################################
####### IMAG SIGMA LAMBDA ######
################################
"""
(σ1+iσ2)(x1+ix2) = σ1*x1-σ2*x2 + i(σ1*x2+σ2*x2)
(x1+ix2)(x1+ix2)(x1+ix2) = (x1*x1-x2*x2 + i(x1*x2 + x2*x2))(x1+ix2)
                         = x1^3 - 2*x1*x2^2 - x2^3 + i(x1^2 * x2 - x2^3 + x1^2 x2 + x1 x2^2)
"""
function a_Imag_Sig_Lmd_real(du,u,param,t)

    p = param.p

    du[1] = -(real(p.σ)*u[1] - imag(p.σ)*u[2]) - (1/6)*p.λ(u[1]*u[1]*u[1] - 2*u[1]*u[2]*u[2] - u[2]*u[2]*u[2])
    du[2] = -(real(p.σ)*u[2] + imag(p.σ)*u[1]) - (1/6)*p.λ(2*u[1]*u[1]*u[2] + u[1]*u[1]*u[2] - u[2]*u[2]*u[2])
end

function b_Imag_Sig_Lmd_real(du,u,param,t)
    du[1,1] = sqrt(2)
end

function get_sde_funcs(::Imag_Sig_Lmd_Param)
    return a_Imag_Sig_Lmd_real, b_Imag_Sig_Lmd_real, nothing
end

####################################
#######  AHO Non-Equilibrium  ######
####################################
"""
    Returns the a term in the SDE for the AHO model
"""
function a_AHO_real_nonEquil(du,u,param,t)

    p = param.p
    a = param.a

    a_R = real(a)
    a_I = imag(a)

    xR, xI, xR_diff_i, xR_diff_pi, xI_diff_i, xI_diff_pi = param.vals
    a_m1, a_p1, _ = param.as

    a_m1_R = real(a_m1)
    a_p1_I = imag(a_p1)

    xR =  @view u[1:div(end,2)]
    xI =  @view u[div(end,2)+1:end]

    xR_diff_i  = xR .- xR[vcat(end,1:end-1)]#circshift(xR,1)
    xR_diff_pi = xR[vcat(2:end,1)] .- xR
    xI_diff_i  = xI .- xI[vcat(end,1:end-1)]#circshift(xI,1)
    xI_diff_pi = xI[vcat(2:end,1)] .- xI

    ####
    # Non-Equilibrium factors
    ####
    #ϕ = 1
    #σ = 1
    #ξ = 1
    nonEqFac = ((p.σ^2 + 1)/ (8*p.ξ^2))
    nonEqFac2 = - ((p.σ^2 - 1)/ (4*p.ξ^2))
    nonEqFac3 = p.η/(2p.ξ)

    #########
    # i=0
    ########

    preFac_1 = 1 / ( abs(a[1]) )

    a_R_1 = a_R[1]
    a_I_1 = a_I[1]

    x_R_1 = xR[1]
    x_I_1 = xI[1]

    x_R_N = xR[end]
    x_I_N = xI[end]

    a2_1 = (a_R_1^2 + a_I_1^2)

    du[1] = preFac_1 * (
            (1/a2_1)*(a_R_1*xI_diff_pi[1]-a_I_1*xR_diff_pi[1])
          + (1/12)*(a_I_1*x_R_1*(6p.m+p.λ*(-3*x_I_1^2+x_R_1^2)) + a_R_1*x_I_1*(6p.m^2+p.λ*(-x_I_1^2+3*x_R_1^2)))
          + (2*nonEqFac * (p.ϕ-x_R_1) + nonEqFac2 * (p.ϕ - x_R_N) + 2nonEqFac3*x_I_1) # Gaussian distribution initial condition
    )

    du[div(end,2)+1] = preFac_1 * (
          - (1/a2_1)*(a_I_1*xI_diff_pi[1]+a_R_1*xR_diff_pi[1])
          - (1/12)*(a_R_1*x_R_1*(6p.m^2+p.λ*(-3*x_I_1^2+x_R_1^2)) - a_I_1*x_I_1*(6p.m^2 + p.λ*(-x_I_1^2+3*x_R_1^2)))
          - (2*nonEqFac * x_I_1 + nonEqFac2 * x_I_N - 2nonEqFac3*(p.ϕ - x_R_1) ) # Gaussian distribution initial condition
    )

    #########
    # i=N
    ########
    preFac_N = 1 / ( abs(a[end]) )

    a_R_Nm1 = a_R[end]
    a_I_Nm1 = a_I[end]

    a2_Nm1 = (a_R_Nm1^2 + a_I_Nm1^2)

    du[div(end,2)] = preFac_N * (
            (1/a2_Nm1)*(a_I_Nm1*xR_diff_i[end]-a_R_Nm1*xI_diff_i[end]) #* - a_R_Nm1*xI_diff_i[end] )
          + (1/12)*(a_I_Nm1*x_R_N*(6p.m + p.λ*(-3*x_I_N^2+x_R_N^2)) + a_R_Nm1*x_I_N*(6p.m^2 + p.λ*(-x_I_N^2+3*x_R_N^2)))
          + (2*nonEqFac * (p.ϕ-x_R_N) + nonEqFac2 * (p.ϕ - x_R_1) - 2nonEqFac3*x_I_N) # Gaussian distribution initial condition
    )

    du[end] = preFac_N * (
            (1/a2_Nm1)*(a_I_Nm1*xI_diff_i[end] + a_R_Nm1*xR_diff_i[end] )
          - (1/12)*(a_R_Nm1*x_R_N*(6p.m + p.λ*(-3*x_I_N^2+x_R_N^2)) - a_I_Nm1*x_I_N*(6p.m^2 + p.λ*(-x_I_N^2+3*x_R_N^2)))
          - (2* nonEqFac * x_I_N + nonEqFac2 * x_I_1 + 2nonEqFac3*(p.ϕ - x_R_N) )# Gaussian distribution initial condition
    )

    ###############
    # i=1 to i=N-1
    ###############
    preFac = 1 ./ (0.5 .* (abs.(a[2:end]) .+ abs.(a_m1[2:end])))

    a_R_j = a_R[2:end]
    a_I_j = a_I[2:end]
    a_R_jm1 = a_R[1:end-1]
    a_I_jm1 = a_I[1:end-1]

    x_R_j = xR[2:end-1]
    x_I_j = xI[2:end-1]
    xR_diff_j  = xR_diff_i[2:end-1]
    xR_diff_pj = xR_diff_pi[2:end-1]
    xI_diff_j  = xI_diff_i[2:end-1]
    xI_diff_pj = xI_diff_pi[2:end-1]

    du[2:div(end,2)-1] .= preFac .* (
        1 .* (a_I_jm1 .* xR_diff_j   .- a_R_jm1 .* xI_diff_j)   ./ (a_I_jm1.^2 .+ a_R_jm1.^2)
     .+ 1 .* (a_R_j   .* xI_diff_pj .- a_I_j   .* xR_diff_pj) ./ (a_I_j.^2   .+ a_R_j.^2  )
     .+ (1. / 12.) .* (
            (a_I_jm1 .+ a_I_j) .* x_R_j .* (6*p.m^2 .+ p.λ .* ( .- 3 .* x_I_j.^2 .+       x_R_j.^2 ))
         .+ (a_R_jm1 .+ a_R_j) .* x_I_j .* (6*p.m^2 .+ p.λ .* ( .-      x_I_j.^2  .+ 3 .* x_R_j.^2 ))
        )
    )

    du[div(end,2)+2:end-1] .=  preFac .* (
        1 .* (a_I_jm1 .* xI_diff_j   .+ a_R_jm1 .* xR_diff_j)   ./ (a_I_jm1.^2 .+ a_R_jm1.^2)
     .- 1 .* (a_I_j   .* xI_diff_pj .+ a_R_j   .* xR_diff_pj) ./ (a_I_j.^2   .+ a_R_j.^2  )
     .+ (1. / 12.) .* (
            (a_I_jm1 .+ a_I_j) .* x_I_j .* (6*p.m^2 .+ p.λ .*( .-      x_I_j.^2  .+ 3 .* x_R_j.^2))
         .- (a_R_jm1 .+ a_R_j) .* x_R_j .* (6*p.m^2 .+ p.λ .*( .- 3 .* x_I_j.^2 .+       x_R_j.^2))
        )
    )

end

function b_AHO_real_nonEquil(du,u,param,t)
    a = param.a

    du[1] = sqrt(2) / sqrt(abs(a[1]))
    du[2:div(end,2)-1] .= sqrt.(2 ./ abs.(a[2:end]))
    du[div(end,2)] = sqrt(2) / sqrt(abs(a[end]))

end

function get_sde_funcs(::NonEquil_AHO_Param)
    return a_AHO_real_nonEquil, b_AHO_real_nonEquil, nothing
end

##################
####### AHO ######
##################
"""
    Returns the a term in the SDE for the AHO model
"""
function a_AHO_real(du,u,param,t)

    p = param.p
    a = param.a

    xR, xI, xR_diff_i, xR_diff_pi, xI_diff_i, xI_diff_pi = param.vals
    a_m1, a_p1, _ = param.as

    xR =  @view u[1:div(end,2)]
    xI =  @view u[div(end,2)+1:end]

    xR_diff_i  = xR .- xR[vcat(end,1:end-1)]#circshift(xR,1)
    xR_diff_pi = xR[vcat(2:end,1)] .- xR
    xI_diff_i  = xI .- xI[vcat(end,1:end-1)]#circshift(xI,1)
    xI_diff_pi = xI[vcat(2:end,1)] .- xI

    preFac = (1.0 ./ (0.5.*(abs.(a) + abs.(a_m1))) ) .* 0.5

    du[1:div(end,2)] .= - preFac .* (
        2. .* p.g .*( real.(a_m1).*xI_diff_i  .- imag.(a_m1).*xR_diff_i ) ./ abs.(a_m1).^2
     .- 2. .* p.g .*( real.(a)  .* xI_diff_pi .- imag.(a)  .* xR_diff_pi )./ abs.(a).^2
     .- real.(a_m1 .+ a) .* (p.m.*xI .+ (1/6)*p.λ .* (-(xI.^3) .+ 3xI.*(xR.^2)))
     .- imag.(a_m1 .+ a) .* (p.m.*xR .+ (1/6)*p.λ .* (xR.^3 .- 3xR.*(xI.^2)))
     #.+ 1e-8
    )

    du[div(end,2)+1:end] .= preFac .* (
        2. .* p.g .* (real.(a_m1).*xR_diff_i .+ imag.(a_m1).*xI_diff_i) ./ abs.(a_m1).^2
     .- 2. .* p.g .* (real.(a).*xR_diff_pi   .+ imag.(a).*xI_diff_pi) ./ abs.(a).^2
     .- real.(a_m1 .+ a) .* (p.m.*xR .+ (1/6)*p.λ .* ((xR.^3) .- 3xR.*(xI.^2)))
     .+ imag.(a_m1 .+ a) .* (p.m.*xI .+ (1/6)*p.λ .* (-xI.^3 .+ 3xI.*(xR.^2)))
    )
end

"""
    Returns the jacobian of the a term in the SDE for the AHO model
"""
function jac_AHO_real(J,u,param,t)

    p = param.p
    a = param.a

    t_steps = div(length(u),2)

    xR =  u[1:t_steps]
    xI =  u[t_steps+1:end]
    a_m1 = circshift(a,1)
    a_p1 = circshift(a,-1)

    for i in 1:t_steps
    preFac = (1.0 / (0.5*(abs(a[i]) + abs(a_m1[i]))) ) * 0.5
    J[i,i] = - preFac * (
       - 2. * p.g * imag(a_m1[i])  / abs(a_m1[i])^2
       - 2. * p.g * imag(a[i])        / abs(a[i])^2
       - real(a_p1[i] + a[i]) * ((1/6)*p.λ * ( 6 * xI[i] * xR[i] ))
       - imag(a_p1[i] + a[i]) * (p.m + (1/6)*p.λ * (3*(xR[i])^2 - 3*(xI[i]).^2))
    )

    J[t_steps+i,t_steps+i] = preFac * (
        2. * p.g  * imag(a_m1[i]) / abs(a_m1[i])^2
      + 2. * p.g  * imag(a[i])    / abs(a[i])^2
      + real(a_p1[i] + a[i]) * (1/6)*p.λ * ( 6*xR[i]*xI[i] )
      + imag(a_p1[i] + a[i]) * (p.m + (1/6)*p.λ * (-3*xI[i]^2 + 3(xR[i]^2)))
    )

    J[i,t_steps+i] = - preFac * (
        2. * p.g * real(a_m1[i]) / abs(a_m1[i])^2
      + 2. * p.g * real(a[i])    / abs(a[i])^2
      - real(a_p1[i] + a[i]) * (p.m + (1/6)*p.λ * (-3(xI[i]^2) + 3(xR[i]^2)))
     .- imag(a_p1[i] + a[i]) * ( (1/6)*p.λ * (6*xR[i]*xI[i]))
    )

    J[t_steps+i,i] = preFac * (
        2. * p.g * real(a_m1[i]) / abs(a_m1[i])^2
      + 2. * p.g * real(a[i])   / abs(a[i])^2
      - real(a_p1[i] + a[i]) * (p.m + (1/6)*p.λ * 3*((xR[i]^2) - 3*(xI[i]^2)))
      + imag(a_p1[i] + a[i]) * (1/6)*p.λ * (6*xI[i]*xR[i])
    )

    J[i,mod1(i+1,t_steps)] = - preFac * (
      + 2. * p.g * imag(a[i]) / abs(a[i])^2
    )

    J[i,t_steps + mod1(i+1,t_steps)] = - preFac * (
      - 2. * p.g *  real(a[i]) / abs(a[i])^2
    )

    J[i,mod1(i-1,t_steps)] = - preFac * (
        2. * p.g * imag(a_m1[i])  / abs(a_m1[i])^2
    )

    J[i,t_steps + mod1(i-1,t_steps)] = - preFac * (
        -2. * p.g * real(a_m1[i])  / abs(a_m1[i])^2
    )

    J[t_steps+i,t_steps + mod1(i+1,t_steps)] = preFac * (
        - 2. * p.g * imag(a[i])/ abs(a[i])^2
    )

    J[t_steps+i, mod1(i+1,t_steps)] = preFac * (
      - 2. * p.g * real(a[i])   / abs(a[i])^2
    )

    J[t_steps+i,t_steps + mod1(i-1,t_steps)] = preFac * (
        -2. * p.g * imag(a_m1[i]) / abs(a_m1[i])^2
    )

    J[t_steps+i, mod1(i-1,t_steps)] = preFac * (
        -2. * p.g * real(a_m1[i]) / abs(a_m1[i])^2
    )

    end
end

"""
    Returns the b term in the SDE for the AHO model
"""
function b_AHO_real(du,u,param,t)
    a = param.a
    a_m1, a_p1, _ = param.as
    du[1:div(end,2)] .= sqrt.(2 ./ (0.5.*(abs.(a) + abs.(a_m1))))
end

function get_sde_funcs(::AHO_Param)
    return a_AHO_real, b_AHO_real, jac_AHO_real
end


text:
using RedPitayaDAQServer
using Statistics
using DataFrames
using CSV

include("../../config.jl")

baseFrequency = 125000000

function configureCluster(rpc::RedPitayaCluster, dec, modulus=4800, periods_per_frame=10)
    samples_per_period = div(modulus, dec)
    decimation(rpc, dec)
    samplesPerPeriod(rpc, samples_per_period)
    periodsPerFrame(rpc, periods_per_frame)
    modeDAC(rpc, "STANDARD")
    triggerMode(rpc, "EXTERNAL")
    ramWriterMode(rpc, "TRIGGERED")

    for rp in rpc.rp
        frequencyDAC(rp, 1, 1, baseFrequency/modulus)
        amplitudeDACNext(rp, 1, 1, 0.5)
        phaseDAC(rp, 1, 1, 0.0 )
        signalTypeDAC(rp, 2 , "SINE")

        frequencyDAC(rp, 2, 1, baseFrequency/modulus)
        amplitudeDACNext(rp, 2, 1, 0.5)
        phaseDAC(rp, 2, 1, 0.0 )
        signalTypeDAC(rp, 2 , "SINE")
    end
end

function readXFullBuffers(rpc::RedPitayaCluster, rpInfo::RPInfo, x)
    adcBufferSize = (1 << 25)
    totalSamples = adcBufferSize * x
    println("ADC Start")
    masterTrigger(rpc, false)
    startADC(rpc)
    masterTrigger(rpc, true)
    readPipelinedSamples(rpc, 0, totalSamples, rpInfo=rpInfo, chunkSize = Int64(ceil(totalSamples/100)))
    stopADC(rpc)
    masterTrigger(rpc, false)
    println("ADC End")
    return totalSamples
end

function managedToReadAllSamples(rpInfo::RPInfo)
    for perf in rpInfo.performances
        for chunkPerf in perf.data
            if chunkPerf.status.overwritten || chunkPerf.status.corrupted
                return false
            end
        end
    end
    return true
end

struct LongTermResults
    readInit::Float64
    readMin::Float64
    readMax::Float64
    readMean::Float64
    readMedian::Float64
    readSTD::Float64
end

function longTermTestChunks(rpc, testTimeInSeconds=60, chunkTimeInSeconds=30)

    secondsRead = 0
    dec = decimation(rpc)
    numberOfChunks = Int64(ceil(testTimeInSeconds / chunkTimeInSeconds))
    println("Continously reading for $testTimeInSeconds sec. in groups of $chunkTimeInSeconds seconds")
    rpInfo = RPInfo(rpc)
    println("ADC Start")
    masterTrigger(rpc, false)
    startADC(rpc)
    masterTrigger(rpc, true)
    # frameStart = enableSlowDAC(rp, true, test_frames, 0.5, 1.0)
    println("Start reading samples")
    samplesRead = 0
    for i = 1:numberOfChunks
        currentChunk = min(testTimeInSeconds - secondsRead, chunkTimeInSeconds)
        samplesInChunk = Int64(ceil(currentChunk * (baseFrequency/dec)))
        println("Chunk $i started")
        @time readPipelinedSamples(rpc, samplesRead,  samplesInChunk, rpInfo = rpInfo, chunkSize = 10066330)
        samplesRead += samplesInChunk
        secondsRead += currentChunk
        println("Chunk $i ended")
    end
    println("End reading samples")
    stopADC(rpc)
    println("ADC Stop")

    println("Computing Stats")
    results = LongTermResults[]
    adcBufferSize = (1 << 25)
    for (rp, perf) in enumerate(rpInfo.performances)
        deltaReads = [perf.data[i].adc.deltaRead / (adcBufferSize) for i = 1:length(perf.data)]
        init = deltaReads[1]
        readMean = mean(deltaReads)
        readMedian = median(deltaReads)
        readVariance = var(deltaReads)
        readStd = std(deltaReads)
        readMin = minimum(deltaReads)
        readMax = maximum(deltaReads)
        result = LongTermResults(init, readMin, readMax, readMean, readMedian, readStd)
        push!(results, result)
        println("DeltaRead Stats for RedPitaya $rp:")
        println("Init: $init")
        println("Min $readMin, Max $readMax")
        println("Mean $readMean, Median $readMedian")
        println("Standard Deviation $readStd")
    end
    return results
end

function persistRPInfo(fileName, rpInfo::RPInfo)
    df = DataFrame(rp = Int[], wpRead = Int[], deltaRead = Int[], deltaSend = Int[])
    for (i, perf) in enumerate(rpInfo.performances)
        # Somewhat inefficient
        for perfData in perf.data
            push!(df, (i, perfData.wpRead, perfData.adc.deltaRead, perfData.adc.deltaSend))
        end
    end
    CSV.write(string(fileName, ".csv"), df)
end

function persistLongTermResults(fileName, results::Array{LongTermResults,1})
    df = DataFrame(rp = Int[], readInit = Float64[], readMin = Float64[], readMax = Float64[], readMean = Float64[], readMedian = Float64[], readSTD = Float64[])
    for (i, result) in enumerate(results)
        push!(df, (i, result.readInit, result.readMin, result.readMax, result.readMean, result.readMedian, result.readSTD))
    end
    CSV.write(string(fileName, ".csv"), df)
end

function fullyTestCluster(connectionStrings, start=1, longTermTimeInSeconds= 60, testName="test"; availableBytes=2147483648)
    # Init Test, see if connections works, first compile of functions
    clusterTestValid = true
    try
        println("Testing if Cluster is valid")
        rpc = RedPitayaCluster(connectionStrings)
        configureCluster(rpc, 64)
        readXFullBuffers(rpc, RPInfo(rpc), 1)
        println()
    catch e
        println("Cluster Valid Test failed")
        println(e)
        clusterTestValid = false
    end

    if clusterTestValid
        # Detailed Test
        println("Detailed Tests")
        decToTest = [64, 32, 16, 8]
        workingCombinations = Pair{Int64, Int64}[]
        for dec in decToTest
            for i = start:length(connectionStrings)
                println("Detailed Test: Dec $dec and $i RedPitayas")
                rpc = RedPitayaCluster(collect(Iterators.take(connectionStrings, i)))
                configureCluster(rpc, dec)
                println("Configured RedPitayas")
                rpInfo = RPInfo(rpc)
                samplesRead = readXFullBuffers(rpc, rpInfo, 3)
                println("Finishing Test")
                persistRPInfo(joinpath(@__DIR__, string(testName, "_detailed_rpc_",i, "_dec_", dec)), rpInfo)
                if managedToReadAllSamples(rpInfo)
                    push!(workingCombinations, Pair(dec, i))
                else
                    println("Samples were lost")
                end
            end
        end

        println("")
        println("Long Term Tests")
        for p in workingCombinations
            dec = p.first
            i = p.second
            println("Long Term Test: Dec $dec and $i RedPitayas")
            rpc = RedPitayaCluster(collect(Iterators.take(connectionStrings, i)))
            configureCluster(rpc, dec)
            println("Configured RedPitayas")

            bytesPerRedPitaya = availableBytes / i
            samplesPerRedPitaya = bytesPerRedPitaya / 4
            timePerChunk =  Int64(ceil(samplesPerRedPitaya / (baseFrequency/dec)))
            timePerChunk = min(longTermTimeInSeconds, timePerChunk)

            results = longTermTestChunks(rpc, longTermTimeInSeconds, timePerChunk)

            println("Finishing Test")
            persistLongTermResults(joinpath(@__DIR__, string(testName, "_long_rpc_",i, "_dec_", dec)), results)
        end

    end
end

ext: jl

max_stars_repo_name: tknopp/RedPitayaDAQServer

max_issues_repo_name: tknopp/RedPitayaDAQServer

max_forks_repo_name: tknopp/RedPitayaDAQServer


text:
#JuliaRL/src/shared.jl

using ReinforcementLearning
using ReinforcementLearningEnvironments
using Flux
using Dates
using Suppressor

"""
    pyenv2env(pyenv::PyObject)

Returns a `ReinforcementLearningEnvironments.GymEnv` environment
from the `gym` environment `pyenv`.
Code is from `ReinforcementLearningEnvironments.GymEnv(::String)` function.
"""
function pyenv2env(pyenv::PyObject)
    obs_space = convert(AbstractSpace, pyenv.observation_space)
    act_space = convert(AbstractSpace, pyenv.action_space)
    obs_type = if obs_space isa Union{MultiContinuousSpace,MultiDiscreteSpace}
        PyArray
    elseif obs_space isa ContinuousSpace
        Float64
    elseif obs_space isa DiscreteSpace
        Int
    elseif obs_space isa VectSpace
        PyVector
    elseif obs_space isa DictSpace
        PyDict
    else
        error("don't know how to get the observation type from observation space of $obs_space")
    end
    env = GymEnv{obs_type,typeof(act_space),typeof(obs_space),typeof(pyenv)}(
        pyenv,
        obs_space,
        act_space,
        PyNULL(),
    )
    RLBase.reset!(env) # reset immediately to init env.state
    env
end

"""
Returns a `LunarLander` environment from the `CustomGym` python module
as a `ReinforcementLearningEnvironments.GymEnv`
"""
function LunarLander(landingFactor=1, legFirstBonus=0)
    gym = pyimport("CustomGym")
    # reload the python module, otherwise changes to the code
    # won't be effective until you restart julia
    gym = pyimport("importlib").reload(gym.lunar_lander)
    gym.LunarLander(landing_factor=landingFactor, leg_first_bonus=legFirstBonus) |>
    pyenv2env |>
    ActionTransformedEnv(a -> a-1;mapping= a -> a+1)

end

function net_model(ns::Int, na::Int, rng)
    Chain(
        Dense(ns, 64, leakyrelu; initW = glorot_uniform(rng)),
        Dense(64, 64, leakyrelu; initW = glorot_uniform(rng)),
        Dense(64, 32, leakyrelu; initW = glorot_uniform(rng)),
        Dense(32, na; initW = glorot_uniform(rng)),
    ) |> cpu
end

function shallow_net_model(ns::Int, na::Int, rng)
    Chain(
        Dense(ns, 128, leakyrelu; initW = glorot_uniform(rng)),
        Dense(128, 64, leakyrelu; initW = glorot_uniform(rng)),
        Dense(64, na; initW = glorot_uniform(rng)),
    ) |> cpu
end

function deep_net_model(ns::Int, na::Int, rng)
    Chain(
        Dense(ns, 64, leakyrelu; initW = glorot_uniform(rng)),
        Dense(64, 64, leakyrelu; initW = glorot_uniform(rng)),
        Dense(64, 64, leakyrelu; initW = glorot_uniform(rng)),
        Dense(64, 32, leakyrelu; initW = glorot_uniform(rng)),
        Dense(32, na; initW = glorot_uniform(rng)),
    ) |> cpu
end

"""
@timeRet expr

Works like the macro `@time` from `Base`,
but returns the printed statistics as a string.
"""
macro timeRet(ex)
    quote
        while false; end # compiler heuristic: compile this block (alter this if the heuristic changes)
        local stats = Base.gc_num()
        local elapsedtime = time_ns()
        local val = $(esc(ex))
        elapsedtime = Base.time_ns() - elapsedtime
        local diff = Base.GC_Diff(Base.gc_num(), stats)
        local output = @capture_out Base.time_print(
            elapsedtime, diff.allocd,
            diff.total_time,
            Base.gc_alloc_count(diff)
        )
        print(output)
        println()
        output
    end
end

"""
    log_training_info(infos, agent, save_dir)

Creates file `training_infos.txt` inside of the `save_dir` directory,
containing training infos (agent structure and elapsed time)
"""
function log_training_info(infos, agent, save_dir)
    file_path = joinpath(save_dir, "training_infos.txt")
    open(f->write(f, infos*"\n", string(agent)), file_path, "w")
end

max_stars_repo_name: FeLusiani/JuliaRL

max_issues_repo_name: FeLusiani/JuliaRL

max_forks_repo_name: FeLusiani/JuliaRL


text:
using Makie, Distributions, Random
import GeometryBasics

mutable struct Bags
    sides
    color
    linewidth
    function Bags(; width, hight, color::Symbol = :green, linewidth = 10)
        sides = [
            Point2f0(0,hight) => Point2f0(0,0);
            Point2f0(0,0) => Point2f0(width,0);
            Point2f0(width,0) => Point2f0(width,hight);
        ]
        new(sides, color, linewidth)
    end
end

mutable struct Hives
    location::Point2f0
    Hives(x,y) = new(Point2f0(x,y))
end

abstract type AbsrtactBee end
struct Workers <: AbsrtactBee end
struct Scouts <: AbsrtactBee end
struct Inactives <: AbsrtactBee end

const worker = Workers()
const scout = Scouts()
const inactive = Inactives()

mutable struct Bees{P}
    node_x::Node
    node_y::Node
    marker::Symbol
    color::Symbol
    Bees(x::T, y::G, marker::Symbol, color::Symbol, type::Type) where {T <: Real, G <: Real} =
        new{type}(Node([x]), Node([y]), marker, color)
end

Bees(data::T) where {T <: Tuple} = Bees(data[1], data[2], data[3])
Bees(x, y, type::T) where {T <: AbsrtactBee} = Bees(T, x, y)
Bees(type::Type{Workers}, x, y) = Bees(x, y, :diamond, :red, type)
Bees(type::Type{Scouts}, x, y) = Bees(x, y, :star4, :purple, type)
Bees(type::Type{Inactives}, x, y) = Bees(x, y, :star4, :black, type)

mutable struct Foods{T <: GeometryBasics.AbstractPoint}
    location::Vector{T}
    # amount
end

mutable struct SimScene
    bag::Bags
    hive::Hives
    bees::Vector{Bees}
    food::Foods
    scene::Scene
    function SimScene(bag, hive, nbees::Tuple, food)
        bees = []
        bee_type = (worker, scout, inactive)
        for ind in 1:length(nbees)
            bees_x_location = fill(Float64(hive.location[1]), nbees[ind])
            bees_y_location = fill(Float64(hive.location[2]), nbees[ind])
            bees_type = fill(bee_type[ind] ,nbees[ind])
            append!(bees, Bees.(zip(bees_x_location, bees_y_location, bees_type)))
        end
        new(bag, hive, bees, food)
    end
end # mutable struct SimScene

"""
    plot_scene(sim_scene::SimScene)

this function plot the simulation scene

inputs:
    sim_scene - the main sim obj

output::
    nothing
"""
function plot_scene!(sim_scene::SimScene)
    bag = sim_scene.bag
    hive = sim_scene.hive
    bees = sim_scene.bees
    food = sim_scene.food
    sim_scene.scene = linesegments(bag.sides , color = bag.color, linewidth = bag.linewidth)
    scatter!(sim_scene.scene, hive.location; color = :orange, markersize = 8 ,marker = :pentagon)
    scatter!(sim_scene.scene, food.location; color = :blue, markersize = 8, marker = :star5)
    for bee in sim_scene.bees
        scatter!(sim_scene.scene, bee.node_x, bee.node_y ; marker = bee.marker, color = bee.color, markersize = 4)
    end
    sim_scene.scene |> display
    return nothing
end # function plot_scene

import AbstractPlotting.linesegments
linesegments(bag::Bags) = linesegments(bag.sides, color = bag.color, linewidth = bag.linewidth)

import Base: convert
convert(t::T, d::Bees{G}) where {T <: AbsrtactBee,G} = Bees(d.node_x.val[1], d.node_y.val[1], t)

"""
    update_plot_bees_location(newsim_scene::SimScene, newx::Vector, new_y::Vect
this function update the bees location in the main scene

inputs:
    sim_scene - the main sim obj
    new_x - a vector that holds all the bees X location
    new_y - a vector that holds all the bees Y location

output:
    nothing
"""
function update_plot_bees_location!(sim_scene::SimScene, new_x::Vector, new_y::Vector)
    sim_scene.bees.node_x[] = new_x
    sim_scene.bees.node_y[] = new_y
    return nothing
end # function update_plot_bees_location

"""
    convert_bee_type!(type::T, sim_scene::SimScene, bee_num::Int) where {T <: AbsrtactBee,G}

this function convert the bee to one of the following {worker, scout, inactive}

inputs:
    bee - the Bees object to convert
    type - the type to convert to
    bee_num - the index in the 'sim_scene.bees' array

output:
    bee - the converted bee object
"""
convert_bee_type!(sim_scene::SimScene, type::T, bee_num::Int) where {T <: AbsrtactBee,G} =
    sim_scene.bees[bee_num] = Bees( sim_scene.bees[bee_num].node_x.val[1],
                                    sim_scene.bees[bee_num].node_y.val[1],
                                    type)

"""
    bees_main_sim(sim_scene::SimScene)

this is the main simulation function

inputs:
    sim_scene - the simulation main obj

output:
    nothing
"""
function bees_main_sim(sim_scene::SimScene)
    body
end # function

#-------------------------- Test Simulation -----------------------------------
width = 100
hight = 100
bag = Bags(width = width, hight = hight)
hive = Hives(width/2,10)
food = Foods(vec(Point2f0[(80,40)]))
sim_scene = SimScene(bag, hive, (3,3,3), food)

plot_scene!(sim_scene)


text:
module QuantumESPRESSOHelpers

include("Inputs/Inputs.jl")
include("Outputs/Outputs.jl")

end # module


text:
function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{L, D},
    ::Type{S},
) where {
    L <: PSY.ControllableLoad,
    D <: AbstractControllablePowerLoadFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(L, sys)

    add_variables!(container, ActivePowerVariable, devices, D())
    add_variables!(container, ReactivePowerVariable, devices, D())

    # Add Variables to expressions
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    add_expressions!(container, ProductionCostExpression, devices, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{L, D},
    ::Type{S},
) where {
    L <: PSY.ControllableLoad,
    D <: AbstractControllablePowerLoadFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(L, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)

    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{L, D},
    ::Type{S},
) where {
    L <: PSY.ControllableLoad,
    D <: AbstractControllablePowerLoadFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(L, sys)

    add_variables!(container, ActivePowerVariable, devices, D())

    # Add Variables to expressions
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    add_expressions!(container, ProductionCostExpression, devices, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{L, D},
    ::Type{S},
) where {
    L <: PSY.ControllableLoad,
    D <: AbstractControllablePowerLoadFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(L, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)

    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{L, InterruptiblePowerLoad},
    ::Type{S},
) where {L <: PSY.ControllableLoad, S <: PM.AbstractPowerModel}
    devices = get_available_components(L, sys)

    add_variables!(container, ActivePowerVariable, devices, InterruptiblePowerLoad())
    add_variables!(container, ReactivePowerVariable, devices, InterruptiblePowerLoad())
    add_variables!(container, OnVariable, devices, InterruptiblePowerLoad())

    # Add Variables to expressions
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    add_expressions!(container, ProductionCostExpression, devices, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{L, InterruptiblePowerLoad},
    ::Type{S},
) where {L <: PSY.ControllableLoad, S <: PM.AbstractPowerModel}
    devices = get_available_components(L, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_constraints!(
        container,
        ReactivePowerVariableLimitsConstraint,
        ReactivePowerVariable,
        devices,
        model,
        S,
    )
    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)

    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{L, InterruptiblePowerLoad},
    ::Type{S},
) where {L <: PSY.ControllableLoad, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(L, sys)

    add_variables!(container, ActivePowerVariable, devices, InterruptiblePowerLoad())
    add_variables!(container, OnVariable, devices, InterruptiblePowerLoad())

    # Add Variables to expressions
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerVariable,
        devices,
        model,
        S,
    )

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    add_expressions!(container, ProductionCostExpression, devices, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{L, InterruptiblePowerLoad},
    ::Type{S},
) where {L <: PSY.ControllableLoad, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(L, sys)

    add_constraints!(
        container,
        ActivePowerVariableLimitsConstraint,
        ActivePowerVariable,
        devices,
        model,
        S,
    )
    add_feedforward_constraints!(container, model, devices)

    objective_function!(container, devices, model, S)

    add_constraint_dual!(container, sys, model)
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{L, StaticPowerLoad},
    ::Type{S},
) where {L <: PSY.ElectricLoad, S <: PM.AbstractPowerModel}
    devices = get_available_components(L, sys)

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)
    add_parameters!(container, ReactivePowerTimeSeriesParameter, devices, model)

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{L, StaticPowerLoad},
    ::Type{S},
) where {L <: PSY.ElectricLoad, S <: PM.AbstractActivePowerModel}
    devices = get_available_components(L, sys)

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ModelConstructStage,
    model::DeviceModel{L, StaticPowerLoad},
    ::Type{S},
) where {L <: PSY.ElectricLoad, S <: PM.AbstractPowerModel}
    # Static PowerLoad doesn't add any constraints to the model. This function covers
    # AbstractPowerModel and AbtractActivePowerModel
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{L, D},
    ::Type{S},
) where {
    L <: PSY.StaticLoad,
    D <: AbstractControllablePowerLoadFormulation,
    S <: PM.AbstractPowerModel,
}
    devices = get_available_components(L, sys)

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)
    add_parameters!(container, ReactivePowerTimeSeriesParameter, devices, model)

    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    add_to_expression!(
        container,
        ReactivePowerBalance,
        ReactivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ::ArgumentConstructStage,
    model::DeviceModel{L, D},
    ::Type{S},
) where {
    L <: PSY.StaticLoad,
    D <: AbstractControllablePowerLoadFormulation,
    S <: PM.AbstractActivePowerModel,
}
    devices = get_available_components(L, sys)

    add_parameters!(container, ActivePowerTimeSeriesParameter, devices, model)
    add_to_expression!(
        container,
        ActivePowerBalance,
        ActivePowerTimeSeriesParameter,
        devices,
        model,
        S,
    )
    return
end

function construct_device!(
    container::OptimizationContainer,
    sys::PSY.System,
    ccs::ModelConstructStage,
    model::DeviceModel{L, D},
    ::Type{S},
) where {
    L <: PSY.StaticLoad,
    D <: AbstractControllablePowerLoadFormulation,
    S <: PM.AbstractPowerModel,
}
    if D != StaticPowerLoad
        @warn(
            "The Formulation $(D) only applies to FormulationControllable Loads, \n Consider Changing the Device Formulation to StaticPowerLoad"
        )
    end

    # Makes a new model with the correct formulation of the type. Needs to recover all the other fields
    # slacks, services and duals are not applicable to StaticPowerLoad so those are ignored
    new_model = DeviceModel(
        L,
        StaticPowerLoad,
        feedforwards=model.feedforwards,
        time_series_names=model.time_series_names,
        attributes=model.attributes,
    )
    construct_device!(container, sys, ccs, new_model, S)
    return
end


text:
@testset "scale" begin
    data = ones(3, 2, 2)
    scales = [1, 2, 3]
    sc_cube = scale(data, scales)
    @test size(sc_cube) == (3, 6, 6)
    @test all(sc_cube[1, 3:4, 3:4] .≈ 1)
    @test all(sc_cube[2, 2:5, 2:5] .≈ 1)
    @test all(sc_cube[3, :, :] .≈ 1)
    @test sum(sc_cube, dims=(2, 3)) ≈ [4, 16, 36]

    for i in axes(data, 1)
        @test scale(data[i, :, :], scales[i], (6, 6)) ≈ sc_cube[i, :, :]
    end
end

@testset "invscale" begin
    data = ones(3, 6, 6)
    scales = [1, 2, 3]
    sc_cube = invscale(data, scales)
    @test size(sc_cube) == (3, 2, 2)
    @test all(sc_cube .≈ 1)

    for i in axes(data, 1)
        @test invscale(data[i, :, :], scales[i], (2, 2)) ≈ sc_cube[i, :, :]
    end
end

@testset "stack and collapse" begin
    data = ones(3, 4, 2, 2)
    scales = [1, 2, 3]
    stack = scale_and_stack(data, scales)
    @test size(stack) == (12, 6, 6)
    @test all(stack[1:3:12, 3:4, 3:4] .≈ 1)
    @test all(stack[2:3:12, 2:5, 2:5] .≈ 1)
    @test all(stack[3:3:12, :, :] .≈ 1)
    @test sum(stack) ≈ 4 * sum((4, 16, 36))

    cube = invscale_and_collapse(stack, scales, (2, 2))
    @test size(cube) == (4, 2, 2)
    @test all(cube .≈ 1)
end

@testset "scale list" begin
    waves = [100, 200, 300]
    @test scale_list(waves) ≈ [3, 1.5, 1]
end


text:
"""
    Prplex{T <: Real} = Perplexification{T}

A perplex number is a perplexification of the real numbers.
"""
const Prplex{T <: Real} = Perplexification{T}

function Prplex(a::T, b::T) where T <: Real
    Prplex{T}(a, b)
end

function Prplex(a::T) where T <: Real
    Prplex{T}(a, zero(T))
end

function Prplex(a::Real, b::Real)
    Prplex(promote(a, b)...)
end

"""
    unreal(z::Prplex)

The unreal part of a perplex number is a 1-vector.
"""
function unreal(z::Prplex{T}) where T <: Real
    T[z.r]
end

function show(io::IO, z::Prplex)
    print(io, "[1: ")
    print(io, z.l)
    print(io, ", R: ")
    print(io, z.r)
    print(io, "]")
end

function random(::Type{Prplex{T}}) where T <: Real
    Prplex{T}(rand(T), rand(T))
end

function random(::Type{Prplex{BigInt}})
    Prplex{BigInt}(randomBigInt(), randomBigInt())
end

function random(::Type{Prplex{BigFloat}})
    Prplex{BigFloat}(randomBigFloat(), randomBigFloat())
end

function random(::Type{Prplex{Rational{BigInt}}})
    Prplex{Rational{BigInt}}(randomBigRational(), randomBigRational())
end

function abs2(z::Prplex)
    abs2(z.l) - abs2(z.r)
end

function barR(z::Prplex)
    Prplex(z.l, -(z.r))
end

function iszerodivisor(z::Prplex)
    z.l == z.r || z.l == -(z.r)
end

function inv(z::Prplex)
    if iszerodivisor(z)
        error(ZeroDivisorInverse)
    end

    barR(z) / abs2(z)
end


text:
"""
This file contains all the functions in ClimateIsca.jl that deal with handling of Isca GCM
parameters, such as:
    - Loading of parameter information
    - Checking if parameter is a surface variable
    - Saving parameter data into NetCDF files
    - Reading parameter data from NetCDF files

"""

function iscaparis(pinfo::Dict,inquiry::AbstractString)

    tf = sum(inquiry==split(pinfo["name"],"_"))
    if !(tf==0); return true; else; return false; end

end

function iscaparfolder(path::AbstractString,parname::AbstractString)
    pfol = joinpath(path,parname)
end

function iscaparfolder(path::AbstractString,parname::AbstractString,level::Integer)
    pfol = joinpath(path,parname,"$(parname)-sig$(sprintf1("%02d",level))")
end

function iscaparload(parameterID::Integer)

    @debug "$(Dates.now()) - Loading information on parameters used in ClimateIsca.jl ..."
    params = readdlm(joinpath(@__DIR__,"parameters.txt"),',',comments=true);
    return params[params[:,1].==parameterID,2:end];

end

function iscaparload(parameterID::Array)

    @debug "$(Dates.now()) - Loading information on parameters used in ClimateIsca.jl ..."
    params = readdlm(joinpath(@__DIR__,"parameters.txt"),',',comments=true);
    cart = findall(params[:,1].==parameterID'); ind = [];
    for ii = 1 : length(cart); append!(ind,cart[ii][1]); end
    return params[ind,2:end];

end

function iscaparameter(parameterID::Integer)

    pinfo  = iscaparload(parameterID);
    return Dict("name"=>pinfo[1],"isca"=>pinfo[2],"full"=>pinfo[3],"unit"=>pinfo[4])

end

function iscapre2lvl(pinfo::Dict,allpre::Array)

    if !iscaparis(pinfo,"sfc")
        pre = [250,300,500,700,850,900]; pind = argmin(abs.(allpre.-pre'),dims=1);
        nlvl = length(bmin); level = zeros(nlvl)
        for ii = 1 : nlvl; level[ii] = pind[ii][1]; end
    else
        @info "$(Dates.now()) - Parameter chosen is a surface variable.  Setting level to 0."
        level = 0;
    end

    return level

end

## Extraction of Parameter Data

function iscaparncread(fnc::AbstractString,pinfo::Dict;timeshape=false)

    data = ncread(fnc,pinfo["isca"]); ndim = numel(size(data));

    if timeshape
        if numel(size(data)) == 4
            @debug "$(Dates.now()) - Parameter chosen has pressure levels and therefore 4 dimensions.  In order to automatically reduce dimensions after time-averaging, we shall permute the array such that the time dimension is the last dimension.";
            data = permutedims(data,[1,2,4,3]);
        end
    end

    if any(iscaparis.(pinfo,["prcp","rcnd","rcnv"])); data = data *24*60*60 end

    return data

end

function iscaparncread(fnc::AbstractString,pinfo::Dict,ncattr::Dict;timeshape=false)

    data = ncread(fnc,pinfo["isca"]); ndim = numel(size(data));

    if !iscaparis(pinfo,"sfc")
        pre = ncattr["fullpre"]; levels = iscapre2lvl(pinfo,pre); data = data[:,:,levels,:];
    end

    if timeshape
        if numel(size(data)) == 4
            @debug "$(Dates.now()) - Parameter chosen has pressure levels and therefore 4 dimensions.  In order to automatically reduce dimensions after time-averaging, we shall permute the array such that the time dimension is the last dimension.";
            data = permutedims(data,[1,2,4,3]);
        end
    end

    if any(iscaparis.(pinfo,["prcp","rcnd","rcnv"])); data = data *24*60*60 end

    return data

end

## Saving of Parameter Data

function iscaparncsave(par::AbstractArray,pinfo::Dict,yr::Integer,iroot::AbstractString,lvl="sfc");

    @load joinpath("$(iroot["data"])","isca_info.jld2") ncattr;
    lon = ncattr["longitude"]; nlon  = length(lon)
    lat = ncattr["latitude"];  nlat  = length(lat);
    pre = ncattr["fullpre"];   issfc = iscaissfc(lvl);
    dim = size(par);           nt    = dim[3];
    tspan = ncattr["timespan"]; pname = pinfo["name"];

    pdir = joinpath(root["data"],pname);

    if issfc
          pfnc = "$(pname)-$(tspan)$(sprintf1("%02d",yr)).nc"
          pdir = joinpath(root["data"],pname);
    else; pfnc = "$(pname)-lvl$(sprintf1("%02d",lvl))-$(tspan)$(sprintf1("%02d",yr)).nc"
          pdir = joinpath(root["data"],pname,"$(pname)-lvl$(sprintf1("%02d",lvl))");
    end

    if !isdir(pdir); mkpath(pdir); end
    if isfile(pfnc)
        @info "$(Dates.now()) - Unfinished netCDF file $(pfnc) detected.  Deleting."
        rm(pfnc);
    end

    var_par = pinfo["name"]; att_par = Dict("units"=>pinfo["unit"]);
    var_lon = "lon";         att_lon = Dict("units"=>"degree");
    var_lat = "lat";         att_lat = Dict("units"=>"degree");
    if issfc; var_pre = "pre"; att_lvl = Dict("units"=>"hPa"); end

    @info "$(Dates.now()) - Creating netCDF file $(pfnc) ..."
    nccreate(pfnc,var_par,"nlon",nlon,"nlat",nlat,"t",nt,atts=att_par,t=NC_FLOAT);
    nccreate(pfnc,var_lon,"nlon",nlon,atts=att_lon,t=NC_FLOAT);
    nccreate(pfnc,var_lat,"nlat",nlat,atts=att_lat,t=NC_FLOAT);
    if issfc; nccreate(pfnc,var_pre,"npre",1,atts=att_lvl,t=NC_FLOAT); end

    @info "$(Dates.now()) - Saving resorted Isca data to netCDF file $(pfnc) ..."
    ncwrite(par,pfnc,var_par);
    ncwrite(lon,pfnc,var_lon);
    ncwrite(lat,pfnc,var_lat);
    if issfc; ncwrite([pre[lvl]],pfnc,var_pre); end

    @debug "$(Dates.now()) - NetCDF.jl's ncread causes memory leakage.  Using ncclose() as a workaround."
    ncclose()

    @info "$(Dates.now()) - Moving $(pfnc) to data directory $(pdir)"
    if isfile(joinpath(pdir,pfnc));
        @info "$(Dates.now()) - An older version of $(pfnc) exists in the $(pdir) directory.  Overwriting."
    end
    mv(pfnc,joinpath(pdir,pfnc),force=true);

end

function iscaparncsave(par::AbstractArray,pinfo::Dict,ncattr::Dict,
                       yr::Integer,iroot::AbstractString,lvl="sfc");

    lon = ncattr["longitude"]; nlon  = length(lon)
    lat = ncattr["latitude"];  nlat  = length(lat);
    pre = ncattr["fullpre"];   issfc = iscaissfc(lvl);
    dim = size(par);           nt    = dim[3];
    tspan = ncattr["timespan"]; pname = pinfo["name"];

    pdir = joinpath(root["data"],pname);

    if issfc
         pfnc = "$(pname)-$(tspan)$(sprintf1("%02d",yr)).nc"
         pdir = joinpath(root["data"],pname);
    else; pfnc = "$(pname)-lvl$(sprintf1("%02d",lvl))-$(tspan)$(sprintf1("%02d",yr)).nc"
         pdir = joinpath(root["data"],pname,"$(pname)-lvl$(sprintf1("%02d",lvl))");
    end

    if !isdir(pdir); mkpath(pdir); end
    if isfile(pfnc)
       @info "$(Dates.now()) - Unfinished netCDF file $(pfnc) detected.  Deleting."
       rm(pfnc);
    end

    var_par = pinfo["name"]; att_par = Dict("units"=>pinfo["unit"]);
    var_lon = "lon";         att_lon = Dict("units"=>"degree");
    var_lat = "lat";         att_lat = Dict("units"=>"degree");
    if issfc; var_pre = "pre"; att_lvl = Dict("units"=>"hPa"); end

    @info "$(Dates.now()) - Creating netCDF file $(pfnc) ..."
    nccreate(pfnc,var_par,"nlon",nlon,"nlat",nlat,"t",nt,atts=att_par,t=NC_FLOAT);
    nccreate(pfnc,var_lon,"nlon",nlon,atts=att_lon,t=NC_FLOAT);
    nccreate(pfnc,var_lat,"nlat",nlat,atts=att_lat,t=NC_FLOAT);
    if issfc; nccreate(pfnc,var_pre,"npre",1,atts=att_lvl,t=NC_FLOAT); end

    @info "$(Dates.now()) - Saving resorted Isca data to netCDF file $(pfnc) ..."
    ncwrite(par,pfnc,var_par);
    ncwrite(lon,pfnc,var_lon);
    ncwrite(lat,pfnc,var_lat);
    if issfc; ncwrite(pre[lvl],pfnc,var_pre); end

    @debug "$(Dates.now()) - NetCDF.jl's ncread causes memory leakage.  Using ncclose() as a workaround."
    ncclose()

    @info "$(Dates.now()) - Moving $(pfnc) to data directory $(pdir)"
    if isfile(joinpath(pdir,pfnc));
       @info "$(Dates.now()) - An older version of $(pfnc) exists in the $(pdir) directory.  Overwriting."
    end
    mv(pfnc,joinpath(pdir,pfnc),force=true);

end


text:
using Test: @test, @testset
import Random
import Statistics
import Zygote
import RestrictedBoltzmannMachines as RBMs
import WeightNormalizedRBMs as wRBMs

@testset "weight_norms" begin
    rbm = RBMs.BinaryRBM(randn(3,2,4), randn(2,5), randn(3,2,4,2,5))
    wn = wRBMs.weight_norms(rbm)
    @test size(wn) == (1,1,1,2,5)
    @test wn ≈ sqrt.(sum(abs2.(rbm.w); dims=1:3))
end

@testset "WeightNormRBM" begin
    rbm = RBMs.BinaryRBM(randn(3,2,4), randn(2,5), randn(3,2,4,2,5))
    un = rand(1,1,1,2,5)
    wrbm = wRBMs.WeightNormRBM(rbm, un)
    rbm1 = RBMs.RBM(wrbm)
    @test rbm1.w ≈ rbm.w
    @test rbm1.visible.θ ≈ rbm.visible.θ
    @test rbm1.hidden.θ ≈ rbm.hidden.θ
end

@testset "∂free_energy" begin
    rbm = RBMs.BinaryRBM(randn(3,2,4), randn(2,5), randn(3,2,4,2,5))
    wrbm = wRBMs.WeightNormRBM(rbm)
    v = Random.bitrand(3,2,4,5,2)
    ∂ = RBMs.∂free_energy(wrbm, v)
    gs, = Zygote.gradient(wrbm) do wrbm
        w = wrbm.g .* wrbm.u ./ sqrt.(sum(abs2, wrbm.u; dims=1:ndims(rbm.visible)))
        return Statistics.mean(RBMs.free_energy(RBMs.RBM(wrbm.visible, wrbm.hidden, w), v))
    end
    @test ∂.g ≈ gs.g
    @test ∂.u ≈ gs.u
    @test ∂.visible.θ ≈ gs.visible.θ
    @test ∂.hidden.θ ≈ gs.hidden.θ
end


text:
using Test

using Catlab.Present

# Schema
########

A, B = Ob(FreeSchema, :A), Ob(FreeSchema, :B)
f, g = Hom(:f, A, B), Hom(:g, B, A)

C = AttrType(FreeSchema.AttrType, :C)
x, y = Attr(:x, A, C), Attr(:y, B, C)

@test dom(x) == A
@test codom(x) == C
@test dom(compose(f,y)) == A
@test dom(compose(g,x)) == B
@test compose(f, compose(g,f)) == compose(compose(f,g), f)
@test compose(f, compose(g, x)) == compose(compose(f, g), x)
@test compose(id(A), f) == f
@test compose(id(A), x) == x


text:
module BackwardKernel

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

using Distributions, SpecialFunctions
import DataStructures: PriorityQueue, peek, enqueue!, dequeue!
import Base: iterate, length, sum, eltype

export EventQueue, sum, length, eltype, iterate

mutable struct EventQueue{ℜ<:Real}
    T::Int
    X::Array{ℜ,2}
    pq::PriorityQueue{Tuple{Int,Int},ℜ,Base.Order.ForwardOrdering}
    function EventQueue(X_new::Array{ℜ,2}) where {ℜ<:Real}
        I::Int = size(X_new,1)
        T_new::Int = Int(round(sum(X_new)))

        pq_new = PriorityQueue([(i₁,i₂) => rand(Dirichlet([1.0,X_new[i₁,i₂]]))[1]
                                           for i₁=1:I, i₂=1:I if X_new[i₁,i₂]>0])
        return new{ℜ}(T_new,copy(X_new),pq_new)
    end
end

Base.sum(L::EventQueue) = sum(L.X)
Base.length(L::EventQueue) = L.T
Base.eltype(L::EventQueue) = Tuple{Int,Int}

function Base.iterate(L::EventQueue{ℜ}, state::Ƶ=1) where {Ƶ<:Int,ℜ<:Real}
    τ::Ƶ = state
    i₁_τ::Ƶ, i₂_τ::Ƶ = 0, 0
    t::ℜ, t_next::ℜ = 0.0, 0.0
    if L.T < τ
        return nothing
    end
    t = peek(L.pq)[2]
    i₁_τ, i₂_τ = dequeue!(L.pq)
    L.X[i₁_τ, i₂_τ] -= 1.0
    if L.X[i₁_τ, i₂_τ] > 0.0
        t_next = t + (1.0 - t)*rand(Dirichlet([1.0, L.X[i₁_τ,i₂_τ]]))[1]
        L.pq[(i₁_τ, i₂_τ)] = t_next
    end
    return (i₁_τ, i₂_τ), state+1
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
# Factorize two different matrices (one real, one complex) and solve for multiple rhs and free memory

using MUMPSjInv
using Test
using LinearAlgebra
using SparseArrays

include("getDivGrad.jl");

A  = getDivGrad(24,23,25);
A2 = getDivGrad(34,32,36);

n  = size(A,1);
n2 = size(A2,1);
A  = A + im*spdiagm(0=>rand(n));

nrhs = 10;
rhs = randn(n,nrhs) + im*randn(n,nrhs);
rhs2 = randn(n2,nrhs);

println("Factorize complex matrix")
@time begin
F1 = factorMUMPS(A,1);
end

println("Factorize real matrix")
@time begin
	F2 = factorMUMPS(A2,1);
end

println("Solve complex system")
@time begin
	x = applyMUMPS(F1,rhs);
end
err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A*x[:,i]-rhs[:,i]) / norm(rhs[:,i]);
end
@test maximum(err) < 1e-14

println("Solve real system")
@time begin
	x2 = applyMUMPS(F2,rhs2);
end
err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A2*x2[:,i]-rhs2[:,i]) / norm(rhs2[:,i]);
end
@test maximum(err) < 1e-14

println("Free memory")
destroyMUMPS(F1);
destroyMUMPS(F2);

println("DONE!")


text:
immutable BinomialRmathSampler <: Sampleable{Univariate,Discrete}
    n::Int
    prob::Float64
end

rand(s::BinomialRmathSampler) =
    round(Int, ccall((:rbinom, "libRmath-julia"), Float64, (Float64, Float64), s.n, s.prob))

# compute probability vector of a Binomial distribution
function binompvec(n::Int, p::Float64)
    pv = Array(Float64, n+1)
    if p == 0.0
        fill!(pv, 0.0)
        pv[1] = 1.0
    elseif p == 1.0
        fill!(pv, 0.0)
        pv[n+1] = 1.0
    else
        q = 1.0 - p
        a = p / q
        @inbounds pv[1] = pk = q ^ n
        for k = 1:n
            @inbounds pv[k+1] = (pk *= ((n - k + 1) / k) * a)
        end
    end
    return pv
end

# Geometric method:
#
#   Devroye. L.
#   "Generating the maximum of independent identically  distributed random variables"
#   Computers and Marhemafics with Applicalions 6, 1960, 305-315.
#
immutable BinomialGeomSampler <: Sampleable{Univariate,Discrete}
    comp::Bool
    n::Int
    scale::Float64
end

BinomialGeomSampler() = BinomialGeomSampler(false, 0, 0.0)

function BinomialGeomSampler(n::Int, prob::Float64)
    if prob <= 0.5
        comp = false
        scale = -1.0/log1p(-prob)
    else
        comp = true
        scale = prob < 1.0 ? -1.0/log(prob) : Inf
    end
    BinomialGeomSampler(comp, n, scale)
end

function rand(s::BinomialGeomSampler)
    y = 0
    x = 0
    n = s.n
    while true
        er = randexp()
        v = er * s.scale
        if v > n  # in case when v is very large or infinity
            break
        end
        y += ceil(Int,v)
        if y > n
            break
        end
        x += 1
    end
    (s.comp ? s.n - x : x)::Int
end

# BTPE algorithm from:
#
#   Kachitvichyanukul, V.; Schmeiser, B. W.
#   "Binomial random variate generation."
#   Comm. ACM 31 (1988), no. 2, 216–222.
#
# Note: only use this sampler when n * min(p, 1-p) is large enough
#       e.g., it is greater than 20.
#
immutable BinomialTPESampler <: Sampleable{Univariate,Discrete}
    comp::Bool
    n::Int
    r::Float64
    q::Float64
    nrq::Float64
    M::Float64
    Mi::Int
    p1::Float64
    p2::Float64
    p3::Float64
    p4::Float64
    xM::Float64
    xL::Float64
    xR::Float64
    c::Float64
    λL::Float64
    λR::Float64
end

BinomialTPESampler() =
    BinomialTPESampler(false, 0, 0., 0., 0., 0., 0,
                       0., 0., 0., 0., 0., 0., 0., 0., 0., 0.)

function BinomialTPESampler(n::Int, prob::Float64)
    if prob <= 0.5
        comp = false
        r = prob
        q = 1.0 - prob
    else
        comp = true
        r = 1.0 - prob
        q = prob
    end

    nrq = n*r*q
    fM = (n+1)*r #
    M = floor(fM)
    Mi = round(Integer, M)
    p1 = floor(2.195*sqrt(nrq)-4.6*q) + 0.5
    xM = M+0.5
    xL = xM-p1
    xR = xM+p1
    c = 0.134 + 20.5/(15.3+M)
    a = (fM-xL)/(fM-xL*r) #
    λL = a*(1.0 + 0.5*a)
    a = (xR-fM)/(xR*q) #
    λR = a*(1.0 + 0.5*a)
    p2 = p1*(1.0 + 2.0*c)
    p3 = p2 + c/λL
    p4 = p3 + c/λR

    BinomialTPESampler(comp,n,r,q,nrq,M,Mi,p1,p2,p3,p4,
                        xM,xL,xR,c,λL,λR)
end

function rand(s::BinomialTPESampler)
    y = 0
    while true
        # Step 1
        u = s.p4*rand()
        v = rand()
        if u <= s.p1
            y = floor(Int,s.xM-s.p1*v+u)
            # Goto 6
            break
        elseif u <= s.p2 # Step 2
            x = s.xL + (u-s.p1)/s.c
            v = v*s.c+1.0-abs(s.M-x+0.5)/s.p1
            if v > 1
                # Goto 1
                continue
            end
            y = floor(Int,x)
            # Goto 5
        elseif u <= s.p3 # Step 3
            y = floor(Int,s.xL + log(v)/s.λL)
            if y < 0
                # Goto 1
                continue
            end
            v *= (u-s.p2)*s.λL
            # Goto 5
        else # Step 4
            y = floor(Int,s.xR-log(v)/s.λR)
            if y > s.n
                # Goto 1
                continue
            end
            v *= (u-s.p3)*s.λR
            # Goto 5
        end

        # Step 5
        # 5.0
        k = abs(y-s.Mi)
        if (k <= 20) || (k >= 0.5*s.nrq-1)
            # 5.1
            S = s.r/s.q
            a = S*(s.n+1)
            F = 1.0
            if s.Mi < y
                for i = (s.Mi+1):y
                    F *= a/i-S
                end
            elseif s.Mi > y
                for i = (y+1):s.Mi
                    F /= a/i-S
                end
            end
            if v > F
                # Goto 1
                continue
            end
            # Goto 6
            break
        else
            # 5.2
            ρ = (k/s.nrq)*((k*(k/3.0+0.625)+1.0/6.0)/s.nrq+0.5)
            t = -k^2/(2.0*s.nrq)
            A = log(v)
            if A < t - ρ
                # Goto 6
                break
            elseif A > t + ρ
                # Goto 1
                continue
            end

            # 5.3
            x1 = Float64(y+1)
            f1 = Float64(s.Mi+1)
            z = Float64(s.n+1-s.Mi)
            w = Float64(s.n-y+1)

            if A > (s.xM*log(f1/x1) + ((s.n-s.Mi)+0.5)*log(z/w) + (y-s.Mi)*log(w*s.r/(x1*s.q)) +
                    lstirling_asym(f1) + lstirling_asym(z) + lstirling_asym(x1) + lstirling_asym(w))
                # Goto 1
                continue
            end

            # Goto 6
            break
        end
    end
    # 6
    (s.comp ? s.n - y : y)::Int
end

# Constructing an alias table by directly computing the probability vector
#
immutable BinomialAliasSampler <: Sampleable{Univariate,Discrete}
    table::AliasTable
end

function BinomialAliasSampler(n::Int, p::Float64)
    pv = binompvec(n, p)
    alias = Array(Int, n+1)
    StatsBase.make_alias_table!(pv, 1.0, pv, alias)
    BinomialAliasSampler(AliasTable(pv, alias, RandIntSampler(n+1)))
end

rand(s::BinomialAliasSampler) = rand(s.table) - 1

# Integrated Polyalgorithm sampler that automatically chooses the proper one
#
# It is important for type-stability
#
type BinomialPolySampler <: Sampleable{Univariate,Discrete}
    use_btpe::Bool
    geom_sampler::BinomialGeomSampler
    btpe_sampler::BinomialTPESampler
end

function BinomialPolySampler(n::Int, p::Float64)
    q = 1.0 - p
    if n * min(p, q) > 20
        use_btpe = true
        geom_sampler = BinomialGeomSampler()
        btpe_sampler = BinomialTPESampler(n, p)
    else
        use_btpe = false
        geom_sampler = BinomialGeomSampler(n, p)
        btpe_sampler = BinomialTPESampler()
    end
    BinomialPolySampler(use_btpe, geom_sampler, btpe_sampler)
end

BinomialPolySampler(n::Real, p::Real) = BinomialPolySampler(round(Int, n), Float64(p))

rand(s::BinomialPolySampler) = s.use_btpe ? rand(s.btpe_sampler) : rand(s.geom_sampler)


text:
struct SplitMeanDiff
    mindiff::Float64
end
SplitMeanDiff() = SplitMeanDiff(0.2)
struct SplitRSS
    minrss::Float64
    diffinstead::Bool
end
SplitRSS(minrss) = SplitRSS(minrss, false)
""" Transform grayscale to binary.
    Use gaussian blur to obtain "background".
    For each pixel find a change point in time.

    Not doing Image alignment/image registration.
    """
function do_split end

"This function is only used for explanation purposes"
function preprocess(A::AbstractArray{T,3}, blurradius=100, ::Type{S}=Float64) where {T,S}
    imgbg = estimatebackground(A, blurradius)
    #Ac = cat(imgbg, imgbg, imgbg, A, dims=3);
    Acc = Array{S, 3}(undef, size(A)[1:2]..., 3+size(A, 3))
    Acc[:,:,1:3] .= oneunit(S)
    Acc[:,:,4:end] .= S.(clamp.(A ./ imgbg, 0.0, 1.0))
    #Acc = cat(ones(T, size(imgbg)..., 3), T.(clamp.(A ./ imgbg, 0.0, 1.0)), dims=3);
    #Acc = T.(Ac ./ imgbg);
    Acc
end

function estimatebackground(A, blurradius)
    imfilter(A[:,:,1], Kernel.gaussian(blurradius));
end
function find_changepoints(m, A::AbstractArray{T,3}, imgbg) where T
    S = Float64
    s3 = size(A, 3) + 3
    Acc = Array{S, 1}(undef, s3)
    Acc[1:3] .= 1.0
    function f(x, y)
        Acc[4:end] .= clamp.(@view(A[x, y, :]) ./ imgbg[x, y], 0.0, 1.0)
        findsplitforpixel(m, Acc)
    end
    [f(x, y) for x in axes(A, 1), y in axes(A, 2)]
end
function do_split(A; blurradius=100, mindiff=0.2)
    do_split(SplitMeanDiff(mindiff), A, blurradius)
end
function do_split(method, A :: AbstractArray{<:Color}, blurradius::Real)
    bg = estimatebackground(A, blurradius)
    Asplit = find_changepoints(method, A, bg)
    do_split(method, A, Asplit)
end
function do_split(method, A, Asplit :: AbstractArray{<:Tuple})
    broadcast(split2binary, Ref(method), Asplit, reshape(1:size(A, 3), (1,1,:)));
end
function find_split(method, A, blurradius=100)
    find_changepoints(method, A, blurradius)
end

# function find_split_threshold(A, blurradius)
#     #Acc = preprocess(A, blurradius)
#     #@time Asplit = mapslices(findsplitforpixel, Acc, dims=3);
#     @time Asplit = find_Asplit(A, blurradius)
#     #A3 = broadcast(dosplit2, A, Asplit, reshape(1:size(A, 3), (1,1,size(A,3))), mindiff);
#     A3 = broadcast(split2threshold, Asplit, reshape(1:size(A, 3), (1,1,size(A,3))));
#     A3
# end

function findsplitrss(x)
    x0 = 0.0
    x1 = Float64(sum(x))
    xx0 = 0.0
    xx1 = Float64(sum(abs2, x))
    n = length(x)
    nullrss = xx1 - x1^2/n

    minss = nullrss
    mini = 0
    meandiff = 0.0
    for i in 1:(n-1)
        xi = Float64(x[i])
        x0 += xi
        x1 -= xi
        xx0 += xi^2
        xx1 -= xi^2
        ss = xx0 - x0^2/i + xx1 - x1^2/(n-i)
        if ss < minss && x0/i > x1/(n-i)
            minss = ss
            meandiff = x0/i - x1/(n-i)
            mini = i
        end
    end
    meandiff, minss, nullrss, mini
end

function findsplitforpixel(m::SplitRSS, x)
    meandiff, rss, nullrss, mini = findsplitrss(x)
    if m.diffinstead
        mini, meandiff
    else
        mini, (nullrss-rss)/(rss/(length(x)-2))
    end
end

function split2binary(m::SplitRSS, split, j)
    i, diff = split
    i -= 3
    if j > i && diff > m.minrss #&& i != 1
        true
    else
        false
    end
end
function split2threshold(split, j)
    i, diff = split
    i -= 3
    if j > i
        diff
    else
        zero(diff)
    end
end

""" Difference in means and sum of variances
"""
function meanandvar(x0, x1, xx0, xx1, i, n)::Tuple{Float64, Float64}
    mean = x0 / i - x1 / (n-i)
    var = (xx0/i - (x0/i)^2) + (xx1/(n-i) - (x1/(n-i))^2)
         # var(x[1:i], corrected=false) + var(x[i+1:n], corrected=false)
    #var = (xx0*i - x0^2)/i^2 + (xx1*(n-i) - x1^2)/(n-i)^2
    #mean = x0 * (n-i) - x1 * i
    #var = (xx0*i - x0^2)*(n-i)^2 + (xx1*(n-i) - x1^2)*i^2
    mean, var
end

function findsplitvar(x)
    x0 = 0.0
    x1 = Float64(sum(x))
    xx0 = 0.0
    xx1 = Float64(sum(abs2, x))
    n = length(x)

    minvar = xx1
    mindiff = x1
    mini = 0
    for i in 1:(n-1)
        xi = Float64(x[i])
        x0 += xi
        x1 -= xi
        xx0 += xi^2
        xx1 -= xi^2
        newdif, newvar = meanandvar(x0, x1, xx0, xx1, i, n)
        if newvar < minvar
            minvar = newvar
            mindiff = newdif
            mini = i
        end
    end
    mini, mindiff, minvar
end

function findsplitforpixel(::SplitMeanDiff, x)
    i, mindiff = findsplitvar(x)
    i, mindiff
end

function split2binary(m::SplitMeanDiff, split, j)
    i, diff = split
    i -= 3
    if j > i && diff > m.mindiff #&& i != 1
        true
    else
        false
    end
end


text:
print("Testing ABCDQR_simula: ")
###################################
nt = 1000
m = ABCDQR_simula(nt)

y = m["y"]
A = m["A"]
C = m["C"]
Q = m["Q"]
R = m["R"]
x10 = 0*A[:,1]
P10 = 0*A

@test sum(y.^2) > 1e-6
####################################
println("OK")


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

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

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


text:
#=
---------------------------------------------------------------------
This function computes the Flux integral contribution to the RHS vector
for a unified CGDG method for the 2D Wave Equation.

Written by F.X. Giraldo on July 14, 2021
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

function create_rhs_flux!(rhs,q,u,face,normals,jac_face,ωq,mapL,mapR,intma,periodicity,DFloat)

    #Get lengths of arrays
    Np=size(intma,1)
    Nface=size(face,2)

    #local arrays
    ql=zeros(DFloat,Np,1)
    qr=zeros(DFloat,Np,1)
    ul=zeros(DFloat,Np,1)
    vl=zeros(DFloat,Np,1)

    #Construct Flux Integral contribution
    for f=1:Nface

        #Store Left Side Variables
        el=face[3,f]
        if (el != -6) #periodic bc
            ilocl=face[1,f]
            for l=1:Np
                #Get Pointers
                il=mapL[1,l,ilocl]
                jl=mapL[2,l,ilocl]
                IL=intma[il,jl,el]

                #Left Element
                ql[l]=q[IL]
                ul[l]=u[1,IL]
                vl[l]=u[2,IL]
            end #l

            #Store Right Side Variables
            er=face[4,f]
            if (er > 0)
                ilocr=face[2,f];
                for l=1:Np
                    #Get Pointers
                    ir=mapR[1,l,ilocr]
                    jr=mapR[2,l,ilocr]
                    IR=intma[ir,jr,er]

                    #Right Element
                    qr[l]=q[IR]
                    #ur[l]=u[IR]; #Not needed since U is continuous
                    #vr[l]=v[IR]; #Not needed since V is continuous
                end #l
            end #if er

            #Do Gauss-Lobatto Integration
            for l=1:Np
                wq=ωq[l]*jac_face[l,f]

                #Store Normal Vectors
                nxl=normals[1,l,f]
                nyl=normals[2,l,f]
                nxr=-nxl
                nyr=-nyl

                #Interpolate onto Quadrature Points
                qlq_k=ql[l]
                qrq_k=qr[l]
                u_k=ul[l]
                v_k=vl[l]
                ul_k=u_k
                vl_k=v_k
                ur_k=u_k
                vr_k=v_k

                #Compute Rusanov flux Constant
                unl=nxl*ul_k + nyl*vl_k
                unr=nxl*ur_k + nyl*vr_k
                claml=abs(unl)
                clamr=abs(unr)
                clam=max(claml,clamr)

                #Flux Variables
                fxl=qlq_k*ul_k
                fyl=qlq_k*vl_k

                fxr=qrq_k*ur_k
                fyr=qrq_k*vr_k

                #Normal Flux Component
                flux_ql=( nxl*fxl + nyl*fyl )
                flux_qr=( nxr*fxr + nyr*fyr )
                flux_q=flux_ql - flux_qr

                #Dissipation Term
                diss_q=clam*(qrq_k - qlq_k)

                #Construct Rusanov Flux
                flux_q=0.5*(flux_q - diss_q)

                #Loop through Side Interpolation Points
                #--------------Left Side------------------#
                il=mapL[1,l,ilocl]
                jl=mapL[2,l,ilocl]
                IL=intma[il,jl,el]
                IL=periodicity[IL]
                rhs[IL] -= wq*flux_q

                #--------------Right Side------------------#
                if (er > 0)
                    ir=mapR[1,l,ilocr]
                    jr=mapR[2,l,ilocr]
                    IR=intma[ir,jr,er]
                    IR=periodicity[IR]
                    rhs[IR] += wq*flux_q
                end #if
            end #l
        end #if el
    end #f

    return(rhs)
end


text:
module Outsource

export outsource

include("connector.jl")

using Distributed: remote_do, workers, myid, remotecall_eval

"""
    outsource(f, id)
    outsource(f, id, T, S)

Remotely issue a processing function `f` to run on worker `id`, optionally
specifying input and output type `T` and `S` respectively.
Returns a `Connector{S,T}` that allows `put!()`ing input data of type `T` and
`take!()`ing output data of type `S`.
The processing function `f` is being called with a `Connector{T,S}` that allows
`take!()`ing input data of type `T` and `put!()`ing output data of type `S`.
Exceptions thrown in `f` are logged as errors in the current task.
"""
outsource(f, id = rand(workers())) = outsource(f, id, Any, Any)
function outsource(f, id, ::Type{T}, ::Type{S}) where {T,S}
    myid() == 1 || @error("need to be called from pid 1")
    id != myid() || @warn("outsourcing to task on same process" *
        " (no parallel execution)", maxlog = 1)

    con = Connector{S, T}()
    rcon = reverse(con)

    # hack to ensure `f` is defined remotely
    # otherwise we get an uncatchable exception that gets hidden
    remotecall_eval(Main, id, f)
    remote_do(id) do
        try
            f(rcon)
        catch e
            # TODO: for non-shared stdout, exceptions need to be
            #       communicated back. Maybe introduce a global pid1
            #       remote-channel for that
            # TODO: prune serialization stuff from stacktrace
            @error(
                "outsourced task failed",
                exception = (e, catch_backtrace()))
        finally
            close(rcon)
        end
    end
    return con
end

end # module


text:
using LibRealSense
using Test

# As the doc says:
# "Note: Due to the USB 3.0 translation layer between native hardware and virtual
#  machine, the librealsense team does not support installation in a VM."
# so we only run tests on MacOS.
get(ENV, "TRAVIS_OS_NAME", "") == "linux" && exit()

err = Ref{Ptr{rs2_error}}(0)

ctx = rs2_create_context(RS2_API_VERSION, err)
@test err[] == C_NULL

device_list = rs2_query_devices(ctx, err)
@test err[] == C_NULL

dev_count = rs2_get_device_count(device_list, err)
@test err[] == C_NULL
println("There are $dev_count connected RealSense devices.")

rs2_delete_device_list(device_list)
rs2_delete_context(ctx)


text:
using Test
using RingLists

a = RingList(1,2,3)
b = RingList(2,3,1)
@test a==b
@test reverse(a)==reverse(b)
@test hash(a) == hash(b)
@test length(a) == 3
@test sort(collect(keys(a))) == [1,2,3]
@test Set(a) == Set([1,2,3])
@test first(a,true) == 1
@test a[1] == 2
@test Vector(b) == [1,2,3]
@test eltype(b) == Int
a = RingList{Int}()
insert!(a,1)
insert!(a,2)
@test Vector(a) == [1,2]

a = RingList(1,2,3,4)
@test next(a,1)==2
@test a[1] == 2
@test previous(a,2) == 1
@test a(2) == 1

@test sort(collect(a)) == [1,2,3,4]
b = RingList(1,2,4)
delete!(a,3)
@test a==b

a = RingList(collect(1:10))
insertafter!(a,99,2)
@test a[2]==99
@test a(99)==2

a = RingList(collect(1:10))
insertbefore!(a,99,2)
@test a(2)==99
@test a[99]==2

a = RingList(collect(1:9))
@test sum(a) == sum(1:9)
v = [x*x for x in a]
@test sum(v) == sum(t^2 for t in 1:9)

a = RingList(3,5,1,2,4)
b = shuffle(a)
@test Set(a)==Set(b)
@test sum(a) == sum(b)


text:
import Base: +, getindex, isempty

export MinkowskiSum, ⊕,
       MinkowskiSumArray,
       MinkowskiSum!,
       array,
       CacheMinkowskiSum,
       forget_sets!

"""
    MinkowskiSum{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}

Type that represents the Minkowski sum of two convex sets.

### Fields

- `X` -- first convex set
- `Y` -- second convex set

### Notes

The `ZeroSet` is the neutral element and the `EmptySet` is the absorbing element
for `MinkowskiSum`.
"""
struct MinkowskiSum{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} <: LazySet{N}
    X::S1
    Y::S2

    # default constructor with dimension check
    function MinkowskiSum{N, S1, S2}(X::S1, Y::S2) where
            {N<:Real, S1<:LazySet{N}, S2<:LazySet{N}}
        @assert dim(X) == dim(Y) "sets in a Minkowski sum must have the " *
            "same dimension"
        return new{N, S1, S2}(X, Y)
    end
end

# convenience constructor without type parameter
MinkowskiSum(X::S1, Y::S2) where {N<:Real, S1<:LazySet{N}, S2<:LazySet{N}} =
    MinkowskiSum{N, S1, S2}(X, Y)

# ZeroSet is the neutral element for MinkowskiSum
@neutral(MinkowskiSum, ZeroSet)

# EmptySet and Universe are the absorbing elements for MinkowskiSum
@absorbing(MinkowskiSum, EmptySet)
# @absorbing(MinkowskiSum, Universe)  # TODO problematic

"""
    X + Y

Convenience constructor for Minkowski sum.

### Input

- `X` -- a convex set
- `Y` -- another convex set

### Output

The symbolic Minkowski sum of ``X`` and ``Y``.
"""
+(X::LazySet, Y::LazySet) = MinkowskiSum(X, Y)

"""
    ⊕(X::LazySet, Y::LazySet)

Unicode alias constructor ⊕ (`oplus`) for the lazy Minkowski sum operator.
"""
⊕(X::LazySet, Y::LazySet) = MinkowskiSum(X, Y)

"""
    dim(ms::MinkowskiSum)::Int

Return the dimension of a Minkowski sum.

### Input

- `ms` -- Minkowski sum

### Output

The ambient dimension of the Minkowski sum.
"""
function dim(ms::MinkowskiSum)::Int
    return dim(ms.X)
end

"""
    σ(d::AbstractVector{N}, ms::MinkowskiSum{N}) where {N<:Real}

Return the support vector of a Minkowski sum.

### Input

- `d`  -- direction
- `ms` -- Minkowski sum

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the summand sets.

### Algorithm

The support vector in direction ``d`` of the Minkowski sum of two sets ``X``
and ``Y`` is the sum of the support vectors of ``X`` and ``Y`` in direction
``d``.
"""
function σ(d::AbstractVector{N}, ms::MinkowskiSum{N}) where {N<:Real}
    return σ(d, ms.X) + σ(d, ms.Y)
end

"""
    ρ(d::AbstractVector{N}, ms::MinkowskiSum{N}) where {N<:Real}

Return the support function of a Minkowski sum.

### Input

- `d`  -- direction
- `ms` -- Minkowski sum

### Output

The support function in the given direction.

### Algorithm

The support function in direction ``d`` of the Minkowski sum of two sets ``X``
and ``Y`` is the sum of the support functions of ``X`` and ``Y`` in direction
``d``.
"""
function ρ(d::AbstractVector{N}, ms::MinkowskiSum{N}) where {N<:Real}
    return ρ(d, ms.X) + ρ(d, ms.Y)
end

"""
    isbounded(ms::MinkowskiSum)::Bool

Determine whether a Minkowski sum is bounded.

### Input

- `ms` -- Minkowski sum

### Output

`true` iff both wrapped sets are bounded.
"""
function isbounded(ms::MinkowskiSum)::Bool
    return isbounded(ms.X) && isbounded(ms.Y)
end

"""
    isempty(ms::MinkowskiSum)::Bool

Return if a Minkowski sum is empty or not.

### Input

- `ms` -- Minkowski sum

### Output

`true` iff any of the wrapped sets are empty.
"""
function isempty(ms::MinkowskiSum)::Bool
    return isempty(ms.X) || isempty(ms.Y)
end

# =================================
# Minkowski sum of an array of sets
# =================================

"""
    MinkowskiSumArray{N<:Real, S<:LazySet{N}} <: LazySet{N}

Type that represents the Minkowski sum of a finite number of convex sets.

### Fields

- `array` -- array of convex sets

### Notes

This type assumes that the dimensions of all elements match.

The `ZeroSet` is the neutral element and the `EmptySet` is the absorbing element
for `MinkowskiSumArray`.

Constructors:

- `MinkowskiSumArray(array::Vector{<:LazySet})` -- default constructor

- `MinkowskiSumArray([n]::Int=0, [N]::Type=Float64)`
  -- constructor for an empty sum with optional size hint and numeric type
"""
struct MinkowskiSumArray{N<:Real, S<:LazySet{N}} <: LazySet{N}
    array::Vector{S}
end

# constructor for an empty sum with optional size hint and numeric type
function MinkowskiSumArray(n::Int=0, N::Type=Float64)::MinkowskiSumArray
    arr = Vector{LazySet{N}}()
    sizehint!(arr, n)
    return MinkowskiSumArray(arr)
end

# ZeroSet is the neutral element for MinkowskiSumArray
@neutral(MinkowskiSumArray, ZeroSet)

# EmptySet and Universe are the absorbing elements for MinkowskiSumArray
@absorbing(MinkowskiSumArray, EmptySet)
# @absorbing(MinkowskiSumArray, Universe)  # TODO problematic

# add functions connecting MinkowskiSum and MinkowskiSumArray
@declare_array_version(MinkowskiSum, MinkowskiSumArray)

"""
    array(msa::MinkowskiSumArray{N, S})::Vector{S} where {N<:Real, S<:LazySet{N}}

Return the array of a Minkowski sum of a finite number of convex sets.

### Input

- `msa` -- Minkowski sum array

### Output

The array of a Minkowski sum of a finite number of convex sets.
"""
function array(msa::MinkowskiSumArray{N, S})::Vector{S} where {N<:Real, S<:LazySet{N}}
    return msa.array
end

"""
    dim(msa::MinkowskiSumArray)::Int

Return the dimension of a Minkowski sum of a finite number of sets.

### Input

- `msa` -- Minkowski sum array

### Output

The ambient dimension of the Minkowski sum of a finite number of sets.
"""
function dim(msa::MinkowskiSumArray)::Int
    return length(msa.array) == 0 ? 0 : dim(msa.array[1])
end

"""
    σ(d::AbstractVector{N}, msa::MinkowskiSumArray{N}) where {N<:Real}

Return the support vector of a Minkowski sum of a finite number of sets in a
given direction.

### Input

- `d`   -- direction
- `msa` -- Minkowski sum array

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the summand sets.
"""
function σ(d::AbstractVector{N}, msa::MinkowskiSumArray{N}) where {N<:Real}
    return σ_helper(d, msa.array)
end

"""
    ρ(d::AbstractVector{N}, msa::MinkowskiSumArray{N}) where {N<:Real}

Return the support function of a Minkowski sum array of a finite number of sets
in a given direction.

### Input

- `d`   -- direction
- `msa` -- Minkowski sum array

### Output

The support function in the given direction.

### Algorithm

The support function of the Minkowski sum of sets is the sum of the support
functions of each set.
"""
function ρ(d::AbstractVector{N}, msa::MinkowskiSumArray{N}) where {N<:Real}
    return sum([ρ(d, Xi) for Xi in msa.array])
end

"""
	isbounded(msa::MinkowskiSumArray)::Bool

Determine whether a Minkowski sum of a finite number of convex sets is bounded.

### Input

- `msa` -- Minkowski sum of a finite number of convex sets

### Output

`true` iff all wrapped sets are bounded.
"""
function isbounded(msa::MinkowskiSumArray)::Bool
    return all(x -> isbounded(x), msa.array)
end

"""
    isempty(msa::MinkowskiSumArray)::Bool

Return if a Minkowski sum array is empty or not.

### Input

- `msa` -- Minkowski sum array

### Output

`true` iff any of the wrapped sets are empty.
"""
function isempty(msa::MinkowskiSumArray)::Bool
    return any(X -> isempty(X), array(msa))
end

# =============================================================
# Minkowski sum of an array of sets with a support vector cache
# =============================================================

"""
    CachedPair{N}

A mutable pair of an index and a vector.

### Fields

- `idx` -- index
- `vec` -- vector
"""
mutable struct CachedPair{N}
    idx::Int
    vec::AbstractVector{N}
end

function getindex(cp::CachedPair, idx::Int)
    if idx == 1
        return cp.idx
    elseif idx == 2
        return cp.vec
    end
    error("invalid index $idx, can only access a pair at index 1 or 2")
end

"""
    CacheMinkowskiSum{N<:Real, S<:LazySet{N}} <: LazySet{N}

Type that represents the Minkowski sum of a finite number of convex sets.
Support vector queries are cached.

### Fields

- `array` -- array of convex sets
- `cache` -- cache of support vector query results

### Notes

This type assumes that the dimensions of all elements match.

The `ZeroSet` is the neutral element and the `EmptySet` is the absorbing element
for `CacheMinkowskiSum`.

The cache (field `cache`) is implemented as dictionary whose keys are directions
and whose values are pairs `(k, s)` where `k` is the number of elements in the
array `array` when the support vector was evaluated last time, and `s` is the
support vector that was obtained.
Thus this type assumes that `array` is not modified except by adding new sets at
the end.

Constructors:

- `CacheMinkowskiSum(array::Vector{<:LazySet})` -- default constructor

- `CacheMinkowskiSum([n]::Int=0, [N]::Type=Float64)`
  -- constructor for an empty sum with optional size hint and numeric type
"""
struct CacheMinkowskiSum{N<:Real, S<:LazySet{N}} <: LazySet{N}
    array::Vector{S}
    cache::Dict{AbstractVector{N}, CachedPair{N}}

    # default constructor that initializes cache
    CacheMinkowskiSum{N, S}(arr::Vector{S}) where {N<:Real, S<:LazySet{N}} =
        new{N, S}(arr, Dict{AbstractVector{N}, CachedPair{N}}())
end

# convenience constructor without type parameter
CacheMinkowskiSum(arr::Vector{S}) where {N<:Real, S<:LazySet{N}} =
    CacheMinkowskiSum{N, S}(arr)

# constructor for an empty sum with optional size hint and numeric type
function CacheMinkowskiSum(n::Int=0, N::Type=Float64)::CacheMinkowskiSum
    arr = Vector{LazySet{N}}()
    sizehint!(arr, n)
    return CacheMinkowskiSum(arr)
end

# ZeroSet is the neutral element for CacheMinkowskiSum
@neutral(CacheMinkowskiSum, ZeroSet)

# EmptySet and Universe are the absorbing element for CacheMinkowskiSum
@absorbing(CacheMinkowskiSum, EmptySet)
# @absorbing(CacheMinkowskiSum, Universe)  # TODO problematic

"""
    array(cms::CacheMinkowskiSum{N, S})::Vector{S} where {N<:Real, S<:LazySet{N}}

Return the array of a caching Minkowski sum.

### Input

- `cms` -- caching Minkowski sum

### Output

The array of a caching Minkowski sum.
"""
function array(cms::CacheMinkowskiSum{N, S}
              )::Vector{S} where {N<:Real, S<:LazySet{N}}
    return cms.array
end

"""
    dim(cms::CacheMinkowskiSum)::Int

Return the dimension of a caching Minkowski sum.

### Input

- `cms` -- caching Minkowski sum

### Output

The ambient dimension of the caching Minkowski sum.
"""
function dim(cms::CacheMinkowskiSum)::Int
    return length(cms.array) == 0 ? 0 : dim(cms.array[1])
end

"""
    σ(d::AbstractVector{N}, cms::CacheMinkowskiSum{N}) where {N<:Real}

Return the support vector of a caching Minkowski sum in a given direction.

### Input

- `d`   -- direction
- `cms` -- caching Minkowski sum

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the summand sets.

### Notes

The result is cached, i.e., any further query with the same direction runs in
constant time.
When sets are added to the caching Minkowski sum, the query is only performed
for the new sets.
"""
function σ(d::AbstractVector{N}, cms::CacheMinkowskiSum{N}) where {N<:Real}
    arr = array(cms)
    l = length(arr)
    cache = cms.cache
    if haskey(cache, d)
        pair = cache[d]
        k = pair.idx
        svec1 = pair.vec
        if k == l
            # has already stored the full support vector
            return svec1
        else
            # has only stored the support vector of the first k sets
            @assert k < l "invalid cache index"
            svec = svec1 + σ_helper(d, @view arr[k+1:l])
        end
    else
        # first-time computation of support vector
        svec = σ_helper(d, arr)
    end
    # NOTE: make a copy of the direction vector (can be modified outside)
    cache[copy(d)] = CachedPair(l, svec)
    return svec
end

"""
	isbounded(cms::CacheMinkowskiSum)::Bool

Determine whether a caching Minkowski sum is bounded.

### Input

- `cms` -- caching Minkowski sum

### Output

`true` iff all wrapped sets are bounded.
"""
function isbounded(cms::CacheMinkowskiSum)::Bool
    return all(x -> isbounded(x), cms.array)
end

"""
    isempty(cms::CacheMinkowskiSum)::Bool

Return if a caching Minkowski sum array is empty or not.

### Input

- `cms` -- caching Minkowski sum

### Output

`true` iff any of the wrapped sets are empty.

### Notes

Forgotten sets cannot be checked anymore.
Usually they have been empty because otherwise the support vector query should
have crashed before.
In that case, the caching Minkowski sum should not be used further.
"""
function isempty(cms::CacheMinkowskiSum)::Bool
    return any(X -> isempty(X), array(cms))
end

"""
    forget_sets!(cms::CacheMinkowskiSum)::Int

Tell a caching Minkowski sum to forget the stored sets (but not the support
vectors).
Only those sets are forgotten such that for each cached direction the support
vector has been computed before.

### Input

- `cms` -- caching Minkowski sum

### Output

The number of sets that have been forgotten.

### Notes

This function should only be used under the assertion that no new directions are
queried in the future; otherwise such support vector results will be incorrect.

This implementation is optimistic and first tries to remove all sets.
However, it also checks that for all cached directions the support vector has
been computed before.
If it finds that this is not the case, the implementation identifies the
biggest index ``k`` such that the above holds for the ``k`` oldest sets, and
then it only removes these.
See the example below.

### Examples

```jldoctest
julia> x1 = BallInf(ones(3), 3.); x2 = Ball1(ones(3), 5.);

julia> cms1 = CacheMinkowskiSum(2); cms2 = CacheMinkowskiSum(2);

julia> d = ones(3);

julia> a1 = array(cms1); a2 = array(cms2);

julia> push!(a1, x1); push!(a2, x1);

julia> σ(d, cms1); σ(d, cms2);

julia> push!(a1, x2); push!(a2, x2);

julia> σ(d, cms1);

julia> idx1 = forget_sets!(cms1) # support vector was computed for both sets
2

julia> idx1 = forget_sets!(cms2) # support vector was only computed for first set
1
```
"""
function forget_sets!(cms::CacheMinkowskiSum)::Int
    len = length(cms.array)
    sets_to_remove = len
    for pair in values(cms.cache)
        if pair.idx != len
            @assert pair.idx < len "invalid cache index"
            sets_to_remove = min(sets_to_remove, pair.idx)
            pair.idx -= len
        else
            pair.idx = 0
        end
    end

    if sets_to_remove != len
        # not all sets shall be forgotten
        offset = len - sets_to_remove
        for pair in values(cms.cache)
            # fix the wrong index assignments from above
            pair.idx += offset
        end
        # forget only the first k sets
        splice!(cms.array, 1:sets_to_remove)
    else
        # forget all sets
        empty!(cms.array)
    end
    return sets_to_remove
end

# ================
# Helper functions
# ================

@inline function σ_helper(d::AbstractVector{N},
                          array::AbstractVector{<:LazySet}
                         )::Vector{N} where {N<:Real}
    svec = zeros(N, length(d))
    for sj in array
        svec += σ(d, sj)
    end
    return svec
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct NumericGroup{T<:AbstractConcept}
    owner::T
    numeric::Number
end

function NumericGroup(res::Proto.NumericGroup)
    owner = instantiate(res.owner)
    numeric = _read_proto_number(res.number)
    return NumericGroup(owner, numeric)
end

function NumericGroup(res::Vector{Proto.NumericGroup})
    return NumericGroup.(res)
end


text:
include("utils.jl")

struct TestVectorSpaceType <: VectorSpaceType end

@testset "Tangent bundle" begin
    M = Sphere(2)

    @testset "FVector" begin
        @test sprint(show, TangentSpace) == "TangentSpace"
        @test sprint(show, CotangentSpace) == "CotangentSpace"
        tvs = ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0])
        fv_tvs = map(v -> FVector(TangentSpace, v), tvs)
        fv1 = fv_tvs[1]
        tv1s = allocate(fv_tvs[1])
        @test isa(tv1s, FVector)
        @test size(tv1s) == (3,)
        @test tv1s.type == TangentSpace
        @test size(tv1s.data) == size(tvs[1])
        @test number_eltype(tv1s) == number_eltype(tvs[1])
        @test number_eltype(tv1s) == number_eltype(typeof(tv1s))
        @test isa(fv1 + fv1, FVector)
        @test (fv1 + fv1).type == TangentSpace
        @test isa(fv1 - fv1, FVector)
        @test (fv1 - fv1).type == TangentSpace
        @test isa(-fv1, FVector)
        @test (-fv1).type == TangentSpace
        @test isa(2 * fv1, FVector)
        @test (2 * fv1).type == TangentSpace

        PM = ProductManifold(Sphere(2), Euclidean(2))
        @test_throws ErrorException flat(
            PM,
            ProductRepr([0.0], [0.0]),
            FVector(CotangentSpace, ProductRepr([0.0], [0.0])),
        )
        @test_throws ErrorException sharp(
            PM,
            ProductRepr([0.0], [0.0]),
            FVector(TangentSpace, ProductRepr([0.0], [0.0])),
        )

        fv2 = FVector(TangentSpace, ProductRepr([1, 0, 0], [1 2]))
        @test submanifold_component(fv2, 1) == [1, 0, 0]
        @test submanifold_component(fv2, 2) == [1 2]
        @test (@inferred submanifold_component(fv2, Val(1))) == [1, 0, 0]
        @test (@inferred submanifold_component(fv2, Val(2))) == [1 2]
        @test submanifold_component(PM, fv2, 1) == [1, 0, 0]
        @test submanifold_component(PM, fv2, 2) == [1 2]
        @test (@inferred submanifold_component(PM, fv2, Val(1))) == [1, 0, 0]
        @test (@inferred submanifold_component(PM, fv2, Val(2))) == [1 2]
        @test submanifold_components(PM, fv2) == ([1.0, 0.0, 0.0], [1.0 2.0])
    end

    types = [Vector{Float64}]
    TEST_FLOAT32 && push!(types, Vector{Float32})
    TEST_STATIC_SIZED && push!(types, MVector{3,Float64})

    for T in types
        x = convert(T, [1.0, 0.0, 0.0])
        TB = TangentBundle(M)
        @test sprint(show, TB) == "TangentBundle(Sphere(2, ℝ))"
        @test base_manifold(TB) == M
        @test manifold_dimension(TB) == 2 * manifold_dimension(M)
        @test representation_size(TB) == (6,)
        CTB = CotangentBundle(M)
        @test sprint(show, CTB) == "CotangentBundle(Sphere(2, ℝ))"
        @test sprint(show, VectorBundle(TestVectorSpaceType(), M)) ==
              "VectorBundle(TestVectorSpaceType(), Sphere(2, ℝ))"
        @testset "Type $T" begin
            pts_tb = [
                ProductRepr(convert(T, [1.0, 0.0, 0.0]), convert(T, [0.0, -1.0, -1.0])),
                ProductRepr(convert(T, [0.0, 1.0, 0.0]), convert(T, [2.0, 0.0, 1.0])),
                ProductRepr(convert(T, [1.0, 0.0, 0.0]), convert(T, [0.0, 2.0, -1.0])),
            ]
            @inferred ProductRepr(
                convert(T, [1.0, 0.0, 0.0]),
                convert(T, [0.0, -1.0, -1.0]),
            )
            for pt in pts_tb
                @test bundle_projection(TB, pt) ≈ pt.parts[1]
            end
            diag_basis = DiagonalizingOrthonormalBasis(log(TB, pts_tb[1], pts_tb[2]))
            basis_types = (
                DefaultOrthonormalBasis(),
                get_basis(TB, pts_tb[1], DefaultOrthonormalBasis()),
                diag_basis,
                get_basis(TB, pts_tb[1], diag_basis),
            )
            test_manifold(
                TB,
                pts_tb,
                test_injectivity_radius = false,
                test_reverse_diff = isa(T, Vector),
                test_forward_diff = isa(T, Vector),
                test_tangent_vector_broadcasting = false,
                test_vee_hat = true,
                test_project_tangent = true,
                test_project_point = true,
                test_default_vector_transport = true,
                basis_types_vecs = basis_types,
                projection_atol_multiplier = 4,
            )
        end
    end

    @test TangentBundle{ℝ,Sphere{2,ℝ}} ==
          VectorBundle{ℝ,Manifolds.TangentSpaceType,Sphere{2,ℝ}}
    @test CotangentBundle{ℝ,Sphere{2,ℝ}} ==
          VectorBundle{ℝ,Manifolds.CotangentSpaceType,Sphere{2,ℝ}}

    @test base_manifold(TangentBundle(M)) == M
    @testset "spaces at point" begin
        x = [1.0, 0.0, 0.0]
        t_x = TangentSpaceAtPoint(M, x)
        t_x2 = TangentSpace(M, x)
        @test t_x == t_x2
        ct_x = CotangentSpaceAtPoint(M, x)
        t_xs = sprint(show, "text/plain", t_x)
        sp = sprint(show, "text/plain", x)
        sp = replace(sp, '\n' => "\n ")
        t_xs_test = "Tangent space to the manifold $(M) at point:\n $(sp)"
        @test t_xs == t_xs_test
        @test base_manifold(t_x) == M
        @test base_manifold(ct_x) == M
        @test t_x.fiber.manifold == M
        @test ct_x.fiber.manifold == M
        @test t_x.fiber.fiber == TangentSpace
        @test ct_x.fiber.fiber == CotangentSpace
        @test t_x.point == x
        @test ct_x.point == x
        # generic vector space at
        fiber = VectorBundleFibers(TestVectorSpaceType(), M)
        v_x = VectorSpaceAtPoint(fiber, x)
        v_xs = sprint(show, "text/plain", v_x)
        fiber_s = sprint(show, "text/plain", fiber)
        v_xs_test = "$(typeof(v_x))\nFiber:\n $(fiber_s)\nBase point:\n $(sp)"
        @test v_xs == v_xs_test
    end

    @testset "tensor product" begin
        TT = Manifolds.TensorProductType(TangentSpace, TangentSpace)
        @test sprint(show, TT) == "TensorProductType(TangentSpace, TangentSpace)"
        @test vector_space_dimension(VectorBundleFibers(TT, Sphere(2))) == 4
        @test vector_space_dimension(VectorBundleFibers(TT, Sphere(3))) == 9
        @test base_manifold(VectorBundleFibers(TT, Sphere(2))) == M
        @test sprint(show, VectorBundleFibers(TT, Sphere(2))) ==
              "VectorBundleFibers(TensorProductType(TangentSpace, TangentSpace), Sphere(2, ℝ))"
    end

    @testset "Error messages" begin
        vbf = VectorBundleFibers(TestVectorSpaceType(), Euclidean(3))
        @test_throws ErrorException inner(vbf, [1, 2, 3], [1, 2, 3], [1, 2, 3])
        @test_throws ErrorException Manifolds.project!(vbf, [1, 2, 3], [1, 2, 3], [1, 2, 3])
        @test_throws ErrorException zero_vector!(vbf, [1, 2, 3], [1, 2, 3])
        @test_throws ErrorException vector_space_dimension(vbf)
        a = fill(0.0, 6)
        @test_throws ErrorException get_coordinates!(
            TangentBundle(M),
            a,
            ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0, 0.0]),
            ProductRepr([1.0, 0.0, 0.0], [0.0, 0.0, 1.0]),
            CachedBasis(DefaultOrthonormalBasis(), []),
        )
    end

    @testset "log and exp on tangent bundle for power and product manifolds" begin
        M = PowerManifold(Circle(ℝ), 2)
        N = TangentBundle(M)
        p1 = ProductRepr([0.0, 0.0], [0.0, 0.0])
        p2 = ProductRepr([-1.047, -1.047], [0.0, 0.0])
        X1 = log(N, p1, p2)
        @test isapprox(N, p2, exp(N, p1, X1))
        @test is_tangent_vector(N, p2, vector_transport_to(N, p1, X1, p2))

        M2 = ProductManifold(Circle(ℝ), Euclidean(2))
        N2 = TangentBundle(M2)
        p1_2 = ProductRepr(ProductRepr([0.0], [0.0, 0.0]), ProductRepr([0.0], [0.0, 0.0]))
        p2_2 = ProductRepr(
            ProductRepr([-1.047], [1.0, 0.0]),
            ProductRepr([-1.047], [0.0, 1.0]),
        )
        @test isapprox(N2, p2_2, exp(N2, p1_2, log(N2, p1_2, p2_2)))

        ppt = PowerVectorTransport(ParallelTransport())
        tbvt = Manifolds.VectorBundleVectorTransport(ppt, ppt)
        @test TangentBundle(M, tbvt).vector_transport === tbvt
        @test CotangentBundle(M, tbvt).vector_transport === tbvt
        @test VectorBundle(TangentSpace, M, tbvt).vector_transport === tbvt
    end
end


text:
mutable struct TouchClassifier <: StackedHourglass.NN
    fb::StackedHourglass.FirstBlock
    r1::StackedHourglass.Residual
    pool1::StackedHourglass.Pool
    r2::StackedHourglass.Residual
    pool2::StackedHourglass.Pool
    r3::StackedHourglass.Residual
    pool3::StackedHourglass.Pool
    w::Union{Param{KnetArray{Float32,2}},Param{Array{Float32,2}}}
    b::StackedHourglass.PType4
end

function TouchClassifier(N,N_input,K,atype=KnetArray)
    fb = StackedHourglass.FirstBlock(N,atype) # 64 x 64 x 64

    fb.c1.w = Knet.Param(convert(atype,Knet.xavier_normal(Float32,7,7,N_input,64)))
    fb.c1.bn_p = Knet.Param(convert(atype{Float32,1},Knet.bnparams(N_input)))
    fb.c1.ms = Knet.bnmoments()

    r1 = StackedHourglass.Residual(N,N,atype)
    p1 = StackedHourglass.Pool()
    r2 = StackedHourglass.Residual(N,N,atype) # 32 x 32 x 64
    p2 = StackedHourglass.Pool()
    r3 = StackedHourglass.Residual(N,N,atype) # 16 x 16 x 64
    p3 = StackedHourglass.Pool() # 8 x 8 x 64
    fc_w = Knet.Param(convert(atype,Knet.xavier_normal(Float32,1,4096)))
    fc_b = Knet.Param(convert(atype,Knet.xavier_normal(Float32,1,1,1,1)))
    TouchClassifier(fb,r1,p1,r2,p2,r3,p3,fc_w,fc_b)
end

function (h::TouchClassifier)(x::StackedHourglass.HGType)

    temp = h.fb(x)
    r1 = h.r1(temp)
    p1 = h.pool1(r1)
    r2 = h.r2(p1)
    p2 = h.pool2(r2)
    r3 = h.r3(p2)
    p3 = h.pool3(r3)
    bmm(h.w,mat(p3)) .+ mat(h.b) #Fully connected layer. bmm does batch multiply, while MAT makes the matrix 2d
end

function (h::TouchClassifier)(x::StackedHourglass.HGType,y)
    bce(vec(h(x)),y)
end

function set_testing(h::TouchClassifier,training=false)
    StackedHourglass.set_testing(h.fb,training)
    StackedHourglass.set_testing(h.r1,training)
    StackedHourglass.set_testing(h.r2,training)
    StackedHourglass.set_testing(h.r3,training)
end


text:
module Exodus

using Base
using Suppressor
using Exodus_jll

include("Constants.jl")
include("Types.jl")
include("Errors.jl")
include("IO.jl")
include("Initialization.jl")
include("Coordinates.jl")
include("Times.jl")
include("Blocks.jl")
include("NodeSets.jl")
include("Meshes.jl")
include("Variables.jl")

end # module