text:
function _build_atsit5_caches(::Type{T}) where {T}

    cs = SVector{6, T}(0.161, 0.327, 0.9, 0.9800255409045097, 1.0, 1.0)

    as = SVector{21, T}(
        #=a21=# convert(T,0.161),
        #=a31=# convert(T,-0.008480655492356989),
        #=a32=# convert(T,0.335480655492357),
        #=a41=# convert(T,2.8971530571054935),
        #=a42=# convert(T,-6.359448489975075),
        #=a43=# convert(T,4.3622954328695815),
        #=a51=# convert(T,5.325864828439257),
        #=a52=# convert(T,-11.748883564062828),
        #=a53=# convert(T,7.4955393428898365),
        #=a54=# convert(T,-0.09249506636175525),
        #=a61=# convert(T,5.86145544294642),
        #=a62=# convert(T,-12.92096931784711),
        #=a63=# convert(T,8.159367898576159),
        #=a64=# convert(T,-0.071584973281401),
        #=a65=# convert(T,-0.028269050394068383),
        #=a71=# convert(T,0.09646076681806523),
        #=a72=# convert(T,0.01),
        #=a73=# convert(T,0.4798896504144996),
        #=a74=# convert(T,1.379008574103742),
        #=a75=# convert(T,-3.290069515436081),
        #=a76=# convert(T,2.324710524099774)
    )

    btildes = SVector{7,T}(
                convert(T,-0.00178001105222577714),
                convert(T,-0.0008164344596567469),
                convert(T,0.007880878010261995),
                convert(T,-0.1447110071732629),
                convert(T,0.5823571654525552),
                convert(T,-0.45808210592918697),
                convert(T,0.015151515151515152)
    )

    rs = SVector{22, T}(
        #=r11=# convert(T,1.0),
        #=r12=# convert(T,-2.763706197274826),
        #=r13=# convert(T,2.9132554618219126),
        #=r14=# convert(T,-1.0530884977290216),
        #=r22=# convert(T,0.13169999999999998),
        #=r23=# convert(T,-0.2234),
        #=r24=# convert(T,0.1017),
        #=r32=# convert(T,3.9302962368947516),
        #=r33=# convert(T,-5.941033872131505),
        #=r34=# convert(T,2.490627285651253),
        #=r42=# convert(T,-12.411077166933676),
        #=r43=# convert(T,30.33818863028232),
        #=r44=# convert(T,-16.548102889244902),
        #=r52=# convert(T,37.50931341651104),
        #=r53=# convert(T,-88.1789048947664),
        #=r54=# convert(T,47.37952196281928),
        #=r62=# convert(T,-27.896526289197286),
        #=r63=# convert(T,65.09189467479366),
        #=r64=# convert(T,-34.87065786149661),
        #=r72=# convert(T,1.5),
        #=r73=# convert(T,-4),
        #=r74=# convert(T,2.5),
    )

    return cs, as, btildes, rs
end


text:
__precompile__()

module FilePathsBase

using Dates
using LinearAlgebra
using Printf

import Base: ==
export
    # Types
    AbstractPath,
    Path,
    PosixPath,
    WindowsPath,
    Mode,
    Status,

    # Methods
    anchor,
    cwd,
    drive,
    home,
    parts,
    root,
    hasparent,
    parents,
    filename,
    extension,
    extensions,
    exists,
    isabs,
    mode,
    created,
    modified,
    relative,
    move,
    remove,
    tmpname,
    tmpdir,
    mktmp,
    mktmpdir,
    chown,
    executable,
    readable,
    writable,
    raw,

    # Macros
    @p_str,
    @__PATH__,
    @__FILEPATH__,

    # Constants
    READ,
    WRITE,
    EXEC

export isexecutable

const PATH_TYPES = DataType[]

function __init__()
    register(PosixPath)
    register(WindowsPath)
end

"""
    AbstractPath

Defines an abstract filesystem path. Subtypes of `AbstractPath` should implement the
following methods:

- `Base.print(io, p)` (default: call base julia's joinpath with drive and path parts)
- `FilePathsBase.parts(p)`
- `FilePathsBase.root(p)`
- `FilePathsBase.drive(p)`
- `FilePathsBase.ispathtype(::Type{MyPath}, x::AbstractString) = true`
"""
abstract type AbstractPath end

function register(T::Type{<:AbstractPath})
    # We add the type to the beginning of our PATH_TYPES,
    # so that they can take precedence over the Posix and
    # Windows paths.
    pushfirst!(PATH_TYPES, T)
end

#=
We only want to print the macro string syntax when compact is true and
we want print to just return the string (this allows `string` to work normally)
=#
Base.print(io::IO, path::AbstractPath) = print(io, joinpath(drive(path), parts(path)...))

function Base.show(io::IO, path::AbstractPath)
    get(io, :compact, false) ? print(io, path) : print(io, "p\"$path\"")
end

Base.parse(::Type{<:AbstractPath}, x::AbstractString) = Path(x)
Base.convert(::Type{AbstractPath}, x::AbstractString) = Path(x)
Base.convert(::Type{String}, x::AbstractPath) = string(x)
Base.promote_rule(::Type{String}, ::Type{<:AbstractPath}) = String
ispathtype(::Type{<:AbstractPath}, x::AbstractString) = false

include("constants.jl")
include("utils.jl")
include("libc.jl")
include("mode.jl")
include("status.jl")
include("posix.jl")
include("windows.jl")
include("path.jl")
include("deprecates.jl")

end # end of module

ext: jl


text:
Base.:(==)(x::nodeTypeDln, y::Real) = ==(Int(x), y)
Base.:(==)(x::Real, y::nodeTypeDln) = ==(x, Int(y))
Base.convert(::Type{nodeTypeDln}, x::Real) = nodeTypeDln(Int(x))
Base.zero(::Type{nodeTypeDln}) = nodeTypeDln(0)

Base.getindex(x::DislocationLoop, i) = i == 1 ? x : throw(BoundsError())
Base.eachindex(x::DislocationLoop) = 1

function Base.zero(::Type{DislocationLoop})
    return DislocationLoop(;
        loopType = loopDln(),
        numSides = convert(Int, 0),
        nodeSide = convert(Int, 0),
        numLoops = convert(Int, 0),
        segLen = convert(Float64, 0),
        slipSystemIdx = convert(Int, 0),
        slipSystem = SlipSystem(nothing, zeros(3), zeros(3)),
        label = zeros(nodeTypeDln, 0),
        buffer = convert(Float64, 0),
        range = zeros(3, 0),
        dist = Zeros(),
    )
end

function Base.zero(::Type{DislocationNetwork})
    return DislocationNetwork(;
        links = zeros(Int, 2, 0),
        slipPlane = zeros(3, 0),
        bVec = zeros(3, 0),
        coord = zeros(3, 0),
        label = zeros(nodeTypeDln, 0),
        nodeVel = zeros(3, 0),
        nodeForce = zeros(3, 0),
        numNode = [0],
        numSeg = [0],
        maxConnect = 0,
        segForce = zeros(3, 2, 0),
        linksConnect = zeros(Int, 2, 0),
        extSeg = zeros(Bool, 0),
        segIdx = zeros(Int, 0, 3),
    )
end
function Base.push!(network::DislocationNetwork, n)
    elemT = eltype(network.coord)
    network = DislocationNetwork(;
        links = hcat(network.links, zeros(Int, 2, n)),
        slipPlane = hcat(network.slipPlane, zeros(elemT, 3, n)),
        bVec = hcat(network.bVec, zeros(elemT, 3, n)),
        coord = hcat(network.coord, zeros(elemT, 3, n)),
        label = vcat(network.label, zeros(nodeTypeDln, n)),
        nodeVel = hcat(network.nodeVel, zeros(elemT, 3, n)),
        nodeForce = hcat(network.nodeForce, zeros(elemT, 3, n)),
        numNode = copy(network.numNode),
        numSeg = copy(network.numSeg),
        maxConnect = copy(network.maxConnect),
        connectivity = hcat(
            network.connectivity,
            zeros(Int, size(network.connectivity, 1), n),
        ),
        linksConnect = hcat(network.linksConnect, zeros(Int, 2, n)),
        extSeg = vcat(network.extSeg, zeros(Bool, n)),
        segIdx = vcat(network.segIdx, zeros(Int, n, 3)),
        segForce = cat(network.segForce, zeros(elemT, 3, 2, n), dims = 3),
    )

    return network
end
function Base.getindex(network::DislocationNetwork, i)
    return @views network.links[:, i],
    network.slipPlane[:, i],
    network.bVec[:, i],
    network.coord[:, i],
    network.label[i],
    network.nodeVel[:, i],
    network.nodeForce[:, i],
    network.connectivity[:, i],
    network.linksConnect[:, i],
    network.extSeg[i],
    network.segIdx[i, :],
    network.segForce[:, :, i]
end
function Base.iszero(network::DislocationNetwork)
    (network.numNode[1] == 0 || network.numSeg[1] == 0) && return true
    return false
end
function Base.getindex(slipSystem::SlipSystem, i)
    return @views slipSystem.slipPlane[:, i], slipSystem.bVec[:, i]
end

"""
```
loopDistribution(::Rand, n, args...; kw...)
```
Returns a `3 × n` matrix whose points follow a random uniform distribution. Used by [`DislocationNetwork`](@ref) when the `dist` parameter of [`DislocationLoop`](@ref) is `Rand()`.
"""
loopDistribution(::Rand, n, args...; kw...) = rand(3, n)
"""
```
loopDistribution(::Zeros, n, args...; kw...)
```
Returns a `3 × n` zeros matrix. Used by [`DislocationNetwork`](@ref) when the `dist` parameter of [`DislocationLoop`](@ref) is `Zeros()`.
"""
loopDistribution(::Zeros, n, args...; kw...) = zeros(3, n)
"""
```
loopDistribution(::Rand, n, args...; kw...)
```
Returns a `3 × n` matrix whose points follow a random normal distribution. Used by [`DislocationNetwork`](@ref) when the `dist` parameter of [`DislocationLoop`](@ref) is `Randn()`.
"""
loopDistribution(::Randn, n, args...; kw...) = randn(3, n)
"""
```
loopDistribution(::Rand, n, args...; kw...)
```
Returns a `3 × n` matrix whose points follow are regularly placed. Used by [`DislocationNetwork`](@ref) when the `dist` parameter of [`DislocationLoop`](@ref) is `Regular()`.

!!! note
    Not yet implemented.
"""
function loopDistribution(::Regular, n, args...; kw...)
    return error("loopDistribution: regular distribution yet not implemented")
end

"""
```
limits!(lims, segLen, range, buffer)
```
Compute the bounding limits within which a [`DislocationLoop`](@ref) will be distributed in a [`DislocationNetwork`](@ref).
"""
function limits!(lims, segLen, range, buffer)
    @inbounds for i in 1:size(lims, 2)
        @simd for j in 1:size(lims, 1)
            lims[j, i] = range[j, i] + buffer * segLen
        end
    end
    return lims
end

"""
```
translatePoints!(coord, lims, disp)
```
Translates coordinates using the limits and displacements calculated by [`limits!`](@ref) and [`loopDistribution`](@ref).
"""
function translatePoints!(coord, lims, disp)
    @inbounds for i in 1:size(coord, 2)
        @simd for j in 1:size(coord, 1)
            coord[j, i] += lims[j, 1] + (lims[j, 2] - lims[j, 1]) * disp[j]
        end
    end
    return coord
end

"""
```
makeSegment(::segEdge, slipPlane, bVec)
```
Return the edge segment `(slipPlane × bVec) / || slipPlane × bVec ||`. Used to create segments for [`DislocationNetwork`](@ref).
"""
function makeSegment(::segEdge, slipPlane, bVec)
    edge = slipPlane × bVec
    return edge / norm(edge)
end
"""
```
makeSegment(::segEdge, slipPlane, bVec)
```
Return the edge segment `slipPlane / || slipPlane ||`. Used to create segments for [`DislocationNetwork`](@ref).
"""
function makeSegment(::segEdgeN, slipPlane, bVec)
    return slipPlane / norm(slipPlane)
end
"""
```
makeSegment(::segScrew, slipPlane, bVec)
```
Return the screw segment `bVec / || bVec ||`. Used to create segments for [`DislocationNetwork`](@ref).
"""
function makeSegment(::segScrew, slipPlane, bVec)
    return bVec / norm(bVec)
end

"""
```
makeConnect(links, maxConnect)
makeConnect!(network::DislocationNetwork)
```
Creates `connectivity` and `linksConnect` matrices for [`DislocationNetwork`](@ref).

* `connectivity` contains the number of other other nodes each node is connected to, up to `maxConnect` other nodes. Every `(2i, j)` entry contains a node connected to node `j`. Every `(2i+1, j)` coordinate contains whether that node is the first or second node in the link.
* `linksConnect` relates connections enabled by a link. Analogous to the connectivity of a link.
"""
function makeConnect(links, maxConnect)
    lenLinks = size(links, 2)
    connectivity = zeros(Int, 1 + 2 * maxConnect, lenLinks)
    linksConnect = zeros(Int, 2, lenLinks)

    _makeConnect!(links, connectivity, linksConnect)

    return connectivity, linksConnect
end
function makeConnect!(network::DislocationNetwork)
    links = network.links
    connectivity = network.connectivity
    linksConnect = network.linksConnect
    connectivity .= 0
    linksConnect .= 0

    _makeConnect!(links, connectivity, linksConnect)

    return nothing
end
function _makeConnect!(links, connectivity, linksConnect)
    lenLinks = size(links, 2)
    # Indices of defined links, undefined links are always at the end so we only need to know the first undefined entry.
    idx = findfirst(x -> x == 0, @view links[1, :])
    isnothing(idx) ? idx = lenLinks : idx -= 1
    # Loop through indices.
    @inbounds @simd for i in 1:idx
        # Node 1, it is the row of the coord matrix.
        n1 = links[1, i]
        n2 = links[2, i]
        # Num of other nodes node n1 is connected to.
        connectivity[1, n1] += 1
        connectivity[1, n2] += 1
        # Next column for n1 in connectivity.
        tmp1 = 2 * connectivity[1, n1]
        tmp2 = 2 * connectivity[1, n2]
        # i = linkID. 1, 2 are the first and second nodes in link with linkID
        connectivity[tmp1:(tmp1 + 1), n1] .= (i, 1)
        connectivity[tmp2:(tmp2 + 1), n2] .= (i, 2)
        # Connectivity of the nodes in link i.
        linksConnect[1, i] = connectivity[1, n1]
        linksConnect[2, i] = connectivity[1, n2]
    end
    return nothing
end

# TODO #15
"""
```
getSegmentIdx(links, label)
```
Finds the indices of a link and corresponding nodes.

## Returns

`n × 3` matrix is of the form `[i, node1, node2]`.

* `i` can be used to find the Burgers vector, slip plane and segment forces of segment `i`, eg `bVec[:, i]`.
* `node1` and `node2` can be used to find the coordinate and velocity of the nodes, eg `l = coord[:, node2] - coord[:, node1]`.
"""
function getSegmentIdx(links, label)
    lenLinks = size(links, 2)
    extSeg = zeros(Bool, lenLinks)
    segIdx = zeros(Int, lenLinks, 3)  # Indexing matrix.

    _getSegmentIdx!(segIdx, extSeg, lenLinks, links, label)

    return segIdx, extSeg
end
"""
```
getSegmentIdx!(network::DislocationNetwork)
```
Mutates the `segIdx` matrix in `network`. Works the same way as [`getSegmentIdx`](@ref).
"""
function getSegmentIdx!(network::DislocationNetwork)
    links = network.links
    label = network.label
    extSeg = network.extSeg
    segIdx = network.segIdx

    segIdx .= 0
    extSeg .= false
    lenLinks = size(links, 2)

    _getSegmentIdx!(segIdx, extSeg, lenLinks, links, label)

    return nothing
end

function _getSegmentIdx!(segIdx, extSeg, lenLinks, links, label)
    idx = findfirst(x -> x == 0, @view links[1, :])
    isnothing(idx) ? idx = lenLinks : idx -= 1

    @inbounds for i in 1:idx
        n1 = links[1, i]
        n2 = links[2, i]
        segIdx[i, :] .= (i, n1, n2)
        (label[n1] == extDln || label[n2] == extDln) ? extSeg[i] = true : continue
    end

    return nothing
end

"""
```
checkNetwork(network::DislocationNetwork)
```
Checks the validity of the dislocation network. It ensures the following conditions are met by the member variables of `network`:
1. `connectivity` and `links` have the same number of non-zero entries;
1. all entries in `bVec` are non-zero;
1. only the trailing columns of `connectivity` are zeros;
1. consistency between `connectivity` and `links`;
1. `bVec` is conserved among connected nodes;
1. entries in `links` are unique;
1. consistency betwen `connectivity` and `linksConnect`

"""
function checkNetwork(network::DislocationNetwork)
    links = network.links
    label = network.label
    connectivity = network.connectivity
    linksConnect = network.linksConnect

    idx = findfirst(x -> x == 0, @view links[1, :])
    isnothing(idx) ? idx = size(links, 2) : idx -= 1

    # Max value of connectivity should be the last link.
    maximum(connectivity) == idx ? nothing :
    error(
        "Non-empty entries of connectivity should be the same as the non-empty entries of links.",
    )

    bVec = network.bVec
    elemT = eltype(network.bVec)
    bSum = zeros(3)
    @inbounds for i in 1:idx
        iLinkBuffer = zeros(Int, 0)
        col = connectivity[1, i]

        # Check for zero Burgers vectors.
        staticBVec = SVector{3, elemT}(bVec[1, i], bVec[2, i], bVec[3, i])
        norm(staticBVec) < eps(elemT) ? error("Burgers vector must be non-zero.
                                       norm(bVec) = $(norm(bVec[i,:]))") : nothing

        # Trailing columns must be zero.
        sum(connectivity[(2 * (col + 1)):end, i]) != 0 ?
        error(
            "Trailing columns of connectivity must be zero.
      connectivity[$(i), $(2 * (col + 1)):end] = $(connectivity[i, (2 * (col + 1)):end]))",
        ) : nothing

        # Burgers vector conservation.
        bSum .= 0
        # No repeated links.
        neighbours = -ones(Int, col)
        for j in 1:col
            j2 = 2 * j
            iLink = connectivity[j2, i]     # Link ID.
            colLink = connectivity[j2 + 1, i] # Link position in links.
            colOppLink = 3 - connectivity[2 * j + 1, i] # Opposite column in links
            staticBVec = SVector{3, elemT}(bVec[1, iLink], bVec[2, iLink], bVec[3, iLink])
            bSum .+= (3 - 2 * colLink) * staticBVec # Running total of burgers vectors.

            neighbours[j] = links[colOppLink, iLink] # Find neighbouring nodes.
            push!(iLinkBuffer, iLink) # Push to link buffer for error printing.

            # Check that node does not appear twice in connectivity.
            @simd for k in (j + 1):col
                connectivity[j2, i] == connectivity[2 * k, i] ?
                error(
                    "Node $(i) cannot appear twice in connectivity, $(connectivity[i, j2])",
                ) : nothing
            end

            # Check links and connectivity are consistent.
            testNode = links[connectivity[j2 + 1, i], connectivity[j2, i]]
            testNode != i ?
            error(
                "Connectivity and links are not consistent.
connectivity[:, $(i)] = $(connectivity[:, i])
links[:, $(connectivity[j2, i])] = [$(links[1, connectivity[j2, i]]), $(links[2, connectivity[j2, i]])]",
            ) : nothing
        end

        # Check for Burgers vector conservation.
        abs(bSum ⋅ bSum) > eps(elemT) ?
        error("Burgers vector is not conserved, bSum = $(bSum).") : nothing

        # Find unique neighbours.
        uniqueNeighbours = unique(neighbours)
        length(uniqueNeighbours) != length(neighbours) ?
        error(
            "There must be no duplicate links. Each entry in links must correspond to a unique pair of adjacent entries in neighbours.
links = $(iLinkBuffer[1:end - 1])
uniqueNeighbours = $(uniqueNeighbours)
neighbours = $(neighbours)",
        ) : nothing

        # Check connectivity and linksConnect are consistent.
        check = (
            connectivity[2 * linksConnect[1, i], links[1, i]],
            connectivity[2 * linksConnect[2, i], links[2, i]],
        )
        i != check[1] | i != check[2] ? error("Inconsistent link.
        links[:, $(i)] = $(links[:, i])
        linksConnect[:, $(i)] = $(linksConnect[:, i])
        connectivity[:, $(links[1, i])] = $(connectivity[:, links[1, i]])
        connectivity[:, $(links[2, i])] = $(connectivity[:, links[2, i]])") : nothing
    end

    return true
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
"
 Setup script for docker

 author: Atsushi Sakai
"

using Pkg
Pkg.add("SciPy")
using SciPy


text:
type PLL
    α # Bandwidth
    β # Phase adjustor
    ƒ # Output frequency
    ϕ # Output phase
end

function PLL( α = 0.002, β = sqrt( α ) )
    PLL( α, β, 0.0, 0.0 )
end

function exec( pll::PLL, x::Number )
    y  = exp( pll.ϕ * im )
    Δϕ = angle( x * conj(y) )
    pll.ƒ += pll.α * Δϕ
    pll.ϕ += pll.β * Δϕ
    pll.ϕ += pll.ƒ

    return y, Δϕ
end

function exec{T}( pll::PLL, x::AbstractVector{T} )
    y = Array( T, length(x) )
    e = Array( typeof(real(x[1])), length(x) )

    for i in 1:length( x )
        ( yi, ei ) = exec( pll, x[i] )
        y[i] = yi
        e[i] = ei
    end

    return y, e
end


text:
export SpaceState, UnRandomized, Randomized

abstract type SpaceState <: AbstractTag end
abstract type UnRandomized <: SpaceState end
abstract type Randomized <: SpaceState end

# Exceptions
export UnRandomizedError
import Base: show

"""
Some operations (like some samplers) only accepts randomized sample space

Throw this error, when the arugment does not meet this condition.
"""
struct UnRandomizedError <: Exception
end

function show(io::IO, e::UnRandomizedError)
    println(io, "Need to randomize first")
end


text:
"""
```
function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matrix{S};
             kwargs...) where {S<:AbstractFloat, U<:Number}
```

### Arguments:
- `loglikelihood::Function`: Log-likelihood function of model being estimated. Takes `parameters`
    and `data` as arguments.
- `parameters::ParameterVector{U}`: Model parameter vector, which stores parameter values,
    prior dists, and bounds.
- `data`: A matrix or dataframe containing the time series of the observables used in
    the calculation of the posterior/loglikelihood

### Keyword Arguments:
- `verbose::Symbol`: Desired frequency of function progress messages printed to standard out.
	- `:none`: No status updates will be reported.
	- `:low`: Status updates for SMC initialization and recursion will be included (default).
	- `:high`: Status updates for every iteration of SMC is output, which includes
    the mean and standard deviation of each parameter draw after each iteration,
    as well as calculated acceptance rate, ESS, and number of times resampled.
- `parallel::Bool = false`: Flag for running algorithm in parallel.
- `n_parts::Int = 5_000`: Number of particles.
- `n_blocks::Int = 1`: Number of parameter blocks in mutation step.
- `n_mh_steps::Int = 1`: Number of Metropolis Hastings steps to attempt during the mutation step.
- `λ::S = 2.1`: The 'bending coefficient' λ in Φ(n) = (n/N(Φ))^λ
- `n_Φ::Int = 300`: Number of stages in the tempering schedule.
- `resampling_method::Symbol`: Which resampling method to use.
    - `:systematic`: Will use systematic resampling (default).
    - `:multinomial`: Will use multinomial resampling.
    - `:polyalgo`: Samples using a polyalgorithm.
- `threshold_ratio::S = 0.5`: Threshold s.t. particles will be resampled when the population
    drops below threshold * N.
- `c::S = 0.5`: Initial scaling factor for covariance of the particles. Controls size of steps in mutation step.
    This value will be adaptively set afterward to reach an accept rate of `target` (see kwarg below).
- `α::S = 1.0`: The mixture proportion for the mutation step's proposal distribution. See `?mvnormal_mixture_draw` for details.
    Note that a value of 0.9 has commonly been used in applications to DSGE models (see citations below).
- `target::S = 0.25`: The target acceptance rate for new particles during mutation.
- `use_fixed_schedule::Bool = true`: Flag for whether or not to use a fixed tempering (ϕ) schedule.
- `tempering_target::S = 0.97`: Coefficient of the sample size metric to be targeted when solving
    for an endogenous ϕ.
- `old_data::Matrix{S} = []`: A matrix containing the time series of observables of previous data
    (with `data` being the new data) for the purposes of a time tempered estimation
    (that is, using the posterior draws from a previous estimation as the initial set
    of draws for an estimation with new data). This matrix is used to compute the
    log-likelihood according to the `old_loglikelihood` function when particles change during tempering.
- `old_cloud::Cloud = Cloud(0, 0)`: associated cloud borne of old data in previous SMC estimation.
    Running a bridge estimation requires `old_data` and `old_cloud`. If no `old_cloud`
    is provided, then we will attempt to load one using `loadpath`.
- `old_loglikelihood::Function = loglikelihood`: the old log-likelihood function when running
    a time tempered estimation. This function should be able to evaluate the log-likelihood
    of `old_data` given the current `parameters`, which may be nontrivial if the current
    `parameters` include additional parameters relative to the old estimation.
    By default, we assume the log-likelihood function has not changed
    and therefore coincides with the current one, so that the log-likelihood of `old_data`
    can be computed as `loglikelihood(parameters, old_data)`
- `old_vintage::String = ""`: String for vintage date of old data
- `smc_iteration::Int = 1`: The iteration index for the number of times SMC has been run on the
     same data vintage. Primarily for numerical accuracy/testing purposes.
- `run_test::Bool = false`: Flag for when testing accuracy of program
- `filestring_addl::Vector{String} = []`: Additional file string extension for loading old cloud.
- `save_intermediate::Bool = false`: Flag for whether one wants to save intermediate Cloud objects
- `intermediate_stage_increment::Int`: Save Clouds at every increment
   (1 = each stage, 10 = every 10th stage, etc.). Useful if you are using a cluster with time
    limits because if you hit the time limit, then you can just
    start from an intermediate stage rather than start over.
- `continue_intermediate::Bool = false`: Flag to indicate whether one is continuing SMC from an
    intermediate stage.
- `intermediate_stage_start::Int = 10`: Intermediate stage at which one wishes to begin the estimation.
- `tempered_update_prior_weight::Float64 = 0.0`: Weight placed on the current priors of parameters
    to construct a convex combination of draws from current priors and the previous estimation's
    cloud. The convex combination serves as the bridge distribution for a time tempered estimation.
- `run_csminwel::Bool = true`: Flag to run the csminwel algorithm to identify the true posterior mode
    (which may not exist) after completing an estimation. The mode identified by SMC is just
    the particle with the highest posterior value, but we do not check it is actually a mode (i.e.
    the Hessian is negative definite).
- `regime_switching::Bool = false`: Flag if there are regime-switching parameters. Otherwise, not all the values of the
    regimes will be used or saved.
- `toggle::Bool = true`: Flag for resetting the fields of parameter values to regime 1 anytime
    the loglikelihood is computed. The regime-switching version of SMC assumes at various points
    that this resetting occurs. If speed is important, then ensure that the fields of parameters
    take their regime 1 values at the end of the loglikelihood computation and set `toggle = false`.
- `debug_assertion::Bool = false`: if true, then when an assertion error is thrown during the estimation,
    output is created in a JLD2 file to help the user debug the problem.
- `log_prob_old_data::Float64 = 0.0`: Log MDD of old data for correct incremental weights when bridging

### Outputs

- `cloud`: The Cloud object containing all of the information about the
    parameter values from the sample, their respective log-likelihoods, the ESS
    schedule, tempering schedule etc., which is saved in the saveroot.

### Overview

Sequential Monte Carlo can be used in lieu of Random Walk Metropolis Hastings to
    generate parameter samples from high-dimensional parameter spaces using
    sequentially constructed proposal densities to be used in iterative importance
    sampling.

This implementation is based on Edward Herbst and Frank Schorfheide's 2014 paper
    'Sequential Monte Carlo Sampling for DSGE Models' and the code accompanying their
    book 'Bayesian Estimation of DSGE Models'. Cai et al. (2021)
    'Online Estimation of DSGE Models' extend the algorithm in Herbst and Schorfheide (2014)
    to include an adaptive schedule and generalized tempering, which are both
    features of this package.

SMC is broken up into three main steps:

- `Correction`: Reweight the particles from stage n-1 by defining incremental weights,
    which gradually "temper in" the loglikelihood function ``p(Y \vert \\theta)^(\\phi_n - \\phi_n-1)`` into the
    normalized particle weights.
- `Selection`: Resample the particles if the distribution of particles begins to
    degenerate, according to a tolerance level for the ESS.
- `Mutation`: Propagate particles ``{\\theta(i), W(n)}`` via a Metropolis
    Hastings algorithm (the number of steps are specified by `n_mh_steps`).
"""
function smc(loglikelihood::Function, parameters::ParameterVector{U}, data::Matrix{S};
             verbose::Symbol = :low,
             testing::Bool   = false,
             data_vintage::String = Dates.format(today(), "yymmdd"),

             parallel::Bool  = false,
             n_parts::Int    = 5_000,
             n_blocks::Int   = 1,
             n_mh_steps::Int = 1,

             λ::S = 2.1,
             n_Φ::Int64 = 300,

             resampling_method::Symbol = :systematic,
             threshold_ratio::S = 0.5,

             # Mutation Settings
             c::S = 0.5,       # step size
             α::S = 1.0,       # mixture proportion
             target::S = 0.25, # target accept rate

             use_fixed_schedule::Bool = true,
             tempering_target::S = 0.97,

             old_data::Matrix{S} = Matrix{S}(undef, size(data, 1), 0),
             old_cloud::Cloud = Cloud(0, 0),
             old_loglikelihood::Function = loglikelihood,
             old_vintage::String = "",
             smc_iteration::Int = 1,

             run_test::Bool = false,
             filestring_addl::Vector{String} = Vector{String}(),
             loadpath::String = "",
             savepath::String = "smc_cloud.jld2",
             particle_store_path::String = "smcsave.h5",
             save_intermediate::Bool = false,
             intermediate_stage_increment::Int = 10,
             continue_intermediate::Bool = false,
             intermediate_stage_start::Int = 0,
             tempered_update_prior_weight::S = 0.0,
             regime_switching::Bool = false,
             toggle::Bool = true,
             debug_assertion::Bool = false,
             log_prob_old_data::Float64 = 0.0) where {S<:AbstractFloat, U<:Number}

    ########################################################################################
    ### Settings
    ########################################################################################

    # Construct closure of mutation function so as to avoid issues with serialization
    # across workers with different Julia system images
    sendto(workers(), parameters = parameters)
    sendto(workers(), data = data)

    function mutation_closure(p::Vector{S}, d_μ::Vector{S}, d_Σ::Matrix{S},
                              n_free_para::Int,
                              blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}},
                              ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0, n_mh_steps::Int = 1,
                              old_data::T = Matrix{S}(undef, size(data, 1), 0)) where {S<:Float64, T<:Matrix}
        return mutation(loglikelihood, parameters, data, p, d_μ, d_Σ, n_free_para, blocks_free, blocks_all,
                        ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data,
                        old_loglikelihood = old_loglikelihood, regime_switching = regime_switching,
                        toggle = toggle)
    end
    @everywhere function mutation_closure(p::Vector{S}, d_μ::Vector{S}, d_Σ::Matrix{S},
                                          blocks_free::Vector{Vector{Int64}}, blocks_all::Vector{Vector{Int64}}, n_free_para::Int,
                                          ϕ_n::S, ϕ_n1::S; c::S = 1.0, α::S = 1.0, n_mh_steps::Int = 1,
                                          old_data::T = Matrix{S}(undef, size(data, 1), 0)) where {S<:Float64, T<:Matrix}
        return mutation(loglikelihood, parameters, data, p, d_μ, d_Σ, blocks_free, blocks_all, n_free_para,
                        ϕ_n, ϕ_n1; c = c, α = α, n_mh_steps = n_mh_steps, old_data = old_data,
                        old_loglikelihood = old_loglikelihood, regime_switching = regime_switching, toggle = toggle)
    end

    # Check that if there's a tempered update, old and current vintages are different
    tempered_update = !isempty(old_data) # Time tempering
    if (verbose != :none) && !(tempered_update & (old_vintage == data_vintage))
        @info "Old & current vintages the same!"
    end

    # General
    i   = 1             # Index tracking the stage of the algorithm
    j   = 2             # Index tracking the fixed_schedule entry ϕ_prop
    ϕ_n = ϕ_prop = 0.   # Instantiate ϕ_n and ϕ_prop variables

    resampled_last_period = false # Ensures proper resetting of ESS_bar after resample
    threshold             = threshold_ratio * n_parts
    #use_fixed_schedule = (tempering_target == 0.0)

    # Now count number of regime switching parameter values (excluding regime 1 values), i.e. if one parameter has 3 regimes, then add 2 to n_para
    n_para = length(parameters)
    for para in parameters
        if !isempty(para.regimes)
            for (ind, val) in para.regimes[:value]
                if ind != 1
                    n_para += 1
                end
            end
        end
    end

    fixed_para_inds = ModelConstructors.get_fixed_para_inds(parameters; regime_switching = regime_switching, toggle = toggle)
    free_para_inds  = ModelConstructors.get_free_para_inds( parameters; regime_switching = regime_switching, toggle = toggle)
    para_symbols    = [θ.key for θ in parameters]
    if regime_switching
        # Concatenate regime symbols for each extra regimes
        reg_switch_symbols = Vector{Symbol}(undef, n_para - length(parameters))
        ind = 0
        for θ in parameters
            if !isempty(θ.regimes)
                for i in 2:length(θ.regimes[:value])
                    ind += 1
                    reg_switch_symbols[ind] = Symbol(θ.key, "_reg$(i)")
                end
            end
        end
        push!(para_symbols, reg_switch_symbols...)
    end

    n_free_para = length(free_para_inds)
    @assert n_free_para > 0 "All model parameters are fixed!"

    #################################################################################
    ### Initialize Algorithm: Draws from prior
    #################################################################################
    println(verbose, :low, "\n\n SMC " * (testing ? "testing " : "") * "starts ....\n\n")

    if tempered_update
        # If user does not input Cloud object themselves, looks for cloud in loadpath.
        cloud = cloud_isempty(old_cloud) ? load(loadpath, "cloud") : old_cloud
        old_n_parts = length(cloud)

        if (tempered_update_prior_weight == 0.0) && (old_n_parts == n_parts)
            # Initialize settings first, with the ESS properly set since starting from an old estimation
            initialize_cloud_settings!(cloud; tempered_update = tempered_update,
                                       n_parts = n_parts, n_Φ = n_Φ, c = c, accept = target)

            # Update the old_loglh column in cloud.particles with the values
            # from the current loglh column in cloud.particles.
            # Then compute log-likelihood of *old* estimation's particles on *new* data and parameters.
            initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel,
                                    toggle = toggle)

        elseif (1. >= tempered_update_prior_weight > 0.) || (old_n_parts != n_parts)
            # Resample from bridge distribution
            n_to_resample = Int(round((1-tempered_update_prior_weight) * n_parts))
            n_from_prior  = n_parts - n_to_resample
            if n_to_resample > 0
                new_inds = resample(get_weights(cloud); n_parts = n_to_resample,
                                    method = resampling_method, parallel = parallel)

            else
                new_inds = Int64[] # empty
            end

            # Add in the samples from the old bridge.
            # Note the update_old_loglh! is commented out b/c
            # it is redundant, as we call initialize_likelihoods! later
            bridge_cloud = Cloud(n_para, n_to_resample)
            update_cloud!(bridge_cloud, cloud.particles[new_inds, :])
            update_loglh!(bridge_cloud, get_loglh(cloud)[new_inds])
            update_logprior!(bridge_cloud, get_logprior(cloud)[new_inds])
            update_old_loglh!(bridge_cloud, get_old_loglh(cloud)[new_inds])

            # Add to the cloud draws from the prior. Note that we are drawing from the current prior,
            # but evaluating the old log-likelihood function on the old data. Thus,
            # you should also make sure that old_loglikelihood can be evaluated on the
            # current set of parameters (which may be nontrivial if new parameters have been added
            # relative to the previous estimation)
            if n_from_prior > 0
                # Evaluate old_loglikelihood on old data and new draws from prior
                prior_cloud = Cloud(n_para, n_from_prior)
                initial_draw!(old_loglikelihood, parameters, old_data, prior_cloud, parallel = parallel,
                              regime_switching = regime_switching, toggle = toggle)

                # Instantiating Cloud object, update draws, loglh, & logprior.
                # Note the update_old_loglh! is commented out b/c
                # it is redundant, as we call initialize_likelihoods! later
                cloud = Cloud(n_para, n_to_resample + n_from_prior)
                update_cloud!(cloud, vcat(bridge_cloud.particles, prior_cloud.particles))
                update_loglh!(cloud, vcat(get_loglh(bridge_cloud), get_loglh(prior_cloud)))
                update_logprior!(cloud, vcat(get_logprior(bridge_cloud), get_logprior(prior_cloud)))
                update_old_loglh!(cloud, vcat(get_old_loglh(bridge_cloud), get_old_loglh(prior_cloud)))
            else
                cloud = bridge_cloud
            end

            # Update the old_loglh column in cloud.particles with the values
            # from the current loglh column in cloud.particles.
            # Then compute current log-likelihood of *old* estimation's particles on *new* data and parameters.
            initialize_likelihoods!(loglikelihood, parameters, data, cloud; parallel = parallel,
                                    toggle = toggle)

            # Ensure no particles yielding -Inf loglhs are kept.
            # These "bad" particles emerge b/c particles from an old estimation
            # won't necessarily imply good loglhs when evaluated on
            # new data and parameters.
            zero_bad_loglh_weights!(cloud)
            normalized_weights = normalize_weights!(cloud) # need to renormalize to ensure weights sum to n_parts

            # Resample weights to remove the -Inf loglh
            new_inds = resample(normalized_weights/n_parts; method = resampling_method,
                                parallel = parallel)
            cloud.particles = [deepcopy(cloud.particles[k,j]) for k in new_inds,
                               j=1:size(cloud.particles, 2)]
            reset_weights!(cloud)

            # Since there was a resampling, set ESS = n_parts
            push!(cloud.ESS, n_parts)

            # Initialize remaining cloud settings
            initialize_cloud_settings!(cloud; tempered_update = tempered_update,
                                       n_parts = n_parts, n_Φ = n_Φ, c = c, accept = target)
        else
            throw(DomainError("The keyword tempered_update_prior_weight must be within the interval [0, 1] but " *
                              "is currently set to $(tempered_update_prior_weight)"))
        end
    elseif continue_intermediate
        cloud = load(loadpath, "cloud")
    else
        # Initialization of Particle Array Cloud
        cloud = Cloud(n_para, n_parts)

        # Instantiating Cloud object, update draws, loglh, & logprior
        initial_draw!(loglikelihood, parameters, data, cloud; parallel = parallel, regime_switching = regime_switching,
                      toggle = toggle)
        initialize_cloud_settings!(cloud; tempered_update = tempered_update,
                                   n_parts = n_parts, n_Φ = n_Φ, c = c, accept = target)
    end

    # Fixed schedule for construction of ϕ_prop
    if use_fixed_schedule
        cloud.tempering_schedule = ((collect(1:n_Φ) .- 1) / (n_Φ-1)) .^ λ
    else
        proposed_fixed_schedule  = ((collect(1:n_Φ) .- 1) / (n_Φ-1)) .^ λ
    end

    # Instantiate incremental and normalized weight matrices for logMDD calculation
    if continue_intermediate
        w_matrix = load(loadpath, "w")
        W_matrix = load(loadpath, "W")
        j        = load(loadpath, "j")
        i        = cloud.stage_index
        c        = cloud.c
        ϕ_prop   = proposed_fixed_schedule[j]
    else
        w_matrix = zeros(n_parts, 1)
        W_matrix = tempered_update ? (sum(get_weights(cloud)) <= 1.0 ?
                                      get_weights(cloud) * n_parts : get_weights(cloud)) :
                                      fill(1,(n_parts, 1))
    end

    # Printing
    init_stage_print(cloud, para_symbols; verbose = verbose,
                     use_fixed_schedule = use_fixed_schedule)
    println(verbose, :low, "\n\n SMC recursion starts... \n\n")

    #################################################################################
    ### Recursion
    #################################################################################
    while ϕ_n < 1.
        start_time = time_ns()
        cloud.stage_index = i += 1

        #############################################################################
        ### Setting ϕ_n (either adaptively or by the fixed schedule)
        #############################################################################
        ϕ_n1 = cloud.tempering_schedule[i-1]

        if use_fixed_schedule
            ϕ_n = cloud.tempering_schedule[i]
        else
            ϕ_n, resampled_last_period, j, ϕ_prop = solve_adaptive_ϕ(cloud,
                                                       proposed_fixed_schedule,
                                                       i, j, ϕ_prop, ϕ_n1,
                                                       tempering_target,
                                                       resampled_last_period)
        end

        #############################################################################
        ### Step 1: Correction
        #############################################################################

        # Calculate incremental weights (if no old data, get_old_loglh(cloud) = 0)
        if tempered_update_prior_weight == 0.0
            incremental_weights = exp.((ϕ_n1 - ϕ_n) * get_old_loglh(cloud) +
                                       (ϕ_n - ϕ_n1) * get_loglh(cloud))
        elseif tempered_update_prior_weight == 1.0
            incremental_weights = exp.((ϕ_n - ϕ_n1) * get_loglh(cloud))
        else
            incremental_weights = exp.((ϕ_n1 - ϕ_n) * log.((exp.(get_old_loglh(cloud) .- log_prob_old_data .+
                                                       log(1-tempered_update_prior_weight)) .+ tempered_update_prior_weight)) +
                                       (ϕ_n - ϕ_n1) * get_loglh(cloud))
        end

        # Update weights
        update_weights!(cloud, incremental_weights)
        mult_weights = get_weights(cloud)

        # Normalize weights
        normalized_weights = normalize_weights!(cloud) #get_weights(cloud)

        w_matrix = hcat(w_matrix, incremental_weights)
        W_matrix = hcat(W_matrix, normalized_weights)

        ##############################################################################
        ### Step 2: Selection
        ##############################################################################

        # Calculate the degeneracy/effective sample size metric
        push!(cloud.ESS, n_parts ^ 2 / sum(normalized_weights .^ 2))

        # Check whether ESS is a NaN and throws an assertion error if it is.
        # In many cases, the problem is that there too few particles.
        check_nan_ess(cloud, i, incremental_weights,
                      normalized_weights, savepath, debug_assertion)

        # Resample if degeneracy/ESS metric falls below the accepted threshold
        if (cloud.ESS[i] < threshold)

            # Resample according to particle weights, uniformly reset weights to 1/n_parts
            new_inds = resample(normalized_weights/n_parts; method = resampling_method,
                                parallel = parallel)
            cloud.particles = [deepcopy(cloud.particles[k,j]) for k in new_inds,
                               j=1:size(cloud.particles, 2)]
            reset_weights!(cloud)
            cloud.resamples += 1
            resampled_last_period = true
            W_matrix[:, i] .= 1
        end

        ##############################################################################
        ### Step 3: Mutation
        ##############################################################################

        # Calculate adaptive c-step for use as scaling coefficient in mutation MH step
        c = c * (0.95 + 0.10 * exp(16.0 * (cloud.accept - target)) /
                 (1.0 + exp(16.0 * (cloud.accept - target))))
        cloud.c = c

        θ_bar = weighted_mean(cloud)
        R     = weighted_cov(cloud)

        # Ensures marix is positive semi-definite symmetric
        # (not off due to numerical error) and values haven't changed
        R_fr = (R[free_para_inds, free_para_inds] + R[free_para_inds, free_para_inds]') / 2.

        # MvNormal centered at ̄θ with var-cov ̄Σ, subsetting out the fixed parameters
        θ_bar_fr = θ_bar[free_para_inds]

        # Generate random parameter blocks
        blocks_free = generate_free_blocks(n_free_para, n_blocks)
        blocks_all  = generate_all_blocks(blocks_free, free_para_inds)

        new_particles = if parallel
            @distributed (hcat) for k in 1:n_parts
                mutation_closure(cloud.particles[k, :], θ_bar_fr, R_fr, n_free_para,
                                 blocks_free, blocks_all, ϕ_n, ϕ_n1; c = c, α = α,
                                 n_mh_steps = n_mh_steps, old_data = old_data)
            end
        else
            hcat([mutation_closure(cloud.particles[k, :], θ_bar_fr, R_fr, n_free_para,
                                   blocks_free, blocks_all, ϕ_n, ϕ_n1; c = c,
                                   α = α, n_mh_steps = n_mh_steps,
                                   old_data = old_data) for k=1:n_parts]...)
        end
        update_cloud!(cloud, new_particles)
        update_acceptance_rate!(cloud)

        ##############################################################################
        ### Timekeeping and Output Generation
        ##############################################################################
        stage_sampling_time = Float64((time_ns() - start_time) * 1e-9)
        cloud.total_sampling_time += stage_sampling_time

        end_stage_print(cloud, para_symbols; verbose = verbose,
                        use_fixed_schedule = use_fixed_schedule)

        if run_test && (i == 3)
            break
        end

        if mod(cloud.stage_index, intermediate_stage_increment) == 0 && save_intermediate
            jldopen(replace(savepath, ".jld2" => "_stage=$(cloud.stage_index).jld2"),
                    true, true, true, IOStream) do file
                write(file, "cloud", cloud)
                write(file, "w", w_matrix)
                write(file, "W", W_matrix)
                write(file, "j", j)
            end
        end
    end

    ##################################################################################
    ### Saving data
    ##################################################################################
    if !testing
        simfile = h5open(particle_store_path, "w")
        particle_store = isdefined(HDF5, :create_dataset) ?
            HDF5.create_dataset(simfile, "smcparams", datatype(Float64), dataspace(n_parts, n_para)) :
            HDF5.d_create(simfile, "smcparams", datatype(Float64), dataspace(n_parts, n_para))

        for k in 1:n_parts; particle_store[k,:] = cloud.particles[k, 1:n_para] end
        close(simfile)
        jldopen(savepath, true, true, true, IOStream) do file
            write(file, "cloud", cloud)
            write(file, "w", w_matrix)
            write(file, "W", W_matrix)
        end
    end
end


text:
# Sample a gas fitness DGP and estimate the parameters for different lenghts of
# time series

using StatsBase
using JLD2

initWorkers()

using DynNets
using StatsFuns
@everywhere using DynNets
@everywhere using StatsFuns
N_est = 300

N = 10 ; Ngroups = 1; GTV =Ngroups; useStartVal = true
NnodesPerGroup = floor(N/Ngroups);groupsInds=zeros(Int64,N)
#Distribuisci i nodi tra i gruppi in parti uguali
for i =0:Ngroups-2 groupsInds[1 + Int64(NnodesPerGroup*i):Int64(NnodesPerGroup*(i+1))] = i+1 end;
groupsInds[find(groupsInds .==0)] = Ngroups
indTvGroups = falses(Ngroups);indTvGroups[1:GTV] = true
Tvals = [1000]
N_ind = length(Tvals)[1]

#the values for the unconidtional means to be used in the DGP are obtained estimating
#an equally spaced degree sequence
tmpDegs = Array{Real,1}(linspace(2,N-2,N))'; tmpdegsT = repmat(tmpDegs,2,1)#ones(N)*N/2
tmpModel =  DynNets.SdErgm1(tmpdegsT,[zeros(N),zeros(GTV),zeros(GTV)],
                                    zeros(Int,N),
                                    falses(1))

#

#Initialize the variable to store
for i =1: size(tmpModel.Par)[1]
    ex1 = "estPar_$(i)  = SharedArray{Float64}(zeros(N_ind,N_est, size(tmpModel.Par[$(i)])[1] ))"
    ex2 = "simPar_$(i) = SharedArray{Float64}(zeros(N_ind,N_est, size(tmpModel.Par[$(i)])[1] ))"
    eval(parse(ex1))
    eval(parse(ex2))
 end
 est_times = SharedArray{Float64}(zeros(N_ind,N_est))
 est_conv_flag = SharedArray{Bool,2}(N_ind,N_est)
 est_flag = SharedArray{Bool,2}(N_ind,N_est)

# Set the value of the parameters common to all replicate of DGP+Estimation
tmp,~ = DynNets.estimate(tmpModel )
realUm = tmp[1]
sim_B = Array(linspace(0.4,0.9,GTV+1)[2:end]) #
sim_A = Array(linspace(0.02,0.2,GTV+1)[2:end]) #
sim_W = realUm.*(1-sim_B)
@progress "T"  for ind = 1:N_ind
 # ind=1
   T = Tvals[ind]
    @sync @parallel   for n = 1:N_est
  #n=1

        simPar_3[ind,n,1:GTV] = sim_A
        simPar_2[ind,n,1:GTV] = sim_B
        simPar_1[ind,n,1:N] = sim_W
        Y_T,fit_T  =  DynNets.sampleDynNetsBin1(T,"Fit-Multi-gas";N=N,
                                                    indTvGroups=indTvGroups,
                                                    groupsInds=groupsInds,
                                                    AgasGroups =sim_A,
                                                    BgasGroups = sim_B,
                                                    WgasNodes = sim_W )
        degsT = squeeze(sum(Y_T,2),2)
        model_data = DynNets.SdErgm1(   degsT,
                                                [zeros(N),zeros(Ngroups),zeros(Ngroups)],
                                                groupsInds,
                                                indTvGroups)
        #estimate
        useStartVal ?   start_val = [sim_W, sim_B, sim_A] : start_val = zeros(3,3)
        est_times[ind,n] = @elapsed  hat_pars , conv_flag =
                            DynNets.estimate(model_data,start_values = start_val )

        indTvNodes =   .!(groupsInds .==0);
        estPar_1[ind,n,:] = hat_pars[1]
        estPar_2[ind,n,:] = hat_pars[2]
        estPar_3[ind,n,:] = hat_pars[3]
        est_conv_flag[ind,n] = conv_flag
        est_flag[ind,n] = true

    end

    # save
    save_fold = "./data/estimatesTest/"
    file_name = "EstTest_$(N)_$(GTV)_$(Tvals)_$(N_est)_$(useStartVal).jld"
    save_path = save_fold*file_name#
        @save(save_path,groupsInds,indTvGroups,simPar_1,simPar_2,simPar_3,estPar_1,estPar_2,estPar_3,est_conv_flag,est_times)
    println([T,sum(est_times[ind,:]),sum(est_conv_flag[ind,:])/sum(est_flag[ind,:])])
end

##

 ##


text:
module ex1 # from Hypre examples
# direct translation from the corresponding C source. See
# src/hypre-2.9.0b/src/examples/ex1.c under the deps directory.
# Run with mpirun -np 2 julia ex1.jl
#
using Hypre
using MPI

function runex1()
  const grid = [Hypre.hypre_StructGrid_struct()]
  const stencil = [Hypre.hypre_StructStencil_struct()]
  const A = [Hypre.hypre_StructMatrix_struct()]
  const b = [Hypre.hypre_StructVector_struct()]
  const x = [Hypre.hypre_StructVector_struct()]
  const solver = [Hypre.hypre_StructSolver_struct()]

  #   /* Initialize MPI */
  MPI.Init()
  const comm = MPI.COMM_WORLD
  const myid = MPI.Comm_rank(comm)
  const num_procs = MPI.Comm_size(comm)

  if num_procs != 2
    if myid == 0
      println("Must run with 2 processors!")
    end
    MPI.Finalize()
    return(0)
  end
  # #   /* 1. Set up a grid. Each processor describes the piece
  # #   of the grid that it owns. */
  # #     /* Create an empty 2D grid object */
  Hypre.HYPRE_StructGridCreate(comm, 2, grid)

  # #     /* Add boxes to the grid */
  if myid == 0
    ilower = Hypre.HYPRE_Int[-3, 1]
    iupper = Hypre.HYPRE_Int[-1, 2]
    Hypre.HYPRE_StructGridSetExtents(grid, ilower, iupper)
  elseif myid == 1
    ilower = Hypre.HYPRE_Int[0, 1]
    iupper = Hypre.HYPRE_Int[2, 4]
    Hypre.HYPRE_StructGridSetExtents(grid, ilower, iupper)
  end
  if myid == 0
    println(grid)
  end
#   # #     /* This is a collective call finalizing the grid assembly.
#   # #     The grid is now ``ready to be used'' */
#   Hypre.HYPRE_StructGridAssemble(gridptr)

#   # #   /* 2. Define the discretization stencil */
#   # #     /* Create an empty 2D, 5-pt stencil object */
#   Hypre.HYPRE_StructStencilCreate(2, 5, stencilptr)

  # #     /* Define the geometry of the stencil. Each represents a
  # #     relative offset (in the index space). */
  offsets = reshape(Hypre.HYPRE_Int[0, 0, -1, 0, 1, 0, 0, -1, 0, 1], 2, 5)

  # #       /* Assign each of the 5 stencil entries */
  for entry in 1:5
    #       Hypre.HYPRE_StructStencilSetElement(stencil, entry, offsets[entry]);
  end
  # #   /* 3. Set up a Struct Matrix */
  # #     /* Create an empty matrix object */
  #     Hypre.HYPRE_StructMatrixCreate(MPI.COMM_WORLD, grid, stencil, &A);

  # #     /* Indicate that the matrix coefficients are ready to be set */
  #     Hypre.HYPRE_StructMatrixInitialize(A);

  # #     /* Set the matrix coefficients.  Each processor assigns coefficients
  # #     for the boxes in the grid that it owns. Note that the coefficients
  # #     associated with each stencil entry may vary from grid point to grid
  # #     point if desired.  Here, we first set the same stencil entries for
  # #     each grid point.  Then we make modifications to grid points near
  # #     the boundary. */
  if myid == 0
    ilower = Hypre.HYPRE_Int[-3, 1]
    iupper = Hypre.HYPRE_Int[-1, 2]
    stencil_indices = Hypre.HYPRE_Int[0, 1, 2, 3, 4] #/* labels for the stencil entries -
    # #         these correspond to the offsets defined above */
    nentries = length(stencil_indices)
    nvalues  = 30 #/* 6 grid points, each with 5 stencil entries */
    values = zeros(nvalues)

    # #       /* We have 6 grid points, each with 5 stencil entries */
    for i in 1:nentries:nvalues
      values[i] = 4.0
      for j in 1:(nentries - 1)
        values[i + j] = -1.0
      end
    end

    #       Hypre.HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, nentries,
    #                                      stencil_indices, values);
  elseif myid == 1
    ilower = Hypre.HYPRE_Int[0, 1]
    iupper = Hypre.HYPRE_Int[2, 4]
    stencil_indices = Hypre.HYPRE_Int[0, 1, 2, 3, 4]
    nentries = length(stencil_indices)
    nvalues  = 60 #/* 12 grid points, each with 5 stencil entries */
    values = zeros(nvalues)
    for i in 1:nentries:nvalues
      values[i] = 4.0
      for j in 1:(nentries - 1)
        values[i + j] = -1.0
      end
    end

    #       Hypre.HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, nentries,
    #                                      stencil_indices, values);
  end

  #     /* Set the coefficients reaching outside of the boundary to 0 */
  if myid == 0
    values = zeros(3)

    # #         /* values below our box */
    ilower = Hypre.HYPRE_Int[-3, 1]
    iupper = Hypre.HYPRE_Int[-1, 1]
    stencil_indices = Hypre.HYPRE_Int[3,]
    #         Hypre.HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
    #                                        stencil_indices, values);

    # #         /* values to the left of our box */
    ilower = Hypre.HYPRE_Int[-3, 1]
    iupper = Hypre.HYPRE_Int[-3, 2]
    stencil_indices = Hypre.HYPRE_Int[1,]
    #         Hypre.HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
    #                                        stencil_indices, values);

    # #         /* values above our box */
    ilower = Hypre.HYPRE_Int[-3, 2]
    iupper = Hypre.HYPRE_Int[-1, 2]
    stencil_indices = Hypre.HYPRE_Int[4,]
    #         Hypre.HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
    #                                        stencil_indices, values);
  elseif myid == 1
    values = zeros(4)
    #       {
    #         /* values below our box */
    ilower = Hypre.HYPRE_Int[0, 1]
    iupper = Hypre.HYPRE_Int[2, 1]
    stencil_indices = Hypre.HYPRE_Int[3,]
    #         Hypre.HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
    #                                        stencil_indices, values);
    #         }
    #       {
    #         /* values to the right of our box */
    ilower = Hypre.HYPRE_Int[2, 1]
    iupper = Hypre.HYPRE_Int[2, 4]
    stencil_indices = Hypre.HYPRE_Int[2,]
    #         Hypre.HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
    #                                        stencil_indices, values);
    #         }
    #       {
    #         /* values above our box */
    ilower = Hypre.HYPRE_Int[0, 4]
    iupper = Hypre.HYPRE_Int[2, 4]
    stencil_indices = Hypre.HYPRE_Int[4,]
    #         Hypre.HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
    #                                        stencil_indices, values);
    #         }
    #       {
    # #         /* values to the left of our box
    # #         (that do not border the other box on proc. 0) */
    ilower = Hypre.HYPRE_Int[0, 3]
    iupper = Hypre.HYPRE_Int[0, 4]
    stencil_indices = Hypre.HYPRE_Int[1,]
    #            Hypre.HYPRE_StructMatrixSetBoxValues(A, ilower, iupper, 1,
    #                                           stencil_indices, values);
    #            }
  end

  # #            /* This is a collective call finalizing the matrix assembly.
  # #            The matrix is now ``ready to be used'' */
  #            Hypre.HYPRE_StructMatrixAssemble(A);
  #            }

  # #            /* 4. Set up Struct Vectors for b and x.  Each processor sets the vectors
  # #              corresponding to its boxes. */
  #              {
  # #                /* Create an empty vector object */
  #                Hypre.HYPRE_StructVectorCreate(MPI.COMM_WORLD, grid, &b);
  #                Hypre.HYPRE_StructVectorCreate(MPI.COMM_WORLD, grid, &x);

  # #                /* Indicate that the vector coefficients are ready to be set */
  #                Hypre.HYPRE_StructVectorInitialize(b);
  #                Hypre.HYPRE_StructVectorInitialize(x);

  # #                /* Set the vector coefficients */
  if myid == 0
    ilower = Hypre.HYPRE_Int[-3, 1]
    iupper = Hypre.HYPRE_Int[-1, 2]
    values = ones(6) #/* 6 grid points */
    #                  Hypre.HYPRE_StructVectorSetBoxValues(b, ilower, iupper, values);

    fill!(values, 0.0)
    #                  Hypre.HYPRE_StructVectorSetBoxValues(x, ilower, iupper, values);
  elseif myid == 1
    ilower = Hypre.HYPRE_Int[0, 1]
    iupper = Hypre.HYPRE_Int[2, 4]
    values = ones(12)
    #                  Hypre.HYPRE_StructVectorSetBoxValues(b, ilower, iupper, values);

    fill!(values, 0.0)
    #                  Hypre.HYPRE_StructVectorSetBoxValues(x, ilower, iupper, values);
  end

  # #                /* This is a collective call finalizing the vector assembly.
  # #                The vectors are now ``ready to be used'' */
  #                  Hypre.HYPRE_StructVectorAssemble(b);
  #                  Hypre.HYPRE_StructVectorAssemble(x);
  #                }

  #              #     5. Set up and use a solver (See the Reference Manual for descriptions of all of the options.)
  # #                /* Create an empty PCG Struct solver */
  #                Hypre.HYPRE_StructPCGCreate(MPI.COMM_WORLD, &solver);

  # #                /* Set some parameters */
  #                Hypre.HYPRE_StructPCGSetTol(solver, 1e-6); #/* convergence tolerance */
  #                Hypre.HYPRE_StructPCGSetPrintLevel(solver, 2); #/* amount of info. printed */

  # #                /* Setup and solve */
  #                Hypre.HYPRE_StructPCGSetup(solver, A, b, x);
  #                Hypre.HYPRE_StructPCGSolve(solver, A, b, x);

  # #              /* Free memory */
  #              Hypre.HYPRE_StructGridDestroy(grid);
  #              Hypre.HYPRE_StructStencilDestroy(stencil);
  #              Hypre.HYPRE_StructMatrixDestroy(A);
  #              Hypre.HYPRE_StructVectorDestroy(b);
  #              Hypre.HYPRE_StructVectorDestroy(x);
  #              Hypre.HYPRE_StructPCGDestroy(solver);

  # #              /* Finalize MPI */
  MPI.Finalize();

  return(0)
end #function

# function runex1mpi()
#   const nprocesses = 2
#   println(readall(`mpirun -np $nprocesses julia $(@__FILE__) go`))
# end
# function __init__()
#   if in("noMPI", ARGS) || in("go", ARGS)
#     gc_disable()
    exitcode = runex1()
#     gc_enable()
    exit(exitcode)
#   else
#     runex1mpi()
#   end
# end

end #module


text:
module Perf
using BenchmarkTools
using BenchmarkTools: Benchmark, TrialEstimate
using Setfield
using Test
using InteractiveUtils
using StaticArrays

struct AB{A,B}
    a::A
    b::B
end

function lens_set_a((obj, val))
    @set obj.a = val
end

function hand_set_a((obj, val))
    AB(val, obj.b)
end

function lens_set_ab((obj, val))
    @set obj.a.b = val
end

function hand_set_ab((obj, val))
    a = AB(obj.a.a, val)
    AB(a, obj.b)
end

function lens_set_a_and_b((obj, val))
    o1 = @set obj.a = val
    o2 = @set o1.b = val
end

function hand_set_a_and_b((obj, val))
    AB(val, val)
end

function lens_set_i((obj, val, i))
    @inbounds (@set obj[i] = val)
end

function hand_set_i((obj, val, i))
    @inbounds Base.setindex(obj, val, i)
end

function benchmark_lens_vs_hand(b_lens::Benchmark, b_hand::Benchmark)

    te_hand = minimum(run(b_lens))
    te_lens = minimum(run(b_hand))
    @show te_lens
    @show te_hand
    @test te_lens.memory == te_hand.memory
    @test te_lens.allocs == te_hand.allocs
    @test te_lens.time <= 2*te_hand.time
end

function uniquecounts(iter)
    ret = Dict{eltype(iter), Int}()
    for x in iter
        ret[x] = get!(ret, x, 0) + 1
    end
    ret
end

function test_ir_lens_vs_hand(info_lens::Core.CodeInfo,
                              info_hand::Core.CodeInfo)

    heads(info) = [ex.head for ex in info.code if ex isa Expr]

    # test no needless kinds of operations
    heads_lens = heads(info_lens)
    heads_hand = heads(info_hand)
    @test Set(heads_lens) == Set(heads_hand)

    # test no intermediate objects or lenses
    @test count(==(:new), heads_lens) == count(==(:new), heads_hand)

    # this test might be too strict
    @test uniquecounts(heads_lens) == uniquecounts(heads_hand)
end

let
    obj = AB(AB(1,2), :b)
    val = (1,2)
    @testset "$(setup.lens)" for setup in [
            (lens=lens_set_a,           hand=hand_set_a,       args=(obj, val)),
            (lens=lens_set_a,           hand=hand_set_a,       args=(obj, val)),
            (lens=lens_set_ab,          hand=hand_set_ab,      args=(obj, val)),
            (lens=lens_set_a_and_b,     hand=hand_set_a_and_b, args=(obj, val)),
            (lens=lens_set_i,           hand=hand_set_i,
             args=(@SVector[1,2], 10, 1))
            ]
        f_lens = setup.lens
        f_hand = setup.hand
        args = setup.args

        @assert f_hand(args) == f_lens(args)

        @testset "IR" begin
            info_lens, _ = @code_typed f_lens(args)
            info_hand, _ = @code_typed f_hand(args)
            test_ir_lens_vs_hand(info_lens, info_hand)
        end

        @testset "benchmark" begin
            b_lens = @benchmarkable $f_lens($args)
            b_hand = @benchmarkable $f_hand($args)
            benchmark_lens_vs_hand(b_lens, b_hand)
        end
    end
end

function compose_left_assoc(obj, val)
    l = @lens ((_.a∘_.b)∘_.c)∘_.d
    set(obj, l, val)
end

function compose_right_assoc(obj, val)
    l = @lens _.a∘(_.b∘(_.c∘_.d))
    set(obj, l, val)
end
function compose_default_assoc(obj, val)
    l = @lens _.a.b.c.d
    set(obj, l, val)
end
@testset "Lens composition compiler prefered associativity" begin

    obj = (a=(b=(c=(d=1,d2=2),c2=2),b2=3), a2=2)
    val = 2.2
    @test compose_left_assoc(obj, val) == compose_default_assoc(obj, val)
    @test compose_right_assoc(obj, val) == compose_default_assoc(obj, val)

    b_default = minimum(@benchmark compose_default_assoc($obj, $val))
    println("Default associative composition: $b_default")
    b_left    = minimum(@benchmark compose_left_assoc($obj, $val)   )
    println("Left associative composition: $b_left")
    b_right   = minimum(@benchmark compose_right_assoc($obj, $val)  )
    println("Right associative composition: $b_right")

    @test b_default.allocs == 0
    @test b_right.allocs == 0
    @test_broken b_left.allocs == 0

    @test b_left.time > 2b_default.time
    @test b_right.time ≈ b_default.time rtol=0.8
end

end


text:
using gazel
using Test

@testset "gazel.jl" begin
    # Write your tests here.
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
    Experimental

!!! warning
    Types, methods, or macros defined in this module are experimental and subject
    to change and will not have deprecations. Caveat emptor.
"""
module Experimental

using Base: Threads, sync_varname

"""
    Const(A::Array)

Mark an Array as constant/read-only. The invariant guaranteed is that you will not
modify an Array (through another reference) within an `@aliasscope` scope.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
struct Const{T,N} <: DenseArray{T,N}
    a::Array{T,N}
end

Base.IndexStyle(::Type{<:Const}) = IndexLinear()
Base.size(C::Const) = size(C.a)
Base.axes(C::Const) = axes(C.a)
@eval Base.getindex(A::Const, i1::Int) =
    (Base.@_inline_meta; Core.const_arrayref($(Expr(:boundscheck)), A.a, i1))
@eval Base.getindex(A::Const, i1::Int, i2::Int, I::Int...) =
  (Base.@_inline_meta; Core.const_arrayref($(Expr(:boundscheck)), A.a, i1, i2, I...))

"""
    @aliasscope expr

Allows the compiler to assume that all `Const`s are not being modified through stores
within this scope, even if the compiler can't prove this to be the case.

!!! warning
    Experimental API. Subject to change without deprecation.
"""
macro aliasscope(body)
    sym = gensym()
    quote
        $(Expr(:aliasscope))
        $sym = $(esc(body))
        $(Expr(:popaliasscope))
        $sym
    end
end

function sync_end(refs)
    local c_ex
    defined = false
    t = current_task()
    cond = Threads.Condition()
    lock(cond)
    nremaining = length(refs)
    for r in refs
        schedule(Task(()->begin
            try
                wait(r)
                lock(cond)
                nremaining -= 1
                nremaining == 0 && notify(cond)
                unlock(cond)
            catch e
                lock(cond)
                notify(cond, e; error=true)
                unlock(cond)
            end
        end))
    end
    wait(cond)
    unlock(cond)
end

"""
    Experimental.@sync

Wait until all lexically-enclosed uses of `@async`, `@spawn`, `@spawnat` and `@distributed`
are complete, or at least one of them has errored. The first exception is immediately
rethrown. It is the responsibility of the user to cancel any still-running operations
during error handling.

!!! Note
    This interface is experimental and subject to change or removal without notice.
"""
macro sync(block)
    var = esc(sync_varname)
    quote
        let $var = Any[]
            v = $(esc(block))
            sync_end($var)
            v
        end
    end
end

end

max_stars_repo_name: fabmazz/julia

max_issues_repo_name: fabmazz/julia

max_forks_repo_name: fabmazz/julia


text:
using RvSpectMLBase
using Test

#=
 This only tests functions from src/utils that don't depend on data structures.
 test/data_structure.jl should test functions that depend on the data structure.
 Functions using multiple data structures should be tested in the more complex data structure's file.
=#

@testset "Simple Physics Functions" begin

    @testset "Util" begin
        @test RvSpectMLBase.calc_doppler_factor(10) ≈ 1.0000000333564094
        @test RvSpectMLBase.predict_intrinsic_stellar_line_width(5780) ≈ 9883.42046054907
    end

    @testset "calc_snr" begin
        @test RvSpectMLBase.calc_snr(8.0,4.0) ≈ 4
        @test RvSpectMLBase.calc_snr(8*ones(100),4*ones(100)) ≈ 40
    end

end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
batchify(resource, batch_size) =
begin
    hm_batches = trunc(Int, length(resource)/batch_size)
    hm_leftover = length(resource)%batch_size
    batches = []
    if batch_size > length(resource)
        push!(batches, resource)
    else
        for i in 1:hm_batches
            push!(batches, resource[(i-1)*batch_size+1:i*batch_size])
        end
        if hm_leftover != 0
            push!(batches, resource[(end-hm_leftover)-1:end])
        end
    end
batches
end

import_data(file) =
    open(file) do f
        println("importing $file")
        data = []
        sample = []
        for line in eachline(f)
            if line == ";"
                if length(sample) != 0
                    push!(data, sample)
                    sample = []
                end
            else
                nrs = split(line, " ")
                arr = []
                for nr in nrs
                    try
                        push!(arr, parse(Float32, nr))
                    catch end
                end
                # arr = [parse(Float32, nr) for nr in nrs]
                push!(sample, reshape(arr, 1, length(arr)))
            end
        end
        if length(sample) != 0
            push!(data, sample)
        end
        println("obtained $(length(data)) samples.")
    data
    end

save_model(model) =
begin
    open("model.txt", "w+") do file
        for mfield in fieldnames(Model)
            layer = getfield(model, mfield)
            for lfield in fieldnames(Layer)
                if lfield != :state
                    items = value(getfield(layer, lfield))
                    for i in items
                        write(file, string(i) * " ")
                    end
                    write(file, "\n")
                end
            end
        end
    end
end

load_model() =
begin

end


text:
#####
##### `Dataset`
#####

struct Dataset
    path::String
    header::Header
    recordings::Dict{UUID,Recording}
end

"""
    Dataset(path; create=false)

Return a `Dataset` instance that contains all metadata necessary to read and
write to the Onda dataset stored at `path`. Note that this constuctor loads all
the `Recording` objects contained in `path/recordings.msgpack.zst`.

If `create` is `true`, then an empty Onda dataset will be created at `path`.
"""
function Dataset(path; create::Bool=false)
    path = rstrip(abspath(path), '/')
    samples_path = joinpath(path, "samples")
    if create
        if isdir(path)
            isempty(readdir(path)) ||
                throw(ArgumentError("cannot create dataset at $path: directory exists and is nonempty"))
        else
            mkdir(path)
        end
        initial_header = Header(ONDA_FORMAT_VERSION, true)
        initial_recordings = Dict{UUID,Recording}()
        write_recordings_file(path, initial_header, initial_recordings)
    elseif !isdir(path)
        throw(ArgumentError("$path is not a valid Onda dataset"))
    end
    !isdir(samples_path) && mkdir(samples_path)
    header, recordings = read_recordings_file(path)
    return Dataset(path, header, recordings)
end

"""
    save_recordings_file(dataset::Dataset)

Overwrite `joinpath(dataset.path, "recordings.msgpack.zst")` with the contents
of `dataset.recordings`.
"""
function save_recordings_file(dataset::Dataset)
    return write_recordings_file(dataset.path, dataset.header, dataset.recordings)
end

#####
##### `merge!`
#####

"""
    merge!(destination::Dataset, datasets::Dataset...; only_recordings::Bool=false)

Write all filesystem content and the `recordings` field of each `Dataset` in
`datasets` to `destination`.

If any filesystem content has a name that conflicts with existing filesystem
content in `destination`, this function will throw an error. An error will also
be thrown if this function encounters multiple recordings with the same UUID.

If `only_recordings` is `true`, then only the `recordings` field of each `Dataset`
is merged, such that no filesystem content is read or written.

NOTE: This function is currently only implemented when `only_recordings = true`.
"""
function Base.merge!(destination::Dataset, datasets::Dataset...;
                     only_recordings::Bool=false)
    only_recordings ||
        error("`merge!(datasets::Dataset...; only_recordings=false)` is not yet implemented")
    for dataset in datasets
        for uuid in keys(dataset.recordings)
            if haskey(destination.recordings, uuid)
                throw(ArgumentError("recording $uuid already exists in the destination dataset"))
            end
        end
        merge!(destination.recordings, dataset.recordings)
    end
    return destination
end

#####
##### `samples_path`
#####

"""
    samples_path(dataset::Dataset, uuid::UUID)

Return the samples subdirectory path corresponding to the recording specified by `uuid`.
"""
function samples_path(dataset::Dataset, uuid::UUID)
    return joinpath(dataset.path, "samples", string(uuid))
end

"""
    samples_path(dataset::Dataset, uuid::UUID, name::Symbol,
                 file_extension=dataset.recordings[uuid].signals[name].file_extension)

Return the samples file path corresponding to the signal named `name` within the
recording specified by `uuid`.
"""
function samples_path(dataset::Dataset, uuid::UUID, name::Symbol,
                      file_extension=dataset.recordings[uuid].signals[name].file_extension)
    file_name = string(name, ".", file_extension)
    return joinpath(samples_path(dataset, uuid), file_name)
end

#####
##### `create_recording!`
#####

"""
    create_recording!(dataset::Dataset, uuid::UUID=uuid4())

Create `uuid::UUID => recording::Recording`, add the pair to `dataset.recordings`,
and return the pair.
"""
function create_recording!(dataset::Dataset, uuid::UUID=uuid4())
    if haskey(dataset.recordings, uuid)
        throw(ArgumentError("recording with UUID $uuid already exists in dataset"))
    end
    recording = Recording(Dict{Symbol,Signal}(), Set{Annotation}())
    dataset.recordings[uuid] = recording
    mkpath(samples_path(dataset, uuid))
    return uuid => recording
end

#####
##### `load`
#####

"""
    load(dataset::Dataset, uuid::UUID, name::Symbol[, span::AbstractTimeSpan])

Load and return the `Samples` object corresponding to the signal named `name`
in the recording specified by `uuid`.

If `span` is provided, this function returns the equivalent of
`load(dataset, uuid, name)[:, span]`, but potentially avoids loading the entire
signal's worth of sample data if the underlying signal file format supports
partial access/random seeks.

See also: [`deserialize_lpcm`](@ref)
"""
function load(dataset::Dataset, uuid::UUID, name::Symbol, span::AbstractTimeSpan...)
    signal = dataset.recordings[uuid].signals[name]
    path = samples_path(dataset, uuid, name, signal.file_extension)
    return load_samples(path, signal, span...)
end

"""
    load(dataset::Dataset, uuid::UUID, names[, span::AbstractTimeSpan])

Return `Dict(name => load(dataset, uuid, name[, span]) for name in names)`.
"""
function load(dataset::Dataset, uuid::UUID, names, span::AbstractTimeSpan...)
    return Dict(name => load(dataset, uuid, name, span...) for name in names)
end

"""
    load(dataset::Dataset, uuid::UUID[, span::AbstractTimeSpan])

Return `load(dataset, uuid, names[, span])` where `names` is a list of all
signal names in the recording specified by `uuid`.
"""
function load(dataset::Dataset, uuid::UUID, span::AbstractTimeSpan...)
    return load(dataset, uuid, keys(dataset.recordings[uuid].signals), span...)
end

#####
##### `store!`
#####

"""
    store!(dataset::Dataset, uuid::UUID, name::Symbol, samples::Samples;
           overwrite::Bool=true)

Add `name => samples.signal` to `dataset.recordings[uuid].signals` and serialize
`samples.data` to the proper file location within `dataset.path`.

If `overwrite` is `false`, an error is thrown if `samples` already exists
in `recording`/`dataset`. Otherwise, existing entries matching `samples.signal`
will be deleted and replaced with `samples`.
"""
function store!(dataset::Dataset, uuid::UUID, name::Symbol, samples::Samples;
                overwrite::Bool=true)
    recording, signal = dataset.recordings[uuid], samples.signal
    if haskey(recording.signals, name) && !overwrite
        throw(ArgumentError("$name already exists in $uuid and `overwrite` is `false`"))
    end
    is_valid(signal) || throw(ArgumentError("signal in `samples` is invalid"))
    if !is_lower_snake_case_alphanumeric(string(name))
        throw(ArgumentError("$name is not lower snake case and alphanumeric"))
    end
    recording.signals[name] = signal
    store_samples!(samples_path(dataset, uuid, name, signal.file_extension), samples;
                   overwrite=overwrite)
    return recording
end

#####
##### `delete!`
#####

"""
    delete!(dataset::Dataset, uuid::UUID)

Delete the recording whose UUID matches `uuid` from `dataset`. This function
removes the matching `Recording` object from `dataset.recordings`, as well as
deletes the corresponding subdirectory in the `dataset`'s `samples` directory.
"""
function Base.delete!(dataset::Dataset, uuid::UUID)
    delete!(dataset.recordings, uuid)
    rm(samples_path(dataset, uuid); recursive=true, force=true)
    return dataset
end

"""
    delete!(dataset::Dataset, uuid::UUID, name::Symbol)

Delete the signal whose name matches `name` from the recording whose UUID matches
`uuid` in `dataset`. This function removes the matching `Signal` object from
`dataset.recordings[uuid]`, as well as deletes the corresponding sample data in
the `dataset`'s `samples` directory.
"""
function Base.delete!(dataset::Dataset, uuid::UUID, name::Symbol)
    rm(samples_path(dataset, uuid, name); force=true)
    delete!(dataset.recordings[uuid].signals, name)
    return dataset
end


text:
@testset "initialization strategies" begin
    n = 100
    d = 5
    dummy_data, labels = generate_mvn_with_outliers(d, n, 42, true, false)
    pools = fill(:U, size(dummy_data, 2))
    model = SVDD.VanillaSVDD(dummy_data)

    @testset "RuleOfThumbSilverman" begin
        # see https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.stats.gaussian_kde.html
        expected = (n * (d + 2) / 4.0)^(-1.0 / (d + 4))
        @test expected == SVDD.calculate_gamma(model, SVDD.RuleOfThumbSilverman())
    end

    @testset "RuleOfThumbScott" begin
        # see https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy.stats.gaussian_kde.html
        expected = n^(-1.0 / (d + 4))
        @test expected == SVDD.calculate_gamma(model, SVDD.RuleOfThumbScott())
    end

    @testset "TaxErrorEstimate" begin
        # see Tax and Duin 2004, Support Vector Data Description
        outlier_percentage = .1
        expected = 1 / (n*outlier_percentage)
        @test expected == SVDD.calculate_C(model, SVDD.TaxErrorEstimate(outlier_percentage))
    end

    @testset "BoundedTaxErrorEstimate-1" begin
        n = 10
        dummy_data_small, _ = generate_mvn_with_outliers(d, n, 42, true, false)
        lower_bound = 0.03
        upper_bound = 0.97
        outlier_percentage = 0.05
        @test upper_bound ≈ SVDD.calculate_C(SVDD.VanillaSVDD(dummy_data_small), SVDD.BoundedTaxErrorEstimate(outlier_percentage, lower_bound, upper_bound))
    end

    @testset "BoundedTaxErrorEstimate-2" begin
        n = 2000
        dummy_data_small, _ = generate_mvn_with_outliers(d, n, 42, true, false)
        lower_bound = 0.03
        upper_bound = 0.97
        outlier_percentage = 0.05
        @test lower_bound ≈ SVDD.calculate_C(SVDD.VanillaSVDD(dummy_data_small), SVDD.BoundedTaxErrorEstimate(outlier_percentage, lower_bound, upper_bound))
    end

    @testset "SimpleCombinedStrategy" begin
        gamma_expected = SVDD.calculate_gamma(model, SVDD.RuleOfThumbScott())
        C_expected = 0.5
        init_strategy = SVDD.SimpleCombinedStrategy(SVDD.RuleOfThumbScott(), SVDD.FixedCStrategy(C_expected))
        C, kernel = SVDD.get_parameters(model, init_strategy)
        @test C == C_expected
        @test kernel == MLKernels.GaussianKernel(gamma_expected)
    end

    @testset "BinarySearchCStrategy" begin
        model = SVDD.VanillaSVDD(dummy_data)

        target_outlier_percentage = 0.1
        C_init = 0.5
        max_iter = 10
        eps = 0.01

        strategy = SVDD.BinarySearchCStrategy(target_outlier_percentage, C_init, max_iter, eps, TEST_SOLVER)
        gamma = 5.0
        C_actual = SVDD.calculate_C(model, strategy, GaussianKernel(gamma))

        init_strategy = SVDD.FixedParameterInitialization(GaussianKernel(gamma), C_actual)

        SVDD.initialize!(model, init_strategy)
        SVDD.fit!(model, TEST_SOLVER)
        actual_outlier_percentage = countmap(SVDD.classify.(SVDD.predict(model, model.data)))[:outlier] / size(model.data, 2)
        @test abs.(actual_outlier_percentage - target_outlier_percentage) < 0.04
    end

    @testset "GammaFirstCombinedStrategy" begin
        model = SVDD.SVDDneg(dummy_data, pools)
        target_outlier_percentage = 0.15

        gamma_strategy = SVDD.RuleOfThumbSilverman()
        C_strategy = SVDD.BinarySearchCStrategy(target_outlier_percentage, 0.5, 10, 0.01, TEST_SOLVER)
        init_strategy = SVDD.GammaFirstCombinedStrategy(gamma_strategy, C_strategy)

        SVDD.initialize!(model, init_strategy)
        SVDD.fit!(model, TEST_SOLVER)
        actual_outlier_percentage = countmap(SVDD.classify.(SVDD.predict(model, model.data)))[:outlier] / size(model.data, 2)
        @test abs.(actual_outlier_percentage - target_outlier_percentage) < 0.04
    end

    @testset "SimpleSubspaceStrategy" begin
        subspaces = [[1, 2], [3, 4, 5]]
        model = SVDD.SubSVDD(dummy_data, subspaces, pools)
        for gamma_strategy in [RuleOfThumbScott(), RuleOfThumbSilverman()]
            @testset "RuleOfThumb-Global-$(typeof(gamma_strategy))" begin
                    strategy = SVDD.SimpleSubspaceStrategy(
                        gamma_strategy,
                        FixedCStrategy(1.0),
                        gamma_scope = Val(:Global))
                    C, kernel = SVDD.get_parameters(model, strategy)
                    @test C ≈ 1.0
                    @test length(kernel) == 2
                    @test all(map(x -> x <: MLKernels.GaussianKernel, typeof.(kernel)))
                    @test kernel[1] == kernel[2]
                end
            @testset "RuleOfThumb-Local-$(typeof(gamma_strategy))" begin
                strategy = SVDD.SimpleSubspaceStrategy(
                    gamma_strategy,
                    FixedCStrategy(1.0),
                    gamma_scope = Val(:Subspace))
                C, kernel = SVDD.get_parameters(model, strategy)
                @test C ≈ 1.0
                @test length(kernel) == 2
                @test all(map(x -> x <: MLKernels.GaussianKernel, typeof.(kernel)))
                @test kernel[1] != kernel[2]
            end
        end
    end

    dummy_data, labels = generate_mvn_with_outliers(2, 50, 42, true, true)
    pools = fill(:Lin, size(dummy_data, 2))

    @testset "WangGammaStrategy" begin
        model = SVDD.VanillaSVDD(dummy_data[:, labels .== :inlier])
        gamma_strategy = SVDD.WangGammaStrategy(TEST_SOLVER, [0.1, 0.5], 1)
        init_strategy = SVDD.SimpleCombinedStrategy(gamma_strategy, SVDD.FixedCStrategy(1))

        SVDD.initialize!(model, init_strategy)
        SVDD.fit!(model, TEST_SOLVER)
        @test all(SVDD.classify.(SVDD.predict(model, dummy_data)) .== labels)
    end

    @testset "WangCombinedInitializationStrategy" begin
        model = SVDD.VanillaSVDD(dummy_data[:, labels .== :inlier])
        fixed_C = 0.95
        init_strategy = SVDD.WangCombinedInitializationStrategy(TEST_SOLVER, [0.1, 0.5], SVDD.FixedCStrategy(fixed_C))

        SVDD.initialize!(model, init_strategy)
        SVDD.fit!(model, TEST_SOLVER)
        @test get_model_params(model)[:C] == fixed_C
        @test all(SVDD.classify.(SVDD.predict(model, dummy_data)) .== labels)
    end
end


text:
struct Model <: UQModel
    func::Function
    name::Symbol
end

function (obj::Model)(df::DataFrame)
    return obj.func(df)
end

function evaluate!(models::Array{T,1} where {T<:UQModel}, df::DataFrame)
    for m in models
        evaluate!(m, df)
    end
    return nothing
end

function evaluate!(m::Model, df::DataFrame)
    df[!, m.name] = m.func(df)
    return nothing
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
using UnicodePlots, ReferenceTests, Test
using ReferenceTests: BeforeAfterFull
using Dates: Date, Day
import Random: seed!
using ColorTypes
using StableRNGs
using StatsBase

include("fixes.jl")

const RNG = StableRNG(1337)
const T_SZ = (24, 80)  # terminal size

# see JuliaTesting/ReferenceTests.jl/pull/91
test_ref(reference, actual) =
    @test_reference(reference, actual, render = BeforeAfterFull(), format = "TXT")

# helpers
macro show_col(p, kv...)
    :(@io2str(
        $(Expr(
            :call,
            :show,
            Expr(:call, :IOContext, :(::IO), :color => true, kv...),
            esc(p),
        ))
    ))
end
macro show_nocol(p, kv...)
    :(@io2str(
        $(Expr(
            :call,
            :show,
            Expr(:call, :IOContext, :(::IO), :color => false, kv...),
            esc(p),
        ))
    ))
end
macro print_col(p, kv...)
    :(@io2str(
        $(Expr(
            :call,
            :print,
            Expr(:call, :IOContext, :(::IO), :color => true, kv...),
            esc(p),
        ))
    ))
end
macro print_nocol(p, kv...)
    :(@io2str(
        $(Expr(
            :call,
            :print,
            Expr(:call, :IOContext, :(::IO), :color => false, kv...),
            esc(p),
        ))
    ))
end

withenv("FORCE_COLOR" => "X") do  # github.com/JuliaPlots/UnicodePlots.jl/issues/134
    for test in (
        "tst_issues.jl",
        "tst_common.jl",
        "tst_graphics.jl",
        "tst_canvas.jl",
        "tst_plot.jl",
        "tst_barplot.jl",
        "tst_histogram.jl",
        "tst_lineplot.jl",
        "tst_scatterplot.jl",
        "tst_contourplot.jl",
        "tst_spy.jl",
        "tst_boxplot.jl",
        "tst_heatmap.jl",
        "tst_deprecated_warns.jl",
    )
        @testset "$test" begin
            include(test)
        end
    end
end


text:
module InterfaceSymbolicUtilsModule

using SymbolicUtils
import ..CoreModule: CONST_TYPE, Node, Options
import ..UtilsModule: isgood, isbad, @return_on_false

const SYMBOLIC_UTILS_TYPES = Union{<:Number,SymbolicUtils.Symbolic{<:Number}}

const SUPPORTED_OPS = (cos, sin, exp, cot, tan, csc, sec, +, -, *, /)

function isgood(x::SymbolicUtils.Symbolic)
    return if SymbolicUtils.istree(x)
        all(isgood.([SymbolicUtils.operation(x); SymbolicUtils.arguments(x)]))
    else
        true
    end
end
subs_bad(x) = isgood(x) ? x : Inf

function parse_tree_to_eqs(tree::Node, options::Options, index_functions::Bool=false)
    if tree.degree == 0
        # Return constant if needed
        tree.constant && return subs_bad(tree.val)
        return SymbolicUtils.Sym{LiteralReal}(Symbol("x$(tree.feature)"))
    end
    # Collect the next children
    children = tree.degree >= 2 ? (tree.l, tree.r) : (tree.l,)
    # Get the operation
    op = tree.degree > 1 ? options.binops[tree.op] : options.unaops[tree.op]
    # Create an N tuple of Numbers for each argument
    dtypes = map(x -> Number, 1:(tree.degree))
    #
    if !(op ∈ SUPPORTED_OPS) && index_functions
        op = SymbolicUtils.Sym{(SymbolicUtils.FnType){Tuple{dtypes...},Number}}(Symbol(op))
    end

    return subs_bad(
        op(map(x -> parse_tree_to_eqs(x, options, index_functions), children)...)
    )
end

# For operators which are indexed, we need to convert them back
# using the string:
function convert_to_function(
    x::SymbolicUtils.Sym{SymbolicUtils.FnType{T,Number}}, options::Options
) where {T<:Tuple}
    degree = length(T.types)
    if degree == 1
        ind = findoperation(x.name, options.unaops)
        return options.unaops[ind]
    elseif degree == 2
        ind = findoperation(x.name, options.binops)
        return options.binops[ind]
    else
        throw(AssertionError("Function $(String(x.name)) has degree > 2 !"))
    end
end

# For normal operators, simply return the function itself:
convert_to_function(x, options::Options) = x

# Split equation
function split_eq(
    op, args, options::Options; varMap::Union{Array{String,1},Nothing}=nothing
)
    !(op ∈ (sum, prod, +, *)) && throw(error("Unsupported operation $op in expression!"))
    if Symbol(op) == Symbol(sum)
        ind = findoperation(+, options.binops)
    elseif Symbol(op) == Symbol(prod)
        ind = findoperation(*, options.binops)
    else
        ind = findoperation(op, options.binops)
    end
    return Node(
        ind,
        convert(Node, args[1], options; varMap=varMap),
        convert(Node, op(args[2:end]...), options; varMap=varMap),
    )
end

function findoperation(op, ops)
    for (i, oi) in enumerate(ops)
        Symbol(oi) == Symbol(op) && return i
    end
    throw(error("Operation $(op) in expression not found in operations $(ops)!"))
end

function Base.convert(
    ::typeof(SymbolicUtils.Symbolic),
    tree::Node,
    options::Options;
    varMap::Union{Array{String,1},Nothing}=nothing,
    index_functions::Bool=false,
)
    return node_to_symbolic(tree, options; varMap=varMap, index_functions=index_functions)
end

function Base.convert(
    ::typeof(Node),
    x::Number,
    options::Options;
    varMap::Union{Array{String,1},Nothing}=nothing,
)
    return Node(CONST_TYPE(x))
end

function Base.convert(
    ::typeof(Node),
    expr::SymbolicUtils.Symbolic,
    options::Options;
    varMap::Union{Array{String,1},Nothing}=nothing,
)
    if !SymbolicUtils.istree(expr)
        varMap === nothing && return Node(String(expr.name))
        return Node(String(expr.name), varMap)
    end

    # First, we remove integer powers:
    y, good_return = multiply_powers(expr)
    if good_return
        expr = y
    end

    op = convert_to_function(SymbolicUtils.operation(expr), options)
    args = SymbolicUtils.arguments(expr)

    length(args) > 2 && return split_eq(op, args, options; varMap=varMap)
    ind = if length(args) == 2
        findoperation(op, options.binops)
    else
        findoperation(op, options.unaops)
    end

    return Node(ind, map(x -> convert(Node, x, options; varMap=varMap), args)...)
end

"""
    node_to_symbolic(tree::Node, options::Options;
                varMap::Union{Array{String, 1}, Nothing}=nothing,
                index_functions::Bool=false)

The interface to SymbolicUtils.jl. Passing a tree to this function
will generate a symbolic equation in SymbolicUtils.jl format.

## Arguments

- `tree::Node`: The equation to convert.
- `options::Options`: Options, which contains the operators used in the equation.
- `varMap::Union{Array{String, 1}, Nothing}=nothing`: What variable names to use for
    each feature. Default is [x1, x2, x3, ...].
- `index_functions::Bool=false`: Whether to generate special names for the
    operators, which then allows one to convert back to a `Node` format
    using `symbolic_to_node`.
    (CURRENTLY UNAVAILABLE - See https://github.com/MilesCranmer/SymbolicRegression.jl/pull/84).
"""
function node_to_symbolic(
    tree::Node,
    options::Options;
    varMap::Union{Array{String,1},Nothing}=nothing,
    index_functions::Bool=false,
)
    expr = subs_bad(parse_tree_to_eqs(tree, options, index_functions))
    # Check for NaN and Inf
    @assert isgood(expr) "The recovered equation contains NaN or Inf."
    # Return if no varMap is given
    varMap === nothing && return expr
    # Create a substitution tuple
    subs = Dict(
        [
            SymbolicUtils.Sym{LiteralReal}(Symbol("x$(i)")) =>
                SymbolicUtils.Sym{LiteralReal}(Symbol(varMap[i])) for i in 1:length(varMap)
        ]...,
    )
    return substitute(expr, subs)
end

function symbolic_to_node(
    eqn::T, options::Options; varMap::Union{Array{String,1},Nothing}=nothing
)::Node where {T<:SymbolicUtils.Symbolic}
    return convert(Node, eqn, options; varMap=varMap)
end

# function Base.convert(::typeof(Node), x::Number, options::Options; varMap::Union{Array{String, 1}, Nothing}=nothing)
# function Base.convert(::typeof(Node), expr::SymbolicUtils.Symbolic, options::Options; varMap::Union{Array{String, 1}, Nothing}=nothing)

function multiply_powers(eqn::Number)::Tuple{SYMBOLIC_UTILS_TYPES,Bool}
    return eqn, true
end

function multiply_powers(eqn::SymbolicUtils.Symbolic)::Tuple{SYMBOLIC_UTILS_TYPES,Bool}
    if !SymbolicUtils.istree(eqn)
        return eqn, true
    end
    op = SymbolicUtils.operation(eqn)
    return multiply_powers(eqn, op)
end

function multiply_powers(
    eqn::SymbolicUtils.Symbolic, op::F
)::Tuple{SYMBOLIC_UTILS_TYPES,Bool} where {F}
    args = SymbolicUtils.arguments(eqn)
    nargs = length(args)
    if nargs == 1
        l, complete = multiply_powers(args[1])
        @return_on_false complete eqn
        @return_on_false isgood(l) eqn
        return op(l), true
    elseif op == ^
        l, complete = multiply_powers(args[1])
        @return_on_false complete eqn
        @return_on_false isgood(l) eqn
        n = args[2]
        if typeof(n) <: Int
            if n == 1
                return l, true
            elseif n == -1
                return 1.0 / l, true
            elseif n > 1
                return reduce(*, [l for i in 1:n]), true
            elseif n < -1
                return reduce(/, vcat([1], [l for i in 1:abs(n)])), true
            else
                return 1.0, true
            end
        else
            r, complete2 = multiply_powers(args[2])
            @return_on_false complete2 eqn
            return l^r, true
        end
    elseif nargs == 2
        l, complete = multiply_powers(args[1])
        @return_on_false complete eqn
        @return_on_false isgood(l) eqn
        r, complete2 = multiply_powers(args[2])
        @return_on_false complete2 eqn
        @return_on_false isgood(r) eqn
        return op(l, r), true
    else
        # return mapreduce(multiply_powers, op, args)
        # ## reduce(op, map(multiply_powers, args))
        out = map(multiply_powers, args) #vector of tuples
        for i in 1:size(out, 1)
            @return_on_false out[i][2] eqn
            @return_on_false isgood(out[i][1]) eqn
        end
        cumulator = out[1][1]
        for i in 2:size(out, 1)
            cumulator = op(cumulator, out[i][1])
            @return_on_false isgood(cumulator) eqn
        end
        return cumulator, true
    end
end

end