text:
"""
    solve(prob::MultiTermsFODEProblem, h, PIIMTrap())

Use implicit product integration trapezoidal type method to solve multi-terms FODE.
"""
struct PIIMTrap <: FractionalDiffEqAlgorithm end

function solve(prob::MultiTermsFODEProblem, h, ::PIIMTrap; abstol=1e-6)
    @unpack parameters, orders, rightfun, u0, tspan = prob
    t0 = tspan[1]; T = tspan[2]
    u0 = u0[:]'

    # Generate the jacobian of the given function
    J_fun(x, y) = ForwardDiff.derivative(x -> rightfun(x, y), x)

    Q = length(orders)
    orders= sort(orders)
    i_al = sortperm(orders, rev=true)
    parameters = parameters[i_al]
    al_Q = orders[end]
    al_i = orders[1:end-1]
    lam_Q = parameters[end]
    lam_rat_i = parameters[1:end-1]/lam_Q
    m_Q = ceil(Int64, al_Q)
    m_i = ceil.(Int64, orders[1:end-1])
    bet = [al_Q .- al_i; al_Q]

    itmax = 100

    gamma_val = zeros(Q, m_Q)
    for i = 1 : Q-1
        k = collect(Int, 0:m_i[i]-1)
        gamma_val[i, k.+1] = gamma.(k.+bet[i].+1)
    end
    k = collect(0:m_Q-1)
    gamma_val[Q, :] = factorial.(k)

    problem_size = size(u0, 1)

    r::Int = 64
    N::Int = ceil(Int, (T-t0)/h)
    Nr::Int = ceil(Int, (N+1)/r)*r
    Qr::Int = ceil(Int, log2((Nr)/r)) - 1
    NNr::Int = 2^(Qr+1)*r

    y = zeros(problem_size, N+1)
    fy = zeros(problem_size, N+1)
    zn = zeros(problem_size, NNr+1, Q)

    nvett = collect(0:NNr+1)
    an = zeros(Q, NNr+1); a0 = zeros(Q, NNr+1)
    for i = 1 : Q
        nbeta = nvett.^bet[i]; nbeta1 = nbeta.*nvett ;
        an[i, :] = [1;  (nbeta1[1:end-2] - 2*nbeta1[2:end-1] + nbeta1[3:end])]*h^bet[i]/gamma(bet[i]+2)
        a0[i, :] = [0;  nbeta1[1:end-2] - nbeta[2:end-1].*(nvett[2:end-1].-bet[i].-1)]*h^bet[i]/gamma(bet[i]+2)
    end
    C = 0
    for i = 1 : Q-1
        C = C + lam_rat_i[i]*an[i, 1]
    end

    t = collect(0:N)*h
    y[:, 1] = u0[:, 1]
    fy[:, 1] .= rightfun(t0, u0[:, 1])
    (y, fy) = PIExTriangolo(1, r-1, t, y, fy, zn, N, a0, an, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val, rightfun, lam_Q, C, abstol, itmax, J_fun)

    ff = zeros(1, 2^(Qr+2)); ff[1:2] = [0 2]; card_ff = 2
    nx0 = 0; nu0 = 0
    for qr = 0 : Qr
        L = 2^qr
        (y, fy) = PIIMTrapDisegnaBlocchi(L, ff, r, Nr, nx0+L*r, nu0, t, y, fy, N, zn, a0, an, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val, rightfun, lam_Q, C, J_fun, itmax, abstol)
        ff[1:2*card_ff] = [ff[1:card_ff] ff[1:card_ff]]
        card_ff = 2*card_ff
        ff[card_ff] = 4*L
    end

    if T<t[N+1]
        c = (T - t[N])/h
        t[N+1] = tfinal
        y[:, N+1] = (1-c)*y[:, N] + c*y[:, N+1]
    end

    t = t[1:N+1] ; y = y[:, 1:N+1]
    return FODESolution(t, y[:])
end

function PIIMTrapDisegnaBlocchi(L, ff, r, Nr, nx0, nu0, t, y, fy, N, zn, a0, an, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val, rightfun, lam_Q, C, J_fun, itmax, abstol)

    nxi::Int = nx0
    nxf::Int = nx0 + L*r - 1
    nyi::Int = nu0
    nyf::Int = nu0 + L*r - 1
    is = 1
    s_nxi = zeros(N)
    s_nxf = zeros(N)
    s_nyi = zeros(N)
    s_nyf = zeros(N)
    s_nxi[1] = nxi
    s_nxf[1] = nxf
    s_nyi[1] = nyi
    s_nyf[1] = nyf

    i_triangolo = 0; stop = false
    while stop == false
        stop = (nxi+r-1 == nx0+L*r-1) || (nxi+r-1 >= Nr-1)

        zn = PIIMTrapQuadrato(nxi, nxf, nyi, nyf, y, fy, zn, an, problem_size, Q)

        (y, fy) = PIExTriangolo(nxi, nxi+r-1, t, y, fy, zn, N, a0, an, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val, rightfun, lam_Q, C, abstol, itmax, J_fun)
        i_triangolo = i_triangolo + 1

        if stop==false
            if nxi+r-1 == nxf
                i_Delta = ff[i_triangolo]
                Delta = i_Delta*r
                nxi = s_nxf[is]+1; nxf = s_nxf[is]  + Delta
                nyi = s_nxf[is] - Delta +1; nyf = s_nxf[is]
                s_nxi[is] = nxi
                s_nxf[is] = nxf
                s_nyi[is] = nyi
                s_nyf[is] = nyf
            else
                nxi = nxi + r; nxf = nxi + r - 1; nyi = nyf + 1; nyf = nyf + r
                is = is + 1
                s_nxi[is] = nxi
                s_nxf[is] = nxf
                s_nyi[is] = nyi
                s_nyf[is] = nyf
            end
        end

    end
    return y, fy
end

function PIIMTrapQuadrato(nxi, nxf, nyi, nyf, y, fy, zn, an, problem_size, Q)
    coef_beg = nxi-nyf; coef_end = nxf-nyi+1
    funz_beg = nyi+1; funz_end = nyf+1

    for i = 1:Q
        vett_coef = an[i, coef_beg:coef_end]
        if nyi == 0
            if i < Q
                vett_funz = [zeros(problem_size, 1)  y[:, funz_beg+1:funz_end] zeros(problem_size, funz_end-funz_beg+1)]
            else
                vett_funz = [zeros(problem_size, 1)  fy[:, funz_beg+1:funz_end] zeros(problem_size, funz_end-funz_beg+1)]
            end
        else
            if i < Q
                vett_funz = [y[:, funz_beg:funz_end] zeros(problem_size, funz_end-funz_beg+1)]
            else
                vett_funz = [fy[:, funz_beg:funz_end] zeros(problem_size, funz_end-funz_beg+1)]
            end
        end
        zzn = real.(FastConv(vett_coef, vett_funz))
        zn[:, nxi+1:nxf+1, i] = zn[:, nxi+1:nxf+1, i] + zzn[:, nxf-nyf+1:end]
    end
    return zn
end

function PIExTriangolo(nxi, nxf, t, y, fy, zn, N, a0, an, t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val, rightfun, lam_Q, C, abstol, itmax, J_fun)

    for n = nxi:min(N, nxf)
        St = PIIMTrapStartingTerm_Multi(t[n+1], t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val)
        Phi_n = copy(St)

        for i = 1:Q-1
            temp = a0[i, n+1]*y[:, 1]+zn[:, n+1, i]
            for j = nxi:n-1
                temp = temp + an[i, n-j+1]*y[:, j+1]
            end
            Phi_n = Phi_n - lam_rat_i[i]*temp
        end
        temp = a0[Q, n+1]*fy[:, 1] + zn[:, n+1, Q]
        for j = nxi : n-1
            temp = temp + an[Q, n-j+1]*fy[:, j+1]
        end
        Phi_n = Phi_n + temp/lam_Q

        yn0 = y[:, n]; fn0 = rightfun(t[n+1], yn0)
        Jfn0 = Jf_vectorfield(t[n+1], yn0, J_fun)
        Gn0 = (1+C)*yn0 .- an[Q, 1]./lam_Q*fn0 .- Phi_n
        stop = false; it = 0

        while ~stop

            JGn0 = (1+C)*(zeros(problem_size, problem_size)+I) .- an[Q, 1]/lam_Q*Jfn0
            global yn1 = yn0 - JGn0\Gn0
            global fn1 = rightfun(t[n+1], yn1)
            Gn1 = (1+C)*yn1 .- an[Q, 1]/lam_Q*fn1 .- Phi_n
            it = it + 1

            stop = (norm(yn1-yn0, Inf) < abstol) || (norm(Gn1, Inf) < abstol)
            if it > itmax && ~stop
                @warn "Non Convergence"
                stop = true
            end

            yn0 = yn1 ; Gn0 = Gn1
            if ~stop
                Jfn0 = Jf_vectorfield(t[n+1], yn0, J_fun)
            end

        end

        y[:, n+1] = yn1
        fy[:, n+1] .= fn1

    end
    return y, fy
end

Jf_vectorfield(t, y, fun)=fun(t, y)

function PIIMTrapStartingTerm_Multi(t,t0, problem_size, u0, Q, m_Q, m_i, bet, lam_rat_i, gamma_val)
    ys = zeros(problem_size)

    for k = 0:m_Q-1
        ys = ys .+ (t-t0)^k./gamma_val[Q, k+1]*u0[:, k+1]
    end
    for i = 1 : Q-1
        temp = zeros(problem_size)
        for k = 0 : m_i[i]-1
            temp = temp .+ (t-t0)^(k+bet[i])/gamma_val[i, k+1]*u0[:, k+1]
        end
        ys = ys + lam_rat_i[i]*temp
    end
    return ys'
end


text:
"""
    module Events

Provides the abstract [`Event`](#) type and concrete event types.

- [`Init`](#) is called once for every callback. Use for initialization code
  that needs access to some `Learner` state.

Events in [`TrainingPhase`](#) and [`ValidationPhase`](#):

- [`EpochBegin`](#) and [`EpochEnd`](#), called at the beginning and end of each
  epoch.
- [`BatchBegin`](#) and [`BatchEnd`](#), called at the beginning and end of each
  batch.
- [`LossBegin`](#), called after the forward pass but before the loss calculation.

[`TrainingPhase`](#) only:

- [`BackwardBegin`](#), called after forward pass and loss calculation but before gradient
  calculation.
- [`BackwardEnd`](#), called after gradient calculation but before parameter update.

"""
module Events

"""
    abstract type Event

Abstract type for events that callbacks can hook into
"""
abstract type Event end

"""
    Init <: Event

Called once when the learner is created/the callback is added.

Hook into this for callback initialization that depends on the
[`Learner`](#)'s state.
"""
struct Init <: Event end

"""
    EpochBegin()

[`Event`](#) called at the beginning of an epoch.
"""
struct EpochBegin <: Event end

"""
    EpochEnd()

[`Event`](#) called at the end of an epoch.
"""
struct EpochEnd <: Event end

"""
    BatchBegin()

[`Event`](#) called at the beginning of a batch.
"""
struct BatchBegin <: Event end

"""
    BatchEnd()

[`Event`](#) called at the end of a batch.
"""
struct BatchEnd <: Event end

"""
    LossBegin()

[`Event`](#) called between calculating `y_pred` and calculating loss
"""
struct LossBegin <: Event end

"""
    BackwardBegin()

[`Event`](#) called between calculating loss and calculating gradients
"""
struct BackwardBegin <: Event end

"""
    BackwardEnd()

[`Event`](#) called between calculating gradients and updating parameters.
"""
struct BackwardEnd <: Event end

export
    # abstract
    Event,
    # concrete
    Init,
    EpochBegin, EpochEnd,
    BatchBegin, BatchEnd,
    LossBegin,
    BackwardBegin, BackwardEnd

end # module

using .Events

ext: jl


text:
@testset "ArchTropAmoeba" begin
    @polyvar x y
    f = x^2+y^2+1

    grid = Grid2D((-2, 2), 30)
    B = amoeba(f, grid=grid, alg=ArchTrop())
    @test any(B.data)
end


text:
## Function for permutation tests
function permtest(Y, res)
    # Permutation of matrix Y
    Y_perm = hcat(shuffle.(eachcol(Y))...)
    # Recalculate BD statistics
    res_p = betadiv(Y_perm)
    # Test if permuted LCBD is greater than original LCBD
    ge = res_p.LCBDi .>= res.LCBDi
    return ge
end

## Function combining BD calculation & permutation tests
function betadiv_perm(Y; nperm=999, distributed=true)
    n = size(Y, 1)
    p = size(Y, 2)
    # Initial BD results
    res = betadiv(Y)
    # Permutations
    if nperm > 0
        nGE_L = ones(Int64, n)
        # Permutation test, clumsy parallelization
        ge = @showprogress pmap(x -> permtest(Y, res), 1:nperm; distributed=distributed)
        # Compile number of permuted LCBDs greater than original LCBD
        for bitarray in ge
            nGE_L[findall(bitarray)] .+= 1
        end
        # Standardize counts
        p_LCBD = nGE_L ./ (nperm + 1)
    else
        p_LCBD = nothing
    end
    # Combine results in tuple
    res_perm = (res..., pLCBD=p_LCBD)
    return res_perm
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
module FringeExplorer

using RadioInterferometry
using Symbolics

export ha2t
export t2ha
export td2wdw
export td2wdw!
export hd2wdw
export hd2wdw!

export telinfo
export antpos
include("telinfo.jl")

export s2c
export c2s
export beamrings
include("spherical.jl")

export MJD0
export radec2hadec
include("erfa_helpers.jl")

"""
Rotaional rate of Earth in radians per second
"""
const Ωe = 7.29211538e-5

"""
    ha2t(ha::Real, Ω=Ωe) -> t

Convert hour angle `ha` in radians to time `t` in seconds. `ha == t == 0` at
transit.  `Ω`, the rate of rotation in radians per second, defaults to `Ωe`, the
rotational rate of Earth.
"""
ha2t(ha::Real, Ω=Ωe) = ha / Ω

"""
    t2ha(t::Real, Ω=Ωe) -> ha

Convert time `t` in seconds to hour angle `ha` in radians. `ha == t == 0` at
transit.  `Ω`, the rate of rotation in radians per second, defaults to `Ωe`, the
rotational rate of Earth.
"""
t2ha(t::Real, Ω=Ωe) = t * Ω

# Symbolic variables
const Ω, t, δ = @variables Ω t δ
const Dt = Differential(t)
const h = Ω * t

# Passing zero for longitude requires rotating XYZ position vectors around Z
# such that positive X intersects the local meridian.
const xyzR3uvw = xyz2uvw(h, δ, 0)

# Substitute value to `Ω`
const xyzR3uvw_e = substitute.(xyzR3uvw, [Ω=>Ωe])

# Delay is w component (use 3:3 to keep as a one row Matrix)
const w = xyzR3uvw_e[3:3,:]
# Delay rate is derivative of w component
const dw = expand_derivatives.(Dt.(w))
# Make 2x3 matrix from w and dw
const wdw_m = [w; dw]

# Create expressions for functions that return wdw for given t and δ
const wdw_expr = build_function(wdw_m, t, δ)

"""
    td2wdw(t, δ) -> Matrix
    tδ2wdw(t, δ) -> Matrix

Returns a 2x3 Matrix that can left multiply a 3 element vector (or a 3 row
matrix) of antenna position(s) to get delay `w` in the first element (or row)
and the delay rate `dw` in the second for a source observed at time `t` and
declination `δ`.  The time `t` should be in seconds, with `t<0` before the
source transits, `t==0` when the source transits, and `t>0` after the source
transits.

# Example
```julia
wdw = td2wdw(0, 0) * antpos
```
"""
td2wdw = eval(wdw_expr[1])
const tδ2wdw = td2wdw

"""
    td2wdw!(dst, t, δ)
    tδ2wdw!(dst, t, δ)

Stores a 2x3 Matrix in `dst` (which should already be a 2x3 Matrix) that can
left multiply a 3 element vector (or a 3 row matrix) of antenna position(s) to
get delay `w` in the first element (or row) and the delay rate `dw` in the
second for a source observed at time `t` and declination `δ`.  The time `t`
should be in seconds, with `t<0` before the source transits, `t==0` when the
source transits, and `t>0` after the source transits.
"""
td2wdw! = eval(wdw_expr[2])
const tδ2wdw! = td2wdw!

"""
    hd2wdw(h, δ, Ω=Ωe) -> Matrix
    hδ2wdw(h, δ, Ω=Ωe) -> Matrix

Returns a 2x3 Matrix that can left multiply a 3 element vector (or a 3 row
matrix) of antenna position(s) to get delay `w` in the first element (or row)
and the delay rate `dw` in the second for a source observed at hour angle `h` and
declination `δ`.  The hour angle `h` should be in radians, with `h<0` before the
source transits, `h==0` when the source transits, and `h>0` after the source
transits.  `Ω` is used to convert `h` to time (see [`ha2t`](@ref)).

# Example
```julia
wdw = hd2wdw(ha, dec) * antpos
```
"""
hd2wdw(h, δ, Ω=Ωe) = td2wdw(ha2t(h, Ω), δ)
hδ2wdw = hd2wdw

"""
    hd2wdw!(dst, t, δ, Ω=Ωe)
    hδ2wdw!(dst, t, δ, Ω=Ωe)

Stores a 2x3 Matrix in `dst` (which should already be a 2x3 Matrix) that can
left multiply a 3 element vector (or a 3 row matrix) of antenna position(s) to
get delay `w` in the first element (or row) and the delay rate `dw` in the
second for a source observed at time `t` and declination `δ`.  The time `t`
should be in seconds, with `t<0` before the source transits, `t==0` when the
source transits, and `t>0` after the source transits.  `Ω` is used to convert
`h` to time (see [`ha2t`](@ref)).
"""
hd2wdw!(dst, h, δ, Ω=Ωe) = td2wdw!(dst, ha2t(h, Ω), δ)
hδ2wdw! = hd2wdw!

end # module FringeExplorer


text:
using Pkg
#Pkg.activate(ENV["HOME"] * "/.julia/environments/" * "covid19")
#Pkg.add(PackageSpec(url="https://github.com/kouketsu/GSW.jl", rev="master"))
Pkg.add("GibbsSeaWater")
Pkg.add(PackageSpec(url="https://github.com/JuliaGeo/GeoDatasets.jl", rev="master"))
Pkg.add("IJulia"); Pkg.build("IJulia"); using IJulia
Pkg.add("Plots")
Pkg.add("PlotlyBase")
Pkg.add("Plotly")
Pkg.add("PlotlyJS")
Pkg.add("ColorSchemes")
Pkg.add("Formatting")
Pkg.add("WebIO")
Pkg.add("PyPlot")
Pkg.add("PyCall")
Pkg.add("Conda")
Pkg.add("CSV")
Pkg.add("XLSX")
Pkg.add("Glob")
Pkg.add("HTTP")
Pkg.add("DataFrames")
Pkg.add("Missings")
Pkg.add("FileIO")
Pkg.add("JLD2")
Pkg.add("Dates")
Pkg.add("Pandas")
Pkg.add("NCDatasets")
Pkg.add("Statistics")
Pkg.add("NaNMath")
Pkg.add("JSON3")
Pkg.add("OnlineStats")
Pkg.add("TimeSeries")
Pkg.add("Measurements")
Pkg.add("Distances")
Pkg.add("DifferentialEquations")
Pkg.add("BenchmarkTools")
Pkg.add("DelimitedFiles")
Pkg.add("Makie")
Pkg.add("GeoMakie")
Pkg.add("GeoJSON")
Pkg.add("Shapefile")
Pkg.add("Genie")

#;conda install matplotlib
#;conda install -c conda-forge cmocean gsw
#;conda install --channel conda-forge erddapy

using Conda
#Conda.add("matplotlib")
#Conda.add("cmocean",channel="conda-forge")
#Conda.add("gsw",channel="conda-forge")


text:
function demangle(name::String)
    return ""
end


text:
#wave = cwts[5]; bc = bcs[3]; β=βs[1]; ave = averagingLengths[1]; n = ns[1]
@testset "Wavelet properties" begin
    bcs = (ContinuousWavelets.PerBoundary(), ContinuousWavelets.SymBoundary(), ContinuousWavelets.ZPBoundary())
    cwts = (cDb2, ContinuousWavelets.dog1, ContinuousWavelets.dog2, ContinuousWavelets.paul1, ContinuousWavelets.morl, ContinuousWavelets.morse)
    βs = (1, 4)
    ns = (128, 1382, 2039)
    averagingLengths = (0, 0.5, 1, 2, 4)
    @testset "length $n, with type $wave, bc $bc, β=$β, ave=$(ave)" for n in ns, wave in cwts, bc in bcs, β in βs, ave in averagingLengths
        wav = ContinuousWavelets.wavelet(wave, β = β, boundary = bc, averagingLength = ave)
        Ŵ, ω = (3, 1)
        with_logger(ConsoleLogger(stderr, Logging.Error)) do
            Ŵ, ω = ContinuousWavelets.computeWavelets(n, wav)
        end
        if wav.waveType isa Union{<:Morlet,<:Paul,<:Morse}
            nFreq, nSpace = ContinuousWavelets.setn(n, wav)
            fullŴ = [Ŵ; zeros(nFreq - 1, size(Ŵ, 2))]
            fullŴ[end-nFreq+1:end, 1] = reverse(Ŵ[:, 1])
            spaceWaves = ifft(fullŴ, 1)
        else
            if bc isa PerBoundary
                nLong = n
            elseif bc isa SymBoundary
                nLong = 2n
            elseif bc isa ZPBoundary
                nLong = 2^(ceil(Int, log2(n)))
            end
            spaceWaves = irfft(Ŵ, nLong, 1)
        end
        # since the signals are centered at the boundary, the area that should be zero is in the middle
        nonSupported = spaceWaves[ceil(Int, n / 2):end-ceil(Int, n / 2), :]
        supported = spaceWaves[[1:floor(Int, n / 2); end-floor(Int, n / 2):end], :]
        # make sure the wavelets don't have significant support out into the mirrored signal
        if !(bc isa PerBoundary) && !(wave isa Paul) && size(Ŵ, 2) > 1
            @test max([norm(nonSupported[:, i], Inf) / norm(supported[:, i], Inf) for i = 2:size(spaceWaves, 2)]...) ≤ 1e-2
        elseif wave isa Paul && size(Ŵ, 2) > 1 # the Paul wavelet decay is just too slow to avoid bleedover, so best to ignore the averaging size
            @test max([norm(nonSupported[:, i], Inf) / norm(supported[:, i], Inf) for i = 2:size(spaceWaves, 2)]...) ≤ 1e-1
        end
        # if the averaging length is 0 and its only averaging at the sizes given, something is wrong
        @test size(Ŵ, 2) > 1 || ave > 0
        if !(wave isa ContOrtho) && size(Ŵ, 2) > 1
            # Guaranteed in a very different way for ContOrtho, and no guarantees made for just averaging
            @test max(abs.(Ŵ[end, :])...) / norm(supported, Inf) ≤ 1e-1
        end
    end
end


text:
using Test
using YaoAPI
using MLStyle
using YaoLocations
using YaoCompiler
using GPUCompiler
using YaoCompiler.Intrinsics
using CompilerPluginTools
using YaoHIR
using ZXCalculus
using YaoArrayRegister

@device function test_a()
    1 => X
    1 => X
end

interp = YaoInterpreter(;options=HardwareFreeOptions(;
        clifford_simplification=true, phase_teleportation=false
    )
)

ir, = code_typed(Intrinsics.apply, (AnyReg, typeof(test_a())); interp)[1]
@test_codeinfo ir begin
    Expr(:invoke, _, apply, _, &(QuoteNode(H)), &(QuoteNode(Locations(1))))::Nothing
    Expr(:invoke, _, apply, _, &(QuoteNode(H)), &(QuoteNode(Locations(1))))::Nothing
end

@device function test_b()
    1 => Rx(2.0)
    1 => Rz(3.0)
end

ir, = code_typed(Intrinsics.apply, (AnyReg, typeof(test_b())); interp)[1]

@test_codeinfo ir begin
    Expr(:invoke, _, apply, _, &(QuoteNode(H)), &(QuoteNode(Locations(1))))::Nothing
    Expr(:invoke, _, apply, _, &(QuoteNode(Rz(2.0))), &(QuoteNode(Locations(1))))::Nothing
    Expr(:invoke, _, apply, _, &(QuoteNode(H)), &(QuoteNode(Locations(1))))::Nothing
    Expr(:invoke, _, apply, _, &(QuoteNode(Rz(3.0))), &(QuoteNode(Locations(1))))::Nothing
end

@device function test_cir()
    5 => H
    @ctrl 4 5 => X
    5 => shift(7 / 4 * π)
    @ctrl 1 5 => X
    5 => T
    @ctrl 4 5 => X
    5 => shift(7 / 4 * π)
    @ctrl 1 5 => X
    4 => T
    5 => T
    @ctrl 1 4 => X
    4 => shift(7 / 4 * π)
    1 => T
    @ctrl 1 4 => X
    @ctrl 4 5 => X
    5 => shift(7 / 4 * π)
    @ctrl 3 5 => X
    5 => T
    @ctrl 4 5 => X
    5 => shift(7 / 4 * π)
    @ctrl 3 5 => X
    4 => T
    5 => T
    @ctrl 3 4 => X
    4 => shift(7 / 4 * π)
    5 => H
    3 => T
    @ctrl 3 4 => X
    @ctrl 4 5 => X
    5 => H
    @ctrl 3 5 => X
    5 => shift(7 / 4 * π)
    @ctrl 2 5 => X
    5 => T
    @ctrl 3 5 => X
    5 => shift(7 / 4 * π)
    @ctrl 2 5 => X
    3 => T
    5 => T
    @ctrl 2 3 => X
    3 => shift(7 / 4 * π)
    5 => H
    2 => T
    @ctrl 2 3 => X
    @ctrl 3 5 => X
    5 => H
    @ctrl 2 5 => X
    5 => shift(7 / 4 * π)
    @ctrl 1 5 => X
    5 => T
    @ctrl 2 5 => X
    5 => shift(7 / 4 * π)
    @ctrl 1 5 => X
    2 => T
    5 => T
    @ctrl 1 2 => X
    2 => shift(7 / 4 * π)
    5 => H
    1 => T
    @ctrl 1 2 => X
    @ctrl 2 5 => X
    @ctrl 1 5 => X
end

ir, = code_typed(Intrinsics.apply, (AnyReg, typeof(test_cir())); interp)[1]

# include("test/emulate.jl")

r = rand_state(5)
cir = test_cir()
target_f = YaoCompiler.compile(
    JLEmulationTarget(), Intrinsics.apply, Tuple{typeof(r),typeof(cir)},
    HardwareFreeOptions(
        clifford_simplification=false, phase_teleportation=false,
    )
)

test_f = YaoCompiler.compile(JLEmulationTarget(), Intrinsics.apply, Tuple{typeof(r),typeof(cir)}, HardwareFreeOptions(
    clifford_simplification=true, phase_teleportation=true,
))

r = rand_state(5)
target_r = copy(r)
test_r = copy(r)

target_f(target_r, cir)
test_f(test_r, cir)

@test test_r ≈ target_r


text:
function naive_filt!(C::SeisChannel;
  fl::Float64=1.0,
  fh::Float64=15.0,
  np::Int=4,
  rp::Int=10,
  rs::Int=30,
  rt::String="Bandpass",
  dm::String="Butterworth"
  )

  T = eltype(C.x)
  fe = 0.5 * C.fs
  low = T(fl / fe)
  high = T(fh / fe)

  # response type
  if rt == "Highpass"
    ff = Highpass(fh, fs=fs)
  elseif rt == "Lowpass"
    ff = Lowpass(fl, fs=fs)
  else
    ff = getfield(DSP.Filters, Symbol(rt))(fl, fh, fs=fs)
  end

  # design method
  if dm == "Elliptic"
    zp = Elliptic(np, rp, rs)
  elseif dm == "Chebyshev1"
    zp = Chebyshev1(np, rp)
  elseif dm == "Chebyshev2"
    zp = Chebyshev2(np, rs)
  else
    zp = Butterworth(np)
  end

  # polynomial ratio
  pr = convert(PolynomialRatio, digitalfilter(ff, zp))

  # zero-phase filter
  C.x[:] = filtfilt(pr, C.x)
  return nothing
end


text:
using Colors
using .Plots
using LinearAlgebra
using RecipesBase
# using ..Client
# using ..Jobs
# using ..Structures

Base.:*(f::Number, r::RGB) = RGB(f * r.r, f * r.b, f * r.g)
Base.:+(r1::RGB, r2::RGB) = RGB(r1.r + r2.r, r1.b + r2.b, r1.g + r2.g)

function blend_color(contribs::Vector, at_colors)
    if any(isnan, contribs)
        return RGB(0.0, 0.0, 0.0)
    end
    if length(contribs) == 1
        return at_colors[1]
    end
    result = weighted_color_mean(normalize([contribs[1], contribs[2]])[1], at_colors[1],
                                 at_colors[2])
    totcontrib = contribs[1] + contribs[2]
    for i in 3:length(contribs)
        result = weighted_color_mean(normalize([totcontrib, contribs[i]])[1], result,
                                     at_colors[i])
        totcontrib += contribs[i]
    end
    return result
end

#@TODO create nice ticks when ks=:relative
# """
# Recipe to plot a Band. If ks is set to relative it calculates the relative offset of the k_points to the middle of the k-point path and puts that on the x-axis.

# Input:    band::Band
# Optional: ks=nothing
# Kwargs:   fermi=0, -> applies fermi level to band eigenvalues before plotting.
# linewidth=2
# """
@recipe function f(band::Band, ks = nothing; fermi = 0, linewidth = 2)
    if ks == :relative_cart
        ks = []
        k_m = band.k_points_cart[div(size(band.k_points_cart)[1] + 1, 2)]
        for k in band.k_points_cart
            push!(ks, norm(k - k_m))
        end
        ks[1:div(length(ks), 2)] = -ks[1:div(length(ks), 2)]
    elseif ks == :relative_cryst
        ks = []
        k_m = band.k_points_cryst[div(size(band.k_points_cryst)[1] + 1, 2)]
        for k in band.k_points_cryst
            push!(ks, norm(k - k_m))
        end
        ks[1:div(length(ks), 2)] = -ks[1:div(length(ks), 2)]
    else
        ks = collect(1:length(band.k_points_cart))
    end
    eigvals = band.eigvals .- fermi

    linewidth --> linewidth
    title --> "Eigenvalues"
    yguide --> "Energy (eV)"
    legend --> false
    return ks, eigvals
end

@recipe function f(bands::Vector{<:Band}, ks = nothing)
    for (i, band) in enumerate(bands)
        @series begin
            band, ks
        end
    end
end

@recipe function f(job::Job, ymin, ymax, occupy_ratio = 0.2; overlap_spin = false, plot_pdos=true)
    palette_ = ismissing(Plots.default(:palette)) ? :default : Plots.default(:palette)
    tc = Plots.plot_color(pop!(plotattributes, :seriescolor,
                                              RGB.(Plots.color_list(Plots.palette(palette_)))))

    bands_label = get(plotattributes, :label, "")
    plt_colors = tc isa Colorant ? repeat([tc], 4) : repeat(tc, 4)
    ylims --> [ymin, ymax]
    gridalpha --> 0.9
    if !any(x -> eltype(x) == QE, job.calculations)
        error("output plotting only implemented for QE jobs.")
    end
    outdat = Client.outputdata(job)
    if !haskey(plotattributes, :fermi)
        frmi = Jobs.readfermi(job, outdat)
    else
        frmi = pop!(plotattributes, :fermi)
    end
    fermi --> frmi
    bands = Jobs.readbands(job, outdat)
    if bands === nothing
        error("No bands found in job $(job.name).")
    end

    # Bands part
    ks = Structures.high_symmetry_kpoints(job.structure)
    tick_vals = Int[]
    tick_syms = String[]
    kpoints = bands isa NamedTuple ? bands.up[1].k_points_cryst : bands[1].k_points_cryst
    for (i, k) in enumerate(kpoints)
        if ks !== nothing
            kpath = ks.kpoints
            for (sym, vec) in kpath
                if vec == k
                    push!(tick_vals, i)
                    push!(tick_syms, " " * string(sym) * " ")
                end
            end
        elseif norm(k) == 0
            push!(tick_vals, i)
            push!(tick_syms, " Γ ")
        end
    end
    if bands isa NamedTuple
        window_ids = map(x -> findall(x) do b
                             min = minimum(b.eigvals) - frmi
                             max = maximum(b.eigvals) - frmi
                             return !((min < ymin && max < ymin) ||
                                      (min > ymax && max > ymax))
                         end, bands)
    else
        window_ids = (findall(bands) do b
                          min = minimum(b.eigvals) - frmi
                          max = maximum(b.eigvals) - frmi
                          return !((min < ymin && max < ymin) || (min > ymax && max > ymax))
                      end,)
    end
    window_ids === nothing && error("No bands inside window")
    # We define a single band plotting series here
    function plot_band(band, color, label, subplot)
        if overlap_spin
            tit = ""
        else
            tit = subplot == 1 ? (length(bands) == 2 ? "Spin up" : "Eigenvalues") :
                  "Spin down"
        end
        @series begin
            xticks --> (tick_vals, tick_syms)
            title --> tit
            yguide := subplot == 1 ? "Energy (eV)" : ""
            serieslabel := label
            subplot := subplot
            seriescolor := color
            legend --> false
            1:length(kpoints), band.eigvals .- frmi
        end
    end

    # PDOS part
    projwfc = Utils.getfirst(x -> Calculations.isprojwfc(x) && haskey(outdat, x.name), job.calculations)
    if projwfc !== nothing && plot_pdos
        if bands isa NamedTuple && !overlap_spin
            doswindow = 3
            layout --> (1, 3)
        else
            doswindow = 2
            layout --> (1, 2)
        end
        states, projbands = outdat[projwfc.name][:states], outdat[projwfc.name][:bands]
        # First we find the amount that all the states appear in the window
        state_occupations = zeros(length(states))
        for ib in 1:(bands isa NamedTuple ? 2 : 1)
            for wid in window_ids[ib]
                b = projbands[wid]
                for ψ in b.extra[:ψ]
                    for is in eachindex(states)
                        state_occupations[is] += ψ[is]
                    end
                end
            end
        end
        # Now we take the most occupied ones, and somehow find out where there's a sudden dropoff
        max_occ = maximum(state_occupations)
        # sorted_occ = sortperm(state_occupations, rev=true)
        # goodids = findall(i -> state_occupations[sorted_occ][i] > occupy_ratio * max_occ, 1:length(state_occupations))
        # ats_orbs = unique(map(x -> (job.structure.atoms[x.atom_id].name, orbital(x.l).name), states[sorted_occ][goodids]))
        goodids = findall(i -> state_occupations[i] > occupy_ratio * max_occ,
                          1:length(state_occupations))
        ats_orbs = unique(map(x -> (job.structure.atoms[x.atom_id].name, Structures.orbital(x.l)),
                              states[goodids]))
        @info "Found $(length(ats_orbs)) atomic orbitals that satisfy the minimum occupation:\n$ats_orbs"

        atom_colors = bands isa NamedTuple ?
                      [plt_colors[2:length(ats_orbs)+1],
                       plt_colors[length(ats_orbs)+2:2*length(ats_orbs)+1]] :
                      [plt_colors[2:length(ats_orbs)+1]]

        bands = bands isa NamedTuple ? bands : [bands]
        band_contribs = [[[zeros(length(ats_orbs)) for i in 1:length(kpoints)]
                          for i1 in 1:length(window_ids[d])] for d in 1:length(bands)]

        @info "Reading pdos files and generating band coloring..."
        for (ia, (atsym, orb)) in enumerate(ats_orbs)
            energies, pd =outdat[projwfc.name][:energies], outdat[projwfc.name][:pdos][atsym][orb]
            #Plots PDOS
            if size(pd, 2) == 2
                @series begin
                    label --> "$(atsym)_$(orb.name)_up"
                    yguide --> ""
                    subplot := doswindow
                    seriescolor := atom_colors[1][ia]
                    title := "DOS"
                    pd[:, 1], energies .- frmi
                end
                @series begin
                    label --> "$(atsym)_$(orb.name)_down"
                    yguide --> ""
                    subplot := doswindow
                    seriescolor := atom_colors[2][ia]
                    title := "DOS"
                    -1 .* pd[:, 2], energies .- frmi
                end
            else
                @series begin
                    label --> "$(atsym)_$(orb.name)"
                    yguide --> ""
                    subplot := doswindow
                    seriescolor := atom_colors[1][ia]
                    title := "DOS"
                    pd, energies .- frmi
                end
            end
            #Calculate band colors
            for (iud, (bnds, contribs)) in enumerate(zip(bands, band_contribs))
                for (ib, b) in enumerate(bnds[window_ids[iud]])
                    for ik in 1:length(kpoints)
                        ibin = findfirst(x -> energies[x] < b.eigvals[ik] <= energies[x+1],
                                         1:length(energies)-1)
                        contribs[ib][ik][ia] += ibin === nothing ? 0.0 : abs(pd[ibin, iud])
                    end
                end
            end
        end
        for contribs in band_contribs
            for ib in contribs
                for ik in 1:length(ib)
                    ib[ik] .= normalize(ib[ik])
                end
            end
        end
        band_colors = [[[blend_color(band_contribs[i][ib][ik], atom_colors[i])
                         for ik in 1:length(kpoints)] for ib in 1:length(window_ids[i])]
                       for i in 1:length(band_contribs)]
        @info "Plotting bands..."
        for (iplt, (bnds, colors)) in enumerate(zip(bands, band_colors))
            if bands_label == ""
                if length(bands) == 2
                    bands_label = iplt == 1 ? "up" : "down"
                else
                    bands_label = ""
                end
            end
            for (ib, (b, c)) in enumerate(zip(bnds[window_ids[iplt]], colors))
                plot_band(b, c, ib == 1 ? bands_label : "", overlap_spin ? 1 : iplt)
            end
        end

        # If no pdos is present
    else
        if bands isa NamedTuple && !overlap_spin
            layout := (1, 2)
        else
            layout := (1, 1)
        end
        @info "Plotting bands..."
        if bands isa NamedTuple
            #loop over up down
            for (iplt, bnds) in enumerate(bands)
                lab = iplt == 1 ? "up" : "down"
                color = iplt == 1 ? :blue : :red
                #loop over bands inside window
                for (ib, b) in enumerate(bnds[window_ids[iplt]])
                    plot_band(b, color, ib == 1 ? lab : "", overlap_spin ? 1 : iplt)
                end
            end
        else
            for (ib, b) in enumerate(bands[window_ids[1]])
                plot_band(b, plt_colors[mod1(ib, length(plt_colors))], "", 1)
            end
        end
    end
end


text:
#=
  Biased coin

  From cplint: http://cplint.eu/example/inference/coin.pl
  """
  Throwing a coin with uncertainty on its fairness, from
  J. Vennekens, S. Verbaeten, and M. Bruynooghe. Logic programs with annotated
  disjunctions. In International Conference on Logic Programming,
  volume 3131 of LNCS, pages 195-209. Springer, 2004.
  """

  Model 1: Probability of tossing a head (no observation)
  Distributions of variable coinHead (num:0)
  1.00000 =>    5049  (0.504900)
  0.00000 =>    4951  (0.495100)

  Model 2: Probability of tossing a head given that we know it's biased
  Distributions of variable coinHead (num:0)
  1.00000 =>    5712  (0.571200)
  0.00000 =>    4288  (0.428800)

  Cf ~/cplint/coin.pl
     ~/blog/biased_coin.blog
     ~/webppl/biased_coin.wppl

=#

using Turing
include("jl_utils.jl")

@model function biased_coin(experiment="observe biased")
    fair = 1
    biased = 2

    coinType ~ Categorical([0.9,0.1]) # fair,biased
    head = 1
    tail = 2
    coinResult ~ coinType == fair ? Categorical([0.5,0.5]) :
            Categorical([0.6,0.4])

    if experiment == "observe biased"
        true ~ Dirac(coinType == biased)
    end

    coinHead ~ Dirac(coinResult == head)
    coinTail ~ Dirac(coinResult == tail)

end

println("Model 1: Probability of tossing a head (no observation)")
model = biased_coin("no observation")

num_chns = 4
# chs = sample(model, Prior(), 1000)
# chs = sample(model, MH(), 10_000)
# chs = sample(model, PG(15), 10_000)
# chs = sample(model, IS(), 10_000)
chs = sample(model, SMC(), 10_000)
# chs = sample(model, SMC(), MCMCThreads(), 10_000, num_chs)

# chs = sample(model, SGLD(), 10_000)
# chs = sample(model,NUTS(), 10_000)
# chs = sample(model,HMC(0.01,5), 10_000)

# display(chs)
# display(plot(chs))

show_var_dist_pct(chs,:coinHead)

println("\nModel 2: Probability of tossing a head given that we know it's biased")
model = biased_coin("observe biased")

num_chns = 4
# chs = sample(model, Prior(), 1000)
# chs = sample(model, MH(), 10_000)
# chs = sample(model, PG(15), 10_000)
# chs = sample(model, IS(), 10_000)
chs = sample(model, SMC(), 10_000)
# chs = sample(model, SMC(), MCMCThreads(), 10_000, num_chs)

# chs = sample(model, SGLD(), 10_000)
# chs = sample(model,NUTS(), 10_000)
# chs = sample(model,HMC(0.01,5), 10_000)

# display(chs)
# display(plot(chs))

show_var_dist_pct(chs,:coinHead)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Test
using JuMP
using DisjunctiveProgramming
using Symbolics

@testset "linear constraints" begin

    function minimal()
        m = Model()
        @variable(m, -1<=x<=10)

        @constraint(m, con1, x<=3)
        @constraint(m, con2, 0<=x)
        @constraint(m, con3, x<=9)
        @constraint(m, con4, 5<=x)

        @disjunction(m,(con1,con2),con3,con4,reformulation=:CHR,name=:y)

        @test true
    end

    function simple_example(reform)
        m = Model()
        @variable(m, -10<=x<=10)
        @constraint(m, con1, x<=-1)
        @constraint(m, con2, 1<=x)

        if (reform == :BMR)
            @disjunction(m, con1, con2, reformulation=:BMR, name=:y)
        elseif (reform == :CHR)
            @disjunction(m, con1, con2, reformulation=:CHR, name=:y)
        end
        return m
    end

    function robustness()
        unreg = m -> unregister(m, :original_model_variables)

        function fresh_model()
            m = Model()
            @variable(m, -10<=x<=10)
            @constraint(m, con1, x<=-1)
            @constraint(m, con2, 1<=x)
            return m, con1, con2
        end

        # not enough constraints
        m = fresh_model()[1]
        # @test_throws DomainError @disjunction(m, (con1, con2), reformulation=:BMR, name=:y)
        unreg(m)
        # @test_throws DomainError @disjunction(m, con1, reformulation=:BMR, name=:y)

        # Big-M reformulation without variable bounds defined
        # should only work if user specifies M value
        m = Model()
        @variable(m, x)
        @constraint(m, con1, x<=-1)
        @constraint(m, con2, 1<=x)
        @test @disjunction(m, con1, con2, reformulation=:BMR, name=:y, M=11) == nothing
        unreg(m)
        @test_throws ErrorException @disjunction(m, con1, con2, reformulation=:BMR, name=:y)

        # CHR reformulation without variable bounds defined should fail
        m = Model()
        @variable(m, x)
        @constraint(m, con1, x<=-1)
        @constraint(m, con2, 1<=x)
        @test_throws AssertionError @disjunction(m, con1, con2, reformulation=:CHR, name=:y)

        # empty constraints on one side of disjunction
        m, con1, con2 = fresh_model()
        @test @disjunction(m, con1, nothing, reformulation=:BMR, name=:y) == nothing
        m, con1, con2 = fresh_model()
        @test @disjunction(m, nothing, con2, reformulation=:CHR, name=:z) == nothing
        m, con1, con2 = fresh_model()
        @test @disjunction(m, (con1, con2), nothing, reformulation=:BMR, name=:y) == nothing

    end

    function model_BMR_valid(m)
        # Expecting to see following constraints:
        # x <= -1 + 11 * (1 - y1)
        # -x <= -1 + 11 * (1 - y2)
        # y1 + y2 = 1
        # -10<=x<=10
        cons = (L -> all_constraints(m, L...)).(list_of_constraint_types(m))
        @test sum(length.(cons)) == 5 + 2 # add 2 for y1/y2 binary statements

        con_strings = replace.(constraints_string(REPLMode, m), "[" => "",
            "]" => "", "+ y" => "+y", " y" => "* y",
            "con1 :" => "", "con2 :" => "")
        cons_actual = Meta.parse.(con_strings[1:5])
        cons_expected = [:(y1 + y2 == 1), :(x + 11 * y1 <= 10),
                         :(-x + 11 * y2 <= 10), :(x >= -10), :(x <= 10)]
        matches = sum(map(i->isequal(i[1], i[2]),
            Base.product(cons_expected, cons_actual)))
        @test matches == length(cons_actual)
    end

    function model_CHR_valid(m)
        # Expecting to see following constraints:
        # x_y1 <= -1 * y1
        # -x_y2 <= -1 * y2
        # y1 + y2 = 1
        # x = x_y1 + x_y2
        # -10 * y1 <= x_y1 <= 10 * y1
        # -10 * y2 <= x_y2 <= 10 * y2
        # -10 <= x <= 10
        # -10 <= x_y1 <= 10
        # -10 <= x_y2 <= 10
        cons = (L -> all_constraints(m, L...)).(list_of_constraint_types(m))
        @test sum(length.(cons)) == 14 + 2 # add 2 for y1/y2 binary statements

        con_strings = replace.(constraints_string(REPLMode, m), "[" => "",
            "]" => "", "+ y" => "+y", " y" => "* y",
            "con1 :" => "", "con2 :" => "")
        println.(con_strings)
        cons_expected = [:(y1 + y2 == 1), :(x - x_1 - x_2 == 0), :(x >= -10),
                     :(x_1 + y1 <= 0), :(-x_2 + y2 <= 0), :(x <= 10),
                     :(x_1 >= -10), :(x_1 <= 10), :(x_2 >= -10), :(x_2 <= 10),
                     :(-10 * y1 - x_1 <= 0), :(-10 * y1 + x_1 <= 0),
                     :(-10 * y2 - x_2 <= 0), :(-10 * y2 + x_2 <= 0)]
        cons_actual = Meta.parse.(con_strings[1:14])
        matches = sum(map(i->isequal(i[1], i[2]),
            Base.product(cons_expected, cons_actual)))
        @test matches == length(cons_actual)
    end

    robustness()

    m1 = simple_example(:BMR)
    model_BMR_valid(m1)

    m2 = simple_example(:CHR)
    model_CHR_valid(m2)

end


text:
function f1()
	 m = 10000
	 for i=1:m
	     k = 6
	     d = k + 2
	     str = string(d,"test")
	 end
end

function f2()
	 m = 100000
	 for i=1:m
	     k = 6
	     d = k + 2
	     str = string(d,"test")
	 end
end

#run section

function runner()
	 local n = 500

	 for i = 1:n
	     f1()
	     f2()
	 end
end

max_stars_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_issues_repo_name: wchmielarz/NobodyDoesItLikeJuliet

max_forks_repo_name: wchmielarz/NobodyDoesItLikeJuliet


text:
using Makie
scene=Scene()
scene=lines(rand(10))

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
# Copying
# =======
#
# Copying methods for biological sequences.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

function Base.copy!{A}(seq::BioSequence{A}, doff::Integer,
                       src::Vector{UInt8},  soff::Integer, len::Integer)
    datalen = seq_data_len(A, len)
    return seq
end

function Base.copy!{A}(dst::BioSequence{A}, src::BioSequence{A})
    return copy!(dst, 1, src, 1)
end

function Base.copy!{A}(dst::BioSequence{A}, doff::Integer,
                       src::BioSequence{A}, soff::Integer)
    return copy!(dst, doff, src, soff, length(src) - soff + 1)
end

function Base.copy!{A}(dst::BioSequence{A}, doff::Integer,
                       src::BioSequence{A}, soff::Integer, len::Integer)
    checkbounds(dst, doff:doff+len-1)
    checkbounds(src, soff:soff+len-1)

    if dst.shared || (dst === src && doff > soff)
        orphan!(dst, length(dst), true)
    end

    id = bitindex(dst, doff)
    is = bitindex(src, soff)
    rest = len * bitsof(A)

    while rest > 0
        # move `k` bits from `src` to `dst`
        x = dst.data[index(id)]
        y = src.data[index(is)]
        if offset(id) < offset(is)
            y >>= offset(is) - offset(id)
            k = min(64 - offset(is), rest)
        else
            y <<= offset(id) - offset(is)
            k = min(64 - offset(id), rest)
        end
        m = mask(k) << offset(id)
        dst.data[index(id)] = y & m | x & ~m

        id += k
        is += k
        rest -= k
    end

    return dst
end

# Actually, users don't need to create a copy of a sequence.
function Base.copy{A}(seq::BioSequence{A})
    newseq = BioSequence{A}(seq, 1:endof(seq))
    orphan!(newseq, length(seq), true)  # force orphan!
    @assert newseq.data !== seq.data
    return newseq
end

# Encoded copying
# ---------------

function encode_copy!{A}(dst::BioSequence{A},
                         src::Union{AbstractVector,AbstractString})
    return encode_copy!(dst, 1, src, 1)
end

function encode_copy!{A}(dst::BioSequence{A},
                         doff::Integer,
                         src::Union{AbstractVector,AbstractString},
                         soff::Integer)
    return encode_copy!(dst, doff, src, soff, length(src) - soff + 1)
end

function encode_copy!{A}(dst::BioSequence{A},
                         doff::Integer,
                         src::Union{AbstractVector,AbstractString},
                         soff::Integer,
                         len::Integer)
    if soff != 1 && isa(src, AbstractString) && !isascii(src)
        throw(ArgumentError("source offset ≠ 1 is not supported for non-ASCII string"))
    end

    checkbounds(dst, doff:doff+len-1)
    if length(src) < soff + len - 1
        throw(ArgumentError("source string does not contain $len elements from $soff"))
    end

    orphan!(dst)
    next = bitindex(dst, doff)
    stop = bitindex(dst, doff + len)
    i = soff
    while next < stop
        x = UInt64(0)
        j = index(next)
        while index(next) == j && next < stop
            char, i = Base.next(src, i)
            x |= enc64(dst, convert(Char, char)) << offset(next)
            next += bitsof(A)
        end
        dst.data[j] = x
    end
    return dst
end

function encode_copy!{A<:Union{DNAAlphabet{4},RNAAlphabet{4}}}(
        dst::BioSequence{A}, doff::Integer,
        src::AbstractVector{UInt8}, soff::Integer, len::Integer)
    checkbounds(dst, doff:doff+len-1)
    if length(src) < soff + len - 1
        throw(ArgumentError("source string does not contain $len elements from $soff"))
    end

    orphan!(dst)
    charmap = A <: DNAAlphabet ? char_to_dna : char_to_rna
    i = soff
    next = bitindex(dst, doff)
    stop = bitindex(dst, doff + len)

    # head
    if offset(next) != 0
        for d in 0:div(64 - offset(next), 4)-1
            dst[doff+d] = charmap[src[i+d]+1]
        end
        i += div(64 - offset(next), 4)
        next += 64 - offset(next)
    end

    # body
    D = 16
    while next < (stop - offset(stop))
        x::UInt64 = 0
        check = 0x00
        @inbounds for d in 0:D-1
            y = reinterpret(UInt8, charmap[src[i+d]+1])
            x |= UInt64(y) << 4d
            check |= y
        end
        if check & 0x80 != 0
            # invalid byte(s) is detected
            for d in 0:D-1
                if !isvalid(charmap[src[i+d]+1])
                    error("cannot encode $(src[i+d])")
                end
            end
        end
        dst.data[index(next)] = x
        i += D
        next += 64
    end

    # tail
    for d in 0:div(stop - next, 4)-1
        dst[doff+i-soff+d] = charmap[src[i+d]+1]
    end

    return dst
end


text:
function sift!(Eav, decomp, d_osf, nsifts=5)
    N = length(decomp)
    e1 = zeros(N)
    e2 = zeros(N)
    avg = zeros(N)
    w = max(div(d_osf-1, 2), 3)
    if iseven(w)
        w += 1
    end
    for j in 1:nsifts
        stream_minmax(e1, e2, decomp, d_osf)
        e1 .= moving_average(e1, d_osf)
        e2 .= moving_average(e2, d_osf)
        avg .= (e1 .+ e2)./2
        avg .= moving_average(avg, d_osf)
        Eav .= Eav .+ avg
        decomp .-= avg
    end
    return nothing
end

function sEMD(signal;maximfs=10,nsifts=2)
    N = length(signal)
    a_new = copy(signal)
    imfs = ElasticArray{Float64}(undef,N,0)
    Eav = zeros(N)
    for i in 1:maximfs
        exs = find_extrema_count(a_new)
        d_osf =  div(N, div(exs,2))
        iseven(d_osf) && (d_osf+=1)

        if d_osf >= div(N,3) || exs < 5
            append!(imfs, a_new)
            break
        end
        fill!(Eav,0)
        sift!(Eav, a_new, d_osf, nsifts)
        append!(imfs, a_new)
        a_new .= Eav
    end
        return convert(Array,imfs)
end


text:
@doc read(joinpath(dirname(@__DIR__), "README.md"), String)
module PkgTemplates

using Base: active_project, contractuser

using Dates: month, today, year
using InteractiveUtils: subtypes
using LibGit2: LibGit2, GitConfig, GitReference, GitRemote, GitRepo, delete_branch
using Pkg: Pkg, TOML, PackageSpec
using REPL.TerminalMenus: MultiSelectMenu, RadioMenu, request
using UUIDs: uuid4

using Mustache: render
using Parameters: @with_kw_noshow

using Mocking

export
    Template,
    AppVeyor,
    BlueStyleBadge,
    CirrusCI,
    Citation,
    Codecov,
    ColPracBadge,
    CompatHelper,
    Coveralls,
    Develop,
    Documenter,
    DroneCI,
    Git,
    GitHubActions,
    GitLabCI,
    License,
    Logo,
    NoDeploy,
    PkgEvalBadge,
    ProjectFile,
    Readme,
    RegisterAction,
    Secret,
    SrcDir,
    TagBot,
    Tests,
    TravisCI

"""
Plugins are PkgTemplates' source of customization and extensibility.
Add plugins to your [`Template`](@ref)s to enable extra pieces of repository setup.

When implementing a new plugin, subtype this type to have full control over its behaviour.
"""
abstract type Plugin end

include("template.jl")
include("plugin.jl")
include("show.jl")
include("interactive.jl")
include("deprecated.jl")

# Run some function with a project activated at the given path.
function with_project(f::Function, path::AbstractString)
    proj = active_project()
    try
        Pkg.activate(path)
        f()
    finally
        proj === nothing ? Pkg.activate() : Pkg.activate(proj)
    end
end

end


text:
const PEDESTRIAN_DEF =  VehicleDef(AgentClass.PEDESTRIAN, 1.0, 1.0)

function AutoViz.render!(rendermodel::RenderModel, env::CrosswalkEnv)
    roadway = gen_straight_roadway(2, env.params.roadway_length)
    render!(rendermodel, roadway)

    curve = env.crosswalk.curve
    n = length(curve)
    pts = Array{Float64}(2, n)
    for (i,pt) in enumerate(curve)
        pts[1,i] = pt.pos.x
        pts[2,i] = pt.pos.y
    end

    add_instruction!(rendermodel, render_dashed_line, (pts, colorant"white", env.crosswalk.width, 1.0, 1.0, 0.0, Cairo.CAIRO_LINE_CAP_BUTT))
    obs = env.obstacles[1]
    for obs in env.obstacles
        pts = Array{Float64}(2, obs.npts)
        for (i, pt) in enumerate(obs.pts)
            pts[1,i] = pt.x
            pts[2,i] = pt.y
        end

        add_instruction!(rendermodel, render_fill_region, (pts, colorant"gray"))
    end

    return rendermodel
end

function animate_record(rec::SceneRecord, overlay::SceneOverlay = EmptyOverlay, save_index::Int64=-1)
    duration =rec.nframes*config.sim_dt
    fps = Int(1/config.sim_dt)
    function render_rec(t, dt)
        frame_index = Int(floor(t/dt)) + 1
        scene = rec[frame_index-nframes(rec)]
        a = render(scene, env, [overlay], cam=cam, car_colors = get_colors(scene))
        if save_index == frame_index
            write_to_png(a,"save$frame_index.png")
        end
        return a
    end
    return duration, fps, render_rec
end

function get_colors(scene::Scene, ego_id::Int64 = 1)
    colors = Dict{Int64, Colorant}()
    for veh in scene
        colors[veh.id] = veh.id == ego_id ?  COLOR_CAR_EGO : COLOR_CAR_OTHER
    end
    return colors
end

function animate_record(rec::SceneRecord)
    duration =rec.nframes*config.sim_dt
    fps = Int(1/config.sim_dt)
    function render_rec(t, dt)
        frame_index = Int(floor(t/dt)) + 1
        return render(rec[frame_index-nframes(rec)], env, cam=cam)
    end
    return duration, fps, render_rec
end

# mutable struct WaitAndGoOverlay <: SceneOverlay
#     Δt::Float64
#     model::WaitAndGoModel
#     ttc::Float64

#     target_id::Int
#     verbosity::Int
#     color::Colorant
#     font_size::Int

#     function WaitAndGoOverlay( Δt::Float64, model::WaitAndGoModel, ttc::Float64 = Inf,
#         target_id::Int=1, verbosity::Int=1;
#         color::Colorant=colorant"white",
#         font_size::Int=20,
#         )

#         new(Δt, model, ttc, target_id, verbosity, color,font_size)
#     end
# end

# function reset_overlay!(overlay::WaitAndGoOverlay, ego::Vehicle)
#     reset_model!(overlay.model, ego)
# end

# function AutoViz.render!(rendermodel::RenderModel, overlay::WaitAndGoOverlay,
#                         scene::Scene, env::CrosswalkEnv)

#     # render overlay
#     font_size = overlay.font_size
#     text_y = font_size
#     text_y_jump = round(Int, font_size*1.2)

#     ego = scene[findfirst(scene, overlay.target_id)]

#     # display ego id
#     add_instruction!( rendermodel, render_text, (@sprintf("id = %d", overlay.target_id), 10, text_y, font_size, overlay.color), incameraframe=false)
#         text_y += text_y_jump
#     # display x position
#     add_instruction!( rendermodel, render_text, (@sprintf("x = %10.3f m", ego.state.posG.x), 10, text_y, font_size, overlay.color), incameraframe=false)
#         text_y += text_y_jump
#     # display ego speed
#     add_instruction!( rendermodel, render_text, (@sprintf("v = %10.3f m/s", ego.state.v), 10, text_y, font_size, overlay.color), incameraframe=false)
#     text_y += text_y_jump

#     # display acceleration
#     add_instruction!( rendermodel, render_text, (@sprintf("a = %10.3f m/s^2", overlay.model.a), 10, text_y, font_size, overlay.color), incameraframe=false)
#     text_y += text_y_jump
#     # display ttc
#     ttc_min = Inf
#     for veh in scene
#         if veh.id != overlay.target_id
#             ttc = (env.params.lane_width - veh.state.posG.y)/veh.state.v
#             if 0 < ttc < ttc_min
#                 ttc_min = ttc
#             end
#         end
#     end
#     overlay.ttc = ttc_min

#     add_instruction!( rendermodel, render_text, (@sprintf("TTC = %10.3f s", overlay.ttc), 10, text_y, font_size, overlay.color), incameraframe=false)
#     text_y += text_y_jump

#     # display state
#     state = "$(overlay.model.policy.reaching)$(overlay.model.policy.wait)$(overlay.model.policy.go)"
#     if overlay.model.policy.reaching
#         state = "reaching"
#     elseif overlay.model.policy.go
#         state = "go"
#     elseif overlay.model.policy.wait
#         state = "wait $(overlay.model.policy.N)"
#     end
#     add_instruction!( rendermodel, render_text, (state, 10, text_y, font_size, overlay.color), incameraframe=false)
#     text_y += text_y_jump

#     #display observation
#     for veh in overlay.model.o
#         if veh.id == overlay.target_id
#             continue
#         end
#         color=RGBA(75./255, 66./255, 244./255, 0.5)
#         x, y, θ = veh.state.posG.x, veh.state.posG.y, veh.state.posG.θ
#         add_instruction!(rendermodel, render_vehicle, (x, y, θ, PEDESTRIAN_DEF.length, PEDESTRIAN_DEF.width, color))
#     end

#     # Update overlay
#     observe!(overlay.model, scene, env.roadway, overlay.target_id)
#     rendermodel
# end

mutable struct QMDPOverlay <: SceneOverlay
    model::AVDriver
    ttc::Float64
    target_id::Int
    verbosity::Int
    color::Colorant
    font_size::Int

    function QMDPOverlay(model::AVDriver, ttc::Float64 = Inf,
        target_id::Int=1, verbosity::Int=1;
        color::Colorant=colorant"white",
        font_size::Int=20,
        )

        new(model, ttc, target_id, verbosity, color,font_size)
    end
end

function reset_overlay!(overlay::QMDPOverlay, ego::Vehicle)
    reset_model!(overlay.model, ego)
end

function AutoViz.render!(rendermodel::RenderModel, overlay::QMDPOverlay, scene::Scene, env::CrosswalkEnv)

    # render overlay
    font_size = overlay.font_size
    text_y = font_size
    text_y_jump = round(Int, font_size*1.2)

    ego = scene[findfirst(overlay.target_id, scene)]

    # display ego id
    add_instruction!( rendermodel, render_text, (@sprintf("Ego car state"), 10, text_y, font_size, overlay.color), incameraframe=false)
        text_y += text_y_jump
    # display x position
    add_instruction!( rendermodel, render_text, (@sprintf("x = %10.3f m", ego.state.posG.x), 10, text_y, font_size, overlay.color), incameraframe=false)
        text_y += text_y_jump
    # display ego speed

    add_instruction!( rendermodel, render_text, (@sprintf("v = %10.3f m/s", ego.state.v), 10, text_y, font_size, overlay.color), incameraframe=false)
    text_y += text_y_jump
    # display acceleration
    add_instruction!( rendermodel, render_text, (@sprintf("a = %10.3f m/s^2", overlay.model.a), 10, text_y, font_size, overlay.color), incameraframe=false)
    text_y += text_y_jump
    # # display ttc
    # ttc_min = Inf
    # for veh in scene
    #     if veh.id != overlay.target_id
    #         ttc = (env.params.lane_width - veh.state.posG.y)/veh.state.v
    #         if 0 < ttc < ttc_min
    #             ttc_min = ttc
    #         end
    #     end
    # end
    # overlay.ttc = ttc_min
    # add_instruction!( rendermodel, render_text, (@sprintf("TTC = %10.3f s", overlay.ttc), 10, text_y, font_size, overlay.color), incameraframe=false)
    # text_y += text_y_jump

    #Display ray
    for oid in keys(overlay.model.o)
        if oid == OFF_KEY; continue; end
        oped = overlay.model.o[oid].ped
        ego_front = get_front(ego)
        line_pts = [oped.posG.x oped.posG.y; ego_front.x ego_front.y]'
        add_instruction!( rendermodel, render_line, (line_pts, COLOR_CAR_EGO, 0.1))
    end

    for veh in scene
        if veh.id == 1 || veh.id ∈ keys(overlay.model.o); continue; end
        ego_front = get_front(ego)
        line_pts = [veh.state.posG.x veh.state.posG.y; ego_front.x ego_front.y]'
        add_instruction!( rendermodel, render_line, (line_pts, RGB(1.0, 0., 0.), 0.1))
    end

    #Display state estimate
    for id in keys(overlay.model.b)
        max, indmax = fargmax(overlay.model.b[id].p)
        for (i,s) in enumerate(overlay.model.b[id].it)
            prob = overlay.model.b[id].p[i]
            if prob > 1e-3
                color=RGBA(75./255, 66./255, 244./255, overlay.model.b[id].p[i])

#             overlay.model.b[id].p[i]
#             println(overlay.model.b[id].p[i])
                p = overlay.model.b[id].it[i].ped.posG
                # if haskey(overlay.model.o, id)
                #     px = overlay.model.o[id].ped.posG.x
                # else
                #     px = p.x
                # end
                px = p.x
                add_instruction!(rendermodel, render_vehicle, (px, p.y, p.θ, PEDESTRIAN_DEF.length, PEDESTRIAN_DEF.width, color))
            end
        end
    end

    # display labels
    x_ego_lab = 2.5
    y_ego_lab = -5.
    θ_ego_lab = 0.
    add_instruction!(rendermodel, render_vehicle, (x_ego_lab, y_ego_lab, θ_ego_lab,
                                                   VehicleDef().length, VehicleDef().width, COLOR_CAR_EGO))
    text_y += 12*text_y_jump
    add_instruction!(rendermodel, render_text, (@sprintf("Ego car"), 120, text_y, font_size, overlay.color), incameraframe=false)
    x_other_lab = 2.5
    y_other_lab = -8
    θ_other_lab = 0.
    add_instruction!(rendermodel, render_vehicle, (x_other_lab, y_other_lab, θ_other_lab,
                                                  PEDESTRIAN_DEF.length, PEDESTRIAN_DEF.width, COLOR_CAR_OTHER))
    text_y += 2.*text_y_jump
    add_instruction!(rendermodel, render_text, (@sprintf("Pedestrian"), 120, text_y, font_size, overlay.color), incameraframe=false)
    text_y += text_y_jump
    add_instruction!(rendermodel, render_text, (@sprintf("ground truth"), 120, text_y, font_size, overlay.color), incameraframe=false)
    x_belief_lab = 2.5
    y_belief_lab = -11
    θ_belief_lab = 0.
    add_instruction!(rendermodel, render_vehicle, (x_belief_lab, y_belief_lab, θ_belief_lab,
                                                 PEDESTRIAN_DEF.length, PEDESTRIAN_DEF.width, RGB(75./255, 66./255, 244./255)))
    text_y += 2.*text_y_jump
    add_instruction!(rendermodel, render_text, (@sprintf("Belief"), 120, text_y, font_size, overlay.color), incameraframe=false)

    add_instruction!(rendermodel, render_text, (@sprintf("Absent state"), 300, 400, font_size, overlay.color), incameraframe=false)

    add_instruction!( rendermodel, render_line, ([35. -5; 36 -5]', RGB(1.0, 0., 0.), 0.1))
    add_instruction!(rendermodel, render_text, (@sprintf("occluded target"), 730, 405, font_size, overlay.color), incameraframe=false)

    add_instruction!( rendermodel, render_line, ([35. -7; 36 -7]', COLOR_CAR_EGO, 0.1))
    add_instruction!(rendermodel, render_text, (@sprintf("visible target"), 730, 445, font_size, overlay.color), incameraframe=false)

    # Display everyones ID
    # for veh in scene
    #     add_instruction!(rendermodel, render_text, ("$(veh.id)", veh.state.posG.x, veh.state.posG.y, font_size, colorant"black"), incameraframe=true)
    # end
    # Update overlay
    observe!(overlay.model, scene, env.roadway, overlay.target_id)
    rendermodel
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
Base.@kwdef struct FakeDSOX4034A <: Oscilloscope
    num_samples = 65104
end

function initialize(model::Type{FakeDSOX4034A})
    return Instr{model}(model(), "", TCPSocket(), true)
end

function initialize(model::FakeDSOX4034A)
    return Instr{typeof(model)}(model, "", TCPSocket(), true)
end

function get_data(instr::Instr{FakeDSOX4034A}, ch::Vector{Int})
    for num in ch
        if num < 1 || num > 4
            error("$num is not a valid channel")
        end
    end
    map(c->get_data(instr, c), ch)
end

num_samples(i::Instr{FakeDSOX4034A}) = i.model.num_samples

function get_data(instr::Instr{FakeDSOX4034A}, ch::Int; scope_stats=false)
    samples = num_samples(instr)
    info = ScopeInfo(
        "8bit",
        "Normal",
        samples,
        7.68e-8,
        -0.0025,
        0.0,
        0.0167364,
        1.28425,
        128.0,
        "",
        "",
        "",
        ch
    )

    volt = if ch == 1
            map(sin, collect(range(0, stop=6pi, length=samples))) .* V
        elseif ch == 2
            map(cos, collect(range(0, stop=6pi, length=samples))) .* V
        elseif ch == 3
            map(x->abs(sin(x)) / -3, collect(range(0, stop=16pi, length=samples))) .* V
        else
            map(x->abs(sin(x)) / 3, collect(range(0, stop=16pi, length=samples))) .* V
    end
    time = collect(range(-0.0025, stop=0.0025, length=samples)) * s
    return ScopeData(info, volt, time)
end


text:
include(joinpath(@__DIR__, "fig1.jl"))
using FoodCachingExperiments, Random, FoodCachingModels, Statistics, DataFrames
# This is a hack, because the labeling in Clayton99C_exp1 is wrong
FoodCachingExperiments.EXPERIMENTS[:Clayton99C_exp1].tests["interactions"].labels[2] = "missing"
import FoodCachingExperiments: EXPERIMENTS, bsave, bload

function plot_sigresults(sigresults;
                         measure = :ratio_same_rejects_mean)
    tmp = combine(groupby(sigresults, [:model, :experiment]),
                  [:samebinarysignificance, :meansamebinarysigA,
                   :meansamebinarysig, :atleasttwo, :ratio_same_rejects] .=> mean,
                 )
    sort!(tmp, :experiment, by = x -> findfirst(==(x), experiment_order))
    sort!(tmp, :model, by = x -> findfirst(==(x), model_order))
    groups = groupby(tmp, :model)
    ys = [getproperty(g, measure)  for g in groups]
    ns = [sum(rejects(EXPERIMENTS[e])) for e in experiment_order
          if e ∉ (:Clayton0103, :Watanabe05)]
    pe = 1:11
    me = 12:22
    se = 23:27
    ce = 28
    reverse(x) = 32 - x
    xtick = reverse.([pe; me .+ 1; se .+ 2; ce .+ 3])
    f1 = @pgf Axis({xtick = xtick, xticklabels = toname.(groups[1].experiment),
               xticklabel_style = {rotate = 60, anchor = "east"},
               width = "16cm", height = "6cm",
               cycle_list = [
                             "{thick, morange, mark = square}",
                             "{thick, mred, mark = square}",
                             "{thick, mblue, mark = square}",
                             "{thick, mgreen, mark = square}",
                             "{thick, gray!60!black, mark = square}"
                            ],
              },
              [Plot(Coordinates(xtick, y)) for y in ys]...)

    f2 = @pgf Axis({cycle_list = [
                                 "{thick, mark = square*, dotted, mark size = 3,
                                 mark options = {solid, fill = white}}",
                                 "{thick, dashed, mark = diamond*,
                                   mark options = {solid, fill =  white}, mark size = 3}",
                                 "{thick, mark = triangle*, mark size = 3,
                                   mark options = {fill = white}}",
                                ],
                   ylabel = "avg. reproducibility [\\%]",
               xticklabel_style = {rotate = 60, anchor = "east"},
               height = "4.8cm", width = "3.8cm",
               ymax = 105,
               ytick = 0:20:100,
               legend_entries = {"satiety experiments", "memory experiments",
                                 "planning experiments"},
               "every axis legend/.append style" = {at = "{(-.4, 1)}",
                                                    anchor = "south west"},
               xtick = 1:5, font = raw"\scriptsize",
               xticklabels = ["\\color{morange}Planning-By-Replay",
                              "\\color{mred}Plastic Caching",
                              "\\color{mblue}No-Plasticity",
                              "\\color{mgreen}No-Plasticity-No-Memory",
                              "\\parbox{3.4cm}{\\color{gray!60!black}\\raggedleft No-Plasticity-No-Memory-\\\\No-Motivational-Control}"],
              },
                  [Plot(Coordinates(1:5, [100*sum(y[x] .* 1/length(x)) for y in ys]))
               for x in Iterators.reverse([pe, me, se])]...)
    (f1, f2)
end

rm_test_locals.(values(EXPERIMENTS))

sigresults1 = run_sigtests(best, Ninner = 1, Nouter = 50)
bsave(joinpath(datapath, "sigresults1"), Dict(:sigresults1 => sigresults1))
sigresults10 = run_sigtests(best, Ninner = 10, Nouter = 50)
bsave(joinpath(datapath, "sigresults10"), Dict(:sigresults10 => sigresults10))

p1indi, p1all = plot_sigresults(sigresults1)
p10indi, p10all = plot_sigresults(sigresults10)

pgfsave(joinpath(figpath, "sig_summary_all1.tikz"), p1all)
pgfsave(joinpath(figpath, "sig_summary_all10.tikz"), p10all)


text:
isdefined(Base, :__precompile__) && __precompile__()

"""
  `SimJulia`

Main module for SimJulia.jl – a combined continuous time / discrete event process oriented simulation framework for Julia.
"""
module SimJulia
  using Base.Collections, Base.Dates

  import Base.show, Base.isless, Base.run, Base.now, Base.schedule, Base.yield, Base.get, Base.interrupt
  import Base.&, Base.|, Base.+, Base.==

  export Event, Process, Simulation
  export Container, Resource, Store
  export run, append_callback, schedule, schedule!
  export timeout
  export yield, interrupt
  export now, value, state
  export (&), (|)
  export get, put, release, request

  include("types.jl")
  include("exceptions.jl")
  include("events.jl")
  include("process.jl")
  include("time.jl")
  include("simulation.jl")
  include("utils.jl")
  include("resources/base.jl")
  include("resources/containers.jl")
  include("resources/stores.jl")
end


text:
# Autogenerated wrapper script for grep_jll for i686-w64-mingw32
export grep

using Libiconv_jll
JLLWrappers.@generate_wrapper_header("grep")
JLLWrappers.@declare_executable_product(grep)
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll)
    JLLWrappers.@init_executable_product(
        grep,
        "bin\\grep.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function HypergeometricDistribution(arg0::jint, arg1::jint, arg2::jint)
    return HypergeometricDistribution((jint, jint, jint), arg0, arg1, arg2)
end

function cumulative_probability(obj::HypergeometricDistribution, arg0::jint)
    return jcall(obj, "cumulativeProbability", jdouble, (jint,), arg0)
end

function get_number_of_successes(obj::HypergeometricDistribution)
    return jcall(obj, "getNumberOfSuccesses", jint, ())
end

function get_numerical_mean(obj::HypergeometricDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::HypergeometricDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_population_size(obj::HypergeometricDistribution)
    return jcall(obj, "getPopulationSize", jint, ())
end

function get_sample_size(obj::HypergeometricDistribution)
    return jcall(obj, "getSampleSize", jint, ())
end

function get_support_lower_bound(obj::HypergeometricDistribution)
    return jcall(obj, "getSupportLowerBound", jint, ())
end

function get_support_upper_bound(obj::HypergeometricDistribution)
    return jcall(obj, "getSupportUpperBound", jint, ())
end

function is_support_connected(obj::HypergeometricDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end

function log_probability(obj::HypergeometricDistribution, arg0::jint)
    return jcall(obj, "logProbability", jdouble, (jint,), arg0)
end

function probability(obj::HypergeometricDistribution, arg0::jint)
    return jcall(obj, "probability", jdouble, (jint,), arg0)
end

function upper_cumulative_probability(obj::HypergeometricDistribution, arg0::jint)
    return jcall(obj, "upperCumulativeProbability", jdouble, (jint,), arg0)
end


text:
# This program is licensed under the MIT License (MIT).
# Copyright (C) 2020 Oliver Schulz <oschulz@mpp.mpg.de>

orig_dir = @__DIR__
@assert isfile(joinpath(orig_dir, "legend-julia-software-tutorial.jmd"))

prj_dir = length(ARGS) >= 1 ? ARGS[1] : orig_dir
@assert isdir(prj_dir)
prj_dir = abspath(prj_dir)
@assert isfile(joinpath(prj_dir, "Project.toml")) && isfile(joinpath(prj_dir, "Manifest.toml"))

scratchdir = mktempdir(prefix="legend_jl_pkgcompile_")
tmp_dir = joinpath(scratchdir, basename(orig_dir))
cp(orig_dir, tmp_dir)
cd(tmp_dir)

import Pkg
Pkg.activate(prj_dir)
Pkg.instantiate()
Pkg.precompile()

include("make.jl")

import PackageCompiler, Libdl

custom_sysimg = joinpath(prj_dir, "JuliaSysimage." * Libdl.dlext)

PackageCompiler.create_sysimage(
    Symbol.(keys(Pkg.project().dependencies)),
    sysimage_path = custom_sysimg,
    precompile_execution_file = [
        "precompile_exec.jl",
        "legend-julia-software-tutorial.jl",
    ],
    cpu_target=PackageCompiler.default_app_cpu_target()
)

default_sysimg = abspath(Sys.BINDIR, "..", "lib", "julia", "sys." * Libdl.dlext)

import Markdown

show(Markdown.parse("
LEGEND Julia system image created.

Default Julia system image is \"$default_sysimg\", to start Julia with the LEGEND environment and system image, use

```shell
julia --project=\"$prj_dir\" --sysimage=\"$custom_sysimg\"
```

Run

```julia
julia> import IJulia
julia> IJulia.installkernel(\"LEGEND Julia\", \"--project=$prj_dir\", \"--sysimage=$custom_sysimg\")
```

to install a Jupyter Julia kernel that will use the LEGEND environment and system image.

To check which system image is active in a Julia session, use

```
julia> Base.JLOptions().image_file |> unsafe_string
```
"))

max_stars_repo_name: fhagemann/legend-julia-tutorial

max_issues_repo_name: fhagemann/legend-julia-tutorial

max_forks_repo_name: fhagemann/legend-julia-tutorial


text:
## extended methods for computing covariance and scatter matrix

# auxiliary functions

function _symmetrize!(a::DenseMatrix)
    m, n = size(a)
    m == n || error("a must be a square matrix.")
    for j = 1:n
        @inbounds for i = j+1:n
            vl = a[i,j]
            vr = a[j,i]
            a[i,j] = a[j,i] = middle(vl, vr)
        end
    end
    return a
end

function _scalevars(x::DenseMatrix, s::DenseVector, dims::Int)
    dims == 1 ? Diagonal(s) * x :
    dims == 2 ? x * Diagonal(s) :
    error("dims should be either 1 or 2.")
end

## scatter matrix

_unscaled_covzm(x::DenseMatrix, dims::Colon)   = unscaled_covzm(x)
_unscaled_covzm(x::DenseMatrix, dims::Integer) = unscaled_covzm(x, dims)

_unscaled_covzm(x::DenseMatrix, wv::AbstractWeights, dims::Colon)   =
    _symmetrize!(unscaled_covzm(x, _scalevars(x, convert(Vector, wv))))
_unscaled_covzm(x::DenseMatrix, wv::AbstractWeights, dims::Integer) =
    _symmetrize!(unscaled_covzm(x, _scalevars(x, convert(Vector, wv), dims), dims))

"""
    scattermat(X, [wv::AbstractWeights]; mean=nothing, dims=1)

Compute the scatter matrix, which is an unnormalized covariance matrix.
A weighting vector `wv` can be specified to weight
the estimate.

# Arguments
* `mean=nothing`: a known mean value. `nothing` indicates that the mean is
  unknown, and the function will compute the mean. Specifying `mean=0` indicates
  that the data are centered and hence there's no need to subtract the mean.
* `dims=1`: the dimension along which the variables are organized.
  When `dims = 1`, the variables are considered columns with observations in rows;
  when `dims = 2`, variables are in rows with observations in columns.
"""
function scattermat end

"""
    cov(X, w::AbstractWeights, vardim=1; mean=nothing,  corrected=false)

Compute the weighted covariance matrix. Similar to `var` and `std` the biased covariance
matrix (`corrected=false`) is computed by multiplying `scattermat(X, w)` by
``\\frac{1}{\\sum{w}}`` to normalize. However, the unbiased covariance matrix
(`corrected=true`) is dependent on the type of weights used:
* `AnalyticWeights`: ``\\frac{1}{\\sum w - \\sum {w^2} / \\sum w}``
* `FrequencyWeights`: ``\\frac{1}{\\sum{w} - 1}``
* `ProbabilityWeights`: ``\\frac{n}{(n - 1) \\sum w}`` where ``n`` equals `count(!iszero, w)`
* `Weights`: `ArgumentError` (bias correction not supported)
"""
cov

"""
    mean_and_cov(x, [wv::AbstractWeights,] vardim=1; corrected=false) -> (mean, cov)

Return the mean and covariance matrix as a tuple. A weighting
vector `wv` can be specified. `vardim` that designates whether
the variables are columns in the matrix (`1`) or rows (`2`).
Finally, bias correction is applied to the covariance calculation if
`corrected=true`. See [`cov`](@ref) documentation for more details.
"""
function mean_and_cov end

scattermat(x::DenseMatrix; mean=nothing, dims::Int=1) =
    _scattermatm(x, mean, dims)
_scattermatm(x::DenseMatrix, ::Nothing, dims::Int) =
    _unscaled_covzm(x .- mean(x, dims=dims), dims)
_scattermatm(x::DenseMatrix, mean, dims::Int=1) =
    _unscaled_covzm(x .- mean, dims)

scattermat(x::DenseMatrix, wv::AbstractWeights; mean=nothing, dims::Int=1) =
    _scattermatm(x, wv, mean, dims)
_scattermatm(x::DenseMatrix, wv::AbstractWeights, ::Nothing, dims::Int) =
    _unscaled_covzm(x .- mean(x, wv, dims=dims), wv, dims)
_scattermatm(x::DenseMatrix, wv::AbstractWeights, mean, dims::Int) =
    _unscaled_covzm(x .- mean, wv, dims)

## weighted cov
covm(x::DenseMatrix, mean, w::AbstractWeights, dims::Int=1;
     corrected::DepBool=nothing) =
    rmul!(scattermat(x, w, mean=mean, dims=dims), varcorrection(w, depcheck(:covm, corrected)))

cov(x::DenseMatrix, w::AbstractWeights, dims::Int=1; corrected::DepBool=nothing) =
    covm(x, mean(x, w, dims=dims), w, dims; corrected=depcheck(:cov, corrected))

function corm(x::DenseMatrix, mean, w::AbstractWeights, vardim::Int=1)
    c = covm(x, mean, w, vardim; corrected=false)
    s = stdm(x, w, mean, vardim; corrected=false)
    cov2cor!(c, s)
end

"""
    cor(X, w::AbstractWeights, dims=1)

Compute the Pearson correlation matrix of `X` along the dimension
`dims` with a weighting `w` .
"""
cor(x::DenseMatrix, w::AbstractWeights, dims::Int=1) =
    corm(x, mean(x, w, dims=dims), w, dims)

function mean_and_cov(x::DenseMatrix, dims::Int=1; corrected::Bool=true)
    m = mean(x, dims=dims)
    return m, covm(x, m, dims, corrected=corrected)
end
function mean_and_cov(x::DenseMatrix, wv::AbstractWeights, dims::Int=1;
                      corrected::DepBool=nothing)
    m = mean(x, wv, dims=dims)
    return m, cov(x, wv, dims; corrected=depcheck(:mean_and_cov, corrected))
end

"""
    cov2cor(C, s)

Compute the correlation matrix from the covariance matrix `C` and a vector of standard
deviations `s`. Use `StatsBase.cov2cor!` for an in-place version.
"""
cov2cor(C::AbstractMatrix, s::AbstractArray) = cov2cor!(copy(C), s)

"""
    cor2cov(C, s)

Compute the covariance matrix from the correlation matrix `C` and a vector of standard
deviations `s`. Use `StatsBase.cor2cov!` for an in-place version.
"""
cor2cov(C::AbstractMatrix, s::AbstractArray) = cor2cov!(copy(C), s)

"""
    cor2cov!(C, s)

Converts the correlation matrix `C` to a covariance matrix in-place using a vector of
standard deviations `s`.
"""
function cor2cov!(C::AbstractMatrix, s::AbstractArray)
    n = length(s)
    size(C) == (n, n) || throw(DimensionMismatch("inconsistent dimensions"))
    for i in CartesianIndices(size(C))
        @inbounds C[i] *= s[i[1]] * s[i[2]]
    end
    return C
end

"""
    CovarianceEstimator

Abstract type for covariance estimators.
"""
abstract type CovarianceEstimator end

"""
    cov(ce::CovarianceEstimator, x::AbstractVector; mean=nothing)

Compute a variance estimate from the observation vector `x` using the  estimator `ce`.
"""
cov(ce::CovarianceEstimator, x::AbstractVector; mean=nothing) =
    error("cov is not defined for $(typeof(ce)) and $(typeof(x))")

"""
    cov(ce::CovarianceEstimator, x::AbstractVector, y::AbstractVector)

Compute the covariance of the vectors `x` and `y` using estimator `ce`.
"""
cov(ce::CovarianceEstimator, x::AbstractVector, y::AbstractVector) =
    error("cov is not defined for $(typeof(ce)), $(typeof(x)) and $(typeof(y))")

"""
    cov(ce::CovarianceEstimator, X::AbstractMatrix, [w::AbstractWeights]; mean=nothing, dims::Int=1)

Compute the covariance matrix of the matrix `X` along dimension `dims`
using estimator `ce`. A weighting vector `w` can be specified.
The keyword argument `mean` can be:

* `nothing` (default) in which case the mean is estimated and subtracted
  from the data `X`,
* a precalculated mean in which case it is subtracted from the data `X`.
  Assuming `size(X)` is `(N,M)`, `mean` can either be:
  * when `dims=1`, an `AbstractMatrix` of size `(1,M)`,
  * when `dims=2`, an `AbstractVector` of length `N` or an `AbstractMatrix`
    of size `(N,1)`.
"""
cov(ce::CovarianceEstimator, X::AbstractMatrix; mean=nothing, dims::Int=1) =
    error("cov is not defined for $(typeof(ce)) and $(typeof(X))")

cov(ce::CovarianceEstimator, X::AbstractMatrix, w::AbstractWeights; mean=nothing, dims::Int=1) =
    error("cov is not defined for $(typeof(ce)), $(typeof(X)) and $(typeof(w))")

"""
    SimpleCovariance(;corrected::Bool=false)

Simple covariance estimator. Estimation calls `cov(x; corrected=corrected)`,
`cov(x, y; corrected=corrected)` or `cov(X, w, dims; corrected=corrected)`
where `x`, `y` are vectors, `X` is a matrix and `w` is a weighting vector.
"""
struct SimpleCovariance <: CovarianceEstimator
    corrected::Bool
    SimpleCovariance(;corrected::Bool=false) = new(corrected)
end

cov(sc::SimpleCovariance, x::AbstractVector) =
    cov(x; corrected=sc.corrected)

cov(sc::SimpleCovariance, x::AbstractVector, y::AbstractVector) =
    cov(x, y; corrected=sc.corrected)

function cov(sc::SimpleCovariance, X::AbstractMatrix; dims::Int=1, mean=nothing)
    dims ∈ (1, 2) || throw(ArgumentError("Argument dims can only be 1 or 2 (given: $dims)"))
    if mean === nothing
        return cov(X; dims=dims, corrected=sc.corrected)
    else
        return covm(X, mean, dims, corrected=sc.corrected)
    end
end

function cov(sc::SimpleCovariance, X::AbstractMatrix, w::AbstractWeights; dims::Int=1, mean=nothing)
    dims ∈ (1, 2) || throw(ArgumentError("Argument dims can only be 1 or 2 (given: $dims)"))
    if mean === nothing
        return cov(X, w, dims, corrected=sc.corrected)
    else
        return covm(X, mean, w, dims, corrected=sc.corrected)
    end
end


text:
using DataFrames
using Test

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = page_model()
    include("../src/components/SLRDamages.jl")

    slrdamages = addslrdamages(m)

    set_param!(m, :SLRDamages, :y_year, [2020, 2030, 2040, 2050, 2075, 2100, 2150, 2200, 2250, 2300])
    set_param!(m, :SLRDamages, :atl_adjustedtolerablelevelofsealevelrise, readpagedata(m, "test/validationdata/atl_adjustedtolerablelevelofsealevelrise.csv"))
    set_param!(m, :SLRDamages, :imp_actualreductionSLR, readpagedata(m, "test/validationdata/$valdir/imp_actualreductionSLR.csv"))
    set_param!(m, :SLRDamages, :y_year_0, 2015.)
    set_param!(m, :SLRDamages, :s_sealevel, readpagedata(m, "test/validationdata/$valdir/s_sealevel.csv"))
    set_param!(m, :SLRDamages, :cons_percap_consumption, readpagedata(m, "test/validationdata/$valdir/cons_percap_consumption.csv"))
    set_param!(m, :SLRDamages, :tct_per_cap_totalcostspercap, readpagedata(m, "test/validationdata/$valdir/tct_per_cap_totalcostspercap.csv"))
    set_param!(m, :SLRDamages, :act_percap_adaptationcosts, readpagedata(m, "test/validationdata/$valdir/act_percap_adaptationcosts.csv"))
    set_param!(m, :SLRDamages, :isatg_impactfxnsaturation, 28.333333333333336)

    ##running Model
    run(m)

    cons_percap_aftercosts = m[:SLRDamages, :cons_percap_aftercosts]
    cons_percap_aftercosts_compare = readpagedata(m, "test/validationdata/$valdir/cons_percap_aftercosts.csv")
    @test ones(10, 8) ≈ cons_percap_aftercosts ./ cons_percap_aftercosts_compare rtol=.001

    i_regionalimpactSLR = m[:SLRDamages, :i_regionalimpactSLR]
    i_regionalimpactSLR_compare = readpagedata(m, "test/validationdata/$valdir/i_regionalimpact_SLRise.csv")
    @test i_regionalimpactSLR ≈ i_regionalimpactSLR_compare rtol=.001

    iref_ImpactatReferenceGDPperCapSLR = m[:SLRDamages, :iref_ImpactatReferenceGDPperCapSLR]
    iref_ImpactatReferenceGDPperCapSLR_compare = readpagedata(m, "test/validationdata/$valdir/iref_ImpactatReferenceGDPperCap_sea.csv")
    @test iref_ImpactatReferenceGDPperCapSLR ≈ iref_ImpactatReferenceGDPperCapSLR_compare rtol=.0001

    rcons_per_cap_SLRRemainConsumption = m[:SLRDamages, :rcons_per_cap_SLRRemainConsumption]
    rcons_per_cap_SLRRemainConsumption_compare = readpagedata(m, "test/validationdata/$valdir/rcons_per_cap_SLRRemainConsumption.csv")
    @test rcons_per_cap_SLRRemainConsumption ≈ rcons_per_cap_SLRRemainConsumption_compare rtol=.01

    rgdp_per_cap_SLRRemainGDP = m[:SLRDamages, :rgdp_per_cap_SLRRemainGDP]
    rgdp_per_cap_SLRRemainGDP_compare = readpagedata(m, "test/validationdata/$valdir/rgdp_per_cap_SLRRemainGDP.csv")
    @test rgdp_per_cap_SLRRemainGDP ≈ rgdp_per_cap_SLRRemainGDP_compare rtol=.01

    igdp_ImpactatActualGDPperCapSLR = m[:SLRDamages, :igdp_ImpactatActualGDPperCapSLR]
    igdp_ImpactatActualGDPperCapSLR_compare = readpagedata(m, "test/validationdata/$valdir/igdp_ImpactatActualGDPperCap_sea.csv")
    @test igdp_ImpactatActualGDPperCapSLR ≈ igdp_ImpactatActualGDPperCapSLR_compare rtol=.01

    isat_ImpactinclSaturationandAdaptationSLR = m[:SLRDamages, :isat_ImpactinclSaturationandAdaptationSLR]
    isat_ImpactinclSaturationandAdaptationSLR_compare = readpagedata(m, "test/validationdata/$valdir/isat_ImpactinclSaturationandAdaptation_SLRise.csv")
    @test isat_ImpactinclSaturationandAdaptationSLR ≈ isat_ImpactinclSaturationandAdaptationSLR_compare rtol=.001
end

max_stars_repo_name: jkikstra/PAGE-VAR

max_issues_repo_name: jkikstra/PAGE-VAR

max_forks_repo_name: jkikstra/PAGE-VAR


text:
push!(LOAD_PATH, "../../src/")

using NESSie
using NESSie.Format

abstract type MeshType end
struct Nodes   <: MeshType end
struct Surface <: MeshType end
struct Volume  <: MeshType end

const formats = Dict{String, String}(
    "hmo"  =>  "hmo",
    "m"    =>  "mcsf",
    ""     =>  "msms",
    "off"  =>  "off",
    "stl"  =>  "stl"
)

const fin = Dict{String, Tuple{DataType, Function, String}}(
    "hmo"   =>  (Surface, readhmo,  "HyperMesh/.hmo"),
    "mcsf"  =>  (Volume,  readmcsf, "GAMer/.m"),
    "msms"  =>  (Surface, readmsms, "MSMS/.{face,vert}"),
    "off"   =>  (Surface, readoff,  "GAMer/.off"),
    "stl"   =>  (Surface, readstl,  "STL/.stl")
)

const fout = Dict{String, Tuple{DataType, Function, String}}(
    "nodes.json"    =>  (Nodes,   writexml3d_json, "XML3D/.json"),
    "nodes.xml"     =>  (Nodes,   writexml3d_xml,  "XML3D/.xml"),
    "surface.hmo"   =>  (Surface, writehmo,        "HyperMesh/.hmo"),
    "surface.json"  =>  (Surface, writexml3d_json, "XML3D/.json"),
    "surface.skel"  =>  (Surface, writeskel,       "SKEL/.skel"),
    "surface.stl"   =>  (Surface, writestl,        "STL/.stl"),
    "surface.vtp"   =>  (Surface, writevtk,        "VTK/.vtp"),
    "volume.skel"   =>  (Volume,  writeskel,       "SKEL/.skel"),
    "volume.vtu"    =>  (Volume,  writevtk,        "VTK/.vtu")
)

iformat = ""; ifname = ""; oformat = ""; ofname = ""

if length(ARGS) == 3
    ifname  = ARGS[1]
    oformat = ARGS[2]
    ofname  = ARGS[3]
    spl = split(split(ifname, "/")[end], ".")
    iformat = length(spl) == 1 ? formats[""] : formats[spl[end]]
elseif length(ARGS) == 4
    iformat = ARGS[1]
    ifname  = ARGS[2]
    oformat = ARGS[3]
    ofname  = ARGS[4]
else
    println("\n\e[1mNESSie.jl file format converter\e[0m")
    println("==================================\n")
    println("\e[1mUsage:\e[0m\n")
    println("\tconverter.jl [<format in>] <file in> <format out> <file out>\n")
    println("If not specified, the input format will be guessed by file extension.")
    println("Multi-file input, as provided by some file formats (e.g., MSMS), has")
    println("to be specified without file extension. For instance,")
    println("\n\tconverter.jl msms mymesh surface.json mymesh.json\n")
    println("will read the MSMS-generated surface mesh from both files mymesh.face")
    println("and mymesh.vert and convert it into a XML3D-compatible JSON format.")
    println("\n\e[1mAvailable input formats:\e[0m\n")
    for (format, (mesh, _, desc)) in fin
        println("\t\e[1m$format\e[0m\t\t$mesh\t\t$desc")
    end
    println("\n\e[1mAvailable output formats:\e[0m\n")
    for (format, (mesh, _, desc)) in fout
        println("\t\e[1m$format\e[0m\t$mesh\t\t$desc")
    end
    println()
    exit(1)
end

mtype, f = fout[oformat][1:2]
model    = fin[iformat][2](ifname)
f(ofname, mtype === Nodes ? model.nodes : model)


text:
# Salsa20 specific tests
# Copyright (C): Jaka Smrekar (vinctux) <vinctux@outlook.com>, 2016

using Krypto
using Base.Test

println("SALSA20 TESTS STARTED...")

# Test values - 1024 octets (random, 0x00, 0xFF, flip)
T = [[csrand(0o1:0o255) for i in 1:256],
     [0o0 for i in 1:256],
     [0o377 for i in 1:256],
     [i % 2 == 0 ? 0o0 : 0o1 for i in 1:256],
     Array{UInt8, 1}("This is Salsa20-Is this Salsa20?")]

IV = csrand(32)

for i in 1:length(T)
    for j in 1:length(T)
        println("Running encryption test #$(i) with key #$(j)... ")
        @time @test decrypt(Salsa20, encrypt(Salsa20, T[i], IV, T[j]), IV, T[j]) == T[i]
    end
end

println("ALL SALSA20 TESTS PASSED.")


text:
#use this to clear a dictionary
function clearcollection(c::AbstractDict)
  for i ∈ 1:length(c)
    pop!(c)
  end
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA