text:
function get_coeffs_composition(g::Vector{Float64})
    n = length(g)
    a = zeros(n+1)
    b = vcat(g, 0)
    a[1] = g[1]/2
    a[n+1] = a[1]
    for j=2:div(n,2)+1
        a[j] = (g[j-1]+g[j])/2
        a[n-j+2] = a[j]
    end
    a,b
end

function gen_interpolation_matrix(x::Vector{Float64}, z::Vector{Float64})
    n = length(x)
    m = length(z)
    L = zeros(n, m)
    for j=1:m
        for i=1:n
            L[i,j] = prod([(z[j]-x[k])/(x[i]-x[k]) for k=1:i-1])*prod([(z[j]-x[k])/(x[i]-x[k]) for k=i+1:n])
        end
    end
    L
end

function gen_interpolation_matrix2(x::Vector{Float64}, z::Vector{Float64})
    L0 = gen_interpolation_matrix(x, [0.0; z])
    L = 0.5*(L0[:,1:end-1]+L0[:,2:end])
end

function gen_extrapolated_wf!(psi::WaveFunction, L::Vector{Float64}, psi_back::Vector{WaveFunction}, ptr::Int;
                              order::Int=length(psi_back)-1)
    n = length(psi_back)
    set!(psi, 0)
    for j=1:order+1
        k = mod(j-order+ptr-2, n)+1
        axpy!(psi, psi_back[k], L[j])
    end
end

function B_step!(psi::WaveFunction, psi_ex::WaveFunction, psi0::WaveFunction, dt::Float64, c::Float64, d::Float64;
                 nonlinear_potential::Bool=false, B_iter::Int=0)
    rhs = psi_ex
    #rhs = wave_function(psi.m)
    if B_iter>0
        copy!(psi0, psi)
    end
    for it=0:B_iter
        set_time!(psi_ex, get_time(psi))
        if nonlinear_potential
            gen_nonlinear_potential!(rhs, psi_ex)
            to_real_space!(psi)
            to_real_space!(rhs)
            u = get_data(psi, true)
            f = get_data(rhs, true)
            u[:] .*= exp(dt*f)
        else
            gen_rhs!(rhs, psi_ex)
            axpy!(psi, rhs, dt)
        end
        if B_iter>0 && it<B_iter
            copy!(psi_ex,psi0)
            scale!(psi_ex, c)
            axpy!(psi_ex, psi, d)
            copy!(psi, psi0)
        end
    end
end

function gen_starting_values(psi::WaveFunction, dt::Number, N::Int, a::Vector{Float64}, b::Vector{Float64};
    psi_ex::WaveFunction=wave_function(psi.m), rhs::WaveFunction=wave_function(psi.m), nonlinear_potential::Bool=false, B_iter::Int=0)
    s = length(a)
    aa = cumsum(a)
    bb = cumsum(b)
    c = (aa-[0.0; bb[1:end-1]])./b
    d = (bb-aa)./b
    psi_back = WaveFunction[wave_function(psi.m) for j=1:N]
    copy!(psi_back[1], psi)
    t0 = get_time(psi)
    for K=2:N
        copy!(psi, psi_back[1])
        for J=2:K
            L = gen_interpolation_matrix(collect(0.0:(K-2.0)), J-2.0+aa)
            #L = gen_interpolation_matrix2(collect(0.0:(K-2.0)), J-2.0+cumsum(b))
            for j = 1:s
                if a[j]!=0.0
                    propagate_A!(psi, a[j]*dt)
                end
                if b[j]!=0.0
                    gen_extrapolated_wf!(psi_ex, L[:,j], psi_back[1:K-1], 2*K-2)
                    B_step!(psi, psi_ex, rhs, b[j]*dt, c[j], d[j], nonlinear_potential=nonlinear_potential, B_iter=B_iter)
                    #set_time!(psi_ex, get_time(psi))
                    #if nonlinear_potential
                    #    gen_nonlinear_potential!(rhs, psi_ex)
                    #    to_real_space!(psi)
                    #    to_real_space!(rhs)
                    #    u = get_data(psi, true)
                    #    f = get_data(rhs, true)
                    #    u[:] .*= exp(b[j]*dt*f)
                    #else
                    #    gen_rhs!(rhs, psi_ex)
                    #    axpy!(psi, rhs, b[j]*dt)
                    #end
                end
            end
            copy!(psi_back[J], psi)
            to_real_space!(psi_back[J])
        end
    end
    psi_back
end

function gen_starting_values3(psi::WaveFunction, dt::Number, N::Int, a::Vector{Float64}, b::Vector{Float64};
    rhs::WaveFunction=wave_function(psi.m), nonlinear_potential::Bool=false)
    s = length(a)
    aa = cumsum(a)
    psi0 = wave_function(psi.m)
    rhs_back = WaveFunction[wave_function(psi.m) for j=1:N]
    copy!(psi0, psi)
    gen_rhs!(rhs_back[1], psi)
    t0 = get_time(psi)
    for K=2:N
        copy!(psi, psi0)
        for J=2:K
            L = gen_interpolation_matrix(collect(0.0:(K-2.0)), J-2.0+aa)
            for j = 1:s
                if a[j]!=0.0
                    propagate_A!(psi, a[j]*dt)
                end
                if b[j]!=0.0
                    gen_extrapolated_wf!(rhs, L[:,j], rhs_back[1:K-1], 2*K-2)
                    if nonlinear_potential
                        to_real_space!(psi)
                        to_real_space!(rhs)
                        u = get_data(psi, true)
                        f = get_data(rhs, true)
                       u[:] .*= exp(b[j]*dt*f)
                    else
                        axpy!(psi, rhs, b[j]*dt)
                    end
                end
            end
            if nonlinear_potential
                gen_nonlinear_potential!(rhs_back[J], psi)
            else
                gen_rhs!(rhs_back[J], psi)
            end
        end
    end
    rhs_back
end

type SplittingWithExtrapolationEquidistantTimeStepperIterator
    psi::WaveFunction
    t0::Real
    tend::Real
    dt::Real
    steps::Int
    a::Vector{Float64}
    b::Vector{Float64}
    L::Matrix{Float64}
    ptr::Int
    psi_back::Vector{WaveFunction} # storage for previous solution values
    psi_ex::WaveFunction  # storage for the extrapolated values
    rhs::WaveFunction     # storage for the righthandside of the B equation
    nonlinear_potential::Bool
    B_iter::Int
    c::Vector{Float64}
    d::Vector{Float64}
end

function splitting_with_extrapolation_equidistant_time_stepper(psi::WaveFunction,
         t0::Real, tend::Real, dt::Real, order::Int, a::Vector{Float64}, b::Vector{Float64}; steps::Int=-1, nonlinear_potential::Bool=false, B_iter::Int=0)
    m = psi.m
    psi_ex = wave_function(m)
    rhs = wave_function(m)
    set_propagate_time_together_with_A!(m, true)
    set_time!(psi, t0)
    psi_back = gen_starting_values(psi, dt, order+1, a, b, psi_ex=psi_ex, rhs=rhs, nonlinear_potential=nonlinear_potential, B_iter=B_iter)
    copy!(psi, psi_back[order+1])
    aa = cumsum(a)
    bb = cumsum(b)
    c = (aa-[0.0; bb[1:end-1]])./b
    d = (bb-aa)./b
    L = gen_interpolation_matrix(collect((-order-0.0):0.0), aa)
    #L = gen_interpolation_matrix2(collect((-order-0.0):0.0), cumsum(b))
    ptr = order+1
    SplittingWithExtrapolationEquidistantTimeStepperIterator(psi, t0, tend, dt, steps, a, b, L, ptr, psi_back, psi_ex, rhs, nonlinear_potential, B_iter, c, d)
end

Base.start(tsi::SplittingWithExtrapolationEquidistantTimeStepperIterator) = (tsi.steps<0 ?
    tsi.t0 + (length(tsi.psi_back)-1)*tsi.dt : length(tsi.psi_back)-1)

Base.done(tsi::SplittingWithExtrapolationEquidistantTimeStepperIterator, t) = (tsi.steps<0 ? t >= tsi.tend : t>=tsi.steps )

function Base.next(tsi::SplittingWithExtrapolationEquidistantTimeStepperIterator, t)
    for j = 1:length(tsi.a)
        if tsi.a[j]!=0.0
            propagate_A!(tsi.psi, tsi.a[j]*tsi.dt)
        end
        if tsi.b[j]!=0.0
            gen_extrapolated_wf!(tsi.psi_ex, tsi.L[:,j], tsi.psi_back, tsi.ptr)
            B_step!(tsi.psi, tsi.psi_ex, tsi.rhs, tsi.b[j]*tsi.dt, tsi.c[j], tsi.d[j],
                    nonlinear_potential=tsi.nonlinear_potential, B_iter=tsi.B_iter)
            #set_time!(tsi.psi_ex, get_time(tsi.psi))
            #if tsi.nonlinear_potential
            #    gen_nonlinear_potential!(tsi.rhs, tsi.psi_ex)
            #    to_real_space!(tsi.psi)
            #    to_real_space!(tsi.rhs)
            #    u = get_data(tsi.psi, true)
            #    f = get_data(tsi.rhs, true)
            #    u[:] .*= exp(tsi.b[j]*tsi.dt*f)
            #else
            #    gen_rhs!(tsi.rhs, tsi.psi_ex)
            #    axpy!(tsi.psi, tsi.rhs, tsi.b[j]*tsi.dt)
            #end
        end
    end
    n = length(tsi.psi_back)
    tsi.ptr = mod(tsi.ptr, n) + 1
    copy!(tsi.psi_back[tsi.ptr], tsi.psi)
    to_real_space!(tsi.psi_back[tsi.ptr])
    if tsi.steps<0
        t1 = t + tsi.dt < tsi.tend ? t + tsi.dt : tsi.tend
    else
        t1 = t+1
    end
    return t1, t1
end

type SplittingWithExtrapolationEquidistantTimeStepperIterator2
    psi::WaveFunction
    t0::Real
    tend::Real
    dt::Real
    steps::Int
    a::Vector{Float64}
    b::Vector{Float64}
    L::Matrix{Float64}
    L1::Matrix{Float64}
    order::Int
    order1::Int
    N::Int
    iter::Int
    ptr::Int
    psi_back::Vector{WaveFunction} # storage for previous solution values
    psi_ex::WaveFunction  # storage for the extrapolated values
    rhs::WaveFunction     # storage for the righthandside of the B equation
    nonlinear_potential::Bool
end

function splitting_with_extrapolation_equidistant_time_stepper2(psi::WaveFunction,
    t0::Real, tend::Real, dt::Real, order::Int, order1::Int, a::Vector{Float64}, b::Vector{Float64}; steps::Int=-1, nonlinear_potential::Bool=false, iter::Int=0, psi_back=nothing)
    m = psi.m
    psi_ex = wave_function(m)
    rhs = wave_function(m)
    set_propagate_time_together_with_A!(m, true)
    set_time!(psi, t0)
    N = max(order, order1)+1
    if psi_back!=nothing
        psi_back1 = psi_back
        @assert length(psi_back)==N
    else
        psi_back1 = gen_starting_values(psi, dt, N, a, b, psi_ex=psi_ex, rhs=rhs, nonlinear_potential=nonlinear_potential)
    end
    copy!(psi, psi_back1[N])
    aa = cumsum(a)
    L = gen_interpolation_matrix(collect((-order-0.0):0.0), aa)
    L1 = gen_interpolation_matrix(collect((-order1-0.0):0.0), aa-1.0)
    ptr = order+1
    SplittingWithExtrapolationEquidistantTimeStepperIterator2(psi, t0, tend, dt, steps, a, b, L, L1, order, order1, N, iter, ptr, psi_back1, psi_ex, rhs, nonlinear_potential)
end

Base.start(tsi::SplittingWithExtrapolationEquidistantTimeStepperIterator2) = (tsi.steps<0 ?
    tsi.t0 + (length(tsi.psi_back)-1)*tsi.dt : length(tsi.psi_back)-1)

Base.done(tsi::SplittingWithExtrapolationEquidistantTimeStepperIterator2, t) = (tsi.steps<0 ? t >= tsi.tend : t>=tsi.steps )

function Base.next(tsi::SplittingWithExtrapolationEquidistantTimeStepperIterator2, t)
    for j = 1:length(tsi.a)
        if tsi.a[j]!=0.0
            propagate_A!(tsi.psi, tsi.a[j]*tsi.dt)
        end
        if tsi.b[j]!=0.0
            gen_extrapolated_wf!(tsi.psi_ex, tsi.L[:,j], tsi.psi_back, tsi.ptr, order=tsi.order)
            set_time!(tsi.psi_ex, get_time(tsi.psi))
            if tsi.nonlinear_potential
                gen_nonlinear_potential!(tsi.rhs, tsi.psi_ex)
                to_real_space!(tsi.psi)
                to_real_space!(tsi.rhs)
                u = get_data(tsi.psi, true)
                f = get_data(tsi.rhs, true)
                u[:] .*= exp(tsi.b[j]*tsi.dt*f)
            else
                gen_rhs!(tsi.rhs, tsi.psi_ex)
                axpy!(tsi.psi, tsi.rhs, tsi.b[j]*tsi.dt)
            end
            #orthonormalize_orbitals!(tsi.psi)
        end
    end
    ptr0 = tsi.ptr
    ptr1 = mod(tsi.ptr, tsi.N) + 1
    copy!(tsi.psi_back[ptr1], tsi.psi)
    to_real_space!(tsi.psi_back[ptr1])
    for iter=1:tsi.iter
        copy!(tsi.psi, tsi.psi_back[ptr0])
        for j = 1:length(tsi.a)
            if tsi.a[j]!=0.0
                propagate_A!(tsi.psi, tsi.a[j]*tsi.dt)
            end
            if tsi.b[j]!=0.0
                gen_extrapolated_wf!(tsi.psi_ex, tsi.L1[:,j], tsi.psi_back, ptr1, order=tsi.order1)
                set_time!(tsi.psi_ex, get_time(tsi.psi))
                if tsi.nonlinear_potential
                    gen_nonlinear_potential!(tsi.rhs, tsi.psi_ex)
                    to_real_space!(tsi.psi)
                    to_real_space!(tsi.rhs)
                    u = get_data(tsi.psi, true)
                    f = get_data(tsi.rhs, true)
                    u[:] .*= exp(tsi.b[j]*tsi.dt*f)
                else
                    gen_rhs!(tsi.rhs, tsi.psi_ex)
                    axpy!(tsi.psi, tsi.rhs, tsi.b[j]*tsi.dt)
                end
                #orthonormalize_orbitals!(tsi.psi)
            end
        end
        copy!(tsi.psi_back[ptr1], tsi.psi)
        to_real_space!(tsi.psi_back[ptr1])
    end
    tsi.ptr = ptr1

    if tsi.steps<0
        t1 = t + tsi.dt < tsi.tend ? t + tsi.dt : tsi.tend
    else
        t1 = t+1
    end
    return t1, t1
end

type SplittingWithExtrapolationEquidistantTimeStepperIterator3
    psi::WaveFunction
    t0::Real
    tend::Real
    dt::Real
    steps::Int
    a::Vector{Float64}
    b::Vector{Float64}
    L::Matrix{Float64}
    L1::Matrix{Float64}
    order::Int
    order1::Int
    N::Int
    iter::Int
    ptr::Int
    rhs_back::Vector{WaveFunction} # storage for previous solution values
    psi0::WaveFunction
    rhs::WaveFunction     # storage for the righthandside of the B equation
    nonlinear_potential::Bool
end

function splitting_with_extrapolation_equidistant_time_stepper3(psi::WaveFunction,
         t0::Real, tend::Real, dt::Real, order::Int, order1::Int, a::Vector{Float64}, b::Vector{Float64}; steps::Int=-1,
         nonlinear_potential::Bool=false, psi_back=nothing, iter::Int=0)
    m = psi.m
    rhs = wave_function(m)
    psi0 = wave_function(m)
    set_propagate_time_together_with_A!(m, true)
    set_time!(psi, t0)
    N = max(order, order1)+1
    if psi_back!=nothing
        rhs_back = psi_back
        @assert length(psi_back)==N
        for k=1:N
            gen_rhs!(rhs, rhs_back[k])
            copy!(rhs_back[k], rhs)
        end
    else
        rhs_back = gen_starting_values3(psi, dt, N, a, b, rhs=rhs, nonlinear_potential=nonlinear_potential)
    end
    aa = cumsum(a)
    L = gen_interpolation_matrix(collect((-order-0.0):0.0), aa)
    L1 = gen_interpolation_matrix(collect((-order1-0.0):0.0), aa-1.0)
    ptr = N
    SplittingWithExtrapolationEquidistantTimeStepperIterator3(psi, t0, tend, dt, steps, a, b, L, L1,
                 order, order1, N, iter, ptr, rhs_back, psi0, rhs, nonlinear_potential)
end

Base.start(tsi::SplittingWithExtrapolationEquidistantTimeStepperIterator3) = (tsi.steps<0 ?
    tsi.t0 + (tsi.N-1)*tsi.dt : tsi.N-1)

Base.done(tsi::SplittingWithExtrapolationEquidistantTimeStepperIterator3, t) = (tsi.steps<0 ? t >= tsi.tend : t>=tsi.steps )

function Base.next(tsi::SplittingWithExtrapolationEquidistantTimeStepperIterator3, t)
    if tsi.iter>=1
        copy!(tsi.psi0, tsi.psi)
    end
    for j = 1:length(tsi.a)
        if tsi.a[j]!=0.0
            propagate_A!(tsi.psi, tsi.a[j]*tsi.dt)
        end
        if tsi.b[j]!=0.0
            gen_extrapolated_wf!(tsi.rhs, tsi.L[:,j], tsi.rhs_back, tsi.ptr, order=tsi.order)
            if tsi.nonlinear_potential
                to_real_space!(tsi.psi)
                to_real_space!(tsi.rhs)
                u = get_data(tsi.psi, true)
                f = get_data(tsi.rhs, true)
                u[:] .*= exp(tsi.b[j]*tsi.dt*f)
            else
                axpy!(tsi.psi, tsi.rhs, tsi.b[j]*tsi.dt)
            end
            #orthonormalize_orbitals!(tsi.psi)
        end
    end

    tsi.ptr = mod(tsi.ptr, tsi.N) + 1
    if tsi.nonlinear_potential
        gen_nonlinear_potential!(tsi.rhs_back[tsi.ptr], tsi.psi)
    else
        gen_rhs!(tsi.rhs_back[tsi.ptr], tsi.psi)
    end

    for iter=1:tsi.iter
        copy!(tsi.psi, tsi.psi0)
        for j = 1:length(tsi.a)
            if tsi.a[j]!=0.0
                propagate_A!(tsi.psi, tsi.a[j]*tsi.dt)
            end
            if tsi.b[j]!=0.0
                gen_extrapolated_wf!(tsi.rhs, tsi.L1[:,j], tsi.rhs_back, tsi.ptr, order=tsi.order1)
                if tsi.nonlinear_potential
                    gen_nonlinear_potential!(tsi.rhs, tsi.psi_ex)
                    to_real_space!(tsi.psi)
                    to_real_space!(tsi.rhs)
                    u = get_data(tsi.psi, true)
                    f = get_data(tsi.rhs, true)
                    u[:] .*= exp(tsi.b[j]*tsi.dt*f)
                else
                    axpy!(tsi.psi, tsi.rhs, tsi.b[j]*tsi.dt)
                end
                #orthonormalize_orbitals!(tsi.psi)
            end
        end
        if tsi.nonlinear_potential
            gen_nonlinear_potential!(tsi.rhs_back[tsi.ptr], tsi.psi)
        else
            gen_rhs!(tsi.rhs_back[tsi.ptr], tsi.psi)
        end
    end

    if tsi.steps<0
        t1 = t + tsi.dt < tsi.tend ? t + tsi.dt : tsi.tend
    else
        t1 = t+1
    end
    return t1, t1
end

function global_orders3(psi::WaveFunction, reference_solution::WaveFunction,
                        t0::Real, tend::Real, dt::Real, order::Int, order1::Int, a::Vector{Float64}, b::Vector{Float64};
                        rows=8, nonlinear_potential::Bool=false, iter::Int=0)
    @assert psi.m==reference_solution.m
    tab = Array(Float64, rows, 3)

    wf_save_initial_value = clone(psi)
    copy!(wf_save_initial_value, psi)

    steps = Int(floor((tend-t0)/dt))
    dt1 = dt
    err_old = 0.0
    println("             dt         err      p")
    println("-----------------------------------")
    for row=1:rows
        for t in splitting_with_extrapolation_equidistant_time_stepper3(psi, t0, tend, dt1, order, order1, a, b, steps=steps,
            nonlinear_potential=nonlinear_potential, iter=iter)
        end
        err = distance(psi, reference_solution)
        if (row==1) then
            @printf("%3i%12.3e%12.3e\n", row, Float64(dt1), Float64(err))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = 0
        else
            p = log(err_old/err)/log(2.0);
            @printf("%3i%12.3e%12.3e%7.2f\n", row, Float64(dt1), Float64(err), Float64(p))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = p
        end
        err_old = err
        dt1 = 0.5*dt1
        steps = 2*steps
        copy!(psi,wf_save_initial_value)
    end
end

function global_orders2(psi::WaveFunction, reference_solution::WaveFunction,
                        t0::Real, tend::Real, dt::Real, order::Int, order1::Int, a::Vector{Float64}, b::Vector{Float64};
                        rows=8, nonlinear_potential::Bool=false, iter::Int=0)
    @assert psi.m==reference_solution.m
    tab = Array(Float64, rows, 3)

    wf_save_initial_value = clone(psi)
    copy!(wf_save_initial_value, psi)

    steps = Int(floor((tend-t0)/dt))
    dt1 = dt
    err_old = 0.0
    println("             dt         err      p")
    println("-----------------------------------")
    for row=1:rows
        for t in splitting_with_extrapolation_equidistant_time_stepper2(psi, t0, tend, dt1, order, order1, a, b, steps=steps,
            nonlinear_potential=nonlinear_potential, iter=iter)
        end
        err = distance(psi, reference_solution)
        if (row==1) then
            @printf("%3i%12.3e%12.3e\n", row, Float64(dt1), Float64(err))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = 0
        else
            p = log(err_old/err)/log(2.0);
            @printf("%3i%12.3e%12.3e%7.2f\n", row, Float64(dt1), Float64(err), Float64(p))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = p
        end
        err_old = err
        dt1 = 0.5*dt1
        steps = 2*steps
        copy!(psi,wf_save_initial_value)
    end
end

function global_orders(psi::WaveFunction, reference_solution::WaveFunction,
                       t0::Real, tend::Real, dt::Real, order::Int, a::Vector{Float64}, b::Vector{Float64};
                        operator_sequence="AB", rows=8, nonlinear_potential::Bool=false, B_iter::Int=0)
    @assert psi.m==reference_solution.m
    tab = Array(Float64, rows, 3)

    wf_save_initial_value = clone(psi)
    copy!(wf_save_initial_value, psi)

    steps = Int(floor((tend-t0)/dt))
    dt1 = dt
    err_old = 0.0
    println("             dt         err      p")
    println("-----------------------------------")
    for row=1:rows
        for t in splitting_with_extrapolation_equidistant_time_stepper(psi, t0, tend, dt1, order, a, b, steps=steps,
            nonlinear_potential=nonlinear_potential, B_iter=B_iter)
        end
        err = distance(psi, reference_solution)
        if (row==1) then
            @printf("%3i%12.3e%12.3e\n", row, Float64(dt1), Float64(err))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = 0
        else
            p = log(err_old/err)/log(2.0);
            @printf("%3i%12.3e%12.3e%7.2f\n", row, Float64(dt1), Float64(err), Float64(p))
            tab[row,1] = dt1
            tab[row,2] = err
            tab[row,3] = p
        end
        err_old = err
        dt1 = 0.5*dt1
        steps = 2*steps
        copy!(psi,wf_save_initial_value)
    end
end

type SplittingWithExtrapolationAdaptiveTimeStepperIterator
    psi::WaveFunction
    t::Float64
    t0::Float64
    tend::Float64
    dt::Float64
    tol::Float64
    dt_old::Float64
    err_old::Float64
    order::Int
    order1::Int
    a1::Vector{Float64}
    b1::Vector{Float64}
    aa1::Vector{Float64}
    order2::Int
    a2::Vector{Float64}
    b2::Vector{Float64}
    aa2::Vector{Float64}
    N::Int
    ptr::Int
    psi_back::Vector{WaveFunction} # storage for previous solution values
    t_back::Vector{Float64}  # previous times
    psi_ex::WaveFunction  # storage for the extrapolated values
    rhs::WaveFunction     # storage for the righthandside of the B equation
    psi2::WaveFunction
    psi0::WaveFunction
    nonlinear_potential::Bool
    B_iter1::Int
    c1::Vector{Float64}
    d1::Vector{Float64}
    B_iter2::Int
    c2::Vector{Float64}
    d2::Vector{Float64}
end

function splitting_with_extrapolation_adaptive_time_stepper(psi::WaveFunction,
        t0::Float64, tend::Float64,  dt::Float64, tol::Float64,order::Int,
        order1::Int, a1::Vector{Float64}, b1::Vector{Float64},
        order2::Int, a2::Vector{Float64}, b2::Vector{Float64}; nonlinear_potential::Bool=false,
        B_iter1::Int=0, B_iter2::Int=0)
    m =psi.m
    psi_ex = wave_function(m)
    rhs = wave_function(m)
    psi0 = wave_function(m)
    psi2 = wave_function(m)

    set_propagate_time_together_with_A!(m, true)
    set_time!(psi, t0)
    N = max(order1, order2)+1
    psi_back = gen_starting_values(psi, dt, N, a1, b1, psi_ex=psi_ex, rhs=rhs, nonlinear_potential=nonlinear_potential, B_iter=B_iter1)
    t_back = [dt*k for k=0:N-1]
    t = t_back[N]
    ptr = order+1
    copy!(psi, psi_back[ptr])
    aa1 = cumsum(a1)
    bb1 = cumsum(b1)
    c1 = (aa1-[0.0; bb1[1:end-1]])./b1
    d1 = (bb1-aa1)./b1
    aa2 = cumsum(a2)
    bb2 = cumsum(b2)
    c2 = (aa2-[0.0; bb2[1:end-1]])./b2
    d2 = (bb2-aa2)./b2
    SplittingWithExtrapolationAdaptiveTimeStepperIterator(psi,t,t0,tend,dt,tol,-1.0,-1.0,order,
        order1,a1,b1,aa1,order2,a2,b2,aa2, N, ptr, psi_back, t_back, psi_ex, rhs, psi2, psi0, nonlinear_potential,
        B_iter1, c1, d1, B_iter2, c2, d2)
end

Base.start(tsi::SplittingWithExtrapolationAdaptiveTimeStepperIterator) = tsi.t

function Base.done(tsi::SplittingWithExtrapolationAdaptiveTimeStepperIterator, state_t::Float64)
   tsi.t >= tsi.tend
end

function Base.next(tsi::SplittingWithExtrapolationAdaptiveTimeStepperIterator, state_t::Float64 )
    const facmin = 0.25
    const facmax = 4.0
    const fac = 0.9
    const beta1=0.25
    const beta2=0.25
    const alpha2=0.25
    dt = tsi.dt
    dt0 = dt
    dt_old=dt
    copy!(tsi.psi0, tsi.psi)
    err = 2.0 #error/tol
    while err>=1.0
        dt = min(dt, tsi.tend-tsi.t)
        dt0 = dt
        tt = Float64[ tsi.t_back[mod(j-tsi.order1+tsi.ptr-2, tsi.N)+1] for j=1:tsi.order1+1]
        L1 = gen_interpolation_matrix(tt, tsi.t + dt*tsi.aa1)
        L2 = gen_interpolation_matrix(tt, tsi.t + dt*tsi.aa2)

        #step_embedded!(tsi.psi, tsi.psi2, dt, tsi.scheme1, tsi.scheme2, tsi.operator_sequence)
        k_start2 = -1
        start2_with_A = false

        for k = 1:length(tsi.a1)
            if k_start2<0.0 && tsi.a1[k]!=tsi.a2[k]
                k_start2 = k
                start2_with_A = true
                copy!(tsi.psi2, tsi.psi)
            end
            if tsi.a1[k]!=0.0
                propagate_A!(tsi.psi, tsi.a1[k]*tsi.dt)
            end
            if k_start2<0.0 && b1[k]!=tsi.b2[k]
                k_start2 = k
                start2_with_A = false
                copy!(tsi.psi2, tsi.psi)
            end
            if tsi.b1[k]!=0.0
                gen_extrapolated_wf!(tsi.psi_ex, L1[:,k], tsi.psi_back, tsi.ptr, order=tsi.order1)
                B_step!(tsi.psi, tsi.psi_ex, tsi.rhs, tsi.b1[j]*tsi.dt, tsi.c1[j], tsi.d1[j], nonlinear_potential=tsi.nonlinear_potential, B_iter=tsi.B_iter1)
                #set_time!(tsi.psi_ex, get_time(tsi.psi))
                #if tsi.nonlinear_potential
                #    gen_nonlinear_potential!(tsi.rhs, tsi.psi_ex)
                #    to_real_space!(tsi.psi)
                #    to_real_space!(tsi.rhs)
                #    u = get_data(tsi.psi, true)
                #    f = get_data(tsi.rhs, true)
                #    u[:] .*= exp(tsi.b1[k]*tsi.dt*f)
                #else
                #    gen_rhs!(tsi.rhs, tsi.psi_ex)
                #    axpy!(tsi.psi, tsi.rhs, tsi.b1[k]*tsi.dt)
                #end
            end
        end

        if k_start2<0
            k_start2 = length(tsi.a1)+1
            start2_with_A = true
            copy!(tsi.psi2, tsi.psi)
        end

        for k = k_start2:length(tsi.a2)
            if tsi.a2[k]!=0.0 &&(k!=k_start2||start2_with_A)
                propagate_A!(tsi.psi2, tsi.a2[k]*tsi.dt)
            end
            if tsi.b2[k]!=0.0
                gen_extrapolated_wf!(tsi.psi_ex, L2[:,k], tsi.psi_back, tsi.ptr, order=tsi.order1)
                B_step!(tsi.psi2, tsi.psi_ex, tsi.rhs, tsi.br2[j]*tsi.dt, tsi.c2[j], tsi.d2[j], nonlinear_potential=tsi.nonlinear_potential, B_iter=tsi.B_iter2)
                #set_time!(tsi.psi_ex, get_time(tsi.psi2))
                #if tsi.nonlinear_potential
                #    gen_nonlinear_potential!(tsi.rhs, tsi.psi_ex)
                #    to_real_space!(tsi.psi2)
                #    to_real_space!(tsi.rhs)
                #    u = get_data(tsi.psi2, true)
                #    f = get_data(tsi.rhs, true)
                #    u[:] .*= exp(tsi.b2[k]*tsi.dt*f)
                #else
                #    gen_rhs!(tsi.rhs, tsi.psi_ex)
                #    axpy!(tsi.psi2, tsi.rhs, tsi.b2[k]*tsi.dt)
                #end
            end
        end

        err = distance(tsi.psi, tsi.psi2)/tsi.tol
        if tsi.dt_old<0.0
           dt = dt*min(facmax, max(facmin, fac*(1.0/err)^(1.0/(tsi.order+1))))
        else
            dt = dt*min(facmax, max(facmin, fac*((1.0/err)^(beta1/(tsi.order+1))*(1.0/tsi.err_old)^(beta2/(tsi.order+1))*(dt/tsi.dt_old)^(-alpha2))))
        end
        if err>=1.0
           copy!(tsi.psi, tsi.psi0)
            @printf("t=%17.9e  err=%17.8e  dt=%17.8e  rejected...\n", Float64(tsi.t), Float64(err), Float64(dt))
        end
    end

    tsi.t += dt0
    tsi.dt = dt
    tsi.dt_old = dt_old
    tsi.err_old = err

    tsi.ptr = mod(tsi.ptr, tsi.N) + 1
    copy!(tsi.psi_back[tsi.ptr], tsi.psi)
    to_real_space!(tsi.psi_back[tsi.ptr])
    tsi.t_back[tsi.ptr] = tsi.t

    tsi.t, tsi.t
end


text:
module Charts

import DataFrames
import Genie, Stipple
import Genie.Renderer.Html: HTMLString, normal_element

using Stipple

export PlotOptions, PlotData, PlotSeries, plot

Genie.Renderer.Html.register_normal_element("apexchart", context = @__MODULE__)

Base.@kwdef mutable struct PlotOptions
  chart_animations_enabled::Bool = false
  chart_animations_easing::Symbol = :easeinout
  chart_animations_speed::Int = 800
  chart_background::String = "#fff"
  chart_font_family::String = "Helvetica, Arial, sans-serif"
  chart_fore_color::String = "#373d3f"
  chart_height::Union{Int,String} = "auto"
  chart_offset_x::Int = 0
  chart_offset_y::Int = 0
  chart_selection_enabled::Bool = true
  chart_selection_type::Symbol = :xy
  chart_sparkline_enabled::Bool = false
  chart_stacked::Bool = false
  chart_stack_type::String = "normal"
  chart_toolbar_show::Bool = false
  chart_type::Symbol = :line
  chart_width::Union{Int,String} = "100%"
  chart_zoom_enabled::Bool = true
  chart_zoom_type::Symbol = :xy
  colors::Vector{String} = ["#2E93fA", "#66DA26", "#546E7A", "#E91E63", "#FF9800"]
  data_labels_enabled::Bool = false
  fill_opacity::Float64 = 1.0
  grid_show::Bool = true
  labels::Vector{String} = String[]
  legend_show::Bool = true
  legend_position::Symbol = :bottom
  legend_font_size::String = "14px"
  legend_font_family::String = "Helvetica, Arial"
  no_data_text::String = ""
  plot_options_bar_horizontal::Bool = false
  plot_options_bar_ending_shape::Symbol = :flat
  plot_options_bar_column_width::String = "100%"
  plot_options_bar_data_labels_position::Symbol = :center
  plot_options_bubble_min_bubble_radius::Union{Int,Symbol} = :undefined
  plot_options_bubble_max_bubble_radius::Union{Int,Symbol} = :undefined
  stroke_curve::Symbol = :smooth
  stroke_show::Bool = true
  stroke_width::Int = 2
  stroke_colors::Vector{String} = String["transparent"]
  subtitle_text::String = ""
  subtitle_align::Symbol = :left
  subtitle_style_font_size::String = "12px"
  theme_mode::Symbol = :light
  theme_palette::Symbol = :palette1
  title_text::String = ""
  title_align::Symbol = :left
  title_margin::Int = 10
  title_style_font_size::String = "14px"
  title_style_font_weight::Union{Int,Symbol} = :bold
  title_style_color::String = "#263238"
  tooltip_enable::Bool = true
  xaxis_type::Symbol = :category
  xaxis_categories::Vector{String} = String[]
  xaxis_tick_amount::Union{Int,Float64,Symbol} = :undefined
  xaxis_max::Union{Int,Float64,Symbol} = :undefined
  xaxis_min::Union{Int,Float64,Symbol} = :undefined
  xaxis_labels_show::Bool = true
  yaxis_tick_amount::Union{Int,Float64,Symbol} = :undefined
  yaxis_max::Union{Int,Float64,Symbol} = :undefined
  yaxis_min::Union{Int,Float64,Symbol} = :undefined
  yaxis_labels_show::Bool = true
end

Base.@kwdef mutable struct PlotData{T<:Vector}
  data::T = T[]
  key::Union{Symbol,Nothing} = :data
end

function PlotData(data::DataFrames.DataFrame; key=:data)
  PlotData(key = key, data = [Array(r) for r in DataFrames.eachrow(data)])
end

function PlotData(data::Vector{T}) where {T}
  PlotData(data = data)
end

Base.@kwdef mutable struct PlotSeries
  name::String = ""
  plotdata::PlotData = PlotData()
end

function plot(fieldname::Symbol;
              options::Union{Symbol,Nothing} = nothing,
              args...) :: String

  k = (Symbol(":series"),)
  v = Any["$fieldname"]

  if options !== nothing
    k = (k..., Symbol("options!"))
    push!(v, options)
  end

  Genie.Renderer.Html.div() do
    apexchart(; args..., NamedTuple{k}(v)...)
  end
end

#===#

function Base.parse(::Type{Vector{PlotSeries}}, d::Vector{Any})
  [PlotSeries(name = x["name"], data = x["data"]) for x in d]
end

function Base.parse(::Type{PlotSeries}, x::Dict{String,Any})
  PlotSeries(name = x["name"], data = x["data"])
end

function Base.parse(::Type{PlotData}, x::Vector{Any})
  PlotData(x)
end

#===#

function Stipple.watch(vue_app_name::String, fieldtype::R{PlotSeries}, fieldname::Symbol, channel::String, model::M)::String where {M<:ReactiveModel}
  string(vue_app_name, raw".\$watch('", fieldname, "', function(newVal, oldVal){

  });\n\n")
end

#===#

function Stipple.render(pd::PlotData{T}, fieldname::Union{Symbol,Nothing} = nothing) where {T<:Vector}
  pd.key === nothing ? pd.data : [Dict(pd.key => pd.data)]
end

function Stipple.render(pdv::Vector{PlotData{T}}, fieldname::Union{Symbol,Nothing} = nothing) where {T<:Vector}
  [Dict(pd.key => pd.data) for pd in pdv]
end

function Stipple.render(ps::PlotSeries, fieldname::Union{Symbol,Nothing} = nothing)
  Dict(:name => ps.name, ps.plotdata.key => ps.plotdata.data)
end

function Stipple.render(psv::Vector{PlotSeries}, fieldname::Union{Symbol,Nothing} = nothing)
  [Dict(:name => ps.name, ps.plotdata.key => ps.plotdata.data) for ps in psv]
end

function Stipple.render(po::PlotOptions, fieldname::Union{Symbol,Nothing} = nothing)
  val = Dict(
    :chart => Dict(
      :animations => Dict(
        :enabled => po.chart_animations_enabled,
        :easing => po.chart_animations_easing,
        :speed => po.chart_animations_speed
      ),
      :background => po.chart_background,
      :fontFamily => po.chart_font_family,
      :foreColor => po.chart_fore_color,
      :height => po.chart_height,
      :offsetX => po.chart_offset_x,
      :offsetY => po.chart_offset_y,
      :selection => Dict(
        :enabled => po.chart_selection_enabled,
        :type => po.chart_selection_type
      ),
      :sparkline => Dict(
        :enabled => po.chart_sparkline_enabled
      ),
      :stacked => po.chart_stacked,
      :stackType => po.chart_stack_type,
      :toolbar => Dict(
        :show => po.chart_toolbar_show
      ),
      :type => po.chart_type,
      :width => po.chart_width,
      :zoom => Dict(
        :enabled => po.chart_zoom_enabled,
        :type => po.chart_zoom_type
      )
    ),
    :colors => po.colors,
    :dataLabels => Dict(
      :enabled => po.data_labels_enabled
    ),
    :fill => Dict(
      :opacity => po.fill_opacity
    ),
    :grid => Dict(
      :show => po.grid_show
    ),
    :labels => po.labels,
    :legend => Dict(
      :position => po.legend_position,
      :fontFamily => po.legend_font_family,
      :fontSize => po.legend_font_size,
      :show => po.legend_show
    ),
    :noData => Dict(
      :text => po.no_data_text
    ),
    :plotOptions => Dict(
      :bar => Dict(
        :horizontal => po.plot_options_bar_horizontal,
        :endingShape => po.plot_options_bar_ending_shape,
        :columnWidth => po.plot_options_bar_column_width,
        :dataLabels => Dict(
          :position => po.plot_options_bar_data_labels_position
        )
      ),
      :bubble => Dict(
        :minBubbleRadius => po.plot_options_bubble_min_bubble_radius,
        :maxBubbleRadius => po.plot_options_bubble_max_bubble_radius
      )
    ),
    :stroke => Dict(
      :curve => po.stroke_curve,
      :show => po.stroke_show,
      :width => po.stroke_width,
      :colors => po.stroke_colors
    ),
    :subtitle => Dict(
      :align => po.subtitle_align,
      :text => po.subtitle_text,
      :style => Dict(
        :fontSize => po.subtitle_style_font_size
      )
    ),
    :theme => Dict(
      :mode => po.theme_mode,
      :palette => po.theme_palette
    ),
    :title => Dict(
      :text => po.title_text,
      :align => po.title_align,
      :margin => po.title_margin,
      :style => Dict(
        :color => po.title_style_color,
        :fontSize => po.title_style_font_size,
        :fontWeight => po.title_style_font_weight
      )
    ),
    :tooltip => Dict(
      :enable => po.tooltip_enable
    ),
    :xaxis => Dict(
      :categories => po.xaxis_categories,
      :tickAmount => po.xaxis_tick_amount,
      :type => po.xaxis_type,
      :max => po.xaxis_max,
      :min => po.xaxis_min,
      :labels => Dict(
        :show => po.xaxis_labels_show
      )
    ),
    :yaxis => Dict(
      :tickAmount => po.yaxis_tick_amount,
      :max => po.yaxis_max,
      :min => po.yaxis_min,
      :labels => Dict(
        :show => po.yaxis_labels_show
      )
    )
  )

  replace(Genie.Renderer.Json.JSONParser.json(val), "\"undefined\""=>"undefined")
end

# #===#

end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 458f9f6e-2c17-11ec-08fe-d92011537eff
begin
	using PlutoUI
end

# ╔═╡ 56ff08ee-dcba-4aae-971a-5eec10a5d8c2
md"""
#### Solution: Coffee Maker Problem
"""

# ╔═╡ 28054bbb-cf38-403f-a563-c8a79551c950
md"""$(PlutoUI.LocalResource("./figs/Fig-CoffeeMaker-Problem-P1.png"))"""

# ╔═╡ 3e6a3f82-70fd-4692-8779-9d43f41947a2
begin

	# setup problem -
	m_3_dot = 0.5 	# units: kg/min
	m_4_dot = 0.36 	# units: kg/min
	α = 0.25 	 	# units: dimensionless
	H₁ = 104.5 		# units: kJ/kg
	H₂ = 283.5 		# units: kJ/kg
	H₃ = 168.5 		# units: kJ/kg
	H₄ = 722.4 		# units: kJ/kg
	m_coffee = 0.5 	# units: kg

	# show -
	nothing
end

# ╔═╡ cfad92dc-25e1-4d45-b94d-9c83888e3695
md"""
a) Determine the flow rate of coffee beans ($\dot{m}_{2}$) and water ($\dot{m}_{1}$) required to make a good cup of coffee (requirement: $\dot{m}_2$ = 1/4$\dot{m}_{1}$).

Start by applying the total mass balance at steady state:

$$\sum_{s=1}^{\mathcal{S}}v_{s}\dot{m}_{s} = 0$$

to the mixer grinder unit:

$$\dot{m}_{1}+\dot{m}_{2} = \dot{m}_{3}$$

In addition to the total mass balance, we also know that $\dot{m}_2$ = 1/4$\dot{m}_{1}$, thus we have two equations and two unknowns which we can solve for $\dot{m}_{1}$ and $\dot{m}_{2}$. You can do this via subsutitution, or by setting a a matrix inversion.

_Substitution approach_: Substitute $\dot{m}_{2}$ = $\alpha\dot{m}_{1}$ into the total mass balance and solve for $\dot{m}_{1}$:

$$\dot{m}_{1} = \frac{\dot{m}_{3}}{1+\alpha}$$

_Matrix vector approach_: We can solve for $\dot{m}_{1}$ and $\dot{m}_{2}$ by setting up the mass balance and coffee requirement as a 2$\times$2 system of equations:

$$\begin{pmatrix}
1 & 1 \\
-\alpha & 1
\end{pmatrix}
\begin{pmatrix}
\dot{m}_{1} \\
\dot{m}_{2}
\end{pmatrix} =
\begin{pmatrix}
\dot{m}_{3} \\
0
\end{pmatrix}$$

which is in the form $\mathbf{A}\mathbf{x}$=$\mathbf{b}$. We can compute the matrix
inverse $\mathbf{A}^{-1}$ to estimate $\mathbf{x}.

"""

# ╔═╡ 36a54cd9-d496-4066-872f-311fb1eba987
begin

	# subsitution method -
	m_1_dot = (m_3_dot)/(1+α)
	m_2_dot = α * m_1_dot

	# show -
	with_terminal() do
		println("m_1_dot = $(m_1_dot) kg/min and m_2_dot = $(m_2_dot) kg/min")
	end
end

# ╔═╡ dbe9f642-feca-4d4e-8cb9-1c851066863d
begin

	# matrix vector solution for m1 and m2 -
	A = [1 1 ; -α 1];
	b = [m_3_dot ; 0];

	# compute -
	x = inv(A)*b
end

# ╔═╡ ea96f2e5-1e4d-43e8-a9ef-a5f0af703349
md"""
b) Because he can’t plug the coffee maker in, he has to generate power to mix/grind the beans and water. What is the rate of work $\dot{W}_{sh}$ that he must supply to the mixer/grinder if this unit is well insulated?

To estimate the rate of shaft work $\dot{W}_{sh}$, we need to use the open energy balance equation at steady-state:

$$\dot{Q}+\dot{W}_{sh}+\sum_{s=1}^{\mathcal{S}}v_{s}\dot{m}_{s}H_{s} = 0$$

We know from the problem that $\dot{Q}$ = 0, thus we can solve the open energy balance to give an expression for the shaft work:

$$\dot{W}_{sh} = \dot{m}_{3}H_{3} - \left(\dot{m}_{1}H_{1} + \dot{m}_{2}H_{2}\right)$$

You already solved for the mass flow rates $\dot{m}_{i}$, and the enthalpy of each stream is given in the problem, thus, we know everything we need to solve for the work required to grind the beans.

"""

# ╔═╡ a3882553-e45f-49a9-8399-9adaaf02a9e7
begin

	W_dot_sh = m_3_dot*H₃ - (m_1_dot*H₁ + m_2_dot*H₂)

	with_terminal() do
		println("Rate of work for grinder = $(W_dot_sh) kJ/min")
	end

end

# ╔═╡ ff263509-d5a8-4227-898b-67491b99bf51
md"""
c) Prof. V loves high enthalpy coffee! In fact, he refuses to drink coffee unless the enthalpy content is between 475 kJ/kg $\leq{H_{5}}\leq$ 500 kJ/kg. What is the minimum and maximum rate of heat input to the brewer/separator that must be supplied to meet this completely reasonable criteria?

To estimate the heat input into the brewer/separator, we need to estimate $\dot{m}_{5}$ (using a total mass balance around the brewer) and then solve for $\dot{Q}$ of the brewer unit using an open energy balance.

Again, we could do solve using algebra and substitution (first find $\dot{m}_{5}$ using the total mass balance, then find $\dot{Q}$ using the open energy balance). However, I'm lazy (and I'm running out of time to post these solutions!) so I'm going to setup the solution as a matrix vector problem, and then invert the matrix:

$$\begin{pmatrix}
-1 & 0 \\
-H_{5} & 1
\end{pmatrix}
\begin{pmatrix}
\dot{m}_{5}\\
\dot{Q}
\end{pmatrix} =
\begin{pmatrix}
\dot{m}_{4} - \dot{m}_{3} \\
\dot{m}_{4}H_{4} - \dot{m}_{3}H_{3}
\end{pmatrix}$$

We can solve the problem twice: once for $H_{5}$ = 475 kJ/kg and again for $H_{5}$ = 500 kJ/kg, this will give us the $\dot{Q}$ range (I'll show the calcultion for the lower bound).

"""

# ╔═╡ 6fff5dba-e6e0-42e8-a06d-519fa582695e
begin

	# Setup -
	H₅ = 475.0 	# units: kJ/kg

	# setup matrix vector form Ax = b, the find x = A^{-1}b
	Ac = [-1 0; -H₅ 1];
	bc = [m_4_dot - m_3_dot ; m_4_dot*H₄ - m_3_dot*H₃];

	# solve -
	xc = inv(Ac)*bc

	# split -
	m_5_dot = xc[1]
	Qdot = xc[2]

	# show -
	with_terminal() do
		println("m_5_dot = $(m_5_dot) kg/min and Q_dot = $(Qdot) kJ/min")
	end
end

# ╔═╡ 2387358c-ee90-4594-9023-b03a51867745
md"""
d) Prof. V has five minutes to make coffee before class starts, did he make it in time (1 cup of coffee is $\simeq$ 0.5 kg)?

To estimate the time, we need to write a total mass balance around the coffee cup, however, the coffee cup will __not__ be at steady-state, this we have:

$$\frac{dm}{dt} = \dot{m}_{5}$$

where $m$ denotes the mass of coffee in the cup, and $t$ denotes time. Let's treat the derivative as an infinitesimal, and rearrange the mass balance as:

$$dm = \dot{m}_{5}dt$$

Next, we can integrate both sides (given that $\dot{m}_{5}$ is constant) between an initial point $t_{1}$ (when the cup is empty, define $t_{1}$ = 0) and someother point in time $t_{2}$:

$$m_{2} = m_{1} + \dot{m}_{5}\left(t_{2} - t_{1}\right)$$

However, at $t_{1}$ = 0, we know that the cup is empty, so $m_{1}$ = 0. Thus, we can rearrange and solve for $t_{2}$:

$$t_{2} = \frac{m_{2}}{\dot{m}_{5}}$$

"""

# ╔═╡ cf7f585a-29a0-4df6-a2a1-dd40dd8a9cf0
begin

	t₂ = m_coffee/(m_5_dot) # units: minutes

	with_terminal() do
		println("Time for coffee to fill cup: $(t₂) min")
	end
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.16"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[HypertextLiteral]]
git-tree-sha1 = "f6532909bf3d40b308a0f360b6a0e626c0e263a8"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.1"

[[IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "a8709b968a1ea6abc2dc1967cb1db6ac9a00dfb6"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.5"

[[PlutoUI]]
deps = ["Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "4c8a7d080daca18545c56f1cac28710c362478f3"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.16"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╟─56ff08ee-dcba-4aae-971a-5eec10a5d8c2
# ╟─28054bbb-cf38-403f-a563-c8a79551c950
# ╠═3e6a3f82-70fd-4692-8779-9d43f41947a2
# ╟─cfad92dc-25e1-4d45-b94d-9c83888e3695
# ╠═36a54cd9-d496-4066-872f-311fb1eba987
# ╠═dbe9f642-feca-4d4e-8cb9-1c851066863d
# ╟─ea96f2e5-1e4d-43e8-a9ef-a5f0af703349
# ╠═a3882553-e45f-49a9-8399-9adaaf02a9e7
# ╟─ff263509-d5a8-4227-898b-67491b99bf51
# ╠═6fff5dba-e6e0-42e8-a06d-519fa582695e
# ╟─2387358c-ee90-4594-9023-b03a51867745
# ╠═cf7f585a-29a0-4df6-a2a1-dd40dd8a9cf0
# ╠═458f9f6e-2c17-11ec-08fe-d92011537eff
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
function isPalindrome(n)
    s = string(n)
    s == reverse(s)
end

max = 0
for i in collect(999:-1:100)
    for j in collect(999:-1:100)
        n = i*j
        if (isPalindrome(n) && n > max)
            max = n
        end
    end
end

println(max)

max_stars_repo_path: 004 - Largest Palindrome Product/004 - Largest_Palindrome_Product.jl

max_stars_repo_name: LanguageAdventures/Project-Euler-Julia

max_issues_repo_name: LanguageAdventures/Project-Euler-Julia

max_forks_repo_name: LanguageAdventures/Project-Euler-Julia


text:
module FiniteStrains

using ..EquationsOfStateOfSolids: _⅔, _⅓, _1½

abstract type FiniteStrain end
struct EulerianStrain <: FiniteStrain end
struct LagrangianStrain <: FiniteStrain end
struct NaturalStrain <: FiniteStrain end
struct InfinitesimalStrain <: FiniteStrain end

"""
    ToEulerianStrain(v0)(v)
    ToLagrangianStrain(v0)(v)
    ToNaturalStrain(v0)(v)
    ToInfinitesimalStrain(v0)(v)

Calculate the finite strain of `v` based on the reference volume `v0`.

!!! info
    See the formulae on the [`Gibbs2` paper](https://www.sciencedirect.com/science/article/pii/S0010465511001470) Table 3.
"""
struct ToStrain{S<:FiniteStrain,T}
    v0::T
end
ToStrain{S}(v0::T) where {S,T} = ToStrain{S,T}(v0)
const ToEulerianStrain = ToStrain{EulerianStrain}
const ToLagrangianStrain = ToStrain{LagrangianStrain}
const ToNaturalStrain = ToStrain{NaturalStrain}
const ToInfinitesimalStrain = ToStrain{InfinitesimalStrain}
(x::ToEulerianStrain)(v) = ((x.v0 / v)^_⅔ - 1) / 2
(x::ToLagrangianStrain)(v) = ((v / x.v0)^_⅔ - 1) / 2
(x::ToNaturalStrain)(v) = log(v / x.v0) / 3
(x::ToInfinitesimalStrain)(v) = 1 - (x.v0 / v)^_⅓

"""
    FromEulerianStrain(v0)(f)
    FromLagrangianStrain(v0)(f)
    FromNaturalStrain(v0)(f)
    FromInfinitesimalStrain(v0)(f)

Calculate the original volume `v` from the finite strain `f` based on the reference volume `v0`.

!!! info
    See the formulae on the [`Gibbs2` paper](https://www.sciencedirect.com/science/article/pii/S0010465511001470) Table 3.
"""
struct FromStrain{S<:FiniteStrain,T}
    v0::T
end
FromStrain{S}(v0::T) where {S,T} = FromStrain{S,T}(v0)
const FromEulerianStrain = FromStrain{EulerianStrain}
const FromLagrangianStrain = FromStrain{LagrangianStrain}
const FromNaturalStrain = FromStrain{NaturalStrain}
const FromInfinitesimalStrain = FromStrain{InfinitesimalStrain}
# Eulerian strain
function (x::FromEulerianStrain)(f)
    v = x.v0 / (2f + 1)^_1½
    return isreal(v) ? real(v) : v
end
# Lagrangian strain
function (x::FromLagrangianStrain)(f)
    v = x.v0 * (2f + 1)^_1½
    return isreal(v) ? real(v) : v
end
# Natural strain
function (x::FromNaturalStrain)(f)
    v = x.v0 * exp(3f)
    return isreal(v) ? real(v) : v
end
# Infinitesimal strain
function (x::FromInfinitesimalStrain)(f)
    v = x.v0 / (1 - f)^3
    return isreal(v) ? real(v) : v
end

"""
    Dⁿᵥf(s::EulerianStrain, deg, v0)
    Dⁿᵥf(s::LagrangianStrain, deg, v0)
    Dⁿᵥf(s::NaturalStrain, deg, v0)
    Dⁿᵥf(s::InfinitesimalStrain, deg, v0)

Return a function of `v` that calculates the `deg`th order derivative of strain wrt volume from `v0`.

!!! info
    See the formulae on Ref. 1 Table 3.
"""
function Dⁿᵥf(s::EulerianStrain, deg, v0)
    function (v)
        if isone(deg)  # Stop recursion
            return -(v0 / v)^_⅔ / 3 / v
        else  # Recursion
            return -(3deg - 1) / 3 / v * Dⁿᵥf(s, deg - 1, v0)(v)
        end
    end
end
function Dⁿᵥf(s::LagrangianStrain, deg, v0)
    function (v)
        if isone(deg)  # Stop recursion
            return -(v / v0)^_⅔ / 3 / v
        else  # Recursion
            return -(3deg - 5) / 3 / v * Dⁿᵥf(s, deg - 1, v0)(v)
        end
    end
end
function Dⁿᵥf(s::NaturalStrain, deg, v0)
    function (v)
        if isone(deg)  # Stop recursion
            return 1 / 3 / v
        else  # Recursion
            return -(deg - 1) / v * Dⁿᵥf(s, deg - 1, v0)(v)
        end
    end
end
function Dⁿᵥf(s::InfinitesimalStrain, deg, v0)
    function (v)
        if isone(deg)  # Stop recursion
            return (1 - ToInfinitesimalStrain(v0)(v))^4 / 3 / v0
        else  # Recursion
            return -(3deg - 2) / 3 / v * Dⁿᵥf(s, deg - 1, v0)(v)
        end
    end
end

function straintype end

end


text:
module Test3D

using Test
using Trixi2Vtk

include("test_trixi2vtk.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

# Create artifacts directory where all files that should be preserved will be stored
artifacts_dir = joinpath(pathof(Trixi2Vtk) |> dirname |> dirname, "artifacts")
if !isdir(artifacts_dir)
  mkdir(artifacts_dir)
end

@testset "3D" begin
  @testset "TreeMesh" begin
    isdir(outdir) && rm(outdir, recursive=true)
    run_trixi(joinpath("tree_3d_dgsem", "elixir_advection_extended.jl"), maxiters=1)

    @timed_testset "uniform mesh" begin
      test_trixi2vtk("solution_000000.h5", outdir,
          hashes=[("solution_000000.vtu", "6ab3aa525851187ee0839e1d670a254a66be4ad7"),
                  ("solution_000000_celldata.vtu", "fdfee2d4200ecdad08067b37908412813016f4e7")])

      # Store output files as artifacts to facilitate debugging of failing tests
      outfiles = ("solution_000000.vtu", "solution_000000_celldata.vtu")
      testname = "3d-tree-mesh-uniform-mesh"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end
  end

  @testset "StructuredMesh" begin
    isdir(outdir) && rm(outdir, recursive=true)
    run_trixi(joinpath("structured_3d_dgsem", "elixir_advection_basic.jl"), maxiters=1)

    @timed_testset "basic" begin
      if Sys.isapple()
        # This file has a different hash on macOS for some reason
        test_trixi2vtk("solution_000000.h5", outdir,
            hashes=[("solution_000000.vtu", "58e07f981fd6c005ea17e47054bd509c2c66d771")])
      else
        test_trixi2vtk("solution_000000.h5", outdir,
            hashes=[("solution_000000.vtu", "58e07f981fd6c005ea17e47054bd509c2c66d771")])
      end

      # Store output files as artifacts to facilitate debugging of failing tests
      outfiles = ("solution_000000.vtu",)
      testname = "3d-curved-mesh-basic"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end
  end

  @testset "P4estMesh" begin
    isdir(outdir) && rm(outdir, recursive=true)
    run_trixi(joinpath("p4est_3d_dgsem", "elixir_advection_amr_unstructured_curved.jl"), maxiters=1)

    @timed_testset "unstructured curved" begin
      if Sys.isapple()
        # This file has a different hash on macOS for some reason
        test_trixi2vtk("solution_000000.h5", outdir,
          hashes=[("solution_000000.vtu", "fe0f45809ef6f3f0c1b7f5331198585f406923c9")])
      else
        test_trixi2vtk("solution_000000.h5", outdir,
          hashes=[("solution_000000.vtu", "0fa5a099378d153aa3a1bb7dcf3559ea5d6bf9c5")])
      end

      # Store output files as artifacts to facilitate debugging of failing tests
      outfiles = ("solution_000000.vtu",)
      testname = "3d-p4est-mesh-unstructured-curved"
      for outfile in outfiles
        println("Copying '", abspath(joinpath(outdir, outfile)),
                "' to '", abspath(joinpath(artifacts_dir, testname * "-" * outfile)),
                "'...")
        cp(joinpath(outdir, outfile), joinpath(artifacts_dir, testname * "-" * outfile), force=true)
      end
    end
  end
end

# Clean up afterwards: delete Trixi output directory
@test_skip rm(outdir, recursive=true)

end


text:
"""
```
fix_upper_boundary(x::T,side) where T
```

Internal function or structure - interface may change.

# Extended help

Move `x` to `x -side` if `x == side`, because we use the convention
that the boundary belongs to the next cell.

"""
@inline fix_upper_boundary(x::T,side) where T = ifelse(x == side, zero(T), x)

"""

```
wrap_cell_fraction(x,unit_cell_matrix)
```

Internal function or structure - interface may change.

# Extended help

Obtaint the coordinates of `x` as a fraction of unit cell vectors, first
positive cell. `x` is a vector of dimension `N` and `cell` a matrix of
dimension `NxN`

## Example

```julia-repl
julia> unit_cell_matrix = [ 10 0
                            0 10 ];

julia> x = [ 15, 13 ];

julia> wrap_cell_fraction(x,unit_cell_matrix)
2-element Vector{Float64}:
 0.5
 0.3
```

"""
@inline function wrap_cell_fraction(x,unit_cell_matrix)
    p = mod.(unit_cell_matrix\x,1)
    p = fix_upper_boundary.(p,1)
    return p
end

"""

```
wrap_to_first(x,unit_cell_matrix)
```

Internal function or structure - interface may change.

# Extended help

Wraps the coordinates of point `x` such that the returning coordinates are in the
first unit cell with all-positive coordinates.

## Example

```julia-repl
julia> unit_cell_matrix = [ 10 0
                            0 10 ];

julia> x = [ 15, 13 ];

julia> wrap_to_first(x,unit_cell_matrix)
2-element Vector{Float64}:
 5.0
 3.0000000000000004
```

"""
@inline function wrap_to_first(x,unit_cell_matrix)
    p = wrap_cell_fraction(x,unit_cell_matrix)
    return unit_cell_matrix*p
end

"""

```
wrap_to_first(x,box::Box)
```

Internal function or structure - interface may change.

# Extended help

Wraps the coordinates of point `x` such that the returning coordinates are in the
first unit cell with all-positive coordinates, given the `Box` structure.

"""
@inline wrap_to_first(x,box::Box) = wrap_to_first(x,box.unit_cell.matrix)

"""

```
wrap_to_first(x,box::Box{OrthorhombicCell,N,T}) where {N,T}
```

Internal function or structure - interface may change.

# Extended help

Wraps the coordinates of point `x` such that the returning coordinates are in the
first unit cell with all-positive coordinates, given an Orthorhombic cell.
This is slightly cheaper than for general cells.

"""
@inline function wrap_to_first(x,box::Box{OrthorhombicCell,N,T}) where {N,T}
    sides = SVector{N,T}(ntuple(i->box.unit_cell.matrix[i,i],N))
    x = mod.(x,sides)
    x = fix_upper_boundary.(x,sides)
    return x
end

"""

```
wrap_relative_to(x, xref, unit_cell_matrix::SMatrix{N,N,T}) where {N,T}
```

Internal function or structure - interface may change.

# Extended help

Wraps the coordinates of point `x` such that it is the minimum image relative to `xref`.

"""
@inline function wrap_relative_to(x, xref, unit_cell_matrix::SMatrix{N,N,T}) where {N,T}
    x_f = wrap_cell_fraction(x,unit_cell_matrix)
    xref_f = wrap_cell_fraction(xref,unit_cell_matrix)
    xw = wrap_relative_to(x_f,xref_f,SVector{N,T}(ntuple(i->1,N)))
    return unit_cell_matrix * (xw - xref_f) + xref
end

"""

```
wrap_relative_to(x,xref,box::Box{UnitCellType,N,T}) where {UnitCellType,N,T}
```

Internal function or structure - interface may change.

# Extended help

Wraps the coordinates of point `x` such that it is the minimum image relative to `xref`,
given a general `Box` structure.

"""
@inline wrap_relative_to(x,xref,box::Box{UnitCellType,N,T}) where {UnitCellType,N,T} =
    wrap_relative_to(x,xref,box.unit_cell.matrix)

"""

```
wrap_relative_to(x,xref,box::Box{OrthorhombicCell,N,T}) where {N,T}
```

Internal function or structure - interface may change.

# Extended help

Wraps the coordinates of point `x` such that it is the minimum image relative to `xref`,
given an Orthorhombic cell. This is slightly cheaper than for general cells.

"""
@inline function wrap_relative_to(x,xref,box::Box{OrthorhombicCell,N,T}) where {N,T}
    sides = SVector{N,T}(ntuple(i->box.unit_cell.matrix[i,i],N))
    return wrap_relative_to(x,xref,sides)
end

"""

```
wrap_relative_to(x,xref,sides::AbstractVector)
```

Internal function or structure - interface may change.

# Extended help

Wraps the coordinates of point `x` such that it is the minimum image relative to `xref`,
for an Orthorhombic cell of which only the side lengths are provided.

"""
@inline function wrap_relative_to(x,xref,sides::AbstractVector)
    xw = mod.(x-xref,sides)
    xw = _wrap_single_coordinate.(xw,sides)
    return xw + xref
end

#
# Wrap a single coordinate
#
@inline function _wrap_single_coordinate(x,s)
    if x >= s/2
        x = x - s
    elseif x < -s/2
        x = x + s
    end
    return x
end

"""

```
translation_image(x::SVector{N,T},unit_cell_matrix,indices) where {N,T}
```

Internal function or structure - interface may change.

# Extended help

Translate vector `x` according to the `unit_cell_matrix` lattice vectors and the `indices`
provided.

"""
@inline translation_image(x::SVector{N,T},unit_cell_matrix,indices) where {N,T} =
    x + unit_cell_matrix*SVector{N,Int}(ntuple(i -> indices[i],N))

"""
```
translation_image(x::AbstractVector{<:AbstractVector},unit_cell_matrix,indices)
```

Translates a complete set of coordinates given a set of indexes of unit-cells. Returns a new
set of coordinates.

## Example

```julia-repl
julia> x = rand(SVector{2,Float64},100);

julia> box = Box([1,1],0.1);

julia> CellListMap.translation_image(x,box.unit_cell.matrix,(1,1))
100-element Vector{SVector{2, Float64}}:
 [1.847791110439223, 1.5989103939725295]
 [1.3493293666090889, 1.4002971843576644]
 [1.4111736701313218, 1.3471780214994182]
 ⋮
 [1.1548437388991908, 1.7034501001177493]
 [1.4066300885242247, 1.2907398318754952]
```

"""
function translation_image(x::AbstractVector{<:AbstractVector},unit_cell_matrix,indices)
    x_new = similar(x)
    for i in eachindex(x)
        x_new[i] = translation_image(x[i],unit_cell_matrix,indices)
    end
    return x_new
end

"""

```
replicate_system!(x::AbstractVector,box::Box,ranges::Tuple)
```

Internal function or structure - interface may change.

# Extended help

Replicate the system (modifying the original array of coordinates) in all
directions defined by the periodic system and by the range of unitary cells
of interest. `x` can be a `(N,M)` matrix, and the unit cell matrix can be
provided instead of the `box`.

## Example

```julia-repl
julia> x = rand(SVector{2,Float64},100);

julia> box = Box([1,1],0.1);

julia> CellListMap.replicate_system!(x,box,(0:0,-1:1))
300-element Vector{SVector{2, Float64}}:
 [0.7119987163255118, 0.6788616154460262]
 [0.6188407316804118, 0.8497116428720384]
 [0.21328895963244354, 0.48932085643862977]
 ⋮
 [0.4114499470191678, 1.1034376619603892]
 [0.6094126258851252, 1.2328989485215263]
```

"""
function replicate_system!(
    x::AbstractVector{SVector{N,T}},
    unit_cell_matrix::AbstractMatrix,
    ranges::Tuple
) where {N,T}
    @assert length(ranges) == N "Tuple of ranges must have the same dimension as the vectors: $N"
    i0 = ntuple(i -> 0, N)
    imgs = Iterators.filter(!isequal(i0),
        Iterators.product(ranges...)
    )
    x0 = copy(x)
    for img in imgs
        x_new = translation_image(x0,unit_cell_matrix,img)
        append!(x,x_new)
    end
    return x
end

replicate_system!(x::AbstractVector,box::Box,ranges::Tuple) =
    replicate_system!(x,box.unit_cell.matrix,ranges)

function replicate_system!(x::AbstractMatrix{T},cell,ranges) where T
    N = size(x,1)
    x_re = [ SVector{N,T}(ntuple(i -> x[i,j], N)) for j in axes(x,2) ]
    replicate_system!(x_re,cell,ranges)
    x = Matrix(reinterpret(reshape,Float64,x_re))
    return x
end

"""

```
neighbor_cells_forward(box::Box{UnitCellType,N}) where UnitCellType
```

Internal function or structure - interface may change.

# Extended help

Function that returns the iterator of the cartesian indices of the cells that must be
evaluated (forward, i. e. to avoid repeated interactions)
if the cells have sides of length `box.cell_size`. `N` can be
`2` or `3`, for two- or three-dimensional systems.

"""
function neighbor_cells_forward(box::Box{UnitCellType,3}) where UnitCellType
    @unpack lcell = box
    nb = Iterators.flatten((
        CartesianIndices((1:lcell,-lcell:lcell,-lcell:lcell)),
        CartesianIndices((0:0,1:lcell,-lcell:lcell)),
        CartesianIndices((0:0,0:0,1:lcell))
    ))
    return nb
end

function neighbor_cells_forward(box::Box{UnitCellType,2}) where UnitCellType
    @unpack lcell = box
    nb = Iterators.flatten((
        CartesianIndices((1:lcell,-lcell:lcell)),
        CartesianIndices((0:0,1:lcell))
    ))
    return nb
end

"""

```
neighbor_cells(box::Box{UnitCellType,N}) where N
```

Internal function or structure - interface may change.

# Extended help

Function that returns the iterator of the cartesian indices of all neighboring
cells of a cell if the cells have sides of `box.cell_size`. `N` can be
`2` or `3`, for two- or three-dimensional systems.

"""
function neighbor_cells(box::Box{UnitCellType,3}) where UnitCellType
    @unpack lcell = box
    return CartesianIndices((-lcell:lcell,-lcell:lcell,-lcell:lcell))
end

function neighbor_cells(box::Box{UnitCellType,2}) where UnitCellType
    @unpack lcell = box
    return CartesianIndices((-lcell:lcell,-lcell:lcell))
end

"""

```
particle_cell(x::SVector{N,T}, box::Box) where {N,T}
```

Internal function or structure - interface may change.

# Extended help

Returns the coordinates of the *computing cell* to which a particle belongs, given its coordinates
and the `cell_size` vector. The computing box is always Orthorhombic, and the first
computing box with positive coordinates has indexes `Box.lcell + 1`.

"""
@inline particle_cell(x::SVector{N}, box::Box) where N =
    CartesianIndex(ntuple(i -> floor(Int,x[i]/box.cell_size[i]) + box.lcell + 1, N))

"""

```
cell_center(c::CartesianIndex{N},box::Box{UnitCellType,N,T}) where {UnitCellType,N,T}
```

Internal function or structure - interface may change.

# Extended help

Computes the geometric center of a computing cell, to be used in the projection
of points. Returns a `SVector{N,T}`

"""
@inline cell_center(c::CartesianIndex{N},box::Box{UnitCellType,N,T}) where {UnitCellType,N,T} =
    SVector{N,T}(ntuple(i -> box.cell_size[i]*(c[i] - box.lcell) - box.cell_size[i]/2, N))

"""

```
cell_cartesian_indices(nc::SVector{N,Int}, i1D) where {N}
```

Internal function or structure - interface may change.

# Extended help

Given the linear index of the cell in the cell list, returns the cartesian indices
of the cell (for arbitrary dimension N).

"""
@inline cell_cartesian_indices(nc::SVector{N,Int}, i1D) where N =
    CartesianIndices(ntuple(i -> nc[i],N))[i1D]

"""

```
cell_linear_index(nc::SVector{N,Int}, indices) where N
```

Internal function or structure - interface may change.

# Extended help

Returns the index of the cell, in the 1D representation, from its cartesian coordinates.

"""
@inline cell_linear_index(nc::SVector{N,Int}, indices) where N =
    LinearIndices(ntuple(i -> nc[i],N))[ntuple(i->indices[i],N)...]

"""

```
out_of_bounding_box(x::SVector{N},box::Box) where N
```

Internal function or structure - interface may change.

# Extended help

Function that evaluates if a particle is outside the computing bounding box,
defined by the maximum and minimum unit cell coordinates.

"""
function out_of_bounding_box(x::SVector{N},box::Box) where N
    @unpack cutoff, unit_cell_max = box
    for i in 1:N
        (x[i] < -cutoff) && return true
        (x[i] >= unit_cell_max[i]+cutoff) && return true
    end
    return false
end
out_of_bounding_box(p::ParticleWithIndex,box::Box) =
    out_of_bounding_box(p.coordinates,box)

"""

```
replicate_particle!(ip,p::SVector{N},box,cl) where N
```

Internal function or structure - interface may change.

# Extended help

Replicates the particle as many times as necessary to fill the computing box.

"""
function replicate_particle!(ip,p::SVector{N},box,cl) where N
    itr = Iterators.product(ntuple(i->box.ranges[i],N)...)
    for indexes in itr
        (count(isequal(0),indexes) == N) && continue
        x = translation_image(p,box.unit_cell.matrix,indexes)
        if ! out_of_bounding_box(x,box)
            cl = add_particle_to_celllist!(ip,x,box,cl;real_particle=false)
        end
    end
    return cl
end

"""

```
check_unit_cell(box::Box)
```

Internal function or structure - interface may change.

# Extended help

Checks if the unit cell satisfies the conditions for using the minimum-image
convention.

"""
check_unit_cell(box::Box) = check_unit_cell(box.unit_cell.matrix,box.cutoff)

function check_unit_cell(unit_cell_matrix::SMatrix{3},cutoff;printerr=true)
    a = @view(unit_cell_matrix[:,1])
    b = @view(unit_cell_matrix[:,2])
    c = @view(unit_cell_matrix[:,3])
    check = true

    if size(unit_cell_matrix) != (3,3)
        printerr && println("UNIT CELL CHECK FAILED: unit cell matrix must have dimenions (3,3).")
        check = false
    end

    if count(el -> el < zero(eltype(unit_cell_matrix)), unit_cell_matrix) != 0
         printerr && println("UNIT CELL CHECK FAILED: unit cell matrix components be strictly positive.")
         check = false
    end

    bc = cross(b,c)
    bc = bc / norm(bc)
    aproj = dot(a,bc)

    ab = cross(a,b)
    ab = ab / norm(ab)
    cproj = dot(c,ab)

    ca = cross(c,a)
    ca = ca / norm(ca)
    bproj = dot(b,ca)

    if (aproj < 2*cutoff) || (bproj < 2*cutoff) || (cproj < 2*cutoff)
        printerr && println("UNIT CELL CHECK FAILED: distance between cell planes too small relative to cutoff.")
        check = false
    end

    return check
end

function check_unit_cell(unit_cell_matrix::SMatrix{2},cutoff;printerr=true)
    a = @view(unit_cell_matrix[:,1])
    b = @view(unit_cell_matrix[:,2])
    check = true

    if size(unit_cell_matrix) != (2,2)
        printerr && println("UNIT CELL CHECK FAILED: unit cell matrix must have dimenions (2,2).")
        check = false
    end

    if count(el -> el < zero(eltype(unit_cell_matrix)), unit_cell_matrix) != 0
         printerr && println("UNIT CELL CHECK FAILED: unit cell matrix components must be strictly positive.")
         check = false
    end

    i = a / norm(a)
    bproj = sqrt(norm_sqr(b) - dot(b,i)^2)

    j = b / norm(b)
    aproj = sqrt(norm_sqr(a) - dot(a,j)^2)

    if (aproj < 2*cutoff) || (bproj < 2*cutoff)
         printerr && println("UNIT CELL CHECK FAILED: distance between cell planes too small relative to cutoff.")
         check = false
    end

    return check
end

#
# Compute the maximum and minimum coordinates of the vectors composing
# the particle sets
#
function _minmax(x::AbstractVector{<:AbstractVector})
    xmin = similar(x[begin])
    xmax = similar(x[begin])
    xmin .= typemax(eltype(xmin))
    xmax .= typemin(eltype(xmax))
    for v in x
        @. xmin = min(xmin,v)
        @. xmax = max(xmax,v)
    end
    return xmin, xmax
end

"""

```
limits(x)
```

Returns the lengths of a orthorhombic box that encompasses all the particles defined in `x`,
to be used to set a box without effective periodic boundary conditions.

"""
function limits(x::AbstractVector{<:AbstractVector})
    xmin, xmax = _minmax(x)
    return Limits(xmax .- xmin)
end

function limits(x::AbstractMatrix)
    N = size(x,1)
    @assert (N == 2 || N == 3) "The first dimension of the matrix must be the dimension (2 or 3)"
    x_re = reinterpret(reshape, SVector{N,eltype(x)}, x)
    return limits(x_re)
end

"""

```
limits(x,y)
```

Returns the lengths of a orthorhombic box that encompasses all the particles defined in `x`
and `y`, to used to set a box without effective periodic boundary conditions.

"""
function limits(x::T,y::T) where T <: AbstractVector{<:AbstractVector}
    xmin, xmax = _minmax(x)
    ymin, ymax = _minmax(y)
    xymin = min.(xmin,ymin)
    xymax = max.(xmax,ymax)
    return Limits(xymax .- xymin)
end

function limits(x::T,y::T) where T <: AbstractMatrix
    N = size(x,1)
    M = size(y,1)
    @assert N == M "The first dimension of the input matrices must be equal. "
    @assert (N == 2 || N == 3) "The first dimension of the matrix must be the dimension (2 or 3)"
    x_re = reinterpret(reshape, SVector{N,eltype(x)}, x)
    y_re = reinterpret(reshape, SVector{N,eltype(y)}, y)
    return limits(x_re,y_re)
end

#
# The following functions are not being used anymore in CellListMap, but may
# be useful if this is split in into a sepparate package for playing with
# periodic systems.
#
"""

```
ranges_of_replicas(cell_size, lcell, nc, unit_cell_matrix::SMatrix{3,3,T}) where T
```

Internal function or structure - interface may change.

# Extended help

Function that sets which is the range of periodic images necessary to fill
the computing box, in 3D.

"""
function ranges_of_replicas(cell_size, lcell, nc, unit_cell_matrix::SMatrix{3,3,T}) where T
    V = SVector{3,T}
    cmin = -lcell*cell_size
    cmax = (nc .- lcell) .* cell_size
    cell_vertices = SVector{8,V}(
        V(   cmin[1],   cmin[2],   cmin[3] ),
        V(   cmin[1],   cmin[2],   cmax[3] ),
        V(   cmin[1],   cmax[2],   cmin[3] ),
        V(   cmin[1],   cmax[2],   cmax[3] ),
        V(   cmax[1],   cmin[2],   cmin[3] ),
        V(   cmax[1],   cmin[2],   cmax[3] ),
        V(   cmax[1],   cmax[2],   cmin[3] ),
        V(   cmax[1],   cmax[2],   cmax[3] )
    )
    r_min, r_max = _ranges_of_replicas(
        SVector{3,Int}(10^6,10^6,10^6), #min
        SVector{3,Int}(-1,-1,-1),       #max
        unit_cell_matrix,
        cell_vertices
    )
    ranges = SVector{3,UnitRange{Int}}(
        r_min[1]:r_max[1],
        r_min[2]:r_max[2],
        r_min[3]:r_max[3]
    )
    return ranges
end

"""

```
ranges_of_replicas(cell_size, lcell, nc, unit_cell_matrix::SMatrix{2,2,T}) where T
```

Internal function or structure - interface may change.

# Extended help

Function that sets which is the range of periodic images necessary to fill
the computing box, in 2D.

"""
function ranges_of_replicas(cell_size, lcell, nc,unit_cell_matrix::SMatrix{2,2,T}) where T
    V = SVector{2,T}
    cmin = -lcell*cell_size
    cmax = (nc .- lcell) .* cell_size
    cell_vertices = SVector{4,V}(
        V(   cmin[1],   cmin[2] ),
        V(   cmin[1],   cmax[2] ),
        V(   cmax[1],   cmin[2] ),
        V(   cmax[1],   cmax[2] ),
    )
    r_min, r_max = _ranges_of_replicas(
        SVector{2,Int}(10^6,10^6), #min
        SVector{2,Int}(-1,-1),     #max
        unit_cell_matrix,
        cell_vertices
    )
    ranges = SVector{2,UnitRange{Int}}(
        r_min[1]:r_max[1],
        r_min[2]:r_max[2]
    )
    return ranges
end

function _ranges_of_replicas(r_min,r_max,unit_cell,cell_vertices)
    for vert in cell_vertices
        r = unit_cell \ vert
        ri = @. ceil(Int,abs(r))
        for (i,el) in pairs(r)
            if el < 0
                @set! ri[i] = -ri[i]
            end
        end
        r_min = min.(ri,r_min)
        r_max = max.(ri,r_max)
    end
    return r_min, r_max
end


text:
function hydrogenerator_dict(;index = 1,
    index_grid= 1,
    name = nothing,
    max_volume=100,
    min_volume=0,
    max_turn = 100,
    min_turn = 0,
    initial_volume = 0.0,
    production_factor = 1.0,
    spill_cost = 0.0,
    downstream_turn= [],
    downstream_spill= [])

    hydrogen = Dict("index" => index,
        "index_grid"=> index_grid,
        "max_volume"=>max_volume,
        "min_volume"=>min_volume,
        "max_turn" => max_turn,
        "min_turn" => min_turn,
        "initial_volume"=>initial_volume,
        "production_factor"=>production_factor,
        "spill_cost"=>spill_cost,
        "downstream_turn"=> downstream_turn,
        "downstream_spill"=> downstream_spill
    )
    if name != nothing
        hydrogen["name"] = name
    end
    return hydrogen
end

function hydro_dict(;hydrogenarray=[])
    hydro_dict = Dict("Hydrogenerators"=>hydrogenarray)
end


text:
type DSSsd <: StepData

	I # total points sampled
	s::Function
	u::Function
	getStep::Function
	stepString::AbstractString

	function DSSsd(s::Function, u::Function, stepString::AbstractString)
		new(0,s, u, DSScomputation, stepString)
	end
end

function DSS(thisUsd::USD)
	DSSsd(thisUsd.u,thisUsd.s, "DSS"*thisUsd.stepString)
end

function DSScomputation(x, k, gnum, sd::DSSsd, problem::Problem; outputLevel = 0)

	sumy=zeros(size(x))
	batchSize = sd.u(k,sd.I)+sd.s(k,sd.I)
	for i in 1:batchSize
		(f,sampleG) = (consume(problem.getNextSampleFunction))(x)
		sumy+=sampleG
	end

	gnum += sd.s(k,sd.I)+sd.u(k,sd.I)

	sd.I = sd.I + sd.u(k,sd.I)

	(sumy/batchSize, gnum)
end


text:
"""
Simplex Method
==============

References
----------

- CLRS's Introduction to Algorithms, Chapter 29.
- https://people.richland.edu/james/ictcm/2006/simplex.html
- http://www.phpsimplex.com/en/simplex_method_theory.htm

"""
module SimplexMethod

using LinearAlgebra

export simplex_method, find_basic_indices, find_optimal_x

#=

The standard form of linear programming:

    max    c^T x
    s.t.   A x = b
           x >= 0

Each intersection point (vertex of the feasible region) can be identified by $n - m$ variables being
0, i.e., $n - m$ non-basic variables. The optimal solution is at an extreme point of the simplex.
The algorithm starts from the origin (a basic feasible solution), iterates through the extreme
points until no improvement can be made.

Let `sB` be the bit vector indicating basic variables, then `sNB = !.sB` indicates the indices of
non-basic variables.

    b = A[:, sB] * x[sB] + A[:, sNB] * x[sNB]

    c' * x = c[sB]' * x[sB] + c[sNB]' * x[sNB]
           = c[sB]' * inv(A[:, sB]) * (b - A[:, sNB] * x[sNB]) + c[sNB]' * x[sNB]
           = c[sB]' * inv(A[:, sB]) * b + (c[sNB] - A[:, sNB]' * inv(A[:, sB])' * c[sB])' * x[sNB]

Basic V.S. Non-basic
--------------------

- Basic variables: nonzero (elementary, foundational, or important)
- Non-basic variables: always zero
- One basic variable for each row of the tableau
- The objective function is always basic in the bottom row
- Entering: non-basic -> basic
- Leaving: basic -> non-basic

- - - -

(Only the slack type implemented)

| Type of inequality |    Type of variable  |
|--------------------|----------------------|
|          ≥         | -surplus +artificial |
|          =         |      +artificial     |
|          ≤         |        +slack        |

=#

function find_basic_indices(S::AbstractMatrix{T})::Vector{Int} where T<:AbstractFloat
    (nr, nc) = size(S)
    indices = zeros(Int, nr)
    for c in 1:nc
        zero_cnt = 0
        one_idx = nothing
        for r = 1:nr
            if S[r, c] ≈ one(T)
                one_idx = r
            elseif S[r, c] ≈ zero(T)
                zero_cnt += 1
            end
        end
        if zero_cnt == nr - 1 && !isnothing(one_idx)
            indices[one_idx] = c
        end
    end
    indices
end

# `S` and `b` should be in the final configuration in the tableau.
function find_optimal_x(
    S::AbstractMatrix{T},
    b::AbstractVector{T},
)::Vector{T} where T<:AbstractFloat
    basic_indices = find_basic_indices(S)
    (m, m_plus_n) = size(S)
    x = zeros(T, m_plus_n)
    x[basic_indices] .= b
    x[1:(m_plus_n - m)]  # Removes slack variables
end

@enum LPResult begin
    # Infinite
    Unbounded
    # Infeasible
end

function simplex_method(
    A::AbstractMatrix{T},
    b::AbstractVector{T},
    c::AbstractVector{T},
)::Union{LPResult, Tuple{T, Vector{T}}} where T<:AbstractFloat
    (m, n) = size(A)
    tableau = zeros(T, (m + 1, n + m + 1))  # introduces m slack variables
    tableau[1:m, 1:n] = A
    tableau[1:m, (n + 1):(end - 1)] = I(m)
    tableau[end, 1:n] = -c'
    tableau[1:m, end] = b

    # Make some aliases
    S = @view tableau[1:(end - 1), 1:(end - 1)]  # the slack form
    coeff = @view tableau[end, 1:n]  # (negative, indeed) coefficients
    bias = @view tableau[1:m, end]

    while true
        (entering_val, entering_idx) = findmin(coeff)
        # There are no more negative coefficients, which means we cannot make the objective larger.
        if entering_val >= zero(T)
            z = tableau[end]
            x = find_optimal_x(S, bias)
            return (z, x)
        end

        # Pivoting. Similar to Gaussian elimination.
        pivot_row = nothing
        pivot_val = typemax(T)
        # Finds the tightest bound
        for r in 1:m
            if S[r, entering_idx] > zero(T)
                t = bias[r] / S[r, entering_idx]
                if t < pivot_val
                    pivot_row = r
                    pivot_val = t
                end
            end
        end

        # If the whole column is negative, the value can grow as large as possible without violating
        # the nonnegative constriaints.
        if isnothing(pivot_row)
            return Unbounded
        end

        # Normalizes row
        tableau[pivot_row, :] .*= 1 / tableau[pivot_row, entering_idx]
        for r in 1:(m + 1)  # including the last row
            if r != pivot_row
                tableau[r, :] .-= tableau[pivot_row, :] * tableau[r, entering_idx]
            end
        end
    end
end

end # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
using ManifoldMeasures
using Test

@testset "ManifoldMeasures.jl" begin
    # Write your tests here.
end


text:
struct LatentDiffEqModel{M,E,D}

    model_type::M
    encoder::E
    decoder::D

    @doc raw"""
        LatentDiffEqModel(model_type, encoder_layers, decoder_layers)
    Construct a LatentDiffEqModel.
    # Arguments
    `model_type`: GOKU() or LatentODE()\
    `encoder_layers`: contains `feature_extractor`, `pattern_extractor` and `latent_in` layers.\
    `decoder_layers`: contains `latent_out`, `diffeq` and `reconstructor` layers.\
    To see an example of how to build these layers, please refer to the `default_layers` function.
    """
    function LatentDiffEqModel(model_type, encoder_layers, decoder_layers)

        encoder = Encoder(model_type, encoder_layers)
        decoder = Decoder(model_type, decoder_layers)
        M, E, D = typeof(model_type), typeof(encoder), typeof(decoder)
        new{M, E, D}(model_type, encoder, decoder)
    end
end

function (model::LatentDiffEqModel)(x, t, variational=false)

    # Get encoded latent initial states and parameters
    μ, logσ² = model.encoder(x)

    # Sample from distributions
    l̃ = variational ? sample(μ, logσ², model) : μ

    # Get predicted output
    X̂ = model.decoder(l̃, t)

    return X̂, μ, logσ²
end

Flux.@functor LatentDiffEqModel

struct Encoder{M,FE,PE,LI}

    model_type::M

    feature_extractor::FE
    pattern_extractor::PE
    latent_in::LI

    @doc raw"""
        Encoder(model_type, encoder_layers)
    Construct an encoder of a LatentDiffEqModel.
    # Arguments
    `model_type`: GOKU() or LatentODE()\
    `encoder_layers`: contains `feature_extractor`, `pattern_extractor` and `latent_in` layers.  Please refer to `default_layers` for an example of how to construct this.
    """
    function Encoder(model_type, encoder_layers)
        M = typeof(model_type)
        FE, PE, LI = typeof.(encoder_layers)
        new{M,FE,PE,LI}(model_type, encoder_layers...)
    end
end

function (encoder::Encoder)(x)

    # Pass every time frame independently through the feature extractor
    fe_out = apply_feature_extractor(encoder, x)

    # Process sequentially with the pattern extractor
    pe_out = apply_pattern_extractor(encoder, fe_out)

    # Pass trough a last layer before sampling
    μ, logσ² = apply_latent_in(encoder, pe_out)

    return μ, logσ²
end

Flux.@functor Encoder

struct Decoder{M,LI,D,R}

    model_type::M

    latent_out::LI
    diffeq::D
    reconstructor::R

    @doc raw"""
        Decoder(model_type, decoder_layers)
    Construct a decoder of a LatentDiffEqModel.
    # Arguments
    `model_type`: GOKU() or LatentODE()\
    `decoder_layers`: contains `latent_out`, `diffeq` and `reconstructor` layers.  Please refer to `default_layers` for an example of how to construct this.
    """
    function Decoder(model_type, decoder_layers)
        M = typeof(model_type)
        LI, D, R = typeof.(decoder_layers)
        new{M,LI,D,R}(model_type, decoder_layers...)
    end
end

function (decoder::Decoder)(l̃, t)

    # Pass sampled latent states through a latent_out layer
    l̂ = apply_latent_out(decoder, l̃)

    # Integrate differential equations
    ẑ = diffeq_layer(decoder, l̂, t)

    # Apply reconstructor independently to each time frame
    x̂ = apply_reconstructor(decoder, ẑ)

    return x̂, ẑ, l̂
end

Flux.@functor Decoder


text:
using DiffEqBase,
    DiffEqOperators, ForwardDiff, LinearAlgebra, SparseDiffTools, Test, LinearSolve
const A = rand(300, 300)
f(du, u) = mul!(du, A, u)
f(u) = A * u
x = rand(300)
v = rand(300)
du = similar(x)

cache1 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x)))))
cache2 = ForwardDiff.Dual{typeof(ForwardDiff.Tag(SparseDiffTools.DeivVecTag(),eltype(x))),eltype(x),1}.(x, ForwardDiff.Partials.(tuple.(reshape(v, size(x)))))
@test num_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f, similar(x), x) * v rtol =
    1e-6
@test num_jacvec!(du, f, x, v, similar(v), similar(v)) ≈
      ForwardDiff.jacobian(f, similar(x), x) * v rtol = 1e-6
@test num_jacvec(f, x, v) ≈ ForwardDiff.jacobian(f, similar(x), x) * v rtol =
    1e-6

@test auto_jacvec!(du, f, x, v) ≈ ForwardDiff.jacobian(f, similar(x), x) * v
@test auto_jacvec!(du, f, x, v, cache1, cache2) ≈
      ForwardDiff.jacobian(f, similar(x), x) * v
@test auto_jacvec(f, x, v) ≈ ForwardDiff.jacobian(f, similar(x), x) * v

f(du, u, p, t) = mul!(du, A, u)
f(u, p, t) = A * u
L = JacVecOperator(f, x)
@test L * x ≈ auto_jacvec(f, x, x)
@test L * v ≈ auto_jacvec(f, x, v)
@test mul!(du, L, v) ≈ auto_jacvec(f, x, v)
DiffEqBase.update_coefficients!(L, v, nothing, nothing)
@test mul!(du, L, v) ≈ auto_jacvec(f, v, v)
# reshape
xx = [x x]
vv = [v v]
duu = [du du]
L = JacVecOperator(f, xx)
@test L * vec(xx) ≈ mul!(vec(duu), L, vec(xx))
@test vec(duu) ≈ vec(A * xx) rtol = 1e-6

L = JacVecOperator(f, x, autodiff = false)
DiffEqBase.update_coefficients!(L, x, nothing, nothing)
@test L * x ≈ num_jacvec(f, x, x)
@test L * v ≈ num_jacvec(f, x, v)
@test mul!(du, L, v) ≈ num_jacvec(f, x, v) rtol = 1e-6
DiffEqBase.update_coefficients!(L, v, nothing, nothing)
@test mul!(du, L, v) ≈ num_jacvec(f, v, v) rtol = 1e-6
L = JacVecOperator(f, xx, autodiff = false)
DiffEqBase.update_coefficients!(L, vv, nothing, nothing)
@test L * vec(xx) ≈ mul!(vec(duu), L, vec(xx)) rtol = 1e-6
@test vec(duu) ≈ vec(A * xx)

L2 = JacVecOperator{Float64}(f)
DiffEqBase.update_coefficients!(L2, x, nothing, nothing)
@test L2 * x ≈ auto_jacvec(f, x, x)
@test L2 * v ≈ auto_jacvec(f, v, v)
@test mul!(du, L2, x) ≈ auto_jacvec(f, x, x)
DiffEqBase.update_coefficients!(L2, xx, nothing, nothing)
@test L2 * vec(xx) ≈ auto_jacvec(f, xx, vec(xx))
@test L2 * vec(vv) ≈ auto_jacvec(f, vv, vec(vv))
@test mul!(vec(duu), L2, vec(xx)) ≈ auto_jacvec(f, xx, vec(xx))

L2 = JacVecOperator{Float64}(f, autodiff = false)
DiffEqBase.update_coefficients!(L2, x, nothing, nothing)
@test L2 * x ≈ num_jacvec(f, x, x)
@test L2 * v ≈ num_jacvec(f, v, v) rtol = 1e-6
DiffEqBase.update_coefficients!(L2, xx, nothing, nothing)
@test L2 * vec(xx) ≈ num_jacvec(f, xx, vec(xx))
@test L2 * vec(vv) ≈ num_jacvec(f, vv, vec(vv)) rtol = 1e-6

using OrdinaryDiffEq
function lorenz(du, u, p, t)
    du[1] = 10.0(u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)
ff1 = ODEFunction(lorenz, jac_prototype = JacVecOperator{Float64}(lorenz, u0))
ff2 = ODEFunction(
    lorenz,
    jac_prototype = JacVecOperator{Float64}(lorenz, u0, autodiff = false),
)

for ff in [ff1, ff2]
    prob = ODEProblem(ff, u0, tspan)
    @test solve(prob, TRBDF2()).retcode == :Success
    @test solve(prob, TRBDF2(linsolve = KrylovJL_GMRES())).retcode == :Success
    @test solve(prob, Exprb32()).retcode == :Success
    @test sol = solve(prob, Rosenbrock23()).retcode == :Success
    @test sol = solve(prob, Rosenbrock23(linsolve = KrylovJL_GMRES())).retcode == :Success
end


text:
# ---
# title: 715. Range Module
# id: problem715
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Segment Tree, Ordered Map
# link: <https://leetcode.com/problems/range-module/description/>
# hidden: true
# ---
#
# A Range Module is a module that tracks ranges of numbers. Your task is to
# design and implement the following interfaces in an efficient manner.
#
# * `addRange(int left, int right)` Adds the half-open interval `[left, right)`, tracking every real number in that interval. Adding an interval that partially overlaps with currently tracked numbers should add any numbers in the interval `[left, right)` that are not already tracked.
#
# * `queryRange(int left, int right)` Returns true if and only if every real number in the interval `[left, right)` is currently being tracked.
#
# * `removeRange(int left, int right)` Stops tracking every real number currently being tracked in the interval `[left, right)`.
#
# **Example 1:**
#
#
#
#     **addRange(10, 20)** : null
#     **removeRange(14, 16)** : null
#     **queryRange(10, 14)** : true (Every number in [10, 14) is being tracked)
#     **queryRange(13, 15)** : false (Numbers like 14, 14.03, 14.17 in [13, 15) are not being tracked)
#     **queryRange(16, 17)** : true (The number 16 in [16, 17) is still being tracked, despite the remove operation)
#
#
# **Note:**
#
# * A half open interval `[left, right)` denotes all real numbers `left <= x < right`.
# * `0 < left < right < 10^9` in all calls to `addRange, queryRange, removeRange`.
# * The total number of calls to `addRange` in a single test case is at most `1000`.
# * The total number of calls to `queryRange` in a single test case is at most `5000`.
# * The total number of calls to `removeRange` in a single test case is at most `1000`.
#
#
## @lc code=start
using LeetCode

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


text:
module DelaunayTessellationFieldEstimator

# Implements the DTFE algoritm as found in
# Schaap, W.E., 2007. The delaunay tessellation field estimator. Ph. D. Thesis.
# https://www.rug.nl/research/portal/files/2816076/c2.pdf
# Author Harvey Devereux (h.devereux@warwick.ac.uk)
# License: MIT
# Any collaboration welcomed

# university wesite https://warwick.ac.uk/fac/sci/mathsys/people/students/2017intake/devereux/
# personal website https://harveydevereux.github.io/

# Works now but is slower than python (possibly due to pyimport overhead?)

# TODO See if GridInterpolations.jl SimplexGrid is faster than Py-importing
# TODO 3D -> N-D
# TODO Gif visualisation

using VoronoiDelaunay
import VoronoiDelaunay: delaunayedges
using GeometricalPredicates
using NearestNeighbors
import NearestNeighbors: knn, KDTree

# for the simplex based interpolation
using PyCall
@pyimport scipy.interpolate as si

export
DTFE,
Coordinates, GeomToDataAreaScalingFactor, BuildDelaunay!, FindTriangles,
ContiguousVoronoiCellArea, DelaunayTessellation2D, DelaunayTessellation,
Point, Point2D, getplotxyd, delaunayedges, UniquePoints, GeomToDataAreaScalingFactor,
DelaunayVertexDensity, DelaunayVertexGeometryDensity, Point2DToArray, DTFEMap,
Grid, meshgrid

const min_lim = VoronoiDelaunay.min_coord
const max_lim = VoronoiDelaunay.max_coord

"""Holds the original data and the rescaled geometric data to
    comform the VoronoiDelaunay\'s range requirements

    NB: GeometricalPredicates requires points ∈ [1,2] but
        VoronoiDelaunay requires points in [1+eps(float),2-2eps(float)]
"""
mutable struct Coordinates{T <: Array{<:Real,2}}
    DataCoordinates::T
    GeometryCoordinates::T
    DataLims::T
    GeomLims::T
end

"""Constructor to deal which automatically calculates the geometry scale
   data and holds the limits for future scaling.
"""
function Coordinates{T <: Array{<:Real,2}}(Points::T)
    Glim = [min_lim max_lim; min_lim max_lim]
    Dlim = [minimum(Points[:,1]) maximum(Points[:,1]); minimum(Points[:,2]) maximum(Points[:,2])]
    GPoints = Scaling(Points, Glim)
    return Coordinates(Points,GPoints,Dlim,Glim)
end

"""Min max scaling to the VoronoiDelaunay range requirement
"""
function Scaling{T <: Array{<:Real,2}}(coords::T,lim::T)
    Rescale = zeros(size(coords))
    for d in 1:2
        min = minimum(coords[:,d])
        max = maximum(coords[:,d])
        a = minimum(lim[d,:])
        b = maximum(lim[d,:])
        Rescale[:,d] = [ ((coords[i,d] - min)/(max-min))*(b-a)+a for i in 1:size(coords,1)]
    end
    return Rescale
end

"""Returns the scaling factor required to calculate the area in the original
   point space.

   Here M,m and the new min and max (assumed the same for all dimensions)
   and min,max are the dimensions maximum and minimum.

   If the scaling is x* = [(x-min)/(max-min)](M-m)+m by writing
   out (a₁-b₁)(a₂-b₂) in the form of x* you\'ll find
   (a₁-b₁)(a₂-b₂)*SF⁻¹
   where SF  = ∏ᵢ(M-m)/(maxᵢ-minᵢ)
"""
function GeomToDataAreaScalingFactor(C::Coordinates)
    prod = 1
    D = C.DataLims
    G = C.GeomLims
    for d in 1:size(D,1)
        prod *= (maximum(D[d,:])-minimum(D[d,:]))/(maximum(G[d,:])-minimum(G[d,:]))
    end
    return prod
end

"""Wrapper to build a Delaunay tessellation from an
   array of coordinates

   returns the points, the tess object is modified so
   no return necessary
"""
# TODO same function but with the coordinates object?
function BuildDelaunay!{A <: Array{<: Real,2}}(Coords::A)
    # doesnt't seem to exist
    # sizehint(tess,size(Coords,1))
    # This because if you try to push
    # to an old tess the function hangs
    # in an infinite loop
    # there is an issue for this in
    # VoronoiDelaunay.jl
    Tess = DelaunayTessellation()
    Points = Point2D[Point(Coords[i,1],Coords[i,2]) for i in 1:size(Coords,1)]
    push!(Tess, Points)
    return Tess,Points
end

"""Usefull for quikcly getting th triangles vertices
all at once, returns a T <: Point2D
"""
function UnpackVertices{T <: VoronoiDelaunay.DelaunayTriangle}(Triangle::T)
    return geta(Triangle), getb(Triangle), getc(Triangle)
end

"""Usefull for quikcly getting th triangles vertices
all at once, returns a T <: Array
"""
function UnpackVerticesAsArray{T <: VoronoiDelaunay.DelaunayTriangle}(Triangle::T)
    a = [getx(geta(Triangle)),gety(geta(Triangle))]
    b = [getx(getb(Triangle)),gety(getb(Triangle))]
    c = [getx(getc(Triangle)),gety(getc(Triangle))]
    return a,b,c
end

"""Finds all the delaunay triangles with the common vertex of \'point\'
   returns an array of these triangles.
"""
function FindTriangles{T <: VoronoiDelaunay.DelaunayTessellation2D,
                       P <: GeometricalPredicates.Point2D}(tess::T,point::P)
    Triangles = Array{VoronoiDelaunay.DelaunayTriangle,1}()
    for triangle in tess
        a, b, c = UnpackVertices(triangle)
        if point == a || point == b || point == c
            push!(Triangles,triangle)
        end
    end
    return Triangles
end

"""The contiguous voronoi cell at a point i (Schaap, 2007 PhD) is the union of
   all delaunay triangles with the coomon vertex i

   returns the Geometry scale area of the union of the triangles
   in \'Triangles\', use FindTriangles to get these for a point.
"""
function ContiguousVoronoiCellArea{T <: VoronoiDelaunay.DelaunayTriangle}(Triangles::Array{T,1})
    Area=0
    for triangle in Triangles
        Area+=GeometricalPredicates.area(triangle)
    end
    return Area
end

"""Finds the set of unique points within a contigous
   Voronoi cell

   returns the array of unique points
"""
function UniquePoints{T <: VoronoiDelaunay.DelaunayTriangle}(Triangles::Array{T,1})
    Uniques = []
    for triangle in Triangles
        for v in UnpackVertices(triangle)
            if !(v in Uniques)
                push!(Uniques,v)
            end
        end
    end
    return Uniques
end

"""Finds the density in the geometry space (no scaling)
   for a single Delaunay vertex point
"""
function DelaunayVertexGeometryDensity{T <: VoronoiDelaunay.DelaunayTessellation2D,
                                       P <: GeometricalPredicates.Point2D}(tess::T,point::P,mass::Float64=1.0)
    Triangles = FindTriangles(tess,point)
    A = ContiguousVoronoiCellArea(Triangles)
    ρ = mass*length(UniquePoints(Triangles))/A
    return ρ
end

"""Wraps around DelaunayVertexGeometryDensity to calculate
   the density in data space for every delaunay vertex point

   Masses can be set to some field value other than mass. The
   default is a uniform unit mass.

   Returns an array of densities scaled to the data space
"""
function DelaunayVertexDensity{T <: VoronoiDelaunay.DelaunayTessellation2D}(tess::T,C::Coordinates,Masses=ones(length(C.GeometryCoordinates)))
    G = C.GeometryCoordinates
    ρ = zeros(size(G,1))
    for i in 1:size(G,1)
        P = Point2D(G[i,1],G[i,2])
        ρ[i] = DelaunayVertexGeometryDensity(tess,P,Masses[i])
    end
    return ρ./GeomToDataAreaScalingFactor(C)
end

meshgrid(x,y) = (repmat(x',length(y),1),repmat(y,1,length(x)))

"""Computes the interpolated DTFE grid from the given points
   with resolution given by Size
"""
function DTFE(Points,Size=500)
    C = Coordinates(Points)
    tess,P = BuildDelaunay!(C.GeometryCoordinates)
    x,y = meshgrid(linspace(C.DataLims[1,:]...,Size),linspace(C.DataLims[2,:]...,Size))

    D = DelaunayTessellationFieldEstimator.DelaunayVertexDensity(tess,C)
    grid = si.griddata(C.DataCoordinates,D,(x,y),method="linear",fill_value=0)
    return grid
end

end # module DelaunayTessellationFieldEstimator


text:
# specialrules.jl - differentiation rules for special functions

# TODO: add a couple of usefule special functions
# See test/special.jl for dummy functions


text:
# run in REPL
using Makie

function deathstar()
    n = 60
    θ = [0; (0.5: n - 0.5) / n; 1]
    φ = [(0: 2n - 2) * 2 / (2n - 1); 2]
    # if x is +0.9 radius units, replace it with the coordinates of sphere surface
    # at (1.2,0,0) center, radius 0.5 units
    x = [(x1 = cospi(φ)*sinpi(θ)) > 0.9 ? 1.2 - x1 * 0.5 : x1 for θ in θ, φ in φ]
    y = [sinpi(φ)*sinpi(θ) for θ in θ, φ in φ]
    z = [cospi(θ) for θ in θ, φ in φ]
    scene = Scene(backgroundcolor=:black)
    surface!(scene, x, y, z, color = rand(RGBAf0, 124, 124), show_axis=false)
end

deathstar()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Distributed
using MixedModelsSim, MixedModels,MixedModelsPermutations
using ProgressMeter
using SharedArrays
function sim_model(f)
    subj_btwn = Dict("age" => ["O", "Y"])

    # there are no between-item factors in this design so you can omit it or set it to nothing
    item_btwn = Dict("stimType" => ["I", "II"])

    # put within-subject/item factors in a Dict
    both_win = Dict("condition" => ["A", "B"])

    # simulate data
    dat = simdat_crossed(
        30,
        30,
        subj_btwn = subj_btwn,
        item_btwn = item_btwn,
        both_win = both_win,
    )

    simMod = MixedModels.fit(MixedModel, f, dat)

    return simMod

end
function run_permutationtest_distributed(n_workers, nRep, simMod,args...)

    if nworkers() < n_workers
        # open as many as necessary
        println("Starting Workers, this might take some time")
        addprocs(
            n_workers - nworkers() + 1,
            exeflags = "--project",
            enable_threaded_blas = true,
        )
    end

    # activate environment
    eval(macroexpand(Distributed, quote
        @everywhere using Pkg
    end))

    @everywhere Pkg.activate(".")
    # load packages on distributed
    eval(
        macroexpand(
            Distributed,
            quote
                @everywhere using MixedModelsSim,
                    Random, MixedModels, MixedModelsPermutations
            end,
        ),
    )
    β_permResult = SharedArray{Float64}(nRep, length(β))
    z_permResult = SharedArray{Float64}(nRep, length(β))

    @everywhere include("test/sim_utilities.jl")

    println("starting @distributed")
    # parallel loop
    @showprogress @distributed for k = 1:nRep
        #println("Thread "*string(Threads.threadid()) * "\t Running "*string(k))
        #res = [1,1.]#
        res = run_permutationtest(MersenneTwister(k), deepcopy(simMod),args...)
        β_permResult[k, :] .= res[1]
        z_permResult[k, :] .= res[2]

    end
    return β_permResult, z_permResult

end
function run_permutationtest(rng, simMod, nPerm, β, σ, θ)
    simMod = simulate!(rng, simMod, β = β, σ = σ, θ = θ)
    refit!(simMod)
    H0 = coef(simMod)
    H0[2] = 0.0
    #H0[1] = 0.0
    perm = permutation(rng, nPerm, simMod, use_threads = false; β = H0)
    p_β = values(permutationtest(perm, simMod; statistic = :β))
    p_z = values(permutationtest(perm, simMod; statistic = :z))
    return (p_β, p_z)
end

function fitsignal(formula, data, signal, contrasts)
    # fit the MixedModel

    fits = Array{Any}(undef, size(signal)[2])
    model = Array{Any}(undef, 1)
    cdata = copy(data)

    for i = 1:(size(signal)[2])
        println(i)
        if i == 1
            cdata[:, formula.lhs.sym] = (signal[:, i])
            model[1] = MixedModels.fit(MixedModel, formula, cdata, contrasts = contrasts)
        else
            model[1] = refit!(model[1], signal[:, i])
        end
        fits[i] = deepcopy(model[1])
    end
    return fits
end

function circulant(x)
    # returns a symmetric matrix where X was circ-shifted.
    lx = length(x)
    ids = [1:1:(lx-1);]
    a = Array{Float64,2}(undef, lx,lx)
    for i = 1:length(x)
        if i==1
            a[i,:] = x
        else
            a[i,:] = vcat(x[i],a[i-1,ids])
        end
    end
    return Symmetric(a)
end

function circulant(x)
    # returns a symmetric matrix where X was circ-shifted.
    lx = length(x)
    ids = [1:1:(lx-1);]
    a = Array{Float64,2}(undef, lx, lx)
    for i = 1:length(x)
        if i == 1
            a[i, :] = x
        else
            a[i, :] = vcat(x[i], a[i-1, ids])
        end
    end
    return Symmetric(a)
#endurn reduce(vcat,a)
end

function exponentialCorrelation(x; nu = 1, length_ratio = 1)
    # generate exponential function
    R = length(x) * length_ratio
    return exp.(-3 * (x / R) .^ nu)
end

function expandgrid(df1, df2)
    # get all combinations of df1&df2

    a = Array{Any}(undef, nrow(df1))
    for i = 1:nrow(df1)
        a[i] = hcat(repeat(df1[[i], :], nrow(df2)), df2)
    end
    return reduce(vcat, a)
end