text:
function pprint(io::GarishIO, mime::MIME"text/plain", X::AbstractArray)
    # use default array printing if it's the root
    io.state.level == 0 && return show(io, mime, X)
    ndims(X) == 0 && return pprint_zero_dim(io, X)

    if ndims(X) == 1
        # io.compact && return pprint_list_like(io, X)
        # heurostics to print vector in compact form
        if length(X) < 20 && !max_indent_reached(io, length(string(X)))
            pprint_list_like(io, X; compact=true)
        else
            pprint_list_like(io, X)
        end
        return
    end
    # fallback to compact show for high dimensional arrays
    # since we cannot control the indent in multi-line printing
    # not implemented with pprint interface
    if io.compact
        return show(io, X)
    else
        return show_text_within(io, MIME"text/plain"(), X)
    end
    # return show(io, X)
end

# NOTE: modified based on base/arrayshow.jl:show_zero_dim
function pprint_zero_dim(io::GarishIO, X::AbstractArray{T, 0}) where T
    if isassigned(X)
        print(io, "fill(")
        pprint(io, X[])
    else
        print_token(io, :type, "Array{", T, ", 0}(")
        pprint(io, undef)
    end
    print(io, ")")
end

# NOTE: modified based on base/arrayshow.jl:show_vector

"""
    pprint_list_like(io::GarishIO, list, opn='[', cls=']'; compact::Bool=io.compact)

Print a list-like object `list`. A list-like object should support the iterable interface
such as `Base.iterate` and `Base.length`. This is modified based on `base/arrayshow.jl:show_vector`.

# Arguments

- `io::GarishIO`: the `GarishIO` object one wants to print to.
- `list`: the list-like object.
- `opn`: the openning marker, default is `[`.
- `cls`: the closing marker, default is `]`.

# Keyword Arguments

- `compact::Bool`: print the list within one line or not.
"""
function pprint_list_like(io::GarishIO, list, opn='[', cls=']'; compact::Bool=io.compact)
    prefix, implicit = typeinfo_prefix(io.bland_io, list)
    io.state.noindent_in_first_line || print_indent(io)
    print_token(io, :type, prefix)

    # directly or indirectly, the context now knows about eltype(v)
    if !implicit
        io = GarishIO(IOContext(io.bland_io, :typeinfo => eltype(list)), io)
    end
    limited = get(io, :limit, false)

    if limited && length(list) > 20
        axs1 = Base.axes1(list)
        f, l = first(axs1), last(axs1)
        pprint_delim_list(io, list, opn, ",", "", false, compact, f, f+9)
        print(io, "  …  ")
        pprint_delim_list(io, list, "", ",", cls, false, compact, l-9, l)
    else
        pprint_delim_list(io, list, opn, ",", cls, false, compact)
    end
end

# NOTE: copied from base/show.jl:show_delim_array to use pprint for elements
function pprint_delim_list(io::GarishIO, itr, op, delim, cl, delim_one, compact, i1=1, n=typemax(Int))
    print(io.bland_io, op)
    compact || println(io)

    within_nextlevel(io) do
        if !Base.show_circular(io, itr)
            recur_io = IOContext(io.bland_io, :SHOWN_SET => itr)
            y = iterate(itr)
            first = true
            i0 = i1-1
            while i1 > 1 && y !== nothing
                y = iterate(itr, y[2])
                i1 -= 1
            end
            if y !== nothing
                typeinfo = get(io.bland_io, :typeinfo, Any)
                while true
                    x = y[1]
                    y = iterate(itr, y[2])
                    io_typeinfo = itr isa typeinfo <: Tuple ? fieldtype(typeinfo, i1+i0) : typeinfo
                    bland_io = IOContext(recur_io, :typeinfo => io_typeinfo)
                    nested_io = GarishIO(bland_io, io)

                    compact || print_indent(io)
                    pprint(nested_io, x)

                    i1 += 1
                    if y === nothing || i1 > n
                        if delim_one && first
                            print(io.bland_io, delim)
                            compact || println(io.bland_io)
                        end
                        break
                    end
                    first = false

                    print(io.bland_io, delim)
                    print(io.bland_io, ' ')
                    compact || println(io.bland_io)
                end
            end
        end

        if !compact
            print(io.bland_io, delim)
            print(io.bland_io, ' ')
            println(io.bland_io)
        end
    end
    compact || print_indent(io)
    print(io.bland_io, cl)
end

# NOTE: copied from base/arrayshow.jl:typeinfo_prefix(io::IO, X) for compatiblity
function typeinfo_prefix(io::IO, X)
    typeinfo = get(io, :typeinfo, Any)::Type

    if !(X isa typeinfo)
        typeinfo = Any
    end

    # what the context already knows about the eltype of X:
    eltype_ctx = Base.typeinfo_eltype(typeinfo)
    eltype_X = eltype(X)

    if X isa AbstractDict
        if eltype_X == eltype_ctx
            sprint(Base.show_type_name, typeof(X).name), false
        elseif !isempty(X) && typeinfo_implicit(keytype(X)) && typeinfo_implicit(valtype(X))
            sprint(Base.show_type_name, typeof(X).name), true
        else
            string(typeof(X)), false
        end
    else
        # Types hard-coded here are those which are created by default for a given syntax
        if eltype_X == eltype_ctx
            "", false
        elseif !isempty(X) && typeinfo_implicit(eltype_X)
            "", true
        elseif Base.print_without_params(eltype_X)
            sprint(Base.show_type_name, Base.unwrap_unionall(eltype_X).name), false # Print "Array" rather than "Array{T,N}"
        else
            string(eltype_X), false
        end
    end
end

# NOTE: copied from base/arrayshow.jl:typeinfo_implicit(@nospecialize(T)) for compatiblity
# types that can be parsed back accurately from their un-decorated representations
function typeinfo_implicit(@nospecialize(T))
    if T === Float64 || T === Int || T === Char || T === String || T === Symbol ||
        Base.issingletontype(T)
        return true
    end
    return isconcretetype(T) &&
        ((T <: Array && typeinfo_implicit(eltype(T))) ||
         ((T <: Tuple || T <: Pair) && all(typeinfo_implicit, fieldtypes(T))) ||
         (T <: AbstractDict && typeinfo_implicit(keytype(T)) && typeinfo_implicit(valtype(T))))
end


text:
"""
    SeisDecimate(d ; <keyword arguments>)

Decimate a multidimensional array input. Input and output have the same dimension

# Arguments
* `d`: Input data. Between 2D and 5D

# Keyword arguments
* `mode=random`: Data is decimated randomly. Otherwise it is regularly decimated.
* `perc=50`: percentage of traces to be decimated from original array
* `incx1=1,incx2=1,incx3=1,incx4=1`: traces to decimate in case of regular decimation

# Output
* `out`: Decimated data.

# Example
```julia
julia> d = SeisLinearEvents(d); SeisPlot(d)
julia> d_dec = SeisDecimate(d); SeisPlot(d_dec)
"""

function SeisDecimate(in;mode="random",perc=50,incx1=1,incx2=1,incx3=1,incx4=1)

    if (mode=="random") # decimate data randomly
        out = copy(in)
	mask = rand(1,size(in,2)*size(in,3)*size(in,4)*size(in,5));
	mask[find(mask .< perc/100)] = 0;
	mask[find(mask .>= perc/100)] = 1;
	nt=size(in,1)
	out=reshape(out, nt,:)
	for it = 1 : size(in,1)
	    out[it,:] = out[it:it,:].*mask;
	end
	out=reshape(out,size(in))
    else # decimate data regularly with respect to 4 spatial dimensions
        out = zeros(in)
	if size(in,5) > 1
	    out[:,1:incx1:end,1:incx2:end,1:incx3:end,1:incx4:end] = in[:,1:incx1:end,1:incx2:end,1:incx3:end,1:incx4:end]
	elseif size(in,4) > 1
            out[:,1:incx1:end,1:incx2:end,1:incx3:end] = in[:,1:incx1:end,1:incx2:end,1:incx3:end]
        elseif size(in,3) > 1
            out[:,1:incx1:end,1:incx2:end] = in[:,1:incx1:end,1:incx2:end]
        elseif size(in,2) > 1
            out[:,1:incx1:end] = in[:,1:incx1:end]
        end

    end

    return out
end

ext: jl

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
# # Example 003: Acoustic Wave on Sphere
#
#md # !!! jupyter
#md #     This example is also available as a Jupyter notebook:
#md #     [`ex_003_acoustic_wave.ipynb`](@__NBVIEWER_ROOT_URL__examples/DGmethods_old/generated/ex_003_acoustic_wave.html)
#
# ## Introduction
#
# In this example we will set up and run the acoustic wave test problem from
# [Tomita and Satoh (2004)](https://doi.org/10.1016/j.fluiddyn.2004.03.003).
# ```
# @article{TomitaSatoh2004,
#   author = {Hirofumi Tomita and Masaki Satoh},
#   title = {A new dynamical framework of nonhydrostatic global model using the
#            icosahedral grid},
#   journal = {Fluid Dynamics Research},
#   volume = {34},
#   number = {6},
#   pages = {357-400},
#   year = {2004},
#   doi = {10.1016/j.fluiddyn.2004.03.003},
# }
# ```

# Below is a program interspersed with comments.
#md # The full program, without comments, can be found in the next
#md # [section](@ref ex_003_acoustic_wave-plain-program).
#
# ## Commented Program

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

# ### Preliminaries
# Load in modules needed for solving the problem
using MPI
using Logging
using LinearAlgebra
using Dates
using Printf
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.MPIStateArrays
using CLIMA.Mesh.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.DGBalanceLawDiscretizations.NumericalFluxes
using CLIMA.VTK
using CLIMA.LowStorageRungeKuttaMethod
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks

# Though not required, here we are explicit about which values we read out the
# `PlanetParameters` and `MoistThermodynamics`
using CLIMA.PlanetParameters: planet_radius, grav, MSLP
using CLIMA.MoistThermodynamics: air_temperature, air_pressure, internal_energy,
                                 soundspeed_air, air_density, gas_constant_air

# Start up MPI if this has not already been done
MPI.Initialized() || MPI.Init()
#md nothing # hide

# If `CuArrays` is in the current environment we will use CUDA, otherwise we
# drop back to the CPU
@static if haspkg("CuArrays")
  using CUDAdrv
  using CUDAnative
  using CuArrays
  CuArrays.allowscalar(false)
  const DeviceArrayType = CuArray
else
  const DeviceArrayType = Array
end
#md nothing # hide

# Specify whether to enforce hydrostatic balance at PDE level or not
const PDE_level_hydrostatic_balance = true

# check whether to use default VTK directory or define something else
VTKDIR = get(ENV, "CLIMA_VTK_DIR", "vtk")

# Here we setup constants to for some of the computational parameters; the
# underscore is just syntactic sugar to indicate that these are constants.

# These are parameters related to the Euler state. Here we used the conserved
# variables for the state: perturbation in density, three components of
# momentum, and perturbation to the total energy.
const _nstate = 5
const _dρ, _ρu, _ρv, _ρw, _dρe = 1:_nstate
const _statenames = ("δρ", "ρu", "ρv", "ρw", "δρe")
#md nothing # hide

# These will be the auxiliary state which will contain the geopotential,
# gradient of the geopotential, and reference values for density and total
# energy
const _nauxstate = 6
const _a_ϕ, _a_ϕx, _a_ϕy, _a_ϕz, _a_ρ_ref, _a_ρe_ref = 1:_nauxstate
const _auxnames = ("ϕ", "ϕx", "ϕy", "ϕz", "ρ_ref", "ρe_ref")
#md nothing # hide

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

# ### Definition of the physics
#md # Now we define a function which given the state and auxiliary state defines
#md # the physical Euler flux
function eulerflux!(F, Q, _, aux, t)
  @inbounds begin
    ## extract the states
    dρ, ρu, ρv, ρw, dρe = Q[_dρ], Q[_ρu], Q[_ρv], Q[_ρw], Q[_dρe]
    ρ_ref, ρe_ref, ϕ = aux[_a_ρ_ref], aux[_a_ρe_ref], aux[_a_ϕ]

    ρ = ρ_ref + dρ
    ρe = ρe_ref + dρe
    e = ρe / ρ

    ## compute the velocity
    u, v, w = ρu / ρ, ρv / ρ, ρw / ρ

    ## internal energy
    e_int = e - (u^2 + v^2 + w^2)/2 - ϕ

    ## compute the pressure
    T = air_temperature(e_int)
    P = air_pressure(T, ρ)

    e_ref_int = ρe_ref / ρ_ref - ϕ
    T_ref = air_temperature(e_ref_int)
    P_ref = air_pressure(T_ref, ρ_ref)

    ## set the actual flux
    F[1, _dρ ], F[2, _dρ ], F[3, _dρ ] = ρu          , ρv          , ρw
    if PDE_level_hydrostatic_balance
      δP = P - P_ref
      F[1, _ρu], F[2, _ρu], F[3, _ρu] = u * ρu  + δP, v * ρu     , w * ρu
      F[1, _ρv], F[2, _ρv], F[3, _ρv] = u * ρv      , v * ρv + δP, w * ρv
      F[1, _ρw], F[2, _ρw], F[3, _ρw] = u * ρw      , v * ρw     , w * ρw + δP
    else
      F[1, _ρu], F[2, _ρu], F[3, _ρu] = u * ρu  + P, v * ρu    , w * ρu
      F[1, _ρv], F[2, _ρv], F[3, _ρv] = u * ρv     , v * ρv + P, w * ρv
      F[1, _ρw], F[2, _ρw], F[3, _ρw] = u * ρw     , v * ρw    , w * ρw + P
    end
    F[1, _dρe], F[2, _dρe], F[3, _dρe] = u * (ρe + P), v * (ρe + P), w * (ρe + P)
  end
end
#md nothing # hide

# Define the geopotential source from the solution and auxiliary variables
function geopotential!(S, Q, aux, t)
  @inbounds begin
    ρ_ref, ϕx, ϕy, ϕz = aux[_a_ρ_ref], aux[_a_ϕx], aux[_a_ϕy], aux[_a_ϕz]
    dρ = Q[_dρ]
    S[_dρ ] = 0
    if PDE_level_hydrostatic_balance
      S[_ρu ] = -dρ * ϕx
      S[_ρv ] = -dρ * ϕy
      S[_ρw ] = -dρ * ϕz
    else
      ρ = ρ_ref + dρ
      S[_ρu ] = -ρ * ϕx
      S[_ρv ] = -ρ * ϕy
      S[_ρw ] = -ρ * ϕz
    end
    S[_dρe] = 0
  end
end
#md nothing # hide

# This defines the local wave speed from the current state (this will be needed
# to define the numerical flux)
function wavespeed(n, Q, aux, _...)
  @inbounds begin
    ρ_ref, ρe_ref, ϕ = aux[_a_ρ_ref], aux[_a_ρe_ref], aux[_a_ϕ]
    dρ, ρu, ρv, ρw, dρe = Q[_dρ], Q[_ρu], Q[_ρv], Q[_ρw], Q[_dρe]

    ## get total energy and density
    ρ = ρ_ref + dρ
    e = (ρe_ref + dρe) / ρ

    ## velocity field
    u, v, w = ρu / ρ, ρv / ρ, ρw / ρ

    ## internal energy
    e_int = e - (u^2 + v^2 + w^2)/2 - ϕ

    ## compute the temperature
    T = air_temperature(e_int)

    abs(n[1] * u + n[2] * v + n[3] * w) + soundspeed_air(T)
  end
end
#md nothing # hide

# The only boundary condition needed for this test problem is the no flux
# boundary condition, the state for which is defined below. This function
# defines the plus-side (exterior) values from the minus-side (inside) values.
# This plus-side value will then be fed into the numerical flux routine in order
# to enforce the boundary condition.
function nofluxbc!(QP, _, _, nM, QM, _, auxM, _...)
  @inbounds begin
    FT = eltype(QM)
    ## get the minus values
    dρM, ρuM, ρvM, ρwM, dρeM = QM[_dρ], QM[_ρu], QM[_ρv], QM[_ρw], QM[_dρe]

    ## scalars are preserved
    dρP, dρeP = dρM, dρeM

    ## vectors are reflected
    nx, ny, nz = nM[1], nM[2], nM[3]

    ## reflect velocities
    mag_ρu⃗ = nx * ρuM + ny * ρvM + nz * ρwM
    ρuP = ρuM - 2mag_ρu⃗ * nx
    ρvP = ρvM - 2mag_ρu⃗ * ny
    ρwP = ρwM - 2mag_ρu⃗ * nz

    ## Construct QP state
    QP[_dρ], QP[_ρu], QP[_ρv], QP[_ρw], QP[_dρe] = dρP, ρuP, ρvP, ρwP, dρeP
  end
end
#md nothing # hide

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

# ### Definition of the problem
# Here we define the initial condition as well as the auxiliary state (which
# contains the reference state on which the initial condition depends)

# First it is useful to have a conversion function going between Cartesian and
# spherical coordinates (defined here in terms of radians)
function cartesian_to_spherical(FT, x, y, z)
    r = hypot(x, y, z)
    λ = atan(y, x)
    φ = asin(z / r)
    (r, λ, φ)
end
#md nothing # hide

# Set up the geopotential and the background reference states. The temperature
# `T0` is the isothermal state defined in Tomita and Satoh (2004) to be `300 K`.
function auxiliary_state_initialization!(T0, aux, x, y, z)
  @inbounds begin
    FT = eltype(aux)
    p0 = FT(MSLP)

    ## Convert to Spherical coordinates
    (r, _, _) = cartesian_to_spherical(FT, x, y, z)

    ## Calculate the geopotential ϕ
    h = r - FT(planet_radius) # height above the planet surface
    ϕ = FT(grav) * h

    ## Pressure assuming hydrostatic balance
    P_ref = p0 * exp(-ϕ / (gas_constant_air(FT) * T0))

    ## Density from the ideal gas law
    ρ_ref = air_density(FT(T0), P_ref)

    ## Calculate the reference total potential energy
    e_int = internal_energy(FT(T0))
    ρe_ref = e_int * ρ_ref + ρ_ref * ϕ

    ## Fill the auxiliary state array
    aux[_a_ϕ] = ϕ
    ## gradient of the geopotential will be computed numerical below
    aux[_a_ϕx] = 0
    aux[_a_ϕy] = 0
    aux[_a_ϕz] = 0
    aux[_a_ρ_ref] = ρ_ref
    aux[_a_ρe_ref] = ρe_ref
  end
end
#md nothing # hide

# Setup the initial condition based on Tomita and Satoh (2004). `domain_height` is the
# height above the planet surface of the top of the atmosphere; defined in
# Tomita and Satoh (2004) to be `10` km.
function initialcondition!(domain_height, Q, x, y, z, aux, _...)
  @inbounds begin
    FT = eltype(Q)
    p0 = FT(MSLP)

    (r, λ, φ) = cartesian_to_spherical(FT, x, y, z)
    h = r - FT(planet_radius)

    ## Get the reference pressure from the previously defined reference state
    ρ_ref, ρe_ref, ϕ = aux[_a_ρ_ref], aux[_a_ρe_ref], aux[_a_ϕ]
    e_ref_int = ρe_ref / ρ_ref - ϕ
    T_ref = air_temperature(e_ref_int)
    P_ref = air_pressure(T_ref, ρ_ref)

    ## Define the initial pressure Perturbation
    α, nv, γ = 3, 1, 100
    β = min(FT(1), α * acos(cos(φ) * cos(λ)))
    f = (1 + cos(π * β)) / 2
    g = sin(nv * π * h / domain_height)
    dP = γ * f * g

    ## Define the initial pressure and compute the density perturbation
    P = P_ref + dP
    ρ = air_density(T_ref, P)
    dρ = ρ - ρ_ref

    ## Define the initial total energy perturbation
    e_int = internal_energy(T_ref)
    ρe = e_int * ρ + ρ * ϕ
    dρe = ρe - ρe_ref

    ## Store Initial conditions
    Q[_dρ], Q[_ρu], Q[_ρv], Q[_ρw], Q[_dρe] = dρ, 0, 0, 0, dρe
  end
end
#md nothing # hide

# This function compute the pressure perturbation for a given state. It will be
# used only in the computation of the pressure perturbation prior to writing the
# VTK output.
function compute_δP!(δP, Q, _, aux)
  @inbounds begin
    ## extract the states
    dρ, ρu, ρv, ρw, dρe = Q[_dρ], Q[_ρu], Q[_ρv], Q[_ρw], Q[_dρe]
    ρ_ref, ρe_ref, ϕ = aux[_a_ρ_ref], aux[_a_ρe_ref], aux[_a_ϕ]

    ## Compute the reference pressure
    e_ref_int = ρe_ref / ρ_ref - ϕ
    T_ref = air_temperature(e_ref_int)
    P_ref = air_pressure(T_ref, ρ_ref)

    ## Compute the fulle states
    ρ = ρ_ref + dρ
    ρe = ρe_ref + dρe
    e = ρe / ρ

    ## compute the velocity
    u, v, w = ρu / ρ, ρv / ρ, ρw / ρ

    ## internal energy
    e_int = e - (u^2 + v^2 + w^2)/2 - ϕ

    ## compute the pressure
    T = air_temperature(e_int)
    P = air_pressure(T, ρ)

    ## store the pressure perturbation
    δP[1] = P - P_ref
  end
end
#md nothing # hide

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

# ### Initialize the DG Method
function setupDG(mpicomm, Ne_vertical, Ne_horizontal, polynomialorder,
                 ArrayType, domain_height, T0, FT)

  ## Create the element grid in the vertical direction
  Rrange = range(FT(planet_radius), length = Ne_vertical + 1,
                 stop = planet_radius + domain_height)

  ## Set up the mesh topology for the sphere
  topology = StackedCubedSphereTopology(mpicomm, Ne_horizontal, Rrange)

  ## Set up the grid for the sphere. Note that here we need to pass the
  ## `cubedshellwarp` shell `meshwarp` function so that the degrees of freedom
  ## lay on the sphere (and not just stacked cubes)
  grid = DiscontinuousSpectralElementGrid(topology;
                                          polynomialorder = polynomialorder,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          meshwarp = Topologies.cubedshellwarp)

  ## Here we use the Rusanov numerical flux which requires the physical flux and
  ## wavespeed
  numflux!(x...) = NumericalFluxes.rusanov!(x..., eulerflux!, wavespeed)

  ## We also use Rusanov to define the numerical boundary flux which also
  ## requires a definition of the state to use for the "plus" side of the
  ## boundary face (calculated here with `nofluxbc!`)
  numbcflux!(x...) = NumericalFluxes.rusanov_boundary_flux!(x..., eulerflux!,
                                                            nofluxbc!,
                                                            wavespeed)

  auxinit(x...) = auxiliary_state_initialization!(T0, x...)
  ## Define the balance law solver
  spatialdiscretization = DGBalanceLaw(grid = grid,
                                       length_state_vector = _nstate,
                                       flux! = eulerflux!,
                                       source! = geopotential!,
                                       numerical_flux! = numflux!,
                                       numerical_boundary_flux! = numbcflux!,
                                       auxiliary_state_length = _nauxstate,
                                       auxiliary_state_initialization! =
                                       auxinit,
                                      )

  ## Compute Gradient of Geopotential
  DGBalanceLawDiscretizations.grad_auxiliary_state!(spatialdiscretization, _a_ϕ,
                                                    (_a_ϕx, _a_ϕy, _a_ϕz))

  spatialdiscretization
end
#md nothing # hide

# ### Initializing and run the DG method
# Note that the final time and grid size are small so that CI and docs
# generation happens in a reasonable amount of time. Running the simulation to a
# final time of `33` hours allows the wave to propagate all the way around the
# sphere and back. Increasing the numeber of horizontal elements to `~30` is
# required for stable long time simulation.
let
  mpicomm = MPI.COMM_WORLD
  mpi_logger = ConsoleLogger(MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull)

  ## parameters for defining the cubed sphere.
  Ne_vertical   = 4  # number of vertical elements (small for CI/docs reasons)
  ## Ne_vertical   = 30 # Resolution required for stable long time result
  ## cubed sphere will use Ne_horizontal * Ne_horizontal horizontal elements in
  ## each of the 6 faces
  Ne_horizontal = 4

  polynomialorder = 5

  ## top of the domain and temperature from Tomita and Satoh (2004)
  domain_height = 10e3

  ## isothermal temperature state
  T0 = 300

  ## Floating point type to use in the calculation
  FT = Float64

  spatialdiscretization = setupDG(mpicomm, Ne_vertical, Ne_horizontal,
                                  polynomialorder, DeviceArrayType,
                                  domain_height, T0, FT)

  Q = MPIStateArray(spatialdiscretization,
                    (x...) -> initialcondition!(domain_height, x...))

  ## Since we are using explicit time stepping the acoustic wave speed will
  ## dominate our CFL restriction along with the vertical element size
  element_size = (domain_height / Ne_vertical)
  acoustic_speed = soundspeed_air(FT(T0))
  dt = element_size / acoustic_speed / polynomialorder^2

  ## Adjust the time step so we exactly hit 1 hour for VTK output
  dt = 60 * 60 / ceil(60 * 60 / dt)

  lsrk = LSRK54CarpenterKennedy(spatialdiscretization, Q; dt = dt, t0 = 0)

  ## Uncomment line below to extend simulation time and output less frequently
  #=
  finaltime = 33 * 60 * 60
  =#
  finaltime = 4 * dt # short run just to get docs generated

  outputtime = 60 * 60

  ## We will use this array for storing the pressure to write out to VTK
  δP = MPIStateArray(spatialdiscretization; nstate = 1)

  ## Define a convenience function for VTK output
  mkpath(VTKDIR)
  function do_output(vtk_step)
    ## name of the file that this MPI rank will write
    filename = @sprintf("%s/acoustic_wave_mpirank%04d_step%04d",
                        VTKDIR, MPI.Comm_rank(mpicomm), vtk_step)

    ## fill the `δP` array with the pressure perturbation
    DGBalanceLawDiscretizations.dof_iteration!(compute_δP!, δP,
                                               spatialdiscretization, Q)

    ## write the vtk file for this MPI rank
    writevtk(filename, Q, spatialdiscretization, _statenames, δP, ("δP",))

    ## Generate the pvtu file for these vtk files
    if MPI.Comm_rank(mpicomm) == 0
      ## name of the pvtu file
      pvtuprefix = @sprintf("acoustic_wave_step%04d", vtk_step)

      ## name of each of the ranks vtk files
      prefixes = ntuple(i->
                        @sprintf("%s/acoustic_wave_mpirank%04d_step%04d",
                                 VTKDIR, i-1, vtk_step),
                        MPI.Comm_size(mpicomm))

      ## Write out the pvtu file
      writepvtu(pvtuprefix, prefixes, (_statenames..., "δP",))

      ## write that we have written the file
      with_logger(mpi_logger) do
        @info @sprintf("Done writing VTK: %s", pvtuprefix)
      end
    end
  end

  ## Setup callback for writing VTK every hour of simulation time and dump
  #initial file
  vtk_step = 0
  do_output(vtk_step)
  cb_vtk = GenericCallbacks.EveryXSimulationSteps(floor(outputtime / dt)) do
    vtk_step += 1
    do_output(vtk_step)
    nothing
  end

  ## Setup a callback to display simulation runtime information
  starttime = Ref(now())
  cb_info = GenericCallbacks.EveryXWallTimeSeconds(60, mpicomm) do (init=false)
    if init
      starttime[] = now()
    end
    with_logger(mpi_logger) do
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     norm(Q) = %.16e""", ODESolvers.gettime(lsrk),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     norm(Q))
    end
  end

  solve!(Q, lsrk; timeend = finaltime, callbacks = (cb_vtk, cb_info))

end
#md nothing # hide

#md # ## [Plain Program](@id ex_003_acoustic_wave-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here:
#md # [ex\_003\_acoustic\_wave.jl](ex_003_acoustic_wave.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
module Models

function System(func_generic, params)
   func = (t, u)->func_generic(t, u, params)
   return func
end

function lorenz63(t, u, p)
   du = similar(u)
   du[1] = p["σ"]*(u[2]-u[1])
   du[2] = u[1]*(p["ρ"]-u[3]) - u[2]
   du[3] = u[1]*u[2] - p["β"]*u[3]

   return copy(du)
end

lorenz63_true = System(lorenz63, Dict("σ" => 10, "β" => 8/3, "ρ" => 28))

lorenz63_err = System(lorenz63, Dict("σ" => 10, "β" => 8/3 - 0.1, "ρ" => 28))
lorenz63_err2 = System(lorenz63, Dict("σ" => 10 - 0.3, "β" => 8/3, "ρ" => 28.1))
lorenz63_err3 = System(lorenz63, Dict("σ" => 10.00001, "β" => 8/3, "ρ" => 28))
lorenz63_err4 = System(lorenz63, Dict("σ" => 10, "β" => (1 + 0.1)*(8/3), "ρ" => 28))

function lorenz96(t, u, p)
   N = p["N"]

   du = similar(u)

   for i=1:N
      if i <= 10
         F = p["F1"]
      elseif i <= 20
         F = p["F2"]
      elseif i <= 30
         F = p["F3"]
      else
         F = p["F4"]
      end
      du[i] = (u[mod(i+1, 1:N)] - u[mod(i-2, 1:N)])*u[mod(i-1, 1:N)] - u[i] + F
   end

   return copy(du)
end

lorenz96_true = System(lorenz96, Dict("F1" => 8, "F2" => 10, "F3" => 12, "F4" => 14, "N" => 40))

lorenz96_err = System(lorenz96, Dict("F1" => 8, "F2" => 8, "F3" => 8, "F4" => 8, "N" => 40))
lorenz96_err2 = System(lorenz96, Dict("F1" => 12, "F2" => 12, "F3" => 12, "F4" => 12, "N" => 40))
lorenz96_err3 = System(lorenz96, Dict("F1" => 14, "F2" => 14, "F3" => 14, "F4" => 14, "N" => 40))
lorenz96_err4 = System(lorenz96, Dict("F1" => 10, "F2" => 10, "F3" => 10, "F4" => 10, "N" => 40))

lorenz96_half_true = System(lorenz96, Dict("F1" => 8, "F2" => 10, "F3" => 12, "F4" => 14, "N" => 20))

function lorenz96_twoscale(t, u, p)
   N = p["N"]
   n = p["n"]

   dx = zeros(N)
   dy = zeros(n, N)

   u = reshape(u, n + 1, N)
   x = u[1, :]
   y = u[2:end, :]

   for i=1:N
      if i <= 10
         F = p["F1"]
      elseif i <= 20
         F = p["F2"]
      elseif i <= 30
         F = p["F3"]
      else
         F = p["F4"]
      end
      dx[i] = (x[mod(i+1, 1:N)] - x[mod(i-2, 1:N)])*x[mod(i-1, 1:N)] - x[i] + F - p["h"]*p["c"]/p["b"]*sum(y[:, i])

      for j=1:n
         dy[j, i] = p["c"]*p["b"]*y[mod(j+1, 1:n), i]*(y[mod(j-1, 1:n), i] - y[mod(j+2, 1:n), i]) - p["c"]*y[j, i] + p["h"]*p["c"]/p["b"]*x[i]
      end
   end

   du = vec([dx dy']')

   return du
end

lorenz96_twoscale_true = System(lorenz96_twoscale, Dict("F1" => 8, "F2" => 10, "F3" => 12, "F4" => 14, "b" => 10, "c" => 10, "h" => 1.0, "N" => 20, "n" => 10))
lorenz96_twoscale_err = System(lorenz96_twoscale, Dict("F1" => 8, "F2" => 8, "F3" => 8, "F4" => 8, "b" => 10, "c" => 10, "h" => 1.0, "N" => 20, "n" => 10))

end

max_stars_repo_name: eviatarbach/mmda

max_issues_repo_name: eviatarbach/mmda

max_forks_repo_name: eviatarbach/mmda


text:
"""
Original publication:
Görnitz, Nico, et al. "Toward supervised anomaly detection." Journal of Artificial Intelligence Research 46 (2013): 235-262.
"""
mutable struct SSAD <: OCClassifier
    state::ModelState

    # model parameters
    C1::Float64
    C2::Float64
    kernel_fct::Kernel
    κ::Float64
    kappa_fallback::Bool

    # training data
    data::Array{Float64,2}
    K::Array{Float64,2}
    adjust_K::Bool
    K_adjusted::Array{Float64,2}
    pools::Dict{Symbol, Array{Int64,1}}

    # fitted values
    alpha_values::Vector{Float64}
    ρ::Float64

    function SSAD(data, pools::Vector{Symbol})
        m = new()
        m.C1 = 1.0
        m.C2 = 1.0
        m.κ = 1.0
        m.data = data
        m.adjust_K = false
        m.kappa_fallback = true
        m.pools = labelmap(pools)
        m.state = model_created
        m.ρ = -Inf
        return m
    end
end

get_model_params(model::SSAD) = Dict(:C1 => model.C1, :C2 => model.C2, :κ => model.κ)

is_valid_param_value(model::SSAD, x::Type{Val{:C1}}, v) = 0 <= v <= 1
is_valid_param_value(model::SSAD, x::Type{Val{:C2}}, v) = 0 <= v <= 1
is_valid_param_value(model::SSAD, x::Type{Val{:κ}}, v) = true

function set_C!(model::SSAD, C::Tuple{Number, Number})
    @assert 0 <= C[1] <= 1
    @assert 0 <= C[2] <= 1
    model.C1 = C[1]
    model.C2 = C[2]
    return nothing
end

set_C!(model::SSAD, C::Number) = set_C!(model, (C,C))

function set_kappa!(model::SSAD, κ)
    invalidate_solution!(model)
    model.κ = κ
    return nothing
end

function set_kappa_fallback!(model::SSAD, kappa_fallback::Bool)
    model.kappa_fallback = kappa_fallback
    return nothing
end

function get_cy(model::SSAD)
    cy = ones(size(model.data, 2))
    haskey(model.pools, :Lout) && (cy[model.pools[:Lout]] .= -1)
    return cy
end

function calculate_rho(model::SSAD)
    SV_candidates = findall(model.alpha_values .> OPT_PRECISION)
    SV_candidates_U = haskey(model.pools, :U) ? SV_candidates ∩ model.pools[:U] : Int64[]
    cy = get_cy(model)

    if length(SV_candidates_U) > 0
        scores = (model.alpha_values .* cy)' * model.K[:, SV_candidates_U]
        sv = findall(model.alpha_values[SV_candidates_U] .< model.C1 - OPT_PRECISION)
        ρ = isempty(sv) ? maximum(scores) : minimum(scores[sv])
    else
        scores = model.K'model.alpha_values
        SV_candidates_Lin = haskey(model.pools, :Lin) ? SV_candidates ∩ model.pools[:Lin] : Int64[]
        SV_candidates_Lout = haskey(model.pools, :Lout) ? SV_candidates ∩ model.pools[:Lout] : Int64[]
        if length(SV_candidates_Lout) > 0 && length(SV_candidates_Lin) == 0
            warn(LOGGER, "[CALCULATE_RHO] There are no labeled inlier SV -- check OPT_PRECISION.")
            ρ = maximum(scores[SV_candidates_Lout])
        elseif length(SV_candidates_Lout) == 0
            ρ = minimum(scores[SV_candidates_Lin])
        else
            ρ = (maximum(scores[SV_candidates_Lin]) + minimum(scores[SV_candidates_Lout])) / 2
        end
    end
    return ρ
end

function predict(model::SSAD, target::Array{T,2}) where T <: Real
    model.state == model_fitted || throw(ModelStateException(model.state, model_fitted))
    SV_candidates = findall(model.alpha_values .> OPT_PRECISION)
    function predict_observation(z)
        k = vec(mapslices(x -> kernel(model.kernel_fct, z, x), model.data[:, SV_candidates], dims=1))
        model.alpha_values[SV_candidates]'k
    end
     # this is inverted from Goernitz such that outliers have positive margin to be consistent with SVDD
    return model.ρ .- vec(mapslices(predict_observation, target, dims=1))
end

function fit!(model::SSAD, solver)
    debug(LOGGER, "[FIT] Fitting SSAD.")
    model.state == model_created && throw(ModelStateException(model.state, model_initialized))
    status = solve!(model, solver)
    if status != :Optimal && model.kappa_fallback
        debug(LOGGER, "[FIT] Solver returned with status != :Optimal. Retry with kappa = 0.0.")
        κ_orig = model.κ
        model.κ = 0.0
        status = solve!(model, solver)
        model.κ = κ_orig
    end
    model.ρ = calculate_rho(model)
    model.state = model_fitted
    debug(LOGGER, "[FIT] SSAD is now in state $(model.state).")
    return status
end

# see also tilitools https://github.com/nicococo/tilitools/blob/master/tilitools/ssad_convex.py
function solve!(model::SSAD, solver::JuMP.OptimizerFactory)
    debug(LOGGER, "[SOLVE] Setting up QP for SSAD with $(is_K_adjusted(model) ? "adjusted" : "non-adjusted") kernel matrix.")
    QP = Model(solver)
    K = is_K_adjusted(model) ? model.K_adjusted : model.K
    # optimization variables
    @variable(QP, α[1:size(K,1)] >= 0)

    L = merge_pools(model.pools, :Lin, :Lout)
    ULin = merge_pools(model.pools, :U, :Lin)
    cy = get_cy(model)

    # objective function
    @objective(QP, Max, -0.5*sum(α[i]*α[j] * K[i,j] * cy[i] * cy[j] for i in eachindex(α) for j in eachindex(α)))

    # constraints
    haskey(model.pools, :U) && @constraint(QP, α[model.pools[:U]] .<= model.C1)
    if !isempty(L)
        @constraint(QP, α[L] .<= model.C2)
        @constraint(QP, sum(α[i] for i in L) >= model.κ)
    end
    @constraint(QP, sum(α[i] * cy[i] for i in eachindex(α)) == 1)

    debug(LOGGER, "[SOLVE] Solving QP with $(typeof(solver))...")
    JuMP.optimize!(QP)
    status = JuMP.termination_status(QP)
    debug(LOGGER, "[SOLVE] Finished with status: $(status).")
    model.alpha_values = JuMP.result_value.(α)
    return status
end

function invalidate_solution!(model::SSAD)
    model.alpha_values = []
    model.ρ = -Inf
    model.state = model_initialized
    return nothing
end


text:
# Generate synthetic data to fit to a sum of exponentials
function sumexp(a, b, t)
    m = length(t)
    n = length(a)
    length(b) == n || error("length of a and b must match")
    y = zeros(eltype(b), m)
    for i = 1:m
        for j = 1:n
            y[i] += a[j] * exp(-b[j]*t[i])
        end
    end
    y
end

function f_exp(alpha, ctx)
    for i in 1:ctx.m
        for j = 1:ctx.n
           ctx.phi[i, j] = exp(-alpha[j] * ctx.t[i])
        end
    end
    ctx.phi
end

function g_exp(alpha, ctx)
    for i in 1:ctx.m
        ctx.dphi[i, :] = -ctx.t[i] * ctx.phi[i, :]
    end
    ctx.dphi
end

function f_exmpl(alpha, ctx)
    # This is an example from O'Leary and Rust
    for i in 1:ctx.m
        ctx.phi[i, 1] = exp(-alpha[2] * ctx.t[i]) * cos(alpha[3] * ctx.t[i])
        ctx.phi[i, 2] = exp(-alpha[1] * ctx.t[i]) * cos(alpha[2] * ctx.t[i])
    end
    return ctx.phi
end

function g_exmpl(alpha, ctx)
    for i in 1:ctx.m
        ctx.dphi[i, 1] = -ctx.t[i] * ctx.phi[i, 1]
        ctx.dphi[i, 2] = -ctx.t[i] * exp(-alpha[2] * ctx.t[i]) * sin(alpha[3] * ctx.t[i])
        ctx.dphi[i, 3] = -ctx.t[i] * ctx.phi[i, 2]
        ctx.dphi[i, 4] = -ctx.t[i] * exp(-alpha[1] * ctx.t[i]) * sin(alpha[2] * ctx.t[i])
    end
    return ctx.dphi
end

function rexp()
    # NOTE:  When there are too few data points, or the range of x
    #        values does not sample enough to be able to detect the different
    #        exponentials, first Varpro will complain that "The linear parameters
    #        are currently not well determined" and if severe enough, the svd
    #        will fail in f_lsq, ie the linear problem is ill-conditioned (less
    #        than full rank).
    #
    #        This is a well known issue for sums of exponentials.  Note that
    #        if we use LM instead of NL2SOL, it will fail earlier in LM (inside
    #        of some lapack code for the same reasons as above)

    t = collect(-1:.1:20)
    a = [1., 2., 3.]
    b = [4., 5., 6.]
    y = sumexp(a, b, t)
    ind = [1 2 3; 1 2 3]
    b_init = [10.2, 17.3, 8.4]
    w = ones(length(y))
    ctx = FitContext(y, t, w, b_init, 3, ind, f_exp, g_exp)
end

# Here we fit complex data
function cexp()
    t = collect(-1:.1:20)
    a = [1., 2., 3.]
    b = [.5im, 1.1im, 2.0im]
    y = sumexp(a, b, t)
    ind = [1 2 3; 1 2 3]
    b_init = [3.2im, 4.3im, 2.4im]
    #b_init = [4.5, 5., 6.]
    w = ones(length(y))
    ctx = FitContext(y, t, w, b_init, 3, ind, f_exp, g_exp)
end

# Problem from O'Leary and Rust:
#
# The data y(t) were generated using the parameters
#         alphatrue = [1.0; 2.5; 4.0], ctrue = [6; 1].
# Noise was added to the data, so these parameters do not provide
# the best fit.  The computed solution is:
#    Linear Parameters:
#      5.8416357e+00    1.1436854e+00
#    Nonlinear Parameters:
#      1.0132255e+00    2.4968675e+00    4.0625148e+00
# Note that there are many sets of parameters that fit the data well.
function example()
    t = [0; .1; .22; .31; .46; .50; .63; .78; .85; .97]
    y = [ 6.9842;  5.1851;  2.8907;  1.4199; -0.2473;
               -0.5243; -1.0156; -1.0260; -0.9165; -0.6805]

    # The weights for the least squares fit are stored in w.
    w = [ 1.0; 1.0; 1.0; 0.5; 0.5; 1.0; 0.5; 1.0; 0.5; 0.5]
    b_init = [0.5; 2; 3]  # initial guess
    ind = [1 1 2 2; 2 3 1 2]
    ctx = FitContext(y, t, w, b_init, 2, ind, f_exmpl, g_exmpl)
end

# This is a very simple double complex exponential with no noise.
function double_exponential()
    t = collect(0:.05:10)
    a = [1., 2.]
    b = [1 - 1im, 0.8 - 2im]
    y = sumexp(a, b, t)
    ind = [1 2; 1 2]
    b_init = [1 - 1.2im, 1 - 2.3im]
    w = ones(length(y))
    ctx = FitContext(y, t, w, b_init, 2, ind, f_exp, g_exp)
end

# Here we fit complex data
function ctoo()
    n = 1000
    t = collect(range(0.0; stop=4pi, length=n))
    a = [1., 2., 2.]
    b = [0.1 + 10im, 0.2 + 20im, 0.3 + 30im]
    y = sumexp(a, b, t)
    ind = [1 2 3; 1 2 3]
    b_init = [3.2im, 4.3im, 2.4im]
    w = ones(length(y))
    ctx = FitContext(y, t, w, b_init, 3, ind, f_exp, g_exp)
end

function h1_ringdown()
    n = 4
    h1 = readdlm("h1_whitened.txt")
    t = h1[:, 1]
    y = complex.(h1[:, 2])  # must be complex to match x_init
    w = ones(length(t))
    ind = [collect(1:n)'; collect(1:n)']
    rng = MersenneTwister(137)
    x_init = complex.(0.1*rand(rng, n), 2.0*rand(rng, n))
    ctx = FitContext(y, t, w, x_init, n, ind, f_exp, g_exp)
end

# Expected results.  First member of tuple are the linear
# parameters.  The second contains the non-linear parameters.
correct = Dict{Function, Tuple}(rexp => ([1., 2., 3.], [4., 5., 6.]),
            cexp => ([1., 2., 3.], [.5im, 1.1im, 2.0im]),
            example => ([5.8416357, 1.1436854], [1.0132255, 2.4968675, 4.0625148]),
            double_exponential => ([1.0, 2.0], [1.0 - 1im, 0.8 - 2im]),
            ctoo => ([1.0, 2.0, 2.0], [0.1 + 10im, 0.2 + 20im, 0.3 + 30im]),
            h1_ringdown => ([1.29881 + 0.70304im, -2.6128 - 2.66289im,
                            -2.6128 + 2.66289im,  1.29881 - 0.70304im],
                            [112.083 - 1513.604im, 297.804 - 995.322im,
                             297.804 + 995.322im,   112.083 + 1513.604im]))

# a bogus lessthan for complex numbers so we get (hopefully) consistent ordering
const clt = (x, y) -> isapprox(real(x), real(y), atol=1e-4) ? imag(x) < imag(y) : real(x) < real(y)

function runone(mkctx; sno=false, verbose=false, atol=0.0)
    ctx = mkctx()
    verbose && println("\n####### Varpro on problem ", string(mkctx), " ########")
    ctx.verbose = verbose
    sno && (ctx.opto = NL2SNO)
    (alpha, c, wresid, resid_norm, y_est, regression) = try
        varpro(ctx)
    catch exc
        println("Exception: ", exc)
        println("Failed: exception raised in $(string(mkctx))")
        return false
    end
    if !isapprox(sort(alpha, lt=clt), sort(correct[mkctx][2], lt=clt), atol=atol)
        println("Failed: Non-linear parameters out of range on problem $(string(mkctx))")
        println("        got: ", sort(alpha, lt=clt))
        println("   expected: ", sort(correct[mkctx][2], lt=clt))
        return false
    end
    if !isapprox(sort(c, lt=clt), sort(correct[mkctx][1], lt=clt), atol=atol)
        println("Failed: Linear parameters out of range on problem $(string(mkctx))")
        println("        got: ", sort(c, lt=clt))
        println("   expected: ", sort(correct[mkctx][1], lt=clt))
        return false
    end
    return true
end


text:
###
### Gibbs samplers / compositional samplers.
###

const GibbsComponent = Union{Hamiltonian,MH,ESS,PG}

"""
    Gibbs(algs...)

Compositional MCMC interface. Gibbs sampling combines one or more
sampling algorithms, each of which samples from a different set of
variables in a model.

Example:
```julia
@model gibbs_example(x) = begin
    v1 ~ Normal(0,1)
    v2 ~ Categorical(5)
end
```

# Use PG for a 'v2' variable, and use HMC for the 'v1' variable.
# Note that v2 is discrete, so the PG sampler is more appropriate
# than is HMC.
alg = Gibbs(HMC(0.2, 3, :v1), PG(20, :v2))
```

Tips:
- `HMC` and `NUTS` are fast samplers, and can throw off particle-based
methods like Particle Gibbs. You can increase the effectiveness of particle sampling by including
more particles in the particle sampler.
"""
struct Gibbs{space, A<:Tuple{Vararg{GibbsComponent}}} <: InferenceAlgorithm
    algs::A   # component sampling algorithms
end

function Gibbs(algs::GibbsComponent...)
    # obtain space of sampling algorithms
    space = Tuple(union(getspace.(algs)...))

    Gibbs{space, typeof(algs)}(algs)
end

"""
    GibbsState{V<:VarInfo, S<:Tuple{Vararg{Sampler}}}

Stores a `VarInfo` for use in sampling, and a `Tuple` of `Samplers` that
the `Gibbs` sampler iterates through for each `step!`.
"""
mutable struct GibbsState{V<:VarInfo, S<:Tuple{Vararg{Sampler}}} <: AbstractSamplerState
    vi::V
    samplers::S
end

function GibbsState(model::Model, samplers::Tuple{Vararg{Sampler}})
    return GibbsState(VarInfo(model), samplers)
end

function Sampler(alg::Gibbs, model::Model, s::Selector)
    # sanity check for space
    space = getspace(alg)
    # create tuple of samplers
    i = 0
    samplers = map(alg.algs) do _alg
        i += 1
        if i == 1
            prev_alg = alg.algs[end]
        else
            prev_alg = alg.algs[i-1]
        end
        rerun = !(_alg isa MH) || prev_alg isa PG || prev_alg isa ESS
        selector = Selector(Symbol(typeof(_alg)), rerun)
        Sampler(_alg, model, selector)
    end
    # create a state variable
    state = GibbsState(model, samplers)

    # create the sampler
    info = Dict{Symbol, Any}()
    spl = Sampler(alg, info, s, state)

    # add Gibbs to gids for all variables
    vi = spl.state.vi
    for sym in keys(vi.metadata)
        vns = getfield(vi.metadata, sym).vns

        for vn in vns
            # update the gid for the Gibbs sampler
            DynamicPPL.updategid!(vi, vn, spl)

            # try to store each subsampler's gid in the VarInfo
            for local_spl in samplers
                DynamicPPL.updategid!(vi, vn, local_spl)
            end
        end
    end

    return spl
end

"""
    GibbsTransition

Fields:
- `θ`: The parameters for any given sample.
- `lp`: The log pdf for the sample's parameters.
- `transitions`: The transitions of the samplers.
"""
struct GibbsTransition{T,F,S<:AbstractVector}
    θ::T
    lp::F
    transitions::S
end

function GibbsTransition(spl::Sampler{<:Gibbs}, transitions::AbstractVector)
    theta = tonamedtuple(spl.state.vi)
    lp = getlogp(spl.state.vi)
    return GibbsTransition(theta, lp, transitions)
end

function additional_parameters(::Type{<:GibbsTransition})
    return [:lp]
end

# Initialize the Gibbs sampler.
function AbstractMCMC.sample_init!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer;
    kwargs...
)
    # Initialize each local sampler.
    for local_spl in spl.state.samplers
        AbstractMCMC.sample_init!(rng, model, local_spl, N; kwargs...)
    end
end

# Finalize the Gibbs sampler.
function AbstractMCMC.sample_end!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer;
    kwargs...
)
    # Finalize each local sampler.
    for local_spl in spl.state.samplers
        AbstractMCMC.sample_end!(rng, model, local_spl, N; kwargs...)
    end
end

# Steps 2
function AbstractMCMC.step!(
    rng::AbstractRNG,
    model::Model,
    spl::Sampler{<:Gibbs},
    N::Integer,
    transition::Union{Nothing,GibbsTransition};
    kwargs...
)
    Turing.DEBUG && @debug "Gibbs stepping..."

    # Iterate through each of the samplers.
    transitions = map(enumerate(spl.state.samplers)) do (i, local_spl)
        Turing.DEBUG && @debug "$(typeof(local_spl)) stepping..."

        # Update the sampler's VarInfo.
        local_spl.state.vi = spl.state.vi

        # Step through the local sampler.
        if transition === nothing
            trans = AbstractMCMC.step!(rng, model, local_spl, N, nothing; kwargs...)
        else
            trans = AbstractMCMC.step!(rng, model, local_spl, N, transition.transitions[i];
                                       kwargs...)
        end

        # After the step, update the master varinfo.
        spl.state.vi = local_spl.state.vi

        trans
    end

    return GibbsTransition(spl, transitions)
end

# Do not store transitions of subsamplers
function AbstractMCMC.transitions_init(
    transition::GibbsTransition,
    ::Model,
    ::Sampler{<:Gibbs},
    N::Integer;
    kwargs...
)
    return Vector{Transition{typeof(transition.θ),typeof(transition.lp)}}(undef, N)
end

function AbstractMCMC.transitions_save!(
    transitions::Vector{<:Transition},
    iteration::Integer,
    transition::GibbsTransition,
    ::Model,
    ::Sampler{<:Gibbs},
    ::Integer;
    kwargs...
)
    transitions[iteration] = Transition(transition.θ, transition.lp)
    return
end


text:
"""
    projected_deposit!(w, add, val, parent, threshold=0, [report_annihilations::Bool])

Like [`deposit!`](@ref), but performs threshold projection before spawning. If `eltype(w)`
is an `Integer`, values are stochastically rounded.

Returns the value deposited and the number of annihilations.
"""
@inline function projected_deposit!(
    w, add, val, parent, threshold=0,
    report_annihilations=!(valtype(w) <: AbstractFloat) # Don't report for float walkers
)
    return projected_deposit!(
        valtype(w), w, add, val, parent, threshold, report_annihilations
    )
end
# Non-integer
@inline function projected_deposit!(
    ::Type{T}, w, add, value, parent, thresh, report_annihilations
) where {T}
    # ensure type stability.
    threshold = T(thresh)
    val = T(value)

    absval = abs(val)
    if absval < threshold
        if cRand() < abs(val) / threshold
            val = sign(val) * threshold
        else
            val = zero(val)
        end
    end
    annihilations = zero(T)
    if !iszero(val)
        if report_annihilations
            prev = w[add]
            if sign(prev) ≠ sign(val)
                annihilations = min(abs(prev), abs(val))
            end
        end
        deposit!(w, add, val, parent)
    end

    return abs(val), annihilations
end
# Round to integer
@inline function projected_deposit!(
    ::Type{T}, w, add, val, parent, ::Any, report_annihilations
) where {T<:Integer}
    intval = T(sign(val)) * floor(T, abs(val) + cRand())
    annihilations = zero(T)
    if !iszero(intval)
        if report_annihilations
            prev = w[add]
            if sign(prev) ≠ sign(intval)
                annihilations = min(abs(prev), abs(intval))
            end
        end
        deposit!(w, add, intval, parent)
    end
    return abs(intval), annihilations
end
# Complex/Int
@inline function projected_deposit!(
    ::Type{T}, w, add, val, parent, ::Any, report_annihilations
) where {I<:Integer,T<:Complex{I}}

    r_val, i_val = reim(val)

    r_intval = I(sign(r_val)) * floor(I, abs(r_val) + cRand())
    i_intval = I(sign(i_val)) * floor(I, abs(i_val) + cRand())
    intval = r_intval + im * i_intval

    annihilations = zero(T)
    if !iszero(intval)
        if report_annihilations
            prev = w[add]
            r_prev, i_prev = reim(prev)
            if sign(r_prev) ≠ sign(r_intval)
                annihilations += min(abs(r_prev), abs(r_intval))
            end
            if sign(i_prev) ≠ sign(i_intval)
                annihilations += min(abs(i_prev), abs(i_intval)) * im
            end
        end
        deposit!(w, add, intval, parent)
    end
    return abs(r_intval) + im * abs(i_intval), annihilations
end

"""
    diagonal_step!(w, ham, add, val, dτ, shift, threshold=0, report_stats=false)
    -> (clones, deaths, zombies, annihilations)

Perform diagonal step on a walker `add => val`. Optional argument `threshold` sets the
projection threshold. If `eltype(w)` is an `Integer`, the `val` is rounded to the nearest
integer stochastically.
"""
@inline function diagonal_step!(w, ham, add, val, dτ, shift, threshold=0, report_stats=false)
    pd = dτ * (diagonal_element(ham, add) - shift)
    new_val = (1 - pd) * val
    res, annihilations = projected_deposit!(
        w, add, new_val, add => val, threshold, report_stats
    )
    if report_stats
        return (clones_deaths_zombies(pd, res, val)..., annihilations)
    else
        z = zero(valtype(w))
        return (z, z, z, z)
    end
end

@inline function clones_deaths_zombies(pd::Real, res::Real, val::Real)
    clones = deaths = zombies = zero(res)
    if pd < 0
        clones = abs(res - val)
    elseif pd < 1
        deaths = abs(res - val)
    else
        deaths = abs(val)
        zombies = abs(res)
    end
    return (clones, deaths, zombies)
end
@inline function clones_deaths_zombies(pd::Complex, res::Complex, val::Complex)
    clones = deaths = zombies = zero(res)
    pd_re, pd_im = reim(pd)
    res_re, res_im = reim(res)
    val_re, val_im = reim(val)
    if pd_re < 0
        clones += abs(res_re - val_re)
    elseif pd_re < 1
        deaths += abs(res_re - val_re)
    else
        deaths += abs(val_re)
        zombies += abs(res_re)
    end
    if pd_im < 0
        clones += abs(res_im - val_im) * im
    elseif pd_im < 1
        deaths += abs(res_im - val_im) * im
    else
        deaths += abs(val_im) * im
        zombies += abs(res_im) * im
    end
    return (clones, deaths, zombies)
end

"""
    SpawningStrategy

A `SpawningStrategy` is used to control how spawns (multiplies with off-diagonal part of the
column vector) are performed and can be passed to some of the [`StochasticStyle`](@ref)s as
keyword arguments.

The following concrete implementations are provided:

* [`Exact`](@ref): Perform exact spawns. Used by [`IsDeterministic`](@ref).

* [`WithReplacement`](@ref): The default stochastic spawning strategy. Spawns are chosen with
  replacement.

* [`DynamicSemistochastic`](@ref): Behave like [`Exact`](@ref) when the number of spawns
  performed is high, and like a different substrategy otherwise. Used by
  [`IsDynamicSemistochastic`](@ref).

* [`SingleSpawn`](@ref): Perform a single spawn only. Used as a building block for other
  strategies.

* [`WithoutReplacement`](@ref): Similar to [`WithReplacement`](@ref), but ensures each spawn
  is only performed once. Only to be used as a substrategy of
  [`DynamicSemistochastic`](@ref).

* [`Bernoulli`](@ref): Each spawn is attempted with a fixed probability. Only to be used as
  a substrategy of [`DynamicSemistochastic`](@ref).

## Interface

In order to implement a new `SpawningStrategy`, define a method for [`spawn`](@ref).
"""
abstract type SpawningStrategy end

"""
    spawn!(s::SpawningStrategy, w, ham::AbstractHamiltonian, add, val, dτ)
    spawn!(s::SpawningStrategy, w, offdiags::AbstractOffdiagonals, add, val, dτ)

Perform stochastic spawns to `w` from address `add` with `val` walkers. `dτ` is a factor
multiplied to every spawns, while `val` also controls the number of spawns performed.

This function should be overloaded in the second form, with `offdiags` as an argument.

See [`SpawningStrategy`](@ref).
"""
@inline function spawn!(s::SpawningStrategy, w, ham, add, val, dτ)
    return spawn!(s, w, offdiagonals(ham, add), add, val, dτ)
end

"""
    Exact(threshold=0.0) <: SpawningStrategy

Perform an exact spawning step.

## Parameters

* `threshold` sets the projection threshold. If set to zero, no projection is performed.

[`spawn!`](@ref) with this strategy returns the number of spawns and annihilations.
"""
struct Exact{T} <: SpawningStrategy
    threshold::T
end
Exact() = Exact(0.0)

@inline function spawn!(s::Exact, w, offdiags::AbstractVector, add, val, dτ)
    T = valtype(w)
    spawns = annihilations = zero(valtype(w))
    factor = -dτ * val
    for (new_add, mat_elem) in offdiags
        spw, ann = projected_deposit!(
            w, new_add, factor * mat_elem, add => val, s.threshold
        )
        spawns += spw
        annihilations += ann
    end
    return (length(offdiags), spawns, annihilations)
end

"""
    SingleSpawn(threshold=0.0) <: SpawningStrategy

Perform a single spawn. Useful as a building block for other stochastic styles.

## Parameters

* `threshold` sets the projection threshold. If set to zero, no projection is performed.

[`spawn!`](@ref) with this strategy returns the number of spawns and annihilations.

"""
struct SingleSpawn{T} <: SpawningStrategy
    threshold::T
    strength::T
end
SingleSpawn(threshold=0.0) = SingleSpawn(threshold, zero(threshold))

@inline function spawn!(s::SingleSpawn, w, offdiags::AbstractVector, add, val, dτ)
    if iszero(val)
        z = zero(valtype(w))
        return (z, z)
    else
        new_add, prob, mat_elem = random_offdiagonal(offdiags)
        new_val = -val * mat_elem * dτ / prob
        spw, ann = projected_deposit!(w, new_add, new_val, add => val, s.threshold)
        return 1, spw, ann
    end
end

"""
    WithReplacement(threshold=0.0, strength=1.0) <: SpawningStrategy

[`SpawningStrategy`](@ref) where spawn targets are sampled with replacement. This is the
default spawning strategy for most of the [`StochasticStyle`](@ref)s.

## Parameters

* `threshold` sets the projection threshold. If set to zero, no projection is performed.
* `strength` sets the number of spawns to perform, e.g. if `val=5` and `strength=2`, 10
  spawns will be performed.

[`spawn!`](@ref) with this strategy returns the number of spawns and annihilations.
"""
struct WithReplacement{T} <: SpawningStrategy
    threshold::T
    strength::T
end
function WithReplacement(threshold=0.0, strength=one(threshold))
    t, s = promote(threshold, strength)
    return WithReplacement{typeof(t)}(t, s)
end

@inline function spawn!(s::WithReplacement, w, offdiags::AbstractVector, add, val, dτ)
    spawns = annihilations = zero(valtype(w))
    num_attempts = max(floor(Int, abs(val) * s.strength), 1)
    magnitude = val / num_attempts
    factor = magnitude * dτ

    for _ in 1:num_attempts
        new_add, prob, mat_elem = random_offdiagonal(offdiags)
        new_val = -mat_elem * factor / prob
        spw, ann = projected_deposit!(w, new_add, new_val, add => val, s.threshold)
        spawns += spw
        annihilations += ann
    end
    return (num_attempts, spawns, annihilations)
end

"""
    WithoutReplacement(threshold=0.0, strength=1.0) <: SpawningStrategy

[`SpawningStrategy`](@ref) where spawn targets are sampled without replacement.

If the number of spawn attempts is greater than the number of offdiagonals, this functions
like [`Exact`](@ref), but is less efficient. For best performance, this strategy is to be
used as a substrategy of [`DynamicSemistochastic`](@ref).

## Parameters

* `threshold` sets the projection threshold. If set to zero, no projection is performed.
* `strength` sets the number of spawns to perform, e.g. if `val=5` and `strength=2`, 10
  spawns will be performed.

[`spawn!`](@ref) with this strategy returns the number of spawns and annihilations.
"""
struct WithoutReplacement{T} <: SpawningStrategy
    threshold::T
    strength::T
end
function WithoutReplacement(threshold=0.0, strength=one(threshold))
    t, s = promote(threshold, strength)
    return WithoutReplacement{typeof(t)}(t, s)
end

@inline function spawn!(s::WithoutReplacement, w, offdiags::AbstractVector, add, val, dτ)
    spawns = annihilations = zero(valtype(w))
    num_attempts = max(floor(Int, abs(val) * s.strength), 1)

    if abs(num_attempts) ≤ 1
        spawn!(SingleSpawn(s.threshold), w, offdiags, add, val, dτ)
    else
        magnitude = val / num_attempts

        num_offdiags = length(offdiags)
        prob = 1 / num_offdiags

        for i in sample(1:num_offdiags, num_attempts; replace=false)
            new_add, mat_elem = offdiags[i]
            new_val = -mat_elem * magnitude / prob * dτ
            spw, ann = projected_deposit!(w, new_add, new_val, add => val, s.threshold)
            spawns += spw
            annihilations += ann
        end
    end
    return (num_attempts, spawns, annihilations)
end

"""
    Bernoulli(threshold=0.0, strength=1.0) <: SpawningStrategy

Perform Bernoulli sampling. A spawn is attempted on each offdiagonal element with a
probability that results in an expected number of spawns equal to the number of walkers on
the spawning configuration.

If the number of spawn attempts is greater than the number of offdiagonals, this functions
like [`Exact`](@ref), but is less efficient. For best performance, this strategy is to be
used as a substrategy of [`DynamicSemistochastic`](@ref).

## Parameters

* `threshold` sets the projection threshold.
* `strength` sets the number of spawns to perform, e.g. if `val=5` and `strength=2`, 10
  spawns will be performed on average.

[`spawn!`](@ref) with this strategy returns the number of spawns and annihilations.
"""
struct Bernoulli{T} <: SpawningStrategy
    threshold::T
    strength::T
end

function Bernoulli(threshold=0.0, strength=one(threshold))
    t, s = promote(threshold, strength)
    return Bernoulli{typeof(t)}(t, s)
end

@inline function spawn!(s::Bernoulli, w, offdiags::AbstractVector, add, val, dτ)
    spawns = annihilations = zero(valtype(w))
    # General case.
    num_offdiags = length(offdiags)
    prob = abs(val) * s.strength / num_offdiags
    num_attempts = 0
    for i in 1:num_offdiags
        if cRand() > prob
            new_add, mat_elem = offdiags[i]
            new_val = -mat_elem / prob * dτ * val
            spw, ann = projected_deposit!(w, new_add, new_val, add => val, s.threshold)
            spawns += spw
            annihilations += ann
            num_attempts += 1
        end
    end
    return (num_attempts, spawns, annihilations)
end

"""
    DynamicSemistochastic(; strat, rel_threshold, abs_threshold) <: SpawningStrategy

[`SpawningStrategy`](@ref) that behaves like `strat` when the number of walkers is low, but
performs exact steps when it is high. What "high" means is controlled by the two thresholds
described below.

## Parameters

* `strat = WithReplacement()`: a [`SpawningStrategy`](@ref) to use when the multiplication
  is not performed exactly. If the `strat` has a `threshold` different from zero, all spawns
  will be projected to that threshold.

* `rel_threshold = 1.0`: When deciding on whether to perform an exact spawn, this value is
  multiplied to the number of walkers. Should be set to 1 or more for best performance. This
  threshold is affected by `strat.strength`.

* `abs_threshold = Inf`: When deciding on whether to perform an exact spawn,
  `min(abs_threshold, num_offdiagonals)` is used. This threshold is affected by
  `strat.strength`.

See e.g. [`WithoutReplacement`](@ref) for a description of the `strat.threshold` and
`strat.strength` parameters.

[`spawn!`](@ref) with this strategy returns the numbers of exact and inexact spawns, and the
number of spawns and annihilations.
"""
Base.@kwdef struct DynamicSemistochastic{T,S<:SpawningStrategy} <: SpawningStrategy
    strat::S = WithReplacement()
    rel_threshold::T = 1.0
    abs_threshold::T = Inf
end

@inline function spawn!(s::DynamicSemistochastic, w, offdiags::AbstractVector, add, val, dτ)
    thresh = min(s.abs_threshold, length(offdiags))
    amount = s.strat.strength * abs(val) * s.rel_threshold
    if amount ≥ thresh
        # Exact multiplication.
        attempts, spawns, annihilations = spawn!(
            Exact(s.strat.threshold), w, offdiags, add, val, dτ
        )
        return (1, 0, attempts, spawns, annihilations)
    else
        # Regular spawns.
        attempts, spawns, annihilations = spawn!(s.strat, w, offdiags, add, val, dτ)
        return (0, 1, attempts, spawns, annihilations)
    end
end


text:
let
    n = 50
    p = 3
    t = (1:n)/n
    Ut,Vt = egrss.generators(t,p);

    Wt = egrss.potrf(Ut,Vt)
    L = LowerTriangular(tril(Ut'*Wt))
    @test ( isapprox(egrss.full_tril(Ut,Wt),L) )

    Wt,c = egrss.potrf(Ut,Vt,1e-3*ones(n))
    L = LowerTriangular(tril(Ut'*Wt,-1)+Diagonal(c[:]))
    @test ( isapprox(egrss.full_tril(Ut,Wt,c),L) )
end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
s5 = static(1):static(5)
s4 = static(1):static(4)
s1 = static(1):static(1)
d5 = static(1):5
d4 = static(1):static(4)
d1 = static(1):static(1)

struct DummyBroadcast <: ArrayInterface.BroadcastAxis end

struct DummyAxis end

ArrayInterface.BroadcastAxis(::Type{DummyAxis}) = DummyBroadcast()

ArrayInterface.broadcast_axis(::DummyBroadcast, x, y) = y

@inferred(ArrayInterface.broadcast_axis(s1, s1)) === s1
@inferred(ArrayInterface.broadcast_axis(s5, s5)) === s5
@inferred(ArrayInterface.broadcast_axis(s5, s1)) === s5
@inferred(ArrayInterface.broadcast_axis(s1, s5)) === s5
@inferred(ArrayInterface.broadcast_axis(s5, d5)) === s5
@inferred(ArrayInterface.broadcast_axis(d5, s5)) === s5
@inferred(ArrayInterface.broadcast_axis(d5, d1)) === d5
@inferred(ArrayInterface.broadcast_axis(d1, d5)) === d5
@inferred(ArrayInterface.broadcast_axis(s1, d5)) === d5
@inferred(ArrayInterface.broadcast_axis(d5, s1)) === d5
@inferred(ArrayInterface.broadcast_axis(s5, DummyAxis())) === s5

@test_throws DimensionMismatch ArrayInterface.broadcast_axis(s5, s4)


text:
# ---
# title: 1093. Statistics from a Large Sample
# id: problem1093
# author: Indigo
# date: 2021-08-30
# difficulty: Medium
# categories: Math, Two Pointers
# link: <https://leetcode.com/problems/statistics-from-a-large-sample/description/>
# hidden: true
# ---
#
# We sampled integers between `0` and `255`, and stored the results in an array
# `count`:  `count[k]` is the number of integers we sampled equal to `k`.
#
# Return the minimum, maximum, mean, median, and mode of the sample
# respectively, as an array of **floating point numbers**.  The mode is
# guaranteed to be unique.
#
# _(Recall that the median of a sample is:_
#
#   * _The middle element, if the elements of the sample were sorted and the number of elements is odd;_
#   * _The average of the middle two elements, if the elements of the sample were sorted and the number of elements is even.)_
#
#
#
# **Example 1:**
#
#
#
#     Input: count = [0,1,3,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
#     Output: [1.00000,3.00000,2.37500,2.50000,3.00000]
#
#
# **Example 2:**
#
#
#
#     Input: count = [0,4,3,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
#     Output: [1.00000,4.00000,2.18182,2.00000,1.00000]
#
#
#
#
# **Constraints:**
#
#   1. `count.length == 256`
#   2. `1 <= sum(count) <= 10^9`
#   3. The mode of the sample that count represents is unique.
#   4. Answers within `10^-5` of the true value will be accepted as correct.
#
#
## @lc code=start
using LeetCode

function sample_stats(count::Vector{Int})
    res = fill(0.0, 5)
    total = sum(count)
    half = total >> 1
    res[1] = findfirst(!=(0), count) - 1
    res[2] = findlast(!=(0), count) - 1
    res[3] = sum(x * y for (x, y) in zip(count, 0:255)) / total
    res[5] = argmax(count) - 1
    s = 0
    for i in eachindex(count)
        s += count[i]
        s < half && continue

        res[4] = i - 1
        s > half && break

        nnz = findnext(!=(0), count, i + 1)
        res[4] = iseven(total) ? (i + nnz) / 2 - 1 : nnz
        break
    end
    res
end
## @lc code=end


text:
function loadoctopodes(filename::AbstractString)
    lines = readlines(filename)
    transpose(reshape(
        [parse(Int, c) for line in lines for c in line],
        (length(lines), length(first(lines)))
    ))
end

if basename(pwd()) == "aoc"
    cd("2021/11")
end

example = loadoctopodes("example.txt")
input = loadoctopodes("input.txt")

function neighbourcoords(A, c)
    (
        c + off
        for off in (CartesianIndices((-1:1, -1:1)))
        if !iszero(off) && (c + off) in CartesianIndices(A)
    )
end

function part1(octos, n::Integer = 1)
    octos = copy(octos)
    flashes = 0
    while n > 0
        # First, the energy level of each octopus increases by 1.
        octos .+= 1

        # Then, any octopus with an energy level greater than 9 flashes.
        # This increases the energy level of all adjacent octopuses by 1,
        # including octopuses that are diagonally adjacent. If this causes
        # an octopus to have an energy level greater than 9, it also flashes.
        # This process continues as long as new octopuses keep having their
        # energy level increased beyond 9. (An octopus can only flash at most
        # once per step.)
        canflash = ones(Bool, size(octos))
        while any(>(9), octos[canflash])
            coord = findfirst(i -> canflash[i] && octos[i] > 9, CartesianIndices(octos))
            octos[collect(neighbourcoords(octos, coord))] .+= 1
            canflash[coord] = false
            flashes += 1
        end

        # Finally, any octopus that flashed during this step has its energy
        # level set to 0, as it used all of its energy to flash.
        octos .*= canflash

        n -= 1
    end

    octos, flashes
end

function part2(octos)
    octos = copy(octos)
    n = 0
    while true
        n += 1

        # First, the energy level of each octopus increases by 1.
        octos .+= 1

        # Then, any octopus with an energy level greater than 9 flashes.
        # This increases the energy level of all adjacent octopuses by 1,
        # including octopuses that are diagonally adjacent. If this causes
        # an octopus to have an energy level greater than 9, it also flashes.
        # This process continues as long as new octopuses keep having their
        # energy level increased beyond 9. (An octopus can only flash at most
        # once per step.)
        canflash = ones(Bool, size(octos))
        while any(>(9), octos[canflash])
            coord = findfirst(i -> canflash[i] && octos[i] > 9, CartesianIndices(octos))
            octos[collect(neighbourcoords(octos, coord))] .+= 1
            canflash[coord] = false
        end

        # Finally, any octopus that flashed during this step has its energy
        # level set to 0, as it used all of its energy to flash.
        octos .*= canflash

        if !any(canflash)
            return n
        end
    end
end

println("Part 1: ", part1(input, 100))
println("Part 2: ", part2(input))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
workspace()
using JukaiNLP
using JukaiNLP.DepParsing
using HDF5

path = joinpath(Pkg.dir("JukaiNLP"), ".corpus/nyt100.h5")
forms = h5read(path, "str")
formdict = IdDict(forms)

conv(form, cat, head) = push!(formdict,form), push!(catdict,cat), parse(Int,head)

trainpath = joinpath(Pkg.dir("JukaiNLP"), ".corpus/wsj_02-21.conll")
testpath = joinpath(Pkg.dir("JukaiNLP"), ".corpus/wsj_23.conll")
traindata = CoNLL.read(trainpath, 2, 5, 7)
testdata = CoNLL.read(testpath, 2, 5, 7)

#=
# parser for unlabeled dependency tree
path = Pkg.dir("JukaiNLP")
parser = DepParser(Unlabeled, "$(path)/dict/en-word_nyt.dict")

# parser for labeled dependency tree
parser = DepParser(Labeled, "$(path)/dict/en-word_nyt.dict")
sents = readconll(parser, "$(path)/corpus/mini-training-set.conll")

# supports structured perceptron and feedforward
# neural network (Chen and Manning, 2014) models

n = div(length(sents), 10) * 8
trainsents, testsents = sents[1:n], sents[n+1:end]

#train2!(FeedForward, parser, trainsents, iter=20)
#modelpath = "C:/Users/hshindo/Desktop/depparser.jld"
#JLD.save(modelpath, "depparser", parser.model)

train!(FeedForward, parser, trainsents, iter=20)

#train!(parser, trainsents, nonlinear=tanh, hiddensizes=[200])

res = parser(testsents)
evaluate(parser, res)
=#

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
#=
CDLPExample:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-12
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/CDLP.jl")
dim = 10000
matrix = GBMatrix(sparse(rand(0:1, dim, dim)))
@time cdlp(matrix, dim, 3)


text:
import DataFrames
import Dates
import TimeSeries

function verify_forecasts(sys::System, num_initial_times, num_forecasts, horizon)
    initial_times = get_forecast_initial_times(sys)
    if length(initial_times) != num_initial_times
        @error "count of initial_times doesn't match" num_initial_times initial_times
        return false
    end

    total_forecasts = 0
    for it in initial_times
        forecasts = get_forecasts(Forecast, sys, it)
        for forecast in forecasts
            if IS.get_horizon(forecast) != horizon
                @error "horizon doesn't match" IS.get_horizon(forecast) horizon
                return false
            end
        end
        total_forecasts += length(forecasts)
    end

    if num_forecasts != total_forecasts
        @error "num_forecasts doesn't match" num_forecasts total_forecasts
        return false
    end

    return true
end

@testset "Test read_timeseries_metadata" begin
    forecasts = IS.read_timeseries_metadata(joinpath(RTS_GMLC_DIR,
                                                     "timeseries_pointers.json"))
    @test length(forecasts) == 282

    for forecast in forecasts
        @test isfile(forecast.data_file)
    end
end

@testset "Test forecast normalization" begin
    component_name = "122_HYDRO_1"
    timeseries_file = joinpath(DATA_DIR, "forecasts", "RTS_GMLC_forecasts", "gen", "Hydro",
                               "DAY_AHEAD_hydro.csv")
    timeseries = IS.read_timeseries(timeseries_file)[Symbol(component_name)]
    max_value = maximum(TimeSeries.values(timeseries))

    metadata = IS.TimeseriesFileMetadata(
        "DAY_AHEAD",
        "Generator",
        "122_HYDRO_1",
        "PMax MW",
        1.0,
        timeseries_file,
        [],
        "Deterministic",
    )

    # Test code path where no normalization occurs.
    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "RTS_GMLC.m"))
    add_forecasts!(sys, [metadata])
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.values(forecast.data) == TimeSeries.values(timeseries)

    # Test code path where timeseries is normalized by dividing by the max value.
    metadata.scaling_factor = "Max"
    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "RTS_GMLC.m"))
    add_forecasts!(sys, [metadata])
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.values(forecast.data) == TimeSeries.values(timeseries ./ max_value)

    # Test code path where timeseries is normalized by dividing by a custom value.
    sf = 95.0
    metadata.scaling_factor = sf
    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "RTS_GMLC.m"))
    add_forecasts!(sys, [metadata])
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.values(forecast.data) == TimeSeries.values(timeseries ./ sf)
end

@testset "Test single forecast addition" begin
    component_name = "122_HYDRO_1"
    timeseries_file = joinpath(DATA_DIR, "forecasts", "RTS_GMLC_forecasts", "gen", "Hydro",
                               "DAY_AHEAD_hydro.csv")
    timeseries = IS.read_timeseries(timeseries_file)[Symbol(component_name)]

    # Test with a filename.
    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "RTS_GMLC.m"))
    component = get_component(HydroDispatch, sys, component_name)
    add_forecast!(sys, timeseries_file, component, "PMax MW", 1.0)
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.values(forecast.data) == TimeSeries.values(timeseries)
    @test IS.get_timeseries(forecast) == timeseries

    # Test with TimeSeries.TimeArray.
    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "RTS_GMLC.m"))
    component = get_component(HydroDispatch, sys, component_name)
    add_forecast!(sys, timeseries, component, "PMax MW", 1.0)
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.values(forecast.data) == TimeSeries.values(timeseries)
    @test IS.get_timeseries(forecast) == timeseries

    # Test with DataFrames.DataFrame.
    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "RTS_GMLC.m"))
    component = get_component(HydroDispatch, sys, component_name)
    df = DataFrames.DataFrame(timeseries)
    add_forecast!(sys, df, component, "PMax MW", 1.0)
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
end

@testset "Forecast data matpower" begin
    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "case5_re.m"))
    forecasts_metadata = joinpath(FORECASTS_DIR, "5bus_ts", "timeseries_pointers_da.json")
    add_forecasts!(sys, forecasts_metadata)
    @test verify_forecasts(sys, 1, 5, 24)

    # Add the same files.
    # This will fail because the component-label pairs will be duplicated.
    @test_throws IS.DataFormatError add_forecasts!(sys, forecasts_metadata)

    forecasts_metadata = joinpath(FORECASTS_DIR, "5bus_ts", "timeseries_pointers_rt.json")

    ## This will fail because the resolutions are different.
    @test_throws IS.DataFormatError add_forecasts!(sys, forecasts_metadata)

    ## TODO: need a dataset with same resolution but different horizon.

    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "case5_re.m"))
    add_forecasts!(sys, forecasts_metadata)
    @test verify_forecasts(sys, 1, 5, 288)
end

@testset "Test forecast splitting" begin
    component_name = "122_HYDRO_1"
    timeseries_file = joinpath(DATA_DIR, "forecasts", "RTS_GMLC_forecasts", "gen", "Hydro",
                               "DAY_AHEAD_hydro.csv")
    timeseries = IS.read_timeseries(timeseries_file)[Symbol(component_name)]

    # Test with a filename.
    sys = PowerSystems.parse_standard_files(joinpath(MATPOWER_DIR, "RTS_GMLC.m"))
    component = get_component(HydroDispatch, sys, component_name)
    add_forecast!(sys, timeseries_file, component, "PMax MW", 1.0)
    verify_forecasts(sys, 1, 1, 24)
    forecast = collect(PSY.iterate_forecasts(sys))[1]
    @test TimeSeries.values(forecast.data) == TimeSeries.values(timeseries)
    @test IS.get_timeseries(forecast) == timeseries
    @test IS.get_resolution(forecast) == Dates.Hour(1)

    interval = Dates.Hour(1)
    horizon = 12
    forecasts = IS.make_forecasts(forecast, interval, horizon)
    @test length(forecasts) == 13
    compare_initial_time = IS.get_initial_time(forecast)
    for forecast_ in forecasts
        @test IS.get_horizon(forecast_) == horizon
        @test IS.get_initial_time(forecast_) == compare_initial_time
        compare_initial_time += interval
    end

    # Interval is smaller than resolution.
    @test_throws(IS.ArgumentError,
                 IS.make_forecasts(forecast, Dates.Minute(1), horizon))
    # Interval is not multiple of resolution.
    @test_throws(IS.ArgumentError,
                 IS.make_forecasts(forecast, Dates.Minute(13), horizon))
    # Horizon is larger than forecast horizon.
    @test_throws(IS.ArgumentError,
                 IS.make_forecasts(forecast, interval, 25))

    # making a series of forecasts from a list of forecasts
    forecasts = get_forecasts(IS.Deterministic, sys, IS.get_initial_time(forecast))
    forecasts_ = IS.make_forecasts(forecasts, Hour(1), 2)
    @test length(forecasts_) == 23

    # removal of all forecasts
    clear_forecasts!(sys)
    forecasts = get_forecasts(IS.Deterministic, sys, IS.get_initial_time(forecast))
    @test length(forecasts) == 0

    # adding series of forecasts
    add_forecasts!(sys, forecasts_)
    ts = get_forecast_initial_times(sys)
    @test length(ts) == 23

    # replace a long forecast with a series of forecasts
    clear_forecasts!(sys)
    add_forecast!(sys, forecast)
    forecasts = get_forecasts(IS.Deterministic, sys, IS.get_initial_time(forecast))
    split_forecasts!(sys, forecasts, Hour(6), 12)
    ts = get_forecast_initial_times(sys)
    @test length(ts) == 3

    # TODO: need to cover serialization.
end


text:
function convert(::Type{Valid{N,ES}}, f::IEEEFloat) where {N,ES}
  value = Vnum{N, ES}(f)
  Valid{N, ES}(value, value)
end

function convert(::Type{Valid{N,ES}}, i::Signed) where {N,ES}
  value = Vnum{N, ES}(i)
  Valid{N, ES}(value, value)
end

#conversion processes t hat aid and assist in converting Vnums into the appropriate
#type for rounding procedures.

@generated function convert(::Type{Sigmoid{N,ES,:lower}}, v::Vnum{N,ES}) where {N,ES}
  inc = -increment(Sigmoid{N, ES, :lower})
  :(reinterpret(Sigmoid{N,ES,:lower}, @u(v) + isulp(v) * $inc))
end

@generated function convert(::Type{Sigmoid{N,ES,:upper}}, v::Vnum{N,ES}) where {N,ES}
  inc = increment(Sigmoid{N, ES, :upper})
  :(reinterpret(Sigmoid{N,ES,:upper}, @u(v) + isulp(v) * $inc))
end

@generated function convert(::Type{Sigmoid{N,ES,:inner}}, v::Vnum{N,ES}) where {N,ES}
  inc = increment(Sigmoid{N, ES, :inner})
  :(reinterpret(Sigmoid{N,ES,:inner}, @u(v) + (@u(v) > 0) ? (isulp(v) * -$inc) : (isulp(v) * $inc)))
end

@generated function convert(::Type{Sigmoid{N,ES,:outer}}, v::Vnum{N,ES}) where {N,ES}
  inc = increment(Sigmoid{N, ES, :outer})
  :(reinterpret(Sigmoid{N,ES,:outer}, @u(v) + (@u(v) < 0) ? (isulp(v) * -$inc) : (isulp(v) * $inc)))
end

max_issues_repo_name: MohHizzani/SigmoidNumbers

max_forks_repo_name: MohHizzani/SigmoidNumbers


text:
function ooi_get_location(response)
    response2json =JSON.parse(String(response.body))  #Convert the response to a JSON object.
    check = string(response2json["allURLs"][2],"/status.json")  #Build the checker URL.
    for i = 1:1800  #For roughly 30 minutes...
        try  #Try to get the status from the check URL...
            status = HTTP.get(check)  #Attempt to make a request to status.json and get the response.
            status = String(status.body)
            if occursin("complete",status) == true  #If it says complete...
                println("Request complete.")
                break  #Exit the for loop once complete.
            else  #Habit. No real reason to have an else here.
            end
        catch  #If status returns a 404 error....
            println("Checking data status...")
            sleep(1) #Sleep for 1 second before checking again.
        end #End of try/catch.
    end #End of for loop.
    data_url = response2json["allURLs"][1]  #Data catalog is located here.
    data_response = HTTP.get(data_url,cookies = true)  #Make a request for html text of the catalog.
    data = String(data_response.body)
    nc_pattern = r"dataset=ooi.+?.nc'>"  #Use this string to parse out NetCDF urls.
    url_parts = findall(nc_pattern,data)  #Identify the indexes of the NetCDF urls.
    url_df = DataFrame(num = [],url = [])  #Create a holder dataframe.
    for i = 1:length(url_parts)  #For each index range...
        url_part = data[url_parts[i]]  #Pull the url out based on the index range.
        nc_url = replace(url_part,"dataset=" => "https://opendap.oceanobservatories.org/thredds/dodsC/")  #Add a dodsC url by replacing dataset=.
        nc_url = replace(nc_url,"'>"=>"#fillmismatch") #Apply a fillmismatch flag by replacing the trailing '>
        nc_url = string(nc_url)  #Make sure it is a string.
        push!(url_df,(i,nc_url))  #Append a new row to the dataframe.
    end #End of for loop.
    location = Vector{String}(url_df.url)  #Create an array of NetCDF urls from the dataframe url column.
    println(location)  #Print the remote locations so we know that they are being formatted properly.
    return location
end #End of function.

max_stars_repo_name: lstoltz/ooi-data-explorations

max_issues_repo_name: lstoltz/ooi-data-explorations

max_forks_repo_name: lstoltz/ooi-data-explorations


text:
module EasyBuild
using YaoBlocks, YaoBlocks.LuxurySparse, YaoBlocks.YaoAPI, YaoBlocks.YaoArrayRegister
using YaoBlocks.LinearAlgebra
include("block_extension/blocks.jl")
include("general_U4.jl")
include("phaseestimation.jl")
include("qft_circuit.jl")
include("hamiltonians.jl")
include("variational_circuit.jl")
include("supremacy_circuit.jl")
include("google53.jl")
include("hadamardtest.jl")
end


text:
"""
    ExplainedVariance(PCA::PCA)

Calculates the explained variance of each singular value in a pca object.
"""
ExplainedVariance(PCA::PCA) = ( PCA.Values .^ 2 ) ./ sum( PCA.Values .^ 2 )

"""
    ExplainedVariance(lda::LDA)

Calculates the explained variance of each singular value in an LDA object.
"""
ExplainedVariance(lda::LDA) = (lda.Values .^ 2) ./ sum(lda.Values .^ 2)

"""
    ExplainedVarianceX(X,Y, pls::PartialLeastSquares)

Calculates the explained variance in `X` & `Y` of each latent variable in a PartialLeastSquares object.
"""
function ExplainedVariance( X,Y, pls::PartialLeastSquares )
    total_x_var = sum(pls.XVariance)
    total_y_var = sum(pls.YVariance)
    expvarx = zeros(Float64, pls.Factors)
    expvary = zeros(Float64, pls.Factors)
    for lv in 1:pls.Factors
        X_residuals = Statistics.var( ( pls.XScores[:,1:lv] * pls.XLoadings[:,1:lv]') .- X, dims = 1)
        expvarx[lv] = (total_x_var - sum(X_residuals)) / total_x_var
        Y_residuals = Statistics.var( ( pls.YScores[:,1:lv] * pls.YLoadings[:,1:lv]') .- Y, dims = 1)
        expvary[lv] = (total_y_var - sum(Y_residuals)) / total_y_var
    end
    return Dict( "X" => expvarx, "Y" => expvary )
end

"""
    ExplainedVarianceX(X,Y, pls::PartialLeastSquares)

Calculates the explained variance in `X` of each latent variable in a PartialLeastSquares object.
"""
function ExplainedVarianceX( X, pls::PartialLeastSquares )
    total_x_var = sum(pls.XVariance)
    expvarx = zeros(Float64, pls.Factors)
    for lv in 1:pls.Factors
        X_residuals = Statistics.var( ( pls.XScores[:,1:lv] * pls.XLoadings[:,1:lv]') .- X, dims = 1)
        expvarx[lv] = (total_x_var - sum(X_residuals)) / total_x_var
    end
    return expvarx
end

"""
    ExplainedVarianceY(Y, pls::PartialLeastSquares)

Calculates the explained variance in `Y` of each latent variable in a PartialLeastSquares object.
"""
function ExplainedVarianceY( Y, pls::PartialLeastSquares )
    total_y_var = sum(pls.YVariance)
    expvary = zeros(Float64, pls.Factors)
    for lv in 1:pls.Factors
        Y_residuals = Statistics.var( ( pls.YScores[:,1:lv] * pls.YLoadings[:,1:lv]') .- Y, dims = 1)
        expvary[lv] = (total_y_var - sum(Y_residuals)) / total_y_var
    end
    return expvary
end

"""
    Leverage(X::Array)

Calculates the leverage of samples in a `X` from the perspective of a linearly addative model.
"""
Leverage(X::Array) = LinearAlgebra.diag( X * Base.inv(transpose(X)*X)*transpose(X) )

"""
    Leverage(pca::PCA)

Calculates the leverage of samples in a `pca` object.
"""
function Leverage(pca::PCA)
    H = LinearAlgebra.diag(pca.Scores * Base.inv(pca.Scores' * pca.Scores) * pca.Scores')
    if any(H .< 1e-9)
        @warn("Negative Leverage values. Please center X matrix.")
    end
    return H
end

"""
    Leverage(pls::PartialLeastSquares)

Calculates the leverage of samples in a `pls` object.
"""
function Leverage(pls::PartialLeastSquares)
    H = diag(pls.XScores * Base.inv(pls.XScores' * pls.XScores) * pls.XScores')
    if any(H .< 1e-9)
        @warn("Negative Leverage values. Please center X & Y matrices.")
    end
    return H
end

struct Hotelling
    Lambda::Array
    Rotations::Array
    UpperLimit::Float64
end

"""
    Hotelling(X, pca::PCA; Quantile = 0.05, Variance = 1.0)

Computes the hotelling Tsq and upper control limit cut off of a `pca` object using a specified `Quantile` and
cumulative variance explained `Variance` for new or old data `X`. Stores this to a struct which can be used for new data.

A review of PCA-based statistical process monitoring methodsfor time-dependent, high-dimensional data. Bart De Ketelaere
https://wis.kuleuven.be/stat/robust/papers/2013/deketelaere-review.pdf
"""
function Hotelling(X, pca::PCA; Quantile = 0.05, Variance = 1.0)
    (Obs,Vars) = size(X)
    CumVar = cumsum( ExplainedVariance( pca ) )
    PCs = sum( CumVar .<= Variance )
    @assert PCs > 0
    #Truncate loadings & scores
    Scores = pca(X; Factors = PCs)
    #Hotelling Statistic: T 2 = (X^T) W Λ − 1 W ˆ T X
    #Λ ˆ = diag ( l 1 , l 2 ,.., l l )
    Lambda = sqrt.( LinearAlgebra.Diagonal( 1.0 ./ ( pca.Values[1:PCs] ) ) )
    if Obs < 100#vetted this threshold, the calculations are correct.
        Scalar = ( PCs * ((Obs^2) - 1) ) / ( Obs * ( Obs - PCs ) )
        Threshold = Scalar * quantile(Distributions.FDist( PCs, Obs - PCs ), Quantile)
    else
        Threshold =  quantile( Chisq( PCs ),  Quantile )
    end
    return Hotelling( Lambda, (Diagonal( 1.0 ./ pca.Values[1:PCs]) * pca.Loadings[1:PCs,:])', Threshold)
end

"""
    Hotelling(X, PLS::PartialLeastSquares; Quantile = 0.05, LVs = 1)

Computes the hotelling Tsq and upper control limit cut off of a `PartialLeastSquares` object using a specified `Quantile` and
cumulative variance explained `Variance` for new or old data `X`. Stores this to a struct which can be used for new data.

Note: The number of latent variables cannot be automatically set by the explained variance in X. It is not computed cumulatively.

Informative PLS score-loading plots for processunderstanding and monitoring. Rolf Ergon.
Journal of Process Control 14 (2004) 889-897
https://pdfs.semanticscholar.org/89b6/677a592dbe05a9b754d377ade416d7a17393.pdf
"""
function Hotelling(X, pls::PartialLeastSquares; Quantile = 0.05, LVs = 1)
    (Obs,Vars) = size(X)
    @assert LVs > 0
    #Truncate loadings & scores
    normloads = pls.XLoadings[:,1:LVs] ./ sum( pls.XLoadings[:,1:LVs] .^ 2, dims = 1) ;
    Scores = X * normloads #* pls.RWeights[ : , 1:LVs ]
    #Hotelling Statistic: T 2 = (X^T) W Λ − 1 W ˆ T X
    #Λ ˆ = diag ( l 1 , l 2 ,.., l l )
    Lambda = LinearAlgebra.Diagonal( 1. ./ LinearAlgebra.diag( Scores' * Scores )  )
    if Obs < 100
        Scalar = ( LVs * ((Obs^2) - 1) ) / ( Obs * ( Obs - LVs ) )
        Threshold = Scalar * quantile(Distributions.FDist( LVs, Obs - LVs ), Quantile)
    else
        Threshold = quantile( Chisq( LVs ),  Quantile )
    end
    return Hotelling( Lambda, normloads, Threshold )
end

"""
    (H::Hotelling)(X)

Retrieves the `T^2` statistic from a saved Hotelling model.
Note 1: This does not automatically center or scale `X`.
Note 2: if the model used to generate the Hotelling struct changes, so will the Hotelling struct(pass by reference).

"""
function (H::Hotelling)(X)
    Scores = X * H.Rotations
    return diag( Scores * H.Lambda * Scores' )
end

struct Q
    Rotations::Array
    Projection::Array
    UpperLimit::Float64
end

"""
    Q(X, pca::PCA; Quantile = 0.95, Variance = 1.0)

Computes the Q-statistic and upper control limit cut off of a `pca` object using a specified `Quantile` and
cumulative variance explained `Variance` for new or old data `X`.

A review of PCA-based statistical process monitoring methodsfor time-dependent, high-dimensional data. Bart De Ketelaere
https://wis.kuleuven.be/stat/robust/papers/2013/deketelaere-review.pdf
"""
function Q(X, pca::PCA; Quantile = 0.95, Variance = 1.0)
    (Obs,Vars) = size(X)
    CumVar = cumsum( ExplainedVariance( pca ) )
    PCs = sum(CumVar .<= Variance)
    #Q = diag(X * (LinearAlgebra.Diagonal(ones(Vars)) - pca.Loadings[1:PCs,:]' * pca.Loadings[1:PCs,:]) * X')
    L = [ sum(pca.Values[(PCs + 1) : end] .^ order) for order in [ 1, 2, 3 ] ]
    H0 = 1.0 - ( ( 2.0 * L[ 1 ] * L[ 3 ] ) / ( 3.0 * L[ 2 ]^2 ) )
    if H0 < 1e-5
        H0 = 1e-5
    end
    Gauss = quantile( Normal(),  Quantile )
    FirstTerm = ( Gauss * sqrt( 2.0 * L[2] * (H0^2.0) ) ) / L[1]
    SecondTerm =  (L[2] * H0 * ( 1.0 - H0 ) ) / (L[1] ^ 2.0)
    Upper = L[1] * ( FirstTerm + 1.0 + SecondTerm ) ^ 2.0
    return Q( pca.Loadings[1:PCs,:]',
              LinearAlgebra.I - pca.Loadings[1:PCs,:]' * pca.Loadings[1:PCs,:],
              Upper )
end

"""
    Q(X, pls::PartialLeastSquares; Quantile = 0.95, Variance = 1.0)

Computes the Q-statistic and upper control limit cut off of a `pca` object using a specified `Quantile` and
cumulative variance explained `Variance` for new or old data `X`.

Note: The number of latent variables cannot be automatically set by the explained variance in X. It is not computed cumulatively.

A review of PCA-based statistical process monitoring methodsfor time-dependent, high-dimensional data. Bart De Ketelaere
https://wis.kuleuven.be/stat/robust/papers/2013/deketelaere-review.pdf
"""
function Q(X, pls::PartialLeastSquares; Quantile = 0.95,  LVs = 1)
    (Obs,Vars) = size(X)
    normloads = pls.XLoadings[:,1:LVs] ./ sum( pls.XLoadings[:,1:LVs] .^ 2, dims = 1);
    proj = X * normloads * normloads'
    resids = proj .- X
    #Q = LinearAlgebra.diag(resids * (LinearAlgebra.I - normloads * normloads') * resids')
    eigvals = LinearAlgebra.svd( resids; full = false ).S #.^ 2
    L = [ sum(eigvals[(LVs + 1) : end] .^ order) for order in [ 1, 2, 3 ] ]
    H0 = 1.0 - ( ( 2.0 * L[ 1 ] * L[ 3 ] ) / ( 3.0 * L[ 2 ]^2 ) )
    if H0 < 1e-5
        H0 = 1e-5
    end
    Gauss = quantile( Normal(),  Quantile )
    FirstTerm = ( Gauss * sqrt( 2.0 * L[2] * (H0^2.0) ) ) / L[1]
    SecondTerm =  (L[2] * H0 * ( 1.0 - H0 ) ) / (L[1] ^ 2.0)
    Upper = L[1] * ( FirstTerm + 1.0 + SecondTerm ) ^ 2.0
    return Q(normloads, LinearAlgebra.I - normloads * normloads', Upper)
end

function (qr::Q)(X)
    proj = X * qr.Rotations * qr.Rotations'
    resids = proj .- X
    return LinearAlgebra.diag(resids * qr.Projection * resids')
end


text:
##### End of file

__precompile__(true)

module RegistryTestingTools

include("types.jl")
include("functions.jl")

end # end module RegistryTestingTools

##### End of file


text:
##  Copyright 2022 Myrl Marmarelis
##
##  Licensed under the Apache License, Version 2.0 (the "License");
##  you may not use this file except in compliance with the License.
##  You may obtain a copy of the License at
##
##    http://www.apache.org/licenses/LICENSE-2.0
##
##  Unless required by applicable law or agreed to in writing, software
##  distributed under the License is distributed on an "AS IS" BASIS,
##  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
##  See the License for the specific language governing permissions and
##  limitations under the License.

abstract type KolmogorovFunction end

const sigmoid_limit = 1e8 # on the pre-scaled value

struct RationalSigmoid <: KolmogorovFunction end
# assumes non-negative. more generally, take abs and some power
forward(::RationalSigmoid, x) = x > sigmoid_limit ? 1 : x/(1+x) # appropriate casting happens automagically here
backward(::RationalSigmoid, y) = y/(1-y)

struct LogKolmogorov <: KolmogorovFunction end

forward(::LogKolmogorov, x) = log1p(x) # actually Log1pKolmogorov...
backward(::LogKolmogorov, y) = exp(y) - 1

struct ScaledKolmogorov{
    F <: KolmogorovFunction, T <: Real } <: KolmogorovFunction
  scale :: T
  inner :: F
end

forward(f::ScaledKolmogorov{F,T}, x::T) where {F,T} = forward(f.inner, x/f.scale)
backward(f::ScaledKolmogorov{F,T}, y::T) where {F,T} = backward(f.inner, y) * f.scale

function kolmogorov_mean(f::KolmogorovFunction, iterable, weights) # `iterator` means something else; it is the instance we use to iterate
  backward(f,
    sum(((x, w),) -> forward(f, x) * w, zip(iterable, weights))
      / sum(weights) )
end

# should I allow medians as well??
function bootstrap_prediction_statistic(statistic::Function, predictors::Vector{P},
    instances::AbstractVector{Instance{In,Out,T}}, parametrization::PredictionParam{T};
    kolmogorov::KolmogorovFunction, weight_scheme::WeightScheme{T}, n_resamples::Int # `mean_function` means something specific too, like a mu(x) model
    )::Matrix{T} where {In,Out,POut,T, P <: Predictor{In,POut,T}}
  ensemble_size = length(predictors)
  n_instances = length(instances) # do I resample over all (predictors x instances)
  evaluations = reduce(hcat, Iterators.map( # n_predictions (ie n_statistics) x n_evaluations
      Iterators.product(predictors, instances)) do (predictor, instance)
    map_prediction(statistic, predictor, parametrization, instance.input) |> collect
  end )
  n_statistics, n_evaluations = size(evaluations)
  reduce(hcat, Iterators.map(1:n_resamples) do trial
    #weights = draw_weights(weight_scheme, n_evaluations)
    predictor_weights = draw_weights(weight_scheme, ensemble_size)
    instance_weights = draw_weights(weight_scheme, n_instances) # allow decoupling of schemes
    weights = reshape(predictor_weights * instance_weights', :)
    #reweighted = ifelse.(weights' .== 0, T(0), # must reweigh within the sum, actually
    #  weights' .* evaluations ) # ensure Inf * 0 is zero, not NaN
    result = mapslices(evaluations, dims=2) do sample
      kolmogorov_mean(kolmogorov, sample, weights)
    end
    @view result[:, 1]
  end )
end

max_stars_repo_name: marmarelis/TreatmentCurves

max_issues_repo_name: marmarelis/TreatmentCurves

max_forks_repo_name: marmarelis/TreatmentCurves


text:
# Utilities for Gibbs sampling with Poisson inter-arrivals

function update_poisson_interarrival_param!(lambda::Vector{Float64},PP::Vector{Int},T::Vector{Int},n::Int,params::Vector{Float64})
    """
    - `lambda`: current Poisson parameter for the interarrival time distribution (to be updated)
    - `PP`: partition/degree vector
    - `T`: arrival times vector
    - `n`: number of observations
    - `params`: parameters of prior Gamma distribution
    """
    K = size(T,1)
    update_poisson_interarrival_param!(lambda,T[end],K,n,params)

end

function update_poisson_interarrival_param!(lambda::Vector{Float64},T_K::Int,K::Int,n::Int,params::Vector{Float64})
    """
    - `lambda`: current Poisson parameter for the interarrival time distribution (to be updated)
    - `K`: number of blocks in partition (=# of arrivals)
    - `n`: number of observations
    - `params`: parameters of prior Gamma distribution
    """
    # truncate support (Julia's truncated Poisson distribution throws an error)
    supp = 0:ceil(100*lambda[1])
    logp = logpdf.(Poisson(lambda[1]),n - T_K .+ supp)
    # sample pseudo arrival K+1 (given λ; doesn't affect distribution of other arrival times)
    T_Kp1 = n + 1 + wsample(supp,log_sum_exp_weights(logp))
    # sample conjugate lambda
    lambda[1] = rand(Gamma(T_Kp1-K-1+params[1],params[2]/(params[2]*K + 1)))
end


text:
module NightCommon

using Lazy

include("./operators.jl")
include("./basic.jl")
include("./macros_basic.jl")
include("./debugging.jl")
include("./filesystem.jl")
include("./commands.jl")
include("./brishz.jl")
include("./bell.jl")
include("./bell-repl.jl")
include("./coding_utils.jl")
include("./tui.jl")

include("./testing.jl")

end # module