text:
"""
    fit_ir(x::AbstractVector, fs, nfft, max_decay)
    fit_ir(x::SampleBuf, nfft, max_decay)

Analyze the given impulse response `x` at sample rate `fs` by curve fitting in
each frequency band. The impulse response audio should be relatively clean (the
impulse should be the time-domain peak, and there should not be too much
nonstationary background noise).

`nfft` gives the number of frequency bands that will be used.

`max_decay` should be the longest time (in seconds) that it might take for
the impulse to decay to the noise floor. Setting this too high may affect
convergence, and setting it too low may effect accuracy.

Returns a vector of named tuples for each band containing:
    - `freq` - the band center frequency (in Hz)
    - `noise_pow` - noise floor (in dB)
    - `peak_snr` - Peak SNR (in dB)
    - `rt60` - RT60 for each band (in seconds)
    - `drr` - Direct-to-Reverberant Ratio for each band (in dB)

If the curve fit does not converge the values may be `missing`.
"""
fit_ir(x::SampleBuf, nfft, max_decay) =
    fit_ir(x.data, samplerate(x), nfft, max_decay)

struct IRFitException <: Exception
    msg::String
end

Base.showerror(io::IO, e::IRFitException) = print(io, e.msg)

function fit_ir(x::AbstractVector, fs, nfft, max_decay)
    maxidx = findmax(abs2.(x))[2]
    # now we take the STFT to generate our frequency bands
    hop = nfft÷2
    if maxidx <= 2hop
        throw(IRFitException("IR Fitting needs at least $(2hop) samples (2 hops) prior to peak (at $maxidx). Provide more signal or decrease your hop size"))
    end
    if maxidx >= length(x) - max_decay*fs
        throw(IRFitException("Signal peak is too close to the end of the signal"))
    end
    # set things up so that the 2nd STFT frame is centered on the peak.
    # this also ensures that no padding is present in that window
    spec = stft2(x[maxidx-2hop:end], nfft, hop; window=hanning)
    specpow = abs2.(spec)
    logspecpow = 10log10.(specpow)

    # TODO: we could probably re-structure this to also allow for nmax
    # starting too low, but then it seems less guaranteed that it will
    # converge
    nmax_init = round(Int, max_decay*fs/hop)

    noise = 10log10.(mean(specpow[:, nmax_init:end], dims=2))

    map(1:size(specpow, 1)) do k
        converged = false
        # The peak index should be 2, but in practice it seems delayed for some
        # frequencies. All the following analysis is performed relative to the
        # peak position
        npeak = findmax(logspecpow[k, :])[2]
        # if the peak is too far delayed from the time-domain peak, it's
        # probably garbage. 10 is pretty arbitrarily-chosen.
        if npeak > 10
            peakdelay = round(Int,(npeak-2)*hop/fs*1000)
            # @warn "Within-band peak $(peakdelay)ms from time-domain peak, skipping"
            @goto cleanup
        end
        # we'll use a fixed endpoint rather then making it relative to the
        # per-band peak so that if the peak gets found very late in the signal
        # it just won't get used.
        nmax = nmax_init
        # from Traer they say the decay starts after about 20ms. We want to make
        # sure we're in the exponential decay portion, so we wayt 50ms  after
        # the peak:
        # nstart = ceil(Int, 0.05 * fs / hop) + npeak
        # unfortunately this process somewhat sensitive to how much of the early
        # IR we include - if I try to skip the first part as in the
        # commented-out code above then it seems to give much more shallow fits,
        # which correspond to much larger RT60s. It's not super clear from Traer
        # whether he included the early IR in the fit.
        nstart = npeak

        # let's try a linear fit from nstart to nmax. We'll then use that to
        # figure out where the decay hits the noise floor. Initially the fitted
        # data will include the noise that's after the decay, so we'll
        # iteratively shorten the region we're fitting over until it converged.
        local L

        while true
            a = [nstart:nmax ones(nmax-nstart+1)]
            b = logspecpow[k, nstart:nmax]
            # aL = b, so a \ b = L
            L = a \ b
            # positive slope - not going to converge
            L[1] >= 0 && break
            # starting below noise floor - not going to converge
            L[2] < noise[k] && break
            nmax_prev = nmax
            nmax = round(Int, (noise[k] - L[2]) / L[1])
            nmax >= nmax_init && break
            # fit region is too small to be reliable
            nmax <= nstart + 10 && break
            noise[k] = 10log10.(mean(specpow[k, nmax:end]))
            # I think it's possible that this could get into an oscillation,
            # so we'll consider any increase in the noise intersect point to be
            # convergence
            if nmax >= nmax_prev
                converged = true
                break
            end
        end

        # we'll estimate the peak power with a quadratic interpolation. The
        # point we interpolate at is based on the relative magnitudes of the
        # neighboring points.
        # details at:
        # https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html

        peakpow, peakpos = if npeak > 1
            α, β, γ = specpow[k, (-1:1).+npeak]
            p = 0.5 * (α-γ)/(α-2β+γ)
            pow = β - 0.25 * (α-γ)*p
            (10log10(pow), p+npeak)
        else
            # this probably means the data isn't good for this band
            (10log10(specpow[k, npeak]), npeak)
        end

        @label cleanup
        freq=(k-1)/nfft*fs
        if converged
            (freq=(k-1)/nfft*fs,
             noise_pow=noise[k], # noise power
             peak_snr=peakpow - noise[k], # peak SNR
             drr=peakpow - (L[2] + peakpos*L[1]), # DRR
             rt60=-60/L[1]*hop/fs) # RT60
        else
            (freq=(k-1)/nfft*fs,
             noise_pow=missing, # noise power
             peak_snr=missing, # peak SNR
             drr=missing, # DRR
             rt60=missing) # RT60
        end
    end
end

max_stars_repo_name: ssfrr/SFRDSP

max_issues_repo_name: ssfrr/SFRDSP

max_forks_repo_name: ssfrr/SFRDSP


text:
doc"""
    Levy(μ, σ)

The *Lévy distribution* with location `μ` and scale `σ` has probability density function

$f(x; \mu, \sigma) = \sqrt{\frac{\sigma}{2 \pi (x - \mu)^3}}
\exp \left( - \frac{\sigma}{2 (x - \mu)} \right), \quad x > \mu$

```julia
Levy()         # Levy distribution with zero location and unit scale, i.e. Levy(0, 1)
Levy(u)        # Levy distribution with location u and unit scale, i.e. Levy(u, 1)
Levy(u, c)     # Levy distribution with location u ans scale c

params(d)      # Get the parameters, i.e. (u, c)
location(d)    # Get the location parameter, i.e. u
```

External links

* [Lévy distribution on Wikipedia](http://en.wikipedia.org/wiki/Lévy_distribution)
"""
immutable Levy{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T

    Levy(μ::T, σ::T) = (@check_args(Levy, σ > zero(σ)); new(μ, σ))
end

Levy{T<:Real}(μ::T, σ::T) = Levy{T}(μ, σ)
Levy(μ::Real, σ::Real) = Levy(promote(μ, σ)...)
Levy(μ::Integer, σ::Integer) = Levy(Float64(μ), Float64(σ))
Levy(μ::Real) = Levy(μ, 1.0)
Levy() = Levy(0.0, 1.0)

@distr_support Levy d.μ Inf

#### Conversions

convert{T <: Real, S <: Real}(::Type{Levy{T}}, μ::S, σ::S) = Levy(T(μ), T(σ))
convert{T <: Real, S <: Real}(::Type{Levy{T}}, d::Levy{S}) = Levy(T(d.μ), T(d.σ))

#### Parameters

location(d::Levy) = d.μ
params(d::Levy) = (d.μ, d.σ)
@inline partype{T<:Real}(d::Levy{T}) = T

#### Statistics

mean{T<:Real}(d::Levy{T}) = T(Inf)
var{T<:Real}(d::Levy{T}) = T(Inf)
skewness{T<:Real}(d::Levy{T}) = T(NaN)
kurtosis{T<:Real}(d::Levy{T}) = T(NaN)

mode(d::Levy) = d.σ / 3 + d.μ

entropy(d::Levy) = (1 - 3digamma(1) + log(16 * d.σ^2 * π)) / 2

median(d::Levy) = d.μ + d.σ / (2 * T(erfcinv(0.5))^2)

#### Evaluation

function pdf(d::Levy, x::Real)
    μ, σ = params(d)
    z = x - μ
    (sqrt(σ) / sqrt2π) * exp((-σ) / (2z)) / z^(3//2)
end

function logpdf(d::Levy, x::Real)
    μ, σ = params(d)
    z = x - μ
    (log(σ) - log2π - σ / z - 3log(z))/2
end

cdf(d::Levy, x::Real) = erfc(sqrt(d.σ / (2(x - d.μ))))
ccdf(d::Levy, x::Real) = erf(sqrt(d.σ / (2(x - d.μ))))

quantile(d::Levy, p::Real) = d.μ + d.σ / (2*erfcinv(p)^2)
cquantile(d::Levy, p::Real) = d.μ + d.σ / (2*erfinv(p)^2)

mgf{T<:Real}(d::Levy{T}, t::Real) = t == zero(t) ? one(T) : T(NaN)

function cf(d::Levy, t::Real)
    μ, σ = params(d)
    exp(im * μ * t - sqrt(-2im * σ * t))
end

#### Sampling

rand(d::Levy) = d.μ + d.σ / randn()^2

ext: jl


text:
# testing subviews

using ArrayViews
using Base.Test

#### Test Subviews

## tools to facilitate aview testing
function _test_arrview_contents(v, r)
    # siz_r = size(r)
    # siz_v = size(v)
    # if siz_r != siz_v
    #     error("Incorrect size: get $(siz_v), but expect $(siz_r)")
    # end

    if length(v) != length(r)
        error("Unmatched length: get $(length(v)), but expect $(length(r))")
    end

    for i = 1 : length(v)
        if v[i] != r[i]
            println("v = ")
            println(v)
            println("r = ")
            println(r)
            error("Incorrect content.")
        end
    end
end

function _test_arrview(a, r, subs...)
    _test_arrview_contents(aview(a, subs...), r)
    _test_arrview_contents(unsafe_aview(a, subs...), r)
end

macro test_arrview(a_, subs...)
    # ArrayViews still allows omitting indices for non-singleton trailing dims,
    # same as  < Julia 0.6.0, pre #23628.
    # For the sake of getting tests to pass on 0.7, test against old behavior
    esc(:(_test_arrview($a_, reshape($(a_),Val($(length(subs))))[$(subs...)], $(subs...))))
end

#### test views from arrays

avparent = copy(reshape(1.:1680., (8, 7, 6, 5)))

# 1D
println("    -- testing 1D views")

@test_arrview(avparent, 3)
@test_arrview(avparent, :)
@test_arrview(avparent, 1:12)
@test_arrview(avparent, 3:2:36)

# 2D
println("    -- testing 2D views")

@test_arrview(avparent, 4, 2)
@test_arrview(avparent, 4, :)
@test_arrview(avparent, 4, 3:10)
@test_arrview(avparent, 4, 2:2:10)

@test_arrview(avparent, :, 2)
@test_arrview(avparent, :, :)
@test_arrview(avparent, :, 3:10)
@test_arrview(avparent, :, 2:2:10)

@test_arrview(avparent, 1:6, 2)
@test_arrview(avparent, 1:6, :)
@test_arrview(avparent, 1:6, 3:10)
@test_arrview(avparent, 1:6, 2:2:10)

@test_arrview(avparent, 1:2:8, 2)
@test_arrview(avparent, 1:2:8, :)
@test_arrview(avparent, 1:2:8, 3:10)
@test_arrview(avparent, 1:2:8, 2:2:10)

# 3D
println("    -- testing 3D views")

@test_arrview(avparent, 4, 3, 2)
@test_arrview(avparent, 4, 3, :)
@test_arrview(avparent, 4, 3, 2:5)
@test_arrview(avparent, 4, 3, 1:2:6)

@test_arrview(avparent, 4, :, 2)
@test_arrview(avparent, 4, :, :)
@test_arrview(avparent, 4, :, 2:5)
@test_arrview(avparent, 4, :, 1:2:6)

@test_arrview(avparent, 4, 3:7, 2)
@test_arrview(avparent, 4, 3:7, :)
@test_arrview(avparent, 4, 3:7, 2:5)
@test_arrview(avparent, 4, 3:7, 1:2:6)

@test_arrview(avparent, 4, 1:2:5, 2)
@test_arrview(avparent, 4, 1:2:5, :)
@test_arrview(avparent, 4, 1:2:5, 2:5)
@test_arrview(avparent, 4, 1:2:5, 1:2:6)

@test_arrview(avparent, :, 3, 2)
@test_arrview(avparent, :, 3, :)
@test_arrview(avparent, :, 3, 2:5)
@test_arrview(avparent, :, 3, 1:2:6)

@test_arrview(avparent, :, :, 2)
@test_arrview(avparent, :, :, :)
@test_arrview(avparent, :, :, 2:5)
@test_arrview(avparent, :, :, 1:2:6)

@test_arrview(avparent, :, 3:7, 2)
@test_arrview(avparent, :, 3:7, :)
@test_arrview(avparent, :, 3:7, 2:5)
@test_arrview(avparent, :, 3:7, 1:2:6)

@test_arrview(avparent, :, 1:2:5, 2)
@test_arrview(avparent, :, 1:2:5, :)
@test_arrview(avparent, :, 1:2:5, 2:5)
@test_arrview(avparent, :, 1:2:5, 1:2:6)

@test_arrview(avparent, 2:7, 3, 2)
@test_arrview(avparent, 2:7, 3, :)
@test_arrview(avparent, 2:7, 3, 2:5)
@test_arrview(avparent, 2:7, 3, 1:2:6)

@test_arrview(avparent, 2:7, :, 2)
@test_arrview(avparent, 2:7, :, :)
@test_arrview(avparent, 2:7, :, 2:5)
@test_arrview(avparent, 2:7, :, 1:2:6)

@test_arrview(avparent, 2:7, 3:7, 2)
@test_arrview(avparent, 2:7, 3:7, :)
@test_arrview(avparent, 2:7, 3:7, 2:5)
@test_arrview(avparent, 2:7, 3:7, 1:2:6)

@test_arrview(avparent, 2:7, 1:2:5, 2)
@test_arrview(avparent, 2:7, 1:2:5, :)
@test_arrview(avparent, 2:7, 1:2:5, 2:5)
@test_arrview(avparent, 2:7, 1:2:5, 1:2:6)

@test_arrview(avparent, 1:2:7, 3, 2)
@test_arrview(avparent, 1:2:7, 3, :)
@test_arrview(avparent, 1:2:7, 3, 2:5)
@test_arrview(avparent, 1:2:7, 3, 1:2:6)

@test_arrview(avparent, 1:2:7, :, 2)
@test_arrview(avparent, 1:2:7, :, :)
@test_arrview(avparent, 1:2:7, :, 2:5)
@test_arrview(avparent, 1:2:7, :, 1:2:6)

@test_arrview(avparent, 1:2:7, 3:7, 2)
@test_arrview(avparent, 1:2:7, 3:7, :)
@test_arrview(avparent, 1:2:7, 3:7, 2:5)
@test_arrview(avparent, 1:2:7, 3:7, 1:2:6)

@test_arrview(avparent, 1:2:7, 1:2:5, 2)
@test_arrview(avparent, 1:2:7, 1:2:5, :)
@test_arrview(avparent, 1:2:7, 1:2:5, 2:5)
@test_arrview(avparent, 1:2:7, 1:2:5, 1:2:6)

# Some 4D Tests
println("    -- testing 4D views")

@test_arrview(avparent, 4, :,     3, 4)
@test_arrview(avparent, 4, :,     :, 4)
@test_arrview(avparent, 4, :,   3:5, 4)
@test_arrview(avparent, 4, :, 1:2:5, 4)

@test_arrview(avparent, :, :,     3, 4)
@test_arrview(avparent, :, :,     :, 4)
@test_arrview(avparent, :, :,   3:5, 4)
@test_arrview(avparent, :, :, 1:2:5, 4)

@test_arrview(avparent, 2:7, :,     3, 4)
@test_arrview(avparent, 2:7, :,     :, 4)
@test_arrview(avparent, 2:7, :,   3:5, 4)
@test_arrview(avparent, 2:7, :, 1:2:5, 4)

@test_arrview(avparent, 1:2:7, :,     3, 4)
@test_arrview(avparent, 1:2:7, :,     :, 4)
@test_arrview(avparent, 1:2:7, :,   3:5, 4)
@test_arrview(avparent, 1:2:7, :, 1:2:5, 4)

@test_arrview(avparent, 4, :,     3, :)
@test_arrview(avparent, 4, :,     :, :)
@test_arrview(avparent, 4, :,   3:5, :)
@test_arrview(avparent, 4, :, 1:2:5, :)

@test_arrview(avparent, :, :,     3, :)
@test_arrview(avparent, :, :,     :, :)
@test_arrview(avparent, :, :,   3:5, :)
@test_arrview(avparent, :, :, 1:2:5, :)

@test_arrview(avparent, 2:7, :,     3, :)
@test_arrview(avparent, 2:7, :,     :, :)
@test_arrview(avparent, 2:7, :,   3:5, :)
@test_arrview(avparent, 2:7, :, 1:2:5, :)

@test_arrview(avparent, 1:2:7, :,     3, :)
@test_arrview(avparent, 1:2:7, :,     :, :)
@test_arrview(avparent, 1:2:7, :,   3:5, :)
@test_arrview(avparent, 1:2:7, :, 1:2:5, :)

@test_arrview(avparent, 4, :,     3, 2:5)
@test_arrview(avparent, 4, :,     :, 2:5)
@test_arrview(avparent, 4, :,   3:5, 2:5)
@test_arrview(avparent, 4, :, 1:2:5, 2:5)

@test_arrview(avparent, :, :,     3, 2:5)
@test_arrview(avparent, :, :,     :, 2:5)
@test_arrview(avparent, :, :,   3:5, 2:5)
@test_arrview(avparent, :, :, 1:2:5, 2:5)

@test_arrview(avparent, 2:7, :,     3, 2:5)
@test_arrview(avparent, 2:7, :,     :, 2:5)
@test_arrview(avparent, 2:7, :,   3:5, 2:5)
@test_arrview(avparent, 2:7, :, 1:2:5, 2:5)

@test_arrview(avparent, 1:2:7, :,     3, 2:5)
@test_arrview(avparent, 1:2:7, :,     :, 2:5)
@test_arrview(avparent, 1:2:7, :,   3:5, 2:5)
@test_arrview(avparent, 1:2:7, :, 1:2:5, 2:5)

@test_arrview(avparent, :, 1, 1, 1)
@test_arrview(avparent, :, 1, 2, 1)
@test_arrview(avparent, :, :, 2, 2)
@test_arrview(avparent, :, :, 1:2, 3:4)

# Some 5D tests
avparent2 = copy(reshape(1.:6720., (8, 7, 6, 5, 4)))

@test_arrview(avparent2, :, 1, 1, 1, 1)
@test_arrview(avparent2, :, 1, 1, 3, 1)
@test_arrview(avparent2, :, 1, 2, 3, 3:4)
@test_arrview(avparent2, 2:7, 1, 1, 1, 3:4)

@test_arrview(avparent2, 1, :, 1, 3, 4)
@test_arrview(avparent2, 2, 3, :, 3, 4)
@test_arrview(avparent2, 3, 4, 5, :, 4)
@test_arrview(avparent2, 2, 3, 4, 4, :)

@test_arrview(avparent2, 1, 2:3, 3, 4, 2)
@test_arrview(avparent2, 2, 1, 3:4, 1, 1)
@test_arrview(avparent2, 3, 1, 1, 4:5, 2)
@test_arrview(avparent2, 4, 1, 1, 3:4, 2:4)

#### Test Subviews of Views

function print_subscripts(subs1, subs2)
    println("Error happens on: ")

    subs1s = join([repr(i) for i in subs1], ", ")
    subs2s = join([repr(i) for i in subs2], ", ")

    println("  subs1 = [$subs1s]")
    println("  subs2 = [$subs2s]")
end

function test_arrview2(a, subs1, subs2)
    v = aview(a, subs1...)
    v2r = aview(copy(v), subs2...)
    _test_arrview_contents(aview(v, subs2...), v2r)
    uv = unsafe_aview(a, subs1...)
    _test_arrview_contents(unsafe_aview(uv, subs2...), v2r)
end

avparent = copy(reshape(1:6912, (12, 12, 8, 6)))

# 1D --> 1D
println("    -- testing 1D sub-views of 1D views")

for sa in Any[Colon(), 1:36, 2:2:36]
    v1 = aview(avparent, sa)
    for sb in Any[4, Colon(), 1:length(v1), 3:2:length(v1)]
        test_arrview2(avparent, (sa,), (sb,))
    end
end

# 2D --> 2D
println("    -- testing 2D sub-views of 2D views")

for sa1 in Any[Colon(), 1:10, 2:2:12], sa2 = Any[Colon(), 1:12, 2:2:16]
    v1 = aview(avparent, sa1, sa2)
    for sb1 in Any[4, Colon(), 2:size(v1,1), 2:2:size(v1,1)],
        sb2 in Any[4, Colon(), 2:size(v1,2), 2:2:size(v1,2)]
        test_arrview2(avparent, (sa1, sa2), (sb1, sb2))
    end
end

# 3D --> 3D
println("    -- testing 3D sub-views of 3D views")

for sa1 in Any[Colon(), 1:10, 2:2:12],
    sa2 in Any[Colon(), 1:10, 2:2:12],
    sa3 in Any[Colon(), 1:7, 2:2:8]
    v1 = aview(avparent, sa1, sa2, sa3)
    (d1, d2, d3) = size(v1)
    for sb1 in Any[Colon(), 2:d1, 2:2:d1],
        sb2 in Any[Colon(), 2:d2, 2:2:d2],
        sb3 in Any[Colon(), 2:d3, 2:2:d3]
        test_arrview2(avparent, (sa1, sa2, sa3), (sb1, sb2, sb3))
    end
end


text:
using Mamba, Distributions

dogs = (Symbol => Any)[:Y => int(readdlm("dogs.wsv")) ];
dogs[:Dogs] = size(dogs[:Y], 1);
dogs[:Trials] = size(dogs[:Y], 2);
dogs[:xa] = mapslices(cumsum, dogs[:Y], 2);
dogs[:xs] = mapslices(x -> [1:25] - x, dogs[:xa], 2);
dogs[:y] = 1 - dogs[:Y][:, 2:25];

model = Model(
  y = Stochastic(2,
      @modelexpr(Dogs, Trials, alpha, xa, beta, xs,
      Distribution[
        begin
          p = exp(alpha * xa[i,j] + beta * xs[i,j])
          Bernoulli(p)
        end
        for i in 1:Dogs, j in 1:Trials-1]),
      false),
  alpha = Stochastic(:(Truncated(Flat(), -Inf, -1e-5))),
  A = Logical(@modelexpr(alpha,exp(alpha))),
  beta = Stochastic(:(Truncated(Flat(), -Inf, -1e-5))),
  B = Logical(@modelexpr(beta,exp(beta))) );

inits = [
  [:y => dogs[:y], :alpha => -1, :beta => -1],
  [:y => dogs[:y], :alpha => -2, :beta => -2]
];

scheme = [Slice([:alpha, :beta], [1.0, 1.0])];
setsamplers!(model, scheme);
sim =
  mcmc(model,dogs,inits,10000,burnin=2500,thin=2, chains=2);

max_stars_repo_path: Module 3/Chapter06/mcmc-dogs.jl

max_issues_repo_path: Module 3/Chapter06/mcmc-dogs.jl

max_forks_repo_path: Module 3/Chapter06/mcmc-dogs.jl


text:
using GLMakie, GeometryBasics, RPRMakie, RadeonProRender
using Colors, FileIO
using Colors: N0f8
RPR = RadeonProRender
GLMakie.set_window_config!(float=true, focus_on_show=false)
cat = load(GLMakie.assetpath("cat.obj"))

begin
    context = RPR.Context()

    matsys = RPR.MaterialSystem(context, 0)
    emissive = RPR.EmissiveMaterial(matsys)
    diffuse = RPR.DiffuseMaterial(matsys)

    fig = Figure(resolution=(1000, 1000))
    ax = LScene(fig[1, 1], scenekw=(show_axis=false,))
    for i in 4:4:12
        n = i + 1
        y = LinRange(0, i, n)
        y2 = (y ./ 2) .- 2
        lines!(ax, fill((i-5) ./ 2, n), y2, sin.(y) .+ 1, linewidth=5)
    end
    mesh!(ax, Rect3f(Vec3f(-3, -3, -0.1), Vec3f(6, 6, 0.1)), color=:white)
    mesh!(ax, Sphere(Point3f(0, 0, 2), 0.1), material=emissive)
    display(fig)
end
# fetch(task)
begin
    context, task = RPRMakie.replace_scene_rpr!(ax.scene, context, matsys)
    emissive.color = Vec3f(4, 2, 2)
end

begin

    uber1.color =Vec4f(1, 1, 1, 1)
    uber1.diffuse_weight =Vec4f(0, 0, 0, 0)
    uber1.diffuse_roughness =Vec4f(1, 1, 1, 1)
    uber1.reflection_color =Vec4f(0.996078, 0.858824, 0.639216, 0)
    uber1.reflection_weight =Vec4f(1, 1, 1, 1)
    uber1.reflection_roughness =Vec4f(0, 0, 0, 0)
    uber1.reflection_anisotropy =Vec4f(0, 0, 0, 0)
    uber1.reflection_anisotropy_rotation =Vec4f(0, 0, 0, 0)
    uber1.reflection_ior =Vec4f(1.36, 1.36, 1.36, 1.36)
    uber1.refraction_color =Vec4f(0.996078, 0.858824, 0.639216, 0)
    uber1.refraction_weight =Vec4f(1, 1, 1, 1)
    uber1.refraction_roughness =Vec4f(0, 0, 0, 0)
    uber1.refraction_ior =Vec4f(1.36, 1.36, 1.36, 1.36)
    uber1.refraction_absorption_color =Vec4f(0.996078, 0.858824, 0.639216, 0)
    uber1.refraction_absorption_distance =Vec4f(0, 0, 0, 0)
    uber1.refraction_caustics =Vec4f(0)
    uber1.coating_color =Vec4f(1, 1, 1, 1)
    uber1.coating_weight =Vec4f(0, 0, 0, 0)
    uber1.coating_roughness =Vec4f(0, 0, 0, 0)
    uber1.coating_ior =Vec4f(3, 3, 3, 3)
    uber1.coating_metalness =Vec4f(0, 0, 0, 0)
    uber1.coating_transmission_color =Vec4f(1, 1, 1, 1)
    uber1.coating_thickness =Vec4f(0, 0, 0, 0)
    uber1.sheen =Vec4f(1, 1, 1, 1)
    uber1.sheen_tint =Vec4f(0, 0, 0, 0)
    uber1.sheen_weight =Vec4f(0, 0, 0, 0)
    uber1.emission_color =Vec4f(1, 1, 1, 1)
    uber1.emission_weight =Vec3f( 0, 0, 0)
    uber1.transparency =Vec4f(0, 0, 0, 0)
    uber1.sss_scatter_color =Vec4f(0, 0, 0, 0)
    uber1.sss_scatter_distance =Vec4f(0, 0, 0, 0)
    uber1.sss_scatter_direction =Vec4f(0, 0, 0, 0)
    uber1.sss_weight =Vec4f(0, 0, 0, 0)
    uber1.backscatter_weight =Vec4f(0, 0, 0, 0)
    uber1.backscatter_color =Vec4f(1, 1, 1, 1)

    uber1.reflection_mode = UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR)
    uber1.emission_mode = UInt(RPR.RPR_UBER_MATERIAL_EMISSION_MODE_SINGLESIDED)
    uber1.coating_mode = UInt(RPR.RPR_UBER_MATERIAL_IOR_MODE_PBR)
    uber1.sss_multiscatter = false
    uber1.refraction_thin_surface = false
end

begin
    uber2.color = Vec4f(0.501961f0,0.0f0,0.0f0,0.0f0)
    uber2.diffuse_weight = Vec4f(1.0f0,1.0f0,1.0f0,1.0f0)
    uber2.diffuse_roughness = Vec4f(0.500000f0,0.500000f0,0.500000f0,0.500000f0)
    uber2.reflection_color = Vec4f(0.490196f0,0.490196f0,0.490196f0,0.0f0)
    uber2.reflection_weight = Vec4f(0.990000f0,0.990000f0,0.990000f0,0.990000f0)
    uber2.reflection_roughness = Vec4f(0.008000f0,0.008000f0,0.008000f0,0.008000f0)
    uber2.reflection_anisotropy = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.reflection_anisotropy_rotation = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.reflection_mode = 1
    uber2.reflection_ior = Vec4f(1.460000f0,1.460000f0,1.460000f0,1.460000f0)
    uber2.refraction_color = Vec4f(1.0f0,1.0f0,1.0f0,1.0f0)
    uber2.refraction_weight = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.refraction_roughness = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.refraction_ior = Vec4f(1.500000f0,1.500000f0,1.500000f0,1.500000f0)
    uber2.refraction_thin_surface = Vec4f(0)
    uber2.refraction_absorption_color = Vec4f(1.0f0,1.0f0,1.0f0,0.0f0)
    uber2.refraction_absorption_distance = Vec4f(1.0f0,1.0f0,1.0f0,0.0f0)
    uber2.refraction_caustics = 1
    uber2.coating_color = Vec4f(0.490196f0,0.490196f0,0.490196f0,0.0f0)
    uber2.coating_weight = Vec4f(1.0f0,1.0f0,1.0f0,1.0f0)
    uber2.coating_roughness = Vec4f(0.008000f0,0.008000f0,0.008000f0,0.008000f0)
    uber2.coating_mode = 1
    uber2.coating_ior = Vec4f(1.460000f0,1.460000f0,1.460000f0,1.460000f0)
    uber2.coating_metalness = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.coating_transmission_color = Vec4f(0.0f0,0.0f0,0.0f0,1.0f0)
    uber2.coating_thickness = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.sheen = Vec4f(1.0f0,1.0f0,1.0f0,1.0f0)
    uber2.sheen_tint = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.sheen_weight = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.emission_color = Vec4f(1.0f0,1.0f0,1.0f0,1.0f0)
    uber2.emission_weight = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.emission_mode = 1
    uber2.transparency = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.sss_scatter_color = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.sss_scatter_distance = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.sss_scatter_direction = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.sss_weight = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.sss_multiscatter = false
    uber2.backscatter_weight = Vec4f(0.0f0,0.0f0,0.0f0,0.0f0)
    uber2.backscatter_color = Vec4f(0.501961f0,0.0f0,0.0f0,0.0f0)
end

begin
    uber3.color = Vec4(0.752941,0.596078,0.443137,0.0)
    uber3.diffuse_weight = Vec4(1.0,1.0,1.0,1.0)
    uber3.diffuse_roughness = Vec4(0.500000,0.500000,0.500000,0.500000)
    uber3.reflection_color = Vec4(0.666667,0.490196,0.313726,0.0)
    uber3.reflection_weight = Vec4(1.0,1.0,1.0,1.0)
    uber3.reflection_roughness = Vec4(0.300000,0.300000,0.300000,0.300000)
    uber3.reflection_anisotropy = Vec4(0.0,0.0,0.0,0.0)
    uber3.reflection_anisotropy_rotation = Vec4(0.0,0.0,0.0,0.0)
    uber3.reflection_mode = 2
    uber3.reflection_ior = Vec4(1.0,1.0,1.0,1.0)
    uber3.refraction_color = Vec4(1.0,1.0,1.0,1.0)
    uber3.refraction_weight = Vec4(0.0,0.0,0.0,0.0)
    uber3.refraction_roughness = Vec4(0.0,0.0,0.0,0.0)
    uber3.refraction_ior = Vec4(1.500000,1.500000,1.500000,1.500000)
    uber3.refraction_thin_surface = 0
    uber3.refraction_absorption_color = Vec4(1.0,1.0,1.0,0.0)
    uber3.refraction_absorption_distance = Vec4(1.0,1.0,1.0,0.0)
    uber3.refraction_caustics = 1
    uber3.coating_color = Vec4(0.752941,0.596078,0.443137,0.0)
    uber3.coating_weight = Vec4(1.0,1.0,1.0,1.0)
    uber3.coating_roughness = Vec4(0.420000,0.420000,0.420000,0.420000)
    uber3.coating_mode = 1
    uber3.coating_ior = Vec4(1.700000,1.700000,1.700000,1.700000)
    uber3.coating_metalness = Vec4(0.0,0.0,0.0,0.0)
    uber3.coating_transmission_color = Vec4(0.0,0.0,0.0,1.0)
    uber3.coating_thickness = Vec4(0.0,0.0,0.0,0.0)
    uber3.sheen = Vec4(1.0,1.0,1.0,1.0)
    uber3.sheen_tint = Vec4(0.0,0.0,0.0,0.0)
    uber3.sheen_weight = Vec4(0.0,0.0,0.0,0.0)
    uber3.emission_color = Vec4(1.0,1.0,1.0,1.0)
    uber3.emission_weight = Vec4(0.0,0.0,0.0,0.0)
    uber3.emission_mode = 1
    uber3.transparency = Vec4(0.0,0.0,0.0,0.0)
    uber3.sss_scatter_color = Vec4(0.0,0.0,0.0,0.0)
    uber3.sss_scatter_distance = Vec4(0.0,0.0,0.0,0.0)
    uber3.sss_scatter_direction = Vec4(0.0,0.0,0.0,0.0)
    uber3.sss_weight = Vec4(0.0,0.0,0.0,0.0)
    uber3.sss_multiscatter = 0
    uber3.backscatter_weight = Vec4(0.0,0.0,0.0,0.0)
    uber3.backscatter_color = Vec4(0.752941,0.596078,0.443137,0.0)
end

begin
    uber4.color = Vec4f(1.0,1.0,1.0,1.0)
    uber4.diffuse_weight = Vec4f(0.0,0.0,0.0,0.0)
    uber4.diffuse_roughness = Vec4f(1.0,1.0,1.0,1.0)
    uber4.reflection_color = Vec4f(0.501961,0.501961,0.501961,0.0)
    uber4.reflection_weight = Vec4f(1.0,1.0,1.0,1.0)
    uber4.reflection_roughness = Vec4f(0.0,0.0,0.0,0.0)
    uber4.reflection_anisotropy = Vec4f(0.0,0.0,0.0,0.0)
    uber4.reflection_anisotropy_rotation = Vec4f(0.0,0.0,0.0,0.0)
    uber4.reflection_mode = 1
    uber4.reflection_ior = Vec4f(1.330000,1.330000,1.330000,1.330000)
    uber4.refraction_color = Vec4f(0.501961,0.898039,0.996078,0.0)
    uber4.refraction_weight = Vec4f(1.0,1.0,1.0,1.0)
    uber4.refraction_roughness = Vec4f(0.0,0.0,0.0,0.0)
    uber4.refraction_ior = Vec4f(1.330000,1.330000,1.330000,1.330000)
    uber4.refraction_thin_surface = 0
    uber4.refraction_absorption_color = Vec4f(0.501961,0.898039,0.996078,0.0)
    uber4.refraction_absorption_distance = Vec4f(0.0,0.0,0.0,0.0)
    uber4.refraction_caustics = 0
    uber4.coating_color = Vec4f(1.0,1.0,1.0,1.0)
    uber4.coating_weight = Vec4f(0.0,0.0,0.0,0.0)
    uber4.coating_roughness = Vec4f(0.0,0.0,0.0,0.0)
    uber4.coating_mode = 1
    uber4.coating_ior = Vec4f(3.0,3.0,3.0,3.0)
    uber4.coating_metalness = Vec4f(0.0,0.0,0.0,0.0)
    uber4.coating_transmission_color = Vec4f(1.0,1.0,1.0,1.0)
    uber4.coating_thickness = Vec4f(0.0,0.0,0.0,0.0)
    uber4.sheen = Vec4f(1.0,1.0,1.0,1.0)
    uber4.sheen_tint = Vec4f(0.0,0.0,0.0,0.0)
    uber4.sheen_weight = Vec4f(0.0,0.0,0.0,0.0)
    uber4.emission_color = Vec4f(1.0,1.0,1.0,1.0)
    uber4.emission_weight = Vec4f(0.0,0.0,0.0,0.0)
    uber4.emission_mode = 1
    uber4.transparency = Vec4f(0.0,0.0,0.0,0.0)
    uber4.sss_scatter_color = Vec4f(0.0,0.0,0.0,0.0)
    uber4.sss_scatter_distance = Vec4f(0.0,0.0,0.0,0.0)
    uber4.sss_scatter_direction = Vec4f(0.0,0.0,0.0,0.0)
    uber4.sss_weight = Vec4f(0.0,0.0,0.0,0.0)
    uber4.sss_multiscatter = 0
    uber4.backscatter_weight = Vec4f(0.0,0.0,0.0,0.0)
    uber4.backscatter_color = Vec4f(1.0,1.0,1.0,1.0)
end


text:
using .Graphs

export inftime_solve, Reorder, InftimeStateResult, state_cumtime, state_prob, mean_time_to_absorption
struct Reorder
    mat2chain::Vector{Int}
    chain2mat::Vector{Int}
    ntransients::Int
    recur_comp_member_count::Vector{Int}
end

function extend!(dst_array, new_array)
    push!(dst_array, new_array...)
end
"""
reorder states into transient states, recurrent components.
"""
function reorder_states(chain::ContMarkovChain)::Reorder
    comps = strongly_connected_components(chain.state_graph)
    mat2chain = Vector{Int}()
    ntrans = 0
    for comp in filter(comp -> !comp.is_bottom, comps)
        extend!(mat2chain, comp.members)
        ntrans += length(comp.members)
    end
    recur_comp_member_count = Vector{Int}()
    for comp in filter(comp -> comp.is_bottom, comps)
        extend!(mat2chain, comp.members)
        push!(recur_comp_member_count, length(comp.members))
    end
    chain2mat = Vector{Int}(undef, length(mat2chain))
    for i in 1:length(mat2chain)
        chain2mat[mat2chain[i]] = i
    end
    return Reorder(mat2chain, chain2mat, ntrans, recur_comp_member_count)
end

"""
    tt_rate_matrix(chain, mat2chain, chain2mat, ntrans)

create transition rate matrix (to-from order) between transient states.
"""
function tt_rate_matrix(chain, order)
    rows = Vector{Int}()
    cols = Vector{Int}()
    vals = Vector{Float64}()
    for src_mat in 1:order.ntransients
        src_chn = order.mat2chain[src_mat]
        for edge in out_edges(chain.state_graph, src_chn)
            dst_chn = dst_node(chain.state_graph, edge)
            dst_mat = order.chain2mat[dst_chn]
            rate = chain.trans_rates[edge]
            if dst_mat <= order.ntransients
                push!(rows, dst_mat)
                push!(cols, src_mat)
                push!(vals, rate)
            end
            push!(rows, src_mat)
            push!(cols, src_mat)
            push!(vals, -rate)
        end
    end
    sparse(rows, cols, vals, order.ntransients, order.ntransients)
end
"""
    ta_rate_matrix(chain, order, abs_start, nabs)

create transition rate matrix (to-from order) from transient states to
recurrent states within range abs_start, abs_start + abs
"""
function ta_rate_matrix(chain, order, abs_start, abs_end)
    rows = Vector{Int}()
    cols = Vector{Int}()
    vals = Vector{Float64}()
    for src_mat in 1:order.ntransients
        src_chn = order.mat2chain[src_mat]
        for edge in out_edges(chain.state_graph, src_chn)
            dst_chn = dst_node(chain.state_graph, edge)
            dst_mat = order.chain2mat[dst_chn]
            rate = chain.trans_rates[edge]
            if dst_mat >= abs_start && dst_mat <= abs_end
                push!(rows, dst_mat - abs_start + 1)
                push!(cols, src_mat)
                push!(vals, rate)
            end
        end
    end
    sparse(rows, cols, vals, abs_end - abs_start + 1, order.ntransients)
end

"""
    aa_rate_matrix(chain, order, abs_start, nabs)

create transition rate matrix (to-from order) from between
recurrent states within range abs_start, abs_start + abs
"""
function aa_rate_matrix(chain, order, abs_start, abs_end)
    rows = Vector{Int}()
    cols = Vector{Int}()
    vals = Vector{Float64}()
    for src_mat in abs_start:abs_end
        src_chn = order.mat2chain[src_mat]
        total_rate = 0.0
        for edge in out_edges(chain.state_graph, src_chn)
            dst_chn = dst_node(chain.state_graph, edge)
            dst_mat = order.chain2mat[dst_chn]
            rate = chain.trans_rates[edge]
            total_rate += rate
            if dst_mat != abs_end
                push!(rows, dst_mat - abs_start + 1)
                push!(cols, src_mat - abs_start + 1)
                push!(vals, rate)
            end
        end
        if src_mat != abs_end
            push!(rows, src_mat - abs_start + 1)
            push!(cols, src_mat - abs_start + 1)
            push!(vals, -total_rate)
        end
        push!(rows, abs_end - abs_start + 1)
        push!(cols, src_mat - abs_start + 1)
        push!(vals, 1.0)
    end
    sparse(rows, cols, vals, abs_end - abs_start + 1, abs_end - abs_start + 1)
end

struct InftimeStateResult
    reorder::Reorder
    solution::Vector{Float64}
end

"""
    state_prob(result, state)

used to retrieve state probablity from result.
"""
function state_prob(res::InftimeStateResult, state)
    state_mat = res.reorder.chain2mat[state]
    if state_mat <= res.reorder.ntransients
        0.0
    else
        res.solution[state_mat]
    end
end

"""
    state_cumtime(result, state)

used to retrieve state cumulative time from result.
"""
function state_cumtime(res::InftimeStateResult, state)
    state_mat = res.reorder.chain2mat[state]
    if state_mat <= res.reorder.ntransients
        res.solution[state_mat]
    else
        Inf
    end
end
"""
    inftime_state(chain, init_prob; spsolve)

compute state cumulative times/probabilites of the markov chain at time infinity.
`state_prob` and `state_cumtime` can be used to retrieve times/probs from the return value.
"""
function inftime_solve(chain::ContMarkovChain, init_prob; spsolve=Base.:\)
    order = reorder_states(chain)
    sol::Vector{Float64} = fill(0.0, state_count(chain))
    for idx in eachindex(init_prob)
        val = init_prob[idx]
        mat_idx = order.chain2mat[idx]
        if mat_idx <= order.ntransients
            sol[mat_idx] = -val
        else
            sol[mat_idx] = val
        end
    end

    #solve Qₜₜ . τ = -π(0)
    if order.ntransients > 0
        QTT = tt_rate_matrix(chain, order)
        sol[1:order.ntransients] = spsolve(QTT, sol[1:order.ntransients])
    end

    abs_start = order.ntransients + 1
    for nabs in order.recur_comp_member_count
        abs_end = abs_start + nabs - 1
        if order.ntransients > 0
            QTA = ta_rate_matrix(chain, order, abs_start, abs_end)
            sol[abs_start:abs_end] += QTA * sol[1:order.ntransients]
        end
        total_prob = sum(sol[abs_start:abs_end])
        QAA = aa_rate_matrix(chain, order, abs_start, abs_end)
        b = fill(0.0, nabs)
        b[nabs] = total_prob
        sol[abs_start:abs_end] = spsolve(QAA, b)
        abs_start = abs_end + 1
    end
    return InftimeStateResult(order, sol)
end

function mean_time_to_absorption(sol::InftimeStateResult)
    g_start = sol.reorder.ntransients + 1
    for g in sol.reorder.recur_comp_member_count
        if g > 1
            for mat_idx in 1:g
                if sol.solution[mat_idx] > 0
                    return Inf
                end
            end
        end
    end
    sum(sol.solution[1:sol.reorder.ntransients])
end

function mean_time_to_absorption(chain::ContMarkovChain, init_prob)
    sol = inftime_solve(chain::ContMarkovChain, init_prob)
    mean_time_to_absorption(sol)
end


text:
import Flamenco

"""
    close_server()

Stop Flamenco server
"""
function close_server()
    Flamenco.close_server()
end

"""
    respond(; headers::Dict{String, String}, status_code::Int16, body::String) :: String

Return response for specified parameters
"""
respond(; headers::Dict{String, String}, status_code::Int16, body::String) = Flamenco.write_response(status_code, headers; body=body)

"""
    start_server(server_host::String, server_port::Int64) :: Nothing

Start Flamenco server and listen for incoming requests and return HTML or JSON body
"""
function start_server(server_host::String, server_port::Int64) :: Nothing
    Flamenco.start_server(server_host, server_port) do request::Flamenco.Server.Request
        render(; request_headers=request.headers, request_method=request.method, request_path=request.target, request_payload=request.body)
    end
end


text:
# Modelo de tanque de armazenamento com tres alimentacoes e perdas
type Tanque3
	Tanque3()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in1}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in2}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in3}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{perda}"
			)),
			[
				:(Fout.F = Fin1.F + Fin2.F + Fin3.F - Fperda.F),
				:(Fout.F * Fout.H = Fin1.F * Fin1.H + Fin2.F * Fin2.H + Fin3.F * Fin3.H - Fperda.F*Fperda.H),
				:([Fout.S,Fout.T] = PP2.propPH(Fout.P,Fout.H)),
				:(Fperda.S = Fout.S),
				:(Fperda.T = Fout.T),
				:(Fperda.P = Fout.P),
				:(Fperda.H = Fout.H),
			],
			[
				"","","","","","","",
			],
			[:PP2,],
			[:Fin1,:Fin2,:Fin3,:Fout,:Fperda,]
		)
	end
	PP2::DanaPlugin
	Fin1::Corrente
	Fin2::Corrente
	Fin3::Corrente
	Fout::Corrente
	Fperda::Corrente
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
end
export Tanque3
function setEquationFlow(in::Tanque3)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
end


text:
handle("ping") do
  "pong"
end

handle("echo") do data
  data
end

handle("cd") do path
  cd(path)
end

handle("evalsimple") do code
  eval(parse(code))
end

handle("exit") do
  exit()
end

handle("clear-workspace") do
  # has to be run by Main otherwise it throws an error
  eval(Main, :(workspace()))
end


text:
# Testing modulation filterbank response
x = [1, zeros(9999)]
mf = compute_modulation_cfs(4,128,8)
mfb = modulation_filterbank(mf, 16000, 2)
y = filt(mfb, x)

mf_matlab = readcsv(open(joinpath(dirname(@__FILE__), "data", "modfilters_freqs.csv")))
y_matlab = readcsv(open(joinpath(dirname(@__FILE__), "data", "modfilters_response.csv")))

@test_approx_eq y y_matlab'
@test_approx_eq mf mf_matlab


text:
abstract type AbstractCurve2D end
abstract type AbstractCircularLinearCurve <: AbstractCurve2D end

abstract type AbstractLinearCurve <: AbstractCircularLinearCurve end
abstract type AbstractLine <: AbstractLinearCurve end
abstract type AbstractRay <: AbstractLinearCurve end
abstract type AbstractLineSegment <: AbstractLinearCurve end

abstract type AbstractCircularCurve <: AbstractCircularLinearCurve end
abstract type AbstractCircle <: AbstractCircularCurve end
abstract type AbstractArc <: AbstractCircularCurve end


text:
using Test
using Distributed
using LinearAlgebra
using NPZ
using Clustering

@everywhere using Random
@everywhere using DPMMSubClustersStreaming
include("multinomial_tests.jl")
include("niw_tests.jl")
include("unitests.jl")
include("module_tests.jl")

# addprocs(2)
# @everywhere using Random
# @everywhere using DPMMSubClustersStreaming

# include("multinomial_tests.jl")
# include("niw_tests.jl")
# include("unitests.jl")
# include("module_tests.jl")


text:
"""
    macro @agent(dims, varargs...)

Basic Macro to create an instance of Agent. It generates an Agent type with the characteristics specified in its arguments.

Check the [documentation](https://dsb-lab.github.io/AgentBasedModels.jl/dev/Usage.html#Defining-Agent-Properties) for details of the possible varargs.
"""
macro agent(dims, varargs...)

    boundaryDeclared = false

    m = Agent()
    m.dims = dims

    if dims == 0
        nothing
    elseif dims == 1
        push!(m.declaredSymbols["Local"],:x)
    elseif dims == 2
        push!(m.declaredSymbols["Local"],:x,:y)
    elseif dims == 3
        push!(m.declaredSymbols["Local"],:x,:y,:z)
    else
        error("Dims has to be an integer between 0 and 3.")
    end
    push!(m.declaredSymbols["Identity"],:id)

    #Add all contributions
    for ii in varargs
        if typeof(ii) != Expr
            error(ii, " should be and expression containing at least a name and a type. Example variable::Variable.")
        end

        if ii.head == :call
            t = eval(ii)
            if typeof(t) != Expr
                error("Function ", ii, " should return a piece of code to be added to the agent.")
            elseif t.head != :block
                t = quote $t end
            end

            t = t.args
        else
            t = [ii]
        end

        for i in t
            if typeof(i) == LineNumberNode
                nothing
            elseif i.head == :(::)
                name = i.args[1]
                type = i.args[2]

                if !(type in VALID_TYPES)
                    error(i, " should be and expression containing at least a name and a type. Example variable::Variable.")
                elseif typeof(name) == Expr
                    if !(name.head == :vect)
                        error(name, " should be a symbol or an array of symbols. Check how to declare ", type, " types.")
                    else
                        names = Symbol[]
                        for j in name.args
                            if typeof(j) != Symbol
                                error(j, " in ", i, " should be a symbol.  Check how to declare ", type, " types.")
                            end
                            push!(names,j)
                        end
                        name = names
                    end
                end

                checkDeclared_(m,name)
                if typeof(name) == Array{Symbol,1}
                    append!(m.declaredSymbols[string(type)],name)
                else
                    push!(m.declaredSymbols[string(type)],name)
                end

            elseif i.head == :(=)
                if i.args[1] in VALID_UPDATES
                    if !(string(i.args[1]) in keys(m.declaredUpdates))
                        if i.args[2].head == :block
                            m.declaredUpdates[string(i.args[1])] = i.args[2]
                        else
                            m.declaredUpdates[string(i.args[1])] = quote $(i.args[2]) end
                        end
                    elseif i.args[2].head == :block
                        append!(m.declaredUpdates[string(i.args[1])].args,i.args[2].args)
                    end
                elseif i.args[1] == :(Boundary)
                    if boundaryDeclared
                        error("Boundary can only be declared once.")
                    end
                    boundary = Main.eval(i.args[2])
                    if !(typeof(boundary)<:Boundary)
                        error("Boundary has to be declared with a type boundary.")
                    elseif dims != boundary.dims
                        error("Boundary must have the same dimensions as the agent.")
                    else
                        m.boundary = boundary
                        boundaryDeclared = true
                    end
                else
                    error(i.args[1], " is not a valid type.")
                end
            else
                error(i, " is not an understood rule or variable declaration. Error in ", ii)
            end
        end
    end

    if !boundaryDeclared && !isempty(m.declaredSymbols["Medium"])
        error("If medium variables are declared, a boundary must be explicitely declared.")
    elseif !boundaryDeclared
        m.boundary = BoundaryFlat(dims)
    end

    return m
end


text:
struct VFE end
const DTC = VFE

struct ApproxPosteriorGP{Tapprox, Tprior, Tdata} <: AbstractGP
    approx::Tapprox
    prior::Tprior
    data::Tdata
end

"""
    approx_posterior(::VFE, fx::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)

Compute the optimal approximate posterior [1] over the process `f`, given observations `y`
of `f` at `x`, and inducing points `u`, where `u = f(z)` for some inducing inputs `z`.

[1] - M. K. Titsias. "Variational learning of inducing variables in sparse Gaussian
processes". In: Proceedings of the Twelfth International Conference on Artificial
Intelligence and Statistics. 2009.
"""
function approx_posterior(::VFE, fx::FiniteGP, y::AbstractVector{<:Real}, u::FiniteGP)
    U_y = cholesky(Symmetric(fx.Σy)).U
    U = cholesky(Symmetric(cov(u))).U

    B_εf = U' \ (U_y' \ cov(fx, u))'

    b_y = U_y' \ (y - mean(fx))

    D = B_εf * B_εf' + I
    Λ_ε = cholesky(Symmetric(D))

    m_ε = Λ_ε \ (B_εf * b_y)

    cache = (
        m_ε=m_ε,
        Λ_ε=Λ_ε,
        U=U,
        α=U \ m_ε,
        z=u.x,
        b_y=b_y,
        B_εf=B_εf,
        x=fx.x,
        Σy=fx.Σy,
    )
    return ApproxPosteriorGP(VFE(), fx.f, cache)
end

"""
    function update_approx_posterior(
        f_post_approx::ApproxPosteriorGP,
        fx::FiniteGP,
        y::AbstractVector{<:Real}
    )

Update the `ApproxPosteriorGP` given a new set of observations. Here, we retain the same
of pseudo-points.
"""
function update_approx_posterior(
    f_post_approx::ApproxPosteriorGP,
    fx::FiniteGP,
    y::AbstractVector{<:Real}
)
    U = f_post_approx.data.U
    z = f_post_approx.data.z

    U_y₂ = cholesky(Symmetric(fx.Σy)).U

    temp = zeros(size(f_post_approx.data.Σy, 1), size(fx.Σy, 2))
    Σy = [f_post_approx.data.Σy temp; temp' fx.Σy]

    b_y = vcat(f_post_approx.data.b_y, U_y₂ \ (y - mean(fx)))

    B_εf₂ = U' \ (U_y₂' \ cov(fx.f, fx.x, z))'
    B_εf = hcat(f_post_approx.data.B_εf, B_εf₂)

    Λ_ε = f_post_approx.data.Λ_ε

    for col in eachcol(B_εf₂)
        lowrankupdate!(Λ_ε, col)
    end

    m_ε = Λ_ε \ (B_εf * b_y)
    α = U \ m_ε
    x = vcat(f_post_approx.data.x, fx.x)

    cache = (
        m_ε=m_ε,
        Λ_ε=Λ_ε,
        U=U,
        α=α,
        z=z,
        b_y=b_y,
        B_εf=B_εf,
        x=x,
        Σy=Σy,
    )
    return ApproxPosteriorGP(VFE(), fx.f, cache)
end

"""
    function update_approx_posterior(
        f_post_approx::ApproxPosteriorGP,
        u::FiniteGP,
    )

Update the `ApproxPosteriorGP` given a new set of pseudo-points to append to the existing
set of pseudo points.
"""
function update_approx_posterior(
    f_post_approx::ApproxPosteriorGP,
    u::FiniteGP,
)
    U11 = f_post_approx.data.U
    C12 = cov(u.f, f_post_approx.data.z, u.x)
    C22 = Symmetric(cov(u))
    U = update_chol(Cholesky(U11,'U', 0), C12, C22).U
    U22 = U[end-length(u)+1:end, end-length(u)+1:end]
    U12 = U[1:length(f_post_approx.data.z), end-length(u)+1:end]

    B_εf₁ = f_post_approx.data.B_εf

    Cu1f = cov(f_post_approx.prior, f_post_approx.data.z, f_post_approx.data.x)
    Cu2f = cov(f_post_approx.prior, u.x, f_post_approx.data.x)

    U_y = cholesky(Symmetric(f_post_approx.data.Σy)).U

    B_εf₂ = U22' \ (Cu2f * inv(U_y)   - U12' * B_εf₁)
    B_εf = vcat(B_εf₁, B_εf₂)

    Λ_ε = update_chol(f_post_approx.data.Λ_ε, B_εf₁ * B_εf₂', B_εf₂ * B_εf₂' + I)

    m_ε = Λ_ε \ (B_εf * f_post_approx.data.b_y)

    α = U \ m_ε

    z = vcat(f_post_approx.data.z, u.x)

    cache = (
        m_ε=m_ε,
        Λ_ε=Λ_ε,
        U=U,
        α=α,
        z=z,
        b_y=f_post_approx.data.b_y,
        B_εf=B_εf,
        x=f_post_approx.data.x,
        Σy=f_post_approx.data.Σy,
    )
    return ApproxPosteriorGP(VFE(), f_post_approx.prior, cache)
end

LinearAlgebra.Symmetric(X::Diagonal) = X

# AbstractGP interface implementation.

function Statistics.mean(f::ApproxPosteriorGP{VFE}, x::AbstractVector)
    return mean(f.prior, x) + cov(f.prior, x, f.data.z) * f.data.α
end

function Statistics.cov(f::ApproxPosteriorGP{VFE}, x::AbstractVector)
    A = f.data.U' \ cov(f.prior, f.data.z, x)
    return cov(f.prior, x) - At_A(A) + Xt_invA_X(f.data.Λ_ε, A)
end

function cov_diag(f::ApproxPosteriorGP{VFE}, x::AbstractVector)
    A = f.data.U' \ cov(f.prior, f.data.z, x)
    return cov_diag(f.prior, x) - diag_At_A(A) + diag_Xt_invA_X(f.data.Λ_ε, A)
end

function Statistics.cov(f::ApproxPosteriorGP{VFE}, x::AbstractVector, y::AbstractVector)
    A_zx = f.data.U' \ cov(f.prior, f.data.z, x)
    A_zy = f.data.U' \ cov(f.prior, f.data.z, y)
    return cov(f.prior, x, y) - A_zx'A_zy + Xt_invA_Y(A_zx, f.data.Λ_ε, A_zy)
end

function mean_and_cov(f::ApproxPosteriorGP{VFE}, x::AbstractVector)
    A = f.data.U' \ cov(f.prior, f.data.z, x)
    m_post = mean(f.prior, x) + A' * f.data.m_ε
    C_post = cov(f.prior, x) - At_A(A) + Xt_invA_X(f.data.Λ_ε, A)
    return m_post, C_post
end

function mean_and_cov_diag(f::ApproxPosteriorGP{VFE}, x::AbstractVector)
    A = f.data.U' \ cov(f.prior, f.data.z, x)
    m_post = mean(f.prior, x) + A' * f.data.m_ε
    c_post = cov_diag(f.prior, x) - diag_At_A(A) + diag_Xt_invA_X(f.data.Λ_ε, A)
    return m_post, c_post
end


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5.3

# Short description of this file: Poisson's 2D Equations with GPU

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Plots
using CUDA

@parameters x y θ
@variables u(..)
@derivatives Dxx''~x
@derivatives Dyy''~y

# 2D PDE
eq  = Dxx(u(x,y,θ)) + Dyy(u(x,y,θ)) ~ sin(2π*x)*cos(2π*y)

# Boundary conditions
bcs = [u(0,y,θ) ~ 0.f0, u(1,y,θ) ~ -sin(pi*1)*sin(pi*y),
       u(x,0,θ) ~ 0.f0, u(x,1,θ) ~ -sin(pi*x)*sin(pi*1)]
# Space and time domains
domains = [x ∈ IntervalDomain(0.0,1.0),
           y ∈ IntervalDomain(0.0,1.0)]

# Neural network
dim = 2 # number of dimensions
chain = Chain(Dense(dim,16,Flux.σ),Dense(16,16,Flux.σ),Dense(16,1))

# Discretization
dx = 0.05
discretization = PhysicsInformedNN(dx,
                                   chain,
                                   strategy = GridTraining())

pde_system = PDESystem(eq,bcs,domains,[x,y],[u])
prob = discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

res = GalacticOptim.solve(prob, BFGS(); cb = cb, maxiters=1000)
phi = discretization.phi

xs,ys = [domain.domain.lower:dx/10:domain.domain.upper for domain in domains]
analytic_sol_func(x,y) = -sin(2π*x)*cos(2π*y)/(8π^2)

u_predict = reshape([first(phi([x,y],res.minimizer)) for x in xs for y in ys],(length(xs),length(ys)))
u_real = reshape([analytic_sol_func(x,y) for x in xs for y in ys], (length(xs),length(ys)))
diff_u = abs.(u_predict .- u_real)

p1 = Plots.plot(xs, ys, u_real, linetype=:contourf,title = "analytic");
p2 = Plots.plot(xs, ys, u_predict, linetype=:contourf,title = "predict");
p3 = Plots.plot(xs, ys, diff_u,linetype=:contourf,title = "error");
plot(p1,p2,p3)

p1 = Plots.plot(xs, ys, u_real, linetype=:surface,title = "analytic");
p2 = Plots.plot(xs, ys, u_predict, linetype=:surface,title = "predict");
p3 = Plots.plot(xs, ys, diff_u,linetype=:surface,title = "error");
plot(p1,p2,p3)
savefig("graph.pdf")


text:
module Test10

import Automa
import Automa.RegExp: @re_str
const re = Automa.RegExp
using Test

@testset "Test10" begin
    machine = Automa.compile(re.primitive(0x61))
    @test Automa.execute(machine, "a")[1] == 0
    @test Automa.execute(machine, "b")[1] < 0

    machine = Automa.compile(re.primitive(0x61:0x62))
    @test Automa.execute(machine, "a")[1] == 0
    @test Automa.execute(machine, "b")[1] == 0
    @test Automa.execute(machine, "c")[1] < 0

    machine = Automa.compile(re.primitive('a'))
    @test Automa.execute(machine, "a")[1] == 0
    @test Automa.execute(machine, "b")[1] < 0

    machine = Automa.compile(re.primitive('樹'))
    @test Automa.execute(machine, "樹")[1] == 0
    @test Automa.execute(machine, "儒")[1] < 0

    machine = Automa.compile(re.primitive("ジュリア"))
    @test Automa.execute(machine, "ジュリア")[1] == 0
    @test Automa.execute(machine, "パイソン")[1] < 0

    machine = Automa.compile(re.primitive([0x61, 0x62, 0x72]))
    @test Automa.execute(machine, "abr")[1] == 0
    @test Automa.execute(machine, "acr")[1] < 0

    machine = Automa.compile(re"[^A-Z]")
    @test Automa.execute(machine, "1")[1] == 0
    @test Automa.execute(machine, "A")[1] < 0
    @test Automa.execute(machine, "a")[1] == 0

    machine = Automa.compile(re"[A-Z]+" & re"FOO?")
    @test Automa.execute(machine, "FO")[1] == 0
    @test Automa.execute(machine, "FOO")[1] == 0
    @test Automa.execute(machine, "foo")[1] < 0

    machine = Automa.compile(re"[A-Z]+" \ re"foo")
    @test Automa.execute(machine, "FOO")[1] == 0
    @test Automa.execute(machine, "foo")[1] < 0

    machine = Automa.compile(!re"foo")
    @test Automa.execute(machine, "bar")[1] == 0
    @test Automa.execute(machine, "foo")[1] < 0

    machine = Automa.compile(re.rep("a"))
    @test Automa.execute(machine, "aaa")[1] == 0
    @test Automa.execute(machine, "aYa")[1] < 0
end

end


text:
#TODO: Fix SMatrix for nonnormal incidence
function smatrix(sim::Simulation{1},
            ω::AbstractVector;
            ky::Real=0,
            kz::Real=0,
            lupack::AbstractLUPACK=DEFAULT_LUPACK)

    spc = SMatrixContainer(sim,ω,lupack,float(ky),float(kz))
    S = SharedArray{ComplexF64,3}((ndims(spc),ndims(spc),length(ω)); init=spc)
    return ScatteringMatrix{ndims(spc),typeof(ω)}(sdata(S),ω)
end

function SMatrixContainer(sim::Simulation{1},ω,lupack,ky,kz)
    N = sum(_which_channels(sim))
    N>0 || throw("no scattering for closed system")
    return SMatrixContainer{N,typeof(sim),typeof(ω),typeof(lupack)}(sim,ω,lupack,ky,kz)
end

"""
    _which_channels(simulation)

Boolean vector indicating which channels are available for scattering.
"""
function _which_channels(sim::Simulation{1})
    bc1,bc2 = sim.boundary.bcs
    bl1,bl2 = sim.boundary.bls
    bch1, bch2 = BCHermiticity(bc1), BCHermiticity(bc2)
    a = falses(2)
    if typeof(bl1)<:AbstractComplexBL || bch1<:NonHermitianBC
        a[1] = true
    end
    if typeof(bl2)<:AbstractComplexBL || bch2<:NonHermitianBC
        a[2] = true
    end
    return a
end

# One-sided, one-dimensional s-matrix core
@inline function dimensional_smatrix!(S::SharedArray, spc::SMatrixContainer{1,TSIM}, indices) where TSIM<:Simulation{1}
    ls = MaxwellLS(spc.sim, mean(spc.ω[indices]), 1, 0)
    a = _which_channels(spc.sim)

    # @fastmath @inbounds @simd
    for i ∈ indices
        if a[1]
            update!(ls,1,0,spc.ω[i])
            alu = lu(ls.maxwell.A,spc.lupack)
            scattering!(ls,alu)
            S[1,1,i] = ls.solution.total.left[2]*sqrt(spc.ω[i])
        elseif a[2]
            update!(ls,0,1,spc.ω[i])
            alu = lu(ls.maxwell.A,spc.lupack)
            scattering!(ls,alu)
            S[1,1,i] = ls.solution.total.right[2]*sqrt(spc.ω[i])
        end
    end
    return nothing
end

# Two-sided, one-dimensional smatrix core
@inline function dimensional_smatrix!(S::SharedArray, spc::SMatrixContainer{2,TSIM}, indices) where TSIM<:Simulation{1}
    ls = MaxwellLS(spc.sim, mean(spc.ω[indices]), 1, 0)
    a = _which_channels(spc.sim)

    # @fastmath @inbounds @simd
    for i ∈ indices
        update!(ls,1,0,spc.ω[i])
        alu = lu(ls.maxwell.A,spc.lupack)
        scattering!(ls,alu)
        S[1,1,i] = ls.solution.scattered.left[2]*sqrt(ls.equivalent_source.channelflux[1])
        S[1,2,i] = ls.solution.total.right[2]*sqrt(ls.equivalent_source.channelflux[2])

        update!(ls,0,1)
        scattering!(ls,alu)
        S[2,1,i] = ls.solution.total.left[2]*sqrt(ls.equivalent_source.channelflux[1])
        S[2,2,i] = ls.solution.scattered.right[2]*sqrt(ls.equivalent_source.channelflux[2])
    end
    return nothing
end


text:
function diffeq(du,u,p,t)
    du[V.M] = p[C.kM]*p[C.S]/(p[C.KS]+p[C.S]) - p[C.dM]*u[V.M]
    du[V.E] = p[C.kE]*(u[V.M]/(p[C.KM]+u[V.M]))*(u[V.E]/(p[C.KE]+u[V.E])) + p[C.kb]*u[V.M]/(p[C.KM]+u[V.M]) +
                p[C.kP1]*u[V.CD]*u[V.RE]/(p[C.KCD]+u[V.RE]) + p[C.kP2]*u[V.CE]*u[V.RE]/(p[C.KCE]+u[V.RE]) - p[C.dE]*u[V.E]-p[C.kRE]*u[V.R]*u[V.E]
    du[V.CD] = p[C.kCD]*u[V.M]/(p[C.KM]+u[V.M]) + p[C.kCDS]*p[C.S]/(p[C.KS]+p[C.S]) - p[C.dCD]*u[V.CD]
    du[V.CE] = p[C.kCE]*u[V.E]/(p[C.KE]+u[V.E]) - p[C.dCE]*u[V.CE]
    du[V.R] = p[C.kR] + p[C.kDP]*u[V.RP]/(p[C.KRP]+u[V.RP]) - p[C.kRE]*u[V.R]*u[V.E] - p[C.kP1]*u[V.CD]*u[V.R]/(p[C.KCD]+u[V.R]) -
                p[C.kP2]*u[V.CE]*u[V.R]/(p[C.KCE]+u[V.R]) - p[C.dR]*u[V.R]
    du[V.RP] = p[C.kP1]*u[V.CD]*u[V.R]/(p[C.KCD]+u[V.R]) + p[C.kP2]*u[V.CE]*u[V.R]/(p[C.KCE]+u[V.R]) + p[C.kP1]*u[V.CD]*u[V.RE]/(p[C.KCD]+u[V.RE]) +
                p[C.kP2]*u[V.CE]*u[V.RE]/(p[C.KCE]+u[V.RE]) - p[C.kDP]*u[V.RP]/(p[C.KRP]+u[V.RP]) - p[C.dRP]*u[V.RP]
    du[V.RE] = p[C.kRE]*u[V.R]*u[V.E]-p[C.kP1]*u[V.CD]*u[V.RE]/(p[C.KCD]+u[V.RE]) - p[C.kP2]*u[V.CE]*u[V.RE]/(p[C.KCE]+u[V.RE]) - p[C.dRE]*u[V.RE]

end

function get_steady_state(p::Vector{Float64})
    tspan::Tuple{Float64,Float64} = (0.0,Inf)
    u0::Vector{Float64} = zeros(V.NUM)

    prob = ODEProblem(diffeq,u0,tspan,p)
    prob = SteadyStateProblem(prob)
    sol = solve(prob,DynamicSS(CVODE_BDF()),dt=1.0)

    return sol.u
end

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
using FileIO
using Test

@testset "Mime save" begin

function Base.show(io::IO, m::MIME"image/svg+xml", data::MimeSaveTestType)
    content = read(joinpath(@__DIR__, "files", "mimesavetest.svg"))
    write(io, content)
end

function Base.show(io::IO, m::MIME"application/pdf", data::MimeSaveTestType)
    content = read(joinpath(@__DIR__, "files", "mimesavetest.pdf"))
    write(io, content)
end

function Base.show(io::IO, m::MIME"application/eps", data::MimeSaveTestType)
    content = read(joinpath(@__DIR__, "files", "mimesavetest.eps"))
    write(io, content)
end

function Base.show(io::IO, m::MIME"image/png", data::MimeSaveTestType)
    content = read(joinpath(@__DIR__, "files", "mimesavetest.png"))
    write(io, content)
end

data = MimeSaveTestType()

output_filename = tempname()

for filetype in [".svg", ".pdf", ".eps", ".png"]

    try
        save(output_filename * filetype, data)

        content_original = read(joinpath(@__DIR__, "files", "mimesavetest$filetype"))
        content_new = read(output_filename * filetype)

        @test content_new == content_original
    finally
        isfile(output_filename * filetype) && rm(output_filename * filetype)
    end

end

end


text:
struct ReferenceRecord
    # TODO: change it to `base` and `bound`
    first::UInt
    last::UInt
    metadata::Any
end

const REFERENCES = Vector{ReferenceRecord}[]

function init_reference()
    if length(REFERENCES) != Threads.nthreads()
        append!(empty!(REFERENCES), (ReferenceRecord[] for _ in 1:Threads.nthreads()))
    end
    clear_reference()
end

function clear_reference()
    foreach(empty!, REFERENCES)
end

function MemoryAccesses.reference(a::AbstractArray, metadata = nothing)
    f = UInt(pointer(a, firstindex(a)))
    l = UInt(pointer(a, lastindex(a)))
    rec = ReferenceRecord(f, l, metadata)
    push!(REFERENCES[Threads.threadid()], rec)
    return
end


text:
using SexSpecificMigrations
using JLD

sttime = Dates.format(now(), "yyyy-m-d-H-M")
hostname = chomp(readall(`hostname`))
filename = "result.$hostname.$sttime.jld"

# One-locus per chromosomal type.
chrs = LocusList((Autosome, 1), (XChromosome, 1), (YChromosome, 1), (Mitochondrion, 1))
popsize = PopulationSizes(1=>(Female=>250, Male=>250), 2=>(Female=>250, Male=>250))
mutrates = MutationRates(Female=>0.000001, Male=>0.0000001)
recrates = RecombinationRates(Female=>0.1, Male=>0.1)
propm = ProportionsOfMigrants(
    1=>(Female=>0.0001, Male=>0.0001),
    2=>(Female=>0.0001, Male=>0.0001))
fitm = FitnessOfMigrants(1=>(Female=>0.8, Male=>0.9), 2=>(Female=>0.8, Male=>0.9))
fitt = FitnessOfForeignTraits(
    1=>(Female=>1.0, Male=>0.3),
    2=>(Female=>1.0, Male=>0.3))
trans = TraitTransmissionRates(
    Mother=>(Daughter=>0.9, Son=>0.6),
    Father=>(Daughter=>0.9, Son=>0.6))

params = ModelParameters(chrs, popsize, mutrates, recrates, propm, fitm, fitt, trans)

pops, genes, t = simulate(params, 5000, 30000, 3, 100)

@save filename pops genes t


text:
@static if VERSION >= v"0.7.0-"
    using Pkg
else
    macro info(x)
        :(info($(esc(x))))
    end
end

if VERSION < v"0.7.0-"
    @info "Pkg.clone(pwd())"
    Pkg.clone(pwd())
end

@info "Pkg.build(IPython)"
if VERSION >= v"1.1.0-rc1"
    Pkg.build("IPython", verbose=true)
else
    Pkg.build("IPython")
end

using IPython

if VERSION >= v"0.7.0-"
    @info "PyCall/deps/build.log:"
    print(read(
        joinpath(dirname(dirname(pathof(IPython.PyCall))), "deps", "build.log"),
        String))
end

@info "show_versions.jl"
include("show_versions.jl")


text:
a_ = zeros(5)
y_ = 0
z_ = -9.4826529886227799e+00
r = zeros(5)

a_[1] = -9.0909090552980257e-02
a_[2] = 9.0909091401531950e-01
a_[3] = -1.0000000051491075e+01
a_[4] = 1.0000000051491075e+01
a_[5] = -1.0000000097329545e+01

r[1] = max(10 - a_[1] - a_[2] + a_[3], 0)
r[2] = max(-a_[2] + a_[2]*y_ - a_[1]*a_[4], 0)
r[3] = a_[3]*a_[2] + a_[1]*a_[4] - a_[5]
r[4] = z_ - sum(a_) + 0.001*sum(a_.^2)
r[5] = a_[3] + a_[4] - 10*y_

for i = 1:5
	println("r[",i,"] = ",r[i])
end


text:
function lin_inter(grid,vals)
    function func(x::Real)
        if x <= grid[1]
            return (vals[2] - vals[1])/(grid[2] - grid[1])*(x - grid[1]) + vals[1]
        elseif x >= grid[end]
            return (vals[end] - vals[end-1])/(grid[end] - grid[end-1])*(x - grid[end]) + vals[end]
        else
            index = searchsortedlast(grid,x)
            return (vals[index+1] - vals[index])/(grid[index+1] - grid[index])*(x - grid[index]) + vals[index]
        end
    end

    function func{T<:Real}(x::AbstractVector{T})
        m = length(x)
        out = Array(Any,m)
        for i in 1:m
            out[i] = func(x[i])
        end
        return out
    end
    return func

end


text:
import .Flux.trainable

### Base Kernels

trainable(k::ConstantKernel) = (k.c,)

trainable(k::GammaExponentialKernel) = (k.γ,)

trainable(k::GammaRationalQuadraticKernel) = (k.α, k.γ)

trainable(k::MaternKernel) = (k.ν,)

trainable(k::LinearKernel) = (k.c,)

trainable(k::PolynomialKernel) = (k.d, k.c)

trainable(k::RationalQuadraticKernel) = (k.α,)

#### Composite kernels

trainable(κ::KernelProduct) = κ.kernels

trainable(κ::KernelSum) = (κ.weights, κ.kernels) #To check

trainable(κ::ScaledKernel) = (κ.σ², κ.kernel)

trainable(κ::TransformedKernel) = (κ.transform, κ.kernel)

### Transforms

trainable(t::ARDTransform) = (t.v,)

trainable(t::ChainTransform) = t.transforms

trainable(t::FunctionTransform) = (t.f,)

trainable(t::LowRankTransform) = (t.proj,)

trainable(t::ScaleTransform) = (t.s,)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "OpenEXR"
version = v"3.1.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/AcademySoftwareFoundation/openexr/archive/refs/tags/v3.1.1.tar.gz", "045254e201c0f87d1d1a4b2b5815c4ae54845af2e6ec0ab88e979b5fdb30a86e")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/openexr*
mkdir build
cd build/
cmake \
    -DBUILD_TESTING=OFF \
    -DOPENEXR_INSTALL_TOOLS=OFF \
    -DOPENEXR_INSTALL_EXAMPLES=OFF \
    -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    ..
make -j${nproc}
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libOpenEXRUtil-3_1", :libOpenEXRUtil),
    LibraryProduct("libOpenEXRCore-3_1", :libOpenEXRCore),
    LibraryProduct("libOpenEXR-3_1", :libOpenEXR),
    LibraryProduct("libIlmThread-3_1", :libIlmThread),
    LibraryProduct("libIex-3_1", :libIex)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("Imath_jll"; compat="=3.1.2"),
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"5.2.0")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using IntegerSmithNormalForm: negateRow, negateCol, swapRows, swapCols, addRow,
addCol, SNF!, SNF, SNFWithoutTransform
using Base.Test

@testset "Negating rows and columns" begin
    A = [[1 2];[3 4]]
    B = [[1 2 3];[4 5 6];[7 8 9]]
    C = [[1 2];[3 4];[5 6]]

    (S,Tmp,T) = negateRow(A,1)
    @test Tmp == [[-1 -2];[3 4]]
    @test S*A*T == Tmp
    @inferred negateRow(A,1)

    (S,Tmp,T) = negateRow(A,2)
    @test Tmp == [[1 2];[-3 -4]]
    @test S*A*T == Tmp
    @inferred negateRow(A,2)

    (S,Tmp,T) = negateRow(B,3)
    @test Tmp == [[1 2 3];[4 5 6];[-7 -8 -9]]
    @test S*B*T == Tmp
    @inferred negateRow(B,3)

    (S,Tmp,T) = negateRow(C,2)
    @test Tmp == [[1 2];[-3 -4];[5 6]]
    @test S*C*T == Tmp
    @inferred negateRow(C,2)

    (S,Tmp,T) = negateCol(A,1)
    @test Tmp == [[-1 2];[-3 4]]
    @test S*A*T == Tmp
    @inferred negateCol(A,1)

    (S,Tmp,T) = negateCol(A,2)
    @test Tmp == [[1 -2];[3 -4]]
    @test S*A*T == Tmp
    @inferred negateCol(A,2)

    (S,Tmp,T) = negateCol(B,3)
    @test Tmp == [[1 2 -3];[4 5 -6];[7 8 -9]]
    @test S*B*T == Tmp
    @inferred negateCol(B,3)

    (S,Tmp,T) = negateCol(C,2)
    @test Tmp == [[1 -2];[3 -4];[5 -6]]
    @test S*C*T == Tmp
    @inferred negateCol(C,2)
end

@testset "Swapping rows and columns" begin
    A = [[1 2];[3 4];[5 6]]

    (S,B,T) = swapRows(A,1,2)
    @test B == [[3 4];[1 2];[5 6]]
    @test S*A*T == B
    @inferred swapRows(A,1,2)

    (S,B,T) = swapRows(A,3,2)
    @test B == [[1 2];[5 6];[3 4]]
    @test S*A*T == B
    @inferred swapRows(A,3,2)

    A = [[1 2 3];[4 5 6]]
    (S,B,T) = swapCols(A,1,2)
    @test B == [[2 1 3];[5 4 6]]
    @test S*A*T == B
    @inferred swapCols(A,1,2)

    (S,B,T) = swapCols(A,3,2)
    @test B == [[1 3 2];[4 6 5]]
    @test S*A*T == B
    @inferred swapCols(A,3,2)
end

@testset "Adding rows and columns" begin
    A = [[1 2];[3 4];[5 6]]

    (S,B,T) = addRow(A,1,2,10)
    @test B == [[31 42];[3 4];[5 6]]
    @test S*A*T == B
    @inferred addRow(A,1,2,10)

    A = [[1 2 3];[4 5 6]]

    (S,B,T) = addCol(A,2,3,10)
    @test B == [[1 32 3];[4 65 6]]
    @test S*A*T == B
    @inferred addCol(A,2,3,10)

end

function getMatrix(
                   k1,k2,k3, a
                  )

  b = 32

  if a == 0
    return [b k1 k2;
            0 b k3;
            0 0 b]
  elseif a == 1
    return [b 0 k2;
            k1 b k3;
            0 0 b]
  elseif a == 2
    return [b k1 0;
            0 b k3;
            k2 0 b]
  elseif a == 3
    return [b k1 k2;
            0 b 0;
            0 k3 b]
  elseif a == 4
    return [b 0 0;
            k1 b k3;
            k2 0 b]
  elseif a == 5
    return [b 0 k2;
            k1 b 0;
            0 k3 b]
  elseif a == 6
    return [b k1 0;
            0 b 0;
            k2 k3 b]
  elseif a == 7
    return [b 0 0;
            k1 b 0;
            k2 k3 b]
  end
end

function testSNF(A)
    B = copy(A)
    @inferred SNF!(B)
    B = copy(A)
    (S,B,T) = SNF!(B)
    @test B == S*A*T
    @test isdiag(B)
    d = diag(B)
    for i in 1:length(d)
        @test d[i] >= 0
        if d[i] != 0
            for j in i:length(d)
                @test mod(d[j],d[i]) == 0
            end
        end
    end
end

@testset "Test Smith Normal Form" begin

    testSNF([16 -16 -16; -16 16 -2; 9 -9 16])
    testSNF([[1 2 3];[4 5 6];[1 4 9]])
    testSNF([[1 2];[4 5];[1 4]])

    aRange = 0:7
    k1Range = [-24:1:24..., 99,100,101]

    simulations = length(aRange)*length(k1Range)^3

    for a in aRange
        for k1 in k1Range
            for k2 in k1Range
                for k3 in k1Range
                    M = getMatrix(k1,k2,k3,a)
                    testSNF(M)
                end
            end
        end
    end

    (S,B,T) = SNF([[1 2 3];[4 5 6];[1 4 9]])
    @inferred SNF([[1 2 3];[4 5 6];[1 4 9]])
    @test B == S*[[1 2 3];[4 5 6];[1 4 9]]*T
    @test isdiag(B)
    d = diag(B)
    for i in 1:length(d)
        @test d[i] >= 0
        for j in i:length(d)
            @test mod(d[j],d[i]) == 0
        end
    end

    B = SNFWithoutTransform([[1 2 3];[4 5 6];[1 4 9]])
    @inferred SNFWithoutTransform([[1 2 3];[4 5 6];[1 4 9]])
    @test isdiag(B)
    d = diag(B)
    for i in 1:length(d)
        @test d[i] >= 0
        for j in i:length(d)
            @test mod(d[j],d[i]) == 0
        end
    end
end


text:
__precompile__()

module PoincareInvariants

    using ApproxFun
    using HDF5
    using ProgressMeter
    using GeometricIntegrators
    using GeometricIntegrators.Utils

    export PoincareInvariant1st, PoincareInvariant1stCanonical,
           PoincareInvariant2nd, PoincareInvariant2ndCanonical,
           PoincareInvariant2ndTrapezoidal,
           evaluate_poincare_invariant, write_to_hdf5

    include("matrix_utils.jl")
    include("poincare_invariant_1st_common.jl")
    include("poincare_invariant_1st.jl")
    include("poincare_invariant_1st_canonical.jl")
    include("poincare_invariant_2nd.jl")
    include("poincare_invariant_2nd_canonical.jl")
    include("poincare_invariant_2nd_trapezoidal.jl")
    include("poincare_invariant_2nd_common.jl")

end


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

using Markdown
using InteractiveUtils

# ╔═╡ e8433caa-1882-4153-bb6b-881b64eda8c5
using StatsPlots,DataFrames,CSV,KernelDensity

# ╔═╡ a524bb50-d918-11eb-18ba-ff416e617daf
md"# Kernel Tests"

# ╔═╡ 2be06875-e8dc-4ac6-bb21-4b094c453fbb
lcf=CSV.File( "/home/graham_s/OU/DD309/2020J/econometrics/data/lcf/lcf_017_18.tab",delim='\t')|>DataFrame

# ╔═╡ 8f3b8c66-9089-467b-a879-f1a8717036aa
atk = KernelDensity.kde((lcf.total_expend, lcf.alcohol_tobacco)) # note this has to be a tuple, for some reason

# ╔═╡ e37800e2-0db4-4e09-8f91-8680e62407fc
StatsPlots.plot(atk)

# ╔═╡ 69ec345f-b97a-40af-97b5-a37600b1d7b2
p1 = @df lcf marginalscatter(:total_expend,:food_and_drink,xlabel="expenditure", ylabel="food and drink")

# ╔═╡ 732de09c-a4e3-47f2-aa81-1152e05b8bec
marginalkde(lcf.total_expend,lcf.food_and_drink,xlabel="expenditure", ylabel="food and drink")

# ╔═╡ Cell order:
# ╠═a524bb50-d918-11eb-18ba-ff416e617daf
# ╠═e8433caa-1882-4153-bb6b-881b64eda8c5
# ╠═2be06875-e8dc-4ac6-bb21-4b094c453fbb
# ╠═8f3b8c66-9089-467b-a879-f1a8717036aa
# ╠═e37800e2-0db4-4e09-8f91-8680e62407fc
# ╠═69ec345f-b97a-40af-97b5-a37600b1d7b2
# ╠═732de09c-a4e3-47f2-aa81-1152e05b8bec

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
## Utility functions for evaluating the performance of classification methods
module PerformanceEval

export auc, rank_enrichment, rank_precision, evaluate_scores, get_enrichment, notHistone, precision_recall

# fill in any blank spots with nearby values
function fill_nans(vals)
    newVals = copy(vals)
    for i in 2:length(vals)
        if isnan(vals[i])
            newVals[i] = newVals[i-1]
        end
    end
    for i in length(vals)-1:-1:1
        if isnan(vals[i])
            newVals[i] = newVals[i+1]
        end
    end
    newVals
end

# compute the rank enrichment curve
function rank_enrichment(truth, pred)
    rank_enrichment(truth[sortperm(pred, rev=true)])
end
function rank_enrichment(rankedTruth, weights=nothing)
    weights = weights == nothing ? ones(length(rankedTruth)) : weights
    x = zeros(Float64, length(rankedTruth))
    y = zeros(Float64, length(rankedTruth))
    randRate = dot(rankedTruth,weights) / sum(weights)
    totalMatchedWeight = 0
    totalSeenWeight = 0
    for i in 1:length(rankedTruth)
        totalMatchedWeight += rankedTruth[i] * weights[i]
        totalSeenWeight += weights[i]
        x[i] = totalSeenWeight
        y[i] = (totalMatchedWeight/totalSeenWeight) / randRate
    end
    x,fill_nans(y)
end
function rank_precision(rankedTruth, weights=nothing)
    weights = weights == nothing ? ones(length(rankedTruth)) : weights
    x,y = rank_enrichment(rankedTruth, weights)
    randRate = dot(rankedTruth,weights) / sum(weights)
    x,y .* randRate
end

function auc(x, y)
    total = 0.0
    lastx = 0.0
    for i in 1:length(x)
        if x[i] != lastx
            total += y[i] * (x[i] - lastx)
            lastx = x[i]
        end
    end
    total
end

function precision_recall(truth, pred)
    precision_recall(truth[sortperm(pred, rev=true)])
end
function precision_recall(rankedTruth)
    x = Int64[]
    y = Float64[]
    total = 0
    for i in 1:length(rankedTruth)
        total += rankedTruth[i]
        push!(x, total)
        push!(y, total/i)
    end

    # prune out the points that fall between other points vertically
    newx = Float64[]
    newy = Float64[]
    push!(newx, x[1])
    push!(newy, y[1])
    for i in 2:length(x)-1
        if x[i-1] == x[i] == x[i+1] && y[i-1] > y[i] > y[i+1]
            # ignore since this point adds no new information
        else
            push!(newx, x[i])
            push!(newy, y[i])
        end
    end
    push!(newx, x[end])
    push!(newy, y[end]);

    newx,newy
end

function get_enrichment(data; quiet=false)
    totalRight = sum(data)
    topFound = sum(data[1:totalRight])
    if !quiet println(STDERR, "totalRight $totalRight, total found $topFound") end
    return (topFound/totalRight)/(totalRight/length(data))
end

end

max_stars_repo_name: slundberg/PerformanceEval

max_issues_repo_name: slundberg/PerformanceEval

max_forks_repo_name: slundberg/PerformanceEval


text:
# There is a misplaced leading comma on line 50, we should get rid of it

dataset_path = joinpath(dirname(@__FILE__), "echocardiogram.data.orig")

f = open(dataset_path)
data = readlines(f)
close(f)

f = open("$dataset_path.custom", "w")
for (i, line) in enumerate(data)
  if i == 50
    line = line[2:end]
  end
  write(f, "$line")
end
close(f)

max_stars_repo_name: CoventryResearch/uci-data

max_issues_repo_name: CoventryResearch/uci-data

max_forks_repo_name: CoventryResearch/uci-data


text:
##############
# GitHubType #
##############
# A `GitHubType` is a Julia type representation of a JSON object defined by the
# GitHub API. Generally:
#
# - The fields of these types should correspond to keys in the JSON object. In
#   the event the JSON object has a "type" key, the corresponding field name
#   used should be `typ` (since `type` is a reserved word in Julia).
#
# - The method `name` should be defined on every GitHubType. This method
#   returns the type's identity in the form used for URI construction. For
#   example, `name` called on an `Owner` will return the owner's login, while
#   `name` called on a `Commit` will return the commit's sha.
#
# - A GitHubType's field types should be Nullables of either concrete types, a
#   Vectors of concrete types, or Dicts.

abstract GitHubType

typealias GitHubString UTF8String

function Base.(:(==))(a::GitHubType, b::GitHubType)
    if typeof(a) != typeof(b)
        return false
    end

    for field in fieldnames(a)
        aval, bval = getfield(a, field), getfield(b, field)
        if isnull(aval) == isnull(bval)
            if !(isnull(aval)) && get(aval) != get(bval)
                return false
            end
        else
            return false
        end
    end

    return true
end

# `namefield` is overloaded by various GitHubTypes to allow for more generic
# input to AP functions that require a name to construct URI paths via `name`
name(val) = val
name(g::GitHubType) = get(namefield(g))

########################################
# Converting JSON Dicts to GitHubTypes #
########################################

function extract_nullable{T}(data::Dict, key, ::Type{T})
    if haskey(data, key)
        val = data[key]
        if !(isa(val, Void))
            if T <: Vector
                V = eltype(T)
                return Nullable{T}(V[prune_github_value(v, V) for v in val])
            else
                return Nullable{T}(prune_github_value(val, T))
            end
        end
    end
    return Nullable{T}()
end

prune_github_value{T}(val, ::Type{T}) = T(val)
prune_github_value(val, ::Type{Dates.DateTime}) = Dates.DateTime(chopz(val))

# ISO 8601 allows for a trailing 'Z' to indicate that the given time is UTC.
# Julia's Dates.DateTime constructor doesn't support this, but GitHub's time
# strings can contain it. This method ensures that a string's trailing 'Z',
# if present, has been removed.
function chopz(str::AbstractString)
    if !(isempty(str)) && last(str) == 'Z'
        return chop(str)
    end
    return str
end

# Calling `json2github(::Type{G<:GitHubType}, data::Dict)` will parse the given
# dictionary into the type `G` with the expectation that the fieldnames of
# `G` are keys of `data`, and the corresponding values can be converted to the
# given field types.
@generated function json2github{G<:GitHubType}(::Type{G}, data::Dict)
    types = G.types
    fields = fieldnames(G)
    args = Vector{Expr}(length(fields))
    for i in eachindex(fields)
        field, T = fields[i], first(types[i].parameters)
        key = field == :typ ? "type" : string(field)
        args[i] = :(extract_nullable(data, $key, $T))
    end
    return :(G($(args...))::G)
end

#############################################
# Converting GitHubType Dicts to JSON Dicts #
#############################################

github2json(val) = val
github2json(uri::HttpCommon.URI) = string(uri)
github2json(dt::Dates.DateTime) = string(dt) * "Z"
github2json(v::Vector) = [github2json(i) for i in v]

function github2json(g::GitHubType)
    results = Dict()
    for field in fieldnames(g)
        val = getfield(g, field)
        if !(isnull(val))
            key = field == :typ ? "type" : string(field)
            results[key] = github2json(get(val))
        end
    end
    return results
end

function github2json{K}(data::Dict{K})
    results = Dict{K,Any}()
    for (key, val) in data
        results[key] = github2json(val)
    end
    return results
end

###################
# Pretty Printing #
###################

function Base.show(io::IO, g::GitHubType)
    print(io, "$(typeof(g)) (all fields are Nullable):")
    for field in fieldnames(g)
        val = getfield(g, field)
        if !(isnull(val))
            gotval = get(val)
            println(io)
            print(io, "  $field: ")
            if isa(gotval, Vector)
                print(io, typeof(gotval))
            else
                showcompact(io, gotval)
            end
        end
    end
end

function Base.showcompact(io::IO, g::GitHubType)
    uri_id = namefield(g)
    if isnull(uri_id)
        print(io, typeof(g), "(…)")
    else
        print(io, typeof(g), "($(repr(get(uri_id))))")
    end
end


text:
# Weave readme
using Pkg
Pkg.activate("./readme-env")

using Weave

weave("README.jmd", out_path = :pwd, doctype = "github")

if false
    tangle("README.jmd")
end


text:
"""
LTESOM - Long Term Energy System Optimization Model
Copyright (C) 2022, Peking University, College of Engineering (COE),
Department of Industrial Engineering and Management
"""

println("Activating the Julia virtual environment")

using Pkg
Pkg.activate("LTESOMEnv")

## Possible solvers for GenX
# Open source solvers
Pkg.add(Pkg.PackageSpec(name="Cbc", version="0.8.0"))
Pkg.add(Pkg.PackageSpec(name="Clp", version="0.8.4"))
Pkg.add(Pkg.PackageSpec(name="GLPK", version="0.14.12"))
Pkg.add(Pkg.PackageSpec(name="Ipopt", version="0.7.0"))
# Pkg.add(Pkg.PackageSpec(name="SCIP"))
# Pkg.build("SCIP")

## Commercial solvers - depending on the host machine
# Pkg.add(Pkg.PackageSpec(name="CPLEX"))
# Pkg.build("CPLEX")
Pkg.add(Pkg.PackageSpec(name="Gurobi"))
Pkg.build("Gurobi")

## Julia interface wrapper
Pkg.add(Pkg.PackageSpec(name="JuMP", version="0.21.8"))
Pkg.add("MathOptInterface")

## Data structures
Pkg.add(Pkg.PackageSpec(name="Dates"))
Pkg.add(Pkg.PackageSpec(name="YAML", version="0.4.7"))
Pkg.add(Pkg.PackageSpec(name="CSV", version="0.8.5"))
Pkg.add(Pkg.PackageSpec(name="Distances", version="0.10.3"))
Pkg.add(Pkg.PackageSpec(name="DataFrames", version="1.0.0"))
Pkg.add(Pkg.PackageSpec(name="DataStructures", version="0.18.9"))
Pkg.add(Pkg.PackageSpec(name="Clustering", version="0.14.2"))
Pkg.add(Pkg.PackageSpec(name="Combinatorics", version="1.0.2"))

## Math base packages
Pkg.add(Pkg.PackageSpec(name="Random"))
Pkg.add(Pkg.PackageSpec(name="MathProgBase"))
Pkg.add(Pkg.PackageSpec(name="StatsBase"))
Pkg.add(Pkg.PackageSpec(name="LinearAlgebra"))
Pkg.add(Pkg.PackageSpec(name="DiffEqSensitivity"))
Pkg.add(Pkg.PackageSpec(name="OrdinaryDiffEq"))
Pkg.add(Pkg.PackageSpec(name="QuasiMonteCarlo"))
Pkg.add(Pkg.PackageSpec(name="Statistics"))
Pkg.add(Pkg.PackageSpec(name="RecursiveArrayTools"))

## Documention tools
Pkg.add(Pkg.PackageSpec(name="Documenter", version="0.27.3"))
Pkg.add(Pkg.PackageSpec(name="DocumenterTools", version="0.1.13"))

## Benchmark tools
Pkg.add(Pkg.PackageSpec(name="BenchmarkTools"))

## Show current environment
Pkg.status()

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
module BijectivePaper

# Columns are odd resolution, even resolution (rmin)
# Columns are quadrants in x and y like you'd graph it.
# 1 2
# 0 3
hilbert_variant_encode_table = Function[
    ((x, y, w) -> (x, y))          ((x, y, w) -> (x, y))
    ((x, y, w) -> (y-w, x))        ((x, y, w) -> (y, x-w))
    ((x, y, w) -> (y-w, x-w))      ((x, y, w) -> (y-w, x-w))
    ((x, y, w) -> (2w-x-1, w-y-1)) ((x, y, w) -> (w-x-1, 2w-y-1))
]

function encode_hilbert_zero(x::Integer, y::Integer)
    z = zero(x)
    if x == 0 && y == 0
        return z
    end
    rmin = convert(typeof(x), floor(log2(max(x, y))) + 1)
    w = 1<<(rmin - 1)
    while rmin > 0
        if rmin&1 == 1  # odd
            quadrant = x < w ? (y < w ? 0 : 1) : (y < w ? 3 : 2)
            parity = 1
        else  # even
            quadrant = x < w ? (y < w ? 0 : 3) : (y < w ? 1 : 2)
            parity = 2
        end
        z = (z<<2) + quadrant
        x, y = hilbert_variant_encode_table[quadrant+1, parity](x, y, w)
        rmin -= 1
        w >>= 1
    end
    z
end

hilbert_variant_decode_table = Function[
    ((x, y, w) -> (x, y))          ((x, y, w) -> (x, y))
    ((x, y, w) -> (y, x+w))        ((x, y, w) -> (y+w, x))
    ((x, y, w) -> (y+w, x+w))      ((x, y, w) -> (y+w, x+w))
    ((x, y, w) -> (2w-x-1, w-y-1)) ((x, y, w) -> (w-x-1, 2w-y-1))
]

function decode_hilbert_zero(z)
    r = z & 3
    x, y = typeof(z).([(0, 0), (0, 1), (1, 1), (1, 0)][r + 1])
    z >>= 2
    rmin = 2
    w = convert(typeof(z), 2)
    while z > 0
        r = z & 3
        parity = 2 - rmin&1
        x, y = hilbert_variant_decode_table[r+1, parity](x, y, w)
        z >>= 2
        rmin += 1
        w <<= 1
    end
    x, y
end

encode_hilbert(x::Integer, y::Integer) = encode_hilbert_zero(x - 1, y - 1) + 1

function decode_hilbert(z::Integer)
    x, y = decode_hilbert_zero(z - 1)
    x + 1, y + 1
end

"""
    hilbert_order(v::AbstractArray, subdivisions)

Calculates the permutation of `v` that would order a 2D array by Hilbert curve,
where `v` is an array of real numbers, dimension (2, N) and subdivisions is
a real number that specifies how many boxes to place the `v` into, per side.

This tells you how to order an array of 2-dimensional numbers so that
they have more memory locality in 1 dimension.
Given an array of real numbers of dimension `(2, n)`, subdivide them in each
dimension by `subdivisions`, and assign each point a Hilbert code. Return
the permutation that would sort the given array by that Hilbert code.

See also: [`encode_hilbert_zero`](@ref).

# Example
```julia
rng = MersenneTwister(984720987)
points_in_space = zeros(2, 100)
rand!(rng, points_in_space)
points_reordered = points_in_space[:, hilbert_order(points_in_space, 50)]
```
"""
function hilbert_order(v::AbstractArray, subdivisions)
    lowx = minimum(v[1, :])
    lowy = minimum(v[2, :])
    highx = maximum(v[1, :])
    highy = maximum(v[2, :])
    iv = zeros(Int64, size(v, 2))
    for i in 1:size(v, 2)
        iv[i] = encode_hilbert_zero(
                    round(Int64, (v[1, i] - lowx) * subdivisions / (highx - lowx)),
                    round(Int64, (v[2, i] - lowy) * subdivisions / (highy - lowy))
                    )
    end
    sortperm(iv)
end

export encode_hilbert_zero, decode_hilbert_zero
export encode_hilbert, decode_hilbert, hilbert_order

end


text:
include("demo_pulsed_solution.jl")
include("demo_steady_solution.jl")

include("arc_solution.jl")

include("act_1_solution.jl")
include("act_2_solution.jl")


text:
module EarlyStopping

using Dates
using Statistics
import Base.+

export StoppingCriterion,
    Never,
    InvalidValue,
    NotANumber, # deprecated
    TimeLimit,
    GL,
    NumberSinceBest,
    Patience,
    UP,
    PQ,
    NumberLimit,
    Threshold,
    Disjunction,
    Warmup,
    criteria,
    stopping_time,
    EarlyStopper,
    done!,
    message,
    needs_training_losses,
    needs_loss

include("api.jl")
include("criteria.jl")
include("disjunction.jl")
include("stopping_time.jl")
include("object_oriented_api.jl")

end # module


text:
using Test

@test 1==1