text:
export switch_map

import Base: show

"""
    switch_map(::Type{R}, mappingFn::F = identity) where { R, F <: Function }

Creates a `switch_map` operator, which returns an Observable that emits items based on applying a function `mappingFn` that you supply to each item
emitted by the source Observable, where that function returns an (so-called "inner") Observable. Each time it observes one of these inner Observables,
the output Observable begins emitting the items emitted by that inner Observable. When a new inner Observable is emitted, `switch_map` stops emitting
items from the earlier-emitted inner Observable and begins emitting items from the new one. It continues to behave like this for
subsequent inner Observables.

# Arguments
- `::Type{R}`: the type of data of output Observables after projection with `mappingFn`
- `mappingFn::F`: projection function with `(data) -> Observable{R}` signature

# Producing

Stream of type `<: Subscribable{R}`

# Examples
```jldoctest
using Rocket

source = from([ of(1), of(2), of(3) ])
subscribe!(source |> switch_map(Int), logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Completed
```

```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> switch_map(Float64, (d) -> of(convert(Float64, d ^ 2))), logger())
;

# output

[LogActor] Data: 1.0
[LogActor] Data: 4.0
[LogActor] Data: 9.0
[LogActor] Completed
```

See also: [`AbstractOperator`](@ref), [`RightTypedOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
switch_map(::Type{R}, mappingFn::F = identity) where { R, F <: Function } = SwitchMapOperator{R, F}(mappingFn)

struct SwitchMapOperator{R, F} <: RightTypedOperator{R}
    mappingFn :: F
end

function on_call!(::Type{L}, ::Type{R}, operator::SwitchMapOperator{R, F}, source) where { L, R, F }
    return proxy(R, source, SwitchMapProxy{L, R, F}(operator.mappingFn))
end

struct SwitchMapProxy{L, R, F} <: ActorSourceProxy
    mappingFn :: F
end

actor_proxy!(::Type{R}, proxy::SwitchMapProxy{L, R, F}, actor::A) where { L, R, F, A } = SwitchMapActor{L, R, F, A}(proxy.mappingFn, actor)

# m - main
# i - inner
mutable struct SwitchMapActorProps
    msubscription :: Teardown
    isubscription :: Teardown
    ismcompleted  :: Bool
    isicompleted  :: Bool
    isdisposed    :: Bool

    SwitchMapActorProps() = new(voidTeardown, voidTeardown, false, true, false)
end

struct SwitchMapActor{L, R, F, A} <: Actor{L}
    mappingFn :: F
    actor     :: A
    props     :: SwitchMapActorProps

    SwitchMapActor{L, R, F, A}(mappingFn::F, actor::A) where { L, R, F, A } = new(mappingFn, actor, SwitchMapActorProps())
end

ismcompleted(actor::SwitchMapActor) = actor.props.ismcompleted
isicompleted(actor::SwitchMapActor) = actor.props.isicompleted

setmcompleted!(actor::SwitchMapActor, value::Bool) = actor.props.ismcompleted = value
seticompleted!(actor::SwitchMapActor, value::Bool) = actor.props.isicompleted = value

isdisposed(actor::SwitchMapActor)   = actor.props.isdisposed

struct SwitchMapInnerActor{R, S} <: Actor{R}
    main :: S
end

on_next!(actor::SwitchMapInnerActor{R}, data::R) where R = next!(actor.main.actor, data)
on_error!(actor::SwitchMapInnerActor,   err)             = error!(actor.main, err)
on_complete!(actor::SwitchMapInnerActor)                 = begin
    seticompleted!(actor.main, true)
    if ismcompleted(actor.main)
        complete!(actor.main)
    end
end

function on_next!(actor::S, data::L) where { L, R, S <: SwitchMapActor{L, R} }
    if !isdisposed(actor)
        unsubscribe!(actor.props.isubscription)
        seticompleted!(actor, false)
        actor.props.isubscription = subscribe!(actor.mappingFn(data), SwitchMapInnerActor{R, S}(actor))
    end
end

function on_error!(actor::SwitchMapActor, err)
    if !isdisposed(actor)
        dispose!(actor)
        error!(actor.actor, err)
    end
end

function on_complete!(actor::SwitchMapActor)
    setmcompleted!(actor, true)
    if !isdisposed(actor) && isicompleted(actor)
        dispose!(actor)
        complete!(actor.actor)
    end
end

function dispose!(actor::SwitchMapActor)
    actor.props.isdisposed = true
    unsubscribe!(actor.props.msubscription)
    unsubscribe!(actor.props.isubscription)
    return nothing
end

struct SwitchMapSource{L, S} <: Subscribable{L}
    source :: S
end

source_proxy!(::Type, proxy::SwitchMapProxy{L}, source::S) where { L, S } = SwitchMapSource{L, S}(source)

function on_subscribe!(source::SwitchMapSource, actor::SwitchMapActor)
    actor.props.msubscription = subscribe!(source.source, actor)
    return SwitchMapSubscription(actor)
end

struct SwitchMapSubscription{A} <: Teardown
    actor :: A
end

as_teardown(::Type{ <: SwitchMapSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::SwitchMapSubscription)
    dispose!(subscription.actor)
    return nothing
end

Base.show(io::IO, ::SwitchMapOperator{R})   where {    R } = print(io, "SwitchMapOperator($R)")
Base.show(io::IO, ::SwitchMapProxy{L, R})   where { L, R } = print(io, "SwitchMapProxy($L, $R)")
Base.show(io::IO, ::SwitchMapActor{L, R})   where { L, R } = print(io, "SwitchMapActor($L -> $R)")
Base.show(io::IO, ::SwitchMapInnerActor{R}) where {    R } = print(io, "SwitchMapInnerActor($R)")
Base.show(io::IO, ::SwitchMapSource{S})     where S        = print(io, "SwitchMapSource($S)")
Base.show(io::IO, ::SwitchMapSubscription)                 = print(io, "SwitchMapSubscription()")


text:
using QuadGK

@testset "Integration and Coupling generation" begin
	# ----------------------------------------
	# Test the numerical integration capabilities
	antiDerX = ChiralWaveguide.__integrate(x -> x, (0.0, 1.0))

	@test antiDerX(1.0) ≈ 1/2	# local accuracy
	@test sum(x -> abs(antiDerX(x) - x^2 / 2), [0.0:0.01:1.0;]) ≤ 1e-12 # global accuracy
end

@testset "Modes with offset" begin
	μ, σ = 10.0, 0.1

	# ----------------------------------------
	# Test normalisation
	compressions = [Algebraic(σ = σ, μ = μ), Exponential(σ = σ, μ = μ), Trigonometric(σ = σ, μ = μ)]
	for compression ∈ compressions
		testMode = Mode(
			t -> abs(t - μ) < σ/2 ? √(1/σ) : 0.0, compression = compression,
			abstol = 1e-8, reltol = 1e-10
		)
		@test abs(testMode.norm(2μ) - 1.0) < 1e-7
	end
end

@testset "Numerical Mode Couplings" begin
	modes = [
		FourierMode(1/√(5), [1/√(5), 1/√(5)], [1/√(5), 1/√(5)], σ = 3.0),
		HardBoxMode(σ = 3.0),    # give a width to help the integrator
		SoftBoxMode(σ = 3.0),    # A little bit of fine tuning is involved here, but
		SoftBoxExpMode(σ = 2.5), # I simple want to test the normalisation conditions.
		GaussMode(),
		ExpMode(γ =  4.0),
		ExpMode(γ = -4.0)
	]
	# ----------------------------------------
	# Test Numerical Derivations
	ts = range(-3.0, 3.0, length = 100)
	for mode ∈ modes
		numericalMode = Mode(t -> mode.modeFunction(t), reltol = 1e-8, abstol = 1e-12)
		@test sum(abs, numericalMode.gₒ.(ts) .- mode.gₒ.(ts)) / length(ts) < 1e-7
	end

	# ----------------------------------------
	# Integration test
	for mode ∈ modes
		@test quadgk(
			t -> mode.modeFunction(t) * mode.gₒ(t),
			-10.0, 10.0
		)[1] ≈ -2.0
	end
end

@testset "FourierMode and HardBoxMode equivalence" begin
	hb = HardBoxMode()
	f  = FourierMode(1., [], [])

	ts = range(0.0, 1.0, length = 101)

	@test sum(abs, hb.modeFunction.(ts) .- f.modeFunction.(ts)) / length(ts) < 1e-9
	@test sum(abs, hb.gᵢ.(ts) .- f.gᵢ.(ts)) / length(ts) < 1e-9
	@test sum(abs, hb.gₒ.(ts) .- f.gₒ.(ts)) / length(ts) < 1e-9
	@test sum(abs, hb.norm.(ts) .- f.norm.(ts)) / length(ts) < 1e-9
end

using QuantumOptics
@testset "Cavity Hopping" begin
	# Idea of this test is, that a cavity for a given mode
	# should have 100% into a cavity for the same mode.
	modes = [
		FourierMode(rand(), rand(2), rand(2), t₀ = -5.0),
		HardBoxMode(t₀ = -5.0),
		SoftBoxMode(),
		SoftBoxExpMode(),
		GaussMode(),
		ExpMode(γ =  3.5), # Needed to fit in the time window
		ExpMode(γ = -3.5)
	]

	for mode ∈ modes
		for n ∈ [1, 2, 3] # Fock state scattering
		  basis = FockBasis(n)
		  a = destroy(basis)

		  gᵢ, gₒ = mode.gᵢ, mode.gₒ

		  H(t) = 0.5im * gᵢ(t) * gₒ(t) * (a' ⊗ a - a ⊗ a')
		  L(t) = gᵢ(t) * a ⊗ one(basis) + gₒ(t) * one(basis) ⊗ a

		  ψ₀ = fockstate(basis, n) ⊗ fockstate(basis, 0)

		  local ts, ρs
		  ts, ρs = timeevolution.master_dynamic(
				[-5:0.01:5.0;], ψ₀, (t,ψ) -> (H(t), [L(t)], [L(t)']), reltol = 1e-7, abstol = 1e-9
			)

		  @test abs(expect(2, a'a, ρs[end]) - n) < 1e-7
		end
	end
end

ext: jl


text:
function read_Species_From_Text_File(file_name)
    return reaction_matrix
end

function calculate_concentration(c_0, λ, Δt)
    Λ = Diagonal(λ)
    R = assemble_transformation_matrix_inverse(λ)
    S = assemble_transformation_matrix(λ)
    A = S * Λ * R * Δt
    return exp(A) * c_0
end

function transformation_matrix_inverse_ij(i, j, λ)
    s_ij::eltype(λ) = 1.0
    # Loop over species ancestor index l for species i with (i-j) generations
    for l = j:(i-1)
        s_ij *= λ[l] / (λ[l] - λ[i])
    end
    return s_ij
end

function transformation_matrix_ij(i, j, λ)
    s_ij = λ[j] / (λ[i] - λ[j])
    # Loop over species ancestor index l for species i with (i-j) generations
    for l = (j+1):(i-1)
        s_ij *= λ[l] / (λ[l] - λ[j])
    end
    return s_ij
end

function assemble_transformation_matrix(λ)
    n = length(λ)
    S = zeros(eltype(λ), length(λ), length(λ))
    S = UnitLowerTriangular(S)
    for i = 2:n
        for j = 1:(i-1)
            S[i, j] = transformation_matrix_ij(i, j, λ)
        end
    end
    return S
end

function assemble_transformation_matrix_inverse(λ)
    n = length(λ)
    R = zeros(eltype(λ), length(λ), length(λ))
    R = UnitLowerTriangular(R)
    for i = 2:n
        for j = 1:i
            R[i, j] = transformation_matrix_inverse_ij(i, j, λ)
        end
    end
    return R
end

function retrieve_generation_number(i, j)
    return (i - j)
end

function transformation_matrix_5x5(λ = [1.0; 2.0; 3.0; 4.0; 5.0])

    s = zeros(5, 5)
    s[1, 1] = 1

    s[2, 1] = λ[1] / (λ[2] - λ[1])
    s[3, 1] = λ[1] / (λ[3] - λ[1]) * λ[2] / (λ[2] - λ[1])
    s[4, 1] = λ[1] / (λ[4] - λ[1]) * λ[2] / (λ[2] - λ[1]) * λ[3] / (λ[3] - λ[1])
    s[5, 1] =
        λ[1] / (λ[5] - λ[1]) * λ[2] / (λ[2] - λ[1]) *
        λ[3] / (λ[3] - λ[1]) * λ[4] / (λ[4] - λ[1])

    s[2, 2] = 1
    s[3, 2] = λ[2] / (λ[3] - λ[2])
    s[4, 2] = λ[2] / (λ[4] - λ[2]) * λ[3] / (λ[3] - λ[2])
    s[5, 2] = λ[2] / (λ[5] - λ[2]) * λ[3] / (λ[3] - λ[2]) * λ[4] / (λ[4] - λ[2])

    s[3, 3] = 1
    s[4, 3] = λ[3] / (λ[4] - λ[3])
    s[5, 3] = λ[3] / (λ[5] - λ[3]) * λ[4] / (λ[4] - λ[3])

    s[4, 4] = 1
    s[5, 4] = λ[4] / (λ[5] - λ[4])

    s[5, 5] = 1

    return s
end

function transformation_matrix_inverse_5x5(λ = [1.0; 2.0; 3.0; 4.0; 5.0])

    s = zeros(5, 5)
    s[1, 1] = 1

    s[2, 1] = λ[1] / (λ[1] - λ[2])
    s[3, 1] = λ[1] / (λ[1] - λ[3]) * λ[2] / (λ[2] - λ[3])
    s[4, 1] = λ[1] / (λ[1] - λ[4]) * λ[2] / (λ[2] - λ[4]) * λ[3] / (λ[3] - λ[4])
    s[5, 1] =
        λ[1] / (λ[1] - λ[5]) * λ[2] / (λ[2] - λ[5]) * λ[3] / (λ[3] - λ[5]) *
        λ[4] / (λ[4] - λ[5])

    s[2, 2] = 1
    s[3, 2] = λ[2] / (λ[2] - λ[3])
    s[4, 2] = λ[2] / (λ[2] - λ[4]) * λ[3] / (λ[3] - λ[4])
    s[5, 2] = λ[2] / (λ[2] - λ[5]) * λ[3] / (λ[3] - λ[5]) * λ[4] / (λ[4] - λ[5])

    s[3, 3] = 1
    s[4, 3] = λ[3] / (λ[3] - λ[4])
    s[5, 3] = λ[3] / (λ[3] - λ[5]) * λ[4] / (λ[4] - λ[5])

    s[4, 4] = 1
    s[5, 4] = λ[4] / (λ[4] - λ[5])

    s[5, 5] = 1

    return s
end

max_stars_repo_name: rannewan/RadioactiveDecayNetwork

max_issues_repo_name: rannewan/RadioactiveDecayNetwork

max_forks_repo_name: rannewan/RadioactiveDecayNetwork


text:
include("utils.jl")

build_so("008")

type Foo
   i::Cint
end

x = 13

y = ccall(
    (:bar, "binary/lib008.so"),
    Foo,
    (Cint, ),
    x
)

@assert y.i == x

y.i = 42

@assert y.i == 42


text:
function process_arg!(srcs_used, dos_ex, ::Type{GroupBy}, ex, index, primary)::Void
    # We assume that any Expr with head :. is just an attribute specification
    # TODO: check that this assumption actually holds
    is_predicate = isa(ex, Symbol) ? false : ifelse(ex.head == :., false, true)
    f_ex, ai = build_f_ex!(srcs_used, ex, index, primary)
    push!(dos_ex.args,
          Expr(:call, GroupBy, is_predicate, f_ex, ai))
    return
end


text:
################################################################################
#                         R MATRIX FUNCTIONS
################################################################################

################################################################################
function rmat_base_inc_b!(mx, θ, zrv, covstr)
        if covstr.repeated.covtype.s == :SI
            rmatp_si!(mx, θ, zrv)
        elseif covstr.repeated.covtype.s == :DIAG
            rmatp_diag!(mx, θ, zrv)
        elseif covstr.repeated.covtype.s == :AR
            rmatp_ar!(mx, θ, zrv)
        elseif covstr.repeated.covtype.s == :ARH
            rmatp_arh!(mx, θ, zrv)
        elseif covstr.repeated.covtype.s == :CSH
            rmatp_csh!(mx, θ, zrv)
        elseif covstr.repeated.covtype.s == :CS
            rmatp_cs!(mx, θ, zrv)
        elseif covstr.repeated.covtype.s == :ARMA
            rmatp_arma!(mx, θ, zrv)
        elseif covstr.repeated.covtype.s == :TOEPP
            rmatp_toepp!(mx, θ, zrv, covstr.repeated.covtype.p)
        elseif covstr.repeated.covtype.s == :TOEPHP
            rmatp_toephp!(mx, θ, zrv, covstr.repeated.covtype.p)
        elseif covstr.repeated.covtype.s == :FUNC
            covstr.repeated.covtype.f.xmat!(mx, θ, zrv, covstr.repeated.covtype.p)
        elseif covstr.repeated.covtype.s == :SPEXP
            rmatp_spexp!(mx, θ, zrv)
        end
end
################################################################################
################################################################################
@inline function rmat_base_inc!(mx, θ, covstr, block, sblock)
    zblock    = view(covstr.rz, block, :)
    @simd for i ∈ axes(sblock[end], 1)
        @inbounds rmat_base_inc_b!(view(mx, sblock[end][i],  sblock[end][i]), θ, view(zblock,  sblock[end][i], :), covstr)
    end
    mx
end
################################################################################
function rmatp_si!(mx, θ, ::AbstractMatrix)
    θsq = θ[1]*θ[1]
    @inbounds @simd for i ∈ axes(mx, 1)
            mx[i, i] += θsq
    end
    nothing
end
function rmatp_diag!(mx, θ, rz)
    #=@turbo=# @inbounds  for i ∈ axes(mx, 1)
        for c ∈ axes(θ, 1)
            mx[i, i] += rz[i, c] * θ[c] * θ[c]
        end
    end
    nothing
end
function rmatp_ar!(mx, θ, ::AbstractMatrix)
    rn  = size(mx, 1)
    de  = θ[1] ^ 2
    @inbounds @simd for m = 1:rn
        mx[m, m] += de
    end
    if rn > 1
        for m = 1:rn - 1
            @inbounds @simd for n = m + 1:rn
                mx[m, n] += de * θ[2] ^ (n - m)
            end
        end
    end
    nothing
end
function rmatp_arh!(mx, θ, rz)
    vec = tmul_unsafe(rz, θ)
    rn    = size(mx, 1)
    if rn > 1
        for m = 1:rn - 1
            @inbounds @simd for n = m + 1:rn
                 mx[m, n] += vec[m] * vec[n] * θ[end] ^ (n - m)
            end
        end
    end
    @inbounds  for m ∈ axes(mx, 1)
        mx[m, m] += vec[m] * vec[m]
    end
    nothing
end
function rmatp_cs!(mx, θ, ::AbstractMatrix)
    rn    = size(mx, 1)
    θsq   =  θ[1]*θ[1]
    θsqp  =  θsq*θ[2]
    @inbounds @simd for i = 1:size(mx, 1)
        mx[i, i] += θsq
    end
    if rn > 1
        for m = 1:rn - 1
            @inbounds @simd for n = m + 1:rn
                mx[m, n] += θsqp
            end
        end
    end
    nothing
end
function rmatp_csh!(mx, θ, rz)
    #vec   = rz * (θ[1:end-1])
    #=
    vec = zeros(T, size(rz, 1))
    #=@turbo=# @inbounds  for r ∈ axes(rz, 1)
        for i ∈ axes(rz, 2)
            vec[r] += rz[r, i] * θ[i]
        end
    end
    =#
    vec = tmul_unsafe(rz, θ)
    rn    = size(mx, 1)
    if rn > 1
        for m = 1:rn - 1
            @inbounds @simd for n = m + 1:rn
                mx[m, n] += vec[m] * vec[n] * θ[end]
            end
        end
    end
    #=@turbo=# @inbounds  for m ∈ axes(mx, 1)
        mx[m, m] += vec[m] * vec[m]
    end
    nothing
end
################################################################################
function rmatp_arma!(mx, θ, ::AbstractMatrix)
    rn  = size(mx, 1)
    de  = θ[1] ^ 2
    @inbounds @simd for m = 1:rn
        mx[m, m] += de
    end
    if rn > 1
        for m = 1:rn - 1
            @inbounds @simd for n = m + 1:rn
                mx[m, n] += de * θ[2] * θ[3] ^ (n - m - 1)
            end
        end
    end
    nothing
end
################################################################################
function rmatp_toepp!(mx, θ, ::AbstractMatrix, p::Int)
    de  = θ[1] ^ 2    #diagonal element
    s   = size(mx, 1) #size
    @inbounds @simd for i = 1:s
        mx[i, i] += de
    end
    if s > 1 && p > 1
        for m = 1:s - 1
            for n = m + 1:(m + p - 1 > s ? s : m + p - 1)
                @inbounds  mx[m, n] += de * θ[n - m + 1]
            end
        end
    end
    nothing
end
################################################################################
function rmatp_toephp!(mx, θ, rz, p::Int)
    l     = size(rz, 2)
    vec   = rz * (θ[1:l])
    s   = size(mx, 1) #size
    if s > 1 && p > 1
        for m = 1:s - 1
            for n = m + 1:(m + p - 1 > s ? s : m + p - 1)
                @inbounds  mx[m, n] += vec[m] * vec[n] * θ[n - m + l]
            end
        end
    end
    @inbounds @simd for m = 1:s
        mx[m, m] += vec[m] * vec[m]
    end
    nothing
end
################################################################################
#=
 Base.@propagate_inbounds function edistance(i::AbstractVector{T1}, j::AbstractVector{T2}) where T1 where T2
    length(i) == length(j) || error("length i not equal j")
    sum = zero(promote_type(T1, T2))
    for c = 1:length(i)
        sum += (i[c]-j[c])^2
    end
    return sqrt(sum)
end
=#
Base.@propagate_inbounds function edistance(mx::AbstractMatrix{T}, i::Int, j::Int) where T
    sum = zero(T)
    for c = 1:size(mx, 2)
        sum += (mx[i,c] - mx[j,c])^2
    end
    return sqrt(sum)
end
################################################################################
 Base.@propagate_inbounds function rmatp_spexp!(mx, θ, rz)
    σ²    = θ[1]^2
    θ     = exp(θ[2])
    rn    = size(mx, 1)
    @simd for i = 1:size(mx, 1)
        mx[i, i] += σ²
    end
    if rn > 1
        for m = 1:rn - 1
            @simd for n = m + 1:rn
                mx[m, n] += σ² * exp(-edistance(rz, m, n) / θ)
            end
        end
    end
    nothing
end


text:
using Yeppp
using Plots
using BenchmarkTools
using DelimitedFiles
using Printf

include("../src/SU_mod.jl")
include("../src/KaiserWindow_mod.jl")
include("../src/LanczosInterpolation_mod.jl")
include("../src/InvSmooth2d_mod.jl")
include("../src/lsq_cgsolver.jl")

using .LSQ_CGSolver
using .SU_mod
using .KaiserWindow_mod
using .LanczosInterpolation_mod
using .InvSmooth2d_mod

using Documenter, AFD2d_mod

makedocs(sitename="AFD2d_docs",modules=[Documenter,AFD2d_mod])


text:
"Holds the tableau of an variational partitioned additive Runge-Kutta method."
struct TableauVPARK{T} <: AbstractTableau{T}
    name::Symbol
    o::Int
    s::Int
    r::Int

    q::CoefficientsARK{T}
    p::CoefficientsARK{T}

    q̃::CoefficientsPRK{T}
    p̃::CoefficientsPRK{T}

    λ::CoefficientsMRK{T}

    d::Vector{T}

    function TableauVPARK{T}(name, o, s, r, q, p, q̃, p̃, λ, d) where {T}
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(r, Integer)
        @assert isa(o, Integer)

        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"

        @assert s==q.s==p.s==q̃.s==p̃.s==length(d)
        @assert r==q.r==p.r==q̃.r==p̃.r==λ.r

        new(name, o, s, r, q, p, q̃, p̃, λ, d)
    end

    function TableauVPARK{T}(name, o, s, r, q, p, q̃, p̃, λ) where {T}
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(r, Integer)
        @assert isa(o, Integer)

        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"

        @assert s==q.s==p.s==q̃.s==p̃.s
        @assert r==q.r==p.r==q̃.r==p̃.r==λ.r

        new(name, o, s, r, q, p, q̃, p̃, λ)
    end
end

function TableauVPARK(name::Symbol, order::Int,
                        a_q::Matrix{T}, a_p::Matrix{T},
                        α_q::Matrix{T}, α_p::Matrix{T},
                        a_q̃::Matrix{T}, a_p̃::Matrix{T},
                        α_q̃::Matrix{T}, α_p̃::Matrix{T},
                        b_q::Vector{T}, b_p::Vector{T},
                        β_q::Vector{T}, β_p::Vector{T},
                        c_q::Vector{T}, c_p::Vector{T},
                        c_λ::Vector{T}, d_λ::Vector{T},
                        d::Vector{T}) where {T <: Real}

    s = length(c_q)
    r = length(c_λ)

    @assert s > 0 "Number of stages s must be > 0"
    @assert r > 0 "Number of stages r must be > 0"

    @assert s==size(a_q,1)==size(a_q,2)==length(b_q)==length(c_q)
    @assert s==size(a_p,1)==size(a_p,2)==length(b_p)==length(c_p)
    @assert s==size(α_q,1)==size(α_p,1)
    @assert r==size(α_q,2)==size(α_p,2)
    @assert s==length(d)
    @assert r==length(c_λ)==length(d_λ)
    @assert r==size(a_q̃,1)==size(a_p̃,1)
    @assert s==size(a_q̃,2)==size(a_p̃,2)
    @assert r==size(α_q̃,1)==size(α_q̃,2)==length(β_q)
    @assert r==size(α_p̃,1)==size(α_p̃,2)==length(β_p)

    q = CoefficientsARK{T}(name, order, s, r, a_q, b_q, c_q, α_q, β_q)
    p = CoefficientsARK{T}(name, order, s, r, a_p, b_p, c_p, α_p, β_p)
    q̃ = CoefficientsPRK{T}(name, order, s, r, a_q̃, c_λ, α_q̃)
    p̃ = CoefficientsPRK{T}(name, order, s, r, a_p̃, c_λ, α_p̃)
    λ = CoefficientsMRK{T}(name, r, d_λ, c_λ)

    TableauVPARK{T}(name, order, s, r, q, p, q̃, p̃, λ, d)
end

function TableauVPARK(name::Symbol, order::Int,
                        a_q::Matrix{T}, a_p::Matrix{T},
                        α_q::Matrix{T}, α_p::Matrix{T},
                        a_q̃::Matrix{T}, a_p̃::Matrix{T},
                        α_q̃::Matrix{T}, α_p̃::Matrix{T},
                        b_q::Vector{T}, b_p::Vector{T},
                        β_q::Vector{T}, β_p::Vector{T},
                        c_q::Vector{T}, c_p::Vector{T},
                        c_λ::Vector{T}, d_λ::Vector{T}) where {T <: Real}

    s = length(c_q)
    r = length(c_λ)

    @assert s > 0 "Number of stages s must be > 0"
    @assert r > 0 "Number of stages r must be > 0"

    @assert s==size(a_q,1)==size(a_q,2)==length(b_q)==length(c_q)
    @assert s==size(a_p,1)==size(a_p,2)==length(b_p)==length(c_p)
    @assert s==size(α_q,1)==size(α_p,1)
    @assert r==size(α_q,2)==size(α_p,2)
    @assert r==length(c_λ)==length(d_λ)
    @assert r==size(a_q̃,1)==size(a_p̃,1)
    @assert s==size(a_q̃,2)==size(a_p̃,2)
    @assert r==size(α_q̃,1)==size(α_q̃,2)==length(β_q)
    @assert r==size(α_p̃,1)==size(α_p̃,2)==length(β_p)

    q = CoefficientsARK{T}(name, order, s, r, a_q, b_q, c_q, α_q, β_q)
    p = CoefficientsARK{T}(name, order, s, r, a_p, b_p, c_p, α_p, β_p)
    q̃ = CoefficientsPRK{T}(name, order, s, r, a_q̃, c_λ, α_q̃)
    p̃ = CoefficientsPRK{T}(name, order, s, r, a_p̃, c_λ, α_p̃)
    λ = CoefficientsMRK{T}(name, r, d_λ, c_λ)

    TableauVPARK{T}(name, order, s, r, q, p, q̃, p̃, λ)
end

# TODO function readTableauVPARKFromFile(dir::AbstractString, name::AbstractString)

"Parameters for right-hand side function of variational partitioned additive Runge-Kutta methods."
mutable struct ParametersVPARK{DT,TT,D,S,R,ϑT,FT,UT,GT,ϕT} <: Parameters{DT,TT}
    f_ϑ::ϑT
    f_f::FT
    f_u::UT
    f_g::GT
    f_ϕ::ϕT

    Δt::TT

    t_q::CoefficientsARK{TT}
    t_p::CoefficientsARK{TT}
    t_q̃::CoefficientsPRK{TT}
    t_p̃::CoefficientsPRK{TT}
    t_λ::CoefficientsMRK{TT}
    d_v::Vector{TT}

    t::TT
    q::Vector{DT}
    p::Vector{DT}
    λ::Vector{DT}

    function ParametersVPARK{DT,TT,D,S,R,ϑT,FT,UT,GT,ϕT}(f_ϑ, f_f, f_u, f_g, f_ϕ, Δt, t_q, t_p, t_q̃, t_p̃, t_λ, d_v) where {DT,TT,D,S,R,ϑT,FT,UT,GT,ϕT}
        # create solution vectors
        q = zeros(DT,D)
        p = zeros(DT,D)
        λ = zeros(DT,D)

        new(f_ϑ, f_f, f_u, f_g, f_ϕ, Δt,
            t_q, t_p, t_q̃, t_p̃, t_λ, d_v,
            zero(TT), q, p, λ)
    end
end

@doc raw"""
Variational partitioned additive Runge-Kutta integrator.

This integrator solves the following system of equations for the internal stages,
```math
\begin{align}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= 0 , & i &= 1, ..., r ,
\end{align}
```
with definitions
```math
\begin{align}
P_{n,i} &= \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
F_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
G_{n,i} &=           -  \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r ,
\end{align}
```
and update rule
```math
\begin{align}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} .
\end{align}
```
"""
struct IntegratorVPARK{DT, TT, ET <: IDAE{DT,TT},
                               PT <: ParametersVPARK{DT,TT},
                               ST <: NonlinearSolver{DT},
                               IT <: InitialGuessPODE{DT,TT}, D, S, R} <: AbstractIntegratorVSPARK{DT, TT}
    equation::ET
    tableau::TableauVPARK{TT}

    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheSPARK{DT, TT, D, S, R}
end

function IntegratorVPARK(equation::IDAE{DT,TT,FT,PT,UT,GT,ϕT,VT},
                         tableau::TableauVPARK{TT}, Δt::TT) where {DT,TT,FT,PT,UT,GT,ϕT,VT}
    D = equation.d
    S = tableau.s
    R = tableau.r

    N = 3*D*S + 3*D*R

    if isdefined(tableau, :d)
        N += D
        d_v = tableau.d
    else
        d_v = DT[]
    end

    # create params
    params = ParametersVPARK{DT,TT,D,S,R,FT,PT,UT,GT,ϕT}(
                                equation.ϑ, equation.f, equation.u, equation.g, equation.ϕ, Δt,
                                tableau.q, tableau.p, tableau.q̃, tableau.p̃, tableau.λ, d_v)

    # create solver
    solver = create_nonlinear_solver(DT, N, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCacheSPARK{DT, TT, D, S, R}()

    # create integrator
    IntegratorVPARK{DT, TT, typeof(equation), typeof(params), typeof(solver), typeof(iguess), D, S, R}(
                                        equation, tableau, params, solver, iguess, cache)
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,TT,D,S,R},
                                        params::ParametersVPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local tpᵢ::TT
    local tλᵢ::TT

    for i in 1:S
        for k in 1:D
            # copy x to Y, Z
            cache.Yi[i][k] = x[3*(D*(i-1)+k-1)+1]
            cache.Zi[i][k] = x[3*(D*(i-1)+k-1)+2]
            cache.Vi[i][k] = x[3*(D*(i-1)+k-1)+3]

            # compute Q and P
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
            cache.Pi[i][k] = params.p[k] + params.Δt * cache.Zi[i][k]
        end

        # compute f(X)
        tpᵢ = params.t + params.Δt * params.t_p.c[i]
        params.f_ϑ(tpᵢ, cache.Qi[i], cache.Vi[i], cache.Φi[i])
        params.f_f(tpᵢ, cache.Qi[i], cache.Vi[i], cache.Fi[i])

        cache.Φi[i] .-= cache.Pi[i]
    end

    for i in 1:R
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+2]
            cache.Λp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+3]

            # compute Q and V
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        tλᵢ = params.t + params.Δt * params.t_λ.c[i]
        params.f_u(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Up[i])
        params.f_g(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Gp[i])
        params.f_ϕ(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Φp[i])
    end

    if length(params.d_v) > 0
        for k in 1:D
            cache.μ[k] = x[3*D*S+3*D*R+k]
        end
    end
end

"Compute stages of variational partitioned additive Runge-Kutta methods."
@generated function function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersVPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    cache = IntegratorCacheSPARK{ST,TT,D,S,R}()

    quote
        compute_stages!(y, $cache, params)

        # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
        for i in 1:S
            for k in 1:D
                b[3*(D*(i-1)+k-1)+1] = - $cache.Yi[i][k]
                b[3*(D*(i-1)+k-1)+2] = - $cache.Zi[i][k]
                b[3*(D*(i-1)+k-1)+3] = - $cache.Φi[i][k]
                for j in 1:S
                    b[3*(D*(i-1)+k-1)+1] += params.t_q.a[i,j] * $cache.Vi[j][k]
                    b[3*(D*(i-1)+k-1)+2] += params.t_p.a[i,j] * $cache.Fi[j][k]
                end
                for j in 1:R
                    b[3*(D*(i-1)+k-1)+1] += params.t_q.α[i,j] * $cache.Up[j][k]
                    b[3*(D*(i-1)+k-1)+2] += params.t_p.α[i,j] * $cache.Gp[j][k]
                end
            end
        end

        # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
        for i in 1:R
            for k in 1:D
                b[3*D*S+3*(D*(i-1)+k-1)+1] = - $cache.Yp[i][k]
                b[3*D*S+3*(D*(i-1)+k-1)+2] = - $cache.Zp[i][k]
                b[3*D*S+3*(D*(i-1)+k-1)+3] = - $cache.Φp[i][k]
                for j in 1:S
                    b[3*D*S+3*(D*(i-1)+k-1)+1] += params.t_q̃.a[i,j] * $cache.Vi[j][k]
                    b[3*D*S+3*(D*(i-1)+k-1)+2] += params.t_p̃.a[i,j] * $cache.Fi[j][k]
                end
                for j in 1:R
                    b[3*D*S+3*(D*(i-1)+k-1)+1] += params.t_q̃.α[i,j] * $cache.Up[j][k]
                    b[3*D*S+3*(D*(i-1)+k-1)+2] += params.t_p̃.α[i,j] * $cache.Gp[j][k]
                end
            end
        end

        # compute b = - [Λ₁-λ]
        if params.t_λ.c[1] == 0
            for k in 1:D
                b[3*D*S+3*(k-1)+3] = - $cache.Λp[1][k] + params.λ[k]
            end
        end

        if length(params.d_v) > 0
            for i in 1:S
                for k in 1:D
                    b[3*(D*(i-1)+k-1)+3] -= $cache.μ[k] * params.d_v[i]
                end
            end

            for k in 1:D
                b[3*D*S+3*D*R+k] = 0
                for i in 1:S
                    b[3*D*S+3*D*R+k] -= $cache.Vi[i][k] * params.d_v[i]
                end
            end
        end
    end
end

function update_params!(params::ParametersVPARK, sol::AtomicSolutionPDAE)
    # set time for nonlinear solver and copy previous solution
    params.t  = sol.t
    params.q .= sol.q
    params.p .= sol.p
    params.λ .= sol.λ
end

function initialize!(int::IntegratorVPARK, sol::AtomicSolutionPDAE)
    sol.t̅ = sol.t - timestep(int)

    equation(int).v(sol.t, sol.q, sol.p, sol.v)
    equation(int).f(sol.t, sol.q, sol.p, sol.f)

    initialize!(int.iguess, sol.t, sol.q, sol.p, sol.v, sol.f,
                            sol.t̅, sol.q̅, sol.p̅, sol.v̅, sol.f̅)
end

function initial_guess!(int::IntegratorVPARK, sol::AtomicSolutionPDAE)
    for i in eachstage(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.p̃, int.cache.ṽ, int.cache.f̃,
                              tableau(int).q.c[i], tableau(int).p.c[i])

        for k in eachdim(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+1] = (int.cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+2] = (int.cache.p̃[k] - sol.p[k])/timestep(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+3] = int.cache.ṽ[k]
        end
    end

    for i in 1:pstages(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.p̃, int.cache.ṽ, int.cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in 1:ndims(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+1] = (int.cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+2] = (int.cache.p̃[k] - sol.p[k])/timestep(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+3] = 0
        end
    end

    if int.params.t_λ.c[1] == 0
        for k in eachdim(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(k-1)+3] = sol.λ[k]
        end
    end

    if isdefined(tableau(int), :d)
        for k in eachdim(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*ndims(int)*pstages(int)+k] = 0
        end
    end
end

function update_solution!(int::IntegratorVPARK{DT,TT}, sol::AtomicSolutionPDAE{DT,TT}) where {DT,TT}
    # compute final update
    update_solution!(sol.q, sol.q̃, int.cache.Vi, int.params.t_q.b, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.Fi, int.params.t_p.b, timestep(int))

    # compute projection
    update_solution!(sol.q, sol.q̃, int.cache.Up, int.params.t_q.β, timestep(int))
    update_solution!(sol.p, sol.p̃, int.cache.Gp, int.params.t_p.β, timestep(int))
    update_multiplier!(sol.λ, int.cache.Λp, int.params.t_λ.b)
end


text:
#========== master evaluator ==========#

"""
    Eval(fwd, rev)
    (e::Eval)(A,B) = fwd(A,B)

This holds the functions `$MAKE` which creates the output array
(before calling `threader($ACT!,...)` to fill it)
and the function `∇$MAKE` which does the reverse pass for differentiation.

It exists so that gradient hooks for various packages can be attached to this,
once. Then `$MAKE` need only be defined in local scope.
"""
struct Eval{F,R}
    fwd::F
    rev::R
end

(ev::Eval)(args...) = ev.fwd(args...)

#========== gradient hooks ==========#
# Macros like @adjoint need to be hidden behind include(), it seems:

using Requires

@init @require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" include("grad/zygote.jl")

@init @require Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" include("grad/tracker.jl")

@init @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" include("grad/reverse.jl")

#=
@init @requite Yota = "cd998857-8626-517d-b929-70ad188a48f0" begin
    using .Yota

#     for (n,A) in enumerate(store.arrays)
#         push!(evalex, quote
#             Yota.@diffrule  $make($(store.arrays...), $(store.scalars...))  $A  getindex($∇make(dy, $(store.arrays...), $(store.scalars...)), $n)
#         end)
#     end

end
=#

#========== not gradients ==========#

getonly(a::AbstractArray) = first(a) # just avoid first(::CuArray)

using Requires

@init @require CuArrays = "3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    using .CuArrays

    Tullio.threader(fun!::F, ::Type{T},
        Z::AbstractArray, As::Tuple, Is::Tuple, Js::Tuple,
        redfun, block=0, keep=nothing) where {F<:Function, T<:CuArray} =
        fun!(T, Z, As..., Is..., Js..., keep)

    Tullio.∇threader(fun!::F, ::Type{T},
        As::Tuple, Is::Tuple, Js::Tuple, block=0) where {F<:Function, T<:CuArray} =
        fun!(T, As..., Is..., Js...,)

    Tullio.getonly(a::CuArray) = sum(a)

end

# CUDA replaces CuArrays, Julia >=1.4 only.
@init @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" begin
    using .CUDA

    Tullio.threader(fun!::F, ::Type{T},
        Z::AbstractArray, As::Tuple, Is::Tuple, Js::Tuple,
        redfun, block=0, keep=nothing) where {F<:Function, T<:CuArray} =
        fun!(T, Z, As..., Is..., Js..., keep)

    Tullio.∇threader(fun!::F, ::Type{T},
        As::Tuple, Is::Tuple, Js::Tuple, block=0) where {F<:Function, T<:CuArray} =
        fun!(T, As..., Is..., Js...,)

    Tullio.getonly(a::CuArray) = sum(a)

end

@init @require LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" begin
    @require ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" begin
        # dual numbers + svec, should live in PaddedMatricesForwardDiff?
        # (And where would the conditional loading go, still here?)
        include("grad/avxdual.jl")
    end
end

#========== storage unwrapper ==========#

"""
    storage_type(adjoint(view(A,...))) == Array{Int,2}
    storage_type(A, B, C) == Array{Int,N} where N

Recursively unwraps wrappers, and combines approximately with `promote_type`.
(Used as the trait to send CuArray to KernelAbstractions and Array{Float or Int}
to LoopVectorization.)
"""
function storage_type(A::AbstractArray)
    P = parent(A)
    typeof(A) === typeof(P) ? typeof(A) : storage_type(P)
end
storage_type(A) = typeof(A)
storage_type(A, Bs...) = promote_storage(storage_type(A), storage_type(Bs...))
storage_type() = AbstractArray

promote_storage(::Type{A}, ::Type{B}) where {A <: Array{T,N}, B <: Array{S,M}} where {T,N,S,M} =
    N==M ? Array{promote_type(T,S), N} : Array{promote_type(T,S)}
promote_storage(::Type{A}, ::Type{B}) where {A <: Array{T,N}, B <: AbstractRange{S}} where {T,N,S} =
    N==1 ? Vector{promote_type(T,S)} : Array{promote_type(T,S)}
promote_storage(::Type{A}, ::Type{B}) where {A <: AbstractRange{T}, B <: Array{S,M}} where {T,S,M} =
    M==1 ? Vector{promote_type(T,S)} : Array{promote_type(T,S)}
promote_storage(T, S) = promote_type(T, S)

#========== fin ==========#


text:
############################################################
# joAbstractLinearOperator subtypes ########################
############################################################

############################################################
# joLinearOperator #########################################

export joLinearOperator, joLinearOperatorException

# type definition
"""
joLinearOperator is glueing type & constructor

!!! Do not use it to create the operators

"""
struct joLinearOperator{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    fop::Function              # forward
    fop_T::Nullable{Function}  # transpose
    fop_A::Nullable{Function}  # adjoint
    fop_C::Nullable{Function}  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
end

# type exception
struct joLinearOperatorException <: Exception
    msg :: String
end

############################################################
# joMatrix #################################################

export joMatrix, joMatrixException

# type definition
"""
joMatrix type

# TYPE PARAMETERS
- DDT::DataType : domain DataType
- RDT::DataType : range DataType

# FIELDS
- name::String : given name
- m::Integer : # of rows
- n::Integer : # of columns
- fop::Function : forward matrix
- fop_T::Function : transpose matrix
- fop_A::Function : adjoint matrix
- fop_C::Function : conj matrix
- iop::Nullable{Function} : inverse for fop
- iop_T::Nullable{Function} : inverse for fop_T
- iop_A::Nullable{Function} : inverse for fop_A
- iop_C::Nullable{Function} : inverse for fop_C

"""
struct joMatrix{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    fop::Function    # forward
    fop_T::Function  # transpose
    fop_A::Function  # adjoint
    fop_C::Function  # conj
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
end

# type exception
struct joMatrixException <: Exception
    msg :: String
end

############################################################
# joLinearFunction #########################################

export joLinearFunction, joLinearFunctionException

# type definition
"""
joLinearFunction type

# TYPE PARAMETERS
- DDT::DataType : domain DataType
- RDT::DataType : range DataType

# FIELDS
- name::String : given name
- m::Integer : # of rows
- n::Integer : # of columns
- fop::Function : forward function
- fop_T::Nullable{Function} : transpose function
- fop_A::Nullable{Function} : adjoint function
- fop_C::Nullable{Function} : conj function
- fMVok : whether fops are rady to handle mvec
- iop::Nullable{Function} : inverse for fop
- iop_T::Nullable{Function} : inverse for fop_T
- iop_A::Nullable{Function} : inverse for fop_A
- iop_C::Nullable{Function} : inverse for fop_C
- iMVok::Bool : whether iops are rady to handle mvec

"""
struct joLinearFunction{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    fop::Function              # forward
    fop_T::Nullable{Function}  # transpose
    fop_A::Nullable{Function}  # adjoint
    fop_C::Nullable{Function}  # conj
    fMVok::Bool                # forward can do mvec
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
    iMVok::Bool
end

# type exception
struct joLinearFunctionException <: Exception
    msg :: String
end

############################################################
# joCoreBlock ##############################################

export joCoreBlock, joCoreBlockException

# type definition
#struct joCoreBlock{DDT,RDT,N} <: joAbstractParallelableLinearOperator{DDT,RDT,N}
struct joCoreBlock{DDT,RDT} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    l::Integer
    ms::Vector{Integer}
    ns::Vector{Integer}
    mo::Vector{Integer}
    no::Vector{Integer}
    ws::Vector{Number}
    fop::Vector{joAbstractLinearOperator}
    fop_T::Vector{joAbstractLinearOperator}
    fop_A::Vector{joAbstractLinearOperator}
    fop_C::Vector{joAbstractLinearOperator}
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
end

# type exception
struct joCoreBlockException <: Exception
    msg :: String
end

############################################################
# joKron ###################################################

export joKron, joKronException

# type definition
#struct joKron{DDT,RDT,N} <: joAbstractParallelableLinearOperator{DDT,RDT,N}
struct joKron{DDT,RDT} <: joAbstractLinearOperator{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    l::Integer
    ms::Vector{Integer}
    ns::Vector{Integer}
    flip::Bool
    fop::Vector{joAbstractLinearOperator}
    fop_T::Vector{joAbstractLinearOperator}
    fop_A::Vector{joAbstractLinearOperator}
    fop_C::Vector{joAbstractLinearOperator}
    iop::Nullable{Function}
    iop_T::Nullable{Function}
    iop_A::Nullable{Function}
    iop_C::Nullable{Function}
end

# type exception
struct joKronException <: Exception
    msg :: String
end


text:
struct Params
    p::Dict{Symbol, String}
end

Params() = Params(Dict{Symbol, String}())

function Base.getproperty(params::Params, prop::Symbol)
    if prop in (:p,)
        getfield(params, :p)
    end
    return getindex(getfield(params, :p), prop)
end

function createparams(request::HTTP.Request, path::String)
    paths = URIs.splitpath(path)
    p_names_str = filter(str -> startswith(str, ':'), paths)
    p_positions = findall(x -> x ∈ p_names_str, paths)
    p_names = map(str -> Symbol(str[2:end]), p_names_str)
    p_values = URIs.splitpath(request.target)[p_positions]
    return Params(Dict(zip(p_names, p_values)))
end

struct Query
    q::Dict{Symbol, String}
end

Query() = Query(Dict{Symbol, String}())

function Base.getproperty(query::Query, prop::Symbol)
    if prop in (:q,)
        getfield(query, :q)
    end
    return getindex(getfield(query, :q), prop)
end

function createquery(request::HTTP.Request)
    dict = queryparams(URI(request.target))
    if isempty(dict)
        return Query()
    else
        return Query(Dict( zip(Symbol.(keys(dict)), values(dict)) ))
    end
end

const BodyTypes = Union{String, JSON3.Object}

struct Req{B<:BodyTypes}
    method::String
    target::String
    params::Params
    query::Query
    body::B
end


text:
int_rules_1_1_2_y = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.2.0.y*P(x)*(c*x)^m*(a+b*x^2)^p =#
    @apply_utils Antiderivative((~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~(p'), ~x) => (1 / 2) * Subst(Antiderivative((~x) ^ ((~m - 1) / 2) * SubstFor((~x) ^ 2, ~Pq, ~x) * (~a + ~b * ~x) ^ ~p, ~x), ~x, (~x) ^ 2) <-- FreeQ([~a, ~b, ~p], ~x) && (PolyQ(~Pq, (~x) ^ 2) && IntegerQ((~m - 1) / 2))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~(p'), ~x) => (1 / ~c) * Antiderivative((~c * ~x) ^ (~m + 1) * PolynomialQuotient(~Pq, ~x, ~x) * (~a + ~b * (~x) ^ 2) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && EqQ(Coeff(~Pq, ~x, 0), 0))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~P2 * (~a + ~(b') * (~x) ^ 2) ^ ~(p'), ~x) => With([f = Coeff(~P2, ~x, 0), g = Coeff(~P2, ~x, 1), h = Coeff(~P2, ~x, 2)], (h * (~c * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ 2) ^ (~p + 1)) / (~b * ~c * (~m + 2 * ~p + 3)) <-- EqQ(g, 0) && EqQ(~a * h * (~m + 1) - ~b * f * (~m + 2 * ~p + 3), 0)) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (PolyQ(~P2, ~x, 2) && NeQ(~m, -1))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~(p'), ~x) => Antiderivative(ExpandIntegrand((~c * ~x) ^ ~m * ~Pq * (~a + ~b * (~x) ^ 2) ^ ~p, ~x), ~x) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (PolyQ(~Pq, ~x) && IGtQ(~p, -2))
    @apply_utils Antiderivative((~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~p, ~x) => With([A = Coeff(~Pq, ~x, 0), Q = PolynomialQuotient(~Pq - Coeff(~Pq, ~x, 0), (~x) ^ 2, ~x)], (A * (~x) ^ (~m + 1) * (~a + ~b * (~x) ^ 2) ^ (~p + 1)) / (~a * (~m + 1)) + (1 / (~a * (~m + 1))) * Antiderivative((~x) ^ (~m + 2) * (~a + ~b * (~x) ^ 2) ^ ~p * (~a * (~m + 1) * Q - A * ~b * (~m + 2 * (~p + 1) + 1)), ~x)) <-- FreeQ([~a, ~b], ~x) && (PolyQ(~Pq, (~x) ^ 2) && (IntegerQ(~m / 2) && (ILtQ((~m + 1) / 2 + ~p, 0) && LtQ(~m + Expon(~Pq, ~x) + 2 * ~p + 1, 0))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~a + ~b * (~x) ^ 2, ~x), f = Coeff(PolynomialRemainder(~Pq, ~a + ~b * (~x) ^ 2, ~x), ~x, 0), g = Coeff(PolynomialRemainder(~Pq, ~a + ~b * (~x) ^ 2, ~x), ~x, 1)], ((~c * ~x) ^ ~m * (~a + ~b * (~x) ^ 2) ^ (~p + 1) * (~a * g - ~b * f * ~x)) / (2 * ~a * ~b * (~p + 1)) + (~c / (2 * ~a * ~b * (~p + 1))) * Antiderivative((~c * ~x) ^ (~m - 1) * (~a + ~b * (~x) ^ 2) ^ (~p + 1) * ExpandToSum((2 * ~a * ~b * (~p + 1) * ~x * Q - ~a * g * ~m) + ~b * f * (~m + 2 * ~p + 3) * ~x, ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (PolyQ(~Pq, ~x) && (LtQ(~p, -1) && GtQ(~m, 0)))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient((~c * ~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ 2, ~x), f = Coeff(PolynomialRemainder((~c * ~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ 2, ~x), ~x, 0), g = Coeff(PolynomialRemainder((~c * ~x) ^ ~m * ~Pq, ~a + ~b * (~x) ^ 2, ~x), ~x, 1)], ((~a * g - ~b * f * ~x) * (~a + ~b * (~x) ^ 2) ^ (~p + 1)) / (2 * ~a * ~b * (~p + 1)) + (1 / (2 * ~a * (~p + 1))) * Antiderivative((~c * ~x) ^ ~m * (~a + ~b * (~x) ^ 2) ^ (~p + 1) * ExpandToSum(2 * ~a * (~p + 1) * (~c * ~x) ^ -(~m) * Q + f * (2 * ~p + 3) * (~c * ~x) ^ -(~m), ~x), ~x)) <-- FreeQ([~a, ~b, ~c], ~x) && (PolyQ(~Pq, ~x) && (LtQ(~p, -1) && ILtQ(~m, 0)))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~a + ~b * (~x) ^ 2, ~x), f = Coeff(PolynomialRemainder(~Pq, ~a + ~b * (~x) ^ 2, ~x), ~x, 0), g = Coeff(PolynomialRemainder(~Pq, ~a + ~b * (~x) ^ 2, ~x), ~x, 1)], (-((~c * ~x) ^ (~m + 1)) * (f + g * ~x) * (~a + ~b * (~x) ^ 2) ^ (~p + 1)) / (2 * ~a * ~c * (~p + 1)) + (1 / (2 * ~a * (~p + 1))) * Antiderivative((~c * ~x) ^ ~m * (~a + ~b * (~x) ^ 2) ^ (~p + 1) * ExpandToSum(2 * ~a * (~p + 1) * Q + f * (~m + 2 * ~p + 3) + g * (~m + 2 * ~p + 4) * ~x, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (PolyQ(~Pq, ~x) && (LtQ(~p, -1) && Not(GtQ(~m, 0))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~m * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~p, ~x) => With([Q = PolynomialQuotient(~Pq, ~c * ~x, ~x), R = PolynomialRemainder(~Pq, ~c * ~x, ~x)], (R * (~c * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ 2) ^ (~p + 1)) / (~a * ~c * (~m + 1)) + (1 / (~a * ~c * (~m + 1))) * Antiderivative((~c * ~x) ^ (~m + 1) * (~a + ~b * (~x) ^ 2) ^ ~p * ExpandToSum(~a * ~c * (~m + 1) * Q - ~b * R * (~m + 2 * ~p + 3) * ~x, ~x), ~x)) <-- FreeQ([~a, ~b, ~c, ~p], ~x) && (PolyQ(~Pq, ~x) && (LtQ(~m, -1) && (IntegerQ(2 * ~p) || NeQ(Expon(~Pq, ~x), 1))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~(p'), ~x) => With([q = Expon(~Pq, ~x)], (Coeff(~Pq, ~x, q) / (~c) ^ q) * Antiderivative((~c * ~x) ^ (~m + q) * (~a + ~b * (~x) ^ 2) ^ ~p, ~x) + (1 / (~c) ^ q) * Antiderivative((~c * ~x) ^ ~m * (~a + ~b * (~x) ^ 2) ^ ~p * ExpandToSum((~c) ^ q * ~Pq - Coeff(~Pq, ~x, q) * (~c * ~x) ^ q, ~x), ~x) <-- EqQ(q, 1) || EqQ(~m + q + 2 * ~p + 1, 0)) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && Not(IGtQ(~m, 0) && ILtQ(~p + 1 / 2, 0)))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * ~Pq * (~a + ~(b') * (~x) ^ 2) ^ ~p, ~x) => With([q = Expon(~Pq, ~x), f = Coeff(~Pq, ~x, Expon(~Pq, ~x))], (f * (~c * ~x) ^ ((~m + q) - 1) * (~a + ~b * (~x) ^ 2) ^ (~p + 1)) / (~b * (~c) ^ (q - 1) * (~m + q + 2 * ~p + 1)) + (1 / (~b * (~m + q + 2 * ~p + 1))) * Antiderivative((~c * ~x) ^ ~m * (~a + ~b * (~x) ^ 2) ^ ~p * ExpandToSum((~b * (~m + q + 2 * ~p + 1) * ~Pq - ~b * f * (~m + q + 2 * ~p + 1) * (~x) ^ q) - ~a * f * ((~m + q) - 1) * (~x) ^ (q - 2), ~x), ~x) <-- GtQ(q, 1) && NeQ(~m + q + 2 * ~p + 1, 0)) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (PolyQ(~Pq, ~x) && (Not(IGtQ(~m, 0)) || IGtQ(~p + 1 / 2, -1)))
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.2 Quadratic/1.1.2.y P(x) (c x)^m (a+b x^2)^p.jl


text:
using Test
using TimeSeriesEcon

@testset "FPConst" begin
    @test 1U == ii(1)
    @test 2000Y == yy(2000)
    @test 1999Q1 == qq(1999, 1)
    @test 1999Q2 == qq(1999, 2)
    @test 1999Q3 == qq(1999, 3)
    @test 1999Q4 == qq(1999, 4)
    @test 1988M12 == mm(1988, 12)
    @test 1988M11 == mm(1988, 11)
    @test 1988M10 == mm(1988, 10)
    @test 1988M9 == mm(1988, 9)
    @test 1988M8 == mm(1988, 8)
    @test 1988M7 == mm(1988, 7)
    @test 1988M6 == mm(1988, 6)
    @test 1988M5 == mm(1988, 5)
    @test 1988M4 == mm(1988, 4)
    @test 1988M3 == mm(1988, 3)
    @test 1988M2 == mm(1988, 2)
    @test 1988M1 == mm(1988, 1)
end

ts_u = TSeries(1:5)
ts_m = TSeries(mm(2018, 1), collect(1.0:12.0))
ts_q = TSeries(qq(2018, 1):qq(2020, 4), collect(1:12))
ts_y = TSeries(yy(2018), collect(1:12))

@testset "TSeries: Construction" begin

    @test ts_u.firstdate == ii(1)
    @test ts_u.values == 1:5

    @test ts_m.firstdate == mm(2018, 1)
    @test ts_m.values == collect(1.0:12.0)

    @test ts_q.firstdate == qq(2018, 1)
    @test ts_q.values == collect(1.0:12.0)

    @test ts_y.firstdate == yy(2018)
    @test ts_y.values == collect(1.0:12.0)

    # Make sure if lengths are different we get an error
    @test_throws ArgumentError TSeries(1U:5U, 1:6)

    let t = TSeries(2U, rand(5))
        @test firstdate(t) == 2U
        @test lastdate(t) == 6U
        @test length(t.values) == 5
        @test length(t) == 5
    end

    let t = TSeries(1991Q1:1992Q4), s = TSeries(2222Y:2225Y, undef), r = TSeries(1006M3:1009M5, 0.3), e = TSeries(2000Y:1995Y, 7)
        @test isempty(e)
        @test length(r) == 10+12+12+5
        @test length(t) == length(t.values) == 8
        @test firstindex(t) == firstdate(t) == 1991Q1
        @test lastindex(t) == lastdate(t) == 1992Q4
    end

end

@testset "Int indexing" begin
    let t = TSeries(4U:8U, rand(5))
        @test t.firstdate == 4U && lastdate(t) == 8U
        # test access
        @test t[1] isa Number
        @test t[1] == t.values[1]
        @test t[2:4] isa TSeries{frequencyof(t),Vector{Float64}}
        @test t[2:4].values == t.values[2:4]
        @test t[[1,3,4]] isa Vector{Float64}
        @test t[[1,3,4]] == t.values[[1,3,4]]
        # test assignment
        @test begin
            t[2:4] = 2.5
            t.values[2:4] == fill(2.5, 3)
        end
        @test 5 == (t[3] = 5)
        @test t.values == [first(t), 2.5, 5.0, 2.5, last(t)]
    end
end

@testset "Views" begin
    let t = TSeries(2010M1, rand(20))
        @test axes(t) == (2010M1 - 1 .+ (1:20),)
        @test Base.axes1(t) == 2010M1 - 1 .+ (1:20)
        z = similar(t)
        @test z isa typeof(t)
        @test z.firstdate == t.firstdate
        @test z != t
        z = copy(t)
        @test z == t
        z[3:5] += 0.2
        z[4:5] = [3,4]
        @test z != t
        z = view(t, 2:5)
        c = view(t, 2010M2:2010M5)
        @test c == z
        @test z == t[2:5]
        z[[1,3]] += [0.5, 0.5]
        z[[2,4]] = [1,1.5]
        @test z == t[2:5]
        @test c == z
    end
end

@testset "MIT arithmetic" begin
    @test 5U < 8U
    @test 5U <= 8U
    @test 5U <= 5U
    @test 5U >= 5U
    @test 5U == 5U
    @test 8U >= 5U
    @test 8U > 5U
    @test 2001Q1 >= 2000Q3
    @test_throws ArgumentError 1M1 > 2Q1
    @test_throws ArgumentError 1M1 <= 2Q1
    @test 1M1 != 2Q1

    @test 2001Y + 5 == 2006Y
    @test 6 + 2001Q3 == 2003Q1
    @test 2003Q1 - 2001Q3 == 6
    @test 2003Q1 - 6 == 2001Q3
    @test_throws ArgumentError 6 - 2003Q1
    @test_throws ArgumentError 2003Q1 + 2003Q1
    @test_throws ArgumentError 2003Q1 + 2003Y
end

@testset "show" begin
    for (nrow, fd) = zip([3, 4, 5, 6, 7, 8, 22, 23, 24, 25, 26, 30], Iterators.cycle((qq(2010, 1), mm(2010, 1), yy(2010), ii(1))))
        let io = IOBuffer()
            t = TSeries(fd, rand(24))
            show(IOContext(io, :displaysize => (nrow, 80)), MIME"text/plain"(), t)
            @test length(readlines(seek(io, 0))) == max(2, min(length(t) + 1, nrow - 3))
        end
    end
end

@testset "frequencyof" begin
    @test frequencyof(qq(2000, 1)) == Quarterly
    @test frequencyof(mm(2000, 1)) == Monthly
    @test frequencyof(yy(2000)) == Yearly
    @test frequencyof(ii(1)) == Unit
    @test frequencyof(qq(2001, 1):qq(2002, 1)) == Quarterly
    @test frequencyof(TSeries(yy(2000), zeros(5))) == Yearly
end

@testset "TSeries: Broadcasting" begin

    allones = TSeries(2020Q1, ones(10))
    @test 5 - allones == 5 .- allones

    # -----------------------------------
    tsbc = TSeries(2020M1, ℯ * ones(12))

    @test Base.BroadcastStyle(typeof(tsbc)) ==  Base.Broadcast.ArrayStyle{TSeries}()

    @test TimeSeriesEcon.find_tseries(Broadcast.Broadcasted(-, (5, tsbc))) == tsbc
    @test TimeSeriesEcon.find_tseries(tsbc) == tsbc
    @test TimeSeriesEcon.find_tseries(Any, tsbc) == tsbc
    @test TimeSeriesEcon.find_tseries(tsbc, nothing) == tsbc

    @test firstdate(similar(tsbc)) == 2020M1
    @test size(similar(tsbc)) == size(tsbc)

    @test log.(tsbc) == TSeries(2020M1, ones(12))
    @test tsbc.firstdate == mm(2020, 1)
    @test exp.(log.(tsbc)) == tsbc
    @test tsbc .* 0 .+ 100 == TSeries(mm(2020, 1), 100 * ones(12))

    tsbc = log.(tsbc) + 99 # 1 + 99
    @test tsbc == TSeries(2020M1, 100 * ones(12))
end

@testset "TSeries: Index using end" begin
    @test ts_m[end] == 12
    @test ts_m[firstdate(ts_m):end] == ts_m
end

@testset "TSeries: Monthly Access" begin
    @test ts_m[mm(2018, 1):mm(2018, 12)] == ts_m
    @test ts_m[mm(2018, 1):mm(2018, 12)].firstdate == mm(2018, 1)

    # access outside of ts boundaries
    @test ts_m[mm(2017, 1):mm(2019, 12)] == ts_m
    @test ts_m[mm(2017, 1):mm(2019, 12)].firstdate == ts_m.firstdate

    # partially out of boundary
    @test ts_m[mm(2017, 1):mm(2018, 6)] == ts_m[mm(2018, 1):mm(2018, 6)]
    @test ts_m[mm(2017, 1):mm(2018, 6)] == ts_m[mm(2018, 1):mm(2018, 6)]

    @test ts_m[mm(2017, 1):mm(2018, 6)] == ts_m[1:6]

    @test ts_m[mm(2018, 6):mm(2019, 12)] == ts_m[mm(2018, 6):mm(2018, 12)]
    @test ts_m[mm(2018, 6):mm(2019, 12)] == ts_m[6:12]

    # fully out of boundary
    @test ts_m[mm(2017, 1)] === nothing
    @test ts_m[mm(2017, 1):mm(2017, 3)] === nothing
end

@testset "TSeries: Quarterly Access" begin
    @test ts_q[qq(2018, 1):qq(2020, 4)] == ts_q

    # access outside of ts boundaries
    @test ts_q[qq(2017, 1):qq(2021, 4)] == ts_q

    # partially out of boundary
    @test ts_q[qq(2017, 1):qq(2018, 4)] == ts_q[qq(2018, 1):qq(2018, 4)]
    @test ts_q[qq(2017, 1):qq(2018, 4)] == ts_q[1:4]

    @test ts_q[qq(2018, 4):qq(2021, 4)] == ts_q[qq(2018, 4):qq(2020, 4)]
    @test ts_q[qq(2018, 4):qq(2021, 4)] == ts_q[4:12]

    # fully out of boundary
    @test ts_q[qq(2017, 1)] == nothing
    @test ts_q[qq(2017, 1):qq(2017, 3)] == nothing
end

@testset "TSeries: Yearly Access" begin
    @test ts_y[yy(2018):yy(2029)] == ts_y

    # access outside of ts boundaries
    @test ts_y[yy(2017):yy(2017) + 100] == ts_y

    # partially out of boundary
    @test ts_y[yy(2017):yy(2018)] == ts_y[yy(2018):yy(2018)]
    @test ts_y[yy(2017):yy(2021)] == ts_y[1:4]

    @test ts_y[yy(2018):yy(2100)] == ts_y[yy(2018):yy(2029)]
    @test ts_y[yy(2021):yy(2100)] == ts_y[4:12]

    # fully out of boundary
    @test ts_y[yy(2017)] == nothing
    @test ts_y[yy(2010):yy(2017)] == nothing
end

# ts_m = TSeries(mm(2018, 1), collect(1.0:12.0))
@testset "TSeries: Monthly Setting" begin
    begin
        ts_m = TSeries(mm(2018, 1):mm(2018, 12), collect(1.0:12.0))
        ts_m[mm(2019, 2):mm(2019, 4)] = 1;
        @test ts_m[mm(2019, 2):mm(2019, 4)].values == [1, 1, 1]
        @test ts_m.firstdate == mm(2018, 1)
        @test isequal(ts_m.values, vcat(collect(1.0:12.0), [NaN], [1, 1, 1]))
    end

    begin
        ts_m = TSeries(mm(2018, 1):mm(2018, 12), collect(1.0:12.0))
        ts_m[mm(2017, 10):mm(2017, 11)] = 1;
        @test ts_m[mm(2017, 10):mm(2017, 11)].values == [1, 1]
        @test ts_m.firstdate == mm(2017, 10)
        @test isequal(ts_m.values, vcat([1, 1], [NaN], collect(1.0:12.0)))
    end

    begin
        ts_m = TSeries(mm(2018, 1):mm(2018, 12), collect(1.0:12.0))
        ts_m[mm(2019, 2):mm(2019, 4)] = [9, 10, 11];
        @test ts_m[mm(2019, 2):mm(2019, 4)].values == [9, 10, 11]
        @test ts_m.firstdate == mm(2018, 1)
        @test isequal(ts_m.values, vcat(collect(1.0:12.0), [NaN], [9, 10, 11]))
    end

    begin
        ts_m = TSeries(mm(2018, 1):mm(2018, 12), collect(1.0:12.0))
        ts_m[mm(2017, 10):mm(2017, 11)] = [9, 10];
        @test ts_m[mm(2017, 10):mm(2017, 11)].values == [9, 10]
        @test ts_m.firstdate == mm(2017, 10)
        @test isequal(ts_m.values, vcat([9, 10], [NaN], collect(1.0:12.0)))
    end
end

@testset "TSeries: Addition" begin
    x = TSeries(ii(1), [7, 7, 7])
    y = TSeries(ii(3), [2, 4, 5])
    @test x + y == TSeries(ii(3), [9])

    x = TSeries(ii(1), [7, 7, 7])
    y = TSeries(ii(2), [2, 4, 5])
    @test x + y == TSeries(ii(2), [9, 11])
end

@testset "TSeries: Iris related" begin
    # IRIS based assignment of values from other TSeries
    x = TSeries(qq(2020, 1), zeros(3));
    y = TSeries(qq(2020, 1), ones(3));
    x[qq(2020, 1):qq(2020, 2)] = y;
    @test x == TSeries(qq(2020, 1), [1, 1, 0])

    # IRIS related: shift
    x = TSeries(qq(2020, 1), zeros(3));
    @test shift(x, 1) == TSeries(qq(2019, 4), zeros(3))

    shift!(x, 1)
    @test x == TSeries(qq(2019, 4), zeros(3))

    # IRIS related: nanrm!
    x = TSeries(qq(2020, 1), [NaN, 123, NaN]);
    nanrm!(x)
    @test x == TSeries(qq(2020, 2), [123])

    # TODO
    # - pct
    # - apct

end

@testset "TSeries: firstdate & lastdate" begin
    x = TSeries(qq(2020, 1), zeros(4));
    @test firstdate(x) == qq(2020, 1)
    @test lastdate(x) == qq(2020, 4)
end

@testset "MIT constructors: mm, qq, yy, ii" begin
    @test mm(2020, 1) == MIT{Monthly}(2020 * 12)
    @test qq(2020, 1) == MIT{Quarterly}(2020 * 4)
    @test yy(2020) == MIT{Yearly}(2020)
    @test ii(2020) == MIT{Unit}(2020)
end

@testset "MIT: year, period" begin
    @test year(mm(2020, 12)) == 2020
    @test period(mm(2020, 12)) == 12
end

@testset "MIT: mitrange" begin
    @test mitrange(TSeries(qq(2020, 1), ones(4))) == qq(2020, 1):qq(2020, 4)
end

@testset "MIT: ppy" begin
    @test ppy(Quarterly) == 4
    @test ppy(qq(2020, 1)) == 4
    @test ppy(TSeries(qq(2020, 1), ones(1))) == 4

end

@testset "axes of range" begin
    @test axes(1U:5U) == axes(1:5)
    @test Base.axes1(2020Y:2030Y) == Base.OneTo(11)
end

@testset "recursive" begin
    ts = TSeries(1U, zeros(0))
    ts[1U] = ts[2U] = 1.0
    @rec 3U:10U ts[t] = ts[t-1]+ts[t-2]
    @test ts.values == [1.0,1,2,3,5,8,13,21,34,55]
    t = zeros(10,7)
    r = rand(1, 7)
    t[1, :] = r
    @rec 2:10 t[s,:] = t[s-1,:] .* s
    @test t ≈ factorial.(1:10) * r
end


text:
using LinearAlgebra
include("poisson_linear_system.jl")

function step_jacobi!(unext, uprev, rhs, stepsize)
    numpts = length(unext)
    for i = 2:numpts-1
        unext[i] = 0.5 * (uprev[i-1] + uprev[i+1] + stepsize^2 * rhs[i])
    end
end

function run_jacobi_iterations(u0, rhs, stepsize, numiter)
    uprev = copy(u0)
    unext = copy(u0)
    for iter = 1:numiter
        step_jacobi!(unext, uprev, rhs, stepsize)
        uprev .= unext
    end
    return unext
end

function run_jacobi_measure_error(u0, rhs, stepsize, numiter, uexact)
    numpts = length(u0)
    uprev = copy(u0)
    unext = copy(u0)
    err = zeros(numiter)
    for iter = 1:numiter
        step_jacobi!(unext, uprev, rhs, stepsize)
        uprev .= unext
        err[iter] = norm(unext - uexact) / numpts
    end
    return unext, err
end

function solve_direct(rhs, bcleft, bcright, stepsize)
    numpts = length(rhs)
    op = -poisson_linear_system(numpts)
    copyrhs = copy(rhs)
    copyrhs[1] = bcleft / stepsize^2
    copyrhs[numpts] = bcright / stepsize^2
    sol = stepsize^2 * (op \ copyrhs)
    return sol
end

function jacobi_error_vs_iterations(numiter, numelements, wavenumber)
    xL = -1.0
    xR = +1.0
    bcleft = bcright = 0.0

    numpts = numelements + 1
    stepsize = (xR - xL) / numelements
    xrange = range(xL, stop = xR, length = numpts)
    rhs = cos.(wavenumber * xrange)

    u0 = zeros(numpts)
    udirect = solve_direct(rhs, bcleft, bcright, stepsize)

    ujacobi, err = run_jacobi_measure_error(u0, rhs, stepsize, numiter, udirect)

    return err
end

# numelements = [4, 8, 16, 32]
# numiter = numelements .^ 2
# wavenumber = 2pi
# err = [
#     jacobi_error_vs_iterations(numiter[i], numelements[i], wavenumber) for
#     i = 1:length(numelements)
# ]
#
#
# using Plots
# plot( err[1], yscale = :log10, label = string(numelements[1]))
# plot!(err[2], yscale = :log10, label = string(numelements[2]))
# plot!(err[3], yscale = :log10, label = string(numelements[3]))
# plot!(err[4], yscale = :log10, label = string(numelements[4]))
# plot!(err[5], yscale = :log10, label = string(numelements[5]))

max_stars_repo_name: ArjunNarayanan/ExploreJacobi

max_issues_repo_name: ArjunNarayanan/ExploreJacobi

max_forks_repo_name: ArjunNarayanan/ExploreJacobi


text:
using CarlemanLinearization

"""
    CARLIN <: AbstractContinuousPost

Implementation of the reachability method using Carleman linearization from [1].

[1] Forets, Marcelo, and Christian Schilling. "Reachability of weakly nonlinear systems using Carleman linearization."
    International Conference on Reachability Problems. Springer, Cham, 2021.
"""
Base.@kwdef struct CARLIN <: AbstractContinuousPost
    N::Int=2 # order of the algorithm
    compress::Bool=false # choose to use compressed Kronecker form
    δ::Float64=0.1 # step size for the linear reachability solver
    bloat::Bool=true # choose to include the error estimate in the result
    resets::Union{Int,Vector{Float64}}=0 # choose the number of resets (equal spacing) or a vector specifying the reset times within the time interval [0, T]
end

include("kronecker.jl")
include("post.jl")
include("reach.jl")


text:
# This file contains functions for the Profit General Direct Approach DEA model
"""
    ProfitGDADEAModel
An data structure representing a profit General Direct Approach DEA model.
"""
struct ProfitGDADEAModel <: AbstractProfitDEAModel
    n::Int64
    m::Int64
    s::Int64
    measure::Symbol
    monetary::Bool
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    techeff::Vector
    alloceff::Vector
    normalization::Vector
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deaprofitgda(X, Y, W, P, measure)
Compute profit efficiency using data envelopment analysis General Direct Approach model for
inputs `X`, outputs `Y`, price of inputs `W`, price of outputs `P`, and efficiency `measure`.

# Measure specification:
- `:ERG`: Enhanced Russell Graph (or Slack Based Measure (SBM)).

# Optional Arguments
- `atol=1e-6`: tolerance for DMU to be considered efficient.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2; 4; 8; 12; 6; 14; 14; 9.412];

julia> Y = [1; 5; 8; 9; 3; 7; 9; 2.353];

julia> W = [1; 1; 1; 1; 1; 1; 1; 1];

julia> P = [2; 2; 2; 2; 2; 2; 2; 2];

julia> deaprofitgda(X, Y, W, P, :ERG)
General Direct Approach Profit DEA Model
DMUs = 8; Inputs = 1; Outputs = 1
Returns to Scale = VRS
Associated efficiency measure = ERG
──────────────────────────────────
     Profit  Technical  Allocative
──────────────────────────────────
1  4.0        0.0         4.0
2  0.5        0.0         0.5
3  0.0        0.0         0.0
4  0.166667   0.0         0.166667
5  0.8        0.6         0.2
6  0.571429   0.52381     0.047619
7  0.285714   0.142857    0.142857
8  0.949449   0.8         0.149449
──────────────────────────────────
```
"""
function deaprofitgda(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector},
    W::Union{Matrix,Vector}, P::Union{Matrix,Vector},
    measure::Symbol;
    monetary::Bool = false, atol::Float64 = 1e-6,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::ProfitGDADEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    nw, mw = size(W, 1), size(W, 2)
    np, sp = size(P, 1), size(P, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nw != nx
        throw(DimensionMismatch("number of rows in W and X ($nw, $nx) are not equal"));
    end
    if np != ny
        throw(DimensionMismatch("number of rows in P and Y ($np, $ny) are not equal"));
    end
    if mw != m
        throw(DimensionMismatch("number of columns in W and X ($mw, $m) are not equal"));
    end
    if sp != s
        throw(DimensionMismatch("number of columns in P and Y ($sp, $s) are not equal"));
    end

    n = nx

    # Get efficiency measure model
    if measure == :ERG
        measuremodel = deaprofiterg(X, Y, W, P, optimizer = optimizer)
        measuretechnicalmodel = deaerg(X, Y, rts = :VRS, optimizer = optimizer)
    else
        throw(ArgumentError("Invalid efficiency `measure`"));
    end

    normmeasure = normfactor(measuremodel)
    pimeasure = efficiency(measuremodel)  .* normmeasure
    timeasure = efficiency(measuremodel, :Technical)
    plambdaeff = peersmatrix(measuremodel)
    Xtarget = targets(measuremodel, :X)
    Ytarget = targets(measuremodel, :Y)
    slackXmeasure = slacks(measuretechnicalmodel, :X)
    slackYmeasure = slacks(measuretechnicalmodel, :Y)
    Xtargetti = targets(measuretechnicalmodel, :X)
    Ytargetti = targets(measuretechnicalmodel, :Y)

    # Calculate slack as the difference between target and observed
    slackXmeasure = X .- Xtargetti
    slackYmeasure = Ytargetti .- Y

    # Get profit inefficiency of projection (target)
    if measure == :ERG
        targetmodel = deaprofiterg(Xtargetti, Ytargetti, W, P, optimizer = optimizer)
    end

    pitarget = efficiency(targetmodel) .* normfactor(targetmodel)

    monetaryti = sum(P .* slackYmeasure, dims=2).+ sum(W .* slackXmeasure, dims = 2)

    # Normalization factor (equation 13.13)
    normalization = zeros(n)
    for i=1:n
        if abs(timeasure[i] ) <= atol
            # For technical efficient DMU's
            normalization[i] = normmeasure[i]
        else
            # For inefficient DMU's
            normalization[i] = monetaryti[i] ./ timeasure[i]
        end
    end

    if monetary
        pefficiency = pimeasure
        techefficiency = timeasure .* normalization
        allocefficiency = pitarget
    else
        pefficiency = pimeasure ./ normalization
        techefficiency = timeasure
        allocefficiency = pitarget ./ normalization
    end

    return ProfitGDADEAModel(n, m, s, measure, monetary, names, pefficiency, plambdaeff, techefficiency, allocefficiency, normalization, Xtarget, Ytarget)
end

function Base.show(io::IO, x::ProfitGDADEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    techeff = efficiency(x, :Technical)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "General Direct Approach Profit DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Returns to Scale = VRS")
        print(io, "\n")
        print(io, "Associated efficiency measure = ", string(x.measure))
        print(io, "\n")
        show(io, CoefTable(hcat(eff, techeff, alloceff), ["Profit", "Technical", "Allocative"], dmunames))
    end

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the first-order derivative of associated Legendre
#   functions.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
# ==============================================================================
#
#   [1] Du, J., Chen, C., Lesur, V., and Wang, L (2015). Non-singular spherical
#       harmonic expressions of geomagnetic vector and gradient tensor fields in
#       the local north-oriented reference frame. Geoscientific Model
#       Development, 8, pp. 1979-1990.
#
#   [2] Ilk, K. H.: Ein Beitrag zur Dynamik ausgedehnter
#       Körper-Gravitationswechselwirkung, Deutsche Geodätische Kommission.
#       Reihe C, Heft Nr. 288, München, 1983.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export dlegendre!, dlegendre
export dlegendre_fully_normalized!, dlegendre_fully_normalized
export dlegendre_schmidt_quasi_normalized!, dlegendre_schmidt_quasi_normalized
export dlegendre_conventional!, dlegendre_conventional

"""
    dlegendre!([N,] dP::AbstractMatrix, ϕ::Number, P::AbstractMatrix, ph_term::Bool = false, n_max::Integer = -1, m_max::Integer = -1)

Compute the first-order derivative of the associated Legendre function
`P_n,m[x]` w.r.t. `ϕ` [rad]:

    dP_n,m[cos(ϕ)]
    --------------
          dϕ

The derivatives will be stored in the matrix `dP`. The maximum degree and order
that will be computed are given by the parameters `n_max` and `m_max`. If they
are negative, then the dimensions of matrix `dP` will be used.

This algorithm needs the matrix `P` with the associated Legendre function. This
can be computed using the function `legendre`. Notice that this matrix must be
computed using the same normalization (see `N`) as the one selected here.

The optional parameter `N` can be used to select the normalization. The
following values are valid:

- `Val(:full)`: Compute the fully normalized associated Legendre function (see
    [`dlegendre_fully_normalized!`](@ref)).
- `Val(:schmidt)`: Compute the Schmidt quasi-normalized associated Legendre
    function (see [`dlegendre_schmidt_quasi_normalized!`](@ref)).
- `Val(:conv)`: Compute the conventional associated Legendre function (see
    [`dlegendre_conventional!`](@ref)).

If `N` is omitted, then the full normalization will be used (`Val(:full)`).

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.
"""
function dlegendre!(
    dP::AbstractMatrix,
    ϕ::Number,
    P::AbstractMatrix,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
)
    return dlegendre_fully_normalized!(dP, ϕ, P, ph_term, n_max, m_max)
end

function dlegendre!(
    ::Val{:full},
    dP::AbstractMatrix,
    ϕ::Number,
    P::AbstractMatrix,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
)
    return dlegendre_fully_normalized!(dP, ϕ, P, ph_term, n_max, m_max)
end

function dlegendre!(
    ::Val{:schmidt},
    dP::AbstractMatrix,
    ϕ::Number,
    P::AbstractMatrix,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
)
    return dlegendre_schmidt_quasi_normalized!(dP, ϕ, P, ph_term, n_max, m_max)
end

function dlegendre!(
    ::Val{:conv},
    dP::AbstractMatrix,
    ϕ::Number,
    P::AbstractMatrix,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
)
    return dlegendre_conventional!(dP, ϕ, P, ph_term, n_max, m_max)
end

"""
    dlegendre([N,] ϕ::Number, n_max::Integer, m_max::Integer = -1, ph_term::Bool = false)

Compute the first-order derivative of the associated Legendre function
`P_n,m[cos(ϕ)]` w.r.t. `ϕ` [rad]:

    dP_n,m[cos(ϕ)]
    --------------
          dϕ

The maximum degree that will be computed is `n_max` and the maximum order is
`m_max`. Notice that if `m_max` is higher than `n_max` or negative, than it is
set to `n_max`.

The optional parameter `N` can be used to select the normalization. The
following values are valid:

- `Val(:full)`: Compute the fully normalized associated Legendre function (see
    [`legendre_fully_normalized`](@ref)).
- `Val(:schmidt)`: Compute the Schmidt quasi-normalized associated Legendre
    function (see [`legendre_schmidt_quasi_normalized`](@ref)).
- `Val(:conv)`: Compute the conventional associated Legendre function (see
    [`dlegendre_conventional!`](@ref)).

If `N` is omitted, then the full normalization will be used (`Val(:full)`).

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Returns

- A matrix with the first-order derivative of the Legendre associated functions
    `P_n,m[cos(ϕ)]`.
- A matrix with the Legendre associated functions `P_n,m[cos(ϕ)]`.
"""
function dlegendre(
    ϕ::Number,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
)
    return dlegendre_fully_normalized(float(ϕ), n_max, m_max, ph_term)
end

function dlegendre(
    ::Val{:full},
    ϕ::Number,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
)
    return dlegendre_fully_normalized(float(ϕ), n_max, m_max, ph_term)
end

function dlegendre(
    ::Val{:schmidt},
    ϕ::Number,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
)
    return dlegendre_schmidt_quasi_normalized(float(ϕ), n_max, m_max, ph_term)
end

function dlegendre(
    ::Val{:conv},
    ϕ::Number,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
)
    return dlegendre_conventional(float(ϕ), n_max, m_max, ph_term)
end

################################################################################
#                Fully Normalized Associated Legendre Functions
################################################################################

"""
    dlegendre_fully_normalized!(dP::AbstractMatrix, ϕ::Number, P::AbstractMatrix, ph_term::Bool = false, n_max::Integer = -1, m_max::Integer = -1)

Compute the first-order derivative of the fully normalized associated Legendre
function `P_n,m[cos(ϕ)]` w.r.t. `ϕ` [rad]:

    dP_n,m[cos(ϕ)]
    --------------
          dϕ

The derivatives will be stored in the matrix `dP`. The maximum degree and order
that will be computed are given by the parameters `n_max` and `m_max`. If they
are negative, then the dimensions of matrix `dP` will be used.

This algorithm needs the matrix `P` with the fully normalized associated
Legendre function. This can be computed using the function
`legendre_fully_normalized`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Remarks

The user is responsible to pass a matrix `P` with the correct values. For
example, if `ph_term` is `true`, then `P` must also be computed with `ph_term`
set to `true`.
"""
function dlegendre_fully_normalized!(
    dP::AbstractMatrix{T1},
    ϕ::Number,
    P::AbstractMatrix{T2},
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
) where {T1<:AbstractFloat, T2<:AbstractFloat}
    T = promote_type(T1, T2)

    # Obtain the maximum degree and order that must be computed.
    n_max, m_max = _get_degree_and_order(dP, P, n_max, m_max)

    # The derivative is compute using the following equation [1, p. 1981]:
    #
    #   dP(n,m)
    #   ------- = a_nm.P(n,m-1) + b_nm.P(n,m+1) ,
    #     dϕ
    #
    #   a_nm = +0.5*sqrt(n+m)*sqrt(n-m+1)*sqrt(C_{m}/C_{m-1})
    #   b_nm = -0.5*sqrt(n+m+1)*sqrt(n-m)*sqrt(C_{m}/C_{m+1})
    #
    #           | 1, m  = 0
    #   C_{m} = |
    #           | 2, m != 0
    #
    # NOTE: The conversion of the coefficients between the full normalization
    # and the Schmidt normalization is performed using:
    #
    #   sqrt(2n+1) ,
    #
    # which depends only on `n`. Since the derivative equation only has terms
    # related to the order `n`, then the same algorithm will work for both full
    # normalization and Schmidt normalization.
    #
    # TODO: This algorithm is based on eq. Z.1.44 of [2]. However, it was
    # verified that it does not provide the correct sign when ϕ ∈ [π, 2π]. This
    # makes sense because the algorithm uses only the values of the
    # coefficients, which are equal for ϕ and -ϕ. However, the derivative w.r.t.
    # ϕ does change. This hack was used so that the values are correct, but
    # further verification is needed.
    #
    # In fact, in [2, p. 119], it is mentioned that `0 <= ϕ <= π`.  However,
    # further clarification is required.

    ϕ    = mod(ϕ, T(2π))
    fact = (ϕ > π) ? -1 : 1

    if ph_term
        fact *= -1
    end

    dP[0+1,0+1] = 0

    m_max < 0 && return nothing

    @inbounds for n in 1:n_max
        for m in 0:n
            if m == 0
                aux  = sqrt(T(n * (n + 1)) / 2)
                a_nm = aux / 2
                b_nm = -a_nm

                # Notice that [1, p. 1985]:
                #
                #                  m
                #   P_(n,-m) = (-1).P_(n,m) ,
                #
                dP[n+1,0+1] = -a_nm * P[n+1,1+1] + b_nm * P[n+1,0+1+1]

            # We should consider the case `m == 1` separately from `n == m`
            # because of the coefficient `C_{m}`.
            elseif m == 1
                a_nm = sqrt(T(2n * (n + 1))) / 2
                dP[n+1,1+1] = a_nm * P[n+1,1-1+1]

                # Only compute `b_nm` if `n > 1`. Otherwise, we could access an
                # invalid memory region if `P` is 2x2.
                if n > 1
                    b_nm = -sqrt(T((n + 2) * (n - 1))) / 2
                    dP[n+1,1+1] += b_nm * P[n+1,1+1+1]
                end

            elseif n != m
                a_nm = +sqrt(T((n + m) * (n - m + 1))) / 2
                b_nm = -sqrt(T((n + m + 1) * (n - m))) / 2

                dP[n+1,m+1] = a_nm * P[n+1,m-1+1] + b_nm * P[n+1,m+1+1]
            else
                a_nm = +sqrt(T((n + m) * (n - m + 1))) / 2

                dP[n+1,m+1] = a_nm * P[n+1,m-1+1]
            end

            dP[n+1,m+1] *= fact

            # Check if the maximum desired order has been reached.
            m >= m_max && break
        end
    end

    return nothing
end

"""
    dlegendre_fully_normalized(ϕ::T, n_max::Integer, m_max::Integer = -1, ph_term::Bool = false) where T<:AbstractFloat

Compute the first-order derivative of the Schmidt fully normalized associated
Legendre function `P_n,m[cos(ϕ)]` w.r.t. `ϕ` [rad]:

    dP_n,m[cos(ϕ)]
    --------------
          dϕ

The maximum degree that will be computed is `n_max` and the maximum order is
`m_max`. Notice that if `m_max` is higher than `n_max` or negative, than it is
set to `n_max`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Returns

A matrix with the first-order derivative of the Legendre associated functions
`P_n,m[cos(ϕ)]`.
"""
function dlegendre_fully_normalized(
    ϕ::T,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
) where T<:AbstractFloat
    if ((m_max < 0) || (m_max > n_max))
        m_max = n_max
    end

    # Check if we need to compute and additional degree in `P` to provide the
    # desire order in `dP`.
    if n_max == m_max
        n_max_P = m_max_P = n_max
    else
        n_max_P = n_max
        m_max_P = m_max + 1
    end

    # First, compute the matrix with the associated Legendre functions.
    P = legendre_fully_normalized(ϕ, n_max_P, m_max_P, ph_term)

    # Now, compute and return the time-derivative of the associated Legendre
    # functions.
    dP = zeros(T, n_max + 1, m_max + 1)
    dlegendre_fully_normalized!(dP, ϕ, P, ph_term)

    return dP, P
end

################################################################################
#            Schmidt Quasi-Normalized Associated Legendre Functions
################################################################################

"""
    dlegendre_schmidt_quasi_normalized!(dP::AbstractMatrix, ϕ::Number, P::AbstractMatrix, ph_term::Bool = false, n_max::Integer = -1, m_max::Integer = -1)

Compute the first-order derivative of the Schmidt quasi-normalized associated
Legendre function `P_n,m[cos(ϕ)]` w.r.t. `ϕ` [rad]:

    dP_n,m[cos(ϕ)]
    --------------
          dϕ

The derivatives will be stored in the matrix `dP`. The maximum degree and order
that will be computed are given by the parameters `n_max` and `m_max`. If they
are negative, then the dimensions of matrix `dP` will be used.

This algorithm needs the matrix `P` with the Schmidt quasi-normalized associated
Legendre function. This can be computed using the function
`legendre_schmidt_quasi_normalized`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Remarks

The user is responsible to pass a matrix `P` with the correct values. For
example, if `ph_term` is `true`, then `P` must also be computed with `ph_term`
set to `true`.
"""
function dlegendre_schmidt_quasi_normalized!(
    dP::AbstractMatrix,
    ϕ::Number,
    P::AbstractMatrix,
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
)
    # The algorithm to compute the first-order derivative using Schmidt
    # normalization is precisely the same as the one that computes using full
    # normalization.
    return dlegendre_fully_normalized!(dP, ϕ, P, ph_term, n_max, m_max)
end

"""
    dlegendre_schmidt_quasi_normalized(ϕ::T, n_max::Integer, m_max::Integer = -1, ph_term::Bool = false) where T<:AbstractFloat

Compute the first-order derivative of the Schmidt quasi-normalized associated
Legendre function `P_n,m[cos(ϕ)]` w.r.t. `ϕ` [rad]:

    dP_n,m[cos(ϕ)]
    --------------
          dϕ

The maximum degree that will be computed is `n_max` and the maximum order is
`m_max`. Notice that if `m_max` is higher than `n_max` or negative, than it is
set to `n_max`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Returns

A matrix with the first-order derivative of the Legendre associated functions
`P_n,m[cos(ϕ)]`.

"""
function dlegendre_schmidt_quasi_normalized(
    ϕ::T,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
) where T<:AbstractFloat
    if ((m_max < 0) || (m_max > n_max))
        m_max = n_max
    end

    # Check if we need to compute and additional degree in `P` to provide the
    # desire order in `dP`.
    if n_max == m_max
        n_max_P = m_max_P = n_max
    else
        n_max_P = n_max
        m_max_P = m_max + 1
    end

    # First, compute the matrix with the associated Legendre functions.
    P = legendre_schmidt_quasi_normalized(ϕ, n_max_P, m_max_P, ph_term)

    # Now, compute and return the time-derivative of the associated Legendre
    # functions.
    dP = zeros(T, n_max + 1, m_max + 1)
    dlegendre_schmidt_quasi_normalized!(dP, ϕ, P, ph_term)

    return dP, P
end

################################################################################
#                  Conventional Associated Legendre Function
################################################################################

"""
    dlegendre_conventional!(dP::AbstractMatrix, ϕ::Number, P::AbstractMatrix, ph_term::Bool = false, n_max::Integer = -1, m_max::Integer = -1)

Compute the first-order derivative of the conventional associated Legendre
function `P_n,m[cos(ϕ)]` w.r.t. `ϕ` [rad]:

    dP_n,m[cos(ϕ)]
    --------------
          dϕ

The derivatives will be stored in the matrix `dP`. The maximum degree and order
that will be computed are given by the parameters `n_max` and `m_max`. If they
are negative, then the dimensions of matrix `dP` will be used.

This algorithm needs the matrix `P` with the conventional associated Legendre
function. This can be computed using the function `legendre_conventional`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Remarks

The user is responsible to pass a matrix `P` with the correct values. For
example, if `ph_term` is `true`, then `P` must also be computed with `ph_term`
set to `true`.
"""
function dlegendre_conventional!(
    dP::AbstractMatrix{T1},
    ϕ::Number,
    P::AbstractMatrix{T2},
    ph_term::Bool = false,
    n_max::Integer = -1,
    m_max::Integer = -1
) where {T1<:AbstractFloat, T2<:AbstractFloat}
    T = promote_type(T1, T2)

    # Obtain the maximum degree and order that must be computed.
    n_max, m_max = _get_degree_and_order(dP, P, n_max, m_max)

    # The derivative is computed using the following equation [1, p. 1981]:
    #
    #   dP(n,m)
    #   ------- = 0.5.( (n+m)(n-m+1)P(n,m-1) - P(n,m+1) ) ,
    #     dθ
    #
    # TODO: This algorithm is based on eq. Z.1.44 of [2]. However, it was
    # verified that it does not provide the correct sign when ϕ ∈ [π, 2π]. This
    # makes sense because the algorithm uses only the values of the
    # coefficients, which are equal for ϕ and -ϕ. However, the derivative w.r.t.
    # ϕ does change. This hack was used so that the values are correct, but
    # further verification is needed.
    #
    # In fact, in [2, p. 119], it is mentioned that `0 <= ϕ <= π`.  However,
    # further clarification is required.

    ϕ    = mod(ϕ, T(2π))
    fact = (ϕ > π) ? -1 : 1

    ph_term && (fact *= -1)

    dP[0+1,0+1] = 0

    m_max < 0 && return nothing

    @inbounds for n = 1:n_max
        for m = 0:n
            if m == 0
                dP[n+1,m+1] = -P[n+1,1+1]
            elseif n != m
                dP[n+1,m+1] = ((n + m) * (n - m + 1) * P[n+1,m-1+1] - P[n+1,m+1+1]) / 2
            else
                dP[n+1,m+1] = ((n + m) * (n - m + 1) * P[n+1,m-1+1]) / 2
            end

            dP[n+1,m+1] *= fact

            # Check if the maximum desired order has been reached.
            m >= m_max && break
        end
    end

    nothing
end

"""
    dlegendre_conventional(ϕ::Number, n_max::Integer, m_max::Integer = -1, ph_term::Bool = false)

Compute the first-order derivative of the conventional associated Legendre
function `P_n,m[cos(ϕ)]` w.r.t. `ϕ` [rad]:

    dP_n,m[cos(ϕ)]
    --------------
          dϕ

The maximum degree that will be computed is `n_max` and the maximum order is
`m_max`. Notice that if `m_max` is higher than `n_max` or negative, than it is
set to `n_max`.

If `ph_term` is set to `true`, then the Condon-Shortley phase term `(-1)ᵐ` will
be included. If `ph_term` is not present, then it defaults to `false`.

# Returns

A matrix with the first-order derivative of the Legendre associated functions
`P_n,m[cos(ϕ)]`.
"""
function dlegendre_conventional(
    ϕ::T,
    n_max::Integer,
    m_max::Integer = -1,
    ph_term::Bool = false
) where T<:AbstractFloat
    if ((m_max < 0) || (m_max > n_max))
        m_max = n_max
    end

    # Check if we need to compute and additional degree in `P` to provide the
    # desire order in `dP`.
    if n_max == m_max
        n_max_P = m_max_P = n_max
    else
        n_max_P = n_max
        m_max_P = m_max + 1
    end

    # First, compute the matrix with the associated Legendre functions.
    P = legendre_conventional(ϕ, n_max_P, m_max_P, ph_term)

    # Now, compute and return the time-derivative of the associated Legendre
    # functions.
    dP = zeros(T, n_max + 1, m_max + 1)
    dlegendre_conventional!(dP, ϕ, P, ph_term)

    return dP, P
end

################################################################################
#                                   Private
################################################################################

"""
    _get_degree_and_order(dP, P, n_max, m_max)

Return the maximum degree and order to compute the Legendre associated functions
given the matrices `dP`, `P`, and the configuration values `n_max` and `m_max`.

"""
@inline function _get_degree_and_order(dP, P, n_max, m_max)
    # Get the size of the matrices.
    Prows,  Pcols  = size(P)
    dProws, dPcols = size(dP)

    rows = min(Prows, dProws)
    cols = min(Pcols, dPcols)

    # If the order or degree is less than 0, then the user wants to use all the
    # available memory.
    if n_max < 0
        n_max = rows - 1
    end

    if m_max < 0
        m_max = (cols <= rows) ? cols - 1 : n_max
    end

    # Make sure that the degree and order fits the matrix.
    if n_max > rows - 1
        n_max = rows - 1
    end

    if ((m_max > cols - 1) || (m_max > n_max))
        m_max = min(cols - 1, n_max)
    end

    return n_max, m_max
end


text:
module Day3
export day3

function parseInput(inputfile="inputs/input03.txt")
    lines = collect(eachline(inputfile))
    data = BitArray{2}(undef, length(lines), length(lines[1]))
    for (i, line) in enumerate(lines)
        data[i,:] .= collect(line) .== '#'
    end
    data
end

function day3(inputfile="inputs/input03.txt")
    parsedInput = parseInput(inputfile)
    slopes = ((1,1), (3,1), (5,1), (7,1), (1,2))
    (countTrees(parsedInput, 3, 1), reduce(*, map(s->countTrees(parsedInput, s...), slopes)))
    #(count(isValidPWD, parsedInput), count(isValidPWDPart2, parsedInput))
end

function countTrees(data, right, down)
    trees = 0
    line = 1
    position = 0
    while line <= size(data, 1)
        if data[line, position+1]
            trees += 1
        end
        line += down
        position = mod(position+right, size(data, 2))
    end
    trees
end

end #module