text:
using Distributions
using ImageFiltering
using LACosmic
using PSFModels: Gaussian
using StableRNGs
using Test

rng = StableRNG(125588)

function make_data(rng, T, N; N_sources=100, N_cosmics=100)
    imdata = fill(200.0, (N, N))

    # Add some fake sources
    for _ in 1:N_sources
        x = rand(rng, Uniform(1, N + 1))
        y = rand(rng, Uniform(1, N + 1))
        brightness = rand(rng, Uniform(1000, 30000)) / (2π * 3.5^2)
		model = Gaussian(T; x=x, y=y, fwhm=3.5, amp=brightness)
        imdata .+= model[axes(imdata)...]
	end

    # Add the poisson noise
    imdata .= rand.(rng, Poisson.(imdata))

    # Add readnoise
    imdata .+= rand(rng, Normal(0, 10), (N, N))

	clean_image = T.(imdata)

    # Add Nc fake cosmic rays
    crmask = falses((N, N))
	for i in 1:N_cosmics
    	cr_x = round(Int, rand(rng, Uniform(6, N - 5)))
    	cr_y = round(Int, rand(rng, Uniform(6, N - 5)))
    	cr_brightnesses = rand(rng, Uniform(1000, 30000))
    	imdata[cr_y, cr_x] += cr_brightnesses
    	crmask[cr_y, cr_x] = true
	end

    # Make a mask where the detected cosmic rays should be
    return (image=T.(imdata), clean_image, mask=crmask)
end

make_data(T, N; kwargs...) = make_data(rng, T, N; kwargs...)
make_data(N; kwargs...) = make_data(Float32, N; kwargs...)

@testset "LACosmic.jl" begin
    @testset "subsample" for T in [Int, Float64, Float32, BigFloat, BigInt]
        data = T.(reshape(1:25, 5, 5))
        sub = @inferred LACosmic.subsample(data)

        @test size(sub) == (10, 10)
        @test parent(sub) === data

        @inbounds for idx in CartesianIndices(data)
            idx_11 = CartesianIndex(idx.I[1] * 2 - 1, idx.I[2] * 2 - 1)
            @test data[idx] === @inferred sub[idx_11]
            sub[idx_11] = 0
            @test data[idx] == 0

            idx_12 = CartesianIndex(idx.I[1] * 2 - 1, idx.I[2] * 2)
            @test data[idx] === @inferred sub[idx_12]
            sub[idx_12] = 0
            @test data[idx] == 0

            idx_21 = CartesianIndex(idx.I[1] * 2, idx.I[2] * 2 - 1)
            @test data[idx] === @inferred sub[idx_21]
            sub[idx_21] = 0
            @test data[idx] == 0

            idx_22 = CartesianIndex(idx.I[1] * 2, idx.I[2] * 2)
            @test data[idx] === @inferred sub[idx_22]
            sub[idx_22] = 0
            @test data[idx] == 0
        end

    end

    @testset "dilation - $k" for k in [3, 5]
        mask = randn(rng, 1001, 1001) .< 0.05
        kernel = ones(k, k)
        expected = mapwindow(any, mask, (k, k), border=Fill(false))
        dilated = @inferred LACosmic.dilate!(copy(mask), k)
        @test dilated == expected
    end

    @testset "rebin" for block_size in [2, 3]
        data = randn(rng, Float32, 512, 512) .+ 100
        sub = @inferred LACosmic.subsample(data, block_size)
        binned = @inferred LACosmic.rebin(sub, block_size)
        @test binned ≈ data

        binned_sum = @inferred LACosmic.rebin(sum, sub, block_size)
        @test binned_sum ≈ block_size^2 .* data
    end

    @testset "cleaning" begin
        input = make_data(1001)
        output, mask = lacosmic(input.image; sigma_clip=6, neighbor_thresh=1, readnoise=10)
        output2, mask2 = lacosmic(output; sigma_clip=6, neighbor_thresh=1, readnoise=10)
        @test iszero(count(mask2))
    end

    @testset "LACosmic algorithm" for T in [Float32, Float64]
        input = make_data(T, 1001)
        output, mask = @inferred lacosmic(input.image; sigma_clip=6, neighbor_thresh=1, readnoise=10)
        @test eltype(output) === T
        @test mask == input.mask
        @test all(output[input.mask] .!= input.image[input.mask])
    end
end


text:
# ---
# title: 1434. Number of Ways to Wear Different Hats to Each Other
# id: problem1434
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming, Bit Manipulation
# link: <https://leetcode.com/problems/number-of-ways-to-wear-different-hats-to-each-other/description/>
# hidden: true
# ---
#
# There are `n` people and 40 types of hats labeled from 1 to 40.
#
# Given a list of list of integers `hats`, where `hats[i]` is a list of all hats
# preferred by the `i-th` person.
#
# Return the number of ways that the n people wear different hats to each other.
#
# Since the answer may be too large, return it modulo `10^9 + 7`.
#
#
#
# **Example 1:**
#
#
#
#     Input: hats = [[3,4],[4,5],[5]]
#     Output: 1
#     Explanation: There is only one way to choose hats given the conditions.
#     First person choose hat 3, Second person choose hat 4 and last one hat 5.
#
# **Example 2:**
#
#
#
#     Input: hats = [[3,5,1],[3,5]]
#     Output: 4
#     Explanation: There are 4 ways to choose hats
#     (3,5), (5,3), (1,3) and (1,5)
#
#
# **Example 3:**
#
#
#
#     Input: hats = [[1,2,3,4],[1,2,3,4],[1,2,3,4],[1,2,3,4]]
#     Output: 24
#     Explanation: Each person can choose hats labeled from 1 to 4.
#     Number of Permutations of (1,2,3,4) = 24.
#
#
# **Example 4:**
#
#
#
#     Input: hats = [[1,2,3],[2,3,5,6],[1,3,7,9],[1,8,9],[2,5,7]]
#     Output: 111
#
#
#
#
# **Constraints:**
#
#   * `n == hats.length`
#   * `1 <= n <= 10`
#   * `1 <= hats[i].length <= 40`
#   * `1 <= hats[i][j] <= 40`
#   * `hats[i]` contains a list of **unique** integers.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
# Global enable of performance counters
#
# See: Section 18.2 (Architectural Performance Monitoring) of Volume 3 of the Software Developer's Guide, as well as in Chapter 35
const IA32_PERF_GLOBAL_CTRL_MSR = IndexZero(0x38f)
const IA32_FIXED_CTR_CTRL_MSR = IndexZero(0x38D)

# MSR for programming programable performance counters.
# Indexed by counter number.
const EVENT_SELECT_MSRS = IndexZero.((
    0x186,
    0x187,
    0x188,
    0x189,
    0x18a,
    0x18b,
    0x18c,
    0x18d,
))

# MSRs for reading the contents of performance counters.
# Indexed by counter number.
const PMC_MSRS = IndexZero.((
    0xc1,
    0xc2,
    0xc3,
    0xc4,
    0xc6,
    0xc7,
    0xc9,
    0xc9,
))

"""
Return the number of programable performance counters on your CPU.
"""
function numcounters()
    # Read from the global control register - use CPU 0 because it always exists
    val = readmsr(IndexZero(0), IA32_PERF_GLOBAL_CTRL_MSR)

    # Clear bits 32, 33, and 34 since these correspond to the fixed function registers.
    val = clearbits(val, (32, 33, 34))

    # The number of set bits is the number of programmable counters.
    return count_ones(val)
end

# TODO: This is a hack at the moment.
#
# Read /dev/cpu.
# There is one entry per cpu, plus a "microcode" directory.
# We subtract 1 for the "microcode" directory to get the number of CPUs.
numcpus() = length(readdir("/dev/cpu")) - 1

# Write 1's to the performance counter locations.
function enablecounters(cpu)
    config = 0x0000000700000000 | (mask(numcounters())-1)
    writemsr(cpu, IA32_PERF_GLOBAL_CTRL_MSR, config)
end
disablecounters(cpu) = writemsr(cpu, IA32_PERF_GLOBAL_CTRL_MSR, zero(UInt64))


text:
"""
```
struct XMMAApprox <: AbstractMMAApprox
    mma::MMAApprox
    approxfg::AbstractVector
    a::AbstractVector
    c::AbstractVector
    d::AbstractVector
    y::AbstractVector
    z::Base.RefValue{<:Real}
    out::AbstractVector
end
```

This is similar to [`MMAApprox`](@ref) but it uses the extended formulation given in the [2002 MMA paper](https://epubs.siam.org/doi/abs/10.1137/S1052623499362822). `y` and `z` are additional variables added and `a`, `c` and `d` are coefficients, all of which are explained in the paper. `mma` is the underlying MMA approximation of the objective and constraints. Calling `mma` returns a vector whose first element is the objective approximation and the rest are the constraint violation approximations.

If the original model is:
 - Variables: x
 - Objective: f(x)
 - Constraint i: g(x) <= 0
 - Bounds: x ∈ X
the transformed model is:
 - Variables: x, y, z
 - Objective: f(x) + a[1] * z + Σ_i c[i] * y[i] + 1/2 * d[i] * y[i]^2
 - Constraint i: g_i(x) - a[i + 1] * z - y[i] <= 0
 - Bounds: x ∈ X, y >= 0, z >= 0
"""
@params struct XMMAApprox <: AbstractMMAApprox
    mma::MMAApprox
    approxfg::AbstractVector
    a::AbstractVector
    c::AbstractVector
    d::AbstractVector
    y::AbstractVector
    z::Base.RefValue{<:Real}
    out::AbstractVector
end

"""
```
XMMAApprox(
    mma::MMAApprox,
    a = fill(1.0, size(mma.p, 1)),
    c = fill(1.1, size(mma.p, 1) - 1),
    d = fill(1.0, size(mma.p, 1) - 1),
    kwargs...,
)
```

Constructs an instance of the extended MMA approximation `XMMAApprox` using the MMA approximation `mma` and coefficients `a`, `b` and `c`. See [`XMMAApprox`](@ref) for an explanation.
"""
function XMMAApprox(
    mma::MMAApprox,
    a = fill(1.0, size(mma.p, 1)),
    c = fill(1.1, size(mma.p, 1) - 1),
    d = fill(1.0, size(mma.p, 1) - 1),
    kwargs...,
)
    x0 = getxk(mma)
    T = eltype(x0)
    approxfg = mma(x0)
    z0 = sqrt(eps(T))
    y0 = max.(approxfg[2:end] .- a[2:end] .* z0, z0)
    out = similar(mma.out)
    return XMMAApprox(mma, approxfg, a, c, d, y0, Ref(z0), out)
end

"""
    (approx::XMMAApprox)(xyz::AbstractVector)

Returns the values of the approximate objective and constraints. `xyz` is `x`, `y` and `z` concatenated in a single vector.
"""
function (approx::XMMAApprox)(xyz::AbstractVector)
    @unpack mma, a, c, d = approx
    T = eltype(xyz)
    z = xyz[end]
    y = xyz[end - length(c) : end - 1]
    x = xyz[begin : end - length(c) - 1]
    approxfg = mma(x)
    ε = 1e-4
    out = approxfg .+ [a[1] * z + ε * z^2 + sum(c .* y .+ d .* y.^2 ./ 2); .-a[2:end] .* z .- y]
    saveout!(approx, out)
    return out
end

"""
    saveout!(approx::XMMAApprox, out)

Saves `out` in the cache `approx.out`.
"""
function saveout!(approx::XMMAApprox, out)
    approx.out .= out
    return approx
end
ChainRulesCore.@non_differentiable saveout!(approx::XMMAApprox, out)

getmmaapprox(approx::XMMAApprox) = approx.mma

getasymptotes(approx::XMMAApprox) = getasymptotes(approx.mma)

getapproxfg(approx::XMMAApprox) = approx.mma.out

getσ(approx::XMMAApprox) = getσ(approx.mma)
function setσ!(approx::XMMAApprox, σ::AbstractVector)
    setσ!(approx.mma, σ)
    return approx
end

getρ(approx::XMMAApprox) = getρ(approx.mma)
function setρ!(approx::XMMAApprox, ρ::AbstractVector)
    setρ!(approx.mma, ρ)
    return approx
end

getfk(approx::XMMAApprox) = getfk(approx.mma)
function setfk!(approx::XMMAApprox, f::AbstractVector)
    setfk!(approx.mma, f)
    return approx
end

getxk(approx::XMMAApprox) = getxk(approx.mma)
function setxk!(approx::XMMAApprox, x::AbstractVector)
    setxk!(approx.mma, x)
    return approx
end

get∇fk(approx::XMMAApprox) = get∇fk(approx.mma)
function set∇fk!(approx::XMMAApprox, ∇f::AbstractVecOrMat)
    set∇fk!(approx.mma, ∇f)
    return approx
end

getparent(approx::XMMAApprox) = getparent(approx.mma)
getdim(approx::XMMAApprox) = getdim(approx.mma)

function updateapprox!(
    approx::XMMAApprox,
    x::AbstractVector,
)
    updateapprox!(approx.mma, x)
    return approx
end
function updateapprox!(
    approx::XMMAApprox,
    x::AbstractVector,
    f::AbstractVector,
    ∇f::AbstractMatrix,
)
    updateapprox!(approx.mma, x, f, ∇f)
    return approx
end
function updateapprox!(approx::XMMAApprox)
    updateapprox!(approx.mma)
    return approx
end


text:
using myJuliaUtils
using Test
using LinearAlgebra
using DelimitedFiles
using Statistics

# Testing both lag functions
@testset "testing lag functions" begin
    # if vector
    T = 100
    P = rand(1:T - 1)
    a = collect(1.0:1.0:T)
    b = lag0(a, P)
    @test b[1:P] == zeros(P)
    @test T - P == b[end]

    # if matrix
    a = reshape(a, T, 1)
    P = round(Int, 1 + rand() * T - 1)
    b = lag0(a, P)
    @test b[1:P,1] == zeros(P)
    @test T - P == b[end,1]

    a = [ones(T, 1) a]
    P = round(Int, 1 + rand() * T - 1)
    b = lag0(a, P)
    @test b[1:P,2] == zeros(P)
    @test T - P == b[end,2]
end

@testset "testing function preparexy()" begin
    T = 100
    P = rand(1:T - 1)
    a = randn(T)
    cnst = rand() > 0.5
    Y, X = preparexy(a, P, cnst)
    # for univariate
    i = rand(1:T - P)
    z = reverse(a[i:i + P - 1])
    if cnst
        push!(z, 1.)
    end
    @test X[i,:] == z
    @test Y[i] == a[i + P]

    # for multivariate
    N = rand(1:20)
    a = randn(T, N)
    Y, X = preparexy(a, P, cnst)
    nT, nX = size(X)

    i = rand(1:T - P)
    z = Vector(vec(reverse(a[i:i + P - 1,:], dims = 1)'))
    if cnst
        push!(z, 1.)
    end
    @test nT == T - P
    @test nX == P * N + cnst
    @test X[i,:] == z
    @test Y[i,:] == a[i + P,:]
end

# testing companion form function
@testset "testing companion form functions" begin
    betta = [2 3;4 5]
    n = 2
    p = 1
    c = false
    @test betta == companionf(betta, n, p, c)

    bettaC = [betta ones(n, 1)]
    c = true
    @test betta == companionf(bettaC, n, p, c)

    c = false
    p = 2
    betta = [2 3 4 5;6 7 8 9]
    betta_comp = [betta;1 0 0 0; 0 1 0 0]
    @test betta_comp == companionf(betta, n, p, c)

    bettaC = [betta ones(n, 1)]
    c = true
    @test betta_comp == companionf(bettaC, n, p, c)

    betta = [1.]
    n = 1
    p = 1
    c = false
    @test reshape(betta, 1, 1) == companionf(betta, n, p, c)

    n = 2
    p = 1
    c = false
    @test_throws ArgumentError companionf(betta, n, p, c)

    n = 1
    p = 2
    c = false
    @test_throws ArgumentError companionf(betta, n, p, c)

    n = 1
    p = 1
    c = true
    @test_throws ArgumentError companionf(betta, n, p, c)
end

# Testing stability function
@testset "testing stability functions" begin
    N = round(Int, 1 + rand() * 100)
# A = randn(N, N)
# A = A'A
    A = genPDmat(N)
    S, Q = eigen(A)
    S[end] = 1.0
    A = Q * Diagonal(S) * Q'
    @test stabcheck(A, N, 1, false)
    @test stabcheckC(A, N, 1, false)

    a = findall(x->x >= 1.0, S)
    b = length(a)
    S[a] = ones(b) - rand(b)
    A = Q * Diagonal(S) * Q'
    @test stabcheck(A, N, 1, false) == false
    @test stabcheckC(A, N, 1, false) == false

    A = [
    0.516090886728283	0.0724233649560897	0.0816329971852510	0.143587089009963	-0.00458107252041236
    0.0724233649560897	0.844836851975328	0.176990845416178	-0.0777482352998716	-0.0921908525982972
    0.0816329971852510	0.176990845416178	0.671194897897464	-0.00975809606339381	-0.0331333053292128
    0.143587089009963	-0.0777482352998717	-0.00975809606339381	0.797277510951764	-0.218043826201017
    -0.00458107252041236	-0.0921908525982972	-0.0331333053292128	-0.218043826201017	0.554361340591299
];

    @test stabcheck(A, size(A, 1), 1, false)
    @test stabcheckC(A, size(A, 1), 1, false)

    A = [
    0.500514390972542	-0.0693316507737557	0.103250184760922	0.157734480886660	0.0498798880488583
    -0.0693316507737557	0.586863186538567	0.00546111896077394	-0.0714652123701029	0.153497763972517
    0.103250184760922	0.00546111896077393	0.487346580847464	-0.0744321935787715	-0.0521727141915581
    0.157734480886660	-0.0714652123701029	-0.0744321935787715	0.566439153014101	0.0799567015086364
    0.0498798880488583	0.153497763972517	-0.0521727141915581	0.0799567015086364	0.603872518537494
];
    @test stabcheck(A, size(A, 1), 1, false) == false
    @test stabcheckC(A, size(A, 1), 1, false) == false

# AR(2) stability check
    Nvar = 1;
    LAGS = 2;
    c = false;
    reps = 50000
    BETTAt = ones(LAGS, reps);
    for i = 1:reps
        chk = 1
        while chk == true
            global BETTAt[:,i] = randn(LAGS, 1)
            chk = stabcheck(BETTAt[:,i]', Nvar, LAGS, c)
        end
    end

    @test all([BETTAt[2,:] .< BETTAt[1,:] .+ 1. BETTAt[2,:] .< -BETTAt[1,:] .+ 1. BETTAt[2,:] .> -1. ])

# some other small tests

    @test stabcheck(1., 1, 1, false) == true
    @test stabcheck([1.], 1, 1, false) == true
    @test stabcheck(ones(1, 1), 1, 1, false) == true
    @test stabcheck([1. randn()], 1, 1, true) == true
    @test stabcheckC(1., 1, 1, false) == true
    @test stabcheckC([1.], 1, 1, false) == true
    @test stabcheckC([1. randn()], 1, 1, true) == true
    @test stabcheck(1, 1, 1, false) == true
    @test stabcheck([1 randn()], 1, 1, true) == true
    @test stabcheckC(1, 1, 1, false) == true
    @test stabcheckC([1 randn()], 1, 1, true) == true
end

# test for normal pdf
x = [-2.,-1.,0.,1.,2.]
y1 = [0.0539909665131881, 0.241970724519143, 0.398942280401433, 0.241970724519143, 0.0539909665131881]

mu = 2.
sigma = 1.
y2 = [0.000133830225764885, 0.00443184841193801, 0.0539909665131881, 0.241970724519143, 0.398942280401433]

@testset "testing normal pdf function" begin
    @test normpdf.(x) ≈ y1
    @test normpdf.(x, mu, sigma) ≈ y2
end

# testing ols function
path = dirname(@__FILE__)
# path = "/Users/zymantas/TVPVARPkg/test/"

data = readdlm("$path/TestOls/data2test.txt", ',');
B2check = readdlm("$path/TestOls/B2check.txt", ',');
stats2check = readdlm("$path/TestOls/stats2check.txt", ',');
covb2check = readdlm("$path/TestOls/covb2check.txt", ',');

y = data[:,1]
x = [ones(length(y)) data[:,2:end]]
r = ols1(y, x)

@testset "testing OLS function" begin
    @test r.bhat ≈ B2check
    @test r.sigbhat ≈ covb2check
    @test r.R2[1] ≈ stats2check[1]
    @test r.sig2hat[1] ≈ stats2check[2]
end

# testing quantileArr function
a = randn(rand(1:20), rand(1:20), rand(1:20))

@testset "testing quantileArr function" begin
    @test quantileArr(a, [0.5], 1) ≈ median(a, dims = 1)
    @test quantileArr(a, [0.5], 2) ≈ median(a, dims = 2)
    @test quantileArr(a, [0.5], 3) ≈ median(a, dims = 3)
end

@testset "testing nan mean" begin
    @test nanmean([1., 2., NaN]) == 1.5
    @test nanmean([1. 2.; NaN 3.]) == 2.0
    @test nanmean([1. 2.; NaN 3.], 2) ≈ [1.5; 3. ]
    @test nanmean([1. 2.; NaN 3.], 1) ≈ [1.0 2.5 ]
    a = [1. 2.; NaN 3.]
    a = cat(a, [1. NaN; 7. 3.], dims = 3)
    @test nanmean(a, 3) ≈ [1. 2.; 7. 3.]

    a = randn(4, 3, 5)
    @test mean(a;dims = 1) ≈ nanmean(a, 1)
    @test mean(a;dims = 2) ≈ nanmean(a, 2)
    @test mean(a;dims = 3) ≈ nanmean(a, 3)

end

@testset "testing inbetween" begin
    b = randn(2, 2, 2)
    b[:,:,1] = [1 2; 3 4]
    b[:,:,2] = b[:,:,1] .+ 2
    a = [3 5; 9 5]
    @test inbetween(a, b) == [true false; false true]

    T = rand(1:1000)
    n = rand(1:1000)

    a = randn(T, n)
    b = randn(T, n, 2)
    int = abs(randn())
    b[:,:,2] = b[:,:,1] .+ int
    a = b[:,:,1] .+ int ./ 2

    @test all(inbetween(a, b))

end

@testset "testing sumsqr" begin
    a = randn(rand(1:100))
    @test sumsqr(a) ≈ a'a
    @test sumsqr(a) ≈ sum(a.^2)

    a = randn(rand(1:100), rand(1:100))
    @test sumsqr(a) ≈ sum(a.^2)
end

@testset "testing getmultdiag!" begin
    n = rand(1:100)
    A = rand() .+ randn(n, n)
    B = rand() .+ randn(n, n)
    v = similar(A, n)

    getmultdiag!(v, A, B)

    @test v ≈ diag(A * B)
end

@testset "testing transf1" begin
    a = transf1To(0.2, 0., 1.)
    @test 0.2 ≈ transf1Back(a, 0., 1.)

    n = rand(1:100)
    a = randn(n)
    aup = a .+ abs.(randn(n))
    adown = a .- abs.(randn(n))

    b = transf1To.(a, adown, aup)
    @test a ≈ transf1Back.(b, adown, aup)
end

@testset "testing regMat2PD!" begin

    for i = 1:100
        a = genPDmat()
        b = copy(a)
        regMat2PD!(a)

        @test a ≈ b
        @test isapprox(a,b; atol=1e-10)
    end

end

max_stars_repo_name: zymbuzz/myJuliaUtils

max_issues_repo_name: zymbuzz/myJuliaUtils

max_forks_repo_name: zymbuzz/myJuliaUtils


text:
# syntax: proto2
using Compat
using ProtoBuf
import ProtoBuf.meta

struct __enum_EventType <: ProtoEnum
    EVENT_CREATE::Int32
    EVENT_CLOSE::Int32
    EVENT_APPEND::Int32
    EVENT_RENAME::Int32
    EVENT_METADATA::Int32
    EVENT_UNLINK::Int32
    __enum_EventType() = new(0,1,2,3,4,5)
end #struct __enum_EventType
const EventType = __enum_EventType()

struct __enum_INodeType <: ProtoEnum
    I_TYPE_FILE::Int32
    I_TYPE_DIRECTORY::Int32
    I_TYPE_SYMLINK::Int32
    __enum_INodeType() = new(0,1,2)
end #struct __enum_INodeType
const INodeType = __enum_INodeType()

struct __enum_MetadataUpdateType <: ProtoEnum
    META_TYPE_TIMES::Int32
    META_TYPE_REPLICATION::Int32
    META_TYPE_OWNER::Int32
    META_TYPE_PERMS::Int32
    META_TYPE_ACLS::Int32
    META_TYPE_XATTRS::Int32
    __enum_MetadataUpdateType() = new(0,1,2,3,4,5)
end #struct __enum_MetadataUpdateType
const MetadataUpdateType = __enum_MetadataUpdateType()

mutable struct EventProto <: ProtoType
    _type::Int32
    contents::Array{UInt8,1}
    EventProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EventProto
const __req_EventProto = Symbol[:_type,:contents]
meta(t::Type{EventProto}) = meta(t, __req_EventProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct EventBatchProto <: ProtoType
    txid::Int64
    events::Base.Vector{EventProto}
    EventBatchProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EventBatchProto
const __req_EventBatchProto = Symbol[:txid]
meta(t::Type{EventBatchProto}) = meta(t, __req_EventBatchProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct CreateEventProto <: ProtoType
    _type::Int32
    path::AbstractString
    ctime::Int64
    ownerName::AbstractString
    groupName::AbstractString
    perms::FsPermissionProto
    replication::Int32
    symlinkTarget::AbstractString
    overwrite::Bool
    defaultBlockSize::Int64
    CreateEventProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CreateEventProto
const __req_CreateEventProto = Symbol[:_type,:path,:ctime,:ownerName,:groupName,:perms]
const __val_CreateEventProto = Dict(:defaultBlockSize => 0)
meta(t::Type{CreateEventProto}) = meta(t, __req_CreateEventProto, ProtoBuf.DEF_FNUM, __val_CreateEventProto, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct CloseEventProto <: ProtoType
    path::AbstractString
    fileSize::Int64
    timestamp::Int64
    CloseEventProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct CloseEventProto
const __req_CloseEventProto = Symbol[:path,:fileSize,:timestamp]
meta(t::Type{CloseEventProto}) = meta(t, __req_CloseEventProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct AppendEventProto <: ProtoType
    path::AbstractString
    newBlock::Bool
    AppendEventProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct AppendEventProto
const __req_AppendEventProto = Symbol[:path]
const __val_AppendEventProto = Dict(:newBlock => false)
meta(t::Type{AppendEventProto}) = meta(t, __req_AppendEventProto, ProtoBuf.DEF_FNUM, __val_AppendEventProto, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct RenameEventProto <: ProtoType
    srcPath::AbstractString
    destPath::AbstractString
    timestamp::Int64
    RenameEventProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct RenameEventProto
const __req_RenameEventProto = Symbol[:srcPath,:destPath,:timestamp]
meta(t::Type{RenameEventProto}) = meta(t, __req_RenameEventProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct MetadataUpdateEventProto <: ProtoType
    path::AbstractString
    _type::Int32
    mtime::Int64
    atime::Int64
    replication::Int32
    ownerName::AbstractString
    groupName::AbstractString
    perms::FsPermissionProto
    acls::Base.Vector{AclEntryProto}
    xAttrs::Base.Vector{XAttrProto}
    xAttrsRemoved::Bool
    MetadataUpdateEventProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct MetadataUpdateEventProto
const __req_MetadataUpdateEventProto = Symbol[:path,:_type]
meta(t::Type{MetadataUpdateEventProto}) = meta(t, __req_MetadataUpdateEventProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct UnlinkEventProto <: ProtoType
    path::AbstractString
    timestamp::Int64
    UnlinkEventProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct UnlinkEventProto
const __req_UnlinkEventProto = Symbol[:path,:timestamp]
meta(t::Type{UnlinkEventProto}) = meta(t, __req_UnlinkEventProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

mutable struct EventsListProto <: ProtoType
    events::Base.Vector{EventProto}
    firstTxid::Int64
    lastTxid::Int64
    syncTxid::Int64
    batch::Base.Vector{EventBatchProto}
    EventsListProto(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct EventsListProto
const __req_EventsListProto = Symbol[:firstTxid,:lastTxid,:syncTxid]
meta(t::Type{EventsListProto}) = meta(t, __req_EventsListProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES, ProtoBuf.DEF_FIELD_TYPES)

export EventType, INodeType, MetadataUpdateType, EventProto, EventBatchProto, CreateEventProto, CloseEventProto, AppendEventProto, RenameEventProto, MetadataUpdateEventProto, UnlinkEventProto, EventsListProto


text:
# plot the distribution at one single point x0
using PyPlot
function plot_distribution_single(distribution_z0, x0, z0_true)
    pdf0 = distribution_z0.pdf
    mean = distribution_z0.mean
    # median = distribution_z0.median
    # CI = distribution_z0.CI
    # CI_left = CI[2]
    # CI_right = CI[3]

    # determine the range of plotting
    kmax = 20
    tol = 1e-3
    b0 = 5.0
    b, iter = zero_finding(pdf0, kmax, tol, b0)

    # plot pdf
    z_grid = range(1e-5, stop = b, step = 0.01)
    p_grid = pdf0.(z_grid);
    PyPlot.plot(z_grid, p_grid, label = "probability density function")
    # plot mean median and true value
    PyPlot.vlines(mean, 0, pdf0(mean), label = "mean", colors = "k")
    # PyPlot.vlines(median, 0, pdf0(median), label = "median",  colors = "b")
    PyPlot.vlines(z0_true, 0, pdf0(z0_true), label = "true value",  colors = "r")
    # plot CI
    # CI_x_range = range(CI_left, stop = CI_right, step = 0.01)
    # CI_y_range = pdf0.(CI_x_range)
    # PyPlot.fill_between(CI_x_range, 0, CI_y_range, alpha = 0.3, label = "95% confidence interval")
    PyPlot.legend()
    PyPlot.grid()
    PyPlot.title("Predicted Distribution at x = $x0")
#     PyPlot.savefig("Figure/Prediction")
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
import Theta: accola_chars_r, accola_chars_sr

@testset "Accola" begin
    @test Set(accola_chars_r([[1,1],[0,1]], [[1,0],[0,0]], [[0,1],[0,0]], [[0,0],[1,0]])) == Set([[[1,1],[0,1]], [[0,1],[0,1]], [[1,0],[0,1]], [[1,1],[1,1]], [[0,0],[0,1]], [[0,1],[1,1]], [[1,0],[1,1]], [[0,0],[1,1]]])
    @test Set(accola_chars_sr([[[1,0,0,0,0],[0,0,0,0,0]], [[0,1,0,0,0],[1,0,0,0,0]], [[0,0,1,0,0],[1,1,0,0,0]], [[0,0,0,1,0],[1,1,1,0,0]], [[0,0,0,0,1],[1,1,1,1,0]], [[1,0,0,0,0],[0,1,1,1,1]], [[0,1,0,0,0],[0,0,1,1,1]], [[0,0,1,0,0],[0,0,0,1,1]], [[0,0,0,1,0],[0,0,0,0,1]], [[1,1,1,1,0],[0,1,0,1,0]], [[1,1,1,1,0],[1,0,1,0,1]]], [[1,2,3,4], [5,6,7,8], [5,6,9,10], [5,7,9,11]])[1]) ==  Set([[[1,0,0,0,0], [0,0,0,0,0]], [[0,1,1,0,1], [1,0,1,0,1]], [[1,1,1,0,1], [1,1,0,1,0]], [[0,0,1,0,1], [0,1,1,0,1]], [[0,0,0,0,0], [0,1,1,1,1]], [[1,1,0,0,0], [1,1,0,0,0]], [[0,1,0,0,0], [1,0,1,1,1]], [[1,0,1,0,1], [0,0,0,1,0]]])
    @test Set(accola_chars()[1][1]) == Set([[[1,0,0,0,0], [0,0,0,0,0]], [[0,1,1,0,1], [1,0,1,0,1]], [[1,1,1,0,1], [1,1,0,1,0]], [[0,0,1,0,1], [0,1,1,0,1]], [[0,0,0,0,0], [0,1,1,1,1]], [[1,1,0,0,0], [1,1,0,0,0]], [[0,1,0,0,0], [1,0,1,1,1]], [[1,0,1,0,1], [0,0,0,1,0]]])
    @test begin
        τ = [0.17012929+3.36021399im 0.42793573+1.65263414im 0.55503983+2.1393021im 0.49182605+2.17345553im 0.93865305+1.81009273im; 0.42793573+1.65263414im 0.87402978+1.33203047im 0.53489381+0.68238397im 0.57601957+1.20509707im 0.50028888+1.11605622im; 0.55503983+2.1393021im 0.53489381+0.68238397im 0.21882206+2.28319153im 0.74421992+1.51352273im 0.41855281+0.87400937im; 0.49182605+2.17345553im 0.57601957+1.20509707im 0.74421992+1.51352273im 0.51070570+1.62305859im 0.52380319+1.15073893im; 0.93865305+1.81009273im 0.50028888+1.11605622im 0.41855281+0.87400937im 0.52380319+1.15073893im 0.79948032+1.38925759im];
        accola(τ) == accola(RiemannMatrix(τ))
        end
end


text:
function savebin(filename,a)
    fileID = open(filename,"w");
    write(fileID,Float64.(a[:]));
    close(fileID);
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
using Test
using Stochatto, MIDI, Statistics, Distributions

@testset "Model" begin

    @testset "Constant" begin
        @test length(NOTE_RANGE) == 97
        @test Int(NOTE_RANGE[1].pitch) == 12
    end

    @testset "Beat" begin
        @test QUARTER * SIXTEENTH == SIXTEENTH
        @test WHOLE == QUARTER * 4
        @test 3*TRIPLET == 2*EIGHTH
        @test 16*THIRTY_SECOND == HALF
    end

    @testset "Key" begin
        cminor = Key(NOTE_RANGE[1], MINOR)
        @test Int(cminor.root.pitch) % 12 == 0
        @test length(cminor.scale) == 7
        show(stdout, cminor)
    end

    @testset "Interval" begin
        @test NOTE_RANGE[1] + ROOT == NOTE_RANGE[1]
        @test pitch_to_name((NOTE_RANGE[1] + OCTAVE * 4).pitch) == "C4"
        @test Note(name_to_pitch("C6"), 127, 0, 0) + PERFECT_FOURTH == Note(name_to_pitch("F6"), 127, 0, 0)
        @test abs(-OCTAVE) == OCTAVE
    end

    @testset "Chord" begin
        note = NOTE_RANGE[1] + OCTAVE*6
        for chord in [MAJOR_CHORD, MINOR_CHORD, POWER_CHORD]
            @test (note .+ chord.notes)[1] == note
        end
    end

    @testset "Scale" begin
        @test length(intersect(String.(keys(SCALES)), String.(keys(MODES)))) == 0
        @test length(unique(values(SCALES))) == length(SCALES)
        @test length(unique(values(MODES))) == length(MODES)
        @test AEOLIAN == MINOR
        @test IONIAN == MAJOR
    end

end

@testset "Compute" begin

    @testset "Utility" begin
        N = 1000
        R = 0.01:0.01:0.99
        X = zeros((N, length(R)))
        for (j, rho) in enumerate(R)
            X[:,j] = acf_series(N, rho)
        end
        A = [cor(X[1:N-1,j], X[2:N,j]) for j in 1:size(X,2)]
        @test cor(A, R) >= 0.99
    end

    @testset "Operation" begin
        note = Note(60, 0, 0, 0)
        major_scale = note .+ Interval.(MAJOR)
        @test major_scale[1] == note
        @test note - MAJOR_CHORD.notes[1] == note
        @test (note - (note - Interval(3))) == Interval(3)
    end

    @testset "Engine" begin
        key = Key(NOTE_RANGE[rand(1:length(NOTE_RANGE))], MINOR)
        note_generator = Normal(0, 9)
        beat_generator = Poisson(2)
        engine = Engine(key, note_generator, beat_generator)
        n = 32
        initial = key.root + OCTAVE*5
        precision = EIGHTH
        notes = generate(engine, n, initial, precision)
        @test length(notes) == n
        @test length(unique([note.pitch for note in notes])) > 1
        @test length(unique([note.duration for note in notes])) > 1
        for i in 2:n
            @test notes[i].position >= notes[i-1].position + notes[i-1].duration
        end
    end

end


text:
using BenchmarkTools
using FiSUDiPPuS

options = joinpath(@__DIR__, "../data/default.jl")
include(options)

x,y = FiSUDiPPuS.get_data(settings) # settings is a dict in the options file

## Convert to correct types
start = settings[:datatype].(settings[:start])
bleach_weight = settings[:datatype](settings[:bleach_weight])

@btime FiSUDiPPuS.model($x[:,:], $start, $settings[:N], $bleach_weight);

x = 3000.0
A = [1.0, -1.0]
Γ = [8.0, 8.0]
ω = [2900.0, 2910.0]
@btime FiSUDiPPuS.sfspec($x, $A, $ω, $Γ)

# 496.331 μs (4505 allocations: 510.06 KiB)
# 336.573 μs (2259 allocations: 194.02 KiB)
# 296.531 μs (22 allocations: 89.84 KiB)
# 298.902 μs (19 allocations: 89.52 KiB)
# 196.338 μs (20 allocations: 89.63 KiB)

# 573.098 ns (11 allocations: 1.19 KiB)
# 447.626 ns (9 allocations: 928 bytes)


text:
using Documenter, SingleFloats

makedocs(
    modules = [SingleFloats],
    sitename = "SingleFloats.jl",
    authors = "Jeffrey Sarnoff",
    pages = Any[
        "Overview" => "index.md",
        "Examples" => "Examples.md",
    ]
)

deploydocs(
    repo = "github.com/JeffreySarnoff/SingleFloats.jl.git",
    target = "build"
)


text:
struct IdMap
    R::AbstractAlgebra.MPolyRing
    map::Dict{Symbol,AbstractAlgebra.MPolyElem}

    function IdMap(R::AbstractAlgebra.MPolyRing)
        vars = gens(R)
        new(R, Dict(zip(map(Symbol ∘ string, vars), vars)))
    end
end

function(m::IdMap)(s::Symbol; force=false)
    !haskey(m.map, s) && force && return zero(m.R)
    m.map[s]
end

(m::IdMap)(s::String; force=false) = m(Symbol(s); force=force)

function(m::IdMap)(p::PolyElem; force=false)
    v = m(var(parent(p)); force=force)
    res = zero(m.R)
    for i in 0:length(p)-1
        c = AbstractAlgebra.coeff(p, i)
        res += m(c; force=force) * v^(i)
    end
    res
end

function(m::IdMap)(p::MPolyElem; force=false)
    vs = [m(s, force=force) for s in AbstractAlgebra.symbols(parent(p))]
    res = zero(m.R)
    for (cc, mm) in zip(AbstractAlgebra.coeffs(p), AbstractAlgebra.monomials(p))
        res += m(cc; force=force) * evaluate(mm, vs)
    end
    res
end

function(m::IdMap)(f::FracElem; force=false)
    m(numerator(f); force=force) // m(denominator(f); force=force)
end

function(m::IdMap)(f::FieldElem; force)
    f
end

function(m::IdMap)(f::RingElem; force)
    f
end

struct IndexMap
    R::AbstractAlgebra.MPolyRing

    function IndexMap(R::AbstractAlgebra.MPolyRing)
        new(R)
    end
end

function(m::IndexMap)(p::MPolyElem; force=false)
    R = parent(p)
    x = length(AbstractAlgebra.symbols(R))
    y = length(AbstractAlgebra.symbols(m.R))

    !force && x != y && error("Number of arguments have to match")
    args = gens(m.R)
    if x > y
        append!(args, zeros(m.R, y-x))
    elseif y > x
        args = args[1:x]
    end
    evaluate(p, args)
end

function(m::Union{IndexMap,IdMap})(I::sideal)
    iszero(I) && return Singular.Ideal(m.R)
    Singular.Ideal(m.R, [m(p; force=true) for p in gens(I)])
end

function imap(I::sideal, R::MPolyRing)
    IdMap(R)(I)
end

function fetch(I::sideal, R::MPolyRing)
    IndexMap(R)(I)
end


text:
## ------------------------------------------------------------------
# References

function reflist(refs::Vector{RBRef})
    ref = RBRefList()
    push!(ref, refs...)
    return ref
end

## ------------------------------------------------------------------
# Accessors
references(r::RBRefList) = getproperty!(() -> Vector{RBRef}(), r, :refs)

## ------------------------------------------------------------------
# Meta

get_sec(q::RBRefList) = get_parent(q)
set_sec!(q::RBRefList, sec::RBSection) = set_parent!(q, sec)

## ------------------------------------------------------------------
# Base
Base.length(r::RBRefList) = length(references(r))
Base.getindex(r::RBRefList, idx) = getindex(references(r), idx)
Base.lastindex(r::RBRefList) = lastindex(references(r))
Base.firstindex(r::RBRefList) = firstindex(references(r))
Base.iterate(r::RBRefList) = iterate(references(r))
Base.iterate(r::RBRefList, state) = iterate(references(r), state)
Base.push!(r::RBRefList, ref::RBRef, refs::RBRef...) = push!(references(r), ref, refs...)

## ------------------------------------------------------------------
# show
function Base.show(io::IO, l::RBRefList)
    nrefs = length(l)
    println(io, "RBRefList(\"", get_label(l), "\") with ", nrefs, " reference(s)")
    # data
    if nrefs > 0
        print(io, "\nreference(s):")
        _show_data_preview(io, l) do ref
            # string("[", get_doi(ref), "] ", get_title(sec))
            string(ref)
        end
    end
    print(io, "\n", _preview(io, "-"^70))
end


text:
using HDF5,JLD,KUnet
isdefined(:xtrn) || (@date @load "zn11oparse.jld")
for y in (:ytrn, :ydev, :ytst); @eval $y=full($y); end # This speeds up the accuracy fn
KUnet.gpu(false) # we don't have full gpu impl for perceptron yet

# xtrn1=xtrn[:,1:100000]
# ytrn1=ytrn[:,1:100000]

for seed=1:10
    srand(seed)
    net = Layer[Perceptron(size(ytrn,1))]
    for epoch=1:20
        @show epoch
        @date train(net, xtrn, ytrn; shuffle=true)
        @date ztrn = predict(net,xtrn)
        @date zdev = predict(net,xdev)
        @date ztst = predict(net,xtst)
        @show (seed, epoch, accuracy(ytrn,ztrn), accuracy(ydev,zdev), accuracy(ytst,ztst))
    end
end


text:
module Distances

include("DistanceTypes.jl")

export
    # Function
    GetDistance,

    # Enum
    DistanceType,
    Euclidean,
    Cityblock,
    TotalVariation,
    Chebyshev,
    Jaccard,
    BrayCurtis,
    SpanNormDist
end


text:
##############################################################################
##
## compare
## compare always return a value between 0 and 1.
##
##############################################################################
"""
    compare(s1::AbstractString, s2::AbstractString, dist::PreMetric)

compare returns a similarity score between the strings `s1` and `s2` based on the distance `dist`
"""
function compare(s1::AbstractString, s2::AbstractString,  dist::Hamming; min_score = 0.0)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len2 == 0 && return 1.0
    1.0 - evaluate(dist, s1, s2) / len2
end

function compare(s1::AbstractString, s2::AbstractString, dist::Union{Jaro, RatcliffObershelp}; min_score = 0.0)
    1.0 - evaluate(dist, s1, s2)
end

function compare(s1::AbstractString, s2::AbstractString, dist::AbstractQGramDistance; min_score = 0.0)
    # When string length < q for qgram distance, returns s1 == s2
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len1 <= dist.q - 1 && return convert(Float64, s1 == s2)
    if typeof(dist) <: QGram
        1 - evaluate(dist, s1, s2) / (len1 + len2 - 2 * dist.q + 2)
    else
        1 - evaluate(dist, s1, s2)
    end
end

function compare(s1::AbstractString, s2::AbstractString,  dist::Union{Levenshtein, DamerauLevenshtein}; min_score = 0.0)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len2 == 0 && return 1.0
    if min_score == 0.0
        return 1.0 - evaluate(dist, s1, s2) / len2
    else
        d = evaluate(dist, s1, s2; max_dist = ceil(Int, len2 * (1 - min_score)))
        out = 1.0 - d / len2
        out < min_score && return 0.0
        return out
    end
end

@deprecate compare(dist::PreMetric, s1::AbstractString, s2::AbstractString) compare(s1, s2, dist)

##############################################################################
##
## Winkler
##
##############################################################################
"""
   Winkler(dist::Premetric, p::Real = 0.1, boosting_threshold::Real = 0.7, l::Integer = 4)

Winkler is a `PreMetric` modifier that boosts the similarity score between two strings by a scale `p` when the strings share a common prefix with lenth lower than `l` (the boost is only applied the similarity score above `boosting_threshold`)
"""
struct Winkler{T1 <: PreMetric, T2 <: Real, T3 <: Real, T4 <: Integer} <: PreMetric
    dist::T1
    p::T2          # scaling factor. Default to 0.1
    boosting_threshold::T3      # boost threshold. Default to 0.7
    l::Integer                  # length of common prefix. Default to 4
    function Winkler(dist::T1, p::T2,  boosting_threshold::T3, l::T4) where {T1, T2, T3, T4}
        p * l >= 1 && throw("scaling factor times length of common prefix must be lower than one")
        new{T1, T2, T3, T4}(dist, p, boosting_threshold, l)
    end
end
Winkler(x) = Winkler(x, 0.1, 0.7, 4)

# hard to use min_score because of whether there is boost or not in the end
function compare(s1::AbstractString, s2::AbstractString, dist::Winkler)
    l = remove_prefix(s1, s2, dist.l)[1]
    # cannot do min_score because of boosting threshold
    score = compare(s1, s2, dist.dist)
    if score >= dist.boosting_threshold
        score += l * dist.p * (1 - score)
    end
    return score
end

JaroWinkler() = Winkler(Jaro(), 0.1, 0.7)

##############################################################################
##
## Partial
## http://chairnerd.seatgeek.com/fuzzywuzzy-fuzzy-string-matching-in-python/
##
##############################################################################
"""
   Partial(dist::Premetric)

Partial is a `PreMetric` modifier that returns the maximal similarity score between the shorter string and substrings of the longer string
"""
struct Partial{T <: PreMetric} <: PreMetric
    dist::T
end

function compare(s1::AbstractString, s2::AbstractString, dist::Partial; min_score = 0.0)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len1 == len2 && return compare(s1, s2, dist.dist; min_score = min_score)
    len1 == 0 && return 1.0
    out = 0.0
    for x in qgram(s2, len1)
        curr = compare(s1, x, dist.dist; min_score = min_score)
        out = max(out, curr)
        min_score = max(out, min_score)
    end
    return out
end

function compare(s1::AbstractString, s2::AbstractString, dist::Partial{RatcliffObershelp};
    min_score = 0.0)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    len1 == len2 && return compare(s1, s2, dist.dist)
    out = 0.0
    for r in matching_blocks(s1, s2)
        # Make sure the substring of s2 has length len1
        s2_start = r[2] - r[1] + 1
        s2_end = s2_start + len1 - 1
        if s2_start <= 0
            s2_end += 1 - s2_start
            s2_start += 1 - s2_start
        elseif s2_end > len2
            s2_start += len2 - s2_end
            s2_end += len2 - s2_end
        end
        i2_start =  nextind(s2, 0, s2_start)
        i2_end = nextind(s2, 0, s2_end)
        curr = compare(s1, SubString(s2, i2_start, i2_end), RatcliffObershelp())
        out = max(out, curr)
    end
    return out
end

##############################################################################
##
## TokenSort
## http://chairnerd.seatgeek.com/fuzzywuzzy-fuzzy-string-matching-in-python/
##
##############################################################################
"""
   TokenSort(dist::Premetric)

TokenSort is a `PreMetric` modifier that adjusts for differences in word orders by reording words alphabetically.
"""
struct TokenSort{T <: PreMetric} <: PreMetric
    dist::T
end

function compare(s1::AbstractString, s2::AbstractString, dist::TokenSort; min_score = 0.0)
    s1 = join(sort!(split(s1)), " ")
    s2 = join(sort!(split(s2)), " ")
    compare(s1, s2, dist.dist; min_score = min_score)
end

##############################################################################
##
## TokenSet
## http://chairnerd.seatgeek.com/fuzzywuzzy-fuzzy-string-matching-in-python/
##
##############################################################################
"""
   TokenSet(dist::Premetric)

TokenSort is a `PreMetric` modifier that adjusts for differences in word orders and word numbers by comparing the intersection of two strings with each string.
"""
struct TokenSet{T <: PreMetric} <: PreMetric
    dist::T
end

function compare(s1::AbstractString, s2::AbstractString, dist::TokenSet; min_score = 0.0)
    v1 = SortedSet(split(s1))
    v2 = SortedSet(split(s2))
    v0 = intersect(v1, v2)
    s0 = join(v0, " ")
    s1 = join(v1, " ")
    s2 = join(v2, " ")
    isempty(s0) && return compare(s1, s2, dist.dist; min_score = min_score)
    dist0 = compare(s0, s1, dist.dist; min_score = min_score)
    min_score = max(min_score, dist0)
    dist1 = compare(s0, s2, dist.dist; min_score = min_score)
    min_score = max(min_score, dist1)
    dist2 = compare(s0, s2, dist.dist; min_score = min_score)
    max(dist0, dist1, dist2)
end

##############################################################################
##
## TokenMax
##
##############################################################################
"""
   TokenMax(dist::Premetric)

TokenSort is a `PreMetric` modifier that combines similarlity scores using the base distance, its Partial, TokenSort and TokenSet modifiers, with penalty terms depending on string lengths.
"""
struct TokenMax{T <: PreMetric} <: PreMetric
    dist::T
end

function compare(s1::AbstractString, s2::AbstractString, dist::TokenMax; min_score = 0.0)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    dist0 = compare(s1, s2, dist.dist; min_score = min_score)
    min_score = max(min_score, dist0)
    unbase_scale = 0.95
    # if one string is much shorter than the other, use partial
    if length(s2) >= 1.5 * length(s1)
        partial_scale = length(s2) > (8 * length(s1)) ? 0.6 : 0.9
        dist1 = partial_scale * compare(s1, s2, Partial(dist.dist);
                                        min_score = min_score / partial_scale)
        min_score = max(min_score, dist1)
        dist2 = unbase_scale * partial_scale *
                compare(s1, s2, TokenSort(Partial(dist.dist));
                            min_score = min_score / (unbase_scale * partial_scale))
        min_score = max(min_score, dist2)
        dist3 = unbase_scale * partial_scale *
                compare(s1, s2, TokenSet(Partial(dist.dist));
                            min_score = min_score / (unbase_scale * partial_scale))
        return max(dist0, dist1, dist2, dist3)
    else
        dist1 = unbase_scale *
                compare(s1, s2, TokenSort(dist.dist);
                            min_score = min_score / unbase_scale)
        min_score = max(min_score, dist1)
        dist2 = unbase_scale *
                compare(s1, s2, TokenSet(dist.dist);
                            min_score = min_score / unbase_scale)
        return max(dist0, dist1, dist2)
    end
end

##############################################################################
##
## Missing Values
##
##############################################################################

function compare(s1::AbstractString, ::Missing, dist::PreMetric; min_score = nothing)
    missing
end
function compare(::Missing, s2::AbstractString, dist::PreMetric; min_score = nothing)
    missing
end
function compare(::Missing, ::Missing, dist::PreMetric; min_score = nothing)
    missing
end


text:
function make_operation(@nospecialize(op), args)
    if op === (*)
        args = filter(!_isone, args)
        if isempty(args)
            return 1
        end
    elseif op === (+)
        args = filter(!_iszero, args)
        if isempty(args)
            return 0
        end
    end
    return op(args...)
end

function detime_dvs(op)
    if !istree(op)
        op
    elseif operation(op) isa Sym
        Sym{Real}(nameof(operation(op)))
    else
        similarterm(op, operation(op),detime_dvs.(arguments(op)))
    end
end

function retime_dvs(op::Sym,dvs,iv)
    Sym{FnType{Tuple{symtype(iv)}, Real}}(nameof(op))(iv)
end

function retime_dvs(op, dvs, iv)
    istree(op) ?
        similarterm(op, operation(op), retime_dvs.(arguments(op),(dvs,),(iv,))) :
        op
end

modified_states!(mstates, e::Equation, statelist=nothing) = get_variables!(mstates, e.lhs, statelist)

macro showarr(x)
    n = string(x)
    quote
        y = $(esc(x))
        println($n, " = ", summary(y))
        Base.print_array(stdout, y)
        println()
        y
    end
end

@deprecate substitute_expr!(expr,s) substitute(expr,s)

function states_to_sym(states::Set)
    function _states_to_sym(O)
        if O isa Equation
            Expr(:(=), _states_to_sym(O.lhs), _states_to_sym(O.rhs))
        elseif istree(O)
            op = operation(O)
            args = arguments(O)
            if op isa Sym
                O in states && return tosymbol(O)
                # dependent variables
                return build_expr(:call, Any[nameof(op); _states_to_sym.(args)])
            else
                canonical, O = canonicalexpr(O)
                return canonical ? O : build_expr(:call, Any[op; _states_to_sym.(args)])
            end
        elseif O isa Num
            return _states_to_sym(value(O))
        else
            return toexpr(O)
        end
    end
end
states_to_sym(states) = states_to_sym(Set(states))

function todict(d)
    eltype(d) <: Pair || throw(ArgumentError("The variable-value mapping must be a Dict."))
    d isa Dict ? d : Dict(d)
end

_merge(d1, d2) = merge(todict(d1), todict(d2))


text:
# Hodgkin-Huxley channel gating functions
αn(V) = V == -55.0 ? 0.1 : (0.01*(V + 55))/(1 - exp(-(V + 55)/10))
βn(V) = 0.125 * exp(-(V + 65)/80)
αm(V) = V == -40.0 ? 1.0 : (0.1*(V + 40))/(1 - exp(-(V + 40)/10))
βm(V) = 4*exp(-(V + 65)/18)
αh(V) = 0.07*exp(-(V+65)/20)
βh(V) = 1/(1 + exp(-(V + 35)/10))

# Steady-state channel gating
m∞(V) = αm(V)/(αm(V) + βm(V))
h∞(V) = αh(V)/(αh(V) + βh(V))
n∞(V) = αn(V)/(αn(V) + βn(V))

function stepm!(m, V, dt)
    @. m = (m + αm(V)*dt)/(1+ (αm(V) + βm(V))*dt)
end

function steph!(h, V, dt)
    @. h = (h + αh(V)*dt)/(1+ (αh(V) + βh(V))*dt)
end

function stepn!(n, V, dt)
    @. n = (n + αn(V)*dt)/(1+ (αn(V) + βn(V))*dt)
end

function update_gates!(V, m, h, n, dt)
    stepm!(m, V, dt)
    steph!(h, V, dt)
    stepn!(n, V, dt)
end

function step!(L, D, U, RHS, a, dt, dx, rl, gl, gna, gk, El, Ena, Ek, m, h, n, N, Iapp)
    # First solve Cable equation
    @. D = 1 + (a*dt)/(rl*cm*dx^2) + (dt/cm)*(gna*m^3*h + gk*n^4 + gl)

    # in first compartment we inject a 1 nA current
    RHS[1] = RHS[1] + (dt/cm)*(gna*m[1]^3*h[1]*Ena + gk*n[1]^4*Ek + gl*El) + (Iapp*dt)/(pi*a*cm*dx)
    RHS[2:end] .= @. RHS[2:end] + (dt/cm)*(gna*m[2:end]^3*h[2:end]*Ena + gk*n[2:end]^4*Ek + gl*El)

    #Backward pass
    for i in N:-1:2
        f = U[i] / D[i]
        D[i-1] -= f * L[i]
        RHS[i-1] -= f * RHS[i]
    end

    # Root solve
    RHS[1] /= D[1]

    # Forward pass
    for j in 2:N
        RHS[j] -= L[j]*RHS[j-1]
        RHS[j] /= D[j]
    end

    # Update dimensionless channel gating vars
    update_gates!(RHS, m, h, n, dt)
end

dt = 0.025 # ms
Iapp = 1.0e-6 # mA (1 nA)
cable_length = 2000.0e-3 # 2000 µm length
N = 2000
dx = cable_length/N
a   = 2.0e-3 # µm in mm
area = 2*pi*a*dx

rl = 1.0e6 # 1 kΩ in milliohms*mm
gl  = 0.003 #mS/mm^2
gna = 1.2 # mS/mm^2 (modified from 1.2 mS)
gk  = .36 # mS/mm^2
El  = -54.4 # mV
Ena = 50.0 # mV
Ek  = -77.0 # mV
cm  = 1.0e-5 # mF/mm^2
Vinit = -55.0 # mV

Li = -(a*dt)/(2*rl*cm*dx^2)
LN = -(a*dt)/(rl*cm*dx^2)
Ui = Li
U1 = LN

L = Vector{Float64}(undef, N)
U = Vector{Float64}(undef, N)
D = Vector{Float64}(undef, N)
RHS = fill(Vinit,(N))
m = fill(m∞(Vinit), (N))
h = fill(h∞(Vinit), (N))
n = fill(n∞(Vinit), (N))
L[1] = 0; L[N] = LN; L[2:end-1] .= Li
U[1] = U1; U[N] = 0; U[2:end-1] .= Ui
history = Vector{Vector{Float64}}()

@time begin
    for i in 1:8000
        step!(L, D, U, RHS, a, dt, dx, rl, gl, gna, gk, El, Ena, Ek, m, h, n, N, Iapp)
    end
end

using Plots

plot(RHS)

max_stars_repo_name: wsphillips/FutureNeuralSimulation

max_issues_repo_name: wsphillips/FutureNeuralSimulation

max_forks_repo_name: wsphillips/FutureNeuralSimulation


text:
using Genie

import AdminController

route("/admin", AdminController.index; method = GET, named = :admin)
route("/admin/:model::String", AdminController.create; method = POST, named = :admin_model_create)
route("/admin/:model::String/:id::Int", AdminController.read; method = GET, named = :admin_model_read)
route("/admin/:model::String/:id::Int", AdminController.update; method = POST, named = :admin_model_update)
route("/admin/:model::String/:id::Int/delete", AdminController.delete; method = GET, named = :admin_model_delete)
route("/admin/:model::String", AdminController.list; method = GET, named = :admin_model_list)


text:
# calculate required values from the solution (for the animation)

derive!(t::TrebuchetState, sol::AbstractODESolution) = derive(t, sol, t.sol)
function derive(t::TrebuchetState, sol::AbstractODESolution, s = Solution())
    stage = t.stage
    p = tuples(sol)
    for (x, y) in p
        derive(t, x, y, stage, s)
    end
    return s
end

function derive(t::TrebuchetState, a::Array, time, ::Union{Val{:Ground},Val{:Hang}}, s=Solution())
    (Aq, Wq, Sq,) = a
    l = t.l
    LAl, LAs, LW, LS, h = l.b, l.c, l.d, l.e, l.a
    LAcg = (LAl - LAs)/2
    SIN = sin
    COS = cos
    push!(s.WeightCG, [LAs*SIN(Aq) + LW*SIN(Aq+Wq), -LAs*COS(Aq) - LW*COS(Aq+Wq)])
    push!(s.WeightArm, [LAs*SIN(Aq), -LAs*COS(Aq)])
    push!(s.ArmSling, [-LAl*SIN(Aq), LAl*COS(Aq)])
    push!(s.Projectile, [-LAl*SIN(Aq) - LS*SIN(Aq+Sq), LAl*COS(Aq) + LS*COS(Aq+Sq)])
    push!(s.SlingEnd,[-LAl*SIN(Aq) - LS*SIN(Aq+Sq), LAl*COS(Aq) + LS*COS(Aq+Sq)])
    push!(s.ArmCG, [-LAcg*SIN(Aq), LAcg*COS(Aq)])
    push!(s.Time, time)
    return s
end

function derive(t::TrebuchetState, a::Array, time, ::Val{:Released}, s=Solution())
    (Px, Py,) = a
    l, an = t.l, t.a
    LAl, LAs, LW, LS, h = l.b, l.c, l.d, l.e, l.a
    Aq, Wq, Sq = an.aq, an.wq, an.sq
    LAcg = (LAl - LAs)/2
    SIN = sin
    COS = cos
    push!(s.WeightCG, [LAs*SIN(Aq) + LW*SIN(Aq+Wq), -LAs*COS(Aq) - LW*COS(Aq+Wq)])
    push!(s.WeightArm, [LAs*SIN(Aq), -LAs*COS(Aq)])
    push!(s.ArmSling, [-LAl*SIN(Aq), LAl*COS(Aq)])
    push!(s.Projectile, [Px, Py])
    push!(s.SlingEnd,[-LAl*SIN(Aq) - LS*SIN(Aq+Sq), LAl*COS(Aq) + LS*COS(Aq+Sq)])
    push!(s.ArmCG, [-LAcg*SIN(Aq), LAcg*COS(Aq)])
    push!(s.Time, time)
    return s
end


text:
import LowLevelJPEG: @code_str

@testset "Constructors" begin
    @test CodeBits() == CodeBits(0, 0x1)
    @test CodeBits(0x0006000D).data === 0x6000D
    @test CodeBits(6, 0b001101).data === 0x6000D
    @test CodeBits(0, 0b001101).data === 0x00000
end

@testset "Traits" begin
    @test LowLevelJPEG.word(code"001101") === UInt16(0b001101)
    @test LowLevelJPEG.isvalid(code"1000100010001000")
    @test !LowLevelJPEG.isvalid(CodeBits(0x10000000))
    @test length(code"") == 0
    @test length(code"001101") == 6
    @test size(code"001101") == (6,)

    @test code"01"[1] == false
end

@testset "Display" begin
    @test shown_string(code"") == "code\"\""
    @test shown_string(code"001101") == "code\"001101\""
    @test shown_string(CodeBits(0x11FFFF)) == "code\"?\""
end

@testset "Functions" begin
    @test LowLevelJPEG.inc(code"001101") === code"001110"
    @test LowLevelJPEG.shl(code"001101") === code"0011010"

    @test BitArray(code"001101") == BitArray((0, 0, 1, 1, 0, 1))
end


text:
field(ca::CA, arg) where {CA<:Union{Caggregate, Carray}} = field(typeof(ca), arg)

field(::Type{CA}, i::Integer) where {T, N, CA<:Carrays{T, N}} = Tuple{T, bitstype(T), integraltype(T), sizeof(Carray{T, i-1}), 0, -1}

field(::Type{CA}, sym::Symbol) where {CA<:Caggregates} = field(fields(bitstype(CA)), Val(sym))
@generated function field(::Type{fields}, ::Val{sym}) where {fields<:Tuple, sym}
	for T in fields
		(n, t) = T
		n === sym && return t
	end
	return nothing
end

function bitmask(::Type{uint}, bits::Integer) where {uint}
	mask = zero(uint)
	for i in 1:bits
		mask = (mask << one(uint)) | one(uint)
	end
	return uint(mask)
end

function loadmasked(::Type{int}, isptr::Bool, bytes::Integer, bits::Integer, size::Integer) where {int<:Integer}
	expr = []

	uint = unsigned(int)
	m = bitmask(uint, size)
	push!(expr, :(load = (load >> $(uint)($(bits))) & $(uint)($(m))))
	int <: Signed && push!(expr, :(load = Bool((load >> $(size-1)) & 1) ? load | (~zero($(uint)) & ~$(uint)($(m))) : load))

	return expr
end

function storemasked(::Type{int}, isptr::Bool, bytes::Integer, bits::Integer, size::Integer) where {int<:Integer}
	expr = []

	uint = unsigned(int)
	mask = bitmask(uint, size) << bits
	append!(expr, loadbytes(int, isptr, bytes, bits, size))
	push!(expr, :(store = ((store << $(bits)) & $(uint)($(mask))) | (load & ~$(uint)($(mask)))))

	return expr
end

function loadbytes(::Type{int}, isptr::Bool, bytes::Integer, bits::Integer, size::Integer) where {int<:Integer}
	expr = []

	uint = unsigned(int)
	size = size < 0 ? sizeof(uint)*8 : size

	push!(expr, :(load = zero($(uint))))
	for i in (ENDIAN_BOM == 0x01020304 ? (1:sizeof(uint)) : reverse(1:sizeof(uint)))
		push!(expr, :(load = (load << 8)))
		if bits÷8 <= i-1 <= (bits+size-1)÷8
			if isptr
				push!(expr, :(load |= $(uint)(unsafe_load(Core.Intrinsics.bitcast(Cptr{UInt8}, ca), $(bytes+i)))))
			else
				push!(expr, :(load |= $(uint)(getfield(getfield(ca, :mem), $(bytes+i)))))
			end
		end
	end

	return expr
end

function storebytes(::Type{int}, casize::Integer, isptr::Bool, bytes::Integer, bits::Integer, size::Integer) where {int<:Integer}
	expr = []

	uint = unsigned(int)
	size = size < 0 ? sizeof(uint)*8 : size

	isptr || push!(expr, :(mem = getfield(ca, :mem)))
	for i in (ENDIAN_BOM == 0x01020304 ? reverse(1:sizeof(uint)) : (1:sizeof(uint)))
		if bits÷8 <= i-1 <= (bits+size-1)÷8
			if isptr
				push!(expr, :(unsafe_store!(Core.Intrinsics.bitcast(Cptr{UInt8}, ca), store & 0xff, $(bytes+i))))
			else
				mem = Expr(:tuple,
					(x == bytes+i ? :(UInt8(store & 0xff)::UInt8) : :(getfield(mem, $(x))) for x in 1:casize)...,
				)
				push!(expr, :(mem = $(mem)))
			end
		end
		push!(expr, :(store = store >> 8))
	end

	return expr
end

function load(ca::CA, ::Type{Tuple{t, u, i, o, b, s}}) where {CA<:Union{Caggregate, Carray, Cptr{<:Caggregates}}, t, u, i, o, b, s}
	return load(ca, Tuple{t, u, i, o, b, s}, Val(!(CA <: Cptr) || s >= 0))
end

@generated function load(ca::CA, ::Type{Tuple{t, u, i, o, b, s}}, ::Val) where {CA<:Union{Caggregate, Carray}, t, u, i<:Nothing, o, b, s}
	u <: Union{Caggregate, Carray} || return error("Unexpected type $(t) in load function")

	expr = quote end

	mem = Expr(:tuple,
		(:(getfield(getfield(ca, :mem), $(x))) for x in o+1:o+sizeof(u))...,
	)
	push!(expr.args, :(return u(undef, $(mem))))

	return expr
end

function load(ca::CA, ::Type{Tuple{t, u, i, o, b, s}}, ::Val{false}) where {CA<:Cptr{<:Caggregates}, t, u, i, o, b, s}
	i <: Integer || u <: Union{Caggregate, Carray} || return error("Unexpected type $(t) in load function")
	x = isqualifiedwith(eltype(CA), Cconst) && !isqualifiedwith(t, Cconst) ? Cconst{t} : t
	return Core.Intrinsics.bitcast(Cptr{x}, Core.Intrinsics.bitcast(Cptr{UInt8}, ca) + o)
end

@generated function load(ca::CA, ::Type{Tuple{t, u, i, o, b, s}}, ::Val{true}) where {CA<:Union{Caggregate, Carray, Cptr{<:Caggregates}}, t, u, i<:Integer, o, b, s}
	expr = quote end

	append!(expr.args, loadbytes(i, CA <: Cptr, o, b, s))
	s >= 0 && append!(expr.args, loadmasked(i, CA <: Cptr, o, b, s))
	push!(expr.args, :(return Core.Intrinsics.bitcast($(u), load)))

	return expr
end

store!(ca::CA, ::Type{Tuple{t, u, i, o, b, s}}, val) where {CA<:Union{Caggregate, Carray, Cptr{<:Caggregates}}, t, u, i, o, b, s} = store!(ca, Tuple{t, u, i, o, b, s}, Val(!(CA <: Cptr) || s >= 0), val)

@generated function store!(ca::CA, ::Type{Tuple{t, u, i, o, b, s}}, ::Val, val) where {CA<:Union{Caggregate, Carray}, t, u, i<:Nothing, o, b, s}
	u <: Union{Caggregate, Carray} || return error("Unexpected type $(t) in store! function")

	expr = quote end

	push!(expr.args, :(v = convert(u, val)))
	mem = Expr(:tuple,
		(x-o in 1:sizeof(u) ? :(getfield(getfield(v, :mem), $(x-o))) : :(getfield(getfield(ca, :mem), $(x))) for x in 1:sizeof(CA))...,
	)
	push!(expr.args, :(return CA(undef, $(mem))))

	return expr
end

function store!(ca::CA, ::Type{Tuple{t, u, i, o, b, s}}, ::Val{false}, val) where {CA<:Cptr{<:Caggregates}, t, u, i, o, b, s}
	u <: Union{Caggregate, Carray} || return error("Unexpected type $(t) in store! function")
	unsafe_store!(Core.Intrinsics.bitcast(Cptr{u}, Core.Intrinsics.bitcast(Cptr{UInt8}, ca) + o), convert(u, val))
	return val
end

@generated function store!(ca::CA, ::Type{Tuple{t, u, i, o, b, s}}, ::Val{true}, val) where {CA<:Union{Caggregate, Carray, Cptr{<:Caggregates}}, t, u, i<:Integer, o, b, s}
	expr = quote end

	uint = unsigned(i)
	push!(expr.args, :(store = Core.Intrinsics.bitcast($(uint), convert($(u), val))))
	s >= 0 && append!(expr.args, storemasked(i, CA <: Cptr, o, b, s))
	append!(expr.args, storebytes(i, sizeof(CA), CA <: Cptr, o, b, s))
	push!(expr.args, CA <: Cptr ? :(return val) : :(return CA(undef, mem)))

	return expr
end


text:
module miscellaneous_examples
using FinEtools

function circular_hole()

    println("""    Meshing  for an elliptical hole in a finite strip.        """)

    t0 = time()

    H = 100. # strip width
    ell = 50. # distance between holes
    a = 10. # radius of the hole
    L = 200. # length of the strip
    nL = 15
    nH = 10
    nR = 50
    fens,fes = Q4elliphole(a, a, L/2, H/2, nL, nH, nR)
    # Note that we have to number the quadrilaterals  in the mirrored mesh
    # differently in order to generate positive areas.
    fens2, fes2 = mirrormesh(fens, fes, [0.0; -1.0],
    [0.0; 0.0]; renumb = c -> c[end:-1:1])
    fens, fes1, fes2 = mergemeshes(fens, fes,  fens2, fes2, a/1000.)
    fes = cat(fes1, fes2)

    geom = NodalField(fens.xyz)

    File =  "elliptical.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.Q4)
    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # circular_hole

function L_shaped_mesh()
    println("""
    Meshing an L-shaped membrane: merging multiple meshes.
    """)

    t0 = time()

    W = 100. # width of the leg
    L = 200. # length of the leg
    nL = 15 # number of elements along the length of the leg
    nW = 10 # number of elements along the width
    tolerance = W / nW / 1.0e5 # tolerance for merging nodes
    Meshes = Array{Tuple{FENodeSet,FESet},1}()
    push!(Meshes, Q4quadrilateral([0.0 0.0; W W], nW, nW))
    push!(Meshes, Q4quadrilateral([-L 0.0; 0.0 W], nL, nW))
    push!(Meshes, Q4quadrilateral([0.0 -L; W 0.0], nW, nL))
    fens, outputfes = mergenmeshes(Meshes, tolerance);
    fes = cat(outputfes[1], cat(outputfes[2], outputfes[3]))

    geom = NodalField(fens.xyz)

    File =  "L_shape.vtk"
    vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.Q4);
    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # L_shaped_mesh

function composite_plate_mesh()
    println(""" Meshing  for a layered (composite) plate. """)

    t0 = time()

    H = 100. # strip width
    L = 200. # length of the strip
    nL = 15
    nH = 10
    xs = collect(linspace(0.0, L, nL+1))
    ys = collect(linspace(0.0, H, nH+1))
    ts = [0.5; 2.0; 0.5; 2.0]
    nts = [1; 2; 2; 1; ]
    fens,fes = H8layeredplatex(xs, ys, ts, nts)

    geom = NodalField(fens.xyz)

    File =  "composite.vtk"
    vtkexportmesh(File, fes.conn, geom.values, FinEtools.MeshExportModule.H8;
    scalars = [( "Layer", fes.label)])
    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # composite_plate_mesh

function boundary_Q4_example()
    t0 = time()

    rho=1.21*1e-9;# mass density
    c =345.0*1000;# millimeters per second
    bulk= c^2*rho;
    Lx=1900.0;# length of the box, millimeters
    Ly=800.0; # length of the box, millimeters

    fens,fes = Q4block(Lx,Ly,3,2); # Mesh
    show(fes.conn)

    bfes = meshboundary(fes)
    File =  "B.vtk"
    vtkexportmesh(File, bfes.conn, fens.xyz, FinEtools.MeshExportModule.L2)
    @async run(`"paraview.exe" $File`)

end # boundary_Q4_example

function allrun()
    println("#####################################################")
    println("# circular_hole ")
    circular_hole()
    println("#####################################################")
    println("# L_shaped_mesh ")
    L_shaped_mesh()
    println("#####################################################")
    println("# composite_plate_mesh ")
    composite_plate_mesh()
    println("#####################################################")
    println("# boundary_Q4_example ")
    boundary_Q4_example()
    return true
end # function allrun

end # module miscellaneous_examples


text:
# License for this file: MIT (expat)
# Copyright 2020-2021, DLR Institute of System Dynamics and Control
#
# This file is part of module ModiaResult (is included by the ModiaPlot_XXMakie packages).

export plot, showFigure, saveFigure, closeFigure, closeAllFigures

#--------------------------- Utility definitions

# color definitions
const colors = ["blue1",
                "green4",
                "red1",
                "cyan1",
                "purple1",
                "orange1",
                "black",
                "blue2",
                "green3",
                "red2",
                "cyan2",
                "purple2",
                "orange2",
                "grey20",
                "blue3",
                "green2",
                "red3",
                "cyan3",
                "purple3",
                "orange3",
                "grey30",
                "blue4",
                "green1",
                "red4",
                "cyan4",
                "purple4",
                "orange4",
                "grey40"]

mutable struct Diagram
    fig
    axis
    row::Int
    col::Int
    curves::Union{AbstractVector, Nothing}
    yLegend::Union{AbstractVector, Nothing}

    function Diagram(fig,row,col,title)
        axis = fig[row,col] = Axis(fig, title=title)
        new(fig, axis, row, col, nothing, nothing)
    end
end

# figures[i] is Figure instance of figureNumber = i
const figures = Dict{Int,Any}()   # Dictionary of MatrixFigures

mutable struct MatrixFigure
    fig
    resolution
    diagrams::Union{Matrix{Diagram},Nothing}

    function MatrixFigure(figureNumber::Int, nrow::Int, ncol::Int, reuse::Bool, resolution)
        if haskey(figures, figureNumber) && reuse
            # Reuse existing figure
            matrixFigure = figures[figureNumber]
            if isnothing(matrixFigure.diagrams) || size(matrixFigure.diagrams) != (nrow,ncol)
                error("From ModiaPlot.plot: reuse=true, but figure=$figureNumber has not $nrow rows and $ncol columns.")
            end
        else
            # Create a new figure
            fig = Figure(resolution=resolution)
            diagrams = Matrix{Diagram}(undef,nrow,ncol)
            matrixFigure = new(fig, resolution, diagrams)
            figures[figureNumber] = matrixFigure
        end
        return matrixFigure
    end
end

getMatrixFigure(figureNumber::Int) = figures[figureNumber]

setTheme() = set_theme!(fontsize = 12)

"""
    getColor(i::Int)

Return color, given an integer (from colors)
"""
function getColor(i::Int)
    j = mod(i, length(colors))
    if j == 0
        j = length(colors)
    end
    return colors[j]
end

"""
    createLegend(diagram::Diagram, maxLegend)

Create a new legend in diagram
"""
function createLegend(diagram::Diagram, maxLegend::Int)::Nothing
    curves  = diagram.curves
    yLegend = diagram.yLegend
    if length(curves) > 1
        nc = min(length(curves), maxLegend)
        curves  = curves[1:nc]
        yLegend = yLegend[1:nc]
    end

    fig  = diagram.fig
    i    = diagram.row
    j    = diagram.col
    fig[i,j] = Legend(fig, curves, yLegend,
                              margin=[0,0,0,0], padding=[5,5,0,0], rowgap=0, halign=:right, valign=:top,
                              linewidth=1, labelsize=10, tellheight=false, tellwidth=false)
    return nothing
end

const reltol = 0.001

function sig_max(xsig, y_min, y_max)
    dy = y_max - y_min
    abstol = typemax(eltype(y_min))
    for (i,v) in enumerate(dy)
        if v > 0.0 && v < abstol
            abstol = v
        end
    end
    @assert(abstol > 0.0)
    y_max2 = similar(y_max)
    for (i,value) = enumerate(y_max)
        y_max2[i] = max(y_max[i], y_min[i] + max(abstol, abs(y_min[i])*reltol))
    end
    return y_max2
end

# Something of the following is required in the file where makie.jl is included:
#   const Makie_Point2f = isdefined(GLMakie, :Point2f) ? Point2f : Point2f0
function fill_between(axis, xsig, ysig_min, ysig_max, color)
    ysig_max2 = sig_max(xsig, ysig_min, ysig_max)
    sig = Makie_Point2f.(xsig,ysig_min)
    append!(sig, reverse(Makie_Point2f.(xsig,ysig_max2)))
    push!(sig, sig[1])
    return poly!(axis, sig, color = color)
end

function plotOneSignal(axis, xsig, ysig, color, ysigType, MonteCarloAsArea)
	if typeof(ysig[1]) <: Measurements.Measurement
		# Plot mean value signal
		xsig_mean = Measurements.value.(xsig)
		ysig_mean = Measurements.value.(ysig)
        curve = lines!(axis, xsig_mean, ysig_mean, color=color)

		# Plot area of uncertainty around mean value signal (use the same color, but transparent)
		ysig_u   = Measurements.uncertainty.(ysig)
		ysig_max = ysig_mean + ysig_u
		ysig_min = ysig_mean - ysig_u
		fill_between(axis, xsig_mean, ysig_min, ysig_max, (color,0.2))

	elseif typeof(ysig[1]) <: MonteCarloMeasurements.StaticParticles ||
           typeof(ysig[1]) <: MonteCarloMeasurements.Particles
        # Plot mean value signal
        pfunctionsDefined = isdefined(MonteCarloMeasurements, :pmean)
        if pfunctionsDefined
            # MonteCarlMeasurements, version >= 1.0
            xsig_mean = MonteCarloMeasurements.pmean.(xsig)
            ysig_mean = MonteCarloMeasurements.pmean.(ysig)
        else
            # MonteCarloMeasurements, version < 1.0
            xsig_mean = MonteCarloMeasurements.mean.(xsig)
            ysig_mean = MonteCarloMeasurements.mean.(ysig)
        end
        xsig_mean = ustrip.(xsig_mean)
        ysig_mean = ustrip.(ysig_mean)
        curve = lines!(axis, xsig_mean, ysig_mean, color=color)

        if MonteCarloAsArea
            # Plot area of uncertainty around mean value signal (use the same color, but transparent)
            if pfunctionsDefined
                # MonteCarlMeasurements, version >= 1.0
                ysig_max = MonteCarloMeasurements.pmaximum.(ysig)
                ysig_min = MonteCarloMeasurements.pminimum.(ysig)
            else
                # MonteCarloMeasurements, version < 1.0
                ysig_max = MonteCarloMeasurements.maximum.(ysig)
                ysig_min = MonteCarloMeasurements.minimum.(ysig)
            end
            ysig_max = ustrip.(ysig_max)
            ysig_min = ustrip.(ysig_min)
            fill_between(axis, xsig_mean, ysig_min, ysig_max, (color,0.2))

        else
            # Plot all particle signals
            value  = ysig[1].particles
            ysig3  = zeros(eltype(value), length(xsig))
            for j in 1:length(value)
                for i in eachindex(ysig)
                    ysig3[i] = ysig[i].particles[j]
                end
                ysig3 = ustrip.(ysig3)
                lines!(axis, xsig, ysig3, color=(color,0.1))
            end
        end

	else
        if typeof(xsig[1]) <: Measurements.Measurement
            xsig = Measurements.value.(xsig)
        elseif typeof(xsig[1]) <: MonteCarloMeasurements.StaticParticles ||
               typeof(xsig[1]) <: MonteCarloMeasurements.Particles
            if isdefined(MonteCarloMeasurements, :pmean)
                # MonteCarlMeasurements, version >= 1.0
                xsig = MonteCarloMeasurements.pmean.(xsig)
            else
                # MonteCarlMeasurements, version < 1.0
                xsig = MonteCarloMeasurements.mean.(xsig)
            end
            xsig = ustrip.(xsig)
        end

        if ysigType == ModiaResult.Clocked
            curve = scatter!(axis, xsig, ysig, color=color, markersize = 2.0)
        else
            curve = lines!(axis, xsig, ysig, color=color)
        end
	end
    return curve
end

function addPlot(names::Tuple, diagram::Diagram, result, grid::Bool, xLabel::Bool, xAxis,
                 prefix::AbstractString, reuse::Bool, maxLegend::Integer, MonteCarloAsArea::Bool)
    xsigLegend = ""
    yLegend = String[]
    curves  = Any[]

    i0 = isnothing(diagram.curves) ? 0 : length(diagram.curves)

    for name in names
        name2 = string(name)
        (xsig, xsigLegend, ysig, ysigLegend, ysigType) = ModiaResult.getPlotSignal(result, name2; xsigName = xAxis)
        if !isnothing(xsig)
            if ndims(ysig) == 1
                push!(yLegend, prefix*ysigLegend[1])
                push!(curves, plotOneSignal(diagram.axis, xsig, ysig, getColor(i0+1), ysigType, MonteCarloAsArea))
                i0 = i0 + 1
            else
                for i = 1:size(ysig,2)
                    curve = plotOneSignal(diagram.axis, xsig, ysig[:,i], getColor(i0+i), ysigType, MonteCarloAsArea)
                    push!(yLegend, prefix*ysigLegend[i])
                    push!(curves, curve)
                end
                i0 = i0 + size(ysig,2)
            end
        end
    end

    #PyPlot.grid(grid)

    if reuse
        diagram.curves  = append!(diagram.curves, curves)
        diagram.yLegend = append!(diagram.yLegend, yLegend)
    else
        diagram.curves  = curves
        diagram.yLegend = yLegend
    end

    createLegend(diagram, maxLegend)

    if xLabel && !reuse && xsigLegend !== nothing
        diagram.axis.xlabel = xsigLegend
    end
end

addPlot(name::AbstractString, args...) = addPlot((name,)        , args...)
addPlot(name::Symbol        , args...) = addPlot((string(name),), args...)

#--------------------------- Plot function

function plot(result, names::AbstractMatrix; heading::AbstractString="", grid::Bool=true, xAxis=nothing,
              figure::Int=1, prefix::AbstractString="", reuse::Bool=false, maxLegend::Integer=10,
              minXaxisTickLabels::Bool=false, MonteCarloAsArea=false)
    if isnothing(result)
        @info "The call of ModiaPlot.plot(result, ...) is ignored, since the first argument is nothing."
        return
    end
    if !reusePossible
        reuse=false
    end

    # Plot a vector or matrix of diagrams
    setTheme()
    (nrow, ncol)  = size(names)
    matrixFigure  = MatrixFigure(figure, nrow, ncol, reuse, (min(ncol*600,1500), min(nrow*350, 900)))
    fig           = matrixFigure.fig
    xAxis2        = isnothing(xAxis) ? xAxis : string(xAxis)
    heading2      = ModiaResult.getHeading(result, heading)
    hasTopHeading = !reuse && ncol > 1 && heading2 != ""

    # Add diagrams
    lastRow = true
    xLabel  = true
    for i = 1:nrow
        lastRow = i == nrow
        for j = 1:ncol
            if reuse
                diagram = matrixFigure.diagrams[i,j]
            else
                if ncol == 1 && i == 1 && !hasTopHeading
                    diagram = Diagram(fig, i, j, heading2)
                else
                    diagram = Diagram(fig, i, j, "")
                end
                matrixFigure.diagrams[i,j] = diagram
            end
            xLabel = !( minXaxisTickLabels && !lastRow )
            addPlot(names[i,j], diagram, result, grid, xLabel, xAxis2, prefix, reuse, maxLegend, MonteCarloAsArea)
        end
    end

    # Add overall heading in case of a matrix of diagrams (ncol > 1) and add a figure label on the top level
    if hasTopHeading
        fig[0,:] = Label(fig, heading2, textsize = 14)
    end
    if showFigureStringInDiagram
        figText = fig[1,1,TopLeft()] = Label(fig, "showFigure(" * string(figure) * ")", textsize=9, color=:blue, halign = :left)
        if hasTopHeading
            figText.padding = (0, 0, 5, 0)
        else
            figText.padding = (0, 0, 30, 0)
        end
    end

    # Update and display fig
    trim!(fig.layout)
    update!(fig.scene)
    if callDisplayFunction
        display(fig)
    end
    return matrixFigure
end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
using Streams
using Test

s1 = Stream(1)
s2 = Stream(2)
s3 = TimeStream{Float64}(sin)
s4 = lift(+, s1, s2, s3)
@info "initial value"
@show s4

@info "changing one addend"
s2[] = 100
@show s4


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

module Enums

import Core.Intrinsics.bitcast
export Enum, @enum

function namemap end

"""
    Enum{T<:Integer}

The abstract supertype of all enumerated types defined with [`@enum`](@ref).
"""
abstract type Enum{T<:Integer} end

basetype(::Type{<:Enum{T}}) where {T<:Integer} = T

(::Type{T})(x::Enum{T2}) where {T<:Integer,T2<:Integer} = T(bitcast(T2, x))::T
Base.cconvert(::Type{T}, x::Enum{T2}) where {T<:Integer,T2<:Integer} = T(x)
Base.write(io::IO, x::Enum{T}) where {T<:Integer} = write(io, T(x))
Base.read(io::IO, ::Type{T}) where {T<:Enum} = T(read(io, basetype(T)))

Base.isless(x::T, y::T) where {T<:Enum} = isless(basetype(T)(x), basetype(T)(y))

Base.Symbol(x::Enum) = namemap(typeof(x))[Integer(x)]::Symbol

Base.print(io::IO, x::Enum) = print(io, Symbol(x))

function Base.show(io::IO, x::Enum)
    sym = Symbol(x)
    if !(get(io, :compact, false)::Bool)
        from = get(io, :module, Main)
        def = typeof(x).name.module
        if from === nothing || !Base.isvisible(sym, def, from)
            show(io, def)
            print(io, ".")
        end
    end
    print(io, sym)
end

function Base.show(io::IO, ::MIME"text/plain", x::Enum)
    print(io, x, "::")
    show(IOContext(io, :compact => true), typeof(x))
    print(io, " = ")
    show(io, Integer(x))
end

function Base.show(io::IO, ::MIME"text/plain", t::Type{<:Enum})
    print(io, "Enum ")
    Base.show_datatype(io, t)
    print(io, ":")
    for x in instances(t)
        print(io, "\n", Symbol(x), " = ")
        show(io, Integer(x))
    end
end

# generate code to test whether expr is in the given set of values
function membershiptest(expr, values)
    lo, hi = extrema(values)
    if length(values) == hi - lo + 1
        :($lo <= $expr <= $hi)
    elseif length(values) < 20
        foldl((x1,x2)->:($x1 || ($expr == $x2)), values[2:end]; init=:($expr == $(values[1])))
    else
        :($expr in $(Set(values)))
    end
end

# give Enum types scalar behavior in broadcasting
Base.broadcastable(x::Enum) = Ref(x)

@noinline enum_argument_error(typename, x) = throw(ArgumentError(string("invalid value for Enum $(typename): $x")))

"""
    @enum EnumName[::BaseType] value1[=x] value2[=y]

Create an `Enum{BaseType}` subtype with name `EnumName` and enum member values of
`value1` and `value2` with optional assigned values of `x` and `y`, respectively.
`EnumName` can be used just like other types and enum member values as regular values, such as

# Examples
```jldoctest fruitenum
julia> @enum Fruit apple=1 orange=2 kiwi=3

julia> f(x::Fruit) = "I'm a Fruit with value: \$(Int(x))"
f (generic function with 1 method)

julia> f(apple)
"I'm a Fruit with value: 1"

julia> Fruit(1)
apple::Fruit = 1
```

Values can also be specified inside a `begin` block, e.g.

```julia
@enum EnumName begin
    value1
    value2
end
```

`BaseType`, which defaults to [`Int32`](@ref), must be a primitive subtype of `Integer`.
Member values can be converted between the enum type and `BaseType`. `read` and `write`
perform these conversions automatically.

To list all the instances of an enum use `instances`, e.g.

```jldoctest fruitenum
julia> instances(Fruit)
(apple, orange, kiwi)
```
"""
macro enum(T, syms...)
    if isempty(syms)
        throw(ArgumentError("no arguments given for Enum $T"))
    end
    basetype = Int32
    typename = T
    if isa(T, Expr) && T.head === :(::) && length(T.args) == 2 && isa(T.args[1], Symbol)
        typename = T.args[1]
        basetype = Core.eval(__module__, T.args[2])
        if !isa(basetype, DataType) || !(basetype <: Integer) || !isbitstype(basetype)
            throw(ArgumentError("invalid base type for Enum $typename, $T=::$basetype; base type must be an integer primitive type"))
        end
    elseif !isa(T, Symbol)
        throw(ArgumentError("invalid type expression for enum $T"))
    end
    values = basetype[]
    seen = Set{Symbol}()
    namemap = Dict{basetype,Symbol}()
    lo = hi = 0
    i = zero(basetype)
    hasexpr = false

    if length(syms) == 1 && syms[1] isa Expr && syms[1].head === :block
        syms = syms[1].args
    end
    for s in syms
        s isa LineNumberNode && continue
        if isa(s, Symbol)
            if i == typemin(basetype) && !isempty(values)
                throw(ArgumentError("overflow in value \"$s\" of Enum $typename"))
            end
        elseif isa(s, Expr) &&
               (s.head === :(=) || s.head === :kw) &&
               length(s.args) == 2 && isa(s.args[1], Symbol)
            i = Core.eval(__module__, s.args[2]) # allow exprs, e.g. uint128"1"
            if !isa(i, Integer)
                throw(ArgumentError("invalid value for Enum $typename, $s; values must be integers"))
            end
            i = convert(basetype, i)
            s = s.args[1]
            hasexpr = true
        else
            throw(ArgumentError(string("invalid argument for Enum ", typename, ": ", s)))
        end
        s = s::Symbol
        if !Base.isidentifier(s)
            throw(ArgumentError("invalid name for Enum $typename; \"$s\" is not a valid identifier"))
        end
        if hasexpr && haskey(namemap, i)
            throw(ArgumentError("both $s and $(namemap[i]) have value $i in Enum $typename; values must be unique"))
        end
        namemap[i] = s
        push!(values, i)
        if s in seen
            throw(ArgumentError("name \"$s\" in Enum $typename is not unique"))
        end
        push!(seen, s)
        if length(values) == 1
            lo = hi = i
        else
            lo = min(lo, i)
            hi = max(hi, i)
        end
        i += oneunit(i)
    end
    blk = quote
        # enum definition
        Base.@__doc__(primitive type $(esc(typename)) <: Enum{$(basetype)} $(sizeof(basetype) * 8) end)
        function $(esc(typename))(x::Integer)
            $(membershiptest(:x, values)) || enum_argument_error($(Expr(:quote, typename)), x)
            return bitcast($(esc(typename)), convert($(basetype), x))
        end
        Enums.namemap(::Type{$(esc(typename))}) = $(esc(namemap))
        Base.typemin(x::Type{$(esc(typename))}) = $(esc(typename))($lo)
        Base.typemax(x::Type{$(esc(typename))}) = $(esc(typename))($hi)
        let insts = (Any[ $(esc(typename))(v) for v in $values ]...,)
            Base.instances(::Type{$(esc(typename))}) = insts
        end
    end
    if isa(typename, Symbol)
        for (i, sym) in namemap
            push!(blk.args, :(const $(esc(sym)) = $(esc(typename))($i)))
        end
    end
    push!(blk.args, :nothing)
    blk.head = :toplevel
    return blk
end

end # module


text:
using Plugins, BenchmarkTools

mutable struct CounterPlugin <: Plugin
    count::UInt
    CounterPlugin() = new(0)
end
Plugins.symbol(::CounterPlugin) = :counter

function tick(me::CounterPlugin, app)
    me.count += 1
    return false
end

mutable struct App
    plugins::PluginStack
    App(plugins, hookfns) = new(PluginStack(plugins, hookfns))
end

function tickerop(app::App)
    tickhook = hooks(app).tick
    tickerop_kern(app, tickhook) # Using a function barrier to get ~7ns per hook activation
end

function tickerop_kern(app::App, tickhook)
    for i = 1:1e6
        tickhook()
    end
end

const app = app = App([CounterPlugin()], [tick])
@btime tickerop(app)
push!(app.plugins, CounterPlugin())
println("Total Tick Count: $(app.plugins[:counter].count)")


text:
if length(ARGS) != 2
    error("Usage: julia --project=. thisfile.jl 127.0.0.1 1234")
end

include("./main.jl")

# start running the server
run(ARGS[1], parse(Int64, ARGS[2]))


text:
x = BigInt(0)
for i = 0:999
    for j = 0:999
        for k = 0:99
            x += BigInt(i) * j * k
        end
    end
end
print(x)

max_stars_repo_name: vagoff/shootout

max_issues_repo_name: vagoff/shootout

max_forks_repo_name: vagoff/shootout


text:
#=
Data for modelling: Structured in Trials (pokes are concatenated to remove gaps)
    - arrival_times = vec : arrival time in the foraging port
    - departure_times = vec : leaving time in the foraging port (moment they start travel poke)
    - average_reward_rate = float64 : average reward rate fixed value
    - block_transitions = vec : block transition from long to short travel or viceversa
    - firts_block = string : travel block state at the session start
    - n_bins = int : total number of 100ms time n_bins
    - n_rewards = int : total number of rewards
    - reward_times = vec : when the reward are given (defined as when the reward available state occurs)
    - time_bin = float : size of the time bin
=#

#= Models Return
return bin by bin vector of
    rewrate(rho)_env
    rewrate_patch
    rew_prediction_error_env
    rew_prediction_error_patch
    p_staying
    p_leaving
=#

#= Reward rate comparison modelling
 5 parameters float randomly initialized 10 times or 20 if didn't stabilise
 alpha_environment = [0,0.1] learning rate of environment
 alpha_patch = [0, 0.2] learning rate of the patch
 beta = [2, 10] noise parameter or inverse temperature
 reset = [0, 0.5] patch env rew rate at the entrance of a trials, because it has decayed in the last trial
 bias = [2, 10] bias towards staying in the softmax

initialization
environment_reward_rate set to the whole session reward rate
patch_reward_rate set to the reset value

bin by bin
    check reward status 1||0 (has a reward)

    compute Env reward prediction error 0||1 - environment reward rate (evolves during travel)
    follow update with relative alpha

    Are we travelling?
        No
            compute Patch reward prediction error 0||1 - patch reward rate (reset after travel, substitute patch_rew_rate with rest param)
            follow update with relative alpha
            calculate probability of staying
            p_staying likelihood = exp(bias + beta*patch_rew_rate)/exp(bias + beta*patch_rew_rate) + exp(beta*env_rew_rate)
            p_leaving = 1 - p_staying
        Yes

        is it the first time bin?
            Yes
            compute Patch reward prediction error 0||1 - patch reward rate (reset after travel, substitute patch_rew_rate with rest param)
            follow update with relative alpha
            calculate probability of leaving
            p_leaving = exp(beta*env)/(exp(beta*env) + exp(bias + beta*patch))
            OR as before 1 - staying
            p_staying likelihood = exp(bias + beta*patch_rew_rate)/exp(bias + beta*patch_rew_rate) + exp(beta*env_rew_rate)
            p_leaving = 1 - p_staying
            1 - p_staying = exp(beta*env)/(exp(beta*env) + exp(bias + beta*patch))
            NO
            reset the patch value to RESET parameter
=#

#= Fixed Treshold
    threshold replacing the environmen_rew_rate with a fixed fitted value
    p_staying = exp(bias+ beta*patch_rew_rate)/(exp(bias+ beta*patch_rew_rate) + exp(beta*threshold)

    Alternatively 2 threshold depending on the travel block
    which is why we have a block transition vector
=#

#= Semiconstant environment / Staggered model: update environment reward rate only at end of a trial
Equivalent to Reward rate comparison model but
    - Store the env_rew_rate at the start of a trial and use as reference throught the trial
    - Keep updating env_rew_rate underhood as before
    if arrival to a new patch is true update the env_rew_rate
    CAREFUL THAT IF YOU UPDATE AT LEAVING FIRST CHECK P_LEAVING THEN CHANGE ENV_RERATE
=#


text:
# This file is part of the FWT.jl package; by Héctor Andrade Loarca

module FWT

using Images
using Wavelets
using PyPlot

export WT,filt_gen, resize_image,load_image, imageplot, rescale,
       plot_wavelet, subsampling, upsampling, cconvol, reverse, clamp,
       perform_wavortho_transform, snr

include("filters.jl")
include("operations.jl")
include("imaging.jl")

#include("")

end # module


text:
abstract type Memory <: LavaAbstraction end

offset(memory::Memory) = 0

abstract type DenseMemory <: Memory end

"""
Memory domains:
- `MEMORY_DOMAIN_HOST` is host-visible memory, ideal for uploads to the GPU. It is preferably coherent and non-cached.
- `MEMORY_DOMAIN_HOST_CACHED` is host-visible, cached memory, ideal for readbacks from the GPU. It is preferably coherent.
- `MEMORY_DOMAIN_DEVICE` is device-local. It may be visible (integrated GPUs).
"""
@enum MemoryDomain::Int8 begin
    MEMORY_DOMAIN_HOST
    MEMORY_DOMAIN_HOST_CACHED
    MEMORY_DOMAIN_DEVICE
    MEMORY_DOMAIN_DEVICE_HOST
end

struct MemoryBlock <: DenseMemory
    handle::Vk.DeviceMemory
    size::Int
    properties::Vk.MemoryPropertyFlag
    domain::MemoryDomain
    is_bound::Ref{Bool}
    ptr::Ref{Ptr{Cvoid}}
end

vk_handle_type(::Type{MemoryBlock}) = Vk.DeviceMemory

Base.size(block::MemoryBlock) = block.size

properties(block::MemoryBlock) = block.properties

ismapped(block::MemoryBlock) = block.ptr[] ≠ C_NULL
isbound(block::MemoryBlock) = block.is_bound[]

function Base.map(memory::DenseMemory, size::Integer = size(memory), offset::Integer = offset(memory))
    if Vk.MEMORY_PROPERTY_HOST_COHERENT_BIT ∉ properties(memory)
        unwrap(Vk.invalidate_mapped_memory_ranges(device(memory), [Vk.MappedMemoryRange(C_NULL, memory, offset, size)]))
    end
    ptr = unwrap(Vk.map_memory(device(memory), memory, offset, size))
    memory.ptr[] = ptr
    ptr
end

function unmap(memory::DenseMemory)
    if Vk.MEMORY_PROPERTY_HOST_COHERENT_BIT ∉ properties(memory)
        unwrap(Vk.flush_mapped_memory_ranges(device(memory), [Vk.MappedMemoryRange(C_NULL, memory, offset(memory), size(memory))]))
    end
    memory.ptr[] = C_NULL
    Vk.unmap_memory(device(memory), memory)
end

function Base.map(f, memory::DenseMemory)
    if ismapped(memory)
        f(memory.ptr[])
    else
        ptr = map(memory)
        ret = f(ptr)
        unmap(memory)
        ret
    end
end

MemoryBlock(device, size::Integer, type::Integer, domain::MemoryDomain) = unwrap(memory_block(device, size, type, domain))

function memory_block(device, size, type, domain)
    prop, i = find_memory_type(physical_device(device), type, domain)
    next = Vk.MemoryAllocateFlagsInfo(0; flags = Vk.MEMORY_ALLOCATE_DEVICE_ADDRESS_BIT)
    @propagate_errors memory = create(MemoryBlock, device, Vk.MemoryAllocateInfo(size, i; next))
    MemoryBlock(memory, size, prop.property_flags, domain, Ref(false), Ref(C_NULL))
end

find_memory_type(physical_device, type_flag, domain::MemoryDomain) = find_memory_type(Base.Fix1(score, domain), physical_device, type_flag)

function score(domain::MemoryDomain, properties)
    @match domain begin
        &MEMORY_DOMAIN_HOST =>
            (10 * (Vk.MEMORY_PROPERTY_HOST_VISIBLE_BIT in properties)) +
            (Vk.MEMORY_PROPERTY_HOST_COHERENT_BIT in properties) -
            (Vk.MEMORY_PROPERTY_HOST_CACHED_BIT in properties)
        &MEMORY_DOMAIN_HOST_CACHED =>
            (10 * (Vk.MEMORY_PROPERTY_HOST_VISIBLE_BIT | Vk.MEMORY_PROPERTY_HOST_CACHED_BIT in properties)) +
            (Vk.MEMORY_PROPERTY_HOST_COHERENT_BIT in properties)
        &MEMORY_DOMAIN_DEVICE =>
            (Vk.MEMORY_PROPERTY_DEVICE_LOCAL_BIT in properties) - 2 * (Vk.MEMORY_PROPERTY_HOST_VISIBLE_BIT in properties)
        &MEMORY_DOMAIN_DEVICE_HOST =>
            (Vk.MEMORY_PROPERTY_DEVICE_LOCAL_BIT in properties) + 2 * (Vk.MEMORY_PROPERTY_HOST_VISIBLE_BIT in properties)
    end
end

function find_memory_type(f, physical_device, type_flag)
    mem_props = Vk.get_physical_device_memory_properties(physical_device)
    memtypes = mem_props.memory_types[1:mem_props.memory_type_count]
    candidate_indices = findall(i -> type_flag & (1 << i) ≠ 0, 0:length(memtypes) - 1)
    index = argmax(i -> f(memtypes[i].property_flags), candidate_indices)
    memtypes[index], index - 1
end

function find_memory_type(physical_device, type_flag, properties::Vk.MemoryPropertyFlag)
    mem_props = Vk.get_physical_device_memory_properties(physical_device)
    indices = findall(x -> properties in x.property_flags, mem_props.memory_types[1:mem_props.memory_type_count]) .- 1
    if isempty(indices)
        error("Could not find memory with properties $properties")
    else
        ind = findfirst(i -> type_flag & (1 << i) ≠ 0, indices)
        if isnothing(ind)
            error("Could not find memory with type $type_flag")
        else
            indices[ind], mem_props[indices[ind]]
        end
    end
end

struct SubMemory{M<:DenseMemory} <: DenseMemory
    memory::M
    offset::Int
    size::Int
end

function SubMemory(memory, size; offset = 0)
    size > 0 || error("Size must be positive")
    offset + size > Lava.size(memory) || error("A SubMemory cannot extend beyond its underlying memory")
    SubMemory(memory, offset, size)
end

handle(sub::SubMemory) = handle(sub.memory)

offset(sub::SubMemory) = sub.offset

Base.size(sub::SubMemory) = sub.size

function Base.view(memory::DenseMemory, range::UnitRange)
    range_check(memory, range)
    SubMemory(memory, range.start, range.stop - range.start)
end

range_check(domain, range) = range.stop ≤ size(domain) || throw(BoundsError(domain, range))

Base.firstindex(memory::Memory) = offset(memory)
Base.lastindex(memory::Memory) = size(memory)

"Memory that can't be accessed by the renderer."
struct OpaqueMemory <: Memory end

function Base.collect(memory::MemoryBlock, size::Integer, device)
    Vk.MEMORY_PROPERTY_HOST_VISIBLE_BIT in properties(memory) && return collect(memory, size)
    device::Device
    src = BufferBlock(device, size; usage = Vk.BUFFER_USAGE_TRANSFER_SRC_BIT)
    bind!(src, memory)

    reqs = Vk.get_buffer_memory_requirements(device, src)
    @assert reqs.size ≤ memory.size
    dst = BufferBlock(device, size; usage = Vk.BUFFER_USAGE_TRANSFER_DST_BIT)
    allocate!(dst, MEMORY_DOMAIN_HOST)
    wait(transfer(device, src, dst; signal_fence = true, free_src = true))
    collect(dst)
end

function Base.collect(memory::MemoryBlock, size::Integer = size(memory))
    @assert Vk.MEMORY_PROPERTY_HOST_VISIBLE_BIT in properties(memory)
    map(memory) do mapped
        ptr = Libc.malloc(size)
        @ccall memmove(ptr::Ptr{Cvoid}, mapped::Ptr{Cvoid}, size::Csize_t)::Ptr{Cvoid}
        Base.unsafe_wrap(Array, Base.unsafe_convert(Ptr{UInt8}, ptr), (size,); own = true)
    end
end

function Base.show(io::IO, block::MemoryBlock)
    print(io, MemoryBlock, "($(Base.format_bytes(size(block))), $(block.domain))")
end


text:
using SimpleDirectMediaLayer.LibSDL2
using Test

include("lib/SDL.jl")
include("lib/SDL_image.jl")
include("lib/SDL_ttf.jl")
include("lib/SDL_mixer.jl")


text:
"""
Test two-stage shooting method
"""

using LinearAlgebra
using DifferentialEquations
using BenchmarkTools
using Plots

pyplot()

include("../../julia-R3BP/R3BP/src/R3BP.jl")
include("../src/ShootingStar.jl")

## Initialize ODE settings
reltol = 1.e-12
abstol = 1.e-12
method = Tsit5()

## define parameters
params = R3BP.get_cr3bp_param(399, 301)
mu = params.mu
println("mu: $mu")
lp = R3BP.lagrangePoints(mu)

# initial condition of halo
X0 = [1.176924090973164, 0.0, -0.060210863312217, 0.0, -0.173836346247689, 0.0]

T = 3.385326412831325

X2 = [1.176924090973164, 0.0, -0.060210863312217, 0.0, -0.173836346247689, 0.0] + 1.e-12*rand(6)

tspan = (0.0, 0.5T)
prob = ODEProblem(R3BP.rhs_cr3bp_sv!, X0, tspan, (mu))
sol = solve(prob, method, reltol=reltol, abstol=abstol)
X1 = sol.u[end] + 1.e-12*rand(6)  #+ [0.0, 0.0,-1.e-5, 0.0, -1.e-6, 0.0]

## prepare for two-stage shooting
n_sv = 6
tofs = [T/2, T/2]
svs = [X0, X1, X2]
n = length(svs)

prob_stm = ODEProblem(R3BP.rhs_cr3bp_svstm!, vcat(X0, reshape(I(6), (36,)))[:], T, (mu))
nodes_conv = ShootingStar.twostage_shooting(n_sv, svs, tofs, prob_stm)

# propagate result
prob = ODEProblem(R3BP.rhs_cr3bp_sv!, X0, T, (mu))
function prob_func(prob, i, repeat)
    remake(prob, u0=nodes_conv[i], tspan=(0.0, tofs[i]))
end

ensemble_prob = EnsembleProblem(prob, prob_func=prob_func)
sim = solve(ensemble_prob, method, EnsembleThreads(), trajectories=n-1, reltol=reltol, abstol=abstol);

display(plot(sim, linealpha=0.4, vars=(1,2), flip=false, aspect_ratio=:equal, size=(800,650), c=:orangered,
     frame_style=:box, gridalpha=0.4, xlabel="x", ylabel="y"))
println("Done!")


text:
readCSV(fn) = DataFrame(CSV.File(fn))

function catalyst_name(m::ReactionSystem;fcall=Catalyst.species,remove_t=true)
    spec_names_t = string.(fcall(m))
    spec_names = remove_t ? replace.(spec_names_t, "(t)"=>"") : spec_names_t
    return spec_names
end

function getTables(csv_fns::NamedTuple)
    dfs = []
    for fn in csv_fns
        push!(dfs, readCSV(fn))
    end
    return NamedTuple{keys(csv_fns)}(tuple(dfs...))
end

"""
Get constructor from a data type
"""
construct(datatype) = typeof(datatype).name.wrapper

function findID(sNames, features; convert_lower = true)
    sNames = convert_lower ? lowercase.(sNames) : sNames
    features = convert_lower ? lowercase.(features) : features
    idxs = Vector{Int}()
    for (ind, i) in enumerate(sNames)
        isEle = true
        for f in features
            if occursin(f, string(i)) == false
                isEle = false
                break
            end
        end
        isEle ? push!(idxs, ind) : nothing
    end
    return idxs
end

"""
Scale maximum element to `capVal`. `v` is a vector
"""
function capVec(v, capVal)
    i = argmax(v)
    scale = capVal/v[i]
    V_cap = v .* scale
    return V_cap
end


text:
using Test
using Graphs

@testset "Incidence Dictionaries" begin

g = Graphs.incdict(Graphs.ExVertex,is_directed=false)
variables = [ExVertex(1, "x1"), ExVertex(2, "x2"), ExVertex(3, "x3"), ExVertex(4, "x4")]
map(v -> add_vertex!(g, v), variables)

for i = 1 : length(variables)
    @test vertices(g)[i] in variables
    @test out_degree(variables[i], g) == 0
end

factors =[ExVertex(5, "x1x2f1"), ExVertex(6, "x2x3f1"), ExVertex(7, "x3x4f1")]
add_vertex!(g, factors[1])
edge = Graphs.make_edge(g, variables[1], factors[1]); Graphs.add_edge!(g, edge)
edge = Graphs.make_edge(g, factors[1], variables[2]); Graphs.add_edge!(g, edge)
add_vertex!(g, factors[2])
edge = Graphs.make_edge(g, variables[2], factors[2]); Graphs.add_edge!(g, edge)
edge = Graphs.make_edge(g, factors[2], variables[3]); Graphs.add_edge!(g, edge)
add_vertex!(g, factors[3])
edge = Graphs.make_edge(g, variables[3], factors[3]); Graphs.add_edge!(g, edge)
edge = Graphs.make_edge(g, factors[3], variables[4]); Graphs.add_edge!(g, edge)

@test num_vertices(g) == 7

@test [out_degree(v, g) for v in variables] == [1, 2, 2, 1]
@test [out_degree(f, g) for f in factors] == [2, 2, 2]

@test setdiff(collect(out_neighbors(variables[1], g)), [factors[1]]) == []
@test setdiff(collect(out_neighbors(variables[2], g)), [factors[1], factors[2]]) == []
@test setdiff(collect(out_neighbors(variables[3], g)), [factors[2], factors[3]]) == []
@test setdiff(collect(out_neighbors(variables[4], g)), [factors[3]]) == []
@test setdiff(collect(in_neighbors(variables[1], g)), [factors[1]]) == []
@test setdiff(collect(in_neighbors(variables[2], g)), [factors[1], factors[2]]) == []
@test setdiff(collect(in_neighbors(variables[3], g)), [factors[2], factors[3]]) == []
@test collect(in_neighbors(variables[4], g)) == [factors[3]]

# Delete testing
delete_vertex!(variables[4], g)
@test num_vertices(g) == 6
@test_throws Exception out_degree(variables[4], g)
@test out_degree(factors[3], g) == 1 # Should be just x3 now

end