text:
const DEFAULT_LINESEARCH_nmax=100
const DEFAULT_LINESEARCH_rmax=100

const DEFAULT_WOLFE_ϵ  = 1E-4

abstract type LinesearchState end

LinesearchState(algorithm, f::Callable, x::Number; kwargs...) = LinesearchState(algorithm, UnivariateObjective(f, x); kwargs...)
LinesearchState(algorithm, f::Callable, x::AbstractVector; kwargs...) = LinesearchState(algorithm, MultivariateObjective(f, x); kwargs...)

solve!(x, δx, ls::LinesearchState) = ls(x, δx)
solve!(x, δx, g, ls::LinesearchState) = ls(x, δx, g)

struct Linesearch{ALG <: LinesearchMethod, OBJ <: AbstractObjective, OPT <: Options, OST <: LinesearchState}
    algorithm::ALG
    objective::OBJ
    config::OPT
    state::OST

    function Linesearch(algorithm, objective, config, state)
        new{typeof(algorithm), typeof(objective), typeof(config), typeof(state)}(algorithm, objective, config, state)
    end
end

function Linesearch(x, objective::AbstractObjective; algorithm = Static(), config = Options(), kwargs...)
    state = LinesearchState(algorithm, objective; kwargs...)
    Linesearch(algorithm, objective, config, state)
end

function Linesearch(x::Number, F::Callable; D = nothing, kwargs...)
    objective = UnivariateObjective(F, D, x)
    Linesearch(x, objective; kwargs...)
end

function Linesearch(x::AbstractVector, F::Callable; D = nothing, kwargs...)
    objective = MultivariateObjective(F, D, x)
    Linesearch(x, objective; kwargs...)
end

(ls::Linesearch)(args...) = ls.state(args...)

solve!(x, δx, ls::Linesearch) = solve!(x, δx, ls.state)


text:
using FirstOrderOptimization

srand(1)

## test asymmetric sketch
m = 100; n = 110; r = 5

X = zeros(m,n)
s = AsymmetricSketch(m,n,r)
for iter=1:r
	Xi = LowRankOperator(rand(m),rand(n)')
	X += Array(Xi)
	additive_update!(s,Xi)
end
@show vecnorm(Array(reconstruct(s)) - X) / vecnorm(X)

## test symmetric sketch
n = 100; r = 5

X = zeros(n,n)
s = SymmetricSketch(n,r)
for iter=1:r
	v = rand(n)
	Xi = LowRankOperator(v,v')
	X += v*v'
	additive_update!(s,Xi)
end
@show vecnorm(Array(reconstruct(s)) - X) / vecnorm(X)

ext: jl

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
print_all(Lexers.TOMLLexer, "toml")

tokentest(
    Lexers.TOMLLexer,
    "# Comment",
    COMMENT => "# Comment",
)

tokentest(
    Lexers.TOMLLexer,
    "x = 1",
    NAME => "x",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    NUMBER_INTEGER => "1",
)

tokentest(
    Lexers.TOMLLexer,
    "[header]",
    PUNCTUATION => "[",
    NAME => "header",
    PUNCTUATION => "]",
)

tokentest(
    Lexers.TOMLLexer,
    "name = \"name\"",
    NAME => "name",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    STRING => "\"name\"",
)

tokentest(
    Lexers.TOMLLexer,
    "dob = 1979-05-27T07:32:00-08:00",
    NAME => "dob",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    LITERAL_DATE => "1979-05-27T07:32:00-08:00",
)

tokentest(
    Lexers.TOMLLexer,
    "ports = [ 8001, 8001, 8002 ]",
    NAME => "ports",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    PUNCTUATION => "[",
    TEXT => " ",
    NUMBER_INTEGER => "8001",
    PUNCTUATION => ",",
    TEXT => " ",
    NUMBER_INTEGER => "8001",
    PUNCTUATION => ",",
    TEXT => " ",
    NUMBER_INTEGER => "8002",
    TEXT => " ",
    PUNCTUATION => "]",
)

tokentest(
    Lexers.TOMLLexer,
    "enabled = true",
    NAME => "enabled",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    KEYWORD_CONSTANT => "true",
)

tokentest(
    Lexers.TOMLLexer,
    "[a.b]",
    PUNCTUATION => "[",
    NAME => "a",
    PUNCTUATION => ".",
    NAME => "b",
    PUNCTUATION => "]",
)

tokentest(
    Lexers.TOMLLexer,
    "\"ʎǝʞ\" = \"value\"",
    STRING => "\"ʎǝʞ\"",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    STRING => "\"value\"",
)

tokentest(
    Lexers.TOMLLexer,
    "'quoted \"value\"' = \"value\"",
    STRING => "'quoted \"value\"'",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    STRING => "\"value\"",
)

tokentest(
    Lexers.TOMLLexer,
    "\" \\n \\t \"",
    STRING => "\" ",
    STRING_ESCAPE => "\\n",
    STRING => " ",
    STRING_ESCAPE => "\\t",
    STRING => " \"",
)

tokentest(
    Lexers.TOMLLexer,
    "'<\\i\\c*\\s*>'",
    STRING => "'<\\i\\c*\\s*>'",
)

tokentest(
    Lexers.TOMLLexer,
    "+99 42 0 -17 1_000 5_349_221 1_2_3_4_5 +1.0 3.1415 -0.01",
    NUMBER_INTEGER => "+99",
    TEXT => " ",
    NUMBER_INTEGER => "42",
    TEXT => " ",
    NUMBER_INTEGER => "0",
    TEXT => " ",
    NUMBER_INTEGER => "-17",
    TEXT => " ",
    NUMBER_INTEGER => "1_000",
    TEXT => " ",
    NUMBER_INTEGER => "5_349_221",
    TEXT => " ",
    NUMBER_INTEGER => "1_2_3_4_5",
    TEXT => " ",
    NUMBER_FLOAT => "+1.0",
    TEXT => " ",
    NUMBER_FLOAT => "3.1415",
    TEXT => " ",
    NUMBER_FLOAT => "-0.01",
)

tokentest(
    Lexers.TOMLLexer,
    "5e+22 1e6 -2E-2 6.626e-34 9_224_617.445_991_228_313",
    NUMBER_FLOAT => "5e+22",
    TEXT => " ",
    NUMBER_FLOAT => "1e6",
    TEXT => " ",
    NUMBER_FLOAT => "-2E-2",
    TEXT => " ",
    NUMBER_FLOAT => "6.626e-34",
    TEXT => " ",
    NUMBER_FLOAT => "9_224_617.445_991_228_313",
)

tokentest(
    Lexers.TOMLLexer,
    "date1 = 1979-05-27T07:32:00Z",
    NAME => "date1",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    LITERAL_DATE => "1979-05-27T07:32:00Z",
)

tokentest(
    Lexers.TOMLLexer,
    "date2 = 1979-05-27T00:32:00-07:00",
    NAME => "date2",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    LITERAL_DATE => "1979-05-27T00:32:00-07:00",
)

tokentest(
    Lexers.TOMLLexer,
    "date3 = 1979-05-27T00:32:00.999999-07:00",
    NAME => "date3",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    LITERAL_DATE => "1979-05-27T00:32:00.999999-07:00",
)

tokentest(
    Lexers.TOMLLexer,
    "1979-05-27T07:32:00",
    LITERAL_DATE => "1979-05-27T07:32:00",
)

tokentest(
    Lexers.TOMLLexer,
    "1979-05-27T00:32:00.999999",
    LITERAL_DATE => "1979-05-27T00:32:00.999999",
)

tokentest(
    Lexers.TOMLLexer,
    "1979-05-27",
    LITERAL_DATE => "1979-05-27",
)

tokentest(
    Lexers.TOMLLexer,
    "point = { x = 1, y = 2 }",
    NAME => "point",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    PUNCTUATION => "{",
    TEXT => " ",
    NAME => "x",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    NUMBER_INTEGER => "1",
    PUNCTUATION => ",",
    TEXT => " ",
    NAME => "y",
    TEXT => " ",
    OPERATOR => "=",
    TEXT => " ",
    NUMBER_INTEGER => "2",
    TEXT => " ",
    PUNCTUATION => "}",
)


text:
module Timer

using ..LIKWID: LibLikwid, timer_initialized

"Initialize LIKWIDs timer module"
function init()
    ret = LibLikwid.timer_init()
    timer_initialized[] = true
    return true
end

"Close and finalize LIKWIDs timer module"
function finalize()
    LibLikwid.timer_finalize()
    timer_initialized[] = false
    return nothing
end

"""
Return the CPU clock determined at [`Timer.init()`](@ref).
"""
function get_cpu_clock()
    if !timer_initialized[]
        init() || error("Couldn't init timer.")
    end
    return Int(LibLikwid.timer_getCpuClock())
end

"""
Return the current CPU clock read from sysfs
"""
function get_cpu_clock_current(cpu_id::Integer)
    if !timer_initialized[]
        init() || error("Couldn't init timer.")
    end
    return Int(LibLikwid.timer_getCpuClockCurrent(cpu_id))
end

"Start the clock and return a `LibLikwid.TimerData` object including the start timestamp."
function start_clock()
    if !timer_initialized[]
        init() || error("Couldn't init timer.")
    end
    timer = LibLikwid.TimerData()
    timer_ref = Ref(timer)
    LibLikwid.timer_start(timer_ref)
    return timer_ref[]
end

"""
    stop_clock(timer::LibLikwid.TimerData) -> newtimer::LibLikwid.TimerData
Stop the clock and return a `LibLikwid.TimerData` object including the start and stop timestamps.
The input `timer` should be the output of [`Timer.start_clock()`](@ref).
"""
function stop_clock(timer::LibLikwid.TimerData)
    if !timer_initialized[]
        init() || error("Couldn't init timer.")
    end
    timer_ref = Ref(timer)
    LibLikwid.timer_stop(timer_ref)
    return timer_ref[]
end

"""
    get_clock(timer::LibLikwid.TimerData)
Return the measured interval in seconds for the given `timer`.
The input `timer` should be the output of [`Timer.stop_clock`](@ref).
"""
function get_clock(timer::LibLikwid.TimerData)
    if !timer_initialized[]
        init() || error("Couldn't init timer.")
    end
    return LibLikwid.timer_print(Ref(timer))
end

"""
    get_clock_cycles(timer::LibLikwid.TimerData)
Return the measured interval in cycles for the given `timer`.
The input `timer` should be the output of [`Timer.stop_clock`](@ref).
"""
function get_clock_cycles(timer::LibLikwid.TimerData)
    if !timer_initialized[]
        init() || error("Couldn't init timer.")
    end
    if iszero(timer.start) || iszero(timer.stop)
        error("Start or stop is zero.")
    end
    return Int(LibLikwid.timer_printCycles(Ref(timer)))
end

"""
    timeit(f)
Time the given function `f` using [`Timer.start_clock`](@ref),
[`Timer.stop_clock`](@ref), etc. under the hood. Automatically
initializes and finalizes the timer module.

# Examples
```julia
julia> LIKWID.Timer.timeit() do
           sleep(1)
       end
(clock = 1.0008815780376372, cycles = 3603224844)
```
"""
function timeit(f)
    init() || error("Couldn't init LIKWIDs timer module.")
    try
        t_start = start_clock()
        f()
        t_stop = stop_clock(t_start)

        return (clock=get_clock(t_stop), cycles=get_clock_cycles(t_stop))
    finally
        finalize()
    end
end

"""
Convenience macro for [`Timer.timeit`](@ref).

# Examples
```julia
julia> LIKWID.Timer.@timeit sleep(1)
(clock = 1.0008815780376372, cycles = 3603224844)
```
"""
macro timeit(expr)
    q = quote
        LIKWID.Timer.timeit() do
            $(expr)
        end
    end
    return esc(q)
end

end # module


text:
using ApproxFunBase, ApproxFunOrthogonalPolynomials, LinearAlgebra, Test
    import ApproxFunBase: testbandedblockbandedoperator, testblockbandedoperator, testraggedbelowoperator, Block, ldiv_coefficients

@testset "PDE" begin
    @testset "Rectangle Laplace/Poisson" begin
        dx = dy = ChebyshevInterval()
        d = dx × dy
        g = Fun((x,y)->exp(x)*cos(y),∂(d))

        B = Dirichlet(d)

        testblockbandedoperator(B)
        testbandedblockbandedoperator(Laplacian(d))
        testbandedblockbandedoperator(Laplacian(d)+0.0I)

        A=[Dirichlet(d);Laplacian(d)]

        @time u=A\[g,0.]
        @test u(.1,.2) ≈ real(exp(0.1+0.2im))

        A=[Dirichlet(d);Laplacian(d)+0.0I]
        @time u=A\[g,0.]

        @test u(.1,.2) ≈ real(exp(0.1+0.2im))

        ## Poisson
        f=Fun((x,y)->exp(-10(x+.2)^2-20(y-.1)^2),ChebyshevInterval()^2,500)  #default is [-1,1]^2
        d=domain(f)
        A=[Dirichlet(d);Laplacian(d)]
        @time  u=\(A,[zeros(∂(d));f];tolerance=1E-7)
        @test ≈(u(.1,.2),-0.04251891975068446;atol=1E-5)
    end

    @testset "Bilaplacian" begin
        dx = dy = ChebyshevInterval()
        d = dx × dy
        Dx = Derivative(dx); Dy = Derivative(dy)
        L = Dx^4⊗I + 2*Dx^2⊗Dy^2 + I⊗Dy^4

        testbandedblockbandedoperator(L)

        B = Dirichlet(dx) ⊗ Operator(I,dy)
        testraggedbelowoperator(B)

        A=[Dirichlet(dx) ⊗ Operator(I,dy);
                Operator(I,dx)  ⊗ Dirichlet(dy);
                Neumann(dx) ⊗ Operator(I,dy);
                Operator(I,dx) ⊗ Neumann(dy);
                 L]

        testraggedbelowoperator(A)

        @time u=\(A,[[1,1],[1,1],[0,0],[0,0],0];tolerance=1E-5)
        @test u(0.1,0.2) ≈ 1.0
    end

    @testset "Schrodinger" begin
        dx=0..1; dt=0.0..0.001
        C=Conversion(Chebyshev(dx)*Ultraspherical(1,dt),Ultraspherical(2,dx)*Ultraspherical(1,dt))
        testbandedblockbandedoperator(C)
        testbandedblockbandedoperator(Operator{ComplexF64}(C))

        d = dx × dt

        x,y = Fun(d)
        @test x(0.1,0.0001) ≈ 0.1
        @test y(0.1,0.0001) ≈ 0.0001

        V = x^2

        Dt=Derivative(d,[0,1]);Dx=Derivative(d,[1,0])

        ϵ = 1.
        u0 = Fun(x->exp(-100*(x-.5)^2)*exp(-1/(5*ϵ)*log(2cosh(5*(x-.5)))),dx)
        L = ϵ*Dt+(.5im*ϵ^2*Dx^2)
        testbandedblockbandedoperator(L)

        @time u = \([timedirichlet(d);L],[u0,[0.,0.],0.];tolerance=1E-5)
        @test u(0.5,0.001) ≈ 0.857215539785593+0.08694948835021317im  # empircal from ≈ schurfact
    end

    @testset "Neumann" begin
        d=ChebyshevInterval()^2

        @time u=\([Neumann(d); Laplacian(d)-100.0I],[[[1,1],[1,1]],0.];tolerance=1E-12)
        @test u(.1,.9) ≈ 0.03679861429138079
    end

    @testset "Transport" begin
        dx=ChebyshevInterval(); dt=Interval(0,2.)
        d=dx × dt
        Dx=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
        x,y=Fun(identity,d)
        @time u=\([I⊗ldirichlet(dt);Dt+x*Dx],[Fun(x->exp(-20x^2),dx);0.];tolerance=1E-12)

        @test u(0.1,0.2) ≈ 0.8745340845783758  # empirical
    end

    @testset "Bilaplacian" begin
        dx=dy=ChebyshevInterval()
        d=dx × dy
        Dx=Derivative(dx);Dy=Derivative(dy)
        L=Dx^4⊗I+2*Dx^2⊗Dy^2+I⊗Dy^4

        testbandedblockbandedoperator(L)

        A=[ldirichlet(dx)⊗Operator(I,dy);
                rdirichlet(dx)⊗Operator(I,dy);
                Operator(I,dx)⊗ldirichlet(dy);
                Operator(I,dx)⊗rdirichlet(dy);
                lneumann(dx)⊗Operator(I,dy);
                rneumann(dx)⊗Operator(I,dy);
                Operator(I,dx)⊗lneumann(dy);
                Operator(I,dx)⊗rneumann(dy);
                L]

        # Checks bug in constructor
        f=Fun((x,y)->real(exp(x+1.0im*y)),component(rangespace(A)[1],1),22)
        @test f(-1.,0.1) ≈ real(exp(-1+0.1im))
        f=Fun((x,y)->real(exp(x+1.0im*y)),component(rangespace(A)[1],1))
        @test f(-1.,0.1) ≈ real(exp(-1+0.1im))

        F=[Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[1]);
            Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[2]);
            Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[3]);
            Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[4]);
            Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[5]);
            Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[6]);
            Fun((x,y)->-imag(exp(x+1.0im*y)),rangespace(A)[7]);
            Fun((x,y)->-imag(exp(x+1.0im*y)),rangespace(A)[8]);
            0]

        @time u=\(A,F;tolerance=1E-10)

        @test u(0.1,0.2)  ≈ exp(0.1)*cos(0.2)

        dx=dy=ChebyshevInterval()
        d=dx×dy
        Dx=Derivative(dx);Dy=Derivative(dy)
        L=Dx^4⊗I+2*Dx^2⊗Dy^2+I⊗Dy^4

        A=[(ldirichlet(dx)+lneumann(dx))⊗Operator(I,dy);
                (rdirichlet(dx)+rneumann(dx))⊗Operator(I,dy);
                Operator(I,dx)⊗(ldirichlet(dy)+lneumann(dy));
                Operator(I,dx)⊗(rdirichlet(dy)+rneumann(dy));
                (ldirichlet(dx)-lneumann(dx))⊗Operator(I,dy);
                (rdirichlet(dx)-rneumann(dx))⊗Operator(I,dy);
                Operator(I,dx)⊗(ldirichlet(dy)-lneumann(dy));
                Operator(I,dx)⊗(rdirichlet(dy)-rneumann(dy));
                L]

        u=\(A,[fill(1.0,8);0];tolerance=1E-5)
        @test u(0.1,0.2) ≈ 1.0

        F=[2Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[1]);
            2Fun((x,y)->real(exp(x+1.0im*y)),rangespace(A)[2]);
            Fun((x,y)->real(exp(x+1.0im*y))-imag(exp(x+1.0im*y)),rangespace(A)[3]);
            Fun((x,y)->real(exp(x+1.0im*y))-imag(exp(x+1.0im*y)),rangespace(A)[4]);
            0;
            0;
            Fun((x,y)->real(exp(x+1.0im*y))+imag(exp(x+1.0im*y)),rangespace(A)[7]);
            Fun((x,y)->real(exp(x+1.0im*y))+imag(exp(x+1.0im*y)),rangespace(A)[8]);
            0]

        u=\(A,F;tolerance=1E-10)

        @test u(0.1,0.2)  ≈ exp(0.1)*cos(0.2)
    end

    @testset "Check resizing" begin
        d=ChebyshevInterval()^2
        A=[Dirichlet(d);Laplacian()+100I]
        QRR = qr(A)
        @time ApproxFunBase.resizedata!(QRR.R_cache,:,2000)
        @test norm(QRR.R_cache.data[1:200,1:200] - A[1:200,1:200]) ==0

        @time ApproxFunBase.resizedata!(QRR,:,200)
        j=56
        v=QRR.R_cache.op[1:100,j]
        @test norm(ldiv_coefficients(QRR.Q,v;maxlength=300)[j+1:end]) < 100eps()

        j=195
        v=QRR.R_cache.op[1:ApproxFunBase.colstop(QRR.R_cache.op,j),j]
        @test norm(ldiv_coefficients(QRR.Q,v;maxlength=1000)[j+1:end]) < 100eps()

        j=300
        v=QRR.R_cache.op[1:ApproxFunBase.colstop(QRR.R_cache.op,j),j]
        @test norm(ldiv_coefficients(QRR.Q,v;maxlength=1000)[j+1:end]) < j*20eps()

        @test ApproxFunBase.colstop(QRR.R_cache.op,195)-194 == ApproxFunBase.colstop(QRR.H,195)

        QR1 = qr(A)
        @time ApproxFunBase.resizedata!(QR1.R_cache,:,1000)
        QR2 = qr([Dirichlet(d);Laplacian()+100I])
        @time ApproxFunBase.resizedata!(QR2.R_cache,:,500)
        n=450;QR1.R_cache.data[1:n,1:n]-QR2.R_cache.data[1:n,1:n]|>norm
        @time ApproxFunBase.resizedata!(QR2.R_cache,:,1000)
        N=450;QR1.R_cache.data[1:N,1:N]-QR2.R_cache.data[1:N,1:N]|>norm
        N=1000;QR1.R_cache.data[1:N,1:N]-QR2.R_cache.data[1:N,1:N]|>norm

        QR1 = qr(A)
        @time ApproxFunBase.resizedata!(QR1,:,1000)
        QR2 = qr([Dirichlet(d);Laplacian()+100I])
        @time ApproxFunBase.resizedata!(QR2,:,500)
        @time ApproxFunBase.resizedata!(QR2,:,1000)

        @test norm(QR1.H[1:225,1:1000]-QR2.H[1:225,1:1000]) ≤ 10eps()

        QR1 = qr(A)
        @time ApproxFunBase.resizedata!(QR1,:,5000)
        @time u=\(QR1,[ones(∂(d));0.];tolerance=1E-7)

        @test norm((Dirichlet(d)*u-ones(∂(d))).coefficients) < 1E-7
        @test norm((A*u-Fun([ones(∂(d));0.])).coefficients) < 1E-7
        @test norm(((A*u)[2]-(Laplacian(space(u))+100I)*u).coefficients) < 1E-10
        @test eltype(ApproxFunBase.promotedomainspace(Laplacian(),space(u))) == Float64
        @test eltype(ApproxFunBase.promotedomainspace(Laplacian()+100I,space(u))) == Float64
        @test norm(((A*u)[2]-(Laplacian()+100I)*u).coefficients) < 1E-10
        @test norm((Laplacian()*u+100*u - (A*u)[2]).coefficients) < 10E-10
        @time v=\(A,[ones(∂(d));0.];tolerance=1E-7)
        @test norm((u-v).coefficients) < 100eps()

        @test u(0.1,1.) ≈ 1.0
        @test u(0.1,-1.) ≈ 1.0
        @test u(1.,0.1) ≈ 1.0
        @test u(-1.,0.1) ≈ 1.0
    end

    @testset "Chebyshev Dirichlet Kronecker" begin
        S=ChebyshevDirichlet()^2
        ff=(x,y)->exp(x)*cos(y)
        u=Fun(ff,S)

        for KO in [Operator(I,factor(S,1))⊗rdirichlet(factor(S,1)),
                    rdirichlet(factor(S,1))⊗Operator(I,factor(S,2))]
            testblockbandedoperator(KO)
            @test norm((KO*u-Fun(ff,rangespace(KO))).coefficients) ≤ 1E-10
        end
    end

    @testset "Small diffusion" begin
        dx=ChebyshevInterval();dt=Interval(0,0.2)
        d=dx×dt
        Dx=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
        x,t=Fun(dx×dt)

        B=0.0
        C=0.0
        V=B+C*x
        ε=0.1
        f=Fun(x->exp(-30x^2),dx)
        u=\([timedirichlet(d);Dt-ε*Dx^2-V*Dx],[f;zeros(3)];tolerance=1E-6)

        @test u(.1,.2) ≈ 0.496524222625512
        B=0.1
        C=0.2
        V=B+C*x
        u=\([timedirichlet(d);Dt-ε*Dx^2-V*Dx],[f;zeros(3)];tolerance=1E-7)
        @test u(.1,.2) ≈ 0.46810331039791464
    end

    @testset "concatenate InterlaceOperator" begin
        a=Fun(x -> 0 ≤ x ≤ 0.5 ? 0.5 : 1, Domain(-1..1) \ [0,0.5])
        @test a(0.1) == 0.5
        @test a(0.7) == 1.0
        s=space(a)
        # Bx=[ldirichlet(s);continuity(s,0)]
        # TODO: this should concat
        dt=Interval(0,2.)
        Dx=Derivative(s);Dt=Derivative(dt)
        Bx=[ldirichlet(s);continuity(s,0)]

        @test ApproxFunBase.rangetype(rangespace(continuity(s,0))) == Vector{Float64}
        @test ApproxFunBase.rangetype(rangespace(Bx)) == Vector{Any}
        @test ApproxFunBase.rangetype(rangespace(Bx⊗Operator(I,Chebyshev()))) == Vector{Any}

        rhs = Fun([0,[0,[0,0]],0],rangespace([I⊗ldirichlet(dt);Bx⊗I;I⊗Dt+(a*Dx)⊗I]))
        @test rhs(-0.5,0.0) == [0,[0,[0,0]],0]

        u=\([I⊗ldirichlet(dt);Bx⊗I;I⊗Dt+(a*Dx)⊗I],
            [Fun(x->exp(-20(x+0.5)^2),s),[[0],[0,0]],0.0];tolerance=1E-2)

        @test u(-0.4,0.1) ≈ u(-0.5,0.0) atol = 0.0001
    end
end


text:
using CMakeWrapper
using BinDeps
using Test

@test isfile(cmake_executable)
run(`$cmake_executable --version`)

ctx = BinDeps.PackageContext(true,
                             pwd(),
                             "CMakeWrapperTest",
                             [])
libdep = BinDeps._library_dependency(ctx, "libfoo")
libdep.helpers = [(BinDeps.NetworkSource(URI("src")), Dict{Symbol, Any}(:filename=>"foo.zip"))]
process = CMakeProcess()
steps = BinDeps.generate_steps(libdep, process, Dict())


text:
using RNGPool
using SMCExamples.LinearGaussian
import SMCExamples.LinearGaussian: LGTheta, simulateLGModel, kalman,
  makeLGModel
import SMCExamples.Particles.Float64Particle

function setupLGModel(n::Int64, obsVariance::Float64)
  theta = LGTheta(1.0, 1.0, 1.0, obsVariance, 0.0, 1.0)
  ys = simulateLGModel(theta, n)
  model = makeLGModel(theta, ys)
  lM = LinearGaussian.makelM(theta)
  ko = kalman(theta, ys)
  return model, lM, ko
end


text:
let
    mpsim = MPsim()

    attrib = Attribute("Attr A")
    addPossibleAttributeValue!( attrib, "Val AA", "Val AB" )
    addSimulationAttribute!( mpsim, attrib )

    node = BaseNode("Test A")
    addNodeRequirement!( node, ("Attr A", "Val AA") )
    addSimulationBaseNode!( mpsim, node )

    recr = Recruitment("EW")
    setRecruitmentTarget!( recr, "Test A" )
    setRecruitmentSchedule!( recr, 12 )
    setRecruitmentFixed!( recr, 10 )
    setRecruitmentAgeFixed!( recr, 18*12 )
    addSimulationRecruitment!( mpsim, recr )

    retir = Retirement()
    setRetirementCareerLength!( retir, 60 )
    setRetirementSchedule!( retir, 12 )
    setSimulationRetirement!( mpsim, retir )

    setSimulationLength!( mpsim, 120 )
    run( mpsim )
end


text:
function simulate(m::Translation, noise::Distribution;
    duration::Number = 6, Δt::Number = 0.01)

    start, stop = m.start, m.stop
    stop <= duration || error("Movement extends beyond simulation event.")
    m.D ∈ ["forwards", "backwards"] || error("Direction not specified.")

    D = m.D == "forwards" ? 1 : -1
    f = 1/(stop - start) # frequency: single sinusiodal head turn

    timesteps = range(0, stop=duration, step=Δt)
    T = length(timesteps)

    # κ1, κ2 = cupula_dynamics(4, Δt)
    α = Array{Float64}(undef, T) # acceleration
    ω = Array{Float64}(undef, T) # velocity
    θ = Array{Float64}(undef, T) # position
    # c = Array{Float64}(undef, T) # cupula dynamics
    y = Array{Float64}(undef, T) # sensory observations

    # α[1]
    ω[1] = 0
    # c[1] = 0
    θ[1] = 0
    y[1] = α[1] + rand(noise)

    for (i, t) ∈ enumerate(timesteps)
        α[i] = t > start && t < stop ? acceleration(D, m.A, f, t, start) : 0

        if i > 1
            ω[i] = ω[i-1] + Δt * α[i-1]
            # c[i] = κ1 * c[i-1] + κ2 * ω[i]
            θ[i] = θ[i-1] + Δt * ω[i-1] + 0.5*Δt^2 * α[i-1]
        end
        y[i] = α[i] + rand(noise)
    end
    return collect(timesteps), y, α, ω, θ, Δt
end

function simulate(m1::HeadTurn, m2::ExternalDisturbance,
    noise::Distribution; duration::Number = 6,
    Δt::Number = 0.01)

    # start, stop = m1.times
    m1.stop <= duration || error("Movement extends beyond simulation event.")
    m2.stop < duration || error("Movement extends beyond simulation event.")
    m1.D ∈ ["left", "right"] || error("Head turn direction not specified.")
    m2.D ∈ ["left", "right"] || error("Disturbance direction not specified.")

    Dh = m1.D == "left" ? 1 : -1
    Dd = m2.D == "left" ? 1 : -1

    f = 1/(m1.stop - m1.start) # frequency: single sinusiodal head turn

    timesteps = range(0, stop=duration, step=Δt)
    T = length(timesteps)

    κ1, κ2 = cupula_dynamics(4, Δt)

    ωᵤ = Array{Float64}(undef, T) # active (controlled)
    ωₑ = Array{Float64}(undef, T) # passive (external)
    ω = Array{Float64}(undef, T) # total velocity
    θ = Array{Float64}(undef, T) # position
    c = Array{Float64}(undef, T) # cupula dynamics
    y = Array{Float64}(undef, T) # sensory observations

    ωᵤ[1] = 0
    ωₑ[1] = 0
    ω[1] = ωᵤ[1] + ωₑ[1]
    θ[1] = 0
    y[1] = ω[1] - c[1] + rand(noise)

    for (i, t) ∈ enumerate(timesteps)
        ωᵤ[i] = t > m1.start && t < m1.stop ? acceleration(Dh, m1.A, f, t, m1.start) : 0
        ωₑ[i] = t > m2.start && t < m2.stop ? Dd * m2.A : 0
        ω[i] = ωᵤ[i] + ωₑ[i]
        if i > 1
            c[i] = κ1 * c[i-1] + κ2 * ω[i]
            θ[i] = θ[i-1] + Δt * ω[i-1]
        end
        y[i] = ω[i] - c[i] + rand(noise)
    end

    p = scatter(timesteps, y,
        # marker = (:circle, 1, 0.5, :black,
        #     stroke(0.5, 0.5, :black)))
            markershape = :circle,
            markersize = 2,
            markeralpha = 0.4,
            markercolor = :black)
    # plot!(p, time, c,
    #     linewidth = 4,
    #     linealpha = 0.5,
    #     linecolor = :steelblue)
    plot!(p, timesteps, ωₑ,
        linewidth = 6,
        linecolor = :orange,
        xlabel = "Time (s)",
        ylabel = "Angular displacement",
        legend = :none)
    plot!(p, timesteps, ωᵤ,
        linewidth = 6,
        linestyle = :dot,
        linecolor = :steelblue1)
    plot!(p, timesteps, ω,
        linewidth = 6,
        linecolor = :steelblue4)

    plot!(p, timesteps, θ,
        linewidth = 6,
        linestyle = :dot,
        linecolor = :purple)

    return collect(timesteps), y, ω, ωᵤ, ωₑ, θ, c, Δt
end


text:
@testset "Events/Element Starts" begin
    evaluate(s) = collect(E.events(L.State(IOBuffer(s))))

    E = Events
    L = Events.Lexical

    ElementEnd   = E.ElementEnd
    ElementStart = E.ElementStart

    DC = E.DataContent
    ME = E.MarkupError

    @testset "Events/Element Starts (Positive)" begin
        @test evaluate("")   == [ ElementStart("a", [ ], L.Location("a buffer", -1)) ]
        @test evaluate("<abc>") == [ ElementStart("abc", [ ], L.Location("a buffer", -1)) ]
        @test evaluate("<é>")   == [ ElementStart("é", [ ], L.Location("a buffer", -1)) ]
        @test evaluate("<aé>")  == [ ElementStart("aé", [ ], L.Location("a buffer", -1)) ]
        @test evaluate("<éb>")  == [ ElementStart("éb", [ ], L.Location("a buffer", -1)) ]
        @test evaluate("<aéb>") == [ ElementStart("aéb", [ ], L.Location("a buffer", -1)) ]

        # White space is allowed after the element name. Just try various combinations, and check that it is discarded.
        #
        @test evaluate("")        == [ ElementStart("a", [ ], L.Location("a buffer", -1)) ]
        @test evaluate("")     == [ ElementStart("a", [ ], L.Location("a buffer", -1)) ]
        @test evaluate("") == [ ElementStart("a", [ ], L.Location("a buffer", -1)) ]

        @test evaluate("") == [ ElementStart("a", [ ], L.Location("a buffer", -1)),
                                      ElementStart("b", [ ], L.Location("a buffer", -1)) ]
    end

    @testset "Events/Empty Elements (Positive)" begin
        @test evaluate("")   == [ ElementStart(true, "a", [ ], L.Location("a buffer", -1)),
                                      ElementEnd("a", L.Location("a buffer", -1)) ]
        @test evaluate("<abc/>") == [ ElementStart(true, "abc", [ ], L.Location("a buffer", -1)),
                                      ElementEnd("abc", L.Location("a buffer", -1)) ]
        @test evaluate("<é/>")   == [ ElementStart(true, "é", [ ], L.Location("a buffer", -1)),
                                      ElementEnd("é", L.Location("a buffer", -1)) ]
        @test evaluate("<aé/>")  == [ ElementStart(true, "aé", [ ], L.Location("a buffer", -1)),
                                      ElementEnd("aé", L.Location("a buffer", -1)) ]
        @test evaluate("<éb/>")  == [ ElementStart(true, "éb", [ ], L.Location("a buffer", -1)),
                                      ElementEnd("éb", L.Location("a buffer", -1)) ]
        @test evaluate("<aéb/>") == [ ElementStart(true, "aéb", [ ], L.Location("a buffer", -1)),
                                      ElementEnd("aéb", L.Location("a buffer", -1)) ]

        # White space is allowed after the element name. Just try various combinations, and check that it is discarded.
        #
        @test evaluate("")        == [ ElementStart(true, "a", [ ], L.Location("a buffer", -1)),
                                            ElementEnd("a", L.Location("a buffer", -1)) ]
        @test evaluate("")     == [ ElementStart(true, "a", [ ], L.Location("a buffer", -1)),
                                            ElementEnd("a", L.Location("a buffer", -1)) ]
        @test evaluate("") == [ ElementStart(true, "a", [ ], L.Location("a buffer", -1)),
                                            ElementEnd("a", L.Location("a buffer", -1)) ]

        @test evaluate("") == [ ElementStart(true, "a", [ ], L.Location("a buffer", -1)),
                                        ElementEnd("a", L.Location("a buffer", -1)),
                                        ElementStart(true, "b", [ ], L.Location("a buffer", -1)),
                                        ElementEnd("b", L.Location("a buffer", -1)) ]
    end

    @testset "Events/Element Starts (Negative ... no terminator)" begin
        @test (evaluate("' to end an element open tag.", L.Location("a buffer", -1)) ])

        @test (evaluate("<abc")
               == [ ElementStart(true, false, "abc", [ ], L.Location("a buffer", -1)),
                    ME("ERROR: Expecting '>' to end an element open tag.", L.Location("a buffer", -1)) ])

        @test (evaluate("' to end an element open tag.", L.Location("a buffer", -1)) ])

        events = evaluate("<a<?")
        @test length(events) == 3
        @test events == [ ElementStart(true, false, "a", [ ], L.Location("a buffer", -1)),
                          ME("ERROR: Expecting '>' to end an element open tag.", L.Location("a buffer", -1)),
                          ME("ERROR: Expecting a PI target.", L.Location("a buffer", -1)) ]

        events = evaluate("<a <?")
        @test length(events) == 3
        @test events == [ ElementStart(true, false, "a", [ ], L.Location("a buffer", -1)),
                          ME("ERROR: Expecting '>' to end an element open tag.", L.Location("a buffer", -1)),
                          ME("ERROR: Expecting a PI target.", L.Location("a buffer", -1)) ]
    end

    @testset "Events/Empty Elements (Negative ... no terminator)" begin
        @test (evaluate("' to end an element open tag.", L.Location("a buffer", -1)),
                    ElementEnd(true, "a", L.Location("a buffer", -1)) ])

        @test (evaluate("<abc/")
               == [ ElementStart(true, true, "abc", [ ], L.Location("a buffer", -1)),
                    ME("ERROR: Expecting '>' to end an element open tag.", L.Location("a buffer", -1)),
                    ElementEnd(true, "abc", L.Location("a buffer", -1)) ])

        @test (evaluate("' to end an element open tag.", L.Location("a buffer", -1)),
                    ElementEnd(true, "a", L.Location("a buffer", -1)) ])

        events = evaluate("<a/<?")
        @test length(events) == 4
        @test events == [ ElementStart(true, false, "a", [ ], L.Location("a buffer", -1)),
                          ME("ERROR: Expecting '>' to end an element open tag.", L.Location("a buffer", -1)),
                          ElementEnd(true, "a", L.Location("a buffer", -1)),
                          ME("ERROR: Expecting a PI target.", L.Location("a buffer", -1)) ]

        events = evaluate("<a /<?")
        @test length(events) == 4
        @test events == [ ElementStart(true, false, "a", [ ], L.Location("a buffer", -1)),
                          ME("ERROR: Expecting '>' to end an element open tag.", L.Location("a buffer", -1)),
                          ElementEnd(true, "a", L.Location("a buffer", -1)),
                          ME("ERROR: Expecting a PI target.", L.Location("a buffer", -1)) ]
    end

    @testset "Events/Element Starts (Negative ... no element name)" begin
        @test (evaluate("<")
               == [ ME("ERROR: Expecting an element name.", L.Location("a buffer", -1)) ])

        events = evaluate("<<?")
        @test length(events) == 2
        @test (events == [ ME("ERROR: Expecting an element name.", L.Location("a buffer", -1)),
                           ME("ERROR: Expecting a PI target.", L.Location("a buffer", -1)) ])

        events = evaluate("< ")
        @test length(events) == 2
        @test (events == [ ME("ERROR: Expecting an element name.", L.Location("a buffer", -1)),
                           DC(" ", true, L.Location("a buffer", -1)) ])
    end
end

max_stars_repo_name: jlegare/NativeXML

max_issues_repo_name: jlegare/NativeXML

max_forks_repo_name: jlegare/NativeXML


text:
@testset "Mergers" begin
    splitname = "splittestname"
    fbop = x -> 0.5 .* x.^2
    combineop = (x, y) -> x + σ.(y)
    m = Merger(splitname, fbop, combineop)
    x = randn(5)
    y = randn(5)
    h = Dict(splitname => y)

    @testset "constructor" begin
        @test m isa Merger{typeof(fbop), typeof(combineop)}
        @test m.splitname == splitname
        @test m.fb(x) ≈ fbop(x)
        @test m.op(x, y) ≈ combineop(x, y)
    end # @testset "constructor"

    @testset "apply" begin
        @test m(x, h) ≈ m.op(x, m.fb(y))
    end # @testset "apply"

    @testset "params" begin
        # Test that a Merger returns the parameters of its internal operations
        fb = Dense(10, 5)
        op = Chain((x, y) -> hcat(x, y), Dense(10, 5))
        m = Merger(splitname, fb, op)
        @test hcat(params(m)...) == hcat(params(fb)..., params(op)...)
    end # @testset "params"

    @testset "name" begin
        @test inputname(m) == splitname
    end # @testset "name"
end # @testset "Mergers"


text:
#!/usr/bin/env julia

"""
feature_counts.jl - Houses functions that deal with GFF feature-BAM alignment overlaps and counts

By: Shaun Adkins (sadkins@som.umaryland.edu)
"""

### Feature Overlap - Miscellaneous information about a feature with respect to its overlaps to alignments

mutable struct FeatureOverlap
    num_alignments::Float32
    feat_counts::Float32
    coords_set::BitSet
end

coordinate_set(f::FeatureOverlap) = f.coords_set
featurecounts(f::FeatureOverlap) = f.feat_counts
totalalignments(f::FeatureOverlap) = f.num_alignments

inc_featurecounts!(f::FeatureOverlap, counts) = f.feat_counts += counts
inc_alignments!(f::FeatureOverlap, num_alignments) = f.num_alignments += num_alignments

### MultimappedAlignment - Alignments that overlap features but are multimapped

struct MultimappedAlignment
    featurename::String
    align_feat_ratio::Float32
    alignmenttype::T where {T<:AbstractAlignment}
end

featurename(m::MultimappedAlignment) = m.featurename
alignment_feature_ratio(m::MultimappedAlignment) = m.align_feat_ratio
alignment_type(m::MultimappedAlignment) = m.alignmenttype

function calc_relative_abundance(feature_counts::Float32, template_totalcounts::Float32)::Float32
    """Calculate relative abundance fraction for the given feature with respect to all features overlapping the template."""
    @assert template_totalcounts > 0 "The template_totalcounts value should be greater than 0"
    # Total relative abundances for all features overlapping this template should = 1.0
    return @fastmath feature_counts / template_totalcounts
end

function calc_subset_totalcounts(feat_overlaps::Dict{String, FeatureOverlap}, feature_counter::Dict{String, UInt})
    """For all features overlapping the alignment template, return the sum of all those feature counts."""
    totalcounts = zero(Float32)
    # Worth noting if a feature overlaps multiple times, it will be reflected in totalcounts
    # This will be appropriately handled when all alignments are iterated through to increment the FeatureOverlap feat_counts for the feature.
    for (featurename, count) in feature_counter
        totalcounts += featurecounts(feat_overlaps[featurename]) * count
    end
    return totalcounts
end

function calc_totalcounts(feat_overlaps::Dict{String, FeatureOverlap})
    """Calculate the sum of all the feature counts."""
    return sum(featurecounts(feat_overlaps[featurename]) for featurename in keys(feat_overlaps))
end

function calc_tpm(len::UInt, totalcounts::Float32, feat_counts::Float32)
    """Calculate TPM score for current feature."""
    return @fastmath(feat_counts * 1000 / len) * 1000000 / totalcounts
end

function compute_alignment_feature_ratio(uniq_feat_coords::BitSet, alignment::Interval)::Float32
    """Calculate the ratio of fragament coordinates that intersect with non-overlapping feature coordinates."""
    # Pertinent alignment info
    alignment_coords =  get_alignment_coords_set(alignment)
    alignment_intersect = intersect(uniq_feat_coords, alignment_coords)
    # Percentage of alignment that aligned with this annotation feature
    # Return as Float32 for performance
    return @fastmath length(alignment_intersect) / length(alignment_coords)
end

function compute_mm_adjusted_counts(total_feat_counts::Float32, template_feat_counts::Float32, template_totalcounts::Float32)::Float32
    """Calculate the adjustment feature counts for a multimapped template to a given feature."""
    relative_abundance = calc_relative_abundance(total_feat_counts, template_totalcounts)
    return template_feat_counts * relative_abundance
end

function compute_mm_counts_by_em(feat_overlaps::Dict{String, FeatureOverlap}, multimapped_dict::Dict{String,StructArray}, args::Dict)
    """Adjust the feature counts of the multimapped overlaps via the Expectation-Maximization algorithm."""
    # Initialize multimapped overlap feature dictionary
    featurenames = collect(keys(feat_overlaps))
    adjusted_mm_overlaps = Dict{String, FeatureOverlap}(featurename => initialize_overlap_info(coordinate_set(feat_overlaps[featurename])) for featurename in featurenames)

    # Iterate through the multimaped alignments for each template
    @inbounds for (tn, templatealignments) in multimapped_dict
        # All multimapped alignments for a record should have the same alignment type
        alignmenttype = alignment_type(first(templatealignments))

        # I do not use "unique" for this because there may be a chance that the template aligns to the same feature multiple times
        overlappingfeatures = [featurename(a) for a in templatealignments]

        template_feature_counter = Dict{String, UInt}(fn => zero(UInt) for fn in unique(overlappingfeatures))
        template_feat_counts = Dict{String, Float32}(fn => zero(Float32) for fn in unique(overlappingfeatures))
        Threads.@threads for fn in unique(overlappingfeatures)
            # Filter all alignments that have this particular template-feature overlap
            templatefeaturealignments = filter(x -> x.featurename == fn, templatealignments)
            template_feature_counter[fn] += length(templatefeaturealignments)
            template_feat_counts[fn] += sum(alignment_feature_ratio, templatefeaturealignments)
        end

        template_totalcounts = calc_subset_totalcounts(feat_overlaps, template_feature_counter)
        # In cases where every feature for this template had no overlaps with any singly-mapped alignments, do not consider them.
        template_totalcounts > zero(Float32) || continue

        # Adjust contribution proportion of multimapped reads, based on estimated relative abundance
        for (featurename, counts) in template_feat_counts
            counts > zero(Float32) || continue
            adjusted_feat_counts = compute_mm_adjusted_counts(featurecounts(feat_overlaps[featurename]), counts, template_totalcounts)
            # If the feature had no contribution to any singly-mapped alignments, ignore it.
            adjusted_feat_counts > zero(Float32) && increment_feature_overlap_information!(adjusted_mm_overlaps[featurename], adjusted_feat_counts, alignmenttype)
        end
    end
    return adjusted_mm_overlaps
end

function create_feat_overlaps_dict(features::Array{GFF3.Record,1}, uniq_coords::Dict{String,Dict}, attr_type::String, stranded_type::String)
    feat_overlaps = Dict{String, FeatureOverlap}()
    for feature in features
        featurename = get_featurename_from_attrs(feature, attr_type)
        uniq_feat_coords = get_feature_nonoverlapping_coords(feature, uniq_coords, isstranded(stranded_type))
        feat_overlaps[featurename] = initialize_overlap_info(uniq_feat_coords)
    end
    return feat_overlaps
end

function filter_features_overlapping_alignments(features::Array{GFF3.Record,1}, alignment::Interval, isstranded::Bool)
    """Filter features to those just that align with the current alignment on the same strand."""
    alignmentstrand = getstrand(alignment, isstranded)
    overlapping_features = filter(x -> isoverlapping(convert(Interval, x), alignment), features)
    filter!(x -> alignmentstrand == getstrand(x, isstranded), overlapping_features)
    return overlapping_features   # Returns an array of GFF3 Records
end

function increment_feature_overlap_information!(feat_overlap::FeatureOverlap, align_feat_ratio::Float32, alignmenttype::T) where {T<:AbstractAlignment}
    """Increment number of alignments and feature count information for feature if alignment overlapped with uniq coords."""
    inc_alignments!(feat_overlap, multiplier(alignmenttype))
    inc_featurecounts!(feat_overlap, align_feat_ratio * multiplier(alignmenttype))
end

function initialize_overlap_info(uniq_feat_coords::BitSet)
    """Initialize overlap diction information."""
    num_alignments = zero(Float32); feat_counts = zero(Float32)
    return FeatureOverlap(num_alignments, feat_counts, uniq_feat_coords)
end

function merge_mm_counts(feat_overlaps::Dict{String, FeatureOverlap}, mm_feat_overlaps::Dict{String, FeatureOverlap}, adjust_num_alignments::Bool=true)
    """Merge EM-computed counts for multimapped reads and singly-mapped feature counts.  Return a new Dict of FeatureOverlaps"""
    adjusted_overlaps = deepcopy(feat_overlaps)
    featurenames = collect(keys(feat_overlaps))
    [inc_featurecounts!(adjusted_overlaps[fn], featurecounts(mm_feat_overlaps[fn])) for fn in featurenames]
    adjust_num_alignments || return adjusted_overlaps

    [inc_alignments!(adjusted_overlaps[fn], totalalignments(mm_feat_overlaps[fn])) for fn in featurenames]
    return adjusted_overlaps
end

function merge_mm_counts!(feat_overlaps::Dict{String, FeatureOverlap}, mm_feat_overlaps::Dict{String, FeatureOverlap}, adjust_num_alignments::Bool=true)
    """Merge EM-computed counts for multimapped reads back into the general feature overlap counts."""
    featurenames = collect(keys(feat_overlaps))
    [inc_featurecounts!(feat_overlaps[fn], featurecounts(mm_feat_overlaps[fn])) for fn in featurenames]
    adjust_num_alignments || return

    [inc_alignments!(feat_overlaps[fn], totalalignments(mm_feat_overlaps[fn])) for fn in featurenames]
end

function process_feature_overlaps!(feat_overlaps::Dict{String, FeatureOverlap}, multimapped_dict::Dict{String,StructArray}, reader::BAM.Reader, feature::GFF3.Record, args)
    """Process all alignment intervals that overlap with the given feature interval."""
    featurename = get_featurename_from_attrs(feature, args["attribute_type"])
    feat_overlap = feat_overlaps[featurename]
    max_frag_size = convert(UInt, args["max_fragment_size"])
    featurestrand = getstrand(feature, isstranded(args["stranded"]))

    # Returning SuperBAMRecord objects instead of BAM.Record objects
    for rec in eachoverlap(reader, feature, args["stranded"], max_frag_size)
        alignmentinterval = interval(rec)
        @assert alignmentinterval !== nothing "All iterated records should be a validated read or fragment"
        strand(rec) == featurestrand || continue

        alignmenttype = alignment_type(rec)
        args["pp_only"] && isa(alignmenttype, ReadAlignment) && continue

        align_feat_ratio =  compute_alignment_feature_ratio(coordinate_set(feat_overlap), alignmentinterval)::Float32
        # Save multimapped records as Interval objects for later, if needed, otherwise increment information for this feature now
        if ismultimapped(record(rec))
            if !args["rm_multimap"]
                templatename = BAM.tempname(record(rec))
                get!(multimapped_dict, templatename, StructArray{MultimappedAlignment}(undef,0))
                push!(multimapped_dict[templatename], MultimappedAlignment(featurename, align_feat_ratio, alignmenttype))
            end
            continue
        end
        align_feat_ratio > 0.0 && increment_feature_overlap_information!(feat_overlap, align_feat_ratio, alignmenttype)
    end

    # NOTE: This function modifies feat_overlaps and multimapped_dict
end

function process_all_feature_overlaps(feat_overlaps::Dict{String,FeatureOverlap}, features::Vector{GFF3.Record}, reader::BAM.Reader, args)
    """Process all features from the annotation file for overlaps with BAM alignments."""
    multimapped_dict = Dict{String,StructArray}()
    [process_feature_overlaps!(feat_overlaps, multimapped_dict, reader, feature, args) for feature in features]
    return multimapped_dict
end

max_stars_repo_name: IGS/FADU

max_issues_repo_name: adkinsrs/FADU

max_forks_repo_name: adkinsrs/FADU


text:
"""
Test script for dirBin0Rec0 model: one parameter for total number of links and one for reciprocity.
Investigate the existence of regions of the θ η plane where the latters cannot be estimated for finite N
"""

using ScoreDrivenERGM
using ScoreDrivenERGM.Utilities

import ScoreDrivenERGM.StaticNets: ErgmDirBin0Rec0(), ergm_par_from_mean_vals,diadProbFromPars , samplSingMatCan, stats_from_mat, get_one_mple, estimate,ErgmDirBin0Rec0, exp_val_stats, sample_ergm, estimate

using PyPlot
pygui(true)
using RCall
using Statistics

ErgmRcall.clean_start_RCall()
ergmTermsString = "edges +  mutual"
model = ErgmDirBin0Rec0()

function sample_ergm_sequence_static_estimate_mle(model::ErgmDirBin0Rec0, N, parDgpSeq, nSample)

    T = size(parDgpSeq)[2]
    parMle = zeros(T, 2, nSample)
    obs = zeros(T, 3, nSample)

    for i =1:T

        A_S = sample_ergm(model, N, parDgpSeq[:,i], nSample )

        obs[i,:,:] = reduce(hcat,[StaticNets.stats_from_mat(model, A) for A in eachslice(A_S, dims=3)])
        parMle[i,:,:] = reduce(hcat,[StaticNets.estimate(model, A) for A in eachslice(A_S, dims=3)])
    end

    return parMle, obs
end

N=30
begin
gix, ax = subplots(2,1)
Nvals = [30, 50, 70, 90]
θ_0 = -5
for N in Nvals
    parDgpSeq = reduce(hcat, [ [θ_0, η ] for η = 2:0.05:7])
    parMle, obs = sample_ergm_sequence_static_estimate_mle(model, N, parDgpSeq, 150 )
    fracZeroRec = mean(obs[:,2,:].==0;dims=2)
    fracUb = mean(obs[:,2,:].==obs[:,1,:]/2;dims=2)
    xplot = parDgpSeq[2,:]
    ax[1].plot(xplot, fracZeroRec)
    ax[2].plot(xplot, fracUb)

end
    ax[1].set_title("Fraction of samples with zero reciprocal pairs θ_dgp = $θ_0 ")
    ax[1].set_xlabel( "η_dgp")
    ax[1].legend(Nvals)
    ax[1].grid()
    ax[2].set_title("Fraction of samples with max reciprocal pairs θ_dgp = $θ_0 ")
    ax[2].set_xlabel( "η_dgp")
    ax[2].legend(Nvals)
    ax[2].grid()
    tight_layout()
end

max_stars_repo_name: domenicodigangi/SDERGM

max_issues_repo_name: domenicodigangi/SDERGM

max_forks_repo_name: domenicodigangi/SDERGM


text:
# sectors.jl
# Explicit realization of sectors to be used for gradation of vector spaces.

# TODO: Add prettier printing for charges.

# Abelian Sectors:
#-----------------
immutable Parity <: Abelian
    charge::Bool
end
Parity(d::Integer)=Parity(iseven(d))

*(c1::Parity,c2::Parity)=Parity(c1.charge != c2.charge)
Base.conj(c::Parity)=Parity(c.charge)
Base.one(c::Parity)=Parity(false)
Base.one(::Type{Parity})=Parity(false)

immutable ZNCharge{N} <: Abelian
    charge::Int
    ZNCharge(charge::Int)=new(mod(charge,N))
end

*{N}(c1::ZNCharge{N},c2::ZNCharge{N})=ZNCharge{N}(c1.charge+c2.charge)
Base.conj{N}(c::ZNCharge{N})=ZNCharge{N}(-c.charge)
Base.one{N}(c::ZNCharge{N})=ZNCharge{N}(0)
Base.one{N}(::Type{ZNCharge{N}})=ZNCharge{N}(0)

immutable U1Charge <: Abelian
    charge::Int
end

*(c1::U1Charge,c2::U1Charge)=U1Charge(c1.charge+c2.charge)
Base.conj(c::U1Charge)=U1Charge(-c.charge)
Base.one(c::U1Charge)=U1Charge(0)
Base.one(::Type{U1Charge})=U1Charge(0)

# Nonabelian Sectors:
#---------------------
# to be done


text:
using DynMultiply
using Base.Test

# first test set
A = rand(1000, 1000)
B = rand(1000, 1000)
v = rand(1000)

@test dynmultiply(A, B, v) ≈ A * B * v
@test dynmultiplystr(A, B, v) == "(A1*(A2*A3))"
@test dynmultiply_plan(A, B, v).data == [0,0,0,1,1,1]

# second test set
A = rand(1000, 1)
B = rand(1, 1000)
C = rand(1000, 50)
D = rand(50, 1)

@test dynmultiply(A, B, C, D) ≈ A * B * C * D
@test dynmultiplystr(A, B, C, D) == "(A1*((A2*A3)*A4))"
@test dynmultiply_plan(A, B, C, D).data == [0,0,0,0,1,1,1,1,2,1]


text:
"""
测试上wuxx的 A V_3 Sb_5
"""

using Plots
using MARY_fRG.Drawers
using MARY_fRG.Fermi.Kagome
using MARY_fRG.Fermi.Surface
using MARY_fRG.FlowEquation

#阻止画图
ENV["GKSwstype"] = "100"

"""
获取极低温的时候的bubble
"""
function get_ult_bubb(Γ4; maxnum=2000000)
    maxi, mini = engpeak_to_surface(Γ4)
    println(maxi, " ", mini)
    Γ4_ult = fliter_away_surface(Γ4, 2.0)
    Γ4_ult = refine_to_surface(Γ4_ult)
    while length(Γ4_ult.ltris) < maxnum
        Γ4_ult = refine_to_surface(Γ4_ult)
    end
    #plt = draw_points([tri.center for tri in Γ4_ult.ltris])
    #@savepng plt "ult_surface"
    maxi, mini = engpeak_to_surface(Γ4_ult)
    println(maxi, " ", mini)
    bubb_pp_u, bubb_fs_u, bubb_ex_u = all_bubble_tf_ult_mt(Γ4_ult, maxi)
    return bubb_pp_u, bubb_fs_u, bubb_ex_u
end

function run_tf()
    model = wuxx_av3sb5_kagome_lattice()
    Γ4 = TFGamma4(
        model, 6.0, 24, 233
    )
    Γ4.V .+= get_wuxx_U_mt(4.0, 2.0, 2.0, Γ4.patches[1], Γ4.patches[2])
    println(
        all(isfinite.(Γ4.V))
    )
    plt = heatmap(Γ4.V[1, 1, 1, 1, :, :, 1])
    png(plt, "Gamma41")
    #return
    lval = 0.
    lstep = 0.01
    bpu, bfu, beu = get_ult_bubb(Γ4)
    #
    #plt = heatmap(1e5*bpu.V[1, 1, 1, 1, 1, :, :])
    #png(plt, "wuxx_bpu1")
    #plt = heatmap(1e5*bpu.V[2, 2, 2, 2, 1, :, :])
    #png(plt, "wuxx_bpu2")
    #plt = heatmap(1e5*bpu.V[1, 2, 2, 1, 1, :, :])
    #png(plt, "wuxx_bpu3")
    #plt = heatmap(1e5*bpu.V[2, 1, 1, 2, 1, :, :])
    #png(plt, "wuxx_bpu4")
    ##
    #plt = heatmap(1e5*bfu.V[1, 1, 1, 1, 1, :, :])
    #png(plt, "wuxx_bfu1")
    #plt = heatmap(1e5*bfu.V[2, 2, 2, 2, 1, :, :])
    #png(plt, "wuxx_bfu2")
    #plt = heatmap(1e5*bfu.V[1, 2, 2, 1, 1, :, :])
    #png(plt, "wuxx_bfu3")
    #plt = heatmap(1e5*bfu.V[2, 1, 1, 2, 1, :, :])
    #png(plt, "wuxx_bfu4")
    ##
    #plt = heatmap(1e5*beu.V[1, 1, 1, 1, 1, :, :])
    #png(plt, "wuxx_beu1")
    #plt = heatmap(1e5*beu.V[2, 2, 2, 2, 1, :, :])
    #png(plt, "wuxx_beu2")
    #plt = heatmap(1e5*beu.V[1, 2, 2, 1, 1, :, :])
    #png(plt, "wuxx_beu3")
    #plt = heatmap(1e5*beu.V[2, 1, 1, 2, 1, :, :])
    #png(plt, "wuxx_beu4")
    @info "ult结束"
    for idx in 1:1:5001
        bubb_pp, bubb_fs, bubb_ex = all_bubble_tf_mt(Γ4, lval)
        dl = dl_tf_mix_ult_mt(Γ4, bubb_pp, bubb_fs, bubb_ex, bpu, bfu, beu)
        Γ4.V .+= dl .* lstep
        if idx % 50 == 1
            plt = heatmap(Γ4.V[1, 1, 1, 1, :, :, 1])
            png(plt, "Gamma41111"*string(idx))
            plt = heatmap(Γ4.V[2, 2, 2, 2, :, :, 1])
            png(plt, "Gamma42222"*string(idx))
            plt = heatmap(Γ4.V[1, 2, 2, 1, :, :, 1])
            png(plt, "Gamma41221"*string(idx))
            plt = heatmap(Γ4.V[2, 1, 1, 2, :, :, 1])
            png(plt, "Gamma42112"*string(idx))
        end
        lval += lstep
    end
end

run_tf()


text:
"""
alert_REPL!(;[duration=2.0], [message="Done!"])

Wraps all code passed to the REPL in an `@alert` macro with the given arguments.
Hence, if anything you run in the REPL takes longer than `duration` seconds, an
alert notification will be displayed. You can set the duration to `Inf` to turn
off the notification.

> NOTE: `onerror` must always be false; therefore no alert will be shown when a
command errors. This is a limitation of REPL design. All REPL statements are
expected to be top-level statements and inserting a `try/catch` block
automatically would prevent the assignment of global variables at the REPL.
"""
function alert_REPL!(;args...)
    if VERSION < v"1.5"
        @error "Requires Julia 1.5 or higher"
        return
    end

    repl_alert_options[] = __at_alert_options__(;args...)[1:2]
    dur = repl_alert_options[][1]
    if dur !== nothing && !isinf(dur)
        secs = @sprintf("%1.2f s", repl_alert_options[][1])
        @info "Alert will be sent if REPL line takes longer than $secs to complete. See "*
            "documentation for `Alert.alert_REPL!`"
    else
        @info "REPL alerts have been turned off."
    end
end

const repl_alert_options = Ref((Inf,"Done!"))
function with_repl_alert(ex)
    dur = repl_alert_options[][1]
    if dur !== nothing && !isinf(dur)
        Expr(:toplevel, quote
            @alert(duration=$(repl_alert_options[][1]),
                message=$(repl_alert_options[][2]),
                onerror=false,
                begin; $ex; end)
        end)
    else
        ex
    end
end

function init_alert_REPL!()
    if VERSION >= v"1.5"
        if isdefined(Base, :active_repl_backend)
            pushfirst!(Base.active_repl_backend.ast_transforms, with_repl_alert)
        else
            pushfirst!(REPL.repl_ast_transforms, with_repl_alert)
        end
    end
end


text:
if Sys.iswindows()
    LIBPATH_env = "PATH"
    LIBPATH_default = ""
    pathsep = ';'
    binlib = "bin"
    shlib_ext = "dll"
elseif Sys.isapple()
    LIBPATH_env = "DYLD_FALLBACK_LIBRARY_PATH"
    LIBPATH_default = "~/lib:/usr/local/lib:/lib:/usr/lib"
    pathsep = ':'
    binlib = "lib"
    shlib_ext = "dylib"
else
    LIBPATH_env = "LD_LIBRARY_PATH"
    LIBPATH_default = ""
    pathsep = ':'
    binlib = "lib"
    shlib_ext = "so"
end

function append_libpath(paths::Vector{<:AbstractString})
    return join(vcat(paths..., get(ENV, LIBPATH_env, String[])), pathsep)
end

function capture_output(cmd::Cmd; verbose::Bool = false)
    out_pipe = Pipe()
    if verbose
        ld_env = filter(e -> startswith(e, "LBT_") || startswith(e, "LD_") || startswith(e, "DYLD_"), something(cmd.env, String[]))
        @info("Running $(basename(cmd.exec[1]))", ld_env)
    end
    p = run(pipeline(ignorestatus(cmd), stdout=out_pipe, stderr=out_pipe), wait=false)
    close(out_pipe.in)
    output = @async read(out_pipe, String)
    wait(p)
    return fetch(output)
end

cc = something(
    get(ENV, "CC", nothing),
    Sys.which("cc"),
    Sys.which("gcc"),
    Sys.which("clang"),
)

# Escape paths on windows by using forward slash so we don't get into escaping hell
pathesc(path::String) = replace(path, "\\" => "/")

@static if Sys.isfreebsd()
    make = "gmake"
else
    make = "make"
end

needs_m32() = startswith(capture_output(`$(cc) -dumpmachine`), "x86_64") && Sys.WORD_SIZE == 32

# Build blastrampoline into a temporary directory, and return that
blastrampoline_build_dir = nothing
function get_blastrampoline_dir()
    if blastrampoline_build_dir !== nothing
        return blastrampoline_build_dir
    end

    cflags_add = needs_m32() ? "-m32" : ""
    dir = mktempdir()
    srcdir = joinpath(dirname(@__DIR__), "src")
    run(`$(make) -sC $(pathesc(srcdir)) CFLAGS_add=$(cflags_add) ARCH=$(Sys.ARCH) clean`)
    run(`$(make) -sC $(pathesc(srcdir)) CFLAGS_add=$(cflags_add) ARCH=$(Sys.ARCH) install builddir=$(pathesc(dir))/build prefix=$(pathesc(dir))/output`)
    global blastrampoline_build_dir = joinpath(dir, "output")
    return blastrampoline_build_dir
end

# Keep these in sync with `src/libblastrampoline_internal.h`
struct lbt_library_info_t
    libname::Cstring
    handle::Ptr{Cvoid}
    suffix::Cstring
    active_forwards::Ptr{UInt8}
    interface::Int32
    f2c::Int32
end
struct LBTLibraryInfo
    libname::String
    handle::Ptr{Cvoid}
    suffix::String
    active_forwards::Vector{UInt8}
    interface::Int32
    f2c::Int32

    LBTLibraryInfo(x::lbt_library_info_t, num_symbols::UInt32) = new(unsafe_string(x.libname), x.handle, unsafe_string(x.suffix), unsafe_wrap(Vector{UInt8}, x.active_forwards, div(num_symbols,8)+1), x.interface, x.f2c)
end
const LBT_INTERFACE_LP64 = 32
const LBT_INTERFACE_ILP64 = 64
const LBT_F2C_PLAIN = 0

struct lbt_config_t
    loaded_libs::Ptr{Ptr{lbt_library_info_t}}
    build_flags::UInt32
    exported_symbols::Ptr{Cstring}
    num_exported_symbols::UInt32
end
const LBT_BUILDFLAGS_F2C_CAPABLE = 0x02

function lbt_forward(handle, path; clear::Bool = false, verbose::Bool = false)
    ccall(dlsym(handle, :lbt_forward), Int32, (Cstring, Int32, Int32), path, clear ? 1 : 0, verbose ? 1 : 0)
end

function lbt_get_config(handle)
    return unsafe_load(ccall(dlsym(handle, :lbt_get_config), Ptr{lbt_config_t}, ()))
end

function lbt_get_num_threads(handle)
    return ccall(dlsym(handle, :lbt_get_num_threads), Int32, ())
end

function lbt_set_num_threads(handle, nthreads)
    return ccall(dlsym(handle, :lbt_set_num_threads), Cvoid, (Int32,), nthreads)
end

function lbt_get_forward(handle, symbol_name, interface, f2c = LBT_F2C_PLAIN)
    return ccall(dlsym(handle, :lbt_get_forward), Ptr{Cvoid}, (Cstring, Int32, Int32), symbol_name, interface, f2c)
end

function lbt_set_forward(handle, symbol_name, addr, interface, f2c = LBT_F2C_PLAIN; verbose::Bool = false)
    return ccall(dlsym(handle, :lbt_set_forward), Int32, (Cstring, Ptr{Cvoid}, Int32, Int32, Int32), symbol_name, addr, interface, f2c, verbose ? 1 : 0)
end

function lbt_set_default_func(handle, addr)
    return ccall(dlsym(handle, :lbt_set_default_func), Cvoid, (Ptr{Cvoid},), addr)
end

function lbt_get_default_func(handle)
    return ccall(dlsym(handle, :lbt_get_default_func), Ptr{Cvoid}, ())
end

max_stars_repo_name: Gnimuc/libblastrampoline

max_issues_repo_name: Gnimuc/libblastrampoline

max_forks_repo_name: Gnimuc/libblastrampoline


text:
# ------------------------------------------------------------------------------------------
# #  Choosing an ODE Algorithm
#
# While the default algorithms, along with `alg_hints = [:stiff]`, will suffice in most
# cases, there are times when you may need to exert more control. The purpose of this part
# of the tutorial is to introduce you to some of the most widely used algorithm choices and
# when they should be used. The corresponding page of the documentation is the [ODE
# Solvers](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html) page which goes into
# more depth.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Diagnosing Stiffness
#
# One of the key things to know for algorithm choices is whether your problem is stiff.
# Let's take for example the driven Van Der Pol equation:
# ------------------------------------------------------------------------------------------

using DifferentialEquations
van! = @ode_def VanDerPol begin
  dy = μ*((1-x^2)*y - x)
  dx = 1*y
end μ

prob = ODEProblem(van!,[0.0,2.0],(0.0,6.3),1e6)

# ------------------------------------------------------------------------------------------
# One indicating factor that should alert you to the fact that this model may be stiff is
# the fact that the parameter is `1e6`: large parameters generally mean stiff models. If we
# try to solve this with the default method:
# ------------------------------------------------------------------------------------------

sol = solve(prob)

# ------------------------------------------------------------------------------------------
# Here it shows that maximum iterations were reached. Another thing that can happen is that
# the solution can return that the solver was unstable (exploded to infinity) or that `dt`
# became too small. If these happen, the first thing to do is to check that your model is
# correct. It could very well be that you made an error that causes the model to be
# unstable!
#
# If the model is the problem, then stiffness could be the reason. We can thus hint to the
# solver to use an appropriate method:
# ------------------------------------------------------------------------------------------

sol = solve(prob,alg_hints = [:stiff])

# ------------------------------------------------------------------------------------------
# Another way to understand stiffness is to look at the solution.
# ------------------------------------------------------------------------------------------

using Plots; gr()
sol = solve(prob,alg_hints = [:stiff],reltol=1e-6)
plot(sol,denseplot=false)

# ------------------------------------------------------------------------------------------
# Let's zoom in on the y-axis to see what's going on:
# ------------------------------------------------------------------------------------------

plot(sol,ylims = (-10.0,10.0))

# ------------------------------------------------------------------------------------------
# Notice how there are some extreme vertical shifts that occur. These vertical shifts are
# places where the derivative term is very large, and this is indicative of stiffness. This
# is an extreme example to highlight the behavior, but this general idea can be carried over
# to your problem. When in doubt, simply try timing using both a stiff solver and a non-
# stiff solver and see which is more efficient.
#
# To try this out, let's use BenchmarkTools, a package that let's us relatively reliably
# time code blocks.
# ------------------------------------------------------------------------------------------

using BenchmarkTools

# ------------------------------------------------------------------------------------------
# Let's compare the performance of non-stiff and stiff solvers on the Lorenz equation, which
# we saw in the last notebook, "ODEIntroduction".
#
# First, let's grab all our information about the Lorenz equation from the last notebook.
# ------------------------------------------------------------------------------------------

function lorenz!(du,u,p,t)
    σ,ρ,β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
end
u0 = [1.0,0.0,0.0]
p = (10,28,8/3)
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)

# ------------------------------------------------------------------------------------------
# And now, let's use the `@btime` macro from benchmark tools to compare the use of non-stiff
# and stiff solvers on this problem.
# ------------------------------------------------------------------------------------------

@btime solve(prob);

@btime solve(prob,alg_hints = [:stiff]);

# ------------------------------------------------------------------------------------------
# In this particular case, we can see that non-stiff solvers get us to the solution much
# more quickly.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## The Recommended Methods
#
# When picking a method, the general rules are as follows:
#
# - Higher order is more efficient at lower tolerances, lower order is more efficient at
# higher tolerances
# - Adaptivity is essential in most real-world scenarios
# - Runge-Kutta methods do well with non-stiff equations, Rosenbrock methods do well with
# small stiff equations, BDF methods do well with large stiff equations
#
# While there are always exceptions to the rule, those are good guiding principles. Based on
# those, a simple way to choose methods is:
#
# - The default is `Tsit5()`, a non-stiff Runge-Kutta method of Order 5
# - If you use low tolerances (`1e-8`), try `Vern7()` or `Vern9()`
# - If you use high tolerances, try `BS3()`
# - If the problem is stiff, try `Rosenbrock23()`, `Rodas4()`, or `CVODE_BDF()`
#
# (This is a simplified version of the default algorithm chooser)
#
# ## Comparison to other Software
#
# If you are familiar with MATLAB, SciPy, or R's DESolve, here's a quick translation start
# to have transfer your knowledge over.
#
# - ode23 –> BS3()
# - ode45/dopri5 –> DP5(), though in most cases Tsit5() is more efficient
# - ode23s –> Rosenbrock23(), though in most cases Rodas4() is more efficient
# - ode113 –> CVODE_Adams(), though in many cases Vern7() is more efficient
# - dop853 –> DP8(), though in most cases Vern7() is more efficient
# - ode15s/vode –> CVODE_BDF(), though in many cases Rodas4() or radau() are more efficient
# - ode23t –> Trapezoid() for efficiency and GenericTrapezoid() for robustness
# - ode23tb –> TRBDF2
# - lsoda –> lsoda() (requires Pkg.add("LSODA"); using LSODA)
# - ode15i –> IDA(), though in many cases Rodas4() can handle the DAE and is significantly
# more efficient
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using NiLangCore
using NiLangCore: compile_ex, dual_ex, precom_ex, memkernel, render_arg, check_args

using Test
import Base: +, -

value(x) = x
NiLangCore.chfield(x::T, ::typeof(value), y::T) where T = y

function add(a!::Number, b::Number)
    a!+b, b
end

function neg(b::Number)
    -b
end
@selfdual neg

@i function add(a!, b)
    add(a! |> value, b |> value)
end

function sub(a!::Number, b::Number)
    a!-b, b
end

@i function sub(a!, b)
    sub(a! |> value, b |> value)
end

@dual add sub

function XOR(a!::Integer, b::Integer)
    xor(a!, b), b
end

@selfdual XOR
#@nograd XOR

@testset "boolean" begin
    x = false
    @instr x ⊻= true
    @test x
    @instr x ⊻= true || false
    @test !x
    @instr x ⊻= true && false
    @instr x ⊻= !false
    @test x
end

@testset "@dual" begin
    @test isreversible(add, Tuple{Any,Any})
    @test isreversible(sub, Tuple{Any,Any})
    @test !isreflexive(add)
    @test ~(add) == sub
    a=2.0
    b=1.0
    @instr add(a, b)
    @test a == 3.0
    args = (1,2)
    @instr add(args...)
    @test args == (3,2)
    @instr sub(a, b)
    @test a == 2.0
    @test check_inv(add, (a, b))
    @test isprimitive(add)
    @test isprimitive(sub)
end

@testset "@selfdual" begin
    @test !isreversible(XOR, Tuple{Any, Any})
    @test !isreversible(~XOR, Tuple{Any, Any})
    @test isreversible(~XOR, Tuple{Integer, Integer})
    @test isreversible(XOR, Tuple{Integer, Integer})
    @test isreflexive(XOR)
    @test isprimitive(XOR)
    @test ~(XOR) == XOR
    a=2
    b=1
    @instr XOR(a, b)
    @test a == 3
    @instr XOR(a, b)
    @test a == 2
end

@testset "+=, -=" begin
    x = 1.0
    y = 1.0
    @instr PlusEq(exp)(y, x)
    @test x ≈ 1
    @test y ≈ 1+exp(1.0)
    @instr (~PlusEq(exp))(y, x)
    @test x ≈ 1
    @test y ≈ 1
end

@testset "+= and const" begin
    x = 0.5
    @instr x += π
    @test x == 0.5+π
    @instr x += log(π)
    @test x == 0.5 + π + log(π)
    @instr x += log(π)/2
    @test x == 0.5 + π + 3*log(π)/2
    @instr x += log(2*π)/2
    @test x == 0.5 + π + 3*log(π)/2 + log(2π)/2
end

@testset "+= keyword functions" begin
    g(x; y=2) = x^y
    z = 0.0
    x = 2.0
    @instr z += g(x; y=4)
    @test z == 16.0
end

@testset "constant value" begin
    @test @const 2 == 2
    @test NiLangCore._isconst(:(@const grad(x)))
end

@testset "+=, -=, *=, /=" begin
    @test compile_ex(@__MODULE__, :(x += y * z), NiLangCore.CompileInfo()).args[1].args[2] == :($PlusEq(*)(x, y, z))
    @test compile_ex(@__MODULE__, dual_ex(@__MODULE__, :(x -= y * z)), NiLangCore.CompileInfo()).args[1].args[2] ==  :($PlusEq(*)(x, y, z))
    @test compile_ex(@__MODULE__, :(x /= y * z), NiLangCore.CompileInfo()).args[1].args[2] == :($DivEq(*)(x, y, z))
    @test compile_ex(@__MODULE__, dual_ex(@__MODULE__, :(x *= y * z)), NiLangCore.CompileInfo()).args[1].args[2] ==  :($DivEq(*)(x, y, z))
    @test ~MulEq(*) == DivEq(*)
    @test ~DivEq(*) == MulEq(*)
    function (g::MulEq)(y, a, b)
        y * g.f(a, b), a, b
    end

    function (g::DivEq)(y, a, b)
        y / g.f(a, b), a, b
    end
    a, b, c = 1.0, 2.0, 3.0
    @instr a *= b + c
    @test a == 5.0
    @instr a /= b + c
    @test a == 1.0
end

@testset "shared read write check" begin
    for (x, y) in [
        (:((-x[3].g' |> NEG).k[5]) , :((x[3]).g.k[5]))
        (:((-(x |> subarray(3)).g' |> NEG).k[5]) , :((x[3]).g.k[5]))
        (:(@skip! x.g) , nothing)
        (:(@const x .|> g) , :x)
        (:(cos.(x[2])) , nothing)
        (:(cos(x[2])) , nothing)
        (:((x |> g)...) , :x)
        (:((x |> g, y.:1)) , [:x, :(y.:1)])
        (:((x |> g, y |> tget(1))) , [:x, :(y[1])])]
        @test memkernel(deepcopy(x)) == y
        @test render_arg(deepcopy(x)) == x
    end
    @test render_arg(:(x.y.[2:3])) == :(x.y |> subarray(2:3))
    @test memkernel(:(x.y |> subarray(2:3))) == (:(x.y[2:3]))
    @test render_arg(:(x.y.[2:3] |> value)) == :(x.y |> subarray(2:3) |> value)
    @test memkernel(:(x.y |> subarray(2:3) |> value)) == :(x.y[2:3])

    @test_throws InvertibilityError check_args([:a, :(a |> grad)])
    @test check_args([:(a.x), :(a.g |> grad)]) isa Nothing
    @test_throws InvertibilityError check_args([:(a.x), :(b[3]), :(b[3])])
    @test_throws InvertibilityError check_args([:(a.x), :((b, a.x))]) isa Nothing
    # TODO: check variable on the same tree, like `a.b` and `a`
end

@testset "tailn" begin
    @test NiLangCore.tailn((1,2,3), Val(1), (1,2,3)) == (2,3)
    @test NiLangCore.tailn((1,2,3), Val(0), (1,2,3)) == (1,2,3)
    @test NiLangCore.tailn(1, Val(1), ()) == ()
    @test NiLangCore.tailn(1, Val(0), (1,)) == (1,)
    @test NiLangCore.tailn((1,), Val(0), ((1,),)) == ((1,),)
end

@testset "dual type" begin
    struct AddX{T}
        x::T
    end
    struct SubX{T}
        x::T
    end
    @dualtype AddX SubX
    @dualtype AddX SubX
    @i function (f::AddX)(x::Real) end
    @test hasmethod(AddX(3), Tuple{Real})
    @test hasmethod(SubX(3), Tuple{Real})
    for (TA, TB) in [(AddX, SubX), (MulEq, DivEq), (XorEq, XorEq), (PlusEq, MinusEq)]
        @test invtype(TA) == TB
        @test invtype(TA{typeof(*)}) == TB{typeof(*)}
        @test invtype(TB) == TA
        @test invtype(TB{typeof(*)}) == TA{typeof(*)}
    end
end


text:
function program_path()
    (isinteractive() || isempty(Base.PROGRAM_FILE)) && return nothing
    return dirname(abspath(Base.PROGRAM_FILE))
end

function app_root_path()
    prog_path = program_path()

    return isnothing(prog_path) ? pwd() : prog_path
end

function pathname_configs(url_base_pathname, requests_pathname_prefix, routes_pathname_prefix)

    raise_error = (s) -> error("""
    $s This is ambiguous.
    To fix this, set `routes_pathname_prefix` instead of `url_base_pathname`.

    Note that `requests_pathname_prefix` is the prefix for the AJAX calls that
    originate from the client (the web browser) and `routes_pathname_prefix` is
    the prefix for the API routes on the backend (as defined within HTTP.jl).
    `url_base_pathname` will set `requests_pathname_prefix` and
    `routes_pathname_prefix` to the same value.
    If you need these to be different values then you should set
    `requests_pathname_prefix` and `routes_pathname_prefix`,
    not `url_base_pathname`.
    """)

    if !isnothing(url_base_pathname) && !isnothing(requests_pathname_prefix)
        raise_error("You supplied `url_base_pathname` and `requests_pathname_prefix`")
    end

    if !isnothing(url_base_pathname) && !isnothing(routes_pathname_prefix)
        raise_error("You supplied `url_base_pathname` and `routes_pathname_prefix`")
    end

    if !isnothing(url_base_pathname) && isnothing(routes_pathname_prefix)
        routes_pathname_prefix = url_base_pathname
    elseif isnothing(routes_pathname_prefix)
        routes_pathname_prefix = "/"
    end

    !startswith(routes_pathname_prefix, "/") && error("routes_pathname_prefix` needs to start with `/`")
    !endswith(routes_pathname_prefix, "/") && error("routes_pathname_prefix` needs to end with `/`")

    app_name = dash_env("APP_NAME")

    if isnothing(requests_pathname_prefix) && !isnothing(app_name)
        requests_pathname_prefix = "/" * app_name * routes_pathname_prefix
    elseif isnothing(requests_pathname_prefix)
        requests_pathname_prefix = routes_pathname_prefix
    end

    !startswith(requests_pathname_prefix, "/") &&
             error("requests_pathname_prefix` needs to start with `/`")
    !endswith(requests_pathname_prefix, routes_pathname_prefix) &&
             error("requests_pathname_prefix` needs to end with `routes_pathname_prefix`")

    return (url_base_pathname, requests_pathname_prefix, routes_pathname_prefix)
end


text:
### convolution.jl
#
# Copyright (C) 2021 Ido Kessler.
#
# Keywords: exact integer convolution
#
# This file is a part of ExactConvolution.jl.
#
# License is MIT.
#
### Commentary:
#
# This file contains implementation of convolutions using bigints multiplication.
#
### Code:

export exact_conv

include("conversion.jl")

"""
    exact_conv([T,] a, b)

Convolve arrays of integers `a` and `b`, and return the resulting array.

```math
exact\\_conv(a,b)[k] = \\sum_{i+j=k} a[i] \\cdot b[j]
```

# Notes:
- The implementation uses GMP's BigInt multiplication as the underling DSP engine (FFT for big enough inputs).
- T must be sufficient to support the result.
- Using negative values in `a` and `b` is supported but discouraged for runtime reasons.

# Examples
```jldoctest
julia> exact_conv(Int32, [1,2], [4,5,6])
4-element Vector{Int32}:
  4
 13
 16
 12
```
"""
function exact_conv end

exact_conv(a, b) = exact_conv(Int64, a, b)

function exact_conv(
    ::Type{T},
    a::AbstractArray{T1},
    b::AbstractArray{T2},
) where {T<:Integer,T1<:Integer,T2<:Integer}
    function make_positive(::Type{T}, v) where {T}
        min_v = min(minimum(v), 0) |> T
        return min_v, (min_v == 0 ? v : (v .- min_v))
    end

    min_a, a_pos = make_positive(T, a)
    min_b, b_pos = make_positive(T, b)

    res = positive_exact_conv(T, a_pos, b_pos)

    # handle the added -min_a * (b[...] - min_b)
    handle_negatives!(T, res, -min_a, b_pos)
    # handle the added -min_b * a[...]
    handle_negatives!(T, res, -min_b, a)

    res
end

# Only accept positive values
function positive_exact_conv(::Type{T}, a, b)::Vector{T} where {T<:Integer}
    nbits(x) = ndigits(x; base = 2)
    max_bits = nbits ∘ maximum

    lg2_max_mult = max_bits(a) + max_bits(b)
    lg2_max_len = [a, b] .|> length |> minimum |> nbits
    nb = lg2_max_mult + lg2_max_len
    nb = 5 * ceil(Int, nb // 5)

    big_a = array2bigint(a, nb)
    big_b = array2bigint(b, nb)
    ab = big_a * big_b
    n = length(a) + length(b) - 1
    bigint2array(T, ab, n, nb)
end

function handle_negatives!(::Type{T}, res, v1_added, v2_pos)::Nothing where {T}
    v1_added == 0 && return
    # min_v1 * v2_[i:j] was added
    n, n2 = length(res), length(v2_pos)

    v1_added = T(v1_added)
    extra = T(0)
    # v2[i] only participate in the conv from index i
    for (i, v) ∈ zip(1:n2, v2_pos)
        extra += v * v1_added
        res[i] -= extra
    end
    res[n2+1:end] .-= extra

    removed_extra = T(0)
    # v2[end-i] only participate in the conv upto index end-i
    for (i, v) ∈ zip(n-n2+1:n, v2_pos)
        res[i] += removed_extra
        removed_extra += v * v1_added
    end
    return
end


text:
module Seregmort

using JSONStat, HTTP, DataStructures, DataFrames, PyCall, PyPlot, Statistics
import JSON
#import Mortchartgen: grpprop

const shpreader = PyNULL()
const ccrs = PyNULL()

function __init__()
	copy!(shpreader, pyimport("cartopy.io.shapereader"))
	copy!(ccrs, pyimport("cartopy.crs"))
end

dfgrp_agemean(df, grpcol, f = mean) =
	combine(groupby(df, grpcol), x -> DataFrame(value = f(x[!, :value])))
dfgrp_sum(df, grpcol, f = sum) =
	combine(groupby(df, grpcol), x -> DataFrame(value = f(x[!, :value]), value_1 = f(x[!, :value_1])))

function grpprop(numframe_sub, denomframe_sub, grpcol, agemean)
	numdenomframe_sub =
		innerjoin(numframe_sub, denomframe_sub, on = grpcol, makeunique = true)
	if agemean
		propfr_agesp = DataFrame()
		propfr_agesp[!, grpcol] = numdenomframe_sub[!, grpcol]
		propfr_agesp[!, :value] = numdenomframe_sub[!, :value]./numdenomframe_sub[!, :value_1]
		return dfgrp_agemean(propfr_agesp, grpcol)
	else
		numdenomgrp = dfgrp_sum(numdenomframe_sub, grpcol)
		propfr_agegr = DataFrame()
		propfr_agegr[!, grpcol] = numdenomgrp[!, grpcol]
		propfr_agegr[!, :value] = numdenomgrp[!, :value]./numdenomgrp[!, :value_1]
		return propfr_agegr
	end
end

export metadata, allregions, unchanged_regions, ndeaths, npop, munis_incounty,
propplotyrs_dict, propplotyrs, propscatsexes_dict, propscatsexes,propmap_dict, propmap,
catot_yrsdict, capop_yrsdict, catot_mapdict, capop_mapdict,
threep, fourp, fivep

MORTURL = "http://api.scb.se/OV0104/v1/doris/sv/ssd/START/HS/HS0301/DodaOrsak"
POPURL = "http://api.scb.se/OV0104/v1/doris/sv/ssd/START/BE/BE0101/BE0101A/BefolkningNy"
MAINPATH = normpath(@__DIR__, "..")
DATAPATH = normpath(MAINPATH, "data")
G_UNITS = JSON.parsefile(normpath(DATAPATH, "g_units.json"))

PyDict(matplotlib["rcParams"])["axes.formatter.use_locale"] = true
matplotlib[:style][:use]("ggplot")

function scb_to_unit(scb)
	scbform = *("SE/", rpad(scb, 9, '0'))
	if scbform in keys(G_UNITS)
		return G_UNITS[scbform]
	else
		return 0
	end
end

function metadata(url = MORTURL)
	req = HTTP.request("GET", url)
	JSON.parse(String(req.body), dicttype = DataStructures.OrderedDict)
end

function causealias(cause, dim)
	if cause == "POP"
		return lowercase(dim["ContentsCode"]["category"]["label"]["BE0101N1"])
	else
		return dim["Dodsorsak"]["category"]["label"][cause]
	end
end

regalias(region, dim) = lstrip(replace(dim["Region"]["category"]["label"][region], region => ""))

sexalias(sex, dim) = dim["Kon"]["category"]["label"][sex]

function agesplitter(age)
	if occursin("-", age)
	    return split(age, "-")
	else
	    return [age]
	end
end

function ageslice(startage, endage, agemean)
	ages = allages()
	startind = indexin([startage], ages)[1]
	endind = indexin([endage], ages)[1]
	if agemean
		agemeanstr = " medel över åldrar"
	else
		agemeanstr = ""
	end
	if startage == endage
		alias = *(replace(startage, "-" => '\u2013'), agemeanstr)
	else
	    alias = *(agesplitter(startage)[1], "\u2013", agesplitter(endage)[end], agemeanstr)
	end
	agelist = ages[startind:endind]
	Dict("agelist" => agelist, "alias" => alias)
end

function allages(ageformat = "mort")
	if ageformat == "mort"
		startint = "0"
		startages = vcat(1, collect(5:5:85))
		endages = vcat(4, collect(9:5:89))
		endint = "90+"
	elseif ageformat == "pop"
		startint = "-4"
		startages = collect(5:5:95)
		endages = collect(9:5:99)
		endint = "100+"
	end
	midints = ["$x-$y" for (x, y) in zip(startages, endages)]
	vcat(startint, midints, endint)
end

function allregions(level, metadict)
	regvalues = metadict["variables"][1]["values"]
	if level == "county"
		return filter(is_county, regvalues)
	elseif level == "municipality"
		return filter(is_municipality, regvalues)
	end
end

unchanged_county(region) = !(region[1:2] in ["11"; "12"; "14"; "15"; "16"])
unchanged_regions(level, metadict) = filter(unchanged_county, allregions(level, metadict))

yearrange(startyear = 1969, endyear = 1996) = ["$(year)" for year in startyear:endyear]

is_county(region) = length(region) == 2 && region != "00"
is_municipality(region) = length(region) == 4

function munis_incounty(county, metadict)
	regvalues = metadict["variables"][1]["values"]
	filter((x) -> (is_municipality(x) && startswith(x, county)), regvalues)
end

function mortreqjson(regvalues, causevalues,  agevalues = allages(),
	sexvalues = ["1"; "2"], yearvalues = yearrange())
	if is_county(regvalues[1])
		regfilter = "vs:RegionLän"
	else
		regfilter = "vs:RegionKommun95"
	end
	if occursin("-", causevalues[1])
		causefilter = "agg:DödsorsakKapitel"
	else
		causefilter = "item"
	end

	Dict("response" => Dict("format" => "json-stat"),
		"query" => [
		Dict("selection" => Dict("filter" => regfilter, "values" => regvalues),
			"code" => "Region");
		Dict("selection" => Dict("filter" => causefilter, "values" => causevalues),
			"code" => "Dodsorsak");
		Dict("selection" => Dict("filter" => "item", "values" => agevalues),
			"code" => "Alder");
		Dict("selection" => Dict("filter" => "item", "values" => sexvalues),
			"code" => "Kon");
		Dict("selection" => Dict("filter" => "item", "values" => yearvalues),
			"code" => "Tid")
		]
		)
end

function popreqjson(regvalues, agevalues = allages("pop"),
	sexvalues = ["1"; "2"], yearvalues = yearrange())
	if is_county(regvalues[1])
		regfilter = "vs:RegionLän07"
	else
		regfilter = "vs:RegionKommun07"
	end

	Dict("response" => Dict("format" => "json-stat"),
		"query" => [
		Dict("selection" => Dict("filter" => regfilter, "values" => regvalues),
			"code" => "Region");
		Dict("selection" => Dict("filter" => "agg:Ålder5år", "values" => agevalues),
			"code" => "Alder");
		Dict("selection" => Dict("filter" => "item", "values" => sexvalues),
			"code" => "Kon");
		Dict("selection" => Dict("filter" => "item", "values" => ["BE0101N1"]),
			"code" => "ContentsCode");
		Dict("selection" => Dict("filter" => "item", "values" => yearvalues),
			"code" => "Tid")
		]
		)
end

function ndeaths(regvalues, causevalues;  agevalues = allages(),
	sexvalues = ["1"; "2"], yearvalues = yearrange())
	qjson = JSON.json(mortreqjson(regvalues, causevalues, agevalues, sexvalues, yearvalues))
	req = HTTP.request("POST", MORTURL, [], qjson)
	reqjsonstat = JSON.parse(String(req.body),
		dicttype = DataStructures.OrderedDict)
	readjsonbundle(reqjsonstat)["dataset"]
end

function npop(regvalues; agevalues = allages("pop"),
	sexvalues = ["1"; "2"], yearvalues = yearrange())
	qjson = JSON.json(popreqjson(regvalues, agevalues, sexvalues, yearvalues))
	req = HTTP.request("POST", POPURL, [], qjson)
	reqjsonstat = JSON.parse(String(req.body),
		dicttype = DataStructures.OrderedDict)
	readjsonbundle(reqjsonstat)["dataset"]
end

dfarrmatch(col, arr) = map((x) -> in(x, arr), Vector(col))

subframe_sray(df, sex, region, agelist, years) = df[((df[!, :Kon].==sex)
	.& (df[!, :Region].==region) .& (dfarrmatch(df[!, :Alder], agelist))
	.& (dfarrmatch(df[!, :Tid], years))), :]

function prop_timegrp(numframe, denomframe, sex, region, agelist, years, agemean)
	numframe_sub = subframe_sray(numframe, sex, region, agelist, years)
	denomframe_sub = subframe_sray(denomframe, sex, region, agelist, years)
	grpprop(numframe_sub, denomframe_sub, :Tid, agemean)
end

propplotyrs_dict(startage, endage, pardict,
	agemean = false, years = yearrange(), sexes = ["2", "1"]) =
propplotyrs(pardict["numframe"], pardict["denomframe"], pardict["numdim"],
	pardict["denomdim"], pardict["numcause"], pardict["denomcause"], pardict["region"],
	startage, endage, agemean, years, sexes)

function propplotyrs(numframe, denomframe, numdim, denomdim, numcause, denomcause,
        region, startage, endage, agemean, years, sexes)
	numcauseal = causealias(numcause, numdim)
	denomcauseal = causealias(denomcause, denomdim)
	regal = regalias(region, numdim)
	ages = ageslice(startage, endage, agemean)
	ageal = ages["alias"]
	agelist = ages["agelist"]
	yrints = map((x)->parse(Int, x), years)
	for sex in sexes
		sexal = sexalias(sex, numdim)
		propframe = prop_timegrp(numframe, denomframe,
			sex, region, agelist, years, agemean)
		plot(yrints, propframe[!, :value], label = sexal, "-*")
	end
	legend(framealpha = 0.5)
	xlim(yrints[1], yrints[end])
	ylim(ymin = 0)
	title("Döda $(numcauseal)/$(denomcauseal)\n$(ageal) $(regal)")
end

subframe_sa(df, sex, agelist) = df[((df[!, :Kon].==sex)
	.& (dfarrmatch(df[!, :Alder], agelist))), :]

function prop_reggrp(numframe, denomframe, sex, agelist, agemean)
	numframe_sub = subframe_sa(numframe, sex, agelist)
	denomframe_sub = subframe_sa(denomframe, sex, agelist)
	grpprop(numframe_sub, denomframe_sub, :Region, agemean)
end

propscatsexes_dict(startage, endage, pardict, agemean = false) =
propscatsexes(pardict["numframe"], pardict["denomframe"], pardict["numdim"],
	pardict["denomdim"], pardict["numcause"], pardict["denomcause"],
	startage, endage, agemean)

function propscatsexes(numframe, denomframe, numdim, denomdim, numcause, denomcause,
        startage, endage, agemean = false)
	numcauseal = causealias(numcause, numdim)
	denomcauseal = causealias(denomcause, denomdim)
	ages = ageslice(startage, endage, agemean)
	ageal = ages["alias"]
	agelist = ages["agelist"]
	yrints = map((x)->parse(Int, x), numframe[!, :Tid])
	startyear = minimum(yrints)
	endyear = maximum(yrints)
	sexframes = Dict()
	for sex in ["2", "1"]
		sexframes[sex] = Dict()
		sexframes[sex]["alias"] = sexalias(sex, numdim)
		sexframes[sex]["propframe"] = prop_reggrp(numframe, denomframe,
			sex, agelist, agemean)
	end
	femprop = sexframes["2"]["propframe"][!, :value]
	maleprop = sexframes["1"]["propframe"][!, :value]
	regcodes = sexframes["2"]["propframe"][!, :Region]
	regals = map((x)->regalias(x, numdim), regcodes)
	scatter(femprop, maleprop)
	for (i, regcode) in enumerate(regcodes)
		annotate(regcode, (femprop[i], maleprop[i]))
	end
	xlabel(sexframes["2"]["alias"])
	ylabel(sexframes["1"]["alias"])
	axminimum = 0.95
	axmaximum = 1.05
	xlim(minimum(femprop)*axminimum, maximum(femprop)*axmaximum)
	ylim(minimum(maleprop)*axminimum, maximum(maleprop)*axmaximum)
	title(*("Döda $(numcauseal)/$(denomcauseal)\n",
		"$(ageal) $(startyear)\u2013$(endyear)"))
	DataFrame(code = regcodes, alias = regals, femprop = femprop,
		maleprop = maleprop)
end

perc_round(value) = replace("$(round(value; digits=4))", "." => ",")

threep(prop) =
	[
	Dict("col" => "lightsalmon", "value" => quantile(prop, 1/3));
	Dict("col" => "tomato", "value" => quantile(prop, 2/3));
	Dict("col" => "red", "value" => quantile(prop, 1))
	]

fourp(prop) =
	[
	Dict("col" => "lightyellow", "value" => quantile(prop, 1/4));
	Dict("col" => "yellow", "value" => quantile(prop, 2/4));
	Dict("col" => "tomato", "value" => quantile(prop, 3/4));
	Dict("col" => "red", "value" => quantile(prop, 1))
	]

fivep(prop) =
	[
	Dict("col" => "lightyellow", "value" => quantile(prop, 1/5));
	Dict("col" => "yellow", "value" => quantile(prop, 2/5));
	Dict("col" => "orange", "value" => quantile(prop, 3/5));
	Dict("col" => "tomato", "value" => quantile(prop, 4/5));
	Dict("col" => "red", "value" => quantile(prop, 1))
	]

propmap_dict(startage, endage, sex, pardict, percfunc = threep, agemean = false) =
propmap(pardict["numframe"], pardict["denomframe"], pardict["numdim"],
	pardict["denomdim"], pardict["numcause"], pardict["denomcause"],
	startage, endage, sex, pardict["shapefname"], percfunc, agemean)

function propmap(numframe, denomframe, numdim, denomdim, numcause, denomcause,
        startage, endage, sex, shapefname, percfunc = threep, agemean = false)
	numcauseal = causealias(numcause, numdim)
	denomcauseal = causealias(denomcause, denomdim)
	sexal = sexalias(sex, numdim)
	ages = ageslice(startage, endage, agemean)
	ageal = ages["alias"]
	agelist = ages["agelist"]
	yrints = map((x)->parse(Int, x), numframe[!, :Tid])
	startyear = minimum(yrints)
	endyear = maximum(yrints)
	region_shp = shpreader[:Reader](shapefname)
	propframe = prop_reggrp(numframe, denomframe, sex, agelist, agemean)
	regcodes = propframe[!, :Region]
	regals = map((x)->regalias(x, numdim), regcodes)
	prop = propframe[!, :value]
	propdict = Dict(zip(regcodes, prop))
	units = map(scb_to_unit, regcodes)
	regdict = Dict(zip(units, regcodes))
	quantiles = percfunc(prop)
	proj = ccrs[:LambertConformal](central_longitude = 10, central_latitude = 52,
		standard_parallels = (35,65), false_easting = 4000000,
		false_northing = 2800000, globe = ccrs[:Globe](ellipse = "GRS80"))
	ax = plt[:axes](projection = proj)
	boundlist = []
	facecolor = "red"
	for region_rec in region_shp[:records]()
		regunit = region_rec[:attributes]["G_UNIT"]
		regend = region_rec[:attributes]["GET_END_YE"]
		if (regunit in keys(regdict) && regend > 1995)
			boundlist = vcat(boundlist, region_rec[:bounds])
			xmean = mean([boundlist[end][1];
				boundlist[end][3]])
			ymean = mean([boundlist[end][2];
				boundlist[end][4]])
			for quantile in quantiles
				if propdict[regdict[regunit]] <= quantile["value"]
					facecolor = quantile["col"]
					break
				end
			end
			ax[:add_geometries]([region_rec[:geometry]], proj,
				edgecolor = "black", facecolor = facecolor)
			ax[:annotate](regdict[regunit], (xmean, ymean),
				ha = "center")
		end
	end
	xminimum = minimum([bound[1] for bound in boundlist])
	xmaximum = maximum([bound[3] for bound in boundlist])
	yminimum = minimum([bound[2] for bound in boundlist])
	ymaximum = maximum([bound[4] for bound in boundlist])
	ax[:set_xlim](xminimum, xmaximum)
	ax[:set_ylim](yminimum, ymaximum)
	percpatches = []
	perclabels = []
	for (i, quantile) in enumerate(quantiles)
		percpatch = matplotlib[:patches][:Rectangle]((0, 0), 1, 1,
			facecolor = quantile["col"])
		percpatches = vcat(percpatches, percpatch)
		if i == 1
			perclabel = *("\u2265", perc_round(minimum(prop)),
				"\n\u2264", perc_round(quantile["value"]))
		else
			perclabel = *("\u2264", perc_round(quantile["value"]))
		end
		perclabels = vcat(perclabels, perclabel)
	end
	legend(percpatches, perclabels, loc = "upper left",
		framealpha = 0.75, bbox_to_anchor=(1,1))
	title(*("Döda $(numcauseal)/$(denomcauseal)\n",
		"$(sexal) $(ageal) $(startyear)\u2013$(endyear)"))
	subplots_adjust(right = 0.8)
	show()
	DataFrame(code = regcodes, alias = regals, prop = prop)
end

function catot_yrsdict(region, cause)
	cadeaths = ndeaths([region], [cause])
	totdeaths = ndeaths([region], ["TOT"])
	Dict(
		"numframe" => cadeaths["datasetframe"],
		"denomframe" => totdeaths["datasetframe"],
		"numdim" => cadeaths["dimension"],
		"denomdim" => totdeaths["dimension"],
		"numcause" => cause,
		"denomcause" => "TOT",
		"region" => region
	)
end

function capop_yrsdict(region, cause)
	cadeaths = ndeaths([region], [cause])
	pop = npop([region])
	Dict(
		"numframe" => cadeaths["datasetframe"],
		"denomframe" => pop["datasetframe"],
		"numdim" => cadeaths["dimension"],
		"denomdim" => pop["dimension"],
		"numcause" => cause,
		"denomcause" => "POP",
		"region" => region
	)
end

function catot_mapdict(regvalues, cause, startyear, endyear,
	shapefname = normpath(DATAPATH, "2504", "__pgsql2shp2504_tmp_table.shp"))
	cadeaths = ndeaths(regvalues, [cause],
		yearvalues = yearrange(startyear, endyear))
	totdeaths = ndeaths(regvalues, ["TOT"],
		yearvalues = yearrange(startyear, endyear))
	Dict(
		"numframe" => cadeaths["datasetframe"],
		"denomframe" => totdeaths["datasetframe"],
		"numdim" => cadeaths["dimension"],
		"denomdim" => totdeaths["dimension"],
		"numcause" => cause,
		"denomcause" => "TOT",
		"shapefname" => shapefname
	)
end

function capop_mapdict(regvalues, cause, startyear, endyear,
	shapefname = normpath(DATAPATH, "2504", "__pgsql2shp2504_tmp_table.shp"))
	cadeaths = ndeaths(regvalues, [cause],
		yearvalues = yearrange(startyear, endyear))
	pop = npop(regvalues, yearvalues = yearrange(startyear, endyear))
	Dict(
		"numframe" => cadeaths["datasetframe"],
		"denomframe" => totdeaths["datasetframe"],
		"numdim" => cadeaths["dimension"],
		"denomdim" => totdeaths["dimension"],
		"numcause" => cause,
		"denomcause" => "POP",
		"shapefname" => shapefname
	)
end

end # module