text:
TESTS = [
    "free" => [
        (
            a,
            [
                (a, [Dict()]),
                (@term(a(b)), []),
                (b, []),
                (f(a), []),
            ]
        ),
        (
            f(a),
            [
                (f(a), [Dict()]),
                (f(), []),
                (f(a, b), []),
                (f(f(a)), []),
                (g(a), []),
                (f(b), []),
                (g(b), []),
            ]
        ),
        (
            f(x),
            [
                (f(a), [Dict(x => a)]),
                (f(b), [Dict(x => b)]),
                (f(f(a)), [Dict(x => f(a))]),
                (f(), []),
                (g(a), []),
                (g(f(a)), []),
            ]
        ),
        (
            g(x, y),
            [
                (g(a, a), [Dict(x => a, y => a)]),
                (g(a, b), [Dict(x => a, y => b)]),
                (g(a, b, c), []),
                (g(), []),
            ]
        ),
        (
            g(x, x),
            [
                (g(a, a), [Dict(x => a)]),
                (g(a, b), []),
                (g(a, b, c), []),
                (g(), []),
            ]
        ),
        (
            g(x, h(x, y, c)),
            [
                (g(a, h(a, a, c)), [Dict(x => a, y => a)]),
                (g(a, h(a, b, c)), [Dict(x => a, y => b)]),
                (g(a, h(a, c, c)), [Dict(x => a, y => c)]),
                (g(a, h(b, b, c)), []),
                (g(a, f(a, b, c)), []),
                (f(a, h(a, b, b)), []),
                (g(a, h(a, b, b)), []),
            ]
        ),
    ],
    "c" => [
        (
            p(a, a),
            [
                (p(a, a), [Dict()]),
                (p(a, b), []),
                (q(a, a), []),
            ]
        ),
        (
            p(a, b),
            [
                (p(a, b), [Dict()]),
                (p(a, a), []),
                (q(a, b), []),
            ]
        ),
        (
            p(x, a),
            [
                (p(a, a), [Dict(x => a)]),
                (p(a, b), [Dict(x => b)]),
                (q(a, a), []),
            ]
        ),
        (
            p(x, x),
            [
                (p(a, a), [Dict(x => a)]),
                (p(a, b), []),
                (q(a, a), []),
            ]
        ),
        (
            p(x, y),
            [
                (p(a, b), [Dict(x => a, y => b), Dict(x => b, y => a)]),
                (p(a, a), [Dict(x => a, y => a)]),
                (q(a, a), []),
            ]
        ),
        (
            p(x, p(a, x)),
            [
                (p(a, p(a, a)), [Dict(x => a)]),
                (p(a, p(a, b)), []),
                (f(a, p(a, a)), []),
                (p(b, p(a, b)), [Dict(x => b)]),
                (p(p(a, b), b), [Dict(x => b)]),
                (p(q(a, b), b), []),
                (p(p(b, b), a), []),
            ]
        ),
        (
            p(x, p(y, x)),
            [
                (p(a, p(a, a)), [Dict(x => a, y => a)]),
                (p(a, p(a, b)), [Dict(x => a, y => b)]),
                (p(b, p(a, c)), []),
            ]
        ),
        (
            p(x, q(y, x)),
            [
                (p(a, q(a, a)), [Dict(x => a, y => a)]),
                (p(q(a, b), a), [Dict(x => a, y => b)]),
                (p(b, q(a, c)), []),
            ]
        ),
        (
            p(x, q(y, z)),
            [
                (p(a, q(a, a)), [Dict(x => a, y => a, z => a)]),
                (p(q(a, b), a), [Dict(x => a, y => a, z => b),
                                 Dict(x => a, y => b, z => a)]),
                (p(q(a, b), q(a, c)), [Dict(x => q(a, b), y => a, z => c),
                                       Dict(x => q(a, b), y => c, z => a),
                                       Dict(x => q(a, c), y => a, z => b),
                                       Dict(x => q(a, c), y => b, z => a)]),
            ]
        ),
        (
            p(q(x, y), q(y, z)),
            [
                (p(q(a, a), q(a, a)), [Dict(x => a, y => a, z => a)]),
                (p(q(a, b), q(c, b)), [Dict(x => a, y => b, z => c),
                                       Dict(x => c, y => b, z => a)]),
            ]
        ),
        (
            p(f(x), q(x, y)),
            [
                (p(f(a), q(a, a)), [Dict(x => a, y => a)]),
                (p(q(a, b), f(a)), [Dict(x => a, y => b)]),
                (p(b, q(a, c)), []),
            ]
        ),
    ],
    "ac" => [
        (
            s(b, x),
            [
                (s(a, b), [Dict(x => a)]),
                (s(a, b, c), [Dict(x => s(a, c))]),
                (s(a, c), []),
            ]
        ),
        (
            s(a, b, c, x),
            [
                (s(a, b, c, b), [Dict(x => b)]),
                (s(a, c, a, b), [Dict(x => a)]),
                (s(c, c, b, a, b), [Dict(x => s(b, c))]),
                (s(a, b, c), []),
                (s(a, b), []),
            ]
        ),
        (
            s(x, y),
            [
                (s(a, b), [Dict(x => a, y => b),
                           Dict(x => b, y => a)]),
                (s(a, b, c), [Dict(x => a, y => s(b, c)),
                              Dict(x => b, y => s(a, c)),
                              Dict(x => c, y => s(a, b)),
                              Dict(y => a, x => s(b, c)),
                              Dict(y => b, x => s(a, c)),
                              Dict(y => c, x => s(a, b))]),
            ]
        )
    ],
    "general" => [
        (
            f(p(a, x), b),
            [
                (f(p(a, a), b), [Dict(x => a)]),
                (f(p(a, b), b), [Dict(x => b)]),
                (f(p(a, b), c), []),
                (f(p(b, b), b), []),
            ]
        ),
        (
            f(p(a, x), g(y)),
            [
                (f(p(a, a), g(b)), [Dict(x => a, y => b)]),
                (f(p(a, b), g(b)), [Dict(x => b, y => b)]),
                (f(p(a, b), h(c)), []),
                (f(p(b, b), g(b)), []),
            ]
        ),
        (
            f(q(p(a, x), g(y)), z),
            [
                (f(q(p(a, a), g(b)), c), [Dict(x => a, y => b, z => c)]),
                (f(q(p(a, b), g(b)), c), [Dict(x => b, y => b, z => c)]),
                (f(q(p(a, b), h(c)), c), []),
                (f(q(p(b, b), g(b)), c), []),
            ]
        ),
        (
            p(x, f(p(y, z), c)),
            [
                (p(a, f(p(b, c), c)), [Dict(x => a, y => b, z => c),
                                               Dict(x => a, y => c, z => b)]),
                (p(a, f(q(b, c), c)), []),
            ]
        ),
        (
            f(x, p(x, y), q(y, z)),
            [
                (f(a, p(a, b), q(b, c)), [Dict(x => a, y => b, z => c)]),
                (f(a, p(b, a), q(b, c)), [Dict(x => a, y => b, z => c)]),
                (f(a, p(b, a), q(c, b)), [Dict(x => a, y => b, z => c)]),
                (f(a, p(a, a), q(c, b)), []),
            ]
        ),
        (
            f(x, q(p(x, y), p(y, z))),
            [
                (f(a, q(p(a, a), p(a, b))), [Dict(x => a, y => a, z => b)]),
                (f(a, q(p(a, b), p(c, b))), [Dict(x => a, y => b, z => c)]),
                (f(b, q(p(a, b), p(c, b))), []),
            ]
        ),
    ],
]

@testset "$set" for (set, tests) ∈ TESTS
    @testset "$pattern" for (pattern, cases) ∈ tests
        @testset "$term" for (term, ms) ∈ cases
            matches = match(pattern, term)
            @test Set(matches) == Set(matches)  # immutability of iterator
            @test Set(matches) == Set(ms)
        end
    end
end


text:
using Base: Bool
function thing_call(f::Function, session::Session, name::AbstractString)
    cm_call(session) do concept_manager
        thing_type_res = TypeDBClient.execute(
            concept_manager, TypeDBClient.ConceptManagerRequestBuilder.get_thing_type_req(name)
        )
        thing_type = TypeDBClient.instantiate(thing_type_res.get_thing_type_res.thing_type)
        remote_thing = TypeDBClient.as_remote(thing_type, concept_manager.transaction)
        f(remote_thing)
    end
end

function cm_call(f::Function, session::Session)
    transaction = Transaction(session, session.sessionID, g.Proto.Transaction_Type.READ, TypeDBOptions())
    concept_manager = ConceptManager(transaction)
    res = f(concept_manager)
    close(transaction)
    return res
end

function _put_attribute_to_db(context, attr_name, type)
    attr_req = g.ConceptManagerRequestBuilder.put_attribute_type_req(attr_name, type)
    return g.execute(context[:transaction], attr_req)
end

@given("put entity type: person") do context
    context[:concept_manager] = ConceptManager(context[:transaction])
    g.put(context[:concept_manager], EntityType, "person")
end

@given("entity(person) set owns attribute type: is-alive") do context
    g.con.entity_set_owns("person", "is-alive", context[:transaction])
end

@given("entity(person) set owns attribute type: age") do context
    g.con.entity_set_owns("person", "age", context[:transaction])
end

@given("entity(person) set owns attribute type: score") do context
    g.con.entity_set_owns("person", "score", context[:transaction])
end

@given("entity(person) set owns attribute type: name") do context
    g.con.entity_set_owns("person", "name", context[:transaction])
end

@given("entity(person) set owns attribute type: email") do context
    g.con.entity_set_owns("person", "email", context[:transaction])
end

@given("entity(person) set owns attribute type: birth-date") do context
    g.con.entity_set_owns("person", "birth-date", context[:transaction])
end

ext: jl


text:
@doc raw"""
    BurresWassertseinMetric <: AbstractMetric

The Bures Wasserstein metric for symmetric positive definite matrices[^MalagoMontruccioPistone2018].

[^MalagoMontruccioPistone2018]:
    > Malagò, L., Montrucchio, L., Pistone, G.:
    > _Wasserstein Riemannian geometry of Gaussian densities_.
    > Information Geometry, 1, pp. 137–179, 2018.
    > doi: [10.1007/s41884-018-0014-4](https://doi.org/10.1007/s41884-018-0014-4)
"""
struct BuresWassersteinMetric <: RiemannianMetric end

@doc raw"""
    change_representer(M::MetricManifold{ℝ,SymmetricPositiveDefinite,BuresWassersteinMetric}, E::EuclideanMetric, p, X)

Given a tangent vector ``X ∈ T_p\mathcal M`` representing a linear function on the tangent
space at `p` with respect to the [`EuclideanMetric`](@ref) `g_E`,
this is turned into the representer with respect to the (default) metric,
the [`BuresWassersteinMetric`](@ref) on the [`SymmetricPositiveDefinite`](@ref) `M`.

To be precise we are looking for ``Z∈T_p\mathcal P(n)`` such that for all ``Y∈T_p\mathcal P(n)```
it holds

```math
⟨X,Y⟩ = \operatorname{tr}(XY) = ⟨Z,Y⟩_{\mathrm{BW}}
```
for all ``Y`` and hence we get
``Z``= 2(A+A^{\mathrm{T}})`` with ``A=Xp``.
"""
change_representer(
    ::MetricManifold{ℝ,SymmetricPositiveDefinite,BuresWassersteinMetric},
    ::EuclideanMetric,
    p,
    X,
)

function change_representer!(
    ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric},
    Y,
    ::EuclideanMetric,
    p,
    X,
)
    A = X * p
    Y .= 2 .* (A .+ A')
    return Y
end

@doc raw"""
    distance(::MatricManifold{SymmetricPositiveDefinite,BuresWassersteinMetric}, p, q)

Compute the distance with respect to the [`BuresWassersteinMetric`](@ref) on [`SymmetricPositiveDefinite`](@ref) matrices, i.e.

```math
d(p,q) =
    \operatorname{tr}(p) + \operatorname{tr}(q) - 2\operatorname{tr}\Bigl( (p^{\frac{1}{2}}qp^{\frac{1}{2}} \bigr)^\frac{1}{2} \Bigr),
```

where the last trace can be simplified (by rotating the matrix products in the trace) to ``\operatorname{tr}(pq)``.
"""
function distance(
    ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric},
    p,
    q,
)
    return sqrt(tr(p) + tr(q) - 2 * tr(sqrt(p * q)))
end

@doc raw"""
    exp(::MatricManifold{ℝ,SymmetricPositiveDefinite,BuresWassersteinMetric}, p, X)

Compute the exponential map on [`SymmetricPositiveDefinite`](@ref) with respect to
the [`BuresWassersteinMetric`](@ref) given by

```math
    \exp_p(X) = p+X+L_p(X)pL_p(X)
```

where ``q=L_p(X)`` denotes the Lyapunov operator, i.e. it solves ``pq + qp = X``.
"""
exp(::MetricManifold{ℝ,SymmetricPositiveDefinite,BuresWassersteinMetric}, p, X)

function exp!(
    ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric},
    q,
    p,
    X,
)
    Y = lyap(p, -X) #lyap solves qp+pq-X=0
    q .= p .+ X .+ Y * p * Y
    return q
end

@doc raw"""
    inner(::MetricManifold{ℝ,SymmetricPositiveDefinite,BuresWassersteinMetric}, p, X, Y)

Compute the inner product [`SymmetricPositiveDefinite`](@ref) with respect to
the [`BuresWassersteinMetric`](@ref) given by

```math
    ⟨X,Y⟩ = \frac{1}{2}\operatorname{tr}(L_p(X)Y)
```

where ``q=L_p(X)`` denotes the Lyapunov operator, i.e. it solves ``pq + qp = X``.
"""
function inner(
    ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric},
    p,
    X,
    Y,
)
    return 1 / 2 * dot(lyap(p, -X), Y)
end

@doc raw"""
    log(::MatricManifold{SymmetricPositiveDefinite,BuresWassersteinMetric}, p, q)

Compute the logarithmic map on [`SymmetricPositiveDefinite`](@ref) with respect to
the [`BuresWassersteinMetric`](@ref) given by

```math
    \log_p(q) = (pq)^{\frac{1}{2}} + (qp)^{\frac{1}{2}} - 2 p
```

where ``q=L_p(X)`` denotes the Lyapunov operator, i.e. it solves ``pq + qp = X``.
"""
log(::MetricManifold{ℝ,SymmetricPositiveDefinite,BuresWassersteinMetric}, p, q)

function log!(
    ::MetricManifold{ℝ,<:SymmetricPositiveDefinite,BuresWassersteinMetric},
    X,
    p,
    q,
)
    X .= sqrt(Symmetric(p * q)) + sqrt(Symmetric(q * p)) - 2 * p
    return X
end


text:
using MarketData

facts("diagnostic tests produce accurate results") do

  context("jaques-berra test") do
      @fact length(cl.values) => 500
  end
end


text: # import WebSockets


text:
##### Precipitation tendencies

#####
##### First order fluxes
#####

function flux(::PrecipitationFlux{Rain}, atmos, args)
    @unpack state, aux = args
    FT = eltype(state)
    u = state.ρu / state.ρ
    q_rai = state.precipitation.ρq_rai / state.ρ

    v_term_rai::FT = FT(0)
    if q_rai > FT(0)
        v_term_rai = terminal_velocity(
            atmos.param_set,
            atmos.param_set.microphys.rai,
            state.ρ,
            q_rai,
        )
    end

    k̂ = vertical_unit_vector(atmos, aux)
    return state.precipitation.ρq_rai * (u - k̂ * v_term_rai)
end

function flux(::PrecipitationFlux{Snow}, atmos, args)
    @unpack state, aux = args
    FT = eltype(state)
    u = state.ρu / state.ρ
    q_sno = state.precipitation.ρq_sno / state.ρ

    v_term_sno::FT = FT(0)
    if q_sno > FT(0)
        v_term_sno = terminal_velocity(
            atmos.param_set,
            atmos.param_set.microphys.sno,
            state.ρ,
            q_sno,
        )
    end

    k̂ = vertical_unit_vector(atmos, aux)
    return state.precipitation.ρq_sno * (u - k̂ * v_term_sno)
end

#####
##### Second order fluxes
#####

function flux(::Diffusion{Rain}, atmos, args)
    @unpack state, aux, t, diffusive = args
    ν, D_t, τ = turbulence_tensors(atmos, state, diffusive, aux, t)
    d_q_rai = (-D_t) .* diffusive.precipitation.∇q_rai
    return d_q_rai * state.ρ
end

function flux(::Diffusion{Snow}, atmos, args)
    @unpack state, aux, t, diffusive = args
    ν, D_t, τ = turbulence_tensors(atmos, state, diffusive, aux, t)
    d_q_sno = (-D_t) .* diffusive.precipitation.∇q_sno
    return d_q_sno * state.ρ
end

#####
##### Sources
#####

function source(s::WarmRain_1M{Rain}, m, args)
    nt = warm_rain_sources(m, args)
    return nt.S_ρ_qr
end

function source(s::RainSnow_1M{Rain}, m, args)
    nt = rain_snow_sources(m, args)
    return nt.S_ρ_qr
end

function source(s::RainSnow_1M{Snow}, m, args)
    nt = rain_snow_sources(m, args)
    return nt.S_ρ_qs
end

max_stars_repo_name: jm-c/CLIMA

max_issues_repo_name: jm-c/CLIMA

max_forks_repo_name: jm-c/CLIMA


text:
@testset "l0_smooth" begin
    @info "Test: L0Smooth"

    @testset "API" begin
        img_gray = imresize(testimage("cameraman"), (64, 64))
        img_rgb = imresize(testimage("lena_color_512"), (64, 64))

        # L0Smooth
        @test L0Smooth() == L0Smooth(λ=2e-2, κ=2.0, βmax=1e5)
        @test L0Smooth(λ=1e-2) == L0Smooth(λ=1e-2, κ=2.0, βmax=1e5)
        @test L0Smooth(κ=1.5) == L0Smooth(λ=2e-2, κ=1.5, βmax=1e5)
        @test L0Smooth(βmax=1e6) == L0Smooth(λ=2e-2, κ=2.0, βmax=1e6)
        @test L0Smooth(λ=1e-2, κ=1.5) == L0Smooth(λ=1e-2, κ=1.5, βmax=1e5)
        @test L0Smooth(βmax=1e6, κ=1.5) == L0Smooth(λ=2e-2, κ=1.5, βmax=1e6)
        @test L0Smooth(βmax=1e6, λ=1e-2) == L0Smooth(λ=1e-2, κ=2.0, βmax=1e6)
        @test L0Smooth(βmax=1e6, λ=1e-2, κ=1.5) == L0Smooth(λ=1e-2, κ=1.5, βmax=1e6)

        # λ should be positive
        @test_throws ArgumentError L0Smooth(λ=-2e-2)
        # κ > 1.0
        @test_throws ArgumentError L0Smooth(κ=0.8)
        # βmax > 1e4
        @test_throws ArgumentError L0Smooth(βmax=1e3)

        # smooth
        f = L0Smooth()
        smoothed_img_gray = smooth(img_gray, f)
        smoothed_img_rgb = smooth(img_rgb, f)
        @test eltype(smoothed_img_gray) == Gray{N0f8}
        @test eltype(smoothed_img_rgb) == RGB{N0f8}
    end

    @testset "ReferenceTests" begin
        # these two reference images are generated from original matlab codes
        # <http://www.cse.cuhk.edu.hk/leojia/projects/L0smoothing/index.html>

        @testset "Gray" begin
        img_gray = imresize(testimage("cameraman"), (64, 64))
        ref = load("algorithms/references/L0_Smooth_Gray.png")
        f = L0Smooth()
        out = smooth(img_gray, f)
        @test assess_psnr(out, eltype(out).(ref)) >= 24
        end

        @testset "RGB" begin
        img_rgb = imresize(testimage("lena_color_512"), (64, 64))
        ref = load("algorithms/references/L0_Smooth_Color3.png")
        f = L0Smooth()
        out = smooth(img_rgb, f)
        @test assess_psnr(out, eltype(out).(ref)) >= 27
        end
    end

    @testset "Offset Arrays" begin
        # Verify that the algorithms give equivalent results on OffsetArrays
        img = imresize(testimage("cameraman"), (64, 64))
        offset_img = OffsetArray(img, OffsetArrays.Origin(0))
        fₛ = L0Smooth()
        imgₛ = smooth(img, fₛ)
        offset_imgₛ = smooth(offset_img, fₛ)
        @test axes(offset_imgₛ) == axes(offset_img)
        @test imgₛ == OffsetArrays.no_offset_view(offset_imgₛ)
    end
end


text:
module rbc
using LinearAlgebra, SymbolicUtils, LaTeXStrings
const max_order = 2
const n_y = 3
const n_x = 2
const n_p = 6
const n_ϵ = 1
const n_z = 2
const η = reshape([0; -1], 2, 1)
const Q = [1.0 0.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0 0.0]
const has_Ω = true
# Display definitions
const x_symbols = [:k, :z]
const y_symbols = [:c, :q, :i]
const u_symbols = [:c, :q, :i, :k, :z]
const p_symbols = [:α, :β, :ρ, :δ, :σ, :Ω_1]
const H_latex = L"\begin{equation}
\left[
\begin{array}{c}
\frac{ - \beta \left( 1 - \delta + \left( k\left( 1 + t \right) \right)^{-1 + \alpha} \alpha e^{z\left( 1 + t \right)} \right)}{c\left( 1 + t \right)} + \frac{1}{c\left( t \right)} \\
 - q\left( t \right) - \left( 1 - \delta \right) k\left( t \right) + c\left( t \right) + k\left( 1 + t \right) \\
 - \left( k\left( t \right) \right)^{\alpha} e^{z\left( t \right)} + q\left( t \right) \\
 - \rho z\left( t \right) + z\left( 1 + t \right) \\
 - k\left( 1 + t \right) + \left( 1 - \delta \right) k\left( t \right) + i\left( t \right) \\
\end{array}
\right]
\end{equation}
"
const steady_states_latex = L"\begin{align}
k\left( \infty \right) =& \left( \frac{-1 + \delta + \frac{1}{\beta}}{\alpha} \right)^{\frac{1}{-1 + \alpha}} \\
z\left( \infty \right) =& 0 \\
c\left( \infty \right) =& \left( \frac{-1 + \delta + \frac{1}{\beta}}{\alpha} \right)^{\frac{\alpha}{-1 + \alpha}} - \left( \frac{-1 + \delta + \frac{1}{\beta}}{\alpha} \right)^{\frac{1}{-1 + \alpha}} \delta \\
q\left( \infty \right) =& \left( \frac{-1 + \delta + \frac{1}{\beta}}{\alpha} \right)^{\frac{\alpha}{-1 + \alpha}} \\
i\left( \infty \right) =& \left( \frac{-1 + \delta + \frac{1}{\beta}}{\alpha} \right)^{\frac{1}{-1 + \alpha}} \delta
\end{align}
"
const steady_states_iv_latex = L"\begin{align}
k\left( \infty \right) =& \left( \frac{-1 + \delta + \frac{1}{\beta}}{\alpha} \right)^{\frac{1}{-1 + \alpha}} \\
z\left( \infty \right) =& 0 \\
c\left( \infty \right) =& \left( \frac{-1 + \delta + \frac{1}{\beta}}{\alpha} \right)^{\frac{\alpha}{-1 + \alpha}} - \left( \frac{-1 + \delta + \frac{1}{\beta}}{\alpha} \right)^{\frac{1}{-1 + \alpha}} \delta \\
q\left( \infty \right) =& \left( \frac{-1 + \delta + \frac{1}{\beta}}{\alpha} \right)^{\frac{\alpha}{-1 + \alpha}} \\
i\left( \infty \right) =& \left( \frac{-1 + \delta + \frac{1}{\beta}}{\alpha} \right)^{\frac{1}{-1 + \alpha}} \delta
\end{align}
"
# Function definitions
include("rbc/zero_order_ip.jl")
include("rbc/first_order_ip.jl")
include("rbc/second_order_ip.jl")
end

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
using DrWatson
@quickactivate "UltrasoundVisualGallery"
DrWatson.greet()


text:
@testset "CartesianIndex" begin
    A = rand(4, 4, 4, 4, 4)
    @test A[.., 1, 2, 3] == A[.., CartesianIndex(1,2,3)]
    @test A[1, .., 2, 3] == A[CartesianIndex(1), .., CartesianIndex(2, 3)]
    @test A[1, 2, .., 3] == A[CartesianIndex(1, 2), .., CartesianIndex(3)]

    @test A[.., 1:1, 2:2, 1:3] == A[.., CartesianIndices((1:1, 2:2, 1:3))]
    @test A[1:1, .., 2:2, 1:3] == A[1:1, .., CartesianIndices((2:2, 1:3))]
    @test A[1:2, 2, .., 1:3] == A[CartesianIndices((1:2,)), 2, .., CartesianIndices((1:3),)]
end


text:
struct TraceObjective{F1, F2} <: Function
    f::F1
    ∇f::F2
    xtrace::Vector{Any}
    ftrace::Vector{Any}
    gtrace::Vector{Any}
end
TraceObjective(f, ∇f) = TraceObjective(f, ∇f, Any[], Any[], Any[])

(to::TraceObjective)(x) = to.f(x)
function ChainRulesCore.rrule(f::TraceObjective, x)
    v, g = f.f(x), f.∇f(x)
    push!(f.xtrace, copy(x))
    push!(f.ftrace, copy(v))
    push!(f.gtrace, copy(g))
    return v, Δ -> (NoTangent(), Δ * g)
end

function maximize_with_trace(f, ∇f, x₀, optimizer::IpoptAlg; power = 2.0, shift = 0.0, deflation_ub = 1000, deflation_radius = 0.0, prev_solutions = Any[], kwargs...)
    _of = TraceObjective(f, ∇f)
    of = x -> -_of(x[1:end-1])
    options = IpoptOptions(; kwargs...)
    model = Nonconvex.Model(of)
    N = length(x₀)
    Nonconvex.addvar!(model, fill(-10000, N), fill(10000, N))
    Nonconvex.addvar!(model, [0.0], [deflation_ub])
    function deflation(x, y)
        if length(prev_solutions) > 0
            d = zero(eltype(x))
            for sol in prev_solutions
                d += (1/max(0, norm(x - sol) - deflation_radius))^power + shift
            end
            return d - y
        else
            return return -one(eltype(x))
        end
    end
    Nonconvex.add_ineq_constraint!(model, x -> deflation(x[1:end-1], x[end]))
    res = Nonconvex.optimize(model, optimizer, [x₀; deflation_ub / 2]; options)
    push!(prev_solutions, res.minimizer[1:end-1])
    return identity.(_of.xtrace), identity.(_of.ftrace), identity.(_of.gtrace)
end

function maximize_with_trace(f, ∇f, x₀, optimizer::Optim.AbstractOptimizer; prev_solutions = [], kwargs...)
    negf(x) = -f(x)
    g!(y, x) = (y .= .-∇f(x))

    function callback(states)
        # terminate if optimization encounters NaNs
        s = states[end]
        md = s.metadata
        return isnan(s.value) || any(isnan, md["x"]) || any(isnan, md["g(x)"])
    end
    options = Optim.Options(;
        store_trace=true, extended_trace=true, callback=callback, kwargs...
    )
    res = Optim.optimize(negf, g!, x₀, optimizer, options)

    xs = Optim.x_trace(res)::Vector{typeof(Optim.minimizer(res))}
    fxs = -Optim.f_trace(res)
    ∇fxs = map(tr -> -tr.metadata["g(x)"], Optim.trace(res))::typeof(xs)

    return xs, fxs, ∇fxs
end


text:
using ResumableFunctions
using Test

@resumable function test_for(a::Int=0; b::Int=a+1) :: Int
  for i in 1:10
    @yield a
    a, b = b, a+b
  end
end

@testset "test_for" begin
@test collect(test_for(4)) == [4, 5, 9, 14, 23, 37, 60, 97, 157, 254]
end

@resumable function test_try(io)
  try
    a = 1
    @yield a
    a = 2
    c = @yield a
    println(io,c)
  catch except
    println(io,except)
  finally
    println(io,"Always")
  end
end

struct SpecialException <: Exception end
@testset "test_try" begin
io = IOBuffer()
try_me = test_try(io)
try_me()
try_me(SpecialException())
@test String(take!(copy(io))) == "SpecialException()\nAlways\n"

io = IOBuffer()
try_me = test_try(io)
try_me()
try_me()
try_me("Hello")
@test String(take!(copy(io))) == "Hello\nAlways\n"

io = IOBuffer()
try_me = test_try(io)
try_me()
try_me()
try_me(SpecialException())
@test_throws ErrorException try_me()
@test String(take!(copy(io))) == "SpecialException()\nAlways\n"
end

@resumable function (test_where1(a::N) :: N) where {N<:Number}
  b = a + one(N)
  for i in 1:10
    @yield a
    a, b = b, a+b
  end
end

@resumable function (test_where2(a::N=4; b::N=a + one(N)) :: N) where N
  for i in 1:10
    @yield a
    a, b = b, a+b
  end
end

@testset "test_where" begin
@test collect(test_where1(4.0)) == [4.0, 5.0, 9.0, 14.0, 23.0, 37.0, 60.0, 97.0, 157.0, 254.0]
@test collect(test_where2(4)) == [4, 5, 9, 14, 23, 37, 60, 97, 157, 254]
end

@resumable function test_varargs(a...)
  for (i, e) in enumerate(a)
    @yield e
  end
end

@testset "test_varargs" begin
@test collect(test_varargs(1, 2, 3)) == [1, 2, 3]
end

@resumable function test_let()
  for u in [[(1,2),(3,4)], [(5,6),(7,8)]]
    for i in 1:2
      let i=i
        val = [a[i] for a in u]
      end
      @yield val
    end
  end
end

@testset "test_let" begin
@test collect(test_let()) == [[1,3],[2,4],[5,7],[6,8]]
end

@resumable function test_nosave()
  for i in 65:74
    @nosave tmp = Char(i)
    @yield tmp
  end
end

@testset "test_nosave" begin
@test collect(test_nosave()) == ['A', 'B', 'C', 'D', 'E',  'F', 'G', 'H', 'I', 'J']
end

@resumable function test_return_value()
  return 1
end

@testset "test_return_value" begin
  @test collect(test_return_value()) == []
end

@resumable function test_continue()
  for i in 1:10
    if i === 2
      continue
    end
    @yield i
  end
end

@resumable function test_continue_double()
  for i in 1:3
    if i === 2
      continue
    end
    for j in 1:3
      if j === 2
        continue
      end
      @yield j
    end
  end
end

@testset "test_continue" begin
  @test collect(test_continue()) == [1, 3, 4, 5, 6, 7, 8, 9, 10]
  @test collect(test_continue_double()) == [1, 3, 1, 3]
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule MAGEMin_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("MAGEMin")
JLLWrappers.@generate_main_file("MAGEMin", UUID("763ebaa8-b0d2-5f6b-90ef-4fc23b5db1c4"))
end  # module MAGEMin_jll


text:
module FinDiff2018

#push!(LOAD_PATH, @__DIR__) # expose all other modules defined in this directory.

#using Plots

export PLCFun, derivRight, samplePLC, evaluate, refine

############################################################
##			Piecewise Linear Continuous Functions

# We approximate functions by a piecewise linear continuous approximation,
# with a regular grid spacing
############################################################

# A piecewise linear continuous function mapping a type T to a type U
#export PLCFun
struct PLCFun{T,U}
    # Domain: [0; 1]
    points::Vector{U}
end

# Functions can be scaled and added
# function Base.length(f::PLCFun{T, U})::Int where {T, U}
#     length(f.points)
# end
function Base. *(a::U, f::PLCFun{T, U})::PLCFun{T, U} where {T, U}
    PLCFun{T, U}(a .* f.points)
end
function Base. +(f::PLCFun{T, U}, g::PLCFun{T, U})::PLCFun{T, U} where {T, U}
    @assert length(f.points) == length(g.points)
    PLCFun{T, U}(f.points .+ g.points)
end

############################################################
# Calculate the x coordinates of the endpoints of the lines
#export xcoord
function xcoord(::Type{T}, nlines::Int, i::Int)::T where {T}
    @assert 1 <= i <= nlines + 1
    dx = 1 / nlines
    x = (i-1) * dx
    x
end

# The inverse of "xcoord": Determine the line segment on which a particular x coordinate lies
#export lineidx
function lineidx(f::PLCFun{T, U}, x::T)::Int where {T, U}
    @assert 0 <= x <= 1
    nlines = length(f.points) - 1
    dx = 1 / nlines
    i = floor(Int, x / dx) + 1
    i = max(1, i)
    i = min(nlines, i)
    i
end

# Convert a general Julia function into a PLCFun. We need to specify the types T and U,
# as well as the number of line segments to use.
#export samplePLC
function samplePLC(::Type{T}, ::Type{U}, nlines::Int, f::Function)::PLCFun{T, U} where {T, U}
    ys = U[f(xcoord(T, nlines, i)) for i in 1:nlines+1]
    PLCFun{T, U}(ys)
end

# Linear interpolation between two points
#export linterp
function linterp(x1::T, y1::U, x2::T, y2::U, x::T)::U where {T, U}
    y1 * (x - x2) / (x1 - x2) + y2 * (x - x1) / (x2 - x1)
end

# Evaluate a PLCFun at point x
#export evaluate
function evaluate(f::PLCFun{T,U}, x::T)::U where {T, U}
    @assert 0 <= x <= 1
    nlines = length(f.points) - 1
    # Find out on which line segment the point lies
    i = lineidx(f, x)
    # Interpolate between the two endpoints of the line segment
    x1 = xcoord(T, nlines, i)
    x2 = xcoord(T, nlines, i+1)
    y1 = f.points[i]
    y2 = f.points[i+1]
    y = linterp(x1, y1, x2, y2, x)
    y
end

# Calculate the derivative of a PLCFun, using a right-biased derivative
#export derivRight
function derivRight(f::PLCFun{T, U})::PLCFun{T, U} where {T, U}
    nlines = length(f.points) - 1
    dx = 1 / nlines
    ys = [[(f.points[i+1] - f.points[i]) / dx for i in 1:nlines];
          (f.points[end] - f.points[end-1]) / dx]
    PLCFun{T, U}(ys)
end

#export refine
function refine(f::PLCFun{T, U})::PLCFun{T, U} where {T, U}
    nlines = length(f.points) - 1
    dx = 1 / nlines
    ys = [evaluate(f, xcoord(T, 2*nlines, i)) for i in 1:2*nlines+1]
    PLCFun{T, U}(ys)
end

#=

Create Plots
#export testplot
#function testplot(arr::Array{Float64,2}, filename::String)
#	heatmap(arr, clim=(-1.0, +1.0), color=:viridis)
#	savefig(filename)
#end

=#

include("AdvectionEq.jl")
include("TimeEvolve.jl")

end # module FinDiff2018


text:
mutable struct MobileUser
	id::Int
	starting_point::Tuple{Int64,Int64}
	canvas::Canvas
	alpha::Int
	y::Float16
end

MobileUser(id::Int,canvas::Canvas,start_point::Tuple{Int64,Int64},alpha::Int,y::Float16) = MobileUser(id,start_point,canvas,alpha,y)

function move(m::MobileUser,iterations::Int64)

	T = DateTime(0)
	logfile = open("log_$(m.id).csv","w+")
	write(logfile,"ID;cluster;node;time\n")

	for iteration in 1:iterations
		#Select cluster proportional to the size
		cluster_index, cluster = select_cluster(m.canvas.clusters)
		println("Select cluster $cluster_index with $(length(cluster.points))")
		#Select y% subset of cluster waypoints
		num_waypoints = ceil(Int64,length(cluster.points) * m.y)
		println("Num waypoints selected $num_waypoints")
		#Extract subpoint list from the cluster
		vertex_list = random_elements_from_array(collect(1:length(cluster.points)),num_waypoints)
		# println("vertex selected: $vertex_list")
		point_list = [cluster.points[i] for i in vertex_list]
		#Add the mobile node to the graph
		push!(point_list,m.starting_point)
		#Generate subgraph
		subgraph = generate_graph(point_list)
		#Visit node with LATP algorithm, starting from the node position
		visiting_order =  latp_algorithm(subgraph, nv(subgraph), m.alpha,T)
		last_node_visited = visiting_order[length(visiting_order)][1]
		m.starting_point =	point_list[last_node_visited]

		for node in visiting_order
			write(logfile,"$(m.id);$cluster_index;$(node[1]);$(node[2])\n")
		end

	end
	close(logfile)
end

#TODO: This would be more efficient if the edge weights are generated first and then the Graph is constructed, according to SimpleWeightedGraph doc
function generate_graph(point_list::Array{Tuple{Int64,Int64},1})
	# println(point_list)
	g = SimpleWeightedGraph()
	for j in 1:length(point_list)
		add_vertex!(g)
		for i in 1:nv(g)-1
			add_edge!( g, i, nv(g), node_distance(point_list[i],point_list[nv(g)]))
		end
	end
	# println(g)
	return g
end

function select_cluster(clusters::Vector{Cluster})
		cluster_ordered = sort([(length(clusters[i].points),i) for i in 1:length(clusters)],rev=true)
		prob_vector = map(d->d[1],cluster_ordered)
		prob_sum = sum( prob_vector )
		prob_vector = map( d -> d / prob_sum, prob_vector)
		#Find node: the probability of selecting a node i is the sum of the probability of all preceiding probs
		next = -1
		p = rand()
		for i in 1:length(cluster_ordered)
			if p <= sum(prob_vector[1:i])
				next = cluster_ordered[i][2]
				break
			end
		end
		return next, clusters[next]
end


text:
function null_error_handler(error_code::Cint,
                  mod::Char,
                  func::Char,
                  eh_data::Ptr{Void})
  nothing
end


text:
using Distributed
addprocs()
using ProgressMeter

@everywhere function get_hundreds_digit(x::Int64)::Int64
    _x = floor(x/100.0)
    return _x % 10
end

@everywhere function power_at_coordinates(x::Int64,y::Int64,serial_number::Int64)::Int64
    rack_id = x+10
    power_level = rack_id*y
    power_level += serial_number
    power_level *= rack_id
    power_level = get_hundreds_digit(power_level)
    power_level -= 5
    return power_level
end

@everywhere function power_cell_grid(serial_number::Int64)::Array{Int64,2}
    grid = zeros(300,300)
    for x in 1:300
        for y in 1:300
            grid[y,x] = power_at_coordinates(x,y,serial_number)
        end
    end
    return grid
end

@everywhere function find_largest_total_power(serial_number::Int64,square_size::Int64)::Pair{Tuple{Int64,Int64},Int64}
    grid = power_cell_grid(serial_number)
    subgrid_power = Dict()
    for x in 1:301-square_size
        for y in 1:301-square_size
            subgrid = grid[y:y+square_size-1,x:x+square_size-1]
            subgrid_power[(x,y)] = sum(subgrid)
        end
    end
    most_power = sort(collect(subgrid_power), by=x->-x[2])
    return most_power[1]
end

function part_2(serial_number::Int64)
    high_scores = Dict()
    futures = []
    @showprogress for square_size in 1:300
        future = @spawn find_largest_total_power(serial_number,square_size)
        push!(futures, (square_size, future))
    end
    @showprogress for f in futures
        high_scores[f[1]] = fetch(f)
    end
    high_score = sort(collect(high_scores), by=x->-x[2][2][2])
    return high_score[1]
end

part_1 = find_largest_total_power(18,3)
@time part_2(18)


text:
const _deprecated_attributes = Dict{Symbol,Symbol}(:orientation => :permute)
const _all_defaults = KW[_series_defaults, _plot_defaults, _subplot_defaults]

const _initial_defaults = deepcopy(_all_defaults)
const _initial_axis_defaults = deepcopy(_axis_defaults)

# add defaults for the letter versions
const _axis_defaults_byletter = KW()

function reset_axis_defaults_byletter!()
    for letter in (:x, :y, :z)
        _axis_defaults_byletter[letter] = KW()
        for (k, v) in _axis_defaults
            _axis_defaults_byletter[letter][k] = v
        end
    end
end
reset_axis_defaults_byletter!()

# to be able to reset font sizes to initial values
const _initial_plt_fontsizes =
    Dict(:plot_titlefontsize => _plot_defaults[:plot_titlefontsize])

const _initial_sp_fontsizes = Dict(
    :titlefontsize => _subplot_defaults[:titlefontsize],
    :legend_font_pointsize => _subplot_defaults[:legend_font_pointsize],
    :legend_title_font_pointsize => _subplot_defaults[:legend_title_font_pointsize],
    :annotationfontsize => _subplot_defaults[:annotationfontsize],
    :colorbar_tickfontsize => _subplot_defaults[:colorbar_tickfontsize],
    :colorbar_titlefontsize => _subplot_defaults[:colorbar_titlefontsize],
)

const _initial_ax_fontsizes = Dict(
    :tickfontsize  => _axis_defaults[:tickfontsize],
    :guidefontsize => _axis_defaults[:guidefontsize],
)

const _initial_fontsizes =
    merge(_initial_plt_fontsizes, _initial_sp_fontsizes, _initial_ax_fontsizes)

const _internal_args =
    [:plot_object, :series_plotindex, :markershape_to_add, :letter, :idxfilter]

const _axis_args = Set(keys(_axis_defaults))
const _series_args = Set(keys(_series_defaults))
const _subplot_args = Set(keys(_subplot_defaults))
const _plot_args = Set(keys(_plot_defaults))

const _magic_axis_args = [:axis, :tickfont, :guidefont, :grid, :minorgrid]
const _magic_subplot_args =
    [:title_font, :legend_font, :legend_title_font, :plot_title_font, :colorbar_titlefont]
const _magic_series_args = [:line, :marker, :fill]
const _all_magic_args =
    Set(union(_magic_axis_args, _magic_series_args, _magic_subplot_args))

const _all_axis_args = union(_axis_args, _magic_axis_args)
const _lettered_all_axis_args =
    Set([Symbol(letter, kw) for letter in (:x, :y, :z) for kw in _all_axis_args])
const _all_subplot_args = union(_subplot_args, _magic_subplot_args)
const _all_series_args = union(_series_args, _magic_series_args)
const _all_plot_args = _plot_args

const _all_args =
    union(_lettered_all_axis_args, _all_subplot_args, _all_series_args, _all_plot_args)

# add all pluralized forms to the _keyAliases dict
for arg in _all_args
    add_aliases(arg, makeplural(arg))
end

# fill symbol cache
for letter in (:x, :y, :z)
    _attrsymbolcache[letter] = Dict{Symbol,Symbol}()
    for k in _axis_args
        # populate attribute cache
        lk = Symbol(letter, k)
        _attrsymbolcache[letter][k] = lk
        # allow the underscore version too: xguide or x_guide
        add_aliases(lk, Symbol(letter, "_", k))
    end
    for k in (_magic_axis_args..., :(_discrete_indices))
        _attrsymbolcache[letter][k] = Symbol(letter, k)
    end
end

# add all non_underscored forms to the _keyAliases
add_non_underscore_aliases!(_keyAliases)


text:
"""
    labeled_accuracy(prediction, gold)

Accuracy score for dependency arcs, including the labels.
"""
function labeled_accuracy end

"""
    unlabeled_accuracy(prediction, gold)

Accuracy score for dependency arcs, not including the labels.
"""
function unlabeled_accuracy end

const labelled_accuracy = labeled_accuracy
const unlabelled_accuracy = unlabeled_accuracy

function labeled_accuracy(tree::DependencyTree, gold::DependencyTree)
    results = map(zip(tokens(tree), tokens(gold))) do (t, g)
        head(t) == head(g) && deprel(t) == deprel(g)
    end
    count(results) / length(results)
end

function labeled_accuracy(trees, goldtrees)
    results = Bool[]
    for (tree, gold) in zip(trees, goldtrees)
        for (t, g) in zip(tokens(tree), tokens(gold))
            append!(results, head(t) == head(g) && deprel(t) == deprel(g))
        end
    end
    count(results) / length(results)
end

function unlabeled_accuracy(tree::DependencyTree, gold::DependencyTree)
    results = map(zip(tokens(tree), tokens(gold))) do (t, g)
        head(t) == head(g)
    end
    count(results) / length(results)
end

function unlabeled_accuracy(trees, goldtrees)
    results = Bool[]
    for (tree, gold) in zip(trees, goldtrees)
        for (t, g) in zip(tokens(tree), tokens(gold))
            push!(results, head(t) == head(g))
        end
    end
    count(results) / length(results)
end


text:
#
# Evaluates s-plane integrand for u (and possibly ux,uy) given params a,b.
#

function integrand!(u::Vector{Complex{Float64}}, ux::Vector{Complex{Float64}}, uy::Vector{Complex{Float64}}, s::Vector{Complex{Float64}}, a::Float64, b::Float64, x::Float64, y::Float64, derivs::Bool, n::Int)
    @simd for i=1:n
        # Precalculations
        @inbounds es = exp(s[i])
        ies = inv(es)
        @inbounds u[i] = cis(a * ies + b * es - es * es * es / 12.0)

        if derivs
            # calculate prefactors for exponentials in order to do derivatives
            c,d = 0.5im*ies,0.5im*es
            # x deriv computation
            @inbounds ux[i] = (-c*x)*u[i]
            # y deriv computation
            @inbounds uy[i] = (-c*y+d)*u[i]
        else
            # make sure that ux, uy aren't going to give us anything funny (or >emach)
            @inbounds ux[i] = ZIM
            @inbounds uy[i] = ZIM
        end
    end
end

# assume if length n not given that vectors are of length 1.

integrand!(u::Vector{Complex{Float64}}, ux::Vector{Complex{Float64}}, uy::Vector{Complex{Float64}}, s::Vector{Complex{Float64}}, a::Float64, b::Float64, x::Float64, y::Float64, derivs::Bool) = integrand!(u, ux, uy, s, a, b, x, y, derivs, 1)


text:
#! /usr/bin/julia

# Rosetta Code, Find largest left truncatable prime in a given base

function addmsdigit{T<:Integer}(p::T, b::T, s::T)
    a = T[]
    q = p
    for i in 1:(b-1)
        q += s
        isprime(q) || continue
        push!(a, q)
    end
    return a
end

function lefttruncprime{T<:Integer}(pbase::T)
    b = convert(BigInt, pbase)
    a = BigInt[]
    append!(a, primes(b-1))
    mlt = zero(BigInt)
    s = one(BigInt)
    while !isempty(a)
        mlt = maximum(a)
        s *= b
        for i in 1:length(a)
            p = shift!(a)
            append!(a, addmsdigit(p, b, s))
        end
    end
    return mlt
end

lo, hi = 3, 17
print("The largest left truncatable primes for bases")
println(@sprintf " %d to %d." lo hi)
for i in lo:hi
    mlt = lefttruncprime(i)
    println(@sprintf "   %3d %d (%s)" i mlt base(i, mlt))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
"""
# Conversion of non-Index data to Matrix

Data in non-index columns of a TS object can be converted into a
`Matrix` type for further numerical analysis using the `Matrix()`
constructor.

# Examples

```jldoctest; setup = :(using TSx, DataFrames, Dates, Random, Statistics)
julia> using Random;
julia> random(x) = rand(MersenneTwister(123), x);
julia> ts = TS([random(10) random(10)])
julia> show(ts)
(10 x 2) TS with Int64 Index

 Index  x1         x2
 Int64  Float64    Float64
─────────────────────────────
     1  0.768448   0.768448
     2  0.940515   0.940515
     3  0.673959   0.673959
     4  0.395453   0.395453
     5  0.313244   0.313244
     6  0.662555   0.662555
     7  0.586022   0.586022
     8  0.0521332  0.0521332
     9  0.26864    0.26864
    10  0.108871   0.108871

julia> Matrix(ts)
10×2 Matrix{Float64}:
 0.768448   0.768448
 0.940515   0.940515
 0.673959   0.673959
 0.395453   0.395453
 0.313244   0.313244
 0.662555   0.662555
 0.586022   0.586022
 0.0521332  0.0521332
 0.26864    0.26864
 0.108871   0.108871
```
"""
function Base.Matrix(ts::TS)
    DataFrames.Matrix(ts.coredata[!, Not(:Index)])
end


text:
"""
    observe(dist::Distribution, x)

Observe sample `x` from distribution `dist` and yield its log-likelihood value.
"""
function observe(dist::Distributions.Distribution, x)
    return Libtask.produce(Distributions.loglikelihood(dist, x))
end


text:
using Sockets
using Distributed: Future

"""
A map from request_id's to pending conditions.
"""
pending_requests = Dict{String, Future}()

"""
    send_request(scope, request[, "key" => "value", ...])

Send a request message for a scope and wait for the response.

Each request has a unique id. This is generated automatically. The client should
(eventually) return a response message with the same id.

**IMPORTANT:** Julia code **CANNOT** synchronously `fetch` (or `wait`, etc.) a
    request that was sent while using IJulia (e.g. Jupyter Notebook/Lab); if you
    exclusively use another provider, it's (probably) fine. The reason for this
    is that the `fetch` makes the current task block, but IJulia will wait for
    the current task to finish before executing the next task. The response will
    be handled in a subsequent task (via a comm message), but since this task
    can never be processed (as the event loop is blocked), this creates a
    deadlock which will freeze the kernel indefinitely.
"""
function send_request(conn, request, data::Pair...)
    request_id = string(UUIDs.uuid1())
    if haskey(pending_requests, request_id)
        error("Cannot register duplicate request id: $(request_id).")
    end
    future = Future()
    pending_requests[request_id] = future
    message = Dict(
      "type" => "request",
      "request" => request,
      "requestId" => request_id,
      data...
    )
    send(conn, message)

    # Run in separate async task to ensure that we delete the pending request
    # future when we get a response (this avoids orphaning responses that we
    # don't actually need).
    return @async begin
        try
            return fetch(future)
        catch exc
            @error "Error fetching request future." exception=exc
            rethrow()
        finally
            delete!(pending_requests, request_id)
        end
    end
end

function Sockets.send(c::AbstractConnection, msg)
    error("No send method for connection of type $(typeof(c))")
end

function logmsg(msg, level="info", data=nothing)
    if level == "error" || level == "warn"
        @warn(msg)
    else
        @info(msg)
    end

    Dict("type"=>"log", "message"=>msg, "level"=>level, data=>data)
end

function log(c::AbstractConnection, msg, level="info", data=nothing)
    send(c, logmsg(msg, level, data))
end

function dispatch(conn::AbstractConnection, data)
    message_type = data["type"]
    if message_type == "request"
        return dispatch_request(conn, data)
    elseif message_type == "response"
        return dispatch_response(conn, data)
    elseif message_type == "command" || haskey(data, "command")
        return dispatch_command(conn, data)
    end
    @error "Unknown WebIO message type: $(message_type)."
end

# function dispatch_command(conn::AbstractConnection, data)
function dispatch_command(conn::AbstractConnection, data)
    # first, check if the message is one of the administrative ones
    cmd = data["command"]
    scope = lookup_scope(data["scope"])
    if cmd == "setup_scope" || cmd == "_setup_scope"
        if cmd == "_setup_scope"
            @warn("Client used deprecated command: _setup_scope.", maxlog=1)
        end
        addconnection!(scope.pool, conn)
    elseif cmd == "update_observable"
        if !haskey(data, "name")
            @error "update_observable message missing \"name\" key."
            return
        elseif !haskey(data, "value")
            @error "update_observable message missing \"value\" key."
            return
        end
        dispatch(scope, data["name"], data["value"])
    else
        @warn "Implicit observable update command is deprecated."
        dispatch(scope, cmd, data["data"])
    end
end

ResponseDict = Dict{String, Any}
request_handlers = Dict{String, Function}()
function register_request_handler(request_type::String, handler::Function)
    if haskey(request_handlers, request_type)
        error("Duplicate request type: $(request_type).")
    end
    request_handlers[request_type] = handler
end

function dispatch_request(conn::AbstractConnection, data)
    request_id = get(data, "requestId", nothing)
    request_type = get(data, "request", nothing)
    if request_id === nothing
        @error("Request message (request=$(repr(request_type))) is missing requestId.")
        return
    end
    try
        handler = get(request_handlers, request_type, nothing)
        if handler === nothing
            error("Unknown request type (request=$(repr(request_type))).")
        end
        # Julia sometimes narrows the type of the returned dict
        # (for example, the handler might return Dict{String, Int64}).
        response = convert(ResponseDict, handler(data))
        response["type"] = "response"
        response["request"] = request_type
        response["requestId"] = request_id
        send(conn, response)
        return
    catch (e)
        send(conn, Dict(
            "type" => "response",
            "request" => get(data, "request", nothing),
            "requestId" => request_id,
            "error" => sprint(showerror, e),
        ))
        return
    end
end

function dispatch_response(conn::AbstractConnection, data)
    request_id = get(data, "requestId", nothing)
    if request_id === nothing
        @error "Response message is missing `requestId` key."
        return
    end

    future = get(pending_requests, request_id, nothing)
    if future === nothing
        @error "Received response message for unknown requestId: $(request_id)."
        return
    end
    put!(future, data)
end


text:
macro meta(ex::Expr)
    @show ex.head ex.args
    ret = Expr(:block)
    if ex.head == :import
        name = string(ex.args[1])
        spec = get_spec(name)
        @show spec
        for name in spec.tagged
            push!(ret.args, :(import $name))
        end
        for pkg in spec.branch
            push!(ret.args, :(import $(pkg.name)))
        end
    end
    esc(ret)
end


text:
include("../src/geom_module.jl")

module dtest
using ..geom
using PyPlot
PyPlot.pygui(true)

# x_ = vec([0 1 0.3 1.5 1.7 -0.9 -1])
# y_ = vec([ 0 0.2 1 1.5 3 2.55 0.65])

x_ = [0.0, 1, 2, 3, 2, 1, 0]
y_ = [0.0, 0, 1, 1, 3, 1, 1]

point_array = map( (x, y) -> geom.MutablePoint(x, y), x_, y_ )

delmesh = geom.build_delaunay_mesh!(point_array)

triidx = geom.triangle_indices(delmesh)

PyPlot.figure()
PyPlot.grid()
geom.plot_triangulation(delmesh; init_triangle = false)

end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using HTTP
using JSON
using WAV
using Images
using SignalAnalysis

function read_audio(filepath)
	val, fs = wavread(filepath)
	lengthneed= 1600
	y = SignalAnalysis.power(tfd(x, Spectrogram(nfft=128+64)))
        y = imresize(Gray.(y), (28, 28));
        y = reshape(convert(Array{Float32},y ), (28,28,1));
	return y
end
y = read_audio("./../mini_speech_commands/down/0a9f9af7_nohash_0.wav")
# y = reshape(convert(Array{Float64},y ), (28,28));
# y = round.(y; digits=5)
# y = [[1.0,2.0,3,4,5],[6,7,8,9,10]]
# y = collect(Iterators.flatten(y))
# println(length(y))
y = string(y)
# println(y)
# data = [1,2,3,4]
# url = "http://localhost:8001/digitreg"
# params = Dict("user"=>"j")

# send the request
# r = HTTP.request("POST", url, [("Content-Type", "application/json")], JSON.json(params))

# query_dict = Dict("retmax" => [1.0,2.0])
# body = HTTP.escapeuri(query_dict)
# HTTP.request("POST", "http://localhost:8001/digitreg", [], body=body)
r = HTTP.post("http://localhost:8001/digitreg", [], y)
println(String(r.body))
# HTTP.post("http://localhost:8001/digitreg", [], body)

max_stars_repo_name: ashwani-rathee/Sardor

max_issues_repo_name: ashwani-rathee/Sardor

max_forks_repo_name: ashwani-rathee/Sardor


text:
module CUDAdrv

using CUDAapi

using CEnum

using Printf

## source code includes

# essential functionality
include("pointer.jl")
const CUdeviceptr = CuPtr{Cvoid}

# low-level wrappers
include("libcuda_common.jl")
include("error.jl")
include("libcuda.jl")
include("libcuda_aliases.jl")

# high-level wrappers
include("types.jl")
include("version.jl")
include("devices.jl")
include("context.jl")
include(joinpath("context", "primary.jl"))
include("stream.jl")
include("memory.jl")
include("module.jl")
include("events.jl")
include("execution.jl")
include("profile.jl")
include("occupancy.jl")

include("deprecated.jl")

## deferred initialization

# CUDA packages require complex initialization (discover CUDA, download artifacts, etc)
# that can't happen at module load time, so defer that to run time upon actual use.

const configured = Ref{Union{Nothing,Bool}}(nothing)

"""
    functional(show_reason=false)

Check if the package has been configured successfully and is ready to use.

This call is intended for packages that support conditionally using an available GPU. If you
fail to check whether CUDA is functional, actual use of functionality might warn and error.
"""
function functional(show_reason::Bool=false)
    if configured[] === nothing
        _functional(show_reason)
    end
    configured[]::Bool
end

const configure_lock = ReentrantLock()
@noinline function _functional(show_reason::Bool=false)
    lock(configure_lock) do
        if configured[] === nothing
            configured[] = false
            if __configure__(show_reason)
                configured[] = true
                try
                    __runtime_init__()
                catch
                    configured[] = false
                    rethrow()
                end
            end
        end
    end
end

# macro to guard code that only can run after the package has successfully initialized
macro after_init(ex)
    quote
        @assert functional(true) "CUDAdrv.jl did not successfully initialize, and is not usable."
        $(esc(ex))
    end
end

## initialization

const __libcuda = Sys.iswindows() ? :nvcuda : :libcuda
libcuda() = @after_init(__libcuda)

function __configure__(show_reason::Bool)
    if haskey(ENV, "_") && basename(ENV["_"]) == "rr"
        show_reason && @error("Running under rr, which is incompatible with CUDA")
        return false
    end

    try
        @debug "Initializing CUDA driver"
        res = @runtime_ccall((:cuInit, __libcuda), CUresult, (UInt32,), 0)
        if res == 0xffffffff
            error("Cannot use the CUDA stub libraries. You either don't have the NVIDIA driver installed, or it is not properly discoverable.")
        elseif res != SUCCESS
            throw_api_error(res)
        end
    catch ex
        show_reason && @error("Could not initialize CUDA", exception=(ex,catch_backtrace()))
        return false
    end

    return true
end

function __runtime_init__()
    if version() < v"9"
        @warn "CUDAdrv.jl only supports NVIDIA drivers for CUDA 9.0 or higher (yours is for CUDA $(version()))"
    end
end

end


text:
module SpecialDirectionCodesModule

export coord_base_indicator_ocean,coord_base_indicator_outflow
export coord_base_indicator_truesink,coord_base_indicator_lake
export dir_based_indicator_ocean,dir_based_indicator_outflow
export dir_based_indicator_truesink,dir_based_indicator_lake

const coord_base_indicator_ocean    = -2
const coord_base_indicator_outflow  = -1
const coord_base_indicator_truesink = -5
const coord_base_indicator_lake     = -4

const dir_based_indicator_ocean =   -1
const dir_based_indicator_outflow =  0
const dir_based_indicator_truesink = 5
const dir_based_indicator_lake =    -2

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
"""
function nJacFD{T}(fcn, n::Int64, lda::Int64, x::Vector{T}, fx::Vector{T},
    yscal::Vector{T}, ajdel::T, ajmin::T, nFcn::Int64,
    a::Array{T,2})

Evaluation of a dense Jacobian matrix using finite difference approximation
adapted for use in nonlinear systems solver.

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable | Description                                                     |
|:---------|:----------------------------------------------------------------|
| fcn      | Function of the form fcn(f, x) to provide right-hand side       |
| n        | Number of rows and columns of the Jacobian                      |
| lda      | Leading dimension of the array "a"                              |
| x[n]     | Current scaled vector                                           |
| fx[n]    | Vector containing fcn(x)                                        |
| yscal[n] | Vector containing scaling factors                               |
| ajdel    | Perturbation of component k: abs(y(k))*ajdel                    |
| ajmin    | Minimum perturbation is ajmin*ajdel                             |
| nFcn*    | fcn evaluation count                                            |

(* marks inout parameters)

## Output parameters
-------------------
| Variable | Description                                                     |
|:---------|:----------------------------------------------------------------|
| a[lda,n] | Array containing the approximated Jacobian                      |
| nFcn*    | fcn evaluation count adjusted                                   |
| iFail    | Return code non-zero if Jacobian could not be computed          |
"""
function nJacFD{T}(fcn, n::Int64, lda::Int64, x::Vector{T}, fx::Vector{T},
    yscal::Vector{T}, ajdel::T, ajmin::T, nFcn::Int64,
    a::Array{T,2})
    # Begin
    # Copy for internal purposes
    xa = x[:]
    iFail = 0
    fu = zero(x)
    for k = 1:n
        w = xa[k]
        su = sign(xa[k])
        if su == 0
            su = 1
        end
        u = max(max(abs(xa[k]),ajmin),yscal[k])*ajdel*su;
        xa[k] = w + u
        try
            fcn(fu,xa)
        catch
            iFail = -1
        end
        nFcn += 1
        if iFail != 0
            break
        end
        xa[k] = w
        a[1:n,k:k] = (fu-fx)./u
    end
    return (nFcn, iFail)
end

"""
function nJacFDb{T}(fcn, n::Int64, lda::Int64, ml::Int64, x::Vector{T},
    fx::Vector{T}, yscal::Vector{T}, ajdel::T, ajmin::T,
    nFcn::Int64, a::Array{T,2})

Evaluation of a banded Jacobian matrix using finite difference approximation
adapted for use in nonlinear systems solver

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable | Description                                                     |
|:---------|:----------------------------------------------------------------|
| fcn      | Function of the form fcn(f, x) to provide right-hand side       |
| n        | Number of rows and columns of the Jacobian                      |
| lda      | Leading dimension of the array "a"                              |
| ml       | Lower bandwidth of the Jacobian matrix                          |
| x[n]     | Current scaled vector                                           |
| fx[n]    | Vector containing fcn(x)                                        |
| yscal[n] | Vector containing scaling factors                               |
| ajdel    | Perturbation of component k: abs(y(k))*ajdel                    |
| ajmin    | Minimum perturbation is ajmin*ajdel                             |
| nFcn*    | fcn evaluation count                                            |

(* marks inout parameters)

## Output parameters
-------------------
| Variable | Description                                                     |
|:---------|:----------------------------------------------------------------|
| a[lda,n] | Array containing the approximated Jacobian                      |
| nFcn*    | fcn evaluation count adjusted                                   |
| iFail    | Return code non-zero if Jacobian could not be computed          |
"""
function nJacFDb{T}(fcn, n::Int64, lda::Int64, ml::Int64, x::Vector{T},
    fx::Vector{T}, yscal::Vector{T}, ajdel::T, ajmin::T,
    nFcn::Int64, a::Array{T,2})
    # Begin
    # Copy for internal purposes
    xa = x[:]
    iFail = 0
    mu = lda - 2*ml -1
    ldab = ml + mu + 1
    fu = zero(x)
    w  = zero(x)
    u  = zero(x)
    for jj = 1:ldab
        for k = jj:ldab:n
            w[k] = xa[k]
            su = sign(xa[k])
            if su == 0
                su = 1
            end
            u[k]  = max(max(abs(xa[k]),ajmin),yscal[k])*ajdel*su
            xa[k] = w[k] + u[k]
        end
        try
            fcn(fu,xa)
        catch
            iFail = -1
        end
        nFcn += 1
        if iFail != 0
            break;
        end
        for k = jj:ldab:n
            xa[k] = w[k]
            i1 = max(1,k-mu)
            i2 = min(n,k+ml)
            mh = mu + 1 - k
            a[mh+i1:mh+i2,k:k] = (fu[i1:i2]-fx[i1:i2])/u[k]
        end
    end
    return (nFcn,iFail)
end

"""
function nJcf{T}(fcn, n::Int64, lda::Int64, x::Vector{T}, fx::Vector{T},
    yscal::Vector{T}, eta::Vector{T}, etamin::T, etamax::T,
    etadif::T, conv::T, nFcn::Int64, a::Array{T,2})

Approximation of dense Jacobian matrix for nonlinear systems solver with
feed-back control of discretization and rounding errors

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable | Description                                                     |
|:---------|:----------------------------------------------------------------|
| fcn      | Function of the form fcn(f, x) to provide right-hand side       |
| n        | Number of rows and columns of the Jacobian                      |
| lda      | Leading dimension of the array "a"                              |
| x[n]     | Current scaled vector                                           |
| fx[n]    | Vector containing fcn(x)                                        |
| yscal[n] | Vector containing scaling factors                               |
| eta[n]*  | Vector of scaled denominator differences                        |
| etamin   | Minimum allowed scaled denominator                              |
| etamax   | Maximum allowed scaled denominator                              |
| etadif   | = sqrt(1.1*epMach)                                              |
| conv     | Maximum norm of last (unrelaxed) Newton correction              |
| nFcn*    | fcn evaluation count                                            |

(* marks inout parameters)

## Output parameters
-------------------
| Variable | Description                                                     |
|:---------|:----------------------------------------------------------------|
| a[lda,n] | Array containing the approximated Jacobian                      |
| eta[n]*  | Vector of scaled denominator differences adjusted               |
| nFcn*    | fcn evaluation count adjusted                                   |
| iFail    | Return code non-zero if Jacobian could not be computed          |
"""
function nJcf{T}(fcn, n::Int64, lda::Int64, x::Vector{T}, fx::Vector{T},
    yscal::Vector{T}, eta::Vector{T}, etamin::T, etamax::T,
    etadif::T, conv::T, nFcn::Int64, a::Array{T,2})
    # Constant
    small2 = 0.1
    # Copy for internal purposes
    xa = x[:]
    # Begin
    fu = zero(x)
    iFail = 0
    for k = 1:n
        is = 0
        qFine = false
        qExit = false
        while !qFine
            w = xa[k]
            su = sign(xa[k])
            if su == 0
                su = 1
            end
            u = eta[k]*yscal[k]*su
            xa[k] = w + u
            try
                fcn(fu,xa)
            catch
                iFail = -1
            end
            nFcn += 1
            if iFail != 0
                qExit = true
                break;
            end
            xa[k] = w
            sumd = 0.0
            for i = 1:n
                hg = max(abs(fx[i]),abs(fu[i]))
                fhi = fu[i] - fx[i]
                if hg != 0.0
                    sumd = sumd + (fhi/hg)^2
                end
                a[i,k] = fhi/u
            end
            sumd = sqrt(sumd/n)
            qFine = true
            if sumd != 0.0 && is == 0
                eta[k] = min(etamax,max(etamin,sqrt(etadif/sumd)*eta[k]))
                is = 1
                qFine = conv < small2 || sumd >= etamin
            end
        end
        if qExit
            break;
        end
    end
    return (nFcn,iFail)
end

"""
function nJcfb{T}(fcn, n::Int64, lda::Int64, ml::Int64, x::Vector{T},
    fx::Vector{T}, yscal::Vector{T}, eta::Vector{T},
    etamin::T, etamax::T, etadif::T, conv::T,
    nFcn::Int64, a::Array{T,2})

Approximation of banded Jacobian matrix for nonlinear systems solver with
feed-back control of discretization and rounding errors

T = Float64 or BigFloat

## Input parameters
-------------------
| Variable | Description                                                     |
|:---------|:----------------------------------------------------------------|
| fcn      | Function of the form fcn(f, x) to provide right-hand side       |
| n        | Number of rows and columns of the Jacobian                      |
| lda      | Leading dimension of the array "a"                              |
| ml       | Lower bandwidth of the Jacobian matrix                          |
| x[n]     | Current scaled vector                                           |
| fx[n]    | Vector containing fcn(x)                                        |
| yscal[n] | Vector containing scaling factors                               |
| eta[n]*  | Vector of scaled denominator differences                        |
| etamin   | Minimum allowed scaled denominator                              |
| etamax   | Maximum allowed scaled denominator                              |
| etadif   | = sqrt(1.1*epMach)                                              |
| conv     | Maximum norm of last (unrelaxed) Newton correction              |
| nFcn*    | fcn evaluation count                                            |

(* marks inout parameters)

## Output parameters
-------------------
| Variable | Description                                                     |
|:---------|:----------------------------------------------------------------|
| a[lda,n] | Array containing the approximated Jacobian                      |
| eta[n]*  | Vector of scaled denominator differences adjusted               |
| nFcn*    | fcn evaluation count adjusted                                   |
| iFail    | Return code non-zero if Jacobian could not be computed          |
"""
function nJcfb{T}(fcn, n::Int64, lda::Int64, ml::Int64, x::Vector{T},
    fx::Vector{T}, yscal::Vector{T}, eta::Vector{T},
    etamin::T, etamax::T, etadif::T, conv::T,
    nFcn::Int64, a::Array{T,2})
    # Constants
    small2 = 0.1
    # Copy for internal purposes
    xa = x[:]
    # Begin
    mu = lda - 2*ml - 1
    ldab = ml + mu + 1
    fu = zero(x)
    w  = zero(x)
    u  = zero(x)
    iFail = 0
    for jj = 1:ldab
        is = 0
        qFine = false
        qExit = false
        while !qFine
            for k = jj:ldab:n
                w[k] = xa[k]
                su = sign(xa[k])
                if su == 0
                    su = 1
                end
                u[k] = eta[k]*yscal[k]*su
                xa[k] = w[k] + u[k]
            end

            try
                fcn(fu,xa)
            catch
                iFail = -1
            end
            nFcn += 1
            if iFail != 0
                qExit = true
                break;
            end

            for k = jj:ldab:n
                xa[k] = w[k]
                sumd = 0.0
                i1 = max(1,k-mu)
                i2 = min(n,k+ml)
                mh = mu + 1 - k
                for i = i1:i2
                    hg = max(abs(fx[i]),abs(fu[i]))
                    fhi = fu[i] - fx[i]
                    if hg != 0.0
                        sumd += (fhi/hg)^2
                    end
                    a[mh+i,k:k] = fhi/u[k]
                end
                sumd = sqrt(sumd/n)
                qFine = true
                if sumd != 0.0 && is == 0
                    eta[k] = min(etamax,max(etamin,sqrt(etadif/sumd)*eta[k]))
                    is = 1
                    qFine = conv < small2 || sumd >= etamin
                end
            end
        end
        if qExit
            break;
        end
    end
    return (nFcn,iFail)
end


text:
module PkgB

using Distributions

randchisq2() = sum(randn(2) .^ 2)

end # module


text:
module Escher

using JSON, Makie

"""
    _bezier(t, p0, p1, p2, p3)

Plot a Bezier curve with `t` in `[0,1]` and fixed points `p0`, `p1`, `p2`, `p3`.
"""
_bezier(t, p0, p1, p2, p3) =
    (1 - t)^3 .* p0 + 3 * (1 - t)^2 * t .* p1 + 3 * (1 - t) * t^2 .* p2 + t^3 .* p3

"""
    _nodes(
        escher;
        primary_node_size = 5,
        secondary_node_size = 3,
        metabolite_identifier = "bigg_id",
        metabolite_node_sizes = Dict{String,Any}(),
        metabolite_node_colors = Dict{String,Any}(),
        metabolite_node_color = :black,
    )

Helper function to extract metabolite node information from `escher`. The kwargs set plot
details.
"""
function _nodes(
    escher;
    primary_node_size = 5,
    secondary_node_size = 3,
    metabolite_identifier = "bigg_id",
    metabolite_node_sizes = Dict{String,Any}(),
    metabolite_node_colors = Dict{String,Any}(),
    metabolite_node_color = :black,
)
    nodexs = Float64[]
    nodeys = Float64[]
    labelxs = Float64[]
    labelys = Float64[]
    metabolite_labels = String[]
    node_pos = Dict()
    markersizes = Float64[]
    markercolors = []
    metabolite_nodes = String[]
    for (node_id, node) in escher["nodes"]
        if haskey(node, "x") && haskey(node, "y")
            x = node["x"]
            y = -node["y"]
            node_pos[node_id] = (x, y)
            # plotted nodes
            if node["node_type"] == "metabolite"
                push!(metabolite_nodes, node_id)
                if haskey(node, metabolite_identifier)
                    push!(metabolite_labels, node[metabolite_identifier])
                    push!(labelxs, node["label_x"])
                    push!(labelys, -node["label_y"])
                end
                push!(nodexs, x)
                push!(nodeys, y)
                if haskey(metabolite_node_colors, node[metabolite_identifier])
                    push!(markercolors, metabolite_node_colors[node[metabolite_identifier]])
                else
                    push!(markercolors, metabolite_node_color) # fallback
                end
                if haskey(metabolite_node_sizes, node[metabolite_identifier])
                    push!(markersizes, metabolite_node_sizes[node[metabolite_identifier]])
                else
                    if node["node_is_primary"] # fallback
                        push!(markersizes, primary_node_size)
                    else
                        push!(markersizes, secondary_node_size)
                    end
                end
            end
        end
    end
    return nodexs,
    nodeys,
    node_pos,
    markersizes,
    labelxs,
    labelys,
    metabolite_labels,
    markercolors,
    metabolite_nodes
end

"""
Plot a metabolic map that is compatible with Escher. The only required argument is the
location of the metabolic map in `json` format.

# Example
```
escherplot("core-map.json"; kwargs...)
```
Here `kwargs` are supported attributes, see the `readme` for more information.

# Attributes
```
metabolite_identifier = "bigg_id",
metabolite_show_text = false,
metabolite_text_size = 4,
metabolite_primary_node_size = 5,
metabolite_secondary_node_size = 3,
metabolite_node_sizes = Dict{String,Any}(),
metabolite_node_colors = Dict{String,Any}(),
metabolite_node_color = :black, # fallback color
metabolite_text_color = :black,
reaction_identifier = "bigg_id",
reaction_show_text = false,
reaction_show_name_instead_of_id = false,
reaction_text_size = 4,
reaction_text_color = :black,
reaction_edge_colors = Dict{String,Any}(),
reaction_edge_color = :black, # fallback color
reaction_edge_widths = Dict{String,Any}(),
reaction_edge_width = 2.0, # fallback width
```
Get or create maps here: `https://escher.github.io/#/`.
"""
@recipe(EscherPlot) do scene
    Attributes(
        metabolite_identifier = "bigg_id",
        metabolite_show_text = false,
        metabolite_text_size = 4,
        metabolite_primary_node_size = 5, # fallback size
        metabolite_secondary_node_size = 3, # fallback size
        metabolite_node_sizes = Dict{String,Any}(),
        metabolite_node_colors = Dict{String,Any}(),
        metabolite_node_color = :black, # fallback color
        metabolite_text_color = :black,
        reaction_identifier = "bigg_id",
        reaction_show_text = false,
        reaction_show_name_instead_of_id = false,
        reaction_text_size = 4,
        reaction_text_color = :black,
        reaction_edge_colors = Dict{String,Any}(), # actual color
        reaction_edge_color = :black, # fallback color
        reaction_edge_widths = Dict{String,Any}(), # actual edge width
        reaction_edge_width = 2.0, # fallback width
        reaction_arrow_size = 6,
        reaction_directions = Dict{String,Any}(), # rid => :f or :b,  relative to model used to construct map
    )
end

function Makie.plot!(ep::EscherPlot{<:Tuple{String}})

    _, escher = JSON.parsefile(to_value(ep[1]))

    metabolite_node_xs,
    metabolite_node_ys,
    metabolite_node_pos,
    metabolite_markersizes,
    metabolite_label_xs,
    metabolite_label_ys,
    metabolite_labels,
    metabolite_marker_colors,
    metabolite_nodes = _nodes(
        escher;
        primary_node_size = to_value(ep.metabolite_primary_node_size),
        secondary_node_size = to_value(ep.metabolite_secondary_node_size),
        metabolite_identifier = to_value(ep.metabolite_identifier),
        metabolite_node_sizes = to_value(ep.metabolite_node_sizes),
        metabolite_node_colors = to_value(ep.metabolite_node_colors),
        metabolite_node_color = to_value(ep.metabolite_node_color),
    )

    reaction_labels = String[]
    reaction_label_xs = Float64[]
    reaction_label_ys = Float64[]
    reaction_id_to_id = Dict(
        rid => rxn[to_value(ep.reaction_identifier)] for (rid, rxn) in escher["reactions"]
    )
    for (rid, rxn) in escher["reactions"]
        # collect reaction label information
        if ep.reaction_show_name_instead_of_id[]
            if haskey(rxn, "name")
                push!(reaction_labels, rxn["name"])
                push!(reaction_label_xs, rxn["label_x"])
                push!(reaction_label_ys, -rxn["label_y"])
            end
        else
            if haskey(rxn, ep.reaction_identifier[])
                push!(reaction_labels, rxn[ep.reaction_identifier[]])
                push!(reaction_label_xs, rxn["label_x"])
                push!(reaction_label_ys, -rxn["label_y"])
            end
        end

        val_rid = reaction_id_to_id[rid]

        lw =
            haskey(to_value(ep.reaction_edge_widths), val_rid) ?
            to_value(ep.reaction_edge_widths)[val_rid] : ep.reaction_edge_width
        c =
            haskey(to_value(ep.reaction_edge_colors), val_rid) ?
            to_value(ep.reaction_edge_colors)[val_rid] : ep.reaction_edge_color

        if isempty(to_value(ep.reaction_edge_widths)) &&
           isempty(to_value(ep.reaction_edge_colors)) # plot everything because no data supplied
            no_reaction_data = false
        elseif haskey(to_value(ep.reaction_edge_widths), val_rid) ||
               haskey(to_value(ep.reaction_edge_colors), val_rid) # data supplied, test if rxn has data
            no_reaction_data = false
        else
            no_reaction_data = true # data supplied, but no rxn data
        end

        # linestyle = :dot if missing information
        for (_, segment) in rxn["segments"]
            if !isnothing(segment["b1"]) && !isnothing(segment["b2"])
                n1 = metabolite_node_pos[segment["from_node_id"]]
                p0 = [n1[1], n1[2]]
                p1 = [segment["b1"]["x"], -segment["b1"]["y"]]
                p2 = [segment["b2"]["x"], -segment["b2"]["y"]]
                n2 = metabolite_node_pos[segment["to_node_id"]]
                p3 = [n2[1], n2[2]]
                ts = range(0.0, 1.0; length = 100) # higher resolution => looks smoother
                ds = _bezier.(ts, Ref(p0), Ref(p1), Ref(p2), Ref(p3))
                xs = [d[1] for d in ds]
                ys = [d[2] for d in ds]
            else
                n1 = metabolite_node_pos[segment["from_node_id"]]
                n2 = metabolite_node_pos[segment["to_node_id"]]
                xs = range(n1[1], n2[1]; length = 100) # higher resolution => looks smoother
                ys = range(n1[2], n2[2]; length = 100) # higher resolution => looks smoother
            end
            # Draw the curves
            if isempty(to_value(ep.reaction_directions)) || !haskey(to_value(ep.reaction_directions), rid)
                # normal
                if segment["to_node_id"] in metabolite_nodes
                    arrow_head_offset = 15
                    arrows!(
                        ep,
                        xs[end-arrow_head_offset:end-arrow_head_offset],
                        ys[end-arrow_head_offset:end-arrow_head_offset],
                        [xs[end-arrow_head_offset+1] - xs[end-arrow_head_offset],],
                        [ys[end-arrow_head_offset+1] - ys[end-arrow_head_offset],],
                        color = c,
                        arrowcolow = c,
                        linewidth = lw,
                        arrowsize = to_value(ep.reaction_arrow_size),
                    )
                end
                if get(rxn, "reversibility", false) && segment["from_node_id"] in metabolite_nodes
                    arrow_head_offset = 85
                    arrows!(
                        ep,
                        xs[end-arrow_head_offset:end-arrow_head_offset],
                        ys[end-arrow_head_offset:end-arrow_head_offset],
                        -[xs[end-arrow_head_offset+1] - xs[end-arrow_head_offset],],
                        -[ys[end-arrow_head_offset+1] - ys[end-arrow_head_offset],],
                        color = c,
                        arrowcolow = c,
                        linewidth = lw,
                        arrowsize = to_value(ep.reaction_arrow_size),
                    )
                end
            elseif to_value(ep.reaction_directions)[rid] == :f # forward
                if segment["to_node_id"] in metabolite_nodes
                    arrow_head_offset = 15
                    arrows!(
                        ep,
                        xs[end-arrow_head_offset:end-arrow_head_offset],
                        ys[end-arrow_head_offset:end-arrow_head_offset],
                        [xs[end-arrow_head_offset+1] - xs[end-arrow_head_offset],],
                        [ys[end-arrow_head_offset+1] - ys[end-arrow_head_offset],],
                        color = c,
                        arrowcolow = c,
                        linewidth = lw,
                        arrowsize = to_value(ep.reaction_arrow_size),
                    )
                end
            elseif to_value(ep.reaction_directions)[rid] == :b # reverse
                if segment["from_node_id"] in metabolite_nodes
                    arrow_head_offset = 85
                    arrows!(
                        ep,
                        xs[end-arrow_head_offset:end-arrow_head_offset],
                        ys[end-arrow_head_offset:end-arrow_head_offset],
                        -[xs[end-arrow_head_offset+1] - xs[end-arrow_head_offset],],
                        -[ys[end-arrow_head_offset+1] - ys[end-arrow_head_offset],],
                        color = c,
                        arrowcolow = c,
                        linewidth = lw,
                        arrowsize = to_value(ep.reaction_arrow_size),
                    )
                end
            end

            if no_reaction_data
                lines!(
                    ep,
                    xs,
                    ys,
                    linewidth = lw,
                    linestyle = :dot,
                    color = c,
                )
            else
                lines!(
                    ep,
                    xs,
                    ys,
                    linewidth = lw,
                    color = c,
                )
            end
        end
    end

    scatter!(
        ep,
        metabolite_node_xs,
        metabolite_node_ys,
        color = metabolite_marker_colors,
        markersize = metabolite_markersizes,
    )

    if ep.metabolite_show_text[]
        positions = [
            (metabolite_label_xs[i], metabolite_label_ys[i]) for
            i in eachindex(metabolite_labels)
        ]
        text!(
            ep,
            metabolite_labels;
            position = positions,
            textsize = ep.metabolite_text_size,
            color = ep.metabolite_text_color,
        )
    end

    if ep.reaction_show_text[]
        positions = [
            (reaction_label_xs[i], reaction_label_ys[i]) for i in eachindex(reaction_labels)
        ]
        text!(
            ep,
            reaction_labels;
            position = positions,
            textsize = ep.reaction_text_size,
            color = ep.reaction_text_color,
        )

    end

    ep
end

export plot!

end # module


text:
## Preloads
using Statistics
using FFTW
using Plots
using BenchmarkTools
using Profile
using LinearAlgebra
using Measures
using ProgressMeter
push!(LOAD_PATH, pwd()*"/main")
using DHC_2DUtils

fhash = fink_filter_hash(1, 8, wd=2, nx=256, Omega=false)

randA = rand(128,128)

no_shift = DHC_compute(randA,fhash)

shift = DHC_compute(circshift(randA,(1,1)),fhash)

maximum(shift - no_shift)

coeff_vec = zeros(256*256,2452)
@showprogress for i=1:256
    @showprogress for j=1:256
        coeff_vec[i+(j-1)*256,:] = DHC_compute(circshift(randA,(i,j)),fhash)
    end
end

using Distributed
addprocs(7);

@everywhere begin
    using Statistics
    using FFTW
    using Plots
    using BenchmarkTools
    using Profile
    using LinearAlgebra
    using Measures
    using ProgressMeter
    push!(LOAD_PATH, pwd()*"/main")
    using DHC_2DUtils

    fhash = fink_filter_hash(1, 8, wd=2, nx=128, Omega=false)

    function calc_shift(x)
        i,j,img = x
        return DHC_compute(circshift(img,(i,j)),fhash)
    end
end

mnist_DHC_out = @showprogress pmap(calc_shift, Iterators.product(1:128,1:128,[randA]))

mnist_DHC_out = hcat(mnist_DHC_out...)

out_test = mnist_DHC_out

maximum(out_test .- out_test[:,1])

max_stars_repo_name: andrew-saydjari/DHC

max_issues_repo_name: andrew-saydjari/DHC

max_forks_repo_name: andrew-saydjari/DHC


text:
using Winston

x = linspace(0., 2pi, 40)
s = sin(x)
c = cos(x)

inset = FramedPlot(title="inset")
setattr(inset.frame, draw_ticks=false)

add(inset, Curve(x, s, kind="dashed"))

p = FramedPlot(aspect_ratio=1)
setattr(p.frame, tickdir=+1, draw_spine=false)

add(p, SymmetricErrorBarsY(x, s, 0.2*ones(length(x))))
add(p, Points(x, s, color="red"))
add(p, PlotInset((.6,.6), (.95,.95), inset))

file(p, "example4.png")


text:
### A Pluto.jl notebook ###
# v0.12.20

using Markdown
using InteractiveUtils

# ╔═╡ b213914c-6fa0-11eb-0ba6-950195900c41
using DifferentialEquations, CairoMakie, StaticArrays

# ╔═╡ 853114f6-6fa0-11eb-39a8-edca713a7c1a
md"### Assignment 5

**Course on Physics of Compact Objects**

Instructor: Prof. P. Ajith

Author: Md Arif Shaikh, ICTS-TIFR, Feb 15, 2021"

# ╔═╡ a5da3e0a-6fa0-11eb-2518-ef262ae935d9
md"### Problem 2
Numerically solve the stellar structure equations

$\frac{dm}{dr} = 4\pi r^2 \rho, \frac{dP}{dr} = - \frac {G m \rho}{r^2}$

assuming the equation of state of non-relativistic degenerate electron gas *including Coulomb correction*s:

$ P(\rho) = K \rho^{5/3}\left[ 1 − (\rho_0 /\rho)^{1/3}\right],$

where $K \simeq 10^{13} \mu_e^{-5/3}$ (cgs units) and $\rho_0 \simeq 0.4 Z^2 \mu_e$ g/cm$^3$. Assume the atomic numeber to be $Z = 6$ (Carbon) an the mean molecular weight per free electron to be $\mu_e = 2$. The obvious boundary conditions
are $m(r = 0) = 0$ and $P(r = 0) = P(\rho_c)$, where $\rho_c$ is the central density. Compute the mass-radius of white dwarfs for central densities ranging $\rho_c = 10^6 − 10^9$ g/cm$^3$ . Compare the mass-radius relation with the same computed earlier
neglecting Coulomb corrections.
"

# ╔═╡ fe34ad36-6fa0-11eb-1187-11d611ecdbb7
begin
	const G = 6.6743 * 1e-8
	const μe = 2.0
	const K = 1e13 * μe^(-5.0 / 3.0)
	const Z = 2.0
	const ρ0 = 0.4 * Z^2 * μe
end

# ╔═╡ 0f31b1e2-6fa1-11eb-3a90-19e6e49800c4
function stellarStructure(u, params, r)
    m, ρ = u
    ρ0, G, K = params
    dmdr = 4.0 * pi * r^2 * abs(ρ)
    dρdr = - (3.0 * G * m * abs(ρ)^(1.0 /3.0 )) / (K * r^2 * (5.0 - 4.0 * abs((ρ0 / ρ))^(1.0 / 3.0)))
    @SVector [dmdr, dρdr]
end

# ╔═╡ 1c2ef3ca-6fa1-11eb-3946-453e872da2ca
begin
	condition(u, r, integrator) =  u[2] - ρ0;
	affect!(integrator) = terminate!(integrator);
	cb = ContinuousCallback(condition, affect!);
end

# ╔═╡ 371db5f2-6fa1-11eb-3546-cf1fa66b48cf
begin
	rspan = [1e-7, 1e10]
	params = [ρ0, G, K]
	ρc = 1e8

	u0 = @SVector [0.0, ρc]
	problem = ODEProblem(stellarStructure, u0, rspan, params)
	solution = solve(problem, callback = cb)
end

# ╔═╡ b115c610-6fa6-11eb-3bf9-eb593e3f7d78
begin
	fig = Figure(resolution = (800, 450), font = "Times New Roman", fontsize = 14)
	ax = fig[1, 1] = Axis(fig, xlabel = "radius [cgs]", ylabel = "density [cgs]")
	lines!(solution.t, solution[2, :], color = :red)
	save("./test.pdf", fig)
	fig
end

# ╔═╡ 83c00ae4-6fa7-11eb-121a-a76e67b869e3
begin
	masses = Float64[]
	radii = Float64[]
	ρcs = LinRange(1.e6, 1.e8, 100)
	for ρc in ρcs
		u0 = @SVector [0.0, ρc]
		problem = ODEProblem(stellarStructure, u0, rspan, params)
		sol = solve(problem, callback = cb)
		push!(radii, sol.t[end])
		push!(masses, sol[1, :][end])
	end
end

# ╔═╡ 7bdf068c-6fab-11eb-31fd-b5b610054ccb
Msun = 1.989 * 1e33

# ╔═╡ eb191c4a-6fab-11eb-18a9-fdad808bddc5
Rsun = 696340 * 1e5

# ╔═╡ 0eb3412c-6fab-11eb-12d0-29d1a35e11ba
begin
	figMR = Figure(resolution = (800, 450), font = "Times New Roman", fontsize = 14)
	axMR = figMR[1, 1] = Axis(figMR,
		ylabel = "radius [km]",
		xlabel = "mass [Msun]")
	lines!(masses / Msun, radii / 1e5, color = :red)
	xlims!(axMR, [masses[1], masses[end]]/Msun)
	ylims!(axMR, [radii[end], radii[1]]/1e5)
	figMR
end

# ╔═╡ Cell order:
# ╟─853114f6-6fa0-11eb-39a8-edca713a7c1a
# ╟─a5da3e0a-6fa0-11eb-2518-ef262ae935d9
# ╠═b213914c-6fa0-11eb-0ba6-950195900c41
# ╠═fe34ad36-6fa0-11eb-1187-11d611ecdbb7
# ╠═0f31b1e2-6fa1-11eb-3a90-19e6e49800c4
# ╠═1c2ef3ca-6fa1-11eb-3946-453e872da2ca
# ╠═371db5f2-6fa1-11eb-3546-cf1fa66b48cf
# ╠═b115c610-6fa6-11eb-3bf9-eb593e3f7d78
# ╠═83c00ae4-6fa7-11eb-121a-a76e67b869e3
# ╠═7bdf068c-6fab-11eb-31fd-b5b610054ccb
# ╠═eb191c4a-6fab-11eb-18a9-fdad808bddc5
# ╠═0eb3412c-6fab-11eb-12d0-29d1a35e11ba


text:
import Pkg # stdlib

struct TomlFile
    filename::String
    function TomlFile(path::String)::TomlFile
        path::String = abspath(strip(path))
        if isfile(path)
            result::TomlFile = new(path)
            return result
        else
            error("File does not exist")
        end
    end
end

function parse_toml_file(x::TomlFile)::Dict{String, Any}
    toml_file_filename::String = x.filename
    toml_file_text::String = read(toml_file_filename, String)
    toml_file_parsed::Dict{String, Any} = Pkg.TOML.parse(toml_file_text)
    return toml_file_parsed
end

function version_string(x::TomlFile)::String
    toml_file_parsed::Dict{String, Any} = parse_toml_file(x)
    version_string::String = toml_file_parsed["version"]
    return version_string
end

function version_string()::String
    predictmd_toml_file::TomlFile = TomlFile(
        package_directory("Project.toml")
        )
    resultversion_string::String = version_string(predictmd_toml_file)
    return resultversion_string
end

function api_version_string()::String
    predictmdapi_toml_file::TomlFile = TomlFile(
        api_package_directory("Project.toml")
        )
    resultversion_string::String = version_string(predictmdapi_toml_file)
    return resultversion_string
end

# function version_string(m::Method)::String
#     m_package_directory::String = package_directory(m)
#     m_toml_file::TomlFile = TomlFile(
#         joinpath(m_package_directory, "Project.toml")
#         )
#     resultversion_string::String = version_string(m_toml_file)
#     return resultversion_string
# end

# function version_string(f::Function)::String
#     m_package_directory::String = package_directory(f)
#     m_toml_file::TomlFile = TomlFile(
#         joinpath(m_package_directory, "Project.toml")
#         )
#     resultversion_string::String = version_string(m_toml_file)
#     return resultversion_string
# end

# function version_string(f::Function, types::Tuple)::String
#     m_package_directory::String = package_directory(f, types)
#     m_toml_file::TomlFile = TomlFile(
#         joinpath(m_package_directory, "Project.toml")
#         )
#     resultversion_string::String = version_string(m_toml_file)
#     return resultversion_string
# end

# function version_string(m::Module)::String
#     m_package_directory::String = package_directory(m)
#     m_toml_file::TomlFile = TomlFile(
#         joinpath(m_package_directory, "Project.toml")
#         )
#     resultversion_string::String = version_string(m_toml_file)
#     return resultversion_string
# end

"""
    version()::VersionNumber

Return the version number of PredictMD.
"""
function version()::VersionNumber
    resultversion_string::String = version_string()
    result_versionnumber::VersionNumber = VersionNumber(resultversion_string)
    return result_versionnumber
end

function api_version()::VersionNumber
    resultversion_string::String = api_version_string()
    result_versionnumber::VersionNumber = VersionNumber(resultversion_string)
    return result_versionnumber
end

"""
    version(m::Method)::VersionNumber

If method `m`
is part of a Julia package, returns the version number of that package.

If method `m`
is not part of a Julia package, throws an error.
"""
# function version(m::Method)::VersionNumber
#     resultversion_string::String = version_string(m)
#     result_versionnumber::VersionNumber = VersionNumber(resultversion_string)
#     return result_versionnumber
# end

"""
    version(f::Function)::VersionNumber

If function `f`
is part of a Julia package, returns the version number of
that package.

If function `f`
is not part of a Julia package, throws an error.
"""
# function version(f::Function)::VersionNumber
#     resultversion_string::String = version_string(f)
#     result_versionnumber::VersionNumber = VersionNumber(resultversion_string)
#     return result_versionnumber
# end

"""
    version(f::Function, types::Tuple)::VersionNumber

If function `f` with type signature `types`
is part of a Julia package, returns the version number of
that package.

If function `f` with type signature `types`
is not part of a Julia package, throws an error.
"""
# function version(f::Function, types::Tuple)::VersionNumber
#     resultversion_string::String = version_string(f, types)
#     result_versionnumber::VersionNumber = VersionNumber(resultversion_string)
#     return result_versionnumber
# end

"""
    version(m::Module)::VersionNumber

If module `m` is part of a Julia package, returns the version number of
that package.

If module `m` is not part of a Julia package, throws an error.
"""
# function version(m::Module)::VersionNumber
#     resultversion_string::String = version_string(m)
#     result_versionnumber::VersionNumber = VersionNumber(resultversion_string)
#     return result_versionnumber
# end


text:
module DemoPackageCT

# Write your package code here.
print_greeting(io::IO = stdout) = print(io, "Hello, world!")

end


text:
using Distributed
addprocs(3)
@everywhere include("reproducibility.jl")
using Random
Random.seed!(42)
##
@everywhere using QSWalk
@everywhere using LightGraphs
@everywhere using LinearAlgebra,SparseArrays
@everywhere using NPZ
data_folder_name = "nonmoralizing_qsw"

## ######################################################################
# #################### premature localization test ######################
# #######################################################################
function gue_random(n::Int)
    @assert n >= 1
    result = zeros(ComplexF64, n, n)
    for i=1:n, j=i:n
        result[i,j] = randn()
    end

    for i=1:n, j=(i+1):n
        result[i,j] += 1im*randn()
        result[j,i] = conj(result[i,j])
    end
    result
end

function goe_random(n::Int)
    @assert n >= 1
    result = zeros(Float64, n, n)
    for i=1:n, j=i:n
        result[i,j] = randn()
        result[j,i] = result[i,j]
    end
    result
end

repeating = 20_000

g = DiGraph(4)
add_edge!(g, 1, 2)
add_edge!(g, 2, 1)
add_edge!(g, 2, 3)
add_edge!(g, 3, 2)
add_edge!(g, 1, 3)
add_edge!(g, 3, 1)
add_edge!(g, 1, 4)

adj = adjacency_matrix(g, Float64, dir=:in)
ll, vset = QSWalk.nm_lind(adj)

data = zeros(repeating)
for m=1:repeating
    block_ham = Dict(v => gue_random(length(v)) for v=vlist(vset))
    rot_ham = nm_loc_ham(vset, block_ham)

    s = evolve_generator(zeros(size(ll)...), [ll], rot_ham)
    data[m] = sort(abs.(eigvals(s)))[2]
end
npzwrite("$data_folder_name/gue_loc_ham.npz", data)

data = zeros(repeating)
for m=1:repeating
    block_ham = Dict(v => goe_random(length(v)) for v=vlist(vset))
    rot_ham = nm_loc_ham(vset, block_ham)

    s = evolve_generator(zeros(size(ll)...), [ll], rot_ham)
    data[m] = sort(abs.(eigvals(s)))[2]
end
npzwrite("$data_folder_name/goe_loc_ham.npz", data)

## ######################################################################
# ######################### symmetry ####################################
# #######################################################################

n = 61
@assert isodd(n)
t = 100

g = PathGraph(n)
a = adjacency_matrix(g)
mid = div(n+1, 2)

l, vset  = nm_lind(a)
h_rot = nm_loc_ham(vset)

s = evolve_generator(spzeros(size(l)...), [l], h_rot)
ϱ = nm_init(vset[[mid]], vset)
ϱ_t = evolve(s, ϱ, t)
npzwrite("$data_folder_name/probdist_non_symmetric.npz", nm_measurement(ϱ_t, vset))

#

linddescription1 = Dict{Int,Matrix{Float64}}(1 => ones(1, 1), 2 => [1 1; 1 -1])
linddescription2 = Dict{Int,Matrix{Float64}}(1 => ones(1, 1), 2 => [1 1; -1 1])
lind1, vset = nm_lind(a, linddescription1)
lind2, _ = nm_lind(a, linddescription2);
h_rot = nm_loc_ham(vset)

s = evolve_generator(spzeros(size(l)...), [lind1, lind2], h_rot)
ϱ = nm_init(vset[[mid]], vset)
ϱ_t = evolve(s, ϱ, t)
npzwrite("$data_folder_name/probdist_symmetric.npz", nm_measurement(ϱ_t, vset))

## ######################################################################
# ######################### propagation ##################################
# #######################################################################

ωs = 0.5:0.1:1.
max_n = 400
t_scales = 0.1:0.1:10
scaling = 1/3
#
αs_global = []
for t_scale = t_scales
    t = t_scale*max_n
    n = max(max_n, ceil(Int, t*scaling))
    n += iseven(n)
    mid = ceil(Int, div(n+1, 2))

    g = PathGraph(n)
    a = adjacency_matrix(g)
    linddescription1 = Dict{Int,Matrix{Float64}}(1 => ones(1, 1), 2 => [1 1; 1 -1])
    linddescription2 = Dict{Int,Matrix{Float64}}(1 => ones(1, 1), 2 => [1 1; -1 1])
    lind1, vset = nm_lind(a, linddescription1)
    lind2, _ = nm_lind(a, linddescription2)
    h_rot = nm_loc_ham(vset)
    h = nm_glob_ham(a)

    s = pmap(ω->evolve_generator(h, [lind1, lind2], h_rot, ω), ωs)
    println(typeof(s))
    initial_state = nm_init(vset[[mid]], vset)
    @time prob_dists = pmap(x-> nm_measurement(evolve(x, initial_state, t), vset), s)

    positions = -(mid-1):(mid-1)
    sec_moments = pmap(pdist -> sum(pdist .* positions.^2), prob_dists)
    push!(αs_global, log.(sec_moments))
    npzwrite("$data_folder_name/scaling_data.npz", hcat(αs_global...))
end