text:
if !(VERSION.major == 0 && VERSION.minor <= 6)
    using Test
else
    using Base.Test
end
using FinancialToolbox

print_colored("Starting Hyper Dual Numbers Test\n", :green)

#Test Parameters
testToll = 1e-14;
spot = 10;
K = 10;
r = 0.02;
T = 2.0;
sigma = 0.2;
d = 0.01;
assert_(value, toll) = @test abs(value) < toll
#EuropeanCall Option
PriceCall = blsprice(spot, K, r, T, sigma, d);
DeltaCall = blsdelta(spot, K, r, T, sigma, d);
ThetaCall = blstheta(spot, K, r, T, sigma, d);
RhoCall = blsrho(spot, K, r, T, sigma, d);
LambdaCall = blslambda(spot, K, r, T, sigma, d);
VannaCall = blsvanna(spot, K, r, T, sigma, d);

#EuropeanPut Option
PricePut = blsprice(spot, K, r, T, sigma, d, false);
DeltaPut = blsdelta(spot, K, r, T, sigma, d, false);
ThetaPut = blstheta(spot, K, r, T, sigma, d, false);
RhoPut = blsrho(spot, K, r, T, sigma, d, false);
LambdaPut = blslambda(spot, K, r, T, sigma, d, false);
VannaPut = blsvanna(spot, K, r, T, sigma, d, false);

#Equals for both Options
Gamma = blsgamma(spot, K, r, T, sigma, d);
Vega = blsvega(spot, K, r, T, sigma, d);

########HYPER DUAL NUMBERS
using HyperDualNumbers;
DerToll = 1e-13;
di = 1e-15;
#Function definition
SpotHyper = hyper(spot, 1.0, 1.0, 0.0);
rHyper = hyper(r, 1.0, 1.0, 0.0);
THyper = hyper(T, 1.0, 1.0, 0.0);
SigmaHyper = hyper(sigma, 1.0, 1.0, 0.0);
KHyper = hyper(K, 1.0, 1.0, 0.0);

#Function definition
#Call
F(spot) = blsprice(spot, K, r, T, sigma, d);
FF(spot) = blsdelta(spot, K, r, T, sigma, d);
G(r) = blsprice(spot, K, r, T, sigma, d);
H(T) = blsprice(spot, K, r, T, sigma, d);
L(sigma) = blsprice(spot, K, r, T, sigma, d);
P(spot) = blsprice(spot, K, r, T, sigma, d) * spot.value / blsprice(spot.value, K, r, T, sigma, d);
V(sigma) = blsdelta(spot, K, r, T, sigma, d);
VV(spot, sigma) = blsprice(spot, K, r, T, sigma, d);
#Put
F1(spot) = blsprice(spot, K, r, T, sigma, d, false);
FF1(spot) = blsdelta(spot, K, r, T, sigma, d, false);
G1(r) = blsprice(spot, K, r, T, sigma, d, false);
H1(T) = blsprice(spot, K, r, T, sigma, d, false);
P1(spot) = blsprice(spot, K, r, T, sigma, d, false) * spot.value / blsprice(spot.value, K, r, T, sigma, d, false);
V1(sigma) = blsdelta(spot, K, r, T, sigma, d, false);
V11(spot, sigma) = blsprice(spot, K, r, T, sigma, d, false);

#TEST
print_colored("--- European Call Sensitivities: HyperDualNumbers\n", :yellow)
print_colored("-----Testing Delta\n", :blue);
assert_(F(SpotHyper).epsilon1 - DeltaCall, DerToll)
print_colored("-----Testing Gamma\n", :blue);
assert_(FF1(SpotHyper).epsilon1 - Gamma, DerToll)
print_colored("-----Testing Gamma  2\n", :blue);
assert_(F(SpotHyper).epsilon12 - Gamma, DerToll)
print_colored("-----Testing Rho\n", :blue);
assert_(G(rHyper).epsilon1 - RhoCall, DerToll)
print_colored("-----Testing Theta\n", :blue);
assert_(-H(THyper).epsilon1 - ThetaCall, DerToll)
print_colored("-----Testing Lambda\n", :blue);
assert_(P(SpotHyper).epsilon1 - LambdaCall, DerToll)
print_colored("-----Testing Vanna\n", :blue);
assert_(V(SigmaHyper).epsilon1 - VannaCall, DerToll)
print_colored("-----Testing Vanna  2\n", :blue);
assert_(VV(hyper(spot, 1, 0, 0), hyper(sigma, 0, 1, 0)).epsilon12 - VannaCall, DerToll)

## Complex Test with Complex Step Approximation for European Put
#TEST
print_colored("--- European Put Sensitivities: HyperDualNumbers\n", :yellow)
print_colored("-----Testing Delta\n", :blue);
assert_(F1(SpotHyper).epsilon1 - DeltaPut, DerToll)
print_colored("-----Testing Rho\n", :blue);
assert_(G1(rHyper).epsilon1 - RhoPut, DerToll)
print_colored("-----Testing Theta\n", :blue);
assert_(-H1(THyper).epsilon1 - ThetaPut, DerToll)
print_colored("-----Testing Lambda\n", :blue);
assert_(P1(SpotHyper).epsilon1 - LambdaPut, DerToll)
print_colored("-----Testing Vanna\n", :blue);
assert_(V1(SigmaHyper).epsilon1 - VannaPut, DerToll)
print_colored("-----Testing Vanna  2\n", :blue);
assert_(V11(hyper(spot, 1, 0, 0), hyper(sigma, 0, 1, 0)).epsilon12 - VannaPut, DerToll)
print_colored("-----Testing Vega\n", :blue);
assert_(L(SigmaHyper).epsilon1 - Vega, DerToll)
print_colored("Hyper Dual Numbers Test Passed\n\n", :green)

#TEST OF INPUT VALIDATION
print_colored("Starting Input Validation Test Hyper Dual Numbers\n", :magenta)
print_colored("----Testing Negative  Spot Price \n", :cyan)
@test_throws(ErrorException, blsprice(-SpotHyper, K, r, T, sigma, d))
@test_throws(ErrorException, blkprice(-SpotHyper, K, r, T, sigma))
@test_throws(ErrorException, blsdelta(-SpotHyper, K, r, T, sigma, d))
@test_throws(ErrorException, blsgamma(-SpotHyper, K, r, T, sigma, d))
@test_throws(ErrorException, blstheta(-SpotHyper, K, r, T, sigma, d))
@test_throws(ErrorException, blsrho(-SpotHyper, K, r, T, sigma, d))
@test_throws(ErrorException, blspsi(-SpotHyper, K, r, T, sigma, d))
@test_throws(ErrorException, blslambda(-SpotHyper, K, r, T, sigma, d))
@test_throws(ErrorException, blsvanna(-SpotHyper, K, r, T, sigma, d))
@test_throws(ErrorException, blsvega(-SpotHyper, K, r, T, sigma, d))

print_colored("----Testing Negative  Strike Price \n", :cyan)
@test_throws(ErrorException, blsprice(spot, -KHyper, r, T, sigma, d))
@test_throws(ErrorException, blkprice(spot, -KHyper, r, T, sigma))
@test_throws(ErrorException, blsdelta(spot, -KHyper, r, T, sigma, d))
@test_throws(ErrorException, blsgamma(spot, -KHyper, r, T, sigma, d))
@test_throws(ErrorException, blstheta(spot, -KHyper, r, T, sigma, d))
@test_throws(ErrorException, blsrho(spot, -KHyper, r, T, sigma, d))
@test_throws(ErrorException, blspsi(spot, -KHyper, r, T, sigma, d))
@test_throws(ErrorException, blslambda(spot, -KHyper, r, T, sigma, d))
@test_throws(ErrorException, blsvanna(spot, -KHyper, r, T, sigma, d))
@test_throws(ErrorException, blsvega(spot, -KHyper, r, T, sigma, d))

print_colored("----Testing Negative  Time to Maturity \n", :cyan)
@test_throws(ErrorException, blsprice(spot, K, r, -THyper, sigma, d))
@test_throws(ErrorException, blkprice(spot, K, r, -THyper, sigma))
@test_throws(ErrorException, blsdelta(spot, K, r, -THyper, sigma, d))
@test_throws(ErrorException, blsgamma(spot, K, r, -THyper, sigma, d))
@test_throws(ErrorException, blstheta(spot, K, r, -THyper, sigma, d))
@test_throws(ErrorException, blsrho(spot, K, r, -THyper, sigma, d))
@test_throws(ErrorException, blspsi(spot, K, r, -THyper, sigma, d))
@test_throws(ErrorException, blslambda(spot, K, r, -THyper, sigma, d))
@test_throws(ErrorException, blsvanna(spot, K, r, -THyper, sigma, d))
@test_throws(ErrorException, blsvega(spot, K, r, -THyper, sigma, d))

print_colored("----Testing Negative  Volatility \n", :cyan)
@test_throws(ErrorException, blsprice(spot, K, r, T, -SigmaHyper, d))
@test_throws(ErrorException, blkprice(spot, K, r, T, -SigmaHyper))
@test_throws(ErrorException, blsdelta(spot, K, r, T, -SigmaHyper, d))
@test_throws(ErrorException, blsgamma(spot, K, r, T, -SigmaHyper, d))
@test_throws(ErrorException, blstheta(spot, K, r, T, -SigmaHyper, d))
@test_throws(ErrorException, blsrho(spot, K, r, T, -SigmaHyper, d))
@test_throws(ErrorException, blspsi(spot, K, r, T, -SigmaHyper, d))
@test_throws(ErrorException, blslambda(spot, K, r, T, -SigmaHyper, d))
@test_throws(ErrorException, blsvanna(spot, K, r, T, -SigmaHyper, d))
@test_throws(ErrorException, blsvega(spot, K, r, T, -SigmaHyper, d))

print_colored("Hyper Dual Input Validation Test Passed\n", :magenta)


text:
function test_transition(pb1::Union{MDP, POMDP}, pb2::Union{SparseTabularMDP, SparseTabularPOMDP})
    @testset "transition" begin
        for s in states(pb1)
            si = stateindex(pb1, s)
            for a in actions(pb1)
                ai = actionindex(pb1, a)
                td1 = transition(pb1, s, a)
                td2 = transition(pb2, si, ai)
                for (sp, p) in weighted_iterator(td1)
                    spi = stateindex(pb1, sp)
                    @test pdf(td2, spi) == p
                end
            end
        end
    end
end

function test_reward(pb1::Union{MDP, POMDP}, pb2::Union{SparseTabularMDP, SparseTabularPOMDP})
    @testset "reward" begin
        for s in states(pb1)
            si = stateindex(pb1, s)
            for a in actions(pb1)
                ai = actionindex(pb1, a)
                if isterminal(pb1, s)
                    @test reward(pb2, si, ai) == 0.0
                else
                    td = transition(pb1, s, a)
                    r_pb1 = 0.0
                    for (sp, p) in weighted_iterator(td)
                        if p > 0.0
                            r_pb1 += p*reward(pb1, s, a, sp)
                        end
                    end
                    @test r_pb1 ≈ reward(pb2, si, ai)
                end
            end
        end
    end
end

function test_observation(pb1::POMDP, pb2::SparseTabularPOMDP)
    @testset "transitions" begin
        for s in states(pb1)
            for a in actions(pb1)
                od1 = observation(pb1, a, s)
                si = stateindex(pb1, s)
                ai = actionindex(pb1, a)
                od2 = observation(pb2, ai, si)
                for (o, p) in weighted_iterator(od1)
                    oi = obsindex(pb1, o)
                    @test pdf(od2, oi) == p
                end
            end
        end
    end
end

## MDP

@testset "Random MDP" begin

mdp = RandomMDP(100, 4, 0.95)

@show_requirements SparseTabularMDP(mdp)

smdp = SparseTabularMDP(mdp)

smdp2 = SparseTabularMDP(smdp)

@test smdp2.T == smdp.T
@test smdp2.R == smdp.R
@test smdp2.discount == smdp.discount

smdp3 = SparseTabularMDP(smdp, reward = zeros(length(states(mdp)), length(actions(mdp))))
@test smdp3.T == smdp.T
@test smdp3.R != smdp.R
smdp3 = SparseTabularMDP(smdp, transition = [sparse(1:length(states(mdp)), 1:length(states(mdp)), 1.0) for a in 1:length(actions(mdp))])
@test smdp3.T != smdp.T
@test smdp3.R == smdp.R
smdp3 = SparseTabularMDP(smdp, discount = 0.5)
@test smdp3.discount != smdp.discount

@test size(transition_matrix(smdp, 1)) == (length(states(smdp)), length(states(smdp)))
@test length(reward_vector(smdp, 1)) == length(states(smdp))

gw = SimpleGridWorld()
sparsegw = SparseTabularMDP(gw)
@test isterminal(sparsegw, 101)
@inferred actions(sparsegw, 101)
@test actions(sparsegw, 101) == collect(actions(sparsegw))

end

## POMDP

function test_pomdp(pomdp::POMDP)
    @show_requirements SparseTabularPOMDP(pomdp)

    rand_mdp = RandomMDP(100, 4, 0.95)

    spomdp = SparseTabularPOMDP(pomdp)

    @testset "generic" begin
        spomdp2 = SparseTabularPOMDP(spomdp)

        @test spomdp2.T == spomdp.T
        @test spomdp2.R == spomdp.R
        @test spomdp2.O == spomdp.O
        @test spomdp2.discount == spomdp.discount

        spomdp3 = SparseTabularPOMDP(spomdp, reward = zeros(length(states(rand_mdp)), length(actions(rand_mdp))))
        @test spomdp3.T == spomdp.T
        @test spomdp3.R != spomdp.R
        spomdp3 = SparseTabularPOMDP(spomdp, transition = [sparse(1:length(states(rand_mdp)), 1:length(states(rand_mdp)), 1.0) for a in 1:length(actions(rand_mdp))])
        @test spomdp3.T != spomdp.T
        @test spomdp3.R == spomdp.R
        spomdp3 = SparseTabularPOMDP(spomdp, discount = 0.5)
        @test spomdp3.discount != spomdp.discount
        @test size(observation_matrix(spomdp, 1)) == (length(states(spomdp)), length(observations(spomdp)))
        @show size(observation_matrix(spomdp, 1)) == (length(states(spomdp)), length(observations(spomdp)))
        @test observation_matrices(spomdp) == spomdp2.O
        @test transition_matrices(spomdp) == spomdp2.T
        @test reward_matrix(spomdp) == spomdp2.R

        @test length(states(pomdp)) == length(states(spomdp))
        @test length(actions(pomdp)) == length(actions(spomdp))
        @test length(observations(pomdp)) == length(observations(spomdp))
        @test statetype(spomdp) == Int64
        @test actiontype(spomdp) == Int64
        @test obstype(spomdp) == Int64
        @test discount(spomdp) == discount(pomdp)
        @test length(spomdp.terminal_states) == length([s for s in states(pomdp) if isterminal(pomdp, s)])
    end

    test_transition(pomdp, spomdp)
    test_reward(pomdp, spomdp)
    test_observation(pomdp, spomdp)
end

@testset "Tiger POMDP" begin
    test_pomdp(TigerPOMDP())
end
@testset "TMaze POMDP" begin
    test_pomdp(TMaze())
end

ext: jl


text:
using InterpretableModels
using Test

@testset "InterpretableModels.jl" begin
    include("test_utils.jl")
    include("test_slim_scoring.jl")
    include("test_lagged_factor_linear_itrm.jl")
    include("test_moving_average_linear_itrm.jl")
end


text:
using Documenter, DemoCards

theme = cardtheme()
simplest_demopage = makedemos("demos/simplest_demopage")
gallery_of_packages = makedemos("demos/gallery_of_packages")

format = Documenter.HTML(edit_branch = "master",
                         assets = [theme])

makedocs(format = format,
         pages = [
            "Home" => "index.md",
            "QuickStart" => "quickstart.md",
            "Concepts" => "concepts.md",
            "Examples" => [
                "Simplest Demopage" => simplest_demopage,
                "Gallery of Packages" => gallery_of_packages,
            ],
            "Package References" => "references.md"
         ],
         sitename = "DemoCards.jl")

deploydocs(repo = "github.com/johnnychen94/DemoCards.jl.git")


text:
using MethodOfLines, Test, ModelingToolkit, SymbolicUtils

@testset "count differentials 1D" begin
    @parameters t x
    @variables u(..)
    Dt = Differential(t)

    Dx = Differential(x)
    eq  = Dt(u(t,x)) ~ -Dx(u(t,x))
    @test first(MethodOfLines.differential_order(eq.rhs, x.val)) == 1
    @test isempty(MethodOfLines.differential_order(eq.rhs, t.val))
    @test first(MethodOfLines.differential_order(eq.lhs, t.val)) == 1
    @test isempty(MethodOfLines.differential_order(eq.lhs, x.val))

    Dxx = Differential(x)^2
    eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
    @test first(MethodOfLines.differential_order(eq.rhs, x.val)) == 2
    @test isempty(MethodOfLines.differential_order(eq.rhs, t.val))
    @test first(MethodOfLines.differential_order(eq.lhs, t.val)) == 1
    @test isempty(MethodOfLines.differential_order(eq.lhs, x.val))

    Dxxxx = Differential(x)^4
    eq  = Dt(u(t,x)) ~ -Dxxxx(u(t,x))
    @test first(MethodOfLines.differential_order(eq.rhs, x.val)) == 4
    @test isempty(MethodOfLines.differential_order(eq.rhs, t.val))
    @test first(MethodOfLines.differential_order(eq.lhs, t.val)) == 1
    @test isempty(MethodOfLines.differential_order(eq.lhs, x.val))
end

@testset "count differentials 2D" begin
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dt = Differential(t)

    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y))
    @test first(MethodOfLines.differential_order(eq.rhs, x.val)) == 2
    @test first(MethodOfLines.differential_order(eq.rhs, y.val)) == 2
    @test isempty(MethodOfLines.differential_order(eq.rhs, t.val))
    @test first(MethodOfLines.differential_order(eq.lhs, t.val)) == 1
    @test isempty(MethodOfLines.differential_order(eq.lhs, x.val))
    @test isempty(MethodOfLines.differential_order(eq.lhs, y.val))
end

@testset "count with mixed terms" begin
    @parameters t x y
    @variables u(..)
    Dxx = Differential(x)^2
    Dyy = Differential(y)^2
    Dx = Differential(x)
    Dy = Differential(y)
    Dt = Differential(t)

    eq  = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + Dx(Dy(u(t,x,y)))
    @test MethodOfLines.differential_order(eq.rhs, x.val) == Set([2, 1])
    @test MethodOfLines.differential_order(eq.rhs, y.val) == Set([2, 1])
end

@testset "Kuramoto–Sivashinsky equation" begin
    @parameters x, t
    @variables u(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dx2 = Differential(x)^2
    Dx3 = Differential(x)^3
    Dx4 = Differential(x)^4

    α = 1
    β = 4
    γ = 1
    eq = Dt(u(x,t)) + u(x,t)*Dx(u(x,t)) + α*Dx2(u(x,t)) + β*Dx3(u(x,t)) + γ*Dx4(u(x,t)) ~ 0
    @test MethodOfLines.differential_order(eq.lhs, x.val) == Set([4, 3, 2, 1])
end

@test_broken begin #@testset "Flatten division" begin
    @parameters x y z t

    #@test_broken isequal(operation(((y^(-1.0+eps(Float64)))*x~0).lhs), *)

    @test_broken isequal(operation(MethodOfLines.flatten_division((x/z~0).lhs)), *)
    @test_broken isequal(operation(MethodOfLines.flatten_division((x*y/z~0).lhs)), *)
    @test_broken isequal(operation(MethodOfLines.flatten_division((x/(y*z)~0).lhs)), *)
    @test_broken isequal(operation(MethodOfLines.flatten_division((x*y/(z*t)~0).lhs)), *)

end

@testset "Split terms" begin
    @parameters t x y r
    @variables u(..) v(..)
    Dt = Differential(t)
    Dx = Differential(x)
    Dr = Differential(r)
    Dxx = Differential(x)^2

    t_min= 0.
    t_max = 2.0
    x_min = 0.
    x_max = 2.
    y_min = 0.
    y_max = 2.

    analytic_sol_func(t,x,y) = exp(x+y)*cos(x+y+4t)

    bcs = [u(0,x) ~ - x * (x-1) * sin(x),
       v(0,x) ~ - x * (x-1) * sin(x),
       u(t,0) ~ 0.0, u(t,1) ~ 0.0,
       v(t,0) ~ 0.0, v(t,1) ~ 0.0,
       u(t_min,x,y) ~ analytic_sol_func(t_min,x,y),
       u(t,x_min,y) ~ analytic_sol_func(t,x_min,y),
       u(t,x_max,y) ~ analytic_sol_func(t,x_max,y),
       u(t,x,y_min) ~ analytic_sol_func(t,x,y_min),
       u(t,x,y_max) ~ analytic_sol_func(t,x,y_max),
       u(0,r) ~ - r * (r-1) * sin(r),
       Dr(u(t,0)) ~ 0.0, u(t,1) ~ sin(1)]

    for bc in bcs
        terms = MethodOfLines.split_terms(bc)
        @test terms isa Vector

    end
end

@testset "Periodic Wraparound" begin
    I = CartesianIndex(2, 5)
    @test MethodOfLines._wrapperiodic(I, 2, 2, 4) == CartesianIndex(2, 2)

    I = CartesianIndex(1, 4)
    @test MethodOfLines._wrapperiodic(I, 2, 1, 4) == CartesianIndex(4, 4)

    I = CartesianIndex(-1, 2)
    @test MethodOfLines._wrapperiodic(I, 2, 1, 4) == CartesianIndex(2, 2)
end


text:
using Base.Test
import Multirate
import DSP
import Multirate.NaiveResamplers

# Disable time and printf macros when not running interactivly ( for travis )
if !isinteractive()
    macro time( ex )
        quote
            $(esc(ex))
        end
    end

    macro printf( args... )
    end
end

#==============================================================================#
#               ____ _ _  _ ____ _    ____    ____ ____ ___ ____               #
#               [__  | |\ | | __ |    |___    |__/ |__|  |  |___               #
#               ___] | | \| |__] |___ |___    |  \ |  |  |  |___               #
#==============================================================================#

function test_singlerate( h, x; atol = 0 )
    xLen       = length( x )
    hLen       = length( h )
    pivotPoint = min( rand(50:150, 1)[1], floor( Int, xLen/4 ))
    x1         = x[ 1 : pivotPoint ]
    x2         = x[ pivotPoint+1 : end ]

    @printf( "\n\n" )
    @printf( "____ _ _  _ ____ _    ____    ____ ____ ___ ____\n" )
    @printf( "[__  | |\\ | | __ |    |___    |__/ |__|  |  |___\n" )
    @printf( "___] | | \\| |__] |___ |___    |  \\ |  |  |  |___\n" )
    @printf( "\nTesting single-rate fitering, h is %s, x is %s. xLen = %d, hLen = %d", string(eltype(h)), string(eltype(x)), xLen, hLen )

    @printf( "\n\tBase.filt\n\t\t")
    @time baseResult = Base.filt( h, 1.0, x )

    if method_exists( DSP.filt, ( typeof(h), typeof(x) ))
        @printf( "\n\tDSP.filt\n\t\t")
        @time dspResult = DSP.filt( h, x )
    end

    @printf( "\n\tMultirate.filt( h, x, 1//1 )\n\t\t" )
    @time statelessResult = Multirate.filt( h, x )

    @printf( "\n\tMultirate.filt. length( x1 ) = %d, length( x2 ) = %d\n\t\t", length( x1 ), length( x2 ) )
    self = Multirate.FIRFilter( h, 1//1 )
    @time begin
        y1 = Multirate.filt( self, x1 )
        y2 = Multirate.filt( self, x2 )
    end
    statefulResult = [ y1; y2 ]

    @printf( "\n\tMultirate.filt filt. Piecewise for first %d inputs\n\t\t", length( x1 ) )
    Multirate.reset( self )
    @time begin
        for i in 1:length(x1)
            y1[i] = Multirate.filt( self, x1[i:i] )[1]
        end
        y2 = Multirate.filt( self, x2 )
    end
    piecewiseResult = [ y1; y2 ]

     if isapprox( baseResult, statelessResult, atol=atol ) && isapprox( baseResult, statefulResult, atol=atol ) && isapprox( baseResult, piecewiseResult, atol=atol )
         return true
     end

     display( [ baseResult statelessResult statefulResult piecewiseResult ] )

     return false

end

#==============================================================================#
#                      ___  ____ ____ _ _  _ ____ ___ ____                     #
#                      |  \ |___ |    | |\/| |__|  |  |___                     #
#                      |__/ |___ |___ | |  | |  |  |  |___                     #
#==============================================================================#

function test_decimation( h, x, decimation; atol=0 )
    xLen       = length( x )
    hLen       = length( h )
    pivotPoint = min( rand(50:150, 1)[1], floor( Int, xLen/4 ))
    x1         = x[ 1 : pivotPoint ]
    x2         = x[ pivotPoint+1 : end ]

    @printf( "\n\n" )
    @printf( "___  ____ ____ _ _  _ ____ ___ _ ____ _  _ \n" )
    @printf( "|  \\ |___ |    | |\\/| |__|  |  | |  | |\\ | \n" )
    @printf( "|__/ |___ |___ | |  | |  |  |  | |__| | \\| \n" )
    @printf( "\nTesting decimation. h::%s, x::%s. xLen = %d, hLen = %d, decimation = %d", string(typeof(h)), string(typeof(h)), xLen, hLen, decimation )

    @printf( "\n\tNaive decimation with Base.filt\n\t\t")
    @time begin
        baseResult   = Base.filt( h, one(eltype(h)), x )
        baseResult   = baseResult[1:decimation:end]
    end

    if method_exists( DSP.filt, ( typeof(h), typeof(x) ))
        @printf( "\n\tNaive decimation with DSP.filt\n\t\t")
        @time begin
            dspResult = DSP.filt( h, x )
            dspResult = dspResult[1:decimation:end]
        end
    end

    @printf( "\n\tMultirate.filt( h, x, 1//%d)\n\t\t", decimation )
    @time statelessResult = Multirate.filt( h, x, 1//decimation )

    @printf( "\n\tMultirate.filt decimation. length( x1 ) = %d, length( x2 ) = %d\n\t\t", length( x1 ), length( x2 ) )
    self = Multirate.FIRFilter( h, 1//decimation )
    @time begin
        y1 = Multirate.filt( self, x1 )
        y2 = Multirate.filt( self, x2 )
    end
    statefulResult = [ y1; y2 ]

    @printf( "\n\tMultirate.filt decimation. Piecewise for first %d inputs.\n\t\t", length( x1 ) )
    Multirate.reset( self )
    y1 = similar( x, 0 )
    @time begin
        for i in 1:length(x1)
            append!( y1, Multirate.filt( self, x1[i:i] ) )
        end
        y2 = Multirate.filt( self, x2 )
    end
    piecewiseResult = [ y1; y2 ]

    if isapprox( baseResult, statelessResult, atol=atol ) && isapprox( baseResult, statefulResult, atol=atol ) && isapprox( baseResult, piecewiseResult, atol=atol )
        return true
    end

    display( [ baseResult statelessResult statefulResult piecewiseResult (baseResult.-statelessResult) (baseResult.-statefulResult) (baseResult.-piecewiseResult)] )
    println()

    return false

end

#==============================================================================#
#               _ _  _ ___ ____ ____ ___  _    ____ ____ ___ ____              #
#               | |\ |  |  |___ |__/ |__] |    |  | |__|  |  |___              #
#               | | \|  |  |___ |  \ |    |___ |__| |  |  |  |___              #
#==============================================================================#

function test_interpolation( h, x, interpolation; atol::Real = 0 )
    xLen       = length( x )
    hLen       = length( h )
    pivotPoint = min( rand(50:150, 1)[1], floor( Int, xLen/4 ))
    x1         = x[ 1 : pivotPoint ]
    x2         = x[ pivotPoint+1 : end ]

    @printf( "\n\n" )
    @printf( "_ _  _ ___ ____ ____ ___  _    ____ ____ ___ _ ____ _  _ \n" )
    @printf( "| |\\ |  |  |___ |__/ |__] |    |  | |__|  |  | |  | |\\ | \n" )
    @printf( "| | \\|  |  |___ |  \\ |    |___ |__| |  |  |  | |__| | \\| \n" )
    @printf( "\nTesting interpolation, h::%s, x::%s. xLen = %d, hLen = %d, interpolation = %d", typeof(h), typeof(x), xLen, hLen, interpolation )

    @printf( "\n\tNaive interpolation with Base.filt\n\t\t")
    @time begin
        xZeroStuffed = zeros( eltype(x), xLen * interpolation )
        for n = 0:xLen-1;
            xZeroStuffed[ n*interpolation+1 ] = x[ n+1 ]
        end
        baseResult = Base.filt( h, one(eltype(h)), xZeroStuffed )
    end

    if method_exists( DSP.filt, ( typeof(h), typeof(x) ))
        @printf( "\n\tNaive interpolation with DSP.filt\n\t\t")
        @time begin
            xZeroStuffed = zeros( eltype(x), xLen * interpolation )
            for n = 0:xLen-1;
                xZeroStuffed[ n*interpolation+1 ] = x[ n+1 ]
            end
            dspResult = DSP.filt( h, xZeroStuffed )
        end
    end

    @printf( "\n\tMultirate.filt( h, x, %d//1 )\n\t\t", interpolation )
    @time statelessResult = Multirate.filt( h, x, interpolation//1 )

    @printf( "\n\tMultirate.filt interpolation. length( x1 ) = %d, length( x2 ) = %d\n\t\t", length( x1 ), length( x2 ) )
    self = Multirate.FIRFilter( h, interpolation//1 )
    @time begin
        y1 = Multirate.filt( self, x1 )
        y2 = Multirate.filt( self, x2 )
    end
    statefulResult = [ y1; y2 ]

    @printf( "\n\tMultirate.filt interpolation. Piecewise for first %d inputs\n\t\t", length( x1 ) )
    Multirate.reset( self )
    y1 = similar( x, 0 )
    @time begin
        for i in 1:length(x1)
            append!( y1, Multirate.filt( self, x1[i:i] ) )
        end
        y2 = Multirate.filt( self, x2 )
    end
    piecewiseResult = [ y1; y2 ]

    if isapprox( baseResult, statelessResult, atol=atol ) && isapprox( baseResult, statefulResult, atol=atol ) && isapprox( piecewiseResult, baseResult, atol=atol )
        return true
    end

    display( [ (1:length(baseResult)) baseResult statelessResult statefulResult piecewiseResult (baseResult.-statelessResult) (baseResult.-statefulResult) (baseResult.-piecewiseResult) ] )

    return false
end

#==============================================================================#
#           ____ ____ ___     ____ ____ ____ ____ _  _ ___  _    ____          #
#           |__/ |__|  |      |__/ |___ [__  |__| |\/| |__] |    |___          #
#           |  \ |  |  |  .   |  \ |___ ___] |  | |  | |    |___ |___          #
#==============================================================================#

function test_rational( h, x, ratio )
    xLen       = length( x )
    hLen       = length( h )
    pivotPoint = min( rand(50:150, 1)[1], floor( Int, xLen/4 ))
    x1         = x[ 1 : pivotPoint ]
    x2         = x[ pivotPoint+1 : end ]
    upfactor   = numerator( ratio )
    downfactor = denominator( ratio )
    resultType = promote_type( eltype(h), eltype(x) )

    @printf( "\n\n" )
    @printf( "      ____ ____ ___ _ ____ _  _ ____ _    \n" )
    @printf( "      |__/ |__|  |  | |  | |\\ | |__| |    \n" )
    @printf( "      |  \\ |  |  |  | |__| | \\| |  | |___ \n" )
    @printf( "                                          \n" )
    @printf( "____ ____ ____ ____ _  _ ___  _    _ _  _ ____\n" )
    @printf( "|__/ |___ [__  |__| |\\/| |__] |    | |\\ | | __\n" )
    @printf( "|  \\ |___ ___] |  | |  | |    |___ | | \\| |__]\n" )
    @printf( "\n\nTesting rational resampling, h::%s, x::%s. xLen = %d, hLen = %d, ratio = %d//%d", string(typeof(h)), string(typeof(x)), xLen, hLen, upfactor, downfactor )

    @printf( "\n\tNaive rational resampling with Base.filt\n\t\t")
    @time begin
        xStuffed   = zeros( resultType, length(x) * upfactor )
        baseResult = Array{resultType}( ceil( Int, length(x) * ratio ))

        for n = 0:length(x)-1;
            xStuffed[ n*upfactor+1 ] = x[ n+1 ]
        end

        baseResult = Base.filt( h, one(eltype(h)), xStuffed );
        baseResult = [ baseResult[n] for n = 1:downfactor:length( baseResult ) ]
    end

    if method_exists( DSP.filt, ( typeof(h), typeof(x) ))
        @printf( "\n\tNaive rational resampling DSP.filt\n\t\t")
        @time begin
            xStuffed  = zeros( resultType, length(x) * upfactor )
            dspResult = Array{resultType}( ceil( Int, length(x) * ratio ))

            for n = 0:length(x)-1;
                xStuffed[ n*upfactor+1 ] = x[ n+1 ]
            end

            dspResult = DSP.filt( h, xStuffed );
            dspResult = [ dspResult[n] for n = 1:downfactor:length( dspResult ) ]
        end
    end

    @printf( "\n\tMultirate.filt( h, x, %d//%d )\n\t\t", upfactor, downfactor )
    @time statelessResult = Multirate.filt( h, x, ratio )

    @printf( "\n\tMultirate.filt rational resampling. length( x1 ) = %d, length( x2 ) = %d\n\t\t", length( x1 ), length( x2 ) )
    self = Multirate.FIRFilter( h, ratio )
    @time begin
        s1 = Multirate.filt( self, x1 )
        s2 = Multirate.filt( self, x2 )
    end
    statefulResult = [ s1; s2 ]

    @printf( "\n\tMultirate.filt rational. Piecewise for all %d inputs\n\t\t", length( x ) )
    self = Multirate.FIRFilter( h, ratio )
    y1 = similar( x, 0 )
    @time begin
        for i in 1:length(x)
            append!( y1, Multirate.filt( self, x[i:i] ) )
        end
    end
    piecewiseResult = y1

    if isapprox( baseResult, statelessResult ) && isapprox( baseResult, statefulResult ) && isapprox( baseResult, piecewiseResult )
        return true
    end

    display( [  baseResult statefulResult statefulResult piecewiseResult ] )

    return false
end

#==============================================================================#
#        ____ ____ ___      ____ ____ ____ ____ _  _ ___  _    _ _  _ ____     #
#        |__| |__/ |__]     |__/ |___ [__  |__| |\/| |__] |    | |\ | | __     #
#        |  | |  \ |__] .   |  \ |___ ___] |  | |  | |    |___ | | \| |__]     #
#==============================================================================#

function test_arbitrary( Th, x, resampleRate, numFilters; atol=0 )
    cutoffFreq      = 0.45
    transitionWidth = 0.05
    (hLen, β)       = Multirate.kaiserlength( transitionWidth, samplerate = numFilters )
    hLen            = ceil(Int,  hLen/numFilters  ) .* numFilters
    h               = Multirate.firdes( hLen, cutoffFreq, DSP.kaiser, samplerate = 32, beta = β ) .* numFilters
    h               = convert( Vector{Th}, h )

    @printf( "____ ____ ___      ____ ____ ____ ____ _  _ ___  _    _ _  _ ____\n" )
    @printf( "|__| |__/ |__]     |__/ |___ [__  |__| |\\/| |__] |    | |\\ | | __\n" )
    @printf( "|  | |  \\ |__] .   |  \\ |___ ___] |  | |  | |    |___ | | \\| |__]\n" )
    @printf( "\n\nh::%s, x::%s, rate = %f, N𝜙 = %d, xLen = %d, ", string(typeof(h)), string(typeof(x)), resampleRate, numFilters, length(x) )

    @printf( "\n\tNaive arbitrary resampling\n\t\t" )
    @time naiveResult = NaiveResamplers.naivefilt( h, x, resampleRate, numFilters )

    @printf( "\n\tStateless arbitrary resampling\n\t\t" )
    @time statelessResult = Multirate.filt( h, x, resampleRate, numFilters )

    @printf( "\n\tPiecewise arbitrary resampling\n\t\t" )
    self           = Multirate.FIRFilter( h, resampleRate, numFilters )
    piecewiseResult = eltype(x)[]
    sizehint!( piecewiseResult, ceil( Int, length(x)*resampleRate ) )
    @time for i in 1:length( x )
        thisY = filt( self, x[i:i] )
        append!( piecewiseResult, thisY )
    end

    commonLen = min( length(naiveResult), length(statelessResult), length(piecewiseResult) )

    resize!( naiveResult, commonLen )
    resize!( statelessResult, commonLen )
    resize!( piecewiseResult, commonLen )

    if isapprox( naiveResult, statelessResult, atol=atol ) && isapprox( naiveResult, piecewiseResult, atol=atol )
        return true
    end

    display( [  (1:commonLen) naiveResult statelessResult  piecewiseResult abs.(naiveResult.-statelessResult) abs.(naiveResult.-piecewiseResult) ] )

    display( [ (1:(commonLen-1)) diff(naiveResult) diff(statelessResult) diff(piecewiseResult) ] )
    return false
end

#==============================================================================#
#                  ____ _  _ _  _    ___ ____ ____ ___ ____                    #
#                  |__/ |  | |\ |     |  |___ [__   |  [__                     #
#                  |  \ |__| | \|     |  |___ ___]  |  ___]                    #
#==============================================================================#

function test_all()
    for interpolation in sort([1; unique(rand(2:32,8))] ),
            decimation in sort([1; unique(rand(2:32,8))] ),
                Th in [Float32, Float64],
                    Tx in [Float32, Float64, Complex64, Complex128]

        h     = rand(Th, rand(16:128,1)[1] )
        xLen  = Int(rand( 200:300, 1 )[1])
        xLen  = xLen-mod( xLen, decimation )
        x     = rand( Tx, xLen )
        ratio = interpolation//decimation

        if ratio == 1
            @test test_singlerate( h, x )
        end

        if decimation != 1
            @test test_decimation( h, x, decimation )
        end

        if interpolation != 1
            @test test_interpolation( h, x, interpolation )
        end

        if numerator(ratio) == interpolation && denominator(ratio) == decimation && numerator(ratio) != 1 && denominator(ratio) != 1
            @test test_rational( h, x, ratio )
            if Tx in [ Float32, Complex64 ]
                @test test_arbitrary( Th, x, Float64(ratio)+rand(), 32 )
            end
        end
    end
end

function test_nextphase()
    for interpolation in 1:8
        for decimation in 1:8
            ratio           = interpolation//decimation
            interpolation   = numerator(ratio)
            decimation      = denominator(ratio)
            x               = repmat( (1:interpolation), decimation )
            reference = [ x[n] for n in 1:decimation:length( x ) ]
            result = [ 1 ]
            for i in 2:interpolation
                append!( result, [ Multirate.nextphase( result[end], ratio ) ] )
            end
            @test isapprox( reference, result )
        end
    end
end

test_nextphase()
test_all()


text:
using Documenter, Aiyagari

PKG_HOME = joinpath(dirname(pathof(Aiyagari)), "..")

using LiterateWeave

EXMPL = joinpath(PKG_HOME, "examples")
TUTORIALS = joinpath(PKG_HOME, "docs", "src", "tutorials")
isdir(TUTORIALS) ? nothing : mkdir(TUTORIALS)

map(Aiyagari.examples) do exmpl
  literateweave(joinpath(EXMPL, exmpl*".jl"), out_path=TUTORIALS, doctype="github")
end

examples_pairs = [exmpl => "tutorials/" * exmpl*".md" for exmpl in Aiyagari.examples]

makedocs(;
    modules=[Aiyagari],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md";
        examples_pairs
    ],
    repo="https://gitlab.com/greimel/Aiyagari.jl/blob/{commit}{path}#L{line}",
    sitename="Aiyagari.jl",
    authors="Fabian Greimel <fabgrei@gmail.com>",
    assets=String[],
 )


text:
function bt02()
#
# Brugnano and Trigiante
#
    p = poly([ones(10);-1;-1;im;-im;2])
    z = [1.0 10; -1 2; im 1; -im 1; 2 1];

    p, PolyZeros(z)
end


text:
using FactorLoadingMatrices
using Random
using Plots
using Turing
# using Memoization and reverse-mode autodiff speeds up MCMC sampling
# *a lot* when fitting high-dimensional models
using Memoization
using ReverseDiff
Turing.setadbackend(:reversediff)

Random.seed!(1)

nx = 50
nfactor = 3
nobs = 300
σ = 0.5

vals = randn(nnz_loading(nx, nfactor))
L = loading_matrix(vals, nx, nfactor)

F = randn(nfactor, nobs)
X = L * F .+ σ*randn()

@model function FactorModel(X, nfactor)
    nx, nobs = size(X)
    # set priors for factors, loading matrix, and observation noise
    F ~ filldist(Normal(0, 1), nfactor, nobs)
    vals ~ filldist(Normal(0, 2), nnz_loading(nx, nfactor))
    σ ~ Exponential(1.0)

    L = loading_matrix(vals, nx, nfactor)
    X ~ arraydist(Normal.(L * F, σ))
end

mod = FactorModel(X, nfactor)
chn = sample(mod, NUTS(), 100)

vals_post = Array(group(chn, :vals))
L_post = [loading_matrix(v, nx, 3) for v in eachrow(vals_post)]

L_post_vm = varimax.(L_post)
plot(mean(L_post_vm), 1:nx, xerror=2std(L_post_vm), markerstrokecolor=1,
    layout=(1, 3), label="Model", title=["L[:, 1]" "L[:, 2]" "L[:, 3]"])
plot!(varimax(L), 1:nx, layout=(1, 3), label="Truth")

savefig("mcmc_factors.png")


text:
module TiledIteration

using OffsetArrays
using Base: tail, Indices, @propagate_inbounds
using Base.IteratorsMD: inc
using ArrayInterface

if VERSION < v"1.6.0-DEV.1174"
    _inc(state, iter) = inc(state, first(iter).I, last(iter).I)
else
    # https://github.com/JuliaLang/julia/pull/37829
    _inc(state, iter) = inc(state, iter.indices)
end

export TileIterator, EdgeIterator, SplitAxis, SplitAxes, padded_tilesize, TileBuffer, RelaxStride, RelaxLastTile

include("tileiterator.jl")

const L1cachesize = 2^15
const cachelinesize = 64

### EdgeIterator ###

struct EdgeIterator{N,UR1,UR2}
    outer::CartesianIndices{N,UR1}
    inner::CartesianIndices{N,UR2}
    function EdgeIterator{N,UR1,UR2}(outer::CartesianIndices{N}, inner::CartesianIndices{N}) where {N,UR1,UR2}
        ((first(inner) ∈ outer) & (last(inner) ∈ outer)) || throw(DimensionMismatch("$inner must be in the interior of $outer"))
        new(outer, inner)
    end
end
EdgeIterator(outer::CartesianIndices{N,UR1}, inner::CartesianIndices{N,UR2}) where {N,UR1,UR2} =
    EdgeIterator{N,UR1,UR2}(outer, inner)
EdgeIterator(outer::Indices{N}, inner::Indices{N}) where N =
    EdgeIterator(promote(CartesianIndices(outer), CartesianIndices(inner))...)

"""
    EdgeIterator(outer, inner)

A Cartesian iterator that efficiently visits sites that are in `outer`
but not in `inner`. This can be useful for calculating edge values
that may require special treatment or boundary conditions.
"""
EdgeIterator

Iterators.IteratorEltype(::Type{<:EdgeIterator}) = Iterators.HasEltype()

Base.eltype(::Type{EdgeIterator{N,UR1,UR2}}) where {N,UR1,UR2} = CartesianIndex{N}
Base.length(iter::EdgeIterator) = length(iter.outer) - length(iter.inner)

function Base.iterate(iter::EdgeIterator)
    iterouter = iterate(iter.outer)
    iterouter === nothing && return nothing
    item = nextedgeitem(iter, iterouter[2])
    item ∉ iter.outer && return nothing
    return item, item
end
function Base.iterate(iter::EdgeIterator, state)
    iterouter = iterate(iter.outer, state)
    iterouter === nothing && return nothing
    item = nextedgeitem(iter, iterouter[2])
    item.I[end] > last(iter.outer.indices[end]) && return nothing
    return item, item
end

@inline function nextedgeitem(iter::EdgeIterator, I::CartesianIndex)
    !(I ∈ iter.inner) && return I
    state = (last(iter.inner)[1], tail(I.I)...)
    newI = CartesianIndex(_inc(state, iter.outer))
    nextedgeitem(iter, newI)
end

Base.show(io::IO, iter::EdgeIterator) = print(io, "EdgeIterator(", iter.outer.indices, ", ", iter.inner.indices, ')')

### SplitAxis and SplitAxes

struct SplitAxis <: AbstractVector{UnitRange{Int}}
    splits::Vector{Int}
end

"""
    SplitAxis(ax::AbstractUnitRange, n::Real)

Split `ax` into `ceil(Int, n)` approximately equal-sized chunks. The first chunk is no larger than any other chunk,
and any fractional "deficit" in `n` will further shrink the first chunk.

This can be useful in splitting work across threads. When the first thread is responsible for assigning work to the others,
it's often useful to assign less work to it to account for the time spent scheduling.

# Examples

```jldoctest; setup=:(using TiledIteration)
julia> collect(SplitAxis(1:16, 4))
4-element Vector{UnitRange{$Int}}:
 1:4
 5:8
 9:12
 13:16

julia> collect(SplitAxis(1:16, 3.5))
4-element Vector{UnitRange{$Int}}:
 1:1
 2:6
 7:11
 12:16
```

In the latter case all the ranges except the first have length 5; consequently, only one element remains for the first chunk.
"""
function SplitAxis(ax::AbstractUnitRange{<:Integer}, n::Real)
    step = ceil(Int, length(ax)/n)
    # Give the smallest amount of work to thread 1, since often it is also scheduling the work for all
    # the other threads.
    SplitAxis(max.(first(ax)-1, collect(reverse(range(last(ax), step=-step, length=ceil(Int, n)+1)))))
end

Base.@propagate_inbounds Base.getindex(sax::SplitAxis, i::Int) = sax.splits[i]+1:sax.splits[i+1]

Base.size(sax::SplitAxis) = (length(sax.splits)-1,)

struct SplitAxes{N} <: AbstractVector{Tuple{UnitRange{Int},Vararg{UnitRange{Int},N}}}
    inner::NTuple{N,UnitRange{Int}}
    splitax::SplitAxis
end

"""
    SplitAxes(axs::NTuple{N,AbstractUnitRange}, n::Real)

Split `axs` into `ceil(Int, n)` approximately equal-sized chunks along the final dimension represented by `axs`.

See [`SplitAxis`](@ref) for further details.

# Examples

```jldoctest; setup=:(using TiledIteration)
julia> A = rand(3, 16);

julia> collect(SplitAxes(axes(A), 4))
4-element Vector{Tuple{UnitRange{$Int}, UnitRange{$Int}}}:
 (1:3, 1:4)
 (1:3, 5:8)
 (1:3, 9:12)
 (1:3, 13:16)
```
"""
SplitAxes(axs::Tuple{AbstractUnitRange,Vararg{AbstractUnitRange}}, n::Real) = SplitAxes{length(axs)-1}(Base.front(axs), SplitAxis(axs[end], n))

Base.@propagate_inbounds Base.getindex(saxs::SplitAxes, i::Int) = (saxs.inner..., saxs.splitax[i])

Base.size(saxs::SplitAxes) = size(saxs.splitax)

### Calculating the size of tiles ###

# If kernelsize-1 is the amount of padding, and s is the extra width of the tile along
# each axis, then the fraction of useful to total work is
# prod(s)/prod(s+kernelsize); this ratio is maximized if s is proportional to
# kernelsize.
"""
    padded_tilesize(T::Type, kernelsize::Dims, [ncache=2]) -> tilesize::Dims

Calculate a suitable tile size to approximately maximize the amount of
productive work, given a stencil of size `kernelsize`. The element
type of the array is `T`. Optionally specify `ncache`, the number of
such arrays that you'd like to have fit simultaneously in L1 cache.

This favors making the first dimension larger, since the first
dimension corresponds to individual cache lines.

# Examples
julia> padded_tilesize(UInt8, (3,3))
(768,18)

julia> padded_tilesize(UInt8, (3,3), 4)
(512,12)

julia> padded_tilesize(Float64, (3,3))
(96,18)

julia> padded_tilesize(Float32, (3,3,3))
(64,6,6)
"""
function padded_tilesize(::Type{T}, kernelsize::Dims, ncache = 2) where T
    nd = max(1, sum(x->x>1, kernelsize))
    # isbits(T) || return map(zero, kernelsize)
    # don't be too minimalist on the cache-friendly dim (use at least 2 cachelines)
    dim1minlen = 2*cachelinesize÷sizeof(T)
    psz = (max(kernelsize[1], dim1minlen), tail(kernelsize)...)
    L = sizeof(T)*prod(psz)
    # try to stay in L1 cache, but in the end we want a reasonably
    # favorable work ratio. f is the constant of proportionality in
    #      s+kernelsize ∝ kernelsize
    f = max(floor(Int, (L1cachesize/(ncache*L))^(1/nd)), 2)
    return _padded_tilesize_scale(f, psz)
end

@noinline _padded_tilesize_scale(f, psz) = map(x->x <= 1 ? x : f*x, psz) # see #15276

### Tile shaping and coordinate transformation

struct TileBuffer{T,N,P} <: DenseArray{T,N}
    view::OffsetArray{T,N,Array{T,N}}  # the currently-active view
    buf::Array{T,P}                    # the original backing buffer
end
viewtype(::Type{TileBuffer{T,N,P}}) where {T,N,P} = OffsetArray{T,N,Array{T,N}}
buftype(::Type{TileBuffer{T,N,P}}) where {T,N,P} = Array{T,P}
parent_type(::Type{OffsetArray{T,N,Array{T,N}}}) where {T,N} = Array{T,N}

"""
    TileBuffer(a, inds::Indices) -> v

Return a buffer-view `v` whose indices match `inds`, using the array
`a` for storage. `inds` does not necessarily have to match the size of
`a` (which allows tiles to be of different sizes, e.g., at the edges).
"""
@inline function TileBuffer(a::Array, inds::Indices)
    l = map(length, inds)
    tilev = OffsetArray(_tileview(a, l), inds)
    TileBuffer(tilev, a)
end

"""
    TileBuffer(T, inds::Indices) -> v

Return a TileBuffer, allocating a new backing array of element type `T`
and size determined by `inds`.
"""
function TileBuffer(::Type{T}, inds::Indices) where T
    l = map(length, inds)
    TileBuffer(Array{T}(undef, l), inds)
end

TileBuffer(tb::TileBuffer, inds::Indices) = TileBuffer(tb.buf, inds)

@inline function _tileview(a::Array, l::Dims)
    if size(a) == l
        return a
    else
        # returning a SubArray would not be type-stable, we must return an Array
        prod(l) > length(a) && throw(DimensionMismatch("array of size $(size(a)) is not adequate for a tile of size $l"))
        return unsafe_wrap(Array, pointer(a), l)
    end
end

Base.axes(tb::TileBuffer) = axes(tb.view)
Base.size(tb::TileBuffer) = size(tb.view)

@inline @propagate_inbounds Base.getindex(tb::TileBuffer{T,N}, I::Vararg{Int,N}) where {T,N} = tb.view[I...]

@inline @propagate_inbounds Base.setindex!(tb::TileBuffer{T,N}, val, I::Vararg{Int,N}) where {T,N} = tb.view[I...] = val

Base.pointer(tb::TileBuffer) = pointer(tb.view)
Base.strides(tb::TileBuffer) = strides(tb.view)
Base.parent(tb::TileBuffer) = tb.view

ArrayInterface.stride_rank(::Type{TB}) where TB<:TileBuffer = ArrayInterface.stride_rank(parent_type(viewtype(TB)))
ArrayInterface.contiguous_axis(::Type{TB}) where TB<:TileBuffer = ArrayInterface.contiguous_axis(parent_type(viewtype(TB)))
ArrayInterface.dense_dims(tb::TileBuffer) = ArrayInterface.dense_dims(parent(tb))
ArrayInterface.offsets(tb::TileBuffer) = ArrayInterface.offsets(parent(tb))

end # module


text:
module LotkaVolterra2dSingularSRK

    const Δt = 0.1
    const nt = 100000

    const PLOT_DIR = "figures"
    const SYMP_DIR = "symplecticity"

    using GeometricIntegrators

    using GeometricProblems.LotkaVolterra2dSingular
    using GeometricProblems.LotkaVolterra2dPlots

    set_config(:nls_atol, 1E-14)
    set_config(:nls_rtol, 1E-14)
    set_config(:nls_nmax, 100);

    include("common.jl")
    include("tableau_lists.jl")

end


text:
module ExtremeLearningMachine

import LinearAlgebra

include("util.jl")
include("activation_functions.jl")
include("neurons.jl")
include("hidden_layer.jl")
include("slfn.jl")
include("elm.jl")

export ELM,
       Linear, ReLU, Sigmoid, Square, Tanh,
       add_neurons!, add_data!, clear!, solve, predict

end # module


text:
function process_dependencies(request::HTTP.Request, state::HandlerState)
    return HTTP.Response(
        200,
        ["Content-Type" => "application/json"],
        body = state.cache.dependencies_json
        )
end


text:
# Much hackier / cleverer solution from reddit, porting to Julia

s = open("input.txt") do file
    readlines(file);
end

function DoTheMath(expression)
    newExpr = replace(expression, "*"=>"-");
    newExpr = replace(newExpr, "+"=>"÷");

    - = *;
    ÷ = +;
    return eval(Meta.parse(newExpr));
end

answer = sum(map(DoTheMath, s))


text:
module Examples

export DenseBNReLU, MLPEncoder, MLPDecoder, MLPAutoEncoder

include("mlp.jl")

end


text:
baremodule TwoBasedIndexing
using Main: Base
export twobased

for n = 1:5
    Base.eval(TwoBasedIndexing, :(getindex(A, $([symbol("i$i") for i = 1:n]...)) =
                                  Base.getindex(A, $([:($(symbol("i$i"))-1) for i = 1:n]...))))
end
getindex(A, inds...) = Base.getindex(A, map(x->x-1, inds)...)

for n = 1:5
    Base.eval(TwoBasedIndexing, :(setindex!(A, X, $([symbol("i$i") for i = 1:n]...)) =
                                  Base.setindex!(A, X, $([:($(symbol("i$i"))-1) for i = 1:n]...))))
end
setindex!(A, X, inds...) = Base.setindex!(A, X, map(x->x-1, inds)...)

function twobased(mod=Base.current_module())
    Base.eval(mod, quote
        const getindex = $getindex
        const setindex! = $setindex!
    end)
    nothing
end
end # module


text:
import SpecialFunctions: digamma, gamma
using Neighborhood: Theiler, KDTree, bulksearch

abstract type NearestNeighborEntropyEstimator <: EntropyEstimator end

function maximum_neighbor_distances(pts, est::EntropyEstimator)
    theiler = Theiler(est.w)
    est.w >= 0 || error("w, the number of neighbors to exclude, must be >= 0")
    tree = KDTree(pts) # Euclidean metric forced
    idxs, dists = bulksearch(tree, pts.data, NeighborNumber(est.k), theiler)
    # Distance to k-th nearest neighbor for each of the points
    return [d[end] for d in dists]
end

"Volume of a unit ball in R^d."
ball_volume(d::Int) = π^(d/2)/gamma((d/2)+1)

include("KozachenkoLeonenko.jl")
include("Kraskov.jl")


text:
module Data

using Random: shuffle!
using Base: @propagate_inbounds

include("dataloader.jl")
export DataLoader

## TODO for v0.13: remove everything below ##############
## Also remove the following deps:
## AbstractTrees, ZipFiles, CodecZLib

import ..Flux
import SHA

deprecation_message() = @warn("Flux's datasets are deprecated, please use the package MLDatasets.jl")

function deps(path...)
  if isnothing(@__DIR__) # sysimages
    joinpath("deps", path...)
  else
    joinpath(@__DIR__, "..", "..", "deps", path...)
  end
end

function download_and_verify(url, path, hash)
    tmppath = tempname()
    download(url, tmppath)
    hash_download = open(tmppath) do f
        bytes2hex(SHA.sha256(f))
    end
    if hash_download !== hash
        msg  = "Hash Mismatch!\n"
        msg *= "  Expected sha256:   $hash\n"
        msg *= "  Calculated sha256: $hash_download"
        error(msg)
    end
    mv(tmppath, path; force=true)
end

function __init__()
  mkpath(deps())
end

include("mnist.jl")
export MNIST

include("fashion-mnist.jl")
export FashionMNIST

include("cmudict.jl")
export CMUDict
using .CMUDict; export cmudict

include("tree.jl")
include("sentiment.jl")
export Sentiment

include("iris.jl")
export Iris

include("housing.jl")
export Housing

#########################################

end#module


text:
using DataFrames
using FileIO
using JSON
using Statistics
using VegaLite

rawdata = JSON.parsefile(joinpath(@__DIR__, "build/results.json"))

potatos = map(rawdata["potatos"]) do info
    result = info["result"]
    ntasks = info["ntasks"]
    impl = Symbol(info["impl"])
    trialid = info["trialid"]
    stats = result["stats"]
    npush = sum(stat["npush"] for stat in stats)
    npop = sum(stat["npop"] for stat in stats)
    t1min = minimum(stat["t1"] for stat in stats)
    t2max = maximum(stat["t2"] for stat in stats)
    duration = (t2max - t1min) / 1e9
    nops = npush + npop
    throughput = nops / duration
    return (; trialid, ntasks, impl, throughput, duration, nops, npush, npop)
end

fais = map(rawdata["fais"]) do info
    result = info["result"]
    ntasks = info["ntasks"]
    impl = :fai
    trialid = info["trialid"]
    stats = result["stats"]
    nops = sum(stat["nfai"] for stat in stats)
    t1min = minimum(stat["t1"] for stat in stats)
    t2max = maximum(stat["t2"] for stat in stats)
    duration = (t2max - t1min) / 1e9
    throughput = nops / duration
    return (; trialid, ntasks, impl, throughput, duration, nops)
end

datadf = vcat(DataFrame(potatos), DataFrame(fais); cols = :union)

impllabels = Dict(
    # :impl -> :Implementation (for plot)
    :base => "Base (Channel)",
    :dlcrq => "Dual LCRQ",
    :fai => "Hardware \"limit\" (FAI)",
)

datadf[!, :Implementation] = getindex.(Ref(impllabels), datadf.impl)

datadf[!, :diff_push_pop] = datadf.npush .- datadf.npop

plt = @vlplot(
    :point,
    x = {:ntasks, title = "Number of Tasks"},
    y = {:throughput, title = "Throughput [#OP/second]"},
    color = :Implementation,
    # title = "Hot Potato Benchmark",
    data = datadf,
)

save(joinpath(@__DIR__, "build/results.png"), plt)

summarydf = let
    idx = datadf.impl .!= :fai
    df = datadf[idx, :]

    sdf = unstack(
        combine(groupby(df, [:impl, :ntasks]), :throughput => median => :throughput),
        :impl,
        :throughput,
    )
    sdf[!, :speedup] = sdf.dlcrq ./ sdf.base
    sdf
end

plt_speedup = @vlplot(
    :point,
    x = {:ntasks, title = "Number of Tasks"},
    y = {:speedup, title = "Speedup [DLCRQ/Base]"},
    # title = "Hot Potato Benchmark (Speedup)",
    data = summarydf,
)

save(joinpath(@__DIR__, "build/speedup.png"), plt_speedup)

plt_trialid = @vlplot(
    :point,
    x = {:ntasks, title = "Number of Tasks"},
    y = {:throughput, title = "Throughput [#OP/second]"},
    color = {
        :trialid,
        scale = {scheme = :magma},
    },
    # title = "",
    data = datadf[datadf.impl .== :base, :],
)

save(joinpath(@__DIR__, "build/trial_dependencies.png"), plt_trialid)

plt_diff_push_pop = @vlplot(
    :point,
    x = :diff_push_pop,
    y = :throughput,
    color = {
        :ntasks,
        # "ntasks:o",
        scale = {scheme = :viridis},
    },
    shape = :impl,
    column = :impl,
    data = datadf,
)

save(joinpath(@__DIR__, "build/diff_push_pop.png"), plt_diff_push_pop)


text:
transD_GP.getchi2forall(optn)
gcf()
##
opt.xall[:] .= zall
transD_GP.plot_posterior(tempest, opt, burninfrac=0.5, figsize=(4,4),
    cmappdf="inferno")
ax = gcf().axes
ax[1].step(log10.(ρ[2:end]), z[2:end], color="w", alpha=0.4)
ax[1].invert_xaxis()
## nuisance histograms
transD_GP.plot_posterior(tempest, optn, burninfrac=0.5, nbins=50)
gcf()


text:
# Find the smallest order of magnitude that is larger than xspan This is a
# little opaque because I want to avoid assuming the log function is defined
# over typeof(xspan)
function bounding_order_of_magnitude(xspan)
    DT = typeof(xspan)
    one_dt = one(DT)

    a = 1
    step = 1
    while xspan < 10.0^a * one_dt
        a -= step
    end

    b = 1
    step = 1
    while xspan > 10.0^b * one_dt
        b += step
    end

    while a + 1 < b
        c = div(a + b, 2)
        if xspan < 10.0^c * one_dt
            b = c
        else
            a = c
        end
    end

    return b
end

# Empty catchall
optimize_ticks() = {}

# Find some reasonable values for tick marks.
#
# This is basically Wilkinson's ad-hoc scoring method that tries to balance
# tight fit around the data, optimal number of ticks, and simple numbers.
#
# Args:
#   x_min: minimum value occuring in the data.
#   x_max: maximum value occuring in the data.
#   Q: tick intervals and scores
#   k_min: minimum number of ticks
#   k_max: maximum number of ticks
#   k_ideal: ideal number of ticks
#   strict_span: true if no ticks should be outside [x_min, x_max].
#
# Returns:
#   A Float64 vector containing tick marks.
#
function optimize_ticks{T}(x_min::T, x_max::T; extend_ticks::Bool=false,
                           Q=[(1.0,1.0), (5.0, 0.9), (2.0, 0.7), (2.5, 0.5), (3.0, 0.2)],
                           k_min=2, k_max=10, k_ideal=5,
                           strict_span=false)

    Qv = [(float64(q[1]), float64(q[2])) for q in Q]
    optimize_ticks_typed(x_min, x_max, extend_ticks, Qv, k_min, k_max, k_ideal, strict_span)
end

function optimize_ticks_typed{T}(x_min::T, x_max::T, extend_ticks, Q::Vector{(Float64,Float64)}, k_min,
                           k_max, k_ideal, strict_span)
    one_t = one(T)
    if x_min == x_max
        R = typeof(1.0 * one_t)
        return R[x_min], x_min - one_t, x_min + one_t
    end

    const n = length(Q)

    # generalizing "order of magnitude"
    xspan = x_max - x_min
    z = bounding_order_of_magnitude(xspan)

    high_score = -Inf
    z_best = 0.0
    k_best = 0.0
    r_best = 0.0
    q_best = 0.0

    while k_max * 10.0^(z+1) * one_t > xspan
        for k in k_min:k_max
            for (q, qscore) in Q
                span = (k - 1) * q * 10.0^z * one_t
                if span < xspan
                    continue
                end

                r = ceil((x_max - span) / (q*10.0^z * one_t))
                while r*q*10.0^z * one_t <= x_min
                    has_zero = r <= 0 && abs(r) < k

                    # simplicity
                    s = has_zero ? 1.0 : 0.0

                    # granularity
                    g = 0 < k < 2k_ideal ? 1 - abs(k - k_ideal) / k_ideal : 0.0

                    # coverage
                    c = 1.5 * xspan/span

                    score = (1/4)g + (1/6)s + (1/3)c + (1/4)qscore

                    # strict limits on coverage
                    if strict_span && span > xspan
                        core -= 10000
                    elseif !strict_span && (span >= 2.0*xspan || span < xspan)
                        score -= 1000
                    end

                    if score > high_score
                        (q_best, r_best, k_best, z_best) = (q, r, k, z)
                        high_score = score
                    end
                    r += 1
                end
            end
        end
        z -= 1
    end

    span = q_best * 10.0^z_best * one_t
    if extend_ticks
        S = Array(typeof(1.0 * one_t), int(3 * k_best))
        for i in 0:(3*k_best - 1)
            S[i+1] = (r_best + i - k_best) * span
        end
        viewmin, viewmax = S[k_best + 1], S[2 * k_best]
    else
        S = Array(typeof(1.0 * one_t), int(k_best))
        for i in 0:(k_best - 1)
            S[i+1] = (r_best + i) * span
        end
        viewmin, viewmax = S[1], S[end]
    end

    S, viewmin, viewmax
end

function optimize_ticks(x_min::Date, x_max::Date; extend_ticks::Bool=false,
                        k_min=nothing, k_max=nothing,
                        scale=:auto)
    # This can be pretty simple. We are choosing ticks on one of three
    # scales: years, months, days.
    if year(x_max) - year(x_min) <= 1 && scale != :year
        if year(x_max) == year(x_min) && month(x_max) - month(x_min) <= 1 && scale != :month
            ticks = Date[]
            if x_max - x_min > Day(7) && scale != :day
                # This will probably need to be smarter
                push!(ticks, x_min)
                while true
                    next_month = Date(year(ticks[end]), month(ticks[end])) + Month(1)
                    while ticks[end] + Week(1) < next_month - Day(2)
                        push!(ticks, ticks[end] + Week(1))
                    end
                    push!(ticks, next_month)
                    if next_month >= x_max
                        break
                    end
                end
            else
                push!(ticks, x_min)
                while ticks[end] < x_max
                    push!(ticks, ticks[end] + day(1))
                end
            end

            viewmin, viewmax = ticks[1], ticks[end]
            ticks, viewmin, viewmax
        else
            ticks = Date[]
            push!(ticks, Date(year(x_min), month(x_min)))
            while ticks[end] < x_max
                push!(ticks, ticks[end] + Month(1))
            end
            viewmin, viewmax = ticks[1], ticks[end]

            ticks, x_min, x_max
        end
    else
        ticks, viewmin, viewmax =
            optimize_ticks(year(x_min), year(x_max), extend_ticks=extend_ticks)
        return Date[Date(round(y)) for y in ticks],
                   Date(round(viewmin)), Date(round(viewmax))
    end
end

# Generate ticks suitable for multiple scales.
function multilevel_ticks{T}(viewmin::T, viewmax::T;
                             scales=[0.5, 5.0, 10.0])

    ticks = Dict()
    for scale in scales
        ticks[scale] = optimize_ticks(viewmin, viewmax,
                                       k_min=2*scale,
                                       k_max=max(3, 10*scale),
                                       k_ideal=max(2, 15*scale))[1]
    end

    return ticks
end

function multilevel_ticks(viewmin::Date, viewmax::Date;
                          scales=[:year, :month, :day])
    span = viewmax - viewmin
    ticks = Dict()
    for scale in scales
        if scale == :year
            s = span / Day(360)
        elseif scale == :month
            s = span / Day(90)
        else
            s = span / Day(1)
        end

        ticks[s/20] = optimize_ticks(viewmin, viewmax, scale=scale)[1]
    end

    return ticks
end


text:
################################################################################
# VisualizationOptions
################################################################################

@with_kw mutable struct PlayerVisualizationOptions{T}
    # Options
	"Path to the material files."
	material_path::String=joinpath(module_dir(), "resources", "material")

	"Path to the object files."
	object_path::String=joinpath(module_dir(), "resources", "object", "car_geometry.obj")

	"Scaling of the car 3D mesh."
	car_scale::T=1.0

	"Car offset that puts the car is the correct orientation/position."
	car_offset::AbstractVector=[
			LinearMap(AngleAxis(-0.012+pi/200, 0, 1, 0)),
			LinearMap(AngleAxis(pi/2, 0, 0, 1)),
			LinearMap(AngleAxis(pi/2, 1, 0, 0))]

	"Height of the collision avoidance cylinder."
	cylinder_height::T=0.02

	"Fading of the cylinder"
	α::T=0.9

	"Colors."
	colors::Dict{Symbol,RGB}=Dict(
		:cornflowerblue => RGB(100/255, 149/255, 237/255),
		:orange         => RGB(255/255, 125/255,   0/255),
		:forestgreen    => RGB( 34/255, 139/255,  34/255),
		:red            => RGB(255/255,   0/255,   0/255),
		:yellow         => RGB(255/255, 255/255,   0/255),
		:gray           => RGB(127/255, 127/255, 127/255),
		:lime           => RGB(  0/255, 255/255,   0/255),
		)
end

@with_kw mutable struct WaypointVisualizationOptions{T}
    # Options
	"Height of the waypoint cylinder."
	waypoint_height::T=0.015

	"Radius of the waypoint cylinder."
	waypoint_radius::T=0.013

	"Fading of the cylinder"
	α::T=0.9

	"Colors."
	colors::Dict{Symbol,RGB}=Dict(
		:cornflowerblue => RGB(100/255, 149/255, 237/255),
		:orange         => RGB(255/255, 125/255,   0/255),
		:forestgreen    => RGB( 34/255, 139/255,  34/255),
		:red            => RGB(255/255,   0/255,   0/255),
		:yellow         => RGB(255/255, 255/255,   0/255),
		:gray           => RGB(127/255, 127/255, 127/255),
		:lime           => RGB(  0/255, 255/255,   0/255),
		)
end

@with_kw mutable struct LineVisualizationOptions{T}
    # Options

	"Line width."
	line_width::T=8.0

	"Fading of the line"
	α::T=0.9

	"Colors."
	colors::Dict{Symbol,RGB}=Dict(
		:cornflowerblue => RGB(100/255, 149/255, 237/255),
		:orange         => RGB(255/255, 125/255,   0/255),
		:forestgreen    => RGB( 34/255, 139/255,  34/255),
		:red            => RGB(255/255,   0/255,   0/255),
		:yellow         => RGB(255/255, 255/255,   0/255),
		:gray           => RGB(127/255, 127/255, 127/255),
		:lime           => RGB(  0/255, 255/255,   0/255),
		)
end

function get_car_offset(vis_opts::PlayerVisualizationOptions)\
	car_scale = vis_opts.car_scale
	car_offset = vis_opts.car_offset
	offset = compose(compose(compose(
			Translation(0., 0., 0.03*car_scale),
			car_offset[1]),
			car_offset[2]),
			car_offset[3])
	return offset
end

################################################################################
# set_player!
################################################################################

function set_player!(vis::Visualizer, player::Player{T};
    vis_opts::PlayerVisualizationOptions=PlayerVisualizationOptions{T}(),
	color=:orange, key::Int=0) where {T}
    clean!(vis)
	mtlpath = joinpath(vis_opts.material_path , "$(string(color)).mtl")
	objpath = vis_opts.object_path
	col = vis_opts.colors[color]
	α = vis_opts.α
	ch = vis_opts.cylinder_height
	player_path = "env/player_key$key/player$(player.id)"
	offset = get_car_offset(vis_opts)

	# Add car object
	car = ModifiedMeshFileObject(objpath, mtlpath, scale=0.010*vis_opts.car_scale)
	setobject!(vis[player_path*"/car"], car)
	settransform!(vis[player_path*"/car"], offset)
	# Add cylinder
	cyl = Cylinder(Point(0.,0,0), Point(0,0,ch), player.r_col)
	mat = MeshPhongMaterial(color=RGBA(col, α))
	setobject!(vis[player_path*"/cyl"], cyl, mat)
    return nothing
end

function set_player!(vis::Visualizer, player::Vector{Player{T}};
    vis_opts::PlayerVisualizationOptions=PlayerVisualizationOptions{T}(),
	color::Vector{Symbol}=Vector{Symbol}(), key::Int=0) where {T}
	length(color) == 0 ? color = collect(keys(vis_opts.colors)) : nothing
	p = length(player)
	c = length(color)
	for i = 1:p
		set_player!(vis, player[i], vis_opts=vis_opts, color=color[1+(i-1)%c], key=key)
	end
	return nothing
end

################################################################################
# build_waypoint!
################################################################################

function build_waypoint!(vis::Visualizer, player::Player{T}, N::Int;
    vis_opts::WaypointVisualizationOptions=WaypointVisualizationOptions{T}(),
	color=:orange, key::Int=0) where {T}
    clean!(vis)
	col = vis_opts.colors[color]
	α = vis_opts.α
	wh = vis_opts.waypoint_height
	wr = vis_opts.waypoint_radius
	waypoint_path = "env/waypoint_key$key/player$(player.id)"

	# Add waypoint
	for k = 1:N
		wpt = Cylinder(Point(0.,0,0), Point(0,0,wh), wr)
		mat = MeshPhongMaterial(color=RGBA(col, α))
		setobject!(vis[waypoint_path*"/$k"], wpt, mat)
	end
    return nothing
end

function build_waypoint!(vis::Visualizer, player::Vector{Player{T}}, N::Int;
    vis_opts::WaypointVisualizationOptions=WaypointVisualizationOptions{T}(),
	color::Vector{Symbol}=Vector{Symbol}(), key::Int=0) where {T}
	length(color) == 0 ? color = collect(keys(vis_opts.colors)) : nothing
	p = length(player)
	c = length(color)
	for i = 1:p
		build_waypoint!(vis, player[i], N, vis_opts=vis_opts, color=color[1+(i-1)%c], key=key)
	end
	return nothing
end

################################################################################
# set_line_traj!
################################################################################

function set_line_traj!(vis::Visualizer, model::AbstractGameModel, player::Vector{Player{T}}, traj;
    vis_opts::LineVisualizationOptions=LineVisualizationOptions{T}(),
	color::Vector{Symbol}=Vector{Symbol}(), key::Int=0) where {T}
	length(color) == 0 ? color = collect(keys(vis_opts.colors)) : nothing
	p = length(player)
	c = length(color)
	for i = 1:p
		set_line_traj!(vis, model, player[i], traj, vis_opts=vis_opts,
			color=color[1+(i-1)%c], key=key)
	end
	return nothing
end

function set_line_traj!(vis::Visualizer, model::AbstractGameModel, player::Player{T}, traj;
    vis_opts::LineVisualizationOptions=LineVisualizationOptions{T}(),
	color=:orange, key::Int=0) where {T}
    clean!(vis)
	col = vis_opts.colors[color]
	α = vis_opts.α
	lw = vis_opts.line_width
	line_path = "env/line_key$key/player$(player.id)"

	# Add line
	mat = LineBasicMaterial(color=color=RGBA(col, α), linewidth=lw)
	N = length(traj)
	vs = [standardize(model, Algames.state(traj[k]))[player.id] for k=1:N]
	points = [[vs[k].x; vs[k].y; 0.01] for k=1:N]
	for k = 1:N
		setobject!(vis[line_path], Object(PointCloud(points), mat, "Line"))
	end
    return nothing
end


text:
using Test
using SeisRequests
using Seis
using Dates: DateTime
import QuakeML

@testset "get_events" begin
    server = "http://service.example.com"
    SUCCESS = SeisRequests.CODE_SUCCESS

    @testset "parse_events" begin
        @testset "No response" begin
            request = FDSNEvent()
            response = HTTP.Messages.Response(request.nodata)
            @test SeisRequests.parse_event_response(Float64, request,
                response, server) == []

            @testset "Element type $T" for T in (Float32, Float64)
                @test SeisRequests.parse_event_response(T, request, response,
                    server) isa Vector{Seis.GeogEvent{T}}
            end
            response.body = [1]
            @test (@test_logs (:warn,
                            "unexpected data in response reporting no data"
                ) SeisRequests.parse_event_response(
                            Float64, request, response, server) == [])
            request = FDSNEvent(nodata=404)
            response = HTTP.Messages.Response(request.nodata)
            @test SeisRequests.parse_event_response(Float64, request, response,
                server) == []
        end

        @testset "'Success' but no data" begin
            request = FDSNEvent()
            response = HTTP.Messages.Response(SUCCESS)
            @test_throws ErrorException SeisRequests.parse_event_response(Float64,
                request, response, server)
        end

        @testset "Content-Type" begin
            @testset "Text" begin
                request = FDSNEvent(format="text")
                body = "usp000jv5f|2012-11-07T16:35:46.930|13.988|-91.895|24|us|us|us|usp000jv5f|mww|7.4|us|offshore Guatemala|earthquake\n"

                @testset "Wrong" begin
                    response_correct = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"text/plain"), body=body)
                    response_type_wrong = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"WRONG"), body=body)
                    @test (@test_logs (:warn, "content type of response is \"WRONG\", " *
                             "not \"text/plain\" as expected"
                        ) SeisRequests.parse_event_response(
                            Float64, request, response_type_wrong, server)) ==
                        SeisRequests.parse_event_response(Float64, request,
                            response_correct, server)
                end

                @testset "Empty" begin
                    response_correct = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"text/plain"), body=body)
                    response_type_empty = HTTP.Messages.Response(SUCCESS,
                        Dict(), body=body)
                    @test SeisRequests.parse_event_response(Float64,
                        request, response_type_empty, server) ==
                        SeisRequests.parse_event_response(Float64,
                            request, response_correct, server)
                end
            end
            @testset "XML" begin
                request = FDSNEvent()
                body = """<?xml version="1.0" encoding="UTF-8"?>
                    <quakeml xmlns="http://quakeml.org/xmlns/quakeml/1.2">
                      <eventParameters publicID="smi:local/b3918029-d67d-4a49-9b40-128a53446e0d">
                        <event publicID="smi:local/bbc44eba-f0d3-40bc-8402-798d7cdfa153">
                          <magnitude publicID="smi:local/fbb0f09d-e23a-4a2d-ac18-caeaa61d671d">
                            <mag>
                              <value>1.0</value>
                            </mag>
                          </magnitude>
                          <origin publicID="smi:local/c61861b7-eb6b-4141-8645-98dc88a71ce7">
                            <time>
                              <value>2020-05-21T10:39:23.469</value>
                            </time>
                            <longitude>
                              <value>0.0</value>
                            </longitude>
                            <latitude>
                              <value>1.0</value>
                            </latitude>
                          </origin>
                        </event>
                      </eventParameters>
                    </quakeml>
                    """
                @testset "Wrong" begin
                    response_correct = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"application/xml"), body=body)
                    response_type_wrong = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"WRONG"), body=body)
                    @test (@test_logs (:warn, "content type of response is \"WRONG\", " *
                                                 "not \"application/xml\" as expected"
                            ) SeisRequests.parse_event_response(
                                Float64, request, response_type_wrong, server)) ==
                        SeisRequests.parse_event_response(Float64, request,
                            response_correct, server)
                end
                @testset "Empty" begin
                    response_correct = HTTP.Messages.Response(SUCCESS,
                        Dict("Content-Type"=>"application/xml"), body=body)
                    response_type_empty = HTTP.Messages.Response(SUCCESS,
                        Dict(), body=body)
                    @test SeisRequests.parse_event_response(
                            Float64, request, response_type_empty, server) ==
                        SeisRequests.parse_event_response(Float64, request,
                            response_correct, server)
                end
            end
        end

        @testset "Text format" begin
            request = FDSNEvent(format="text")
            body = """
                usp000jv5f|2012-11-07T16:35:46.930|13.988|-91.895|24|us|us|us|usp000jv5f|mww|7.4|us|offshore Guatemala|earthquake
                usp000juhz|2012-10-28T03:04:08.820|52.788|-132.101|14|us|us|us|usp000juhz|mww|7.8|us|Haida Gwaii, Canada|earthquake
                """
            @testset "Eltype $T" for T in (Float32, Float64)
                response = HTTP.Messages.Response(SUCCESS, body)
                out = SeisRequests.parse_event_response(T, request, response, server)
                @test out isa Vector{Seis.GeogEvent{T}}
                @test length(out) == 2
                @test out[1] == Seis.GeogEvent{T}(id="usp000jv5f",
                    time=DateTime("2012-11-07T16:35:46.930"), lat=13.988,
                    lon=-91.895, dep=24, meta=Dict(
                        :author=>"us", :catalog=>"us", :contributor=>"us",
                        :contributor_id=>"usp000jv5f", :mag_type=>"mww",
                        :mag=>7.4, :mag_author=>"us",
                        :location_name=>"offshore Guatemala",
                        :type=>"earthquake", :server=>server))
                @test out[2] == Seis.GeogEvent{T}(id="usp000juhz",
                    time=DateTime("2012-10-28T03:04:08.820"), lat=52.788,
                    lon=-132.101, dep=14, meta=Dict(
                        :author=>"us", :catalog=>"us", :contributor=>"us",
                        :contributor_id=>"usp000juhz", :mag_type=>"mww",
                        :mag=>7.8, :mag_author=>"us",
                        :location_name=>"Haida Gwaii, Canada",
                        :type=>"earthquake", :server=>server))
            end
        end

        @testset "QuakeML format" begin
            request = FDSNEvent()
            # Fake mashup of IRIS and ISC response.
            # First event is from IRIS; second from ISC
            body = """
                <?xml version="1.0" encoding="UTF-8"?>
                <q:quakeml xmlns:q="http://quakeml.org/xmlns/quakeml/1.2"
                    xmlns:iris="http://service.iris.edu/fdsnws/event/1/"
                    xmlns="http://quakeml.org/xmlns/bed/1.2"
                    xmlns:xsi="http://www.w3.org/2000/10/XMLSchema-instance"
                    xsi:schemaLocation="http://quakeml.org/schema/xsd http://quakeml.org/schema/xsd/QuakeML-1.2.xsd">
                <eventParameters publicID="smi:service.iris.edu/fdsnws/event/1/query">

                    <!-- IRIS event -->
                    <event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=11041250">
                        <type>earthquake</type>
                        <description xmlns:iris="http://service.iris.edu/fdsnws/event/1/" iris:FEcode="111">
                            <type>Flinn-Engdahl region</type>
                            <text>NORTHERN PERU</text>
                        </description>
                        <preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=193627740</preferredMagnitudeID>
                        <preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=38683486</preferredOriginID>
                        <!-- Fake non-preferred origin -->
                        <origin publicID="smi:example.com/fake/event">
                            <time>
                                <value>3000-01-01T00:00:00.000</value>
                            </time>
                            <creationInfo>
                                <author>Me</author>
                            </creationInfo>
                            <latitude>
                                <value>0</value>
                            </latitude>
                            <longitude>
                                <value>0</value>
                            </longitude>
                            <depth>
                                <value>700000.0</value>
                            </depth>
                        </origin>

                        <!-- End of fake origin -->
                        <origin xmlns:iris="http://service.iris.edu/fdsnws/event/1/"
                                publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=38683486"
                                iris:contributorOriginId="pt19146001" iris:contributor="us"
                                iris:contributorEventId="us60003sc0,at00ps3pco,pt19146001"
                                iris:catalog="NEIC PDE">
                            <time>
                                <value>2019-05-26T07:41:15.058</value>
                            </time>
                            <creationInfo>
                                <author>at,pt,us</author>
                            </creationInfo>
                            <latitude>
                                <value>-5.8132</value>
                            </latitude>
                            <longitude>
                                <value>-75.2775</value>
                            </longitude>
                            <depth>
                                <value>122400.0</value>
                            </depth>
                        </origin>

                        <!-- Fake non-preferred magnitude -->
                        <magnitude publicID="smi:example.com/fake/magnitude">
                            <mag>
                                <value>15</value>
                            </mag>
                            <type>mb</type>
                            <creationInfo>
                                <author>Me</author>
                            </creationInfo>
                        </magnitude>
                        <!-- End of fake magnitude -->

                        <magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=193627740">
                            <mag>
                                <value>8.0</value>
                            </mag>
                            <type>Mww</type>
                            <creationInfo>
                                <author>us</author>
                            </creationInfo>
                        </magnitude>
                    </event>

                    <!-- ISC event -->
                    <event publicID="smi:ISC/evid=618237507">
                      <preferredOriginID>smi:ISC/origid=613854533</preferredOriginID>
                      <description>
                        <text>Mozambique Channel</text>
                        <type>Flinn-Engdahl region</type>
                      </description>
                      <type>earthquake</type>
                      <typeCertainty>known</typeCertainty>
                      <origin publicID="smi:ISC/origid=613854533">
                        <time>
                          <value>2018-02-16T12:02:45.80Z</value>
                          <uncertainty>5.66</uncertainty>
                        </time>
                        <latitude>
                          <value>-25.4270</value>
                        </latitude>
                        <longitude>
                          <value>38.7170</value>
                        </longitude>
                        <depth>
                          <value>5000.0</value>
                        </depth>
                        <quality>
                          <associatedStationCount>5</associatedStationCount>
                          <standardError>18.0000</standardError>
                          <azimuthalGap>349.000</azimuthalGap>
                        </quality>
                        <creationInfo>
                          <author>NAM</author>
                        </creationInfo>
                        <originUncertainty>
                          <preferredDescription>horizontal uncertainty</preferredDescription>
                          <minHorizontalUncertainty>999900</minHorizontalUncertainty>
                          <maxHorizontalUncertainty>999900</maxHorizontalUncertainty>
                          <azimuthMaxHorizontalUncertainty>90</azimuthMaxHorizontalUncertainty>
                        </originUncertainty>
                      </origin>
                      <magnitude publicID="smi:ISC/magid=619919995">
                        <mag>
                          <value>9.90</value>
                        </mag>
                        <type>MD</type>
                        <originID>smi:ISC/origid=613854533</originID>
                        <creationInfo>
                          <author>NAM</author>
                        </creationInfo>
                      </magnitude>
                    </event>
                </eventParameters>
                </q:quakeml>
                """
            @testset "Eltype $T" for T in (Float32, Float64)
                response = HTTP.Messages.Response(SUCCESS, body)
                out = SeisRequests.parse_event_response(T, request, response, server)
                @test out isa Vector{Seis.GeogEvent{T}}
                @test length(out) == 2
                @test out[1] == Seis.GeogEvent{T}(lon=-75.2775, lat=-5.8132,
                    dep=122.4, time=DateTime(2019, 05, 26, 07, 41, 15, 058),
                    id="smi:service.iris.edu/fdsnws/event/1/query?originid=38683486",
                    meta=Dict(
                        :author=>"at,pt,us",
                        :description=>"NORTHERN PERU (Flinn-Engdahl region)",
                        :mag_type=>"Mww", :mag=>8.0, :mag_author=>"us",
                        :type=>"earthquake",
                        :quakeml=>QuakeML.readstring(body).event[1],
                        :server=>server))
                @test out[2] == Seis.GeogEvent{T}(lon=38.7170, lat=-25.427,
                    dep=5.0, time=DateTime(2018, 02, 16, 12, 02, 45, 800),
                    id="smi:ISC/origid=613854533",
                    meta=Dict(
                        :author=>"NAM",
                        :description=>"Mozambique Channel (Flinn-Engdahl region)",
                        :mag_type=>"MD", :mag=>9.9, :mag_author=>"NAM",
                        :type=>"earthquake",
                        :quakeml=>QuakeML.readstring(body).event[2],
                        :server=>server))
            end
        end
    end

    @testset "Real requests" begin
        @testset "Verbosity" begin
            @test_logs (:info, "Request status: Successful request, results follow"
                ) get_events(minmagnitude=8, starttime="2018-01-01", verbose=true)
            @test_logs get_events(minmagnitude=8, starttime="2018-01-01", verbose=false)
        end

        @testset "Eltype" begin
            @testset "$T" for T in (Float32, Float64)
                @test get_events(starttime=1990-01-01, endtime=1991-01-01,
                    minmagnitude=7.5, T=T, verbose=false) isa Vector{Seis.GeogEvent{T}}
            end
            @testset "Default eltype" begin
                @test get_events(starttime="2018-01-01", endtime="2018-01-02",
                    minmagnitude=4.7, verbose=false) isa Vector{Seis.GeogEvent{Float64}}
            end
        end

        @testset "Server" begin
            @test first(get_events(starttime="2000-01-01", endtime="2001-01-01",
                    minmagnitude=7, verbose=false, server="ISC")
                ).meta.server == "ISC"
        end

        @testset "Simple" begin
            evts = get_events(starttime="2019-01-01", endtime="2020-01-01",
                minmagnitude=7.5, format="text", verbose=false)
            @test length(evts) == 3
            @test evts[1] == Event(lon=-75.2775, lat=-5.8132, dep=122.4,
                time=DateTime(2019, 05, 26, 07, 41, 15),
                id="11041250",
                meta=Dict(
                    :catalog=>"NEIC PDE",
                    :contributor_id=>"us60003sc0,at00ps3pco,pt19146001",
                    :contributor=>"us", :author=>"at,pt,us",
                    :mag_type=>"Mww", :mag_author=>"us", :mag=>8.0,
                    :location_name=>"NORTHERN PERU", :server=>"IRIS"))
        end

        @testset "Station" begin
            date1, date2 = DateTime(1990), DateTime(1990, 2)
            sta = Station(lon=0, lat=0, meta=Dict(:startdate=>date1, :enddate=>date2))
            evts = get_events(sta, maxradius=90, minmagnitude=5, format="text",
                verbose=false)
            @test all(e -> distance_deg(e, sta) <= 90, evts)
            @test all(e -> date1 <= e.time <= date2, evts)
        end
    end
end


text:
using Bridge, ConjugatePriors, Plots, Distributions, LaTeXStrings

srand(5);
pgfplots()

simname = :gammachain

#joinpath(ENV["BRIDGE_OUTDIR"],  "$simname$(simid)paths.jld2")

Nf = 800001 # divisible by 25, 50, 100
step = [25, 50, 100, 200][4] # subsampling stepsize
H = [20,40,80,160,320][4]

skipp = 100

s1(t) = 3/2 + sin(2*(4t-2)) + 2exp(-16((4t-2)).^2) # called s_1 in the paper
b(x) = -10x + 20. # called b_1 in the paper

s = s1
dtf = 1/(Nf-1)
ttf = 0:dtf:1

N = (Nf-1) ÷ step + 1 # number of observations: n
td = (0:dtf:1)[1:step:end]
dt = td[2] - td[1]

M = (length(td)-1) ÷ H # observed intervals per bin
pri = InverseGamma(0.1, 0.1)

r = length(td)-1 - H*M
assert(H*M == length(td)-1)

a1 = 1/10
a, aζ = 1000.0/sqrt(H), 1000.0/sqrt(H)

import Bridge.piecewise

function piecewise(Y::SamplePath, tend = Y.tt[end])
    tt = [Y.tt[1]]
    n = length(Y.yy)
    append!(tt, repeat(Y.tt[2:n], inner=2))
    push!(tt, tend)
    tt, repeat(Y.yy, inner=2)
end
function piecewise(X::Bridge.EstSamplePath, tend = X.tt[end])
    Bridge.EstSamplePath(Bridge.piecewise(SamplePath(X.tt, X.yy), tend)...,
        Bridge.piecewise(SamplePath(X.tt, X.vv), tend)[2])
end

@recipe function plot(X::Bridge.EstSamplePath{Float64}; quantile = 0.95)
    q = Distributions.quantile(Normal(), 1 - (1 - quantile)/2)
    if Plots.backend() != Plots.PGFPlotsBackend()
        ribbon --> q*sqrt.(X.vv)
        X.tt, X.yy
    else
        ribbon --> (NaN,NaN)
        linewidth --> 0
        vcat(X.tt, reverse(X.tt)), vcat(X.yy + q*sqrt.(X.vv), reverse(X.yy - q*sqrt.(X.vv)))
    end
end

@recipe function plot(X::Bridge.SamplePathBand)
    if Plots.backend() != Plots.PGFPlotsBackend()
        error("todo")
    else
        ribbon --> (NaN,NaN)
        linewidth --> 0
        vcat(X.tt, reverse(X.tt)), vcat(X.uu, reverse(X.ll))
    end
end

"""
White noise on the grid tt
"""
dW(tt) = randn(length(tt)-1).*sqrt.(diff(tt))

"""
Finite difference operator
"""
function diff0(x)
    y = similar(x)
    y[1] = zero(x)
    for i in 2:length(x)
        y[i] = x[i]-x[i-1]
    end
    y
end

"""
Euler scheme with start in y1, space dependent drift b and time dependent dispersion coefficient
"""
function euler1(tt, y1, b, sigma, yy = zeros(tt), dW = dW(tt))
    yy[1] = y1
    for i in 2:length(yy)
        yy[i] = yy[i-1] + b(yy[i-1]) * (tt[i]-tt[i-1]) + sigma(tt[i-1])*dW[i-1]
    end
    yy
end


text:
module SPF

using HTTP
using XLSX
using UrlDownload
using DataFrames
using Dates

export
	get_additional_cpie10,
	get_anxious_index,
	get_dispersion,
	get_individual,
	get_inflation,
	get_mean,
	get_median,
	get_microdata,
	get_prob

include("url.jl")
include("get_additional_cpie10.jl")
include("get_anxious_index.jl")
include("get_dispersion.jl")
include("get_individual.jl")
include("get_inflation.jl")
include("get_mean.jl")
include("get_median.jl")
include("get_microdata.jl")
include("get_prob.jl")

#TODO Quarter(1) is introduced after julia v1.5.2 - change in future versions
#TODO add checks wether survey exists
#TODO write proper documentation

end


text:
"""
    @with_kw mutable struct ECRRepo
        repositoryArn::Union{Missing, String} = missing
        registryId::Union{Missing, String} = missing
        repositoryName::Union{Missing, String} = missing
        repositoryUri::Union{Missing, String} = missing
        createdAt::Union{Missing, String} = missing
        imageTagMutability::Union{Missing, String} = missing
        imageScanningConfiguration::Union{Missing, Any} = missing
        encryptionConfiguration::Union{Missing, Any} = missing
        exists::Bool = true
    end

Represents an AWS ECR (Elastic Container Registry) Repo. Should not be instantiated directly.
If `exists` is `true`, then the image is assumed to exit and so should be visible from utilities
such as `docker image ls`.
"""
@with_kw mutable struct ECRRepo
  repositoryArn::Union{Missing, String} = missing
  registryId::Union{Missing, String} = missing
  repositoryName::Union{Missing, String} = missing
  repositoryUri::Union{Missing, String} = missing
  createdAt::Union{Missing, String} = missing
  imageTagMutability::Union{Missing, String} = missing
  imageScanningConfiguration::Union{Missing, Any} = missing
  encryptionConfiguration::Union{Missing, Any} = missing
  exists::Bool = true
end
StructTypes.StructType(::Type{ECRRepo}) = StructTypes.Mutable()
Base.:(==)(a::ECRRepo, b::ECRRepo) = a.repositoryUri == b.repositoryUri

"""
    get_ecr_repo(image::LocalImage)::Union{Nothing, ECRRepo}

Queries AWS and returns an `ECRRepo` instance that is associated with the passed `local_image`.
Returns `nothing` if one cannot be found.
"""
function get_ecr_repo(local_image::LocalImage)::Union{Nothing, ECRRepo}
  all_repos = get_all_ecr_repos()
  image_full_name = get_image_full_name(local_image)
  index = findfirst(repo -> repo.repositoryUri == image_full_name, all_repos)
  isnothing(index) ? nothing : all_repos[index]
end

"""
    get_all_ecr_repos(jot_generated_only::Bool = true)::Vector{ECRRepo}

Returns a vector of `ECRRepo`s, representing all AWS-hosted ECR Repositories.

`jot_generated_only` specifies whether to filter for jot-generated repos only.
"""
function get_all_ecr_repos(jot_generated_only::Bool = false)::Vector{ECRRepo}
  all_repos_json = readchomp(`aws ecr describe-repositories`)
  all = JSON3.read(all_repos_json, Dict{String, Vector{ECRRepo}})
  all_repos = all["repositories"]
  jot_generated_only ? filter(is_jot_generated, all_repos) : all_repos
end

"""
    get_ecr_repo(repo_name::String)::Union{Nothing, ECRRepo}

Queries AWS and returns an `ECRRepo` instance with the passed `repo_name`.
Returns `nothing` if one cannot be found.
"""
function get_ecr_repo(repo_name::String)::Union{Nothing, ECRRepo}
  all_repos = get_all_ecr_repos()
  index = findfirst(repo -> repo.repositoryName == repo_name, all_repos)
  isnothing(index) ? nothing : all_repos[index]
end

function create_ecr_repo(image::LocalImage)::ECRRepo
  image.exists || error("Image does not exist")
  labels = get_labels(image)
  create_script = get_create_ecr_repo_script(
                                             get_lambda_name(image),
                                             get_aws_region(image),
                                             labels,
                                            )
  repo_json = readchomp(`bash -c $create_script`)
  @debug repo_json
  JSON3.read(repo_json, Dict{String, ECRRepo})["repository"]
end

"""
    delete!(repo::ECRRepo)

Removes the passed `ECRRepo` instance from AWS ECR, and sets the `exists` attribute to `false` to
indicate it no longer exists.
"""
function delete!(repo::ECRRepo)
  repo.exists || error("Repo does not exist")
  delete_script = get_delete_ecr_repo_script(repo.repositoryName)
  run(`bash -c $delete_script`)
  repo.exists = false
end


text:
using Combinatorics: with_replacement_combinations
using StatsBase: countmap

function compute(n::Int, m::Int)::Int
    factorials = map(factorial, 0:9)
    chains = Dict{Int, Int}()
    result = 0
    for combination ∈ with_replacement_combinations(9:-1:0, n)
        for _ ∈ 1:count(isequal(0), combination) + 1
            number = parse(Int, join(combination))
            chain = [number]
            term = sum(factorials[parse.(Int, collect(string(number))) .+ 1])
            while term ∉ chain
                push!(chain, term)
                if haskey(chains, term) break end
                term = sum(factorials[parse.(Int, collect(string(term))) .+ 1])
            end
            index = findfirst(isequal(term), chain)
            loop_len = haskey(chains, term) ? chains[term] : length(chain) - index + 1
            for (i, j) ∈ enumerate(chain) chains[j] = i ≥ index ? loop_len : index - i + loop_len end
            if chains[number] == m
                digits = length(combination)
                result += (digits - count(isequal(0), combination)) * factorials[digits] ÷ prod(factorials[values(countmap(combination)) .+ 1])
            end
            if 0 ∈ combination
                deleteat!(combination, findfirst(isequal(0), combination))
                if length(combination) == 0 break end
            end
        end
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler