text:
using GeometricIntegrators.CommonFunctions
using GeometricIntegrators.Solutions
using GeometricIntegrators.Solutions: createHDF5
using GeometricIntegrators.TestProblems.HarmonicOscillatorProblem
using HDF5: HDF5File
using Test

nt    = 10
Δt    = .1
ni    = 5

t0    = 0.
x0    = rand(2)
y0    = rand(2)
q0    = rand(1)
p0    = rand(1)
z0    = rand(3)
λ0    = rand(1)
μ0    = rand(2)

t1    = 1.
x1    = rand(2, ni)
y1    = rand(2, ni)
q1    = rand(1, ni)
p1    = rand(1, ni)
z1    = rand(3, ni)
λ1    = rand(1, ni)
μ1    = rand(2, ni)

t2    = t1 + (t1-t0)
x2    = rand(2, ni)
y2    = rand(2, ni)
q2    = rand(1, ni)
p2    = rand(1, ni)
z2    = rand(3, ni)
λ2    = rand(1, ni)
μ2    = rand(2, ni)

tx    = zero(x0)
ty    = zero(y0)
tq    = zero(q0)
tp    = zero(p0)
tz    = zero(z0)
tλ    = zero(λ0)
tμ    = zero(μ0)

xs    = rand(2,nt)
ys    = rand(2,nt)
qs    = rand(1,nt)
ps    = rand(1,nt)
zs    = rand(3,nt)
λs    = rand(1,nt)
μs    = rand(2,nt)

Xs    = rand(2,nt,ni)
Ys    = rand(2,nt,ni)
Qs    = rand(1,nt,ni)
Ps    = rand(1,nt,ni)
Zs    = rand(3,nt,ni)
Λs    = rand(1,nt,ni)
Ms    = rand(2,nt,ni)

ode   = harmonic_oscillator_ode()
dae   = harmonic_oscillator_dae()
pode  = harmonic_oscillator_pode()
pdae  = harmonic_oscillator_pdae()

x100  = rand(ndims(ode),  100)
y100  = rand(ndims(ode),  100)
z100  = rand(ndims(dae),  100)
q100  = rand(ndims(pode), 100)
p100  = rand(ndims(pode), 100)
λ100  = rand(dae.m,  100)
μ100  = rand(pdae.m, 100)

h5file = "test.hdf5"

@testset "$(rpad("ODE Solution",80))" begin
    asol = AtomicSolution(ode)

    # test constructors and general functionality
    sol = Solution(ode, Δt, nt)
    @test typeof(sol) <: SolutionODE

    sol0 = Solution(similar(ode, x0), Δt, nt)
    @test typeof(sol0) <: SolutionODE

    sol1 = Solution(similar(ode, x1), Δt, nt)
    @test typeof(sol1) <: SolutionODE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, x0)
    get_initial_conditions!(sol, tx, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tx == x0
    @test asol.t == t0
    @test asol.q == x0

    δt, δx = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δx == x0

    set_initial_conditions!(sol1, similar(ode, t1, x2))
    for i in 1:ni
        get_initial_conditions!(sol1, tx, i)
        @test tx == x2[:,i]
    end

    # test set/get solution
    sol1 = Solution(similar(ode, x0), Δt, nt)
    sol2 = Solution(similar(ode, x0), Δt, nt)
    for i in 1:nt
        tx .= xs[:,i]
        asol.q .= xs[:,i]
        set_solution!(sol1, tx, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[:,1:nt] == xs
    @test sol2.q[:,1:nt] == xs

    sol1 = Solution(similar(ode, x1), Δt, nt)
    sol2 = Solution(similar(ode, x1), Δt, nt)
    for i in 1:nt
        for k in 1:ni
            tx .= Xs[:,i,k]
            asol.q .= Xs[:,i,k]
            set_solution!(sol1, tx, i, k)
            set_solution!(sol2, asol, i, k)
        end
    end
    @test sol1.q[:,1:nt,:] == Xs
    @test sol2.q[:,1:nt,:] == Xs

    # test general hdf5 functions
    sol1 = Solution(ode, Δt, nt)
    h5 = createHDF5(sol1, h5file)
    @test typeof(h5) == HDF5File
    close(h5)
    @test isfile(h5file)

    sol1 = Solution(ode, Δt, nt)
    h5 = createHDF5(sol1, h5file; overwrite=false)
    @test typeof(h5) == HDF5File
    close(h5)
    @test isfile(h5file)
    rm(h5file)

    sol1 = Solution(ode, Δt, nt)
    write_to_hdf5(sol1, h5file)
    @test isfile(h5file)
    rm(h5file)

    # test hdf5 in- and output
    sol1 = Solution(harmonic_oscillator_ode(x0), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionODE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_ode(x1), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionODE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(ode, Δt, 20, nsave=2)
    @test sol.nt == 10

    sol = Solution(ode, Δt, 20, nsave=2, nwrite=10)
    @test sol.nt == 5

    # test reset
    sol = Solution(ode, Δt, nt)
    reset!(sol)
    @test sol.t[0]   == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt

    sol1 = Solution(ode, Δt, 100, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            set_solution!(sol1, x100[:,(i-1)*10+j], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(ode.q₀, (ndims(ode),1)), x100)
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(ode, Δt, 100, nsave=2, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            set_solution!(sol1, x100[:,(i-1)*10+j], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(ode.q₀, (ndims(ode),1)), x100)[:,1:2:end]
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end

@testset "$(rpad("PODE Solution",80))" begin
    asol = AtomicSolution(pode)

    # test constructors and general functionality
    sol = Solution(pode, Δt, nt)
    @test typeof(sol) <: SolutionPODE

    sol0 = Solution(similar(pode, q0, p0), Δt, nt)
    @test typeof(sol0) <: SolutionPODE

    sol1 = Solution(similar(pode, q1, p1), Δt, nt)
    @test typeof(sol1) <: SolutionPODE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, q0, p0)
    get_initial_conditions!(sol, tq, tp, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tq == q0
    @test tp == p0
    @test asol.t == t0
    @test asol.q == q0
    @test asol.p == p0

    δt, δq, δp = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δq == q0
    @test δp == p0

    set_initial_conditions!(sol1, similar(pode, t1, q2, p2))
    for i in 1:ni
        get_initial_conditions!(sol1, tq, tp, i)
        @test tq == q2[:,i]
        @test tp == p2[:,i]
    end

    # test set/get solution
    sol1 = Solution(similar(pode, q0, p0), Δt, nt)
    sol2 = Solution(similar(pode, q0, p0), Δt, nt)
    for i in 1:nt
        tq .= qs[:,i]
        tp .= ps[:,i]
        asol.q .= qs[:,i]
        asol.p .= ps[:,i]
        set_solution!(sol1, tq, tp, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[:,1:nt] == qs
    @test sol1.p[:,1:nt] == ps
    @test sol2.q[:,1:nt] == qs
    @test sol2.p[:,1:nt] == ps

    sol1 = Solution(similar(pode, q1, p1), Δt, nt)
    sol2 = Solution(similar(pode, q1, p1), Δt, nt)
    for i in 1:nt
        for k in 1:ni
            tq .= Qs[:,i,k]
            tp .= Ps[:,i,k]
            asol.q .= Qs[:,i,k]
            asol.p .= Ps[:,i,k]
            set_solution!(sol1, tq, tp, i, k)
            set_solution!(sol2, asol, i, k)
        end
    end
    @test sol1.q[:,1:nt,:] == Qs
    @test sol1.p[:,1:nt,:] == Ps
    @test sol2.q[:,1:nt,:] == Qs
    @test sol2.p[:,1:nt,:] == Ps

    # test hdf5 in- and output
    sol1 = Solution(harmonic_oscillator_pode(q0, p0), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPODE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_pode(q1, p1), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPODE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(pode, Δt, 20, nsave=2)
    @test sol.nt == 10

    sol = Solution(pode, Δt, 20, nsave=2, nwrite=10)
    @test sol.nt == 5

    # test reset
    sol = Solution(pode, Δt, nt)
    reset!(sol)
    @test sol.t[0]   == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt

    sol1 = Solution(pode, Δt, 100, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, q100[:,oj], p100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionPODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(pode.q₀, (ndims(pode),1)), q100)
    @test sol2.p.d == hcat(reshape(pode.p₀, (ndims(pode),1)), p100)
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(pode, Δt, 100, nsave=2, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, q100[:,oj], p100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionPODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(pode.q₀, (ndims(pode),1)), q100)[:,1:2:end]
    @test sol2.p.d == hcat(reshape(pode.p₀, (ndims(pode),1)), p100)[:,1:2:end]
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end

@testset "$(rpad("DAE Solution",80))" begin
    asol = AtomicSolution(dae)

    # test constructors and general functionality
    sol = Solution(dae, Δt, nt)
    @test typeof(sol) <: SolutionDAE

    sol0 = Solution(similar(dae, z0, λ0), Δt, nt)
    @test typeof(sol0) <: SolutionDAE

    sol1 = Solution(similar(dae, z1, λ1), Δt, nt)
    @test typeof(sol1) <: SolutionDAE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, z0, λ0)
    get_initial_conditions!(sol, tz, tλ, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tz == z0
    @test tλ == λ0
    @test asol.t == t0
    @test asol.q == z0
    @test asol.λ == λ0

    δt, δz, δλ = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δz == z0
    @test δλ == λ0

    set_initial_conditions!(sol1, similar(dae, t1, z2, λ2))
    for i in 1:ni
        get_initial_conditions!(sol1, tz, tλ, i)
        @test tz == z2[:,i]
        @test tλ == λ2[:,i]
    end

    # test set/get solution
    sol1 = Solution(similar(dae, z0, λ0), Δt, nt)
    sol2 = Solution(similar(dae, z0, λ0), Δt, nt)
    for i in 1:nt
        tz .= zs[:,i]
        tλ .= λs[:,i]
        asol.q .= zs[:,i]
        asol.λ .= λs[:,i]
        set_solution!(sol1, tz, tλ, i)
        set_solution!(sol2, asol, i)
    end
    @test sol1.q[:,1:nt] == zs
    @test sol1.λ[:,1:nt] == λs
    @test sol2.q[:,1:nt] == zs
    @test sol2.λ[:,1:nt] == λs

    sol1 = Solution(similar(dae, z1, λ1), Δt, nt)
    sol2 = Solution(similar(dae, z1, λ1), Δt, nt)
    for i in 1:nt
        for k in 1:ni
            tz .= Zs[:,i,k]
            tλ .= Λs[:,i,k]
            asol.q .= Zs[:,i,k]
            asol.λ .= Λs[:,i,k]
            set_solution!(sol1, tz, tλ, i, k)
            set_solution!(sol2, asol, i, k)
        end
    end
    @test sol1.q[:,1:nt,:] == Zs
    @test sol1.λ[:,1:nt,:] == Λs
    @test sol2.q[:,1:nt,:] == Zs
    @test sol2.λ[:,1:nt,:] == Λs

    # test hdf5 in- and output
    sol1 = Solution(harmonic_oscillator_dae(z0, λ0), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionDAE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_dae(z1, λ1), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionDAE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(dae, Δt, 20, nsave=2)
    @test sol.nt == 10

    sol = Solution(dae, Δt, 20, nsave=2, nwrite=10)
    @test sol.nt == 5

    # test reset
    sol = Solution(dae, Δt, nt)
    reset!(sol)
    @test sol.t[0]   == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt

    sol1 = Solution(dae, Δt, 100, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, z100[:,oj], λ100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(dae.q₀, (ndims(dae),1)), z100)
    @test sol2.λ.d == hcat(reshape(dae.λ₀, (dae.m,1)), λ100)
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(dae, Δt, 100, nsave=2, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, z100[:,oj], λ100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(dae.q₀, (ndims(dae),1)), z100)[:,1:2:end]
    @test sol2.λ.d == hcat(reshape(dae.λ₀, (dae.m,1)), λ100)[:,1:2:end]
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end

@testset "$(rpad("PDAE Solution",80))" begin
    asol = AtomicSolution(pdae)

    # test constructors and general functionality
    sol = Solution(pdae, Δt, nt)
    @test typeof(sol) <: SolutionPDAE

    sol0 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    @test typeof(sol0) <: SolutionPDAE

    sol1 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    @test typeof(sol1) <: SolutionPDAE

    @test sol != sol0
    @test sol != sol1

    @test nsave(sol) == 1
    @test ntime(sol) == nt
    @test timesteps(sol) == Δt .* collect(0:nt)

    # test initial conditions
    set_initial_conditions!(sol, t0, x0, y0, μ0)
    get_initial_conditions!(sol, tx, ty, tμ, 1)
    get_initial_conditions!(sol, asol, 1)

    @test tx == x0
    @test ty == y0
    @test tμ == μ0
    @test asol.t == t0
    @test asol.q == x0
    @test asol.p == y0
    @test asol.λ == μ0

    tx .= 0
    ty .= 0
    get_initial_conditions!(sol, tx, ty, 1)
    @test tx == x0
    @test ty == y0

    δt, δx, δy, δμ = get_initial_conditions(sol, 1)
    @test δt == t0
    @test δx == x0
    @test δy == y0
    @test δμ == μ0

    set_initial_conditions!(sol1, similar(pdae, t1, x2, y2, μ2))
    for i in 1:ni
        get_initial_conditions!(sol1, tx, ty, tμ, i)
        @test tx == x2[:,i]
        @test ty == y2[:,i]
        @test tμ == μ2[:,i]
    end

    # test set/get solution
    sol1 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    sol2 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    sol3 = Solution(similar(pdae, x0, y0, μ0), Δt, nt)
    for i in 1:nt
        tx .= xs[:,i]
        ty .= ys[:,i]
        tμ .= μs[:,i]
        asol.q .= xs[:,i]
        asol.p .= ys[:,i]
        asol.λ .= μs[:,i]
        set_solution!(sol1, tx, ty, tμ, i)
        set_solution!(sol2, tx, ty, i)
        set_solution!(sol3, asol, i)
    end
    @test sol1.q[:,1:nt] == xs
    @test sol1.p[:,1:nt] == ys
    @test sol1.λ[:,1:nt] == μs
    @test sol2.q[:,1:nt] == xs
    @test sol2.p[:,1:nt] == ys
    @test sol3.q[:,1:nt] == xs
    @test sol3.p[:,1:nt] == ys
    @test sol3.λ[:,1:nt] == μs

    sol1 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    sol2 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    sol3 = Solution(similar(pdae, x1, y1, μ1), Δt, nt)
    for i in 1:nt
        for k in 1:ni
            tx .= Xs[:,i,k]
            ty .= Ys[:,i,k]
            tμ .= Ms[:,i,k]
            asol.q .= Xs[:,i,k]
            asol.p .= Ys[:,i,k]
            asol.λ .= Ms[:,i,k]
            set_solution!(sol1, tx, ty, tμ, i, k)
            set_solution!(sol2, tx, ty, i, k)
            set_solution!(sol3, asol, i, k)
        end
    end
    @test sol1.q[:,1:nt,:] == Xs
    @test sol1.p[:,1:nt,:] == Ys
    @test sol1.λ[:,1:nt,:] == Ms
    @test sol2.q[:,1:nt,:] == Xs
    @test sol2.p[:,1:nt,:] == Ys
    @test sol3.q[:,1:nt,:] == Xs
    @test sol3.p[:,1:nt,:] == Ys
    @test sol3.λ[:,1:nt,:] == Ms

    # test hdf5 in- and output
    sol1 = Solution(harmonic_oscillator_pdae(x0, y0, μ0), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPDAE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_pdae(x1, y1, μ1), Δt, nt)
    create_hdf5!(sol1, h5file)
    write_to_hdf5(sol1)
    close(sol1)
    @test isfile(h5file)
    sol2 = SSolutionPDAE(h5file)
    @test sol1   != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    # test nsave and nwrite parameters
    sol = Solution(pdae, Δt, 20, nsave=2)
    @test sol.nt == 10

    sol = Solution(pdae, Δt, 20, nsave=2, nwrite=10)
    @test sol.nt == 5

    # test reset
    sol = Solution(pdae, Δt, nt)
    reset!(sol)
    @test sol.t[0]   == t1
    @test sol.t[end] == t2
    @test offset(sol) == nt

    sol1 = Solution(pdae, Δt, 100, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, x100[:,oj], y100[:,oj], μ100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionPDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(pdae.q₀, (ndims(pdae),1)), x100)
    @test sol2.p.d == hcat(reshape(pdae.p₀, (ndims(pdae),1)), y100)
    @test sol2.λ.d == hcat(reshape(pdae.λ₀, (pdae.m,1)), μ100)
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(pdae, Δt, 100, nsave=2, nwrite=10)
    create_hdf5!(sol1, h5file)
    for i in 1:10
        for j in 1:10
            oj = (i-1)*10+j
            set_solution!(sol1, x100[:,oj], y100[:,oj], μ100[:,oj], j)
        end
        write_to_hdf5(sol1)
        reset!(sol1)
    end
    close(sol1)
    sol2 = SSolutionPDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol=1E-14
    @test sol2.q.d == hcat(reshape(pdae.q₀, (ndims(pdae),1)), x100)[:,1:2:end]
    @test sol2.p.d == hcat(reshape(pdae.p₀, (ndims(pdae),1)), y100)[:,1:2:end]
    @test sol2.λ.d == hcat(reshape(pdae.λ₀, (pdae.m,1)), μ100)[:,1:2:end]
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end


text:
using Test

using PlutoHooks

include("./helpers.jl")

#=
We run the tests without Pluto twice:
  1. Without Main.PlutoRunner defined (No Pluto in sight)
  2. With Main.PlutoRunner defined (Pluto is defined but the macro is not run in Pluto)
=#

include("./without_pluto.jl")
include("./with_pluto.jl")
include("./without_pluto.jl")

ext: jl


text:
@testset "Jacobi identity" begin
  @test A × (B × C) ≈ -(B × (C × A) + C × (A × B))
end


text:
import Random
import Dates
import Statistics

function read_noise(path="../noise.dat")
    # Read binary quantum randomness from noise.dat; return it as an Int8 array.
    open(path) do f
        raw_data = read(f)
        as_int8 = reinterpret(Int8, raw_data)
        return(as_int8)
    end
end

function alter_noise(raw_noise)
    # Take an array of Int8 noise, slice off a bit of it, and return it.
    n = 10000 + rand(1:10000) # the pseudorandom number
    new_noise = raw_noise[n:length(raw_noise)-n] # slice from n to the end
    return(new_noise)
end

function create_multiplier(noise::Int8)::Float64
    # Take an Int8 number; convert it into a 'noise_multiplier' for ease of use.
    noise_multiplier = 1 + (noise/256)/4
    return(noise_multiplier)
end

function add_defer_date(workout::String, date::Dates.Date)::String
    workout = "$(workout) @defer($(date) 5am)"
    return(workout)
end

function add_due_date(workout::String, date::Dates.Date)::String
    workout = "$(workout) @due($(date) 5pm)"
    return(workout)
end

function distance_to_sprints(distance_km, sprint_length_m)
    # Takes a distance number, rounds it, squashes it, and converts it to int of
    # n repeats of the specified sprint_length. Pretty ugly/stupid, but works.
    denom = sprint_length_m / 75 # so that longer sprints have fewer repeats
    n_sprints = ((distance_km*1000) / denom) / sprint_length_m
    return(Int(round(n_sprints)))
end

function distance_to_200m(distance)
    # Wrapper for distance_to_sprints, with a cap, for 200m sprints.
    n_sprints = distance_to_sprints(distance, 200)
    return(maximum([minimum([n_sprints, 20]),5]))
end

function distance_to_400m(distance)
    # Wrapper to get n of 400m sprints. Pretty much as above.
    n_sprints = distance_to_sprints(distance, 400)
    return(maximum([minimum([n_sprints, 12]),5]))
end

function distance_to_800m(distance)
    # Wrapper to get 800m sprint n with a cap. Again, as above.
    n_sprints = distance_to_sprints(distance, 800)
    return(maximum([minimum([n_sprints, 8]),4]))
end

function get_run_workout_options(distance::Float64)
    # Takes a workout distance and returns an array of workout options. The
    # proportion of hard vs easy workouts changes as the distance increases.
    hard_w = ["hill_run",
              "sprints_200m",
              "sprints_400m",
              "sprints_800m"]
    easy_w = ["fartlek",
              "tempo"]

    if distance < 15
        options = easy_w
    else
        options = vcat(repeat(easy_w, outer=5),
                       repeat(hard_w, outer=Int(div(distance,10))))
    end

    return(options)
end

function choose_a_run_workout_type(distance::Float64)
    # Take distance, generate the options list, and then choose a workout type.
    options = get_run_workout_options(distance)
    return(options[rand(1:length(options))])
end

function add_estimate(workout::String, distance::Float64)::String
    # Take a workout and its distance and return workout with estimate added.
    workout = "$(workout) @estimate($(10 + Int(round(distance))*6)m)"
    return(workout)
end

function build_run_workout(distance::Float64, kind::String)::String
    # Take distance and workout kind and return full task (with estimate).
    if kind == "sprints_200m"
        workout = "- Go for a run: $(distance_to_200m(distance)) repeats of 200m"
    elseif kind == "sprints_400m"
        workout = "- Go for a run: $(distance_to_400m(distance)) repeats of 400m"
    elseif kind == "sprints_800m"
        workout = "- Go for a run: $(distance_to_800m(distance)) repeats of 800m"
    elseif kind == "tempo"
        workout = "- Go for a tempo run: $(minimum([round((0.7*distance); digits = 1), 20.0])) km"
    elseif kind == "fartlek"
        workout = "- Go for a fartlek run: $(minimum([round((0.55*distance); digits = 1), 20.0])) km"
    elseif kind == "hill_run"
        workout = "- Go for a hill run: $(minimum([round((0.35*distance); digits = 1), 15.0])) km"
    elseif kind == "long_run"
        workout = "- Go for a long run: $(round(distance; digits = 1)) km"
    else
        workout = "- Go for a run"
    end

    return(workout)
end

function build_cross_workout()::String
    options = ["Do a kettlebell workout @estimate(90m)",
               "Do the /r/bodyweightfitness routine @estimate(90m)"]
    workout = "- $(options[rand(1:length(options))])"
    return(workout)
end

function is_run_workout(workout::String)::Bool
    # Takes a workout string and returns whether or not it's a run workout.
    return(occursin("run: ", workout) && occursin("- Go for a ", workout))
end

function is_bodyweight_workout(workout::String)::Bool
    return(occursin("/r/bodyweightfitness routine", workout))
end

function is_kettlebell_workout(workout::String)::Bool
    return(occursin("kettlebell workout", workout))
end

function expand_run_workout(workout::String, pushups::Int64, pullups::Int64)::String
    short_workout = replace(workout, r"@(defer|due).+" => "")

    expanded_workout = "$(workout) @autodone(true) @parallel(false)\n"
    expanded_workout = "$(expanded_workout)    $(short_workout)\n"

    pushups_sets = div(pushups, 25)
    pullups_sets = div(pullups, 10)
    n_sets = maximum([minimum([pushups_sets, pullups_sets]), 1])

    pushups_options = ["wide","diamond","regular"]

    for i in 1:n_sets
        expanded_workout = "$(expanded_workout)    - Do 25 $(pushups_options[rand(1:length(pushups_options))]) pushups @estimate(5m)\n"
        expanded_workout = "$(expanded_workout)    - Do 10 pullups @estimate(5m)\n"
    end

    expanded_workout = "$(expanded_workout)    - Stretch to touch toes for 100 seconds @estimate(2m)\n"
    expanded_workout = "$(expanded_workout)    - Log $(n_sets*25) total pushups to exist.io @estimate(5m)\n"
    expanded_workout = "$(expanded_workout)    - Log $(n_sets*10) total pullups to exist.io @estimate(5m)"

    return(expanded_workout)
end

function get_descending_reps_sets(n_max::Int64)
    all_sets = [1]
    while sum(all_sets) <= n_max
        push!(all_sets, all_sets[length(all_sets)] + 1)
    end

    return(reverse(all_sets))
end

function expand_kettlebell_workout(workout::String, pushups::Int64, pullups::Int64, kettlebells::Int64)::String
    expanded_workout = "$(workout) @autodone(true) @parallel(false)\n"
    sets = [get_descending_reps_sets(x) for x in [Int64(pushups), Int64(pullups), Int64(round(kettlebells/20))]]
    min_set = sets[findmin([length(x) for x in sets])[2]]

    for i in min_set
        expanded_workout = "$(expanded_workout)    - Do $(i*10) kettlebell goblet squats @estimate(5m)\n    - Do $(i) pullups @estimate(5m)\n    - Do $(i*10) kettlebell swings @estimate(5m)\n    - Do $(i) pushups @estimate(5m)\n"
    end

    expanded_workout = "$(expanded_workout)    - Log $(sum(min_set)*10) total kettlebell goblet squats to exist.io @estimate(5m)\n"
    expanded_workout = "$(expanded_workout)    - Log $(sum(min_set)*10) total kettlebell swings to exist.io @estimate(5m)\n"
    expanded_workout = "$(expanded_workout)    - Log $(sum(min_set)) total pullups to exist.io @estimate(5m)\n"
    expanded_workout = "$(expanded_workout)    - Log $(sum(min_set)) total pushups to exist.io @estimate(5m)"

    return(expanded_workout)

end

function expand_bodyweight_workout(workout::String)::String
    expanded_workout = "$(workout) @autodone(true) @parallel(false)\n    - Do 10 resistance band shoulder dislocates @estimate(5m)\n    - Do 10 resistance band chest flies @estimate(5m)\n    - Do 10 resistance band lateral raises @estimate(5m)\n    - Do 10 resistance band front raises @estimate(5m)\n    - Do 10 squat sky reaches @estimate(5m)\n    - Do 10+ wrist warm-up stretches @estimate(5m)\n    - Do 10 cat-camels @estimate(5m)\n    - Do 10 shoulder shrugs @estimate(5m)\n    - Do 60 second front plank @estimate(5m)\n    - Do 60 second left-side plank @estimate(5m)\n    - Do 60 second right-side plank @estimate(5m)\n    - Do 60 second reverse plank @estimate(5m)\n    - Do 60 second hollow hold @estimate(5m)\n    - Do 60 second arch hold @estimate(5m)\n    - Do 5 pull-ups @estimate(5m)\n    - Do 5 dips @estimate(5m)\n    - Do 5 pull-ups @estimate(5m)\n    - Do 5 dips @estimate(5m)\n    - Do 5 pull-ups @estimate(5m)\n    - Do 5 dips @estimate(5m)\n    - Do 30 second L-Sit @estimate(5m)\n    - Do 5 squats @estimate(5m)\n    - Do 30 second L-Sit @estimate(5m)\n    - Do 5 squats @estimate(5m)\n    - Do 30 second L-Sit @estimate(5m)\n    - Do 5 squats @estimate(5m)\n    - Do 5 push-ups @estimate(5m)\n    - Do 5 rows @estimate(5m)\n    - Do 5 push-ups @estimate(5m)\n    - Do 5 rows @estimate(5m)\n    - Do 5 push-ups @estimate(5m)\n    - Do 5 rows @estimate(5m)\n    - Do hip-flexors mobility stretches @estimate(5m)\n    - Hold pigeon-pose for 60 seconds each side @estimate(5m)\n    - Stretch to touch toes for 100 seconds @estimate(5m)\n    - Log 15 total pull-ups to exist.io @estimate(5m)\n    - Log 15 total push-ups to exist.io @estimate(5m)"
    return(expanded_workout)
end

function expand_cross_workout(workout::String, pushups::Int64, pullups::Int64, kettlebells::Int64)::String
    if is_kettlebell_workout(workout)
        expanded_workout = expand_kettlebell_workout(workout, pushups, pullups, kettlebells)
    elseif is_bodyweight_workout(workout)
        expanded_workout = expand_bodyweight_workout(workout)
    else
        expanded_workout = workout
    end
    return(expanded_workout)
end

function detect_distance(workout::String)::Float64
    # Take a raw workout string (possibly expanded) and return estimate distance
    headline = string(match(r"^.+", workout).match)
    estimate = string(match(r"@estimate\([0-9]+m\)", headline).match)

    if occursin("repeats of", headline)
        return(parse(Float64, string(match(r"[0-9]+", estimate).match)) / 6)
    elseif is_run_workout(headline)
        distance = match(r"[0-9]+\.[0-9]+ km", headline).match
        return(parse(Float64, distance[1:end-3]))
    else
        return(parse(Float64, string(match(r"[0-9]+", estimate).match)) / 12)
    end
end

function generate_workout(raw_distance::Float64,
                          noise_multiplier::Float64,
                          n_runs_per_week::Int64,
                          push_up_max::Int64,
                          pull_up_max::Int64,
                          kettlebell_max::Int64,
                          current_day::Dates.Date,
                          expand_task = true,
                          bias = 0.01,
                          rest_day = Dates.Sunday,
                          must_run = false)

    distance = raw_distance * (noise_multiplier + bias)

    if Dates.dayofweek(current_day + Dates.Day(1)) == rest_day
        workout = build_run_workout(distance, "long_run")
        workout = add_estimate(workout, distance)
    elseif Dates.dayofweek(current_day) == rest_day
        workout = "- Do a session of yoga with Adriene @estimate(90m)"
    else
        if rand(1:5) <= n_runs_per_week-1 || must_run
            workout = build_run_workout(distance,
                                        choose_a_run_workout_type(distance))
            workout = add_estimate(workout, distance)
        else
            workout = build_cross_workout()
        end
    end

    workout = add_defer_date(workout, current_day)
    workout = add_due_date(workout, current_day)

    if expand_task
        if is_run_workout(workout)
            workout = expand_run_workout(workout,
                                         push_up_max,
                                         pull_up_max)
        else
            workout = expand_cross_workout(workout,
                                           push_up_max,
                                           pull_up_max,
                                           kettlebell_max)
        end
    end

    return(workout)
end

function generate_workout_sequence(last_run_distance::Float64,
                                   goal_run_distance::Float64,
                                   raw_noise_array,
                                   n_runs_per_week::Int64,
                                   current_pushups::Int64,
                                   pushups_max::Int64,
                                   current_pullups::Int64,
                                   pullups_max::Int64,
                                   current_kettlebells::Int64,
                                   kettlebells_max::Int64,
                                   program_start_date::Dates.Date,
                                   expand_task = true,
                                   bias = 0.05,
                                   rest_day = Dates.Sunday)
    pushups = current_pushups
    pullups = current_pullups
    kettlebells = current_kettlebells
    workouts = []
    days_since_last_run = 0
    workout_distances = [last_run_distance]
    workout_date = program_start_date

    while workout_distances[length(workout_distances)] < goal_run_distance

        target_distance = Statistics.mean([5 + length(workouts)/21,
                                           5 + length(workouts)/14,
                                           Statistics.mean(workout_distances),
                                           Statisitics.median(workout_distances)])

        workout = generate_workout(target_distance,
                                   create_multiplier(pop!(raw_noise_array)),
                                   n_runs_per_week,
                                   pushups,
                                   pullups,
                                   kettlebells,
                                   workout_date,
                                   expand_task,
                                   bias,
                                   rest_day,
                                   days_since_last_run > 5)

        if pushups < pushups_max
            pushups += rand(0:1)
        end

        if pullups < pullups_max && Dates.dayofweek(workout_date) == rest_day
            pullups += 1
        end

        if kettlebells < kettlebells_max
            kettlebells += 5
        end

        push!(workouts, workout)

        if is_run_workout(workout)
            append!(workout_distances, detect_distance(workout))
            days_since_last_run = 0
        else
            days_since_last_run += 1
        end

        workout_date = workout_date + Dates.Day(1)
    end
    return(workouts)
end


text:
## runtests.jl  Unit tests for ROC
## (c) 2014--2015 David A. van Leeuwen
##
## Licensed under the MIT software license, see LICENSE.md

using LinearAlgebra
using ROCAnalysis
using Random

include("extremes.jl")
include("big.jl")
include("dcf.jl")
include("delong.jl")


text:
import Base.<=

export AbstractPolygon,
       tohrep,
       tovrep

"""
    AbstractPolygon{N<:Real} <: AbstractPolytope{N}

Abstract type for polygons (i.e., 2D polytopes).

### Notes

Every concrete `AbstractPolygon` must define the following functions:
- `tovrep(::AbstractPolygon{N})::VPolygon{N}`         -- transform into
    V-representation
- `tohrep(::AbstractPolygon{N})::S where {S<:AbstractHPolygon{N}}` -- transform
    into H-representation

```jldoctest
julia> subtypes(AbstractPolygon)
2-element Array{Any,1}:
 AbstractHPolygon
 VPolygon
```
"""
abstract type AbstractPolygon{N<:Real} <: AbstractPolytope{N} end

# --- LazySet interface functions ---

"""
    dim(P::AbstractPolygon)::Int

Return the ambient dimension of a polygon.

### Input

- `P` -- polygon

### Output

The ambient dimension of the polygon, which is 2.
"""
@inline function dim(P::AbstractPolygon)::Int
    return 2
end

"""
    jump2pi(x::N)::N where {N<:AbstractFloat}

Return ``x + 2π`` if ``x`` is negative, otherwise return ``x``.

### Input

- `x` -- real scalar

### Output

``x + 2π`` if ``x`` is negative, ``x`` otherwise.

### Examples

```jldoctest
julia> using LazySets: jump2pi

julia> jump2pi(0.0)
0.0

julia> jump2pi(-0.5)
5.783185307179586

julia> jump2pi(0.5)
0.5
```
"""
@inline function jump2pi(x::N)::N where {N<:AbstractFloat}
    x < zero(N) ? 2 * pi + x : x
end

"""
    quadrant(w::AbstractVector{N})::Int where {N<:Real}

Compute the quadrant where the direction `w` belongs.

### Input

- `w` --  direction

### Output

An integer from 0 to 3, with the following convention:

```
     ^
   1 | 0
  ---+-->
   2 | 3
```

### Algorithm

The idea is to encode the following logic function:
``11 ↦ 0, 01 ↦ 1, 00 ↦ 2, 10 ↦ 3``, according to the convention of above.

This function is inspired from AGPX's answer in:
[Sort points in clockwise order?](https://stackoverflow.com/a/46635372)
"""
@inline function quadrant(w::AbstractVector{N})::Int where {N<:Real}
    dwx = w[1] >= zero(N) ? 1 : 0
    dwy = w[2] >= zero(N) ? 1 : 0
    return (1 - dwx) + (1 - dwy) + ((dwx & (1 - dwy)) << 1)
end

"""
    is_right_turn(u::AbstractVector{N},
                  v::AbstractVector{N})::Bool where {N<:Real}

Determine if the acute angle defined by two 2D vectors is a right turn (< 180°
counter-clockwise).

### Input

- `u` -- first 2D direction
- `v` -- second 2D direction

### Output

`true` iff the two vectors constitute a right turn.
"""
@inline function is_right_turn(u::AbstractVector{N},
                               v::AbstractVector{N})::Bool where {N<:Real}
    return u[1] * v[2] - v[1] * u[2] >= zero(N)
end

"""
    <=(u::AbstractVector{N}, v::AbstractVector{N})::Bool where {N<:Real}

Compare two 2D vectors by their direction.

### Input

- `u` -- first 2D direction
- `v` -- second 2D direction

### Output

`true` iff ``\\arg(u) [2π] ≤ \\arg(v) [2π]``.

### Notes

The argument is measured in counter-clockwise fashion, with the 0 being the
direction (1, 0).

### Algorithm

The implementation checks the quadrant of each direction, and compares
directions using the right-hand rule (see [`is_right_turn`](@ref)).
In particular, this method does not use the arctangent.
"""
function <=(u::AbstractVector{N}, v::AbstractVector{N})::Bool where {N<:Real}
    qu, qv = quadrant(u), quadrant(v)
    if qu == qv
        # same quadrant, check right-turn with center 0
        return is_right_turn(u, v)
    end
    # different quadrant
    return qu < qv
end

"""
    <=(u::AbstractVector{N}, v::AbstractVector{N})::Bool where {N<:AbstractFloat}

Compares two 2D vectors by their direction.

### Input

- `u` --  first 2D direction
- `v` --  second 2D direction

### Output

True iff ``\\arg(u) [2π] ≤ \\arg(v) [2π]``

### Notes

The argument is measured in counter-clockwise fashion, with the 0 being the
direction (1, 0).

### Algorithm

The implementation uses the arctangent function with sign, `atan`, which for two
arguments implements the
[`atan2` function](https://en.wikipedia.org/wiki/Atan2).
"""
function <=(u::AbstractVector{N},
            v::AbstractVector{N})::Bool where {N<:AbstractFloat}
    return jump2pi(atan(u[2], u[1])) <= jump2pi(atan(v[2], v[1]))
end


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 3f3dbc34-3593-11eb-2250-13fc496e403e
begin
	import Pkg; Pkg.add.(["LaTeXStrings", "PlutoUI", "Distributions", "ColorSchemes", "Measures"])
	using Plots, LaTeXStrings, PlutoUI, Measures
	import Distributions: Beta, pdf
	import ColorSchemes.Paired_12;
	colors = Paired_12[[1,2,7,8,3,4,5,6,9,10]]
	Plots.default()
	gr(legend = nothing, grid = false, color = colors[2], lw = 10, legendfontsize=20,
    		xtickfontsize=20, ytickfontsize=20, xguidefontsize=20, yguidefontsize=20,
		titlefontsize = 20)
	nothing
end

# ╔═╡ 2ea253b6-358a-11eb-10c9-6d257cbca066
md"### Bayesian analysis of the Bernoulli model"

# ╔═╡ 5d6396b0-358a-11eb-2995-eb061a2476aa
md" **Bernoulli model**:

$$x_1,\ldots,x_n |\theta \overset{\mathrm{iid}}{\sim} \mathrm{Bernoulli}(\theta)$$"

# ╔═╡ b7b9303e-358a-11eb-067d-fd0020d579a3
md"**Beta prior**:

$$\theta \sim \mathrm{Beta}(\alpha,\beta)$$"

# ╔═╡ da5f2ce0-358a-11eb-2522-81165870626e
md" **Beta posterior**:

$$\theta | x_1,\ldots,x_n \sim \mathrm{Beta}(\alpha+s,\beta+f),$$

where $s=\sum_{i=1}^n x_i$ is the number of successes and $f = n-s$ is the number of failures."

# ╔═╡ d39a7b72-358b-11eb-1d11-dfdcf238a7f4
md"""
**Input data and prior hyperparameters**:

n =  $(@bind n Slider(1:100; default=10, show_value=true))

s/n =  $(@bind sFrac Slider(0:0.01:1; default=0.5, show_value=true))

\$\alpha = \$ $(@bind α Slider(eps():100; default=3, show_value=true))

\$\beta = \$ $(@bind β Slider(eps():100; default=3, show_value=true))

"""

# ╔═╡ 698e86b0-3590-11eb-185b-df426ce40664
begin
	# Function to compute the Beta posterior and moments
	function BernBeta(s, f, α, β)
		priorMean = α/(α+β)
		priorStd = √(α*β/((α+β)^2*(α+β+1)))
		αPost = α + s
		βPost = β + f
		postMean = αPost/(αPost+βPost)
		postStd = √(αPost*βPost/((αPost+βPost)^2*(αPost+βPost+1)))
		post = Beta(αPost, βPost)
		return post, αPost, βPost, priorMean, priorStd, postMean, postStd
	end
	s = sFrac*n
	post, αPost, βPost, priorMean, priorStd, postMean, postStd = BernBeta(s, n-s, α, β)
end;

# ╔═╡ ca4114a4-3591-11eb-1777-4593b8760a68
begin
	f = n - s
	θGrid = range(eps(),1-eps(), length = 1000)
	binWidth = θGrid[2]-θGrid[1]
	logℓ = s*log.(θGrid) .+ f*log.(1 .- θGrid)
	like = exp.(logℓ .- maximum(logℓ))
	normLike = like/(sum(like)*binWidth)
	plot(size = (1200,800), θGrid, pdf.(Beta(α,β),θGrid), lw = 10, color = colors[6],
		xlabel = L"\theta", label = "Prior", yaxis = false, margin = 5mm)
	plot!(θGrid, normLike, lw = 10, color = colors[2], label = "Likelihood")
	plot!(θGrid, pdf.(Beta(αPost,βPost),θGrid), lw = 10, color = colors[4], label = "Posterior")
end

# ╔═╡ 7cdb3f2c-358d-11eb-03e6-df23873fb7b4
md"""
**Prior mean and standard deviation**:

\$\mathbb{E}(\theta)=\$ $(priorMean)

\$\mathbb{S}(\theta)=\$ $(priorStd)
"""

# ╔═╡ 25b33df2-3598-11eb-2364-a540b649f80d
md"""
**Posterior mean and standard deviation**:

\$\mathbb{E}(\theta|\mathbf{x})=\$  $(postMean)

\$\mathbb{S}(\theta|\mathbf{x})=\$ $(postStd)
"""

# ╔═╡ 3d25b63c-359c-11eb-2aa0-c558e2eb0a6b

# ╔═╡ 1f03b924-359c-11eb-0e7a-41c84b5dc313
md"""
---
Notebook for the book: Villani, M. (2020). [Bayesian Learning]()."""

# ╔═╡ Cell order:
# ╟─2ea253b6-358a-11eb-10c9-6d257cbca066
# ╟─3f3dbc34-3593-11eb-2250-13fc496e403e
# ╟─5d6396b0-358a-11eb-2995-eb061a2476aa
# ╟─b7b9303e-358a-11eb-067d-fd0020d579a3
# ╟─da5f2ce0-358a-11eb-2522-81165870626e
# ╠═d39a7b72-358b-11eb-1d11-dfdcf238a7f4
# ╠═698e86b0-3590-11eb-185b-df426ce40664
# ╟─ca4114a4-3591-11eb-1777-4593b8760a68
# ╟─7cdb3f2c-358d-11eb-03e6-df23873fb7b4
# ╟─25b33df2-3598-11eb-2364-a540b649f80d
# ╟─3d25b63c-359c-11eb-2aa0-c558e2eb0a6b
# ╠═1f03b924-359c-11eb-0e7a-41c84b5dc313

max_stars_repo_name: mattiasvillani/BLnotebooks

max_issues_repo_name: mattiasvillani/BLnotebooks

max_forks_repo_name: mattiasvillani/BLnotebooks


text:
# Advection on the Sphere test case.
#
# This is a case similar to the one presented in:
# @article{WILLIAMSON1992211,
#   title = {A standard test set for numerical approximations to the
#            shallow water equations in spherical geometry},
#   journal = {Journal of Computational Physics},
#   volume = {102},
#   number = {1},
#   pages = {211 - 224},
#   year = {1992},
#   doi = {10.1016/S0021-9991(05)80016-6},
#   url = {https://doi.org/10.1016/S0021-9991(05)80016-6},
# }
#
# This version runs the advection on the sphere as a stand alone test (no
# dependence on CLIMA moist thermodynamics)
#--------------------------------#
#--------------------------------#
# Can be run with:
# Integration Testing: JULIA_CLIMA_INTEGRATION_TESTING=true mpirun -n 1 julia --project=@. advection_sphere_lsrk.jl
# No Integration Testing: JULIA_CLIMA_INTEGRATION_TESTING=false mpirun -n 2 julia --project=@. advection_sphere_lsrk.jl
#--------------------------------#
#--------------------------------#
# Can be run on the GPU with:
# Integration Testing: JULIA_CLIMA_INTEGRATION_TESTING=true mpirun -n 1 julia --project=/home/fxgiraldo/CLIMA/env/gpu advection_sphere_lsrk.jl
# No Integration Testing: JULIA_CLIMA_INTEGRATION_TESTING=false mpirun -n 1 julia --project=/home/fxgiraldo/CLIMA/env/gpu advection_sphere_lsrk.jl
#--------------------------------#
#--------------------------------#

using MPI
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGBalanceLawDiscretizations
using CLIMA.DGBalanceLawDiscretizations.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.VTK
using LinearAlgebra
using StaticArrays
using Logging, Printf, Dates

const uid, vid, wid = 1:3
const radians = true
const γ_exact = 7 // 5

if !@isdefined integration_testing
  const integration_testing =
  parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

@inline function preflux(Q)
  γ::eltype(Q) = γ_exact
  @inbounds ρ = Q[1]
  ρinv = 1 / ρ
  u, v, w, P = 0, 0, 0, 0
  (P, u, v, w, ρinv)
end

#{{{ advectionflux
@inline function advectionflux!(F,Q,QV,aux,t)
  P,u,v,w,ρinv = preflux(Q)
  @inbounds begin
    u,v,w = aux[uid], aux[vid], aux[wid]
    ρ=Q[1]
    F[1,1], F[2,1], F[3,1]=u*ρ, v*ρ, w*ρ
  end
end
#}}} advectionflux

#{{{ wavespeed
@inline function wavespeed(n, Q, vel, t...)
  # abs(dot(vel,n)) # FIXME: Nice to do this, but breaks with cassette
  abs(vel[1] * n[1] + vel[2] * n[2] + vel[3] * n[3])
end
#}}} wavespeed

#{{{ Cartesian->Spherical
@inline function cartesian_to_spherical(x,y,z,radians)
  # Conversion Constants
  if radians
    c=1.0
  else
    c=180/π
  end
  λ_max=2π*c
  λ_min=0*c
  ϕ_max=+0.5*π*c
  ϕ_min=-0.5*π*c

  # Conversion functions
  r = hypot(x,y,z)
  λ=atan(y,x)*c
  ϕ=asin(z/r)*c
  return (r, λ, ϕ)
end
#}}} Cartesian->Spherical

#{{{ velocity initial condition
@inline function velocity_init!(vel, x, y, z)
  @inbounds begin
    FT = eltype(vel)
    (r, λ, ϕ) = cartesian_to_spherical(x,y,z,radians)
    # w = 2 * FT(π) * cos(ϕ) # Case 1 -> shear flow
    w = 2 * FT(π) * cos(ϕ) * r # Case 2 -> solid body flow
    uλ, uϕ = w, 0
    vel[uid] = -uλ*sin(λ) - uϕ*cos(λ)*sin(ϕ)
    vel[vid] = +uλ*cos(λ) - uϕ*sin(λ)*sin(ϕ)
    vel[wid] = +uϕ*cos(ϕ)
  end
end
#}}} velocity initial condition

# initial condition
function advection_sphere!(Q, t, x, y, z, vel)
  FT = eltype(Q)
  rc=1.5
  (r, λ, ϕ) = cartesian_to_spherical(x,y,z,radians)
  ρ = exp(-((3λ)^2 + (3ϕ)^2))

  @inbounds Q[1] = ρ
end

#{{{ Main
function main(mpicomm, FT, topl, N, timeend, ArrayType, dt, ti_method)
  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                          meshwarp = Topologies.cubedshellwarp)

  #{{{ numerical fluxex
  numflux!(x...) = NumericalFluxes.rusanov!(x..., advectionflux!,
                                            wavespeed)

  # zero flux at the boundaries: not entirely accurate but good enough
  function numbcflux!(F, _...)
    for s = 1:length(F)
      F[s] = 0
    end
  end

  # Define Spatial Discretization
  spacedisc=DGBalanceLaw(grid = grid,
                         length_state_vector = 1,
                         flux! = advectionflux!,
                         numerical_flux! = numflux!,
                         numerical_boundary_flux! = numbcflux!,
                         auxiliary_state_length = 3,
                         auxiliary_state_initialization! = velocity_init!)

  # Initial Condition
  initialcondition(Q, x...) = advection_sphere!(Q, 0, x...)
  Q = MPIStateArray(spacedisc, initialcondition)

  # Store Initial Condition as Exact Solution
  Qe = copy(Q)

  # Define Time-Integration Method
  if ti_method == "LSRK"
    TimeIntegrator = LSRK54CarpenterKennedy(spacedisc, Q; dt = dt, t0 = 0)
  elseif ti_method == "SSP33"
    TimeIntegrator = SSPRK33ShuOsher(spacedisc, Q; dt = dt, t0 = 0)
  elseif ti_method == "SSP34"
    TimeIntegrator = SSPRK34SpiteriRuuth(spacedisc, Q; dt = dt, t0 = 0)
  end

  #------------Set Callback Info--------------------------------#
  # Set up the information callback
  starttime = Ref(now())
  cbinfo = GenericCallbacks.EveryXWallTimeSeconds(10, mpicomm) do (s=false)
    if s
      starttime[] = now()
    else
      @info @sprintf("""Update
                     simtime = %.16e
                     runtime = %s
                     Δmass   = %.16e""",
                     ODESolvers.gettime(TimeIntegrator),
                     Dates.format(convert(Dates.DateTime,
                                          Dates.now()-starttime[]),
                                  Dates.dateformat"HH:MM:SS"),
                     abs(weightedsum(Q) - weightedsum(Qe)) / weightedsum(Qe))
    end
    nothing
  end

  # Set up the information callback
  cbmass = GenericCallbacks.EveryXSimulationSteps(1000) do
    @info @sprintf """Update
    Δmass   = %.16e""" abs(weightedsum(Q) - weightedsum(Qe)) / weightedsum(Qe)
  end

  step = [0]
  mkpath("vtk")
  cbvtk = GenericCallbacks.EveryXSimulationSteps(1000) do (init=false)
    outprefix = @sprintf("vtk/advection_sphere_%dD_mpirank%04d_step%04d",
                         3, MPI.Comm_rank(mpicomm), step[1])
    @debug "doing VTK output" outprefix
    writevtk(outprefix, Q, spacedisc, ("ρ", ))
    step[1] += 1
    nothing
  end
  #------------Set Callback Info--------------------------------#

  # Perform Time-Integration
  solve!(Q, TimeIntegrator; timeend=timeend, callbacks=(cbinfo, cbmass, cbvtk))

  # Print some end of the simulation information
  error = euclidean_distance(Q, Qe) / norm(Qe)
  Δmass = abs(weightedsum(Q) - weightedsum(Qe)) / weightedsum(Qe)
  @info @sprintf """Finished
  error = %.16e
  Δmass = %.16e
  """ error Δmass

  # return diagnostics
  return (error, Δmass)
end
#}}} Main

#{{{ Run Script
function run(mpicomm, Nhorizontal, Nvertical, N, timeend, FT, dt, ti_method,
             ArrayType)
  Rrange=range(FT(1); length=Nvertical+1, stop=2)
  topl = StackedCubedSphereTopology(mpicomm,Nhorizontal,Rrange; boundary=(1,1))
  (error, Δmass) = main(mpicomm, FT, topl, N, timeend, ArrayType, dt,
                        ti_method)
end
#}}} Run Script

#{{{ Run Program
using Test
let
  CLIMA.init()
  ArrayTypes = (CLIMA.array_type(),)

  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
  ll == "WARN"  ? Logging.Warn  :
  ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  # Perform Integration Testing for three different grid resolutions
  ti_method = "LSRK" #LSRK or SSP
  timeend = 1
  numelem = (2, 2) #(Nhorizontal,Nvertical)
  N = 4

  expected_error = Array{Float64}(undef, 3) # h-refinement levels lvl
  expected_error[1] = 1.7501966878661540e-01 # Ne=2
  expected_error[2] = 9.8077688360037538e-03 # Ne=4
  expected_error[3] = 1.7739284567936140e-04 # Ne=8
  expected_mass = Array{Float64}(undef, 3) # h-refinement levels lvl
  expected_mass[1] = 2.8041360758393849e-16 # Ne=2
  expected_mass[2] = 1.8219438767875646e-15 # Ne=4
  expected_mass[3] = 6.3067022934564928e-15 # Ne=8
  lvls = integration_testing ? length(expected_error) : 1

  @testset "$(@__FILE__)" for ArrayType in ArrayTypes
    dt=1e-2*5 # stable dt for N=4 and Ne=5
    for FT in (Float64,) # Float32)
      err = zeros(FT, lvls)
      mass= zeros(FT, lvls)
      for l = 1:lvls
        Nhorizontal = 2^(l-1) * numelem[1]
        Nvertical   = 2^(l-1) * numelem[2]
        dt=dt/Nhorizontal
        nsteps = ceil(Int64, timeend / dt)
        dt = timeend / nsteps
        @info @sprintf """Run Configuration
        Nhorizontal = %d
        Nvertical   = %d
        N           = %d
        dt          = %.16e
        nsteps       = %d
        """ Nhorizontal Nvertical N dt nsteps
        @info (ArrayType, FT)
        (err[l], mass[l]) = run(mpicomm, Nhorizontal, Nvertical, N, timeend,
                                FT, dt, ti_method, ArrayType)
        @test err[l]  ≈ FT(expected_error[l])
        #                @test mass[l] ≈ FT(expected_mass[l])
      end
      @info begin
        msg = ""
        for l = 1:lvls-1
          rate = log2(err[l]) - log2(err[l+1])
          msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
        end
        msg
      end
    end
  end

  #=
  # This snippet of code allows one to run just one instance/configuration.
  # Before running this, Comment the Integration Testing block above
  FT = Float64
  N=4
  ArrayType = Array
  dt=1e-2*5 # stable dt for N=4 and Ne=5
  ti_method = "SSP34" #LSRK or SSP
  timeend=1.0
  Nhorizontal = 2 # number of horizontal elements per face of cubed-sphere grid
  Nvertical = 2 # number of horizontal elements per face of cubed-sphere grid
  dt=dt/Nhorizontal
  nsteps = ceil(Int64, timeend / dt)
  dt = timeend / nsteps
  (error, Δmass) = run(mpicomm, Nhorizontal, Nvertical, N, timeend, FT, dt,
                       ti_method, ArrayType)
  =#

end # Test

# nothing

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
#################### ChainProgress Meter ####################

#################### Types and Constructors ####################

struct ChainProgressFrame
  verbose::Bool

  function ChainProgressFrame(title::AbstractString, verbose::Bool)
    verbose && print(title * "...\n\n")
    new(verbose)
  end
end

mutable struct ChainProgress
  frame::ChainProgressFrame
  chain::Int
  iters::Int
  counter::Int
  runin::Int
  threshold::Float64
  t0::Float64

  function ChainProgress(frame::ChainProgressFrame, chain::Integer,
                         iters::Integer)
    new(frame, chain, iters, 0, max(1, min(10, round(Int, 0.01 * iters))),
        0.0, time())
  end
end

#################### Base Methods ####################

function reset!(p::ChainProgress)
  p.counter = 0
  p.threshold = 0.0
  p.t0 = time()
  p
end

function next!(p::ChainProgress)
  p.counter += 1
  if p.counter / p.iters >= p.threshold && p.counter >= p.runin
    p.threshold += 0.10
    p.frame.verbose && print(stdout, p)
  end
  p
end

function Base.print(io::IO, p::ChainProgress)
  elapsed = time() - p.t0
  remaining = elapsed * (p.iters / p.counter - 1.0)
  str = @sprintf("Chain %u: %3u%% [%s of %s remaining]%s",
                 p.chain,
                 100.0 * p.counter / p.iters,
                 strfsec(remaining),
                 strfsec(elapsed + remaining),
                 "\n"^(p.counter == p.iters))
  println(io, str)
end

function strfsec(sec::Real)
  hr, sec = divrem(sec, 3600)
  min, sec = divrem(sec, 60)
  @sprintf "%u:%02u:%02u" hr min sec
end


text:
using Random, Distributions, KernelDensity, Plots, LaTeXStrings; pyplot()
Random.seed!(1)

function pval(mu0,mu,sig,n)
    sample = rand(Normal(mu,sig),n)
    xBar   = mean(sample)
    s      = std(sample)
    tStat  = (xBar-mu0) / (s/sqrt(n))
    ccdf(TDist(n-1), tStat)
end

mu0, mu1A, mu1B  = 20, 22, 24
sig, n, N = 7, 5, 10^6

pValsH0 = [pval(mu0,mu0,sig,n) for _ in 1:N]
pValsH1A = [pval(mu0,mu1A,sig,n) for _ in 1:N]
pValsH1B = [pval(mu0,mu1B,sig,n) for _ in 1:N]

alpha = 0.05
estPwr(pVals) = sum(pVals .< alpha)/N

println("Power under H0:  ", estPwr(pValsH0))
println("Power under H1A: ", estPwr(pValsH1A))
println("Power under H1B (mu's farther apart): ", estPwr(pValsH1B))

stephist(pValsH0, bins=100,
	normed=true, c=:blue, label="Under H0")
stephist!(pValsH1A, bins=100,
	normed=true, c=:red, label="Under H1A")
stephist!(pValsH1B, bins=100,
	normed=true, c=:green, label="Under H1B",
    xlims=(0,1), ylims=(0,6), xlabel = "p-value", ylabel = "Density")

plot!([alpha,alpha],[0,6],
	c=:black, ls=:dash, label=L"\alpha")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
mutable struct RvrpComputedData
    # instance_id::String
    # pickup_ids::Vector{String}
    # delivery_ids::Vector{String}
    # depot_ids::Vector{String}
    # recharging_ids::Vector{String}
    location_id_2_index::Dict{String, Int}
    location_group_id_2_index::Dict{String, Int}
    product_specification_class_id_2_index::Dict{String, Int}
    request_id_2_index::Dict{String, Int}
    vehicle_category_id_2_index::Dict{String, Int}
    vehicle_set_id_2_index::Dict{String, Int}
    capacity_id_2_index::Dict{String, Int}
    property_id_2_index::Dict{String, Int}
    travel_specification_id_2_index::Dict{String, Int}
    features::BitSet
    uses_default_vehicle_category::Bool
end


text:
#Test code
#-------------------------------------------------------------------------------

using InspectDR
using Base.Test

#Simply test for exceptions on sample plots for now:
include("../sample/runsamples.jl")

:Test_Complete


text:
using ITensors
include("2d_classical_ising.jl")
include("ctmrg.jl")

function main()
  # Unit cell size
  ny,nx = 2,2

  # Make the site indices
  sh,sv = site_inds(ny,nx,2)

  # Make the Ising partition function
  β = 1.1*βc()
  T = ising_partition(sh,sv,β)

  # Make the initial environment
  (Clu,Cru,Cld,Crd),(Al,Ar,Au,Ad) = ctmrg_environment((sh,sv))

  # Check the initialize environment
  check_environment(T,(Clu,Cru,Cld,Crd),(Al,Ar,Au,Ad))

  # Run ctmrg
  @show κave = ctmrg(T,(Clu,Cru,Cld,Crd),(Al,Ar,Au,Ad))

  return
end


text:
# Author: Sean Wu
# First Written: December 15

"""
    getdata_ppdstepsize([ppd_stepsizestotest, seed]) -> DataFrame

Run controlled ppd experiments at different ppd step sizes (input in an array).

Other simulation controls include:
    set_PPD_ON() -> IMPORTANT: Make sure PPD_ON is true before running.
    set_MAX_STEPS_PER_WALK()
    set_NUMBER_OF_WALKS()
Lab data shows no ppd experiments (ppd stepsize = 0.1) sees double the yield of ppd.
Target when the sensor receives half as many peroxides with smaller ppd stepsizes.
"""
function getdata_ppdstepsize(ppd_stepsizestotest::Array{Float64, 1}=[0.001], seed::Int=randseed())
    ppd_print(seed)
    ss_arr, sensor_arr, escape_arr, unresv_arr = ([], [], [], [])
    set_PPD_ON(true)
    for ss in ppd_stepsizestotest
        set_PPD_STEP_SIZE(ss)
        append!(ss_arr, ss)

        data = Dict{String,Integer}()
        data["sensor"]  = 0
        data["escape"] = 0
        data["particles unresolved"] = 0
        data = run_sim!(data, seed)
        append!(sensor_arr, data["sensor"])
        append!(escape_arr, data["escape"])
        append!(unresv_arr, data["particles unresolved"])
    end
    df = DataFrame(ppd_step_size = ss_arr,
        sensor_yield = sensor_arr,
        escaped = escape_arr,
        unresolved = unresv_arr)
    current_df(df)
    current_seed(seed)
    current_path("out/ppddata/")
    return df
end

"Helper function to print parameters, not exported in PPD module"
function ppd_print(seed::Int)
    println("############################")
    println("   Compare PPD Step Sizes   ")
    println("############################")
    println("_________Parameters_________")
    println("# of trials:\t\t", string(NUMBER_OF_WALKS))
    println("# of steps (max):\t", string(MAX_STEPS_PER_WALK))
    println("Step size, water:\t", string(WATER_STEP_SIZE))
    println("Random seed:\t\t$seed")
    nothing
end

"""
    ppd_sim([seed])

Run one simulation of the ppd experiment with one sensor in a well.

Simulation controls:
    set_PPD_ON() -> if PPD_ON === false, then PPD_STEP_SIZE = WATER_STEP_SIZE
    set_MAX_STEPS_PER_WALK()
    set_NUMBER_OF_WALKS()
    set_PPD_STEP_SIZE()
"""
function ppd_sim(seed::Int=randseed())
    println("############################")
    println("       PPD Experiment       ")
    println("############################")
    println("_________Parameters_________")
    println("PPD on sensor:\t\t", string(PPD_ON))
    println("# of trials:\t\t", string(NUMBER_OF_WALKS))
    println("# of steps (max):\t", string(MAX_STEPS_PER_WALK))
    println("Step size, water:\t", string(WATER_STEP_SIZE))
    if PPD_ON
        println("Step size, PPD:\t\t", string(PPD_STEP_SIZE))
    end
    println("Random seed:\t\t$seed")

    data = Dict{String,Integer}()
    data["sensor"] = 0
    data["escape"] = 0
    data["particles unresolved"] = 0
    output = run_sim!(data, seed)

    println("_________Results____________")
    println("# in sensor:\t\t", string(data["sensor"]))
    println("# of escaped:\t\t", string(data["escape"]))
    println("# unresolved:\t\t", string(data["particles unresolved"]))
    nothing
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
function is_palindrome(str)
  str == reverse(str)
end

sum = 0
for i = 1:1000000
  if is_palindrome(dec(i)) && is_palindrome(bin(i))
    sum = i + sum
  end
end

println(sum)

max_stars_repo_name: vmchale/julia-problem-solving

max_issues_repo_name: vmchale/julia-problem-solving

max_forks_repo_name: vmchale/julia-problem-solving


text:
using Bridge, Distributions
using Test, LinearAlgebra

h = 1e-7
n, m = 50, 10000
T1 = 1.
T2 = 2.
T = T2-T1
ss = range(T1, stop=T2, length=n)
tt = Bridge.tofs(ss, T1, T2)

u = 0.5
v = 0.3
a = .7
a2 = 0.4
P = LinPro(-0.8, 0.0, sqrt(a))
P2 = LinPro(-0.8, 0.0, sqrt(a2))

la = 1
cs = Bridge.CSpline(T1, T2, la*P.B*u, la*P.B*v)
Po = BridgeProp(P, tt, (u, v), a, cs)
Pt = Bridge.ptilde(Po)

cs2 = Bridge.CSpline(T1, T2, P2.B*u, P2.B*v)
Po2 = BridgeProp(P2, tt, (u, v), a2, cs)
Pt2 = Bridge.ptilde(Po2)

s = ss[div(n,2)]
t = tt[div(n,2)]
@test Bridge.V(t, T2, v, Pt) ≈ Bridge.Vs(s, T1, T2, v, Pt)
@test Bridge.dotV(t, T2, v, Pt) ≈ Bridge.dotVs(s, T1, T2, v, Pt)
@test abs((Bridge.V(t+h, T2, v, Pt) - Bridge.V(t, T2, v, Pt))/h - Bridge.dotV(t, T2, v, Pt)) < h

@test tt[1] == T1
@test tt[end] == T2
@test (Bridge.V(T1, T2, v, Pt) - u)/T ≈ Bridge.uofx(T1, Po.v[1], T1, T2, v, Pt) #
@test [T1, u] ≈ [Bridge.txofsu(T1, Bridge.uofx(T1, Po.v[1], T1, T2, v, Pt), T1, T2, v, Pt)...]
@test norm(Bridge.soft(Bridge.tofs(1:0.1:2, 1, 2), 1,2 ) .- (1:0.1:2)) < sqrt(eps())

if la == 1
    @test norm(Bridge.b(T1, u, P) - Bridge.b(T1, u, Pt)) + norm(Bridge.b(T2, v, P) - Bridge.b(T2, v,Pt)) < sqrt(eps())
end

X = Bridge.solve(EulerMaruyama(), sample(ss, Wiener{Float64}()), Po)
Y = Bridge.solve(EulerMaruyama(), Bridge.innovations(EulerMaruyama(), X, Po), Po)
@test Y.tt ≈ X.tt
@test Y.yy ≈ X.yy

X = ubridge(sample(ss, Wiener{Float64}()), Po)
@test X.tt ≈ tt
Y = ubridge(Bridge.uinnovations(X, Po), Po)
@test Y.tt ≈ X.tt
@test Y.yy ≈ X.yy

X = Bridge.solve(Bridge.Mdb(), sample(ss, Wiener{Float64}()), Po)
Y = Bridge.solve(Bridge.Mdb(), innovations(Bridge.Mdb(), X, Po), Po)
@test Y.tt ≈ X.tt
@test Y.yy ≈ X.yy

p = pdf(transitionprob(T1, u, T2, P), v)
pt = exp(lptilde(Po))

C = []
Co = []
Cnames = []
push!(Cnames, "Euler")
z = Float64[
    let
        X = solve(EulerMaruyama(), sample(tt, Wiener{Float64}()), Po)
         Bridge.llikelihoodleft(X, Po)
    end
    for i in 1:m]

o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

push!(Cnames, "Euler + Trapez")
z = Float64[
    let
        X = solve(EulerMaruyama(), sample(tt, Wiener{Float64}()), Po)
         Bridge.llikelihoodtrapez(X, Po)
    end
    for i in 1:m]

o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

push!(Cnames, "MDGP+Left")
z = Float64[
    let
        X = solve(Bridge.Mdb(), sample(tt, Wiener{Float64}()), Po)
        Bridge.llikelihoodleft(X, Po)
    end
    for i in 1:m]

o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

push!(Cnames, "MDGP+Trapez")
z = Float64[
    let
        X = solve(Bridge.Mdb(), sample(tt, Wiener{Float64}()), Po)
        Bridge.llikelihoodtrapez(X, Po)
    end
    for i in 1:m]

o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

push!(Cnames, "TCSGP")
z = Float64[
    let
        X = ubridge(sample(ss, Wiener{Float64}()), Po)
        Bridge.llikelihoodleft(X, Po)
    end
    for i in 1:m]

o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

push!(Cnames, "TCSGP + Trapez")
z = Float64[
    let
        X = ubridge(sample(ss, Wiener{Float64}()), Po)
        Bridge.llikelihoodtrapez(X, Po)
    end
    for i in 1:m]

o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

push!(Cnames, "TCSGPLL")
z = Float64[
    let
        X = ubridge(sample(ss, Wiener{Float64}()), Po)
        ullikelihood(X, Po)
    end
    for i in 1:m]
o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

push!(Cnames, "TCSGPLL + Trapez")
z = Float64[
    let
        X = ubridge(sample(ss, Wiener{Float64}()), Po)
        Bridge.ullikelihoodtrapez(X, Po)
    end
    for i in 1:m]
o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

push!(Cnames, "TCSGP + Inno + Update")
z = Float64[
                  let
                      X = ubridge(sample(ss, Wiener{Float64}()), Po2)
                      Z = Bridge.uinnovations(X, Po2)
                      Z2 = sample(Z.tt, Wiener{Float64}())
                      Z.yy[:] = sqrt(.8)*Z.yy + sqrt(0.2)*Z2.yy
                      Z.yy[end]
                      X = ubridge(Z, Po)
                      Bridge.llikelihoodleft(X, Po)
                  end
                  for i in 1:m]

o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

push!(Cnames, "TCSGPLL + Inno Mix")
z = Float64[
           let
               W = sample(ss, Wiener{Float64}())
               X = ubridge(W, Po2)
               Z = Bridge.innovations(EulerMaruyama(), X, Po2)
               X = solve(EulerMaruyama(), Z, Po)
               Bridge.llikelihoodtrapez(X, Po)
           end
           for i in 1:m]
o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

push!(Cnames, "MDGP + Trapez + Inno")
z = Float64[
    let
        X = solve(Bridge.Mdb(), sample(tt, Wiener{Float64}()), Po2)
        W = innovations(EulerMaruyama(), X, Po2)
        X = solve(EulerMaruyama(), W, Po)
        Bridge.llikelihoodtrapez(X, Po)
    end
    for i in 1:m]

o = mean(exp.(z)*pt/p); push!(Co, o); push!(C, abs(o - 1)*sqrt(m)/std(exp.(z)*pt/p))

display([Cnames C Co])
println("\n  Method                   rel error  mean")
println("Remark: Change of measure is successfull if mean is close to 1")


text:
struct SpatialGrid{T <: AbstractFloat }
    x          :: OffsetVector{T}
    global_DoF :: Vector{Vector{Int64}}
    nfree      :: Int64
    nfixed     :: Int64
    max_r      :: Int64
    pt         :: Vector{T} # Quadrature points on (-1,1)
    wt         :: Vector{T} # Quadrature weights
end

"""
   Φ = shape_funcs(r, ξ)

Returns shape functions `Φ[i,j] = Φ_(i-1)(ξ[j])`.
"""
function shape_funcs(r::Int64, ξ::AbstractVector{T}
                    ) where { T <: AbstractFloat }
    J = length(ξ)
    Φ = zeros(T, r, J)
    for j = 1:J
        Φ[1,j] = ( 1 - ξ[j] ) / 2
        Φ[2,j] = ( 1 + ξ[j] ) / 2
    end
    P = Array{T}(undef, r, J)
    legendre_polys!(P, ξ) # P[i,j] = P_(i-1)(ξ[j])
    for j = 1:J, i = 3:r
        a = convert(T, 2(2i-3))
        Φ[i,j] = ( P[i,j] - P[i-2,j] ) / sqrt(a)
    end
    return Φ
end

function SpatialGrid(x::OffsetVector{T}, r::Vector{Int64},
                     bc::Tuple{String,String},
                     num_x_Gauss_pts::Int64) where { T <: AbstractFloat }
    P = length(x) - 1
    @argcheck length(r) == P DimensionMismatch
    global_DoF = Vector{Vector{Integer}}(undef, P)
    total =  0
    max_r = 0
    for p = 1:P
        @argcheck r[p] ≥ 2
        global_DoF[p] = Array{Integer}(undef, r[p])
        total += r[p]
        if r[p] ≥ max_r
            max_r = r[p]
        end
    end
    total -= P-1
    if bc == ("essential", "essential")
        nfixed = 2
        nfree = total - nfixed
        global_DoF[1][1] = nfree + 1
        m = 0
        if P > 1
            for j = 2:r[1]
                m += 1
                global_DoF[1][j] = m
            end
            for p = 2:P-1
                global_DoF[p][1] = global_DoF[p-1][2]
                for j = 2:r[p]
                    m += 1
                    global_DoF[p][j] = m
                end
            end
            global_DoF[P][1] = global_DoF[P-1][2]
            global_DoF[P][2] = nfree + 2
            for j = 3:r[P]
                m += 1
                global_DoF[P][j] = m
            end
        else
            global_DoF[1][2] = nfree + 2
            for j = 3:r[1]
                global_DoF[1][j] = j-2
            end
        end
    elseif bc == ("essential", "natural")
        nfixed = 1
        nfree = total - nfixed
        global_DoF[1][1] = nfree + 1
        m = 0
        for j = 2:r[1]
            m += 1
            global_DoF[1][j] = m
        end
        for p = 2:P
            global_DoF[p][1] = global_DoF[p-1][2]
            for j = 2:r[p]
                m += 1
                global_DoF[p][j] = m
            end
        end
    elseif bc == ("natural", "essential")
        nfixed = 1
        nfree = total - nfixed
        m = 0
        for j = 1:r[1]
            m += 1
            global_DoF[1][j] = m
        end
        for p = 2:P-1
            global_DoF[p][1] = global_DoF[p-1][2]
            for j = 2:r[p]
                m += 1
                global_DoF[p][j] = m
            end
        end
        global_DoF[P][1] = global_DoF[P-1][2]
        global_DoF[P][2] = nfree + 1
        for j = 3:r[P]
            m += 1
            global_DoF[P][j] = m
        end
    elseif bc == ("natural", "natural")
        nfixed = 0
        nfree = total - nfixed
        m = 0
        for j = 1:r[1]
            m += 1
            global_DoF[1][j] = m
        end
        for p = 2:P
            global_DoF[p][1] = global_DoF[p-1][2]
            for j = 2:r[p]
                m += 1
                global_DoF[p][j] = m
            end
        end
    else
        throw(ArgumentError("unrecognised boundary conditions"))
    end
    pt, wt = GaussQuadrature.legendre(T, num_x_Gauss_pts)
    return SpatialGrid{T}(x, global_DoF, nfree, nfixed, max_r, pt, wt)
end

function SpatialGrid(x::OffsetVector{T}, r::Int64,
                     bc::Tuple{String,String},
                     num_Gauss_pts) where { T <: AbstractFloat }
    P = length(x) - 1
    vecr = fill(r, P)
    return SpatialGrid(x, vecr, bc, num_Gauss_pts)
end

function ref_element_matrices!(M::Matrix{T}, S::Matrix{T},
                               κ::T) where { T <: AbstractFloat }
    sz = size(M, 1)
    @argcheck size(M, 2) == sz
    @argcheck size(S) == (sz, sz)
    @argcheck sz >= 2
    fill!(M, zero(T))
    fill!(S, zero(T))
    two = convert(T, 2)
    three = convert(T, 3)
    M[1, 1] = M[2, 2] =  2 / three
    M[2, 1] = M[1, 2] =  1 / three
    S[1, 1] = S[2, 2] =  κ / 2
    S[2, 1] = S[1, 2] = -κ / 2
    if sz == 2
        return
    end
    six = convert(T, 6)
    if sz >= 3
        M[1, 3] = M[2, 3] = -1 /sqrt(six)
        M[3, 1] = M[1, 3]
        M[3, 2] = M[2, 3]
    end
    if sz >= 4
        ten = convert(T, 10)
        M[1, 4] = 1 / ( 3 * sqrt(ten) )
        M[2, 4] = -M[1, 4]
        M[4, 1] = M[1, 4]
        M[4, 2] = M[2, 4]
    end
    for j = 3:sz
        M[j, j] = two / ((2j-5)*(2j-1))
        S[j, j] = κ
    end
    for j = 3:sz-2
        a = convert(T, (2j-3)*(2j+1))
        M[j, j+2] = M[j+2, j] = -1 / ( (2j-1) * sqrt(a) )
    end
end

function assembled_matrices(grid::SpatialGrid,
                            κ::T) where { T <: AbstractFloat }
    x, gdof = grid.x, grid.global_DoF
    nfree, nfixed, max_r = grid.nfree, grid.nfixed, grid.max_r
    P = length(gdof)
    Mref = zeros(T, max_r, max_r)
    Sref = zeros(T, max_r, max_r)
    ref_element_matrices!(Mref, Sref, κ)
    IM, JM, VM = Int64[], Int64[], Float64[]
    IS, JS, VS = Int64[], Int64[], Float64[]
    for p = 1:P
        hp = x[p] - x[p-1]
        rp = length(gdof[p])
        for i = 1:rp
            gi = gdof[p][i]
            if gi ≤ nfree
                for j = 1:rp
                    gj = gdof[p][j]
                    push!(IM, gi)
                    push!(JM, gj)
                    push!(VM, (hp/2)*Mref[i,j])
                    push!(IS, gi)
                    push!(JS, gj)
                    push!(VS, (2/hp)*Sref[i,j])
                end
            end
        end
    end
    M = SparseArrays.sparse(IM, JM, VM, nfree, nfree+nfixed)
    S = SparseArrays.sparse(IS, JS, VS, nfree, nfree+nfixed)
    return M, S
end

function element_load_vector!(Fp::AbstractArray{T},
                              elt::Tuple{T,T}, f::Function,
                              Φ::Matrix{T}, pt::Vector{T}, wt::Vector{T}
                             ) where { T <: AbstractFloat }
    r = length(Fp)
    J = length(pt)
    @argcheck size(Φ,1) ≤ r  DimensionMismatch
    @argcheck size(Φ,2) == J DimensionMismatch
    Jacobian = ( elt[2] - elt[1] ) / 2
    fill!(Fp, zero(T))
    for j = 1:J
        xj = ( ( 1 - pt[j] ) * elt[1] + ( 1 + pt[j] ) * elt[2] ) / 2
        for i = 1:r
            Fp[i] += wt[j] * f(xj) * Φ[i,j] * Jacobian
        end
    end
end

function assembled_load_vector!(F::Vector{T}, grid::SpatialGrid{T},
                                f::Function) where { T <: AbstractFloat }
    x, gdof, pt, wt = grid.x, grid.global_DoF, grid.pt, grid.wt
    nfree, nfixed, max_r = grid.nfree, grid.nfixed, grid.max_r
    @argcheck length(F) == nfree
    P = length(x) - 1
    Φ = shape_funcs(max_r, pt)
    Fp = zeros(T, max_r)
    fill!(F, zero(T))
    for p = 1:P
        rp = length(gdof[p])
        element_load_vector!(view(Fp, 1:rp), (x[p-1],x[p]), f, Φ, pt, wt)
        for i = 1:length(gdof[p])
            gi = gdof[p][i]
            if gi ≤ nfree
                F[gi] += Fp[i]
            end
        end
    end
end

function assembled_load_vector(grid::SpatialGrid{T},
                               f::Function) where { T <: AbstractFloat }
    nfree = grid.nfree
    F = zeros(T, nfree)
    assembled_load_vector!(F, grid, f)
    return F
end

function spatial_points!(xvals::Vector{T}, grid::SpatialGrid{T},
                         ξ::AbstractVector{T}) where { T <: AbstractFloat }
    x, gdof = grid.x, grid.global_DoF
    Nx = length(xvals)
    P = length(x) - 1
    pts_per_interval = length(ξ)
    n = 0
    if Nx == 1 + P * ( pts_per_interval - 1 )
        for p = 1:P
            for j = 1:pts_per_interval - 1
                xpj = ( ( 1 - ξ[j] ) * x[p-1] + (1 + ξ[j] ) * x[p] ) / 2
                n += 1
                xvals[n] = xpj
            end
        end
        n += 1
        xvals[n] = x[P]
    elseif Nx == P * pts_per_interval
        for p = 1:P
            for j = 1:pts_per_interval
                xpj = ( ( 1 - ξ[j] ) * x[p-1] + (1 + ξ[j] ) * x[p] ) / 2
                n += 1
                xvals[n] = xpj
            end
        end
    else
        throw(ArgumentError("xvals length does not match grid"))
    end
end

function evaluate_fem1d_soln!(Uvals::Vector{T}, U::Vector{T},
                              grid::SpatialGrid{T},
                              ξ::AbstractVector{T}
                             ) where { T <: AbstractFloat }
    x, gdof = grid.x, grid.global_DoF
    nfree, nfixed, max_r = grid.nfree, grid.nfixed, grid.max_r
    @argcheck length(U) == nfree DimensionMismatch
    Nx = length(Uvals)
    P = length(x) - 1
    pts_per_interval = length(ξ)
    Φ = shape_funcs(max_r, ξ)
    n = 0
    if Nx == 1 + P * ( pts_per_interval - 1 ) # ξ[1] = -1, ξ[end] = +1
        for p = 1:P
            rp = length(gdof[p])
            for j = 1:pts_per_interval-1
                n += 1
                Uvals[n] = zero(T)
                for i = 1:rp
                    gpi = gdof[p][i]
                    if gpi ≤ nfree
                        Uvals[n] += U[gpi] * Φ[i,j]
                    end
                end
            end
        end
        n += 1
        gpi = gdof[P][2]
        if gpi ≤ nfree
            Uvals[n] = U[gpi]
        end
    elseif Nx == P * pts_per_interval  # ξ[1] > -1 or ξ[end] < +1
        for p = 1:P
            rp = length(gdof[p])
            for j = 1:pts_per_interval
                n += 1
                Uvals[n] = zero(T)
                for i = 1:rp
                    gpi = gdof[p][i]
                    if gpi ≤ nfree
                        Uvals[n] += U[gpi] * Φ[i,j]
                    end
                end
            end
        end
    else
        throw(ArgumentError("Uvals length does not match grid"))
    end
end

function evaluate_fem1d_soln(U::Vector{T},
                             grid::SpatialGrid{T},
                             pts_per_interval::Int64
                            ) where { T <: AbstractFloat }
    x, nfree = grid.x, grid.nfree
    @argcheck length(U) == nfree DimensionMismatch
    P = ubound(x, 1)
    Nx = 1 + P * ( pts_per_interval - 1 )
    ξ = collect( range(-one(T), stop=one(T), length=pts_per_interval) )
    xvals = zeros(T, Nx)
    spatial_points!(xvals, grid, ξ)
    Uvals = zeros(T, Nx)
    evaluate_fem1d_soln!(Uvals, U, grid, ξ)
    return xvals, Uvals
end

function evaluate_pcwise_poly(U::Array{T}, t::OffsetVector{T},
                              grid::SpatialGrid{T},
                              pts_per_time_interval::Int64,
                              pts_per_space_interval::Int64
                             ) where { T <: AbstractFloat }
    τ = range(-one(T), stop=one(T), length=pts_per_time_interval)
    ξ = range(-one(T), stop=one(T), length=pts_per_space_interval)
    xvals, pcwise_t, pcwise_U = evaluate_pcwise_poly(U, t, grid, τ, ξ)
    return xvals, pcwise_t, pcwise_U
end

function evaluate_pcwise_poly(U::Array{T}, t::OffsetVector{T},
                              grid::SpatialGrid{T},
                              τ::AbstractVector{T},
                              ξ::AbstractVector{T}
                             ) where { T <: AbstractFloat }
    x = grid.x
    P = length(x) - 1
    nfree, rt, N = size(U)
    pts_per_space_interval = length(ξ)
    pts_per_time_interval = length(τ)
    @argcheck axes(t, 1) == 0:N
    if ξ[1] > -1.0 || ξ[end] < 1.0
        Nx = P * pts_per_space_interval
    else
        Nx = 1 + P * ( pts_per_space_interval - 1 )
    end
    xvals = zeros(T, Nx)
    spatial_points!(xvals, grid, ξ)
    Ψ = Array{Float64}(undef, rt, pts_per_time_interval)
    legendre_polys!(Ψ, τ)
    pcwise_t = zeros(T, pts_per_time_interval, N)
    pcwise_U = zeros(T, Nx, pts_per_time_interval, N)
    Uvals = zeros(T, Nx)
    for n = 1:N
        for i = 1:pts_per_time_interval
            pcwise_t[i,n] = ( (1-τ[i])*t[n-1] + (1+τ[i])*t[n] ) / 2
            for k = 1:rt
                evaluate_fem1d_soln!(Uvals, U[:,k,n], grid, ξ)
                pcwise_U[:,i,n] .+= Uvals * Ψ[k,i]
            end
        end
    end
    return xvals, pcwise_t, pcwise_U
end

function evaluate_pcwise_error(pcwise_U::Array{T}, u::Function,
                               pcwise_t::Array{T}, xvals::Vector{T}
                              ) where { T <: AbstractFloat }
    Nx = length(xvals)
    N = size(pcwise_t, 2)
    pts_per_time_interval = size(pcwise_t, 1)
    pcwise_err = zeros(Nx, pts_per_time_interval, N)
    for n = 1:N
        for i = 1:pts_per_time_interval
            for p = 1:Nx
                pcwise_err[p,i,n] = ( pcwise_U[p,i,n]
                                     - u(xvals[p], pcwise_t[i,n]) )
            end
        end
    end
    return pcwise_err
end

function dynamic_load_vector!(Fn::Matrix{T}, In::Vector{Float64},
                              grid::SpatialGrid{T}, f::Function,
                              τ::Vector{T}, w::Vector{T},
                              Ψ::AbstractMatrix{T}) where { T <: AbstractFloat }
    nfree, rt = size(Fn)
    @argcheck grid.nfree == nfree
    @argcheck length(In) == 2
    fill!(Fn, 0.0)
    Fni = zeros(nfree)
    Δt = In[2] - In[1]
    for i = 1:length(τ)
        ti = ( ( 1 - τ[i] ) * In[1] + ( 1 + τ[i] ) * In[2] ) / 2
        assembled_load_vector!(Fni, grid, x -> f(x,ti))
        for j = 1:rt
            Fn[:,j] += (Δt/2) * w[i] * Fni * Ψ[j,i]
        end
    end
end

function PDEDG(κ::Float64, f::Function, U0::Vector{Float64},
               grid::SpatialGrid{Float64}, t::OffsetVector{Float64},
               rt::Int64, Mt::Integer) where T <: AbstractFloat
    x, gdof = grid.x, grid.global_DoF
    nfree, nfixed = grid.nfree, grid.nfixed

    G = coef_G(Float64, rt)
    K = coef_K(Float64, rt, rt)
    H = zeros(rt, rt)
    for j = 1:rt
        H[j,j] = 1 / (2j-1)
    end

    N = length(t) - 1
    U = zeros(nfree*rt, N)
    Fn = zeros(nfree, rt)
    M, S = assembled_matrices(grid, κ)
    M = M[1:nfree,1:nfree]
    S = S[1:nfree,1:nfree]
    Ψ_at_minus1 = Vector{Float64}(undef, rt)
    legendre_polys!(Ψ_at_minus1, -1.0)
    τ, w = GaussQuadrature.legendre(Mt)
    Ψ = Array{Float64}(undef, rt, Mt)
    legendre_polys!(Ψ, τ)
    dynamic_load_vector!(Fn, t[0:1], grid, f, τ, w, Ψ)
    b = Fn[1:nfree*rt] + kron(Ψ_at_minus1, M*U0)
    Δt1 = t[1] - t[0]
    G_kron_M = kron(G, M)
    H_kron_S = kron(H, S)
    A = G_kron_M + Δt1 * H_kron_S
    U[:,1] = A \ b
    K_kron_M = kron(K, M)
    for n = 2:N
        Δtn = t[n] - t[n-1]
        dynamic_load_vector!(Fn, t[n-1:n], grid, f, τ, w, Ψ)
        b = Fn[1:nfree*rt] + K_kron_M * U[:,n-1]
        A = G_kron_M + Δtn * H_kron_S
        U[:,n] = A \ b
    end
    return reshape(U, nfree, rt, N)
end

"""
    t, U = PDEDG(κ, f, U0, grid, N, rt, num_Gauss_pts)

Faster version for uniform time steps.
"""
function PDEDG(κ::Float64, f::Function, U0::Vector{Float64},
               grid::SpatialGrid{Float64}, max_t::Float64, N::Int64,
               rt::Int64, Mt::Integer)
    x, gdof = grid.x, grid.global_DoF
    nfree, nfixed = grid.nfree, grid.nfixed

    G = coef_G(Float64, rt)
    K = coef_K(Float64, rt, rt)
    H = zeros(rt, rt)
    for j = 1:rt
        H[j,j] = 1 / (2j-1)
    end

    t_ = collect(range(0, max_t, length=N+1))
    t = OffsetVector(t_, 0:N)
    U = zeros(nfree*rt, N)
    Fn = zeros(nfree, rt)
    M, S = assembled_matrices(grid, κ)
    M = M[1:nfree,1:nfree]
    S = S[1:nfree,1:nfree]
    Ψ_at_minus1 = Vector{Float64}(undef, rt)
    legendre_polys!(Ψ_at_minus1, -1.0)
    τ, w = GaussQuadrature.legendre(Mt)
    Ψ = Array{Float64}(undef, rt, Mt)
    legendre_polys!(Ψ, τ)
    dynamic_load_vector!(Fn, t[0:1], grid, f, τ, w, Ψ)
    b = Fn[1:nfree*rt] + kron(Ψ_at_minus1, M*U0)
    Δt = max_t / N
    G_kron_M = kron(G, M)
    H_kron_S = kron(H, S)
    A = G_kron_M + Δt * H_kron_S
    LU = SparseArrays.lu(A)
    U[:,1] = LU \ b
    K_kron_M = kron(K, M)
    for n = 2:N
        dynamic_load_vector!(Fn, t[n-1:n], grid, f, τ, w, Ψ)
        b = Fn[1:nfree*rt] + K_kron_M * U[:,n-1]
        U[:,n] = LU \ b
    end
    return t, reshape(U, nfree, rt, N)
end

function spatial_L2_projection(f::Function, grid::SpatialGrid{Float64})
    M, S = assembled_matrices(grid, 0.0)
    nfree = grid.nfree
    Pf = M[1:nfree,1:nfree] \ assembled_load_vector(grid, f)
    return Pf
end

"""
    err_τ = maxerr(τ, U, t, u, ξ, grid, cutoff=0.0)

Returns err_τ[i] = max |U(x,t)-u(x,t)| over x values in each spatial
interval corresponding to ξ the values in [-1,1], and over the t value
in each time subinterval corresponding to τ[i] in [-1,1], except for time
intervals with t[n] ≤ cutoff.
"""
function maxerr(τ:: AbstractVector{Float64}, U::Array{Float64},
                t::OffsetVector{Float64}, u::Function,
                ξ::AbstractVector{Float64}, grid::SpatialGrid,
                cutoff=0.0)
    xvals, pcwise_t, pcwise_U = evaluate_pcwise_poly(U, t, grid, τ, ξ)
    pcwise_err = evaluate_pcwise_error(pcwise_U, u, pcwise_t, xvals)
    N = size(pcwise_t, 2)
    pts_per_time_interval = size(pcwise_t, 1)
    Nx = size(pcwise_err, 1)
    err_τ = zeros(pts_per_time_interval)
    for n = 1:N
        if t[n] ≤ cutoff
            continue
        end
        for i = 1:pts_per_time_interval
            for p = 1:Nx
                err_τ[i] = max(abs(pcwise_err[p,i,n]), err_τ[i])
            end
        end
    end
    return err_τ
end

"""
    U = FPDEDG(κ, f, U0, grid, t, rt, num_t_Gauss_pts)

DG solver for fractional PDE.
"""
function FPDEDG(κ::Float64, f::Function, U0::Vector{Float64},
               grid::SpatialGrid{Float64}, t::OffsetVector{Float64},
               rt::Int64, Mt::Int64, store::Store{Float64})
    α = store.α
    x, gdof = grid.x, grid.global_DoF
    nfree, nfixed = grid.nfree, grid.nfixed

    G = coef_G(Float64, rt)
    K = coef_K(Float64, rt, rt)
    Nt = length(t) - 1
    Hn = OffsetArray{Matrix{Float64}}(undef, 0:Nt-1)
    for n = 0:Nt-1
        Hn[n] = Array{Float64}(undef, rt, rt)
    end
    τ, wτ = rule(store.legendre[Mt])
    Ψ = view(store.Ψ, 1:rt, 1:Mt)
    legendre_polys!(Ψ, τ)

    U = zeros(nfree*rt, Nt)
    Fn = zeros(nfree, rt)
    M, S = assembled_matrices(grid, κ)
    M = M[1:nfree,1:nfree]
    S = S[1:nfree,1:nfree]
    Ψ_at_minus1 = Array{Float64}(undef, rt)
    legendre_polys!(Ψ_at_minus1, -1.0)
    dynamic_load_vector!(Fn, t[0:1], grid, f, τ, wτ, Ψ)
    b = Fn[1:nfree*rt] + kron(Ψ_at_minus1, M*U0)
    coef_Hn!(Hn, 1, 0, t, Mt, store)
    G_kron_M = kron(G, M)
    Hnn_kron_S = kron(Hn[0], S)
    A = G_kron_M + Hnn_kron_S
    U[:,1] = A \ b
    K_kron_M = kron(K, M)
    for n = 2:Nt
        coef_Hn!(Hn, n, n-1, t, Mt, store)
        Hnn_kron_S = kron(Hn[0], S)
        A = G_kron_M + Hnn_kron_S
        fill!(b, 0.0)
        for l = 1:n-1
            Hnl_kron_S = kron(Hn[n-l], S)
            b .= b .+ Hnl_kron_S * U[:,l]
        end
        dynamic_load_vector!(Fn, t[n-1:n], grid, f, τ, wτ, Ψ)
        b .= Fn[1:nfree*rt] .+ K_kron_M * U[:,n-1] .- b
        U[:,n] = A \ b
    end
    return reshape(U, nfree, rt, Nt)
end

"""
    t, U = FPDEDG(κ, f, U0, grid, max_t, N, rt, num_t_Gauss_pts)

Specialised version for uniform time steps.
"""
function FPDEDG(κ::Float64, f::Function, U0::Vector{Float64},
                grid::SpatialGrid{Float64}, max_t::Float64,
                N::Int64, rt::Int64, Mt::Int64, store::Store{Float64})
    α = store.α
    x, gdof = grid.x, grid.global_DoF
    nfree, nfixed = grid.nfree, grid.nfixed

    Δt = max_t / N
    G = coef_G(Float64, rt)
    K = coef_K(Float64, rt, rt)
    H = coef_H_uniform!(rt, N, Mt, store)

    t_ = collect(range(0, max_t, length=N+1))
    t = OffsetVector(t_, 0:N)

    U = zeros(nfree*rt, N)
    Fn = zeros(nfree, rt)
    M, S = assembled_matrices(grid, κ)
    M = M[1:nfree,1:nfree]
    S = S[1:nfree,1:nfree]
    Ψ_at_minus1 = Vector{Float64}(undef, rt)
    legendre_polys!(Ψ_at_minus1, -1.0)
    τ, w = GaussQuadrature.legendre(Mt)
    Ψ = Array{Float64}(undef, rt, Mt)
    legendre_polys!(Ψ, τ)
    dynamic_load_vector!(Fn, t[0:1], grid, f, τ, w, Ψ)
    b = Fn[1:nfree*rt] + kron(Ψ_at_minus1, M*U0)
    G_kron_M = kron(G, M)
    Hnn_kron_S = kron(H[0], S)
    A = G_kron_M + Δt^α * Hnn_kron_S
    Fact = SparseArrays.lu(A)
    U[:,1] = Fact \ b
    K_kron_M = kron(K, M)
    for n = 2:N
        fill!(b, 0.0)
        for l = 1:n-1
            Hnl_kron_S = kron(H[n-l], S)
            b .= b .+ Δt^α * Hnl_kron_S * U[:,l]
        end
        dynamic_load_vector!(Fn, t[n-1:n], grid, f, τ, w, Ψ)
        b .= Fn[1:nfree*rt] .+ K_kron_M * U[:,n-1] .- b
        U[:,n] = Fact \ b
    end
    return t, reshape(U, nfree, rt, N)
end

function functional(pcwise_U::Array{Float64}, pcwise_t::Array{Float64},
                    func::Function, grid::SpatialGrid,
                    ξ::Vector{Float64}, w::Vector{Float64})
    x = grid.x
    pts_per_time_interval = size(pcwise_t,1)
    N = size(pcwise_t, 2)
    P = ubound(x, 1)
    pts_per_space_interval = length(ξ)
    Fnctl = zeros(pts_per_time_interval, N)

end


text:
#####
##### two-dimensional fluxes
#####
"
function euler_fluxes_2D(rhoL,uL,vL,betaL,rhoR,uR,vR,betaR,
                         rhologL,betalogL,rhologR,betalogR)
assumes primitive variables ordering: UL = (rhoL,uL,...,betaL),
                                      UR = (rhoR,uR,...,betaR)
"
function euler_fluxes_2D(rhoL,uL,vL,betaL,rhologL,betalogL,
                         rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))
    vavg   = (@. .5*(vL+vR))

    unorm = (@. uL*uR + vL*vR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. FxS1*vavg)
    FxS4 = (@. f4aux*uavg)

    FyS1 = (@. rholog*vavg)
    FyS2 = (@. FxS3)
    FyS3 = (@. FyS1*vavg + pa)
    FyS4 = (@. f4aux*vavg)
    return (FxS1,FxS2,FxS3,FxS4),(FyS1,FyS2,FyS3,FyS4)
end

function euler_fluxes_2D_x(rhoL,uL,vL,betaL,rhoR,uR,vR,betaR)
    rhologL,betalogL,rhologR,betalogR = map(x->log.(x),(rhoL,betaL,rhoR,betaR))
    return euler_fluxes_2D_x(rhoL,uL,vL,betaL,rhologL,betalogL,
                             rhoR,uR,vR,betaR,rhologR,betalogR)
end
function euler_fluxes_2D_y(rhoL,uL,vL,betaL,rhoR,uR,vR,betaR)
    rhologL,betalogL,rhologR,betalogR = map(x->log.(x),(rhoL,betaL,rhoR,betaR))
    return euler_fluxes_2D_y(rhoL,uL,vL,betaL,rhologL,betalogL,
                             rhoR,uR,vR,betaR,rhologR,betalogR)
end

function euler_fluxes_2D_x(rhoL,uL,vL,betaL,rhologL,betalogL,
                           rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))
    vavg   = (@. .5*(vL+vR))

    unorm = (@. uL*uR + vL*vR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. FxS1*vavg)
    FxS4 = (@. f4aux*uavg)

    return (FxS1,FxS2,FxS3,FxS4)
end

function euler_fluxes_2D_y(rhoL,uL,vL,betaL,rhologL,betalogL,
                           rhoR,uR,vR,betaR,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))
    vavg   = (@. .5*(vL+vR))

    unorm = (@. uL*uR + vL*vR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FyS1 = (@. rholog*vavg)
    FyS2 = (@. FyS1*uavg)
    FyS3 = (@. FyS1*vavg + pa)
    FyS4 = (@. f4aux*vavg)
    return (FyS1,FyS2,FyS3,FyS4)
end

#####
##### shim functions
#####

# dispatch to n-dimensional constitutive routines, with optional entropy scaling
function primitive_to_conservative(rho,u,v,p)
   rho,rhoU,E = primitive_to_conservative_nd(rho,(u,v),p)
   return rho,rhoU...,E
end
function v_ufun(rho,rhou,rhov,E)
    v1,vU,vE = v_ufun_nd(rho,(rhou,rhov),E)
    return scale_entropy_output(v1,vU...,vE)
end
function u_vfun(v1,vU1,vU2,vE)
    v1,vU1,vU2,vE = scale_entropy_input(v1,vU1,vU2,vE)
    rho,rhoU,E = u_vfun_nd(v1,(vU1,vU2),vE)
    return rho,rhoU...,E
end
function conservative_to_primitive_beta(rho,rhou,rhov,E)
    rho,U,beta = conservative_to_primitive_beta_nd(rho,(rhou,rhov),E)
    return rho,U...,beta
end

Sfun(rho,rhou,rhov,E) = Sfun_nd(rho,(rhou,rhov),E)
pfun(rho,rhou,rhov,E) = pfun_nd(rho,(rhou,rhov),E)
betafun(rho,rhou,rhov,E) = betafun_nd(rho,(rhou,rhov),E)

"function euler_fluxes_2D(rhoL,uL,vL,betaL,rhoR,uR,vR,betaR)
assumes primitive variables ordering: UL = (rhoL,uL,...,betaL),
                                      UR = (rhoR,uR,...,betaR)"
function euler_fluxes(rhoL,uL,vL,betaL,rhoR,uR,vR,betaR)
    rhologL,betalogL,rhologR,betalogR = map(x->log.(x),(rhoL,betaL,rhoR,betaR))
    return euler_fluxes_2D(rhoL,uL,vL,betaL,rhologL,betalogL,
                           rhoR,uR,vR,betaR,rhologR,betalogR)
end

euler_fluxes(rhoL,uL,vL,betaL,rhologL,betalogL,
             rhoR,uR,vR,betaR,rhologR,betalogR) =
     euler_fluxes_2D(rhoL,uL,vL,betaL,rhologL,betalogL,
                     rhoR,uR,vR,betaR,rhologR,betalogR)

euler_fluxes(varL,varR) = euler_fluxes(varL...,varR...)
euler_fluxes(varL,varR,logL,logR) = euler_fluxes(varL...,logL...,varR...,logR...)

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
@doc raw"""
Axisymmetric tokamak equilibrium in (R,Z,ϕ) coordinates with covariant
components of the vector potential given by
```math
A (R, Z, \phi) = \frac{B_0}{2} \, \bigg( R_0 \, \frac{Z}{R} , \, - R_0 \, \ln \bigg( \frac{R}{R_0} \bigg) , \, - \frac{r^2}{q_0} \bigg)^T ,
```
resulting in the magnetic field with covariant components
```math
B (R, Z, \phi) = \frac{B_0}{q_0} \, \bigg( - \frac{Z}{R} , \, \frac{R - R_0}{R} , \, - q_0 R_0 \bigg)^T ,
```
where $r = \sqrt{ (R - R_0)^2 + Z^2 }$.

Parameters:
 * `R₀`: position of magnetic axis
 * `B₀`: B-field at magnetic axis
 * `q₀`: safety factor at magnetic axis
"""
module AxisymmetricTokamakCylindrical

    using RecipesBase

    import ..ElectromagneticFields
    import ..ElectromagneticFields: AnalyticEquilibrium, code

    export  AxisymmetricTokamakCylindricalEquilibrium

    const DEFAULT_R₀ = 1.0
    const DEFAULT_B₀ = 1.0
    const DEFAULT_q₀ = 2.0

    const ITER_R₀ = 6.2
    const ITER_B₀ = 5.3
    const ITER_q₀ = √2

    struct AxisymmetricTokamakCylindricalEquilibrium{T <: Number} <: AnalyticEquilibrium
        name::String
        R₀::T
        B₀::T
        q₀::T

        function AxisymmetricTokamakCylindricalEquilibrium{T}(R₀::T, B₀::T, q₀::T) where T <: Number
            new("AxisymmetricTokamakCylindricalEquilibrium", R₀, B₀, q₀)
        end
    end

    AxisymmetricTokamakCylindricalEquilibrium(R₀::T=DEFAULT_R₀, B₀::T=DEFAULT_B₀, q₀::T=DEFAULT_q₀) where T <: Number = AxisymmetricTokamakCylindricalEquilibrium{T}(R₀, B₀, q₀)

    function init(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀)
        AxisymmetricTokamakCylindricalEquilibrium(R₀, B₀, q₀)
    end

    function ITER()
        AxisymmetricTokamakCylindricalEquilibrium(ITER_R₀, ITER_B₀, ITER_q₀)
    end

    macro code(R₀=DEFAULT_R₀, B₀=DEFAULT_B₀, q₀=DEFAULT_q₀)
        code(init(R₀, B₀, q₀); escape=true)
    end

    macro code_iter()
        code(ITER(); escape=true)
    end

    function Base.show(io::IO, equ::AxisymmetricTokamakCylindricalEquilibrium)
        print(io, "Axisymmetric Tokamak Equilibrium in (R,Z,ϕ) Coordinates with\n")
        print(io, "  R₀ = ", equ.R₀, "\n")
        print(io, "  B₀ = ", equ.B₀, "\n")
        print(io, "  q₀ = ", equ.q₀)
    end

    R(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[1]
    Z(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[2]
    ϕ(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[3]
    r²(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2
    r(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = sqrt(r²(x, equ))
    X(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x,equ) * cos(ϕ(x,equ))
    Y(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x,equ) * sin(ϕ(x,equ))
    θ(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = atan(Z(x,equ), R(x,equ) - equ.R₀)

    ElectromagneticFields.J(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x,equ)

    ElectromagneticFields.A₁(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = + equ.B₀ * equ.R₀ * Z(x,equ) / R(x,equ) / 2
    ElectromagneticFields.A₂(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2
    ElectromagneticFields.A₃(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = + equ.B₀ * r²(x,equ) / equ.q₀ / 2

    ElectromagneticFields.x¹(ξ::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = X(ξ,equ)
    ElectromagneticFields.x²(ξ::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = Y(ξ,equ)
    ElectromagneticFields.x³(ξ::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = Z(ξ,equ)

    ElectromagneticFields.ξ¹(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = sqrt(x[1]^2 + x[2]^2)
    ElectromagneticFields.ξ²(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = x[3]
    ElectromagneticFields.ξ³(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = atan(x[2], x[1])

    ElectromagneticFields.g₁₁(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = one(eltype(x))
    ElectromagneticFields.g₂₂(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = one(eltype(x))
    ElectromagneticFields.g₃₃(x::AbstractVector, equ::AxisymmetricTokamakCylindricalEquilibrium) = R(x, equ)^2

    ElectromagneticFields.get_functions(::AxisymmetricTokamakCylindricalEquilibrium) = (X=X, Y=Y, Z=Z, R=R, r=r, θ=θ, ϕ=ϕ, r²=r²)

    function ElectromagneticFields.periodicity(x::AbstractVector, ::AxisymmetricTokamakCylindricalEquilibrium)
        p = zero(x)
        p[3] = 2π
        return p
    end

    @recipe function f(equ::AxisymmetricTokamakCylindricalEquilibrium;
                       nx = 100, ny = 120, levels = 50, size = (400,400),
                       xlims = (  0.5 * equ.R₀,   1.5 * equ.R₀),
                       ylims = (- 0.5 * equ.R₀, + 0.5 * equ.R₀))

        xgrid = LinRange(xlims[1], xlims[2], nx)
        zgrid = LinRange(ylims[1], ylims[2], ny)
        pot   = [ElectromagneticFields.A₃([xgrid[i], zgrid[j], 0.0], equ) / xgrid[i] for i in eachindex(xgrid), j in eachindex(zgrid)]

        seriestype   := :contour
        aspect_ratio := :equal
        size   := size
        xlims  := xlims
        ylims  := ylims
        levels := levels
        legend := :none

        (xgrid, zgrid, pot')
    end

end