text:
@safetestset hi_g_definition = "g is the direction of gray code" begin
using BijectiveHilbert: brgc, hi_g
for i in 0:1000
    @test brgc(i) ⊻ brgc(i + 1) == 1<<hi_g(i)
    @test brgc(i) == brgc(i + 1) ⊻ 1<<hi_g(i)  # xor. It's how it works.
end
end

@safetestset binary_reflection_of_gc = "gray code is binary reflected" begin
using BijectiveHilbert: brgc
n = 10
for i in 0:(1<<n - 1)
    @test brgc(1<<n - 1 - i) == brgc(i) ⊻ 1<<(n-1)
end
end

@safetestset hi_g_matches_long_form = "hi_g is the same as its definition" begin
using BijectiveHilbert: hi_g, hi_g_orig
for i = 0:100
    @test hi_g(i) == hi_g_orig(i)
end
end

@safetestset hi_g_symmetry = "hi_g symmetry invariant" begin
using BijectiveHilbert: brgc, hi_g, bshow
# n = 3
# for i in 0:(1<<n)
#     println(
#         bshow(brgc(i)), " ",
#         bshow(brgc(i + 1)), " ",
#         bshow(brgc(i) ⊻ brgc(i + 1)), " ",
#         bshow(1<<hi_g(i)), " ",
#         bshow(i)
#     )
# end

# Test for the hi_g function.
n = 5
for i in 0:(1<<n - 2)
    # println(hi_g(i), " ", hi_g(1<<n - 2 - i))
    @test(hi_g(i) == hi_g(1<<n - 2 - i))
end
end

@safetestset efdg_match_table = "directions match figure 5 in hamilton report" begin
using BijectiveHilbert: hi_e, hi_f, hi_d, hi_g
n = 2
trials = [
    #i     e     f  d  g
    [0, 0b00, 0b01, 0, 0],
    [1, 0b00, 0b10, 1, 1],
    [2, 0b00, 0b10, 1, 0],
    [3, 0b11, 0b10, 0, 2]
]
for (i, e, f, d, g) in trials
    @test hi_e(i) == Int(e)
    @test hi_f(i, n) == Int(f)
    @test hi_d(i, n) == d
    if g != 42
        @test hi_g(i) == g
    end
end
end

# The 2d example in the paper has a simpler 1d equivalent,
# which is two segments in a row. d=0 for both. f is 1.
@safetestset efdg_match_1d = "directions match a one-dimensional equivalent" begin
using BijectiveHilbert: hi_e, hi_f, hi_d, hi_g
n = 1
trials = [
    #i  e  f  d  g
    [0, 0, 1, 0, 0],  # the first segment should match n=2
    [1, 0, 1, 0, 1]  # the second segment will be redirected for n=2
]
for (i, e, f, d, g) in trials
    @test hi_e(i) == e
    @test hi_f(i, n) == f
    @test hi_d(i, n) == d
    @test hi_g(i) == g
end
end

# A 3d example.
@safetestset efdg_match_3d = "directions match a three-dimensional equivalent" begin
using BijectiveHilbert: hi_e, hi_f, hi_d, hi_g
n = 3
trials = [
    #i      e      f  d  g
    [0, 0b000, 0b001, 0, 0], # segments 0-2 match n=2
    [1, 0b000, 0b010, 1, 1],
    [2, 0b000, 0b010, 1, 0],
    [3, 0b011, 0b111, 2, 2], # segment 3 gets redirected in z-direction
    [4, 0b011, 0b111, 2, 0], # segment 4 also in z direction
    [5, 0b110, 0b100, 1, 1], # then we repeat the 2d plane in reverse.
    [6, 0b110, 0b100, 1, 0],
    [7, 0b101, 0b100, 0, 3]
]
for (i, e, f, d, g) in trials
    @test hi_e(i) == Int(e)
    @test hi_f(i, n) == Int(f)
    @test hi_d(i, n) == d
    @test hi_g(i) == g
end
end

@safetestset hi_e_matches_float_version = "hi_e matches its floating point version" begin
using BijectiveHilbert: hi_e, hi_e_orig
# Show our implementation for hi_e matches the specification.
for i = 1:100
    @test(hi_e(i) == hi_e_orig(i))
end
end

@safetestset cubes_neighbor_along_g_coord = "cubes are neighbors along gth coordinate" begin
# From page 12, just below the figure. The definition of the directions.
using BijectiveHilbert: hi_e, hi_f, hi_g
for n in 2:5
    # Goes to 2^n-2, not 2^n-1 because there is no neighbor for the last point.
    for i in 0:(1<<n - 2)
        fside = hi_f(i, n) ⊻ (1<<hi_g(i))
        eside = hi_e(i + 1)
        @test fside == eside
        if fside != eside
            @show n, i, fside, eside
            break
        end
    end
end
end

@safetestset hi_d_symmetry_invariant = "hi_d symmetry invariant works corollary 2.7" begin
using BijectiveHilbert: hi_d, bshow
# n = 3
# for i = 0:(1<<n)
#     println(i, " ", bshow(hi_d(i, n)))
# end

# invariant for d on page 12. Corollary 2.7.
# Does not hold true for i=0 and i=1<<n - 1.
n = 5
for i = 0:(1<<n - 1)
    #println(hi_d(i, n), " ", hi_d((1<<n)-1-i, n))
    d = hi_d(i, n)
    dflip = hi_d((1<<n) - 1 - i, n)
    @test(d == dflip)
end
end

@safetestset hi_d_compare_g = "hi_d equals hi_g at 2n-1 lemma 2.8" begin
using BijectiveHilbert: hi_d, hi_g
for n = 1:5
    d = hi_d(1<<n - 1, n)
    g = hi_g(1<<n - 1)
    @test d % n == g % n
end
end

@safetestset hi_d_compares_with_g = "hi_d equals g for mod 2" begin
using BijectiveHilbert: hi_d, hi_g
for n = 2:5
    for i = 0:(1<<n - 1)
        d = hi_d(i, n)
        if i == 0
            @test d == 0
        elseif i % 2 == 0
            g = hi_g(i - 1)
            @test d == (g % n)
        else
            g = hi_g(i)
            @test d == (g % n)
        end
    end
end
end

@safetestset hi_edg_invariant = "hi e, d, g invariant is consistent" begin
using BijectiveHilbert: hi_d, hi_e, hi_g
# invariant for e, d, and g. pg. 11. Eq. 1.
# Contrary to the tech report, this is NOT true for 2^n-1.
n = 5
for i = 0:(1<<n - 2)
    @test(hi_e(i + 1) == hi_e(i) ⊻ (1 << hi_d(i, n)) ⊻ (1 << hi_g(i)))
end
end

@safetestset hi_e_reflection = "hi_e reflects to f" begin
using BijectiveHilbert: hi_d, hi_e, hi_f
n = 5
for i = 0:(1<<n - 1)
    e = hi_e(i)
    d = hi_d(i, n)
    f = hi_f(i, n)
    # Corollary 2.7, pg 12, states that e ⊻ d = f, but that's not true.
    # Page 11, first paragraph, says e(i) \xor f(i) = 2^d(i). That's this.
    @test e ⊻ (1<<d) == f
end
end

@safetestset rotation_invariant_zero = "rotation invariant on hi_e is zero" begin
using BijectiveHilbert: hi_T, hi_e, hi_d
# lemma 2.11, page 15
# Assert T_{e,d}(e) == 0
n = 3
for i = 0:(1<<n - 1)
    # println(hi_T(hi_e(i), hi_d(i, n), hi_e(i), n))
    @test(0 == hi_T(hi_e(i), hi_d(i, n), hi_e(i), n))
end
end

# "rotation invariant on hi_f is power of two" begin
# lemma 2.11, page 15
# Assert T_{e,d}(f) == 2^(n-1)
# Contrary to the tech report, this is NOT true.
# @test(0x1<<(n-1) == hi_T(hi_f(i, n), hi_d(i, n), hi_e(i), n))

# Top of page 16, stated invariant:
# (T_(e,d)(a) rotleft (d + 1)) ⊻ e == a
# Also NOT true

@safetestset ef_invariant = "ef relationship holds" begin
using BijectiveHilbert: hi_e, hi_f
# invariant bottom of pg. 11. Lemma 2.6.
# fails for i=0 and i=2^n
n = 5
for i = 0b1:(0b1<<(n - 1))
    e = hi_e(i)
    ff = hi_f((0b1<<n) - 0b1 - i, n) ⊻ (0b1<<(n-1))
    # println(e, " ", ff)
    @test(e == ff)
    @test typeof(e) == UInt8
    @test typeof(ff) == UInt8
end
end

@safetestset inverse_t_itself = "transform is its own inverse" begin
using BijectiveHilbert: hi_d, hi_e, hi_T, hi_T_inv
# Check that the inverse of T is an inverse.
n = 3
for i = 0x0:(0x1<<n - 0x1)
    for b = 0x0:(0x1<<n - 0x1)
        @test(typeof(i) == UInt8)
        d = UInt8(hi_d(i, n))
        e = UInt8(hi_e(i))
        a = hi_T(b, d, e, n)
        b1 = hi_T_inv(a, d, e, n)
        # println(join(string.(typeof.((i, b, a, b1))), " "))
        # println("b ", bitstring(b), " a ", bitstring(a), " b1 ", bitstring(b1))
        @test(b == b1)
    end
end
end

@safetestset get_ith_bit_trials = "trials show can get ith bit of vectors" begin
using BijectiveHilbert: ith_bit_of_indices, get_location

v = [0b01110100, 0b01101010, 0b01011001]
trials = [
    [0, "00000100"],
    [1, "00000010"],
    [2, "00000001"],
    [3, "00000110"],
    [4, "00000101"],
    [5, "00000011"],
    [6, "00000111"],
    [7, "00000000"]
]
for (idx, t0) in trials
    @test bitstring(ith_bit_of_indices(3, v, idx)) == t0
    # Look, Hamilton's version matches.
    @test bitstring(get_location(v, idx)) == t0
end
end

@safetestset ith_bit_types = "ith_bit ops use all types" begin
using BijectiveHilbert: ith_bit_of_indices, set_indices_bits!
using Random
rng = MersenneTwister(4972433234)
tt = [Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128]
for T in tt
    for j in 1:10
    n = 1
    p = zeros(T, 1)
    p[1] = T(rand(rng, 0:typemax(T)))
    q = zeros(T, 1)
    for i in 0:(8*sizeof(T) - 1)
        b = ith_bit_of_indices(n, p, i)
        set_indices_bits!(q, b, n, i)
    end
    @test q[1] == p[1]
    end
end
end

@safetestset hilbert_index_complete4d = "hilbert index is a complete set for 4d" begin
using BijectiveHilbert: SpaceGray
using ..HilbertTestSuite: check_complete_set
b = 3
n = 4
gg = SpaceGray(b, n)
@test check_complete_set(gg, b, n)
end

# @safetestset hilbert_index_paper_tr = "paper and techreport agree" begin
# using BijectiveHilbert: hilbert_index, hilbert_index_paper
# They do NOT agree.
# m = 0x4  # One larger won't fit into a byte and must fail.
# for i in 0:(1<<m - 1)
#     for j in 0:(1<<m - 1)
#         h_tr = hilbert_index(0x2, m, UInt8[i, j])
#         h_p = hilbert_index_paper(0x2, m, UInt8[i, j])
#         @test typeof(h_tr) == typeof(h_p)
#         @test h_p == h_tr
#         if h_p != h_tr
#             break
#         end
#     end
# end
# end

@safetestset libhilbert_matches = "libhilbert output matches" begin
using BijectiveHilbert: SpaceGray, encode_hilbert_zero
n = 4
b = 5
gg = SpaceGray(b, n)
trials = [
    [10, 1, 1, 1, 1],
    [11, 1, 0, 1, 1],
    [12, 1, 0, 1, 0],
    [13, 1, 1, 1, 0],
    [14, 1, 1, 0, 0],
    [15, 1, 0, 0, 0],
    [16, 2, 0, 0, 0],
    [17, 2, 0, 1, 0],
    [18, 2, 0, 1, 1],
    [19, 2, 0, 0, 1]
]
for trial in trials
    v = convert(Vector{UInt8}, trial[2:end])
    result = Int(encode_hilbert_zero(gg, v))
    @test result == trial[1]
end
end

@safetestset correct_bits_set_in_indices = "indices bits match" begin
using BijectiveHilbert: set_indices_bits!
m = 3
i = 4
l = 0b110
p = zeros(Int, 3)
set_indices_bits!(p, l, m, i)
@test p == [0, 1<<i, 1<<i]
end

@safetestset paper_inv = "paper hilbert is own inverse" begin
using BijectiveHilbert: SpaceGray
using ..HilbertTestSuite: check_own_inverse
for n in 2:5
    for b in 2:4
        gg = SpaceGray(b, n)
        @test check_own_inverse(gg, b, n)
    end
end
end

@safetestset brgc_mask_equal = "mask calculation matches paper" begin
using BijectiveHilbert: extract_mask, extract_mask_paper, hi_d
using Random
rng = MersenneTwister(974073242)
for trial in 1:1000
    n = rand(rng, 2:7)
    m = rand(rng, 2:5, n)
    i = rand(rng, 0:(n - 1))
    d = hi_d(i, n)
    a = extract_mask(m, n, d, i)
    b = extract_mask_paper(m, n, d, i)
    @test a == b
end
end

@safetestset libhilbert_hc_matches = "libhilbert compressed output matches" begin
using BijectiveHilbert: encode_hilbert_zero, Compact
n = 4
b = 5
trials = [
    [10, 1, 1, 1, 1],
    [11, 1, 0, 1, 1],
    [12, 1, 0, 1, 0],
    [13, 1, 1, 1, 0],
    [14, 1, 1, 0, 0],
    [15, 1, 0, 0, 0],
    [16, 2, 0, 0, 0],
    [17, 2, 0, 1, 0],
    [18, 2, 0, 1, 1],
    [19, 2, 0, 0, 1]
]
for trial in trials
    v = convert(Vector{UInt8}, trial[2:end])
    ms = [b, b, b, b]
    gg = Compact(ms)
    result = Int(encode_hilbert_zero(gg, v))
    @test result == trial[1]
end
end

@safetestset compact_index_is_its_inverse = "compact index is its inverse" begin
using BijectiveHilbert: Compact
using ..HilbertTestSuite: check_own_inverse
using Random
rng = MersenneTwister(432479874)
for n = [5]
    for i in 1:1
        ms = rand(rng, 2:5, n)
        gg = Compact(ms)
        @test check_own_inverse(gg, ms, n)
    end
end
end

@safetestset compact_index_matches_examples = "compact index matches examples" begin
include(joinpath(dirname(@__FILE__), "check_hamilton.jl"))
using BijectiveHilbert: Compact, decode_hilbert_zero!, index_type, axis_type
fns = ["test/3_2_5_7.txt", "test/2_5_3.txt", "test/3_4_2_3.txt",
       "test/3_5_2_3.txt", "test/4_5_5_5.txt"]
for fn in fns
    if !isfile(fn)
        fn = basename(fn)
    end
    if isfile(fn)
        h, X, ms = hamilton_example(fn)
        n = size(X, 1)
        gg = Compact(ms)
        p = zeros(axis_type(gg), n)
        H = index_type(gg)
        for i in eachindex(h)
            decode_hilbert_zero!(gg, p, H(h[i]))
            @test convert(Vector{UInt8}, p) == X[:, i]
        end
    end
end
end

@safetestset compact_index_inv_matches_examples = "inverse compact index matches examples" begin
include(joinpath(dirname(@__FILE__), "check_hamilton.jl"))
using BijectiveHilbert: Compact, encode_hilbert_zero, axis_type, index_type
fns = ["test/3_2_5_7.txt", "test/2_5_3.txt", "test/3_4_2_3.txt",
       "test/3_5_2_3.txt", "test/4_5_5_5.txt"]
for fn in fns
    if !isfile(fn)
        fn = basename(fn)
    end
    if isfile(fn)
        h, X, ms = hamilton_example(fn)
        n = size(X, 1)
        gg = Compact(ms)
        H = index_type(gg)
        for i in eachindex(h)
            XX = convert(Vector{axis_type(gg)}, X[:, i])
            hc = encode_hilbert_zero(gg, XX)
            @test hc == H(h[i])
            if hc != H(h[i])
                @show n, ms, i, h[i], hc, X[:, i]
                break
            end
        end
    end
end
end

@safetestset compact_for_signed = "compact hilbert for signed integers" begin
using BijectiveHilbert
using Random
rng = MersenneTwister(9742439)
for i in 1:1000
    n = rand(rng, 3:5)
    ms = rand(rng, 2:7, n)
    gg = Compact(Int, ms)
    gg2 = Compact(ms)
    x = zeros(Int, n)
    A = axis_type(gg2)
    y = zeros(A, n)
    for i in 1:n
        x[i] = rand(rng, 0:(1<<ms[i] - 1))
        y[i] = A(x[i])
    end
    a = encode_hilbert_zero(gg, x)
    b = Int(encode_hilbert_zero(gg2, y))
    @test a == b
end
end

@safetestset hamilton_zero_based = "hamilton works for zero based" begin
using BijectiveHilbert
gg = Compact(UInt, [4, 3, 3])
X = [9, 4, 5]
h = encode_hilbert(gg, X)
Y = zeros(Int, 3)
decode_hilbert!(gg, Y, h)
@test X == Y
end

@safetestset spacegray_type_interactions = "SpaceGray type interactions" begin
    using BijectiveHilbert
    using UnitTestDesign
    using Random
    rng = Random.MersenneTwister(9790323)
    for retrial in 1:5
        AxisTypes = shuffle(rng, [Int8, Int, UInt, Int128, UInt8, UInt128])
        IndexTypes = shuffle(rng, [Union{}, Int8, UInt8, Int, UInt, Int128, UInt128])
        Count= shuffle(rng, [0, 1])
        Dims = shuffle(rng, [2, 3, 4])
        Bits = shuffle(rng, [2, 3, 4, 5])
        test_set = all_pairs(
            AxisTypes, IndexTypes, Count, Dims, Bits;
        )
        for (A, I, C, D, B) in test_set
            if I == Union{}
                gg = SpaceGray(B, D)
                I = index_type(gg)
            else
                gg = SpaceGray(I, B, D)
            end
            if B * D > log2(typemax(I))
                continue
            end
            last = (one(I) << (B * D)) - one(I) + I(C)
            mid = one(I) << (B * D - 1)
            few = 5
            X = zeros(A, D)
            hlarr = vcat(C:min(mid, few), max(mid + 1, last - few):last)
            for hl in hlarr
                hli = I(hl)
                if C == 0
                    decode_hilbert_zero!(gg, X, hli)
                    hl2 = encode_hilbert_zero(gg, X)
                    if hl2 != hli
                        @show A, I, C, D, B, X
                        @test hl2 == hli
                    end
                    @test typeof(hl2) == typeof(hli)
                else
                    decode_hilbert!(gg, X, hli)
                    hl3 = encode_hilbert(gg, X)
                    @test hl3 == hli
                    @test typeof(hl3) == typeof(hli)
                end
            end
        end
    end
end

@safetestset spacegray_random = "spacegray random" begin
using BijectiveHilbert
aas = [Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128]
tts = [UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128]
ggbase = SpaceGray(UInt32, 4, 4)
hli = UInt32(2)
Xbase = zeros(UInt8, 4)
decode_hilbert_zero!(ggbase, Xbase, hli)

for idx in Base.IteratorsMD.CartesianIndices((length(aas), length(tts)))
    i, j = Tuple(idx)
    A = aas[i]
    T = tts[j]
    gg = SpaceGray(T, 4, 4)
    X = zeros(A, 4)
    hli = T(2)
    decode_hilbert_zero!(gg, X, hli)
    @test X == Xbase
    hl2 = encode_hilbert_zero(gg, X)
    @test hl2 == hli
end
end

@safetestset compact_type_interactions = "Compact type interactions" begin
    using BijectiveHilbert
    using UnitTestDesign
    using Random
    rng = Random.MersenneTwister(9790323)
    for retrial in 1:5
        AxisTypes = shuffle(rng, [Int8, Int, UInt, Int128, UInt8, UInt128])
        IndexTypes = shuffle(rng, [Union{}, Int8, UInt8, Int, UInt, Int128, UInt128])
        Count= shuffle(rng, [0, 1])
        Dims = shuffle(rng, [2, 3, 4])
        Bits = shuffle(rng, [2, 3, 4, 5])
        test_set = all_pairs(
            AxisTypes, IndexTypes, Count, Dims, Bits;
        )
        for (A, I, C, D, B) in test_set
            if I == Union{}
                gg = Compact(fill(B, D))
                I = index_type(gg)
            else
                gg = Compact(I, fill(B, D))
            end
            if B * D > log2(typemax(I))
                continue
            end
            last = (one(I) << (B * D)) - one(I) + I(C)
            mid = one(I) << (B * D - 1)
            few = 5
            X = zeros(A, D)
            hlarr = vcat(C:min(mid, few), max(mid + 1, last - few):last)
            for hl in hlarr
                hli = I(hl)
                if C == 0
                    decode_hilbert_zero!(gg, X, hli)
                    hl2 = encode_hilbert_zero(gg, X)
                    if hl2 != hli
                        @show A, I, C, D, B, X
                        @test hl2 == hli
                    end
                    @test typeof(hl2) == typeof(hli)
                else
                    decode_hilbert!(gg, X, hli)
                    hl3 = encode_hilbert(gg, X)
                    @test hl3 == hli
                    @test typeof(hl3) == typeof(hli)
                end
            end
        end
    end
end


text:
require("Ito")

using Calendar
using Ito.Time
using Ito.TermStructure

@test_approx_eq 1.7 compound_factor(0.7, :Compounded, Ito.TermStructure.Annual, Actual360(), today(), today()+days(360))
@test_approx_eq 1.35^2 compound_factor(0.7, :Compounded, Ito.TermStructure.Semiannual, Actual360(), today(), today()+days(360))

@test_approx_eq 1/1.7 discount_factor(0.7, :Compounded, Ito.TermStructure.Annual, Actual360(), today(), today()+days(360))

ext: jl


text:
module HNScraper

    export scrape_topstories
    export get_topstories

    using JSON
    using HTTP
    using DataFrames
    using Dates
    using Feather

    include("functions.jl")

end

max_stars_repo_name: JohannesNakayama/hacker-news-analysis

max_issues_repo_name: JohannesNakayama/hacker-news-analysis

max_forks_repo_name: JohannesNakayama/hacker-news-analysis


text:
module SpectralEnvelope

using FFTW
using Polynomials
using LinearAlgebra
using Statistics

function detrend(x,y; deg = 1)
    fit = polyfit(x,y,deg)
    detrended = y - fit(x)
    return detrended
end

"""
smoothens the given series by mixen each points with it's neighboors.

    input:
    - data: the series you want to smooth
    - m : the numbers of points to be involved in the mixing process.

    output:
    - smoothened series, shorter of m points than the original.

"""

function smooth(data::Array{Float64,1}; m = 5)
    if m != 0
        smoothed_data = Float64[]
        weigths = Int[]
        for w in collect(-m:m)
            p = abs(2*w)+1
            if p >= 10
                append!(weigths,p)
            else
                append!(weigths,10)
            end
        end
        p_tot = sum([1/w for w in weigths])
        for i in 1:length(data)
            if i > m && i < length(data) - m
                smoothed_Iw = 0
                for (l,k) in enumerate(collect(-m:m))
                    smoothed_Iw += (1/weigths[l])*data[i+k]
                end
                append!(smoothed_data, smoothed_Iw/p_tot)
            end
        end
        return smoothed_data
    elseif m == 0
        return data
    end
end

"""
One-hot encodes the time-series. For k categories,
each category gets a k-1 unit vector, except the last one which
is mapped to zeros(k - 1).

    Example:

    vectorize([1,2,3,1,2]) returns :
    [[1,0],[0,1],[0,0],[1,0]]

"""

function vectorize(data)
    categories = unique(data)
    deleteat!(categories, length(categories))
    sorted_data = zeros(length(categories),length(data))
    for (t,d) in enumerate(data)
        for (i,j) in enumerate(categories)
            if d == j
                sorted_data[i,t] = 1
            end
        end
    end
    return sorted_data
end

"""
Splits the time-series into overlapping equal-lengthed chunks of given size.

    input :
    - x : the time-series
    - window : the size of the equal-lengthed chunks
    - step : size of the step by which the window is sliding along the time-series
    the smaller the value of size, the bigger the overlap between the different chunks.

    output :
    -  overlapping equal-lengthed windowed time-series
"""

function partitioning(x,window::Int,step::Int)
    return [x[i:i+window] for i in 1:step:length(x) if  i + window <= length(x)]
end

function partitioning(x::Array{Float64,2},window::Int,step::Int)
    return [x[:,i:i+window-1] for i in 1:step:length(x[1,:]) if  i + window <= length(x[1,:])]
end

"""
returns the power spectrum of a given time serie,

    input :
    - time serie
    - window for averaging (default = length(time serie)/10)
    - step for the window's sliding

    output :
    - value of the power spectrum

You will have to provide the frequencies yourself
"""

function power_spectrum(x::Array{Float64,1}, window::Int, step::Int)
    ts = partitioning(x,window,step)
    ps = [real(fft(i .- mean(i)).*conj(fft(i .- mean(i)))) for i in ts]
    pxx = mean(ps)
    return [pxx[i] for i in 1:div(window,2)]
end

"""
    periodogram_matrix(ts::Array{Float64,2}, demean = false; m=2)

Computes the k x k periodogram matrix (k number of categories).
If 'average' true, segments data and computes the average periodogram matrix over all segments.
m is the smoothing parameter.
"""
function periodogram_matrix(ts::Array{Float64,2}, average = false; m=2)
    if average == false
        periodo = zeros(length(ts[:,1]), length(ts[:,1]), length(ts[1,:]))
        DFT = zeros(Complex{Float64}, length(ts[:,1]), length(ts[1,:]))
        for i in 1:length(ts[:,1])
            DFT[i,:] = fft(ts[i,:] .- mean(ts[i,:]))
        end
        for i in 1:length(ts[:,1])
            for j in 1:length(ts[:,1])
                tmp = smooth(real((1/length(ts[1,:]))*DFT[i,:].*conj(DFT[j,:])); m = m)
                for w in 1:length(ts[1,:])-(2*m+1)
                    periodo[i,j,w] = tmp[w]
                end
            end
        end
        return periodo
    elseif average == true
        window = Int(floor(length(ts[1,:])/2))
        step = Int(floor(length(ts[1,:])/4))
        dims = length(ts[:,1])
        periodo = zeros(length(ts[:,1]), length(ts[:,1]), window)
        seg_ts = partitioning(ts, window, step)
        DFT = []
        for el in seg_ts
            tmp = zeros(Complex, dims, window)
            for i in 1:length(el[:,1])
                tmp[i,:] = fft(el[i,:] .- mean(el[i,:]) )
            end
            push!(DFT, tmp)
        end
        for el in DFT
            for i in 1:dims
                for j in 1:dims
                    tmp = smooth(real((1/window)*el[i,:].*conj(el[j,:])); m = m)
                    for w in 1:window-(2*m+1)
                        periodo[i,j,w] += tmp[w]
                    end
                end
            end
        end
        return periodo./length(seg_ts)
    end
end

"""
Computes the covariance-variance matrix of a given time-series.
"""
function varcov(ts::Array{Float64,2})
    cov_matrix = zeros(length(ts[:,1]), length(ts[:,1]))
    for i in 1:length(ts[:,1])
        for j in 1:length(ts[:,1])
            cov_matrix[i,j] = sum((1/(length(ts[1,:])-1))*(ts[i,:].-mean(ts[i,:])).*(ts[j,:].-mean(ts[j,:])))
        end
    end
    return cov_matrix
end

"""
Computes the spectral envelope of the given time-series.

    input:

    - ts: the time series to analyse
    - m : the degree of smoothing wished.
          m corresponds to the number of neighbooring points that are mixed
          with given point to realize the smoothing.

    output:

    - Frequencies : list of points corresponding to the involved frequencies.
                    contained in [0,0.5]
    - amplitude : values of the spectral envelope for each given frequency point.
    - eigenvectors : the optimal scaling for the different categories, for each frequency point.
    - categories : the list of categories present in the data.

"""
function spectral_envelope(ts; m = 3)
    result = Float64[]
    vec = vectorize(ts)
    eigvec = zeros(length(vec[:,1]), trunc(Int,length(ts)/2))
    S = inv(sqrt(varcov(vec)))
    # alternative definition of S
    # S = inv(varcov(vec))
    p = periodogram_matrix(vec; m = m)
    f_len = trunc(Int,length(p[1,1,:])/2)
    for i in 1:f_len
        ev = findmax(real.(eigvals(S*p[:,:,i]*S/f_len)))
        # ev = findmax(real.(eigvals(S*p[:,:,i]/f_len)))
        append!(result, ev[1])
        b = S*eigvecs(S*p[:,:,i])[ev[2],:]
        b = b./sqrt(sum(b.^2))
        eigvec[:,i] = real.(b)
    end
    return collect(1:length(result))/length(p[1,1,:]), result[1:end], eigvec
end

"""
     get_mappings(data, freq; m = 3)

Returns the optimal mappings corresponding to the frequency 'freq'.
Prints the position and strengh of peak at 'f' for control purposes.
input :
    - data : input categorical time series
    - freq : the frequency where one wants the optimal mappings
    - m : smoothing parameters (see 'spectral_envelope')
returns :
    - mappings : the mapping at the peak associated with the desired frequency.

example:
    # we have some data, plotting the spectral envelope, we see a peak at 0.33
    m = get_mappings(data, 0.33)
"""
function get_mappings(data, freq; m = 3)
    categories = unique(data)
    f, se, eigvecs = spectral_envelope(data; m = m)
    delta = f[2] - f[1]
    window = Int(div(0.04*length(f),1))
    peak_pos = findmin([abs(freq - delta*i) for i in 1:length(f)])[2]
    true_pos = findmax(se[peak_pos-window:peak_pos+window])[2] + peak_pos - window - 1
    mappings = Dict(categories[i] => round(eigvecs[i,true_pos]; digits = 2) for i in 1:length(eigvecs[:,1]))
    if length(categories) != length(eigvecs[:,1])
        push!(mappings, unique(data)[end] => 0.0)
    end
    println("position of peak: ",round(f[true_pos],digits = 2)," strengh of peak: ",round(se[true_pos], digits = 2))
    return mappings
end

function findmax_in(xserie,yserie,xlim)
    range = []
    real_pos = []
    for (index,value) in enumerate(xserie)
        if value>= xlim[1] && value <= xlim[2]
            append!(range,yserie[index])
            append!(real_pos,index)
        end
    end
    max,pos = findmax(range)
    return max, xserie[real_pos[pos]], real_pos[pos]
end

"""
    apply_mapping(ts, mapping)

Transforms the input time-series 'ts' according to the mapping 'mapping'.
Typically, 'mapping' would be the output of the 'get_mappings' function, although you can provide your own.
'mapping' should be a Dict specifying each value and what it should be mappied to. Example:
{"A" : 0.54, "G" : 0.62, "T" : -0.57, "C" : 0.0}

Returns an array which correspond to the mapped series according to 'mapping'
"""
function apply_mapping(ts, mapping)
    mapped_series = deepcopy(ts)
    for (index, value) in enumerate(mapped_series)
        mapped_series[index] = mapping[value]
    end
    return mapped_series
end

export spectral_envelope, get_mappings, detrend, smooth, power_spectrum, apply_mapping

end


text:
module test_8
using AdventOfCode2018.Day8
using Test

test_input = "2 3 0 3 10 11 12 1 1 0 1 99 2 1 1 2\n"
numbers  = reverse([2, 3, 0, 3, 10, 11, 12, 1, 1, 0, 1, 99, 2, 1, 1, 2])
tree = Node(metadata=[1, 1, 2], children=[
    Node(metadata=[10, 11, 12]),
    Node(metadata=[2], children=[Node(metadata=[99])])
])

@testset "day 8, part 1" begin
	@test parse_input(test_input) == numbers
    @test build_tree(numbers) == tree
    @test part_1(tree) == 138
    @test part_2(tree) == 66

end

end


text:
using Toolips, Documenter
import Toolips: Servable, ServerExtension, ToolipsServer
Documenter.makedocs(root = "../",
       source = "src",
       build = "build",
       clean = true,
       doctest = true,
       modules = [Toolips],
       repo = "https://github.com/ChifiSource/Toolips.jl",
       highlightsig = true,
       sitename = "toolips",
       expandfirst = [],
       pages = [
                "toolips" => "index.md",
                "projects" => "projects.md",
                "core" => "core.md",
                "servables" => "servables.md",
                "developer api" => "developer_api.md",
               ]
       )


text:
################################################################################
#  Copyright 2022, Tom Van Acker (BASF)                                        #
################################################################################
# TFT.jl                                                                       #
# A Julia package for Taylor-Fourier Transform.                                #
# See http://github.com/timmyfaraday/TFT.jl                                    #
################################################################################

module TFT

    # import pkg
    import DSP
    import Polynomials
    import FFTW

    # pkg constants
    const _DSP  = DSP
    const _POL  = Polynomials
    const _FFTW = FFTW

    # paths
    const BASE_DIR = dirname(@__DIR__)

    # include
    include("types/dtft.jl")

    include("core/estimator.jl")
    include("core/ospline.jl")

    include("prob/tft.jl")
    include("prob/ftft.jl")

    include("util/util.jl")

    # export
    export  tft, ftft

    export  amplitude, phase, ar_phase, frequency, rocof, phasor, ar_phasor,
            signal, error
    export  a, ϕ, φ, f, r, ξ, ψ

end


text:
# finite mixture models

####
#
#  All subtypes of AbstractMixtureModel should implement the following methods:
#
#  - ncomponents(d): the number of components
#
#  - component(d, k):  return the k-th component
#
#  - probs(d):       return a vector of prior probabilities over components.
#

abstract AbstractMixtureModel{VF<:VariateForm,VS<:ValueSupport,C<:Distribution} <: Distribution{VF, VS}

immutable MixtureModel{VF<:VariateForm,VS<:ValueSupport,C<:Distribution} <: AbstractMixtureModel{VF,VS,C}
    components::Vector{C}
    prior::Categorical

    function MixtureModel(cs::Vector{C}, pri::Categorical)
        length(cs) == ncategories(pri) ||
            error("The number of components does not match the length of prior.")
        new(cs, pri)
    end
end

typealias UnivariateMixture{S<:ValueSupport,   C<:Distribution} AbstractMixtureModel{Univariate,S,C}
typealias MultivariateMixture{S<:ValueSupport, C<:Distribution} AbstractMixtureModel{Multivariate,S,C}
typealias MatrixvariateMixture{S<:ValueSupport,C<:Distribution} AbstractMixtureModel{Matrixvariate,S,C}

component_type{VF,VS,C}(d::AbstractMixtureModel{VF,VS,C}) = C

#### Constructors

function MixtureModel{C<:Distribution}(components::Vector{C}, prior::Categorical)
    VF = variate_form(C)
    VS = value_support(C)
    MixtureModel{VF,VS,C}(components, prior)
end

MixtureModel{C<:Distribution}(components::Vector{C}, p::Vector{Float64}) =
    MixtureModel(components, Categorical(p))

# all components have the same prior probabilities
MixtureModel{C<:Distribution}(components::Vector{C}) =
    MixtureModel(components, Categorical(length(components)))

_construct_component{C<:Distribution}(::Type{C}, arg) = C(arg)
_construct_component{C<:Distribution}(::Type{C}, args::Tuple) = C(args...)

function MixtureModel{C<:Distribution}(::Type{C}, params::AbstractArray, p::Vector{Float64})
    components = C[_construct_component(C, a) for a in params]
    MixtureModel(components, p)
end

function MixtureModel{C<:Distribution}(::Type{C}, params::AbstractArray)
    components = C[_construct_component(C, a) for a in params]
    MixtureModel(components)
end

#### Basic properties

length(d::MultivariateMixture) = length(d.components[1])
size(d::MatrixvariateMixture) = size(d.components[1])

ncomponents(d::MixtureModel) = length(d.components)
components(d::MixtureModel) = d.components
component(d::MixtureModel, k::Int) = d.components[k]

probs(d::MixtureModel) = probs(d.prior)
params(d::MixtureModel) = ([params(c) for c in d.components], params(d.prior)[1])

function mean(d::UnivariateMixture)
    K = ncomponents(d)
    p = probs(d)
    m = 0.0
    for i = 1:K
        pi = p[i]
        if pi > 0.0
            c = component(d, i)
            m += mean(c) * pi
        end
    end
    return m
end

function mean(d::MultivariateMixture)
    K = ncomponents(d)
    p = probs(d)
    m = zeros(length(d))
    for i = 1:K
        pi = p[i]
        if pi > 0.0
            c = component(d, i)
            BLAS.axpy!(pi, mean(c), m)
        end
    end
    return m
end

function var(d::UnivariateMixture)
    K = ncomponents(d)
    p = probs(d)
    means = Array(Float64, K)
    m = 0.0
    v = 0.0
    for i = 1:K
        pi = p[i]
        if pi > 0.0
            ci = component(d, i)
            means[i] = mi = mean(ci)
            m += pi * mi
            v += pi * var(ci)
        end
    end
    for i = 1:K
        pi = p[i]
        if pi > 0.0
            v += pi * abs2(means[i] - m)
        end
    end
    return v
end

function cov(d::MultivariateMixture)
    K = ncomponents(d)
    p = probs(d)
    m = zeros(length(d))
    md = zeros(length(d))
    V = zeros(length(d),length(d))

    for i = 1:K
        pi = p[i]
        if pi > 0.0
            c = component(d, i)
            BLAS.axpy!(pi, mean(c), m)
            BLAS.axpy!(pi, cov(c), V)
        end
    end
    for i = 1:K
        pi = p[i]
        if pi > 0.0
            c = component(d, i)
            # todo: use more in-place operations
            md = mean(c) - m
            BLAS.axpy!(pi, md*md', V)
        end
    end
    return V
end

#### show

function show(io::IO, d::MixtureModel)
    K = ncomponents(d)
    pr = probs(d)
    println(io, "MixtureModel{$(component_type(d))}(K = $K)")
    Ks = min(K, 8)
    for i = 1:Ks
        @printf(io, "components[%d] (prior = %.4f): ", i, pr[i])
        println(io, component(d, i))
    end
    if Ks < K
        println(io, "The rest are omitted ...")
    end
end

#### Evaluation

function _cdf(d::UnivariateMixture, x::Real)
    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K
    r = 0.0
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0
            c = component(d, i)
            r += pi * cdf(c, x)
        end
    end
    return r
end

cdf(d::UnivariateMixture{Continuous}, x::Float64) = _cdf(d, x)
cdf(d::UnivariateMixture{Discrete}, x::Int) = _cdf(d, x)

function _mixpdf1(d::AbstractMixtureModel, x)
    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K
    v = 0.0
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0
            c = component(d, i)
            v += pdf(c, x) * pi
        end
    end
    return v
end

function _mixpdf!(r::AbstractArray, d::AbstractMixtureModel, x)
    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K
    fill!(r, 0.0)
    t = Array(Float64, size(r))
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0
            pdf!(t, component(d, i), x)
            BLAS.axpy!(pi, t, r)
        end
    end
    return r
end

function _mixlogpdf1(d::AbstractMixtureModel, x)
    # using the formula below for numerical stability
    #
    # logpdf(d, x) = log(sum_i pri[i] * pdf(cs[i], x))
    #              = log(sum_i pri[i] * exp(logpdf(cs[i], x)))
    #              = log(sum_i exp(logpri[i] + logpdf(cs[i], x)))
    #              = m + log(sum_i exp(logpri[i] + logpdf(cs[i], x) - m))
    #
    #  m is chosen to be the maximum of logpri[i] + logpdf(cs[i], x)
    #  such that the argument of exp is in a reasonable range
    #

    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K

    lp = Array(Float64, K)
    m = -Inf   # m <- the maximum of log(p(cs[i], x)) + log(pri[i])
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0
            # lp[i] <- log(p(cs[i], x)) + log(pri[i])
            lp_i = logpdf(component(d, i), x) + log(pi)
            @inbounds lp[i] = lp_i
            if lp_i > m
                m = lp_i
            end
        end
    end
    v = 0.0
    @inbounds for i = 1:K
        if p[i] > 0.0
            v += exp(lp[i] - m)
        end
    end
    return m + log(v)
end

function _mixlogpdf!(r::AbstractArray, d::AbstractMixtureModel, x)
    K = ncomponents(d)
    p = probs(d)
    @assert length(p) == K
    n = length(r)
    Lp = Array(Float64, n, K)
    m = fill(-Inf, n)
    for i = 1:K
        @inbounds pi = p[i]
        if pi > 0.0
            lpri = log(pi)
            lp_i = slice(Lp, :, i)
            # compute logpdf in batch and store
            logpdf!(lp_i, component(d, i), x)

            # in the mean time, add log(prior) to lp and
            # update the maximum for each sample
            for j = 1:n
                lp_i[j] += lpri
                if lp_i[j] > m[j]
                    m[j] = lp_i[j]
                end
            end
        end
    end

    fill!(r, 0.0)
    @inbounds for i = 1:K
        if p[i] > 0.0
            lp_i = slice(Lp, :, i)
            for j = 1:n
                r[j] += exp(lp_i[j] - m[j])
            end
        end
    end

    @inbounds for j = 1:n
        r[j] = log(r[j]) + m[j]
    end
    return r
end

pdf(d::UnivariateMixture{Continuous}, x::Real) = _mixpdf1(d, x)
pdf(d::UnivariateMixture{Discrete}, x::Int) = _mixpdf1(d, x)
logpdf(d::UnivariateMixture{Continuous}, x::Real) = _mixlogpdf1(d, x)
logpdf(d::UnivariateMixture{Discrete}, x::Int) = _mixlogpdf1(d, x)

_pdf!(r::AbstractArray, d::UnivariateMixture{Discrete}, x::UnitRange) = _mixpdf!(r, d, x)
_pdf!(r::AbstractArray, d::UnivariateMixture, x::AbstractArray) = _mixpdf!(r, d, x)
_logpdf!(r::AbstractArray, d::UnivariateMixture, x::AbstractArray) = _mixlogpdf!(r, d, x)

_pdf(d::MultivariateMixture, x::AbstractVector) = _mixpdf1(d, x)
_logpdf(d::MultivariateMixture, x::AbstractVector) = _mixlogpdf1(d, x)
_pdf!(r::AbstractArray, d::MultivariateMixture, x::AbstractMatrix) = _mixpdf!(r, d, x)
_lodpdf!(r::AbstractArray, d::MultivariateMixture, x::AbstractMatrix) = _mixlogpdf!(r, d, x)

## component-wise pdf and logpdf

function _cwise_pdf1!(r::AbstractVector, d::AbstractMixtureModel, x)
    K = ncomponents(d)
    length(r) == K || error("The length of r should match the number of components.")
    for i = 1:K
        r[i] = pdf(component(d, i), x)
    end
    r
end

function _cwise_logpdf1!(r::AbstractVector, d::AbstractMixtureModel, x)
    K = ncomponents(d)
    length(r) == K || error("The length of r should match the number of components.")
    for i = 1:K
        r[i] = logpdf(component(d, i), x)
    end
    r
end

function _cwise_pdf!(r::AbstractMatrix, d::AbstractMixtureModel, X)
    K = ncomponents(d)
    n = size(X, ndims(X))
    size(r) == (n, K) || error("The size of r is incorrect.")
    for i = 1:K
        pdf!(slice(r,:,i), component(d, i), X)
    end
    r
end

function _cwise_logpdf!(r::AbstractMatrix, d::AbstractMixtureModel, X)
    K = ncomponents(d)
    n = size(X, ndims(X))
    size(r) == (n, K) || error("The size of r is incorrect.")
    for i = 1:K
        logpdf!(slice(r,:,i), component(d, i), X)
    end
    r
end

componentwise_pdf!(r::AbstractVector, d::UnivariateMixture, x::Real) = _cwise_pdf1!(r, d, x)
componentwise_pdf!(r::AbstractVector, d::MultivariateMixture, x::AbstractVector) = _cwise_pdf1!(r, d, x)
componentwise_pdf!(r::AbstractMatrix, d::UnivariateMixture, x::AbstractVector) = _cwise_pdf!(r, d, x)
componentwise_pdf!(r::AbstractMatrix, d::MultivariateMixture, x::AbstractMatrix) = _cwise_pdf!(r, d, x)

componentwise_logpdf!(r::AbstractVector, d::UnivariateMixture, x::Real) = _cwise_logpdf1!(r, d, x)
componentwise_logpdf!(r::AbstractVector, d::MultivariateMixture, x::AbstractVector) = _cwise_logpdf1!(r, d, x)
componentwise_logpdf!(r::AbstractMatrix, d::UnivariateMixture, x::AbstractVector) = _cwise_logpdf!(r, d, x)
componentwise_logpdf!(r::AbstractMatrix, d::MultivariateMixture, x::AbstractMatrix) = _cwise_logpdf!(r, d, x)

componentwise_pdf(d::UnivariateMixture, x::Real) = componentwise_pdf!(Array(Float64, ncomponents(d)), d, x)
componentwise_pdf(d::UnivariateMixture, x::AbstractVector) = componentwise_pdf!(Array(Float64, length(x), ncomponents(d)), d, x)
componentwise_pdf(d::MultivariateMixture, x::AbstractVector) = componentwise_pdf!(Array(Float64, ncomponents(d)), d, x)
componentwise_pdf(d::MultivariateMixture, x::AbstractMatrix) = componentwise_pdf!(Array(Float64, size(x,2), ncomponents(d)), d, x)

componentwise_logpdf(d::UnivariateMixture, x::Real) = componentwise_logpdf!(Array(Float64, ncomponents(d)), d, x)
componentwise_logpdf(d::UnivariateMixture, x::AbstractVector) = componentwise_logpdf!(Array(Float64, length(x), ncomponents(d)), d, x)
componentwise_logpdf(d::MultivariateMixture, x::AbstractVector) = componentwise_logpdf!(Array(Float64, ncomponents(d)), d, x)
componentwise_logpdf(d::MultivariateMixture, x::AbstractMatrix) = componentwise_logpdf!(Array(Float64, size(x,2), ncomponents(d)), d, x)

## Sampling

immutable MixtureSampler{VF,VS,Sampler} <: Sampleable{VF,VS}
    csamplers::Vector{Sampler}
    psampler::AliasTable
end

function MixtureSampler{VF,VS}(d::MixtureModel{VF,VS})
    csamplers = map(sampler, d.components)
    psampler = sampler(d.prior)
    MixtureSampler{VF,VS,eltype(csamplers)}(csamplers, psampler)
end

rand(d::MixtureModel) = rand(component(d, rand(d.prior)))

rand(s::MixtureSampler) = rand(s.csamplers[rand(s.psampler)])
_rand!(s::MixtureSampler{Multivariate}, x::AbstractVector) = _rand!(s.csamplers[rand(s.psampler)], x)

sampler(d::MixtureModel) = MixtureSampler(d)


text:
using BinaryProvider

JLCXX_DIR = get(ENV, "JLCXX_DIR", "")
if(isempty(JLCXX_DIR))
    JLCXX_DIR =  get(ENV, "JlCxx_DIR", "")
end

const verbose = "--verbose" in ARGS
const prefix = Prefix(JLCXX_DIR == "" ? get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")) : JLCXX_DIR)

products = Product[
    LibraryProduct(prefix, "libcxxwrap_julia", :libcxxwrap_julia),
    LibraryProduct(prefix, "libcxxwrap_julia_stl", :libcxxwrap_julia_stl)
]

supported = true

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaBinaryWrappers/libcxxwrap_julia_jll.jl/releases/download/libcxxwrap_julia-v0.6.6%2B0/"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Windows(:i686, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/libcxxwrap_julia.v0.6.6.i686-w64-mingw32-cxx11.tar.gz", "cc5503ad6bad6e7dc3e17ab6e95f08a807d3959dcb634fb5e33be99162fa3245"),
    MacOS(:x86_64) => ("$bin_prefix/libcxxwrap_julia.v0.6.6.x86_64-apple-darwin14-cxx11.tar.gz", "d62b2065e4c0dcaefc0612d7c00d6d2c1dd96601dc87dd24730df71c06c0fcc2"),
    Linux(:x86_64, libc=:glibc, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/libcxxwrap_julia.v0.6.6.x86_64-linux-gnu-cxx11.tar.gz", "bbd1e5fe3fcb9b5c9ee26091b18abe4eb38c0875a96c63a63feef0c09b66cd54"),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7, :cxx11)) => ("$bin_prefix/libcxxwrap_julia.v0.6.6.x86_64-w64-mingw32-cxx11.tar.gz", "c764ba83621a12196e5824414d5205d62b3803674fa8be4107edcaaa10065973"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
transform_platform(platform) = typeof(platform)(platform.arch;libc=platform.libc,call_abi=platform.call_abi,compiler_abi=CompilerABI(:gcc7,:cxx11))
transform_platform(platform::MacOS) = MacOS(:x86_64)
if JLCXX_DIR == ""
    platform = transform_platform(platform_key_abi())
    if haskey(download_info, platform)
        if !supported
            error("Julia version $VERSION is not supported for binary download. Please build libcxxwrap-julia from source and set the JLCXX_DIR environment variable to the build dir or installation prefix.")
        end
        url, tarball_hash = download_info[platform]
        if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
            # Download and install binaries
            install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose,ignore_platform=true)
        end
    elseif unsatisfied
        # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
        # Alternatively, you could attempt to install from a separate provider,
        # build from source or something even more ambitious here.
        error("Your platform $(triplet(platform)) is not supported by this package!")
    end
else
    if unsatisfied
        error("The required libraries were not found in the provided JLCXX_DIR directory $JLCXX_DIR")
    end
end

write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
# An implementation of the Playfair cipher in Julia

# Loops an index if it is greater than 5
function loop5(i)
    if i > 5
        i -= 5
    end
    return i
end

# The main function
function playfair()

    # A keyword is needed to create the square
    println("Please enter a keyword")
    keyword = readline()
    keyword = replace(keyword, " " => "")
    keyword = uppercase(keyword)

    alphabet = "ABCDEFGHIKLMNOPQRSTUVWXYZ" # No j

    # Substitute all j's for i's
    keyword = replace(keyword, "J" => "I")

    # Create a string representation of the square
    keyString = ""
    for c in keyword
        if !in(c, keyString)
            keyString = string(keyString, string(c))
        end
    end
    for c in alphabet
        if !in(c, keyString)
            keyString = string(keyString, string(c))
        end
    end

    # Fill 5x5 grid
    str1 = split(SubString(keyString, 1, 5), "")
    str2 = split(SubString(keyString, 6, 10), "")
    str3 = split(SubString(keyString, 11, 15), "")
    str4 = split(SubString(keyString, 16, 20), "")
    str5 = split(SubString(keyString, 21, 25), "")
    keygrid = [str1 str2 str3 str4 str5]

    # Needs rotating the right way
    grid = Array{String}(undef, 5, 5)
    for c in keygrid
        i = findfirst(isequal(c), keygrid)
        x = i[2]
        y = i[1]
        grid[x,y] = c
    end
    #show(stdout, "text/plain", grid)

    # Take the plaintext
    println("Please enter the plaintext")
    plaintext = readline()
    plaintext = uppercase(plaintext)

    # Replace J's with I's
    plaintext = replace(plaintext, "J" => "I")

    # Remove spaces and turn into bigrams
    plaintext = replace(plaintext, " " => "")

    # Add X on the end if the text is of odd length
    if isodd(length(plaintext))
        plaintext = plaintext * "X"
    end

    # Turn into bigrams
    bigrams = []
    token = 0
    for i in 1:length(plaintext)
        if iseven(i + token)
            bigram = plaintext[i-1] * plaintext[i]
            if bigram[1] == bigram[2]
                bigram = bigram[1] * "X"
                token += 1
            end
            push!(bigrams, bigram)
        end
    end

    # Replace repeated letters with X
    for i in eachindex(bigrams)
        if bigrams[i][1] == bigrams[i][2]
            bigrams[i] = bigrams[i][1] * "X"
        end
    end

    # Use the grid to encrypt the plaintext
    ciphertext = ""
    for b in bigrams
        pos1 = findfirst(isequal(string(b[1])), grid)
        pos2 = findfirst(isequal(string(b[2])), grid)
        if pos1[2] == pos2[2] # Same X
            ciphertext *= grid[loop5(pos1[1] + 1), pos1[2]] * grid[loop5(pos2[1] + 1), pos2[2]] * " "
        elseif pos1[1] == pos2[1] # Same Y
            ciphertext *= grid[pos1[1], loop5(pos1[2] + 1)] * grid[pos2[1], loop5(pos2[2] + 1)] * " "
        else # Must form a square
            ciphertext *= grid[pos1[1], pos2[2]] * grid[pos2[1], pos1[2]] * " "
        end
    end
    println(ciphertext)
end

playfair()


text:
#!/usr/bin/env julia6
using ReadGlobal, MyStats, JLD2

function dbkhistindices(input::AbstractString,nbins::Integer)
    field = readfield(input)
    minf,maxf = min_max(field)
    ind = dbkhist_indices(field,nbins)
    centers =  zeros(nbins)
    @inbounds for i in 1:nbins
        centers[i] = mean(view(field,ind[i]))
    end
    bins = centers
    @save "pdf_dbkindices_$(input)_nbins$(nbins).jld2" ind bins
end

function main(f,ns)
    n = parse(Int,ns)
    dbkhistindices(f,n)
end

function main(f)
    n = 36
    dbkhistindices(f,n)
end

main(ARGS...)


text:
module Regularization

using Flux.Losses: logitcrossentropy
using Statistics

"""
	spectral_decoupling(ŷ, y, λ)

Logitcrossentropy extended with spectral decoupling regularization as proposed by
[Pezeshki et al.](https://arxiv.org/abs/2011.09468). The implementation here somewhat differs from the formula presented
in the paper but the methodology is consistent with the python code provided by the authors.
"""
function spectral_decoupling(ŷ, y, λ)
	return logitcrossentropy(ŷ, y) + λ * mean(ŷ .^ 2)
end

end


text:
function obs_lh_given_hmm(observations, hmm; linear=true)
    linear ? (return linear_likelihood(observations, hmm)) : (return bw_likelihood(Matrix(transpose(observations)), hmm))
end

function bw_likelihood(observations, hmm)
    obs_lengths = [findfirst(iszero,observations[:,o])-1 for o in 1:size(observations,2)]

    lls = bw_llhs(hmm, observations)
    log_α = messages_forwards_log(hmm.a, hmm.A, lls, obs_lengths)
    log_β = messages_backwards_log(hmm.A, lls, obs_lengths)

    O = size(lls)[3] #the last T value is the 0 end marker of the longest T
    log_pobs = zeros(O)
    Threads.@threads for o in 1:O
        log_pobs[o] = logsumexp(lps.(log_α[:,1,o], log_β[:,1,o]))
    end

    return sum(log_pobs)
end

function linear_likelihood(observations,hmm)
    O= size(observations,1);
    obs_lengths = [findfirst(iszero,observations[o,:])-1 for o in 1:size(observations,1)] #mask calculations here rather than mle_step to prevent recalculation every iterate

    a = transpose(log.(hmm.a)); A = log.(hmm.A)
    N = length(hmm.B)
    mask=observations.!=0
    #INITIALIZATION
    βoi_T = zeros(O,N); βoi_t = zeros(O,N) #log betas at T initialised as zeros

    #RECURRENCE
    βoi_T = backwards_lh_sweep!(hmm, A, N, βoi_T, βoi_t, observations, mask, obs_lengths)

    #TERMINATION
    lls = c_llhs(hmm,observations[:,1])
    α1om = lls .+ a #first position forward msgs

    return lps([logsumexp(lps.(α1om[o,:], βoi_T[o,:])) for o in 1:O])
end
                #LINEAR_STEP SUBFUNCS
                function backwards_lh_sweep!(hmm, A, N, βoi_T, βoi_t, observations, mask, obs_lengths)
                    for t in maximum(obs_lengths)-1:-1:1
                        last_β=copy(βoi_T)
                        lls = c_llhs(hmm,observations[:,t+1])
                        omask = findall(mask[:,t+1])
                        βoi_T[omask,:] .+= view(lls,omask,:)
                        Threads.@threads for m in 1:N
                            βoi_t[omask,m] = logsumexp.(eachrow(view(βoi_T,omask,:).+transpose(view(A,m,:))))
                        end
                        βoi_T=copy(βoi_t)
                    end
                    return βoi_T
                end


text:
module InnerLoaders

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export innerloader

register_normal_element("q__inner__loading", context = @__MODULE__)

function innerloader(label::String = "",
                          state::Union{Symbol,Nothing} = nothing,
                          args...;
                          wrap::Function = StippleUI.DEFAULT_WRAPPER,
                          kwargs...)

  wrap() do
    q__inner__loading(  args...;
              attributes(
                        [ Symbol(":state") => state,
                          :label => label,
                          kwargs...
                        ], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

end


text:
#=
Suppose you are given a table of currency exchange rates, represented as a 2D array. Determine whether there is a possible arbitrage: that is, whether there is some sequence of trades you can make, starting with some amount A of any currency, so that you can end up with some amount greater than A of that currency.

There are no transaction costs and you can trade fractional quantities.
#=

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
function getelemface(dim,elemtype)

if dim== 1
    face=[1;2];
    nfe = 2;
    nvf = 1;
elseif dim==2
    if elemtype==0
        face=[[2,3];[3,1];[1,2]];
        nfe = 3;
        nvf = 2;
    elseif elemtype==1
        face=[[1,2];[2,3];[3,4];[4,1]];
        nfe = 4;
        nvf = 2;
    end
elseif dim==3
    if elemtype==0
        face=[[2,3,4];[1,4,3];[1,2,4];[1,3,2]];
        nfe = 4;
        nvf = 3;
    elseif elemtype==1
        face=[[1,4,3,2];[5,6,7,8];[1,2,6,5];[3,4,8,7];[2,3,7,6];[4,1,5,8]];
        nfe = 6;
        nvf = 4;
    end
else
    error("Only can handle dim=1, dim=2 or dim=3");
end

face = reshape(face,nvf,nfe);

return face;

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
###################################################
## geo/timings/fromnetwork.jl
## extends NetworkTimings object to continuous timings, using a projector
###################################################

"""
    `TimingsFromNetwork`, uses projector to extend timings from network
"""
mutable struct TimingsFromNetwork <: GeoTimings
    trips::GeoData # compulsory attribute

    "timings in network"
    timings::NetworkTimings
    "Projector that links network timings to continuous timings"
    proj::NetworkProjector

    function TimingsFromNetwork(nt::NetworkTimings, proj::NetworkProjector)
        gt = new()
        gt.trips = proj.trips
        gt.timings = nt
        gt.proj  = proj
        return gt
    end
end

estimateTime(gt::TimingsFromNetwork, pLon::Float32, pLat::Float32, dLon::Float32, dLat::Float32)=
    getTripTiming(gt.proj, gt.timings, pLon, pLat, dLon, dLat)

estimateTime(gt::TimingsFromNetwork, tId::Int) =
    getTripTiming(gt.proj, gt.timings, tId)


text:
# These are extra tests + some tyre kicking with the ordering
# following CommonMark.jl's blocks and some extras
#   ✅ means there's some tests,
#   🚫 means it's not supported because ambiguous with something else
#
# 0 paragraphs ✅
# 1 atxheading ✅
# 2 blockquote ✅
# 3 fencedcodeblock ✅
# x htmlblock 🚫
# x indentedcodeblock 🚫
# 4 list ✅ (⚠ validation done in Franklin)
# x setextheading 🚫
# 5 hrule ✅
# 6 paragraphs --> 0 ✅
# 7 emphasis *, **, ***, _, __, ___  ✅
# 8 autolink <..> ✅ (⚠ normalisation via URIs is done in Franklin)
# 9 htmlentity ✅ (they're left as is)
# x htmlinline 🚫
# 10 inlinecode ✅
# 11 image, links, footnotes ✅ (⚠ no check that ref exists)
#
# x hard line breaks 🚫
# 12 comments ✅
# 13 backslash escapes ✅
#
# 14 table blocks ✅ (⚠ validation done in Franklin)
#
# -- Franklin
#
# f0 raw
# f1 inline math ✅ (including switchoff)
# f2 block math ✅
# f3 code block ✅
# f4 code block with lang ✅
# f5 code block eval ✅ (see 3)
# f6 newcom ✅ (⚠ assembly done in Franklin, needs the def)
# f7 com ✅
# f7i internal coms
# f8 newenv
# f9 env
# f9i internal envs (e.g. eqs)
# f10 cu brackets  ✅ (see f7 etc)
# f11 dbb ✅
# f12 emojis
# f13 def line   @def ...
# f14 def block  +++...+++
# f15 div block
# f16 html block

@testset "0>paragraphs" begin
    p = """
        abc

        def

        ghi
        """ |> grouper
    @test ct(p[1]) // "abc"
    @test ct(p[2]) // "def"
    @test ct(p[3]) // "ghi"
    @test all(p_i.role == :PARAGRAPH for p_i in p)

    p = """
        abc
        ```
        def
        ```
        ghi
        """ |> grouper
    @test ct(p[1]) // "abc"
    @test p[2].role == :CODE_BLOCK
    @test ctf(p[2]) // "def"
    @test p[3].role == :PARAGRAPH
    @test ct(p[3]) // "ghi"

    p = raw"""
        abc `def` ghi $jkl$ & 123 @@c foo@@ end
        """ |> grouper
    @test ct(p[1]) // raw"abc `def` ghi $jkl$ & 123"
    @test ctf(p[2]) // "foo"
    @test ct(p[3]) // "end"

    # nesting should not matter as inner blocks are disabled
    p = raw"""
        abc @@A aa @@B bb @@ @@ def
        """ |> grouper
    @test ct(p[1]) // "abc"
    @test ctf(p[2]) // "aa @@B bb @@"
    @test ct(p[3]) // "def"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "1>atxheading" begin
    p = """
        # a
        ## b
        ### c
        #### d
        ##### e
        ###### f
        """ |> grouper
    @test length(p) == 6
    @test ct(p[1]) // "# a"
    @test ct(p[2]) // "## b"
    @test ct(p[3]) // "### c"
    @test ct(p[4]) // "#### d"
    @test ct(p[5]) // "##### e"
    @test ct(p[6]) // "###### f"
    @test ctf(p[1]) // "a"
    @test ctf(p[2]) // "b"
    # spaces don't matter
    p = """
        # a
            ## b
        """ |> grouper
    @test ctf(p[1]) // "a"
    @test ctf(p[2]) // "b"

    # has to be at start of line
    p = """
        a # bc
        """ |> grouper
    @test ctf(p[1]) // "a # bc"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "2>blockquotes" begin
    p = """
        > abc
        """ |> grouper
    @test p[1] // "> abc"
    p = """
        > abc
        > def
        > ghi
        """ |> grouper
    @test p[1] // "> abc\n> def\n> ghi"
    p = """
        > abc
        def

        ghi
        """ |> grouper
    @test p[1] // "> abc\ndef"
    @test p[2] // "ghi"
    p = """
        > abc
        > def

        >ghi
        > jkl
        """ |> grouper
    @test p[1] // "> abc\n> def"
    @test p[2] // ">ghi\n> jkl"

    p = """
        > abc
        > > def
        > > ghi
        """ |> grouper
    @test p[1] // "> abc\n> > def\n> > ghi"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "3+f3+f4+f5>codeblock" begin
    p = """
        abc
        ```
        def
        ```
        ghi
        """ |> grouper
    @test p[1].role == :PARAGRAPH
    @test ctf(p[2]) // "def"
    @test p[3].role == :PARAGRAPH

    p = """
        abc
        ````
        def
        ```
        ghi
        ```
        jkl
        ````
        mno
        """ |> grouper
    @test p[1] // "abc"
    @test ctf(p[2]) // "def\n```\nghi\n```\njkl"
    @test p[3] // "mno"

    p = """
        abc
        ```def
        ghi
        ```
        jkl
        """ |> grouper
    @test ctf(p[2]) // "def\nghi"

    p = """
        abc
        ```!
        def
        ```
        ghi
        """ |> grouper
    @test ctf(p[2]) // "!\ndef"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "4>list" begin
    p = """
        abc
        * i
        * j
            * k
        """ |> grouper
    @test p[1] // "abc"
    @test p[2] // "* i\n* j\n    * k"

    p = """
        abc
        * i
          * j
        still part of j
          * k
        """ |> grouper
    @test p[1] // "abc"
    @test p[2] // "* i\n  * j\nstill part of j\n  * k"

    p = """
        abc
        1. i
        1) j
        1 k
        """ |> grouper
    @test p[2] // "1. i\n1) j\n1 k"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "5>hrule" begin
    p = """
        abc
        ---
        def
        """ |> grouper
    @test isp(p[1])
    @test p[2] // "---"
    @test isp(p[3])

    p = """
        > abc
        ***
        > def
        """ |> grouper
    @test p[1] // "> abc"
    @test p[2] // "***"
    @test p[3] // "> def"

    p = """
        > abc
        xxx
        ___
        > def
        """ |> grouper
    @test p[1] // "> abc\nxxx"
    @test p[2] // "___"
    @test p[3] // "> def"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "6>paragraph" begin
    # see 0>
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "7>emph" begin
    p = """
        a *b* _c_ **d** __e__ ***f*** ___g___

        **b _c_ d**
        """ |> md_blockifier
    @test p[1].name == :EMPH_EM
    @test ct(p[1]) == "b"
    @test p[2].name == :EMPH_EM
    @test ct(p[2]) == "c"
    @test p[3].name == :EMPH_STRONG
    @test ct(p[3]) == "d"
    @test p[4].name == :EMPH_STRONG
    @test ct(p[4]) == "e"
    @test p[5].name == :EMPH_EM_STRONG
    @test ct(p[5]) == "f"
    @test p[6].name == :EMPH_EM_STRONG
    @test ct(p[6]) == "g"
    @test p[end-1].name == :EMPH_STRONG
    @test ct(p[end-1]) == "b _c_ d"

    s = "a*b*c*"
    t = s |> toks
    b = s |> md_blockifier
    @test ct(b[1]) == "b"
    @test length(b) == 1
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "8>autolink" begin
    p = """
        a <bc> def <http://example.com> and < done >>.
        """ |> md_blockifier
    @test p[1].name == :AUTOLINK
    @test ct(p[1]) == "bc"
    @test p[2].name == :AUTOLINK
    @test ct(p[2]) == "http://example.com"
    @test p[3].name == :P_BREAK
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

# entities are just left 'as is'.
@testset "9>entity" begin
    p = """
        abc & def & but &amp. *
        """ |> grouper
    @test ctf(p[1]) == "abc & def & but &amp. *"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "10>inline code" begin
    s = """
        abc `def` and `` ` `` and *`fo*o`*.
        """
    p = s |> slice
    @test ct(p[1]) // "abc"
    @test ct(p[2]) // "def"  # `def`
    @test ct(p[3]) // "and"
    @test ct(p[4]) // " ` "
    @test ct(p[5]) // "and"
    @test ct(p[6]) // "`fo*o`"
    @test ct(p[7]) // "."
    g = s |> grouper
    @test length(g) == 1  # all inline blocks
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "11>img, links, footnotes" begin
    s = """
        [abc] [def](ghi) ![jkl] ![mno](pqr)
        [ref]: aaa
        """
    b = s |> md_blockifier
    @test b[1].name == :LINK_A
    @test b[1] // "[abc]"
    @test b[2].name == :LINK_AB
    @test b[2] // "[def](ghi)"
    @test b[3].name == :IMG_A
    @test b[3] // "![jkl]"
    @test b[4].name == :IMG_AB
    @test b[4] // "![mno](pqr)"
    @test b[5].name == :REF
    @test b[5] // "[ref]:"

    # not ok because not at the start of a line bar spaces
    s = "abc [def]: hello" |> md_blockifier
    @test length(s) == 0

    # aggregation over multiline for ref
    s = """
        abc
        [def]: foo
        bar

        baz
        """
    p = s |> grouper
    @test ctf(p[1]) // "abc"
    @test length(p[2].blocks) == 1
    @test p[2].blocks[1] // "[def]: foo\nbar"
    @test ctf(p[3]) // "baz"

    # AR links
    p = """
        [A][B] ![C][D] ![E][F]:
        """ |> grouper
    @test p[1].blocks[1] // "[A][B]"
    @test p[1].blocks[3] // "![C][D]"
    @test p[1].blocks[5] // "![E][F]" # precedence of the img
    @test p[1].blocks[6] // ":"

    # final img
    p = """
        [A][B] ![C]
        """ |> grouper
    @test p[1].blocks[1] // "[A][B]"
    @test p[1].blocks[3] // "![C]"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "12>comment" begin
    p = "Hello <!--bar--> baz foo `<!--aa-->`" |> grouper
    @test length(p) == 1
    @test p[1].blocks[1] // "Hello"
    @test p[1].blocks[2] // "<!--bar-->"
    @test p[1].blocks[3] // "baz foo"
    @test p[1].blocks[4] // "`<!--aa-->`"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "13>escapes" begin
    p = raw"abc \_ foo" |> slice
    @test Int('_') == 95
    @test text(p[1]) // "abc _ foo"
    for c in raw"""*_`@#{}$~!"%&'+,-./:;<=>?^|"""
        p = "a \\" * c * " b" |> slice
        ic = Int(c)
        @test text(p[1]) // "a &#$(ic); b"
    end
    p = "a \\ b" |> slice
    @test text(p[1]) // "a \ b"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "14>table" begin
    p = """
        abc
        | a | b | c |
        | - | - | - |
        | 1 | 2 | 3 |
        def
        """ |> grouper
    @test p[1] // "abc"
    @test p[3] // "def"
    p = "abc | def" |> md_blockifier
    @test length(p) == 0
    p = "| abc | def |" |> md_blockifier
    @test p[1].name == :TABLE_ROW_CAND
    p = "|abc|def|" |> md_blockifier
    @test p[1].name == :TABLE_ROW_CAND
    p = "|abc|def" |> md_blockifier   # doesn't end with '|'
    @test length(p) == 0
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# ===========================================================================
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# ===========================================================================
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "f0>raw" begin
    # doesn't fail because we don't inspect what goes on in `???`.
    p = """
        foo bar ??? <!-- etc __ ??? baz
        """ |> grouper
    @test p[1].blocks[1] // "foo bar"
    @test p[1].blocks[2] // "??? <!-- etc __ ???"
    @test p[1].blocks[3] // "baz"
end

@testset "f1>inline math" begin
    s = raw"abc $ghi$ mkl"
    b = s |> md_blockifier
    @test ct(b[1]) == "ghi"
    p = s |> grouper
    @test p[1] // s
    @test p[1].blocks[2] // raw"$ghi$"

    # disable math
    s = raw"foo $800"
    p = FP.md_partition(s, disable=[:MATH_INLINE])
    @test p[1].name == :TEXT
    @test ct(p[1]) == s
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "f2>block math" begin
    p = raw"""
        abc $$x = 1+1$$ end
        """ |> grouper
    @test length(p) == 3
    @test ctf(p[2]) // "x = 1+1"
    p = raw"""
        abc \[x = 1+1\] end
        """ |> grouper
    @test length(p) == 3
    @test ctf(p[2]) // "x = 1+1"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "f6+f7>newcoms,coms" begin
    s = raw"""
        \newcommand{\foo}{abc}
        \newcommand{\bar}[1]{abc#1}
        \foo \bar{def}
        """
    p = s |> grouper
    @test p[1] // s
    @test p[1].blocks[1] // raw"\newcommand"
    @test p[1].blocks[2] // raw"{\foo}"
    @test p[1].blocks[3] // raw"{abc}"
    @test p[1].blocks[4] // "\n"
    @test p[1].blocks[5] // raw"\newcommand"
    @test p[1].blocks[6] // raw"{\bar}"
    @test p[1].blocks[7] // raw"[1]"
    @test p[1].blocks[8] // raw"{abc#1}"

    p = raw"\newcommand{\foo}  [1 ] {abc}" |> grouper
    @test p[1].blocks[2] // raw"{\foo}"
    @test p[1].blocks[3].ss == "  "
    @test p[1].blocks[4] // "[1 ]"
    @test p[1].blocks[5].ss == " "
    @test p[1].blocks[6] // "{abc}"
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "f10>cu_brackets" begin
    # see f6,f7
end

# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX

@testset "f11>dbb" begin
    p = "abc {{def}} ghi" |> grouper
    @test length(p) == 1
    @test p[1].blocks[2].name == :DBB
    @test ct(p[1].blocks[2]) // "def"
end

# ////////////////////////////////////////////////////////////////////////////
# ////////////////////////////////////////////////////////////////////////////
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# ////////////////////////////////////////////////////////////////////////////
# XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX XXX
# ////////////////////////////////////////////////////////////////////////////
# ////////////////////////////////////////////////////////////////////////////

@testset "xx corner cases" begin
    s = "[`]`]"
    b = s |> md_blockifier
    @test b[1].ss == "[`]`]"

    s = "*abc<!--d*-->"
    b = s |> md_blockifier
    @test ct(b[1]) == "d*"
    g = s |> grouper
    @test g[1].ss == s
end


text:
using ConservationLawsParticles
using SpecialFunctions, RecursiveArrayTools, DifferentialEquations, Plots
using Random

# reference solution
boxsol(t, x; a, b) = (erf((x-a)/2sqrt(t)) - erf((x-b)/2sqrt(t))) / (2*(b-a))
# refsol(t, x) = sum(boxsol(t, x; a=x0[i-1], b=x0[i]) for i in eachindex(x0)[2:end]) / (length(x0)-1)
# refsol(t, x) = 1 / (pi * (1 + x^2))

# model
V(t, x) = -2x / (1 + x^2) # 2 + 0.5sin(2x) + 0.3sin(2sqrt(2)*x) + 0.2cos(2sqrt(7)*x)
W(t, x) = 0.0 # x^2/2 - abs(x)
W′(t, x) = 0.0 # x - sign(x)
mob(ρ) = 1.0 # max((ρ-1)^2 * (2ρ+1), 0)
A(ρ) = ρ^2
simple_model = DiffusiveSampledModel((V,), ((W′,),), (mob,), (SimpleDiffusion(A),))
mean_model = DiffusiveSampledModel((V,), ((W′,),), (mob,), (MeanDiffusion(A),))
min_model = DiffusiveSampledModel((V,), ((W′,),), (mob,), (MinDiffusion(A),))

# initial condition
Random.seed!(42)
n = 60
x0 = ArrayPartition(sort(vcat(randn(n) .- 1, .5randn(n) .+ 1)))

# time span
tspan = (0., 15.)

# ODE system for the particles
simple_prob = ODEProblem(velocities_diff!, x0, tspan, simple_model)
mean_prob = ODEProblem(velocities_diff!, x0, tspan, mean_model)
min_prob = ODEProblem(velocities_diff!, x0, tspan, min_model)

# solve it
abstol = reltol = 2e-7
@time simple_sol = solve(simple_prob, BS5(); abstol=abstol, reltol=reltol);
@time mean_sol = solve(mean_prob, BS5(); abstol=abstol, reltol=reltol);
@time min_sol = solve(min_prob, BS5(); abstol=abstol, reltol=reltol);

# plot the particle trajectories
plot(title="Trajectories", xlabel="time", ylabel="position", legend=false, ylims=(-5,5))
plot!(simple_sol, vars=eachindex(x0)[1:3:end], color=1)
plot!(mean_sol, vars=eachindex(x0)[1:3:end], color=2)
plot!(min_sol, vars=eachindex(x0)[1:3:end], color=3)
savefig("diffusion9.png")

# plot an animation of the density
anim = @animate for t in range(tspan...; length=100)
    plot(title="Density", xlabel="position", ylabel="density",
        xlims=(-2,2), ylims=(0,.6), legend=:topleft, size=(1200, 800))
    # plot!(x -> refsol(t, x); color=:black, label="analytic")
    densityplot!(simple_sol(t); color=1, label="simple")
    densityplot!(mean_sol(t); color=2, label="mean")
    densityplot!(min_sol(t); color=3, label="min")
end
gif(anim, "diffusion9.gif")


text:
module Minidolfin

using PyCall

function __init__()
    global FIAT = pyimport("FIAT")
    global ffc = pyimport("ffc")
    global ufl = pyimport("ufl")
end

include("mesh.jl")
include("dofmap.jl")
include("assemble.jl")
include("plot.jl")
include("io.jl")

end # module


text:
ismutable(::Type{<:StaticArrays.StaticArray}) = false
ismutable(::Type{<:StaticArrays.MArray}) = true


text:
module GraphPlot

importall BrmlFunctions
using PyPlot

export VariablePlot
type VariablePlot
    text
    x
    y
    nodetype
    nodecolor
    nodesize
    function VariablePlot(text;x=0,y=0,nodetype="circle",nodecolor="black",nodesize=0.05)
        new(text,x,y,nodetype,nodecolor,nodesize)
    end
end

export PlotGraph
function PlotGraph{I<:Number}(D::Dict{Integer,VariablePlot},LL::Dict{Integer,Array{I,}};arrowsize=0.1,arrowwidthscale=0.25,arrowcolor="black",PlotBidirectionalEdgeAsUndirected=true) # D is the Variable Dictionary and LL is the Dictionary List of directed connections (parent=>child)

    if false
    L=Dict{Integer,Array{Integer}}()
    # first switch to child=>parent list:
    for p=keys(LL)
        ch=LL[p]
        for c=ch
            if haskey(L,c)
                L[c]=union(L[c],p)
            else
                L[c]=union([],p)
            end
        end
    end
        end
    L=deepcopy(LL)

    an=linspace(0,2*pi,100)
    for i=1:length(D)
        xx=D[i].x; yy=D[i].y
        annotate(D[i].text,xy=[xx;yy],horizontalalignment="center",verticalalignment="center")
        if D[i].nodetype=="circle"
            r=D[i].nodesize
            plot(xx+r*cos(an),yy+r*sin(an),color=D[i].nodecolor)
        end
    end
    axis("equal")

    ky=collect(keys(L))
    lambda=linspace(0,1,100)
    for i=ky
        for j=L[i]
            xstart=D[i].x; xend=D[j].x
            ystart=D[i].y; yend=D[j].y

            if D[j].nodetype=="circle" && D[i].nodetype=="circle"
                rstart=D[i].nodesize
                rend=D[j].nodesize
                e=zeros(2,1)
                e[1]=xend-xstart
                e[2]=yend-ystart
                xvec=zeros(2,1)
                xvec[1]=xstart
                xvec[2]=ystart
                yvec=zeros(2,1)
                yvec[1]=xend
                yvec[2]=yend
                e=e./sqrt(e'*e)

                xxvec=xvec+rstart*e
                yyvec=yvec-rend*e
                xxstart=xxvec[1]
                yystart=xxvec[2]
                xxend=yyvec[1]
                yyend=yyvec[2]

                plot(lambda*xxstart+(1-lambda)*xxend,lambda*yystart+(1-lambda)*yyend,color=arrowcolor)

                if haskey(L,j) && any(L[j].==i)
                    bidirectionaledge=true
                else
                    bidirectionaledge=false
                end

                if !bidirectionaledge | (bidirectionaledge && !PlotBidirectionalEdgeAsUndirected)
                    # add an arrow:
                    eb=zeros(2,1); eb[1]=e[2]; eb[2]=-e[1];
                    avec=yvec-(rend+arrowsize)*e;
                    Bvec=avec+arrowwidthscale*arrowsize*eb;
                    Bx=Bvec[1]; By=Bvec[2];
                    Cvec=avec-arrowwidthscale*arrowsize*eb;
                    Cx=Cvec[1]; Cy=Cvec[2];
                    Ax=xxend; Ay=yyend;

                    plot(lambda*Ax+(1-lambda)*Bx,lambda*Ay+(1-lambda)*By,color=arrowcolor)
                    plot(lambda*Ax+(1-lambda)*Cx,lambda*Ay+(1-lambda)*Cy,color=arrowcolor)
                    plot(lambda*Bx+(1-lambda)*Cx,lambda*By+(1-lambda)*Cy,color=arrowcolor)
                end
            end
        end
    end

end

include("make_layout.jl")
include("poset.jl")
include("toposort.jl")

export PlaceVertices
function PlaceVertices(D,L;UseEnergyApproach=false,scale=4,crosspenalty=0.5,proximitypenalty=2,linelengthpenalty=0.5,attempts=200)

    if UseEnergyApproach
        ## This is just a quick attempt to places vertices such that they prefer not to be too close and also for connecting lines not to cross. Needs more work!

        Ebest=1000000000
        DDD=deepcopy(D)
        for goes=1:attempts
            V=length(D) # number of variables
            an=linspace(0,2pi,V)
            DD=deepcopy(D)
            for i=1:V
                #DD[i].x=cos(an[i]);     DD[i].y=sin(an[i])
                DD[i].x=scale*rand();     DD[i].y=scale*rand()
            end

            E=DAGenergy(DD,L,scale,crosspenalty,proximitypenalty,linelengthpenalty)
            for loop=1:10
                for v=randperm(V)
                    Dcand=deepcopy(DD)
                    Dcand[v].x=scale*rand()
                    Dcand[v].y=scale*rand()
                    Ecand=DAGenergy(DD,L,scale,crosspenalty,proximitypenalty,linelengthpenalty)
                    if Ecand<E
                        E=Ecand
                        DD[v]=deepcopy(Dcand[v])
                    end
                end
            end
            if E<Ebest
                Ebest=E
                DDD=deepcopy(DD)
            #println(Ebest)
            end

        end

        return DDD
    else # use Taylan's method:
        DD=deepcopy(D)
        V=length(D)
        A=zeros(V,V)
        for i=collect(keys(L))
            for j=L[i]
                A[i,j]=1
            end
        end
        xx,yy=make_layout(A)
        for i=1:V
            DD[i].x=scale*xx[i]
            DD[i].y=scale*yy[i]
        end
        return DD

    end

end

export DAGenergy
function DAGenergy(points,L,scale=4,crosspenalty=0.5,proximitypenalty=2,linelengthpenalty=0.5)

    E=0.0;
    if true
    for i=1:length(points)
        vi=vec([points[i].x,points[i].y])
        for j=setdiff(1:length(points),i)
            vj=vec([points[j].x,points[j].y])
            d=vi-vj;
            dist=sqrt(dot(d,d))
                if dist<0.1*scale
                    E=E+1
                elseif dist>0.5*scale
                    E=E+1
                end
        end
    end
    E=proximitypenalty*E
    end
    if true
        for a=collect(keys(L))
            for b=L[a]
                aa=vec([points[a].x points[a].y])
                bb=vec([points[b].x points[b].y])
                E=E+linelengthpenalty*sqrt(dot( bb-aa,bb-aa)) # encourage lines to not be too long

                for c=setdiff(collect(keys(L)),[a b])
                    for d=setdiff(L[c],[a b])
                        E=E+crosspenalty*linescross(points[a],points[b],points[c],points[d])
                    end
                end
            end
        end
    end
    return E
end

export linescross
function linescross(line1start,line1end,line2start,line2end)
    # returns true is line1 crosses line2
    a=vec([line1start.x line1start.y])
    b=vec([line1end.x line1end.y])
    c=vec([line2start.x line2start.y])
    d=vec([line2end.x line2end.y])

    g=a-c;e=b-a;f=d-c;

    alpha=dot(e,e); beta=-dot(f,e); theta=-dot(g,e);
    gamma=dot(e,f); delta=-dot(f,f); phi=-dot(g,f);
    mu=(alpha*phi-theta*gamma)/(delta*alpha-beta*gamma)
    lambda=(theta-mu*beta)/alpha
    v=g+lambda*e-mu*f
    E=dot(v,v)

    if 0<= lambda <=1 && 0<=mu<=1 && E<0.0001
        return true
    else
        return false
    end
end

end


text:
using ProjectEuler
using Test

@testset "ProjectEuler.jl" begin
    # Primality test
    @test ProjectEuler.isprime(37)
    @test ProjectEuler.isprime(7919)
    @test !ProjectEuler.isprime(4500)

    # Get smallest divisor
    @test ProjectEuler.get_smallest_divisor(155077777) == 7919
    @test ProjectEuler.get_smallest_divisor(100000000) == 2

    # Reverse digits of an integer
    @test ProjectEuler.reverse_int(155077777) == 777770551

    # Euler's totient
    @test ProjectEuler.totient(110) ≈ 40
    @test ProjectEuler.totient(10) ≈ 4
    @test ProjectEuler.totient(9) ≈ 6
end


text:
using SchumakerSpline
using Dates
tol = 10*eps()

StartDate = Date(2018, 7, 21)
x = Array{Date}(undef,1000)
for i in 1:1000
    x[i] = StartDate +Dates.Day(2* (i-1))
end

function f(x::Date)
    days_between = Dates.days(x - StartDate)
    return log(days_between+1) + sqrt(days_between)
end
y = f.(x)

spline = Schumaker(x,y)
for i in 1:length(x)
    abs(evaluate(spline, x[i]) - y[2]) < tol
end
# Evaluation with a Float64.
evaluate(spline, 11.5)

# Testing second derivatives
second_derivatives = evaluate.(spline, x,2)
maximum(second_derivatives) < tol

# Testing Integrals
function analytic_integral(lhs,rhs)
    lhs_in_days = Dates.days(lhs - StartDate)
    rhs_in_days = Dates.days(rhs - StartDate)
    return (rhs_in_days+1)*log(rhs_in_days+1)-rhs_in_days + (2/3)*rhs_in_days^(3/2) - ((lhs_in_days+1)*log(lhs_in_days+1) - lhs_in_days + (2/3)*lhs_in_days^(3/2))
end

lhs = StartDate
rhs = StartDate + Dates.Month(16)
numerical_integral = evaluate_integral(spline, lhs,rhs)
analytical = analytic_integral(lhs,rhs)
abs(  analytical - numerical_integral  ) < 1

## Testing with only one date provided.
x = Array{Date}(undef, 1)
x[1] = Date(2018, 7, 21)
y = Array{Float64}(undef, 1)
y[1] = 0.0
spline = Schumaker(x,y)
abs(evaluate(spline,  Date(2018, 7, 21))) < tol
abs(evaluate(spline,  Date(2019, 7, 21))) < tol
abs(evaluate(spline,  Date(2000, 7, 21))) < tol

## Testing with two dates provided.
x = Array{Date}(undef,2)
x[1] = Date(2018, 7, 21)
x[2] = Date(2018, 8, 21)
y = Array{Float64}(undef,2)
y[1] = 0.0
y[2] = 1.0
spline = Schumaker(x,y)
abs(evaluate(spline,  Date(2018, 7, 21))) < tol
abs(evaluate(spline,  Date(2018, 7, 30))) > tol
abs(evaluate(spline,  Date(2018, 8, 21)) - y[2]) < tol
abs(evaluate(spline,  Date(2019, 8, 21)) - y[2]) > tol

spline = Schumaker(x, y , extrapolation =  (Constant,Constant))
abs(evaluate(spline,  Date(2018, 8, 21)) - y[2]) < tol
abs(evaluate(spline,  Date(2019, 8, 21)) - y[2]) < tol

## Testing with three dates provided.
x = Array{Date}(undef,3)
x[1] = Date(2018, 7, 21)
x[2] = Date(2018, 8, 21)
x[3] = Date(2018, 9, 21)
y = Array{Float64}(undef,3)
y[1] = 0.0
y[2] = 1.0
y[3] = 1.3
spline = Schumaker(x,y)
abs(evaluate(spline,  x[2]) - y[2]) < tol


text:
## Miscelaneous methods for GASS
##

## init population ...
function init_pop(obs , sub , wei , ga )
    arrcfg = CSV.read(obs.Array_Configuration_File, datarow=4 , header=["X" , "Y", "Z" , "diam" , "name"] , delim= " ")
    pop = population(arrcfg , obs , sub , wei , ga)

    return(pop)
end

max_stars_repo_name: bosscha/GASS

max_issues_repo_name: bosscha/GASS

max_forks_repo_name: bosscha/GASS


text:
# ---
# title: 921. Minimum Add to Make Parentheses Valid
# id: problem921
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack, Greedy
# link: <https://leetcode.com/problems/minimum-add-to-make-parentheses-valid/description/>
# hidden: true
# ---
#
# Given a string `S` of `'('` and `')'` parentheses, we add the minimum number
# of parentheses ( `'('` or `')'`, and in any positions ) so that the resulting
# parentheses string is valid.
#
# Formally, a parentheses string is valid if and only if:
#
#   * It is the empty string, or
#   * It can be written as `AB` (`A` concatenated with `B`), where `A` and `B` are valid strings, or
#   * It can be written as `(A)`, where `A` is a valid string.
#
# Given a parentheses string, return the minimum number of parentheses we must
# add to make the resulting string valid.
#
#
#
# **Example 1:**
#
#
#
#     Input: "())"
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: "((("
#     Output: 3
#
#
# **Example 3:**
#
#
#
#     Input: "()"
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: "()))(("
#     Output: 4
#
#
#
# **Note:**
#
#   1. `S.length <= 1000`
#   2. `S` only consists of `'('` and `')'` characters.
#
#
#
#
## @lc code=start
using LeetCode

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


text:
doc"""
`multistep_rule(a_jm1,a_j,a_jp1,y_jm1,y_j,h)`

computes the multistep rule in Numerov's method:

$y_{j+1} = { \left(2-5h^2 a_j/6\right)y_j - \left(1+h^2 a_{j-1}/12\right)y_{j-1}
    \over 1+h^2 a_{j+1}/12 }$
"""
@fastmath function multistep_rule(a_jm1::Float64, a_j::Float64, a_jp1::Float64,
                                    y_jm1::Float64, y_j::Float64, h::Float64)
    ((2 - 5*h*h*a_j/6)*y_j - (1 + h*h*a_jm1/12)*y_jm1) / (1 + h*h*a_jp1/12)
end

doc"""
`multistep_integrator!(y,a,h,y1,y2)` is the same as `multistep_integrator`,
but overwrites `y`.
"""
function multistep_integrator!(y::AbstractVector{Float64}, a::AbstractVector{Float64},
                                    h::Float64, y1::Float64, y2::Float64)
    y[1] = y1
    y[2] = y2
    @inbounds for j in 2:length(a)-1
        y[j+1] = multistep_rule(a[j-1], a[j], a[j+1], y[j-1], y[j], h)
    end
end

doc"""
`multistep_integrator(a,h,y1,y2)` solves

$y''(x) + a(x) y(x) = 0$

on a uniform grid via the 4th order
[Numerov's method](https://en.wikipedia.org/wiki/Numerov%27s_method).

Assuming a uniform grid $[x_1, \dots, x_n]$ with spacing `h`

`a` - discretization of $a(x)$, $[a(x_1), \dots, a(x_n)]$

`y1,y2` - initial boundary condition, specifying $y(x_1)$ and $y(x_2)$

Returns the solution $[y(x_1), \dots, y(x_n)]$.
"""
function multistep_integrator(a::AbstractVector{Float64}, h::Float64, y1::Float64, y2::Float64)
    y = similar(a)
    multistep_integrator!(y,a,h,y1,y2)
    y
end

doc"""
`multistep_integrator_endpoint(y,a,h,y1,y2)`

is the same as `multistep_integrator`,
but returns only $y(x_n)$.

`y` is used as temporary storage.
"""
function multistep_integrator_endpoint(y::AbstractVector{Float64}, a::AbstractVector{Float64},
                                    h::Float64, y1::Float64, y2::Float64)
    multistep_integrator!(y,a,h,y1,y2)
    y[end]
end

doc"""
`multistep_integrator_node_count(y,a,h,y1,y2)`

is the same as `multistep_integrator`,
but returns the number of nodes in $y(x)$ (i.e. the number of sign changes).

`y` is used as temporary storage.
"""
function multistep_integrator_node_count(y::AbstractVector{Float64}, a::AbstractVector{Float64},
                                    h::Float64, y1::Float64, y2::Float64)
    y[1] = y1
    y[2] = y2
    p = 0
    @inbounds for j in 2:length(a)-1
        y[j+1] = multistep_rule(a[j-1], a[j], a[j+1], y[j-1], y[j], h)
        if sign(y[j+1]) != sign(y[j])
            p += 1
        end
    end
    p
end

function expandrange(q::Function, qmin::Int, qmax::Int, Emin::Float64, Emax::Float64, Elim::Float64)
    while q(Emin)::Int > qmin
        Emin -= Emax - Emin
        # @assert q(Emin) <= q(Emax) "node count is not monotonically increasing"
        @assert abs(Emin) < Elim "Emin exceeded Elim. decrease grid spacing, or increase Elim"
    end

    while q(Emax)::Int <= qmax
        Emax += Emax - Emin
        @assert abs(Emax) < Elim "Emax exceeded Elim. decrease grid spacing, or increase Elim"
    end

    Emin,Emax
end

doc"""
`eig(aE,dx,qmin,qmax;Emin=0.0,Emax=1.0,Elim=10000.)`

solves the eigenvalue equation

$y''(x) + a_E(x) y(x) = 0$

via the shooting method and Numerov's multistep integrator (see `multistep_integrator`),
on a uniform grid $[x_1, \dots, x_n]$ with spacing `dx`
and boundary conditions $y(x_1) = 0$ and $y(x_n) = 0$,

`aE` - the function $a_E(x)$, which returns a vector $[a_E(x_1), \dots, a_E(x_n)]$
for a given value of `E`.

`qmin,qmax` - solve for eigenfunctions with `qmin < q < qmax` nodes.

`Emin,Emax` - initial `E` interval to search. this will be expanded automatically
 (up to `±Elim`) until the range includes all requested eigenvalues.

returns `(u,v,q)`

where `u` are the eigenvalues, `v` are the eigenfunctions (the jth
eigenfunction given by `v[:,j]`), and `q` are the respective number of nodes in
the eigenfunctions.

e.g. *particle in a box*
```julia
x = lingrid(0,1,100)
dx = x[2]-x[1]
v = zeros(x)
evs,efs,q = eig(E -> 2*(E-v),dx,0,4)
```
"""
function eig(aE::Function, dx::Float64, qmin::Integer, qmax::Integer;
        Emin=0., Emax=1., Elim=10000.)
    y = similar(aE(Emin))::Vector{Float64}

    integrator = E -> multistep_integrator(aE(E), dx, 0., 1.)
    F = E -> multistep_integrator_endpoint(y, aE(E), dx, 0., 1.)  # matching function
    q = E -> multistep_integrator_node_count(y, aE(E), dx, 0., 1.)  # node count

    Emin,Emax = expandrange(q,qmin,qmax,Emin,Emax,Elim)

    I = findsteps(q, Emin, Emax, qmin, qmax)  # intervals
    u = Float64[fzero(F, a, b) for (a,b) in I]  # compute eigenvalues
    v = hcat([integrator(E) for E in u]...)  # compute eigenvectors
    q = Int64[q(a) for (a,b) in I]  # node counts

    u, v, q
end

doc"""
`eig!(u,v,aE,dx,qmin,qmax;Emin=0.0,Emax=1.0,Elim=10000.)`

same as `eig` except that the eigenvalues and eigenfunctions are appended to `u`
and `v`, respectively.

returns `q`, the number of nodes in the eigenfunctions.
"""
function eig!(u::Vector{Float64}, v::Vector{Vector{Float64}},
        aE::Function, dx::Float64, qmin::Integer, qmax::Integer;
        Emin=0., Emax=1., Elim=10000.)
    y = similar(aE(Emin))::Vector{Float64}

    integrator = E -> multistep_integrator(aE(E), dx, 0., 1.)
    F = E -> multistep_integrator_endpoint(y, aE(E), dx, 0., 1.)  # matching function
    q = E -> multistep_integrator_node_count(y, aE(E), dx, 0., 1.)  # node count

    Emin,Emax = expandrange(q,qmin,qmax,Emin,Emax,Elim)

    I = findsteps(q, Emin, Emax, qmin, qmax)  # intervals
    u_ = Float64[fzero(F, a, b) for (a,b) in I]  # compute eigenvalues
    q = Int64[q(a) for (a,b) in I]  # node counts

    append!(u, u_)
    append!(v, [integrator(E) for E in u_])  # compute eigenvectors

    q
end

max_stars_repo_name: johncsnyder/Element


text:
# .................................. read-io (header sector) .......................................................

function _fits_read_IO(filename::String)

    Base.Filesystem.isfile(filename) || error("Error: $filename: file not found in current directory")

    a = _validate_FITS_name(filename)

    o = IOBuffer()

    nbytes = Base.write(o,Base.read(filename))            # number of bytes
    nblock = nbytes ÷ 2880                                # number of blocks (2880 bytes/block)
    remain = nbytes % 2880                                # remainder (incomplete block)

    remain > 0 && return error("FitsError: format requires integer number of blocks of 2880 bytes")

    return o

end

function _read_header(o::IO, hduindex::Int)               # reads non-blank records from header

    ptr = _hdu_pointers(o)

    header::Array{String,1} = []
    record = " "
    Base.seek(o, ptr[hduindex])

    while record ≠ "END" * Base.repeat(" ",77)
        record = String(Base.read(o,80))
        Base.push!(header,record)
    end

    return FITS_header = _cast_header(header, hduindex)

end

# .................................. read-io (data sector) .......................................................

function _read_data(o::IO, hduindex::Int)                   # read all data using header information

    FITS_header = _read_header(o, hduindex)

      dicts = FITS_header.dict
    hdutype = Base.get(dicts,"XTENSION", "'PRIMARY '")
    hdutype = Base.strip(hdutype[2:9])

    hdutype == "PRIMARY"  && return _read_PRIMARY_data(o, hduindex)
    hdutype == "IMAGE"    && return _read_IMAGE_data(o, hduindex)
    hdutype == "TABLE"    && return _read_TABLE_data(o, hduindex)
    hdutype == "BINTABLE" && return _read_BINTABLE_data(o, hduindex)

    return error("FitsError: '$hdutype': not a 'FITS standard extension'")

end

function _read_PRIMARY_data(o::IO, hduindex::Int)             # read all data using header information

    ptr = _data_pointers(o)                                 # ptrd: start-of-data pointers

    FITS_header = _read_header(o, hduindex)

    Base.seek(o,ptr[hduindex])

    dicts = FITS_header.dict
    ndims = Base.get(dicts,"NAXIS", 0)

    if ndims > 0
         dims = Core.tuple([Base.get(dicts,"NAXIS$n", 0) for n=1:ndims[1]]...)      # e.g. dims[1]=(512,512,1)
        ndata = Base.prod(dims)                                                     # number of data points
        nbits = Base.get(dicts,"BITPIX", 0)
        bzero = Base.get(dicts,"BZERO", 0.0)                # default 0.0 by convention (also for unsigned Ints)
            E = _fits_eltype(nbits, bzero)
         data = [Base.read(o,E) for n=1:ndata]
         data = Base.ntoh.(data)                            # change from network to host ordering
         data = data .+ E(bzero)                            # offset from Int to UInt
         data = Base.reshape(data, dims)
    else
        data = []
    end

    return FITS_data = _cast_data(hduindex, "PRIMARY", data)

end

function _read_IMAGE_data(o::IO, hduindex::Int)             # read all data using header information

    ptr = _data_pointers(o)                                 # ptrd: start-of-data pointers

    FITS_header = _read_header(o, hduindex)

    Base.seek(o,ptr[hduindex])

    dicts = FITS_header.dict
    ndims = Base.get(dicts,"NAXIS", 0)

    if ndims > 0
         dims = Core.tuple([Base.get(dicts,"NAXIS$n", 0) for n=1:ndims[1]]...)      # e.g. dims[1]=(512,512,1)
        ndata = Base.prod(dims)                                                     # number of data points
        nbits = Base.get(dicts,"BITPIX", 0)
        bzero = Base.get(dicts,"BZERO", 0.0)                # default 0.0 by convention (also for unsigned Ints)
            E = _fits_eltype(nbits, bzero)
         data = [Base.read(o,E) for n=1:ndata]
         data = Base.ntoh.(data)                            # change from network to host ordering
         data = data .+ E(bzero)                            # offset from Int to UInt
         data = Base.reshape(data, dims)
    else
        data = []
    end

    return FITS_data = _cast_data(hduindex, "IMAGE", data)

end

function _read_TABLE_data(o::IO, hduindex::Int)

    ptr = _data_pointers(o)

    FITS_header = _read_header(o, hduindex)

    dicts = FITS_header.dict
    lrecs = Base.get(dicts,"NAXIS1", 0)
    nrecs = Base.get(dicts,"NAXIS2", 0)

    Base.seek(o,ptr[hduindex])

    data = [String(Base.read(o,lrecs)) for i=1:nrecs]

    return FITS_data = _cast_data(hduindex, "TABLE", data)

end


text:
using CollegeStratCollege
using Test, Random
using ModelObjectsLH

@testset "All" begin
    include("dropout_rule_test.jl");
    include("grad_rule_test.jl");
    include("logistic_test.jl");
    include("course_grid_test.jl");
    include("study_time_grid_test.jl");
    include("hproduction_test.jl");
    include("h_shocks_test.jl");
    include("learning_test.jl");
    include("tuition_test.jl");
    include("college_test.jl");
end

# ------------


text:
@testset "Basic Structural With Explanatory Model" begin
    @test has_fit_methods(BasicStructuralExplanatory)
    y = CSV.File(StateSpaceModels.AIR_PASSENGERS) |> DataFrame
    logap = log.(y.passengers)
    X = rand(length(logap), 2)
    model = BasicStructuralExplanatory(logap, 12, X)
    fit!(model)
    model.results
    # forecasting
    # For a fixed forecasting explanatory the variance must not decrease
    forec = forecast(model, ones(10, 2))
    @test monotone_forecast_variance(forec)
    kf = kalman_filter(model)
    ks = kalman_smoother(model)
    @test_throws AssertionError simulate_scenarios(model, 10, 1000, ones(5, 2))
    scenarios = simulate_scenarios(model, 10, 1000, ones(10, 2))
    test_scenarios_adequacy_with_forecast(forec, scenarios)
    scenarios = simulate_scenarios(model, 10, 500, ones(10, 2, 500))
    test_scenarios_adequacy_with_forecast(forec, scenarios)
end


text:
# params(f) Based on deepcopy_internal:

params(f) = (ps=Param[]; params_internal(f,ps,IdDict()); ps)

# Tapes can only have params at the top level, so here is a more efficient implementation for tapes:
params(t::Tape) = (n.Value for n in t.list if n.Value isa Param)

params_internal(p::Param, ps::Vector{Param}, d::IdDict) = if !haskey(d,p); d[p]=true; push!(ps,p); end

params_internal(x::Union{Symbol,Core.MethodInstance,Method,GlobalRef,DataType,Union,UnionAll,Task,Regex},
                ps::Vector{Param}, stackdict::IdDict) = return
params_internal(x::Tuple, ps::Vector{Param}, stackdict::IdDict) =
    for p in x; params_internal(p, ps, stackdict); end

params_internal(x::Module, ps::Vector{Param}, stackdict::IdDict) = return

params_internal(x::String, ps::Vector{Param}, stackdict::IdDict) = return

function params_internal(x::Core.SimpleVector, ps::Vector{Param}, stackdict::IdDict)
    if haskey(stackdict, x)
        return
    end
    stackdict[x] = true
    for p in x; params_internal(p, ps, stackdict); end
end

function params_internal(@nospecialize(x), ps::Vector{Param}, stackdict::IdDict)
    T = typeof(x)::DataType
    nf = nfields(x)
    (isbitstype(T) || nf == 0) && return
    if haskey(stackdict, x)
        return
    end
    if T.mutable
        stackdict[x] = true
    end
    for i in 1:nf
        if isdefined(x,i)
            params_internal(getfield(x,i), ps, stackdict)
        end
    end
end

function params_internal(x::Array, ps::Vector{Param}, stackdict::IdDict)
    if haskey(stackdict, x)
        return
    end
    _params_array_t(x, eltype(x), ps, stackdict)
end

function _params_array_t(@nospecialize(x), T, ps::Vector{Param}, stackdict::IdDict)
    stackdict[x] = true
    if isbitstype(T)
        return
    end
    for i = 1:(length(x)::Int)
        if ccall(:jl_array_isassigned, Cint, (Any, Csize_t), x, i-1) != 0
            xi = ccall(:jl_arrayref, Any, (Any, Csize_t), x, i-1)
            if !isbits(xi)
                params_internal(xi, ps, stackdict)
            end
        end
    end
end

function params_internal(x::Union{Dict,IdDict}, ps::Vector{Param}, stackdict::IdDict)
    if haskey(stackdict, x)
        return
    end
    stackdict[x] = true
    for (k, v) in x
        params_internal(k, ps, stackdict)
        params_internal(v, ps, stackdict)
    end
end


text:
immutable Timeseries{T} <: AbstractArray{T,1}
    n::Int
    t::Vector{T}

    Δt::T
    step::Int

    function Timeseries(n, Δt, step)
        @assert T <: Real
        @assert n > 0
        @assert step > 0
        @assert n ≥ step

        t = zeros(T, n+1)
        new(n, t, Δt, step)
    end
end

function Timeseries{T}(n::Int, Δt::T, step::Int=1)
    return Timeseries{T}(n, Δt, step)
end

Base.eltype{T}(t::Timeseries{T}) = T
Base.ndims(t::Timeseries) = 1
Base.size(t::Timeseries) = size(t.t)
Base.length(t::Timeseries) = length(t.t)
Base.endof(t::Timeseries) = length(t)
Base.indices(t::Timeseries) = (0:t.n,)
Base.indices(t::Timeseries, d) = indices(t)[d]
Base.strides(t::Timeseries) = (1,)
Base.stride(t::Timeseries, d) = strides(t)[d]

@inline function Base.setindex!(ts::Timeseries, t, i::Int)
    @boundscheck checkbounds(ts.t, i+1)
    @inbounds setindex!(ts.t, t, i+1)
end

@inline function Base.getindex(ts::Timeseries, i::Int)
    @boundscheck checkbounds(ts.t, i+1)
    @inbounds t = getindex(ts.t, i+1)
    return t
end

function compute_timeseries!(ts::Timeseries)
    for n in 1:ts.n
        ts[n] = ts[0] + n * ts.step * ts.Δt
    end
end


text:
module WTExperiments

import AbstractActuators as Move

end