text:
module TestImages
using FileIO, AxisArrays, OffsetArrays
using Pkg.Artifacts
using StringDistances
using ColorTypes
const artifacts_toml = abspath(joinpath(@__DIR__, "..", "Artifacts.toml"))

export testimage

# TODO: export shepp_logan when we remove Images.shepp_logan
# https://github.com/JuliaImages/Images.jl/pull/904
# export shepp_logan

remotefiles = [
    "autumn_leaves.png" ,
    "blobs.gif" ,
    "cameraman.tif" ,
    "earth_apollo17.jpg" ,
    "fabio_color_256.png" ,
    "fabio_color_512.png" ,
    "fabio_gray_256.png" ,
    "fabio_gray_512.png" ,
    "hela-cells.tif" ,
    "house.tif" ,
    "jetplane.tif" ,
    "lake_color.tif" ,
    "lake_gray.tif" ,
    "lena_color_256.tif" ,
    "lena_color_512.tif" ,
    "lena_gray_256.tif" ,
    "lena_gray_512.tif" ,
    "lena_gray_16bit.png" ,
    "livingroom.tif" ,
    "lighthouse.png" ,
    "mandril_color.tif" ,
    "mandril_gray.tif" ,
    "mandrill.tiff" ,
    "m51.tif" ,
    "moonsurface.tiff" ,
    "mountainstream.png" ,
    "mri-stack.tif" ,
    "multi-channel-time-series.ome.tif" ,
    "peppers_color.tif" ,
    "peppers_gray.tif" ,
    "pirate.tif" ,
    "resolution_test_512.tif" ,
    "resolution_test_1920.tif" ,
    "simple_3d_ball" ,
    "simple_3d_psf" ,
    "toucan.png" ,
    "walkbridge.tif" ,
    "woman_blonde.tif" ,
    "woman_darkhair.tif" ,
]

"""
    img = testimage(filename; download_only=false, [ops...])

Load test image that partially matches `filename`, the first match is used if there're more
than one.

If `download_only=true`, the full filepath is returned.
Any other keyword arguments `ops` will be passed to image IO backend through `FileIO.load`.

# Example

```julia
julia> using TestImages
julia> img = testimage("cameraman.tif"); # fullname
julia> img = testimage("cameraman"); # without extension works
julia> img = testimage("c"); # with only partial name also works
```

# Extended help

The following is a complete list of testimages, you can also check them at
https://testimages.juliaimages.org/

$(reduce((x, y)->join([x, "\n - \`\"", splitext(y)[1], "\"\`"]), sort(remotefiles); init=""))
"""
function testimage(filename; download_only::Bool = false, ops...)
    imagefile = image_path(full_imagename(filename))

    if startswith(basename(imagefile), "lena_")
        @warn "Usage of \"lena\" is not recommended, and the image may be removed in a later release. See https://womenlovetech.com/losing-lena-why-we-need-to-remove-one-image-and-end-techs-original-sin/ for more information." maxlog=1
    end

    download_only && return imagefile

    img = load(imagefile; ops...)
    if basename(imagefile) == "mri-stack.tif"
        # orientation is posterior-right-superior,
        # see http://www.grahamwideman.com/gw/brain/orientation/orientterms.htm
        return AxisArray(img, (:P, :R, :S), (1, 1, 5))
    elseif basename(imagefile) == "simple_3d_psf"
        # zero frequency is at (0, 0, 0)
        return OffsetArray(img, (-33, -33, -33))
    end
    img
end

"""
    fullname = full_imagename(shortname)

Get the first match of `shortname` in `TestImages.remotefiles`
"""
function full_imagename(filename)
    idx = findfirst(remotefiles) do x
        startswith(x, filename)
    end
    if idx === nothing
        warn_msg = "\"$filename\" not found in `TestImages.remotefiles`."

        best_match = _findmax(filename)
        if isnothing(best_match[2])
            similar_matches = remotefiles[_findall(filename)]
            if !isempty(similar_matches)
                similar_matches_msg = "  * \"" * join(similar_matches, "\"\n  * \"") * "\""
                warn_msg = "$(warn_msg) Do you mean one of the following?\n$(similar_matches_msg)"
            end
            throw(ArgumentError(warn_msg))
        else
            idx = best_match[2]
            @warn "$(warn_msg) Load \"$(remotefiles[idx])\" instead."
        end
    end
    return remotefiles[idx]
end

function image_path(imagename)
    ensure_artifact_installed("images", artifacts_toml)

    image_dir = artifact_path(artifact_hash("images", artifacts_toml))
    return joinpath(image_dir, imagename)
end

_findall(name; min_score=0.6) = findall(name, remotefiles, Winkler(Jaro()), min_score=min_score)
_findmax(name; min_score=0.8) = findmax(name, remotefiles, Winkler(Jaro()), min_score=min_score)

"""
    phantom = shepp_logan(N, [M=N]; high_contrast=true)

Output the NxM Shepp-Logan phantom, which is a standard test image usually used
for comparing image reconstruction algorithms in the field of computed
tomography (CT) and magnetic resonance imaging (MRI).

If the argument `M` is omitted, the phantom is of size `NxN`. When setting the keyword argument
`high_constrast` to false, the CT version [1] of the phantom is created. Otherwise, the high contrast
MRI version [2] is calculated.

# References

[1] Shepp, Lawrence A., and Benjamin F. Logan. "The Fourier reconstruction of a head section." _IEEE Transactions on nuclear science_ 21.3 (1974): 21-43.

[2] Toft, Peter Aundal. "The Radon transform-theory and implementation." (1996): 201.

[3] Jain, Anil K. Fundamentals of digital image processing. _Prentice-Hall, Inc._, (1989): 439.
"""
function shepp_logan(N::Int, M::Int; high_contrast::Bool=true, highContrast=nothing)
    if !isnothing(highContrast)
        # compatbitity to Images.shepp_logan
        # remove this when we remove Images.shepp_logan
        Base.depwarn("keyword `highContrast` is deprecated, use `high_contrast` instead.", :shepp_logan)
    end
    x = Array(range(-1, stop=1, length=M)')
    y = Array(range(1, stop=-1, length=N))

    # follow the notation in [2]
    A = high_contrast ?
          # high contrast (MRI) version of the phantom -- [2] p.201
          (1.0 , -0.8   , -0.2  , -0.2  , 0.1 , 0.1  ,  0.1  ,  0.1  ,  0.1  ,  0.1  ) :
          # original (CT) version of the phantom -- [1]
          (2.0 , -0.98  , -0.02 , -0.02 , 0.01, 0.01 ,  0.01 ,  0.01 ,  0.01 ,  0.01 )
          # [3] p.439 uses the following setting for the CT version
          # and is used by MATLAB's built-in `phantom` with method `Shepp-Logan`
        # (1.0 , -0.98  , -0.02 , -0.02 , 0.01, 0.01 ,  0.01 ,  0.01 ,  0.01 ,  0.01 )
    x₀ =  (0.0 ,  0.0   ,  0.22 , -0.22 , 0.0 , 0.0  ,  0.0  , -0.08 ,  0.0  ,  0.06 )
    y₀ =  (0.0 , -0.0184,  0.0  ,  0.0  , 0.35, 0.1  , -0.1  , -0.605, -0.605, -0.605)
    a  =  (0.69,  0.6624,  0.11 ,  0.16 , 0.21, 0.046,  0.046,  0.046,  0.023,  0.023)
    b  =  (0.92,  0.874 ,  0.31 ,  0.41 , 0.25, 0.046,  0.046,  0.023,  0.023,  0.046)
    ϕ  =  (0.0 ,  0.0   , -18.0 ,  18.0 , 0.0 , 0.0  ,  0.0  ,  0.0  ,  0.0  ,  0.0  )

    function _ellipse(dx, dy, a, b, sin_ϕ, cos_ϕ)
        tx = cos_ϕ * dx + sin_ϕ * dy
        ty = sin_ϕ * dx - cos_ϕ * dy
        abs2(tx * b) + abs2(ty * a) < (a * b)^2
    end
    function _ellipse(dx, dy, a, b)
        # a faster case when ϕ == 0.0
        abs2(dx * b) + abs2(dy * a) < (a * b)^2
    end

    P = zeros(Gray{Float64}, N, M)
    for l = 1:length(ϕ)
        if ϕ[l] == 0.0
            @. P = gray(P) + A[l] * _ellipse(x - x₀[l], y - y₀[l], a[l], b[l])
        else
            sin_ϕ, cos_ϕ = sincosd(ϕ[l])
            @. P = gray(P) + A[l] * _ellipse(x - x₀[l], y - y₀[l], a[l], b[l], sin_ϕ, cos_ϕ)
        end
    end

    return P
end
shepp_logan(N::Integer, M::Integer; kwargs...) = shepp_logan(Int(N), Int(M); kwargs...)
shepp_logan(N::Integer; kwargs...) = shepp_logan(Int(N), Int(N); kwargs...)

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    @assert precompile(image_path, (String,))
    @assert precompile(testimage, (String,))
    @assert precompile(shepp_logan, (Int,Int))
end
VERSION >= v"1.4.2" && _precompile_() # https://github.com/JuliaLang/julia/pull/35378

end # module


text:
## References
# - [1] Zhan, Zhang, and Chung. Adaptive particle swarm optimization, IEEE Transactions on Systems, Man, and Cybernetics, Part B: CyberneticsVolume 39, Issue 6 (2009): 1362-1381
struct ParticleSwarm{Tn,T}
    n_particles::Tn
    limit_search_space::Bool
    elitist_learning::Bool
    c₁::T
    c₂::T
    σmin::T
    σmax::T
end
# defaults for c's are from the flowchart of page 1370, σ's are from the bottom of 1370
ParticleSwarm(;
    n_particles = nothing,
    limit_search_space = false,
    elitist_learning = true,
    c₁ = 2.0,
    c₂ = 2.0,
    σmin = 0.1,
    σmax = 1.0,
) = ParticleSwarm(n_particles, limit_search_space, elitist_learning, c₁, c₂, σmin, σmax)

summary(::ParticleSwarm) = "Adaptive Particle Swarm"
function solve(
    problem::OptimizationProblem,
    x0,
    method::ParticleSwarm,
    options::OptimizationOptions,
)
    if !(mstyle(problem) === InPlace())
        throw(
            ErrorException("solve() not defined for OutOfPlace() with $(summary(method))"),
        )
    end
    t0 = time()
    n_particles = method.n_particles isa Nothing ? max(length(x0), 5) : method.n_particles
    lower, upper = lowerbounds(problem), upperbounds(problem)

    T, n = eltype(x0), length(x0)

    c₁, c₂, σmin, σmax = T(method.c₁), T(method.c₂), T(method.σmin), T(method.σmax)
    ω = T(9) / 10 # see page 1368

    X, V = [copy(x0) for i = 1:n_particles], [x0 * T(0) for i = 1:n_particles]
    X_best = [x0 * T(0) for i = 1:n_particles]

    Fs, Fs_best = zeros(T, n_particles), zeros(T, n_particles)

    x = copy(x0)
    if method.elitist_learning
        x_learn = copy(x0)
    end
    current_state = 0

    # spread the initial population uniformly over the whole search space
    width = upper .- lower
    for i = 1:n_particles
        X[i] .= lower .+ width .* rand(T)
        X_best[i] .= X[i]
    end

    best_f = T(0)
    swarm_f = best_f
    f0 = swarm_f
    X_best[1] .= x0
    X[1] .= x0
    iter = 0
    while iter <= options.maxiter
        iter += 1
        limit_X!(X, lower, upper)
        Fs = batched_value(problem, Fs, X)

        if iter == 1
            copyto!(Fs_best, Fs)
            best_f = Base.minimum(Fs)
        end
        best_f = housekeeping!(Fs, Fs_best, X, X_best, x, best_f)
        if iter == 1
            f0 = best_f
        end
        if method.elitist_learning # try to avoid non-global minima
            x_learn .= x
            # Perturb a random dimension and replace the current worst
            # solution in X with x_learn if x_learn presents the new
            # best solution. Else, discard x_learn.
            worst_Fs, i_worst = findmax(Fs)
            random_index = rand(1:n)

            # Eqn (14) - See Section VI-C for a discussion of max and min
            # Could allow for an annealing scheme. Maybe a SAMIN style?
            σlearn = σmax - (σmax - σmin) * iter / options.maxiter

            x_learn[random_index] =
                x_learn[random_index] + width[random_index] * σlearn * randn()
            x_learn[random_index] =
                max(lower[random_index], min(upper[random_index], x_learn[random_index]))

            Fs_learn = value(problem, x_learn)
            if Fs_learn < best_f
                X_best[i_worst] .= x_learn
                X[i_worst] .= x_learn
                x .= x_learn

                Fs_best[i_worst] = Fs_learn
                Fs[i_worst] = Fs_learn
                best_f = Fs_learn
            end
        end
        # TODO find a better name for _f (look inthe paper, it might be called f there)
        current_state, swarm_f = get_swarm_state(X, Fs, x, current_state)
        ω, c₁, c₂ = update_swarm_params!(c₁, c₂, ω, current_state, swarm_f)
        update_swarm!(X, X_best, x, n, V, ω, c₁, c₂)
    end
    best_f, x
    ConvergenceInfo(
        method,
        (
            swarm_f = swarm_f,
            X = X,
            Fs = Fs,
            solution = x,
            minimum = best_f,
            f0 = f0,
            iter = iter,
            time = time() - t0,
        ),
        options,
    )
end

function update_swarm!(X, X_best, best_point, n, V, ω, c₁, c₂)

    Tx = eltype(first(X))
    # compute new positions for the swarm particles
    # FIXME vmax should be 20%
    for i in eachindex(X, X_best)
        for j = 1:n
            r1 = rand(Tx)
            r2 = rand(Tx)
            vx = X_best[i][j] - X[i][j]
            vg = best_point[j] - X[i][j]
            V[i][j] = V[i][j] * ω + c₁ * r1 * vx + c₂ * r2 * vg
            X[i][j] = X[i][j] + V[i][j]
        end
    end
end

function get_mu_1(f::Tx) where {Tx}
    if Tx(0) <= f <= Tx(4) / 10
        return Tx(0)
    elseif Tx(4) / 10 < f <= Tx(6) / 10
        return Tx(5) * f - Tx(2)
    elseif Tx(6) / 10 < f <= Tx(7) / 10
        return Tx(1)
    elseif Tx(7) / 10 < f <= Tx(8) / 10
        return -Tx(10) * f + Tx(8)
    else
        return Tx(0)
    end
end

function get_mu_2(f::Tx) where {Tx}
    if Tx(0) <= f <= Tx(2) / 10
        return Tx(0)
    elseif Tx(2) / 10 < f <= Tx(3) / 10
        return Tx(10) * f - Tx(2)
    elseif Tx(3) / 10 < f <= Tx(4) / 10
        return Tx(1)
    elseif Tx(4) / 10 < f <= Tx(6) / 10
        return -Tx(5) * f + Tx(3)
    else
        return Tx(0)
    end
end

function get_mu_3(f::Tx) where {Tx}
    if Tx(0) <= f <= Tx(1) / 10
        return Tx(1)
    elseif Tx(1) / 10 < f <= Tx(3) / 10
        return -Tx(5) * f + Tx(3) / 2
    else
        return Tx(0)
    end
end

function get_mu_4(f::Tx) where {Tx}
    if Tx(0) <= f <= Tx(7) / 10
        return Tx(0)
    elseif Tx(7) / 10 < f <= Tx(9) / 10
        return Tx(5) * f - Tx(7) / 2
    else
        return Tx(1)
    end
end

function get_swarm_state(X, Fs, best_point, previous_state)
    # swarm can be in 4 different states, depending on which
    # the weighing factors c₁ and c₂ are adapted.
    # New state is not only depending on the current swarm state,
    # but also from the previous
    n_particles = length(X)
    n = length(first(X))
    Tx = eltype(first(X))
    f_best, i_best = findmin(Fs)
    d = zeros(Tx, n_particles)
    for i = 1:n_particles
        dd = Tx(0)
        for k = 1:n_particles
            for dim = 1:n
                @inbounds ddd = (X[i][dim] - X[k][dim])
                dd += ddd * ddd
            end
        end
        d[i] = sqrt(dd)
    end
    dg = d[i_best]
    dmin = Base.minimum(d)
    dmax = Base.maximum(d)

    f = (dg - dmin) / max(dmax - dmin, sqrt(eps(Tx)))

    mu = zeros(Tx, 4)
    mu[1] = get_mu_1(f)
    mu[2] = get_mu_2(f)
    mu[3] = get_mu_3(f)
    mu[4] = get_mu_4(f)
    best_mu, i_best_mu = findmax(mu)
    current_state = 0

    if previous_state == 0
        current_state = i_best_mu
    elseif previous_state == 1
        if mu[1] > 0
            current_state = 1
        else
            if mu[2] > 0
                current_state = 2
            elseif mu[4] > 0
                current_state = 4
            else
                current_state = 3
            end
        end
    elseif previous_state == 2
        if mu[2] > 0
            current_state = 2
        else
            if mu[3] > 0
                current_state = 3
            elseif mu[1] > 0
                current_state = 1
            else
                current_state = 4
            end
        end
    elseif previous_state == 3
        if mu[3] > 0
            current_state = 3
        else
            if mu[4] > 0
                current_state = 4
            elseif mu[2] > 0
                current_state = 2
            else
                current_state = 1
            end
        end
    elseif previous_state == 4
        if mu[4] > 0
            current_state = 4
        else
            if mu[1] > 0
                current_state = 1
            elseif mu[2] > 0
                current_state = 2
            else
                current_state = 3
            end
        end
    end
    return current_state, f
end

function update_swarm_params!(c₁, c₂, ω, current_state, f::T) where {T}

    Δc₁ = T(5) / 100 + rand(T) / T(20)
    Δc₂ = T(5) / 100 + rand(T) / T(20)

    if current_state == 1
        c₁ += Δc₁
        c₂ -= Δc₂
    elseif current_state == 2
        c₁ += Δc₁ / 2
        c₂ -= Δc₂ / 2
    elseif current_state == 3
        c₁ += Δc₁ / 2
        c₂ += Δc₂ / 2
    elseif current_state == 4
        c₁ -= Δc₁
        c₂ -= Δc₂
    end

    if c₁ < T(3) / 2
        c₁ = T(3) / 2
    elseif c₁ > T(5) / 2
        c₁ = T(5) / 2
    end

    if c₂ < T(3) / 2
        c₂ = T(5) / 2
    elseif c₂ > T(5) / 2
        c₂ = T(5) / 2
    end

    if c₁ + c₂ > T(4)
        c_total = c₁ + c₂
        c₁ = c₁ / c_total * 4
        c₂ = c₂ / c_total * 4
    end

    ω = 1 / (1 + T(3) / 2 * exp(-T(26) / 10 * f))
    return ω, c₁, c₂
end

function housekeeping!(Fs, Fs_best, X, X_best, best_point, F)
    for i in eachindex(X, X_best)
        if Fs[i] <= Fs_best[i]
            Fs_best[i] = Fs[i]
            X_best[i] .= X[i]

            if Fs[i] <= F
                best_point .= X[i]
                F = Fs[i]
            end
        end
    end
    return F
end

function limit_X!(X, lower, upper)
    # limit X values to boundaries
    for x in X
        x .= min.(max.(x, lower), upper)
    end
    X
end

ext: jl

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
function one_away(str1 :: ASCIIString, str2 :: ASCIIString)
  """
  Check if the shorter string is only one
  edit distance form the longer string
  """
  function is_one_away(shorter :: ASCIIString, longer :: ASCIIString)
    idx_s, idx_l = 1, 1
    len_s, len_l = length(shorter), length(longer)
    find_diff = false
    while idx_s <= len_s && idx_l <= len_l
      if shorter[idx_s] == longer[idx_l]
        idx_s += 1
        idx_l += 1
      else
        if find_diff return false; else find_diff = true; end
        if len_s == len_l idx_s += 1; end
        idx_l += 1
      end
    end
    return true
  end

  if abs(length(str1) - length(str2)) > 1
    return false
  else
    if length(str1) < length(str2)        # pass two strings in the right order
      return is_one_away(str1, str2)
    else
      return is_one_away(str2, str1)
    end
  end
end

# Time  : O(n)
# Space : O(1)

using Base: Test

@test one_away("abc", "abl") == true
@test one_away("ab", "abl") == true
@test one_away("abc", "ab") == true
@test one_away("abc", "abc") == true
@test one_away("abc", "bca") == false

max_stars_repo_name: Kiandr/crackingcodinginterview

max_issues_repo_name: Kiandr/crackingcodinginterview

max_forks_repo_name: Kiandr/crackingcodinginterview


text:
"""
SmetsWouters{T} <: AbstractModel{T}

The SmetsWouters type defines the structure of the FRBNY DSGE
model. We can then concisely pass around a Model object to the
remaining steps of the model (solve, estimate, and forecast).

### Fields

#### Parameters and Steady-States
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model parameters.

* `steady_state::Vector`: Model steady-state values, computed as a
function of elements of `parameters`.

* `keys::Dict{Symbol,Int}`: Maps human-readable names for all model
parameters and steady-states to their indices in `parameters` and
`steady_state`.

#### Inputs to Measurement and Equilibrium Condition Equations

The following fields are dictionaries that map human-readible names to
row and column indices in the matrix representations of of the
measurement equation and equilibrium conditions.

* `endogenous_states::Dict{Symbol,Int}`: Maps each state to a column
in the measurement and equilibrium condition matrices.

* `exogenous_shocks::Dict{Symbol,Int}`: Maps each shock to a column in
the measurement and equilibrium condition matrices.

* `expected_shocks::Dict{Symbol,Int}`: Maps each expected shock to a
column in the measurement and equilibrium condition matrices.

* `equilibrium_conditions::Dict{Symbol,Int}`: Maps each equlibrium
condition to a row in the model's equilibrium condition matrices.

* `endogenous_states_augmented::Dict{Symbol,Int}`: Maps lagged states
to their columns in the measurement and equilibrium condition
equations. These are added after Gensys solves the model.

* `observables::Dict{Symbol,Int}`: Maps each observable to a row in
the model's measurement equation matrices.

#### Model Specifications and Settings

* `spec::AbstractString`: The model specification identifier, "m990",
cached here for filepath computation.

* `subspec::AbstractString`: The model subspecification number,
indicating that some parameters from the original model spec ("ss0")
are initialized differently. Cached here for filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect
computation without changing the economic or mathematical setup of
the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. By default, it is
seeded to ensure replicability in algorithms that involve randomness
(such as Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If
`true`, settings from `m.test_settings` are used in place of those in
`m.settings`.

"""
type SmetsWouters{T} <: AbstractModel{T}
    parameters::ParameterVector{T}                  # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                # model steady-state values
    keys::Dict{Symbol,Int}                          # human-readable names for all the model
                                                    # parameters and steady-n_states

    endogenous_states::Dict{Symbol,Int}             # these fields used to create matrices in the
    exogenous_shocks::Dict{Symbol,Int}              # measurement and equilibrium condition equations.
    expected_shocks::Dict{Symbol,Int}               #
    equilibrium_conditions::Dict{Symbol,Int}        #
    endogenous_states_augmented::Dict{Symbol,Int}   #
    observables::Dict{Symbol,Int}                   #

    spec::ASCIIString                               # Model specification number
    subspec::ASCIIString                            # Model subspecification
    settings::Dict{Symbol,Setting}                  # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}             # Settings/flags for testing mode
    rng::MersenneTwister                            # Random number generator
    testing::Bool                                   # Whether we are in testing mode or not
end

description(m::SmetsWouters) = "Smets-Wouters Model"

"""
`init_model_indices!(m::SmetsWouters)`

Arguments:
`m:: SmetsWouters`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::SmetsWouters)
    # Endogenous states
    endogenous_states = [[
        :y_t, :c_t, :i_t, :qk_t, :k_t, :kbar_t, :u_t, :rk_t, :mc_t,
        :π_t, :μ_ω_t, :w_t, :L_t, :R_t, :g_t, :b_t, :μ_t, :z_t,
        :λ_f_t, :λ_f_t1, :λ_w_t, :λ_w_t1, :rm_t, :Ec_t, :Eqk_t, :Ei_t,
        :Eπ_t, :EL_t, :Erk_t, :Ew_t, :y_f_t, :c_f_t,
        :i_f_t, :qk_f_t, :k_f_t, :kbar_f_t, :u_f_t, :rk_f_t, :w_f_t,
        :L_f_t, :r_f_t, :Ec_f_t, :Eqk_f_t, :Ei_f_t, :EL_f_t, :Erk_f_t, :ztil_t];
        [symbol("rm_tl$i") for i = 1:n_anticipated_shocks(m)]]

    # Exogenous shocks
    exogenous_shocks = [[
        :g_sh, :b_sh, :μ_sh, :z_sh, :λ_f_sh, :λ_w_sh, :rm_sh];
        [symbol("rm_shl$i") for i = 1:n_anticipated_shocks(m)]]

    # Expectations shocks
    expected_shocks = [
        :Ec_sh, :Eqk_sh, :Ei_sh, :Eπ_sh, :EL_sh, :Erk_sh, :Ew_sh, :Ec_f_sh,
        :Eqk_f_sh, :Ei_f_sh, :EL_f_sh, :Erk_f_sh]

    # Equilibrium conditions
    equilibrium_conditions = [[
        :eq_euler, :eq_inv, :eq_capval, :eq_output, :eq_caputl, :eq_capsrv, :eq_capev,
        :eq_mkupp, :eq_phlps, :eq_caprnt, :eq_msub, :eq_wage, :eq_mp, :eq_res, :eq_g, :eq_b, :eq_μ, :eq_z,
        :eq_λ_f, :eq_λ_w, :eq_rm, :eq_λ_f1, :eq_λ_w1, :eq_Ec,
        :eq_Eqk, :eq_Ei, :eq_Eπ, :eq_EL, :eq_Erk, :eq_Ew, :eq_euler_f, :eq_inv_f,
        :eq_capval_f, :eq_output_f, :eq_caputl_f, :eq_capsrv_f, :eq_capev_f, :eq_mkupp_f, :eq_caprnt_f, :eq_msub_f,
        :eq_res_f, :eq_Ec_f, :eq_Eqk_f, :eq_Ei_f, :eq_EL_f, :eq_Erk_f, :eq_ztil];
        [symbol("eq_rml$i") for i=1:n_anticipated_shocks(m)]]

    # Additional states added after solving model
    # Lagged states and observables measurement error
    endogenous_states_augmented = [
        :y_t1, :c_t1, :i_t1, :w_t1, :π_t1, :L_t1, :Et_π_t]

    # Measurement equation observables
    observables = [[
        :obs_gdp,              # quarterly output growth
        :obs_hours,            # aggregate hours growth
        :obs_wages,            # real wage growth
        :obs_gdpdeflator,      # inflation (GDP deflator)
        :obs_nominalrate,      # nominal interest rate
        :obs_consumption,      # consumption growth
        :obs_investment];       # investment growth
        [symbol("obs_nominalrate$i") for i=1:n_anticipated_shocks(m)]] # compounded nominal rates

    for (i,k) in enumerate(endogenous_states);            m.endogenous_states[k]            = i end
    for (i,k) in enumerate(exogenous_shocks);             m.exogenous_shocks[k]             = i end
    for (i,k) in enumerate(expected_shocks);              m.expected_shocks[k]              = i end
    for (i,k) in enumerate(equilibrium_conditions);       m.equilibrium_conditions[k]       = i end
    for (i,k) in enumerate(endogenous_states);            m.endogenous_states[k]            = i end
    for (i,k) in enumerate(endogenous_states_augmented); m.endogenous_states_augmented[k] = i+length(endogenous_states) end
    for (i,k) in enumerate(observables);                  m.observables[k]                  = i end
end

function SmetsWouters(subspec::AbstractString="ss0")

    # Model-specific specifications
    spec               = split(basename(@__FILE__),'.')[1]
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister()        # Random Number Generator
    testing            = false

    # initialize empty model
    m = SmetsWouters{Float64}(
            # model parameters and steady state values
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(), Dict{Symbol,Int}(),

            # model indices
            Dict{Symbol,Int}(), Dict{Symbol,Int}(), Dict{Symbol,Int}(), Dict{Symbol,Int}(), Dict{Symbol,Int}(), Dict{Symbol,Int}(),

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing)

    # Set settings
    settings_smets_wouters!(m)
    default_test_settings!(m)

    # Initialize parameters
    m <= parameter(:α,      0.24, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     Normal(0.30, 0.05),         fixed=false,
                   description="α: Capital elasticity in the intermediate goods sector's Cobb-Douglas production function.",
                   tex_label="\\alpha")

    m <= parameter(:ζ_p,   0.7813, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     BetaAlt(0.5, 0.1),          fixed=false,
                   description="ζ_p: The Calvo parameter. In every period, intermediate goods producers optimize prices with probability (1-ζ_p). With probability ζ_p, prices are adjusted according to a weighted average of the previous period's inflation (π_t1) and steady-state inflation (π_star).",
                   tex_label="\\zeta_p")

    m <= parameter(:ι_p,   0.3291, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     BetaAlt(0.5, 0.15),         fixed=false,
                   description="ι_p: The weight attributed to last period's inflation in price indexation. (1-ι_p) is the weight attributed to steady-state inflation.",

                   tex_label="\\iota_p")

    m <= parameter(:δ,      0.025,  fixed=true,
                   description="δ: The capital depreciation rate.", tex_label="\\delta" )

    m <= parameter(:Upsilon,  1.000,  (0., 10.),     (1e-5, 0.),      DSGE.Exponential(),    GammaAlt(1., 0.5),          fixed=true,
                   description="Υ: The trend evolution of the price of investment goods relative to consumption goods. Set equal to 1.",
                   tex_label="\\mathcal{\\Upsilon}")

    m <= parameter(:Φ,   1.4672, (1., 10.),     (1.00, 10.00),   DSGE.Exponential(),    Normal(1.25, 0.12),         fixed=false,
                   description="Φ: Fixed costs.",
                   tex_label="\\Phi")

    m <= parameter(:S′′,       6.3325, (-15., 15.),   (-15., 15.),     Untransformed(),  Normal(4., 1.5),            fixed=false,
                   description="S'': The second derivative of households' cost of adjusting investment.", tex_label="S\\prime\\prime")

    m <= parameter(:h,        0.7205, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     BetaAlt(0.7, 0.1),          fixed=false,
                   description="h: Consumption habit persistence.", tex_label="h")

    m <= parameter(:ppsi,     0.2648, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     BetaAlt(0.5, 0.15),         fixed=false,
                   description="ppsi: Utilization costs.", tex_label="ppsi")

    m <= parameter(:ν_l,     2.8401, (1e-5, 10.),   (1e-5, 10.),     DSGE.Exponential(),    Normal(2, 0.75),            fixed=false,
                   description="ν_l: The coefficient of relative risk aversion on the labor term of households' utility function.", tex_label="\\nu_l")

    m <= parameter(:ζ_w,   0.7937, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     BetaAlt(0.5, 0.1),          fixed=false,
                   description="ζ_w: (1-ζ_w) is the probability with which households can freely choose wages in each period. With probability ζ_w, wages increase at a geometrically weighted average of the steady state rate of wage increases and last period's productivity times last period's inflation.",
                   tex_label="\\zeta_w")

    m <= parameter(:ι_w,   0.4425, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     BetaAlt(0.5, 0.15),         fixed=false,
                   description="ι_w: No description available.",
                   tex_label="\\iota_w")

    m <= parameter(:λ_w,      1.5000,                                                                               fixed=true,
                   description="λ_w: The wage markup, which affects the elasticity of substitution between differentiated labor services.",
                   tex_label="\\lambda_w")

    m <= parameter(:β, 0.7420, (1e-5, 10.),   (1e-5, 10.),     DSGE.Exponential(),    GammaAlt(0.25, 0.1),        fixed=false,  scaling = x -> 1/(1 + x/100),
                   description="β: Discount rate.",
                   tex_label="\\beta ")

    m <= parameter(:ψ1,  1.7985, (1e-5, 10.),   (1e-5, 10.00),   DSGE.Exponential(),    Normal(1.5, 0.25),          fixed=false,
                   description="ψ₁: Weight on inflation gap in monetary policy rule.",
                   tex_label="\\psi_1")

    m <= parameter(:ψ2,  0.0893, (-0.5, 0.5),   (-0.5, 0.5),     Untransformed(),  Normal(0.12, 0.05),         fixed=false,
                   description="ψ₂: Weight on output gap in monetary policy rule.",
                   tex_label="\\psi_2")

    m <= parameter(:ψ3, 0.2239, (-0.5, 0.5),   (-0.5, 0.5),     Untransformed(),  Normal(0.12, 0.05),         fixed=false,
                   description="ψ₃: Weight on rate of change of output gap in the monetary policy rule.",
                   tex_label="\\psi_3")

    m <= parameter(:π_star,   0.7000, (1e-5, 10.),   (1e-5, 10.),     DSGE.Exponential(),    GammaAlt(0.62, 0.1),        fixed=false,  scaling = x -> 1 + x/100,
                   description="π_star: The steady-state rate of inflation.",
                   tex_label="\\pi_*")

    m <= parameter(:σ_c, 1.2312, (1e-5, 10.),   (1e-5, 10.),     DSGE.Exponential(),    Normal(1.5, 0.37),          fixed=false,
                   description="σ_c: No description available.",
                   tex_label="\\sigma_{c}")

    m <= parameter(:ρ,      .8258, (1e-5, 0.999), (1e-5, 0.999),   SquareRoot(),     BetaAlt(0.75, 0.10),        fixed=false,
                   description="ρ: The degree of inertia in the monetary policy rule.",
                   tex_label="\\rho")

    m <= parameter(:ϵ_p,     10.000,                                                                               fixed=true,
                   description="ϵ_p: No description available.",
                   tex_label="\\varepsilon_{p}")

    m <= parameter(:ϵ_w,     10.000,                                                                               fixed=true,
                   description="ϵ_w: No description available.",
                   tex_label="\\varepsilon_{w}")

    # exogenous processes - level
    m <= parameter(:γ,      0.3982, (-5.0, 5.0),     (-5., 5.),     Untransformed(), Normal(0.4, 0.1),            fixed=false, scaling = x -> x/100,
                   description="γ: The log of the steady-state growth rate of technology.",
                   tex_label="\\gamma")

    m <= parameter(:Lmean,  875., (-1000., 1000.), (-1e3, 1e3),   Untransformed(), Normal(-45, 5),   fixed=false,
                   description="Lmean: No description available.",
                   tex_label="Lmean")

    m <= parameter(:g_star,    0.1800,                                                                               fixed=true,
                   description="g_star: No description available.",
                   tex_label="g_*")

    # exogenous processes - autocorrelation
    m <= parameter(:ρ_g,      0.9930, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_g: AR(1) coefficient in the government spending process.",
                   tex_label="\\rho_g")

    m <= parameter(:ρ_b,      0.2703, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_b: AR(1) coefficient in the intertemporal preference shifter process.",
                   tex_label="\\rho_b")

    m <= parameter(:ρ_μ,     0.5724, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_μ: AR(1) coefficient in capital adjustment cost process.",
                   tex_label="\\rho_{\\mu}")

    m <= parameter(:ρ_z,      0.9676, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")

    m <= parameter(:ρ_λ_f,    0.8692, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_λ_f: AR(1) coefficient in the price mark-up shock process.",
                   tex_label="\\rho_{\\lambda_f}")

    m <= parameter(:ρ_λ_w,    0.9546, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_λ_w: AR(1) coefficient in the wage mark-up shock process.", # CHECK THIS
                   tex_label="\\rho_{\\lambda_f}")

    #monetary policy shock - see eqcond
    m <= parameter(:ρ_rm,     0.3000, (1e-5, 0.999),   (1e-5, 0.999), SquareRoot(),    BetaAlt(0.5, 0.2),           fixed=false,
                   description="ρ_rm: AR(1) coefficient in the monetary policy shock process.", # CHECK THIS
                   tex_label="\\rho_{rm}")

    # exogenous processes - standard deviation
    m <= parameter(:σ_g,      0.6090, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_g: The standard deviation of the government spending process.",
                   tex_label="\\sigma_{g}")

    m <= parameter(:σ_b,      0.1818, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_b: No description available.",
                   tex_label="\\sigma_{b}")

    m <= parameter(:σ_μ,     0.4601, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_μ: The standard deviation of the exogenous marginal efficiency of investment shock process.",
                   tex_label="\\sigma_{\\mu}")

    m <= parameter(:σ_z,      0.4618, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_z: No description available.",
                   tex_label="\\sigma_{z}")

    m <= parameter(:σ_λ_f,    0.1455, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_λ_f: The mean of the process that generates the price elasticity of the composite good.  Specifically, the elasticity is (1+λ_{f,t})/(λ_{f_t}).",
                   tex_label="\\sigma_{\\lambda_f}")

    m <= parameter(:σ_λ_w,    0.2089, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_λ_w: No description available.",
                   tex_label="\\sigma_{\\lambda_w}")

    m <= parameter(:σ_rm,     0.2397, (1e-8, 5.),      (1e-8, 5.),    DSGE.Exponential(),   RootInverseGamma(2., 0.10),  fixed=false,
                   description="σ_rm: No description available.",
                   tex_label="\\sigma_{rm}")

    m <= parameter(:η_gz,       0.0500, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.50, 0.20), fixed=false,
                   description="η_gz: Correlate g and z shocks.",
                   tex_label="\\eta_{gz}")

    m <= parameter(:η_λ_f,      0.7652, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.50, 0.20),         fixed=false,
                   description="η_λ_f: No description available.",
                   tex_label="\\eta_{\\lambda_f}")

    m <= parameter(:η_λ_w,      0.8936, (1e-5, 0.999), (1e-5, 0.999), SquareRoot(),
                   BetaAlt(0.50, 0.20),         fixed=false,
                   description="η_λ_w: AR(2) coefficient on wage markup shock process.",
                   tex_label="\\eta_{\\lambda_w}")

    # steady states
    m <= SteadyStateParameter(:zstar,  NaN, description="steady-state growth rate of productivity", tex_label="\\z_*")
    m <= SteadyStateParameter(:rstar,   NaN, description="steady-state something something", tex_label="\\r_*")
    m <= SteadyStateParameter(:Rstarn,  NaN, description="steady-state something something", tex_label="\\R_*_n")
    m <= SteadyStateParameter(:rkstar,  NaN, description="steady-state something something", tex_label="\\BLAH")
    m <= SteadyStateParameter(:wstar,   NaN, description="steady-state something something", tex_label="\\w_*")
    m <= SteadyStateParameter(:Lstar,   NaN, description="steady-state something something", tex_label="\\L_*")
    m <= SteadyStateParameter(:kstar,   NaN, description="Effective capital that households rent to firms in the steady state.", tex_label="\\k_*")
    m <= SteadyStateParameter(:kbarstar, NaN, description="Total capital owned by households in the steady state.", tex_label="\\bar{k}_*")
    m <= SteadyStateParameter(:istar,  NaN, description="Detrended steady-state investment", tex_label="\\i_*")
    m <= SteadyStateParameter(:ystar,  NaN, description="steady-state something something", tex_label="\\y_*")
    m <= SteadyStateParameter(:cstar,  NaN, description="steady-state something something", tex_label="\\c_*")
    m <= SteadyStateParameter(:wl_c,   NaN, description="steady-state something something", tex_label="\\wl_c")

    init_model_indices!(m)
    init_subspec!(m)
    steadystate!(m)
    return m
end

"""
```
steadystate!(m::SmetsWouters)
```

Calculates the model's steady-state values. `steadystate!(m)` must be called whenever the parameters of `m` are updated.
"""
function steadystate!(m::SmetsWouters)
    m[:zstar]    = log(1+m[:γ]) + m[:α]/(1-m[:α])*log(m[:Upsilon])
    m[:rstar]    = exp(m[:σ_c]*m[:zstar]) / m[:β]
    m[:Rstarn]   = 100*(m[:rstar]*m[:π_star] - 1)
    m[:rkstar]   = m[:rstar]*m[:Upsilon] - (1-m[:δ])
    m[:wstar]    = (m[:α]^m[:α] * (1-m[:α])^(1-m[:α]) * m[:rkstar]^(-m[:α]) / m[:Φ])^(1/(1-m[:α]))
    m[:Lstar]    = 1.
    m[:kstar]    = (m[:α]/(1-m[:α])) * m[:wstar] * m[:Lstar] / m[:rkstar]
    m[:kbarstar] = m[:kstar] * (1+m[:γ]) * m[:Upsilon]^(1 / (1-m[:α]))
    m[:istar]    = m[:kbarstar] * (1-((1-m[:δ])/((1+m[:γ]) * m[:Upsilon]^(1/(1-m[:α])))))
    m[:ystar]    = (m[:kstar]^m[:α]) * (m[:Lstar]^(1-m[:α])) / m[:Φ]
    m[:cstar]    = (1-m[:g_star])*m[:ystar] - m[:istar]
    m[:wl_c]     = (m[:wstar]*m[:Lstar])/(m[:cstar]*m[:λ_w])

    return m
end

function settings_smets_wouters!(m::SmetsWouters)

    default_settings!(m)

    # Anticipated shocks
    m <= Setting(:n_anticipated_shocks, 0)
    m <= Setting(:n_anticipated_shocks_padding, 20)
    m <= Setting(:n_anticipated_lags, 26)

    # Estimation
    m <= Setting(:reoptimize, true)
    m <= Setting(:recalculate_hessian, true)

    # Data vintage
    m <= Setting(:data_vintage, "150827")

    m.settings
end


text:
@testset "SizedString" begin
    @testset "Constructors" begin
        @test SizedString{0}() == SizedString{0}("")
        @test SizedString{4}(NTuple{4, UInt8}(('w', 'h', 'o', 'm'))) == SizedString{4}("whom")
        @test SizedString{5}('B', 'u', 'z', 'z', '!') == SizedString{5}("Buzz!")
        @test SizedString{6}(UInt8['u', 'v', 'w', 'x', 'y', 'z']) == SizedString{6}("uvwxyz")
        @test SizedString{7}(b"abcdefg") == SizedString{7}("abcdefg")

        s = SizedString{3}("foo")
        @test SizedString{3}(s) === s

        @test_throws ArgumentError SizedString{3}("ab")
        @test_throws ArgumentError SizedString{4}(NTuple{5, UInt8}(('!', '(', '/', '%', ')')))
        @test_throws ArgumentError SizedString{5}('B', 'u', 'z', 'z')
        @test_throws ArgumentError SizedString{6}(UInt8['t', 'u', 'v', 'w', 'x', 'y', 'z'])
        @test_throws ArgumentError SizedString{7}(b"abcdef")

        @test_throws InexactError SizedString{3}("Σi")

        @test sized"bar" == SizedString{3}("bar")
    end

    @testset "AbstractString Interface" begin
        N = rand(3:9)
        s = rand(Char(32):Char(126), N)

        ss = SizedString{N}(s)

        @test ncodeunits(ss) == N

        @test codeunit(ss) == UInt8

        @test isvalid(ss, rand(eachindex(s))) == true
        @test isvalid(ss, N+1) == false

        @test sizeof(ss) == N

        @test firstindex(ss) == 1
        @test lastindex(ss) == N

        s0 = SizedString{0}("")
        @test isempty(s0)
        @test !isempty(ss)
    end

    @testset "Iteration" begin
        N = 9
        s = rand(Char(32):Char(126), N)

        ss = SizedString{N}(s)

        for c in ss
            @test c isa Char
            @test c == popfirst!(s)
        end

        @test length(s) == 0

        @test Base.IteratorSize(typeof(ss)) == Base.HasLength()
        @test length(ss) == length(typeof(ss)) == N

        @test Base.IteratorEltype(typeof(ss)) == Base.HasEltype()
        @test eltype(ss) == eltype(typeof(ss)) == Char
    end

    @testset "Indexing" begin
        N = 7
        s = String(rand(Char(32):Char(126), N))

        ss = SizedString{N}(s)

        # Not sure how to test that inbounds actually works
        inbound_codeunit(i) = (local s = ss; @inbounds codeunit(s, i))
        @test_broken inbound_codeunit(0) isa UInt8

        inbound_getindex(i) = (local s = ss; @inbounds s[i])
        @test_broken inbound_getindex(0) isa Char

        for i in eachindex(s)
            @test codeunit(ss, i) === codeunit(s, i)
            @test inbound_codeunit(i) === codeunit(s, i)

            @test ss[i] === s[i]
            @test inbound_getindex(i) === s[i]

            @test ss[1:i] == s[1:i]
            @test ss[collect(i:N)] == s[collect(i:N)]
        end
    end

    let e = SizedString{0}("")
        @test e == ""
        @test length(e) == 0
        @test convert(String, e) === ""
        @test isa(convert(String, e), String)
        @test convert(SizedString{0}, "") == e
    end

    let s = SizedString{3}("xyZ")
        @test s == "xyZ"
        @test isless(s, "yyZ")
        @test isless("xyY", s)
        @test length(s) == sizeof(s) == 3
        @test collect(s) == ['x', 'y', 'Z']
        @test convert(String, s) == "xyZ"
        @test convert(SizedString{3}, "xyZ") == s
        @test isa(convert(SizedString{3}, "xyZ"), SizedString)
    end

    let b = IOBuffer()
        write(b, "a tEst str1ng")
        seekstart(b)
        @test read(b, SizedString{4}) == "a tE"
        @test read(b, SizedString{2}) == "st"
        b = IOBuffer()
        data = "\0Te\$t\0_"
        write(b, SizedString{length(data)}(data))
        @test String(take!(b)) == data
    end
end


text:
include("./imports.jl")

include("initialization.jl")


text:
import Base: push!
import SimilaritySearch: search
export InvIndex, prune, search
using SimilaritySearch

"""
Inverted index search structure

"""
mutable struct InvIndex <: Index
    lists::Dict{Symbol, Vector{SparseVectorEntry}}
    n::Int
    InvIndex() = new(Dict{Symbol, Vector{SparseVectorEntry}}(), 0)
    InvIndex(lists, n) = new(lists, n)
end

# useful constant for searching
const EMPTY_SPARSE_VECTOR = SparseVectorEntry[]

"""
    update!(a::InvIndex, b::InvIndex)

Updates inverted index `a` with `b`.
"""
function update!(a::InvIndex, b::InvIndex)
    for (sym, list) in b.lists
        if haskey(a.lists, sym)
            append!(a.lists[sym], list)
        else
            a.lists[sym] = list
        end
    end

    a.n += b.n
    a
end

"""
    save(f::IO, invindex::InvIndex)

Writes `invindex` to the given stream
"""
function save(f::IO, invindex::InvIndex)
    write(f, invindex.n)
    write(f, length(invindex.lists))
    for (sym, lst) in invindex.lists
        println(f, string(sym))
        write(f, length(lst))
        for p in lst
            write(f, p.id, p.weight)
        end
    end
end

"""
    load(f::IO, ::Type{InvIndex})

Reads an inverted index from the given file
"""
function load(f::IO, ::Type{InvIndex})
    invindex = InvIndex()
    invindex.n = read(f, Int)
    m = read(f, Int)
    for i in 1:m
        sym = Symbol(readline(f))
        l = read(f, Int)
        lst = Vector{SparseVectorEntry}(undef, l)
        invindex.lists[sym] = lst
        for j in 1:l
            id = read(f, Int)
            weight = read(f, Float64)
            lst[j] = SparseVectorEntry(id, weight)
        end
    end

    invindex
end

"""
    push!(index::InvIndex, objID::UInt64, bow::Dict{Symbol,Float64})

Inserts a weighted bag of words (BOW) into the index.

See [compute_bow](@ref) to compute a BOW from a text
"""
function push!(index::InvIndex, objID::Integer, bow::Dict{Symbol,Float64})
    index.n += 1
    for (sym, weight) in bow
        if haskey(index.lists, sym)
            push!(index.lists[sym], SparseVectorEntry(objID, weight))
        else
            index.lists[sym] = [SparseVectorEntry(objID, weight)]
        end
    end
end

function fit(::Type{InvIndex}, db::AbstractVector{Dict{Symbol,Float64}})
    invindex = InvIndex()
    for (i, bow) in enumerate(db)
        push!(invindex, i, bow)
    end

    invindex
end

"""
    prune(invindex::InvIndex, k)

Creates a new inverted index using the given `invindex` discarding many entries with low weight.
It keeps at most `k` entries for each posting list; it keeps those entries with more wight values.
"""
function prune(invindex::InvIndex, k)
    I = InvIndex()
    I.n = invindex.n
    sizehint!(I.lists, length(invindex.lists))

    for (t, list) in invindex.lists
        I.lists[t] = l = deepcopy(list)
        sort!(l, by=x -> -x.weight)
        if length(list) > k
            resize!(l, k)
        end
    end

    # normalizing prunned vectors
    _norm_pruned!(I)
end

function _norm_pruned!(I::InvIndex)
    D = Dict{Int,Float64}()

    for (t, list) in I.lists
        for p in list
            D[p.id] = get(D, p.id, 0.0) + p.weight * p.weight
        end
    end

    for (k, v) in D
        D[k] = 1.0 / sqrt(v)
    end

    for (t, list) in I.lists
        for p in list
            p.weight *= D[p.id]
        end
    end

    I
end

"""
    search(invindex::InvIndex, dist::Function, q::Dict{Symbol, R}, res::KnnResult) where R <: Real

Seaches for the k-nearest neighbors of `q` inside the index `invindex`. The number of nearest
neighbors is specified in `res`; it is also used to collect the results. Returns the object `res`.
If `dist` is set to `angle_distance` then the angle is reported; otherwise the
`cosine_distance` (i.e., 1 - cos) is computed.
"""
function search(invindex::InvIndex, dist::Function, q::Dict{Symbol, R}, res::KnnResult) where R <: Real
    D = Dict{Int, Float64}()
    # normalize!(q) # we expect a normalized q
    for (sym, weight) in q
        lst = get(invindex.lists, sym, EMPTY_SPARSE_VECTOR)
        if length(lst) > 0
            for e in lst
                D[e.id] = get(D, e.id, 0.0) + weight * e.weight
            end
        end
    end

    for (i, w) in D
        if dist == angle_distance
            w = max(-1.0, w)
            w = min(1.0, w)
            w = acos(w)
            push!(res, i, w)
        else
            push!(res, i, 1.0 - w)  # cosine distance
        end
    end

    res
end


text:
@polly function kernel_atax(A, x, y, tmp)
    m,n = size(A)

    for i = 1:n
        y[i] = 0
    end

    for i = 1:m
        tmp[i] = zero(eltype(tmp))

        for j = 1:n
            tmp[i] = tmp[i] + A[i,j] * x[j]
        end

        for j = 1:n
            y[j] = y[j] + A[i,j] * tmp[i]
        end
    end
end

let
    m = 1900
    n = 2100

    A = zeros(Float32, m, n)
    x = zeros(Float32, n)
    y = zeros(Float32, n)
    tmp = zeros(Float32, m)

    for i = 1:n
        x[i] = 1 + (i / n)
    end
    for i = 1:m, j = 1:n
        A[i,j] = ((i+j) % n) / (5*m)
    end

    SUITE["atax"] = @benchmarkable kernel_atax(A, x, y, tmp) setup = (A = copy($A); x = copy($x); y = copy($y); tmp = copy($tmp))
end


text:
arr = [1,2,3]
# arr[0] #=> causes a program stop with: ERROR: BoundsError()
# sqrt(-3) #=>  causes an ERROR: DomainError: sqrt will only
sqrt(complex(-3)) #=> 0.0 + 1.7320508075688772im

type CustomException <: Exception end

# throw:
codes = ["AO", "ZD", "SG", "EZ"]
code = "AR"
if in(code, codes)
    println("This is an acceptable code")
else
    throw(DomainError())
end

# error:
function sqroot(n::Int)
    n >= 0 || error("n must be non-negative")
    n == 0 && return 0
    sqrt(n)
end
sqroot(4) #=> 2.0
sqroot(0) #=> 0.0
# sqroot(-6) #=> ERROR: n must be non-negative

warn("Something is not right here")
info("Did you know this?")

type CustomException <: Exception
end

# try / catch
a = []
try
    pop!(a)
catch ex
    println(typeof(ex))  #=> ErrorException
    backtrace()
    showerror(STDOUT, ex) #=> array must be non-empty
end

try
   # try this code
catch ex
   if isa(ex, DomainError)
# do this
   elseif isa(ex, BoundsError)
# do this
   end
end

ret = try
         a = 4 * 2
      catch ex
      end
println(ret) #=> 8

# finally
f = open("file1.txt") #=> IOStream(<file file1.txt>)
try
    # operate on file f
finally
    close(f)
end

max_stars_repo_path: Module 1/Chapter04/Old/errors.jl

max_issues_repo_path: Module 1/Chapter04/Old/errors.jl

max_forks_repo_path: Module 1/Chapter04/Old/errors.jl


text:
#=
The code here is mainly deals with propagating derivative information between input and
output values. Usually, this means incrementing/decrementing the input's derivative(s) by an
amount scaled by the output's derivative(s). Sometimes, extra partial information (in the
form of normal scalars/arrays, or `Dual` numbers) needs to accounted for as well.

Often, a function's input and output values are not similarly shaped. To account for these
cases, `broadcast` and `reduce` versions of the propagation functions have been implemented.
These special cases should be decently fast, but some might be prone to numerical error if
the accumulated derivative becomes too large compared to the individual terms being added to
it. This can be overcome by using the divide-and-conquer strategy from Base.mapreducedim,
but that strategy is less cache efficient and more complicated to implement.

A lot of the code here is pretty repetitive, and only covers the use cases that actually
arise from the derivative definitions implemented elsewhere in ReverseDiff. At some point,
we should figure out a cleaner, more general implementation pattern that doesn't sacrifice
efficiency.
=#

#############
# utilities #
#############

index_bound(x::Any, ::AbstractArray{T,N}) where {T,N} = nothing

index_bound(x::AbstractArray, ::AbstractArray{T,N}) where {T,N} = CartesianIndex{N}(ntuple(i -> size(x, i), Val(N)))

###################
# increment_deriv #
###################

@inline increment_deriv!(t::TrackedArray, x::AbstractArray, i) = (t.deriv[i] += x[i]; nothing)
@inline increment_deriv!(t::TrackedArray, x::Real, i) = (t.deriv[i] += x; nothing)
@inline increment_deriv!(t::AbstractArray, x::AbstractArray, i) = increment_deriv!(t[i], x[i])
@inline increment_deriv!(t::AbstractArray, x::Real, i) = increment_deriv!(t[i], x)

function increment_deriv!(t::AbstractArray, x)
    for i in eachindex(t)
        increment_deriv!(t, x, i)
    end
    return nothing
end

function increment_deriv!(t::TrackedReal, x::Real)
    pull_deriv!(t)
    t.deriv += x
    push_deriv!(t)
    return nothing
end

###################
# decrement_deriv #
###################

@inline decrement_deriv!(t::TrackedArray, x::AbstractArray, i) = (t.deriv[i] -= x[i]; nothing)
@inline decrement_deriv!(t::TrackedArray, x::Real, i) = (t.deriv[i] -= x; nothing)
@inline decrement_deriv!(t::AbstractArray, x::AbstractArray, i) = decrement_deriv!(t[i], x[i])
@inline decrement_deriv!(t::AbstractArray, x::Real, i) = decrement_deriv!(t[i], x)

function decrement_deriv!(t::AbstractArray, x)
    for i in eachindex(t)
        decrement_deriv!(t, x, i)
    end
    return nothing
end

function decrement_deriv!(t::TrackedReal, x::Real)
    pull_deriv!(t)
    t.deriv -= x
    push_deriv!(t)
    return nothing
end

###############################
# diffresult_increment_deriv! #
###############################

@inline getpartial(r::DiffResults.ImmutableDiffResult{1,V,Tuple{D}}, p) where {V,D<:AbstractArray} = DiffResults.derivative(r)[p]
@inline getpartial(r::DiffResults.ImmutableDiffResult{1,V,Tuple{D}}, p) where {V,D<:Number} = DiffResults.derivative(r)

function diffresult_increment_deriv!(input::AbstractArray, x::AbstractArray,
                                     results, p::Int)
    for i in eachindex(x)
        increment_deriv!(input, x[i] * getpartial(results[i], p), i)
    end
    return nothing
end

function diffresult_increment_deriv!(input::AbstractArray, x::AbstractArray,
                                     results, p::Int, bound::CartesianIndex)
    for i in CartesianIndices(size(x))
        increment_deriv!(input, x[i] * getpartial(results[i], p), min(bound, i))
    end
    return nothing
end

function diffresult_increment_deriv!(input::TrackedReal, x::AbstractArray,
                                     results, p::Int, ::Nothing)
    pull_deriv!(input)
    input_deriv = input.deriv
    for i in eachindex(results)
        input_deriv += x[i] * getpartial(results[i], p)
    end
    input.deriv = input_deriv
    push_deriv!(input)
    return nothing
end

##############################
# broadcast_increment_deriv! #
##############################

# without partials #
#------------------#

function broadcast_increment_deriv!(input::AbstractArray, x::AbstractArray,
                                    bound::CartesianIndex)
    for i in CartesianIndices(size(x))
        increment_deriv!(input, x[i], min(bound, i))
    end
    return nothing
end

function broadcast_increment_deriv!(input::TrackedReal, x::AbstractArray, ::Nothing)
    pull_deriv!(input)
    input_deriv = input.deriv
    for i in x
        input_deriv += i
    end
    input.deriv = input_deriv
    push_deriv!(input)
    return nothing
end

# with partial array #
#--------------------#

function broadcast_increment_deriv!(input::AbstractArray, x::AbstractArray,
                                    partials::AbstractArray,
                                    input_bound::CartesianIndex,
                                    partials_bound::CartesianIndex)
    for i in CartesianIndices(size(x))
        current_deriv = x[i] * partials[min(partials_bound, i)]
        increment_deriv!(input, current_deriv, min(input_bound, i))
    end
    return nothing
end

function broadcast_increment_deriv!(input::TrackedReal, x::AbstractArray,
                                    partials::AbstractArray, ::Nothing, ::CartesianIndex)
    pull_deriv!(input)
    input_deriv = input.deriv
    for i in eachindex(x)
        input_deriv += x[i] * partials[i]
    end
    input.deriv = input_deriv
    push_deriv!(input)
    return nothing
end

# with partial scalar #
#---------------------#

@inline function broadcast_increment_deriv!(input, x, partial::Ref, i, j)
    return broadcast_increment_deriv!(input, x, partial[], i, j)
end

function broadcast_increment_deriv!(input::AbstractArray, x::AbstractArray,
                                    partial::Real, input_bound::CartesianIndex,
                                    ::Nothing)
    for i in CartesianIndices(size(x))
        increment_deriv!(input, x[i] * partial, min(input_bound, i))
    end
    return nothing
end

function broadcast_increment_deriv!(input::TrackedReal, x::AbstractArray,
                                    partial::Real, ::Nothing, ::Nothing)
    pull_deriv!(input)
    input_deriv = input.deriv
    for i in eachindex(x)
        input_deriv += x[i] * partial
    end
    input.deriv = input_deriv
    push_deriv!(input)
    return nothing
end

##############################
# broadcast_decrement_deriv! #
##############################

function broadcast_decrement_deriv!(input::AbstractArray, x::AbstractArray,
                                    bound::CartesianIndex)
    for i in CartesianIndices(size(x))
        decrement_deriv!(input, x[i], min(bound, i))
    end
    return nothing
end

function broadcast_decrement_deriv!(input::TrackedReal, x::AbstractArray, ::Nothing)
    pull_deriv!(input)
    input_deriv = input.deriv
    for i in x
        input_deriv -= i
    end
    input.deriv = input_deriv
    push_deriv!(input)
    return nothing
end

##############################
# reduction_increment_deriv! #
##############################

function reduction_increment_deriv!(input::AbstractArray, x::AbstractArray,
                                    bound::CartesianIndex)
    for i in CartesianIndices(size(input))
        increment_deriv!(input, x[min(bound, i)], i)
    end
    return nothing
end


text:
function  exasim(pde,mesh)

# search compilers and set options
pde = Gencode.setcompilers(pde);

# generate input files and store them in datain folder
pde, mesh, master, dmd = Preprocessing.preprocessing(pde,mesh);

# generate source codes and store them in app folder
Gencode.gencode(pde);

if pde.usecmake==1
    runstr = Main.cmakecompile(pde); # use cmake to compile source codes
    compilerstr = "";
else
    # compile source codes to build an executable file and store it in app folder
    compilerstr = Gencode.compilecode(pde);

    # run executable file to compute solution and store it in dataout folder
    runstr = Gencode.runcode(pde);
end

# get solution from output files in dataout folder
sol = Postprocessing.fetchsolution(pde,master,dmd);

return sol,pde,mesh,master,dmd,compilerstr,runstr

end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
module Day6

export solvePart1, solvePart2

import DelimitedFiles
import Distances

"""Get a sorted array of tuples containing xy coordinates."""
function getInput()
    open("day6input.txt") do f
        input = DelimitedFiles.readdlm(f, ',', Int64, '\n')
        sorted = sortslices(input, dims=1)
        mapslices(row -> (row[1], row[2]), sorted, dims=2)[:]
    end
end

""""Create a grid of coordinates to consider, 1 unit larger than the bounding
box of coordinates."""
function createGrid(input)
    xmin = minimum(map(first, input))
    ymin = minimum(map(tuple -> tuple[2], input))
    xmax = maximum(map(first, input))
    ymax = maximum(map(tuple -> tuple[2], input))
    xs = xmin-1:xmax+1
    ys = ymin-1:ymax+1
    return collect(Base.product(xs, ys))
end

"""Create a dictionary mapping from (x, y) to the distances to each point."""
function findDistances(grid, points)
    allDistances = Dict()
    for square in grid
        distances = map(point -> Distances.cityblock([x for x in point],
                                                     [x for x in square]),
                        points)
        allDistances[square] = distances
    end
    return allDistances
end

"""Using the distance map, create a mapping to the closest point."""
function findClosestPoints(distancesMap, points)
    closest = Dict()
    for square in keys(distancesMap)
        distances = distancesMap[square]
        min_index = argmin(distances)
        min_point = points[min_index]
        min_distance = distances[min_index]
        if count(d -> d == min_distance, distances) == 1
            closest[square] = min_point
        else
            closest[square] = "tie"
        end
    end
    return closest
end

"""Find all points that are not closest to an edge square."""
function getEligiblePointsSet(grid, points, closest)
    eligible = Set(points)

    (m, n) = size(grid)
    # Edges along the x direction
    for i in 1:m
        if closest[grid[i, 1]] in eligible
            pop!(eligible, closest[grid[i, 1]])
        end
        if closest[grid[i, n]] in eligible
            pop!(eligible, closest[grid[i, n]])
        end
    end
    # Edges along the y direction
    for i in 1:n
        if closest[grid[1, i]] in eligible
            pop!(eligible, closest[grid[1, i]])
        end
        if closest[grid[m, i]] in eligible
            pop!(eligible, closest[grid[m, i]])
        end
    end
    return eligible
end

"""Find the size of the regions around each point."""
function findRegionSizes(closest, eligible)
    regionSizes = Dict()
    for square in keys(closest)
        p = closest[square]
        if p in eligible
            regionSizes[p] = get(regionSizes, p, 0) + 1
        end
    end
    return regionSizes
end

function solvePart1()
    points = getInput()
    grid = createGrid(points)
    allDistances = findDistances(grid, points)
    closest = findClosestPoints(allDistances, points)
    eligiblePoints = getEligiblePointsSet(grid, points, closest)
    regionSizes = findRegionSizes(closest, eligiblePoints)
    return max(values(regionSizes)...)
end

function solvePart2()
    points = getInput()
    grid = createGrid(points)
    allDistances = findDistances(grid, points)

    totalDistances = Dict()
    for (square, distances) in allDistances
        totalDistances[square] = sum(distances)
    end

    return length(filter(d -> d < 10000, collect(values(totalDistances))))
end
end


text:
# Test file for 5-d Chebyshev polynomials

using ChebyshevApprox

n1 = 10
n2 = 15
n3 = 7
n4 = 16
n5 = 31

dom_1 = [2.0, -3.0]
dom_2 = [1.5,  0.5]
dom_3 = [0.5, -0.5]
dom_4 = [-0.66, -1.79]
dom_5 = [1.0, -1.0]
dom = [dom_1 dom_2 dom_3 dom_4 dom_5]

nodesr_1  = chebyshev_nodes(n1,dom_1)
nodesr_2  = chebyshev_nodes(n2,dom_2)
nodesr_3  = chebyshev_nodes(n3,dom_3)
nodesr_4  = chebyshev_nodes(n4,dom_4)
nodesr_5  = chebyshev_nodes(n5,dom_5)
nodese_1  = chebyshev_extrema(n1,dom_1)
nodese_2  = chebyshev_extrema(n2,dom_2)
nodese_3  = chebyshev_extrema(n3,dom_3)
nodese_4  = chebyshev_extrema(n4,dom_4)
nodese_5  = chebyshev_extrema(n5,dom_5)
nodesex_1 = chebyshev_extended(n1,dom_1)
nodesex_2 = chebyshev_extended(n2,dom_2)
nodesex_3 = chebyshev_extended(n3,dom_3)
nodesex_4 = chebyshev_extended(n4,dom_4)
nodesex_5 = chebyshev_extended(n5,dom_5)
nodesv_1  = vertesi_nodes(n1,dom_1)
nodesv_2  = vertesi_nodes(n2,dom_2)
nodesv_3  = vertesi_nodes(n3,dom_3)
nodesv_4  = vertesi_nodes(n4,dom_4)
nodesv_5  = vertesi_nodes(n5,dom_5)

yr  = zeros(n1,n2,n3,n4,n5)
ye  = zeros(n1,n2,n3,n4,n5)
yex = zeros(n1,n2,n3,n4,n5)
yv  = zeros(n1,n2,n3,n4,n5)

for i = 1:n1
  for j = 1:n2
    for k = 1:n3
      for l = 1:n4
        for m = 1:n5

          yr[i,j,k,l,m]  = (nodesr_1[i]+4.0)^0.5+nodesr_1[i]*sqrt(nodesr_2[j])+exp(nodesr_3[k])*nodesr_4[l]-(1.0+nodesr_5[m])^2
          ye[i,j,k,l,m]  = (nodese_1[i]+4.0)^0.5+nodese_1[i]*sqrt(nodese_2[j])+exp(nodese_3[k])*nodese_4[l]-(1.0+nodese_5[m])^2
          yex[i,j,k,l,m] = (nodesex_1[i]+4.0)^0.5+nodesex_1[i]*sqrt(nodesex_2[j])+exp(nodesex_3[k])*nodesex_4[l]-(1.0+nodesex_5[m])^2
          yv[i,j,k,l,m]  = (nodesv_1[i]+4.0)^0.5+nodesv_1[i]*sqrt(nodesv_2[j])+exp(nodesv_3[k])*nodesv_4[l]-(1.0+nodesv_5[m])^2

        end
      end
    end
  end
end

order_tensor = [6, 6, 6, 6, 6]
order_complete = 6

wr_tensor    = chebyshev_weights_threaded(yr,[nodesr_1,nodesr_2,nodesr_3,nodesr_4,nodesr_5],order_tensor,dom)
wr_complete  = chebyshev_weights_threaded(yr,(nodesr_1,nodesr_2,nodesr_3,nodesr_4,nodesr_5),order_complete,dom)
we_tensor    = chebyshev_weights_extrema_threaded(ye,(nodese_1,nodese_2,nodese_3,nodese_4,nodese_5),order_tensor,dom)
we_complete  = chebyshev_weights_extrema_threaded(ye,(nodese_1,nodese_2,nodese_3,nodese_4,nodese_5),order_complete,dom)
wex_tensor   = chebyshev_weights_extended_threaded(yex,(nodesex_1,nodesex_2,nodesex_3,nodesex_4,nodesex_5),order_tensor,dom)
wex_complete = chebyshev_weights_extended_threaded(yex,(nodesex_1,nodesex_2,nodesex_3,nodesex_4,nodesex_5),order_complete,dom)
wv_tensor    = chebyshev_weights_vertesi_threaded(yv,(nodesv_1,nodesv_2,nodesv_3,nodesv_4,nodesv_5),order_tensor,dom)
wv_complete  = chebyshev_weights_vertesi_threaded(yv,(nodesv_1,nodesv_2,nodesv_3,nodesv_4,nodesv_5),order_complete,dom)

point = [1.748, 0.753, 0.119, -0.947, -0.23]

yr_tensor    = chebyshev_evaluate(wr_tensor,point,order_tensor,dom)
yr_complete  = chebyshev_evaluate(wr_complete,point,order_complete,dom)
ye_tensor    = chebyshev_evaluate(we_tensor,point,order_tensor,dom)
ye_complete  = chebyshev_evaluate(we_complete,point,order_complete,dom)
yex_tensor   = chebyshev_evaluate(wex_tensor,point,order_tensor,dom)
yex_complete = chebyshev_evaluate(wex_complete,point,order_complete,dom)
yv_tensor    = chebyshev_evaluate(wv_tensor,point,order_tensor,dom)
yv_complete  = chebyshev_evaluate(wv_complete,point,order_complete,dom)

y_actual = (point[1]+4.0)^0.5+point[1]*sqrt(point[2])+exp(point[3])*point[4]-(1.0+point[5])^2

println([y_actual yr_tensor yr_complete ye_tensor ye_complete yex_tensor yex_complete yv_tensor yv_complete])

yr_grad_tensor       = chebyshev_gradient(wr_tensor,point,order_tensor,dom)
yr_grad_complete     = chebyshev_gradient(wr_complete,point,order_complete,dom)
yr_deriv_tensor_1    = chebyshev_derivative(wr_tensor,point,1,order_tensor,dom)
yr_deriv_complete_1  = chebyshev_derivative(wr_tensor,point,1,order_complete,dom)
ye_grad_tensor       = chebyshev_gradient(we_tensor,point,order_tensor,dom)
ye_grad_complete     = chebyshev_gradient(we_complete,point,order_complete,dom)
ye_deriv_tensor_1    = chebyshev_derivative(we_tensor,point,1,order_tensor,dom)
ye_deriv_complete_1  = chebyshev_derivative(we_complete,point,1,order_complete,dom)
yex_grad_tensor      = chebyshev_gradient(wex_tensor,point,order_tensor,dom)
yex_grad_complete    = chebyshev_gradient(wex_complete,point,order_complete,dom)
yex_deriv_tensor_1   = chebyshev_derivative(wex_tensor,point,1,order_tensor,dom)
yex_deriv_complete_1 = chebyshev_derivative(wex_complete,point,1,order_complete,dom)
yv_grad_tensor       = chebyshev_gradient(wv_tensor,point,order_tensor,dom)
yv_grad_complete     = chebyshev_gradient(wv_complete,point,order_complete,dom)
yv_deriv_tensor_1    = chebyshev_derivative(wv_tensor,point,1,order_tensor,dom)
yv_deriv_complete_1  = chebyshev_derivative(wv_complete,point,1,order_complete,dom)

println(yr_grad_tensor, yr_grad_complete, yr_deriv_tensor_1, yr_deriv_complete_1)
println(ye_grad_tensor, ye_grad_complete, ye_deriv_tensor_1, ye_deriv_complete_1)
println(yex_grad_tensor, yex_grad_complete, yex_deriv_tensor_1, yex_deriv_complete_1)
println(yv_grad_tensor, yv_grad_complete, yv_deriv_tensor_1, yv_deriv_complete_1)

# Test the structures

chebpoly = ChebPoly(wr_complete,order_complete,dom)
chebyshev_evaluate(chebpoly,point)

chebr  = ChebInterpRoots(yr,(nodesr_1,nodesr_2,nodesr_3,nodesr_4,nodesr_5),order_complete,dom)
chebe  = ChebInterpExtrema(ye,(nodese_1,nodese_2,nodese_3,nodese_4,nodese_5),order_complete,dom)
chebex = ChebInterpExtended(yex,(nodesex_1,nodesex_2,nodesex_3,nodesex_4,nodesex_5),order_complete,dom)
chebv  = ChebInterpVertesi(yv,(nodesv_1,nodesv_2,nodesv_3,nodesv_4,nodesv_5),order_complete,dom)
cheb_interp(chebr,point)
cheb_interp(chebe,point)
cheb_interp(chebex,point)
cheb_interp(chebv,point)

chebr_fn  = cheb_interp(chebr)
chebe_fn  = cheb_interp(chebe)
chebex_fn = cheb_interp(chebex)
chebv_fn  = cheb_interp(chebv)

chebr_fn(point)
chebe_fn(point)
chebex_fn(point)
chebv_fn(point)

max_issues_repo_name: RJDennis/ChebyshevApprox

max_forks_repo_name: RJDennis/ChebyshevApprox


text:
using Distributions
using PSIS
using Random
using Test

@testset "Generalized Pareto distribution" begin
    @testset "fit_gpd" begin
        @testset for μ in (-1, 5),
            σ in (0.5, 1.0, 2.0),
            k in (-1.0, 0.0, 0.3, 1.0, 2.0),
            improved in (true, false)

            d = Distributions.GeneralizedPareto(μ, σ, k)
            rng = MersenneTwister(42)
            x = rand(rng, d, 200_000)
            dhat = PSIS.fit_gpd(x; μ=μ, min_points=80, improved=improved)
            @test dhat.μ == μ
            @test dhat.σ ≈ σ atol = 0.01
            @test dhat.k ≈ k atol = 0.01
        end
        @testset "nearly uniform" begin
            x = ones(200_000)
            dhat = PSIS.fit_gpd(x; μ=1.0, min_points=80)
            @test dhat.μ == 1.0
            @test dhat.σ == 0.0
            @test dhat.k == -1
        end
    end
    @testset "convert to Distributions type" begin
        d = PSIS.GeneralizedPareto(1.0, 2.0, 3.0)
        @test @inferred(convert(Distributions.GeneralizedPareto, d)) ===
            Distributions.GeneralizedPareto(1.0, 2.0, 3.0)
        d = PSIS.GeneralizedPareto(0.0, 0.0, -1.0)
        @test @inferred(convert(Distributions.GeneralizedPareto{Float32}, d)) ===
            Distributions.GeneralizedPareto{Float32}(0.0f0, 0.0f0, -1.0f0)
    end
end


text:
export
    SetRule,
    dnaSize,
    resultSize,
    initialize!,
    decode,
    asValue

immutable SetRule{T} <: SynthRule
    count::Int
    items::Vector{T}
    dnaSize::Int
end

function SetRule{T}(items::Vector{T}, count = 1)
    dnaSize = length(items) * count

    SetRule{T}(count, items, dnaSize)
end

dnaSize(rule::SetRule) = rule.dnaSize

resultSize(rule::SetRule) = rule.count

function initialize!(rule::SetRule, ctx, state) end

function decode(rule::SetRule, pars)
    scope!(pars.ctx) do this
        const depth = length(rule.items)
        dist = moddims(pars.dnaFragment, rule.count, depth)
        rnd = randu(pars.ctx, Float32, rule.count, depth)
        test = rnd ./ maxOf(dist, 0.000001f0)
        max, idx = imax(test, 1)
        pars.result[pars.resultSeq] = idx
    end
end

function asValues(rule::SetRule, values, startIndex)
    result = Vector(resultSize(rule))
    for i in 1:resultSize(rule)
        itemIndex = values[startIndex - 1 + i]
        result[i] = rule.items[Int(round(itemIndex + 1))]
    end
    result
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
abstract type Group end

Base.broadcasted(f::Function, x::Group, y::AbstractVector{<:Integer}) = f.((x for i in 1:length(y)), y)

Base.broadcasted(::typeof(*), x::G, y::Vector{G}) where G <: Group = (x for i in 1:length(y)) .* y
Base.broadcasted(::typeof(*), x::Vector{G}, y::G) where G <: Group =  x .* (y for i in 1:length(x))

order(x::G) where G <: Group = order(G)

import Base./
/(x::G, y::G) where G <: Group = x * inv(y)

name(x::G) where G <: Group = name(G)

struct ECGroup{P<:ECPoint} <: Group
    x::P
end

<|(::Type{ECGroup{P}}, x) where P <: ECPoint = ECGroup{P}(P <| x)

Base.:*(x::G, y::G) where G <: ECGroup = G(x.x + y.x)

function Base.:^(x::G, n::Integer) where G <: ECGroup
    n_mod = mod(n, order(G))

    @assert n_mod != 0 "A bad exponent"

    return G(n_mod * x.x)
end

Base.inv(g::G) where G <: ECGroup = g^(order(G) - 1)

order(::Type{ECGroup{P}}) where P = order(P)

validate(x::ECGroup) = validate(x.x)

Base.:(==)(x::G, y::G) where G <: ECGroup = x.x == y.x

gx(p::ECGroup) = gx(p.x)
gy(p::ECGroup) = gy(p.x)

modulus(::Type{ECGroup{P}}) where P <: ECPoint = modulus(P)

name(::Type{ECGroup{P}}) where P <: ECPoint = name(P)

Base.isless(x::G, y::G) where G <: ECGroup = gx(x) == gx(y) ? gx(x) < gx(y) : gy(x) < gy(y)

function Base.show(io::IO, g::G) where G <: ECGroup
    show(io, G)
    print(io, " <| (")
    show(io, gx(g))
    print(io, ", ")
    show(io, gy(g))
    print(io, ")")
end

# Static fields for PGroup
struct static_PGroup{N} ### Type parameter is essential to ensure it to be bitstype
    p::StaticBigInt{N}
    q::StaticBigInt{N} # for the cases where order is unkown I could have a 0 here.
    name::UInt128

    static_PGroup(p::StaticBigInt{N}, q::StaticBigInt{N}, name::Symbol) where N = new{N}(p, q, string2uint(string(name)))
    function static_PGroup(p::Integer, q::Integer, name::Symbol)

        n = bitlength(p)

        _p = StaticBigInt(p; n)
        _q = StaticBigInt(q; n)
        static_PGroup(_p, _q, name)
    end

    static_PGroup(p::Integer, q::Integer) where N = static_PGroup(p, q, Symbol(""))
end

struct PGroup{S} <: Group
    g::BigInt

    function PGroup{S}(x::BigInt) where S
        #@assert 1 < x < S.p "Not in range"
        # Relaxing a little so that group products could happen propeerly
        # But is a presence of 1 in multiplicative products a vulnerability?
        @assert 0 < x < S.p "Not in range"
        new{S}(x)
    end

    PGroup{S}(x::Integer) where S = PGroup{S}(BigInt(x))
end

specialize(::Type{PGroup}, p::Integer, q::Integer) = PGroup{static_PGroup(p, q)}
specialize(::Type{PGroup}, q::Integer) = specialize(PGroup, 2*q + 1, q) # Safe prime constructor

specialize(::Type{PGroup}, p::Integer, q::Integer, name::Symbol) = PGroup{static_PGroup(p, q, name)}
specialize(::Type{PGroup}, q::Integer, name::Symbol) = specialize(PGroup, 2*q + 1, q, name) # Safe prime constructor

modulus(::Type{PGroup{S}}) where S = BigInt(S.p)
modulus(::G) where G <: PGroup = modulus(G)

order(::Type{PGroup{S}}) where S = BigInt(S.q)

name(::Type{PGroup}) = nothing

function name(::Type{PGroup{S}}) where S

    str = uint2string(S.name)

    if str == ""
        return nothing
    else
        return Symbol(str)
    end
end

Base.show(io::IO, g::PGroup) = print(io, value(g))

Base.show(io::IO, ::Type{PGroup}) = print(io, "PGroup")

groupstr(m) = "𝐙/($(trimnumber(m)))"

function Base.show(io::IO, ::Type{G}) where G <: PGroup

    if name(G) == nothing
        print(io, groupstr(modulus(G)))
    else
        print(io, name(G))
    end
end

function Base.display(::Type{G}) where G <: PGroup
    show(G)
    print(" (order = $(order(G)))")
end

value(g::PGroup) = g.g

<|(::Type{P}, x::Integer) where P <: PGroup = P(BigInt(x))

validate(g::G) where G <: PGroup = value(g) != 1 && powermod(value(g), order(G), modulus(G)) == 1

import Base.*
*(x::G, y::G) where G <: PGroup = G(mod(value(x) * value(y), modulus(G)))

import Base.^
function ^(x::G, n::Integer) where G <: PGroup

    n_mod = mod(n, order(G))
    @assert n_mod != 0 "A bad exponent"
    @assert value(x) != 1 "A value 1 is not part of prime group"

    return G(powermod(value(x), n_mod, modulus(G)))
end

import Base.inv
inv(x::G) where G <: PGroup = G(modinv(value(x), modulus(G)))

Base.:(==)(x::G, y::G) where G <: PGroup = x.g == y.g

Base.isless(x::G, y::G) where G <: PGroup = value(x) < value(y)

Base.convert(::Type{G}, x::Integer) where G <: PGroup = G(x)

# function Base.prod(x::Vector{G}) where G <: PGroup

#     p = modulus(G)
#     s = value(x[1])

#     for i in x[2:end]
#         s = mod(s * value(i), p)
#     end

#     return G(s)
# end

max_stars_repo_name: PeaceFounder/DiffieHellmanGroups

max_issues_repo_name: PeaceFounder/DiffieHellmanGroups

max_forks_repo_name: PeaceFounder/DiffieHellmanGroups


text:
struct WireNet{P}
    specified_connections::Vector{Vector{Observable{P}}}
    connections::Observable{Vector{Vector{P}}}
    function WireNet(connections::Vector{Vector{Q}};P=Point2f) where Q
        obs = filter(x->isa(x,Observable), collect(flatten(connections)))

        new_connections = Observable([
            [zero(P) for i in 1:(length(connection)*2-1)] for connection in connections
        ])

        function update!(args...)
            for (connection,spec_connection) in zip(new_connections[],connections)
                end_horizontal = if length(spec_connection) >= 3
                    n1,n2,n3 = to_value.(spec_connection[end-2:end])
                    !(n1[1] ≤ n2[1] ≤ n3[1] || n1[1] ≥ n2[1] ≥ n3[1])
                else
                    false
                end

                for i in reverse(eachindex(spec_connection[1:end-1]))
                    n1,n2 = to_value.(spec_connection[i:i+1])
                    intermediate = if n1[1] == n2[1]
                        end_horizontal = true
                        (n1 + n2)/2
                    elseif n1[2] == n2[2]
                        end_horizontal = false
                        (n1 + n2)/2
                    elseif end_horizontal
                        end_horizontal = false
                        P(n1[1],n2[2])
                    else
                        end_horizontal = true
                        P(n2[1],n1[2])
                    end
                    connection[2i-1] = n1
                    connection[2i] = intermediate
                end
                connection[end] = to_value(spec_connection[end])
            end
            notify(new_connections)
        end

        onany(update!, obs...)
        update!()

        new{P}(connections, new_connections)
    end
end

##
##
function intersectSegments(A::P, B::P, C::P, D::P)::Union{Nothing, P} where P
    ΔAB = B-A
    ΔAC = C-A
    ΔAD = D-A

    # fast check of bounding box
    lu1 = min.(A,B)
    ro1 = max.(A,B)
    lu2 = min.(C,D)
    ro2 = max.(C,D)
    if any(ro2 .< lu1) || any(ro1 .< lu2)
        # bounding boxes don't intersect -> no chance for overlap
        return nothing
    end

    ΔAB′ = ΔAB/(ΔAB⋅ΔAB)

    # get first projection
    ρ1 = (ΔAB′⋅ΔAC)
    δ1 = ΔAC-ρ1 * ΔAB
    # get second projection
    ρ2 = (ΔAB′⋅ΔAD)
    δ2 = ΔAD-ρ2 * ΔAB

    # check if both δ point in the same direction
    if δ1⋅δ2 > 0
        # they point in the same direction -> no intersection
        return nothing
    else
        # they point in opposite directions -> compute intersection
        d1 = √(δ1⋅δ1)
        d2 = √(δ2⋅δ2)

        d12 = d1+d2
        if d12 ≈ 0
            # CD  must lie on the line AB -> Error
            AssertionError("Line segments $((A,B)) and $((C,D)) are colinear!")
        end

        ρ3 = d1 / d12 * (ρ2-ρ1) + ρ1
        if 0.0 ≤ ρ3 ≤ 1.0
            return A + ρ3 * ΔAB
        else
            return nothing
        end
    end
end

function splitSegment(p_start,p_end,p_split,r_split)
    Δ1 = p_split-p_start
    Δ2 = p_split-p_end
    n1 = norm(Δ1)
    n2 = norm(Δ2)

    ret1 = n1 < r_split ? nothing : (p_start, p_split-r_split*Δ1/n1)
    ret2 = n2 < r_split ? nothing : (p_split-r_split*Δ2/n2, p_end)
    return ret1, ret2
end

##

(anglepoint(center::P, angle, radius) where {P}) = center + P(cos(angle),sin(angle))*radius

@recipe(WirePlot, nets) do scene
    Attributes(
        bridge_radius=0.1,
        linewidth=3,
        csegs=20,
        net_colors=DefaultDict{Symbol,Union{Symbol,RGBf,RGBAf}}(:black)
    )
end

function Makie.plot!(wireplot::WirePlot)
    nets = wireplot[1][]
    net_names = keys(nets)
    all_wires = Dict((k=>Observable(Point2f[]) for k in net_names)...)
    function update_wires!(bridge_radius, csegs, nets...)
        all_segments=[[Tuple{Point2f,Point2f}[] for connection in net] for net in nets]
        for (i,(key,net)) in enumerate(zip(net_names, nets))
            line = all_wires[key][]
            empty!(line)
            for (j,connection) in enumerate(net)
                for k in 1:length(connection)-1
                    segments = Union{Nothing, Tuple{Point2f,Point2f}}[(connection[k],connection[k+1])]
                    while !isempty(segments)
                        segment = pop!(segments)
                        split=false
                        for other_net_connections in all_segments[1:i-1]
                            for other_net_connection in other_net_connections
                                for other_segment in other_net_connection
                                    pt = intersectSegments(segment[1],segment[2],other_segment[1],other_segment[2])
                                    if !isnothing(pt)
                                        (seg1,seg2) = splitSegment(segment[1], segment[2], pt, bridge_radius)
                                        if !isnothing(seg2)
                                            push!(segments, seg2)
                                        end
                                        if !isnothing(seg1)
                                            push!(segments, seg1)
                                        end
                                        split=true
                                        break
                                    end
                                end
                                if split
                                    break
                                end
                            end
                            if split
                                break
                            end
                        end

                        if !split
                            push!(all_segments[i][j], segment)

                            if isempty(line) || isnan(line[end])
                                append!(line,segment)
                            elseif segment[1] ≈ line[end]
                                push!(line, segment[end])
                            else
                                Δ = segment[1]-line[end]
                                slope = Δ[1] ≈ 0.0 ? sign(Δ[2])*π/2 : atan(Δ[2]/Δ[1])
                                append!(line,anglepoint.(Ref(0.5*line[end] + 0.5*segment[1]), LinRange(slope+π, slope, csegs), 0.5*norm(Δ)))
                                push!(line, segment[end])
                            end
                        end
                    end
                end
                push!(line, Point2f(NaN,NaN))
            end
            notify(all_wires[key])
        end
    end
    onany(update_wires!, wireplot[:bridge_radius], wireplot[:csegs], (net.connections for net in values(nets))...)
    notify(first(values(nets)).connections)

    col = get(wireplot.attributes,:color,Observable(:black))
    Dict(name=>to_value(col) for name in keys(all_wires))

    merged_colors = lift(wireplot[:net_colors], values(all_wires)...) do cols, wires...
        colors = []
        for (key,wire) in zip(keys(all_wires),wires)
            append!(colors, fill(cols[key], length(wire)))
        end
        colors
    end

    merged_wires = lift(values(all_wires)...) do wires...
        [wires...;]
    end
    lines!(wireplot, merged_wires; linewidth=wireplot[:linewidth], color=merged_colors)

    wireplot
end


text:
using SnoopCompile

SnoopCompile.@snoop "compiles.csv" begin
    using TableReader
    readcsv("test/test.csv")
    readtsv("test/test.tsv")
end

data = SnoopCompile.read("compiles.csv")
pc = SnoopCompile.parcel(reverse!(data[2]))
SnoopCompile.write("precompile", pc)